張?zhí)旌?，季華益，曾德國

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

0 引言

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

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

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

1 數(shù)學(xué)模型

1.1 陣列數(shù)據(jù)DOA估計(jì)數(shù)學(xué)模型

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

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

其中第i個(gè)導(dǎo)向矢量為：

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

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

1.2 基于稀疏表示的DOA估計(jì)模型

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

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

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

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

式中，Nf為寬帶信號(hào)頻點(diǎn)個(gè)數(shù)。當(dāng)角度空間劃分為N個(gè)網(wǎng)格時(shí)，Ψ為M×（Nf×N）維矩陣，隨著頻段的展寬，稀疏分解的計(jì)算量成指數(shù)倍增長。若fisinθi＝fjsinθj，則字典中α（fi，θi）＝α（fj，θj），此時(shí)字典的相關(guān)度為1，算法的性能急劇下降。

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

2.1 空間頻率模型的建立

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

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

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

考察矩陣X（f）中的任一行，假設(shè)在頻點(diǎn)fi處有峰值，則取出各通道在該頻點(diǎn)處對(duì)應(yīng)的頻域值組成列向量Y，構(gòu)造稀疏表達(dá)式：

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

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

寬帶接收機(jī)所接收到的一幀數(shù)據(jù)中可能包含多個(gè)頻率的信號(hào)（包括窄帶和寬帶）。如果將整個(gè)帶寬劃分為若干子頻帶，再分別對(duì)每個(gè)子頻帶進(jìn)行處理，則會(huì)對(duì)沒有信號(hào)的子頻段進(jìn)行處理，浪費(fèi)計(jì)算資源。本文提出的方法根據(jù)信號(hào)在頻域內(nèi)的分布特點(diǎn)，自適應(yīng)地進(jìn)行信號(hào)分離，然后在有信號(hào)的頻點(diǎn)處進(jìn)行DOA估計(jì)。方法如下：

1）對(duì)陣列接收到的數(shù)據(jù)進(jìn)行傅立葉變換得到信號(hào)的頻域譜。任意選取一個(gè)陣列單元的頻域譜進(jìn)行譜峰搜索，記錄下峰值及其對(duì)應(yīng)的頻點(diǎn)fi。

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

3）對(duì)上式進(jìn)行稀疏分解，得到頻率為fi處的空間頻率γ＝［γ1，γ2，…，γK］，K為信號(hào)個(gè)數(shù)。則信號(hào)的DOA為：

重復(fù)以上步驟，即可得到寬頻段內(nèi)所有窄帶信號(hào)的DOA。

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

2.3 稀疏表達(dá)式的求解

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

最小l0范數(shù)是一個(gè)NP難題，而文獻(xiàn)［11］指出，采用l1范數(shù)代替l0范數(shù)可以得到如下凸優(yōu)化問題：

文獻(xiàn)［12］證明了在滿足一定條件時(shí)式（12）和式（13）的等價(jià)性，因此，信號(hào)重構(gòu)問題可以轉(zhuǎn)換為一個(gè)線性規(guī)劃問題加以求解。如果考慮存在噪聲的情況，上述問題可以轉(zhuǎn)化為如下最小l1范數(shù)問題：

式中，σ是與噪聲能量有關(guān)的參數(shù)。

本文采用基于梯度投影的BP算法求解最小l1范數(shù)優(yōu)化問題來重構(gòu)稀疏信號(hào)。得到稀疏信號(hào)后，進(jìn)而根據(jù)非零元素的位置求得信源DOA的估計(jì)值。

3 算法仿真

3.1 空間頻率法與直接恢復(fù)法和傳統(tǒng)MUSIC法的比較

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

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

可以看到，空間頻率法的估計(jì)誤差與傳統(tǒng)的MUSIC法相當(dāng)，但遠(yuǎn)優(yōu)于直接恢復(fù)法；直接恢復(fù)法所要求解的稀疏矩陣是1801×200維的（該矩陣的每一列都有相同的稀疏結(jié)構(gòu)），相比而言空間頻率法求解的稀疏矩陣則為1001×1維，且無需MUSIC法的空間譜峰搜索，因而空間頻率法比另兩種算法計(jì)算量更小。

3.2 相干信號(hào)DOA估計(jì)

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

圖2 相干信號(hào)DOA估計(jì)結(jié)果

對(duì)相干信號(hào)而言，傳統(tǒng)的MUSIC算法已經(jīng)基本失效，空間頻率法與空間平滑MUSIC法均能對(duì)目標(biāo)角度進(jìn)行正確估計(jì)。為比較空間平滑MUSIC法和本文算法的估計(jì)精度，以1dB為步長，考察信噪比為0～20dB時(shí)兩種算法的性能，其它條件不變，每個(gè)信噪比做500次蒙特卡洛實(shí)驗(yàn)。仿真結(jié)果如圖3所示。

圖3 相干信號(hào)估計(jì)均方根誤差隨信噪比變化

在低信噪比情況下，本文算法的估計(jì)精度要優(yōu)于平滑MUSIC法，在SNR＞10dB時(shí)兩種算法性能相當(dāng)。

3.3 寬帶相干信號(hào)DOA估計(jì)

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

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

由圖4可以看出，本文方法的均方根誤差與CSM法相當(dāng)?shù)∮贗SM法，在信噪比大于5dB的情況下，均方根誤差小于0.5°。由于CSM法需要預(yù)估角度，而預(yù)估角的選擇對(duì)最終的估計(jì)結(jié)果有較大影響，因此綜合來看，本文提出的空間頻率法較優(yōu)。

4 結(jié)束語

適合寬帶接收機(jī)的高精度、高分辨率寬帶DOA估計(jì)算法具有重要的理論研究價(jià)值和軍事應(yīng)用前景。本文提出了一種基于空間頻率系數(shù)分解的寬帶DOA估計(jì)方法。仿真實(shí)驗(yàn)表明，該方法對(duì)寬頻帶內(nèi)的多個(gè)窄帶和寬帶信號(hào)、相干和非相干信號(hào)的DOA進(jìn)行準(zhǔn)確估計(jì)，其精度、對(duì)噪聲的魯棒性均優(yōu)于傳統(tǒng)的MUSIC法，且明顯降低了計(jì)算量，是一種適應(yīng)性很強(qiáng)的超分辨算法。■

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

［2］ 張小飛，汪飛，徐大專.陣列信號(hào)處理的理論和應(yīng)用［M］.北京：國防工業(yè)出版社，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.