基于声压振速联合处理的稀疏协方差DOA估计

Sparse Covariance DOA Estimation Based on Combined Information Processing of Sound Pressure and Particle Velocity

摘要

为充分利用矢量水听器中声压振速信息之间的关系来提高DOA估计精度， 本文提出了基于声压振速联合处理的稀疏协方差DOA（Direction of Arrival）估计方法。该方法首先利用声压振速之间的相关性， 构造阵列协方差矩阵； 其次， 将空间入射角度集合进行等角度划分， 构造超完备冗余字典； 然后， 在过完备基上寻找阵列协方差矩阵的最稀疏系数， 利用系数向量中的非零行所对应的行号得到DOA估计值。将该算法与CBF算法及L1-SVD算法进行对比仿真实验， 结果表明， 在信号源数分别为3， 4， 5的情形下， 本文所提算法在低信噪比和小快拍数情形时， 具有更低的均方根误差， DOA估计性能优势明显。

Abstract

In order to make full use of the relationship between sound pressure and particle velocity in vector hydrophone， and improve the accuracy of DOA estimation， sparse covariance DOA （Direction of Arrival） estimation method based on combined information processing of pressure and particle velocity was proposed. Firstly， the correlation between sound pressure and particle vibration velocity was used to construct the array of covariance matrix. Secondly， the spatial incident angle set was divided into equal angles and the super-complete redundant dictionary was built. Then looking for the most sparse coefficient of array covariance matrix on the over-complete basis and DOA estimates were obtained by using the row numbers corresponding to the non-zero rows in the coefficient vector. Finally， the algorithm was compared with CBF algorithm and L1-SVD algorithm， simulation results show that the proposed algorithm has lower root-mean-square error （RMSE） error when the number of signal sources is 3，4，5， respectively， the method has good DOA estimation performance compared to other algorithms under the conditions of small snapshot number and low signal-to-noise ratio （SNR）.

Graphical abstract

关键词

DOA估计 / 稀疏表示 / 阵列协方差矩阵 / 矢量线性阵

Key words

DOA estimation / sparse representation / arrray covariance matrix / vector sensor line array

0 引　言

波达方向估计（Direction of Arrival， DOA）是阵列信号处理的主要组成部分， 已被广泛应用到雷达预警、 移动通信、 电子对抗、 天文观测等领域。常规波束形成算法^［1］（CBF）为最早应用于阵列信号处理的算法， 康春玉等^［2］在CBF的基础上进行改进， 用平均时间延迟协方差矩阵来代替传感器接收数据的协方差矩阵进行奇异值分解， 改进后的算法在仿真实验和海上实测数据中都能实现良好的DOA估计性能。随后在CBF算法的基础上涌现出了大量的研究成果， 其中较为突出的有子空间分解类算法和子空间拟合类算法， 子空间分解类算法中的典型代表算法为以多重信号分类^［3］（MUltiple SIgnal Classification， MUSIC）算法为代表的噪声子空间类算法和以旋转不变子空间^［4］ （Estimation of Signal Parameters via Rotation Invariant Technique， ESPRIT）为代表的信号子空间类算法， 此类算法的共同特点为对接收到的信号进行特征分解， 将其分解为噪声子空间和信号子空间， 利用两个空间的正交性来进行求解， 此方法大大提高了分辨率， 子空间拟合类算法中最具代表性的为极大似然估计类算法^［5］， 但是最大似然估计算法由于在求解DOA问题时， 需要进行多维搜索， 运算量巨大。

近年来， 出现了基于信号稀疏表示（Sparse Signal Representation， SSR）的DOA估计方法，该方法利用信源在空域中的稀疏性， 在采样率远低于奈奎斯特采样的前提下， 对数据进行优化， 提高了DOA估计精度。Malioutov等^［6］提出了L1-SVD（L1-Singular Value Decomposition）方法， 用左奇异向量稀疏表示信号， 并基于L1范数对信号重构， 从而获得信源角度。窦慧晶等^［7］对L1-SVD算法进行改进， 提出了基于加权l_{\mathrm{1}}范数稀疏信号表示的DOA估计方法， 利用子空间的正交性得到一个加权矩阵， 然后基于L1范数对信号重构， 这种方法在低信噪比下能准确稳定地估计出信源角度。

由于声矢量传感器（AVS）由声压传感器和振速传感器复合而成， 使其相对于声压传感器可接收更多的信息， 因此， 声矢量传感器近年来成为众多学者研究的重点。考虑到声矢量传感器的优势， 研究者们基于矢量阵提出了许多高分辨率DOA估计算法^［8‐10］，但是这些算法都没有利用声压和振速之间相干的性能来求解DOA， 忽略了由此带来的抗各向同性噪声优势。对此， 许多学者开始研究声压振速联合处理的方法， 白兴宇等^［11］将AVS的抗各向同性噪声能力与相干信号子空间（CSS）方法的高性能宽带测向能力进行结合，对信噪比低的信号实现了高精度信源数检测与高分辨方位估计。李海涛等^［12］利用(p+v_c)v_c构造矩阵协方差， 提出了一种基于声压振速联合处理的矢量阵时域解析MVDR波束形成方法， 能将矢量阵时域解析MVDR的高分辨能力与声矢量阵的抗噪能力有机结合， 与常规矢量阵的时域解析方法相比， 降低了运算量， 可以更好地抑制各向同性噪声， 实现有效的DOA估计。马林等^［13］提出了一种基于矢量传感器声压振速时反卷积联合处理输出的被动宽带信号检测方法，有效提高了低信噪比下DOA估计的精度。

本文利用声压振速之间相干的关系， 提出了一种基于声压振速联合处理的阵列协方差矩阵的DOA估计方法， 通过对振速进行投影， 与声压信息进行联合来求解协方差矩阵， 构造过完备基， 利用阵列协方差向量的稀疏性进行DOA估计。仿真实验表明， 相比现有的几种DOA估计算法相比， 本文方法的性能更好。

1 二维声矢量阵列信号模型

考虑浅海远场一般情况， 本文不考虑振速的z分量。假设由M个矢量水听器组成的均匀线阵， 阵元间距为d， 接收D(M>D)个远场窄带信号， 第D个信号的到达角为{\theta }_D， 则在t时刻， 第m个矢量水听器的输出为

x_m\left(t\right)={\displaystyle \sum _{i=\mathrm{1}}^D}{a}_m\left({\theta }_i\right){\left[\mathrm{1},\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_i,\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i\right]}^{\mathrm{T}}{s}_i\left(t\right)+n_m\left(t\right),

m=\mathrm{1,2},\cdots ,M， t=\mathrm{1,2},\cdots ,L,（1）

式中： a_m\left({\theta }_i\right)={\mathrm{e}}^{-\mathrm{j}\mathrm{2}\mathrm{\pi }d\left(m-\mathrm{1}\right)\frac{\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i}{\lambda }}为第m个阵元对于第i(i=\mathrm{1,2},\cdots ,D)个信号的声压响应系数； s_i为第i个信号的波前； n_m\left(t\right)为第m个阵元接收到的均值为\mathrm{0}， 方差为{\sigma }^{\mathrm{2}}的高斯白噪声； L为数据样本量， 将上述阵列输出模型写成矩阵形式， 则

\mathit{X}(t)={\displaystyle \sum _{i=\mathrm{1}}^D}{\mathit{a}}_I\otimes \mathit{a}({\theta }_i)s_i(t)+{\mathit{N}}_v(t)=

{\displaystyle \sum _{i=\mathrm{1}}^D}{\mathit{a}}_p({\theta }_i)\otimes \mathit{a}({\theta }_i)s_i(t)+{\mathit{N}}_v(t)=

{\mathit{A}}_v(\theta )\mathit{S}(t)+{\mathit{N}}_v(t),（2）

式中： {\mathit{a}}_p({\theta }_i)={\left[\mathrm{1},{\mathrm{e}}^{-j\mathrm{2}\pi d\frac{\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i}{\lambda }},\cdots ,{\mathrm{e}}^{-j\mathrm{2}\pi d(M-\mathrm{1})\frac{\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i}{\lambda }}\right]}^{\mathrm{T}}； \mathit{a}({\theta }_i)={\left[\mathrm{1},\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_i,\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i\right]}^{\mathrm{T}}； {\mathit{a}}_I={[\mathrm{1,1},\cdots ,\mathrm{1}]}^{\mathrm{T}}为M维列向量； \mathrm{T}为矩阵的转置。

{\mathit{A}}_v\left(\theta \right)=\left[{\mathit{a}}_p\left({\theta }_{\mathrm{1}}\right)\otimes \mathit{a}\left({\theta }_{\mathrm{1}}\right),\right.{\mathit{a}}_p\left({\theta }_{\mathrm{2}}\right)\otimes \mathit{a}\left({\theta }_{\mathrm{2}}\right),\cdots ,\left.{\mathit{a}}_p\left({\theta }_D\right)\otimes \mathit{a}\left({\theta }_D\right)\right]， {\mathit{A}}_v\left(\theta \right)\in C^{\mathrm{3}M\times D}称为声矢量水听器的信号方向矩阵。符号\otimes表示矩阵克罗内克积或直积； \mathit{S}\left(t\right)={\left[s_{\mathrm{1}}\left(t\right),s_{\mathrm{2}}\left(t\right),\cdots ,s_D\left(t\right)\right]}^{\mathrm{T}}，\mathit{S}\left(t\right)∈C^{D\times L}为信号源矢量， 符号x_{mp}(t),x_{m{v}_x}(t),x_{m{v}_y}(t)分别为t时刻第m个阵元接收到的声压信息和振速在x、 y路2个分量的时域波形。该阵列接收到的数据为

\mathit{X}(t)=[x_{\mathrm{1}p}(t),x_{\mathrm{1}{v}_x}(t),x_{\mathrm{1}{v}_y}(t),\cdots ,x_{Mp},x_{M{v}_x}(t),x_{M{v}_y}{(t)]}^{\mathrm{T}}， \mathit{X}(t)\in C^{\mathrm{3}M\times L}，

{\mathit{N}}_v(t)=[n_{\mathrm{1}}(t),n_{\mathrm{2}}(t),\cdots ,n_{\mathrm{3}M}{(t)]}^{\mathrm{T}}，

{\mathit{N}}_v(t)\in C^{\mathrm{3}M\times L}为阵列接收到的噪声矢量。

2 声压振速联合处理模型

本文仅考虑二维情形， 在声矢量阵列信号输出模型的基础上， 将接收到的声压振速信息分别表示为

{\mathit{X}}_p(t)={\displaystyle \sum _{i=\mathrm{1}}^D}{\mathit{a}}_p(\theta )s_i(t)+{\mathit{N}}_p(t)=

{\mathit{A}}_p(\theta )\mathit{S}(t)+{\mathit{N}}_p(t),（3）

{\mathit{X}}_{v_x}(t)={\displaystyle \sum _{i=\mathrm{1}}^D}{\mathit{a}}_p(\theta )s_i(t)\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_i+{\mathit{N}}_{v_x}(t)=

{\mathit{A}}_p(\theta ){\theta }_{v_x}\mathit{S}(t)+{\mathit{N}}_{v_x}(t),（4）

{\mathit{X}}_{v_y}(t)={\displaystyle \sum _{i=\mathrm{1}}^D}{\mathit{a}}_p(\theta )s_i(t)\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i+{\mathit{N}}_{v_y}(t)=

{\mathit{A}}_p(\theta ){\theta }_{v_y}\mathit{S}(t)+{\mathit{N}}_{v_y}(t),（5）

其中，

{\mathit{X}}_p(t)=[x_{\mathrm{1}p}(t),x_{\mathrm{2}p}(t),\cdots ,x_{Mp}{(t)]}^{\mathrm{T}}，

{\mathit{X}}_{v_x}(t)=[x_{\mathrm{1}{v}_x}(t),x_{\mathrm{2}{v}_x}(t),\cdots ,x_{M{v}_x}{(t)]}^{\mathrm{T}}，

{\mathit{X}}_{v_y}(t)=[x_{\mathrm{1}{v}_y}(t),x_{\mathrm{2}{v}_y}(t),\cdots ,x_{M{v}_y}{(t)]}^{\mathrm{T}}，

A_p(\theta )=[a_p({\theta }_{\mathrm{1}}),a_p({\theta }_{\mathrm{2}}),\cdots ,a_p({\theta }_D)]，

式中： {\mathit{X}}_p(t),{\mathit{X}}_{v_x}(t),{\mathit{X}}_{v_y}(t)\in C^{M\times L}分别为矢量水听器的声压信息和振速在x、 y路的输出； {\mathit{N}}_p(t)， {\mathit{N}}_{v_x}(t)， {\mathit{N}}_{v_y}(t)\in C^{M\times L}分别为矢量水听器的声压、 振速x、 y路输出的噪声向量； {\theta }_{v_x}， {\theta }_{v_y}\in R^{D\times D}是对应于振速在x、 y路输出的系数。

{\theta }_{v_x}=diag[\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{1}},\cdots ,\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_D],（6）

{\theta }_{vy}=diag[\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{1}},\cdots ,\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_D]\mathrm{。}（7）

矢量阵两轴向振速输出信息可以进行组合， 使用电子旋转法得到振速在某观测方位角\psi上的投影， 则

{\mathit{X}}_{vr}(t)={\mathit{X}}_{v_x}(t)\mathrm{c}\mathrm{o}\mathrm{s}\psi +X_{v_y}(t)\mathrm{s}\mathrm{i}\mathrm{n}\psi =

{\displaystyle \sum _{i=\mathrm{1}}^D}{\mathit{a}}_p(\theta )s_i(t)\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_i\mathrm{c}\mathrm{o}\mathrm{s}\psi +

{\displaystyle \sum _{i=\mathrm{1}}^D}{\mathit{a}}_p(\theta )s_i(t)\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_i\mathrm{s}\mathrm{i}\mathrm{n}\psi +{\mathit{N}}_{vr}(t)=

{\displaystyle \sum _{i=\mathrm{1}}^D}{\mathit{a}}_p(\theta )s_i(t)\mathrm{c}\mathrm{o}\mathrm{s}(\psi -{\theta }_i)+{\mathit{N}}_{vr}(t)={\mathit{A}}_p(\theta ){\theta }_{vr}S(t)+{\mathit{N}}_{vr}(t),（8）

式中： {\mathit{X}}_{vr}(t)=[x_{\mathrm{1}}(t),x_{\mathrm{2}}(t),\cdots ,x_M{(t)]}^{\mathrm{T}}； {\mathit{N}}_{vr}(t)={\mathit{N}}_{v_x}(t)\mathrm{c}\mathrm{o}\mathrm{s}\psi +{\mathit{N}}_{v_y}(t)\mathrm{s}\mathrm{i}\mathrm{n}\psi； {\theta }_{vr}=diag[\mathrm{c}\mathrm{o}\mathrm{s}(\psi-{\theta }_{\mathrm{1}})，\cdots ,\mathrm{c}\mathrm{o}\mathrm{s}(\psi -{\theta }_D)]，{\mathit{X}}_{vr}(t)，{\mathit{N}}_{vr}(t)\in C^{M\times L}，{\theta }_{vr}∈R^{D\times D}。

基于以上声压振速联合处理模型， 构造协方差矩阵

{\mathit{B}}_r=E[（{\mathit{X}}_p(t)+{\mathit{X}}_{vr}）{\mathit{X}}_{vr}^{\mathrm{H}}(t)]=

E[{\mathit{A}}_p(\theta )\mathit{S}(t){\mathrm{S}}^{\mathrm{H}}(t){\theta }_{vr}{\mathit{A}}_p^{\mathrm{H}}(\theta )]+E[{\mathit{N}}_p(t){\mathit{N}}_{vr}^{\mathrm{H}}(t)]+

E[{\mathit{A}}_p(\theta ){\theta }_{vr}\mathit{S}(t){\mathit{S}}^{\mathrm{H}}(t){\theta }_{vr}{\mathit{A}}_p^{\mathrm{H}}(\theta )]+

E[{\mathit{N}}_{vr}(t){\mathit{N}}_{vr}^{\mathrm{H}}(t)]\mathrm{。}（9）

在各向同性噪声场景中， 阵元间距为半波长时， 不同阵元接收到的声压振速不相关， 即

E[{\mathit{N}}_p(t){\mathit{N}}_{vr}^{\mathrm{H}}(t)]=\mathrm{0}\mathrm{。}（10）

此时， 互协方差阵为

{\mathit{B}}_r={\mathit{A}}_p(\theta )(\mathrm{1}+{\theta }_{vr})E[\mathit{S}(t){\mathit{S}}^{\mathrm{H}}(t)]{\theta }_{vr}{\mathit{A}}_p^{\mathrm{H}}(\theta )+

E[{\mathit{N}}_{vr}{\mathit{N}}_{vr}^{\mathrm{H}}]\mathrm{。}（11）

令{\mathit{B}}_s=E[\mathit{S}(t){\mathit{S}}^{\mathrm{H}}(t)]{\theta }_{vr}，{\mathit{B}}_N=E[{\mathit{N}}_{vr}{\mathit{N}}_{vr}^{\mathrm{H}}]，则

{\mathit{B}}_r={\mathit{A}}_p(\theta )(\mathrm{1}+{\theta }_{vr}){\mathit{B}}_s{}_p^{\mathrm{H}}\mathit{A}(\theta )+{\mathit{B}}_N,

{\mathit{B}}_r\in C^{M\times M}\mathrm{。}（12）

3 DOA估计方法

将矩阵{\mathit{A}}_p(\theta )进行扩展， 形成过完备字典{\mathit{A}}_{\Omega }(\tilde{\theta })， 假设空间入射角度集合\Omega =\left\{{\tilde{\theta }}_{\mathrm{1}},{\tilde{\theta }}_{\mathrm{2}},\cdots ,{\tilde{\theta }}_W\right\}(W>D)中包含了所有可能的入射角度， {\left\{{\tilde{\theta }}_w\right\}}_{w=\mathrm{1}}^W是在空间域中采样的特定方向，例如从0°~180°以1°为间隔划分。可构造过完备冗余字典^［14］为

{\mathit{A}}_{\Omega }(\tilde{\theta })=[a_{\Omega }({\tilde{\theta }}_{\mathrm{1}}),\cdots ,a_{\Omega }({\tilde{\theta }}_w),\cdots ,a_{\Omega }({\tilde{\theta }}_W)]；

将{\mathit{B}}_r的第m列重构为

r_m=E[（{\mathit{X}}_p(t)+{\mathit{X}}_{vr}）X_{vrm}^{*}(t)]={\mathit{A}}_{\Omega }(\tilde{\theta }){\mathit{b}}_m+{\sigma }^{\mathrm{2}}{\mathit{e}}_m,

m=\mathrm{1,2},\cdots ,M,（13）

式中： {()}^{*}是复共轭运算； {\mathit{X}}_{vrm}(t)为{\mathit{X}}_{vr}(t)的第m列； {\mathit{A}}_{\Omega }(\tilde{\theta })是M\times W维阵列流形矩阵； {\mathit{b}}_m是W\times \mathrm{1}维用上述超完备基表示的系数向量； 误差项{\mathit{e}}_m为M\times \mathrm{1}维向量， 第m项为1， 其余为0。由于{\mathit{A}}_{\Omega }(\tilde{\theta })有一个非平凡的零空间， 因此， 有不唯一的{\mathit{b}}_m满足式（13）。如果{\left\{{\tilde{\theta }}_w\right\}}_{w=\mathrm{1}}^W足够密集， 那么{\left\{a_{\Omega }({\tilde{\theta }}_w)\right\}}_{w=\mathrm{1}}^W中的D个向量可以非常接近（甚至等于）{\left\{a_{\Omega }({\tilde{\theta }}_i)\right\}}_{i=\mathrm{1}}^D{}_{}， 所以， 一个理想的{\mathit{b}}_m应该是除了D个与这D个基向量相关的元素之外， 所有元素都为零的向量。将式（13）中的模型写成矩阵形式为

\mathit{R}={\mathit{A}}_{\Omega }(\tilde{\theta })\mathit{B}+{\sigma }^{\mathrm{2}}{\mathit{I}}_M\mathrm{。}（14）

显然， 理想的{\left\{{\mathit{b}}_m\right\}}_{m=\mathrm{1}}^M应该与\mathit{B}=[b_{\mathrm{1}},b_{\mathrm{2}},\dots ,b_M]具有相同的稀疏结构， 即每个理想的b_m中非零元素的位置序号相同， 出现在\mathit{B}的同一行。然后， 引入一个向量{\mathit{b}}^{\circ }=[b_{\mathrm{1}}^{\circ },b_{\mathrm{2}}^{\circ },\dots ,b_W^{\circ }{]}^{\mathrm{T}}， 第w个元素b_w^{\circ }等于\mathit{B}的第w行L_{\mathrm{2}}范数， 即：

b_w^{\circ }={‖{(B)}_{w•}‖}_{\mathrm{2}}={\left({\displaystyle \sum _{m=\mathrm{1}}^M}{b}_{wm}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}}\mathrm{。}（15）

由于每个理想的{\left\{{\mathit{b}}_m\right\}}_{m=\mathrm{1}}^M的稀疏结构相同， 故用一个稀疏结构相同的{\mathit{b}}^{\circ }来描述{\left\{{\mathit{b}}_m\right\}}_{m=\mathrm{1}}^M的稀疏性， 且可寻找一个足够稀疏的{\mathit{b}}^{\circ }使得{\left\{{\mathit{b}}_m\right\}}_{m=\mathrm{1}}^M尽可能拟合{\left\{r_m\right\}}_{m=\mathrm{1}}^M， 从而使\mathit{B}的一行中的所有元素同时趋于零或非零。因此， 关于式（14）的DOA估计可以转化为找到一个足够稀疏的{\mathit{b}}^{\circ }， 通过利用{\mathit{b}}^{\circ }的稀疏结构来确定DOA值。注意， 估计精度将受限于网格划分大小， 这种限制是空间连续体的离散逼近所固有的， 但可以通过网格细化来缓解^［15］。

将上述引入的{\mathit{b}}^{\circ }称为协方差向量的稀疏表示。找到恰当的稀疏度度量准则来约束{\mathit{b}}^{\circ }对DOA估计的精度是必要的。由于需要通过寻找最稀疏的{\mathit{b}}^{\circ }来恢复信号， 但L_{\mathrm{0}}范数最小化问题为NP-hard问题， 无法直接求解， 常用L_{\mathrm{1}}范数问题逼近， 这里用L_{\mathrm{1}}范数^［16］来衡量{\mathit{b}}^{\circ }的稀疏度， 其定义为{‖\mathit{z}‖}_{\mathrm{1}}=\sum _i\left|z_i\right|， 其中， z_i是\mathit{z}的第i个元素的模， 表示向量\mathit{z}中非零元素的绝对值之和。向量越稀疏， 其L_{\mathrm{1}}范数越小， 故可用向量的L_{\mathrm{1}}范数的这一特性来寻找一个足够稀疏的{\mathit{b}}^{\circ }。

基于式（14）的DOA估计可以表示为约束优化形式

\underset{B}{\mathrm{m}\mathrm{i}\mathrm{n}}{‖{\mathit{b}}^{\circ }‖}_{\mathrm{1}} s.t. \mathit{R}={\mathit{A}}_{\Omega }(\tilde{\theta })\mathit{B}+{\sigma }^{\mathrm{2}}{\mathit{I}}_M\mathrm{。}（16）

在实际中， 未知\mathit{R}可以通过\widehat{\mathit{R}}={\mathit{B}}_r=\mathit{R}+\Delta \mathit{R}来估计，\Delta R=\widehat{\mathit{R}}-\mathit{R}是估计误差， 在大样本条件下\Delta R的分布特性表明， 其向量化形式满足

vec(\Delta \mathit{R})\sim AsN\left(\mathrm{0},\frac{\mathrm{1}}{L}{\mathit{R}}^{\mathrm{T}}\otimes \mathit{R}\right)\mathrm{。}（17）

式中： vec()表示堆栈操作， 按顺序排列， 矩阵的列一列一列地排列； AsN(\mu ,{\sigma }^{\mathrm{2}})表示均值为\mu、 方差为{\sigma }^{\mathrm{2}}的渐近正态分布； \otimes为克罗内克积。用\widehat{\mathit{R}}代替\mathit{R}的直接后果是式（16）中的等式约束由于\Delta \mathit{R}而崩溃。在这种情况下， 将“数据”vec(\widehat{\mathit{R}})拟合到其数据模型最佳的加权最小二乘框架^［17］中， 即

\underset{B}{\mathrm{m}\mathrm{i}\mathrm{n}}vec[\widehat{\mathit{R}}-{\mathit{A}}_{\Omega }(\tilde{\theta })\mathit{B}-{\sigma }^{\mathrm{2}}{\mathit{I}}_M{]}^H{\mathit{G}}^{-\mathrm{1}}\times vec[\widehat{\mathit{R}}-{\mathit{A}}_{\Omega }(\tilde{\theta })\mathit{B}-{\sigma }^{\mathrm{2}}{\mathit{I}}_M]+\lambda {‖{\mathit{b}}^{\circ }‖}_{\mathrm{1}},（18）

等价于

\underset{B}{\mathrm{m}\mathrm{i}\mathrm{n}}{‖{\mathit{G}}^{-\frac{\mathrm{1}}{\mathrm{2}}}vec[\widehat{\mathit{R}}-{\mathit{A}}_{\Omega }(\tilde{\theta })\mathit{B}-{\sigma }^{\mathrm{2}}{\mathit{I}}_M]‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda {‖{\mathit{b}}^{\circ }‖}_{\mathrm{1}}\mathrm{。}（19）

权矩阵{\mathit{G}}^{-\mathrm{1}}是渐近协方差矩阵vec(\Delta \mathit{R})的逆， {\mathit{G}}^{-\frac{\mathrm{1}}{\mathrm{2}}}=\sqrt[]{L}{\mathit{R}}^{-\frac{T}{\mathrm{2}}}\otimes {\mathit{R}}^{-\frac{\mathrm{1}}{\mathrm{2}}}是{\mathit{G}}^{-\mathrm{1}}的厄米特平方根； 正则化参数\lambda用来平衡{\mathit{b}}^{\circ }的稀疏性， 一般情况下， 用与式（17）或式（18）相同或类似的公式很难直接确定\lambda， 下面将推导一个更易于处理的DOA估计公式， 给出在L足够大的情况下拟合误差\Delta \mathit{R}的判别依据。

从（16）式可以推断出

{\mathit{G}}^{-\frac{\mathrm{1}}{\mathrm{2}}}vec[\widehat{\mathit{R}}-{\mathit{A}}_{\Omega }(\tilde{\theta })\mathit{B}-{\sigma }^{\mathrm{2}}{\mathit{I}}_M]\sim AsN(\mathrm{0},{\mathit{I}}_{M^{\mathrm{2}}}),（20）

可以直接得到

{‖{\mathit{G}}^{-\frac{\mathrm{1}}{\mathrm{2}}}vec[\widehat{\mathit{R}}-{\mathit{A}}_{\Omega }(\tilde{\theta })\mathit{B}-{\sigma }^{\mathrm{2}}{\mathit{I}}_M]‖}_{\mathrm{2}}^{\mathrm{2}}\sim As{\chi }^{\mathrm{2}}(M^{\mathrm{2}})\mathrm{。}（21）

As{\chi }^{\mathrm{2}}(M^{\mathrm{2}})表示自由度为M^{\mathrm{2}}的渐进卡方分布， 在式（20）中引入参数\eta， 有

{‖{\mathit{G}}^{-\frac{\mathrm{1}}{\mathrm{2}}}vec[\widehat{\mathit{R}}-{\mathit{A}}_{\Omega }(\tilde{\theta })\mathit{B}-{\sigma }^{\mathrm{2}}{I}_M]‖}_{\mathrm{2}}^{\mathrm{2}}\le \eta \mathrm{。}（22）

p为概率， 通常来说， 设置p=\mathrm{0.999}来确定\eta的值是足够的， 这使得上述不等式对任意数组成立， 概率为0.999（几乎是一个确定的事件）。

结合上述推导和分析， 可以得到一个统计上鲁棒性更好且更易处理的DOA估计， 公式为

\underset{\mathit{B}}{\mathrm{m}\mathrm{i}\mathrm{n}}{‖{\mathit{b}}^{\circ }‖}_{\mathrm{1}} s.t.{‖{\mathit{G}}^{-\frac{\mathrm{1}}{\mathrm{2}}}vec(\widehat{\mathit{R}}-{\sigma }^{\mathrm{2}}{\mathit{I}}_M)-{\mathit{G}}^{-\frac{\mathrm{1}}{\mathrm{2}}}[{\mathit{I}}_M\otimes {\mathit{A}}_{\Omega }(\tilde{\theta })]vec(\mathit{B})‖}_{\mathrm{2}}\le

\sqrt[]{\eta }\mathrm{。}（23）

式（23）中的误差抑制准则是式（22）的平方根版本（变体）^［18］，{\mathit{G}}^{-\frac{\mathrm{1}}{\mathrm{2}}}=\sqrt[]{L}{\mathit{R}}^{-\frac{T}{\mathrm{2}}}\otimes {\mathit{R}}^{-\frac{\mathrm{1}}{\mathrm{2}}}， {\mathit{G}}^{-\frac{\mathrm{1}}{\mathrm{2}}}vec(\widehat{\mathit{R}}-{\sigma }^{\mathrm{2}}{\mathit{I}}_M)， {\mathit{G}}^{-\frac{\mathrm{1}}{\mathrm{2}}}[{\mathit{I}}_M\otimes {\mathit{A}}_{\Omega }(\tilde{\theta })]可以表示为vec[\sqrt[]{L}{\widehat{\mathit{R}}}^{-\frac{\mathrm{1}}{\mathrm{2}}}(\widehat{\mathit{R}}-{\sigma }^{\mathrm{2}}{\mathit{I}}_M){\widehat{\mathit{R}}}^{-\frac{\mathrm{1}}{\mathrm{2}}}]和\sqrt[]{L}{\widehat{\mathit{R}}}^{-\frac{T}{\mathrm{2}}}⊗[{\widehat{\mathit{R}}}^{-\frac{\mathrm{1}}{\mathrm{2}}}{\mathit{A}}_{\Omega }(\tilde{\theta })]， {\sigma }^{\mathrm{2}}的最大似然估计（用{\widehat{\sigma }}^{\mathrm{2}}表示）是对\mathit{R}的\mathrm{3}(M-D)个较小特征值求平均值得到的。

为了用标准SOC编程形式表述式（23）， 引入辅助变量c使{‖{\mathit{b}}^{\circ }‖}_{\mathrm{1}}\le c， 不等式（23）进一步等价地转化为一组SOC约束模型， 即

\underset{\mathit{B},\gamma ,c}{\mathrm{m}\mathrm{i}\mathrm{n}}c s.t. {\mathrm{1}}^T\mathit{\gamma }\le c，{‖{(\mathit{B})}_w‖}_{\mathrm{2}}\le {\gamma }_w，

{‖y-\varphi vec(\mathit{B})‖}_{\mathrm{2}}\le \sqrt[]{\eta }， w=\mathrm{1,2},\cdots ,W,（24）

式中： 1是W\times \mathrm{1}维的单位向量； \mathit{\gamma }是第q个元素为{\gamma }_q的向量； {‖{(\mathit{B})}_w‖}_{\mathrm{2}}为{\mathit{b}}^{\circ }的第w个元素； y={\widehat{\mathit{G}}}^{-\frac{\mathrm{1}}{\mathrm{2}}}vec(\widehat{\mathit{R}}-{\widehat{\sigma }}^{\mathrm{2}}{\mathit{I}}_M)和\varphi ={\widehat{\mathit{G}}}^{-\frac{\mathrm{1}}{\mathrm{2}}}[{\mathit{I}}_M\otimes {\mathit{A}}_{\Omega }(\tilde{\theta })]，注意， 在计算式（14）之前， 需要对\psi的每一列进行重正化处理， 相对于vec(\mathit{B})的非归一化过完备基对计算其正确的稀疏结构有不利影响。式（24）为二次约束优化问题， 可用CVX数值软件包进行求解。

4 仿真分析

为验证本文方法的可行性， 研究信源数、 信噪比、 快拍数对算法DOA估计精度的影响。本文采用计算机仿真方式分析验证， 将本文算法与L1-SVD算法、 CBF（波束形成）算法进行对比， 在以下实验中设置矢量均匀线阵阵元数为8， 阵元间距设置为半波长（0.5 m）， 空间入射角度的集合\Omega的范围设定为0°~180°， 各算法设置字典离散网格为{\mathrm{1}}^{\circ }， 网格数W=\mathrm{181}， 概率p=\mathrm{0.999}， 将算法的估计均方根误差（RMSE）定义为

E_{\mathrm{R}\mathrm{M}\mathrm{S}}(i)\underline{\Delta }\frac{\mathrm{1}}{ND}\sqrt[]{{\displaystyle \sum _{n=\mathrm{1}}^N}{\displaystyle \sum _{i=\mathrm{1}}^D}[({\widehat{\theta }}_i^n-{\theta }_i{)}^{\mathrm{2}}]},（25）

式中： {\widehat{\theta }}_i^n为第n次实验对第i个信号的DOA估计值； N为蒙特卡洛实验次数。

4.1 多信源下各算法空间谱图对比

为了初步探究信号来源的个数对算法求解DOA的影响， 本实验中设置各信源的角度间隔均为20°， 信号源输入为［30°，50°，70°］、 ［30°，50°，70°，90°］和［30°，50°，70°，90°，110°］， 信号频率分别为［1 000 Hz，1 200 Hz，1 400 Hz］、 ［1 000 Hz，1 200 Hz，1 400 Hz，1 600 Hz］和［1 000 Hz，1 200 Hz，1 400 Hz，1 600 Hz， 1 800 Hz］， 快拍数均为100， 信噪比为0 dB，在此设置下比较L1-SVD算法、 CBF算法与本文算法在3个信号源、 4个信号源和5个信号源的空间谱图， 如图 1 所示。

由图 1（a） 可以看出， 在3个信源情况下， 该算法可以准确求出30°， 50°， 70°， CBF算法会出现伪峰或分辨不清相邻角度， 将30°， 50°估计成40°， 仅能求出70°， 而L1-SVD算法正确估计到的角度只有30°。随着信源数的增加， 从图 1（b）、 图 1（c） 中可以看出， 在4个信源和5个信源时， 本文算法可以完成DOA估计， 求得角度30°， 50°， 70°， 90°和330°， 50°， 70°， 90°， 110°。在图 1（b） 中， CBF算法在信源数为4时， 30°， 50°依然分辨不清， 仅可以求出70°， 90°。L1-SVD算法无法准确估计角度， 在图 1（c） 中， CBF算法可以求出分辨出5个信号源， 且可以求出70°， 90°， 110°， 而L1-SVD算法在此设置下并不能正确完成DOA估计。综上， 对比以上两种算法， 本文所提算法在0 dB，快拍数为100的情况下， DOA估计性能良好。

4.2 不同信噪比下各算法估计误差的比较

本节探究在信源数为3， 4， 5时信噪比的大小对算法DOA估计性能的影响。由于4.1节实验表明， 在快拍数为100时， CBF算法和L1-SVD算法的DOA估计性能较差， 因此， 在该实验中设置快拍数为300， 信噪比以5　dB为步长从-5 dB增加到20 dB， 信号源输入与4.1节实验相同， 在每个信噪比下进行100次蒙特卡洛实验。图 2 展示了3个信号源、 4个信号源和5个信号源情况下不同算法的DOA估计性能随信噪比变化的曲线。

由图 2 可以看出， 在信源数为3和4的情况下， 信噪比为-5 dB时， 本文算法的均方根误差为10°， 而CBF算法和L1-SVD算法的均方根误差已大于{\mathrm{10}}^{}， 随着信噪比的增大， L1-SVD算法和CBF算法由于间隔较小， DOA估计误差并未明显减小， 而本文所提算法的DOA估计误差一直呈现下降的趋势。

当信噪比达到0 dB时， 本文所提算法以较低的均方根误差完成了DOA估计。在信源数为5 时，随着信源数变多， 信噪比为-5时， 精度比信源数为3和4时有所降低， 但均方根误差随着信噪比降低， 与CBF算法和L1-SVD算法相比， DOA估计误差较低， 性能更好。

4.3 不同快拍数下各算法估计误差的比较

本节探究在信源数为3， 4， 5时快拍数的大小对算法DOA估计性能的影响。实验中设置信噪比为0 dB， 快拍数以50个快拍为步长从50增加到500， 信号源输入与4.1节实验相同， 每个信噪比下进行100次蒙特卡洛实验。图 3 展示了3个信号源， 4个信号源和5个信号源情况下， 不同算法的DOA估计性能随快拍数的变化曲线。

由图 3 可以看出， 3个信源时， 快拍数从50增加到100时， 本文算法的均方根误差急剧下降， 随后趋于稳定， 其大小远低于CBF算法和L1-SVD的误差。在信源数为4和5时， 随着快拍数的增加， CBF算法和L1-SVD算法的均方根误差没有下降的趋势， 而本文所提算法的均方根误差随着快拍数的增加逐渐减小。综上， 多信源情况下， 在信噪比为0时， 当快拍数增加时， DOA估计精度提高， 且本文算法与L1-SVD算法和CBF算法相比具有较低的均方根误差。

5 结　论

本文主要讨论了在矢量水听器中基于声压振速联合处理的DOA估计方法。该方法将接收到的振速信息在某一角度上进行投影， 与声压信息相加后再与投影处理后的振速相加求和来构建协方差矩阵； 然后， 将阵列流形矩阵进行扩展， 得到稀疏矩阵， 与矩阵运算结合构造出一个稀疏向量， 利用协方差误差服从渐进正态分布特性， 运用概率论的知识来约束协方差矩阵的误差， 在二次约束规划问题的框架下， 用CVX工具箱进行求解， 完成DOA估计。仿真实验表明， 该方法在多信源情形下具有良好的DOA估计性能。

参考文献

[1] THAKUR D， BAGHEL V， TALLURI S R. Proximal gradient method based robust Capon beamforming against large DOA mismatch［J］. Frequenz， 2021， 75（7/8）： 259-266.

[2] 康春玉， 章新华， 吴清华.一种实用的盲波束形成和DOA估计方法［J］.数据采集与处理， 2009， 24（S1）： 79-83.

[3] KANG Chunyu， ZHANG Xinhua， WU Qinghua. Improved blind beamforming and direction of arrival（DOA） estimation method［J］. Journal of Data Acquisition and Processing， 2009， 24（S1）： 79-83. （in Chinese）

[4] WAGNER M， PARK Y， GERSTOFT P. Gridless DOA estimation and Root-MUSIC for non-uniform linear arrays［J］.IEEE Transactions on Signal Processing， 2021， 69： 2144-2157.

[5] XIAO M， DUAN Z， YANG Z. A weighted forward-backward spatial smoothing DOA estimation algorithm based on TLS-ESPRIT［J］. IEICE Transactions on Information and Systems， 2021， 104（6）： 881-884.

[6] JIANG Z M， ZHANG P， RIHAN M， et al. Maximum likelihood approach to DOA estimation using lens antenna array［J］. EURASIP Journal on Wireless Communications and Networking， 2019， 2019（1）： 1-7.

[7] MALIOUTOV D， CETIN M， WILLSKY A S.A sparse signal reconstruction perspective for source localization with sensor arrays［J］.IEEE Transactions on Signal Processing， 2005， 53（8）： 3010-3022.

[8] 窦慧晶， 高立菁， 朱子云.基于加权l₁范数稀疏信号表示的DOA估计［J］.北京工业大学学报， 2018， 44（10）： 1297-1302.

[9] DOU Huijing， GAO Lijing， ZHU Ziyun. DOA estimation based on weighted l₁ norm sparse signal representation［J］. Journal of Beijing University of Technology， 2018， 44（10）： 1297-1302.（in Chinese）

[10] WONG K T， ZOLTOWSKI M D. Self-initiating MUSIC-based direction finding in underwater acoustic particle velocity-field beamspace［J］. IEEE Journal. of Oceanic Engineering， 2000， 25（2）： 262-273.

[11] WONG K T， ZOLTOWSKI M D. Near-field azimuth and elevation angle estimation using a single vector hydrophone［J］.IEEE Trans Signal Procesing，2001， 49（11）： 2498-2510.

[12] CHEN H W， ZHAO J W. Coherent signal-subspace processing of acoustic vector sensor array for DOA estimation of wideband sources［J］.Signal Processing， 2005， 85（1）： 837-847.

[13] 白兴宇， 杨德森， 赵春晖.基于声压振速联合信息处理的声矢量阵相干信号子空间方法［J］.声学学报，2006（5）： 410-417.

[14] BAI Xingyu， YANG Desen， ZHAO Chunhui.The coherent signal-subspace method based on combined information processing of pressure and particle velocity using the acoustic vector sensor array［J］. Acta Acustica， 2006（5）： 410-417. （in Chinese）

[15] 李海涛， 王余， 王易川.基于声压振速联合处理的矢量线阵时域解析MVDR方法研究［J］.舰船科学技术， 2015， 37（12）： 110-113.

[16] LI Haitao， WANG Yu， WANG Yichuan， Acoustic vector sensor line array time-domain analysis MVDR method based on the combined information processing of pressure and particle velocity［J］. Ship Science and Technology， 2015， 37（12）： 110-113. （in Chinese）

[17] 马林， 王彪， 赵安邦， 等.基于声压振速时反卷积联合处理的宽带信号检测［J］.江苏科技大学学报（自然科学版）， 2021， 35（6）： 64-71.

[18] MA Lin， WANG Biao， ZHAO Anbang， et al Broadband source detection based on sound pressure and vibration velocity time reversal convolution joint processing［J］. Journal of Jiangsu University of Science and Technology （Natural Science Edition）， 2021， 35（6）： 64-71. （in Chinese）

[19] SHI S， LI， Y， YANG D.et al. DOA Estimation of coherent signals based on the sparse representation for acoustic vector-sensor arrays［J］.Circuits Syst Signal Process， 2020， 39： 3553-3573.

[20] 蒋留兵， 荣书伟， 车俐.基于网格部分细化的DOA估计方法［J］.太赫兹科学与电子信息学报， 2020， 18（5）： 786-792.

[21] JIANG Liubing， RONG Shuwei， CHE Li. Partially refined grid for DOA estimation［J］. Journal of Terahertz Science and Electronic Information Technology， 2020， 18（5）： 786-792. （in Chinese）

[22] 孙转转， 宋儒瑛. L1范数下有噪信号恢复的误差估计［J］.成都大学学报（自然科学版）， 2016， 35（4）： 345-347.

[23] SUN Zhuanzhuan， SONG Ruying. Error estimate of noisy signal recovery by L_1-norm［J］.Journal of Chengdu University（Natural Science Edition）， 2016， 35（4）： 345-347. （in Chinese）

[24] OTTERSTEN B， STOICA P， ROY R AND， Covariance matching estimation techniques for array signal processing applications［J］.Digit Signal Process， 1998（8）： 185-210.

[25] 苏龙， 谷绍湖.基于阵列协方差矩阵的DOA估计方法［J］.计量与测试技术， 2022， 49（12）： 7-10.

[26] SU Long， GU Shaohu. DOA estimation method based on array covariance matrix［J］. Metrology & Measurement Technique， 2022， 49（12）： 7-10.（in Chinese）

基金资助

国家自然科学基金项目(61774137)

山西省基础研究计划资助项目(202103021224195)

山西省基础研究计划资助项目(202103021224212)

山西省基础研究计划资助项目(202103021223189)

山西省基础研究计划资助项目(20210302123019)

山西省回国留学人员科研项目(2020-104)