基于DML－ESPRIT算法的多維參數(shù)快速聯(lián)合估計

羅 爭，張 旻，劉 圓

(1.中國人民解放軍61123部隊，北京 101309;2.合肥電子工程學院，安徽合肥 230037)

二維DOA(Direction－of－Arrival)估計一直是通信、雷達、導航等領域研究的熱點問題［1－4］。在現(xiàn)代密集復雜的信號環(huán)境下，電子偵察中的信號處理任務越來越難，對無線電測向提出了高精度、高分辨率、實時快速等要求。特別是對遠場目標信號的方位角、俯仰角、功率與頻率［5－6］等參數(shù)的快速聯(lián)合估計，是電子戰(zhàn)與雷達領域的一個新興研究方向，目前在國內外尚處于初步研究階段。

眾所周知，貝葉斯方法是基于統(tǒng)計理論的一種經典方法，適用于參數(shù)估計問題。最大似然估計方法就是貝葉斯方法的一種特例，是在已知白噪聲情況下的貝葉斯最優(yōu)估計［7－8］。在最大似然估計算法中，觀測所得信號的似然函數(shù)被定義為含有未知參數(shù)的條件概率密度函數(shù)，目的是選定未知的參數(shù)以使得該似然函數(shù)盡可能大，通過最大化似然函數(shù)求出的解都被認為是未知參數(shù)的一個估計［9］。Ziskind L和 Wax M于1988年將最大似然參數(shù)估計方法應用于DOA估計問題［10］。目前，最大似然算法在空間譜估計領域中已有大量高價值的研究成果，與 MUSIC(Multiple Signal Classification，MUSIC)［7］等子空間分解類 DOA 估計算法相比，最大似然算法的估計精度高，具有良好的魯棒性和穩(wěn)定性，尤其是低信噪比、小快拍數(shù)據(jù)情況下，最大似然算法比子空間分解類算法性能更好。最大似然算法的思想簡單、估計性能優(yōu)越，但其算法的實現(xiàn)過程復雜，求解涉及到一個多維的非線性優(yōu)化問題，計算量大。而二維ESPRIT(Estimation of Signal Parameters via Rational Invariance Techniques，ESPRIT)算法［11］在進行二維DOA估計時無需進行二維譜峰搜索，因此大幅降低了計算量，能滿足強實時性的要求，具有廣泛的應用前景。

針對以上問題，論文在雙L陣列的基礎上，探索了一種二維DOA－功率－頻率的快速聯(lián)合估計算法。通過引入空間錐角表示方法，將方位角與俯仰角的估計轉換為了獨立的空間錐角估計，并利用確定性最大似然方法實現(xiàn)了多輻射源的二維DOA－功率－頻率的聯(lián)合估計;圍繞最大似然估計模型計算復雜的問題，通過對陣列進行擴展，進一步利用TLS－ESPRIT算法實現(xiàn)了模型的快速求解，從而達到二維DOA－功率－頻率快速聯(lián)合估計的目的。

1 陣列結構與信號模型

如圖1所示的3M+1元的雙L型陣列，其中子陣X1、X2位于 x 軸上，子陣 Y1、Y2位于 y 軸上，子陣 Z1、Z2位于z軸上，x軸、y軸與z軸互成90°，為便于描述，本文用X、Y、Z分別表示x軸、y軸、z軸上的子陣，子陣 X、Y、Z的結構相同，均為等距線陣，且陣元數(shù)均為M+1。

圖1 雙L陣列的結構示意圖

假設p個遠場信號入射上述雙L型陣列，第k(k=1，2，…，p)個信號源的二維 DOA 信息為(θk，φk)，其中θk、φk分別表示其方位角和俯仰角，則在某一時刻t，子陣 X、Y、Z接收的快拍數(shù)據(jù)矢量 x(t)、y(t)、z(t)可分別表示如下

式(1)中，Ax(θ，φ，f)、Ay(θ，φ，f)、Az(θ，φ，f)分別表示子陣 X、Y、Z 的陣列流形矩陣;nx(t)、ny(t)、nz(t)分別表示子陣X、Y、Z的陣列噪聲數(shù)據(jù)矢量;s(t)=［s1(t)，s2(t)，…，sp(t)］表示 t時刻的 p×1維的空間信號矢量。

且空間信號矢量中的元素sk(t)可表示為

式(2)中，ui(t)、Ф(t)和w0分別表示目標信號k的幅度、相位和頻率信息。

2 二維DOA－功率－頻率快速聯(lián)合估計

2.1 二維DOA的空間錐角降維表示

受降維思想的啟發(fā)，提出三維空間錐角概念:在圖2所示的空間圖形中，以三維坐標系的原點0為頂點，以坐標軸為中心線可獲得各個子陣的半圓錐曲面，將半圓錐面上的母線與坐標軸正方向的夾角稱為空間錐角，并分別用α、β、γ表示X軸、Y軸、Z軸上的空間錐角，且{α，β，γ}∈［0，180°］。顯然，雙 L 陣中 X 軸、Y軸上的半圓錐曲面的交線就是輻射源的入射角。

圖2 二維DOA的空間錐角降維表示

由幾何關系可知，空間錐角 α、β、γ與信號二維DOA信息θ、φ存在以下的轉換關系

對于多輻射源的二維DOA估計問題，如果能夠分別獲取與雙L陣列的入射信號配對的空間錐角，就能準確獲取各輻射源的二維DOA信息，而且空間錐角定義的測向方法，每一維角度均是獨立估計，不是聯(lián)合估計，也就不會產生聯(lián)合估計中的因為仰角估計精度不高而導致方位角估計失敗的問題。

2.2 空間錐角－功率－頻率估計模型

根據(jù)上文所描述地陣列信號結構以及噪聲模型，在此以子陣X為例進行說明，其他子陣同理類推。對于未知的確定性目標信號源，觀測數(shù)據(jù)的一階矩、二階矩參數(shù)滿足如下條件

上式中，式(4)表示確定性最大似然(DML)準則的均值;式(5)表示其方差，則觀測矢量N次快拍數(shù)據(jù)的聯(lián)合概率密度函數(shù)可表示為

對式兩邊同時取負對數(shù)，可得

由式可知，聯(lián)合概率密度函數(shù)F(x1，x2，…，xN)是一個關于未知參量α、σ2、si的函數(shù)，要得到參數(shù)的最大似然估計，即求得一組參量使得估計準則式(8)最小，故參量α、σ2、S的最大似然估計公式分別為

式中，A*表示陣列流型A的偽逆;PA表示陣列流型A的投影矩陣;P⊥A表示PA的正交矩陣;R表示陣列快拍數(shù)據(jù)協(xié)方差矩陣R的估計值，即協(xié)方差矩陣，inv{·}為求跡算子;arg{·}為取相位算子。

對確定性最大似然估計原則進行擴展，根據(jù)式可得信號功率與頻率估計

DML算法的角度估計過程是一個復雜的多維搜索問題，計算量隨著目標個數(shù)的增加呈指數(shù)增長。因此，文獻［12～13］提出利用交替投影算法(AP)來實現(xiàn)角度估計，該算法將復雜的多維網(wǎng)格搜索轉化為簡單的多個一維搜索，在一定程度上降低了DML算法的計算量，但當信源數(shù)較多時，算法的收斂速度相當慢。

因此，文中提出利用快速求解算法TLS－ESPRIT進行空間錐角估計，優(yōu)點在于:(1)ESPRIT算法無需譜峰搜索，算法估計時效性強，算法運算時間不受信源數(shù)影響，適用于多信源的二維DOA估計。(2)LSESPRIT、TLS－ESPRIT、TAM 及實值空間的 ESPRIT算法的估計性能接近，但在低信噪比情況下的TLSESPRIT算法估計性能最優(yōu)，可滿足復雜環(huán)境下的定位探測需求。

2.3 多維參數(shù)聯(lián)合估計模型的求解

假設子陣X1、X2的接收數(shù)據(jù)分別為x1、x2，并構造如下數(shù)據(jù)矩陣

對上式的協(xié)方差矩陣R=E［xxH］進行特征值分解可得

式中，∑S為大特征值組成的對角陣;∑N為小特征值組成的對角陣;US為大特征值對應的特征矢量張成的子空間也即信號子空間;UN為小特征值對應的特征矢量張成的子空間也即噪聲子空間，可知

顯然，子陣X1的大特征值對應的特征矢量張成的子空間US1、子陣X2的大特征值對應的特征矢量張成的子空間US2與陣列流型A張成的子空間三者相等，即

此時，有且只有一個非奇異矩陣T，使得

總體最小二乘的基本思想就是在約束條件下，同時使得校正項ΔUS1，ΔUS2盡可能小。因此，TLS的解等價于

其中與特征值對角矩陣∑對應的特征向量矩陣E又可寫成

由文獻［14］可知，構造一個矩陣Ψ就可以得到來波方位信息，構造方法可通過式得到

然后對ΨTLS進行特征值便可得到目標的空間錐角信息［9］。

2.4 空間錐角配對與角度轉換

利用TLS－ESPRIT估計算法獲得的只是獨立的空間錐角信息，最終需要配對算法對各個子陣的空間錐角參數(shù)進行配對組合，從而獲得目標信號源的二維DOA 信息。從空間幾何知識可知，(cosαk，cosβk，cosγk)為一個信號源來波方向向量的空間余弦，其間存在下列數(shù)學關系

定義三維空間錐角組合(αi，βj，γk)的數(shù)據(jù)關聯(lián)殘差

式中，{i，j，k}∈［1，ρ］。如果空間錐角組合(αi，βj，γk)為同一目標的來波方向時，應滿足

此外，還可利用估計所得的信號功率、頻率信息對空間錐角進行配對，故空間錐角的配對問題在本文中可得到較好的解決。

根據(jù)式(26)對 αi，βj，γk空間錐角進行配對之后，需要將空間錐角參數(shù)轉換為信號源的二維DOA數(shù)據(jù)。依據(jù)式(3)，令根據(jù)幾何關系可知信號源的方位角與俯仰角為

以上內容以子陣X為例進行了推導，對子陣Y、Z同樣適用。

2.5 DML－ESPRIT算法

綜上可知，利用最大似然估計理論可得到α、σ2、S等參數(shù)的聯(lián)合估計模型，具體見式(9)～式(11)，然后利用式(12)和式(13)可實現(xiàn)輻射源功率、頻率的估計，但現(xiàn)有的DML求解方法計算量過大，而利用TLSESPRIT算法恰好解決了該問題，文中充分利用了最大似然估計理論與ESPRIT算法的各自獨特優(yōu)點，實現(xiàn)了多目標二維DOA－功率－頻率的快速聯(lián)合估計，故文中稱該算法為DML－ESPRIT算法。

表1歸納了基于DML－ESPRIT算法的二維DOA－功率－頻率快速聯(lián)合估計步驟。

表1 DML－ESPRIT算法

3 仿真結果

為驗證DML－ESPRIT算法的性能，特進行如下實驗。實驗仿真條件:考慮如圖1所示的陣列結構，各子陣為6元均勻線陣，陣元間隔d=λmin/2，噪聲模型為加性高斯白噪聲，采樣快拍數(shù)為512。實驗環(huán)境:Intel Core(TM)2 Duo CPU 2.19 GHz;2.0 GB內存;Matlab 7.8.0仿真平臺;Windows XP SP2操作系統(tǒng)。

3.1 二維DOA－功率－頻率聯(lián)合估計實驗

假設3個不同功率的遠場獨立信號入射到如前所述雙L型陣列上，各信號的方位角、俯仰角、中心頻率與功率參數(shù)組合(θk，φk，Pk，fk)分別為:(29°，332°，0.35，20 MHz)、(44°，113°，1，25 MHz)、(78°，52°，0.5，30 MHz);利用公式組(3)可知，各信源的方位及俯仰角的空間錐角(α，β，γ)對應表示為:(39.44°，114.24°，61°)、(106.32°，48.54°，46°)、(82.65°，80.57°，12°)。

為了對算法的估計性能進行比較，分別使用APDML算法、AP－SSF算法、DML－ESPRIT算法與MUSIC算法對信號源的參數(shù)進行估計。圖3顯示了在SNR=10 dB情況下，DML－ESPRIT算法1 000次Monte－Carlo實驗的α、β、γ及功率估計結果直方圖;圖4為4種算法二維DOA估計與真實二維DOA星座圖;圖5為基于DML－ESPRIT算法頻率與功率聯(lián)合估計隨信噪比變化曲線，信噪比范圍為－10～15 dB，步進為3 dB。

圖3 空間錐角與功率聯(lián)合估計

圖4 二維DOA估計星座圖

圖5 DML－ESPRIT算法的頻率與功率聯(lián)合估計結果隨SNR的變化曲線

如圖3所示，DML－ESPRIT算法能準確地實現(xiàn)各個信號空間錐角和功率的聯(lián)合估計，估計結果接近真實值;圖4的仿真結果表明，根據(jù)角度配對方法正確實現(xiàn)了空間錐角α、β的精確配對，轉換后的二維DOA估計結果基本與真實值重合;另外，在低信噪比條件下，SNR=－10 dB時，DML－ESPRIT算法仍能精確估計出各個目標信號的功率、頻率參數(shù)，如圖5所示，SNR=5 dB時，DML－ESPRIT算法的二維DOA－功率－頻率聯(lián)合估計結果，真實值與估計值具體對比情況如表2所示，可見本文所提算法各個參數(shù)的估計結果均接近3個目標的真實值。

表2 DML－ESPRIT算法的二維DOA－功率－頻率聯(lián)合估計結果

3.2 算法估計精度實驗

采用上述天線陣列結構，假設在(0°～360°，0°～90°)范圍內隨機產生的3個信號源，歸一化信號功率在0.1～1之間內隨機產生，頻率分別為10 MHz、15 MHz、20 MHz，信噪比范圍為－10 ～15 dB，步進3 dB，每個信噪比下進行1 000次Monte－Carlo實驗。圖6為不同信噪比下 MUSIC、AP－DML、AP－SSF與 DML－ESPRIT算法的角度估計誤差隨信噪比變化曲線，圖7為4種算法的功率估計誤差隨信噪比變化曲線;圖8為4種算法的頻率估計結果隨信噪比變化曲線。為了量化算法估計精度，特將二維DOA與功率的估計誤差定義如下

式中，L表示Monte－Carlo實驗次數(shù)。

圖6 二維DOA估計RMSE隨SNR變化曲線

圖7 功率估計RMSE隨SNR變化曲線

圖8 頻率估計結果隨SNR變化曲線

由圖6～圖8顯示的仿真結果可知，DML－ESPRIT算法能實現(xiàn)對目標信號的二維DOA－功率－頻率的高精度估計，且二維DOA估計精度高于其他3種算法。

3.3 算法時效性驗證實驗

實驗目的是對實驗2中所涉及4種算法的時效性進行比較。天線陣列結構及信號模型與實驗1相同，假設兩個來波信號的參數(shù)(θk，φk，fk，Pk)分別為:(30°，120°，15，0.5)、(60°，240°，20，1)，在 SNR=10 dB條件下，進行1 000次Monte－Carlo實驗，表3所示MUSIC、AP－DML、AP－SSF與 DML－ESPRIT算法每次運行的平均耗時。

表3 算法運行時間比較

表3統(tǒng)計結果表明DML－ESPRIT算法的運行時間遠小于其他3種算法，平均耗時為MUSIC算法的1/26、AP－DML算法的1/36。另外，由于 MUSIC、APSSF、AP－DML算法需進行譜峰搜索，因此計算量會隨著空間分辨率的增大成指數(shù)增長，而DML－ESPRIT因無需譜峰搜索，故不存在該問題。

以上3個實驗的仿真結果證實:本文提出的DMLESPRIT算法不但實現(xiàn)了多輻射源的二維DOA－功率－頻率聯(lián)合估計，且在保持較高估計精度的前提下，算法時效性遠優(yōu)于其他算法，平均耗時約為35 ms，可滿足實際工程應用。

4 結束語

文中提出的DML－ESPRIT算法主要依托雙L陣列的空間特性，將方位角與俯仰角的估計問題轉換為空間錐角估計，并利用最大似然估計理論推導了輻射源的空間錐角、功率和頻率的聯(lián)合估計數(shù)學模型;然后通過應用TLS－ESPRIT算法實現(xiàn)了參數(shù)的快速估計。目標輻射源多維參數(shù)快速聯(lián)合估計的實現(xiàn)為多站多目標測向定位［15］、方位數(shù)據(jù)關聯(lián)等問題［16－17］的解決提供了新思路。

［1］Mashud H，Kaushik M.Direction－of－arrival estimation using a mixed L2.0 norm approximation ［J］.IEEE Transactions on Signal Processing，2010，58(9):4646 －4655.

［2］Chen F J，Sam K，Chaiwah K.ESPRIT － like two － dimensional DOA estimation for coherent signals［J］.IEEE Transactions on Aerospace and Electronic Systems，2010，46(3):1477－1484.

［3］Liang J L，Liu D.Joint elevation and azimuth direction finding using L － shaped array［J］.IEEE Transactions on Antennas and Propagation，2010，58(6):2136 －2141.

［4］Shannon D B，Tszping C，Karl G.Robust DOA estimation:the reiterative super resolution(RISR)algorithm ［J］.IEEE Transactions on Aerospace and Electronic Systems，2011，47(1):332－346.

［5］杜剛，王永良，張永順，等.空間相干信號的頻率和二維到達角估計算法［J］.系統(tǒng)工程與電子技術，2008，30(6):1050－1053.

［6］劉聰鋒，廖桂生.寬帶接收機的窄帶信號頻率和二維角度估計新方法［J］.電子學報，2009，27(3):523－528.

［7］Petre S，Arye N.MUSIC，maximum likelihood，and cramerrao bound［J］.IEEE Transactions on Acoustic，Speech and Signal Processing，1989，37(5):720 －740.

［8］Petre S，Arye N.MUSIC，maximum likelihood，and cramerrao bound:further results and comparisons［J］.IEEE Transactions on Acoustic，Speech，and Signal Processing，1990，38(12):2140－2150.

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

［10］Liskind I，Wax M.Maximum likelihood localization of multiple sources by alternating projection［J］.IEEE Transactions on Acoustic，Speech，and Signal Processing，1988，36(10):1553－1559.

［11］Roy R K T.ESPRIT－estimation of signal parameters via rotational invariance techniques［J］.IEEE Transactions on A-coustic，Speech，and Signal Processing，1986，37(7):984－995.

［12］Wong K M，Reilly J P，Wu Q，et al.Estimation of the directions of arrival of signals in unknown correlated noise，Part I:the MAP approach and its implementation［J］.IEEE Transactions on Signal Processing，1992，40(8):2007 －2017.

［13］Wong K M，Reilly J P，Wu Q，et al.Estimation of the directions of arrival of signals in unknown correlated noise，Part II:Asymptotic behavior and performance of the MAP approach［J］.IEEE Transactions on Signal Processing，1992，40(8):2018－2028.

［14］楊磊，趙擁軍，王志剛.基于功率和相位聯(lián)合估計TLSESPRIT算法的極化干涉SAR數(shù)據(jù)分析［J］.測繪學報，2007，36(2):163 －168.

［15］周益明.一種基于RSS的環(huán)境自適應目標定位算法［J］.通信對抗，2011(2):9－11.

［16］Chen SL，Xu Y B.A new joint possibility data association algorithm avoiding track coalescence［J］.International Journal of Intelligent Systems and Applications，2011，3(2):45 －51.

［17］Oh S，Russell S，Sastry S.Markov chain monte carlo data association for multi－ target tracking［J］.IEEE Transactions on Automatic Control，2009，54(3):481 －497.