基于迂回式匹配追踪算法的DOA估计方法

黄嘉俊，张靖奇，何伟嘉，王 鹏

(中北大学 理学院，山西 太原 030051)

0 引 言

水声信号的波达方向(Direction of arrival，DOA)估计一直是水声信号处理的重要研究领域.早期的DOA估计算法主要围绕波束形成技术开展，以Capon算法[1]为代表的波束形成算法，通过对传感器阵列接收到的各路信号进行加权处理，抑制非目标方向的干扰信号，增强目标方向的期望信号来获得空间谱.然而，该技术无法在低信噪比和高混响环境中准确估计波达方向.随着子空间类算法的出现，波达方向估计的分辨率得到了极大的提高.一类是子空间分解类算法，包括多重信号分类(Multiple signal classification,MUSIC)算法[2]和旋转不变子空间(Estimation of signal parameters via rotation invariant technique,ESPRIT)算法[3].这类算法利用信号子空间与噪声子空间的正交特性来构造空间谱，突破了DOA估计中分辨率的瑞利限，实现了超分辨率的DOA估计.但是子空间分解类算法由于依赖阵列接收信号的统计特性，需要大量的快拍数据，难以实现小快拍数据下精确的空间谱估计，此外，这类方法无法有效估计相干信号的空间谱.另一类是子空间拟合类DOA估计算法，其中具有代表性的是最大似然(Maximum likelihood，ML)算法[4]，该方法通过求解非线性方向估计似然函数来获得空间谱，由于需要进行多维搜索来找寻最优解，其运算量较大.

近年来，Donoho D、Candes E和Tao T提出了一种新理论—压缩感知(Compressive sensing，CS)[5].该方法可以在确保信息完整的情况下，以小于奈奎斯特采样率的频率对信号进行采样.其过程主要是：首先利用信号的稀疏性，以远小于奈奎斯特采样率的随机采样率对稀疏信号进行采样，然后通过重构算法较好地重构原始信号.在实际空域中，实际目标信号数相对于整个空间潜在的辐射信号数来说是极少的，这使得DOA估计中的空域信号本身即满足稀疏性.因此，基于压缩感知原理的DOA估计算法是一种利用信号稀疏性来实现高分辨率估计的方法[6].由于其在单快拍以及信源信号相干的条件下仍能获得良好的DOA估计而受到广泛的关注.

学者们在此基础上提出了许多信号重构的方法，主要包括凸优化算法以及贪婪算法[7-8]等.基于凸优化类方法的信号重构算法，包括奇异值分解(Singular value decomposition，SVD)[9]、基追踪降噪(Basis pursuit de-noising，BPDN)[10]和硬阈值迭代(Iterative hard thresholding，IHT)[11]等算法.该类算法是一种全局优化算法，通过将非凸问题转化为凸问题求解找到待重构信号的近似值.文献[12]给出了一种基于同伦思想的基追踪问题解法，该方法通过迭代得出最有效的规整因子和重构结果，适合数据量较少、无法确定误差界限的场景，联系宽带信号的联合稀疏性，该方法被扩展到单个频域快拍时的宽带信号测向中.文献[13]提出使用基追踪降噪算法进行DOA估计，基于压缩感知的DOA估计模型的求解为最小l0范数问题.在一定条件下，l1范数与l0范数最小化问题是等价的[5]，该算法的思路是将l0范数最小化问题转换成l1范数最优化问题，以便使用最优化方法求解.基追踪降噪算法可以对相干信号进行DOA估计，具有较高的稳定性和重构精度，但是其最终求解的无约束优化问题涉及到正则化参数的选择，正则化参数的大小影响着最终重构效果对噪声的适应性，需要选择合适的正则化参数来保证算法性能.迂回式匹配追踪(Detouring matching pursuit，DMP)算法[14]是由裴廷睿等提出的一种贪婪重构算法.该算法严格证明了子内积逆和系数矩阵递增递减核心式，利用子内积逆和系数矩阵减少计算残差误差变化量的计算量，降低了计算复杂度，具有重构准确率高、计算复杂度低、对传感矩阵列的相关性要求低等特点.

综合上述分析，为了实现矢量水听器阵列在小快拍及低信噪比下的波达方向估计，本文将迂回式匹配追踪算法与波达方向估计相结合，提出一种基于迂回式匹配追踪的波达方向估计算法.在信号重构过程中，采用先逐个最优缩减、后逐个最优扩增的方式更新支撑集，提高重构准确率，扩大重构稀疏信号的稀疏度范围.在存在相干信号、相同的低信噪比、小快拍数条件下的仿真实验中，本方法可以得到有效的DOA估计，与传统超分辨率算法相比更有优势.在真实的水声实验环境中，应用本文方法仍可以得到有效的DOA估计.

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

考虑到水声环境的特殊性，假设声信号为远场窄带信号；水下声速为1 500 m/s；水下传播介质为各向同性的非色散均匀线性介质；仿真实验中使用的噪声为高斯白噪声；阵列的各传感器间不存在互耦情形.

假设t时刻有K个空间信号s1(t),s2(t),…,sK(t)入射到由M个矢量水听器按间距d组成的均匀线阵中，第k个信号的到达角为θk,k=1,2,…,K.以阵列中的第一个矢量水听器为参考，t时刻第m个阵元的输出为

m=1,2,…,M,

(1)

(2)

(3)

式中：k=1,2,…,K，将该阵列的输出写成矩阵形式，即

X(t)=A(θ)S(t)+N(t),

(4)

式中：X(t)=[x1(t),x2(t),…,xM(t)]T，X(t)∈CM×1为阵列输出向量；A(θ)=[a(θ1)⊗u1,a(θ2)⊗u2,…,a(θK)⊗uK]为矢量水听器阵列的信号方向矩阵，也叫流形矩阵；符号⊗为克罗内克积；S(t)=[s1(t),s2(t),…,sK(t)]T，S(t)∈CK×1称为信号源矢量；N(t)=[n1(t),n2(t),…,nM(t)]T，N(t)∈CM×1为阵列接收噪声矢量.

2 基于迂回式匹配追踪算法的DOA估计

2.1 DOA估计模型的稀疏化描述

G(θ)=[a(θ1),a(θ2),…,a(θN)]，

(5)

故DOA估计的数学模型可以转化为

X(t)=G(θ)S(t)+N(t).

(6)

该式为典型的稀疏表示模型，因为S(t)为空域内的稀疏信号，不需要额外的稀疏基来对S(t)稀疏表示，所以扩展后的超完备冗余字典G(θ)既是压缩感知模型中的测量矩阵也是传感矩阵.同时，测量矩阵是通过对空间进行等角度网格划分得到的，由这种稀疏表示方式得到的测量矩阵满足等距约束性质(Restricted isometry property,RIP)，因此可以通过求解式(6)来解决稀疏信号重构问题，从而得到DOA估计值.

然而，式(6)是一个欠定方程，无法直接求解，信号的重构为求解l0范数最小化问题，即通过

s.t.X(t)=G(θ)S(t)+N(t),

(7)

进行求解.

然而，l0范数最小化为NP-hard问题，无法直接求解，具有代表性的解法主要有凸优化算法和贪婪类重构算法.凸优化算法是一种全局优化算法，如基追踪(BP)算法.凸优化算法通常将式(7) 转化为带有不等式约束的优化问题，即最小l1范数优化问题

(8)

式中：γ为误差.一般地，将式(8)转化成

(9)

式中：λ是一个与噪声有关的正则化参数，即求解无约束优化问题.

然而，最小l1范数优化问题是一个广义线性规划问题，需要求解二阶锥规划，以牺牲计算复杂度来实现准确重构，计算复杂度较高，难以适用于大规模信号的快速DOA估计.

2.2 DMP算法

输入：M×N维传感矩阵G(θ)，M维观测信号X(t)，稀疏信号维度N，稀疏度K，阈值ε.

步骤1(初始化)：

步骤2(获取初始值支撑集)：

步骤3：

令b=b+1.判定是否满足b

步骤5(扩增支撑集)：

νt=

Δ(Λ/{ω1,ω2,…,ωb})∪{νb,…,νt}).

选择νt加入支撑集Λ，判定{ωb,…,ωt}∪{νb,…,νt}和{ωb,…,ωt}所包含的元素是否相同，若不相同，则t=t-1，计算子内积逆矩阵P(Λ/{ω1,…,ωb})∪{νb,…,νt}和系数矩阵C(Λ/{ω1,…,ωb})∪{νb,…,νt}，重复步骤5，直至t=0；若相同，则跳转至步骤3.

步骤6：

步骤7(更新支撑集)：

2.3 计算DOA估计值

3 算法复杂度分析

4 仿真实验

4.1 DMP算法、MUSIC算法在小快拍数、低信噪比条件下的估计性能

在该实验中，分别设置单信源与多信源的环境.其中阵元数目为10，阵元间距为0.5 m，DMP算法为单快拍，MUSIC算法快拍数为10，信噪比(Signal-to-noise ratio，SNR)均为-10 dB，单信源实验的入射角为-15°，信号频率为 1 500 Hz.多信源实验的3个入射角分别为-15°，30°，50°，信号频率分别为1 000 Hz，1 500 Hz，2 000 Hz.实验结果如图1 所示.

由图1 可以看出，在快拍数为10，SNR为-10 dB 的环境中，无论是单信源还是多信源实验，MUSIC算法和ESPRIT算法已经无法估计出正确的DOA值，MUSIC算法还产生了大量的伪峰.在图1(a)中标记极值点的位置也可以看出，MUSIC算法在单信源数下搜寻到的谱峰为-5°，ESPRIT算法得出的估计角度为-6.262 25°，与实际的DOA值误差较大.而本文所提出的DMP算法在单快拍条件下，估计信号入射角度为-14.5°，与实际信号入射角误差为0.5°，依旧准确估计出了入射信号方位角.在图1(b)中，MUSIC算法在多信源数下搜寻到的谱峰为-17.5°，8.5°和46°，ESPRIT算法得出的估计角度为-11.861 4°，44.388 8°和65.119 8°.DMP算法在单快拍条件下，估计信号入射角度为-15.5°，30.5°和48.5°，与实际信号入射角的误差介于0.5°～1.5°之间，表现出较好的估计效果.

(a) 单信源实验

(a) 1个信号源

为了说明DMP算法在单快拍、低信噪比环境中的优势，分别将MUSIC算法和ESPRIT算法的快拍数取为10，20，50和100，其他条件与上述实验一样，表1 展示了不同快拍数下MUSIC算法和ESPRIT算法的单信源实验DOA估计结果.

表1 不同快拍数下MUSIC算法和ESPRIT算法的单信源实验DOA估计结果Tab.1 Single source experiment DOA estimation results of MUSIC and ESPRIT under different snapshots

由表1 可知，当快拍数增加到100时，MUSIC和ESPRIT算法均可估计出正确的DOA值，但是在实际应用中，要得到持续不断的高质量采样往往是困难的，DMP算法可以在小快拍情形下得到正确的DOA估计值，有效减小了MEMS矢量水听器阵列在数据采样时的压力，说明DMP算法具有实际应用的价值.

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

设置矢量均匀线阵阵元数为10，快拍数为10，信噪比以5 dB为步长从-10 dB增加到20 dB，每个信噪比下进行100次蒙特卡洛实验.图2 展示了1个信号源、2个信号源和3个信号源条件下，DMP算法、MUSIC算法和ESPRIT算法各自的均方根误差随信噪比变化的曲线.

由图2 可知，随着信噪比增大，3种算法的均方根误差都在逐渐减小，估计精度都在逐渐增大.在SNR处于-10 dB～-5 dB区间时，DMP算法的估计性能与MUSIC和ESPRIT算法接近，但考虑到DMP算法只取了单快拍下的数据进行DOA估计，因此在信噪比较低时DMP算法只需要一次采样便可达到多快拍下MUSIC和ESPRIT算法的估计性能.-5 dB以上时DMP算法的性能明显优于MUSIC和ESPRIT算法.从整体来看，DMP算法的均方根误差都处于比较低的位置，估计精度较高.总之，DMP算法具有良好的抗噪能力，相同信噪比下只需单快拍数据即可取得相对较高的估计精度，在信噪比较高的环境下具有更显著的优势.

4.3 相干信号源条件下的DOA估计性能对比

设置矢量均匀线阵阵元数为10，间距为0.5 m，3个声源信号的入射角分别为-50°，10°，60°，其中10°与60°方向信号为相干信号，频率分别为1 000 Hz，1 500 Hz和1 500 Hz.为了观察相干信号对MUSIC算法的影响，将MUSIC算法的快拍数设置为100，DMP算法为单快拍，信噪比为0 dB，实验结果见图3.

图3 相干信号的DOA估计实验Fig.3 DOA estimation experiment of coherent signal

由图3 可以看出，当存在相干信源时，MUSIC算法出现了一些伪峰，并且其谱峰较为平滑，几乎无法分辨信源角度.ESPRIT算法则无法得出正确的DOA估计值.DMP算法在单快拍、低信噪比情形下仍能得到有效的DOA估计值，这是由于DMP算法仅仅利用信号源的稀疏性来求解DOA估计值，不受入射信号频率的影响，所以相干信号不会对DOA估计的结果产生较大影响.

5 结 论

针对单快拍、信噪比较低条件下的波达方向估计问题，本文提出一种基于迂回式匹配追踪算法的DOA估计算法.该算法的优势和特点为：

1) 利用压缩感知理论求解DOA问题，通过单快拍即可求解出声源信号的入射角度，尤其是在低信噪比情形下，单快拍即可达到小快拍数下传统的子空间分解类方法的估计精度.

2) 仿真实验表明，基于迂回式匹配追踪算法的DOA估计算法求解的过程不涉及矩阵的特征分解，克服了相干信号对估计精度造成的影响，在处理小块拍、低信噪比的相干信号时仍然有效.