基于改进稀疏度自适应匹配追踪算法的压缩感知DOA估计

2021-11-17窦慧晶肖子恒

北京工业大学学报 2021年11期

窦慧晶,肖子恒,杨帆

(北京工业大学信息学部，北京 100124)

作为阵列信号处理领域的关键问题，波达方向(direction of arrival,DOA)估计在雷达、通信、地震等众多领域得到广泛应用[1-2].目前，国内外学者提出的DOA估计算法主要有两大类：第一类是以Capon算法为代表的波束形成算法[3]，通过抑制无用信号、提升有用信号来获得信号的空域谱；另一类是以多重信号分类(multiple signal classification,MUSIC)算法[4]以及旋转不变子空间(estimation of signal parameters via rotational invariance technique,ESPRIT)算法[5]为代表的子空间类算法，它对接收数据的协方差矩阵进行特征值分解得到信号子空间和噪声子空间，再根据子空间的正交性构造出空间谱.然而，基于波束形成的算法虽克服了阵列孔径限制，但估计分辨率不高且抗噪能力弱；子空间类算法需要大量独立同分布的采样数据使阵列协方差矩阵满足统计特性，因此计算量较大，难以满足现代信号处理的实时性要求.

近年来Donoho等[6]提出压缩感知理论受到广泛关注，它利用信号的稀疏性，在远小于奈奎斯特频率的采样率下，通过求解l0最小范数优化问题得到高精度的重构信号.在实际空域中目标信号相比于整个空间角度是极少的，本身具有稀疏性，因此非常适合在压缩感知模型下完成DOA估计的求解.2003年，文献[7]首次提出l1-奇异值分解(singular value decomposition，SVD)算法进行DOA估计，该算法利用奇异值分解和l1范数的凸优化，求解优化函数得到稀疏向量.l1-SVD算法可直接用于相干信号的DOA估计，且估计精度高于其他传统算法.但由于使用子空间类算法的SVD，不精确的信源数会导致估计性能的下降或失效.为解决这一问题，文献[8]提出l1-阵列协方差稀疏表示(sparse representation of array covariance vectors，SRACV)算法，它利用矢量化样本协方差矩阵的统计特性，仿照l1-SVD算法构造基追踪问题并求解得出DOA估计.该算法在具有无须已知信源数以及可以处理相干信号的优点的同时，还保证低信噪比(signal noise ratio，SNR)下的鲁棒性，但是需要大量运算来计算协方差矩阵及其特征值分解，计算复杂度远大于l1-SVD算法以及子空间算法.文献[9]提出一种基于单快拍数的动态参数l1正则化算法，该算法能有效解决真实方位角与离散的空间网格的不一致问题，同时DOA估计的准确度也比较高.

l1范数凸优化类算法重构精度虽高，但由于涉及到二阶锥的求解使得算法运行时间较长，难以应用于大规模信号，实用性较差.因此，基于压缩感知贪婪算法的DOA估计近年来受到研究学者们的关注.文献[10]提出单快拍数下压缩感知DOA估计的算法，并将正交匹配追踪(orthogonal matching pursuit，OMP)算法[11]应用于DOA估计.OMP算法是基于(matching pursuit，MP)算法改进的，它在筛选测量矩阵中的列向量原子时，不仅保证所选取的列向量与信号相关系数值最大，而且还使得与残差矢量正交，从而能够得到最优解.文献[12]在OMP算法的基础上提出一种结合局部优化处理的正交匹配追踪DOA估计方法，只需要少量计算复杂度的代价就能提高DOA的估计精度.文献[13]提出一种使用正则化方法筛选原子集的估计算法，达到减少迭代次数和提升重构精度的目的.但是以上贪婪匹配追踪类算法都需要以信源数已知为先验条件，而在现实情况里往往无法得知确切的信源数.针对这一问题，Do等[14]首次提出在每次迭代中增加固定步长以逐步逼近信号真实稀疏度值的稀疏度自适应匹配追踪(sparsity adaptive matching pursuit，SAMP)算法，然而该算法在噪声环境下无法收敛，从而引入大量错误原子导致估计算法失效，并且由于使用固定步长，灵活度不高，选择小的初始步长会大大增加迭代次数，导致运行时间变长，选择大的初始步长容易造成稀疏度过估计，影响估计精度.

针对上述问题，本文将稀疏度自适应匹配追踪算法与DOA估计相结合，在此基础上根据噪声环境下迭代残差的变化规律设计一种合适的迭代终止条件，并通过迭代状态控制步长自适应变化，提出一种改进稀疏度自适应匹配追踪的DOA估计算法.本文算法能够在单快拍、不预设信源数的条件下有效估计相干信号，相比改进前的算法能提升噪声环境下的DOA估计性能，且通过选择大初始步长，既降低运行时间，又保证良好的估计精度.

1 阵列信号模型

假设在t时刻有K个远场窄带信号s1(t),s2(t),…,sK(t)入射到阵元数为M、阵元间距为d的均匀线阵中，其来波方向分别为θ1,θ2,…,θK，以第1个阵元为参考阵，则信号入射到天线的示意图如图1所示.

图1 接受天线阵列图Fig.1 Receiving antenna array

因为K个入射信号都是远场信源，所以各信号被阵元接收到的角度方向都相等.则各阵元的输出信号xm(t)可表示为所有入射信号延时相加的形式

(1)

式中：nm(t)为第m个阵元上的加性高斯白噪声；τm,k为第m个阵元接收到信号sk(t)时相对参考阵元的延迟时间，在均匀线阵的条件下可表示为τm,k=2πd(m-1)sinθk/λ，其中λ为信号波长.因此，sk(t-τm,k)可以表示为sk(t)e-j2πd(m-1)sin θk/λ.

将sk(t)到达各阵元的相位差组成向量形式记为

α(θk)=[1,e-j2πdsin θk/λ,…,e-j2πd(M-1)sin θk/λ]T=
[α1(θk),α2(θk),…,αM(θk)]T

(2)

则xm(t)可以写成

(3)

将K个信号的导向矢量构成矩阵A(θ)=[a(θ1),a(θ2),…,a(θK)]，则接收信号模型(3)可描述为

x(t)=A(θ)s(t)+n(t)

(4)

式中：x(t)=[x1(t),x2(t),…,xM(t)]T为阵列输出向量；s(t)=[s1(t),s2(t)，…，sK(t)]T为信源向量；n(t)=[n1(t),n2(t),…,nM(t)]T为加性噪声向量；A(θ)为M×K阶阵列流型矩阵，矩阵第m行k列的元素表示第m个天线阵元接收到第k个远场信号的增益和相位延迟信息.

2 算法原理

2.1 压缩感知理论基础

压缩感知理论表明，当信号具有稀疏性或可压缩性时，通过少量投影数据就可实现信号的精确重构.设复信号z=[z1,z2,…,zn]T可以用某稀疏基矩阵ψ=[ψ1,ψ2,…,ψP]的线性组合表示，即z=ψs，其中s是一个p×1维的具有k(k≪p)稀疏的向量.存在一个与稀疏基ψ不相关的m×n维的测量矩阵φ，z在投影测量矩阵φ下的线性测量可表示为

y=φ·z=φψ·s=Θ·s

(5)

式中Θ=φψ为m×p维矩阵，表示推广后的测量矩阵，称为感知矩阵.则y可看成是信号s在感知矩阵Θ下的线性投影.如果Θ满足有限等距性[15]等稀疏重构条件，就可以求解以下l0范数问题以极高的概率重构稀疏信号s：

min‖s‖0s.t.y=Θ·s

(6)

2.2 基于压缩感知的DOA估计模型

(7)

式中：

2.3 基于改进稀疏度自适应匹配追踪的DOA估计算法

(8)

式中：‖h‖0为h的0范数，表示h中非零元素的数量；ε为与噪声有关的常量，式(8)在本质上是具有非凸稀疏约束的优化问题，并且已经被证明是计算复杂度很高的NP-hard问题[17].目前常见重构算法分为凸优化算法和贪婪算法.凸优化算法是在一定条件下将l0范数的非凸优化问题转为l1范数的凸优化问题进行求解，最具代表性的方法是基追踪算法.此类算法虽具有较高的重构精度，但计算复杂度高，难以满足现代信号处理的实时性要求，实用性不强[18].

但是前述的所有算法都需要以信号真实的稀疏度K为先验条件，在算法中以K来控制迭代次数，这与现实情况相矛盾.所以在稀疏度未知情况下，SAMP算法可以克服前述的贪婪算法缺陷实现信号精确重构.SAMP算法通过分阶段的思想实现稀疏度自适应，其具体过程如下：首先设置一个初始步长s，首个阶段在稀疏度I=s下进行循环迭代重构信号.如果稀疏度是合适的，则保持支撑集大小为该稀疏度进行迭代直到满足迭代终止条件.否则迭代后的残差会增大，表明稀疏度不够，需要增大支撑集，令稀疏度I=I+s继续迭代.如此进行下去，直到满足迭代停止条件输出结果.然而如何选择初始步长是个关键问题，它的取值对信号的重建效果有很大影响.一般s≤K，如果s取1，则一定能准确达到信号的真实稀疏度，但会增加算法的迭代次数，从而导致运行时间变长.如果s取较大值，可以减少重构时间，但容易出现估计的稀疏度超过真实稀疏度K的过估计现象，从而影响重建效果.因此SAMP算法无法兼顾重构速度和重构精度.

鉴于这一情况，本文选择大步长作为初始步长快速逼近，当支撑集还未达到K时，迭代的残差是不满足算法终止条件的，当迭代过程中第1次达到迭代终止条件时，说明出现精确估计或过估计这两者情况中的一种.对此，可以将迭代状态回溯至上一次，再使用小步长(s=1)迭代来逼近K，以解决SAMP算法中速度与精度的矛盾.

在SAMP算法中通常以残差的l2范数来判定迭代是否终止，即当‖r‖2≤ε时终止迭代，其中ε取值为较小的固定常数.对于信噪比较大的信号，这种方法对信号的重构精度影响不大.然而，如果原始信号中夹杂大量噪声，则噪声干扰会导致迭代无法达到该终止条件并超过真实的稀疏度，在后续迭代中还会加入多余的错误原子，干扰信号的精确重构.算法在每次迭代过程中，都会优先选择与观测信号最相关的原子加入支撑集，从而使得残差逐渐减小.当稀疏信号中所有非零元素对应的原子都加入后，剩下的都是和观测信号相关度低的原子.如果后面继续保持迭代过程，只能选出低相关度的原子，残差的能量变化很小.因此可以利用这一特点来判断迭代是否达到准确的稀疏度，改进后的迭代终止条件可表示为

(9)

式中：‖rt‖2为迭代t次的残差值；It为第t次迭代的稀疏度.式(9)左侧表示的意义是相邻2次迭代残差变化的比例，分子分母都除以迭代的稀疏度之差是为了在步长变化后也能进行终止条件的判断.如果满足式(9)，则说明在第t-1次迭代后残差下降减缓，出现“拐点”，那么It-1即为真实的信源稀疏度，直接输出第t-1次的估计结果.在大部分情况下这种判断方式都是正确有效的，但如果在信噪比极低的情况下，噪声的能量远大于信号的能量.即使迭代次数达到真实稀疏度，后续迭代还会选出与噪声匹配的原子，使残差继续减小，“拐点”不明显，从而导致式(9)的判定条件失效.这时候可以选择迭代残差首次小于噪声的能量作为迭代终止的判定依据，用公式可表示为

‖rt‖2

(10)

式中εn为噪声能量，但实际迭代过程中会因为噪声干扰使得算法的重构误差变大，因此还要考虑加入误差平衡因子c，使其能够迭代到正确的稀疏度再输出结果.通过大量实验数据测得，c取0.7最为合适.通过式(10)使得算法在噪声干扰大的情况下也能保证一定的估计效果.

改进的稀疏度自适应匹配追踪算法(improved SAMP,ISAMP)步骤如下：

输出N×1维稀疏信号，稀疏度I.

步骤1初始化.残差r0=X，支撑集Λ0=∅，稀疏度I0=0，迭代次数t=1，稀疏度I1=s.

步骤6判断是否初次满足式(9)(10)的迭代终止条件，若满足则t=t-1，It=It-s，s=1转入步骤8；否则转入步骤7.

步骤7判断是否满足式(9)(10)的迭代终止条件，若满足则停止迭代并转入步骤9；否则转入步骤8.

步骤8t=t+1,It=It-1+s转入步骤2.

如果信号的稀疏度不是初始步长的整数倍，那么SAMP算法无论迭代多少次也无法准确估计稀疏度，但ISAMP算法通过步骤6将迭代状态回溯至上一次并采用小步长迭代，避免了因为大步长迭代造成的过估计.同时利用式(9)(10)作为迭代终止条件摆脱了对K的依赖，保证算法在低信噪比条件下依然能够有效估计.

3 仿真结果分析

仿真条件说明：本节实验采用的阵列都为均匀线阵，阵元间距为半波长，在网格划分时角度空间[-90°,90°]都均分为361个角度，即角度网格为0.5°.实验中评价估计性能的指标为均方根误差和估计成功率，均方根误差的定义为

(11)

式中：K为信源数；N设为100是蒙特卡罗实验次数；kn为第k个信源第n次试验的角度估计值；θk为第k个信源真实的波达角.估计成功率定义为

(12)

式中sk表示在同一组空间信源的N次DOA估计实验中第k个信源估计成功的次数.信源估计成功的判定条件是角度偏差不大于1°，实验中蒙特卡洛实验次数N为100，即重复100次相同条件下的DOA估计，统计估计成功的实验次数，并以此计算估计成功率.(处理器为Intel(R)Core(TM)i5-8400 CPU@2.80 GHz；软件平台为Matlab 2008b)

实验一：仿真比较信源相干条件下各算法的DOA估计性能.

仿真设置阵元数M=21，信噪比为20 dB，3个远场窄带信号的入射角度分别为-10°、10°和25°，在相干信源和非相干信源的条件下，比较单快拍下SAMP算法、ISAMP算法以及快拍数为100的MUSIC算法的DOA估计结果，实验结果见图2、3.

图2 各算法对非相干信号源的DOA估计Fig.2 DOA estimation of incoherent signal sources by various algorithms

图3 各算法对相干信号源的DOA估计Fig.3 DOA estimation of coherent signal sources by various algorithms

由图2、3可知，在估计非相干信号源时，3种算法都准确地估计出波达方向角度.对于相干信源的估计，MUSIC算法中出现多数伪峰，且谱峰不够尖锐，几乎无法分辨信源角度.而SAMP算法和ISAMP算法仍然准确地估计出3个信源角度，实现对于相干信号源的DOA估计.这是因为基于压缩感知的DOA估计不用对矩阵进行特征分解，因此具有解相干的天然优势，解决了子空间类算法无法估计相干信源的问题.并且在实验中贪婪类算法只用单快拍就完成DOA的精确估计，证明压缩感知DOA估计算法对采样信息量要求低的特点，减少数据存储及运算的成本.

实验二：仿真比较不同信噪比环境下残差能量随迭代次数变化的趋势.

仿真设置阵元数为21，3个远场窄带信号的入射角度分别为-10°、10°和25°.实验对比信噪比分别为10、0、-10 dB时单快拍ISAMP算法的残差大小随迭代次数的变化曲线，结果见图4.

图4 不同信噪比下残差能量随迭代次数变化曲线Fig.4 Variation curve of residual energy with iteration times under different SNR

在图4中可以看到，10、0 dB条件下的ISAMP算法在迭代次数没到真实稀疏度3时，每次迭代残差能量都会快速减小；当迭代次数达到3之后，残差在后续的迭代过程中几乎不变.所以真实的信源稀疏度就是图中拐点出现时的迭代次数.但在信噪比为-10 dB时，即使达到真实稀疏度后残差依然保持着一定的速度减小，上述拐点并不明显.因此需要依靠式(10)作为迭代终止条件.

实验三：比较SAMP算法、ISAMP算法、MUSIC算法、ESPRIT算法在噪声环境下的估计性能.

1)仿真设置阵元数为21，3个远场窄带信号的入射角度分别为-10°、10°和25°.仿真中加入信噪比为0 dB的高斯白噪声干扰，比较单快拍下SAMP算法和ISAMP算法在低信噪比下的估计结果.

从图5可以看出，SAMP算法在0 dB的信噪比条件下已经无法正确估计出DOA，产生大量虚假信号，这是因为信号中混杂的噪声使迭代无法达到预设的阈值，并在后续迭代过程中不断匹配错误的原子，直到达到最大迭代次数，输出与真实信号无关的结果，算法已经失效，因此在后续的一系列实验中将不再比较SAMP算法.从图6可以看出，本文所提出的ISAMP算法准确估计出3个方位角，证明本文算法在噪声环境下依然能够正确估计波达方向角度.

图5 信噪比为0 dB时SAMP算法的DOA估计Fig.5 DOA estimation of SAMP algorithm when SNR is 0 dB

图6 信噪比为0 dB时ISAMP算法的DOA估计Fig.6 DOA estimation of ISAMP algorithm when SNR is 0 dB

2)为进一步考察ISAMP算法在不同信噪比条件下的估计精度，本次实验采用均方根误差作为衡量DOA估计算法性能好坏的指标.仿真中设置阵元数为12，快拍数为15，加入高斯白噪声干扰，信噪比从-10 dB以2 dB为间隔提升到10 dB，比较ISAMP算法、MUSIC算法、ESPRIT算法的均方根误差随信噪比变化的曲线，实验结果如图7所示.

图7 信噪比变化时不同算法的均方根误差比较Fig.7 Comparison of RMSE of different algorithms with SNR

由图7可知，随着信噪比增大，3种算法的均方根误差逐渐减小，估计精度变高.在信噪比处于-10～0 dB区间时ISAMP算法的均方根误差低于另外2种算法，0 dB以上时3种算法的均方根误差几乎相同.从整体来看，本文算法的均方根误差都处于较低的位置，估计精度较高.改进后的算法具有良好的抗噪能力，同样信噪比条件下估计精度更高，且在低信噪比的情况下具有更明显的优势.

实验四：比较ISAMP算法在不同初始步长下的估计性能.

由于实验中ISAMP算法在某些特殊的阵元下进行DOA估计时可能出现对信源数目的误判，或是个别错误估计的角度与真实角度有很大偏差的情况，导致算法的均方估计误差波动较大，无法真实体现算法性能，因此，本实验使用成功估计概率作为评价指标.仿真设置阵元数为12，信噪比为-5 dB，快拍数为20.图8是不同初始步长的ISAMP算法、MUSIC算法、ESPRIT算法的估计成功率随阵元数变化的实验结果，图9是阵元数为40时各算法估计时间随信号稀疏度变化的实验结果.

图8 不同步长下估计成功率比较Fig.8 Comparison of estimated success rate under different step sizes

图9 不同步长下运行时间比较Fig.9 Comparison of running time under different step sizes

从图8可以看出，随着阵元数M的增加各算法的估计成功率逐渐提高,在阵元数为6～15时，ISAMP算法的估计成功率高于另外2种算法.并且ISAMP算法在初始步长为1、3时的估计成功率差别很小，这是因为当选择大步长进行迭代时，ISAMP算法使用了过估回调的方法控制步长变化，本质上还是小步长的迭代，所以同样能够准确估计DOA.

从图9可以看出，MUSIC算法和ESPRIT算法在不同稀疏度条件下的运行时间变化不大.而ISAMP算法随着稀疏度的增加估计时间也越长，且大步长的运行时间明显低于小步长的.这是因为稀疏度越大，ISAMP算法需要迭代的次数就越多，而大步长能够以更少的迭代次数达到迭代终止条件，提高估计速度.从本实验可以看出，ISAMP算法在选择大的初始步长时具有明显优势，可以兼顾估计精度以及运行时间.