采用信号子空间稀疏表示的DOA估计方法

解 虎，冯大政，魏倩茹

（西安电子科技大学雷达信号处理国家重点实验室，陕西 西 安710071）

0 引 言

波达方向（direction of arrival，DOA）估计是阵列信号处理的重要研究内容之一，被广泛应用于雷达、无线通信、电磁场、声纳、地震勘探等诸多领域。DOA估计方法的主要目标是在噪声环境下，分辨两个在方位向非常接近的目标，为此，利用信号源个数小于阵元数的特点，人们提出了大量高分辨DOA估计方法。例如，经典的多重信号分类（multiple signal classification，MUSIC）法，其关键是信号子空间维数的确定。本文也采用类似的方法，首先确定信号子空间的维数，然后将DOA估计问题构造成一个稀疏表示问题以获得较高的角度分辨率。

基于稀疏表示的DOA估计方法与已有的方法虽有相同之处，但是在本质上却有不同。常用的DOA估计方法分为两类，即参数化估计和非参数化估计。对于非参数估计方法，主要有波束形成法，基于子空间方法的MUSIC和基于最小方差无畸变（minimum variance distortionless response，MVDR）的高分辨谱估计法等。其中波束扫描方法虽然与信噪比（signal-to-noise ratio，SNR）无关，但是受到瑞利限的制约，无法分辨一个波束宽度内的两个信号。而MUSIC和MVDR方法虽然能够分辨一个瑞利单元内的两个目标，但是却需要较高的SNR和大量的样本，而且对信号源之间的相干性有一定要求，不能对相干信号源直接进行有效分辨。基于最大似然的参数化估计方法分为：确定最大似然和统计最大似然［1］。这类方法虽然有很好的统计特性，但是需要准确的初始值才能收敛到全局最优点。而对于基于稀疏表示的DOA估计方法，不需要一个好的初始值和大量样本，对信号相关性也没有特别限定，依然能获得高分辨率。

近几十年来，稀疏表示取得了重大进展，已经被广泛应用于图像重建与理解［2-3］、模式识别［4］和高分辨谱估计［5-6］等诸多领域。Gorodnitsky等提出了一种加权迭代最小2-范数（focal underdetermined system solver，FOCUSS）［7］方法来求解稀疏表示问题，并在DOA估计上取得较好的结果，但是仅适用于单次快拍。然而单快拍信号容易受到信号幅度变化的影响，而且在低信噪比的情况下，正确估计和恢复来波信号的概率较小。实际中，用于DOA估计的样本个数比较充足，因此人们提出了多样本下的基于稀疏表示的DOA估计方法，以达到高分辨率。文献［8］提出了一种基于奇异值分解的多测量矢量欠定系统正则化聚焦求解算法实现DOA超分辨估计。该算法计算量快，但是对正则化参数比较敏感。文献［9］采用l1，2混合范数对多快拍接收信号的协方差矩阵进行联合稀疏表示，达到DOA估计的目的。由于需要对协方差矩阵求逆，当样本较少时，会对导致算法性能下降，但是该算法不需要确定参数，适用范围更广。文献［10］中，作者采用加权l1范数确保解的稀疏性，并将其应用于DOA估计。然而其需要Capon估计加权系数，因此当存在相关信号时，算法性能会下降。文献［11］提出一种将基于协方差的多快拍信号联合稀疏表示转化成单一向量的稀疏表示模型，该方法计算量小，但是要求信号之间互不相关。文献［12］进一步将稀疏表示方法应用于宽带DOA估计中，在低信噪比下取得了良好的估计精度。目前基于多快拍信号稀疏表示的模型主要分为两类：对信号子空间进行稀疏表示和对信号协方差矩阵进行稀疏表示。同等条件下（样本数目相同，参数选择合理），基于信号子空间的算法性能一般要优于基于协方差矩阵的算法，但是后者可以自动确定误差界，适用性更强。

本文提出一种基于信号子空间的迭代加权最小方差算法，该方法首先对多快拍信号做奇异值分解提取信号子空间，得到低维的接收信号矩阵，再采用迭代加权最小方差方法求解欠定的稀疏表示问题。本文算法本质上是对文献［8］算法的改进，考虑了噪声对稀疏解的影响，对代价函数进行改进，使得算法对参数不敏感。所提方法角度分辨性能好、DOA估计精度高，且对具有任意相关性的信号能直接进行处理，不需要常规的去相关处理，另一方面该方法虽然需要已知信源数目，但是在信源数估计不准的情况下依然具有很好的性能。

1 信号模型

为简化模型，假设一均匀线阵有M个阵元，阵元间距d＝λ／2，其中λ表示雷达的工作波长。有随机分布在远场的P个信号源si（t），i＝1，2，…，P，分别以方向θi入射到M个阵元上。则阵列接收到的信号可以表示为

式中，y（t）为阵列在t时刻接收到的M×1维信号矢量；n为噪声；S（t）为P×1维的信号源矢量，S（t）＝［s1（t），s2（t），…，sP（t）］T；A（θ）为M×P维阵列流形矩阵，A（θ）＝［a（θ1），a（θ2），…，a（θP）］，其中

对观测数据在t＝tj，j＝1，2，…，L时刻进行采样，将所得的L个快拍进行列累积构成新的接收数据矩阵Y＝［y（t1），y（t2），…，y（tL）］。

图1 基于稀疏表示的DOA估计方法示意图

根据稀疏表示理论，对于某一单快拍接收信号矢量y＝y（ti），在不考虑噪声的情况下，可以表示为

式中，A∈CM×N（M≪N）为超完备的阵列流形矩阵；Θ＝（θ1，θ2，…，θN）为待选方向角度组成的集合，基于稀疏表示的DOA估计模型如图1所示。如果没有关于目标信息的先验知识，一般将整个空间平均分成N份。一般来说，Θ的选择应当满足θ＝｛θ1，θ2，…，θP｝∈Θ或者θ中的所有元素在Θ中都存在对应的近似值。因此，根据稀疏表示理论，求解目标角度等价于搜索最稀疏的组合矢量x，即

式中，‖x‖0表示x中非零元素的个数，由于求解是组合优化问题，属于非确定性多项式（nondeterministic polynomialhard，NP-hard）难题［13］，计算量极大。

本文提出一种基于迭代最小方差的方法来近似逼近式（2）的解，该方法基于一种新的稀疏性测度定义，即：x含有最多均匀微小元素，而不是最多的零元素。其本质上是在传统的迭代l2-范数基础上，加入了对噪声的约束，达到提高算法性能的目的。

2 单快拍稀疏表示

假设只有单快拍数据y用于估计目标源方位，在加性高斯噪声情况下将稀疏解分解为两部分，得到新的线性约束方程：

式中，ys为理想情况下的接收信号；xs为对应的稀疏解（A（Θ）·xs＝ys）；n为噪声，xn为n的对应解（xn在xs非零元素的位置上为0），满足A（Θ）xn＝n。为减少xn中出现较大值对最终结果x的影响（例如伪峰），期望xn中的非零元素尽可能的小且均匀。而实际中，也可以假设噪声在整个空间中都是均匀或近似均匀的，那么用超完备基A（Θ）均匀的表示噪声n是可行的。因此，稀疏表示问题即可转换为求解下述无约束优化问题：

式中，λ为正则化参数，对于平衡加权方差var（xn）和残留‖Ax-y‖22大小至关重要，λ越小，var（xn）也就越大，反之亦然；‖xn‖22表示xn的2-范数，以避免xn的值出现大而均匀的情况，β为对应的权值；‖·‖2代表向量的2-范数；var（xn）为xn的方差，var（xn）越小，xn中的元素也就越均匀，其可以写作：

式中，1为元素全为1的列向量，C＝I-（1／N）11H是计算向量方差的N×N系数矩阵。然而实际中xs和xn是不可分离的，因此假设存在一个加权矩阵使得xn＝Wx，从而隔离两者之间的影响，其中W＝diag（w），w为仅由｛0，1｝组成的向量，0位于xs的非零元素位置上，其余位置上为1。假设W已知，即假设信号方位已知，但是信号的幅度未知，根据式（4）及新的稀疏性的定义，稀疏表示问题可以表示如下：

求解式（6）可以得到：

可以看出权矩阵W对于求解x至关重要，而实际中这样的W是不可能已知的。因此，采用加权迭代的方法来逼近W，式（7）可以改写成迭代的形式：

式中，wk＝｜xk｜／max（｜xk｜），wint为初始值，｜xk｜表示对xk逐项求绝对值所得的向量，此时wk中的所有元素位于区间［1，0）。由于方差矩阵为半正定矩阵，对C加载一个对角矩阵相当于不仅约束微小值的均匀性，而且对微小值大小也进行了约束，使得xn小且均匀，达到减少伪峰的目的。为了简化参数选择，一般固定β／λ＝β，则式（8）可写作：

与经典的FOUCSS相比，本文算法使用了矩阵（λC＋βI），而不是单位矩阵I，相当于引入了对噪声项的约束，因此算法性能比FOUCSS较好。而本文算法具有与FOUCSS具有相同的收敛性且证明方法相似，文献［7］证明了其收敛性，用同样的方法也可以证明本文算法的收敛性，此处不再证明。

3 多快拍联合稀疏表示

当多快拍数据可获的时，其模型可扩展为

式中，Y＝［y（t1），y（t2），…，y（tL）］∈CM×L为数据矩阵；X＝［x1，x2，…，xL］∈CN×L为对应的稀疏解矩阵；Ψ为噪声。当有多快拍数据可用时，多快拍联合稀疏处理自然比单快拍更为稳健，这里，联合稀疏指各信号源的入射角度在不同快拍数据间保持不变，变化的仅仅是幅度。当快拍数据个数L较多时，传统的联合稀疏表示所需的计算量随着快拍数据个数超线性增加，不适用于实时化DOA估计要求。

本文采用基于奇异值分解分解（singular value decomposition，SVD）的迭代最小方差方法对多快拍信号进行联合处理，首先对多快拍信号进行SVD并将观测数据转化到信号子空间上，只对信号子空间进行处理，由于目标源个数远小于快拍数，所以大大提高了计算效率。首先对数据矩阵进行SVD分解

式中，对角矩阵Λ代表Y的奇异值；矩阵U和V分表代表左右酉矩阵。假设已知实际信号源的个数为P（P≪N），则向信号子空间投影后的观测数据矩阵YSVD＝UPΛP，此时的数据规模为M×P，其中UP为U的前P列，ΛP为对角矩阵Λ前P个元素组成的P×P对角矩阵。由上所述可知新的观测模型为

对于新生成的低维子空间XSVD是否满足联合稀疏性，给出如下证明：

命题1 已知Y＝A（Θ）X，X满足列向量之间联合稀疏，则经过降维处理后的新模型YSVD＝AXSVD，XSVD也满足列向量联合稀疏，且＝［XSVD，X］也满足联合稀疏。

证明 令＝YV＝UΛ，则

式中，XSVD＝XVJP，已知X列向量满足之间联合稀疏，对其右乘一个非零矩阵，不影响其列稀疏性及稀疏结构，因此命题得证。

通过SVD降维后的模型保持了原模型的稀疏性及稀疏结构，因此不影响其对信号方向的估计，且经过SVD估计得到的信号子空间中，多快拍数据中信号能量得到累计，信噪比大，能更加有效地估计DOA。对降维后的数据模型进行联合稀疏表示，结合式（9）得到

其中wk为加权矩阵Wk的对角线元素，令

4 初始值与参数选择

初始加权矩阵W对与式（9）和式（11）能否正确收敛以及收敛速度都是至关重要的，因此如何选择一个合适的初始值对于所提算法也同样是关键。理论上，当初始值越接近真实值，收敛速度也就越快，结果也就更加的逼近真实值。而实际中一个合适的初始值一般很难获得。本文算法一个优点在于，当没有合适的初始值时，令Wint＝I，也能获得非常好的结果。在下面的实验中，均采用Wint＝I，实验结果表明该初始值能很好地保证算法收敛。

另外，正则化参数对本文算法的性能也起着至关重要的作用，这里针对两种情况，给出相应的参数选择的方法。

（1）当剩余残差的上界ε2（‖YSVD-A（Θ）XSVD‖F≤ε2）已知时，根据‖YSVD-A（Θ）Xreg‖F≤ε2来确定正则化参数［14］，其中

（2）当剩余残差未知时，正则化参数选择需要合理的权衡稀疏解的稀疏性和剩余残差的大小，然而如何定义什么是合理，目前还没有公认的准则。针对单快拍稀疏表示问题，文献［15］提出了L-curve法来确定正则化参数。在多快拍联合处理时，文献［16］对L-curve法进行推广和改进，提出了修正的L-curve法来确定β。修正L-curve需要已知信噪比的最小可能范围，才能确定准确的参数。

以上两种方法可以分别适用于本文算法，而本文算法另一个优点在于对于参数的不敏感性，即已知信噪比情况下，根据上述方法确定一个合适的参数，在以后的迭代中保持不变。当信噪比大于-10dB时，实验表明，β＝10-3能取得良好的性能。在信噪比发生变化时，对性能影响不大，因此参数对信噪比变化相对稳健。

5 算法分析及计算复杂度

本文算法正则化参数在很大信噪比范围内皆能取得很好的DOA估计性能，从理论上说有两个原因：①对每一步迭代中的权值矩阵进行归一化并且限定权值向量中的元素值域在［1，μ］区间内。每步迭代前对权值向量的取值范围进行限定，避免了因噪声和计算误差引起的极小和极大值对迭代结果的影响，提高了算法的稳健性。②对噪声部分的均匀性进行约束，使得微小值更加平滑。对噪声的约束，使得在迭代过程中wk中的微小值更加均匀，减少极小值的出现。应当注意到，本文提出的迭代加权最小方差算法与文献［8］算法相比，在于矩阵C的选取不同。矩阵C与单位阵的实质在于是否考虑噪声对稀疏结果的影响。本文算法不但可以对信号子空间进行稀疏恢复，也可以对最近提出的样本协方差矩阵进行稀疏表示，两种模型本质上是一样的。

计算复杂度是衡量算法可行性的一个重要指标，对其进行分析十分必要。本文算法的计算量主要来自与两方面：①对数据矩阵进行SVD分解，估计信号子空间，假设信号源个数已知为P，那么估计信号源子空间所需计算量为O（PML）；②应用加权迭代算法求解稀疏表示问题，估计目标源方位。对于迭代算法只需讨论单次迭代所需的计算量，其计算量主要在于计算一个N×N矩阵的逆，计算量大约为O（N3）。而传统的l1-SVD方法所需的计算量为O（P3N3）＋O（PML），当目标数较多时，l1-SVD所需的计算量明显大于所提算法。

6 仿真实验

实验中，采用M＝16元均匀线阵，阵元间距为半波长。假设有6个远场目标源，信号到达角为：［-20.1，-16.2，-5.3，0，17，22］。设接收机噪声为零均值复高斯噪声，各信号源的辐射功率相同。信噪比设为SNR＝0dB，平稳快拍观测数L＝300，将方位角按-90°～90°等间隔的分为181份，角度间隔为1°来构造超完备阵列流形矩阵A（Θ）。

实验1 信号相关性对所提算法的影响

图2和图3为本文算法分别在相关信号和不相关信号时随所用子空间维数变化所得空域谱，由图可知，不论信号是否相关本文算法均能有效估计出信号源方位，而且与传统高分辨方法不同，对信号子空间维数估计不准不敏感。这是由于本文算法对噪声也进行约束，当子空间中包括噪声时，噪声被完备空间平均表示了，进而减少了伪峰的出现，提高了基于稀疏表示的DOA估计性能。图4为在正则化参数固定时，本文算法随信噪比变化时的输出空域谱。由图可知，本文算法参数在信噪比变化时，可以保持不变，在DOA估计中可以预先设定，不需要利用L-curve法在每次迭代中分别确定。

图2 性能随信号空间维数P变化（信号之间不相关）

图3 性能随信号空间维数P变化（信号相关（-20°，-16°）；（17°，22°）方向信号分别相关）

图4 随信噪比变化时所得DOA谱（参数β＝10-3）

实验2 与传统算法比较

图5给出了不同算法在信号相关和不相关时的空域谱，其中样本数L＝300，SNR＝0dB，信号方位角为［-20°，-16°］。当信号相关时，基于子空间的高分辨方法MUSIC和CAPON法均不能正确检测目标方位，而基于稀疏表示的本文算法、M-FOCUSS［8］和l1-SVD［14］算法均能有效估计目标方位。但是M-FOCUSS对噪声比较敏感，容易产生伪峰，对目标个数的估计造成困难。当信号不相关时，MUSIC算法能正确的估计出目标源方位（目标个数已知），而由于目标相距过近，CAPON法未能有效分辨两个目标；而基于稀疏表示的本文算法、l1-SVD和M-FOCUSS算法性能不受信号相关性的影响，验证了上述推论。

图6和图7中，从统计意义上衡量本文算法的性能，并与l1-SVD进行比较。图中每个点都是由200次独立的蒙特卡罗实验所得。当估计所得的信号方位角度与实际角度相差的绝对值和小于2°时，即为一次正确检测。那么，正确检测的次数与实验总次数的比值即定义为检测概率。图6给出了4种算法随信噪比变化时的检测概率，由图可以看出，本文算法正确检测来波方向的概率高于l1-SVD，而且在低信噪比情况下也能以较高概率估计目标方位。另一方面l1-SVD需要已知信噪比来确定正则化参数，文献［14］提出一种确定l1-SVD正则化参数的方法，实验表明该方法在高信噪比和大样本下效果比较好，在低信噪比下性能较差。

图5 不同算法所得空域谱

图6 检测概率随信噪比变化（相关系数为0.98）

图7为检测概率随角度间隔的变化曲线，其中SNR＝0dB，L＝300，两目标信号相关（相关系数为0.98）。显然本文算法能检测相距较近的目标，而此时l1-SVD算法性能容易受到正则化参数的影响，很难正确的估计出目标信号方向；当目标信号相距较远时，两种算法性能相当；MFOUCSS算法虽然能检测方位上较近的目标，但是由于伪峰的影响，稳健性较差；而传统的MUSIC和CAPON算法对于强相关信号，检测性能很差。图8所示为检测概率随样本变化曲线，可以看出本文算法在样本较少时可以较好的估计出目标方位，总体上要优于l1-SVD和M-FOUCSS算法。

图7 检测概率随角度间隔变化（相关系数为0.98）

图8 检测概率随样本个数变化（SNR＝0dB）

7 结 论

本文提出一种基于多快拍稀疏表示的DOA估计方法。该算法给出对稀疏性测度的一种新定义，并应用迭代加权最小方差方法来求解该稀疏解，以达到估计目标源方位的目的。由于该方法没有利用样本的统计信息，因而与传统DOA方法不同，能估计具有任意相关性目标信号，且不需已知相关信号个数。相比于M-FOUCSS算法，本文方法引入了对噪声均匀性的约束，极大地降低了伪峰的出现，提高了目标的检测概率。该方法参数选则策略简单，且参数选择对信噪比变化不敏感，具有很强的适应性和稳健性。

［1］Krimand H，Viberg M.Two decades of array signal processing research：the parametric approach［J］.IEEESignalProcessing，1996，13（4）：67-94.

［2］He Z S，Cichocki A，Zdunek R，et al.Improved FOCUSS method with conjugate gradient iterations［J］.IEEETrans.onSignal Processing，2009，57（1）：309-404.

［3］Ptucha R，Savakis A.Manifold based sparse representation for facial understanding in natural images［J］.ImageandVision Computing，2013，31（5）：365-378.

［4］Cheng H，Liu Z，Yang L，et al.Sparse representation and learning in visual recognition：theory and application［J］.Signal Processing，2012，93（6）：1408-1425.

［5］Wang B，Liu J J，Sun X Y.Mixed sources localization based on sparse signal reconstruction［J］.IEEESignalProcessingLetters，2012，19（8）：487-490.

［6］Dai J S，Xu X，Zhao D.Direction-of-arrival estimation via realvalued sparse representation［J］.IEEEAntennasandWireless PropagationLetters，2013，12：376-379.

［7］Gorodnitsky I F，Rao B D.Sparse signal reconstruction from limited data using FOCUSS：a re-weighted minimum norm algorithm［J］.IEEETrans.onSignalProcessing，1997，45（3）：600-616.

［8］Cotter S F，Rao B D，Engan K.Sparse solutions to linear inverse problems with multiple measurement vectors［J］.IEEETrans.onSignalProcessing，2005，53（7）：2477-2488.

［9］Yin J H，Chen T Q.Direction-of-arrival estimation using a sparse representation of array covariance vectors［J］.IEEE Trans.onSignalProcessing，2011，59（9）：4489-4493.

［10］Xu X，Wei X H，Ye Z F.DOA estimation based on sparse signal recovery utilizing weighted-norm penalty［J］.IEEESignal Processingletters，2012，19（3）：155-158.

［11］He Z Q，Liu Q H，Jin L N，et al.Low complexity method for DOA estimation using array covariance matrix sparse representation［J］.ElectronicsLetters，2013，49（3）：228-230.

［12］Li P F，Zhang M，Zhong Z F，et al.Broadband DOA estimation based on sparse representation in spatial frequency domain［J］.JournalofElectronics＆InformationTechnology，2012，34（2）：404-409.（李鹏飞，张曼，钟子发，等.基于空频域稀疏表示的宽频段 DOA 估计［J］.电子与信息学报，2012，34（2）：404-409.）

［13］Richard G B.Compressive sensing［J］.IEEESignalProcessing Magazine，2007：118-124.

［14］Malioutov D，Cetin M，Willsky A S.A sparse signal reconstruction perspective for source localization with sensor arrays［J］.IEEETrans.onSignalProcessing，2005，53（8）：3010-3022.

［15］Hansen P C，O’Leary D P.The use of the L-curve in the regularization of discrete ill-posed problems［J］.SIAMJournalof ScienceComputationation，1993，14（6）：1487-1503.

［16］Rao B D，Engan K，Cotter S F.Subset selection in noise based on diversity measure minimization［J］.IEEETrans.onSignal Processing，2003，51（3）：760-770.