羅 爭，張 旻，劉 圓

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

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

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

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

1 陣列結(jié)構(gòu)與信號(hào)模型

如圖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的結(jié)構(gòu)相同，均為等距線陣，且陣元數(shù)均為M+1。

圖1 雙L陣列的結(jié)構(gòu)示意圖

假設(shè)p個(gè)遠(yuǎn)場信號(hào)入射上述雙L型陣列，第k(k=1，2，…，p)個(gè)信號(hào)源的二維 DOA 信息為(θk，φk)，其中θk、φk分別表示其方位角和俯仰角，則在某一時(shí)刻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時(shí)刻的 p×1維的空間信號(hào)矢量。

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

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

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

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

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

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

由幾何關(guān)系可知，空間錐角 α、β、γ與信號(hào)二維DOA信息θ、φ存在以下的轉(zhuǎn)換關(guān)系

對(duì)于多輻射源的二維DOA估計(jì)問題，如果能夠分別獲取與雙L陣列的入射信號(hào)配對(duì)的空間錐角，就能準(zhǔn)確獲取各輻射源的二維DOA信息，而且空間錐角定義的測向方法，每一維角度均是獨(dú)立估計(jì)，不是聯(lián)合估計(jì)，也就不會(huì)產(chǎn)生聯(lián)合估計(jì)中的因?yàn)檠鼋枪烙?jì)精度不高而導(dǎo)致方位角估計(jì)失敗的問題。

2.2 空間錐角－功率－頻率估計(jì)模型

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

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

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

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

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

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

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

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

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

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

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

式中，∑S為大特征值組成的對(duì)角陣;∑N為小特征值組成的對(duì)角陣;US為大特征值對(duì)應(yīng)的特征矢量張成的子空間也即信號(hào)子空間;UN為小特征值對(duì)應(yīng)的特征矢量張成的子空間也即噪聲子空間，可知

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

此時(shí)，有且只有一個(gè)非奇異矩陣T，使得

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

其中與特征值對(duì)角矩陣∑對(duì)應(yīng)的特征向量矩陣E又可寫成

由文獻(xiàn)［14］可知，構(gòu)造一個(gè)矩陣Ψ就可以得到來波方位信息，構(gòu)造方法可通過式得到

然后對(duì)ΨTLS進(jìn)行特征值便可得到目標(biāo)的空間錐角信息［9］。

2.4 空間錐角配對(duì)與角度轉(zhuǎn)換

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

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

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

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

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

以上內(nèi)容以子陣X為例進(jìn)行了推導(dǎo)，對(duì)子陣Y、Z同樣適用。

2.5 DML－ESPRIT算法

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

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

表1 DML－ESPRIT算法

3 仿真結(jié)果

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

3.1 二維DOA－功率－頻率聯(lián)合估計(jì)實(shí)驗(yàn)

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

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

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

圖4 二維DOA估計(jì)星座圖

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

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

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

3.2 算法估計(jì)精度實(shí)驗(yàn)

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

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

圖6 二維DOA估計(jì)RMSE隨SNR變化曲線

圖7 功率估計(jì)RMSE隨SNR變化曲線

圖8 頻率估計(jì)結(jié)果隨SNR變化曲線

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

3.3 算法時(shí)效性驗(yàn)證實(shí)驗(yàn)

實(shí)驗(yàn)?zāi)康氖菍?duì)實(shí)驗(yàn)2中所涉及4種算法的時(shí)效性進(jìn)行比較。天線陣列結(jié)構(gòu)及信號(hào)模型與實(shí)驗(yàn)1相同，假設(shè)兩個(gè)來波信號(hào)的參數(shù)(θk，φk，fk，Pk)分別為:(30°，120°，15，0.5)、(60°，240°，20，1)，在 SNR=10 dB條件下，進(jìn)行1 000次Monte－Carlo實(shí)驗(yàn)，表3所示MUSIC、AP－DML、AP－SSF與 DML－ESPRIT算法每次運(yùn)行的平均耗時(shí)。

表3 算法運(yùn)行時(shí)間比較

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

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

4 結(jié)束語

文中提出的DML－ESPRIT算法主要依托雙L陣列的空間特性，將方位角與俯仰角的估計(jì)問題轉(zhuǎn)換為空間錐角估計(jì)，并利用最大似然估計(jì)理論推導(dǎo)了輻射源的空間錐角、功率和頻率的聯(lián)合估計(jì)數(shù)學(xué)模型;然后通過應(yīng)用TLS－ESPRIT算法實(shí)現(xiàn)了參數(shù)的快速估計(jì)。目標(biāo)輻射源多維參數(shù)快速聯(lián)合估計(jì)的實(shí)現(xiàn)為多站多目標(biāo)測向定位［15］、方位數(shù)據(jù)關(guān)聯(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］杜剛，王永良，張永順，等.空間相干信號(hào)的頻率和二維到達(dá)角估計(jì)算法［J］.系統(tǒng)工程與電子技術(shù)，2008，30(6):1050－1053.

［6］劉聰鋒，廖桂生.寬帶接收機(jī)的窄帶信號(hào)頻率和二維角度估計(jì)新方法［J］.電子學(xué)報(bào)，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］王永良，陳輝，彭應(yīng)寧，等.空間譜估計(jì)理論與算法［M］.北京:清華大學(xué)出版社，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)合估計(jì)TLSESPRIT算法的極化干涉SAR數(shù)據(jù)分析［J］.測繪學(xué)報(bào)，2007，36(2):163 －168.

［15］周益明.一種基于RSS的環(huán)境自適應(yīng)目標(biāo)定位算法［J］.通信對(duì)抗，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.