矢量水听器由传统的无指向性标量水听器(声压水听器)和具有与频率无关的偶极子自然指向性的质点振速水听器复合而成，它可同步测量声场中的声压和质点振速信息，并利用这些信息实现对目标的无模糊定向^[1-2]、多个不同频率目标的分辨^[3]，具有旁瓣低、检测性能和参数估计精度高的优势^[4]。因此，矢量水听器相较于标量水听器具有更广泛的实用性。

近年来，有关矢量水听器及其信号处理技术的研究成为相关研究人员关注的热点。文献[5]针对远场平面波和高斯白噪声，提出了时域矢量水听器阵列信号处理数据模型，研究表明利用矢量水听器阵可以实现无左右舷模糊定向。文献[1-2]把矢量水听器阵列应用到常规波束形成(conventional beamforming, CBF)，最小方差无畸变响应(minimum variance distortionless response, MVDR)，多重信号分类(multiple signal classfication method, MUSIC)方位估计中，提高了目标波达方向(direction of arrival, DOA)估计的精度。文献[6]研究了矢量水听器阵列的空间分辨能力，从理论上推导了矢量水听器阵列所能分辨的最大目标个数，表明单个矢量水听器最多能分辨2个相干目标。目前，有关矢量水听器阵列处理非相干信号的算法层出不穷：文献[7]针对色噪声背景下矢量水听器方位估计性能下降问题，提出了不需要加权参数设置的噪声协方差矩阵重构算法，提高了矢量水听器在色噪声背景下的方位估计性能；文献[8]基于克罗内克积，文献[9]基于嵌套矢量水听器阵列，分别提出了目标数目大于阵元数目欠定条件下矢量水听器阵列的估计方法；文献[10-11]利用解析振速与声压-振速抑制各向同性噪声能力，研究了声压振速联合处理的方位估计方法，但是该方法无法估计目标个数大于阵元个数的情况，损失了矢量阵的自由度。矢量水听器阵列处理非相干信号的研究为矢量水听器定向的工程应用奠定了基础。

然而，有关矢量水听器阵列处理相干信号的研究鲜少。针对标量水听器而言，空间平滑算法是处理相干信号常用的方法之一，该方法是以牺牲部分阵列孔径为代价。当矢量水听器个数较少时，该类算法失效。相较于标量水听器阵列而言，矢量水听器自身具有单边指向性，基于这一原理文献[12]提出了一种振速域平滑算法，文献[13]提出了一种振速域差分平滑算法，振速信息用于相干项抑制，其单边指向特性损失，无法实现全空间无模糊定向。文献[14]在文献[12-13]的基础上提出了一种改进算法，解决了振速域平滑算法与振速域差分平滑算法存在左右舷模糊问题。文献[15]将压缩感知理论应用于矢量水听器阵的DOA估计中，提出了基于压缩感知技术的时空联合滤波高分辨方位估计方法，对快拍数与信号的相关性具有较低的敏感度，有效地提高了相干信号方位估计性能。然而，该方法受噪声功率影响较大，不合理的噪声功率约束，将无法得到正确的DOA估计。另外，该方法是利用CVX从观测数据域求解DOA，随着快拍数的增加计算量较大。此外，文献[12-15]的方法都无法估计声源个数大于等于矢量传感器个数的情况。文献[16]利用向量化的思想先把协方差矩阵向量化再利用稀疏表示的方法求解DOA。该方法解决了相干源的难题。然而，在求解稀疏信号的过程中该方法需要利用CVX工具恢复稀疏信号。凸优化的方法面临计算量较大及需要合理的正则化参数设置才可获得高精度的DOA估计。然而，在复杂的水声环境中，正则化参数难以做到恰当选择，将会导致方位估计性能急剧下降。针对文献[16]正则化参数难以选择以及凸优化方法计算量较大的难题，文献[17]提出了一种无需参数设置的迭代自适应方法(iterative adaptive approach, IAA)。相较于文献[16]利用凸优化方法恢复稀疏信号，IAA方法通过迭代的方式恢复稀疏信号减少了计算量并且不需要正则化参数的选择，被广泛应用于雷达系统中。然而，IAA方法面临分辨率及估计精度不高的难题。此外，文献[16-17]都是基于标量传感器阵列的DOA估计性能进行的研究。

为此，本文把文献[17]的方法扩展到矢量水听器阵列中，有效解决了现有方法在解决相干源问题时损失自由度的难题。并且针对IAA方法分辨率较低及DOA估计精度不高的问题，提出了迭代稀疏协方差矩阵拟合算法(iterative sparse covariance matrix fitting algorithm, ISCMFA)。基于加权协方差矩阵拟合准则，构建了关于稀疏信号功率的目标函数，利用Frobenius范数的性质推导了稀疏信号功率迭代更新的递推式。相较于文献[17]采用迭代方式恢复稀疏信号，所提方法采用Frobenius范数从样本协方差矩阵中直接估计稀疏信号的功率减少了多快拍带来的计算负担。理论分析结果表明所提算法对信号相关性具有较低的敏感性。仿真实验结果表明，在相干源与非相干源情况下，所提方法相较于CBF、MVDR、MUSIC及IAA方法提高了DOA估计的精度。

1 问题描述 1.1 矢量水听器阵列数据模型

考虑K个来自不同方向(θ_k, k=1, …, K)的远场等功率窄带信号源入射在有M个矢量水听器组成的均匀线列阵上, 矢量阵列配置如图 1所示。

每个矢量水听器由2个振速水听器与1个标量水听器组成, 振速水听器与标量水听器在空间位于同一点。第k个声源入射在矢量水听器上时3×1的阵列流形定义为

(1)

式中, θ_k∈(-180°, 180°)为第k个声源的方位角。u_xk, u_yk分别表示第k个声源沿x轴、y轴方向的振速分量。由(1)式可以看出单个矢量水听器阵列流形不包含时延与相位信息, 这是由于矢量水听器的阵列流形与入射信号的频率是相互独立的, 并且假设振速水听器和标量水听器在同一点, 两者没有时延或相位差^[12]。在t时刻, 第m个矢量水听器接收到3×1的矢量数据为

(2)

式中, s_k(t)表示第k个信号的复包络, w_m(t)表示第m个矢量水听器接收到的均值为零、方差为σ²的3×1维平稳复高斯白噪声。q_m(θ_k)=e^{-j2π(m-1)dcosθ_k/λ}为阵元m对第k个声源的声压响应系数, d表示相邻矢量水听器阵元间距, λ表示波长。

M个矢量水听器阵列构成3M×K的阵列流形矩阵可表示为A(θ)=[a(θ₁), …, a(θ_K)], 其中

(3)

式中: q(θ_k)=[q₁(θ_k), q₂(θ_k), …, q_M(θ_k)]^T; “⊗”表示Kronecker积。

在t时刻, 矢量水听器阵列接收到3M×1的数据矢量为

(4)

式中: A(θ)=Q(θ)⊙H(θ); “⊙”表示Khatri-Rao积; Q(θ)=[q(θ₁), q(θ₂), …, q(θ_K)]; H(θ)=[h(θ₁), h(θ₂), …, h(θ_K)]; x(t)=[x₁^T(t), x₂^T(t), …, x_M^T(t)]^T表示3M×1的观测数据; s(t)=[s₁(t), s₂(t), …, s_K(t)]^T为t时刻各声源的波形; w(t)=[w₁^T(t), w₂^T(t), …, w_M^T(t)]^T表示3M×1的噪声矢量。阵列的协方差矩阵可表示为

(5)

式中, P, C分别表示K×K的信号协方差矩阵与3M×3M的噪声协方差矩阵。

(6)

式中, σ²表示噪声功率, I表示3M×3M的单位矩阵。在实际应用中, L次快拍下阵列的协方差矩阵可表示为

(7)

式中: X=A(θ)S+W; X=[x(1), x(2), …, x(L)]为3M×L的观测数据; S=[s(1), s(2), …, s(L)]为K×L的信号数据, 并且W=[w(1), w(2), …, w(L)]为3M×L的噪声数据。

1.2 矢量水听器阵列网格数据模型

在实际应用中, 声源个数相对于整个空间角度范围是有限的, 也就是说, 来波方向的声源具有稀疏分布的特性。因此, DOA估计问题就可以转化为稀疏信号表示问题。假设, 声源在空间角度范围内所有可能来波方向的集合为, 其中N是在空间角度范围-180°~180°之间按照一定的角度间隔划分的网格数, 并且需要满足N≫K。假设声源的真实来波方向落在所划分的网格角度上, 划分网格后的离散阵列流形为

(8)

因此, 稀疏信号的矢量阵列流行矩阵可表示为, 并且是已知的与来波方向无关, 这与(4)式中的A(θ)是有所不同的。于是, 在不考虑噪声的情况下, t时刻矢量水听器阵列稀疏信号矢量可表示为

(9)

式中: L表示快拍数; , 是列矢量, 并且仅有K个非零行。的每一个非零行就对应着一个目标。若通过(9)式, 可估计出, 并且估计值用表示。定义

(10)

式中, 表示第n个网格上信号功率的估计值。因此, 所有网格点上信号的功率。

在实际应用中, 受环境噪声的影响, 不能保证仅有K个非零行。也就是说, 在中不可能仅有K个非零元素。因此, 可以对谱峰搜索来获得DOA的估计值。通过以上分析, 若想求出DOA的估计值, 求解或者是至关重要的。下文将讨论如何求解。

2 迭代稀疏协方差矩阵拟合DOA估计方法

在求解的过程中, 为了降低计算复杂度以及高斯白噪声的影响, 在假设信源数目已知的条件下, 对样本协方差矩阵进行特征值分解, 得

(11)

式中, λ₁, λ₂, …, λ_3M表示3M各特征值, 与之对应的特征向量分别为e₁, e₂, …, e_3M。Σ_S是由K个大特征值组成的对角阵, U_S是由与大特征值相对应的特征向量张成的子空间, 称作信号子空间, 称作信号协方差矩阵; Σ_N为3M-K个小特征值组成的对角阵, U_W是由与小特征值相对应的特征向量张成的子空间, 称作噪声子空间, 称作噪声协方差矩阵。为此, 基于加权协方差矩阵拟合准则, 提出一种从信号协方差矩中求解稀疏信号功率的方法, 获得后对所有网格上的信号功率进行谱峰搜索即可获得DOA估计值。

2.1 基于Frobenius范数的稀疏信号功率迭代估计方法

根据(9)式, 我们可以获得第n个网格上信号的协方差矩阵可表示为

(12)

式中, A_n表示的第n列。并且, 定义第n个网格上干扰加噪声的协方差矩阵为

(13)

并且

(14)

式中, 为N×N的对角阵, 对角线上的元素为相应网格上的信号功率。根据矩阵求逆定理, 对(13)式两边分别求逆可得

(15)

式中, 。为了求解, 基于加权协方差矩阵拟合准则, 构造如下函数

(16)

令, 利用Frobenius范数的性质对(16)式化简可得

(17)

根据(17)式, 可看出是与无关的项, 为此, 对关于求一阶导数, 可得

(18)

令, 可得

(19)

不难发现, 并且, ≥0(因为表示功率, 所以必须满足此条件), 这也就意味着是唯一的最小值, 即是唯一的最小值。

为了减少求解Q_n逆矩阵的计算量, 把(15)式代入(19)式, 得

(20)

式中, 。根据(15)式可知, 1+ α_n≠0。因此式, (20)式可变为

(21)

根据(21)式可以看出, 为求解需要R′的先验知识, 根据(14)式可知, 求解R′需要P的先验知识, 为此所提算法需要进行迭代求解。

不妨假设第j+1次迭代中第n个网格点上的信号功率表示为, 则

(22)

(23)

可通过CBF算法^[1]估计的信号功率进行初始化

(24)

2.2 所提算法相关性性能分析

根据(22)式, 对第n个网格点上第j次迭代获得的信号功率进一步转化为

(25)

式中

(26)

(27)

根据(14)公式, 可知为N×N的对角阵, 因此, R′为Toeplitz矩阵。从而可得

(28)

不妨假设, 有2个相关信号源的来波方向分别落在第f个网格上与第z个网格上, 其中f≠z, 在x中的2个相关信号源可表示为与。根据(27)式, 可知

(29)

观察(29)式, 不难发现w″_n^(j)近似等效为滤波器, 与A_n相对应的信号会无失真地通过该滤波器, 与此同时, w″_n^(j)还会衰减A_f≠A_n, A_z≠A_n。当A_n=A_f时, w″_n^(j)会衰减其他信号(包含与A_z相对应的信号), 与A_f相对应的信号会无失真通过。反之, 当A_n=A_z时, w″_n^(j)会衰减其他信号(包含与A_f相对应的信号), 与A_z相对应的信号会无失真通过。因此, 所提算法对信号的相关性具有较低的敏感度, 避免了传统算法依赖于样本协方差矩阵秩的弊端。

从(28)式可看出, 去掉w″_n^(j)上的迭代上标即可变为w″_n, w″n与MVDR算法的最优权矢量有相似之处:所提算法利用第j次获得的稀疏信号功率重构的协方差矩阵R′替换了MVDR算法最优权值中的样本协方差矩阵, 利用重构的协方差矩阵R′对权值进行更新, 使R′更逼近一个理想的协方差矩阵, 收敛时获得更精确的信号功率。因此, 对网格点上的信号功率进行谱峰搜索时会获得更精确的DOA估计。

2.3 计算复杂度分析

本小节对所提ISCMFA算法及IAA方法的计算复杂度进行分析。首先假设在角度扫描范围[-180°:180°]之间总共需要搜索的次数为G=表示向下取整数, ξ表示扫描角度间隔)。本文中所研究的声矢量水听器是由一个全向的声压水听器与2个振速传感器组合而成的, 因此观测数据的维数为声矢量水听器个数的3倍, 即M′=3M, 其中, M表示矢量传感器的个数, M′表示观测数据的维数。CBF、MVDR和MUSIC算法不需要进行迭代运算具有较低的计算复杂度, 为此, 下文仅对ISCMFA方法与IAA方法的计算复杂度进行分析。

在本文中, 所提ISCMFA方法与IAA方法都是根据(24)式进行初始化, 因此在2种迭代方法在初始化阶段具有相同的计算复杂度, 初始化共需要M′G+(2L+1)G次乘法与M′G+(2L-1)G次加法。所提ISCMFA方法与IAA方法在求解稀疏信号功率的过程中都需要进行协方差矩阵R′的逆运算。在计算每个网格点上的稀疏信号功率时并不需要每次都对R′进行逆运算, 而是在所有网格点上的信号功率都计算完毕后计算一次R′的逆矩阵, 以此在循环中不断更新协方差矩阵R′, 其R′的逆矩阵需要M′²G+M′G+15M′³次乘法与M′²G+15M′³次加法。2种迭代算法计算复杂度的不同主要来源于每一次迭代更新的稀疏信号功率。所提算法每次迭代更新稀疏信号功率需要7M′²G+3M′G+G次乘法与6M′²G+2M′G+G次加法, IAA算法需要7M′²L²G+3M′GL+2G次乘法与6M′²L²G+2M′GL+3G次加法。考虑J次迭代, 2种迭代算法总共需要的计算复杂度见表 1。

从表 1可以看出, 在阵元个数、扫描间隔的总数与迭代次数一定的情况下, 快拍数对IAA方法计算量的影响要高于对所提方法的影响, 主要归因于所提方法是直接从样本协方差矩阵中来恢复稀疏信号的功率, 减少了快拍数对所提算法计算复杂度的影响。

2.4 所提算法流程

采用ISCMFA算法进行方位估计的流程如下:

第1步  根据(7)式, 计算样本协方差矩阵, 对样本协方差矩阵进行特征值分解, 得。

第2步  按照特征值的大小排序, 把与信号个数K个相等的特征值和对应的特征矢量看作信号子空间, 有。

第3步  根据(24)式初始化。其中, n=1, 2, …, N。

第4步  根据(23)式更新R′^(j), 其中, j表示迭代次数。

第5步  根据(22)式, 计算, 其中, 。

第6步  收敛性判断, 如果 < ε或者j≥β, 成立, 则终止迭代, 其中, ε表示误差容忍度, β表示最大迭代次数。否则, j=j+1, 返回第五步。

第7步  通过对峰值搜索得到波达方向的估计值。

3 仿真与性能分析

本小节从收敛速度、空间谱、分辨概率与估计精度4个方面评估所提ISCMFA方法的有效性。在仿真中, 信号的中心频率f=700 Hz, 并且满足远场等功率窄带的条件, 采样频率f_s=50 kHz, 噪声为均值为零方差为1的高斯白噪声, 矢量水听器个数M=2, 阵元间距d=λ, 角度在-180°到180°之间以1°为步长进行扫描, ε=10^-4。

仿真中, 正确分辨需要满足的条件为

(30)

式中, θ_k表示第k个目标的真实方位, 表示第k个目标方位的估计值, B_W表示矢量阵的波束宽度。然而, 鲜少有文献给出矢量传感器阵列B_W的解析式。在这里基于2个矢量水听器按照等波长配置采用CBF算法^[1]在无噪声环境下得到的波束下降3 dB两方向对应的夹角作为2个矢量水听器阵列的B_W, 仿真结果如图 2所示, 图 2可知B_W=20°。

本文以均方根误差作为估计精度的衡量标准, 定义均方根误差为

(31)

式中, θ_k表示第k个目标的真实方位, 表示第j次蒙特卡罗实验获得的第k个目标方位的估计值, J表示独立蒙特卡罗实验次数, 其中, 蒙特卡罗实验次数为1 000。

3.1 收敛性分析

图 3比较了方位分别位于θ₁=-15°, θ₂=(-15+Δ)°的2个非相干源, Δ表示2个声源的角度间隔, 信噪比为10 dB, 快拍数为500时, 所提ISCMFA方法与IAA方法的均方根误差与迭代次数的关系。由图 3可以看出, 所提ISCMFA方法相较于IAA方法具有较快的收敛速度; 当角度间隔相对较小时, 通过迭代次数未能使IAA方法的均方根误差减少, 这主要归因于IAA方法分辨率不高导致的; 随着角度间隔的增加, IAA方法与ISCMFA方法的均方根误差都会逐渐减少并收敛; 在角度间隔较大或者较小时, 2种方法都会在迭代次数为15附近收敛, 为此在下文的仿真中, 最大迭代次数β=15。

3.2 空间谱分析

以所提ISCMFA算法为例, 图 4比较了3个相干声源的来波方向为θ₁=-25°, θ₂=0°, θ₃=25°, 信噪比为20 dB, 快拍数为5 000时, 矢量水听器阵(在空间谱图中标记为V-ISCMFA)与标量水听器阵(在空间谱图中标记为P-ISCMFA)对方位估计的性能差异, 其中, 2个矢量水听器组成线阵, 阵元间距d=λ, 在空间谱图中标记为2V-ISCMFA, 4个标量水听器、5个标量水听器及6个标量水听器分别组成均匀线阵, 阵元间距d=0.5λ, 在空间谱图中分别标记为4P-ISCMFA, 5P-ISCMFA, 6P-ISCMFA。图 4给出了不同阵列配置下的空间谱图。

从图 4可以看出, 矢量水听器具有单边指向性, 其在全空间扫描范围内无模糊现象; 无方向性的标量水听器在真实方位的镜像方位上出现伪峰; 标量阵中, 阵元个数越多谱峰越尖锐; 2个矢量水听器阵列的主瓣宽度明显窄于4个标量水听器阵列的主瓣宽度, 近似于5个标量水听器阵列的主瓣宽度, 因此, 在具有相同的主瓣宽度下矢量水听器减少了阵列的孔径, 为低频信号探测的进一步研究奠定了基础。在下文的方位估计性能分析中, 以2个矢量水听器分辨2个目标为例, 对所提ISCMFA算法与CBF算法^[1]、MVDR算法^[2]、MUSIC算法^[2]、IAA算法^[17]的方位估计性能进行研究。

3.3 分辨概率分析

在声源的来波方向、信噪比、快拍数与图 3相同的条件下, 图 5给出了分辨概率与角度间隔的关系曲线。从图 5可以看出, CBF算法无法分辨角度小于30°的2个非相干目标, MUSIC算法与本文所提算法都具有较高的分辨概率, MUSIC算法通过信号子空间与噪声子空间的正交性实现了高分辨, 而本文所提算法利用最小二乘使信号协方差矩阵与估计的信号协方差矩阵的残差最小, 减少了噪声协方差矩阵对稀疏信号功率估计值的影响, 获得较高的分辨概率。

3.4 估计精度分析

从图 5可以看出, 在角度间隔Δ=30°时, 所有算法都可达到100%的分辨率。为此, 固定θ₁=-15°, θ₂=15°, 研究几种算法的均方根误差与相关系数的关系曲线。图 6给出了信噪比为15 dB, 快拍数为500时, 均方根误差与相关系数的关系曲线。从图 6可以看出, 当相关系数为0时, 所提方法估计的平均误差为MUSIC方法的39.4%;当相关系数为1时, 所提方法估计的平均误差为IAA方法的12.9%;相较于相干系数为0时, 所提方法在相干系数为1时的估计误差有所增加, 但是相较于CBF、MVDR、MUSIC与IAA方法, 所提方法的估计精度仍是最高的, 这归因于所提算法通过每次迭代估计的稀疏信号功率对信号协方差矩阵重构使其逼近理想协方差矩阵, 使得估计的信号功率更精确, 对网格点上的信号功率进行谱峰搜索时会获得更加精确的DOA估计。为此, 所提算法在估计具有由高度相关性目标的DOA时具有显著的优势。

从图 6可以看出, 当相关系数为1时, MUSIC算法方位估计精度急剧下降。为了研究几种算法与信噪比的关系, 图 7给出了在快拍数为500, 相关系数为0时, 均方根误差与信噪比的关系曲线。从图 7可以看出, 在所有信噪比范围内, CBF算法、IAA算法的方位估计精度都存在较大误差, MVDR算法与MUSIC算法在较高信噪比下才可获得较好的方位估计结果。所提算法在0~7dB范围内方位估计性能较于MVDR算法与MUSIC算法具有明显的优势; 在信噪比大于7 dB时, 所提算法较于MVDR算法和MUSIC算法方位估计性能约分别提高了9和3 dB。

图 8给出了信噪比为10 dB, 快拍数与均方根误差的关系曲线。仿真条件与图 7相同。从图 8可以看出, 在快拍数较小时, CBF、MVDR、MUSIC方法的均方根误差变化率较大; IAA方法与所提方法的均方根误差的变化率较小, 并且所提方法的均方根误差要小于IAA方法; 当快拍数较大时, CBF、MVDR及IAA方法的均方根误差变化极小, 而MUSIC方法与所提方法的均方根误差会随着快拍数的增加均方根误差逐渐减少, 并且MUSIC方法在快拍数较大时均方根误差的变化率要明显小于快拍数较少时均方根误差的变化率, 这主要是由于在快拍数较小时, 使得样本协方差矩阵的大特征值与小特征值不能明显区分开来, 错误的子空间划分导致方位估计出现较大的误差。所提算法相较于其它算法在所有快拍数范围内具有更精确的方位估计, 这是由于所提算法对利用Frobenius范数求解的稀疏信号功率进行迭代求解, 使其估计的稀疏信号功率更精确从而获得精确的DOA估计。

4 结论

为提高矢量水听器阵列相干信号方位估计性能, 本文提出了一种稀疏协方差矩阵拟合波达方向估计方法。该方法首先利用特征值分解获得信号协方差矩阵, 然后基于加权协方差矩阵拟合准则, 构建了关于信号协方差矩阵与稀疏信号功率的目标函数, 利用Frobenius范数的性质推导了稀疏信号功率迭代更新的递推式, 通过迭代估计得到稀疏信号功率谱, 最终通过峰搜索来获得DOA的估计值。理论分析表明, 所提算法对信号的相关性具有较低的敏感度。仿真实验结果证明, 在相干源情况下, 所提方法相较于IAA方法具有较高的估计精度, 其估计误差是IAA方法的12.9%;在非相干情况下, 所提算法估计性能优越于CBF、MVDR、MUSIC及IAA方法, 其估计误差分别是高分辨MUSIC算法与IAA算法的39.4%与73.7%。因此, 所提方法为实际应用中面临目标具有高度相关性的研究提供了理论基础, 下一步将通过水池与湖试实验数据验证算法的有效性。