一种改进广义正交匹配追踪的DOA 估计方法

刁弘扬，胡洲勇，禹永植

1.哈尔滨工程大学 信息与通信工程学院，黑龙江 哈尔滨 150001

2.北京遥感设备研究所，北京 100854

信号的波达方向(direction of arrival,DOA)估计在阵列信号处理领域占据重要地位，广泛应用于雷达、医学诊疗等领域，近年来学者们进行了大量的研究，取得了非常显著的成果[1−3]。以MUSIC[4]和ESPRIT[5]算法为代表的经典子空间分解类算法，实现简单且分辨率高，并且在此基础上发展出了多种改进算法[6−8]，但是这类算法运算量大，需要在信噪比高且快拍数多的情况下才能得到较好的估计性能。压缩感知[9−10]理论是一个新兴起的领域，它能够突破奈奎斯特采样定理的限制，压缩感知理论中的很多重构算法都被用于DOA 估计中，其中平滑l0范数重构算法[11−12]、

l1−SVD算法[13−14]和贪婪算法是广泛应用的方法。平滑l0范数重构算法采用最速下降法[15]和梯度投影的原理，构造一个函数来近似l0范数，但是此算法不能保证重构的精确度。l1−SVD算法利用l1范数重构算法对信号进行重构，估计出信号的DOA，能够在低快拍下实现DOA 估计；但是该算法重构过程复杂，会导致较高的计算复杂度，并且在求解过程中，需要人工去设置所需参数，参数选取不好将会对最终的结果产生较大的影响。常用的贪婪算法有匹配追踪(matching pursuit,MP)、正交匹配追踪(orthogonal matching pursuit,OMP)[16]、正则OMP[17]等，该类算法需要找到测量矩阵中能够提供最大能量的某列原子，然后一直进行迭代，直到满足迭代条件，得到最优解。OMP 算法是一种常用的重构算法，研究者们提出了多种改进算法[18−19]，但是OMP 算法运算时间长，DOA 估计精度有待进一步提高。

针对稀疏阵列中贪婪算法DOA 估计求解问题，本文引入广义正交匹配追踪(generalized orthogonal matching pursuit,GOMP）用于DOA 估计，仿真表明直接将GOMP 算法用于DOA 估计较OMP算法性能有所提高；为了进一步提高DOA 估计的精度，引入最优化理论中的最速下降法，提出基于最速下降法的改进GOMP 算法。

1 基于稀疏表示的DOA 估计模型

假设有一个由M个阵元组成的均匀线阵，阵元间距为d，有K个方向分别为θi(i=1,2,···,K)的远场窄带信号xk(t),k∈{1,2,···,K}。通过网格划分将空间信号稀疏化，整个空间范围等角度均匀划分为N份，则所有可能的方向可以表示为θi(i=1,2,···,N)。将空间信号进行稀疏化，只有K个方向上存在真实的目标信号源，其中，N远远大于实际信号源的数量K。则空间中所有的方向集合都与一个导向矢量相对应，导向矢量为A=[a(θ1),a(θ2),···,a(θN)]，导向矢量也称为阵列流型矩阵，则t时刻阵列接收的信号为

式中：x(t)=[x1(t),x2(t),···,xN(t)]T为空间入射信号矢量，其稀疏度为K(K≪N)，x(t)中只有K个非零元素，其余N−K个为零元素，非零元素所在的位置即为真实信号源的目标方向；n(t)=[n1(t),n2(t),···,nM(t)]T是均值为0，方差为 σ2的高斯白噪声；y(t)=[y1(t),y2(t),···yM(t)]T为阵列接收信号矢量。

稀疏表示的DOA 估计就是在已知阵列接收矢量y以及阵列流型矩阵A的条件下，找出x(t)中所有非零元素所在的位置集合，该位置集合即是θi(i=1,2,···,N)中实际信号源所在的角度。

2 基于改进GOMP 算法的DOA 估计

2.1 OMP 算法

由式(1)、(2)可知，x在空间上是K阶稀疏的，只有K个非零的元素，而且非零元素对应的位置即为θi(i=1,2,···,N)中实际信号源的角度。

为了找出x中非零元素所在的位置，需要确定阵列流型矩阵A中的有效矢量。OMP 算法通过贪婪搜索的方式来搜索有效列。在每次迭代过程中，找出A中 与阵列接收信号y中剩余部分最相关的一列，剔除该列对于y的贡献，标记该列并更新残差。当完成K次迭代后，即可找出参与对y的线性组合中列的标记集合，最后用该标记集合，在集合θi(i=1,2,···,N)中找到信号源所在的真实角度位置。

2.2 GOMP 算法

基于OMP 算法的DOA 估计方法，耗时比较长且估计的精度不高，为了提高DOA 估计精度，采用广义正交匹配追踪算法来求解DOA 估计。

OMP 算法和GOMP 算法均是通过贪婪迭代从阵列流型矩阵A中找到能够提供最大能量的某些列原子，代价函数选用y残差最小。2 种算法的区别在于进行重构时，OMP 算法每次迭代过程中只选择与y中剩余部分最相关的一列，而GOMP算法在每次迭代过程中会选择与y中剩余部分最相关的K列，然后在所有选出的可能的坐标集合中进行搜寻，在索引坐标集合中加入使估计能量值最大的坐标。相对于OMP 算法，GOMP 算法的搜索选择方式能够更好地保证从信号能量相近或者不完全和不精确的阵列接收信号中恢复信号，能更准确地找出信号源所在的真实角度位置。

2.3 最速下降法的相关理论

最速下降法基于最优化理论[20]，从梯度的角度来求解最优化问题。函数f(x)的泰勒展开式为

式中：o(x)是高阶无穷小量，其在xk的梯度为gk=∇f(xk)。记x−xk=akdk，其中，ak是步长，是一个大于零的数值，dk为梯度的方向，则泰勒展开式为f(x)=f(xk)++o(x)。在＜0时f(x)＜f(xk)，此时的函数值是下降的，dk则是下降的方向。

由柯西−施瓦兹不等式可知：

当且仅当dk=gk时，等号成立，最大。所以当且仅当dk=−gk时，最小，以此为搜索方向，此时f(x)是下降最快的。

最速下降法在每次迭代时选择一个合适的步长，使得目标函数每次减小的程度最大，即每次得到的值和目标函数值之间的差距最小。用f(x)代表目标函数，迭代步长为

梯度的迭代为

因此计算步长的公式就是为找到最好的点xk+1，函数f(x)在这点可以取得最小值，此时f(xk−agk)的极小值就是此次迭代的步长。

最速下降法的过程可以表示为在每次迭代中，从xk出发，按照梯度的负方向去寻找最好的步长值，根据寻找到的步长值来确定下一次迭代过程的出发点xk+1。

将函数的类型扩展至二次型，设函数为

式中：Q是 对称正定阵；b∈Rn；x∈Rn；令

迭代表示为

步长为

在gk=0时迭代过程停止，在不为零时，设φk(a)=f(xk−agk)，由局部最小点的一阶必要条件得到：

2.4 改进GOMP 算法

将最速下降法和GOMP 相结合得到改进GOMP算法，改进算法在信号重构时不使用计算量大的最小二乘法，而是沿着负梯度方向进行搜索而得到最优解。

改进算法中迭代选择原子及更新索引集的过程和GOMP 算法相一致，在恢复信号过程中，GOMP等贪婪算法均使用最小二乘法进行重构，其矩阵公式为x=(ATA)−1ATy，由于需要进行矩阵的求逆运算，导致运算量很大，会消耗很长的时间；而最速下降法虽然需要进行多次运算，但都是基本运算，相对而言计算量较小，所要耗费的时间也会更少。因此改进算法使用最速下降法代替最小二乘法来恢复信号。

用Φ={φk}表示测量矩阵，此测量矩阵即为DOA 估计模型中的阵列流型矩阵，y表示要恢复的信号，xi表示第i次循环得到的近似信号，ri表示残差，Λi表示索引集合，Ψ 表示选择原子的列向量集合，先得到矩阵的负梯度方向为gi=ΨiTri，然后计算每次迭代的步长为

式中ci=Ψigi。从而得到恢复信号：

最后更新残差：

完成循环，判断是否满足迭代停止条件。

改进后的GOMP 算法步骤可总结如下：

1)初始化残差r0=y，初始化索引集 Λ为空集，初始化迭代次数i=1；

2)计算残差和测量矩阵 Φ的每一列的内积，找出若干个绝对值最大的元素，并将这些值对应的列序号构成集合P；

3)更新索引集Λi=Λi−1∪P，更新列向量集合Ψi=Ψi−1∪φP；

5)利用式(4)更新残差；

6)若满足迭代停止条件，根据索引集 Λ在角度θi(i=1,2,···,N)中找到信号源所在的真实角度位置；若不满足，则跳转到步骤2)进行下一次循环。

3 仿真结果分析

本文通过MATLAB 仿真验证所提算法的可行性，并与l1−SVD算法、OMP 算法、GOMP 算法进行实验对比验证其性能。在实验仿真中，采用均匀直线阵列，噪声为高斯白噪声，均匀线阵的阵元间距为半波长，按照等角度划分方式将空间信号等角度划分为181 份，每份角度为1°。均方根误差表示为

式中：M表示蒙特卡洛实验次数；K表示信号源的总个数；θkm表示第k个空间信号源的DOA 真实值；表示第k个空间信号源的DOA 估计值。

实验1残差二范数对比实验。设2 个信号源以入射角10°、30°入射到均匀直线阵列上，对结合最速下降法的改进GOMP 算法和原始GOMP算法进行性能对比，在单快拍条件下，比较2 种情况下的残差二范数，进行500 次蒙特卡洛实验，仿真结果如图1 所示。

图1 改进GOMP 算法与GOMP 算法残差比较

通过图1 可以看出，改进GOMP 算法比GOMP算法的残差值小，随着信噪比的提升，改进算法的残差值下降更快，可见对GOMP 算法的优化是有效的。

实验2不同信噪比下DOA 估计的性能对比实验。设2 个信号源以入射角10°、30°入射到均匀直线阵列上，均匀线阵阵元数为12，快拍数为8，信噪比以2 dB 为步进从−8～8 dB 变化，在相同的实验条件下进行500 次蒙特卡洛实验。当估计偏差小于或等于1°时，认为估计成功，在不同信噪比条件下，统计各个条件下的成功概率和均方根误差，仿真结果如图2、3 所示。

图2 不同信噪比下DOA 估计成功概率

图3 不同信噪比下DOA 估计均方根误差

在不同信噪比条件下，由图2 可以看出，与其他3 种算法相比，改进GOMP 算法具有更高的估计成功概率，在信噪比为−2 dB 时估计成功概率接近100%。由图3 可以看出，在高信噪比情况下，4 种算法均趋于稳定，改进GOMP 算法具有更小的均方根误差，特别是在低信噪比情况下，改进GOMP 算法具有明显的优势。

实验3不同快拍数下的DOA 估计的性能对比实验。设2 个信号源以入射角10°、30°入射到均匀直线阵列上，均匀线阵阵元数为12，信噪比为−4 dB，在相同的实验条件下进行500 次蒙特卡洛实验。当估计偏差小于或等于1°时，认为估计成功，在快拍数不同的情况下，统计各个条件下的成功概率和均方根误差，仿真结果如图4、5 所示。

图4 不同快拍数下DOA 估计成功率

图5 不同快拍数下DOA 估计均方根误差

在不同快拍数情况下，由图4 可以看出，改进GOMP 算法具有更高的估计成功概率，且具有很高的稳定性。在快拍数为12 时，已经达到以100%的概率估计出信号的DOA。随着快拍数的增加，其他3 种算法的成功率缓慢增加至趋于平稳。由图5 可以看出，4 种算法的均方根误差均随着快拍数的增加逐渐减小，改进GOMP 算法一直保持着较小的均方根误差，特别是在低快拍数情况下，改进GOMP 算法具有明显的优势。

实验4不同快拍数情况下算法运行时间对比实验。由于l1-SVD 算法运行时间相对于贪婪类算法运行时间长，所以只比较另外3 种算法的运行时间。由图6 可以看出，随着快拍数的增加，3 种算法的运行时间均逐步上升，且改进GOMP算法相较于其他2 种算法的运行时间相对较短；在相同的条件下，同其他2 种算法相比，改进GOMP 算法具有更短的运算时间。因此从整体上看改进GOMP 算法的运行时间相对偏低，相比于其他2 种算法具有运行时间上的优势。

图6 不同快拍数下算法运行时间比较

4 结论

本文提出了一种新的贪婪DOA 估计方法，该方法可有效提高DOA 估计的精度。GOMP 是一种改进的贪婪恢复算法，其利用更好的原子选择方式，提高了贪婪算法的稳定性和角度估计精度。

1)利用GOMP 算法求解DOA 问题，相较于原有贪婪DOA 估计方法提高了估计精度和成功率。

2)引入最速下降法对GOMP 算法进行优化，并通过仿真验证了改进方法具有更高的估计精度和成功率，降低了算法的运行时间，相对于原有贪婪DOA 估计方法的性能有一定的提升。