2.5 水声信号定位方法

2.5.1 波达方位（DOA）估计方法

声呐的重要任务之一是使用水听器阵列进行DOA估计。常见的DOA估计方法可以分为谱估计类算法和参数化算法两类方法。谱估计类算法主要有CBF算法、Capon波束形成方法、子空间类方法（如MUSIC方法、ESPRIT方法等）；参数化算法主要有确定性最大似然方法和随机最大似然方法等。

以下针对均匀线阵列说明 DOA 估计方法原理。假设 K 个中心频率为 f 的信号s （ t）=[^s₁（t），…，s_K（t）]^T以角度θ=[θ₁，…，θ_k]^T入射到一个水听器数量为M（K＜M）、阵元间距为d的均匀线阵列上。假设远场目标信号为窄带信号，信号包络在信号通过阵列期间不会发生改变。阵列接收信号模型可以表示为

式中， A=[^a（θ₁），…，a（θ_k）] 是阵列流形矩阵，是对应于θ_k方向的阵列流形向量， c是声速，ω=2πf 是信号圆频率，n（t）=[ⁿ₁ （t），…，n_M（t）]^T是阵列接收到的噪声向量。

1）CBF算法

在远场下信号以平面波形式入射到各个阵元上，各阵元接收到信号的相位差只与来波方向和相邻阵元的间隔有关。对不同的阵元，按照特定方向进行时延补偿可以使得各个阵元输出端的信号进行同相叠加，从而使得阵列在此方向产生一个主瓣波束。若扫描的方向与来波方向相同，则此方向波束的能量较强；反之，则较弱，这便是常规波束形成进行方位估计的基本原理。

常规波束形成利用对阵元接收信号进行复加权实现对阵列的时延补偿，假设加权向量为w（θ）=[^w₁（θ），…，w_M（θ）]^T，则波束形成后阵列的输出为

波束形成器的输出功率可以表示为

式中，R=E[x（t） x（t）^H]是阵列接收数据的协方差矩阵，实际中无法求得R，通常使用采样协方差矩阵 ˆR来代替

式中，L代表数据快拍数。常规波束形成器的加权向量为

其值与θ方向对应的阵列流形向量a（θ）相同，因此常规波束形成器的输出功率可以表示为

当目标方向为 _sθ时，常规波束形成定位半功率点波束宽度为

从（2-51）中可以看出，常规波束形成器的方位分辨率随阵列的孔径增大而提高。对于两个频率相同的信号，若其入射角的间隔小于波束宽度，则常规波束的形成将无法进行区分，只能通过增加阵列孔径的方法提高方位分辨率。

2）Capon波束形成方法

Capon 波束形成方法是一种自适应波束形成算法，与常规波束形成相比，具有高方位分辨率、低旁瓣的特点。Capon波束形成方法的目标函数为

使用Lagrange乘子法，式（2-52）可以转化为无约束优化问题，即

与式（2-50）相比，Capon 波束形成方法最小化波束形成器的功率输出并约束扫描方向的能量为 1，它抑制了从扫描方向以外入射到阵列的干扰和噪声。式（2-53）对w_MV求偏导数可得

因此

将约束条件代入式（2-54）可得

将式（2-56）代入式（2-54）可得

Capon波束形成器的输出功率为

实际中使用采样协方差矩阵代替R，在目标信号相干、快拍数较小等情况下，Capon波束形成器性能会严重下降。为了提高Capon波束形成器的稳健性，通常在求逆前对的对角线进行对角加载，即

式中，ρ为加载量，I为单位矩阵。

3）MUSIC方法

MUSIC方法是最经典的子空间类方法，子空间类方法的基本思想：在信号个数少于传感器个数的情况下，对阵列接收数据的协方差矩阵进行特征分解，将特征向量划分为信号子空间和噪声子空间两个相互正交的子空间。若子空间被正确地划分，则真实来波方向对应的阵列流形向量将位于信号子空间且与噪声子空间正交。利用此正交性可以构造一个尖锐的目标函数，从而达到超分辨的目的。

假设噪声为白噪声且服从均值为0、方差为σ²的高斯分布，噪声与信号不相关，根据式（2-45），协方差矩阵可以写为

式中， R₁=AR_sA^H ， R_s=E[^ss^H]是源信号协方差矩阵，假设信号互不相关，那么有

因此R_s是秩为K的满秩矩阵，由于M＞K，rank （ AR_sA^H ）=K，因此R₁中有M-K个为零的特征值，令其非零特征值为，因此R的特征值可以分为两组，即

令s_k代表 R 的特征值_kλ对应的特征向量， S=[s₁，s₂，…，s_K]代表信号子空间， G=[s_K+1，…， s_M]代表噪声子空间。由于G中所有列向量对应的特征值为σ² ，则

根据式（2-61），有

因此有

对式（2-67）等号两端同乘G^H ，可得

由于R_s为正定矩阵，式（2-68）成立的充要条件为

式（2-69）说明阵列流形矩阵中的各个列向量与噪声子空间G正交，因此可以构造方位估计目标函数：

当θ扫描到真实目标方位时， P_MUSIC （θ）将会出现尖锐的峰值。

2.5.2 匹配场信号处理

匹配场处理器是对平面波常规波束形成器在海洋信道中的推广，其基本思想是对阵列实际的测量信号与信号所有可能出现位置产生的拷贝场信号进行匹配。拷贝场信号可以通过2.2节中的各种传播模型进行数值计算得到。与2.5.1节中各种方位估计算法不同，匹配场处理器的性能不仅与定位目标函数的具体形式有关，还强烈依赖具体的海洋环境。匹配场处理器利用声场的干涉模式进行定位，具体过程包括对扫描网格位置放置测试声源、使用传播模型计算测试声源在阵列各个阵元产生的声场并与阵列实际测量的到达声场数据进行相关联，当测试声源放置到真实声源位置时，相关取得最大值。

在频域上，类似于常规波束形成器，常规匹配滤波器（巴特雷匹配滤波器）是对阵列数据按照假设的海洋环境和声源位置进行加权，即

式中，a是扫描位置参数（通常包括深度和距离）；a_T是真实的声源位置；p （ a_T ）为频域上阵列的测量数据；K=E[^p（a_T ） p^H （a_T ）]是互谱密度矩阵；w（a）是对应扫描位置参数a的加权向量，通过对波动方程进行数值求解得到，通常w（a）须进行归一化，即 w（a）||₂ =1。式（2-71）中最大值对应的a即||a_T的估计值，但与常规波束形成器类似，常规匹配滤波器也存在低分辨率、高旁瓣等缺陷。为了提高空间分辨率和抑制旁瓣，可以将高分辨方位估计算法的思想推广到匹配场处理中，从而导出多种自适应匹配场处理器。

类似 Capon 算法，自适应匹配场处理器^[133]的权值由拷贝场和接收数据协方差矩阵共同决定，能够有效地抑制旁瓣和干扰。Capon匹配场处理器的目标函数为

由于高分辨处理器对所有的参数都较敏感，环境参数的细微不确定性通常会严重降低估计的精确程度，因此环境容忍能力对匹配场处理器至关重要。在Capon匹配场处理器的基础上又发展出邻域位置多约束匹配场处理器^[134]和环境扰动约束匹配场处理器^[98]。邻域位置多约束匹配场处理器认为环境扰动会使匹配场定位结果在一定的领域位置内变化，因此，在Capon匹配场处理器的基础上对多个领域位置加入定值约束，适用于海深等参数失配的情况。环境扰动约束匹配场处理器首先使用一定范围内环境参数扰动后的拷贝场作为约束矩阵，利用约束矩阵构造拷贝信号的相关矩阵，对相关矩阵进行特征分解以提取出环境参数扰动的一、二阶统计特性；最后利用环境参数的一、二阶统计特性构造自适应匹配场处理器的权值向量。与邻域位置多约束匹配场处理器相比，环境扰动约束匹配场处理器具有更大的环境参数宽容性。

上述匹配场处理都是针对单频情况的，宽带信号情况下可以直接对常规匹配滤波器在不同的频率上输出的功率进行求和

式中， N _f为频点个数。式（2-73）是一种声压数据的宽带非相干处理方法，也可以使用声压数据的二阶统计量进行匹配。假设M_p（f，a）和M_q（f，a）分别代表声源位置为a时通过声场模型计算得到的第p个和第q个水听器接收到信号的频谱，M_p，q（f，a）=M_p（f，a）M_q（f，a）^*代表拷贝互功率谱函数；D_p（ f，a_T）和D_q（ f，a_T）是第 p 个和第 q 个水听器实际测量到信号的频谱， D_p，q（ f， a_T ）=D_p（ f， a_T ）D_q（ f， a_T ）^*代表实测互功率谱，则 Westwood 宽带匹配场处理器可以表示为

式中，

在频域对测量场和拷贝场的互谱矩阵进行相关并进行宽带相干求和。

使用声压场进行匹配场定位时，一般需要使用垂直/水平线阵列对声场的干涉模式进行采样。除此之外，还可以利用单阵元多途到达时延的差异进行匹配定位。在声线模型中，某一点的声压可以看作从声源发射的多条声线在该点叠加的效果，假设z_s ， z_r和r分别代表声源深度、接收器深度及声源与接收器之间的水平距离，则声源到接收器之间的信道响应可以表示为

式中，N是本征声线的条数，g_i代表参数为（z_s，z_r，r）时第i条特征声线的幅度，n_i代表参数为（ z_s ， z_r ， r）时第i条特征声线的到达时延。当声源信号为s（n）时，水听器收到的信号x（n）可以表示为

x（n）具有N-1个时延，其时延差向量可以表示为

式中，ΔT_i代表实际接收信号中第i+1个时延和第i个时延之差。类似声压匹配的流程，使用射线模型计算得到的拷贝时延差向量可以表示为

式中，Δ_iτ代表拷贝时延向量中第i+1个时延和第i个时延之差。对实测数据时延差和拷贝时延差进行匹配，其目标函数为

式（2-80）可以很容易地扩展到阵列处理上，假设阵列有 J 个水听器，则时延差匹配的阵列处理目标函数可以表示为

式中，d_j和O_j分别表示第j个水听器对应的实测数据到达时延差向量和拷贝到达时延差向量。

使用到达时延差进行匹配具有计算量小、物理意义直观的优点，但目标函数式（2-80）对模型误差非常敏感，在接收信号中若缺少或多出一个到达时延，则各个时延差的对应关系变乱， F （a）会发生剧烈变化。

当只有单水听器且接收信号未知时，可以采用自相关函数来提取耦合的时延关系，然后与拷贝场的自相关结构进行匹配。假设声源信号与噪声不相关，根据式（2-77）， x （ n）的自相关函数可以表示为

式中， n_i=n _j=0表示第一个相对时延，为声源信号的自相关函数，为噪声的自相关函数。假设s（n）为宽带信号，噪声为高斯白噪声，则当m=0或m=n_i-n _j时， R_xx（m）会出现峰值。式（2-82）中的峰值结构对应了多途到达时延之间的差值，一般是相当复杂的，很难通过式（2-82）反推信号的到达时延。由于N条本征声线会造成数量为N（N-1）/2个自相关峰值出现，即使环境参数剧烈变化，信道响应也相对稳定，因此特定声源信号经信道传输后，接收信号的自相关函数也相对稳定，并且是声源位置、接收器位置及水平距离等因素的函数。对信道响应求自相关函数，可得

可以看出，R_hh（m）和R_xx（m）出现峰值的时间相同，类似于匹配场，使用不同扫描点计算得到的信道函数自相关函数来匹配由实际接收信号的自相关函数得到声源的位置估计，可以避免对接收信号的到达时延进行估计。假设信道响应自相关函数和数据自相关函数均截短为2N+1点，使用如下目标函数进行单水听器匹配场定位：

2.5.3 基于稀疏性表示的DOA方法

假设一个阵元间距为d的M元均匀水听器阵列用于水声信号的定位，如图2-36所示。在平面波假设条件下，K个角频率为ω的声源信号以来波方向向量θ入射到水听器阵列：

θ=[θ₁,θ₂,…,θ_K]^T

此时阵列流形矩阵为

A(θ)=[a(₁θ), a(θ₂ ),…, a(θ_K)]

式中，是导向向量，c是水中的声速。t时刻阵列接收的信号可以表示为

式中， x（t）=[x₁（t），x₂（t），…，x_M（t）]^T是阵列接收信号，u（t）=[u₁（t），u₂（t），…，u_N（t）]^T是源信号向量， n（t）=[n₁（t），n₂（t），…，n_M（t）]^T是噪声向量。

为了将声源定位问题转化为一个信号重构问题，构造新的阵列流形矩阵：

式中，包含了声源所有可能的入射方向，假设所有来波方位离散点数量N一般远大于声源实际个数K和阵元数M。此时，阵列模型可以表示为

式中，s（t）=[s₁（t），s₂（t），…，s_N（t）]^T是假设所有来波方位的信号向量。若第k个实际声源u_k（t）对应的方位θ_k和第n个假设来波s_n（t）对应的方位相同，则u_k（t）=s_n（t）。因此，K个实际声源的情况下，s（t）中只有K个不为零分量，s（t）是K稀疏的。式（2-87）只针对单个观测向量（Single Measurement Vectors，SMV），在离散的多个观测向量（Multiple Measurement Vectors，MMV）情况下，接收信号的模型可以表示为

式中， X=[x（t₁），x（t₂），…，x（t_T）]， S=[s（t₁），s（t₂），…，s（t_T）]， N=[n（t₁），n（t₂），…，n（t_T）]。式（2-87）也被称为x（t）的过完备表示，中的每一个列向量被称为过完备基，被称为过完备基矩阵。式（2-87）将信号用假设来波达到方向的函数表示，因此，信号s（t）的能量谱是稀疏的。由于A的列数远大于行数，使用x（t）对s（t）进行重构是一个欠定问题，有无穷多个解。为了得到一个唯一解，通常需要采用正则化技术。正则化技术通过对优化目标函数增加一个正则化项，来约束优化变量的可行域；从贝叶斯的观点看，正则化技术实际上是对优化变量施加了一定的先验信息。式（2-88）的正则化通用框架为

式中，函数 f （s）代表了测量数据拟合的失真度等代价函数；R（s）为正则项，代表了对s的惩罚度；λ是正则化参数，平衡 f （s）和R（s）的相对重要程度。常见的正则化模型见表2-1。

正则化技术已经广泛地应用于阵列信号处理中。例如，自适应波束形成对角加载的优化目标函数为

式中，w为加权向量，Σ=E「_└^xx^H┐_」^为协方差矩阵，a （θ_s ）为指向 _sθ的导向向量，式（2-90）等价于

因此，对角加载技术是Tikhonov正则化的一种。表2-1中的正则化模型除了Tikhonov模型，均具有促进向量稀疏的作用且都涉及l₁范数优化问题。下面将介绍几种广泛应用于DOA估计的l₁范数优化模型。

1.l₁-SVD

具有促稀疏性的l₁范数正则化模型被广泛地应用于 DOA 估计中，其目标函数可以表示为

式中，λ为正则化参数，其取值与噪声能量有关。l₁范数正则化技术是一种利用稀疏先验信息来提高角测向分辨率的DOA方法，与流行的高分辨方法（如Capon方法、MUSIC算法）相比，具有更好的测向分辨率，能够在相干声源存在时工作，因此受到了广泛关注^[135]。

当信号源在时域采样点处于平稳时，可以通过联合时间处理对多个快拍的数据进行相干融合，得到更精确的方位估计。定义S的l_p，q混合范数|| S||_p，q为

式（2-88）中S在空域具有稀疏性但在时域无稀疏性，可以通过构造以下的优化目标函数来实现联合时间处理

式中，||·||_F表示矩阵的 Frobenius 范数。|| S||_2，1对特定方向上不同时间的数据使用 l₂范数进行融合，对特定时间不同方向的数据使用l₁范数进行融合。因此，只对空域进行稀疏增强。式（2-94）代表的是一个复值凸优化问题，可以通过二阶锥优化等方法求解，但是计算量随T的增大而超线性地增大，当T值很大时，计算量过大，计算时间过长。为了同时降低计算复杂度和噪声的影响，使用 SVD 对多快拍数据进行相干融合，对接收信号矩阵S进行奇异值分解，即

式中，X的维度为M×T；U的维度为M×M，是X的左奇异矩阵；Y的维度为M×T，是X的奇异值矩阵；V的维度为T×T，是X的右奇异矩阵。假设Y对角线上的元素按照从大到小排列，令V=[V_KV_T-K]。将X分解为信号子空间和噪声子空间，信号子空间对应前K个奇异值，噪声子空间对应后T-K个奇异值。对阵列模型即式（2-88）等号左右两边均右乘V_K：

此时，信号矩阵的维度降低到了M× K。因此，式（2-94）可改写为

此方法称为l₁-SVD方法。与MUSIC算法类似，进行信号子空间划分时，l₁-SVD需要知道声源个数K，但是L1-SVD算法并未利用信号子空间与噪声子空间的正交性，因此K的误差不会对定位结果造成剧烈影响。式（2-97）也被称为问题^[136]，它是一个凸优化问题，其等价于

式（2-98）是一个二阶锥优化问题，可以采用内点法来求解。参数λ的选取与噪声能量有关，实际上，给定一个λ值，当噪声能量上界δ满足一定条件时，式（2-92）的解与以下问题的解相同，即

式（2-99）也被称为基追踪去噪问题，同样可以使用二阶锥规划等凸优化方法进行求解，当阵元数M值较大时，正则化参数λ通常可以选取为^[137]

2.稀疏行范数重构

l₁-SVD算法利用奇异值分解，对多个观测数据进行噪声子空间和信号子空间的划分，通过舍弃噪声子空间部分来达到降低MMV问题重构的复杂度。除了l₁-SVD，还存在另一类能够有效计算MMV问题的算法，即稀疏行范数重构算法。与l₁-SVD算法不同的是，它无须知道全部观测数据，只须观测数据协方差矩阵即可对信号进行定位。稀疏行范数重构算法适用于某些由于存储空间和传输数据率的限制而不能传输原始观测信号的应用场景^[137]。

对式（2-94）进行变形，则有

可以证明，式（2-101）等价于

式中，是采样协方差矩阵， D₊代表非负实对角矩阵的集合，注意：，W 中的特定元素代表了此方向多个时间接收的来波平均功率。假设一个维度为T× T埃尔米特矩阵U_T满足以下条件

矩阵M≥0，表示M半正定，由于X^H （ AWA^H+λI_M）^-1X半正定，所以U_T也半正定，则

所以

式（2-102）等价于

因为（ AWA^H+λI_M）^-1和U_T都是埃尔米特矩阵，根据 Schur 补定理^[138] ，矩阵正定等价于U_T-X^H （ AWA^H+λI_M）^-1X≥0^[138]，式（2-106）可以转化为一个半正定规划的问题，即

类似地，式（2-106）还可以转化为

式（2-107）中U_T的维度为T，式（2-108）中U_M的维度为M，当阵元数M大于等于观测数量T时，选用式（2-107）进行优化；反之，则使用式（2-108）进行优化。由于式（2-108）不需要知道观测数据矩阵X，只与采样协方差矩阵 ˆR有关，计算复杂度与 T无关，因此适用于T值很大的情况。

3.稀疏谱拟合

稀疏行范数重构方法的出发点仍然在数据域的l_2，1混合范数优化问题的求解上，仍是利用S的块稀疏性。考虑到协方差矩阵R=E[xx^H]，R的维度为N×N，使用采样协方差矩阵代替R，即

式中，L为多观测的快拍数，若令为源信号的采样协方差矩阵，则式（2-109）可以改写为

R_s中的对角元素表示特定方向信号的功率，在K 个信号情况下， R_s中的非零元素个数最多为K²（此时信号全部相干），最少为K（此时信号全部不相干）。因此，源信号采样协方差矩阵R_s是稀疏的，可以通过对测量信号协方差矩阵进行稀疏重构来得到R_s，从而得到信号的DOA估计。在信号相干情况下，重构目标函数为

在信号非相干情况下，可以表示为

式中，p_n代表了θ_n方向入射的信号功率。若令p=[p₁， p₂，…， p_n]^T ，a_v（θ_n）=vec[^a（θ_n）a（θ_n）^H]， A_N=[a_v（₁θ），…a_v（θ_N）]，则非相干情况下的重构目标函数为

式（2-111）与式（2-113）都可以通过二阶锥优化来进行求解，此方法也被称为稀疏谱拟合^[139]（Sparse Spectrum Fitting，SSF）。SSF方法利用了源信号协方差矩阵的稀疏性，在处理MMV问题时计算量远小于l₁-SVD方法。