一種基于稀疏表示的DOA估計新方法

張天恒，季華益，曾德國

（中國航天科工集團8511研究所，江蘇南京210007）

0 引言

對輻射源進行精確的測向定位是現代戰爭中的重要作戰環節之一，因此，波達方向（DOA）估計理論和算法的研究具有極大的軍事應用價值［1］。實際的測向裝備面臨的是一個很寬的頻段，同時空間中本身包含著諸如線性調頻、擴頻等寬帶信號，如何對這些信號的DOA進行高精度的估計是一個巨大的挑戰。目前經典的DOA估計方法是基于陣列的空間譜估計法。以MUSIC為代表的信號子空間算法使測向定位技術突破了分辨率的限制，但該算法有計算量大、在低信噪比條件下性能不佳等缺點，且該算法不能直接對相干信號進行處理，而利用空間平滑的方法解相干則會損失一定的陣列孔徑［2］。

與經典的譜估計法相比，基于稀疏表示的DOA估計法具有很高的估計精度，無需任何預處理，并可直接應用于相干信號，因而得到了國內外學者的廣泛關注。Malioutov等［3］提出了利用頻率和角度構建聯合字典的方法。但該方法隨著頻段的展寬，字典的長度隨之增加。這一方面增加了計算量和存儲量；另一方面，在寬頻段范圍內很難保證字典內沒有相同的元素，當存在相近的元素時，字典的相干性不能滿足，稀疏分解的精度急劇下降。

針對以上問題，本文引入空間頻率的概念，提出了一種基于空間頻率的DOA估計方法。該方法的字典在整個帶寬內恒定不變。利用陣列接收信號的頻域數據建立稀疏模型，可以同時對帶寬范圍內的多個窄帶和寬帶信號進行高精度DOA估計。

1 數學模型

1.1 陣列數據DOA估計數學模型

假設K個遠場信號入射到陣元間距為d的M個均勻線陣上。第i（i＝1，2，…，K）個信號的方位角為θi，頻率為fi。M個陣元接收到的數據可以表示為一個矢量：

式中，A（f，θ）是M×K維陣列流形矩陣，S（t）＝［sf1，θ1（t），sf2，θ2（t），…，sfK，θK（t）］T為K×1維空間信號矩陣，噪聲矩陣N（t）是與信號無關的M×1維高斯白噪聲。陣列流形矩陣為：

其中第i個導向矢量為：

式中，c為電磁波傳播速度。

寬帶DOA估計算法大都由窄帶算法發展而來。現有的基于信號子空間的寬帶DOA估計方法主要分為兩類：非相干子空間法（ISM）和相干子空間法（CSM）。ISM［4］將整個寬帶頻段劃分為若干個窄帶子頻段，根據每個子頻段的中心頻率分別建立方向矢量，這大大增加了存儲空間和計算量；CSM［5］則需要構造聚焦矩陣，將不同頻率的信號子空間映射到同一個參考頻率上，然后將所有頻率成分的信號功率譜密度矩陣做平均。該方法需要對角度進行預估計，當預估值與真實值相差較大時算法性能較差。這兩類算法都需要對空間譜峰值進行搜索，且不能直接用于相干信號，而利用空間平滑的方法解相干則會損失陣列孔徑。

1.2 基于稀疏表示的DOA估計模型

由于空間輻射源的個數相對于整個角度空間來說是有限的，因此輻射源在角度空間的分布具有稀疏性［6］。Mallat和Zhang首次提出了基于過完備原子庫的信號表示［7］，利用其進行DOA估計的主要思想是將信號的陣列流形矩陣擴展成一個過完備的字典Φ，它包含了輻射源所有可能的角度。令θ＝｛θ1，θ2，…，θN｝，N為整個角度空間劃分的網格數，則陣列流形矩陣A可以擴展成如下形式：

從Φ中找到具有最佳線性組合的m項原子來表示原信號，由于m?N，因此該過程稱為信號的稀疏表示。

定義N×1維矢量h＝［h1，h2，…，hN］，當且僅當信號源位于θn時，hn不等于零。于是DOA估計的稀疏模型為：

利用匹配追蹤（MP）或基追蹤［8］（BP）等稀疏分解算法求得h，再根據h中非零元素的位置即可得到信號的DOA。在進行寬帶DOA估計時，Malioutov提出將寬帶范圍內的若干頻點分別建立字典，再將若干個字典組合成一個大字典的方法，即：

式中，Nf為寬帶信號頻點個數。當角度空間劃分為N個網格時，Ψ為M×（Nf×N）維矩陣，隨著頻段的展寬，稀疏分解的計算量成指數倍增長。若fisinθi＝fjsinθj，則字典中α（fi，θi）＝α（fj，θj），此時字典的相關度為1，算法的性能急劇下降。

2 基于空間頻率稀疏表示的寬帶DOA估計

2.1 空間頻率模型的建立

為解決寬帶信號字典隨帶寬展寬而變大的問題，令γi＝fidsinθi／c，i＝1，2，…，N，則：

為了保證無模糊測向，要求陣元間距小于波長的一半，即γ的取值范圍是－0.5≤γ≤0.5。根據空間頻率的取值范圍劃分N個網格，則過完備字典可表示為：通過求解式（5）即可得到輻射源的空間頻率。

由于空間頻率是頻率和角度的二維函數，因此必須在頻率已知時才能求出DOA。對矩陣X和S做快速傅里葉變換后的結果分別為X（f）和S（f），則有：

考察矩陣X（f）中的任一行，假設在頻點fi處有峰值，則取出各通道在該頻點處對應的頻域值組成列向量Y，構造稀疏表達式：

通過求出向量h中非零值所在的位置，即可得到頻率為fi的所有信號的DOA。

2.2 基于稀疏表示的空間頻率DOA估計方法

寬帶接收機所接收到的一幀數據中可能包含多個頻率的信號（包括窄帶和寬帶）。如果將整個帶寬劃分為若干子頻帶，再分別對每個子頻帶進行處理，則會對沒有信號的子頻段進行處理，浪費計算資源。本文提出的方法根據信號在頻域內的分布特點，自適應地進行信號分離，然后在有信號的頻點處進行DOA估計。方法如下：

1）對陣列接收到的數據進行傅立葉變換得到信號的頻域譜。任意選取一個陣列單元的頻域譜進行譜峰搜索，記錄下峰值及其對應的頻點fi。

2）取出所有陣列單元在頻點fi處的值組成一個列向量Y，按2.1節中的方法構造空間頻率矩陣Φ，得到模型Y＝Φh。

3）對上式進行稀疏分解，得到頻率為fi處的空間頻率γ＝［γ1，γ2，…，γK］，K為信號個數。則信號的DOA為：

重復以上步驟，即可得到寬頻段內所有窄帶信號的DOA。

對于有一定帶寬的信號，假設根據頻譜判斷其帶寬為B，將帶寬分為P段，每一段的中心頻率為fp，p＝1，2，…，P。按照上面的方法構建稀疏模型，得到每個子頻段的DOA估計值（fp）。對所有的估計值求平均，即得到寬帶信號的DOA的最終結果

2.3 稀疏表達式的求解

求解式（10）的本質是從M個方程中求解出N個未知量（M?N），從表面上看，這是個無法直接求解的欠定方程，然而文獻［9］指出，由于信號矢量h是稀疏的，若矩陣Φ滿足有限等距準則（RIP）等稀疏重構條件，則該信號時可以以很高的概率被重構出來。Candès等證明［10］，信號重構問題可以通過求解以下最小l0范數問題加以解決。

最小l0范數是一個NP難題，而文獻［11］指出，采用l1范數代替l0范數可以得到如下凸優化問題：

文獻［12］證明了在滿足一定條件時式（12）和式（13）的等價性，因此，信號重構問題可以轉換為一個線性規劃問題加以求解。如果考慮存在噪聲的情況，上述問題可以轉化為如下最小l1范數問題：

式中，σ是與噪聲能量有關的參數。

本文采用基于梯度投影的BP算法求解最小l1范數優化問題來重構稀疏信號。得到稀疏信號后，進而根據非零元素的位置求得信源DOA的估計值。

3 算法仿真

3.1 空間頻率法與直接恢復法和傳統MUSIC法的比較

考慮6陣元均勻線陣，假設遠場入射信號入射角為－30.3°，頻率800MHz，采樣率1GHz，陣元間距為波長的一半，時域快拍數K＝200，對時域數據做1024點FFT得到頻域快拍。文獻［1］采用的直接恢復法取所有的200個時域快拍構建稀疏表達式，空間域范圍－90°～90°，網格劃分為1801個，以多測量矢量法（MMV）求解；空間頻率法的空間頻率范圍－0.5～0.5，網格劃分為1001個，取信號頻率所在的頻點構建稀疏表達式，以單測量矢量法（SMV）求解。以為步長，考察信噪比為0～20dB時兩種算法的性能，每個信噪比做500次蒙特卡洛實驗。仿真結果如圖1所示。

圖1 均方根誤差隨信噪比變化

可以看到，空間頻率法的估計誤差與傳統的MUSIC法相當，但遠優于直接恢復法；直接恢復法所要求解的稀疏矩陣是1801×200維的（該矩陣的每一列都有相同的稀疏結構），相比而言空間頻率法求解的稀疏矩陣則為1001×1維，且無需MUSIC法的空間譜峰搜索，因而空間頻率法比另兩種算法計算量更小。

3.2 相干信號DOA估計

假設有兩個遠場相干窄帶入射信號，其入射角度分別－20.6°和7.8°，頻率為800MHz。其它條件如前所述。SNR＝10dB，分別采用MUSIC法、空間平滑MUSIC法和空間頻率法估計DOA。仿真如圖2所示。

圖2 相干信號DOA估計結果

對相干信號而言，傳統的MUSIC算法已經基本失效，空間頻率法與空間平滑MUSIC法均能對目標角度進行正確估計。為比較空間平滑MUSIC法和本文算法的估計精度，以1dB為步長，考察信噪比為0～20dB時兩種算法的性能，其它條件不變，每個信噪比做500次蒙特卡洛實驗。仿真結果如圖3所示。

圖3 相干信號估計均方根誤差隨信噪比變化

在低信噪比情況下，本文算法的估計精度要優于平滑MUSIC法，在SNR＞10dB時兩種算法性能相當。

3.3 寬帶相干信號DOA估計

假設兩個寬帶相干LFM信號，中心頻率750MHz，帶寬100MHz，入射角分別為－22.8°和6.2°。ISM法和空間頻率法將帶寬分為5段，在每一段都求出一個估計值，將它們的平均值作為最終結果；CSM法構造聚焦矩陣時參考頻率選為750MHz，角度預估值為－5°。SNR＝0～20dB，步長1dB，每個信噪比處做500次蒙特卡洛實驗。均方根誤差隨信噪比的變化如圖4所示。

圖4 均方根誤差隨信噪比變化

由圖4可以看出，本文方法的均方根誤差與CSM法相當但小于ISM法，在信噪比大于5dB的情況下，均方根誤差小于0.5°。由于CSM法需要預估角度，而預估角的選擇對最終的估計結果有較大影響，因此綜合來看，本文提出的空間頻率法較優。

4 結束語

適合寬帶接收機的高精度、高分辨率寬帶DOA估計算法具有重要的理論研究價值和軍事應用前景。本文提出了一種基于空間頻率系數分解的寬帶DOA估計方法。仿真實驗表明，該方法對寬頻帶內的多個窄帶和寬帶信號、相干和非相干信號的DOA進行準確估計，其精度、對噪聲的魯棒性均優于傳統的MUSIC法，且明顯降低了計算量，是一種適應性很強的超分辨算法。■

［1］ 林波.基于壓縮感知的輻射源DOA估計［D］.長沙：國防科學技術大學，2010.

［2］ 張小飛，汪飛，徐大專.陣列信號處理的理論和應用［M］.北京：國防工業出版社，2010.

［3］ Malioutov D，Cetin M，Willsky AS.A sparse signal reconstruction perspective for source localization with sensor arrays［J］.IEEE Transactions on Signal Processing，2005，53（8）：3010－3022.

［4］ Wang Y，Leus G.Direction estimation using compressive sampling array processing［C］.15th IEEE Workshop on Statistical Signal Processing，2009：626－629.

［5］ Wang H，Kaveh M.Coherent signal－subspace processing for the detection and estimation of angles of multiple wideband sources［J］.IEEE Transactions Acoustics，Speech，and Signal Processing，1985，33（4）：823－831.

［6］ Cotter SF，Rao BD，Engan K，et al.Sparse solutions to linear inverse problems with multiple measurement vectors［J］.IEEE Transactions on Signal Processing，2005，53（7）：2477－2488.

［7］ Mallat S，Zhang Z.Matching pursuits with time－frequency dictionaries［J］.IEEE Transactions Signal Process，1993，41（12）：3397－3415.

［8］ Chen HB，Donoho DL，Saunders MA.Atomic decom－position by basis pursuit［J］.SIAM Joumal on Scientific Compting，2001，20（1）：33－66.

［9］ Donoho DL.Compressed sensing［J］.IEEE Transactions on Information Theory，2006，52（4）：1289－1306.

［10］Candès E.Compressive sampling.In：Proceedings of International Congress of Mathematicians［M］.Madrid，Spain：European Mathematical Society Publishing House，2006：1433－1452.

［11］Candès E，Romberg J，Tao T.Robust uncertainty principles：Exact signal reconstruction from highly incomplete frequency information［J］.IEEE Transactions on Information Theory，2006，52（2）：489－509.

［12］Donoho DL，Elad M，Temlyakov VN.Stable recovery of sparse overcomplete representations in the presence of noise［J］.IEEE Transactions on Information Theory，2006，52（1）：6－18.