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

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

0 引言

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

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

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

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

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

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

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

p(xn)=ejk·xn

(1)

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

(2)

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

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

(3)

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

(4)

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

(5)

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

(6)

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

[xlm,ylm]T=[rlcosφl(shuí)m,rlsinφl(shuí)m)]T

(7)

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

(8)

定義波束旁瓣區(qū)域?yàn)棣?，則最大旁瓣級(jí)峰值(PPSLL)是旁瓣區(qū)域的最大值為

(9)

波束主瓣寬度(PMLW)為

(10)

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

(11)

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

(12)

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

J1(kdmnsinφmin)sinφmin

(13)

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

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

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

由式(8)可知，fsinφ可看做原始波束圖的綜合自變量，而sinφ與掃描區(qū)域的范圍相關(guān)。綜合自變量不變時(shí)，原始波束圖不變，由此可以推斷，低頻時(shí)掃描區(qū)域的原始波束圖與高頻時(shí)縮小掃描區(qū)域的原始波束圖相似。因此，可將寬帶范圍最低旁瓣水平約束等效為最高頻率點(diǎn)局部掃描區(qū)域最低旁瓣水平約束。以需求頻率范圍上限確定目標(biāo)函數(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分別為兩個(gè)函數(shù)值的權(quán)重。

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

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

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

圖2 半徑余量示意圖

(15)

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

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

(16)

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

(i=1,2,…,I)

(18)

(i=1,2,…,I)

(19)

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

(20)

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

(21)

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

(22)

(23)

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

(24)

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

圖3 粒子群算法流程圖

3 仿真實(shí)驗(yàn)

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

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

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

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

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

3.2 兩種適應(yīng)度函數(shù)對(duì)比實(shí)驗(yàn)

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

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

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

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

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

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

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

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

表4 兩種陣型不同頻率的最高旁瓣級(jí)對(duì)比

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

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

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

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

3.4 兩種優(yōu)化方式對(duì)比實(shí)驗(yàn)

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

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

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

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

4 結(jié)束語(yǔ)

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