張志成，林 君，石要武，孫曉東

（1.吉林大學(xué) 通信工程學(xué)院，長(zhǎng)春 130012；2.吉林大學(xué) 儀器科學(xué)與電氣工程學(xué)院，長(zhǎng)春 130061）

0 引 言

近年來(lái)，波達(dá)方向（DOA）和多普勒頻率聯(lián)合估計(jì)技術(shù)在多個(gè)領(lǐng)域引起了廣泛的關(guān)注和研究，但估計(jì)精度、運(yùn)算量和參數(shù)配對(duì)等問(wèn)題一直是DOA和多普勒頻率聯(lián)合估計(jì)算法的主要制約因素?，F(xiàn)有的DOA和多普勒頻率聯(lián)合估計(jì)方法大多由一維的子空間分解類(lèi)算法發(fā)展而來(lái)。文獻(xiàn)［1－5］基于ESPRIT算法實(shí)現(xiàn)了對(duì)DOA和多普勒頻率的聯(lián)合估計(jì)，雖然其中一些算法實(shí)現(xiàn)了參數(shù)的自動(dòng)配對(duì)，但其無(wú)法達(dá)到最佳的估計(jì)精度，且只適用于均勻線陣的情況，在低信噪比或采樣快拍數(shù)較小的條件下，其參數(shù)估計(jì)性能較差。文獻(xiàn)［6］提出了一種簡(jiǎn)單有效的頻率－空間－頻率（FSF）的樹(shù)狀結(jié)構(gòu)，利用三次一維 MUSIC／ESPRIT算法實(shí)現(xiàn)了對(duì)DOA和多普勒頻率的聯(lián)合估計(jì)，且參數(shù)可以自動(dòng)配對(duì)，但該方法得到的估計(jì)結(jié)果依然不是最優(yōu)估計(jì)。最大似然方法是一個(gè)性能優(yōu)良的最優(yōu)估計(jì)方法，但其求解過(guò)程需要進(jìn)行非線性的多維搜索，大大增加了算法的復(fù)雜度。文獻(xiàn)［7］利用最大似然方法實(shí)現(xiàn)了DOA和多普勒頻率的聯(lián)合估計(jì)，并將多維最優(yōu)化問(wèn)題轉(zhuǎn)換成多個(gè)一維最優(yōu)化問(wèn)題，簡(jiǎn)化了計(jì)算，但其估計(jì)誤差較大且需要額外的參數(shù)配對(duì)運(yùn)算。文獻(xiàn)［8－9］利用重要性采樣技術(shù)降低了求解似然函數(shù)最大值的計(jì)算量，并實(shí)現(xiàn)了DOA和多普勒頻率的參數(shù)自動(dòng)配對(duì)，但其需要根據(jù)信噪比的不同而不斷調(diào)整重要性采樣函數(shù)，且算法所需要的控制參數(shù)較多，不利于實(shí)際應(yīng)用。

本文將狀態(tài)空間模型和最大似然方法相結(jié)合，提出了一種基于最大似然方法的DOA和多普勒頻率聯(lián)合估計(jì)算法。利用狀態(tài)空間模型能夠分析和處理時(shí)變系統(tǒng)的特性，將待估計(jì)的參數(shù)轉(zhuǎn)換到狀態(tài)空間模型的系統(tǒng)矩陣中，利用Hankel矩陣構(gòu)造數(shù)據(jù)協(xié)方差矩陣來(lái)抑制噪聲，進(jìn)而通過(guò)最大似然方法得到信號(hào)角度和多普勒頻率的估計(jì)值，并利用人工蜂群算法對(duì)似然函數(shù)求解進(jìn)行優(yōu)化。該方法可得到漸近無(wú)偏的參數(shù)估計(jì)，參數(shù)能夠自動(dòng)配對(duì)，且大大減小了最大似然方法的計(jì)算量。

1 角度和多普勒頻率聯(lián)合估計(jì)

1.1 接收數(shù)據(jù)模型

假設(shè)接收陣列為具有M個(gè)陣元的線列陣，以陣列的第一個(gè)陣元作為基準(zhǔn)，陣元間距為（0，d1，…，dM－1），共有P個(gè)遠(yuǎn)場(chǎng)窄帶信號(hào)源，其中M ＞P，設(shè)波達(dá)方向?yàn)?（θ1，θ2，…，θP），且假設(shè)所有的信號(hào)都具有相同的載頻，則陣列接收到的信號(hào)可寫(xiě)為

式中:v（t）為復(fù)噪聲向量，假設(shè)為高斯白噪聲；sp（t）為信號(hào)的復(fù)包絡(luò)，可寫(xiě)成sp（t）＝αpej2πfpt，其中αp為復(fù)幅度，fp為第p個(gè)信號(hào)源的多普勒頻率；θp為第p個(gè)信號(hào)源的到達(dá)角度；a（θp）表示θp的方向向量，其定義為

式中:fc為信號(hào)的載頻；c為信號(hào)傳播速度。

利用歸一化頻率來(lái)描述多普勒頻率，式（1）可寫(xiě)成下列矩陣形式

由此可以得到如下形式的狀態(tài)空間模型

式中:ωp＝2πfp。

1.2 信號(hào)參數(shù)的估計(jì)

為了對(duì)噪聲進(jìn)行抑制，并能夠準(zhǔn)確地對(duì)波達(dá)方向和多普勒頻率進(jìn)行估計(jì)，本文方法從構(gòu)造Hankel矩陣開(kāi)始。

定義陣列接收數(shù)據(jù)和噪聲數(shù)據(jù)Hankel矩陣如下

式中:N、K分別為構(gòu)成Hankel矩陣的行數(shù)和列數(shù)；N＋K為采樣快拍數(shù)。

定義信號(hào)采樣向量為

則式（4）（5）的狀態(tài)空間模型可以被寫(xiě)成

式中:ΓK為廣義可觀測(cè)矩陣，其具有如下的形式

式（10）具有與陣列接收數(shù)據(jù)矩陣相同的形式，其可被定義為一種經(jīng)過(guò)擴(kuò)展的陣列接收數(shù)據(jù)矩陣。Hankel矩陣Y可表示為采樣快拍數(shù)為N的陣列接收數(shù)據(jù)矩陣，其每一列可被視作一次采樣，則式（10）可寫(xiě)作

式中:ΓK已包含了待估計(jì)的角度和多普勒頻率信息。

取觀測(cè)數(shù)據(jù)Y（n）的對(duì)數(shù)似然函數(shù):

式（13）是一個(gè)關(guān)于未知參量θ、ω、σ2和S（n）的非線性函數(shù)，其中未知參量σ2和S（n）的確定性最大似然估計(jì)分別為

式（16）的求解需要通過(guò)非線性多維搜索實(shí)現(xiàn)，由于ΓK是包含DOA和多普勒頻率信息的已知結(jié)構(gòu)矩陣，因此在搜索過(guò)程中可實(shí)現(xiàn)DOA和多普勒頻率的自動(dòng)配對(duì)。

2 似然函數(shù)的全局優(yōu)化

2.1 人工蜂群算法

人工蜂群算法（Artificial bee colony alogrithm，ABCA）是一種基于群體智能的啟發(fā)式搜索算法［10－11］，它模仿自然界中的蜂群總能尋找到最好的蜜源這種行為，具有控制參數(shù)少、易于實(shí)現(xiàn)、計(jì)算簡(jiǎn)潔等優(yōu)點(diǎn)，在求解組合優(yōu)化、多峰函數(shù)求極值等問(wèn)題中顯現(xiàn)出了優(yōu)良的特性。

人工蜂群算法將蜜蜂種群分為三部分:采蜜蜂、待工蜂和偵察蜂。對(duì)于函數(shù)優(yōu)化問(wèn)題，蜂群中的每個(gè)采蜜蜂對(duì)應(yīng)著一組解。開(kāi)始搜索時(shí)，由偵察蜂隨機(jī)搜索解空間內(nèi)的解，然后偵察蜂變成采蜜蜂，將所搜索到的解按照由好到壞（根據(jù)適應(yīng)度函數(shù)的高低）排序，并對(duì)其鄰域進(jìn)行搜索。如果鄰域搜索到的解優(yōu)于現(xiàn)有的最優(yōu)解，則替換現(xiàn)有的解，反之則保持現(xiàn)有的解不變。在所有的采蜜蜂都完成鄰域搜索后，采蜜蜂將所搜索到的解信息與待工蜂共享。待工蜂根據(jù)一定的概率選擇采蜜蜂進(jìn)行鄰域搜索，適應(yīng)度高的采蜜蜂吸引更多的待工蜂。待工蜂搜索到的鄰域解與采蜜蜂的現(xiàn)有最優(yōu)解進(jìn)行比較，如果優(yōu)于現(xiàn)有最優(yōu)解則替換該解；反之則保持現(xiàn)有解不變并繼續(xù)進(jìn)行鄰域搜索。鄰域搜索不能無(wú)限制地進(jìn)行，當(dāng)搜索達(dá)到一定次數(shù)但搜到的解仍沒(méi)有改善時(shí)，這組解將被放棄，所對(duì)應(yīng)的采蜜蜂變成偵察蜂重新開(kāi)始隨機(jī)搜索一組新的解。

在人工蜂群算法中，三種不同蜜蜂扮演了各自不同的角色。采蜜蜂用來(lái)保持優(yōu)良的解集；待工蜂用來(lái)加快全局最優(yōu)化的收斂速度；偵察蜂有助于跳出局部最優(yōu)解。整個(gè)蜂群分工協(xié)作，直到找到全局最優(yōu)解。

2.2 基于人工蜂群算法的似然函數(shù)優(yōu)化

假設(shè)蜂群總數(shù)為L(zhǎng)S，其中采蜜蜂和待工蜂數(shù)量分別為L(zhǎng)E和LO（令LE＝LO），搜索維數(shù)為2P，角度和多普勒頻率的搜索范圍分別為（－90°，90°） 和 （－π，π）。每一組解（i＝1，2，…，LE）為一個(gè)2P維的向量，所有的初始化解均由偵察蜂按照下面的公式進(jìn)行隨機(jī)搜索

式中:Xmin和Xmax分別為搜索范圍的最小值和最大值；β為在（0，1）范圍內(nèi)的隨機(jī)數(shù)。由此可以得到初始解集X2P。

偵察蜂在完成初始化搜索后轉(zhuǎn)變?yōu)椴擅鄯?，并進(jìn)行鄰域搜索。鄰域搜索可按照下式進(jìn)行:

式中:ca和cb為在（0，1）范圍內(nèi)的概率系數(shù)，其取值大小將決定待工蜂選擇變?yōu)椴擅鄯涞母怕?。適應(yīng)度高的解將吸引更多的待工蜂，適應(yīng)度低的解將吸引少量的待工蜂或者根本吸引不到待工蜂。

如果某一組解經(jīng)過(guò)一定次數(shù)的鄰域搜索仍然無(wú)法提高其適應(yīng)度，該組解將被放棄，其所對(duì)應(yīng)的采蜜蜂將轉(zhuǎn)變?yōu)閭刹旆渲匦码S機(jī)搜索一組新的解。

按照以上的步驟進(jìn)行迭代，最終將得到適應(yīng)度最高的一組解，該組解即為使似然函數(shù)最大的一組角度和多普勒頻率的估計(jì)值。

2.3 計(jì)算量分析

最大似然算法需要通過(guò)多維網(wǎng)格搜索實(shí)現(xiàn)，其計(jì)算復(fù)雜度為O｛［（θmax－ θmin）（ωmax－。其中，（θmin，θmax）和Δθ分別為角度搜索范圍和搜索步長(zhǎng)；（ωmin，ωmax）和Δω分別為頻率搜索范圍和搜索步長(zhǎng)。利用人工蜂群算法實(shí)現(xiàn)參數(shù)估計(jì)的計(jì)算復(fù)雜度約為O（2PLSIN＋LS）。其中，LS和IN分別為蜂群數(shù)量和迭代次數(shù)。

表1給出了兩種算法在不同信源數(shù)情況下的計(jì)算量比較。其中，角度搜索范圍和搜索步長(zhǎng)分別選取 （－90°，90°） 和0.5°；頻率搜索范圍和搜索步長(zhǎng)分別選取 （－π，π）和0.01rad；蜂群數(shù)量和迭代次數(shù)分別選取100和200。

表1 網(wǎng)格搜索與人工蜂群算法計(jì)算量比較Table 1 Computation load comparison of grid search and artificial bee colony algorithm

由表1可以看出，人工蜂群算法的計(jì)算量遠(yuǎn)小于網(wǎng)格搜索，特別是當(dāng)信源數(shù)增加時(shí)，網(wǎng)格搜索的計(jì)算量呈指數(shù)增長(zhǎng)，而人工蜂群算法的計(jì)算量呈線性增長(zhǎng)。當(dāng)信源數(shù)較多時(shí)，人工蜂群算法的優(yōu)勢(shì)更加明顯。

3 仿真驗(yàn)證

仿真實(shí)驗(yàn)采用線列陣，陣元數(shù)M＝10，采樣快拍數(shù)為100。在高斯白噪聲背景下，采用100次獨(dú)立Monte Carlo實(shí)驗(yàn)，以如下均方根誤差（RMSE）作為衡量算法性能的標(biāo)準(zhǔn):

實(shí)驗(yàn)1 分別采用人工蜂群算法、微分進(jìn)化算法和遺傳算法對(duì)似然函數(shù)進(jìn)行優(yōu)化求解。由于種群數(shù)量和迭代次數(shù)是決定群體智能優(yōu)化算法的關(guān)鍵參數(shù)，分別選取相同的種群數(shù)量和迭代次數(shù)可使以上三種優(yōu)化算法的計(jì)算量大致相等。因此，選取種群數(shù)量為50，迭代次數(shù)為100次。仿真實(shí)驗(yàn)取入射信號(hào)數(shù)P＝2，入射信號(hào)的入射角度分別為（20°，30°），入射信號(hào)的歸一化多普勒頻率分別為（0.7，0.9）。

圖1和圖2分別給出了三種算法在信噪比變化時(shí)DOA估計(jì)和多普勒頻率估計(jì)的均方根誤差的變化。從圖1、圖2中可看出，人工蜂群算法能夠很好地實(shí)現(xiàn)對(duì)似然函數(shù)的優(yōu)化，每次優(yōu)化均可以找到似然函數(shù)的最大值，具有很好的魯棒性。對(duì)于微分進(jìn)化算法和遺傳算法，由于種群數(shù)量和迭代次數(shù)的限制，算法無(wú)法保證每次都能找到全局最優(yōu)解，要想獲得更好的函數(shù)優(yōu)化性能，必須增加種群數(shù)量和迭代次數(shù)，這將大大增加算法的計(jì)算量。由此可見(jiàn)，人工蜂群算法可以用較小的計(jì)算量實(shí)現(xiàn)對(duì)似然函數(shù)較精確的求解。

圖1 不同優(yōu)化算法DOA估計(jì)均方根誤差比較Fig.1 RMSE of DOA estimation compared with other optimization methods

圖2 不同優(yōu)化算法頻率估計(jì)均方根誤差比較Fig.2 RMSE of frequency estimation compared with other optimization methods

實(shí)驗(yàn)2 選取兩個(gè)較具代表性的DOA和多普勒頻率聯(lián)合估計(jì)算法:JAFE［5］算法和 FSFMUSIC［6］算法，與本文算法進(jìn)行比較。由于以上兩種算法只適用于均勻線陣，故采用均勻線陣對(duì)算法進(jìn)行仿真?？紤]入射信號(hào)數(shù)P＝3的情況，信號(hào)源的波達(dá)方向分別為（20°，30°，40°），信號(hào)源的歸一化多普勒頻率分別為（0.7，0.9，1）。圖3和圖4分別給出了三種算法在信噪比變化時(shí)DOA估計(jì)和多普勒頻率估計(jì)的均方根誤差的變化。從圖3、圖4中可以看出，本文算法的估計(jì)性能優(yōu)于其他兩種算法，尤其在低信噪比下，本文算法能夠更早地貼近克拉美羅界。這主要得益于本文所使用的最大似然方法可得到漸近無(wú)偏的參數(shù)估計(jì)。

圖3 不同聯(lián)合估計(jì)算法DOA估計(jì)均方根誤差比較Fig.3 RMSE of DOA estimation compared with other joint estimation methods

圖4 不同聯(lián)合估計(jì)算法頻率估計(jì)均方根誤差比較Fig.4 RMSE of frequency estimation compared with other joint estimation methods

實(shí)驗(yàn)3 考慮三個(gè)入射信號(hào)源中的兩個(gè)信源參數(shù)比較接近的情況。信源參數(shù)選擇如表2所示，取SNR＝0dB，對(duì)JAFE算法、FSF－MUSIC算法和本文算法的估計(jì)性能進(jìn)行比較。由表2可見(jiàn)，本文算法能夠比較準(zhǔn)確地區(qū)分出各個(gè)信源的波達(dá)方向和多普勒頻率信息，而其他兩種算法的參數(shù)估計(jì)性能均受到了不同程度的影響，不能夠準(zhǔn)確分辨出參數(shù)比較接近的信號(hào)源。

表2 信號(hào)參數(shù)較接近時(shí)的RMSETable 2 RMSE when parameters are close

4 結(jié)束語(yǔ)

針對(duì)波達(dá)方向和多普勒頻率聯(lián)合估計(jì)問(wèn)題，將人工蜂群算法與最大似然估計(jì)相結(jié)合，提出了一種基于人工蜂群算法的DOA和多普勒頻率聯(lián)合估計(jì)方法。該方法利用狀態(tài)空間模型將所要估計(jì)的參數(shù)轉(zhuǎn)換到廣義可觀測(cè)矩陣中，再通過(guò)最大似然方法得到信號(hào)到達(dá)角和多普勒頻率的估計(jì)。由于求解似然函數(shù)計(jì)算量巨大，本文將人工蜂群算法應(yīng)用于似然函數(shù)的優(yōu)化求解中，大大減小了算法的計(jì)算量。與其他DOA和多普勒頻率聯(lián)合估計(jì)算法相比，本文方法具有估計(jì)精度高、控制參數(shù)少、魯棒性好、參數(shù)自動(dòng)配對(duì)等特點(diǎn)，且在低信噪比和小樣本下依然能夠得到較滿意的估計(jì)結(jié)果。

［1］Haardt M，Nossek J A.3－D unitary ESPRIT for joint 2－D angle and carrier estimation［C］∥IEEE International Conference on Acoustics，Speech and Signal Processing，Munich，1997.

［2］Lemma A N，van der Veen A J，Deprettere E F.Joint angle－frequency estimation using multi－resolution ESPRIT［C］∥IEEE International Conference on Acoustics，Speech and Signal Processing，Seattle，1998.

［3］Viberg M，Stoica P.A computationally efficient method for joint direction finding and frequency estimation in colored noise［C］∥The Thirty－Second Asilomar Conference on Signal，Systems and Computers，Pacific Grove，1998.

［4］Lemma A N，van der Veen A J，Deprettere E F.A－nalysis of ESPRIT based joint angle－frequency estimation［C］∥IEEE International Conference on A－coustics，Speech and Signal Processing，Istanbul，2000.

［5］Lemma A N，van der Veen A J，Deprettere E F.Analysis of joint angle－frequency estimation using ESPRIT［J］.IEEE Transactions on Signal Processing，2003，51（5）:1264－1283.

［6］Lin Jen－der，F(xiàn)ang Wen－h(huán)sien，Wang Yung－yi，et al.FSF MUSIC for joint DOA and frequency estimation and its performance analysis［J］.IEEE Transactions on Signal Processing，2006，54（12）:4529－4542.

［7］Atheley F.Asymptotically decoupled angle－frequency estimation with sensor arrays［C］∥The Thirtythird Asilomar Conference on Signal，Systems and Computers，Pacific Grove，1999.

［8］王惠剛，劉強(qiáng).多普勒方位聯(lián)合估計(jì)的蒙特卡洛算法［J］.電子學(xué)報(bào)，2009，37（9）:1965－1970.Wang Hui－gang，Liu Qiang.A Monte Carlo method for joint estimation of Dopplers and DOAs［J］.Acta Electronica Sinica，2009，37（9）:1965－1970.

［9］Wang Hui－gang，Steven Key.Maximum likelihood angle－Doppler estimator using importance sampling［J］.IEEE Transactions on Aerospace and Electronics Systems，2010，46（2）:610－622.

［10］Karaboga D，Basturk B.A powerful and efficient algorithm for numerical function optimization:artificial bee colony （ABC）algorithm［J］.Journal of Global Optimization，2007，39（3）:459－471.

［11］Karaboga D，Basturk B.On the performance of artificial bee colony（ABC）algorithm［J］.Applied Soft Computing，2008，8（1）:687－697.