基于深度学习的原子范数最小化无网格DOA估计方法

A Gridless DOA Estimation Method for Atomic Norm Minimization Based on Deep Learning

摘要

针对传统子空间及现有的一些稀疏重构算法在小快拍、 多信源、 相干信源条件下阵列信号DOA估计精度低等问题， 提出了一种基于深度学习的原子范数最小化无网格DOA估计方法。首先在Toeplitz矩阵范德蒙德分解定理的前提下， 利用矢量水听器声压振速相关性构造Toeplitz矩阵， 然后利用深度学习交替乘子算法重构Toeplitz矩阵， 最后利用Rootmusic算法分解Toeplitz矩阵求得DOA估计值。仿真实验表明该算法相较于其他传统算法和一些稀疏重构算法在小快拍、 多信源、 相干信源条件下具有更低的均方根误差， DOA估计性能优势明显。

Abstract

Addressing the issues of low DOA estimation accuracy of array signals under conditions of small snapshots， multiple sources， and coherent sources in traditional subspace and existing sparse reconstruction algorithms， a grid-less DOA estimation method is proposed based on deep learning and atomic norm minimization. Firstly， under the premise of the Toeplitz matrix Vandermonde decomposition theorem， the Toeplitz matrix is constructed by the correlation between pressure and velocity from the vector hydrophone. Then， the Toeplitz matrix is reconstructed using the deep learning-based alternating direction method of the multipliers algorithm. Finally， the Root-music algorithm is used to decompose the Toeplitz matrix and obtain the DOA estimates. Simulation experiments show that this algorithm has a lower root mean square error compared to other traditional algorithms and some sparse reconstruction algorithms under the above conditions， demonstrating a clear advantage in DOA estimation performance.

Graphical abstract

关键词

波达方向 / 范德蒙德分解 / 原子范数最小化 / 深度学习交替乘子算法

Key words

direction of arrival / Vandermonde decomposition / atomic norm minimization / depth alternating multiplier algorithm

0 引　言

波达方向（Direction of Arrival，DOA）估计广泛应用于航海、 军事、 航空等领域。早期的DOA估计以常规波束形成（Conventional Beamforming， CBF）算法^［1］为代表， 但是该算法受到相关瑞利限的约束， 因此算法估计精度受阵列孔径及阵元数的影响。随着子空间类算法的涌现， DOA估计精度得到了极大的提升， 其中一类是子空间分解算法， 典型代表有多重信号分类（Multiple Signal Classification， MUSIC）算法^［2］和旋转不变子空间（Estimation of Signal Parameters via Rotation Invariant Technique， ESPRIT）算法^［3］。其中， MUSIC算法是一种基于特征结构的DOA估计方法， 其通过谱峰搜索对应的空域网格划分角度间隔求得DOA估计值， 而ESPRIT算法则避免了MUSIC算法的谱峰搜索步骤， 直接从旋转因子中求解DOA， 算法的运算量进一步减小。这些算法充分利用了信号子空间与噪声子空间正交的特性， 突破了瑞利限的约束， 但其估计精度完全依赖于阵列接收信号的统计特性， 需要大量的快拍数。另一类是子空间拟合算法， 具有代表性的是最大似然（Maximum Likelihood， ML）算法^［4‐5］， 该类算法通过求解一组参变量来使得多次联合采集快拍数据的概率密度函数最小以达到DOA估计的目的， 但由于计算量大， 限制了实际应用。

近年来， 由D.Donoho和Candes等^［6‐7］提出的压缩感知（Compressed Sensing， CS）理论对信号的压缩采样提供了新的方法， 打破了奈奎斯特采样定理中信号采样速率要低于两倍信号带宽的条件限制， 在实现对信号采样的同时完成压缩的过程中并不直接测量信号本身， 它使用非自适应的线性投影来获得信号的整体构造， 从而直接得到重要信息， 忽略那些在有损压缩中会被丢弃的信息。这些研究学者的主要工作可以分为三部分： 首先是信号的稀疏表示， 该过程是压缩感知理论的基础和前提； 其次是投影测量， 现已证明大部分具有一致分布的随机矩阵都可以作为观测矩阵； 最后是重构算法， 该过程采用全局非自适应测量方法， 观测数量远少于信号长度， 从而使得数据量的采集大大减少。压缩感知理论指出， 只要信号是稀疏的或在某个变换域是稀疏的， 那么就可以通过一个稀疏阵将所得的高维信号投影成一个低维信号， 然后从低维信号中以较高的准确率重构出原始信号^［8］。

之后， 压缩感知在DOA估计中的应用蓬勃兴起， 其中压缩感知算法分为贪婪算法及凸优化算法。贪婪算法的代表是正交匹配追踪（Orthogonal Matching Pursuit，OMP）算法^［9］， 其通过构造完备字典集， 利用最小二乘思想不断更新原子集， 从而选取与阵列信号最相关的原子集来达到重构信号的目的， 这类算法运算量小， 重构效率高， 但其属于网格划分类DOA估计方法， 存在分辨力不足及收敛性无法保证等缺点。凸优化算法将非凸问题转化为凸问题， 从而得到原始信号的逼近解， 比较典型的算法有索套回归（Least Absolute Shrinkage and Selection Operator， LASSO）算法^［10］， 但这类算法往往计算速度慢， 而且求解的无约束优化问题需要合适的正则化参数来保证算法的重建性能^［11］。刘永豪等^［12］等提出了交替方向乘子法的矢量水听器DOA估计方案， 提高了算法效率； 郑文康等^［13］提出了基于可分离替代函数算法的DOA估计方法， 并应用于矢量水听器阵列信号估计， 相比于同类算法估计优势明显。矢量水听器由声压传感器及振速传感器组成， 只考虑二维情况下的矢量水听器， 其可以联合拾取声压和x\mathrm{、}y路的振速分量信息， 相较于声压传感器可以获得更多的信息， 将上算法应用于矢量水听器进行DOA估计时直接利用了矢量水听器的接收数据而未利用矢量水听器的声压振速相关性， 忽略了由此带来的各向同性噪声优势。为了进一步提高矢量水听器DOA估计精度， 许多学者开始研究有关声压振速的DOA高精度估计方法。孙国仓等^［14］ 提出了以声压振速互协方差矩阵为基础的ESPRIT频率估计算法， 在不损失估计性能的前提下， 计算量减小为现有算法的1/27； 姚直象等^［15］提出了一种基于声压振速联合处理的矢量阵旋转不变子空间估计方法， 该方法充分利用声压振速组合指向性抗干扰能力， 可以更好地抑制各向同性干扰， 提高阵列的处理增益， 从而得到更好的方位估计性能； 张旭等^［16］提出一种使用协方差矩阵分解（Covariance Matrix Decompose， CMD）的声矢量阵声压振速联合处理方法， 该方法使用重构的协方差矩阵能够获得更高的阵处理增益， 较传统声压振速联合处理方法具有更高的DOA估计精度。

综上所述， 为了进一步提高矢量水听器的阵列信号DOA估计精度， 本文在托普利兹矩阵（Toeplitz Matrix）\mathit{T}(\mathit{u})范德蒙德（Vandermonde ）分解定理的前提下， 提出了一种基于深度学习的原子范数最小化无网格DOA估计方法（ANM-ADMM-NET）， 通过矢量水听器声压振速相关性构建Toeplitz 矩阵\mathit{T}(\mathit{u})， 将原子范数最小化（Atomic Norm Minimization，ANM）问题转化为半正定（Semi Definite Programming， SDP）优化问题， 然后利用深度学习交替方向乘子（Alternating Driection Me-thod of Multipliers， ADMM）算法重构Toeplitz矩阵\mathit{T}(\mathit{u})， 最后通过求根多重信号分类（Root Multiple Signal Classification， Rootmusic）算法对Toeplitz矩阵\mathit{T}(\mathit{u})进行Vandermonde分解来求得DOA估计值。仿真实验表明该算法在小快拍、 多信源、 相干信源条件下均体现了较稳定的DOA估计性能， 而且始终具有较强的稳健性。

1 矢量水听器阵列声压振速处理模型

假设t时刻有K个远场窄带信号s_{\mathrm{1}}(t),s_{\mathrm{2}}(t),\cdots ,s_K(t)入射到由M个矢量水听器按间距d组成的均匀线阵中。并且第k个入射信号的水平方位角为{\theta }_k,k=\mathrm{1,2},\cdots ,K。如果仅考虑二维情况， 每个矢量水听器由1个声压传感器及2个振速传感器组成， 它可以联合拾取声压振速信息， 如图 1 所示， 其中v_x， v_y为入射信源声压部分沿x\mathrm{、}y轴方向的振速分量。

以均匀线阵阵列中第一个阵元为参考点， t时刻第m个声压水听器输出为

p_m(t)={\displaystyle \sum _{k=\mathrm{1}}^K}{a}_m({\theta }_k)s_k(t)+n_m(t),（1）

式中： a_m({\theta }_k)={\mathrm{e}}^{-\mathrm{j}\mathrm{2}\mathrm{\pi }d(m-\mathrm{1})\frac{\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k}{\lambda }}为t时刻第m个阵元对第k个入射信号的声压响应系数， \lambda为入射信号的波长； s_k(t)为t时刻第k个入射信号的波前； n_m(t)为t时刻第m个阵元接收的噪声。阵元间距d为半波长。

式（1）写成矢量形式为

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

式中： \mathit{P}(t)=[p_{\mathrm{1}}(t),p_{\mathrm{2}}(t),\cdots ,p_M{(t)]}^{\mathrm{T}}为t时刻声压阵列输出矢量， \mathrm{T}表示转置； \mathit{S}(t)=[s_{\mathrm{1}}(t),s_{\mathrm{2}}(t),\cdots ,s_K{(t)]}^{\mathrm{T}}为信号源矢量； {\mathit{A}}_p(\theta )=[{\mathit{a}}_p({\theta }_{\mathrm{1}}),{\mathit{a}}_p({\theta }_{\mathrm{2}}),\cdots ,{\mathit{a}}_p({\theta }_K)]为声压阵列信号方向矩阵， {\mathit{a}}_p({\theta }_k)=[a_{\mathrm{1}}({\theta }_k),a_{\mathrm{2}}({\theta }_k),\cdots ,a_M({\theta }_k{)]}^{\mathrm{T}}为声压阵对第k个信号的方向向量。

假设p_m(t)为t时刻入射信源声压部分， 则声压与振速分量有如下关系

\left\{\begin{array}{l}{p}_m(t)=x_m(t)+n_{pm}(t),\\ v_{xm}(t)=x_m(t)\mathrm{c}\mathrm{o}\mathrm{s}\theta +n_{vxm}(t),\\ v_{ym}(t)=x_m(t)\mathrm{s}\mathrm{i}\mathrm{n}\theta +n_{vym}(t),\end{array}\right.（3）

式中： n_p(t)\mathrm{、}{n}_{vxm}(t)\mathrm{、}{n}_{vym}(t)分别为噪声沿声压和振速方向的分量。于是t时刻第m个阵元的输出为

y_m(t)={\displaystyle \sum _{k=\mathrm{1}}^K}{a}_m({\theta }_k)[\mathrm{1},cos{\theta }_k,sin{\theta }_k{]}^T{s}_k(t)+n_{vm}(t),（4）

式中： n_{vm}(t)为第m个阵元接收到的噪声。

n_{vm}(t)=[n_{pm}(t),n_{vxm}(t),n_{vym}{(t)]}^{\mathrm{H}},（5）

式中： \mathrm{H}表示共轭转置。令u_k=[\mathrm{1},\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k,\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k{]}^{\mathrm{T}}， 则式（4）可写为

y_m(t)={\displaystyle \sum _{k=\mathrm{1}}^K}{a}_m({\theta }_k)u_k{s}_k(t)+n_{vm}(t)\mathrm{。}（6）

将该阵列的输出写成矩阵形式\mathit{Y}(t)=[y_{\mathrm{1}}(t),y_{\mathrm{2}}(t),\cdots ,y_M{(t)]}^{\mathrm{T}}， 则

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

式中： {\mathit{A}}_v(\theta )=[{\mathit{a}}_p({\theta }_{\mathrm{1}})\otimes u_{\mathrm{1}}，{\mathit{a}}_p({\theta }_{\mathrm{2}})\otimes u_{\mathrm{2}}，\cdots，{\mathit{a}}_p({\theta }_K)⊗u_K]为声矢量水听器阵列的信号方向矩阵， 符号\otimes为直积或克罗内克积； \mathit{S}(t)=[s_{\mathrm{1}}(t),s_{\mathrm{2}}(t),\cdots ,s_K{(t)]}^{\mathrm{T}}为信号源矢量； {\mathit{N}}_v(t)=[n_{v\mathrm{1}}(t),n_{v\mathrm{2}}(t),\cdots ,n_{v\mathrm{3}M}{(t)]}^{\mathrm{T}}为阵列接收噪声矢量， 其中n_{v\mathrm{3}M}(t)为t时刻第M个矢量水听器接收到的第3路声压振速信息。

一般处理以上矢量阵列信号的DOA估计都是基于Nehorai处理框架， 其基本思想是对矢量阵列输出的声压振速信息分别独立处理， 没有利用声压振速处理模型对声压振速信息进行联合处理。

下面将矢量阵列输出信号模型式（7）分别写成声压、 振速形式，

\left\{\begin{array}{l}{\mathit{Y}}_p(t)={\mathit{A}}_p(\theta )\mathit{S}(t)+{\mathit{N}}_p(t),\\ {\mathit{Y}}_{vx}(t)={\mathit{A}}_p(\theta ){\Phi }_{vx}\mathit{S}(t)+{\mathit{N}}_{vx}(t),\\ {\mathit{Y}}_{vy}(t)={\mathit{A}}_p(\theta ){\Phi }_{vy}\mathit{S}(t)+{\mathit{N}}_{vy}(t),\end{array}\right.（8）

式中： {\mathit{Y}}_p(t)\mathrm{、}{\mathit{Y}}_{vx}(t)和{\mathit{Y}}_{vy}(t)分别为矢量水听器的声压和振速x\mathrm{、}y路的输出； {\mathit{N}}_p(t)、 {\mathit{N}}_{vx}(t)、 {\mathit{N}}_{vy}(t)分别为矢量水听器的声压和振速x、 y输出的噪声向量； {\Phi }_{vx}、 {\Phi }_{vy}为x\mathrm{、}y路振速输出系数，

\begin{array}{c}{\Phi }_{\mathrm{v}\mathrm{x}}=diag[\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{\mathrm{1}},\cdots ,\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_K],\\ {\Phi }_{\mathrm{v}\mathrm{y}}=diag[\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{\mathrm{1}},\cdots ,\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_K],\end{array}（9）

式中： diag[·]表示以某一元素为对角线构成的对角矩阵。将上式中的振速投影到某观测方向\phi上， 可得

\begin{array}{c}{\mathit{Y}}_{vr}(t)={\mathit{Y}}_{vx}(t)\mathrm{c}\mathrm{o}\mathrm{s}\phi +{\mathit{Y}}_{vy}(t)\mathrm{s}\mathrm{i}\mathrm{n}\phi =\\ {\mathit{A}}_p(\theta ){\Phi }_{vr}\mathit{S}(t)+{\mathit{N}}_{vr}(t),\end{array}（10）

其中，

\begin{array}{c}{\Phi }_{vr}=diag[\mathrm{c}\mathrm{o}\mathrm{s}(\phi -{\theta }_{\mathrm{1}}),\cdots ,\mathrm{c}\mathrm{o}\mathrm{s}(\phi -{\theta }_K)],\\ {\mathit{N}}_{vr}(t)={\mathit{N}}_{vx}(t)\mathrm{c}\mathrm{o}\mathrm{s}\phi +{\mathit{N}}_{vy}(t)\mathrm{s}\mathrm{i}\mathrm{n}\phi \mathrm{。}\end{array}（11）

对于观测方向\phi的选择， 可以依据电子旋转方向确定， 或者指定为某一感兴趣方向。接着构造声压振速互协方差矩阵，

{\mathit{R}}_r=E[{\mathit{Y}}_p(t){\mathit{Y}}_{vr}^{\mathrm{H}}(t)]=

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

E[{\mathit{N}}_p(t){\mathit{N}}_{vr}^{\mathrm{H}}(t)],（12）

式中： E\left\{\cdot \right\}表示期望运算。在各向同性噪声场中， 共点同步测量的声压和振速是互不相关的。同时， 当阵元间距为半波长时， 空间上不同阵元所测得的声压和振速也互不相关， 即

\left\{\begin{array}{l}E[{\mathit{N}}_p(t){\mathit{N}}_{vx}^{\mathrm{H}}(t)]=\mathrm{0},\\ E[{\mathit{N}}_p(t){\mathit{N}}_{vy}^{\mathrm{H}}(t)]=\mathrm{0},\end{array}\right.（13）

所以，

\begin{array}{c}E[{\mathit{n}}_p(t){\mathit{n}}_{vr}^{\mathrm{H}}(t)]=\mathrm{c}\mathrm{o}\mathrm{s}\phi E[{\mathit{N}}_p(t){\mathit{N}}_{vr}^{\mathrm{H}}(t)]+\\ \mathrm{s}\mathrm{i}\mathrm{n}\phi E[{\mathit{N}}_p(t){\mathit{N}}_{vy}^{\mathrm{H}}(t)]=\mathrm{0}\mathrm{。}\end{array}（14）

将式（14）代入式（12）可得

{\mathit{R}}_r=E[{\mathit{Y}}_p(t){\mathit{Y}}_{vr}^{\mathrm{H}}(t)]={\mathit{A}}_p(\theta ){\mathit{R}}_S{\Phi }_{vr}{\mathit{A}}_p^{\mathrm{H}}(\theta ),（15）

式中： {\mathit{R}}_S=E[\mathit{s}(t){\mathit{s}}^{\mathrm{H}}(t)]为信源协方差矩阵。由于入射信号间互不相关， 故信源协方差矩阵为

{\mathit{R}}_S=diag[{\sigma }_{\mathrm{1}}^{\mathrm{2}},\cdots ,{\sigma }_K^{\mathrm{2}}],（16）

式中： {\sigma }_k^{\mathrm{2}}=E[{\mathit{S}}_k(t){\mathit{S}}_k^{\mathrm{H}}(t)]为入射信号功率。

令{\mathit{R}}_S^{\mathrm{\text{'}}}={\mathit{R}}_S{\Phi }_{vr}为等效信源协方差矩阵， 则

{\mathit{R}}_S^{\text{'}}=diag[{\sigma }_{\mathrm{1}}^{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}(\phi -{\theta }_{\mathrm{1}}),\cdots ,{\sigma }_N^{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}(\phi -{\theta }_K)]\mathrm{。}（17）

此时， 声压振速互协方差矩阵为

{\mathit{R}}_r={\mathit{A}}_p(\theta ){\mathit{R}}_S^{\mathrm{\text{'}}}{\mathit{A}}_p^{\mathrm{H}}(\theta ),（18）

{\mathit{R}}_r具有与Toeplitz矩阵^{［17‐18］}类似的Vandermonde 分解形式。

2 原子范数模型

由于上述矢量水听器沿振速x\mathrm{、}y路方向输出已具备声压方向信息， 故取输出信息{\mathit{Y}}_{vr}(t)进行以下理论分析。

\begin{array}{c}{\mathit{Y}}_{vr}(t)={\mathit{A}}_p(\theta ){\Phi }_{vr}\mathit{S}(t)+{\mathit{N}}_{vr}(t)=\\ {\displaystyle \sum _{k=\mathrm{1}}^K}{\mathit{a}}_p({\theta }_k)\mathrm{c}\mathrm{o}\mathrm{s}(\phi -{\theta }_k){\mathit{s}}_k(t)+{\mathit{N}}_{vr}(t)=\\ {\displaystyle \sum _{k=\mathrm{1}}^K}{c}_k{\mathit{a}}_p({\theta }_k)\mathrm{c}\mathrm{o}\mathrm{s}(\phi -{\theta }_k){\mathit{b}}_k+{\mathit{N}}_{vr}(t)=\\ \mathit{X}+{\mathit{N}}_{vr}(t),\end{array}（19）

式中： {\mathit{s}}_k(t)=[s_{k\mathrm{1}}(t),s_{k\mathrm{2}}(t),\cdots ,s_{kL}(t)]\in {\mathbb{C}}^{\mathrm{1}\times L}， {\mathit{b}}_k=c_k^{-\mathrm{1}}{s}_k(t)， c_k={‖{\mathit{s}}_k(t)‖}_{\mathrm{2}}>\mathrm{0}并且{‖{\mathit{b}}_k‖}_{\mathrm{2}}=\mathrm{1}， 其中， {‖·‖}_{\mathrm{2}}表示向量的2-范数， L为快拍数。依据原子范数理论^［19］， 信号\mathit{X}可以看作集合{\mathit{A}}_M中的K个原子的线性组合， 集合{\mathit{A}}_M可以被定义为

{\mathit{A}}_M\triangleq \left\{\mathit{A}({\mathit{a}}_p({\theta }_k)\mathrm{c}\mathrm{o}\mathrm{s}(\phi -{\theta }_k),{\mathit{b}}_k)=\right.

\left.{\mathit{a}}_p({\theta }_k)\mathrm{c}\mathrm{o}\mathrm{s}(\phi -{\theta }_k){\mathit{b}}_k\right\},（20）

则关于\mathit{X}的l_{\mathrm{0}}-原子范数被定义为

{‖\mathit{X}‖}_{{\mathit{A}}_{M,\mathrm{0}}}=\underset{K}{\mathrm{i}\mathrm{n}\mathrm{f}}\left\{K\left|\mathit{X}={\displaystyle \sum _{k=\mathrm{1}}^K}{c}_k\mathit{A}({\mathit{a}}_p({\theta }_k)\mathrm{c}\mathrm{o}\mathrm{s}(\phi -{\theta }_k),{\mathit{b}}_k),{\mathit{a}}_p({\theta }_k)\mathrm{c}\mathrm{o}\mathrm{s}(\phi -{\theta }_k)\in {\mathit{A}}_M,c_k\ge \mathrm{0}\right.\right\},（21）

因此， 式（21）等价于处理以下优化问题：

\underset{\mathit{X}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}{‖\mathit{X}‖}_{{\mathit{A}}_M,\mathrm{0}}},\mathrm{s}.\mathrm{t}.{‖{\mathit{Y}}_{vr}(t)-\mathit{X}‖}_{\mathrm{2}}^{\mathrm{2}}\le \epsilon \mathrm{。}（22）

通常上式问题的解可以由其对偶问题的解得到^{［20‐21］}， 利用声压振速互协方差矩阵{\mathit{R}}_r\triangleq \mathit{T}(\mathit{u})^［19］的半正定托普利兹性和低秩性， 且T(\cdot )表示从向量\mathit{u}\in {\mathbb{C}}^{M\times \mathrm{1}}到赫尔米特矩阵\mathit{T}(\mathit{u})的映射， 也可以将其表示为如下的SDP问题

\underset{\mathit{W},\mathit{u},\mathit{X}}{\mathrm{m}\mathrm{i}\mathrm{n}\frac{\mathrm{1}}{\mathrm{2}}}\left[\mathrm{T}\mathrm{r}(\mathit{W})+\mathrm{T}\mathrm{r}(\mathit{T}(u))\right],

\mathrm{s}.\mathrm{t}.\left[\begin{array}{l}\mathit{T}(u) \mathit{X}\\ {\mathit{X}}^{\mathrm{H}} \mathit{W}\end{array}\right]\ge \mathrm{0},{‖{\mathit{Y}}_{vr}(t)-\mathit{X}‖}_{\mathrm{F}}^{\mathrm{2}}\le \epsilon ,（23）

式中： {‖·‖}_{\mathrm{F}}表示矩阵的 F 范数； \mathrm{T}\mathrm{r}(\cdot )表示迹； \mathit{W}\in {\mathbb{C}}^{L\times L}为辅助变量矩阵。对于上述问题， CVX求解器SDPT3可以有效地保证求解最优\mathit{T}(\mathit{u})^［22］， 然后利用Vandermonde分解^{［20‐21］}\mathit{T}(\mathit{u})求得DOA， 但由于其计算量巨大， 不能满足实时性要求。

故将上述问题改写为

[\mathit{X},\mathit{u}]=\underset{\mathit{X},\mathit{W},\mathit{u},\mathit{\Theta }}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\tau }{\mathrm{2}}[\mathrm{T}\mathrm{r}(\mathit{W})+\mathrm{T}\mathrm{r}(\mathit{T}(\mathit{u}))]+\frac{\mathrm{1}}{\mathrm{2}}{‖{\mathit{Y}}_{vr}(t)-\mathit{X}‖}_{\mathrm{F}}^{\mathrm{2}},\mathrm{s}.\mathrm{t}.\mathit{\Theta }=\left[\begin{array}{l}\mathit{T}(\mathit{u}) \mathit{X}\\ {\mathit{X}}^{\mathrm{H}} \mathit{W}\end{array}\right]\ge \mathrm{0},（24）

式中： \tau为正则化因子， 上式的增广拉格朗日函数为^［23］

\underset{\mathit{X},\mathit{W},\mathit{u},\mathbf{\Theta }\ge \mathrm{0},\mathbf{\Lambda }}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\tau }{\mathrm{2}}\left[\mathrm{T}\mathrm{r}(\mathit{W})+\mathrm{T}\mathrm{r}(\mathit{T}(\mathit{u}))\right]+\frac{\mathrm{1}}{\mathrm{2}}{‖{\mathit{Y}}_{vr}(t)-\mathit{X}‖}_{\mathrm{F}}^{\mathrm{2}}+〈\mathit{\Lambda },\mathit{\Theta }-\left[\begin{array}{l}\mathit{T}(\mathit{u}) \mathit{X}\\ {\mathit{X}}^{\mathrm{H}} \mathit{W}\end{array}\right]〉+

\frac{\rho }{\mathrm{2}}{‖\mathit{\Theta }-\left[\begin{array}{l}\mathit{T}(\mathit{u}) \mathit{X}\\ {\mathit{X}}^{\mathrm{H}} \mathit{W}\end{array}\right]‖}_{\mathrm{F}}^{\mathrm{2}}=\underset{\mathit{X},\mathit{W},\mathit{u},\mathbf{\Theta }\ge \mathrm{0},\mathbf{\Lambda }}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\tau }{\mathrm{2}}[\mathrm{T}\mathrm{r}(\mathit{W})+\mathrm{T}\mathrm{r}(T(\mathit{u}))]+\frac{\mathrm{1}}{\mathrm{2}}{‖{\mathit{Y}}_{vr}(t)-\mathit{X}‖}_{\mathrm{2}}^{\mathrm{2}}+\Delta ,（25）

式中： \mathit{\Lambda }\in {\mathbb{C}}^{(M+L)\times (M+L)}为朗格拉日乘子； \rho >\mathrm{0}为惩罚因子。上式是一个无约束优化问题， 在获得接收信号{\mathit{Y}}_{vr}(t)的情况下， 仅利用{\mathit{Y}}_{vr}(t)所含的声压振速信息拟合协方差矩阵{\mathit{R}}_r，然后通过ADMM算法交替迭代最小化代价函数\mathit{X},\mathit{W},\mathit{u},\mathit{\Theta }和\mathit{\Lambda }来得到最优估计\mathit{u}^{［23‐24］}， 从而通过Rootmusic算法分解\mathit{T}(\mathit{u})求解DOA。

3 基于原子范数最小化的深度网络无网格DOA估计方法

使用ADMM算法求解DOA需要提前选择合适的参数， 包括惩罚因子\rho和正则化因子\tau， 这对实际应用是一个挑战。同时， 参数设置不当会降低ADMM算法的收敛速度和精度， 从而增加计算复杂度， 降低DOA的估计性能。即使通过理论分析和交叉验证方法可以选择合适的参数^［23］， 固定的参数设置也不能保证最优ADMM算法的收敛性。故本文提出一种深度神经网络ANM-ADMM-NET算法， 其在ADMM算法的迭代步骤上新增学习参数\eta， 最终通过训练构造的样本数据集选择最优参数， 算法具体步骤如下。

其中， z_i=\mathrm{e}\mathrm{x}\mathrm{p}(j{w}_i)， w_i=\mathrm{2}\mathrm{\pi }{f}_i， f_i为信号的中心频率， z_i为由Rootmusic依据Vandermonde分解\mathit{T}({\mathit{u}}^{(K)})^{［25‐26］}构造的多项式的根， \mathit{a}=T\mathrm{*}(\mathit{A})为从矩阵\mathit{A}\in {\mathbb{C}}^{M\times M}到向量\mathit{a}\in {\mathbb{C}}^{M\times \mathrm{1}}的映射， 且a_i=sum({\mathit{A}}_{p,q}\left|q-p+\mathrm{1}=i)\right.， {\mathit{e}}_{\mathrm{1}}\in {\mathbb{C}}^{M\times \mathrm{1}}为单位矩阵的第一列， g=\mathrm{1,2}\cdots ,M+N,{\left\{{\delta }_g\right\}}_{+}为所有小于0的元素都被设置成0， \mathit{\Gamma }=diag\left(\frac{\mathrm{1}}{M},\frac{\mathrm{1}}{M-\mathrm{1}}\right.，…，\left.\frac{\mathrm{1}}{M-(M-\mathrm{1})}\right)， \mathit{G}为一个正交矩阵， 满足\mathit{G}{\mathit{G}}^{\mathrm{T}}={\mathit{I}}_{M+L}。以下详细介绍ANM-ADMM-NET算法的4个组成部分： 网络结构、 数据集构建、 网络初始化和训练过程。

3.1 网络结构

根据ANM-ADMM-NET算法迭代的步骤， 对应于如图 2 所示的K层网络ANM-ADMM-NET。其输入为{\mathit{Y}}_{vr}、 {\mathit{\Theta }}^{(\mathrm{0})}={\mathrm{0}}_{M+L}、 {\mathit{\Lambda }}^{(\mathrm{0})}={\mathrm{0}}_{M+L}， 可学习的参数为\mathit{\Omega }={\left\{{\mathit{\Omega }}^{(k+\mathrm{1})}\right\}}_{k=o}^{K-\mathrm{1}}=\{{\rho }_{k+\mathrm{1}}，{\tau }_{k+\mathrm{1}}，{\eta }_{k+\mathrm{1}}{\}}_{k=\mathrm{0}}^{K-\mathrm{1}}， 其中， 新参数{\eta }_{k+\mathrm{1}}是为了进一步增强ANM-ADMM-NET的学习能力及性能所添加的。通过构建以下网络层迭代得到\mathit{T}({\mathit{u}}^{(K)})， 进而得到DOA的估计值。其中， k+\mathrm{1}层ANM-ADMM-NET的网络参数为

{\mathit{u}}^{(k+\mathrm{1})}=F_{k+\mathrm{1}}\{{\mathit{Y}}_{vr},{\mathit{X}}^{(k)},{\mathit{W}}^{(k)},{\mathit{u}}^{(k)},{\mathit{\Theta }}^{(k)},{\mathit{\Lambda }}^{(k)},{\mathit{\Omega }}^{(k+\mathrm{1})}\},（26）

式中： F_{k+\mathrm{1}}\{\cdot \}由5个迭代子层组成， 包括重建子层A^{(k+\mathrm{1})}， 辅助变量更新子层B^{(k+\mathrm{1})}， Toeplitz变换子层C^{(k+\mathrm{1})}， 非线性子层D^{(k+\mathrm{1})}， 乘数更新子层E^{(k+\mathrm{1})}； 分别对应于\mathit{X}\mathrm{、}\mathit{W}\mathrm{、}\mathit{u}\mathrm{、}\mathit{\Theta }\mathrm{、}\mathit{\Lambda }的更新迭代， 具体描述如图 2 所示。

A^{(k+\mathrm{1})}\mathrm{、}{B}^{(k+\mathrm{1})}\mathrm{、}{C}^{(k+\mathrm{1})}\mathrm{、}{D}^{(k+\mathrm{1})}\mathrm{、}{E}^{(k+\mathrm{1})}层的输出分别为{\mathit{X}}^{(k+\mathrm{1})}\mathrm{、}{\mathit{W}}^{(k+\mathrm{1})}\mathrm{、}{\mathit{u}}^{(k+\mathrm{1})}\mathrm{、}{\mathit{\Theta }}^{(k+\mathrm{1})}\mathrm{、}{\mathit{\Lambda }}^{(k+\mathrm{1})}， 各层的输入如下：

1） 重建子层A^{(k+\mathrm{1})}： 输入{\mathit{\Theta }}_{\mathit{X}}^{(k)}、 {\mathit{\Lambda }}_{\mathit{X}}^{(k)}及{\mathit{Y}}_{vr}；

2） 辅助变量更新子层B^{(k+\mathrm{1})}： 输入{\mathit{\Theta }}_{\mathit{W}}^{(k)}、 {\mathit{\Lambda }}_{\mathit{W}}^{(k)}；

3） Toeplitz变换子层C^{(k+\mathrm{1})}： 输入{\mathit{\Theta }}_{T(\mathit{u})}^{(k)}、 {\mathit{\Lambda }}_{T(\mathit{u})}^{(k)}；

4） 非线性子层D^{(k+\mathrm{1})}： 输入{\mathit{\Lambda }}^{(k)}、 {\mathit{X}}^{(k+\mathrm{1})}、 {\mathit{W}}^{(k+\mathrm{1})}、 {\mathit{u}}^{(k+\mathrm{1})}；

5） 乘法器更新子层E^{(k+\mathrm{1})}： 输入{\mathbf{\Lambda }}^{(k)}、 {\mathit{X}}^{(k+\mathrm{1})}、 {\mathit{W}}^{(k+\mathrm{1})}、 {\mathit{u}}^{(k+\mathrm{1})}、 {\mathit{\Theta }}^{(k+\mathrm{1})}。

考虑到每个子层的参数都是学习和调优的， 最后将产生K层的\mathrm{3}K个最优参数， 分别为\{{\rho }_{\mathrm{1}},{\rho }_{\mathrm{2}},\cdots ,{\rho }_k\}、 \{{\tau }_{\mathrm{1}},{\tau }_{\mathrm{2}},\cdots ,{\tau }_k\}、 \{{\eta }_{\mathrm{1}},{\eta }_{\mathrm{2}},\cdots ,{\eta }_k\}， 与ADMM算法中固定参数{\rho }_{\mathrm{1}},{\rho }_{\mathrm{2}},\cdots ,{\rho }_k， {\tau }_{\mathrm{1}},{\tau }_{\mathrm{2}},\cdots ,{\tau }_k的输入方式相比， 这种ANM-ADMM-NET参数学习策略具有优越的灵活性以及出色的非线性拟合能力^［27］。

3.2 数据集构造

本文提出的ANM-ADMM-NET是一种结合模型驱动和数据驱动的稀疏恢复方法， 其有效性关键在于构建一个良好泛化能力的数据集， 通过建立一个充分且完整的数据集， ANM-ADMM-NET在训练过程中能够减少过拟合现象， 从而更准确地进行DOA估计。因此， 本文生成训练信号数据集{\mathit{Y}}_{vr}和理想标签向量\mathit{u}， 具体为：

1） 给定矢量水听器数量M， 快拍数L；

2） 给定信源数最大值K， 并随机生成源数k；

3） 对每一个信号k(\mathrm{1,2},\cdots K)， 其频率服从均匀分布U(f_{\mathrm{m}\mathrm{i}\mathrm{n}},f_{\mathrm{m}\mathrm{a}\mathrm{x}})。在多源情况下任意两个源之间的频率间隔需要满足\underset{i\ne j}{\mathrm{m}\mathrm{i}\mathrm{n}}\left|f_i-f_j\right|>\mathrm{1}/M^［22］， 接收信号由式（7）产生， 训练验证数据由式（10）产生， 噪声为高斯随机白噪声；

4） 重复以上步骤可得到D组接收信号\{{\mathit{Y}}_d{\}}_{d=\mathrm{1}}^D， 理想标签集向量\{{\mathit{u}}_d{\}}_{d=\mathrm{1}}^D经对偶原子范数最小化通过Vandermonde分解得到^{［20‐21］}；

5） 将上述集合随机划分为训练数据集\{{\mathit{Y}}_q^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}},{\mathit{u}}_q^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}{\}}_{q=\mathrm{1}}^{Q=\mathrm{0.8}D}和验证数据集\{{\mathit{Y}}_o^{\mathrm{v}\mathrm{a}\mathrm{l}},{\mathit{u}}_o^{\mathrm{v}\mathrm{a}\mathrm{l}}{\}}_{o=\mathrm{1}}^{O=\mathrm{0.2}D}。

3.3 网络初始化与训练

3.3.1 网络初始化

考虑在式（23）中T(\cdot )映射的难度以及算法迭代中T^{\mathrm{*}}(\cdot )和特征分解的难度， 正确初始化参数\mathit{\Omega }=\{{\rho }_{k+\mathrm{1}},{\tau }_{k+\mathrm{1}},{\eta }_{k+\mathrm{1}}{\}}_{k=\mathrm{0}}^{K-\mathrm{1}}可以使其更容易收敛， 而且在一定程度上避免陷入局部最优解。每层参数的初始值设为{\rho }_{\mathrm{1}:K}={\rho }_{\mathrm{0}}、 {\tau }_{\mathrm{1}:K}={\tau }_{\mathrm{0}}、 {\eta }_{\mathrm{1}:K}={\eta }_{\mathrm{0}}， 在理论分析的基础上增强所提出方法的灵活性^［23］。

3.3.2 网络训练

采用Adam算法对参数进行学习和调整， 初始学习率为0.1， 以便快速实现可能的全局最优。基于3.2节中所构建的训练数据集\{{\mathit{Y}}_q^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}},{\mathit{u}}_q^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}{\}}_{q=\mathrm{1}}^{Q=\mathrm{0.8}D}， 网络层K的最优参数{\mathit{\Omega }}^{\mathrm{*}}=\{{\rho }_{k+\mathrm{1}}^{\mathrm{*}},{{\tau }_{k+\mathrm{1}}}^{\mathrm{*}},{\eta }_{k+\mathrm{1}}^{\mathrm{*}}{\}}_{k=\mathrm{0}}^{K-\mathrm{1}}可以通过使用反向传播（Back Proagation，BP）原理最小化以下归一化均方误差损失函数^［28］， 即

{\mathit{\Omega }}^{\mathrm{*}}=arg\mathrm{m}\mathrm{i}\mathrm{n}\left({\displaystyle \sum _{q=\mathrm{1}}^Q}{L}_{\mathit{u}}^q\right)/Q,（27）

式中： L_u^q={‖{\mathit{u}}^{(K)}(\mathit{\Omega },{\mathit{\Theta }}^{(\mathrm{0})},{\mathit{\Lambda }}^{(\mathrm{0})},{\mathit{Y}}_q^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}})-{\mathit{u}}_q^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}‖}_{\mathrm{F}}^{\mathrm{2}}/{‖{\mathit{u}}_q^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}‖}_{\mathrm{F}}^{\mathrm{2}}， {\mathit{u}}^{(K)}(\mathit{\Omega }，{\mathit{\Theta }}^{(\mathrm{0})}，{\mathit{\Lambda }}^{(\mathrm{0})}，{\mathit{Y}}_q^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}})为网络参数\mathit{\Omega }的第K层迭代Toeplitz变换子层输出的估计分量； {\mathit{\Theta }}^{(\mathrm{0})}={\mathrm{0}}_{M+L}、 {\mathit{\Lambda }}^{(\mathrm{0})}={\mathrm{0}}_{M+L}、 {\mathit{Y}}_q^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}作为其输入。根据验证数据集\{{\mathit{Y}}_o^{\mathrm{v}\mathrm{a}\mathrm{l}},u_o^{\mathrm{v}\mathrm{a}\mathrm{l}}{\}}_{o=\mathrm{1}}^{O=\mathrm{0.2}D}得到最优参数{\mathit{\Omega }}^{\mathrm{*}}后， 将其代入后续实验测试数据{\mathit{Y}}^{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}中可得最优{\mathit{u}}^{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}估计为

{\mathit{u}}^{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}={\mathit{u}}^{(K)}\{{\mathit{\Omega }}^{\mathrm{*}},{\mathit{\Theta }}^{(\mathrm{0})},{\mathit{\Lambda }}^{(\mathrm{0})},{\mathit{Y}}^{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}\},（28）

最后通过Rootmusic算法Vandermonde分解Toeplitz矩阵\mathit{T}({\mathit{u}}^{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}})求得DOA的估计值。

4 仿真实验与结果分析

4.1 不同网络层和快拍数的模型训练与验证损失

为了选取本文所提算法的最佳实验参数， 并验证DOA估计性能， 数据集构造中选取矢量水听器均匀线阵阵元数为12， 分别设计以下两个实验验证所构建网络层的训练验证损失性能： ① 固定快拍数为100， 网络层数K分别选取5、 10、 20； ② 固定网络层数K=\mathrm{20}， 并分别选取快拍数L=\mathrm{10}\mathrm{、}\mathrm{50}\mathrm{、}\mathrm{100}。在以上实验条件下分别随机产生1、 2、 3个信源角度， 每个信源角度产生300个样本数据集D， 在每个样本数据集中随机加入\left[-\mathrm{5,20}\right]之间的随机整数噪声， 训练周期Epoch为180， 不同的网络层数与不同的快拍数产生的训练验证损失如图 3、 图 4 所示。

从图 3 可以看出， 在固定样本数据集L=\mathrm{100}的情况下， 随着网络层数K的增加， 算法的训练损失与验证损失逐渐减小， 说明随着网络层的加深， 算法的参数学习能力也随之加强。在图4固定网络层数K的情况下， 随着训练样本集快拍数的增加， ANM-ADMM-NET的训练损失与验证损失达到了最低。进一步说明了在样本数据充足与网络层足够深的前提下， 可以获得理想的最优参数。为了测试在上述样本数据集情况下所选的最优参数在估计DOA时的性能优势， 选取L=\mathrm{100},K=\mathrm{20}的网络层的最优参数代入所提算法并进行后续实验， 以验证本文所提方法的可行性。

4.2 不同信噪比下各算法均方根误差对比实验

为了验证所提算法的性能， 选取4.1节训练所得最优参数进行以下实验。在本实验中设置30个矢量水听器的均匀线性阵列， 阵列中每个水听器的间距为0.5 m， 各算法快拍数设置为10， 实验信噪比范围为［-5 dB， 20 dB］， 步长为2 dB。4个信源角度为-\mathrm{50}°、 -\mathrm{40}°、 -\mathrm{20}°、 \mathrm{10}°， 频率为1 000、 1 200、 1 400、 1 600 Hz； 5个信源角度为-\mathrm{50}°、 -\mathrm{40}°、 -\mathrm{20}°、 \mathrm{10}°、 \mathrm{20}°， 频率为1 000、 1 200、 1 400、 1 600、 1 800 Hz； 6个信源角度为-\mathrm{50}°、 -\mathrm{40}°、 -\mathrm{20}°、 \mathrm{10}°、 \mathrm{20}°、 \mathrm{30}°， 频率为1 000、 1 200、 1 400、 1 600、 1 800、 2 000 Hz。均方根误差（Root Mean Square Error， RMSE）作为衡量各算法的DOA估计性能指标之一， 在每个信噪比下进行100次Monte Carlo实验。RMSE定义为

RMSE=\sqrt[]{\frac{\mathrm{1}}{KJ}{\displaystyle \sum _{k=\mathrm{1}}^K}{\displaystyle \sum _{j=\mathrm{1}}^J}[{\stackrel{\wedge }{\theta }}_k(j)-{\theta }_k{]}^{\mathrm{2}}} ,（29）

式中： K为信源总数； J为Monte Carlo实验次数； {\stackrel{\wedge }{\theta }}_k(j)为第j次Monte Carlo实验中第k个信源的DOA估计值； {\theta }_k为第k个入射信源的真实值。实验针对OMP、 MUSIC、 ESPRIT、 CBF、 LASSO、 ANM算法进行了比较， 图 4（a）~图 4（c） 分别展示了不同信源情况下均方根误差与信噪比变化的曲线。

从图 5 可以看出在多信源的条件下， 随着信噪比的增加， 除ESPRIT算法外， 其余算法的均方根误差均出现了递减趋势。

SNR在小于0 dB时各算法RMSE均出现了较大的递减趋势， 在0~5 dB之间也出现了一定的递减趋势， 在5 dB之后各算法的RMSE处于相对稳定的状态。而ESPRIT算法在4、 5信源的情况下， 随着信噪比的增加算法的均方根误差逐渐减小， 特别在4信源的情况下， 该算法在10 dB之后出现了较大的递减趋势， 但随着信源数的增加该算法的均方根误差递减趋势逐渐变缓， 在6信源的条件下出现了较大的估计误差。对于ANM方法求解DOA， 虽然该算法相较于其他算法具有一定的性能优势， 其在5 dB后相较于对比实验中的大部分算法均具有较低的均方误差， 但由于该算法是在固定参数的情况下进行DOA估计， 所以对于接收到的大量声压振速数据不具有泛化性， 因此实际应用还存在一定的局限性。相比之下， 本文所提算法在多信源条件下均表现出较小的均方误差， 即使在0 dB以下的低信噪比环境下也具有良好的估计精度。此外， 随着信源数量的增加， 该算法的估计误差并未受到显著影响， 这表明所构建的网络层在一定程度上对样本数据集中加入的随机噪声具有抗噪作用， 进一步验证了训练数据集构造的合理性。所选的最优参数因此具有一定的实用价值， 能够在实际应用中发挥良好的性能。

综上所述， 本文所提算法在多信源、 小快拍条件下展现出了优良的性能， 不仅提高了ANM方法在求解DOA时的估计精度， 还克服了ESPRIT算法在处理多信源问题时的局限性。

4.3 不同快拍下各算法的均方根误差对比实验

设信噪比（Signal-to-Noise Ratio，SNR）为0 dB， 快拍数以10为间隔从10增加到110， 在每个快拍下进行100次Monte Carlo实验， 其余实验条件与4.2节相同。实验结果如图 6 所示。

从图 6 的均方根误差随快拍数变化曲线中可以看出， 随着快拍数的增加， 所有算法达到了较低的均方误差， 在多信源的条件下各算法在20快拍之前均出现了较大递减趋势， ESPRIT算法的变化趋势较为明显。同时也进一步说明了该算法在快拍数较大的情况下才具有较高的估计精度， 但在实际环境中并未能同步拾取较大的快拍数， 因而该算法具有一定的局限性。相较于实验中的大部分对比算法， 利用ANM在求解DOA时展现出了一定的性能优势。然而， 由于参数选择的局限性， 它在与本文所提的深度学习算法竞争最优参数选择以进行DOA估计时， 仍呈现出一定的劣势。相比之下， 本文所提算法在均方误差方面始终保持着最低水平。

从图6可以清晰观察到， 无论是在小快拍还是多快拍的条件下， 本文算法相较于OMP、 LASSO等压缩感知算法均具备显著优势。同时， 对于MUSIC、 ESPRIT、 CBF等传统算法在小快拍条件下估计精度不足的问题， 本文算法也实现了进一步的提升。此外， 研究结果还表明了最优参数的选取与样本数据集大小之间的关联： 样本数据集越大， 所构建的网络层参数对数据的学习能力越强， 从而能够更准确地选择最优参数。这些通过训练得到的参数对于不同的快拍数均表现出良好的普适性， 在给定条件下， 所选最优参数始终能够保持较高的DOA估计精度。

综上所述， 本文所提算法在多信源环境下， 随着快拍数的增加， 相较于其他对比算法， 均展现出了更低的均方根误差。这充分证明了通过模型训练获得的最优参数具有良好的泛化能力， 能够显著提升DOA估计精度。

4.4 不同信噪比下各算法成功率对比实验

设置与4.2节相同的实验前提来验证各算法在不同信噪比下的成功率， 实验结果如图 7 所示。实验选取4个信源角度进行成功率的估计， 在每个信噪比下进行100次Monte Carlo实验， 每次Monte Carlo实验成功的条件为{\displaystyle \sum _{i=\mathrm{1}}^K}[{\stackrel{\wedge }{\theta }}_k(j)-{\theta }_k]\le \mathrm{0.25}K，其中， K为信源总数， {\stackrel{\wedge }{\theta }}_k(j)为第j次Monte Carlo实验中第k个信源的估计值， {\theta }_k为第k个信源的真实值。

从图 7 可以看出， 随着信噪比的增加， 各算法的成功率也随之增加。其中变化较明显的是ESPRIT算法， 由于实验是在10快拍的环境下进行， 并且由图4、 5可以看出只有在信噪比与快拍数均足够大的前提下该算法才具有较高的估计精度。在信噪比达到5 dB后， 所有算法均达到较高的估计精度， 而ESPRIT算法在15 dB后才达到较高的估计精度， 但是本文所提算法始终具有较高的估计精度， 比较明显的是在0 dB以前该算法就已达到较高的成功率， 这说明在训练数据集中加入的噪声使得该算法相较于其他算法具有更强的抗噪能力， 在低信噪比的条件下仍保持较高的成功率， 同时也说明训练所选的最优参数对接收到的不同信号数据集具有泛化性， 受噪声影响较小， 与其余算法相比性能更优， 可适用于复杂环境。

综上所述， 与其余对比算法相比， 本文所提算法在DOA估计性能方面展现出了显著的优势， 尤其在与ESPRIT算法的比较中更为突出。此外， 本文算法对ANM方法进行了进一步的优化， 使得参数的选取更加合理， 从而有效提升了矢量水听器阵列信号的DOA估计精度。

4.5 相干信源误差对比实验

由于一些传统算法在实际环境中易受到相干信源的影响， 为了进一步验证本文算法的实验性能， 实验设置水听器均匀线阵阵元数为8， 快拍数为5， 信噪比为-5 dB，入射信源角度为-\mathrm{48}°、 -\mathrm{13}°、 \mathrm{23}°、 \mathrm{50}°， 频率为1 000 Hz， 实验结果如图 8 所示。

从图 8 中可以看出在多个相干源存在的条件下， OMP算法在-\mathrm{48}°和\mathrm{50}°附近估计出现了较大偏差， 在其余两个相干信源的条件下基本能准确估计角度。MUSIC算法除在-\mathrm{13}°处准确估计相干信源外， 在其余3个相干信源的环境下出现了较大的估计误差， 这也是传统子空间算法在DOA估计所存在的不足之处。ESPRIT算法只可以准确估计\mathrm{23}°处的相干信源， 在相干信源-{\mathrm{48}}^{\circ }处出现了较大的估计偏差， 而其余相干信源也已无法准确估计。LASSO算法对最后一个相干信源角度\mathrm{50}°的估计出现了一定的偏差， 其余角度的估计较为准确， 而ANM方法在-\mathrm{48}°处的估计出现了较大的误差， 其余角度均能较准估计相干信源的值， 可以看出ANM方法相较于其他传统算法可以较好的消除相干信源的影响， 但依然受到相干信源的一定制约。本文算法以及传统CBF算法在所有相干信源角度上都实现了较为精确的估计。这一优势主要源于本文算法的独特设计： 它并未直接利用信号的接收数据进行DOA估计， 而是利用矢量水听器的声压振速相关性来重构Toeplitz协方差矩阵。在这个过程中， 本文算法充分利用矢量水听器所接收的声压和振速信息， 有效去除了接收数据中的冗余部分， 最后通过直接对Toeplitz矩阵进行Rootmusic算法分解， 进而得到了DOA的估计值， 这种处理方式使得算法在相干信源的影响下表现出较小的误差。

综上所述， 本文算法在相干信源环境下展现出了优良的DOA估计性能， 其精度远超其他算法。此外， 在大量样本数据训练的基础上， 该算法对相干信源还具有一定的泛化能力。这些特点进一步凸显了本文所提算法在解决相干信源DOA估计问题上的优良性能。

4.6 各算法分辨率对比实验

实验设置矢量水听器均匀线阵阵元数为10， 各算法快拍数为10， 信噪比为0 dB。假设入射信号为\theta =\{-\mathrm{30}°，(-\mathrm{30}+\mathrm{\Delta }s)°\}， 角度差\mathrm{\Delta }s由\mathrm{2}°以步长2°变化到\mathrm{20}°， 在每个角度下进行100次Monte Carlo实验， 实验结果如图 9 所示。

从图 9 可以看出， 随着信源间角度差的增大， 各算法的均方根误差逐渐减小， 即估计精度在逐步提升。其中MUSIC、 ESPRIT、 CBF等传统算法在角度间隔差达10°后才展现出较高的估计性能， 这反映了传统算法在分辨率方面的局限性。相比之下， OMP、 LASSO、 ANM等重构算法在较小的角度间隔下就表现出明显的优势， 其分辨率有所提升， 并在角度间隔8°下达到了较高的水平。本文所提算法在更小的角度间隔差4°下就实现了较低的均方误差， 即该算法相对于其它算法具有较高的分辨率。这得益于该算法通过构建网络模型来训练数据集， 从而获取一组优化的参数解。在随机选取训练数据集角度的情况下， 所得最优参数具有一定的泛化能力， 进一步提升了算法的分辨率和DOA估计精度。

综上所述， 本文所提算法在分辨力方面相较于其他对比算法有了显著提高， 这充分验证了模型训练的合理性和有效性。

5 结　语

本文针对传统及现有稀疏重构算法在小快拍、 多信源、 相干信源下阵列信号DOA估计的不足做了一些针对性改善， 而且相较于其他对比算法具有较高的分辨率。本文所提出的基于深度学习的交替迭代乘子算法避免了固定参数选择的困境， 且通过仿真实验验证了所提算法的DOA估计性能， 使得该算法的DOA估计精度进一步提高， 该算法随着信噪比、 快拍数的增加始终具有较稳定的DOA估计性能， 而且该过程不需要进行网格划分。在原子范数理论的前提下， 借助深度交替乘子算法可以有效求解SDP问题， 为今后无网格DOA估计提供了一定的研究思路。

参考文献

[1] ZHONG S， XU L， YANG H， et al. DOA estimation of underwater object based on cross-spectrum deconvolution algorithm［C］//Proceedings of Chinese Intelligent Automation Conference. Singapore： Springer， 2022： 627-636.

[2] PENG S， CHEN B， YANG M. Joint sparse recovery for direction of arrival based on the generalized MUSIC criterion［J］. Digital Signal Processing， 2022， 122： 103382.

[3] BENSALEM M， BARKAT O. DOA estimation of linear dipole array with known mutual coupling based on ESPRIT and MUSIC［J］. Radio Science， 2022， 57（2）： 1-15.

[4] WANG P， HUANG J， HE W， et al. Maximum likelihood DOA estimation based on improved invasive weed optimization algorithm and application of MEMS vector hydrophone array［J］. AIMS Mathematics， 2022， 7（7）： 12342-12363.

[5] 何伟嘉， 黄嘉俊， 张靖奇， 等. 基于金枪鱼算法的最大似然信号方位与频率联合估计［J］. 测试技术学报， 2023， 37（1）： 60-65.

[6] HE Weijia， HUANG Jiajun， ZHANG Jingqi， et al. Joint angle-frequency estimation with maximum likelihood based on tuna swarm optimization algorithm［J］. Journal of Test and Measurement Technology， 2023， 37（1）： 60-65. （in Chinese）

[7] DONOHO D L. Compressed sensing［J］. IEEE Transactions on Information Theory， 2006， 52（4）： 1289-1306.

[8] CANDES E J， ROMBERG J， TAO T. Robust uncertainty principles： exact signal reconstruction from highly incomplete frequency information［J］. IEEE Transactions on Information Theory， 2006， 52（2）： 489-509.

[9] ZHU S， CHEN S， PENG X， et al. A signal reconstruction method of wireless sensor network based on compressed sensing［J］. EURASIP Journal on Wireless Communications and Networking， 2020（1）： 106.

[10] 刁弘扬， 胡洲勇， 禹永植. 一种改进广义正交匹配追踪的DOA估计方法［J］. 应用科技， 2020， 47（4）： 54-58.

[11] DIAO Hongyang， HU Zhouyong， YU Yongzhi. An improved generalized orthogonal matching pursuit method used in the direction of arrival estimation［J］. Applied Science and Technology， 2020， 47（4）： 54-58. （in Chinese）

[12] KUKREJA S L， LÖFBERG J， BRENNER M J. A least absolute shrinkage and selection operator （LASSO） for nonlinear system identification［J］. IFAC Proceedings Volumes， 2006， 39（1）： 814-819.

[13] 王勇， 李韬， 项建弘. 一种基于加权改进平滑l0范数的DOA估计方法［J］.应用科技， 2022， 49（4）： 38-43.

[14] WANG Yong， LI Tao， XIANG Jianhong. A DOA estimation method based on weighted improved smoothing l0 norm［J］. Applied technology， 2022， 49（4）： 38-43. （in Chinese）

[15] 刘永豪， 徐明. 交替方向乘子法的矢量水听器DOA估计方案［J］. 舰船科学技术， 2024， 46（12）： 140-143.

[16] LIU Yonghao， XU Ming. DOA estimation scheme of vector hydrophone by alternating direction multiplier method［J］. Ship Science and Technology， 2024， 46（12）： 140-143. （in Chinese）

[17] 郑文康， 魏志晴， 白艳萍， 等. 基于可分离替代函数算法的DOA估计方法［J］. 陕西科技大学学报， 2024， 42（1）： 197-205.

[18] ZHENG Wenkang， WEI Zhiqing， BAI Yanping， et al. DOA estimation based on separable surrogate function algorithm［J］. Journal of Shaanxi University of Science & Technology， 2024， 42（1）： 197-205. （in Chinese）

[19] 孙国仓， 惠俊英， 郭龙祥. 声压振速联合频率估计［J］. 计算机工程与应用， 2008（10）： 154-156.

[20] SUN Goucang， HUI Junying， GUO Longxiang. Combined frequency estimation of sound pressure and velocity［J］. Computer Engineering and Applications， 2008（10）： 154-156. （in Chinese）

[21] 姚直象， 姜可宇， 郭瑞， 等. 基于声压振速联合处理的矢量阵旋转不变子空间方位估计方法［J］. 北京理工大学学报， 2012， 32（5）： 513-516.

[22] YAO Zhixiang， JIANG Keyu， GUO Rui， et al. An improved bearing estimation algorithm using acoustic vector sensor array based on rotational invariance technique［J］. Transactions of Beijing Institute of Technology， 2012， 32（5）： 513-516. （in Chinese）

[23] 张旭， 时胜国， 朱晓春， 等. 使用协方差矩阵分解的矢量阵声压振速联合处理波达方向估计［J］. 声学学报， 2023， 48（5）： 937-949.

[24] ZHANG Xu， SHI Shengguo， ZHU Xiaochun， et al. Direction of arrival estimation by joint processing of vector array sound pressure and vibration velocity using covariance matrix decomposition［J］. Acta Acustica， 2023， 48（5）： 937-949. （in Chinese）

[25] 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.

[26] YANG Z， LI J， STOICA P， et al. Sparse methods for direction-of-arrival estimation［M］. Amsterdam： Elsevier， 2018： 509-581.

[27] BHASKAR B N， RECHT B. Atomic norm denoising with applications to line spectral estimation［C］//2011 49th Annual Allerton Conference on Communication， Control， and Computing （Allerton）， 2011： 261-268.

[28] XIAO P， LIAO B， DELIGIANNIS N. DeepFPC： a deep unfolded network for sparse signal recovery from 1-Bit measurements with application to DOA estimation［J］. Signal Processing， 2020， 176： 107699.

[29] LI Y， CHI Y. Off-the-grid line spectrum denoising and estimation with multiple measurement vectors［J］. IEEE Transactions on Signal Processing， 2016， 64（5）： 1257-1269.

[30] TOH K C， TODD M J， TÜTÜNCÜ R H. On the implementation and usage of SDPT3–A Matlab software package for semidefinite-quadratic-linear programming， version 4.0［M］. New York： Springer， 2012： 715-754.

[31] BOYD S. Distributed optimization and statistical learning via the alternating direction method of multipliers ［J］. Foundations and Trends in Machine Learning， 2010， 3（1）： 1-122.

[32] LIU S， MAO Z， ZHANG Y D， et al. Rank minimization-based toeplitz reconstruction for DoA estimation using coprime array ［J］. IEEE Communications Letters， 2021， 25（7）： 2265-2269.

[33] BARABELL A. Improving the resolution performance of eigenstructure-based direction-finding algorithms［C］//ICASSP IEEE International Conference on Acoustics， Speech， and Signal Processing， 1983： 336-339.

[34] RAO B D， HARI K V S. Performance analysis of root-music［J］. IEEE Transactions on Acoustics， Speech， and Signal Processing， 1989， 37（12）： 1939-1949.

[35] ZHU H， FENG C， FENG W， et al. A deep learning approach for sparse single snapshot DOA estimation［J］. Journal of Signal Processing， 2022， 38： 2114-2123.

[36] RUMELHART D E， HINTON G E， WILLIAMS R J. Learning representations by back-propagating errors［J］. Nature， 1986， 323： 533-536.

基金资助

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

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

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

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

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

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

山西省回国留学人员科研资助项目(2021-108)

山西省回国留学人员科研资助项目(2022-149)