一種基于相位平均的旋轉(zhuǎn)聲源高分辨率定位方法

初 寧， 黃 乾， 余 亮， 寧 岳， 徐建鋒， 吳大轉(zhuǎn)

(1. 浙江大學(xué) 能源工程學(xué)院， 杭州 310027； 2. 浙江上風(fēng)高科專(zhuān)風(fēng)實(shí)業(yè)有限公司， 浙江 紹興 312300;3. 上海交通大學(xué) 機(jī)械與動(dòng)力工程學(xué)院， 上海 200240)

風(fēng)機(jī)等通風(fēng)設(shè)備在工業(yè)中應(yīng)用廣泛，發(fā)揮著送風(fēng)、換氣和排煙等重要作用，但是設(shè)備在運(yùn)行過(guò)程中和潛在故障發(fā)生時(shí)都不可避免地產(chǎn)生噪聲。利用麥克風(fēng)陣列采集風(fēng)機(jī)聲源信號(hào)，對(duì)旋轉(zhuǎn)葉片進(jìn)行聲源定位、繪制聲源能量的空間分布云圖，能夠有效定位異常聲源的位置，對(duì)風(fēng)機(jī)的噪聲品質(zhì)評(píng)價(jià)、故障預(yù)警和結(jié)構(gòu)優(yōu)化提供新的技術(shù)支撐。

旋轉(zhuǎn)聲源定位的研究目前主要集中于提高算法分辨率和定位精度。Lowis等[1]將管道內(nèi)旋轉(zhuǎn)聲源的聲壓表示成旋轉(zhuǎn)框架下的頻域格林函數(shù)與聲源源強(qiáng)之間的乘積，再利用波束形成識(shí)別聲源，但分辨率有待提高。Oerlemans等[2-3]提出了通過(guò)時(shí)域插值法去除多普勒效應(yīng)的影響，從而實(shí)現(xiàn)旋轉(zhuǎn)聲源定位，但是該算法在插值運(yùn)算時(shí)破壞了原始采樣序列，會(huì)產(chǎn)生較大定位誤差。王梟等[4]通過(guò)在靜止聲源框架中加入轉(zhuǎn)速，推導(dǎo)得到DAMAS2 (extensions of DAMAS) 的修正算法，有效提升了單極子聲源的定位分辨率，但是定位較復(fù)雜聲源時(shí)存在一定的誤差。Zhang等[5]提出了一種基于時(shí)域等效源方法的時(shí)域反演技術(shù)，不僅可以準(zhǔn)確定位旋轉(zhuǎn)單極子，而且可以定量地預(yù)測(cè)聲源的輻射聲場(chǎng)特征，但是該方法僅適用于近場(chǎng)測(cè)試。Yu等[6]提出了一種降階循環(huán)維納濾波器，可以從轉(zhuǎn)子葉片的風(fēng)洞測(cè)試信號(hào)中提取出轉(zhuǎn)子葉片信號(hào)以及其他干擾，從而純化波束形成的定位結(jié)果。Ralph等[7-8]分別基于虛擬陣列旋轉(zhuǎn)法實(shí)現(xiàn)了旋轉(zhuǎn)風(fēng)機(jī)和四葉片風(fēng)扇的聲源定位。

為了提高旋轉(zhuǎn)聲源定位算法的分辨率，不少學(xué)者做過(guò)大量研究工作。由于常規(guī)波束形成給出了一個(gè)模糊的低分辨率定位成像結(jié)果，附帶有較多的旁瓣影響，所以可以考慮對(duì)常規(guī)波束形成的結(jié)果進(jìn)行解卷積得到高分辨率結(jié)果。經(jīng)典的解卷積方法有DAMAS(deconvolution approach for the mapping of acoustic sources)[9]，CLEAN[10]，DAMAS2[11]以及SC-DAMAS(sparsity constrained-DAMAS)[12]等。將解卷積方法應(yīng)用到旋轉(zhuǎn)聲源定位中，即可實(shí)現(xiàn)旋轉(zhuǎn)聲源高分辨率定位。莫品西[13]基于去多普勒效應(yīng)給出了旋轉(zhuǎn)聲源的高分辨率成像算法，通過(guò)直升機(jī)旋翼試驗(yàn)驗(yàn)證了方法的可靠性。許丹[14]提出了一種旋轉(zhuǎn)框架技術(shù)來(lái)消除多普勒效應(yīng)的影響，再采用改進(jìn)的DAMAS解卷積算法識(shí)別旋轉(zhuǎn)聲源，通過(guò)仿真和試驗(yàn)驗(yàn)證了該方法的有效性。Ma等[15]提出了一種通過(guò)小波壓縮計(jì)算網(wǎng)格來(lái)提高DAMAS效率的新方法，同時(shí)還能保證DAMAS算法的高分辨率。

但是，采用時(shí)域插值法對(duì)旋轉(zhuǎn)聲源進(jìn)行去多普勒效應(yīng)時(shí)，會(huì)破壞原始采樣序列的結(jié)構(gòu)和相位的連續(xù)性，從而在插值時(shí)產(chǎn)生無(wú)法預(yù)知的誤差，而且時(shí)域方法的分辨率普遍有待提高。而基于頻域去多普勒效應(yīng)的格林函數(shù)只在聲源運(yùn)動(dòng)過(guò)程中的某一個(gè)瞬時(shí)成立，因此只能適應(yīng)低速移動(dòng)目標(biāo)、遠(yuǎn)距離探測(cè)下旋轉(zhuǎn)聲源定位。

為此，本文提出一種基于相位平均法[16-17]的旋轉(zhuǎn)聲源定位方法。相位平均主要被應(yīng)用于消除干擾噪聲的影響[18-20]，針對(duì)干擾噪聲的統(tǒng)計(jì)平穩(wěn)特性，采用平均互功率譜消除功率譜中由干擾噪聲引起的成分。針對(duì)旋轉(zhuǎn)聲源定位時(shí)，本文借用相位平均法的思想將旋轉(zhuǎn)聲信號(hào)劃分成足夠多個(gè)時(shí)間序列切片，在每一切片內(nèi)，微小位移的旋轉(zhuǎn)聲源可近似等效為靜止聲源進(jìn)行定位，并給出了相位平均技術(shù)的應(yīng)用條件。為進(jìn)一步提高定位分辨率，運(yùn)用卷積模型改造基于波束形成的能量傳播模型，利用解卷積算法消除定位模糊和提高定位精度。其優(yōu)勢(shì)在于，卷積模型中的點(diǎn)擴(kuò)散函數(shù)(point spread function ，PSF)能夠建立起陣列測(cè)量與聲源旋轉(zhuǎn)的定量聯(lián)系，由此通過(guò)PSF來(lái)表征麥克風(fēng)陣列定位性能、分析頻率對(duì)定位分辨率的影響。最后，通過(guò)仿真信號(hào)、旋轉(zhuǎn)單極子試驗(yàn)、風(fēng)扇試驗(yàn)驗(yàn)證了方法的有效性和魯棒性。

1 理論分析

1.1 旋轉(zhuǎn)聲源信號(hào)傳播模型

如圖1所示，在介質(zhì)運(yùn)動(dòng)的聲場(chǎng)中，聲場(chǎng)聲壓滿足下面的非齊次聲學(xué)波動(dòng)方程

圖1 旋轉(zhuǎn)單極子在有流聲場(chǎng)下的傳播模型

(1)

根據(jù)聲場(chǎng)疊加原理，在某一位置x所測(cè)量的聲壓可以表示成所有等效源對(duì)于該處聲壓的貢獻(xiàn)值疊加，那么式(1)的解為

(2)

式中：dVy為位置y處聲源的體積微元；G(x,t|y,t′)為格林函數(shù),表示位置y處聲源在t′時(shí)刻發(fā)出的脈沖信號(hào)在接收端x處和t時(shí)刻的聲場(chǎng)響應(yīng)；δ(·)為狄拉克函數(shù);聲源按照運(yùn)動(dòng)軌跡xs(t′)旋轉(zhuǎn);s(t′)為聲源在t′時(shí)刻的聲強(qiáng)。

將式(2)中的3個(gè)等式聯(lián)立，再利用狄拉克函數(shù)的篩分性質(zhì)可以解得單極子旋轉(zhuǎn)聲源在自由聲場(chǎng)下的聲壓解為

(3)

在有流聲場(chǎng)中，運(yùn)動(dòng)聲源的頻域傳播模型表示為

(4)

式(4)給出了運(yùn)動(dòng)聲源在頻域上的連續(xù)傳播模型，基于此可以得出麥克風(fēng)陣列采集聲信號(hào)的離散傳播模型。假定麥克風(fēng)陣列由M個(gè)麥克風(fēng)組成，聲源平面被劃分成N個(gè)網(wǎng)格點(diǎn)，其中有K個(gè)不相關(guān)單極子聲源(K

P=GS+E

(5)

式中：P∈M×Z為M個(gè)通道麥克風(fēng)信號(hào)做Z點(diǎn)傅里葉變換后的頻域矩陣；S=[S1,S2,…,SK]T∈K×Z為K個(gè)聲源信號(hào)的頻域形式,(·)T為轉(zhuǎn)置；E∈M×Z為干擾噪聲；G∈M×K為頻域格林函數(shù)，由式(4)可以得出

(6)

1.2 相位平均旋轉(zhuǎn)聲源定位算法

式(5)給出了麥克風(fēng)陣列采集聲信號(hào)在頻域上的正向傳播模型，但是這是一個(gè)非線性系統(tǒng)，因?yàn)槁曉葱盘?hào)和聲源位置均是未知量，因此求解該反問(wèn)題時(shí)，在聲源所在平面劃出一塊掃描區(qū)域并劃分成多個(gè)網(wǎng)格點(diǎn)，每一個(gè)網(wǎng)格點(diǎn)可以看成一個(gè)等效聲源，這里假定有N個(gè)網(wǎng)格點(diǎn)，由此將非線性系統(tǒng)問(wèn)題的求解轉(zhuǎn)化成線性系統(tǒng)求解。常規(guī)波束形成(conventional beamforming, CBF)對(duì)于該線性反問(wèn)題的求解為

(7)

(8)

式中：f為需要分析的目標(biāo)頻率；rm,n為第m個(gè)麥克風(fēng)到第n個(gè)網(wǎng)格點(diǎn)的距離，n=1,2,…,N;m=1,2,…,M。RCSM(f)為麥克風(fēng)陣列頻域信號(hào)的互譜矩陣，具體表示為

RCSM(f)=[P(f)P(f)H]

(9)

式中：P(f)∈M×Z為陣列頻域信號(hào)；[·]為數(shù)學(xué)期望，由數(shù)理統(tǒng)計(jì)的原理可知，樣本較大時(shí)，均值近似等于數(shù)學(xué)期望，式(9)可近似表達(dá)為

(10)

式中，L為快拍數(shù)，表示將原始信號(hào)分成L段，相鄰兩個(gè)子段信號(hào)重疊率為50%。將式(10)的互譜矩陣估計(jì)值代入式(7)即可得到常規(guī)波束形成的輸出。常規(guī)波束形成的算法原理示意圖，如圖2上半部分[21]所示。對(duì)于靜止聲源的求解，CBF是一種魯棒性很高的方法，但是當(dāng)聲源旋轉(zhuǎn)時(shí)，由圖2可知CBF是對(duì)整個(gè)采樣序列求互譜矩陣，所以無(wú)法準(zhǔn)確給出第n個(gè)網(wǎng)格點(diǎn)的功率輸出。

而相位平均法則是將聲源旋轉(zhuǎn)一周的采樣序列劃分成足夠多個(gè)切片，每一個(gè)切片聲源位移微小，可以近似等效為靜止聲源求解。為了提高分辨率，可以利用多個(gè)周期內(nèi)聲源經(jīng)過(guò)同一弧度的數(shù)據(jù)進(jìn)行多次平均。相位平均的原理示意圖，如圖2下半部分所示。第z個(gè)切片下聲信號(hào)經(jīng)過(guò)相位平均處理的互譜矩陣表示為

r{z}CSM(ma,mb)=1QLp∑Q-1q=0[({pma(I{z}+qD)})× ({pmb(I{z}+qD)})H]I{z}=Im0,n0+(z-1)Lp+II= -Lp2,-Lp2+1,…,Lp-12 Im0,n0=round(fsΔtm0,n0)

(11)

-Lp2,

Lp-12;·

圖2 常規(guī)波束形成(上)與相位平均波束形成(下)的對(duì)比

(12)

式中:fr為聲源轉(zhuǎn)速，r/min;R為旋轉(zhuǎn)半徑，m，fs為采樣頻率；grid_div為掃描網(wǎng)格點(diǎn)的邊長(zhǎng)，m；δCBF為波束形成分辨率，m；c為聲速常量，m/s；r為陣列到聲源的平均距離,m;f為分析頻率，Hz；Da為陣列孔徑，m。式(12)給出了聲源轉(zhuǎn)速fr，r/min，采樣頻率fs,Hz,分析頻率f,Hz，旋轉(zhuǎn)半徑R,Hz，單個(gè)切片采樣長(zhǎng)度Lp，samples，這幾個(gè)主要變量之間的關(guān)系，可以看出這些變量之間是相互制約的。例如，取c=340 m/s,r=2 m,Da=1 m,f=2 500 Hz,fs=50 kHz,fr=1 200 r/min,R=0.5 m，求得Lp≤217。

式(12)第一個(gè)不等式限制了每個(gè)切片下聲源旋轉(zhuǎn)過(guò)的弧長(zhǎng)不超過(guò)一個(gè)掃描網(wǎng)格點(diǎn)，這是因?yàn)閷?duì)每一個(gè)切片內(nèi)的旋轉(zhuǎn)聲源定位時(shí)將其近似看作是靜止聲源處理，將劃分的每一個(gè)網(wǎng)格點(diǎn)看作是一個(gè)等效點(diǎn)聲源，因此為了滿足點(diǎn)聲源的物理分辨率，需要保證聲源在單個(gè)切片內(nèi)運(yùn)動(dòng)過(guò)的軌跡不超過(guò)一個(gè)網(wǎng)格點(diǎn)的大??；第二個(gè)不等式限制了掃描網(wǎng)格點(diǎn)邊長(zhǎng)不應(yīng)超過(guò)常規(guī)波束形成(conventional beamformer，CBF)算法的理論分辨率(即劃分的網(wǎng)格物理分辨率尺寸應(yīng)當(dāng)小于算法的理論分辨率)，這是因?yàn)樗惴ńo出了理論分辨率的極限值，即理論上可以分辨出的兩個(gè)聲源的最小距離，而我們假定劃分的每一個(gè)掃描網(wǎng)格點(diǎn)是一個(gè)等效點(diǎn)聲源，在不犧牲算法的理論分辨率要求下，應(yīng)當(dāng)限制單個(gè)掃描網(wǎng)格點(diǎn)的大小不超過(guò)該理論分辨率。

上述分析給出了其他變量給定條件下Lp取值的上限，需要注意的是，采樣長(zhǎng)度Lp過(guò)小會(huì)導(dǎo)致傅里葉變換的頻率分辨率過(guò)低，能量不夠集中；除此之外，采樣長(zhǎng)度過(guò)小還會(huì)導(dǎo)致樣本量太小，誤差較大。因此，應(yīng)用相位平均法時(shí)還需要對(duì)采樣長(zhǎng)度Lp的下界進(jìn)行約束，根據(jù)大量的仿真和試驗(yàn)研究，我們發(fā)現(xiàn)當(dāng)Lp<100時(shí)，能量幾乎無(wú)法聚焦，難以實(shí)現(xiàn)成像，所以本文所述的基于相位平均的旋轉(zhuǎn)聲源定位方法選取的單個(gè)切片采樣長(zhǎng)度均為100或者200。當(dāng)Lp=100時(shí)，可以通過(guò)式(12)確定另外幾個(gè)變量的取值范圍，即相位平均法的適用范圍。例如，取c=340 m/s,r=2 m,Da=1 m,f=2 500 Hz,fs=50 kHz,Lp=100 samples,R=0.5 m，求得fr<2 597 r/min。當(dāng)然，根據(jù)式(12)確定的只是轉(zhuǎn)速的理論上限，在實(shí)際情況下，轉(zhuǎn)速如果略高于理論上限，采用本方法也可以定位，只是會(huì)犧牲一定的空間分辨率。同理，當(dāng)聲源轉(zhuǎn)速fr給定后，也可以通過(guò)式(12)確定出分析頻率、聲源旋轉(zhuǎn)半徑等其他變量的適用條件，在此不一一贅述。

將式(11)基于相位平均的互譜矩陣代入式(7)即可獲得切片z中，第n個(gè)網(wǎng)格點(diǎn)基于相位平均處理的波束形成輸出功率，具體表示

(13)

為了提高空間分辨率，Capon在波束形成基礎(chǔ)上，提出了最小方差無(wú)失真響應(yīng)方位估計(jì)波束器(minimum variance distortionless response，MVDR)，極大地抑制了波束形成的旁瓣[21]，基于相位平均處理的MVDR算法的方向向量為

(14)

將式(14)代入式(13)可得旋轉(zhuǎn)框架下的MVDR空間功率譜表達(dá)式

(15)

式中，(·)-1為矩陣求逆。為了區(qū)別常規(guī)MVDR算法，本文中將旋轉(zhuǎn)框架下基于相位平均處理的MVDR算法稱(chēng)作RMVDR(rotating MVDR)。

1.3 點(diǎn)擴(kuò)散函數(shù)和RDAMAS算法

假定各個(gè)單極子聲源信號(hào)互不相關(guān)，將式(5)代入式(9)可得

RCSM=G[SSH]GH+

(16)

式中，xn為[SSH]的對(duì)角線一個(gè)元素，其代表聲源能量，將式(16)代入式(7)可得波束形成輸出為

(17)

式中，εn1為由噪聲引起的誤差項(xiàng)。由1.2節(jié)中的分析可知，相位平均可以將每一個(gè)切片內(nèi)產(chǎn)生微小位移的旋轉(zhuǎn)聲源近似等效為靜止聲源求解。基于以上近似，各個(gè)切片內(nèi)信號(hào)進(jìn)行處理時(shí)不考慮多普勒效應(yīng)的影響，運(yùn)動(dòng)馬赫數(shù)Mam,k=0，相當(dāng)于式(6)中的gm,k可以由式(8)來(lái)表示，即gn等價(jià)于en，則式(17)可以表示為

(18)

式中：n1,n2=1,…,N；yn1為第n1個(gè)網(wǎng)格點(diǎn)的波束形成輸出；en2為第n2個(gè)網(wǎng)格點(diǎn)的波束形成方向向量。將式(18)表示為矩陣形式

(19)

式中：ε=Diag([EEH])∈N×1為模型誤差；Diag(·)為將矩陣對(duì)角線元素以向量形式表示；H∈N×N稱(chēng)之為傳播矩陣；x∈N×1為掃描區(qū)域的所有等效聲源。

利用式(19)求解x可以看成一個(gè)解卷積問(wèn)題，在這里采用DAMAS算法求解該問(wèn)題。用DAMAS求解時(shí)，假定ε=0，y可以看作是基于相位平均處理后的波束形成輸出或者M(jìn)VDR輸出，即式(13)或式(15)的結(jié)果，采用高斯-賽德?tīng)柕ㄇ蠼?，迭代步驟

(20)

式(19)中，由于傳播矩陣H是近似的對(duì)稱(chēng)塊托普利茲(symmetric block toeplitz，SBT)矩陣，可以通過(guò)這個(gè)矩陣獲得一個(gè)不變卷積核[22]，將Hx以二維卷積形式表示為

(21)

式中：x0為向量x的矩陣形式，此處即為掃描區(qū)域的網(wǎng)格點(diǎn)矩陣;h*∈A×B為不變卷積核，也叫做點(diǎn)擴(kuò)散函數(shù)(point spread function, PSF)；*為同維卷積操作，即卷積結(jié)果依然是一個(gè)A×B矩陣；[·]′為將矩陣按照列轉(zhuǎn)化為一個(gè)列向量。

按照式(19)可以將傳播矩陣H中的元素hn1,n2表示

(22)

(23)

通常，式(21)中的PSF尺寸選為A=B=Nm=min {row(x0),column(x0)}，即取掃描網(wǎng)格矩陣的行數(shù)和列數(shù)中較小的值，PSF表示

h*a,b=h～n1,n2n1=N+12n2=n1+Nm+12-a Nm+Nm+12-b

(24)

2 仿真模擬

2.1 仿真條件

如圖1(a)所示，聲源掃描平面設(shè)置為2 m×2 m的正方形平面，掃描網(wǎng)格大小0.05 m×0.05 m，設(shè)置兩個(gè)仿真聲源，聲源的初始坐標(biāo)為(-0.7,0)，(0.7,0)，旋轉(zhuǎn)半徑為0.7 m，圓心為坐標(biāo)原點(diǎn)；右側(cè)的色標(biāo)顯示了聲源的聲壓級(jí)大小，單位為dB，動(dòng)態(tài)顯示范圍設(shè)置為16 dB。仿真陣列采用56通道螺旋陣列，示意圖如圖1(b)所示。陣列參數(shù)如表1所示。陣列平面與聲源掃描平面相距1 m。聲源轉(zhuǎn)速設(shè)置為10 Hz，逆時(shí)針?lè)较蛐D(zhuǎn)，為了模擬實(shí)際旋轉(zhuǎn)工況，仿真在10 Hz的基本轉(zhuǎn)速上又添加了±1 Hz(即±10%)的頻率波動(dòng)作為隨機(jī)干擾。采樣頻率設(shè)置為50 kHz，單個(gè)切片的采樣點(diǎn)數(shù)為100個(gè)采樣點(diǎn)，RDAMAS迭代次數(shù)300次，均選取RBF的結(jié)果作為RDAMAS的輸入。

表1 麥克風(fēng)陣列參數(shù)

2.2 仿真結(jié)果

本仿真模擬了1 000 Hz中心頻率下的聲信號(hào)定位仿真，計(jì)算了當(dāng)信噪比為5 dB時(shí)的定位結(jié)果，聲源信號(hào)分別采用單頻信號(hào)和高斯白噪聲信號(hào)進(jìn)行模擬，噪聲為高斯白噪聲，將噪聲添加到麥克風(fēng)陣列接收的聲信號(hào)上；定位結(jié)果云圖中環(huán)形虛線表示聲源旋轉(zhuǎn)軌跡，虛線上的*表示的是該切片的中心時(shí)刻聲源真實(shí)所在位置，如圖3所示。

(a) 仿真聲源初始位置

圖4給出了當(dāng)源信號(hào)為1 000 Hz單頻信號(hào)時(shí)的定位結(jié)果，分別展示了未經(jīng)相位平均的常規(guī)波束形成算法以及通過(guò)相位平均處理的3種算法的定位結(jié)果，同時(shí)還給出了點(diǎn)擴(kuò)散函數(shù)的圖像?？梢钥闯觯珻BF將聲源的整個(gè)旋轉(zhuǎn)軌跡定位出來(lái)，無(wú)法準(zhǔn)確定位某一處位置的聲源。經(jīng)過(guò)相位平均處理之后，RBF,RMVDR和RDAMAS均可以實(shí)現(xiàn)準(zhǔn)確定位，其中，RDAMAS算法在RBF的基礎(chǔ)進(jìn)行了解卷積操作，有效提升了分辨率。但是RMVDR算法分辨率不高，這是由于源信號(hào)為單頻信號(hào)，導(dǎo)致麥克風(fēng)陣列各通道信號(hào)相關(guān)性較強(qiáng)，使得互譜矩陣的秩較小，在求逆時(shí)會(huì)出現(xiàn)病態(tài)解。

圖5展示了信號(hào)源改為高斯白噪聲信號(hào)，其余各項(xiàng)參數(shù)和圖4仿真保持一致的定位結(jié)果。對(duì)比圖4可以看出，RBF和RDAMAS并沒(méi)有發(fā)生明顯變化，但是RMVDR的性能得到了明顯提升，說(shuō)明RMVDR比較適合于高斯白噪聲信號(hào)定位。因?yàn)樾盘?hào)源為寬帶高斯信號(hào)時(shí)，信號(hào)之間相關(guān)性小，互譜矩陣幾乎是滿秩矩陣，求逆不會(huì)產(chǎn)生病態(tài)解，而對(duì)于單頻等窄帶信號(hào)由于式(15)中互譜矩陣求逆產(chǎn)生的誤差會(huì)影響定位精度。

(a) CBF

(a) CBF

除此之外，從仿真的結(jié)果可以看出轉(zhuǎn)速波動(dòng)對(duì)本文所述的相位平均法沒(méi)有造成影響。這是因?yàn)橄辔黄骄ǖ幕驹硎菍⒃夹盘?hào)劃分成多個(gè)切片分塊處理，每一個(gè)切片內(nèi)的聲源位移微小，可以近似等效為靜止聲源求解，所以不必關(guān)心每個(gè)時(shí)刻聲源的具體轉(zhuǎn)速與運(yùn)動(dòng)軌跡，只需采集到每個(gè)切片內(nèi)的聲信號(hào)就可進(jìn)行運(yùn)算，這也是相位平均法與時(shí)域去多普勒效應(yīng)法的本質(zhì)區(qū)別。時(shí)域去多普勒效應(yīng)是在已知聲源運(yùn)動(dòng)軌跡和轉(zhuǎn)速的條件下，通過(guò)時(shí)域插值法由麥克風(fēng)信號(hào)重建聲源位置信號(hào)，所以當(dāng)轉(zhuǎn)速存在隨機(jī)波動(dòng)時(shí)，無(wú)法準(zhǔn)確判斷掃描區(qū)域的網(wǎng)格點(diǎn)(等效源)在各個(gè)采樣時(shí)刻的位置，插值時(shí)會(huì)產(chǎn)生較大誤差。這也說(shuō)明相位平均法比時(shí)域去多普勒效應(yīng)法的魯棒性更高，能夠適應(yīng)轉(zhuǎn)速波動(dòng)在±10%以內(nèi)的旋轉(zhuǎn)聲源定位。

3 試驗(yàn)驗(yàn)證

3.1 試驗(yàn)參數(shù)

為了驗(yàn)證本文提出的方法的有效性，我們進(jìn)行了旋轉(zhuǎn)藍(lán)牙音箱試驗(yàn)和工業(yè)風(fēng)扇試驗(yàn)，兩組試驗(yàn)均在消聲室中進(jìn)行。其中，藍(lán)牙音箱和工業(yè)風(fēng)扇的各項(xiàng)試驗(yàn)參數(shù)如表2所示。麥克風(fēng)陣列采用和仿真相同的56通道螺旋陣列，麥克風(fēng)采用1/4英寸駐極體麥克風(fēng)，頻率響應(yīng)20 Hz～20 kHz。旋轉(zhuǎn)藍(lán)牙音箱正對(duì)麥克風(fēng)陣列平面發(fā)聲，測(cè)試距離1.3 m，位置關(guān)系如圖6所示。工業(yè)風(fēng)扇由于旋轉(zhuǎn)時(shí)產(chǎn)生的氣流太大，防止麥克風(fēng)陣列采集信號(hào)被過(guò)多的流滯噪聲干擾，將風(fēng)扇背對(duì)陣列面，距離1.52 m，位置關(guān)系如圖7所示。

表2 試驗(yàn)參數(shù)

圖6 旋轉(zhuǎn)藍(lán)牙音箱位置圖

圖7 工業(yè)風(fēng)扇位置圖

3.2 試驗(yàn)結(jié)果

3.2.1 藍(lán)牙音箱旋轉(zhuǎn)試驗(yàn)結(jié)果

如圖8所示為旋轉(zhuǎn)藍(lán)牙音箱在600 r/min轉(zhuǎn)速，2 500 Hz單頻信號(hào)下的定位結(jié)果。采樣頻率為50 kHz，因此聲源旋轉(zhuǎn)一周總采樣點(diǎn)數(shù)為5 000，將每一圈的麥克風(fēng)陣列采集到的原始信號(hào)劃分成25個(gè)切片，單個(gè)切片的采樣點(diǎn)數(shù)為200，單個(gè)切片聲源旋轉(zhuǎn)過(guò)的角度為14.4°。圖8分別展示了CBF的定位結(jié)果，以及基于相位平均處理的RBF，RMVDR和RDAMAS定位結(jié)果，采用RBF的結(jié)果作為式(20)中RDAMAS迭代算法的輸入。CBF，RBF和RDAMAS動(dòng)態(tài)范圍16 dB, RMVDR動(dòng)態(tài)范圍8 dB。虛線大圈表示聲源旋轉(zhuǎn)圓周，每一個(gè)切片下聲源實(shí)際位置用虛線小圈表示(注：由于測(cè)量和對(duì)準(zhǔn)誤差，標(biāo)示位置可能和實(shí)際位置略有微小誤差)。

(a) CBF

圖8(a)的CBF因?yàn)閷?duì)整段采樣信號(hào)進(jìn)行波束形成，所以無(wú)法準(zhǔn)確定位某一處位置聲源，只能得到聲源運(yùn)動(dòng)軌跡。RBF，RMVDR和RDAMAS均可以準(zhǔn)確定位對(duì)應(yīng)切片下的聲源位置，由于4個(gè)藍(lán)牙音箱一致性不太好，所以展示出來(lái)的聲壓幅值大小不一致，在定位云圖上表現(xiàn)為對(duì)應(yīng)聲源顏色深淺不一樣，圖8(c)和圖8(d)尤為明顯。另外對(duì)比圖8(c)和圖8(d)聲源旋轉(zhuǎn)過(guò)的角度基本等于4個(gè)切片長(zhǎng)度下聲源真實(shí)旋轉(zhuǎn)過(guò)的角度(14.4°×4=57.6°)。對(duì)比圖8(c)和圖8(g)，RDAMAS在RBF基礎(chǔ)上采用迭代法極好的抑制了旁瓣，提高了定位分辨率。圖8(e)和圖8(f)展示的RMVDR結(jié)果雖然也能夠定位出聲源位置，但成像較模糊，分析原因可能是單頻信號(hào)的相關(guān)性較大，導(dǎo)致互譜矩陣秩較小，是一個(gè)病態(tài)矩陣，在求逆時(shí)容易產(chǎn)生較大誤差，這與圖4(e)和圖4(f)仿真結(jié)果一致。

其余試驗(yàn)參數(shù)不變，僅將聲源信號(hào)頻率提高到4 000 Hz，結(jié)果如圖9所示。對(duì)比圖8(b)和圖9(b)可以看出，隨著頻率升高，PSF變得尖銳，分辨率提高，這與理論分析相吻合。同時(shí)，各種算法的空間分辨率也隨著頻率提高而升高，尤其是RBF的旁瓣明顯削弱，性能提升顯著。但是，隨著頻率升高，波長(zhǎng)變短，可能會(huì)超出陣列的空間采樣極限(麥克風(fēng)間距不得大于半個(gè)波長(zhǎng))，因此定位結(jié)果中出現(xiàn)了較多的柵瓣影響。

(a) CBF

將圖8的聲源信號(hào)改為高斯白噪聲信號(hào)，定位中心頻率依然為2 500 Hz，其余各項(xiàng)參數(shù)保持不變，結(jié)果如圖10所示。對(duì)比圖10和圖8，RBF和RDAMAS的結(jié)果變化不大，但是RMVDR性能提升明顯，分辨率顯著提高，說(shuō)明RMVDR算法對(duì)高斯白噪聲信號(hào)適應(yīng)性更好，這與仿真的結(jié)果相一致。

(a) CBF

3.2.2 工業(yè)風(fēng)扇試驗(yàn)結(jié)果

工業(yè)風(fēng)扇為三葉片軸流風(fēng)扇，分別采集了1 130 r/min轉(zhuǎn)速和1 278 r/min轉(zhuǎn)速下的信號(hào)，測(cè)試距離1.52 m。圖12為分析頻率2 500 Hz的結(jié)果，單個(gè)切片采樣點(diǎn)數(shù)為100，1個(gè)大圓虛線框表示葉尖的運(yùn)行軌跡，4個(gè)小圓虛線框表示對(duì)應(yīng)切片下的扇葉實(shí)際位置(注：由于測(cè)量和對(duì)準(zhǔn)誤差，標(biāo)示位置可能和實(shí)際位置略有微小誤差)。CBF依然無(wú)法準(zhǔn)確定位聲源位置，基于相位平均處理的RBF，RMVDR和RDAMAS均可以準(zhǔn)確定位該采樣切片對(duì)應(yīng)的聲源位置。其中，RBF的分辨率較差，這是因?yàn)楸绕鹚{(lán)牙音箱旋轉(zhuǎn)，工業(yè)風(fēng)扇的氣動(dòng)噪聲源聲場(chǎng)更加復(fù)雜，干擾大，而且聲源不是呈現(xiàn)離散點(diǎn)聲源分布，而是呈現(xiàn)分布式聲源的分布特征(即不滿足式(5)中的K?N)。相對(duì)而言，RMVDR的結(jié)果比較理想，這是因?yàn)闅鈩?dòng)噪聲源的頻譜呈現(xiàn)寬帶連續(xù)譜+低頻離散線譜的特征，比較適用于RMVDR算法，具體表現(xiàn)形式如圖11所示。低頻相鄰線譜頻率差值約為55.3 Hz，約等于風(fēng)扇的葉頻(1 130/60×3=56.5 Hz)，中高頻為寬帶連續(xù)譜，分析頻率選為中高頻率。將RMVDR算法的結(jié)果作為RDAMAS迭代的輸入，對(duì)比圖12(e)和圖12(g)，經(jīng)過(guò)迭代運(yùn)算之后，定位分辨率顯著提升。

(a) 時(shí)域波形

RBF，RMVDR和RDAMAS算法定位的聲源位置均分布在扇葉尖部，這是因?yàn)閷?duì)于旋轉(zhuǎn)葉片，氣流流過(guò)葉尖引起的不規(guī)則流動(dòng)產(chǎn)生的渦流噪聲是氣動(dòng)噪聲源的主要組成部分，而這一部分噪聲主要為中高頻噪聲[24]。隨機(jī)渦從尾緣脫落，引起與其相反方向的局部環(huán)量變化，從而引起槳盤(pán)上誘導(dǎo)速度的變換，導(dǎo)致葉尖產(chǎn)生隨機(jī)升力的脈動(dòng)，最終產(chǎn)生頻帶較寬的渦流脫落噪聲。渦流脫落噪聲位于葉尖，圖12的工業(yè)風(fēng)扇試驗(yàn)定位的最大噪聲源就位于葉輪尖部，與理論分析相吻合。

(a) CBF

將風(fēng)扇轉(zhuǎn)速調(diào)為高檔(1 278 r/min)，分析頻率為3 500 Hz，其余試驗(yàn)條件不變，定位結(jié)果如圖13所示。RBF，RMVDR和RDAMAS算法依舊保持很好的性能，RDAMAS依舊展現(xiàn)出極高的分辨率。隨著分析頻率增加，從圖13(b)可以看出，PSF變得更加尖銳，分辨率提高；定位結(jié)果的分辨率也明顯提高，尤其是RBF對(duì)旁瓣的抑制作用變好。

(a) CBF

4 結(jié) 論

本文基于相位平均原理提出了針對(duì)旋轉(zhuǎn)聲源的定位方法，并給出了相位平均技術(shù)所適用的條件。其基本思路是將旋轉(zhuǎn)聲信號(hào)進(jìn)行切片平均，在每一個(gè)切片內(nèi)近似等效于靜止聲源處理，并在此基礎(chǔ)上提出了基于解卷積算法的RDAMAS算法，在波束形成基礎(chǔ)上顯著提高了定位的空間分辨率，通過(guò)仿真、旋轉(zhuǎn)藍(lán)牙音箱試驗(yàn)以及工業(yè)風(fēng)扇試驗(yàn)驗(yàn)證了算法的有效性。由于相位平均只對(duì)原始采樣信號(hào)進(jìn)行切片重排序，不像插值法去多普勒效應(yīng)那樣會(huì)改變信號(hào)自身原始結(jié)構(gòu)，所以該方法具有很高的魯棒性，幾乎所有的靜止高分辨率方法均可以基于相位平均來(lái)實(shí)現(xiàn)對(duì)旋轉(zhuǎn)聲源的定位。除此之外，本文還引入了PSF來(lái)表征陣列性能以及頻率對(duì)于定位分辨率的影響，通過(guò)PSF圖像的尖銳程度反映定位的分辨率大小。

仿真和試驗(yàn)表明，RBF能保持較高的魯棒性，但分辨率不足；而RMVDR更適用于高斯白噪聲信號(hào)定位，空間分辨率高于RBF。RDAMAS在RBF或RMVDR基礎(chǔ)上進(jìn)行迭代運(yùn)算，大幅提高了定位的分辨率。綜合1.2節(jié)針對(duì)式(12)的分析討論，再結(jié)合大量仿真和實(shí)驗(yàn)結(jié)果，可以大致給出本文所述方法的適用范圍：當(dāng)采樣頻率保持50 kHz,測(cè)試距離2 m時(shí)，本方法適用于聲源轉(zhuǎn)速≤3 000 r/min，分析頻率≥1 500 Hz,旋轉(zhuǎn)半徑0.3 m≤R≤1.0 m，聲源轉(zhuǎn)速波動(dòng)在基本轉(zhuǎn)速的±10%的旋轉(zhuǎn)單極子定位成像。針對(duì)三葉片工業(yè)風(fēng)扇，本方法可以實(shí)現(xiàn)對(duì)轉(zhuǎn)速1 278 r/min風(fēng)扇葉片的葉尖氣動(dòng)噪聲源定位，由此本方法可以推廣應(yīng)用到中低轉(zhuǎn)速(≤1 500 r/min)下，葉片數(shù)量較少(不超過(guò)5片)的風(fēng)機(jī)氣動(dòng)噪聲源定位。由于氣動(dòng)噪聲呈現(xiàn)出中高頻(≥1 500 Hz)寬帶譜的特性，所以能夠采用本文所述的算法進(jìn)行定位，而且定位的能量最大的聲源位置位于葉輪尖部。