基于粒子群優(yōu)化的寬帶同心圓環(huán)陣型優(yōu)化方法

康雅聰，魏明洋，田巳睿，丁林寧

(1.南京理工大學(xué) 電子工程與光電技術(shù)學(xué)院 江蘇 南京 210094；2.中國科學(xué)院 聲學(xué)研究所，北京 100190)

0 引言

聲學(xué)成像技術(shù)在航空氣動聲學(xué)和機械設(shè)備噪聲檢測[1-2]等領(lǐng)域具有廣泛的應(yīng)用。它結(jié)合了攝像頭獲取的實時畫面和平面?zhèn)髀暺麝嚵械玫降穆曉茨芰枯椛浞植迹苤庇^地確定目標聲源的位置。陣列點擴散函數(shù)(PSF)[3-4]是指聲源位于畫面正中心時的陣列波束響應(yīng)，它指示了陣列對空間聲源的響應(yīng)情況，且與陣元位置有關(guān)。此外，衡量聲成像性能的一個重要指標是傳聲器陣列的抗干擾能力，即期望陣列處理后形成的波束圖具有很低的旁瓣水平。自適應(yīng)波束形成、反卷積波束形成等陣列算法能有效消除波束圖中的旁瓣，但是這些算法的效果受限于原始波束圖的質(zhì)量，即點擴散函數(shù)的旁瓣水平。

在實際工程應(yīng)用中，陣列的布放范圍、陣元數(shù)量、陣元間距等都有一定的要求。為了在這些實際限定的因素下得到更低的旁瓣水平，需要對陣列的陣元分布進行優(yōu)化設(shè)計。陣型優(yōu)化是指以陣元空間位置為自變量，點擴散函數(shù)的旁瓣水平為因變量，獲取最低旁瓣水平時陣元位置的優(yōu)化過程?？紤]到陣型優(yōu)化是一個高維度、非線性、不連續(xù)、不可微的優(yōu)化問題，當前的解決方案是利用全局優(yōu)化的智能算法進行求解，如遺傳算法[5]、粒子群算法[6]等。同心圓環(huán)陣列是一簇以坐標原點為圓心，不同圓半徑的圓環(huán)陣組成的陣列。其圓周對稱特性使陣列的方向圖形狀相對穩(wěn)定，且在俯仰角方向具有理想的輻射特性[7]，因而受到越來越多研究者的關(guān)注。具體的相關(guān)研究包括但不限于：Haupt[8]利用混合遺傳算法，以圓環(huán)半徑和圓環(huán)上陣元數(shù)量作為優(yōu)化對象對同心圓環(huán)陣列進行綜合，取得了較好的效果。陳客松等[9]通過約束同一圓環(huán)上陣元間距相等,利用修正遺傳算法優(yōu)化圓環(huán)的半徑,獲得最小的峰值旁瓣電平。欒曉明等[7]提出了一種基于差分進化算法的稀布同心圓環(huán)陣列半徑和陣元間距的聯(lián)合優(yōu)化方法，獲得了陣列響應(yīng)低峰值旁瓣電平的陣元分布。這些研究都是基于單一頻率的研究，目前關(guān)于寬帶的陣列優(yōu)化方法的研究還很少。

本文提出了一種基于平面同心圓環(huán)陣列的陣型優(yōu)化設(shè)計方法，構(gòu)造了一種度量寬帶點擴散函數(shù)旁瓣水平的目標函數(shù)，采用粒子群全局優(yōu)化算法，有效降低了波束掃描區(qū)域的整體旁瓣水平。同時，目標函數(shù)的低維度流形表示以及粒子群的高效搜索使得優(yōu)化算法快速收斂，有效地降低了優(yōu)化設(shè)計耗時。

1 同心圓環(huán)陣列波束掃描模型

如圖1所示，考慮一個由不同半徑的圓環(huán)陣列構(gòu)成的同心圓環(huán)傳聲器陣列，以同心圓環(huán)的圓心建立三維坐標系。陣列陣元位于xOy平面內(nèi)，假設(shè)同心圓環(huán)陣列共有L個圓環(huán)，每個圓環(huán)的半徑分別為r1、r2、…、rL、每個環(huán)上的陣元數(shù)量分別為H1、H2、…、HL，總的傳聲器個數(shù)為N。

圖1 同心圓環(huán)陣列掃描模型

假設(shè)聲源分布在距傳聲器陣列中心點為z0的掃描平面上，聲源滿足陣列遠場條件，則第n個傳聲器接收到掃描點(θ,φ)處的聲壓信號表示為

p(xn)=ejk·xn

(1)

式中：xn=[xn,yn](n=1,2,…,N)是第n個陣元的位置向量；k為波導(dǎo)矢量，即

(2)

式中：c為聲音在介質(zhì)中的傳播速度；f為入射聲波頻率；θ為入射波方位角；φ為入射波俯仰角。

假設(shè)掃描平面上有一單位幅度的聲源位于(0,0,z0)，其波矢量為k0，則傳聲器陣列接收到的聲壓為

(3)

一個掃描點沿著方向k的導(dǎo)向矢量為

(4)

則基于常規(guī)波束形成得到的原始波束圖為

(5)

由于傳聲器陣列位于z=0平面，則式(5)改寫為

(6)

具體到同心圓環(huán)陣列，陣元的位置坐標為

[xlm,ylm]T=[rlcosφlm,rlsinφlm)]T

(7)

式中：rl為第l個圓環(huán)的半徑；φlm表示第l個圓環(huán)上的第m個陣元與x軸的夾角。則同心圓環(huán)陣列的掃描聲源的原始波束圖為

(8)

定義波束旁瓣區(qū)域為Θ，則最大旁瓣級峰值(PPSLL)是旁瓣區(qū)域的最大值為

(9)

波束主瓣寬度(PMLW)為

(10)

式中D為圓環(huán)陣列最大直徑。文獻[10]定義了以下表示旁瓣能量水平的函數(shù)(即M函數(shù))：

(11)

M函數(shù)反映了圖1中環(huán)形掃描區(qū)域輻射聲源能量水平，Ω包含所有感興趣的掃描方向，k=|k|，φmin和φmax分別代表環(huán)形掃描區(qū)域的最小和最大俯仰角。為了最小化旁瓣，選擇的區(qū)域應(yīng)包含所有感興趣的掃描方向但不包含主瓣，故式(11)改寫為

(12)

Amn=J1(kdmnsinφmax)sinφmax-

J1(kdmnsinφmin)sinφmin

(13)

式中：dmn是第m個陣元與第n個陣元之間的距離；J1是第一類的一階貝塞爾函數(shù)。

2 基于粒子群算法的陣型優(yōu)化

2.1 構(gòu)造適應(yīng)度函數(shù)

由式(8)可知，fsinφ可看做原始波束圖的綜合自變量，而sinφ與掃描區(qū)域的范圍相關(guān)。綜合自變量不變時，原始波束圖不變，由此可以推斷，低頻時掃描區(qū)域的原始波束圖與高頻時縮小掃描區(qū)域的原始波束圖相似。因此，可將寬帶范圍最低旁瓣水平約束等效為最高頻率點局部掃描區(qū)域最低旁瓣水平約束。以需求頻率范圍上限確定目標函數(shù)，構(gòu)造如下適應(yīng)度函數(shù)：

J(x1…,xN)=w1M1(x1…,xN)+

w2M2(x1…,xN)

(14)

式中：M1(x1，…，xN)，M2(x1，…，xN)分別為不同積分區(qū)域內(nèi)的M函數(shù)；w1,w2分別為兩個函數(shù)值的權(quán)重。

另外，由于第一旁瓣存在于整個設(shè)計頻帶范圍內(nèi)，我們的設(shè)計目標是在減小第一旁瓣的同時最小化整體旁瓣水平，因此需要確定合適的權(quán)重以及M1、M2函數(shù)合適的積分邊界φmin,1、φmax,1和φmin,2、φmax,2。對于M1函數(shù)，φmin,1為主瓣第一次過零處的位置，φmax,1為第一旁瓣第一次過零處的位置；對于M2函數(shù)，φmin,2為M1函數(shù)所確定的φmax,1，φmax,2為實際掃描區(qū)域的最大俯仰角。

2.2 陣形設(shè)計與優(yōu)化

實際工程應(yīng)用中，陣列的布局區(qū)域有一定限制。為了降低優(yōu)化變量維度，本文假設(shè)環(huán)狀區(qū)域中內(nèi)圈圓半徑為rmin，外圈圓半徑為rmax，在約束區(qū)域內(nèi)設(shè)計L個圓環(huán)組成同心圓環(huán)陣列，各圓環(huán)的半徑分別為r1、r2、…、rL，每個環(huán)上的陣元數(shù)目分別為H1、H2、…、HL，并且每個環(huán)上的陣元等角度分布，區(qū)域內(nèi)陣元數(shù)量為N，約束所有陣元的最小間距為dmin。由于陣元的最小間距限制，則相鄰圓環(huán)之間的間隔至少為dmin。故采取文獻[11]中提出的余量思想：選取各個圓環(huán)的半徑余量Δr1、Δr2、…、ΔrL-1作為半徑的優(yōu)化參量以減少無效空間的搜索。半徑余量為相鄰兩個圓環(huán)之間除最小間距dmin外的其他可優(yōu)化的間距，如圖2所示。圖中陰影部分表示半徑余量，空白部分表示最小陣元間距,則有：

圖2 半徑余量示意圖

(15)

式中：rl為第l個圓環(huán)的半徑；Δrl為第l環(huán)和第l+1環(huán)之間的可優(yōu)化距離。

定義每個圓環(huán)上第一個陣元與水平方向的偏轉(zhuǎn)角φ1、φ2、…、φL，并將其作為待優(yōu)化變量。第m個陣元是指坐標系中從x軸起始，沿逆時針方向第m個陣元，則第l個圓環(huán)上的第m個陣元的偏轉(zhuǎn)角可表示為

(16)

綜合以上約束，對于粒子群優(yōu)化算法，定義待優(yōu)化粒子為X=(Δr1,Δr2，…，ΔrL-1,φ1,φ2…φL)，構(gòu)建以下粒子群優(yōu)化模型：

minJ(X)

s.t. J(X)=w1M1(X)+w2M2(X)

0≤Δrl≤rmax-rmin-(L-1)dmin

(17)

粒子群優(yōu)化算法(PSO)過程中每個優(yōu)化問題的解都是搜索空間中的一只鳥，稱為“粒子”。所有的粒子都有一個由被優(yōu)化的適應(yīng)度函數(shù)決定的適應(yīng)度值(Fitness Value)，每個粒子根據(jù)飛行速度決定其飛行的方向和距離，群體中的所有粒子追隨當前的最優(yōu)粒子在解空間中搜索。對于本文的陣型優(yōu)化，粒子的解就是能代表陣元位置的一組參數(shù)向量Xopt,適應(yīng)度函數(shù)就是2.1節(jié)所提出的函數(shù)。本文需要根據(jù)選取的(2L-1)個待優(yōu)化參數(shù)布置陣元，用一個(2L-1)維的向量表示一個粒子，隨機生成I個粒子，其中第t次迭代第i個粒子的位置Xi和速度vi分別表示為

(i=1,2,…,I)

(18)

(i=1,2,…,I)

(19)

第i個粒子當前t次迭代為止搜索到的最優(yōu)布陣方式即為個體極值，記為

(20)

整個粒子群當前t次迭代為止搜索到的最優(yōu)布陣方式為全局極值，記為

(21)

根據(jù)下式更新粒子的速度和位置：

(22)

(23)

式中：w為慣性權(quán)重；c1,c2為加速常數(shù)；s1,s2為兩個在[0,1]范圍內(nèi)變化的隨機數(shù)。本文w采用隨迭代次數(shù)不斷減小的慣性權(quán)重：

(24)

式中：wmax為最大的慣性權(quán)重；wmin為最小的慣性權(quán)重；T為總迭代次數(shù)；t為當前迭代次數(shù)。當w較大時，粒子搜索傾向于全局搜索；當w較小時，粒子可以在局部區(qū)域獲得更精確的極值。因此，采用隨迭代次數(shù)不斷減少的慣性權(quán)重使得搜索過程首先進行全局搜索，然后在逐步減小的搜索區(qū)域進行細化搜索以獲得更精確的解。整個算法流程如圖3所示。

圖3 粒子群算法流程圖

3 仿真實驗

3.1 寬帶轉(zhuǎn)最高頻點可行性分析實驗

為驗證2.1節(jié)提出的將寬帶范圍最低旁瓣水平約束等效為最高頻率點局部掃描區(qū)域最低旁瓣水平約束方法的可行性，本文設(shè)計的同心圓環(huán)陣參數(shù)如表1所示。每個環(huán)的半徑相差0.007 2 m，陣列平面距離掃描平面2 m，掃描平面為2.4 m×2.4 m的方形網(wǎng)格點，聲源頻率范圍為8～24 kHz，得到的原始波束圖沿x軸的頻率切片圖如圖4所示。

表1 同心圓環(huán)參數(shù)

圖4 原始波束圖沿x軸的頻率切片圖

由圖4可以看出，隨著頻率的增大，主瓣寬度逐漸減小，旁瓣向主瓣靠攏，頻率越高，則出現(xiàn)在視野范圍內(nèi)的旁瓣越多。因此，低頻時掃描區(qū)域的原始波束圖與高頻時縮小掃描區(qū)域的原始波束圖相似的推論是正確的。

3.2 兩種適應(yīng)度函數(shù)對比實驗

利用粒子群工具包，分別以本文所提出的適應(yīng)度函數(shù)和最高旁瓣級(PSLL)為優(yōu)化目標，對同心圓環(huán)陣列進行優(yōu)化設(shè)計，同心圓環(huán)陣的相關(guān)參數(shù)如表2所示。粒子群算法的參數(shù)如表3所示。兩種適應(yīng)度函數(shù)都在相同條件下進行仿真，陣列平面距離掃描平面2 m，掃描平面為2.4 m×2.4 m的方形網(wǎng)格點，網(wǎng)格分辨率為0.01 m。根據(jù)環(huán)形的尺寸可以求得主瓣邊界對應(yīng)的俯仰角φmin,M1選擇最大俯仰角為8°，M2選擇最小俯仰角為8°，最大俯仰角為33°。設(shè)w1=5，w2=0.1，聲源頻率為48 kHz，聲源強度為30 dB。

表2 同心圓環(huán)參數(shù)

表3 粒子群算法參數(shù)

針對兩種適應(yīng)度函數(shù)分別進行10次優(yōu)化求解實驗，取仿真中最好的結(jié)果，算法最優(yōu)值隨迭代次數(shù)的變化曲線如圖5所示。圖中黑色曲線表示全局最優(yōu)值隨迭代次數(shù)的變化情況，藍色點線代表種群的平均適應(yīng)度值隨迭代次數(shù)的變化情況，該曲線很好地反映了粒子群算法根據(jù)全局最優(yōu)值和局部最優(yōu)值不斷迭代，并最終收斂到最優(yōu)解的過程。

圖5 兩種適應(yīng)度函數(shù)迭代對比圖

將本文提出的適應(yīng)度函數(shù)稱為第一種適應(yīng)度函數(shù)，最高旁瓣級函數(shù)(PSLL)稱為第二種適應(yīng)度函數(shù)。圖5(a)是第一種適應(yīng)度函數(shù)的迭代曲線，由圖5(a)可看出迭代約200次后開始收斂。圖5(b)是第二種適應(yīng)度的迭代曲線，迭代300次開始收斂。對于第一種適應(yīng)度函數(shù)，平均每次迭代耗時1.436 1 s；而第二種適應(yīng)度函數(shù)平均每次迭代耗時17.876 s。這是由于第一種適應(yīng)度函數(shù)的計算量只與傳聲器的數(shù)量有關(guān)，第二種適應(yīng)度的計算量與傳聲器的數(shù)量、掃描的網(wǎng)格點數(shù)有關(guān)，而實際應(yīng)用中往往需要高精度的波束形成圖，掃描的網(wǎng)格點數(shù)遠大于傳聲器數(shù)量，故第一種適應(yīng)度函數(shù)平均每次迭代的運算量遠小于第二種適應(yīng)度函數(shù)。由此可見，本文所提出的適應(yīng)度函數(shù)對寬帶信號可以有效降低旁瓣，說明優(yōu)化方案可行。選取本文提出適應(yīng)度函數(shù)進行尋優(yōu)可以減少運算量，提高計算速度。

這兩種適應(yīng)度函數(shù)得到的優(yōu)化陣型和原始波束形成圖如圖6所示。表4為兩種陣型的參數(shù)對比。

圖6 兩種適應(yīng)度函數(shù)得到的優(yōu)化陣型比較

表4 兩種陣型不同頻率的最高旁瓣級對比

通過圖6(b)、(d)和表4可以看出，以第一種適應(yīng)度函數(shù)優(yōu)化得到的陣列一和以第二種適應(yīng)度函數(shù)優(yōu)化得到的陣列二擁有相似的主瓣寬度，48 kHz時均為0.32 m。在48 kHz時，雖然陣列二的最高旁瓣級為-19.55 dB，低于陣列一最高旁瓣級，但是該最高旁瓣出現(xiàn)在離主瓣較近的第一旁瓣處，而陣列一的最高旁瓣出現(xiàn)在遠離主瓣視野的四角處。另外，陣列一的第一旁瓣級比陣列二低10.71 dB，這表明第一種適應(yīng)度函數(shù)能有效降低主瓣附近的第一旁瓣高度。由表4可見，陣列一低頻部分的最高旁瓣級總是低于陣列二，這說明對于寬頻波束圖，陣列一的整體旁瓣水平低于陣列二。以上對比分析證實了本文適應(yīng)度函數(shù)的有效性。

3.3 權(quán)重系數(shù)對第一種適應(yīng)度函數(shù)優(yōu)化結(jié)果的影響

設(shè)置不同的權(quán)重系數(shù)會對實驗結(jié)果產(chǎn)生不同的影響。為了說明實驗第3.2節(jié)設(shè)置的權(quán)重的必要性，本文設(shè)計了另外兩組不同權(quán)重的仿真實驗，其他條件同仿真實驗第3.2節(jié)，得到的陣型原始波束圖如圖7所示。對比圖7(a)和圖6(b)可以看出，相比于兩個權(quán)重系數(shù)相等的情況，w1大于w2時可使靠近主瓣的第一旁瓣級更低，而其他范圍的旁瓣水平相似；對比圖7(b)和圖6(a)可以看出，相比于兩個權(quán)重系數(shù)相等的情況下，w1

圖7 兩組權(quán)重得到的原始波束圖

3.4 兩種優(yōu)化方式對比實驗

為了說明本文所提出的加入偏轉(zhuǎn)角的優(yōu)化方式能有效降低陣列在寬頻信號的旁瓣水平，本節(jié)設(shè)計了僅優(yōu)化半徑的仿真實驗。采用本文提出的適應(yīng)度函數(shù)，其他條件同仿真實驗第3.2節(jié)，得到僅優(yōu)化半徑的優(yōu)化陣型和原始波束圖如圖8所示，所得陣型參數(shù)與前文陣型一對比如表5所示。

圖8 僅優(yōu)化半徑得到的陣型

表5 兩種陣型參數(shù)對比

由圖8和表5可以看出，對于僅優(yōu)化半徑得到的陣列波束形成圖，其主瓣寬度和整體旁瓣能量水平都遠不及優(yōu)化半徑聯(lián)合偏轉(zhuǎn)角度的優(yōu)化方式，這表明本文提出優(yōu)化方式的有效性。

4 結(jié)束語

針對實際工程中約束陣列陣元總數(shù)、最大孔徑、最小陣元間距等條件下的陣形優(yōu)化設(shè)計問題，提出了一種基于粒子群優(yōu)化算法且適用于寬帶信號的同心圓環(huán)陣列陣型優(yōu)化方法。該方法以約束的最高頻率構(gòu)造了反映寬帶信號在掃描區(qū)域內(nèi)旁瓣水平的適應(yīng)度函數(shù)，同時以圓環(huán)半徑和陣元偏轉(zhuǎn)角作為聯(lián)合優(yōu)化變量，基于粒子群優(yōu)化算法對陣型優(yōu)化問題進行求解。數(shù)值結(jié)果仿真表明，本文所提出的適應(yīng)度函數(shù)在優(yōu)化效率上優(yōu)于傳統(tǒng)的以最高旁瓣級作為適應(yīng)度函數(shù)的方法，優(yōu)化后的陣型在掃描區(qū)域內(nèi)旁瓣能量水平更低，更適合于寬帶信號場景。下一步的工作重點為聲像儀信號采集系統(tǒng)的實現(xiàn)，將本文設(shè)計的同心圓環(huán)傳聲器陣型應(yīng)用到該系統(tǒng)中，對優(yōu)化后陣型的性能進行實驗驗證。