李 磊，李國林

（海軍航空工程學(xué)院，山東 煙臺(tái) 264001）

0 引言

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

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

1 信號(hào)接收模型

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

圖1 陣列模型Fig.1 Array model

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

以子陣X 編號(hào)為0的陣元作為參考陣元，則子陣X 第i個(gè)陣元的輸出信號(hào)可表示為：

式中，

同理，子陣Y 第i個(gè)陣元的輸出信號(hào)可表示為：

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

2 基于矩陣重構(gòu)的解相干算法

常規(guī)子空間類譜估計(jì)算法都需要進(jìn)行時(shí)域快拍累積（快拍次數(shù)大多數(shù)通常在數(shù)百次以上），然后計(jì)算陣列接收數(shù)據(jù)協(xié)方差矩陣，需要進(jìn)行大量自相關(guān)及互相關(guān)運(yùn)算，并且協(xié)方差矩陣在相干入射信源背景會(huì)導(dǎo)致虧秩現(xiàn)象［10］。本節(jié)利用信號(hào)一階統(tǒng)計(jì)特性構(gòu)造偽協(xié)方差矩陣，無需復(fù)雜的相關(guān)運(yùn)算，計(jì)算量小，且保證重構(gòu)的協(xié)方差矩陣在相干源背景下不會(huì)出現(xiàn)虧秩，從而實(shí)現(xiàn)解相干。

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

子陣Y 的信號(hào)輸出矢量表示為：

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

子陣X 和子陣Y 接收數(shù)據(jù)的一階矩定義為：

其中，

表示入射信號(hào)的均值。利用子陣X1接收數(shù)據(jù)的一階 矩mxi，i=- （M -1），…，0，…，M -1 構(gòu) 造Toeplitz形式的偽協(xié)方差矩陣Q1

則易推得

第M 行至第2 M-1行

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

同理，分別利用子陣X2和子陣Y 接收數(shù)據(jù)的一階矩構(gòu)造偽協(xié)方差矩陣Q2和Q3

利用上述三個(gè)具有旋轉(zhuǎn)不變關(guān)系的偽協(xié)方差矩陣，進(jìn)一步構(gòu)造3 M×M 維的擴(kuò)展協(xié)方差矩陣Q

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

其中，US1=US（1：M，：），US2=US（M +1：2 M，：），US3=US（2 M+1：3 M，：），則信號(hào)子空間之間滿足下列關(guān)系式

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

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

理論上，對(duì)Ψ1和Ψ2進(jìn)行特征分解可得到同樣的特征向量矩陣T，但實(shí)際運(yùn)算中兩次特征分解是獨(dú)立進(jìn)行的，導(dǎo)致兩次特征分解的特征向量矩陣T1和T2的排序可能不同，但同一信號(hào)在兩次特征分解中的特征向量強(qiáng)相關(guān)，可構(gòu)造式（28）所示排序矩陣R，依據(jù)矩陣R 每列的最大值來調(diào)整β 的順序，即可實(shí)現(xiàn)二維角度參數(shù)匹配。

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

步驟1 獲取子陣接收數(shù)據(jù)X（t）和Y（t）；

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

步驟3 根據(jù)式（16）、式（20）和式（21）構(gòu)造偽協(xié)方差矩陣Q1、Q2和Q3；

步驟4 根據(jù)式（22）構(gòu)造擴(kuò)展協(xié)方差矩陣Q，對(duì)其奇異值分解并取左奇異向量矩陣U 的前N 列為信號(hào)子空間US；

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

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

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

3 數(shù)值仿真

為驗(yàn)證算法的正確性和有效性，設(shè)計(jì)如下matlab仿真?？紤]DOA 矩陣類算法［8］在分辨力和計(jì)算復(fù)雜度方面相較其他二維DOA 估計(jì)算法的優(yōu)勢(shì)，將本文算法與DOA 矩陣法及空域平滑DOA 矩陣法進(jìn)行對(duì)比，分析算法優(yōu)勢(shì)和局限性。

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

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

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

仿真2：信噪比對(duì)算法估計(jì)性能的影響?？紤]陣列和信號(hào)模型不變，快拍次數(shù)不變，使信噪比從0 dB到20dB 區(qū)間以1dB 的步長遞增，每個(gè)信噪比下分別進(jìn)行100次獨(dú)立Monte-Carlo統(tǒng)計(jì)仿真，考察本文算法和空域平滑DOA 矩陣法的均方根誤差（root-mean-square error，RMSE）和分辨成功概率隨SNR 變 化 關(guān) 系。DOA 估 計(jì) 的RMSE 定 義［11］為：

式中，（αk，βk）為實(shí)際值，（^αk，^βk）為估計(jì)值。

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

圖3（b）為分辨成功概率隨信噪比變化曲線。設(shè)定信號(hào)DOA 估計(jì)的RMSE小于1時(shí)視為分辨成功?？梢钥闯?，兩種算法對(duì)每個(gè)信源的DOA 估計(jì)分辨概率隨信噪比的增大而提高，且本文算法優(yōu)于空域平滑DOA 矩陣法，與之前的理論分析保持一致。

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

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

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

4 結(jié)論

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

［1］蔣柏峰，呂曉德.一種基于導(dǎo)向矢量變換的DOA 估計(jì)預(yù)處理方法［J］.電子與信息學(xué)報(bào)，2012，34（7）：1552-1557.

［2］宋愛民，李堰，劉劍，等.非圓信號(hào)多級(jí)維納濾波DOA估計(jì)求根算法［J］.電子科技大學(xué)學(xué)報(bào)，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］景小榮，隋偉偉，周圍.基于四階累積量和時(shí)間平滑的相干信號(hào)DOA 估計(jì)［J］.系統(tǒng)工程與電子技術(shù)，2012，34（4）：789-795.

［5］張聰，胡謀法，盧煥章.基于虛擬陣列空間平滑的相干信號(hào)DOA 估計(jì)［J］.電子學(xué)報(bào)，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］韓勇，喬曉林，金銘，等.空間相關(guān)信號(hào)源高分辨處理的Toep-MUSIC 改進(jìn)算法［J］.系統(tǒng)工程與電子技術(shù)，2009，31（7）：1544-1547.

［8］殷勤業(yè)，鄒理和，Newwcomb W R.一種高分辨率二維信號(hào)參量估計(jì)方法-波達(dá)方向矩陣法［J］.通信學(xué)報(bào)，1991，12（4）：1-7.

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

［10］張賢達(dá).矩陣分析與應(yīng)用［M］.北京：清華大學(xué)出版社，2009.

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

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