基于一階統計量的子空間旋轉不變解相干算法

李 磊，李國林

（海軍航空工程學院，山東 煙臺 264001）

0 引言

在現代戰爭中，精確制導武器已經成為信息化局部戰爭中物理殺傷的主要手段，這就需要引信探測系統具備高精度、快速目標空間定向能力，實時精確的估計目標的二維波達方向（direction of arrival，DOA）。同時，由于敵方轉發式干擾或多徑效應，探測空間存在大量相干或強相關信源，使得常規空間譜估計算法估計性能下降甚至失效，因此，研究相干信源和低信噪比背景下精確快速的二維解相干算法顯得尤為重要。

近年來，陣列信號處理在DOA 估計領域涌現出了許多子空間類算法［1-3］，但相干或強相關信源會導致陣元接收數據的協方差矩陣虧秩，信號子空間向噪聲子空間擴散，算法估計錯誤。針對相干信源DOA 估計問題，常用處理方式為空間平滑及其改進算法。文獻［4］提出一種基于四階累積量和時間平滑的聯合處理算法，擴展了陣列的自由度，且具有抑制高斯色噪聲的能力。文獻［5］基于共軛循環平穩信號特點，提出了一種虛擬陣列空間平滑算法，提高了DOA 估計的精度和穩健性。另一種常用的解相干算法是矩陣重構類算法，文獻［6］研究了一種基于矩陣重構的波束形成自適應算法，利用Toeplitz特性恢復矩陣的秩實現信源完全解相干。文獻［7］提出了一種改進的Hermitian型矩陣構造方法，有效降低了信源相關度，但上述算法均基于數據的二階統計特性，需要在時域上進行大快拍數據累計及相關運算，計算量大。殷勤業等提出了DOA 矩陣算法［8］，具有計算量小，參數自動配對等優點，同時針對相干信源入射背景提出空域平滑DOA 矩陣算法，但空域平滑處理過程增加了算法的計算復雜度。單次快拍類算法［9］通過一次快拍數據構造解相干矩陣，以一次快拍攜帶信息代替整體數據的統計特性，雖然計算量小，但在信噪比較低時會導致DOA估計精度大幅下降，魯棒性差，無法滿足武器探測系統的估計精度指標。針對以上問題，本文提出基于一階統計量的子空間旋轉不變（ESPRIT）解相干算法。

1 信號接收模型

假設N 個中心波長為λ 且均值不為零的遠場窄帶信號源si（t），i=0，1，…，N 入射到圖1所示雙平行線陣系統。陣列由兩個各向同性且均勻分布的線陣子陣構成，其中子陣X 共2 M 個陣元，陣元編號分別為-（M-1），…，0，…，M ，子陣Y 共2 M-1個陣元，編號分別為-（M-1），…，0，…，M-1。兩子陣間距為D，陣元間距為d。

圖1 陣列模型Fig.1 Array model

圖中，αi、βi 和γi分別表示第i個信源si（t）分別與x、y 和z 軸之間的夾角，則由幾何關系可知cos2αi+cos2βi+cos2γi=1成立，即αi、βi 和γi中只有兩個獨立元素，此處定義αi和βi 分別為入射信源sk（t）的方位角和俯仰角，表述信源的二維波達方向。假設噪聲為零均值加性高斯白噪聲，陣元間噪聲彼此獨立，且與信號不相關。

以子陣X 編號為0的陣元作為參考陣元，則子陣X 第i個陣元的輸出信號可表示為：

式中，

同理，子陣Y 第i個陣元的輸出信號可表示為：

式中，nxi（t）和nyi（t）分別表示子陣X 和子陣Y 第i個陣元的零均值加性高斯白噪聲，滿足

2 基于矩陣重構的解相干算法

常規子空間類譜估計算法都需要進行時域快拍累積（快拍次數大多數通常在數百次以上），然后計算陣列接收數據協方差矩陣，需要進行大量自相關及互相關運算，并且協方差矩陣在相干入射信源背景會導致虧秩現象［10］。本節利用信號一階統計特性構造偽協方差矩陣，無需復雜的相關運算，計算量小，且保證重構的協方差矩陣在相干源背景下不會出現虧秩，從而實現解相干。

將子陣X 進一步分為兩個子陣，前2 M-1個陣元（編號為-（M-1），…，0，…，M-1）組成子陣X1，后2 M-1個陣元（編號-（M-2），…，0，…，M ）組成子陣X2。則子陣的信號輸出矢量分別表示為：

子陣Y 的信號輸出矢量表示為：

式中，A（α）是（2 M-1）×N 維陣列流型矩陣。

子陣X 和子陣Y 接收數據的一階矩定義為：

其中，

表示入射信號的均值。利用子陣X1接收數據的一階 矩mxi，i=- （M -1），…，0，…，M -1 構 造Toeplitz形式的偽協方差矩陣Q1

則易推得

第M 行至第2 M-1行

可以看出，矩陣Rm為對角陣，則Q1的秩只與信源個數相同，而與信源之間的相關性無關，即可實現完全解相干。另一方面，由式（16）可知本文構造的偽協方差矩陣能夠有效濾除噪聲分量，保證了算法在低信噪比下的估計性能。

同理，分別利用子陣X2和子陣Y 接收數據的一階矩構造偽協方差矩陣Q2和Q3

利用上述三個具有旋轉不變關系的偽協方差矩陣，進一步構造3 M×M 維的擴展協方差矩陣Q

對矩陣Q 進行奇異值分解（Q=UΣVH），由奇異值分解的意義可知，若矩陣Q 的秩為N ，則矩陣U 的前N 列組成Q 的列空間的標 準正交基［11］。因此，奇異值分解得到的信號子空間可表示為US=U（：，1：N），而奇異值分解的信號子空間與陣列流型張成的信號子空間相同，即滿足如下關系：

其中，US1=US（1：M，：），US2=US（M +1：2 M，：），US3=US（2 M+1：3 M，：），則信號子空間之間滿足下列關系式

由式（26）及式（27）求最小二乘解，可得

由Φ1和Φ2的形式可知，Φ1只和方位角α 有關，Φ2只和俯仰角β有關。分別對Ψ1和Ψ2進行特征分解，則Ψ1和Ψ2的特征值組成的對角陣與Φ1和Φ2相同，結合式（2）和式（4）即可得到信源入射的二維DOA。

理論上，對Ψ1和Ψ2進行特征分解可得到同樣的特征向量矩陣T，但實際運算中兩次特征分解是獨立進行的，導致兩次特征分解的特征向量矩陣T1和T2的排序可能不同，但同一信號在兩次特征分解中的特征向量強相關，可構造式（28）所示排序矩陣R，依據矩陣R 每列的最大值來調整β 的順序，即可實現二維角度參數匹配。

下面給出本文算法的具體步驟：

步驟1 獲取子陣接收數據X（t）和Y（t）；

步驟2 根據式（13）和式（14）計算子陣X 和子陣Y 接收數據的一階距mxi和myi；

步驟3 根據式（16）、式（20）和式（21）構造偽協方差矩陣Q1、Q2和Q3；

步驟4 根據式（22）構造擴展協方差矩陣Q，對其奇異值分解并取左奇異向量矩陣U 的前N 列為信號子空間US；

步驟5將US按式（23）劃分3個子空間，并按式（26）和式（27）求得Ψ1和Ψ2；

步驟6將Ψ1和Ψ2進行特征值分解，得到Φ1和Φ2。根據式（2）和式（4）的求得αi和βi，并按照式（28）進行配對。

由算法步驟可知，本文算法基于陣列接收數據的一階統計量，步驟1-步驟3對數據預處理過程不存在復乘運算，算法的復乘運算主要集中在后三個實現步驟，約為O（M2N+29 M3）。常規平滑類算法需要進行數據平滑處理，則設定平滑次數為K ，子陣陣元數為M ，快拍次數為L。空域平滑DOA 矩 陣 法［8］需 要 的 復 乘 次 數 約 為O（M2NL +3 M3），由上文可知，基于二階統計量的DOA 算法所需時域快拍累計次數多為數百次以上，則本文算法比空域平滑DOA 矩陣法的計算量低。另外，傳統二維ESPRIT 平滑算法［12］所需復乘次數約為O（8 M2NL+65 M3），則本文算法的計算復雜度較之下降更為明顯。

3 數值仿真

為驗證算法的正確性和有效性，設計如下matlab仿真。考慮DOA 矩陣類算法［8］在分辨力和計算復雜度方面相較其他二維DOA 估計算法的優勢，將本文算法與DOA 矩陣法及空域平滑DOA 矩陣法進行對比，分析算法優勢和局限性。

仿真1：相干信源的二維DOA 估計。假設三個均值非零遠場窄帶相干信源分別從方向（30°，80°），（55°，35°）和（85°，60°）入 射 至 圖1 所 示 陣列，設定陣元參數M=9，信噪比SNR=10dB，陣列噪聲為零均值加性高斯白噪聲。為避免陣列產生測向模糊，設定陣元間距d=λ／2，快拍數均為50次，每種算法分別進行100 次的Monte-Carlo統計仿真，其估計星座圖如圖2所示。

圖2 相干信源下算法DOA 估計星座圖Fig.2 DOA estimation constellation of coherent signal sources

由圖2（a）可知，本文算法可實現相關信源的完全解相干，這是由于算法構造的偽協方差矩陣的秩僅與信源個數有關，保證了算法的二維解相干能力。圖2（b）所示空域平滑DOA 矩陣法通過空域平滑處理也能實現二維DOA 估計，但其解相干能力是以空域平滑所產生的計算量為代價，且其估計結果散布較大，估計性能不及本文算法。分析其原因，一方面，本文算法在構造偽協方差矩陣的過程中有效濾除了噪聲分量，受噪聲影響較小。另一方面，本文算法基于陣列數據的一階統計特性，受快拍次數影響較小。而DOA 矩陣法建立在陣列數據的二階統計特性，對數據長度的要求很高，只有在時域快拍累積非常大時才能達到較高的估計性能。圖2（c）所示常規DOA 矩陣法由于協方差矩陣降秩，DOA 估計已經失效。同時必須看到，本文算法利用信號一階矩特性直接構造偽協方差矩陣，避免了常規子空間類算法的復雜相關運算，運算量小，更適用于武器引信探測系統的應用背景。

仿真2：信噪比對算法估計性能的影響。考慮陣列和信號模型不變，快拍次數不變，使信噪比從0 dB到20dB 區間以1dB 的步長遞增，每個信噪比下分別進行100次獨立Monte-Carlo統計仿真，考察本文算法和空域平滑DOA 矩陣法的均方根誤差（root-mean-square error，RMSE）和分辨成功概率隨SNR 變 化 關 系。DOA 估 計 的RMSE 定 義［11］為：

式中，（αk，βk）為實際值，（^αk，^βk）為估計值。

圖3（a）為RMSE 隨信噪比變化曲線。隨著信噪比的增大，兩種算法的DOA 估計的RMSE 都逐漸減小，且本文算法RMSE 始終比空域平滑DOA矩陣法小，尤其當信噪比較低時。分析其原因，當陣列噪聲為零均值高斯白噪聲時，本文構造的偽協方差矩陣能夠有效濾除噪聲分量。雖然空域平滑DOA 矩陣法也進行了相應的去噪處理，但由于快拍次數較小，其去噪能力受到一定限制。當SNR=0 dB時，空域平滑DOA 矩陣法的三個信源RMSE 分別約為1.8、1.9 和2.1，而本文算法三個信源RMSE相對較小，分別約為0.6、1和1.3，估計精度更高。

圖3（b）為分辨成功概率隨信噪比變化曲線。設定信號DOA 估計的RMSE小于1時視為分辨成功。可以看出，兩種算法對每個信源的DOA 估計分辨概率隨信噪比的增大而提高，且本文算法優于空域平滑DOA 矩陣法，與之前的理論分析保持一致。

圖3 估計性能隨信噪比變化曲線Fig.3 Performance curves vs.SNR

仿真3：快拍次數對估計性能的影響。考慮陣列和信號模型不變，設定信噪比SNR=10dB，快拍數從50到1000以50為步長遞增，每個快拍數下分別進行100 次獨立Monte-Carlo統計仿真，得到本文算法和空域平滑DOA 矩陣法的RMSE隨快拍次數的變化關系如圖4所示。顯然，隨著快拍數的增加，兩種算法的3個信源DOA 估計的RMSE 都逐漸減小，且本文算法DOA 估計的RMSE始終比相同快拍數時的空域平滑DOA 矩陣法小，與上文分析一致。在快拍數為50次時，本文算法將三個相干信源的估計誤差均控制在1以內，更符合少快拍數下引信探測系統對目標二維方位估計精度的需求。

圖4 RMSE隨快拍次數變化曲線Fig.4 RMSE curves vs.snapshots

4 結論

本文提出基于一階統計量的二維ESPRIT 解相干算法。算法利用陣元接收數據的一階統計量構造三個Toeplitz形式的偽協方差矩陣，實現了相干信源的解相干。然后通過構造擴展協方差矩陣并對其進行一次奇異值分解即實現了相干信源的二維DOA 估計。仿真結果表明：該算法在低信噪比和少快拍數下具有比空域平滑DOA 矩陣法更高的估計性能，且避免了傳統算法對大量時域快拍累計的協方差矩陣的處理過程，計算復雜度低，實時性高。

［1］蔣柏峰，呂曉德.一種基于導向矢量變換的DOA 估計預處理方法［J］.電子與信息學報，2012，34（7）：1552-1557.

［2］宋愛民，李堰，劉劍，等.非圓信號多級維納濾波DOA估計求根算法［J］.電子科技大學學報，2013，42（1）：53-57.

［3］Mendoza Montoya F，Covarrubias R D H，Lopez-Miranda C A.DOA estimation in mobile communication system usingsubspace tracking methods［J］.IEEE Latin America Transactions，2008，6（2）：123-129.

［4］景小榮，隋偉偉，周圍.基于四階累積量和時間平滑的相干信號DOA 估計［J］.系統工程與電子技術，2012，34（4）：789-795.

［5］張聰，胡謀法，盧煥章.基于虛擬陣列空間平滑的相干信號DOA 估計［J］.電子學報，2010，38（4）：929-933.

［6］Gao Jian-yong，Gao Yong.The adaptive beamforming algorithm for coherent sources non-uniform linear array based on improved Toeplitz approximation［J］.Journal of Astronautics，2007，28（6）：1715-1718.

［7］韓勇，喬曉林，金銘，等.空間相關信號源高分辨處理的Toep-MUSIC 改進算法［J］.系統工程與電子技術，2009，31（7）：1544-1547.

［8］殷勤業，鄒理和，Newwcomb W R.一種高分辨率二維信號參量估計方法-波達方向矩陣法［J］.通信學報，1991，12（4）：1-7.

［9］洪升，萬顯榮，易建新，等.基于單次快拍的雙基地MIMO雷達多目標角度估計方法［J］.電子與信息學報，2013，35（5）：1149-1156.

［10］張賢達.矩陣分析與應用［M］.北京：清華大學出版社，2009.

［11］王永良，陳輝，彭應寧，等.空間譜估計理論與算法［M］.北京：清華大學出版社，2009.

［12］董軼，吳云韜，廖桂生.一種二維到達方向估計的ESPRIT 新方法［J］.西安電子科大學報，2003，30（5）：569-575.