基于复近似信息传递算法和K近邻算法的DOA估计方法

DOA Estimation Method Based on Complex Approximate Message Passing Algorithm and K⁃Nearest Neighbour Algorithm

摘要

针对传统算法在低信噪比、 小快拍、 多信源等情况下波达方向（Direction of Arrival， DOA）估计精度低的问题， 提出了一种基于多测量向量模型的复近似信息传递算法（MMV Complex Approximate Message Passing， MCAMP）和K近邻算法（K-Nearest Neighbour， KNN）的矢量水听器阵列DOA估计方法。首先， 对空域进行等角度划分， 构造出超完备阵列流形矩阵， 建立基于稀疏表示的多快拍DOA估计模型。然后， 采用MCAMP算法进行初步估计， 保存估计结果的峰值数据， 使用KNN算法对此数据进行聚类。最后， 使用内积匹配准则选择每类信号值最大的原子， 从而得到DOA估计值。仿真实验结果表明， 与传统算法相比， 该方法具有抗噪能力强， 估计精度高等优点。

Abstract

Aiming at the problem of low estimation accuracy of direction of arrival （DOA） of traditional algorithms in low signal-to-noise ratio， small snapshots， and multiple sources， a method of DOA estimation for vector hydrophone arrays based on complex approximate message passing （MCAMP） and K-nearest neighbour （KNN） algorithm of the multiple measurement vector model was proposed. Firstly， based on the equal angle division of airspace， an overcomplete array manifold matrix was constructed， and a multi-snapshot DOA estimation model based on sparse representation was established. Then， the MCAMP algorithm was used for preliminary estimation， and the peak data of the estimation results were saved and clustered using the KNN algorithm. Finally， the inner product matching criterion was used to select the atom with the largest signal value in each class to obtain the DOA estimation. The results of simulation experiments show that compared with the traditional algorithm， the method has the advantages of strong anti-noise capability and high estimation accuracy.

Graphical abstract

关键词

DOA估计 / 压缩感知 / 复近似信息传递算法 / K近邻算法 / 矢量水听器

Key words

DOA estimation / compressed sensing / complex approximate message passing algorithm / K-nearest neighbour algorithm / vector hydrophone

0 引　言

波达方向估计（Direction of Arrival， DOA）是阵列信号的研究热点之一， 它在5G和物联网通信^［1］、 地震监测^［2］和水声信号处理^［3］等众多领域中的应用极大地方便了人们的生产和生活。早期的DOA估计算法， 例如常规波束形成（Conventional Beamforming， CBF）^［4］， 直接对应于传统的傅里叶谱估计方法， 使得阵列角度分辨率受到瑞利限的限制， 不适合一些精度要求较高的应用场景。为了提高算法的估计精度， 学者们提出了两大类算法。一类是子空间分解类算法， 例如： 多重信号分类（Multiple Signal Classification， MUSIC）^［5］， MUSIC算法利用信号子空间与噪声子空间的正交性， 构建空间谱函数， 通过谱峰搜索实现DOA估计； 而旋转不变子空间算法（Estimation of Signal Parameters using Rotational Invariance Techniques， ESPRIT）^［6］利用各子阵的信号子空间的旋转不变性实现了无需谱峰搜索的DOA估计。但是， 这类算法在小快拍数和相干信号环境下的估计性能较差。另一类是子空间拟合类算法， 例如极大似然算法（Maximum Likelihood， ML）^［7］， 需要对所有信号来向进行联合搜索， 严重影响了工作效率。

Donoho^［8］提出了压缩感知理论， 该理论利用信号的稀疏性， 通过重构算法即可获得高精度的重构信号。与子空间算法相比， 该类算法不仅适用于相干信号， 在小快拍和低信噪比情况下仍具有良好的估计性能。空间信号的自然稀疏性使得DOA估计可以作为压缩感知理论框架下的稀疏重构问题来解决。

典型的稀疏重构算法有贪婪算法、 凸优化算法以及迭代阈值算法。其中， 贪婪算法主要是通过迭代选择原子和计算残差来实现信号重构， 如正交匹配追踪算法（Orthogonal Matching Pursuit， OMP）属于贪婪算法^［9］， 这类算法迭代次数少， 计算复杂度低， 但是重构精度不高。凸优化算法则是将非凸优化问题转化为凸优化问题进而求解目标信号， 如L1-SVD方法（L1-Singular Value Decomposition）^［10］。相比贪婪算法， 凸优化算法重构时产生的误差较小。但是， 凸优化类方法的时间复杂度较大。迭代阈值算法相比凸优化算法有更快的重构速度， 这类算法通过硬阈值和软阈值重构^{［11‐13］}， 然而在性能上不如凸优化算法。近似信息传递算法（Approximate Message Passing， AMP）^［14］是由Donoho等提出的新型稀疏重构算法， 它不仅具有迭代阈值类算法在计算复杂度上的优势， 还具有基追踪类算法强大的信号重构能力。然而， AMP算法主要用于实值信号处理。复近似信息传递算法（Complex Approximate Message Passing， CAMP）^［15］将AMP算法从实值域扩展到复值域， CAMP算法在求解复值最小绝对收缩和选择算子（Complex-Value Least Absolute Shrinkage and Selection Operator， c-LASSO）问题时， 避免了分别处理实部和虚部， 相比AMP算法具有更好的重构性能。

目前， 基于压缩感知算法的DOA估计主要针对SMV（Single Measurement Vector， SMV）模型， 其在低信噪比下的鲁棒性较差， 重构效率较低。相比之下， 多测量向量（Multiple Measurement Vector， MMV）模型更符合实际DOA估计的应用场景。此外， 压缩感知算法会受到感知矩阵相干性影响， 而DOA估计问题中的导向矩阵是高度相干矩阵， 直接使用CAMP算法不能很好地实现DOA估计。

因此， 为了提升小快拍数以及低信噪比条件下的DOA估计精度， 同时鉴于DOA估计处理的对象为复数信号， 本文提出了一种基于MMV模型的复近似信息传递算法（MMV Complex Approximate Message Passing， MCAMP）和K近邻算法（K-Nearest Neighbour Algorithm， KNN）相结合的DOA估计方法。该方法首先将CAMP算法推广至MMV模型下的复近似信息传递算法， 并进行初步估计。接着， 在MCAMP算法重构结果的支撑集中进一步寻找更精确的解。在相同的信噪比、 快拍数的仿真实验中， 该算法不需要预先确定信源数就可以得到高精度的DOA估计值， 比传统的DOA方法更具优势。

1 阵列信号模型

1.1 矢量水听器阵列信号模型

考虑t时刻在M个间距为d的均匀线性矢量水听器阵列上， 接收K个远场窄带信号， 这些信号的入射角度为{\theta }_k， k=\mathrm{1,2},\cdots ,K， 如图 1 所示。

图 1 中每个矢量水听器由两个振速通道和一个声压通道组成， t时刻声压和振速的关系为

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

式中： \mathit{\xi }， \mathit{\eta }为相互正交的单位向量， 则声压与两个振速通道的接收数据模型为

\left\{\begin{array}{l}p(t)=x(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，\end{array}\right.（2）

式中： p(t)、 v_x(t)、 v_y(t)分别为t时刻的声压通道以及位于x轴和y轴方向的振速通道接收的数据。

综上可得， 单个矢量水听器输出的声压与振速三路信号为

\begin{array}{c}\mathit{y}(t)={\left[p(t),v_x(t),v_y(t)\right]}^{\mathrm{T}}=\\ x(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{Y}(t)=\mathit{A}(\theta )\mathit{S}(t)+\mathit{N}(t)，（4）

式中： \mathit{Y}(t)={\left[y_{\mathrm{1}}(t),y_{\mathrm{2}}(t),\cdots ,y_M(t)\right]}^{\mathrm{T}}， \mathit{Y}(t)∈{\mathbb{C}}^{\mathrm{3}M\times \mathrm{1}}， 表示矢量水听器阵列输出数据； \mathit{A}(\theta )=\left[\begin{array}{c}\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)\end{array}\right],\mathit{A}(\theta )\in {\mathbb{C}}^{\mathrm{3}M\times K}， 表示矢量水听器阵列的信号方向矩阵； \mathit{a}({\theta }_k)=[\mathrm{1},{\mathrm{e}}^{-\mathrm{j}\mathrm{2}\mathrm{\pi }d\frac{\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k}{\lambda }},\cdots ,{\mathrm{e}}^{-\mathrm{j}\mathrm{2}\mathrm{\pi }d\frac{\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k}{\lambda }}{]}^{\mathrm{T}}表示第k个声源的声压响应系数， \lambda表示波长； \mathit{u}({\theta }_k)=[\mathrm{1},\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k,\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k{]}^{\mathrm{T}}， k=\mathrm{1,2},\cdots ,K， 表示阵列流形矩阵扩展向量； \otimes表示克罗内克积； \mathit{S}(t)={\left[\begin{array}{c}{s}_{\mathrm{1}}(t),s_{\mathrm{2}}(t),\cdots ,s_K(t)\end{array}\right]}^{\mathrm{T}}， \mathit{S}(t)\in {\mathbb{C}}^{K\times \mathrm{1}}， 表示t时刻的信源矢量； \mathit{N}(t)={\left[n_{\mathrm{1}}(t),n_{\mathrm{2}}(t),\cdots ,n_M(t)\right]}^{\mathrm{T}}， \mathit{N}(t)\in {\mathbb{C}}^{\mathrm{3}M\times \mathrm{1}}， 表示加性噪声矢量。

1.2 基于压缩感知的多快拍DOA估计模型

下面考虑多个快拍的情况， L个快拍且含噪情况下阵列的输出模型为

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

式中： \overline{\mathit{Y}(t)}=[{\mathit{Y}}_{\mathrm{1}}(t),{\mathit{Y}}_{\mathrm{2}}(t),\cdots ,{\mathit{Y}}_L(t)]\in {\mathbb{C}}^{\mathrm{3}M\times L}， 表示L个快拍下阵列的输出数据； \overline{\mathit{S}(t)}={\left[\begin{array}{c}{\mathit{S}}_{\mathrm{1}}(t),{\mathit{S}}_{\mathrm{2}}(t),\cdots ,{\mathit{S}}_L(t)\end{array}\right]}^{\mathrm{T}}\in {\mathbb{C}}^{K\times L}， 表示K\times L维信源矢量； \overline{\mathit{N}(t)}=[{\mathit{N}}_{\mathrm{1}}(t)，{\mathit{N}}_{\mathrm{2}}(t)，\cdots，{\mathit{N}}_L(t)]∈{\mathbb{C}}^{\mathrm{3}M\times L}。

阵列流形矩阵\mathit{A}的每一列对应一个入射声源的空间角度， 直接对信号的可能空间角度［0°，180°］进行网格划分， 可以得到扩展后的阵列流形矩阵， 进而构建DOA估计的稀疏表示模型。划分后的空间稀疏模型如图 2 所示。

在图 2 中， 整个空域角被划分为N份。为了体现目标声源的稀疏性， 离散网格的个数会比实际信号源的个数大很多， 也就是N\gg K。在这样构造的空间稀疏信号中， 在\stackrel{̑}{\mathit{S}(t)}中只有K个实际存在信号的位置的元素是非零的。在空间网格划分后， 导向矢量矩阵由\mathrm{3}M\times K阶扩展为\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 ,\right.

\left.\mathit{a}({\theta }_N)\otimes \mathit{u}({\theta }_N)\right]\mathrm{。}（6）

与此同时， 基于压缩感知的多快拍阵列输出模型变为

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

通过上述估计模型的构建， 求解式（7）的稀疏重构问题就可以得到DOA估计， 即给定观测信号\stackrel{̑}{\mathit{Y}(t)}和感知矩阵\mathit{G}(\theta )， 重构目标稀疏信号\stackrel{̑}{\mathit{S}(t)}。可以将其转化为最小化{\mathrm{𝓁}}_{\mathrm{0}}范数优化问题进行求解， 即

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

式中： {\epsilon }^{\mathrm{2}}为与噪声有关的常量。

式（8）在本质上是具有非凸稀疏约束的优化问题， 并且已经被证明是计算复杂度很高的NP-hard问题， 通常将其转化为{\mathrm{𝓁}}_{\mathrm{1}}范数优化问题， 即

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

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

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

式中： \lambda为正则化参数， \lambda >\mathrm{0}， 当噪声强度不同时， 其值也不同。

最小化{\mathrm{𝓁}}_{\mathrm{1}}范数优化问题可以将其转化为二阶锥模型^［16］进行求解， 利用Matlab中的CVX工具箱， 求解出\widehat{\mathit{S}(t)}之后进行谱峰搜索， 所得峰值位置就是DOA估计值。这一过程虽然能够实现精确重构但影响了计算效率， 计算复杂度较高， 因此不适用于大规模信号DOA估计。针对这一问题， 本文提出了一种迭代阈值算法和K近邻算法相结合的方法， 以解决基于压缩感知模型的DOA估计问题。

2 基于MCAMP算法和KNN算法的DOA估计方法

2.1 MCAMP算法

MCAMP算法的关键步骤主要有两个： 1） 通过对式（10）进行向量化处理， 将DOA估计的无约束优化问题的求解转化为c-LASSO问题； 2） 基于矩阵乘积的向量化与Kronecker积的关系， 将CAMP算法推广至求解多快拍情况下的DOA估计的基本迭代过程。

令\stackrel{̑}{\mathit{y}}=vec(\stackrel{̑}{\mathit{Y}})， \stackrel{̑}{\mathit{G}}=（{\mathit{I}}_N\otimes \mathit{G})， \stackrel{̑}{\mathit{s}}=vec(\stackrel{̑}{\mathit{S}})， 则式（10）中的无约束优化问题可以转化为如式（11）的c-LASSO问题。

\widehat{\mathit{s}}=\underset{\stackrel{̑}{\mathit{s}}}{\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖\stackrel{̑}{\mathit{y}}-\stackrel{̑}{\mathit{G}}\stackrel{̑}{\mathit{s}}‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda {‖\stackrel{̑}{\mathit{s}}‖}_{\mathrm{1}}，（11）

式中： vec(·)为矩阵的向量化操作， 把\stackrel{̑}{\mathit{S}}按列堆叠成向量\stackrel{̑}{\mathit{s}}； \otimes为Kronecker积； {\mathit{I}}_N\in {\mathbb{R}}^{N\times N}为单位矩阵。

上述c-LASSO问题即可由CAMP算法进行求解。CAMP算法迭代公式为

\left\{\begin{array}{l}{\tilde{\mathit{s}}}^t={\stackrel{̑}{\mathit{G}}}^H{\mathit{z}}^{t-\mathrm{1}}+{\widehat{\mathit{s}}}^{t-\mathrm{1}},\\ {\widehat{\mathit{s}}}^t=\eta ({\tilde{\mathit{s}}}^t;\gamma {\sigma }^t),\\ \begin{array}{l}{\mathit{z}}^t=\stackrel{̑}{\mathit{y}}-\stackrel{̑}{\mathit{G}}{\widehat{\mathit{s}}}^t+{\mathit{z}}^{t-\mathrm{1}}·\\ \frac{\mathrm{1}}{\mathrm{2}\delta }\left[〈\partial {\eta }_u^{\mathrm{R}\mathrm{e}}({\tilde{\mathit{s}}}^t;\gamma {\sigma }^t〉+〈\partial {\eta }_v^{\mathrm{I}\mathrm{m}}({\tilde{\mathit{s}}}^t;\gamma {\sigma }^t)〉\right],\end{array}\end{array}\right.（12）

其中，

\eta (u+vi;\tau )=\left\{\begin{array}{ll}u+vi-\frac{\tau (u+vi)}{\sqrt[]{u^{\mathrm{2}}+v^{\mathrm{2}}}},& u^{\mathrm{2}}+v^{\mathrm{2}}>{\tau }^{\mathrm{2}},\\ \mathrm{0},& u^{\mathrm{2}}+v^{\mathrm{2}}\le {\tau }^{\mathrm{2}}，\end{array}\right.（13）

式中： \eta (u+vi;\tau )为复数软阈值函数； {\tilde{\mathit{s}}}^t为第t次迭代信号的噪声估计信号； {\widehat{\mathit{s}}}^t为第t次迭代后的重构信号； {\mathit{z}}^t为第t次迭代后的残差； \delta =M/N为测量速率； {\eta }_{}^{\mathrm{R}\mathrm{e}}和{\eta }^{\mathrm{I}\mathrm{m}}分别为\eta的实部和虚部， \partial {\eta }_u^{\mathrm{R}\mathrm{e}}为{\eta }^{\mathrm{R}\mathrm{e}}关于实部u的偏导数， \partial {\eta }_v^{\mathrm{I}\mathrm{m}}为{\eta }^{\mathrm{I}\mathrm{m}}关于虚部v的偏导数； 〈·〉为向量的平均值； 参数\gamma与正则化参数\lambda等效^［17］， 这种关系使得CAMP算法能够准确解决c-LASSO问题， 噪声标准差{\sigma }^t可使用{\sigma }^t=\sqrt[]{\mathrm{1}/\mathrm{l}\mathrm{n}\mathrm{2}}\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{n}(\left|{\tilde{\mathit{s}}}^t\right|)进行估计。

由于考虑的是多快拍情况下的DOA估计问题， 故本文根据矩阵乘积的向量化与Kronecker积的关系， 依据式（14）将迭代过程进行了矩阵化处理。

\left\{\begin{array}{l}\begin{array}{l}{\stackrel{̑}{\mathit{G}}}^H{\mathit{z}}^t={\left({\mathit{I}}_N\otimes \mathit{G}\right)}^H{\mathit{z}}^t=\left({\mathit{I}}_N\otimes {\mathit{G}}^H\right)vec\left({\mathit{Z}}^t\right)=\\ vec\left({\mathit{G}}^H{\mathit{Z}}^t{\left({\mathit{I}}_N\right)}^{\mathrm{T}}\right)=vec\left({\mathit{G}}^H{\mathit{Z}}^t\right),\end{array}\\ \begin{array}{l}\stackrel{̑}{\mathit{G}}{\widehat{\mathit{s}}}^t=({\mathit{I}}_N\otimes \mathit{G}){\widehat{\mathit{s}}}^t=({\mathit{I}}_N\otimes \mathit{G})vec\left({\widehat{\mathit{S}}}^t\right)=\\ vec\left(\mathit{G}{\widehat{\mathit{S}}}^t{\left({\mathit{I}}_N\right)}^{\mathrm{T}}\right)=vec\left(\mathit{G}{\widehat{\mathit{S}}}^t\right)\mathrm{。}\end{array}\end{array}\right.（14）

可以得到矩阵化形式

\left\{\begin{array}{l}{\tilde{\mathit{S}}}^t={\mathit{G}}^H{\mathit{Z}}^{t-\mathrm{1}}+{\widehat{\mathit{S}}}^{t-\mathrm{1}},\\ {\widehat{\mathit{S}}}^t=\eta ({\tilde{\mathit{S}}}^t;\gamma {\sigma }^t),\\ \begin{array}{l}{\mathit{Z}}^t=\stackrel{̑}{\mathit{Y}}-\mathit{G}{\widehat{\mathit{S}}}^t+{\mathit{Z}}^{t-\mathrm{1}}\frac{\mathrm{1}}{\mathrm{2}\delta }·\\ \left[〈\partial {{\eta }_u}^{\mathrm{R}\mathrm{e}}({\tilde{\mathit{S}}}^t;\gamma {\sigma }^t〉+〈\partial {{\eta }_v}^{\mathrm{I}\mathrm{m}}({\tilde{\mathit{S}}}^t;\gamma {\sigma }^t)〉\right]\mathrm{。}\end{array}\end{array}\right.（15）

综上， 输入观测矩阵\stackrel{̑}{\mathit{Y}(t)}和感知矩阵\mathit{G}(\theta )， 通过式（15）的迭代过程， 可以得到式（10）的解， 即得到重构信号\widehat{\mathit{S}(t)}。

2.2 基于K近邻的聚类方法

压缩感知算法的性能会受到感知矩阵相干性和噪声环境的影响， 而基于压缩感知的DOA估计算法能够大概率地在真实角度邻近区域得到估计结果， 这使得在真实角度邻近区域能够收集到更多的峰值信息^［18］。因此， 通过对这些峰值数据的综合分析， 能够明确基于压缩感知所得的DOA估计值主要集中于哪些特定角度范围， 从而有效确定信源数量并精确估计DOA角度。考虑利用KNN聚类方法^［19］对MCAMP-DOA得到的峰值数据进行处理， 通过将聚类分析理论与MCAMP算法结合， 对MCAMP算法DOA估计的结果进行进一步的聚类， 将大大提高MCAMP算法的DOA估计精度。

基于KNN的聚类方法具体流程如下： 首先， 初始化一个空的聚类集合\mathit{C}， 输入样本集合\mathit{J}和最大距离d_max。接着， 遍历每个原子， 计算各类中心欧氏距离d， 如果d>d_max， 创建新类， 否则将当前原子加入到最近邻类。最后 ， 得到聚类簇\mathit{C}=\left\{C_{\mathrm{1}},C_{\mathrm{2}},\cdots ,C_K\right\}。

2.3 基于MCAMP算法与KNN算法的DOA估计方法

本文MCAMP-KNN算法是将MCAMP算法和KNN算法进行了结合。具体而言， 首先， MCAMP算法通过阵列阵元接收到的信号数据， 并运用稀疏重构理论， 确定信号源最可能存在的角度； 同时， MCAMP算法可被视为一个初步的DOA估计的过程， 用以确定大概的到达方向。然后， 将MCAMP算法得到的解支撑集作为输入的数据集， 利用KNN算法进行聚类， 确定聚类簇的数量， 即信源数。最后， 使用OMP算法内积匹配准则选择每类对应信号值最大的原子， 得到重构信号， 将重构信号中的非零元素所在的位置转换为对应的信号源入射角度， 即得到DOA估计值。使用MCAMP-KNN算法， 不仅可以准确估计真实信号源的数量， 还可以在低信噪比和小快拍条件下明显提高DOA估计的准确度。综上， 本文MCAMP-KNN算法流程如图 3 所示。

算法具体如下：

算法1　基于MCAMP算法与KNN算法的DOA估计方法

输入： M\times L维观测矩阵\stackrel{̑}{\mathit{Y}(t)}， M\times N感知矩阵\mathit{G}(\theta )， 最大迭代次数itermax， 最大距离dmax， 参数P（P可根据数据量适当调整）；

步骤1： 使用式（15）进行初步估计， 保存估计结果的前P个最大的峰值数据；

步骤2： 对步骤1中保存的P个峰值数据进行基于K近邻的聚类， 得到聚类簇数目K， 即信源数；

步骤3： 对每类数据进行去噪处理， 筛选出部分噪声数据；

步骤4： 使用OMP算法的内积匹配准则， 选择每类中对应信号值最大的原子， 逐步扩充支撑集， 得到重构信号\widehat{\mathit{S}(t)}；

步骤5： 将重构信号\widehat{\mathit{S}(t)}中的非零元素所在的位置转换为对应的信号源入射角度， 即得到正确DOA估计值。

输出： DOA估计值。

3 仿真结果及分析

为了验证与评估本文所提算法的估计性能， 实验同时执行了本文所提MCAMP-KNN算法、 MCAMP、 OMP以及MUSIC和CBF算法。在以下实验中， 矢量均匀线性阵列阵元数为20， 阵元间距为0.5 m， 角度搜索范围为［0°，180°］， 噪声为高斯白噪声。选用均方根误差（Root Mean Square Error， RMSE）来评估DOA估计性能， 定义为

E_{\mathrm{R}\mathrm{M}\mathrm{S}}=\sqrt[]{\frac{\mathrm{1}}{K\cdot T}{\displaystyle \sum _{c=\mathrm{1}}^T}{\displaystyle \sum _{r=\mathrm{1}}^K}{\left[{\widehat{\theta }}_r(c)-{\theta }_r\right]}^{\mathrm{2}}}，（16）

式中： T为实验的次数； K为信源的数量； {\widehat{\theta }}_r(c)为第c次蒙特卡洛实验中对第r个目标的DOA估计结果； {\theta }_r为第r个目标的实际方向角。仿真实验在Matlab操作平台上进行。

3.1 各算法在低信噪比和小快拍下的空间谱图

在本实验中， 设置快拍数均为10， 信噪比均为-10 dB， 信号源入射角度为30.4°， 60.2°和120.8°， 信号频率分别为1 000， 1 200和1 400 Hz。在此设置下比较MCAMP算法、 OMP算法、 MUSIC算法、 CBF算法与本文所提算法的空间谱图， 实验结果如图 4 所示。

由图 4 可以看出， 在低信噪比和小快拍条件下， MUSIC算法和CBF算法存在大量伪峰， 估计性能最低。MCAMP算法虽然也存在伪峰， 但其主峰峰值都远大于伪峰的峰值， 相比OMP算法估计精度较高。本文所提算法相比MCAMP算法无大量伪峰， 估计精度高。综上所述， 本文所提算法具有较高的估计精度， 能够保证在低信噪比和小快拍数条件下的DOA估计性能。

为了直观地说明本文所提算法在低信噪比和小快拍环境中的优势， 通过表 1 展现了不同算法的DOA估计值。

由表 1 可以看出， 只有MCAMP-KNN算法准确估计出信号2的入射角度， 其估计值为60.2°， 与真实值60.2°相同。

综上可知， 本文所提算法在低信噪比和小快拍情况下估计的准确度最高， 其次是MCAMP算法。其余的算法相比本文所提算法估计精度低， 误差较大。实验结果充分验证了本文所提算法在DOA估计上具有更高的精度。

3.2 各算法在不同信噪比下的DOA估计误差

在本实验中， 设置快拍数为10， 信噪比的变化范围为-10~10 dB， 步进大小为2 dB， 针对每个信噪比， 执行100次蒙特卡洛实验。图 5 展示了1~3个信源时， 本文所提MCAMP-KNN算法与MCAMP算法、 OMP算法、 MUSIC算法和CBF算法的DOA估计均方根误差随着信噪比变化的趋势。

在1个信号源情况下， 入射角度设置为9.4°， 信号频率为1 000 Hz； 在2个信号源情况下， 入射角度设置为9.4°， 29.8°， 信号频率分别为1 000， 1 200 Hz； 在3个信号源情况下， 入射角度设置为9.4°， 29.8°和60.2°， 信号频率分别为1 000， 1 200和1 400 Hz。

由图 5 可以看出， 在信噪比变化范围内， 随着信噪比的增加， 不同算法的均方根误差都不断减小， 即DOA估计性能越来越强。在信噪比为-10 dB时， 即使信源数增加， 本文提出的MCAMP-KNN算法相较于MCAMP算法， 其均方根误差依然保持较低水平， 这归因于本文所提算法是在MCAMP算法的解支撑集中进一步寻找更精确的解， 因此， 相比MCAMP， 本算法具有更低的均方根误差。同时， MCAMP算法的均方根误差也低于OMP算法、 MUSIC算法和CBF算法。另外， 信噪比在-10 dB至10 dB的变化范围内， 本文所提算法的均方根误差都是最低的。综上所述， 本文所提算法无论是在单信源环境还是多信源环境， 在DOA估计过程中都展现出更强的鲁棒性， 能够有效地减小噪声对DOA估计结果的影响， 显著提升DOA估计的准确性。

3.3 各算法在不同快拍数下的DOA估计误差

在本实验中， 设置与上节相同的实验条件， 所不同的是设置信噪比为-5 dB， 快拍数从2增加到102， 步进大小为10。图 6 展示了在1~3个信源情况下， 本文所提算法与MCAMP算法、 OMP算法、 MUSIC算法和CBF算法在DOA估计中的均方根误差随快拍数变化的关系曲线。

由图 6 可以看出， 不同算法的均方根误差都随快拍数的增加而减小。值得注意的是， 在快拍数的变化范围内， 即使信源数增加， 传统CBF算法、 MUSIC算法和OMP算法的DOA估计精度都存在较大误差， 而MCAMP-KNN算法的DOA估计均方根误差明显小于MCAMP算法， 这主要是由于本文所提算法是在MCAMP算法的支撑集中寻找更精确的解。综上所述， 无论是单信源还是多信源实验， 本文所提算法的估计性能都优于其它算法， 并且在小快拍的情况下也能够保持良好的估计性能。

3.4 各算法的估计成功概率随信噪比和快拍数变化的情况

在本实验中， 设3个信号分别从20.3°， 50.2°和90.8°入射到阵元数为20的矢量水听器上。 1） 信噪比仿真实验条件： 信噪比从-10~10 dB变化， 步进大小为2 dB， 快拍数为10。2） 快拍数仿真实验条件： 快拍数从2增加到102， 步进大小为10， 信噪比为-5 dB。图 7 分别展示了在不同信噪比和快拍数下， 本文所提算法与MCAMP算法、 OMP算法、 MUSIC算法和CBF算法的成功率对比情况。在本实验中， 估计误差小于或等于0.5°时则视为估计成功。

由图 7 可以看出， 随着信噪比和快拍数的增大， 各算法的DOA估计成功率都逐渐增大， 并且本文所提算法的成功率都是率先达到100%。由图 7（a） 可知， 在信噪比为-10 dB时， 本文所提算法成功率最高， 在整个信噪比的变化区间内， 本文算法都比其它算法成功率高。由图 7（b） 可知， 在快拍数为2时， 本文所提算法DOA估计成功率超过其它算法， 在整个快拍数的变化区间内， 本文算法都比其它算法成功率高。综上所述， 在低信噪比和小快拍数情况下， 本文所提算法都优于其余4种算法， 在DOA估计中更具优势。

3.5 各算法的均方根误差和估计成功概率随阵元数变化的情况

在本实验中， 设置与上节相同的实验条件， 所不同的是阵元数从10~60变化， 步进大小为5， 信噪比为0 dB， 快拍数为20。图 8 展示了在不同阵元数情况下， 本文所提算法与MCAMP算法、 OMP算法、 MUSIC算法和CBF算法在DOA估计中的均方根误差和成功率的对比情况。

由 8（a） 可知， 随着阵元数的增加， 各算法的均方根误差都出现了较大的递减趋势， 这主要是由于阵元数增加， 阵列接收矩阵中包含的有效信息也随之增加， 从而减小了DOA估计的均方根误差。同时， 在阵元数大于15之后， 本文所提算法的均方根误差都是最低的。由图 8（b） 可知， 不同算法的估计成功率随着阵元数的增加不断增加， 并且本文所提算法率先达到100%。在阵元数为10时， 本文所提算法的DOA估计成功率略低于MCAMP算法， 在阵元数超过15之后， 本文算法的成功率都高于其它算法。综上所述， 对于更大规模的阵列， 本文所提MCAMP-KNN算法也能保持很高的成功率以及较低的均方根误差， 这为其在实际大规模阵列系统中的应用提供了有力支撑。

4 结　论

本文基于压缩感知的DOA估计模型， 提出了一种结合基于MMV模型的复近似信息传递算法和K近邻算法的矢量水听器DOA估计方法， 该算法有以下特点和优势：

1） 本文借助矩阵乘积的向量化与Kronecker积的关系， 将CAMP算法推广至求解多快拍情况下的DOA估计值， 相较于传统算法， 其估计性能有一定的提升；

2） 由于基于压缩感知的DOA估计算法能够大概率地在真实角度邻近区域得到估计结果， 且本文在初始估计结果的支撑集中进一步寻找更精确的解， 因此相较于原本的CAMP算法估计性能上有一定的提升；

3） 利用压缩感知进行DOA估计， 在单信源和多信源条件下， 矢量水听器在小快拍、 低信噪比下也可估计出信号源的方向角度；

4） 该算法不需要预先确定信源数目， 就可以准确估计， 具有良好的实用价值。

参考文献

[1] FABREGAT G， BELLOCH J A， BADÍA J M， et al. Design and implementation of acoustic source localization on a low-cost IoT edge platform［J］. IEEE Transactions on Circuits and Systems II： Express Briefs， 2020， 67（12）： 3547-3551.

[2] KATZ N Z， KHATIB M， BEN-HORIN Y， et al. Geometry design for DOA estimation in seismic 2D-arrays： Simulation study［J］. IEEE Access， 2024， 12： 35827-35843.

[3] YAN H， CHEN T， WANG P， et al. A direction-of-arrival estimation algorithm based on compressed sensing and density-based spatial clustering and its application in signal processing of MEMS vector hydrophone［J］. Sensors， 2021， 21（6）： 2191.

[4] XU Z， LI H， YANG K. A modified differential beamforming and its application for DOA estimation of low frequency underwater signal［J］. IEEE Sensors Journal， 2020， 20（16）： 8890-8902.

[5] 王绪虎， 侯玉君， 金序， 等. 基于协方差矩阵Toeplitz重构的MUSIC算法的波达方位估计方法［J］. 声学技术， 2023， 42（6）： 825-831.

[6] WANG Xuhu， HOU Yujun， JIN Xu， et al. DOA estimation method using MUSIC algorithm based on Toeplitz reconstruction of covariance matrix［J］. Technical Acoustics， 2023， 42（6）： 825-831. （in Chinese）

[7] VEERENDRA D， BALAMURUGAN K S， VILLAGÓMEZ-GALINDO M， et al. Optimizing sensor array DOA estimation with the manifold reconstruction unitary ESPRIT algorithm［J］. IEEE Sensors Letters， 2023， 7（12）： 7006804.

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

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

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

[11] 窦慧晶， 肖子恒， 杨帆. 基于改进稀疏度自适应匹配追踪算法的压缩感知DOA估计［J］. 北京工业大学学报， 2021， 47（11）： 1239-1246.

[12] DOU Huijing， XIAO Ziheng， YANG Fan. Compressed sensing DOA estimation based on improved sparsity adaptive matching pursuit algorithm［J］. Journal of Beijing University of Technology， 2021， 47（11）： 1239-1246. （in Chinese）

[13] 赵季红， 马兆恬， 曲桦， 等. 基于加权1范数稀疏信号重建的DoA估计［J］. 北京邮电大学学报， 2016， 39（5）： 33-36.

[14] ZHAO Jihong， MA Zhaotian， QU Hua， et al. DoA estimation based on sparse signal recovery utilizing weighted 1 norm［J］. Journal of Beijing University of Posts and Telecommunications， 2016， 39（5）： 33-36. （in Chinese）

[15] BLUMENSATH T， DAVIES M E. Iterative hard thresholding for compressed sensing［J］. Applied and Computational Harmonic Analysis， 2009， 27（3）： 265-274.

[16] DAUBECHIES I， DEFRISE M， DE MOL C. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint［J］. Communications on Pure and Applied Mathematics， 2004， 57（11）： 1413-1457.

[17] BECK A， TEBOULLE M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems［J］. SIAM Journal on Imaging Sciences， 2009， 2（1）： 183-202.

[18] DONOHO D L， MALEKI A， MONTANARI A. Message-passing algorithms for compressed sensing［J］. Proceedings of the National Academy of Sciences of the United States of America， 2009， 106（45）： 18914-18919.

[19] LUO X， YU L， LI M， et al. Complex approximate message passing equivalent source method for sparse acoustic source reconstruction［J］. Mechanical Systems and Signal Processing， 2024， 217： 111476.

[20] LOBO M S， VANDENBERGHE L， BOYD S， et al. Applications of second-order cone programming［J］. Linear Algebra and Its Applications， 1998， 284（1/2/3）： 193-228.

[21] MALEKI A. Approximate message passing algorithms for compressed sensing［D］. Stanford： Stanford University， 2010.

[22] 朱宝灼. 基于FOCUSS算法的阵列分析综合及DOA估计研究［D］. 成都： 电子科技大学， 2022.

[23] 宋柳. 基于信号稀疏表示的轴承故障诊断研究［D］. 南京： 东南大学， 2020.

基金资助

山西省留学回国人员科技活动择优资助项目(20240011)

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

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

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