朱 任，莊銘杰

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

0 引 言

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

近年來，以遺傳算法(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ì)中，經(jīng)過仿真證實(shí)了這類算法能有效降低譜峰搜索的計(jì)算量，大幅提升DOA估計(jì)的精度和實(shí)時(shí)性，在實(shí)際工程應(yīng)用中具有潛在的應(yīng)用價(jià)值。

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

1 基于MUSIC算法的DOA估計(jì)

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

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

圖1 陣列模型

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

(1)

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

(2)

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

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

(3)

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

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

(4)

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

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

(5)

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

(6)

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

aH(α)UN=0。

(7)

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

(8)

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

(9)

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

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

(10)

當(dāng)有k個(gè)入射信號(hào)時(shí)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)的極大值點(diǎn)，就可以估算出達(dá)波信號(hào)的俯仰角α和方位角β。

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

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

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

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

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

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

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

以遺傳算法、粒子群算法為代表的經(jīng)典進(jìn)化算法在單信號(hào)DOA估計(jì)中表現(xiàn)出優(yōu)異的性能，但卻難以估計(jì)多信號(hào)的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°處的極大值點(diǎn)，而卻無法發(fā)現(xiàn)α=-60°和α=20°兩處的極大值點(diǎn)。所有粒子都趨向于最優(yōu)粒子，降低了粒子的多樣性，從而失去發(fā)現(xiàn)其他極大值點(diǎn)的機(jī)會(huì)。為了降低最優(yōu)粒子對(duì)全局粒子的吸引，本文引入小生境思想將解空間劃分成多個(gè)區(qū)域，避免PSO在多峰的譜峰搜索中出現(xiàn)上述問題。

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

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

2.1 中點(diǎn)順序聚類算法

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

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

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

(15)

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

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

(16)

中點(diǎn)順序聚類算法的偽代碼如下：

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

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

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

遍歷：從i=2到m

找到Pi最近的小生境nicheN，并計(jì)算中點(diǎn)midi。

如果：滿足式(15)：

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

否則：

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

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

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

(17)

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

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

(18)

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

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

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

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

循環(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對(duì)產(chǎn)生的粒子聚類，劃分小生境。

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

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

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

3 仿真與分析

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

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

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

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

(19)

為了確定NPSO中的預(yù)設(shè)年齡參數(shù)，選取粒子數(shù)為200，預(yù)設(shè)年齡為1～20，并選擇表3中的樣例C，運(yùn)行500次蒙特卡洛仿真并統(tǒng)計(jì)達(dá)到預(yù)設(shè)年齡時(shí)的PR，結(jié)果如表4所示。因?yàn)镹PSO在達(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將在多信號(hào)DOA估計(jì)中與它們進(jìn)行比較。

表5 算法參數(shù)

3.2 算法在低計(jì)算資源下的性能

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

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

表7 精度級(jí)別ε=0.5時(shí)各算法的峰識(shí)別率

表8 精度級(jí)別ε=0.1時(shí)各算法的峰識(shí)別率

表9 精度級(jí)別ε=0.01時(shí)各算法的峰識(shí)別率

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

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

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

3.3 算法在高計(jì)算資源下的性能

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

(20)

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

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

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

3.4 SNR對(duì)算法性能的影響

為了檢驗(yàn)算法在不同信噪比(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較低時(shí)空間譜會(huì)出現(xiàn)許多假峰、干擾峰，且峰值的坐標(biāo)與DOA之間會(huì)存在較大誤差，因此為了準(zhǔn)確比較結(jié)果，圖中給出了網(wǎng)格搜索的曲線，與該線重合度越高則說明結(jié)果越準(zhǔn)確。

圖8 RMSE隨SNR變化曲線

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

3.5 不同高度譜峰對(duì)算法性能的影響

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

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

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

表10 不同精度級(jí)別時(shí)各算法的峰識(shí)別率

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

4 結(jié)束語

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

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

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

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

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

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