0 引　言

水下目标声源定位是水声学研究的重点和热点问题，匹配场定位是目前水下声源定位的主要手段之一。首先将声源可能存在的区域在深度和距离上网格化处理，然后以网格点为假想声源位置，利用声场传播模型计算拷贝场，最后利用实际声场与拷贝场的相关性估计声源位置。

1976年Bucker^[1]首次提出了匹配场定位的概念，他从Hinich的结论推导出线性处理器，但这种处理器对多普勒失配非常敏感，在采用长时间积分时更明显。后来Bucker便构造了一种称为检测因子的估计器，它可以有效抑制这种敏感，使用真实环境模型，并引入模糊度平面的概念^[2]。在此之前，所有研究都是基于垂直阵进行匹配场处理。1977年，Hinich^[3]使用水平阵对声源定深，并得出结论：水平列阵比垂直列阵能够更大程度的利用费雪信息来估计平稳均匀水平波导的源深度。到了1981年，Klemmn等^[4]又提出了一种高分辨的最大熵波束形成器，他从列阵处理器的角度进行研究，发现这种高分辨波束形成器能很大程度提高匹配场定位性能。紧接着，Fizell^[5]和Baggeroer^[6]等系统地研究用高分辨的波束形成器对匹配场处理器的性能进行改进，使用这种思想的典型代表就是由Capon^[7]提出的最小方差无失真响应波束形成器。1984年，尚尔昌^[8]采用模态滤波器将垂直列阵的场采样数据转化为由各模态特征函数和模态相位实现了声源深度估计，使匹配场定位产生了一个新的研究方向。1987年，Yang^[9]使用一种特征向量分解技术提取了一个长1 km垂直列阵接收到的数据的模态振幅，对模拟数据和北冰洋上实验数据进行了成功的估计。1990年，Yang在文献[9]的基础上指出，匹配域处理相当于一种不同模式分解算法下的模态波束形成方法。实验表明采用模态滤波可抑制上部声速不匹配引起的侧偏^[10]。2000年，Collison和Dosso^[11]将测量数据分解得到了传播模式的激励，然后将这些激励与网格计算的模型复制激励相匹配实现定位，该算法可以适应环境的不匹配。同年，Orris和Nicholas等^[12]采用匹配相位相干匹配场处理器研究了相干多频匹配场处理，这种处理器在声源水平方向的灵敏度较高，表现出更优越的距离分辨力。2003年，Soares和Jesus^[13]提出了一个由确定的信道乘以一个扰动随机因子加上噪声组成观测信号的一致性模型，该模型的模糊面最大时具有与匹配相位处理器相同的性能，无需估计相位项，有很低的计算代价。2004年，杨坤德等^[14]提出了一种在不确定环境和强干扰背景下检测水下弱目标的稳健自适应匹配场方案，该方案可以在环境失配条件下有效抑制水面强干扰，还可以为水下微弱目标检测提供较大范围的空时相干积累，有利于提高目标的定位精度和输出信干比。2005年，Touze等^[15]利用实际模型对模态振幅进行时间-频域表征，声源深度估计令人满意，同时对噪声也有很好的鲁棒性。2006年，Nicolas等^[16]利用水平列阵提出一种利用传感器HLA记录压力场的模态分解估计声源深度的有效方法，结果表明通过合并模态振幅，源深度估计得到了很大改善。2007年，Debever和Kuperman^[17]使用相干宽带白噪声约束处理器进行窄带处理，结果表明该处理器可以增强对环境不匹配和快照缺陷的鲁棒性。

本文将基于矩阵空间特征分解的算法应用于水下声源定位，将信号处理的观测空间分解为信号子空间和噪声子空间，构造网格化的空间谱实现声源定位。首先简要阐述基本声学理论，给出了2种传统匹配场处理算法的相关结论，推导给出基于矩阵空间特征分解的匹配场算法的定位表达式。经过单声源、双声源和环境失配条件下的仿真，比较了线性处理器、最小方差处理器和基于矩阵空间特征分解处理器的性能。最后通过实验数据进行验证，结果表明基于特征分解的处理器可以实现声源准确定位，尤其在双声源定位时较最小方差处理器有更大优势。

1 匹配场处理的原理 1.1 硬底均匀浅海声场模型

根据简正波理论，在海洋分层波导环境下，声场将表示为简正波展开或叠加的形式。在柱坐标系下，位于 $ {z_s} $ 深处的声源在 $ \left( {r,z} \right) $ 处的远场声场表示为：

$ p\left( {r,z} \right) = \frac{i}{{\rho \left( {{z_s}} \right)\sqrt {8\text{π} r} }}{e^{ - i\frac{\text{π} }{4}}}\sum\limits_{n = 1}^\infty {{Z_n}\left( {{z_s}} \right)} {Z_n}\left( z \right)\frac{{{e^{i{k_n}r}}}}{{\sqrt {{k_n}} }}\text{。} $

(1)

其中： $ i $ 为虚数， $ r $ 为接收点水平距离， $ z $ 为接收点深度， $\rho \left( {{z_s}} \right)$ 为海水密度， $ {z_s} $ 为声源深度， $ s = {\rm{1,2}}, \cdots ,N $ ， $ {Z_n} $ 和 $ {k_n} $ 为第 $ n $ 阶简正波的波数和正交归一化的本征函数。

1.2 Bartlett处理器

对于 $ M $ 元垂直阵接收到的数据在频域表示为：

$ {\boldsymbol{x}}\left( {{\omega _n}} \right) = {\boldsymbol{A}}\left( {{\omega _n},{r_s},{z_s},\ell } \right){\boldsymbol{s}}\left( {{\omega _n}} \right) + {\boldsymbol{n}}\left( {{\omega _n}} \right)\text{。} $

(2)

其中： $ {\boldsymbol{A}}\left( {{\omega _n},{r_s},{z_s},\ell } \right) $ 为水听器在频率 $ {\omega _n} $ 处采集到的信号强度， ${\boldsymbol{s}}\left( {{\omega _n}} \right)$ 为声源在点 $ \left( {{r_s},{z_s}} \right) $ 处环境参数为 $ \ell $ 时传播模型计算得到的拷贝场向量， ${\boldsymbol{n}}\left( {{\omega _n}} \right)$ 为噪声向量。

得到测量数据之后，再将声源可能存在的区域进行网格化处理，即将指定区域作深度和距离上的间隔划分得到一个网格区域，然后通过传播模型计算每个网格点处的拷贝场，最后把测量数据与所有网格点的拷贝场作相关匹配。

Bartlett处理器的模糊度表达式可表示为：

$ {B_{B{\rm{artlett}}}} = {\boldsymbol{p}}_{\boldsymbol{c}}^H{\boldsymbol{R}}{{\boldsymbol{p}}_{\boldsymbol{c}}} \text{。}$

(3)

其中： $ {\boldsymbol{R}} = \displaystyle\frac{1}{L}\sum\limits_{l = 1}^L {{{\boldsymbol{x}}_{\boldsymbol{l}}}\left( {{\omega _n}} \right)} {\boldsymbol{x}}_{\boldsymbol{l}}^H\left( {{\omega _n}} \right) $ ，表示在一个时间周期内采集了 $ L $ 组统计独立的数据 $\left[ {{{\boldsymbol{x}}_{\boldsymbol{l}}}\left( {{\omega _n}} \right)} \right]_{l = 1}^L$ ； $ {{\boldsymbol{p}}_{\boldsymbol{c}}} $ 为传播模型计算得到的网格点处拷贝场。

1.3 最小方差处理器

自适应匹配场的功率输出为：

$ {B_{MV}} = {{\boldsymbol{w}}^H}\left( {{\omega _n}} \right){\boldsymbol{Rw}}\left( {{\omega _n}} \right) ,$

(4)

权向量 $ {\boldsymbol{w}}\left( {{\omega _n}} \right) $ 通过拷贝场向量 $ {{\boldsymbol{p}}_{\boldsymbol{c}}} $ 进行归一化：

$ {{\boldsymbol{w}}^H}\left( {{\omega _n}} \right){{\boldsymbol{p}}_{\boldsymbol{c}}} = 1 ,$

(5)

由此得到最小方差处理器的自适应权向量：

$ {\boldsymbol{w}}\left( {{\omega _n}} \right) = \frac{{{{\boldsymbol{R}}^{ - 1}}{{\boldsymbol{p}}_{\boldsymbol{c}}}}}{{{\boldsymbol{p}}_{\boldsymbol{c}}^H{{\boldsymbol{R}}^{ - 1}}{{\boldsymbol{p}}_{\boldsymbol{c}}}}} ,$

(6)

将上式代入式（4）中，得到简化的功率输出：

$ {B_{MV}} = \frac{1}{{{\boldsymbol{p}}_{\boldsymbol{c}}^H{{\boldsymbol{R}}^{ - 1}}{{\boldsymbol{p}}_{\boldsymbol{c}}}}} \text{。}$

(7)

1.4 基于特征分解的匹配场算法

特征分解是线性代数中的一种矩阵分解方式，本文将矩阵空间特征分解引入到匹配场声源定位中。依据特征分解的有关知识，将阵元接收数据的协方差矩阵进行特征分解求出特征值和特征向量，然后可以根据特征值的大小找出噪声向量对应的特征值，再根据特征值找到噪声向量构成噪声子空间，利用噪声子空间和信号子空间相互正交这一关系将噪声子空间与每一个网格点处的拷贝场向量作相关，得到与噪声子空间最为接近正交的拷贝场向量坐标，即为声源位置。

依据式（2）所示的信号模型，此算法假定输入信号与噪声不相关，在理想条件下，接收数据协方差矩阵可表示为：

$ \begin{split} {\boldsymbol{R}} = \,& E\left[ {{\boldsymbol{x}}\left( {{\omega _n}} \right){\boldsymbol{x}}{{\left( {{\omega _n}} \right)}^H}} \right] \hfill = {\boldsymbol{A}}E\left[ {{\boldsymbol{s}}\left( {{\omega _n}} \right){\boldsymbol{s}}{{\left( {{\omega _n}} \right)}^H}} \right]{{\boldsymbol{A}}^H} + \\ &E\left[ {{\boldsymbol{n}}\left( {{\omega _n}} \right){\boldsymbol{n}}{{\left( {{\omega _n}} \right)}^H}} \right] \hfill = {\boldsymbol{A}}{{\boldsymbol{R}}_{\boldsymbol{s}}}{{\boldsymbol{A}}^H} + \sigma _n^2{\boldsymbol{I}} \hfill \text{。} \end{split} $

(8)

其中： $ {{\boldsymbol{R}}_{\boldsymbol{s}}} = E\left[ {{\boldsymbol{s}}\left( {{\omega _n}} \right){\boldsymbol{s}}{{\left( {{\omega _n}} \right)}^H}} \right] $ ，为输入信号的协方差矩阵； $ \sigma _n^2 $ 为噪声方差。

根据线性代数知识有 $ \left| {{\boldsymbol{R}} - {\lambda _j}{\boldsymbol{I}}} \right| = 0 $ ，其中 $ {\lambda _j} $ 为接收数据协方差矩阵 $ {\boldsymbol{R}} $ 的特征值，将式（8）代入得：

$ \left| {{\boldsymbol{A}}{{\boldsymbol{R}}_{\boldsymbol{s}}}{{\boldsymbol{A}}^H} + \sigma _n^2{\boldsymbol{I}} - {\lambda _j}{\boldsymbol{I}}} \right| = \left| {{\boldsymbol{A}}{{\boldsymbol{R}}_{\boldsymbol{s}}}{{\boldsymbol{A}}^H} - ({\lambda _j} - \sigma _n^2){\boldsymbol{I}}} \right| = 0 ,$

(9)

式中， $ {\lambda _j} - \sigma _n^2 $ 为矩阵 $ {\boldsymbol{A}}{{\boldsymbol{R}}_{\boldsymbol{s}}}{{\boldsymbol{A}}^H} $ 的特征值。假定输入信号数目为 $ K $ ，当输入信号数目小于阵元数 $ M $ 时，得到矩阵 ${\boldsymbol{A}}{{\boldsymbol{R}}_{\boldsymbol{s}}}{{\boldsymbol{A}}^H}$ 的秩为 $ K $ ，等价于此矩阵有 $ M - K $ 个零特征值。

样本数在无穷大情况下，接收数据协方差矩阵 $ {\boldsymbol{R}} $ 就会有 $ M - K $ 个特征值 $ {\lambda _j} $ 等于噪声方差 $ \sigma _n^2 $ 。然而现实情况下样本数不可能无穷大，所以只能得到 $ M - K $ 个非常接近 $ \sigma _n^2 $ 的特征值。将 $ M - K $ 个近似接近 $ \sigma _n^2 $ 特征值对应的特征向量组成的向量矩阵称为噪声子空间，其余特征向量组成的向量矩阵为信号子空间，下面具体分析它们之间的关系。

设 $ {q_j} $ 为小特征值对应的特征向量，则有 $ \left( {{\boldsymbol{R}} - {\lambda _j}{\boldsymbol{I}}} \right){{\boldsymbol{q}}_{\boldsymbol{j}}} = $ $ 0 $ ，其中 $ {\lambda _j} = \sigma _n^2 $ ，化简后得：

$ {\boldsymbol{A}}{{\boldsymbol{R}}_{\boldsymbol{s}}}{{\boldsymbol{A}}^H}{{\boldsymbol{q}}_{\boldsymbol{j}}} = 0 \text{。}$

(10)

因为 $ {\boldsymbol{A}} $ 是满秩， $ {{\boldsymbol{R}}_{\boldsymbol{s}}} $ 非奇异，所以有 $ {{\boldsymbol{A}}^H}{{\boldsymbol{q}}_{\boldsymbol{j}}} = 0 $ 。上式表明 $M - K$ 个小特征值对应的特征向量与拷贝场向量相互正交，即噪声子空间与信号子空间相互正交，借助这个关系将其应用到水下声源匹配场定位中。先将实际声场数据进行特征分解得到噪声子空间，然后利用信号子空间与噪声子空间相互正交关系将噪声子空间和网格点处的拷贝场取相关，此拷贝场由声场传播模型计算得到，取相关之后，式（12）的最大值所对应的拷贝场坐标就是声源实际位置。

噪声子空间用矩阵表示为：

$ {\boldsymbol N}_{\boldsymbol n}={\left[{\boldsymbol q}_{\boldsymbol k},{\boldsymbol q}_{\boldsymbol {k+1}},\cdots ,{\boldsymbol q}_{\boldsymbol {M+1}}\right]}^{{\rm{T}}}, $

(11)

当真实声源距离和深度为 $\left( {r,z} \right)$ 时，应当满足 $ {\boldsymbol{p}}_{\boldsymbol{c}}^H\left( {r,z} \right){{\boldsymbol{N}}_{\boldsymbol{n}}}{\boldsymbol{N}}_{\boldsymbol{n}}^H{{\boldsymbol{p}}_{\boldsymbol{c}}}\left( {r,z} \right) = 0 $ ，即 $ {{\boldsymbol{p}}_{\boldsymbol{c}}}\left( {r,z} \right) $ 和 $ {{\boldsymbol{N}}_{\boldsymbol{n}}} $ 的正交性使得 $ {\boldsymbol{p}}_{\boldsymbol{c}}^H{{\boldsymbol{N}}_{\boldsymbol{n}}}{\boldsymbol{N}}_{\boldsymbol{n}}^H{{\boldsymbol{p}}_{\boldsymbol{c}}} $ 有最小值0。

由此可得基于矩阵空间特征分解的匹配场声源定位算法表达式为：

$ B = \frac{1}{{{\boldsymbol{p}}_{\boldsymbol{c}}^H{{\boldsymbol{N}}_{\boldsymbol{n}}}{\boldsymbol{N}}_{\boldsymbol{n}}^H{{\boldsymbol{p}}_{\boldsymbol{c}}}}} \text{。}$

(12)

2 仿真与分析 2.1 仿真计算的环境模型

采用如图1所示的海洋环境进行仿真。假设海域水深为60.5 m，沉积层为10 m，上表面声速为1590 m/s，下表面声速为1598 m/s，沉积层密度为1.65 g/cm³，衰减系数为 $0.35\;{\rm{dB}}/\lambda $ ；海底位于水下70.5 m处，声速为1650 m/s，密度为1.8 g/cm³，衰减系数为 $0.2\;{\rm{dB}}/\lambda $ 。

2.2 仿真及性能评价 2.2.1 单声源定位仿真

用图1海洋环境参数编辑环境文件，然后用kraken声场计算模型对于500 Hz位于5822 m距离、30 m深度的声源计算数据场，8阵元垂直接收阵均匀分布在11.5 ～60.5 m水深范围内。类似地采用网格搜索，搜索水平距离3800～8000 m，步长1 m，搜索深度0 ～60.5 m，步长0.5 m计算拷贝场。然后分别用3种处理器仿真得到模糊度平面。

从图2（a）看，Bartlett处理器定位出现了旁瓣，且真实声源附近存在模糊性，无法准确确定声源深度和距离；由图2（b）看出，与Bartlett处理器相比最小方差处理器有效抑制了旁瓣，但在真实声源附近也有一定的模糊性，仍然无法非常准确定位声源；从图2（c）可以看出，基于特征分解的处理器实现了声源的准确定位。选择一定的信噪比，通过50组的仿真数据统计得出了3种处理器在不同信噪比情况下的声源定位准确度。根据文献[18]，如果模糊度平面的最大值位于真实位置附近（距离方向±400 m，深度方向±10 m）说明定位成功，声源定位准确度定义为测量数据中定位成功的比例。

从图3看，基于特征分解的处理器定位准确度整体上高于最小方差处理器，最小方差处理器的定位准确度又整体上高于Bartlett处理器。

2.2.2 双声源定位仿真

仿真条件不变，假定有2个声源：声源1深度30 m，距离5822 m；声源2深度34 m，距离5822 m。

图4中的两声源在深度上相距4 m，而距离相同。声源相距很近时，图4（a）的Bartlett处理器结果变得很差，无法确定声源的数目和位置；由图4（b）可看出，最小方差处理器有一定的分辨力，但它在深度方向分辨力仍然不足，无法将2个声源准确分开；从图4（c）看出，基于特征分解的处理器有很好的辨别能力，可以实现声源的准确定位。图5为3种处理器在不同信噪比下对应的定位准确度，可以看出，基于特征分解的处理器仍然性能最好，其次是最小方差处理器，最后是Bartlett处理器。值得注意的是随着信噪比增大，Bartlett处理器的定位准确度没有明显变化，且准确度都很低，从而说明Bartlett处理器缺乏对双声源的定位能力；最小方差处理器相比在单声源定位时，它的定位准确度有所降低；基于特征分解的处理器随着信噪比的增大，准确度不断提高，信噪比在10 dB时准确度达到了单声源定位时的大小。综上所述，基于特征分解的处理器是一种高分辨算法，在双声源距离较近的情况下仍然可以实现声源的定位，且比最小方差处理器性能更优越。

2.2.3 环境失配对定位结果的影响

1）海水深度失配

由于潮汐与海面波浪等可能会导致海平面上升，传播模型计算得到的拷贝场出现失配，从而使匹配场定位结果产生误差，甚至定位错误。图6给出了3种处理器定位准确度与海水深度增量的关系。

图6横坐标是海水增长量，从图中观察到随着海深误差的增大，3种处理器的定位准确度都呈下降趋势。但可以看出最小方差处理器和基于特征分解的处理器的定位准确度依然要高于Bartlett处理器。这2种算法在海洋环境深度失配不严重时的定位性能仍然要优于Bartlett处理器。此现象说明，对于海水深度受潮水影响进而影响垂直阵阵元深度问题的情况，3种匹配场定位算法的定位能力没有发生质变。比较2个高分辨算法定位准确度的变化情况，从图中看出，基于特征分解处理器性能仍优于最小方差处理器。

2）基阵位置偏移

在实际海洋环境定位中，经常会由于洋流缘故使接收阵元位置偏移，从而导致定位的不准确。下面从最简单的偏移情况出发仿真分析。基阵最下方阵元在海底固定，其他阵元分别向右偏移一定位置，使接收阵所有阵元处于同一条直线上，我们将这条直线与没发生偏移时阵元所在垂直线的夹角称为偏移角。

如图7所示，Bartlett处理器从一开始定位准确度就没有另外2种处理器高，但是从整体上看它的变化幅度较小，所以Bartlett处理器对环境失配的敏感性较低。基于特征分解的处理器和最小方差处理器随着角度的增大变化幅度相似，在角度较小时，前者的准确度略好于后者。当偏移角大于6°时，这2种处理器的定位准确度均低于Bartlett处理器，这说明在失配较严重时，这种高分辨算法表现出了较差的定位性能，它们对此失配情况较为敏感。

3）海水声速失配

海水声速作为传播模型中一个重要参数，它测量的准确与否会直接影响到拷贝声场的精确程度。由于匹配场定位非常依赖于拷贝场，从而可能会导致定位结果的不准确。

由图8看出，Bartlett处理器的定位准确度随着声速的失配变化很小，说明它对海水声速失配有更好的稳健性。但基于特征分解的处理器和最小方差处理器的定位准确度比它高。从这3种处理器的定位准确度来看，海水声速的失配相比较海水深度和接收阵偏移对它的影响很小，有很大的失配量，声源定位准确度变化幅却较小。

4）海底声速失配

在匹配场定位过程中，由于海底介质很难实现采样分析，所以海底声速的测量准确与否就会对定位结果产生影响，从而带来定位误差。

从图9可以看出，海底声速虽然变换范围很大，增量达到了−100～100，但是其定位准确度却变化的非常小，只在0.02的范围内变化，这就说明海底声速的变化对定位影响非常小，表现在模糊度平面上因海底声速失配而产生的变化基本可以忽略。

综合以上4种情况的分析，可以得出结论：因洋流导致基阵发生偏移对声源定位结果的影响最大，其次是海水深度的失配，再者是海水声速的失配，海底声速的失配对定位准确度影响最小。而且基于特征分解的处理器和最小方差处理器在以上情况下定位准确度的变化趋势比较相似，但前者的定位准确度又略高于后者。从稳健性来看，这两种高分辨算法对环境失配的敏感性要高于Bartlett处理器。总体来看，基于特征分解的处理器在环境发生失配时有较好的稳健性，可以实现声源的准确定位。

3 实验验证

采用陵水海域某次垂直阵定位实验数据对定位算法进行验证。实验环境即为图1所示的海洋环境，其中垂直阵含有32个阵元，阵元间距为1.5 m，均匀分布在12.5～58 m水深中。采用频率相近的双声源作为目标，声源1深度7 m，距离5300 m，频率为710 Hz。声源2深度40 m，距离5300，频率为690 Hz。第1组实验时，两声源信噪比分别为170 dB，160 dB。第2组实验时，两声源信噪比都为170 dB。

利用kraken声场传播模型生成拷贝场，对真实声源位置附近作网格化处理，搜索水平距离3800 ～8000 m，步长1 m，搜索深度0 ～60.5 m，步长0.5 m。然后分别用3种处理器对实验数据进行定位得出模糊度平面。

比较图10和图11的模糊度平面，首先观察图10（a）和图11（a）可以看出，无论2个声源是相同噪声比还是不同噪声比，Bartlett处理器只能分辨出一个声源，且在2个声源信噪比不同时，它会定位出信噪比较高的那个声源，缺乏双声源的定位能力。比较图10（b）和图11（b）可以看出，两声源信噪比不同时，最小方差处理器也无法定位双声源，它只定位出了信噪比为170 dB的声源，但当2个声源都为170 dB时，从图11（b）看出，最小方差处理器可以定位得到2个声源的位置，但对40 m深处声源的定位很模糊，没有7 m处的定位效果清晰。然而基于特征分解的处理器在这2种情况下定位都很准确，既可以定位出相同信噪比的2个声源，也可以定位出不同信噪比的2个声源。此结果说明，基于特征分解的处理器能够实现存在海面强干扰条件下的水下弱目标定位。

4 结　语

本文将基于矩阵空间特征分解的算法应用于水下声源匹配场定位，与Bartlett处理器相比，可以有效抑制旁瓣，提高了声源定位精度。与最小方差处理器相比，更适用于双声源定位，尤其是在声源的定深方面，能分辨出深度相近的2个声源。其次，本文还分析了3种处理器在不同环境失配条件下的定位能力，得出结论：Bartlett处理器对环境变化有很好的宽容性，但对深度分辨力不足，定位精度不高；最小方差处理器相比于Bartlett处理器对旁瓣有很好的抑制作用；基于特征分解的处理器与最小方差处理器在环境失配时性能相近，但定位准确度又略高于最小方差处理器，尤其是双声源定位时，最小方差处理器可以辨别声源数目，但无法准确定位它们的深度，而基于特征分解的处理器可以在辨别声源数目的情况下准确定位声源位置。海上实验数据表明，基于特征分解的处理器能够实现存在海面强干扰条件下的水下弱目标定位。