基于粒子群優化的寬帶同心圓環陣型優化方法

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

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

0 引言

聲學成像技術在航空氣動聲學和機械設備噪聲檢測[1-2]等領域具有廣泛的應用。它結合了攝像頭獲取的實時畫面和平面傳聲器陣列得到的聲源能量輻射分布，能直觀地確定目標聲源的位置。陣列點擴散函數(PSF)[3-4]是指聲源位于畫面正中心時的陣列波束響應，它指示了陣列對空間聲源的響應情況，且與陣元位置有關。此外，衡量聲成像性能的一個重要指標是傳聲器陣列的抗干擾能力，即期望陣列處理后形成的波束圖具有很低的旁瓣水平。自適應波束形成、反卷積波束形成等陣列算法能有效消除波束圖中的旁瓣，但是這些算法的效果受限于原始波束圖的質量，即點擴散函數的旁瓣水平。

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

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

1 同心圓環陣列波束掃描模型

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

圖1 同心圓環陣列掃描模型

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

p(xn)=ejk·xn

(1)

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

(2)

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

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

(3)

一個掃描點沿著方向k的導向矢量為

(4)

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

(5)

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

(6)

具體到同心圓環陣列，陣元的位置坐標為

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

(7)

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

(8)

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

(9)

波束主瓣寬度(PMLW)為

(10)

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

(11)

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

(12)

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

J1(kdmnsinφmin)sinφmin

(13)

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

2 基于粒子群算法的陣型優化

2.1 構造適應度函數

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

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

w2M2(x1…,xN)

(14)

式中：M1(x1，…，xN)，M2(x1，…，xN)分別為不同積分區域內的M函數；w1,w2分別為兩個函數值的權重。

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

2.2 陣形設計與優化

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

圖2 半徑余量示意圖

(15)

式中：rl為第l個圓環的半徑；Δrl為第l環和第l+1環之間的可優化距離。

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

(16)

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

minJ(X)

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

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

(17)

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

(i=1,2,…,I)

(18)

(i=1,2,…,I)

(19)

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

(20)

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

(21)

根據下式更新粒子的速度和位置：

(22)

(23)

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

(24)

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

圖3 粒子群算法流程圖

3 仿真實驗

3.1 寬帶轉最高頻點可行性分析實驗

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

表1 同心圓環參數

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

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

3.2 兩種適應度函數對比實驗

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

表2 同心圓環參數

表3 粒子群算法參數

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

圖5 兩種適應度函數迭代對比圖

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

這兩種適應度函數得到的優化陣型和原始波束形成圖如圖6所示。表4為兩種陣型的參數對比。

圖6 兩種適應度函數得到的優化陣型比較

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

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

3.3 權重系數對第一種適應度函數優化結果的影響

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

圖7 兩組權重得到的原始波束圖

3.4 兩種優化方式對比實驗

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

圖8 僅優化半徑得到的陣型

表5 兩種陣型參數對比

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

4 結束語

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