水声匹配场定位原理及实验

qianpeng897
0 阅读5分钟

1. 原理

匹配场定位是通过数值模拟海洋声场，进行拷贝场计算，并通过拷贝场和测量场相关匹配估计目标距离深度的一种方法。匹配场的可看作为非平面条件下的波束形成，核心思路在“匹配”二字，在海洋环境无失配条件下，匹配场的定位效果较好。 本文代码开源，地址 UWSL-MFP: 水下声源定位的匹配场算法

image.png

在海洋中的接收阵接收频域信号可以表示为

X(ω)⏟接收信号=ps(ω)S(ω)⏟声场信号+N(ω)⏟阵噪声=[p1(ω)⋮pM(ω)]⏟声场向量S(ω)⏟源信号+[n1(ω)⋮nM(ω)]⏟噪声向量\underbrace{\boldsymbol{X}(\omega)}_{\text{接收信号}} = \underbrace{\boldsymbol{p}_s(\omega)S(\omega)}_{\text{声场信号}} + \underbrace{\boldsymbol{N}(\omega)}_{\text{阵噪声}} = \underbrace{\begin{bmatrix}p_1(\omega)\\\vdots\\p_M(\omega)\end{bmatrix}}_{\text{声场向量}} \underbrace{S(\omega)}_{\text{源信号}} + \underbrace{\begin{bmatrix}n_1(\omega)\\\vdots\\n_M(\omega)\end{bmatrix}}_{\text{噪声向量}}

这里使用频域表达式的原因是表达式更加简洁，在实际匹配场中，也是先对接收时域信号做傅里叶变换转为频域后，再做匹配场处理。

对接收信号求协方差矩阵，得到

RX=E{X(ω)XH(ω)}=ps(ω)σS2psH(ω)+RN\mathbf{R}_X = \mathbb{E}\left\{\boldsymbol{X}(\omega)\boldsymbol{X}^\mathrm{H}(\omega)\right\} = \boldsymbol{p}_s(\omega)\sigma_S^2\boldsymbol{p}_s^\mathrm{H}(\omega) + \boldsymbol{R}_N

其中：

  • RX\mathbf{R}_X阵列接收信号：协方差矩阵
  • σS2=E{∣S(ω)∣2}\sigma_S^2 = \mathbb{E}\{|S(\omega)|^2\}：源信号功率
  • RN=E{N(ω)NH(ω)}\mathbf{R}_N = \mathbb{E}\{\boldsymbol{N}(\omega)\boldsymbol{N}^\mathrm{H}(\omega)\}：噪声协方差矩阵

默认假设：信号与噪声不相关。

一个简单的设想是，通过假定声源位置，计算出声场向量，并与协方差矩阵做相关处理，如果假设位置刚好为声源位置的话，使相关峰值最大，这样就得到了估计的声源位置。

按照这样的设想，设计权重向量 W\boldsymbol{W} 。设 V\boldsymbol{V} 为拷贝场向量，设定归一化约束条件

WH(ϕ)V(ϕ)=1\boldsymbol{W}^{H}\left(\phi\right)\boldsymbol{V}\left(\phi\right)=1

在匹配场定位中，模糊平面 PP 由互谱密度矩阵 K\mathbf{K} 和权向量表示为

P(ϕ)=WH(ϕ)RXW(ϕ)P(\phi)=\boldsymbol{W}^H(\phi)\mathbf{R}_X\boldsymbol{W}(\phi)

式中，ϕ\phi 为声源位置参数，需根据定位精度预先划分搜索网格，精度要求越高，网格越密集。对每个网格点假设声源位于该处，通过匹配场处理得到各点输出功率，最终模糊平面的峰值位置即为声源估计位置。

显然，最简单且直观的权重向量就是拷贝场向量本身（向量自相关性最强），这种处理器称为 Bartlett 处理器，又叫常规处理器、线性处理器。

WBart(ϕ)=V(ϕ)VH(ϕ)V(ϕ)\boldsymbol{W}_{Bart}(\phi)=\frac{\boldsymbol{V}(\phi)}{\boldsymbol{V}^{H}(\phi)\boldsymbol{V}(\phi)}

MVDR 匹配处理器基于最小均方误差准则，在约束目标方向增益不变的同时使总输出能量最小，因此可以减小 Bartlett 处理器的旁瓣，抑制干扰，增强失配稳健性。 其核心是在约束条件下最小化输出能量：

min⁡WH(ϕ)RXW(ϕ)s.t.WH(ϕ)V(ϕ)=1\min \boldsymbol{W}^H(\phi) \mathbf{R}_X \boldsymbol{W}(\phi) \quad \text{s.t.} \quad \boldsymbol{W}^H(\phi) \boldsymbol{V}(\phi) = 1

通过拉格朗日乘子法，可求得权向量：

WMV(ϕ)=RX−1V(ϕ)VH(ϕ)RX−1V(ϕ)\boldsymbol{W}_{MV}(\phi) = \frac{\mathbf{R}_X^{-1} \boldsymbol{V}(\phi)}{\boldsymbol{V}^H(\phi) \mathbf{R}_X^{-1} \boldsymbol{V}(\phi)}

代入后得到输出模糊度平面：

PMV(ϕ)=1VH(ϕ)RX−1V(ϕ)P_{MV}(\phi) = \frac{1}{\boldsymbol{V}^H(\phi) \mathbf{R}_X^{-1} \boldsymbol{V}(\phi)}

2. 仿真实现

为了更直观的展现匹配场的算法逻辑，下面基于简正波模型仿真实现单频的频域匹配场算法。

简正波模型的声场计算简单，声场表示为：

p(r,z)=iρ(zs)8πre−iπ/4∑m=1∞Ψm(zs)Ψm(z)eikrmrkrmp(r,z) = \frac{i}{\rho(z_s)\sqrt{8\pi r}} e^{-i\pi/4} \sum\limits_{m=1}^{\infty} \Psi_m(z_s) \Psi_m(z) \frac{e^{ik_{rm}r}}{\sqrt{k_{rm}}}

式里各变量的物理含义如下：

  • rr：水平距离，一般指声源到接收点的水平距离（单位：m）。
  • zsz_s：声源深度，即声源在水中的垂直位置（单位：m）。
  • zz：接收点深度，即水听器/接收器在水中的垂直位置（单位：m）。
  • krmk_{rm}：第 mm 阶简正波的水平波数（本征值），描述该阶模式在水平方向上的传播特性（单位：1/m）。
  • Ψm(z)\Psi_m(z)：第 mm 阶简正波的深度本征函数，描述该阶模式在深度 zz 处的振幅分布。
  • ρ(zs)\rho(z_s)：声源深度 zsz_s 处的海水密度（单位：kg/m³），在浅海近似中常取常数。
  • ii：虚数单位，满足 i2=−1i^2 = -1。
  • p(r,z)p(r,z)：在深度 zz、水平距离 rr 处的声压，是声场的核心物理量。

在典型的Pekeris浅海环境下进行匹配场定位，Pekeris 波导为典型的浅海波导环境，其海洋环境参数下图所示。海洋声场模型分为真空层上半空间、海水层和海底层下半空间，其中海水层为声速 cw=1500 m/sc_w = 1500\ \text{m/s} 、密度 ρw=1000 kg/m3\rho_w = 1000\ \text{kg/m}^3 、海深 D=200 mD = 200\ \text{m} 的均匀介质；海底声速 cb=1600 m/sc_b = 1600\ \text{m/s} ，海底密度 ρb=1800 kg/m3\rho_b = 1800\ \text{kg/m}^3 ，海底衰减 αb=0.2 dB/λ\alpha_b = 0.2\ \text{dB/}\lambda 。声源深度 zs=25 mz_s = 25\ \text{m} ，接收器深度 zr=1∼91 mz_r = 1 \sim 91\ \text{m} ，间隔5m，阵元个数为19个，接收器距离 rr=2 kmr_r = 2\ \text{km} 。

image.png

首先进行环境设定

% 系统参数
f = 100;              % 声源频率
fs = 3e4;             % 采样频率
rho = 1.03;           % 海水密度

% 声源参数
source_depth = 25;    % 实际声源深度
source_range = 2000;  % 实际声源距离
w = 2*pi*f;           % 角频率

% 接收基阵参数
array_start_depth = 1;    % 接收水听器基阵的第一个阵元深度
array_end_depth = 91;     % 接收水听器基阵的最后一个阵元深度
array_element_spacing = 5; % 阵元间隔
num_elements = 19;        % 阵元数目
SNR = 5;                  % 信噪比

% 扫描参数
scan_range_start = 100;   % 扫描距离的近点
scan_range_end = 10000;   % 扫描距离的远点
scan_range_step = 10;
scan_depth_start = 1;     % 扫描深度的初始点
scan_depth_end = 99;      % 扫描深度的最深点
scan_depth_step = 1;      % 扫描深度的精度

% 时间向量
T = 1;                    % 信号持续时间
t = 0:1/fs:(T-1/fs);
fk = (0:length(t)-1)/length(t)*fs;

% 匹配场类型
type = 'Bart'; % 'MV': MVDR 匹配处理器  'Bart': Bartlett 处理器

env 设定时，将接收器1m一个放置，后续计算声场会比较方便，直接把接收器深度（整数）当作索引带入到代码进行计算即可

'pekeris'		! 标题
100			! 频率
1			! NMEDIA，介质层数
'CVW'		! TOP OPTION，海面选项。第一个字母代表声速剖面的插值方法，第二个字母代表上半边界条件，第三个字母代表衰减系数单位。
0  0.0  200.0		! 海水竖直网格层个数，界面粗糙度，海水深度
    0.0  1500.00  0.0  1.0  0.0  0.0 	! 深度，纵波速度，横波速度，密度，纵波衰减系数、横波衰减系数
    200.0 1500.00  /		! 深度，纵波速度，'/'代表后面几项与上面保持一致
'A' 0.0              	! BOTTOM OPTION，海底选项，A：海底半空间条件，0.0海底粗糙度
    200.0  1600.00 0.0 1.8 0.2 0.0 / ! 海底深度，纵波速度度，横波速度，密度，纵波衰减，横波衰减
0 20000		! CLOW，CHIGH，相速度计算的范围
0			! RMAX (km) 
1			! NSD 规定声源个数
25/			! SD(1:NSD) (m) 规定声源深度
200			! NRD  规定接收器深度取值个数
1 200/		! RD(1:NRD) (m) 规定接收器深度

直接调用kraken.exe生成.mod，使用.mod中的简正波模态值和模态函数信息，无需计算.shd声场

%% 调用kraken程序
!kraken.exe pekerisK
clear read_modes_bin
[modes] = read_modes_bin_once('pekerisK.mod', 0);
horizontal_wavenumber = modes.k.';  % 简正波模态值
mode_matrix = modes.phi;          % 简正波模态函数
num_modes = modes.M;              % 简正波数量

匹配场计算核心代码如下，在计算signal时，先将其转为了频域，取100Hz中心3个频率点作为信号，计算接收信号协方差矩阵。

%% 匹配场计算
measured_data = mes_vector(...
    horizontal_wavenumber, mode_matrix, num_modes, num_elements, ...
    source_depth, source_range, array_start_depth, array_end_depth, ...
    array_element_spacing, rho);

% 生成多采样数据
sampled_data = zeros(num_elements, length(t));
for n = 1:num_elements
    sampled_data(n,:) = measured_data(n) * exp(1i*w*t);
end

% 添加噪声
signal_with_noise = awgn(sampled_data, SNR, 'measured'); % 信噪比5dB

% 声源做傅里叶变换
source_fourier = fft(sampled_data, [], 2);
[~, index] = min(abs(fk - f));
% 取中心3个频率点作为信号
index_l = min(index-3, 1);
index_r = max(index+3, length(t));
source_fourier_v = source_fourier(:, index_l:index_r);

% 计算网格点处的拷贝场
copy_field = rep_vector(...
    horizontal_wavenumber, mode_matrix, num_elements, scan_range_start, ...
    scan_range_end, scan_range_step, num_modes, scan_depth_start, scan_depth_end, ...
    scan_depth_step, rho, array_start_depth, array_end_depth, ...
    array_element_spacing);

% 计算模糊度平面
ambiguity_surface = calculate_ambiguity_surface(copy_field, source_fourier_v, num_elements, type);

% 生成深度和距离向量
depth_vector = scan_depth_start:scan_depth_step:scan_depth_end;
range_vector = scan_range_start:scan_range_step:scan_range_end;

首先计算测量场，直接使用简正波理论，将简正波本征值和本征函数代入公式计算。代码中horizontal_wavenumber代表m号本征值krmk_{rm}，mode_matrix代表本征函数Ψm(z)\Psi_m(z)，计算1到num_modes阶简正波的声压贡献，直接通过矩阵点乘得到所有阵元上的接收声场。

% 根据点源声场的远场解表达式算出水听器的实测数据
function measured_data = mes_vector(...
    horizontal_wavenumber, mode_matrix, num_modes, num_elements, ...
    source_depth, source_range, array_start_depth, array_end_depth, ...
    array_element_spacing, rho)
% horizontal_wavenumber - 水平波数
% mode_matrix - 深度分布函数
% num_modes - 所取简正波号数
% num_elements - 阵元数目
% source_depth - 实际声源深度
% source_range - 实际声源距离
% array_start_depth - 实际接收水听器基阵的第一个阵元深度
% array_end_depth - 实际接收水听器基阵的最后一个阵元深度
% array_element_spacing - 实际的阵元间隔
% rho - 海水密度


% 初始化变量
measured_data = zeros(num_elements,1);
element_index = 1;

% 计算每个阵元的接收信号
for depth = array_start_depth:array_element_spacing:array_end_depth
    % 简正波叠加和，存入数组
    measured_data(element_index) = -1i * sum(sqrt(2*pi./(horizontal_wavenumber(1:num_modes)*source_range*rho)) .* ...
            mode_matrix(depth, 1:num_modes) .* mode_matrix(source_depth, 1:num_modes) .* ...
            exp(-1i*(horizontal_wavenumber(1:num_modes)*source_range - pi/4)));
    element_index = element_index + 1;
end

拷贝场计算通过扫描声源网格，将简正波本征值和本征函数代入公式计算。

% 网格点处的拷贝场
% 根据点源声场的远场解表达式算出网格点处的拷贝场
% 最后输出的copy_field是每个网格点处的拷贝场
function copy_field = rep_vector(...
    horizontal_wavenumber, mode_matrix, num_elements, scan_range_start, scan_range_end, scan_range_step, ...
    num_modes, scan_depth_start, scan_depth_end, scan_depth_step, ...
    rho, array_start_depth, array_end_depth, array_element_spacing)
% horizontal_wavenumber - 水平波数
% mode_matrix - 简正波深度分布矩阵
% num_elements - 阵元数目
% scan_range_start - 扫描距离的近点
% scan_range_end - 扫描距离的远点
% num_modes - 所取简正波号数
% scan_depth_start - 扫描深度的初始点
% scan_depth_end - 扫描距离的最深点
% scan_depth_step - 扫描深度的精度
% rho - 海水密度（本程序未使用）
% array_start_depth - 实际接收水听器基阵的第一个阵元深度
% array_end_depth - 实际接收水听器基阵的最后一个阵元深度
% array_element_spacing - 实际的阵元间隔

% 计算网格大小
num_range_steps = (scan_range_end - scan_range_start)/scan_range_step + 1;
num_depth_steps = (scan_depth_end - scan_depth_start)/scan_depth_step + 1;

% 初始化拷贝场矩阵
copy_field = zeros(num_elements, num_range_steps, num_depth_steps);

depth_index = 1;

% 遍历每个扫描深度
for scan_depth = scan_depth_start:scan_depth_step:scan_depth_end
    range_index = 1;
    
    % 遍历每个扫描距离
    for scan_range = scan_range_start:scan_range_step:scan_range_end
        element_index = 1;
        element_pressure = zeros(num_elements,1);
        
        % 遍历每个阵元深度
        for element_depth = array_start_depth:array_element_spacing:array_end_depth           
            % 简正波叠加和，存入数组
            element_pressure(element_index) = -1i * sum(sqrt(2*pi./(horizontal_wavenumber(1:num_modes)*scan_range*rho)) .* ...
            mode_matrix(element_depth, 1:num_modes) .* mode_matrix(scan_depth, 1:num_modes) .* ...
            exp(-1i*(horizontal_wavenumber(1:num_modes)*scan_range - pi/4)));
            element_index = element_index + 1;
        end
        
        % 存储当前距离和深度的拷贝场
        copy_field(:, range_index, depth_index) = element_pressure;
        range_index = range_index + 1;
    end
    
    depth_index = depth_index + 1;
end
end

最后通过相关处理，得到模糊平面。注意MV处理器中，协方差矩阵进行了对角加载，防止协方差矩阵不满秩导致求逆的不稳定。

% 对拷贝场和实测数据进行相关处理，得到模糊度平面
function ambiguity_surface = calculate_ambiguity_surface(copy_field, signal, num_elements, type)
% angular_frequency - 2*pi*f角频率
% copy_field - 每个网格点处的拷贝场
% time_vector - 时间向量
% signal - 信号
% num_elements - 阵元数目


% 计算信号协方差矩阵
covariance_matrix = signal * signal';

ambiguity_surface = zeros(size(copy_field, 2), size(copy_field, 3));

% 遍历每个深度
for depth_idx = 1:size(copy_field, 3)
    % 遍历每个距离
    for range_idx = 1:size(copy_field, 2)
        % 提取某一距离和深度的拷贝场
        copy_vector = copy_field(:, range_idx, depth_idx);

        % 归一化拷贝向量
        copy_vector = copy_vector / sqrt(copy_vector' * copy_vector);

        % 计算最小方差处理器输出
        if strcmp(type, 'MV')
            denominator = copy_vector' / (covariance_matrix + 0.0001 * sum(diag(covariance_matrix)+1e-9) * eye(num_elements)) * copy_vector;
            ambiguity_value = 1 / denominator;
        elseif strcmp(type, 'Bart')
            ambiguity_value = copy_vector' * covariance_matrix * copy_vector;
        end
        % 存储结果
        ambiguity_surface(range_idx, depth_idx) = ambiguity_value;
    end
end
end

绘图和寻找最大值位置程序如下

% 绘制模糊度平面
figure();
pcolor(range_vector, depth_vector, abs(ambiguity_surface)');
set(gca, 'ydir', 'reverse');
shading interp;
xlabel('水平距离/m');
ylabel('深度/m');

% 寻找最大值位置
[max_values, depth_indices] = max(ambiguity_surface');
[~, range_index] = max(max_values);
depth_index = depth_indices(range_index);

% 显示结果
disp(['处理结果：深度 ', num2str(depth_vector(depth_index)), 'm   距离 ', num2str(range_vector(range_index)), 'm;']);
toc

Bartlet处理器模糊度图结果如下，处理结果：深度 25m 距离 2000m，估计结果准确。

veticalArrayResultBart.png

MVDR处理器处理结果：深度 25m 距离 2000m，其模糊度图比较尖锐。

veticalArrayResult.png

avatar
文章
阅读
粉丝