小生境粒子群算法在多信號DOA估計中的應(yīng)用*

朱 任，莊銘杰

(華僑大學(xué) 工學(xué)院，福建 泉州 362021)

0 引 言

空間譜估計是智能天線中陣列信號處理的一個重要研究方向，在現(xiàn)代通信、聲吶、雷達(dá)和導(dǎo)航等領(lǐng)域中被廣泛應(yīng)用。波達(dá)方向(Direction of Arrival,DOA)是空間譜估計中的一項重要參數(shù)，快速求解高精度DOA對智能天線系統(tǒng)的應(yīng)用具有重要意義。多重信號分類(Multiple Signal Classification,MUSIC)算法[1]具有高精度和高分辨的特點，是估計DOA參數(shù)的經(jīng)典算法之一。但是MUSIC算法需要對空間譜函數(shù)進(jìn)行網(wǎng)格搜索來尋找極大值點，即譜峰搜索。譜峰搜索通常是一個計算量大、耗時長的窮舉最優(yōu)化問題，在高精度搜索時會大幅降低算法的實時性，因此保持MUSIC算法的譜峰搜索高精度同時降低計算量是一項很有價值的研究。

近年來，以遺傳算法(Genetic Algorithm,GA)、粒子群算法(Particle Swarm Optimization,PSO)等為代表的經(jīng)典進(jìn)化算法在許多領(lǐng)域的問題優(yōu)化、求解方面表現(xiàn)出優(yōu)異的性能[2]，有許多研究者將其改進(jìn)后應(yīng)用在譜峰搜索中：文獻(xiàn)[3-6]分別將雞群算法、GA、PSO和差分進(jìn)化算法應(yīng)用到DOA估計中，經(jīng)過仿真證實了這類算法能有效降低譜峰搜索的計算量，大幅提升DOA估計的精度和實時性，在實際工程應(yīng)用中具有潛在的應(yīng)用價值。

但是，許多情況下空間譜函數(shù)并不只有一個“峰”。在多用戶環(huán)境下，陣列同時接收到多個入射信號，此時空間譜函數(shù)的圖形會出現(xiàn)多個與DOA相對應(yīng)的極大值點，因此這是一個多峰函數(shù)優(yōu)化問題。顯然，經(jīng)典進(jìn)化算法求解此類多峰值優(yōu)化問題卻無能為力，它們只能搜索一個最優(yōu)的極大值點，而無法找到其他的極大值點。對于這類問題，引入小生境機(jī)制是一種有效的解決辦法，它模擬了自然環(huán)境中不同種群占據(jù)不同類型資源的生存方式，每個小生境都代表著解空間中包含一個局部解的區(qū)域以及該區(qū)域內(nèi)的粒子。這些粒子只會在所屬小生境內(nèi)尋找最優(yōu)解，這樣有利于保持種群的多樣性來防止收斂于全局最優(yōu)解，非常適合用來解決多峰函數(shù)優(yōu)化問題。文獻(xiàn)[7-9]在幾種進(jìn)化算法的基礎(chǔ)上，分別提出了基于并行免疫、排擠、隔離策略的小生境技術(shù)，并通過仿真驗證了算法的有效性。文獻(xiàn)[10]分析了PSO的拓?fù)浣Y(jié)構(gòu)并提出多種不同鄰域結(jié)構(gòu)的小生境PSO。適應(yīng)值與歐氏距離比PSO(Fitness Euclidean Distance Ratio PSO,FER-PSO)[11]利用適應(yīng)值與歐氏距離的比值來調(diào)整粒子位置，實驗結(jié)果表明算法能可靠地定位多個局部最優(yōu)解。鄰域擁擠差分進(jìn)化算法(Neighborhood based Crowding DE,NCDE)[12]是一種鄰域變異策略的差分進(jìn)化算法，變異是在每個基于歐氏距離的鄰域內(nèi)進(jìn)行的，因而能保持進(jìn)化過程中發(fā)現(xiàn)的最優(yōu)解，實驗證明能很好解決多峰搜索問題。文獻(xiàn)[13]模擬了蒲公英種子的傳播方式，提出了小生境蒲公英算法并應(yīng)用于DOA估計中。雖然實驗結(jié)果顯示有效，但是該算法的運(yùn)行前提是設(shè)置合適的小生境半徑參數(shù)，而這個參數(shù)的選取則依賴于待求解問題的先驗知識，其中包括解的分布間隔、疏密程度等。在沒有相關(guān)先驗知識的情況下，算法很難順利求解問題，這也是許多進(jìn)化計算算法普遍存在的缺點。在諸如雷達(dá)跟蹤、聲吶探測等許多工程應(yīng)用中，盡管適用于多峰函數(shù)求解的小生境優(yōu)化算法很多，但仍難以滿足MUSIC算法要求的實時處理能力和多信號DOA估計能力，同時還要盡可能減少對先驗知識的依賴等?？傊?，目前的小生境算法可以解決多峰函數(shù)搜索的優(yōu)化問題，但難以完全滿足應(yīng)用中對DOA估計的全部要求。

1 基于MUSIC算法的DOA估計

1.1 均勻直線陣列和L形陣列

如圖1示，yoz平面的y軸上放置M根天線，按均勻直線陣列(Uniform Linear Array,ULA)方式排列。

圖1 陣列模型

令相鄰兩陣元的間距d，達(dá)波方向與z軸的夾角為α，S(t)表示t時刻入射到陣列上的帶限遠(yuǎn)場信號，則該達(dá)波信號到達(dá)第m個陣元時相對第1個陣元所產(chǎn)生的時延大小如式(1)所示，其中c為光速。

(1)

由時延易得陣元間的相位差。因此對于ULA，達(dá)波信號的導(dǎo)向向量為

(2)

ULA的陣列流形由關(guān)于各個入射信號的導(dǎo)向向量組成，當(dāng)有k個入射信號時，ULA的陣列流形為

A=[a(α1),a(α2),…,a(αk)]，

(3)

則t時刻經(jīng)ULA接收處理后信號X(t)的矩陣表達(dá)式為

X(t)=AS(t)+N(t)。

(4)

式中：N(t)是均值為0、方差為σ2的高斯白噪聲，與入射信號不相關(guān)。陣列接收數(shù)據(jù)的協(xié)方差矩陣Rx可改寫為

Rx=E[XXH]=AE[SSH]AH+σ2I=ARSAH+σ2I。

(5)

式中：RS是信號的相關(guān)矩陣，I是單位矩陣。由于信號與噪聲相互獨立，Rx可經(jīng)特征值分解為信號和噪聲兩部分：

(6)

其中，US和UN分別是由較大特征值和較小特征值所對應(yīng)的特征向量組成的信號子空間和噪聲子空間。由于這兩個子空間是正交的，因此導(dǎo)向向量aH(α)和噪聲子空間UN同樣也是正交的，即

aH(α)UN=0。

(7)

由于在實際應(yīng)用中采集到的數(shù)據(jù)長度是有限的，無法直接得到接收數(shù)據(jù)的協(xié)方差矩陣，所以通過L次有限快拍樣本來估算：

(8)

通過樣本估算的協(xié)方差矩陣無法滿足信號子空間與噪聲子空間完全正交，即式(7)為趨近于0的數(shù)值。應(yīng)用其倒數(shù)為極大值的原理，MUSIC算法以最小優(yōu)化搜索的方式建立式(9)的功率譜函數(shù)來估計達(dá)波方向：

(9)

通過搜索式(9)函數(shù)的極大值點，就可以估算信號的入射角度α。

ULA無法同時檢測三維空間中的俯仰角和方位角，因此本文采用如圖1所示的L形陣列結(jié)構(gòu)。令x軸和y軸各有M個間距均為d的均勻分布陣元，達(dá)波信號S與z軸夾角為α，x軸與信號在xoy平面上投影夾角為β，則L形陣列的導(dǎo)向向量為

(10)

當(dāng)有k個入射信號時L形陣列的流形為

(11)

式中：x軸上的陣列流形Ax和y軸上的陣列流形Ay分別為

Ax=[ax(α1,β1),ax(α2,β2),…,ax(αk,βk)]，

(12)

Ay=[ay(α1,β1),ay(α2,β2),…,ay(αk,βk)]。

(13)

所以L型陣列的功率譜函數(shù)為

(14)

同樣地，通過搜索式(14)的極大值點，就可以估算出達(dá)波信號的俯仰角α和方位角β。

1.2 譜峰搜索及復(fù)雜度

假設(shè)有一圖1中yoz平面上的ULA，陣元數(shù)M=6，陣元間距d=λ/2；有三個用戶信號的DOA分別為α=-60°、α=20°和α=40°，SNR=20 dB，快拍數(shù)為500。將式(9)在α∈[-90°,90°]范圍內(nèi)繪制得到圖2(a)的空間譜圖形。另假設(shè)有一圖1中xoy平面上的L形陣列，在x軸和y軸各有陣元數(shù)M=6，陣元間距d=λ/2;有兩個用戶信號分別從α=30°、β=60°和α=60°、β=120°入射，按照式(14)的譜函數(shù)在α∈[0°,90°]、β∈[0°,180°]范圍內(nèi)繪制得出圖2(b)的空間譜圖。其中，所有用戶接收信號的SNR=20 dB，快拍數(shù)均為500。圖2中，空間譜圖形出現(xiàn)了與信號源數(shù)量相同、高度相近的極大值點(即譜峰)，其對應(yīng)的坐標(biāo)角度就是用戶信號的DOA，因此MUSIC算法的最后一個步驟就是使用網(wǎng)格搜索的方式尋找空間譜函數(shù)的峰值。網(wǎng)格的大小將決定DOA估計的復(fù)雜度和精度：一元和二元網(wǎng)格搜索的時間復(fù)雜度分別為O(n)和O(mn)，其中m和n是α和β軸的網(wǎng)格數(shù)量。網(wǎng)格搜索的誤差小于網(wǎng)格的步長，例如步長為1°的網(wǎng)格搜索得出的結(jié)果誤差在1°以內(nèi)。

(a)一元MUSIC算法空間譜圖形

為了展示網(wǎng)格搜索的性能，仿真實驗選擇圖2的譜函數(shù)圖形，仿真環(huán)境:CPU為Ryzen 5 3600，平臺為Matlab2018。表1和表2給出了網(wǎng)格搜索的結(jié)果和用時，結(jié)果顯示搜索用時隨著網(wǎng)格步長縮小而大幅增長，導(dǎo)致MUSIC算法不適用于實時性要求高的場合中，而若選擇大步長則使得結(jié)果得誤差較大。因此，改進(jìn)和優(yōu)化譜峰搜索，使其在保證高精度的前提下降低計算量，是本文的主要研究內(nèi)容。

表1 一元等步長網(wǎng)格搜索性能

表2 二元等步長網(wǎng)格搜索性能

1.3 經(jīng)典進(jìn)化算法在譜峰搜索中的問題

以遺傳算法、粒子群算法為代表的經(jīng)典進(jìn)化算法在單信號DOA估計中表現(xiàn)出優(yōu)異的性能，但卻難以估計多信號的DOA。以圖2(a)中的譜函數(shù)為例，使用經(jīng)典PSO進(jìn)行優(yōu)化，算法參數(shù)選?。篊1=C2=2.05，iwt=0.729 8，粒子數(shù)50。搜索結(jié)果如圖3所示，其中黑色圓圈表示粒子，紅色三角表示最優(yōu)粒子。

圖3 經(jīng)典PSO優(yōu)化圖2(a)的效果圖

從圖3可以看出，經(jīng)典PSO僅能找出α=40°處的極大值點，而卻無法發(fā)現(xiàn)α=-60°和α=20°兩處的極大值點。所有粒子都趨向于最優(yōu)粒子，降低了粒子的多樣性，從而失去發(fā)現(xiàn)其他極大值點的機(jī)會。為了降低最優(yōu)粒子對全局粒子的吸引，本文引入小生境思想將解空間劃分成多個區(qū)域，避免PSO在多峰的譜峰搜索中出現(xiàn)上述問題。

2 一種改進(jìn)的粒子群算法

針對網(wǎng)格搜索計算量大而經(jīng)典進(jìn)化計算無法優(yōu)化多峰譜函數(shù)的情況，本文提出了小生境粒子群算法。算法主要由中點順序聚類和基于迭代數(shù)的搜索策略構(gòu)成。本節(jié)將描述算法的原理并給出偽代碼。

2.1 中點順序聚類算法

順序聚類(Sequential Algorithmic Scheme,SAS)是一種簡單有效的聚類方式，可以將對象的集合劃分成由相似對象組成的數(shù)個類。因此，順序聚類結(jié)合PSO應(yīng)用時，可以根據(jù)粒子間距和函數(shù)值將隨機(jī)分布的粒子劃分成數(shù)個小粒子群，形成多個小生境。本文將粒子的坐標(biāo)中點信息用于改進(jìn)順序聚類，提出了一種中點順序聚類(Midpoint SAS,MSAS)算法。

在MSAS中，f(x)表示粒子x的適應(yīng)值，將所有粒子按照適應(yīng)值大小進(jìn)行降序排列，隨后從大到小遍歷每個粒子。如圖4(a)所示，計算粒子Pi和距離它最近的小生境nicheN中最優(yōu)秀粒子headN的歐氏距離中點midi，若midi的適應(yīng)值小于Pi和headN，則認(rèn)為兩者之間存在一個“山谷”，Pi不屬于nicheN。算法會建立一個新的小生境nichenew，并將Pi歸入其中，如式(15)所示；否則，如圖4(b)、(c)所示，Pi屬于離它最近的小生境nicheN，算法連同中點midi一起歸入nicheN，并比較這三個點的適應(yīng)值，更新headN。

(a)中點適應(yīng)值最小的情況

(15)

圖5的黑色曲線是式(16)在0.1≤x<3.4取值范圍內(nèi)的圖形，圖中的“*”和“o”是運(yùn)行MSAS將30個隨機(jī)分布粒子聚類的結(jié)果，其中不同小生境的成員以不同顏色區(qū)分，“*”代表一個小生境中的最優(yōu)粒子，“o”代表其他粒子。通過觀察圖5可知式(16)具有5個高度不同的“山峰”，且每個“山峰”的范圍不同。經(jīng)過MSAS對隨機(jī)粒子聚類后，搜索空間被劃分為5個不同大小的小生境，可見MSAS具有無需相關(guān)參數(shù)、自適應(yīng)調(diào)整小生境大小的能力，對待這類大小不一、分布不均勻的多峰函數(shù)十分有效。

圖5 中點順序聚類優(yōu)化式(16)的結(jié)果

(16)

中點順序聚類算法的偽代碼如下：

輸入:粒子總數(shù)m,目標(biāo)函數(shù)function

根據(jù)適應(yīng)值大小將所有粒子降序排列。

P1∈niche1,head1←P1(第一個粒子屬于第一個小生境，并記錄最優(yōu)值)

遍歷：從i=2到m

找到Pi最近的小生境nicheN，并計算中點midi。

如果：滿足式(15)：

Pi∈nichenew,headnew←Pi(Pi屬于新的小生境，并記錄最優(yōu)值)

否則：

Pi∈nicheN,比較Pi、headN和midi適應(yīng)值大小，更新headN。

2.2 粗略搜索與精細(xì)搜索

經(jīng)過MSAS的處理后NPSO可以獲取到大部分小生境的模糊位置和大小，并將每個已知小生境的年齡設(shè)置為1。之后，由于初始的隨機(jī)粒子沒能落到個別隱蔽的“山峰”中，導(dǎo)致這些小生境未能被建立，所以算法開始執(zhí)行粗略搜索來對解空間進(jìn)行充分探索。粗略搜索是發(fā)生在小生境之間的搜索，目的是模糊定位所有局部最優(yōu)解。如式(17)所示，粗略搜索會遍歷每個小生境，計算其最優(yōu)粒子head與相鄰小生境的最優(yōu)粒子headN位置之間的歐氏距離，以該最優(yōu)粒子為中心、距離一半作為半徑，在目標(biāo)函數(shù)上隨機(jī)產(chǎn)生一個粒子，最后將產(chǎn)生的所有粒子通過MSAS劃分到小生境中。

(17)

為了增強(qiáng)粗略搜索發(fā)現(xiàn)較隱蔽“山峰”的能力，每個小生境會在整個解空間中隨機(jī)產(chǎn)生一個粒子。每經(jīng)過一次粗搜的小生境年齡會增加1。

經(jīng)過數(shù)輪的粗略搜索，NPSO已經(jīng)能發(fā)現(xiàn)所有局部最優(yōu)解的模糊位置，建立了對應(yīng)數(shù)量的小生境，但是還無法精確定位每個小生境的最優(yōu)粒子位置。此時，達(dá)到預(yù)設(shè)年齡參數(shù)的小生境會轉(zhuǎn)為精細(xì)搜索策略：精搜是發(fā)生在小生境內(nèi)部的搜索，目的是找出最優(yōu)解的精確位置。如式(18)所示，隨機(jī)選取小生境內(nèi)一點P，以P至本小生境內(nèi)最優(yōu)粒子head的距離一半為半徑，以head為中心，在目標(biāo)函數(shù)上隨機(jī)生成一個粒子，并和head比較適應(yīng)值，判斷是否需要更新最優(yōu)粒子。

(18)

在算法初始階段，小生境還未被完全發(fā)現(xiàn)或者其中的最優(yōu)解還不能模糊定位時，粒子和已知的小生境頭部之間的關(guān)系十分復(fù)雜。例如，當(dāng)粒子和已知最近的小生境之間隔著一個未被發(fā)現(xiàn)的小生境時，MSAS會錯誤判斷。因此，NPSO在每輪迭代的結(jié)尾都會判斷兩個相鄰的小生境的最優(yōu)粒子是否滿足MSAS的規(guī)則，對錯誤的小生境進(jìn)行合并。小生境粒子群算法(NPSO)偽代碼如下：

輸入:最大粒子數(shù)m，搜索范圍，預(yù)設(shè)年齡，最大迭代數(shù)，目標(biāo)函數(shù)。

在搜索范圍內(nèi)隨機(jī)生成m個粒子。

運(yùn)行MSAS對粒子聚類，劃分小生境。

循環(huán)：從迭代數(shù)=1到最大迭代數(shù)。

循環(huán)：從i=1 到小生境總數(shù)。

如果：年齡<=預(yù)設(shè)年齡

根據(jù)式(17)產(chǎn)生粒子。

在搜索范圍內(nèi)隨機(jī)產(chǎn)生粒子。

否則：

根據(jù)式(18)產(chǎn)生粒子。

MSAS對產(chǎn)生的粒子聚類，劃分小生境。

判斷兩個相鄰的小生境的最優(yōu)粒子是否滿足MSAS中的規(guī)則并合并誤判的小生境。

保持粒子總數(shù)不超過最大粒子數(shù)m。

返回每個小生境的最優(yōu)粒子。

3 仿真與分析

3.1 仿真環(huán)境與參數(shù)

假設(shè)有一圖1中xoy平面的L形陣列，其中x和y軸上各有6個間距半波長的均勻分布陣元。陣列接收信號受到零均值、方差為σ2的加性高斯噪聲信號干擾，信噪比為20 dB，快拍數(shù)500。仿真所采用的空間譜參數(shù)如表3所示，三個樣例分別有2、3和4個譜峰，且非均勻分布，搜索范圍均為α∈[0°,90°],β∈[0°,180°]。選擇數(shù)量不同、分布不均的樣例是為了能更好地考察算法對小生境范圍的自適應(yīng)能力。由于干擾噪聲的影響，譜峰坐標(biāo)和DOA之間存在一定誤差，因此譜峰坐標(biāo)由0.001°邊長的網(wǎng)格搜索得出。本節(jié)中的誤差均是指算法輸出結(jié)果與譜峰坐標(biāo)之間的誤差。

表3 空間譜樣例的參數(shù)

為了比較算法在相同計算資源下的性能表現(xiàn)，實驗限制了目標(biāo)函數(shù)(譜函數(shù))的最大求值次數(shù)(Max Function Evaluation,MaxFEs)為10 000次。同時使用峰識別率(Peak Ratio,PR)來統(tǒng)計算法輸出結(jié)果的準(zhǔn)確度。式(19)是PR的定義式，由算法發(fā)現(xiàn)的解數(shù)量除以總的解數(shù)量得出，若輸出結(jié)果和實際坐標(biāo)間的歐氏距離達(dá)到精度要求，則認(rèn)為算法成功獲取到該解。

(19)

為了確定NPSO中的預(yù)設(shè)年齡參數(shù)，選取粒子數(shù)為200，預(yù)設(shè)年齡為1～20，并選擇表3中的樣例C，運(yùn)行500次蒙特卡洛仿真并統(tǒng)計達(dá)到預(yù)設(shè)年齡時的PR，結(jié)果如表4所示。因為NPSO在達(dá)到預(yù)設(shè)年齡前的搜索策略是粗搜，更注重搜索的廣度，旨在發(fā)現(xiàn)解的大致位置，故此處精度要求設(shè)置為5。

表4 NPSO不同預(yù)設(shè)年齡的仿真結(jié)果

表4中PR隨著預(yù)設(shè)年齡增大而逐漸增大并趨近于100%。觀察發(fā)現(xiàn)當(dāng)預(yù)設(shè)年齡大于10之后PR的增量逐漸減少，呈現(xiàn)邊際收益遞減的規(guī)律。為了平衡搜索的廣度和深度，預(yù)設(shè)年齡取10較為適宜。為全面體現(xiàn)NPSO算法是一種高效的多峰問題優(yōu)化算法，表5給出了6種同類算法及參數(shù)，它們是一組性能比較優(yōu)秀的算法及其參數(shù)，NPSO將在多信號DOA估計中與它們進(jìn)行比較。

表5 算法參數(shù)

3.2 算法在低計算資源下的性能

仿真實驗采用的譜峰搜索對象是表3中給出的樣例A、B和C，粒子數(shù)均取200，每種算法各運(yùn)行500次蒙特卡洛仿真。表6統(tǒng)計了算法的運(yùn)行時間,而表7～9統(tǒng)計了MaxFEs為10 000時各算法在三個不同精度(0.5、0.1和0.01)要求下的PR。圖6是各算法輸出結(jié)果與譜峰坐標(biāo)之間誤差的累積分布函數(shù)(Cumulative Distribution Function,CDF)圖。

表6 所有算法的運(yùn)行時間

表7 精度級別ε=0.5時各算法的峰識別率

表8 精度級別ε=0.1時各算法的峰識別率

表9 精度級別ε=0.01時各算法的峰識別率

圖6 各算法輸出結(jié)果誤差的CDF圖

表6～9的結(jié)果表明，r2pso、r2psolhc、r3pso和r3psolhc四種算法因為策略相似，所以在耗時及各個精度要求下的PR都接近；NPSO以外的其他算法中，F(xiàn)ER-PSO在0.5和0.1精度范圍內(nèi)下具有較高的PR，但所有算法在0.01精度條件下PR都非常低，只有NPSO算法在各個精度條件下均能保持高PR，這意味著其獲取到了更精確的譜峰坐標(biāo)值。圖6中CDF曲線越往左移表明DOA的估計誤差越小，其中NPSO的線條具有最大的斜率，表明絕大部分誤差(大約95%)都分布在極小的范圍內(nèi)，估值結(jié)果的方差小、數(shù)值穩(wěn)定；其余算法的曲線符合表7～9的數(shù)據(jù)，F(xiàn)ER-PSO在線條起始階段與其余算法重合，表明其低誤差結(jié)果的數(shù)量與其他算法相近，之后斜率增大，意味著較低誤差結(jié)果的數(shù)量會比其他算法多，而其他五種算法之間的誤差辨識度不高。

通過上述實驗可以看出，NPSO算法具有計算速度快的特點，相較于表2中0.01°網(wǎng)格算法所需的2 177 s,NPSO用時僅約其1/7 000。在保證高精度的前提下，新算法可大大提升譜峰搜索的速度。

3.3 算法在高計算資源下的性能

上節(jié)仿真實驗均運(yùn)行在低算力資源的環(huán)境中，且最大函數(shù)求值次數(shù)都選取10 000次，可以看出除NPSO外，在低計算資源條件下其余算法均表現(xiàn)出很低的峰識別率，這可能是因為算法受到算力的限制無法發(fā)揮出真實的性能。為此，本節(jié)實驗中放寬算力資源的限制，使算法能對解空間充分地進(jìn)行探索和挖掘，以此來檢驗算法的最大性能。選取表3樣例C，算法參數(shù)不變，最大函數(shù)求值次數(shù)取值從1 000次開始，之后為2 500次，然后每隔2 500次取一個數(shù)值，最大為50 000次。在每個MaxFEs取值下都對7種算法運(yùn)行500次蒙特卡洛仿真。圖7統(tǒng)計了結(jié)果的均方根誤差(Root Mean Square Error,RMSE)，它衡量了輸出結(jié)果(觀測值)與譜峰坐標(biāo)(真值)之間的偏差，其解析式為

(20)

式中:N表示進(jìn)行的蒙特卡洛仿真的輪數(shù)，S是入射信號總數(shù)，αn(i)表示第i輪仿真中算法輸出的第n個入射信號的俯仰角α與真值之間的誤差，βn(i)表示第i輪仿真中算法輸出的第n個信號的方位角β與真值之間的誤差。

觀察圖7可以看出，在較低算力資源下，NPSO較其余算法能取得更高的精度。隨著MaxFEs的增大，算法獲得了足夠的計算資源，RMSE變化逐漸穩(wěn)定，此時NPSO、FER-PSO和NCDE的RMSE相較其他算法能領(lǐng)先一個數(shù)量級。實驗證明，無論算力資源是否充足，NPSO算法較其他算法均能保持精度的優(yōu)勢。

圖7 RMSE隨函數(shù)求值次數(shù)變化曲線

3.4 SNR對算法性能的影響

為了檢驗算法在不同信噪比(Signal-to-Noise Ratio,SNR)下的性能，SNR分別設(shè)置為-20 dB、-10 dB、0 dB、10 dB和20 dB，其余參數(shù)不變，運(yùn)行500次蒙特卡洛仿真。圖8給出了8種方法的RMSE隨SNR變化曲線。當(dāng)SNR較低時空間譜會出現(xiàn)許多假峰、干擾峰，且峰值的坐標(biāo)與DOA之間會存在較大誤差，因此為了準(zhǔn)確比較結(jié)果，圖中給出了網(wǎng)格搜索的曲線，與該線重合度越高則說明結(jié)果越準(zhǔn)確。

圖8 RMSE隨SNR變化曲線

從圖8中可以看出，NPSO的曲線與網(wǎng)格搜索的曲線高度重合，表明無論是在信噪比高還是低的情況下，本文算法不受“假”峰、“干擾”峰等的影響，均能準(zhǔn)確地識別出譜函數(shù)的真正譜峰，有助于保持空MUSIC算法在低信噪比條件下性能維持穩(wěn)定。

3.5 不同高度譜峰對算法性能的影響

在實際應(yīng)用中，各信號的強(qiáng)度存在強(qiáng)弱之分，由此產(chǎn)生的譜峰高度就會有較大反差。為了檢驗算法在信號強(qiáng)度反差較大情況下的性能，本節(jié)仿真在樣例C的基礎(chǔ)上，PMUSIC采用功率表示，其余參數(shù)不變，得到的空間譜圖形如圖9所示。圖中四個信號的功率強(qiáng)度反差較大，表現(xiàn)為譜峰高度大小不一。

圖9 不同譜峰高度的空間譜圖形

實驗使用表5的算法參數(shù)，MaxFEs設(shè)置為10 000次，運(yùn)行500次蒙特卡洛仿真。統(tǒng)計算法在三個精度級別下的PR，結(jié)果如表10所示。

表10 不同精度級別時各算法的峰識別率

仿真結(jié)果顯示，在信號強(qiáng)度反差較大情況下，本文提出的NPSO不受譜峰之間落差大的影響，依舊能夠保持極高的精度。

4 結(jié)束語

譜峰搜索耗時長、實時性低是制約MUSIC算法應(yīng)用因素之一。在單譜峰中，研究者將進(jìn)化計算算法應(yīng)用于該領(lǐng)域已取得不錯的效果，但是面對多譜峰尚未有十分滿意的解決方案。本文提出了一種改進(jìn)的小生境粒子群算法(NPSO)，并將其應(yīng)用于MUSIC算法的多信號DOA估計中，通過與同類小生境算法和網(wǎng)格搜索進(jìn)行比較，得到如下研究結(jié)論：

(1)精度高。NPSO能選擇不同的搜索策略，在提高搜索廣度的同時，也兼顧了搜索精度。實驗結(jié)果顯示，在低算力資源情況下可以取得較高精度，在算力資源豐富的情況下精度能達(dá)到10-3，滿足了高精度應(yīng)用場景的需求，并且不易受低信噪比以及信號強(qiáng)度反差大的影響。

(2)速度快。小生境的并行性使得在10-2精度級別下，NPSO僅需網(wǎng)格搜索用時的1/7 000，大大提高了DOA估計速度，對提升通信、雷達(dá)等的性能具有十分重要的意義。

(3)本文所提出的小生境粒子群算法NPSO適用于求解多峰函數(shù)優(yōu)化問題，不僅僅是應(yīng)用于多信號的MUSIC算法譜峰搜索當(dāng)中，也可推廣到其他具有相似問題的工程中，具有潛在的應(yīng)用價值。

(4)本文僅僅提升了譜峰搜索性能，并未優(yōu)化MUSIC算法的其他環(huán)節(jié)，而DOA估計的整體誤差由所有環(huán)節(jié)的誤差組成。因此，下一步的研究內(nèi)容是將PSO用于提升MUSIC算法在接收數(shù)據(jù)處理等環(huán)節(jié)的性能。

綜上所述，本文提出的方法不僅可以很好地應(yīng)用于MUSIC算法的多信號DOA估計中，而且也可推廣到其他需要求解多峰優(yōu)化問題的工程應(yīng)用中。