張義元，張志成，石要武，劉昱兵

(吉林大學(xué)通信工程學(xué)院，長春130012)

0 引言

陣列信號處理是信號處理重點(diǎn)研究方向之一，近些年來得到迅速發(fā)展。信號的波達(dá)方向(DOA:Direction of Arrival)估計是陣列信號處理的重要研究內(nèi)容和持續(xù)研究熱點(diǎn)之一，在雷達(dá)、聲納、通信和地震學(xué)等領(lǐng)域應(yīng)用前景很大。最大似然(ML:Maximum Likelihood)參數(shù)估計方法是參數(shù)估計理論的典型估計方法，Ziskind首先將最大似然參數(shù)估計應(yīng)用于信號的 DOA估計［1，2］。與 MUSIC(Multiple Single Classification)算法和ESPRIT(Estimating Single Parameters via Rotational Invariance Techniques)算法相比，ML算法的估計性能良好，具有很好的穩(wěn)定性，尤其是低信噪比情況下，ML算法比MUSIC及其他子空間分解類算法性能強(qiáng)很多。但由于最大似然DOA估計的似然函數(shù)是多維非線性的極值問題，需要全局極值的多維搜索，運(yùn)算量非常大。

20世紀(jì)后期，仿生智能算法脫穎而出，如粒子群優(yōu)化算法(PSO:Particle Swarm Optimization)、遺傳算法(GA:Genetic Algorithm)和人工蜂群算法(ABC:Artificial Bee Colony)等，它們通過模擬自然界中生物的行為解決大規(guī)模復(fù)雜優(yōu)化問題［3，4］。針對最大似然DOA估計問題的復(fù)雜性，一些研究學(xué)者將仿生智能算法應(yīng)用于信號的DOA估計，取得了一些效果。Li等［5］將GA算法應(yīng)用于最大似然DOA估計，首次實現(xiàn)了仿生算法的最大似然DOA估計。但GA算法容易出現(xiàn)過早收斂等問題。李俊武等［6］利用PSO算法做最大似然DOA估計，性能良好，但由于PSO算法收斂速度較慢，容易陷入局部最優(yōu)解，在求解多個方向角同時估計時仍存在不足。

布谷鳥搜索算法(CS:Cuckoo Search)是YANG等［7］提出的仿生智能算法。該算法主要基于布谷鳥寄生繁殖機(jī)理和萊維飛行搜索兩方面。基于CS算法方法簡捷，求解復(fù)雜問題時優(yōu)化速度快。筆者將改進(jìn)后的CS算法應(yīng)用于最大似然算法DOA估計，解決最大似然DOA估計存在的計算量大，容易陷入局部最優(yōu)等問題，對于多個方向角的聯(lián)合估計，力圖獲得較好的估計精度和收斂性能。

1 信號的DOA估計

1.1 陣列信號模型

假設(shè)均勻線陣信號模型的陣元是均勻且各向同性的，陣元數(shù)為M，相鄰陣元間距為d，有P個窄帶遠(yuǎn)場信號源以平面波入射，入射角度為θk，入射波長為λ［8］，則陣元接收數(shù)據(jù)可寫成矢量形式

其中 X(t)=［x1(t)，x2(t)，…，xM(t)］T為 M 維的接收數(shù)據(jù)矢量;S(t)=［s1(t)，s2(t)，…，sP(t)］T為 P維的信號源矢量;N(t)=［n1(t)，n2(t)，…，nM(t)］T為M維的噪聲矢量，并假定噪聲為高斯白噪聲，在空間和時間上均獨(dú)立;A=［a(θ1)，a(θ2)，…，a(θP)］為M×P維的陣列導(dǎo)向矩陣，其第k列a(θk)為第 k個信號源 sk的導(dǎo)向矢量為a(θk)=［1，exp(jφ(θk))，…，exp(j(M －1)φ(θk))］T，φ(θk)=為第k個信號源矢量與陣列法向量的夾角。

通過以上陣列信號DOA估計的數(shù)學(xué)模型，可以得到陣列接收數(shù)據(jù)的協(xié)方差矩陣

其中RS為信號的協(xié)方差矩陣，RN為噪聲的協(xié)方差矩陣。

1.2 最大似然算法測向原理

假設(shè)噪聲為平穩(wěn)、空間和時間均不相關(guān)的高斯隨機(jī)過程，均值為零，方差為σ2;源信號為未知的確定性信號，則有

由概率論的理論［9］，幾個獨(dú)立且同分布高斯隨機(jī)過程的概率密度函數(shù)

對式(5)取對數(shù)，有

整理可得，最大似然算法的DOA估計譜函數(shù)

經(jīng)數(shù)學(xué)推導(dǎo)后可得代價函數(shù)

目標(biāo)函數(shù)式(8)是關(guān)于入射方向角θ1，θ2，…，θP的多維非線性函數(shù)，直接求解最優(yōu)值，運(yùn)算量非常大，且由于其通常具有復(fù)雜的局部極值結(jié)構(gòu)，很難保證收斂到全局最優(yōu)解;另一方面，求其導(dǎo)數(shù)較難。因此，按照傳統(tǒng)的優(yōu)化算法通過目標(biāo)函數(shù)值及其目標(biāo)函數(shù)導(dǎo)數(shù)值求解，難以解決全局極值問題。

2 最大似然DOA估計的函數(shù)優(yōu)化

2.1 布谷鳥搜索算法

布谷鳥搜索算法是隨機(jī)全局優(yōu)化算法，其靈感來源于布谷鳥將它育雛一些小布谷鳥的任務(wù)寄生在其他鳥類的巢中的現(xiàn)象。在自然界中，布谷鳥是通過隨機(jī)的或是類似隨機(jī)的方式尋找適合自己產(chǎn)蛋的鳥窩位置，為了能更好模擬布谷鳥尋窩的方式，首先需要設(shè)定以下3個理想的狀態(tài):

1)一只布谷鳥每次只產(chǎn)一個蛋，并隨機(jī)將這個蛋放到某個鳥窩中;

2)最佳鳥窩中高性能的蛋將會被保留到下一代;

3)鳥窩的數(shù)量是固定的，蛋放到鳥窩中被宿主發(fā)現(xiàn)的概率為Pa∈［0，1］。

在被宿主發(fā)現(xiàn)時，宿主鳥可選擇拋出異樣鳥蛋，或直接放棄該鳥巢，重新建造一個新的鳥巢［10］。

在CS搜索算法中，1個鳥巢表示1個候選解。布谷鳥搜索算法采用萊維隨機(jī)行走方式更新鳥窩位置

其中Xg，t表示第g代第t個解;α＞0表示步長的大小，用于控制隨機(jī)搜索的范圍。⊕是點(diǎn)乘積，L(β)服從Levy概率分布

基于上述假設(shè)，布谷鳥搜索算法的主要步驟描述如下:

Step1 選取目標(biāo)函數(shù)f(X)，X=(x1，…，xd)T，隨機(jī)產(chǎn)生一組鳥窩的初始位置Xi(i=1，2，…，n)，并設(shè)置算法參數(shù);

Step2 計算每個鳥窩的目標(biāo)函數(shù)值，并記錄當(dāng)前的最好解;

Step3 保留上一代中最好的鳥窩位置，對其他的鳥窩位置按式(9)進(jìn)行更新;

Step4 將新得到一組鳥窩位置與上一代進(jìn)行對比，若較好，則將該組最優(yōu)解作為當(dāng)前的最優(yōu)位置;

Step5 用服從均勻分布的隨機(jī)數(shù)R∈［0，1］作為鳥窩宿主發(fā)現(xiàn)外來鳥蛋的可能性與Pa比較，若R＞Pa，則隨機(jī)改變鳥窩位置，得到一組新位置;

Step6 若未滿足結(jié)束條件，則返回Step2;

Step7 輸出全局最優(yōu)位置。

2.2 基于布谷鳥搜索算法的最大似然DOA估計

將布谷鳥搜索算法應(yīng)用于空間信號的DOA估計問題。最大似然算法的譜函數(shù)P(θ)作為待搜索的目標(biāo)函數(shù)，每組鳥窩位置Xi(i=1，2，…，n)代表待估計方向角在第t次迭代中的n個解θi(t)(i=1，2，…，n)。經(jīng)過每次迭代，按照式(9)得到一組新解θi(t+1)，并對基本的布谷鳥搜索算法作如下改進(jìn)。

1)將初始化過程產(chǎn)生N個標(biāo)準(zhǔn)正態(tài)分布的初始化種群，改為服從(－π/2，π/2)的均勻分布的N個初始化種群。

2)對式(9)更改步長控制方式，具體處理如下。

①為能從當(dāng)前最優(yōu)解中獲得更多有用的步長信息［11］，計算步長

并加入慣性權(quán)重α0=0.1(1－k/M)，其中M為最大迭代次數(shù)，k為當(dāng)前的迭代次數(shù)，Xbest表示當(dāng)前最優(yōu)解。

②為了便于計算，采用式(12)計算Levy隨機(jī)數(shù)

其中u，v服從標(biāo)準(zhǔn)正態(tài)分布，β=1.5。

生成新的解。

如被宿主發(fā)現(xiàn)，在按一定的概率(發(fā)現(xiàn)概率Pa)丟棄部分解之后，算法采用偏好隨機(jī)游動重新生成相同數(shù)量的解

其中r為縮放因子，為服從(0，1)均勻分布的隨機(jī)數(shù);Xg，t和Xg，k表示g代的兩個隨機(jī)解。

因此，按照布谷鳥搜索算法的程序和式(9)進(jìn)行算法的迭代，經(jīng)過一定次數(shù)的迭代，可得到相應(yīng)性能良好的解，作為最大似然算法DOA估計的方向角。

顯然，綜合式(11)～式(13)，在Levy飛行中，CS搜索算法采用下面等效的位置更新公式

3 實驗與結(jié)果

在Matlab仿真實驗平臺上，選用其他幾種較受歡迎的元啟發(fā)式仿生智能算法，包括微分進(jìn)化算法(DE:Differential Evolution)，粒子群優(yōu)化算法(PSO:Particle Swarm Optimization)和克隆選擇算法(CLONALG:Clonal Genetic)，分別做最大似然DOA估計，比較布谷鳥搜索算法與上述仿生智能算法的收斂性及其精度。

實驗1 收斂實驗。陣列模型選用均勻線陣，陣元數(shù)為10。陣元之間的間距設(shè)為信號載波波長的一半，快拍數(shù)為100，噪聲為0的高斯白噪聲。搜索線陣范圍為(－90°，90°);4種仿生優(yōu)化算法的參數(shù)設(shè)置如表1所示。其中Q為種群數(shù)量，I為最大迭代次數(shù)，F(xiàn)為變異因子，c為交叉概率，c1為加速因子1，c2為加速因子2，m為變異概率，r為抗體個數(shù)，p為發(fā)現(xiàn)概率，a為搜索步長。圖1是入射角為(20°，50°)時最優(yōu)目標(biāo)函數(shù)值隨迭代次數(shù)的變化情況，圖2是入射角為(20°，50°，－10°)時最優(yōu)目標(biāo)函數(shù)值隨迭代次數(shù)的變化情況。從圖1可以看出，有兩個入射信號時，布谷鳥搜索算法較其他的仿生優(yōu)化算法具有較快的收斂特性;而圖2中，在同時有3個入射信號時，布谷鳥搜索算法的收斂特性也較好，搜索到的似然函數(shù)值更接近于全局最優(yōu)解，所以在入射信號數(shù)目增多、待優(yōu)化目標(biāo)函數(shù)復(fù)雜的情況下，布谷鳥搜索算法的收斂速度和全局優(yōu)化能力較為突出。

表1 不同優(yōu)化算法的參數(shù)設(shè)置Tab.1 Parameter settings of different algorithm

圖1 兩個入射信號時布谷鳥搜索算法和其他算法的收斂性比較Fig.1 Comparison of the convergence for cuckoo search algorithm and other algorithms when the two incident signals

圖2 3個入射信號時布谷鳥搜索算法和其他算法的收斂性比較Fig.2 Comparison of the convergence for cuckoo search algorithm and other algorithms when the three incident signals

實驗2 算法的精度實驗。所有條件均與實驗1相同，用均方根誤差

作為衡量算法精度的標(biāo)準(zhǔn)。其中N為信號源的個數(shù)，NMC為實驗次數(shù)，θi表示第i個信源的真實到達(dá)角，^θi(k)為第k次估計的到達(dá)角的估計值。實驗選取信噪比的變化范圍在－15～15 dB之間，在高斯白噪聲背景下，采用100次獨(dú)立Monte Carlo實驗，得到均方根誤差隨信噪比的變化曲線比較如圖3所示。從圖3可以看出，幾種優(yōu)化算法在低信噪比時誤差較為接近，而隨著信噪比的增加，布谷鳥搜索算法優(yōu)化的DOA譜估計的估計誤差較低，估計的穩(wěn)定性較好。

圖3 與其他優(yōu)化算法比較的DOA估計均方根誤差Fig.3 RMSE of DOA estimation compared with other optimization methods

4 算法的計算量分析

最大似然DOA估計算法需要通過譜峰搜索實現(xiàn)參數(shù)估計，計算量主要由信號參數(shù)的維數(shù)和搜索步長決定。以估計兩個入射角為例，計算量為 O((θmax1－θmin1)(θmax2－ θmin2)/Δθ1Δθ2)，其中(θmin，θmax)和 Δθ分別為DOA的搜索范圍和搜索步長。而基于布谷鳥搜索算法的最大似然DOA估計中譜峰搜索的計算量主要由種群數(shù)量和迭代次數(shù)決定，可表示為O(NIN+N)。其中N和IN分別為種群數(shù)量和迭代次數(shù)。

取DOA搜索范圍為(－90°，90°)，傳統(tǒng)網(wǎng)格搜索方法和筆者方法的計算量比較如表2所示。傳統(tǒng)網(wǎng)格搜索方法的步長Δθ=0.1;筆者方法種群數(shù)量為100，迭代次數(shù)為100。從表2中可以清楚地看到，在最大似然DOA估計中，同傳統(tǒng)的網(wǎng)格搜索方法相比，布谷鳥搜索算法可大幅度減少DOA估計的計算量。

表2 算法的計算量比較Tab.2 Conparison of calculation amount

5 結(jié) 語

筆者將最大似然DOA譜估計和布谷鳥搜索算法相結(jié)合，提出了基于布谷鳥搜索算法的最大似然DOA估計方法。該方法在保留最大似然DOA估計良好性能的同時，加入了改進(jìn)的布谷鳥搜索策略，實現(xiàn)了對最大似然譜函數(shù)的優(yōu)化搜索，大幅度降低了DOA估計的計算量。仿真實驗表明，在相同的實驗條件下，與其他方法相比，該算法具有搜索精度高、優(yōu)化速度快、全局搜索能力強(qiáng)等特點(diǎn)。而在待估計信號源個數(shù)增多，最大似然DOA譜函數(shù)更加復(fù)雜的情況下，基于CS優(yōu)化的最大似然DOA估計也能實現(xiàn)較快、較好的估計。所以筆者提出的改進(jìn)CS優(yōu)化策略，是解決此類復(fù)雜DOA譜峰搜索問題的有效方法。

［1］LEE J Y.Acoustic DOA Estimation:An Approximate Maximum Likelihood Approach ［J］.IEEE Systems Journal，2013，8(7):131－141.

［2］王永德，陳旗.基于最大似然估計的空間譜測向技術(shù)［J］.計算機(jī)與數(shù)字工程，2010，38(9):123－127.WANG Yongde，CHEN Qi.Direction Finding of Spatial Spectrum Based on Maximum Likelihood Estimation［J］.Computer＆Digital Engineering，2010，38(9):123－127.

［3］CIVICIOGLU P.A Conception Comparison of the Cuckoo Search，Particle Swarm Optimization Differential Evolution and Artificial Bee Colonial Gorithms［J］.Artificial Intelligence Review，2011，39(4):315－346.

［4］YANG X S，DEB S.Multi－Objective Cuckoo Search for Design Optimization［J］.Computers ＆ Operations Research，2013，40(6):1616－1624.

［5］LI M，LU Y.Genetic Algorithm Based Maximum Likelihood DOA Estimation［C］∥Proceedings of Radar Conference 2002.Edinburgh，UK:［s.n.］，2002，10:502－506.

［6］李俊武，俞志富.改進(jìn)粒子群算法在DOA估計中的應(yīng)用［J］.計算機(jī)工程與應(yīng)用，2013，49(9):203－206.LI Junwu，YU Zhifu.Improved PSO and Its Application to DOA Estimation ［J］.Computer Engineering and Applications，2013，49(9):203－206.

［7］YANG X S，DEB S.Cuckoo Search via Lévy Flights［C］∥Proc of the World Congress on Nature and Biologically Inspired Computing.Piscataway:IEEE Publications，2009:210－214.

［8］ZHANG Zhicheng，SHI Yaowu.Application of Artificial Bee Colony Algorithm to Maximum Likelihood DOA Estimation［J］.Journal of Bionic Engineering，2013，10(1):100－109.

［9］WANG H，KAY S.Maximum Likelihood Angle－Doppler Estimator Using Importance Sampling ［J］.IEEE Transactions on Aerospace and Electronic Systems，2010，46(2):610－622.

［10］李煜，馬良.新型元啟發(fā)式布谷鳥搜索算法［J］.系統(tǒng)工程，2012，30(8):64－69.LI Yu，MA Liang.The New Meta－Heuristic Cuckoo Search Algorithm ［J］.Systems Engineering，2012，30(8):64－69.

［11］王李進(jìn)，尹義龍.逐維改進(jìn)的布谷鳥搜索算法［J］.軟件學(xué)報，2013，24(11):2687－2698.WANG Lijin，YIN Yilong.Cuckoo Search Algorithm with Dimension by Dimension Improvement［J］.Journal of Software，2013，24(11):2687－2698.