基于迭代收缩阈值算法的DOA估计方法

DOA Estimation Method Based on Iterative Shrinkage Threshold Algorithm

摘要

针对传统波达方向（Direction of Arrival， DOA）估计算法在低信噪比、 小快拍的条件下估计精度不高的问题， 提出了一种基于迭代收缩阈值算法的矢量水听器阵列多快拍DOA估计方法。首先， 对空域进行等角度划分， 构造超完备冗余字典， 建立基于信号多快拍条件下的DOA估计模型， 然后， 采用迭代收缩阈值算法解决稀疏重构问题， 求解出信号的稀疏系数矩阵， 最后， 将稀疏矩阵中行向量的范数映射到划分好的网格上， 得到DOA估计值。仿真实验结果表明： 该方法在低信噪比、 小快拍条件下比OMP、 MUSIC和CBF等传统算法拥有更高的DOA估计精度和更强的鲁棒性。

Abstract

Aiming at the problem that the traditional direction of arrival （DOA） estimation algorithm does not have high estimation accuracy under the conditions of low signal-to-noise ratio and small snapshots， a multi-snapshot DOA estimation method based on the iterative shrinkage threshold algorithm for vector hydrophone arrays is proposed. Firstly， the airspace domain is divided into equal angles， and an ultra-complete redundant dictionary is constructed to establish a DOA estimation model based on the multi-fast-beat condition of the signal.Then，the iterative shrinkage threshold algorithm is used to solve the sparse reconstruction problem， and the sparse coefficient matrix of the signal is solved. Finally， the paradigms of the row vectors in the sparse matrix are mapped onto the well-demarcated mesh， and the DOA estimation value is obtained. Simulation experimental results show that the method has higher DOA estimation accuracy and stronger robustness than traditional algorithms such as OMP， MUSIC and CBF algorithms under low signal-to-noise ratio and small snap conditions.

Graphical abstract

关键词

波达方向（DOA）估计 / 压缩感知 / 稀疏重构 / 迭代收缩阈值算法 / 矢量水听器

Key words

direction of arrival（DOA） estimation / compressive sensing / sparse reconstruction / iterative shrinkage threshold algorithm / vector hydrophones

0 引　言

波达方向（Direction of Arrival， DOA）估计在雷达定位、 水声目标识别以及军事等领域发挥着重要作用， 并促进了这些领域的快速发展^［1‐3］。早期的DOA估计以常规波束形成（Conventional Beam Forming， CBF）^［4］为代表， 该算法的分辨率受到瑞利限约束， 导致无法准确分辨位于一个波束宽度内的声源目标。之后针对DOA估计分辨率低的问题， 提出了两大类算法： 经典子空间分解类算法^［5］包括多重信号分类（Multiple Signal Classification， MUSIC）算法^［6］和旋转不变子空间法（Estimation of Signal Parameters using Rotational Invariance Techniques，ESPRIT）^［7］。其中， MUSIC算法是在已知信源个数前提下， 首先通过特征分解划分出信号和噪声两个子空间， 然后利用谱峰搜寻的方式寻找目标声源， 计算量相对较大且只能用于不相关信源。而ESPRIT算法利用信号的相位信息进行频率估计， 然后利用阵列接收数据的旋转不变特性来求解DOA， 相较MUSIC 算法， 计算量减小并且不需要谱峰搜索。但该类算法主要依据接收到的信号的统计特性来进行DOA估计， 这往往需要大量快拍数， 也就意味着需要进行较长时间运算， 并且在小快拍下其算法的性能会大大下降。另一类是子空间拟合类算法^［8］， 其中最典型的是最大似然（Maximum Likelihood， ML）算法^［9］， 该算法需要搜索多维谱峰， 具有很好的统计特性， 但计算很复杂， 而且还需设置合适的初始值。因此， 这些传统的DOA估计方法并不适用于实际应用。

压缩感知（Compressive Sensing， CS）理论自2006年以来在数学界得到了迅速发展。Donoho等^［10］将压缩感知中稀疏重构思想引入到信号的DOA估计中， 这类方法比传统空谱估计方法具有更好的性能^［11］。

近年来众多DOA估计方法中， 稀疏信号重构（Sparse Signal Reconstruction， SSR）技术^［12］是一种重要的方法， 其基本原理是利用信号的稀疏特性求解欠定方程。在数学中， l_{\mathrm{0}}范数可以用来表征信号的稀疏性， 但是计算l_{\mathrm{0}}范数是一个NP-hard问题。对此， 学者们提出了许多方法来解决稀疏信号的重构问题， 主要包括贪婪类算法、 凸优化方法以及稀疏贝叶斯方法^［13］。其中， 采用贪婪算法求解欠定方程最为常用。而贪婪类算法是稀疏重构算法中的一类， 目前已经提出了一些经典算法， 如匹配追踪（Matching Pursuit， MP）算法^［14］、 正交匹配追踪（Orthogonal-Matching-Pursuit， OMP）算法^［15］等， 这类算法的实现较为简单并且计算复杂度较低， 信号在信噪比较高的情况下重构效率高， 但是在低信噪比下重构精度并不佳。而且， 已有许多凸优化算法被成功应用于信号重构问题， 如基追踪降噪（Basis Pursuit Denoising， BPDN）算法^［16］等， 然而， 这些算法在DOA估计性能方面会受到方向网格的限制。若网格划分过粗， 则DOA估计精度会大大下降； 而网格划分过细， 计算量必然增大， 从而导致耗时间太长， 因此都不适合实际应用。

迭代收缩算法是一种新兴的用于解决l_{\mathrm{1}}优化问题的方法， 通过在迭代过程中对解进行收缩操作来提高求解效率。迭代收缩阈值算法（Iterative Shrinkage Thresholding Algorithm， ISTA）^［17］在每次迭代过程中， 用收缩阈值的方式来更新待恢复的信号。为了提高矢量水听器阵列信号在低信噪比和小快拍条件下的DOA估计精度， 本文提出了一种基于迭代收缩阈值算法的DOA估计方法， 并将其成功应用于矢量水听器阵列信号DOA估计中， 该算法不仅计算简单而且收敛性得到了一定的保证。在信号具有空域稀疏性的前提下， 该算法利用信号所对应的稀疏解相较于其他非稀疏解的局部最大性来重构原始信号， 可获得高精度的DOA估计值。在实际的水声实验中， 应用本文算法可以有效进行DOA估计， 在水声信号处理领域具有一定的研究价值。

1 矢量水听器阵列信号输出模型

现在考虑矢量水听器阵列的接收信号模型， 假设t时刻有K个信源s_{\mathrm{1}}(t),s_{\mathrm{2}}(t), \cdots ,s_K(t)入射到由M个矢量水听器按间距d组成的均匀线阵中， 在t时刻第k个信号的水平方位角为{\theta }_k, k=\mathrm{1,2},\cdots ,K。在水声系统中， 远场信号可近似当作平面波处理， 得到由矢量水听器阵列组成的均匀线阵模型， 示意图如图 1 所示。

图 1 中， v_x,v_y为声压部分沿x轴和y轴方向的振速分量。在t时刻， 矢量水听器的声压和2个振速分量有如下关系

\mathit{v}(t)=[\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k\cdot \mathit{\xi }+\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k\cdot \mathit{\eta }]p(t),（1）

式中： p(t)为声压标量； \mathit{\xi }， \mathit{\eta }为正交的单位向量， 则p(t)与v_x,v_y的关系可以表示为

\left\{\begin{array}{l}p(t)=z(t),\\ v_x(t)=p(t)\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k,\\ v_y(t)=p(t)\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k\mathrm{。}\end{array}\right.（2）

可见， 矢量水听器可输出声压与振速三路信号， 即

\begin{array}{c}\mathit{x}(t)={\left[p(t),v_x(t),v_y(t)\right]}^{\mathrm{T}}=\\ z(t){\left[\mathrm{1},\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k,\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k\right]}^{\mathrm{T}}\mathrm{。}\end{array}（3）

将矢量水听器阵列输出改写成t时刻单快拍的形式， 即

\mathit{X}(t)=\mathit{A}(\theta )\mathit{S}(t)+\mathit{N}(t)，（4）

式中： \mathit{X}(t)={\left[x_{\mathrm{1}}(t),x_{\mathrm{2}}(t),\cdots ,x_M(t)\right]}^{\mathrm{T}}为在t时刻阵列的输出数据，\mathit{X}(t)\in {\mathbb{C}}^{\mathrm{3}M\times \mathrm{1}}； \mathit{A}(\theta )=[\mathit{a}({\theta }_{\mathrm{1}})\otimes \mathit{u}({\theta }_{\mathrm{1}}),\mathit{a}({\theta }_{\mathrm{2}})\otimes \mathit{u}({\theta }_{\mathrm{2}}),\cdots ,\mathit{a}({\theta }_K)\otimes \mathit{u}({\theta }_K)]为\mathrm{3}M\times K维阵列流形矩阵， \otimes为克罗内克积， \mathit{a}({\theta }_k)={\left[\mathrm{1},{\mathrm{e}}^{-j\mathrm{2}\pi d\frac{\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k}{\tau }},\cdots ,{\mathrm{e}}^{-j\mathrm{2}\pi d(M-\mathrm{1})\frac{\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k}{\tau }}\right]}^{\mathrm{T}}为第k个信源的声压响应系数， \tau为入射信源的波长， \mathit{u}({\theta }_k)={\left[\mathrm{1},\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k,\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k\right]}^{\mathrm{T}}， k=\mathrm{1,2},\cdots ,K； \mathit{S}(t)=\left[s\mathrm{1}(t)\right.， {\left.s\mathrm{2}(t),\cdots ,{\mathit{s}}_K(t)\right]}^{\mathrm{T}}\in {\mathbb{C}}^{K\times \mathrm{1}}为t时刻信号源矢量； \mathit{N}(t)={\left[{\mathit{n}}_{\mathrm{1}}(t),{\mathit{n}}_{\mathrm{2}}(t),\cdots ,{\mathit{n}}_M(t)\right]}^{\mathrm{T}}为t时刻阵列接收的噪声矢量， \mathit{N}(t)\in {\mathbb{C}}^{\mathrm{3}M\times \mathrm{1}}。

2 基于迭代阈值收缩算法的DOA估计方法

2.1 稀疏信号DOA估计模型

下面考虑多快拍的情况。由阵列的输出模型可知， 阵列接收信号的数据矢量可表示为

\overline{\mathit{X}(t)}=\mathit{A}(\theta )\overline{\mathit{S}(t)}+\overline{\mathit{N}(t)}，（5）

式中： \overline{\mathit{X}(t)}=\left[{\mathit{X}}_{\mathrm{1}}(t),{\mathit{X}}_{\mathrm{2}}(t),\cdots ,{\mathit{X}}_L(t)\right]\in {\mathbb{C}}^{\mathrm{3}M\times L}； \overline{\mathit{S}(t)}={\left[{\mathit{S}}_{\mathrm{1}}(t),{\mathit{S}}_{\mathrm{2}}(t),\cdots ,{\mathit{S}}_L(t)\right]}^{\mathrm{T}}为K\times L维稀疏信号； \overline{\mathit{N}(t)}={\left[{\mathit{N}}_{\mathrm{1}}(t),{\mathit{N}}_{\mathrm{2}}(t),\cdots ,{\mathit{N}}_L(t)\right]}^{\mathrm{T}}∈{\mathbb{C}}^{\mathrm{3}M\times L}。

实际入射信号在空域上是稀疏的， 由于空间中有K个入射信号， 将空间进行等间隔划分为N份， 且声源个数K\ll N， 此时就构造出一个N\times L维的稀疏信号\stackrel{̑}{\mathit{S}（t)}， \stackrel{̑}{\mathit{S}（t)}中只有实际存在目标声源的K个位置的元素非零， 其它N-K个位置的元素为零。同时， 将阵列流形矩阵\mathit{A}(\theta )∈{\mathbb{C}}^{\mathrm{3}M\times K}扩充为划分好的全部方向角的超完备冗余字典\mathit{G}(\theta )\in {\mathbb{C}}^{\mathrm{3}M\times N}， 即

\mathit{G}(\theta )=\left[\mathit{a}({\theta }_{\mathrm{1}})\otimes \mathit{u}({\theta }_{\mathrm{1}}),\mathit{a}({\theta }_{\mathrm{2}})\otimes \mathit{u}({\theta }_{\mathrm{2}}),\cdots ,\mathit{a}({\theta }_N)\otimes \mathit{u}({\theta }_N)\right]\mathrm{。}（6）

故基于压缩感知理论的多快拍DOA估计模型可以转化为

\mathit{Y}(t)=\mathit{G}(\theta )\stackrel{̑}{\mathit{S}（t)}+\stackrel{̑}{\mathit{N}(t)}。（7）

实际DOA估计就是对稀疏信号的重建过程， 其本质是通过阵列接收数据\mathit{Y}(t)和字典\mathit{G}(\theta )来重构稀疏信号‖\stackrel{̑}{\mathit{S}（t)}‖， \mathit{G}(\theta )的列之间需要满足有限等距性质（ Restrited Isomely Property， RIP）^［17］。显然， 式（7）是欠定方程， 因此， 需要将其转化为求解l_{\mathrm{0}}范数的优化问题， 即

\stackrel{̑}{\mathit{S}（t)}=\mathrm{a} \mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}{‖\stackrel{̑}{\mathit{S}（t)}‖}_{\mathrm{0}}, s.t. {‖\mathit{Y}(t)-\mathit{G}(\theta )\stackrel{̑}{\mathit{S}（t)}‖}_{\mathrm{2}}^{\mathrm{2}}\le {\epsilon }^{\mathrm{2}}\mathrm{。}（8）

式（8）中， 寻找非零元素最少的解也就是最稀疏的解。{‖\stackrel{̑}{\mathit{S}（t)}‖}_{\mathrm{0}}为稀疏信号中非零元素的个数， {\epsilon }^{\mathrm{2}}为受噪声能量影响的容噪参数。由于式（8）是一个NP-hard问题^［18］， 无法直接求解。为确保信号的稀疏性和信号重建的精确性， 学者们提出将式（11）转化为l_{\mathrm{1}}范数的优化问题来求解， 即

\stackrel{̑}{\mathit{S}（t)}=\mathrm{a} \mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}{‖\stackrel{̑}{\mathit{S}（t)}‖}_{\mathrm{1}}, s.t. {‖\mathit{Y}(t)-\mathit{G}(\theta )\stackrel{̑}{\mathit{S}（t)}‖}_{\mathrm{2}}^{\mathrm{2}}\le {\epsilon }^{\mathrm{2}}\mathrm{。}（9）

一般地， 将式（9）转换成求解无约束优化问题， 即

\stackrel{̑}{\mathit{S}（t)}=\underset{\stackrel{̑}{\mathit{S}(t)}}{\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}} {‖\mathit{Y}(t)-\mathit{G}(\theta )\stackrel{̑}{\mathit{S}（t)}‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda {‖\stackrel{̑}{\mathit{S}（t)}‖}_{\mathrm{1}}\mathrm{。}（10）

式（10）是一个l_{\mathrm{1}}范数正则化模型， 在DOA估计中得到广泛运用。其中， \lambda为正则化参数， 在凸优化类方法中， 正则化参数的选择取决于背景噪声， 而正则化参数大小的选择会影响到信号重构精度， 因此， 如何选择合适的正则化参数极其困难^［19］。

最小化l_{\mathrm{1}}范数问题是一个线性规划问题， 需要转换为二阶锥规划问题^［20-21］， 然后使用内点法求解上述问题。但大多数情况下， 由于数据维数过大， 算法的复杂性很高， 所以求解起来非常耗时。因此， 对于大规模信号快速DOA估计来说， 这种方法并不适用。目前比较常用的方法是利用 Matlab中CVX工具箱来解决最小化l_{\mathrm{1}}范数优化问题。但是由于凸优化工具箱CVX求解速度较慢， 不适合实际应用， 许多学者利用梯度法求解最小化l_{\mathrm{1}}范数问题。梯度法是一种新的计算方法， 它的计算量主要在矩阵与向量的乘积上， 具有较低的运算量， 并且构造简单， 易于实现。针对上述问题， 本文在此基础上提出了一种新的基于梯度的方法， 即迭代收缩阈值算法， 然后采用该算法来解决压缩感知模型中的信号重构问题。

2.2 基于迭代收缩阈值算法的稀疏重构算法

迭代收缩阈值算法（Iterative Shrinkage Thresholding Algorithm， ISTA）是求解线性逆问题的经典方法， 是经典梯度法的扩展， 也是一种迭代优化算法， 用最小化替代函数来确保原目标函数的迭代下降。同时， ISTA算法也被认为是一种有效的稀疏重构算法， 该算法相较于传统算法求解过程简单便捷， 重构精度更好。

本文首先考虑如下问题

\underset{\stackrel{̑}{S}}{\mathrm{m}\mathrm{i}\mathrm{n}}\mathit{F}\left(\stackrel{̑}{\mathit{S}}\right)=\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{Y}-\mathit{G}\stackrel{̑}{\mathit{S}}‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda {‖\stackrel{̑}{\mathit{S}}‖}_{\mathrm{1}}，（11）

式中： \stackrel{̑}{\mathit{S}}为待恢复信号； \mathit{G}为给定的矩阵； \mathit{Y}为测量信号； \lambda >\mathrm{0}为正则化参数。由于目标函数\mathit{F}\left(\stackrel{̑}{\mathit{S}}\right)并不容易优化， 根据Majorization-Minimization优化框架， 可以优化替代目标函数

\mathit{G}\left(\stackrel{̑}{\mathit{S}},\mathit{Z}\right)=\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{Y}-\mathit{G}\stackrel{̑}{\mathit{S}}‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda {‖\stackrel{̑}{\mathit{S}}‖}_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{G}\stackrel{̑}{\mathit{S}}-\mathit{G}\mathit{Z}‖}_{\mathrm{2}}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}{‖\stackrel{̑}{\mathit{S}}-\mathit{Z}‖}_{\mathrm{2}}^{\mathrm{2}}\mathrm{。}（12）

式（12）中， 需满足如下关系

\mathit{G}\left(\stackrel{̑}{\mathit{S}},\mathit{Z}\right)\ge \mathit{F}\left(\stackrel{̑}{\mathit{S}}\right)，\mathit{G}\left(\mathit{Z},\mathit{Z}\right)=\mathit{F}\left(\stackrel{̑}{\mathit{S}}\right)\mathrm{。}（13）

下面对替代目标函数\mathit{G}\left(\stackrel{̑}{\mathit{S}},\mathit{Z}\right)进行化简， 如式（14）所示。

\begin{array}{c}\mathit{G}\left(\stackrel{̑}{\mathit{S}},\mathit{Z}\right)=\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{Y}-\mathit{G}\stackrel{̑}{\mathit{S}}‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda {‖\stackrel{̑}{\mathit{S}}‖}_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{G}\stackrel{̑}{\mathit{S}}-\mathit{G}\mathit{Z}‖}_{\mathrm{2}}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}{‖\stackrel{̑}{\mathit{S}}-\mathit{Z}‖}_{\mathrm{2}}^{\mathrm{2}}=\\ \frac{\mathrm{1}}{\mathrm{2}}\left({‖\mathit{Y}‖}_{\mathrm{2}}^{\mathrm{2}}+\mathrm{2}\mathrm{t}\mathrm{r}({\stackrel{̑}{\mathit{S}}}^{\mathrm{{\rm T}}}{\mathit{G}}^{\mathrm{{\rm T}}}\mathit{Y})+{‖\mathit{G}\stackrel{̑}{\mathit{S}}‖}_{\mathrm{2}}^{\mathrm{2}}\right)-\frac{\mathrm{1}}{\mathrm{2}}\left({‖\mathit{G}\stackrel{̑}{\mathit{S}}‖}_{\mathrm{2}}^{\mathrm{2}}-\mathrm{2}\mathrm{t}\mathrm{r}({\stackrel{̑}{\mathit{S}}}^{\mathrm{{\rm T}}}{\mathit{G}}^{\mathrm{{\rm T}}}\mathit{G}\mathit{Z})+{‖\mathit{G}\mathit{Z}‖}_{\mathrm{2}}^{\mathrm{2}}\right)+\\ \frac{\mathrm{1}}{\mathrm{2}}\left({‖\stackrel{̑}{\mathit{S}}‖}_{\mathrm{2}}^{\mathrm{2}}-\mathrm{2}\mathrm{t}\mathrm{r}({\stackrel{̑}{\mathit{S}}}^{\mathrm{{\rm T}}}\mathit{Z})+{‖\mathit{Z}‖}_{\mathrm{2}}^{\mathrm{2}}\right)+\lambda {‖\stackrel{̑}{\mathit{S}}‖}_{\mathrm{1}}=\\ \frac{\mathrm{1}}{\mathrm{2}}\left({‖\stackrel{̑}{\mathit{S}}‖}_{\mathrm{2}}^{\mathrm{2}}-\mathrm{2}\mathrm{t}\mathrm{r}({\stackrel{̑}{\mathit{S}}}^{\mathrm{{\rm T}}}{\mathit{G}}^{\mathrm{{\rm T}}}\mathit{Y})+\mathrm{2}\mathrm{t}\mathrm{r}({\stackrel{̑}{\mathit{S}}}^{\mathrm{{\rm T}}}{\mathit{G}}^{\mathrm{{\rm T}}}\mathit{G}\mathit{Z})-\mathrm{2}\mathrm{t}\mathrm{r}({\stackrel{̑}{\mathit{S}}}^{\mathrm{{\rm T}}}\mathit{Z})\right)+\frac{\mathrm{1}}{\mathrm{2}}\left({‖\mathit{Y}‖}_{\mathrm{2}}^{\mathrm{2}}-{‖\mathit{G}\mathit{Z}‖}_{\mathrm{2}}^{\mathrm{2}}+{‖\mathit{Z}‖}_{\mathrm{2}}^{\mathrm{2}}\right)+\lambda {‖\stackrel{̑}{\mathit{S}}‖}_{\mathrm{1}}=\\ \frac{\mathrm{1}}{\mathrm{2}}\left({‖\stackrel{̑}{\mathit{S}}‖}_{\mathrm{2}}^{\mathrm{2}}-\mathrm{2}{\stackrel{̑}{\mathit{S}}}^{\mathrm{{\rm T}}}\mathrm{t}\mathrm{r}\left({\mathit{G}}^{\mathrm{{\rm T}}}\mathit{Y}-{\mathit{G}}^{\mathrm{{\rm T}}}\mathit{G}\mathit{Z}+\mathit{Z}\right)\right)+\frac{\mathrm{1}}{\mathrm{2}}\left({‖\mathit{Y}‖}_{\mathrm{2}}^{\mathrm{2}}-{‖\mathit{G}\mathit{Z}‖}_{\mathrm{2}}^{\mathrm{2}}+{‖\mathit{Z}‖}_{\mathrm{2}}^{\mathrm{2}}\right)+\lambda {‖\stackrel{̑}{\mathit{S}}‖}_{\mathrm{1}}=\\ \frac{\mathrm{1}}{\mathrm{2}}\left({‖\stackrel{̑}{\mathit{S}}‖}_{\mathrm{2}}^{\mathrm{2}}-\mathrm{2}{\stackrel{̑}{\mathit{S}}}^{\mathrm{{\rm T}}}\mathrm{t}\mathrm{r}\left({\mathit{G}}^{\mathrm{{\rm T}}}\mathit{Y}-{\mathit{G}}^{\mathrm{{\rm T}}}\mathit{G}\mathit{Z}+\mathit{Z}\right)+{‖{\mathit{G}}^{\mathrm{{\rm T}}}\mathit{Y}-{\mathit{G}}^{\mathrm{{\rm T}}}\mathit{G}\mathit{Z}+\mathit{Z}‖}_{\mathrm{2}}^{\mathrm{2}}\right)+\\ \frac{\mathrm{1}}{\mathrm{2}}\left({‖\mathit{Y}‖}_{\mathrm{2}}^{\mathrm{2}}-{‖\mathit{G}\mathit{Z}‖}_{\mathrm{2}}^{\mathrm{2}}+{‖\mathit{Z}‖}_{\mathrm{2}}^{\mathrm{2}}-{‖{\mathit{G}}^{\mathrm{{\rm T}}}\mathit{Y}-{\mathit{G}}^{\mathrm{{\rm T}}}\mathit{G}\mathit{Z}+\mathit{Z}‖}_{\mathrm{2}}^{\mathrm{2}}\right)+\lambda {‖\stackrel{̑}{\mathit{S}}‖}_{\mathrm{1}}=\\ \frac{\mathrm{1}}{\mathrm{2}}\left({‖\stackrel{̑}{\mathit{S}}-\left(\left({\mathit{G}}^{\mathrm{{\rm T}}}(\mathit{Y}-\mathit{G}\mathit{Z})+\mathit{Z}\right)\right)‖}_{\mathrm{2}}^{\mathrm{2}}\right)+\lambda {‖\stackrel{̑}{\mathit{S}}‖}_{\mathrm{1}}+k，\end{array}（14）

式中： k为一个与\stackrel{̑}{\mathit{S}}无关的项，且k={‖\mathit{Y}‖}_{\mathrm{2}}^{\mathrm{2}}-{‖\mathit{G}\mathit{Z}‖}_{\mathrm{2}}^{\mathrm{2}}+{‖\mathit{Z}‖}_{\mathrm{2}}^{\mathrm{2}}-{‖{\mathit{G}}^{\mathrm{{\rm T}}}\mathit{Y}-{\mathit{G}}^{\mathrm{{\rm T}}}\mathit{G}\mathit{Z}+\mathit{Z}‖}_{\mathrm{2}}^{\mathrm{2}}。令{\mathit{U}}_{\mathrm{i}}^{(k)}={\mathit{G}}^{\mathrm{{\rm T}}}(\mathit{Y}-\mathit{G}{\mathit{Z}}^{(k)})+{\mathit{Z}}^{(k)}。所以优化替代函数时， 可以等价于优化式（15）， 即

\sum _i{\mathit{G}}^{\mathrm{*}}\left(\stackrel{̑}{{\mathit{S}}_i},{\mathit{Z}}^{\mathrm{*}}\right)=\sum _i\frac{\mathrm{1}}{\mathrm{2}}{\left(\stackrel{̑}{{\mathit{S}}_i}-{\mathit{U}}_i^{(k)}\right)}^{\mathrm{2}}+\lambda \left|\stackrel{̑}{{\mathit{S}}_i}\right|\mathrm{。}（15）

当{\mathit{G}}^{\mathrm{*}}\left(\stackrel{̑}{{\mathit{S}}_i},{\mathit{Z}}^{\mathrm{*}}\right)的偏导为0时能取到最小值， 即

\frac{\partial {\mathit{G}}^{\mathrm{*}}\left(\stackrel{̑}{{\mathit{S}}_i},{\mathit{Z}}^{\mathrm{*}}\right)}{\partial \stackrel{̑}{{\mathit{S}}_i}}\in \left({\stackrel{̑}{\mathit{S}}}_i^{\mathrm{*}}-{\mathit{U}}_i^{(k)}\right)+\lambda \mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(\stackrel{̑}{{\mathit{S}}_i^{\mathrm{*}}})=\mathrm{0},（16）

移项得到

{\mathit{U}}_i^{(k)}\in {\stackrel{̑}{\mathit{S}}}_i^{\mathrm{*}}+\lambda \mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left({\stackrel{̑}{\mathit{S}}}_i^{\mathrm{*}}\right),（17）

则可得梯度下降公式为

\stackrel{̑}{S^{(k)}}={\mathit{G}}^{\mathrm{{\rm T}}}(Y-\mathit{G}{\mathit{Z}}^{(k)})+{\mathit{Z}}^{(k)}\mathrm{。}（18）

在复数域上， 软阈值函数为

T_{\omega }(\stackrel{̑}{\mathit{S}})=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(\stackrel{̑}{\mathit{S}})\mathrm{m}\mathrm{a}\mathrm{x}(|\stackrel{̑}{\mathit{S}}|-\omega ,\mathrm{0})\mathrm{。}（19）

则可得到式（7）的迭代解为：

\stackrel{̑}{S^{(k+\mathrm{1})}}=T_{\omega }(\stackrel{̑}{S^{(k)}})\mathrm{。}（20）

式（20）中的收缩算子T_{\omega }(\cdot )旨在近似满足正则项， 它的作用是将输入值通过收缩函数映射到所期待的输出值。本文的收缩算子T_{\omega }(\cdot )是把原点附近的值映射为零， 也就是进行收缩操作。以上就 是ISTA算法的求解过程以及核心思想。

因此， 通过迭代实现\mathit{G}(\stackrel{̑}{\mathit{S}},\mathit{Z})的最小化， 其中ISTA求解的迭代值\stackrel{̑}{S^{(k+\mathrm{1})}}仅依赖于前一次的迭代值\stackrel{̑}{S^{(k)}}， 然后解序列最终会收敛到替换函数\mathit{G}(\stackrel{̑}{\mathit{S}},\mathit{Z})的最小点， 即原函数\mathit{F}\left(\stackrel{̑}{\mathit{S}}\right)的最小点。

根据上述一系列推导， 能够找到式（10）中所提问题的最优解。也就是在基于压缩感知的多快拍DOA估计模型中， 可以使用ISTA算法来解决上述问题， 该算法的原理图如图 2 所示。

本文提出的迭代收缩阈值算法的DOA估计方法的具体实现步骤描述如下：

输入： 稀疏信号维度N， 稀疏度K， 快拍数L， 字典矩阵\mathit{G}(\theta )\in {\mathbb{C}}^{\mathrm{3}M\times N}， 矢量水听器阵列的接收信号\mathit{Y}(t)\in {\mathbb{C}}^{\mathrm{3}M\times L}。

第1步： 设置初始解\stackrel{̑}{{\mathit{S}}_{\mathrm{0}}(t)}=0， 初始化残差R_{\mathrm{0}}=\mathit{Y}(t)-\mathit{G}(\theta ){\stackrel{̑}{\mathit{S}(t)}}^{(k)}=\mathit{Y}(t)， 步长为1。

迭代次数： k=\mathrm{1}。

第2步： 计算第k次梯度迭代的结果\stackrel{̑}{\mathit{S}{(t)}^{(\mathrm{k})}}， 其中\stackrel{̑}{\mathit{S}{(t)}^{(\mathrm{k})}}=\mathit{G}{(\theta )}^{\mathrm{{\rm T}}}(\mathit{Y}(t)-\mathit{G}(\theta )Z{(t)}^{(k)})+Z{(t)}^{(k)}。

第3步： 用公式\stackrel{̑}{\mathit{S}{(t)}^{(\mathrm{k}+\mathrm{1})}}=T_{\omega }\left(\stackrel{̑}{\mathit{S}{(t)}^{(\mathrm{k})}}\right)对上一步所求的迭代结果进行收缩运算。

第4步： 更新残差{\mathit{R}}_k=\mathit{Y}(t)-\mathit{G}(\theta )\stackrel{̑}{\mathit{S}{(t)}^{(\mathrm{k}+\mathrm{1})}}。

第5步： 判断收敛性， 如果\frac{{‖\stackrel{̑}{\mathit{S}{(t)}^{(\mathrm{k}+\mathrm{1})}}-\stackrel{̑}{\mathit{S}{(t)}^{(\mathrm{k})}}‖}_{\mathrm{2}}^{\mathrm{2}}}{{‖\stackrel{̑}{\mathit{S}{(t)}^{(\mathrm{k}+\mathrm{1})}}‖}_{\mathrm{2}}^{\mathrm{2}}}小于预先设置的阈值， 则结束迭代； 否则， k=k+\mathrm{1}， 返回第二步。

输出： N\times L维重构信号\stackrel{̑}{\mathit{S}(t)}。

2.3 计算DOA的估计值

在上节输出的稀疏重构信号矩阵\stackrel{̑}{\mathit{S}(t)}\in {\mathbb{C}}^{N\times L}中只有少数行有非零值。因此， 假设存在N个入射信号， 在L个快拍下恢复的信号中， 只需要计算出矩阵\stackrel{̑}{\mathit{S}(t)}每一行的l_{\mathrm{2}}范数， 找到所有的局部最大值， 然后对这些值进行降序排列， 选取前K个， 选择的这K个值对应的行所代表的角度即为信号源的入射角度， 即DOA估计值。

3 仿真实验及结果分析

3.1 不同信噪比下各算法DOA估计误差比较

为了验证本文提出的ISTA算法的DOA估计效果， 与OMP、 CBF和MUSIC经典算法进行对比。仿真实验在Matlab平台测试。用均方根误差（Root Mean Square Error，RMSE）来评价各算法的DOA估计性能， RMSE定义如下

RMSE=\sqrt[]{\frac{\mathrm{1}}{KC}{\displaystyle \sum _{k=\mathrm{1}}^K}{\displaystyle \sum _{c=\mathrm{1}}^C}\left({\widehat{\theta }}_k\left(c\right)-{\theta }_k\right)},（21）

式中： K为信源数； C为蒙特卡洛实验次数； {\widehat{\theta }}_k(c)为第c次实验中第k个信源的DOA估计值； {\theta }_k为第k个入射信源的真实角度值。

矢量水听器采用均匀线阵， 共有30个阵元， 阵元间距为0.5 m， 快拍数为10个。入射角度范围为［0°，180°］， 信噪比范围为-10~20 dB， 以2 dB为步长递增。图 3 展示了在单个信源、 3个信源和5个信源条件下， 随着信噪比变化， ISTA、 OMP、 CBF和MUSIC算法的 RMSE曲线变化情况。其中， 单个信源的入射角度为\mathrm{10.1}°， 3个信源的入射角度为10.1°、 30.5°、 80.4°， 5个信源的入射角度为10.1°、 30.5°、 60.3°、 80.4°、 130.5°。

从图 3 可以看出， 对于单个信号源， ISTA算法在信噪比［-10 dB， 20 dB］内的RMSE均小于OMP、 CBF、 MUSIC算法。在多信源条件下， CBF算法的RMSE明显高于其他3种算法。在信源数量增多的情况下， OMP和MUSIC算法的均方根误差显著增大。综上所述， 信噪比为-10~20 dB范围内， 在小快拍的情况下， 无论在单信源还是多信源， 本文算法的RMSE递减幅度最为显著， DOA估计性能都比其他算法的估计精度高， 随着信噪比的增加， 该算法的性能优势进一步体现。

3.2 不同快拍数下各算法DOA估计误差比较

设置与上述相同的实验前提， 在信噪比为0 dB的情况下， 以步长为30从5到395设置快拍数。实验分别在单个信源、 3个信源和5个信源条件下， ISTA、 OMP、 CBF和MUSIC算法的RMSE随快拍数变化的实验结果如图 4 所示。

由图 4 可以看到， 在特定信噪比条件下， 本文算法相较于其他算法具有更低的RMSE， CBF算法的误差最大。随着快拍数的增加， 4种算法的RMSE在100快拍之前普遍呈现较大的递减趋势， 但本文算法的下降幅度更为明显。随着快拍数增加至100后， 4种算法的RMSE逐渐趋于稳定。实验结果表明， 在单信源和多信源条件下， OMP算法的RMSE始终仅次于本文算法。传统算法的估计误差相对较大， 尤其是随着信源个数的增加， DOA估计误差变得更大， 而本文算法则不受信源个数的影响。这也说明在理想的实验条件下， 拥有足够大的采样数时， 该算法仍能获得较高的DOA估计精度。总体来说， 固定信噪比为0 dB的情况下， 本文所提ISTA算法的DOA估计精度都是最高的， 且随着快拍数的增大， DOA估计性能增强， 进一步体现了该算法的稳健性， 表明它可以适用于复杂环境的DOA估计。

3.3 不同信噪比下各算法成功率比较

在本节实验中， 设置单个信源的入射角度为 20°， 3个信源的入射角度为20°、 50°和90°， 5个信源的入射角度为20°、 50°、 90°、 130°和 150°， 设置快拍数为10， 阵元数为30。信噪比为-10~20 dB， 以2 dB为步长递增。图 5 为不同信源情况下4种算法的成功概率。为了保证实验的可靠性， 在每个信噪比条件下进行100次蒙特卡罗实验。在本实验中， 成功率的判别准则为

{\displaystyle \sum _{k=\mathrm{1}}^K}\left|{\widehat{\theta }}_k\left(c\right)-{\theta }_k\right|\le \mathrm{0.25}K\mathrm{。}（22）

从图 5 可以看出， 当信噪比为-10 dB时， 在单信源及多信源中， ISTA算法的成功率都是最高的。无论是在单信源还是多信源的条件下， 随着信噪比增大， 上述所有算法的DOA估计成功率都在不断提升， 但ISTA算法的成功率总能最早达到100%。因此， ISTA算法的DOA估计准确率要远比OMP、 MUSIC和CBF算法高。由图 5（a）~图 5（c） 可以看出， 所提算法的成功率始终最高， 特别在低信噪比情况下表现更为明显。总之， 本文算法与其余3种算法相比表现更优， 验证了本文算法的优良性能。

4 结　论

本文基于压缩感知和稀疏重构理论， 提出了基于迭代阈值收缩算法的DOA估计方法， 该方法在解决矢量水听器在低信噪比、 小快拍条件下的DOA估计问题上有以下优势和特点：

1） 本文将ISTA算法成功应用于矢量水听器阵列信号DOA估计， 通过针对性的改善获得了更高的DOA估计精度。

2） 经仿真实验， 无论是在单信源还是多信源的条件下， 在小快拍、 低信噪比下DOA估计具有较高的准确率， 而且在低信噪比下， 本文算法的成功率远高于OMP、 MUSIC和CBF算法。

3） 该方法不需要对协方差矩阵进行分解， 因此具有更强的鲁棒性。

参考文献

[1] GANGULY S， GHOSH J， KUMAR P K， et al. An efficient DOA estimation and jammer mitigation method by means of a single snapshot compressive sensing based sparse coprime array［J］. Wireless Personal Communications， 2022， 123（3）： 2737-2757.

[2] ZHANG X， TAO H， FANG Z， et al. Efficient DOA estimation for wideband sources in multipath environment［J］. Remote Sensing， 2022， 14（16）： 3951.

[3] 周锋， 张文博， 张宝胜， 等. 基于空时处理的水声扩频通信［J］. 电子与信息学报， 2022， 44（6）： 2006-2013.

[4] ZHOU Feng， ZHANG Wenbo， ZHANG Baosheng， et al. Underwater acoustic spread spectrum communications based on space-time cluster processing［J］. Journal of Electronics & Information Technology， 2022， 44（6）： 2006-2013. （in Chinese）

[5] JIANG L Z. A beam domain selection method based on the conventional beam forming［J］. Applied Mechanics and Materials， 2014， 644-650： 4663-4667.

[6] SCHMIDT R. Multiple emitter location and signal parameter estimation［J］. IEEE Transactions on Antennas and Propagation， 1986， 34（3）： 276-280.

[7] ZHANG X， FENG D. An efficient MUSIC algorithm enhanced by iteratively estimating signal subspace and its applications in spatial colored noise［J］. Remote Sensing， 2022， 14（17）： 4260.

[8] YU J， LI J， SUN B， et al. Multiple RFI sources location method combining two-dimensional ESPRIT DOA estimation and particle swarm optimization for spaceborne SAR［J］. Remote Sensing， 2021， 13（6）： 1207.

[9] STOICA P， SHARMAN K C. Maximum likelihood methods for direction-of-arrival estimation［J］. IEEE Transactions on Acoustics， Speech， and Signal Processing， 1990， 38（7）： 1132-1143.

[10] WRIGHT S E. A note on computing maximum likelihood estimates for the three-parameter asymmetric Laplace distribution［J］. Applied Mathematics and Computation， 2024， 464： 128381.

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

[12] 任东升， 侯进， 党辉， 等. 压缩感知算法在阵列测向中的相干信号分辨［J］. 计算机应用研究， 2023， 40（4）： 1165-1171.

[13] REN Dongsheng， HOU Jin， DANG Hui， et al. Coherent signal resolution in array direction finding using compressed sensing algorithm［J］. Application Research of Computers， 2023， 40（4）： 1165-1171. （in Chinese）

[14] YIN Y， WANG Y， DAI T， et al. DOA estimation based on smoothed sparse reconstruction with time-modulated linear arrays［J］. Signal Processing， 2024， 214： 109229.

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

[16] 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）

[17] ŚLESICKA A， KAWALEC A. An application of the orthogonal matching pursuit algorithm in space-time adaptive processing［J］. Sensors， 2020， 20（12）： 3468.

[18] 张骄， 王敏. 基于正交匹配追踪算法的虚拟内插空间平滑DOA估计［J］. 测试技术学报， 2022， 36（5）： 369-375.

[19] ZHANG Jiao， WANG Min. DOA estimation of virtual interpolation and spatial smoothing based on orthogonal matching pursuit algorithm［J］. Journal of Test and Measurement Technology， 2022， 36（5）： 369-375. （in Chinese）

[20] GILL P R， WANG A， MOLNAR A. The In-crowd algorithm for fast basis pursuit denoising［J］. IEEE Transactions on Signal Processing， 2011， 59（10）： 4595-4605.

[21] BECK A， TEBOULLE M. A fast Iterative Shrinkage-Thresholding algorithm with application to wavelet-based image deblurring［C］//IEEE International Conference on Acoustics， Speech and Signal Processing， 2009： 693-696.

[22] CANDES E J， TAO T. Decoding by linear programming［J］. IEEE Transactions on Information Theory， 2005， 51（12）： 4203-4215.

[23] DING X， CHEN J， QIU J. A fast and differentiable optimization model for finding the minimum dominating set in large-scale graphs［J］. Journal of Computational Science， 2023， 73： 102146.

[24] 刘天娇. 稀疏贝叶斯学习理论及其在水声阵列信号处理中的应用［D］. 哈尔滨： 哈尔滨工程大学， 2020.

[25] 虞飞， 宋俊， 余赟， 等. 基于模型稀疏表示的单快拍波达方向鲁棒估计［J］. 声学技术， 2020， 39（5）： 627-631.

[26] YU Fei， SONG Jun， YU Yun， et al. Robust single-snapshot DOA estimation based on sparse representation of array output model［J］. Technical Acoustics， 2020， 39（5）： 627-631. （in Chinese）

基金资助

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

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

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

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

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

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

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

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