基于廣義近似消息傳遞的快速DOA估計方法

張 俊，張新禹，姜衛(wèi)東，劉永祥，黎 湘

（國防科技大學電子科學學院，湖南 長沙410073）

0 引 言

波達方向（direction of arrival，DOA）估計是陣列信號處理研究的一類重點問題，其廣泛應用于雷達、聲納、無線通信、地震及水下探測等眾多領域。經過近幾十年的發(fā)展，研究人員提出了大量高分辨DOA估計算法。但是這些傳統的DOA估計算法在實際應用中受到諸多限制。例如，以多重信號分類方法（multiple signal classification，MUSIC）為代表的子空間分解類算法通常需要較高的信噪比（signal to noise ratio，SNR）、足夠的快拍數量以及已知信源數量。子空間擬合類算法，如最大似然算法和加權子空間擬合算法等，雖然在低SNR和小快拍數的情況下表現比子空間分解類算法好，但因其有計算量大的缺點而難以實際應用。

近年來，稀疏恢復類DOA算法由于在低SNR、有限快拍數和相關信號等場景中的良好性能，成為研究熱點。此類算法的基本思想是在滿足一定約束條件下，利用過完備的基函數求解可以表示原始信號的稀疏系數向量。稀疏恢復類算法主要包括貪婪追蹤、l （≥0）范數優(yōu)化和稀疏貝葉斯學習（sparse Bayesian learning，SBL）三類算法。貪婪追蹤類算法是在每次迭代中根據相關性選擇與當前殘余信號特征結構匹配的最優(yōu)原子，實現對原始信號的逼近，典型算法有匹配追蹤算法、正交匹配追蹤算法等。雖然此類算法不需要信源數量的先驗信息，也可用于處理相關信號，但是存在計算復雜度高、收斂速度慢等問題。l 范數優(yōu)化算法是將l 范數懲罰項引入目標函數，保證解的稀疏性，防止出現過擬合，常用的算法有LASSO（least absolute shrinkage and selection operator）算法、-奇異值分解（singular value decomposition，SVD）算法等。此類算法通常也具有較高的計算復雜度，同時正則化參數的選擇對性能的影響比較明顯。SBL類方法是從貝葉斯觀點出發(fā)，利用信號的稀疏先驗信息對信號進行稀疏重構。相對于上述兩類方法，SBL算法能夠保證達到全局最優(yōu)且計算效率也有顯著提高。同時，SBL算法的在信號分量之間存在強相關性時仍然具有良好的重構性能。文獻［6-7］提出了基于SBL的雙基地無源雷達DOA估計方法，其思想主要是通過空間幾何關系建立雙站接收數據之間的關系，聯合進行DOA估計，在一定程度上緩解了快拍數的限制，性能要優(yōu)于只使用單站數據的情況，但是在計算后驗概率密度時涉及到階矩陣求逆問題，需要在計算效率和估計精度之間做出平衡。為避免大規(guī)模矩陣求逆問題，文獻［8］提出了一種基于酉變換的近似消息傳遞（approximate message passing，AMP）的SBL-DOA估計算法。但是AMP算法只能解決線性問題，不適用于任意的信號分布，且要求字典矩陣滿足獨立同分布的零均值高斯分布。文獻［12］將AMP算法進行了拓展，提出廣義AMP（generalized AMP，GAMP）算法。GAMP算法對輸入、輸出信號的分布沒有要求，且可以解決非線性問題，適用于具有精確漸近性能保證的非高斯估計問題。文獻［13］進一步證明了在獨立高斯先驗和似然函數的條件下，當字典矩陣奇異值的峰—均值比較小時，可以保證GAMP算法收斂，并給出了一種阻尼GAMP算法，使得在上述條件不滿足時也可以保證算法收斂，極大地拓展了算法的適用范圍。

本文提出一種SBL框架下的基于GAMP的雙基地聯合DOA估計方法。首先，構建了雙站快拍數據下的信號模型。通過數據間的獨立性，利用GAMP網絡將多維聯合后驗概率密度的計算簡化為標量形式的邊緣概率密度計算，并推導了參數更新的閉式表達式。然后，針對離網目標，提出了一種基于梯度的網格節(jié)點更新策略。仿真實驗表明，在低SNR和有限快拍數的情況下，本算法相比于其他算法具有較高的DOA估計精度和較低的計算復雜度。

1 雙基地無源雷達接收信號模型

1.1 系統模型

圖1 雙基地無源雷達系統空間幾何關系Fig.1 Spatial geometry of the bi-static passive radar system

根據空間幾何關系，容易得到第個目標在兩個天線陣列上的入射角滿足：

化簡可得

因此，通過式（2）就可以建立兩個雷達之間接收數據的關系。兩個雷達站多快拍模型的接收數據可以寫成

1.2 貝葉斯框架下的稀疏信號模型

對于信號的先驗概率，根據目標的獨立性假設可以認為、中的元素是相互獨立的，且每一行的元素是獨立同分布的，服從均值為零、方差為γ≥0（＝1，2，…，）的復高斯分布。這里假設入射信號在兩個陣列處的方差相同，這樣設定主要是為了簡化模型復雜度。實際應用中，可以根據雷達方程求解出兩個陣列接收到同一目標信號的功率比值，故可以通過補償使兩處的接收信號功率一致。故信號幅值的先驗概率可寫為

式中：x 表示X （＝1，2）中第行列的元素；（·）為沖擊函數。那么接收信號的先驗概率密度可表示為

在SBL框架下，本文假設這些超參數服從伽馬分布：

2 基于SBL的快速DOA估計算法

2.1 利用GAMP算法求解信號幅值的后驗概率

根據模型中的獨立性假設及傳輸信道的湍流效應，不同陣列之間、不同信源之間信號幅值都是相互獨立的，故入射信號的后驗概率可以分解為

式中：、（＝1，2，…，）分別表示、的第行。由于、的先驗分布為復高斯分布，根據共軛先驗的性質，其后驗分布也應為復高斯分布。故只需求出均值和協方差矩陣即可獲得后驗概率密度。

根據各樣本之間的獨立性，有

式中：、 （＝1，2，…，）分別是、中的第個元素。

通過式（11）可以看出，矢量形式的后驗概率密度轉化成了一系列標量后驗概率密度的乘積。這里可以采用文獻［12］提出的GAMP算法計算、的后驗概率密度，具體步驟如算法1所示。

2.2 估計γ

按照標準EM算法的操作，的最大后驗估計為

對于式（19）的求解一般是構建ln（（，，））的緊湊下界（），根據Jensen不等式可得

式中：［·］表示在概率密度（）上求期望。

將式（20）代入式（19）并寫成標量形式，忽略（）中與γ（＝1，2，…，）無關的項，可以得到

式中：目標函數結合式可以寫為

在無信息先驗的情況中，是一個很小的值，忽略并令式（23）等于零，可以得到γ的迭代更新表達式：

經實驗驗證，式（25）比式（24）收斂更快。

2.3 估計目標數量K和噪聲功率 和

一般情況下，目標數量可通過搜索的峰值來估計，但是在樣本數較少、SNR較低的條件下往往會出現多個虛假尖峰。參考文獻［15］的方法，本文采用如下表達式對進行估計：

準確的噪聲功率估計可以加快算法的收斂速度，然而通過類似式（19）的方式求解也會涉及階矩陣求逆的運算，不但計算復雜度較高，而且估計結果比實際值要小。因此，本算法采用下列實證方法來估計噪聲功率：

式中：和可能在每次迭代中會改變，但是在角度估計正確的條件下，上式是對噪聲功率的一致無偏漸近估計。當→∞時，該估計的方差趨近于克拉美羅下限。

2.4 算法初始值的選擇

雖然本算法對初始值不敏感，但是一個好的初始值可以加快算法的收斂速度。在已知SNR的情況下，噪聲功率的初始值可設為

如果缺少SNR的信息，噪聲功率也可初始化為一個較小的值。因為當角度估計不準時，通過式（29）對噪聲功率的估計值會偏高。

信號功率可通過數據協方差矩陣來進行初始化。初始時假設所有方向的信號功率均為，即＝1，其中

2.5 對離網目標的處理

從模型的構建可以看出，角度網格的密度決定著本算法DOA估計的精度。對于離網目標，文獻［19］給出了一個簡單的網格更新策略，但是會增加網格的數量和密度，從而使計算復雜度提高；文獻［10］根據陣列流行矩陣的范德蒙德結構提出了Root-SBL方法，但此方法只適用于均勻線陣，不具有普遍性。因此，本算法在噪聲功率估計完成后，考慮使用梯度方法求解目標所處網格的偏移量來實現對網格位置的動態(tài)調整。

最優(yōu)的網格分布應滿足如下條件：

式中：上標為網格更新步驟中的迭代次數；sign（·）為符號函數；為更新步長。使用單站數據分別計算對同一目標的網格更新，然后換算成一個站的角度，取平均值作為最終的網格更新。

2.6 算法復雜度分析

最后，對本文所提算法的具體實現步驟總結如下。

3 仿真分析

本節(jié)通過數值仿真的方式對所提出的基于GAMP的快速SBL算法（GAMP-SBL）進行有效性驗證，并從空間譜、不同SNR和快拍數下的均方根誤差以及計算時間等4個方面與經典的Lasso算法、BM-SBL算法、M-SBL算法進行對比分析。其中，LASSO算法采用CVX工具箱求解，M-SBL算法默認已知目標數量，故該算法沒有估計這一步驟。

（1）仿真1：算法空間譜分析

圖2為各算法在陣列1上的空間譜。實驗條件為SNR＝-7 dB，SNR＝-10 dB，圖2（a）快拍數為40，圖2（b）快拍數為10。從圖2（a）中可以看出，各算法均能正確識別出兩個目標，本算法與BM-SBL算法的波束寬度基本一致，M-SBL算法波束稍寬。當快拍數進一步降低，如圖2（b）所示，BM-SBL算法副瓣起伏較大，M-SBL算法只估計出了10°方向的目標，LASSO算法性能最好，但這是調整正則化參數后的結果，本算法雖然也出現了峰值較低的副瓣，但不影響目標方位的估計，說明了本算法在小快拍數條件下的穩(wěn)健性。

圖2 算法空間譜分析Fig.2 Analysis of algorithm spatial spectrum

（2）仿真2：算法DOA估計精度分析

圖3為不同SNR下各算法的估計精度。各算法所用快拍數為40，站1的SNR從-10 dB變化到2 dB。從圖3中可以看出，各算法的估計誤差均隨SNR增大而減小，本算法與M-SBL算法的誤差基本一致，均明顯小于BM-SBL和LASSO算法。

圖3 SNR對不同算法估計精度的影響Fig.3 Influence of SNR on the estimation accuracy of different algorithms

圖4為快拍數對算法估計精度的影響。SNR為-10 dB，快拍數為10到150，其中M-SBL算法在10個快拍下不能正確估計，故只統計了快拍數大于20的情況。由圖4可知，隨快拍數的增多，各算法的估計誤差呈下降趨勢，BM-SBL的誤差下降最快，本算法在小快拍時的估計誤差最小，當快拍數較大時，估計誤差與BM-SBL算法接近。

圖4 快拍數對不同算法估計精度的影響Fig.4 Influence of snapshot number on the estimation accuracy of different algorithms

（3）仿真3：算法計算時間分析

平均計算時間可以體現算法的計算復雜度，本實驗中各算法均不采用角度網格更新策略。

圖5為不同SNR下各算法的平均計算時間對比，實驗條件同圖3實驗?？梢钥闯?，除BM-SBL算法外，其他算法的平均計算時間基本不隨SNR變化。另外，本算法的平均計算時間遠小于其他算法，一定程度上說明本算的計算復雜度要小于其他算法。

圖5 不同SNR下各算法的平均計算時間Fig.5 Average calculation time of each method with different SNR

圖6為不同快拍下各算法的平均計算時間對比，實驗條件同圖5實驗。從圖6中可以看出，各算法的平均計算時間均隨快拍數的增加而變大，但是M-SBL、BM-SBL和本算法上升的比較平緩，LASSO算法的增加速度最快，受快拍數的影響最大。本算法平均計算時間仍然遠小于其他算法。對比圖5和圖6可以看出，快拍數對計算時間的影響較大。

圖6 不同快拍下各算法的平均計算時間Fig.6 Average calculation time of each method with different number of snapshots

（4）仿真4：角度空間網格更新分析

為更好地衡量角度更新策略的有效性，實驗設置空間目標源在陣列1上的入射角度為-5°和5.4°，其他實驗條件同圖3實驗，此時無網格更新算法的理論均方根誤差為0.28°。

圖7為不同SNR下算法采用網格更新策略前后的估計誤差。從圖7中可以看出，在SNR較小時，使用更新策略后會導致誤差增大，但是當SNR較大時，更新策略會明顯降低誤差，突破0.28°的理論誤差，說明本算法提出的網格更新策略是有效的。對比圖7與圖3，雖然在小SNR時誤差有所上升，但依然比其他算法要低。

圖7 不同SNR下網格更新前后算法的估計誤差Fig.7 RMSE of the proposed algorithm versus SNR before and after the refinement

圖8為在上述實驗條件下，算法采用網格更新策略前后的平均計算時間。從圖8中可知，算法的平均計算時間對SNR的變化不敏感，增加網格更新策略后，算法計算時間有所提高，但與圖5對比可知仍然顯著低于其他算法。

圖8 不同SNR下網格更新前后算法的平均計算時間Fig.8 Average calculation time of the proposed algorithm versus different SNR before and after the refinement

4 結 論

本文提出一種基于GAMP算法改進的雙基地無源雷達SBL的DOA估計方法。該算法根據雙站雷達接收數據之間的相互獨立性，構建了適用于GAMP算法的信號模型，將SBL中后驗概率密度的計算轉變?yōu)闃肆窟吘壐怕拭芏鹊某朔e，避免了大規(guī)模矩陣求逆的計算。同時針對離網目標，推導了基于梯度下降的角度網格空間更新方法。仿真結果表明，與經典的LASSO、M-SBL以及BM-SBL算法相比，本算法在低SNR和有限快拍數條件下估計誤差低、穩(wěn)健性好，尤其是計算復雜度上有顯著的優(yōu)勢。但本算法主要針對的是相互獨立的目標信號，后續(xù)工作重點將是相關性目標信號下的算法改進及應用研究。