采用信號子空間稀疏表示的DOA估計方法

解 虎，馮大政，魏倩茹

（西安電子科技大學雷達信號處理國家重點實驗室，陜西 西 安710071）

0 引 言

波達方向（direction of arrival，DOA）估計是陣列信號處理的重要研究內容之一，被廣泛應用于雷達、無線通信、電磁場、聲納、地震勘探等諸多領域。DOA估計方法的主要目標是在噪聲環(huán)境下，分辨兩個在方位向非常接近的目標，為此，利用信號源個數(shù)小于陣元數(shù)的特點，人們提出了大量高分辨DOA估計方法。例如，經(jīng)典的多重信號分類（multiple signal classification，MUSIC）法，其關鍵是信號子空間維數(shù)的確定。本文也采用類似的方法，首先確定信號子空間的維數(shù)，然后將DOA估計問題構造成一個稀疏表示問題以獲得較高的角度分辨率。

基于稀疏表示的DOA估計方法與已有的方法雖有相同之處，但是在本質上卻有不同。常用的DOA估計方法分為兩類，即參數(shù)化估計和非參數(shù)化估計。對于非參數(shù)估計方法，主要有波束形成法，基于子空間方法的MUSIC和基于最小方差無畸變（minimum variance distortionless response，MVDR）的高分辨譜估計法等。其中波束掃描方法雖然與信噪比（signal-to-noise ratio，SNR）無關，但是受到瑞利限的制約，無法分辨一個波束寬度內的兩個信號。而MUSIC和MVDR方法雖然能夠分辨一個瑞利單元內的兩個目標，但是卻需要較高的SNR和大量的樣本，而且對信號源之間的相干性有一定要求，不能對相干信號源直接進行有效分辨?；谧畲笏迫坏膮?shù)化估計方法分為：確定最大似然和統(tǒng)計最大似然［1］。這類方法雖然有很好的統(tǒng)計特性，但是需要準確的初始值才能收斂到全局最優(yōu)點。而對于基于稀疏表示的DOA估計方法，不需要一個好的初始值和大量樣本，對信號相關性也沒有特別限定，依然能獲得高分辨率。

近幾十年來，稀疏表示取得了重大進展，已經(jīng)被廣泛應用于圖像重建與理解［2-3］、模式識別［4］和高分辨譜估計［5-6］等諸多領域。Gorodnitsky等提出了一種加權迭代最小2-范數(shù)（focal underdetermined system solver，F(xiàn)OCUSS）［7］方法來求解稀疏表示問題，并在DOA估計上取得較好的結果，但是僅適用于單次快拍。然而單快拍信號容易受到信號幅度變化的影響，而且在低信噪比的情況下，正確估計和恢復來波信號的概率較小。實際中，用于DOA估計的樣本個數(shù)比較充足，因此人們提出了多樣本下的基于稀疏表示的DOA估計方法，以達到高分辨率。文獻［8］提出了一種基于奇異值分解的多測量矢量欠定系統(tǒng)正則化聚焦求解算法實現(xiàn)DOA超分辨估計。該算法計算量快，但是對正則化參數(shù)比較敏感。文獻［9］采用l1，2混合范數(shù)對多快拍接收信號的協(xié)方差矩陣進行聯(lián)合稀疏表示，達到DOA估計的目的。由于需要對協(xié)方差矩陣求逆，當樣本較少時，會對導致算法性能下降，但是該算法不需要確定參數(shù)，適用范圍更廣。文獻［10］中，作者采用加權l(xiāng)1范數(shù)確保解的稀疏性，并將其應用于DOA估計。然而其需要Capon估計加權系數(shù)，因此當存在相關信號時，算法性能會下降。文獻［11］提出一種將基于協(xié)方差的多快拍信號聯(lián)合稀疏表示轉化成單一向量的稀疏表示模型，該方法計算量小，但是要求信號之間互不相關。文獻［12］進一步將稀疏表示方法應用于寬帶DOA估計中，在低信噪比下取得了良好的估計精度。目前基于多快拍信號稀疏表示的模型主要分為兩類：對信號子空間進行稀疏表示和對信號協(xié)方差矩陣進行稀疏表示。同等條件下（樣本數(shù)目相同，參數(shù)選擇合理），基于信號子空間的算法性能一般要優(yōu)于基于協(xié)方差矩陣的算法，但是后者可以自動確定誤差界，適用性更強。

本文提出一種基于信號子空間的迭代加權最小方差算法，該方法首先對多快拍信號做奇異值分解提取信號子空間，得到低維的接收信號矩陣，再采用迭代加權最小方差方法求解欠定的稀疏表示問題。本文算法本質上是對文獻［8］算法的改進，考慮了噪聲對稀疏解的影響，對代價函數(shù)進行改進，使得算法對參數(shù)不敏感。所提方法角度分辨性能好、DOA估計精度高，且對具有任意相關性的信號能直接進行處理，不需要常規(guī)的去相關處理，另一方面該方法雖然需要已知信源數(shù)目，但是在信源數(shù)估計不準的情況下依然具有很好的性能。

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）］，其中

對觀測數(shù)據(jù)在t＝tj，j＝1，2，…，L時刻進行采樣，將所得的L個快拍進行列累積構成新的接收數(shù)據(jù)矩陣Y＝［y（t1），y（t2），…，y（tL）］。

圖1 基于稀疏表示的DOA估計方法示意圖

根據(jù)稀疏表示理論，對于某一單快拍接收信號矢量y＝y(tǒng)（ti），在不考慮噪聲的情況下，可以表示為

式中，A∈CM×N（M?N）為超完備的陣列流形矩陣；Θ＝（θ1，θ2，…，θN）為待選方向角度組成的集合，基于稀疏表示的DOA估計模型如圖1所示。如果沒有關于目標信息的先驗知識，一般將整個空間平均分成N份。一般來說，Θ的選擇應當滿足θ＝｛θ1，θ2，…，θP｝∈Θ或者θ中的所有元素在Θ中都存在對應的近似值。因此，根據(jù)稀疏表示理論，求解目標角度等價于搜索最稀疏的組合矢量x，即

式中，‖x‖0表示x中非零元素的個數(shù)，由于求解是組合優(yōu)化問題，屬于非確定性多項式（nondeterministic polynomialhard，NP-hard）難題［13］，計算量極大。

本文提出一種基于迭代最小方差的方法來近似逼近式（2）的解，該方法基于一種新的稀疏性測度定義，即：x含有最多均勻微小元素，而不是最多的零元素。其本質上是在傳統(tǒng)的迭代l2-范數(shù)基礎上，加入了對噪聲的約束，達到提高算法性能的目的。

2 單快拍稀疏表示

假設只有單快拍數(shù)據(jù)y用于估計目標源方位，在加性高斯噪聲情況下將稀疏解分解為兩部分，得到新的線性約束方程：

式中，ys為理想情況下的接收信號；xs為對應的稀疏解（A（Θ）·xs＝y(tǒng)s）；n為噪聲，xn為n的對應解（xn在xs非零元素的位置上為0），滿足A（Θ）xn＝n。為減少xn中出現(xiàn)較大值對最終結果x的影響（例如偽峰），期望xn中的非零元素盡可能的小且均勻。而實際中，也可以假設噪聲在整個空間中都是均勻或近似均勻的，那么用超完備基A（Θ）均勻的表示噪聲n是可行的。因此，稀疏表示問題即可轉換為求解下述無約束優(yōu)化問題：

式中，λ為正則化參數(shù)，對于平衡加權方差var（xn）和殘留‖Ax-y‖22大小至關重要，λ越小，var（xn）也就越大，反之亦然；‖xn‖22表示xn的2-范數(shù)，以避免xn的值出現(xiàn)大而均勻的情況，β為對應的權值；‖·‖2代表向量的2-范數(shù)；var（xn）為xn的方差，var（xn）越小，xn中的元素也就越均勻，其可以寫作：

式中，1為元素全為1的列向量，C＝I-（1／N）11H是計算向量方差的N×N系數(shù)矩陣。然而實際中xs和xn是不可分離的，因此假設存在一個加權矩陣使得xn＝Wx，從而隔離兩者之間的影響，其中W＝diag（w），w為僅由｛0，1｝組成的向量，0位于xs的非零元素位置上，其余位置上為1。假設W已知，即假設信號方位已知，但是信號的幅度未知，根據(jù)式（4）及新的稀疏性的定義，稀疏表示問題可以表示如下：

求解式（6）可以得到：

可以看出權矩陣W對于求解x至關重要，而實際中這樣的W是不可能已知的。因此，采用加權迭代的方法來逼近W，式（7）可以改寫成迭代的形式：

式中，wk＝｜xk｜／max（｜xk｜），wint為初始值，｜xk｜表示對xk逐項求絕對值所得的向量，此時wk中的所有元素位于區(qū)間［1，0）。由于方差矩陣為半正定矩陣，對C加載一個對角矩陣相當于不僅約束微小值的均勻性，而且對微小值大小也進行了約束，使得xn小且均勻，達到減少偽峰的目的。為了簡化參數(shù)選擇，一般固定β／λ＝β，則式（8）可寫作：

與經(jīng)典的FOUCSS相比，本文算法使用了矩陣（λC＋βI），而不是單位矩陣I，相當于引入了對噪聲項的約束，因此算法性能比FOUCSS較好。而本文算法具有與FOUCSS具有相同的收斂性且證明方法相似，文獻［7］證明了其收斂性，用同樣的方法也可以證明本文算法的收斂性，此處不再證明。

3 多快拍聯(lián)合稀疏表示

當多快拍數(shù)據(jù)可獲的時，其模型可擴展為

式中，Y＝［y（t1），y（t2），…，y（tL）］∈CM×L為數(shù)據(jù)矩陣；X＝［x1，x2，…，xL］∈CN×L為對應的稀疏解矩陣；Ψ為噪聲。當有多快拍數(shù)據(jù)可用時，多快拍聯(lián)合稀疏處理自然比單快拍更為穩(wěn)健，這里，聯(lián)合稀疏指各信號源的入射角度在不同快拍數(shù)據(jù)間保持不變，變化的僅僅是幅度。當快拍數(shù)據(jù)個數(shù)L較多時，傳統(tǒng)的聯(lián)合稀疏表示所需的計算量隨著快拍數(shù)據(jù)個數(shù)超線性增加，不適用于實時化DOA估計要求。

本文采用基于奇異值分解分解（singular value decomposition，SVD）的迭代最小方差方法對多快拍信號進行聯(lián)合處理，首先對多快拍信號進行SVD并將觀測數(shù)據(jù)轉化到信號子空間上，只對信號子空間進行處理，由于目標源個數(shù)遠小于快拍數(shù)，所以大大提高了計算效率。首先對數(shù)據(jù)矩陣進行SVD分解

式中，對角矩陣Λ代表Y的奇異值；矩陣U和V分表代表左右酉矩陣。假設已知實際信號源的個數(shù)為P（P?N），則向信號子空間投影后的觀測數(shù)據(jù)矩陣YSVD＝UPΛP，此時的數(shù)據(jù)規(guī)模為M×P，其中UP為U的前P列，ΛP為對角矩陣Λ前P個元素組成的P×P對角矩陣。由上所述可知新的觀測模型為

對于新生成的低維子空間XSVD是否滿足聯(lián)合稀疏性，給出如下證明：

命題1 已知Y＝A（Θ）X，X滿足列向量之間聯(lián)合稀疏，則經(jīng)過降維處理后的新模型YSVD＝AXSVD，XSVD也滿足列向量聯(lián)合稀疏，且＝［XSVD，X］也滿足聯(lián)合稀疏。

證明 令＝YV＝UΛ，則

式中，XSVD＝XVJP，已知X列向量滿足之間聯(lián)合稀疏，對其右乘一個非零矩陣，不影響其列稀疏性及稀疏結構，因此命題得證。

通過SVD降維后的模型保持了原模型的稀疏性及稀疏結構，因此不影響其對信號方向的估計，且經(jīng)過SVD估計得到的信號子空間中，多快拍數(shù)據(jù)中信號能量得到累計，信噪比大，能更加有效地估計DOA。對降維后的數(shù)據(jù)模型進行聯(lián)合稀疏表示，結合式（9）得到

其中wk為加權矩陣Wk的對角線元素，令

4 初始值與參數(shù)選擇

初始加權矩陣W對與式（9）和式（11）能否正確收斂以及收斂速度都是至關重要的，因此如何選擇一個合適的初始值對于所提算法也同樣是關鍵。理論上，當初始值越接近真實值，收斂速度也就越快，結果也就更加的逼近真實值。而實際中一個合適的初始值一般很難獲得。本文算法一個優(yōu)點在于，當沒有合適的初始值時，令Wint＝I，也能獲得非常好的結果。在下面的實驗中，均采用Wint＝I，實驗結果表明該初始值能很好地保證算法收斂。

另外，正則化參數(shù)對本文算法的性能也起著至關重要的作用，這里針對兩種情況，給出相應的參數(shù)選擇的方法。

（1）當剩余殘差的上界ε2（‖YSVD-A（Θ）XSVD‖F(xiàn)≤ε2）已知時，根據(jù)‖YSVD-A（Θ）Xreg‖F(xiàn)≤ε2來確定正則化參數(shù)［14］，其中

（2）當剩余殘差未知時，正則化參數(shù)選擇需要合理的權衡稀疏解的稀疏性和剩余殘差的大小，然而如何定義什么是合理，目前還沒有公認的準則。針對單快拍稀疏表示問題，文獻［15］提出了L-curve法來確定正則化參數(shù)。在多快拍聯(lián)合處理時，文獻［16］對L-curve法進行推廣和改進，提出了修正的L-curve法來確定β。修正L-curve需要已知信噪比的最小可能范圍，才能確定準確的參數(shù)。

以上兩種方法可以分別適用于本文算法，而本文算法另一個優(yōu)點在于對于參數(shù)的不敏感性，即已知信噪比情況下，根據(jù)上述方法確定一個合適的參數(shù)，在以后的迭代中保持不變。當信噪比大于-10dB時，實驗表明，β＝10-3能取得良好的性能。在信噪比發(fā)生變化時，對性能影響不大，因此參數(shù)對信噪比變化相對穩(wěn)健。

5 算法分析及計算復雜度

本文算法正則化參數(shù)在很大信噪比范圍內皆能取得很好的DOA估計性能，從理論上說有兩個原因：①對每一步迭代中的權值矩陣進行歸一化并且限定權值向量中的元素值域在［1，μ］區(qū)間內。每步迭代前對權值向量的取值范圍進行限定，避免了因噪聲和計算誤差引起的極小和極大值對迭代結果的影響，提高了算法的穩(wěn)健性。②對噪聲部分的均勻性進行約束，使得微小值更加平滑。對噪聲的約束，使得在迭代過程中wk中的微小值更加均勻，減少極小值的出現(xiàn)。應當注意到，本文提出的迭代加權最小方差算法與文獻［8］算法相比，在于矩陣C的選取不同。矩陣C與單位陣的實質在于是否考慮噪聲對稀疏結果的影響。本文算法不但可以對信號子空間進行稀疏恢復，也可以對最近提出的樣本協(xié)方差矩陣進行稀疏表示，兩種模型本質上是一樣的。

計算復雜度是衡量算法可行性的一個重要指標，對其進行分析十分必要。本文算法的計算量主要來自與兩方面：①對數(shù)據(jù)矩陣進行SVD分解，估計信號子空間，假設信號源個數(shù)已知為P，那么估計信號源子空間所需計算量為O（PML）；②應用加權迭代算法求解稀疏表示問題，估計目標源方位。對于迭代算法只需討論單次迭代所需的計算量，其計算量主要在于計算一個N×N矩陣的逆，計算量大約為O（N3）。而傳統(tǒng)的l1-SVD方法所需的計算量為O（P3N3）＋O（PML），當目標數(shù)較多時，l1-SVD所需的計算量明顯大于所提算法。

6 仿真實驗

實驗中，采用M＝16元均勻線陣，陣元間距為半波長。假設有6個遠場目標源，信號到達角為：［-20.1，-16.2，-5.3，0，17，22］。設接收機噪聲為零均值復高斯噪聲，各信號源的輻射功率相同。信噪比設為SNR＝0dB，平穩(wěn)快拍觀測數(shù)L＝300，將方位角按-90°～90°等間隔的分為181份，角度間隔為1°來構造超完備陣列流形矩陣A（Θ）。

實驗1 信號相關性對所提算法的影響

圖2和圖3為本文算法分別在相關信號和不相關信號時隨所用子空間維數(shù)變化所得空域譜，由圖可知，不論信號是否相關本文算法均能有效估計出信號源方位，而且與傳統(tǒng)高分辨方法不同，對信號子空間維數(shù)估計不準不敏感。這是由于本文算法對噪聲也進行約束，當子空間中包括噪聲時，噪聲被完備空間平均表示了，進而減少了偽峰的出現(xiàn)，提高了基于稀疏表示的DOA估計性能。圖4為在正則化參數(shù)固定時，本文算法隨信噪比變化時的輸出空域譜。由圖可知，本文算法參數(shù)在信噪比變化時，可以保持不變，在DOA估計中可以預先設定，不需要利用L-curve法在每次迭代中分別確定。

圖2 性能隨信號空間維數(shù)P變化（信號之間不相關）

圖3 性能隨信號空間維數(shù)P變化（信號相關（-20°，-16°）；（17°，22°）方向信號分別相關）

圖4 隨信噪比變化時所得DOA譜（參數(shù)β＝10-3）

實驗2 與傳統(tǒng)算法比較

圖5給出了不同算法在信號相關和不相關時的空域譜，其中樣本數(shù)L＝300，SNR＝0dB，信號方位角為［-20°，-16°］。當信號相關時，基于子空間的高分辨方法MUSIC和CAPON法均不能正確檢測目標方位，而基于稀疏表示的本文算法、M-FOCUSS［8］和l1-SVD［14］算法均能有效估計目標方位。但是M-FOCUSS對噪聲比較敏感，容易產(chǎn)生偽峰，對目標個數(shù)的估計造成困難。當信號不相關時，MUSIC算法能正確的估計出目標源方位（目標個數(shù)已知），而由于目標相距過近，CAPON法未能有效分辨兩個目標；而基于稀疏表示的本文算法、l1-SVD和M-FOCUSS算法性能不受信號相關性的影響，驗證了上述推論。

圖6和圖7中，從統(tǒng)計意義上衡量本文算法的性能，并與l1-SVD進行比較。圖中每個點都是由200次獨立的蒙特卡羅實驗所得。當估計所得的信號方位角度與實際角度相差的絕對值和小于2°時，即為一次正確檢測。那么，正確檢測的次數(shù)與實驗總次數(shù)的比值即定義為檢測概率。圖6給出了4種算法隨信噪比變化時的檢測概率，由圖可以看出，本文算法正確檢測來波方向的概率高于l1-SVD，而且在低信噪比情況下也能以較高概率估計目標方位。另一方面l1-SVD需要已知信噪比來確定正則化參數(shù)，文獻［14］提出一種確定l1-SVD正則化參數(shù)的方法，實驗表明該方法在高信噪比和大樣本下效果比較好，在低信噪比下性能較差。

圖5 不同算法所得空域譜

圖6 檢測概率隨信噪比變化（相關系數(shù)為0.98）

圖7為檢測概率隨角度間隔的變化曲線，其中SNR＝0dB，L＝300，兩目標信號相關（相關系數(shù)為0.98）。顯然本文算法能檢測相距較近的目標，而此時l1-SVD算法性能容易受到正則化參數(shù)的影響，很難正確的估計出目標信號方向；當目標信號相距較遠時，兩種算法性能相當；MFOUCSS算法雖然能檢測方位上較近的目標，但是由于偽峰的影響，穩(wěn)健性較差；而傳統(tǒng)的MUSIC和CAPON算法對于強相關信號，檢測性能很差。圖8所示為檢測概率隨樣本變化曲線，可以看出本文算法在樣本較少時可以較好的估計出目標方位，總體上要優(yōu)于l1-SVD和M-FOUCSS算法。

圖7 檢測概率隨角度間隔變化（相關系數(shù)為0.98）

圖8 檢測概率隨樣本個數(shù)變化（SNR＝0dB）

7 結 論

本文提出一種基于多快拍稀疏表示的DOA估計方法。該算法給出對稀疏性測度的一種新定義，并應用迭代加權最小方差方法來求解該稀疏解，以達到估計目標源方位的目的。由于該方法沒有利用樣本的統(tǒng)計信息，因而與傳統(tǒng)DOA方法不同，能估計具有任意相關性目標信號，且不需已知相關信號個數(shù)。相比于M-FOUCSS算法，本文方法引入了對噪聲均勻性的約束，極大地降低了偽峰的出現(xiàn)，提高了目標的檢測概率。該方法參數(shù)選則策略簡單，且參數(shù)選擇對信噪比變化不敏感，具有很強的適應性和穩(wěn)健性。

［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.（李鵬飛，張曼，鐘子發(fā)，等.基于空頻域稀疏表示的寬頻段 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.