基于矩形單元波函數(shù)疊加法的平面近場聲全息

范少杰 向宇 石梓玉

摘 要：波疊加法和等效源法都是平面近場聲全息中的常用方法。傳統(tǒng)波疊加法需要對所有單元計算數(shù)值積分，效率較低；等效源法直接將單元積分簡化為單極子點源，計算效率獲得提高，但其過度簡化導(dǎo)致精度損失。針對以上問題，在基于波疊加法的平面近場聲全息基礎(chǔ)上，提出一種替代矩形單元積分的波函數(shù)疊加法。該方法將矩形單元的積分聲場表示成Helmholtz方程球面波譜的波函數(shù)形式，并結(jié)合球面波譜的性質(zhì)，將單元的遠(yuǎn)場球面波譜內(nèi)推至近場，由此獲得了與矩形單元積分聲場等價的外部聲場顯式波函數(shù)解析表達(dá)式，有效避免了傳統(tǒng)波疊加法的數(shù)值積分以及傳統(tǒng)等效源法的過度簡化。利用矩形簡支板的近場聲全息算例對比了所提方法、傳統(tǒng)波疊加法以及傳統(tǒng)等效源法的重建效果。結(jié)果表明：所提方法在全頻段的聲場重建精度均高于傳統(tǒng)等效源法，在高頻段的重建精度高于傳統(tǒng)波疊加法，其計算效率顯著高于傳統(tǒng)波疊加法。

關(guān)鍵詞：近場聲全息；波疊加法；等效源法；波函數(shù)

中圖分類號：TB52 DOI：10.16375/j.cnki.cn45-1395/t.2023.02.001

0 引言

近場聲全息是噪聲檢測與聲源識別定位中的常用技術(shù)，該技術(shù)最早由Williams等[1]在20世紀(jì)80年代應(yīng)用于聲學(xué)計算領(lǐng)域。經(jīng)過各國學(xué)者數(shù)十年研究，由最初的基于空間傅里葉變換算法（SFT）的近場聲全息技術(shù)[2]發(fā)展出了更多適用范圍更廣的近場聲全息算法[3]。為了克服空間傅里葉變換算法只能適用于規(guī)則形狀聲源的局限性，Veronesi等[4]在1989年將邊界元法引入聲場近場聲全息技術(shù)中。Bai[5]對該方法進(jìn)行了改進(jìn)，隨后被大量應(yīng)用于工程實踐中[6-7]。雖然邊界元法適用于任意形狀的聲源，但其復(fù)雜的奇異積分處理也相應(yīng)地提高了計算成本。Koopmann等[8]于1989年將波疊加的思想應(yīng)用到聲場計算中，通過在輻射體內(nèi)部布置虛擬源來模擬聲源外部聲場。由于虛擬源與真實邊界不重合，從而避免了邊界元法中的奇異積分，進(jìn)一步提高了計算效率。

在實際工程及試驗測量中，平面?zhèn)髀暺麝嚵胁贾孟鄬唵?，被廣泛應(yīng)用于聲源近場的聲場重建中[9-10]。在基于傳統(tǒng)波疊加法的平面近場聲全息中，需要計算虛擬邊界上離散單元的數(shù)值積分，其精度較高，但單元積分運算將耗費大量時間成本。若將單元積分簡化為單極子點源，則得到傳統(tǒng)的單極子等效源法。該方法無需計算單元積分，計算效率得到很大提高[11]，但從單元積分簡化為點源的過程中，由于存在較大的邊界離散誤差和積分近似誤差[12-13]，其計算精度低于波疊加法。

為提高傳統(tǒng)波疊加法的計算效率并兼顧其計算精度，本文以工程中常用的矩形單元為例，提出了一種利用具有顯式解析表達(dá)式的波函數(shù)代替單元積分的方法。文中詳細(xì)闡述了波函數(shù)的推導(dǎo)過程，并簡要分析了波函數(shù)的性質(zhì)，最后利用簡支板的算例對比了所提方法與傳統(tǒng)方法在聲場重建中的效果。

1 基于波疊加法的平面近場聲全息原理

在基于波疊加法的平面近場聲全息中，通常在求解域之外的虛擬邊界上布置虛擬源來模擬真實聲源產(chǎn)生的聲場，如圖1所示。

對于常用矩形采樣陣列，為了保證波疊加法的計算精度以及穩(wěn)定性，一般設(shè)置虛擬面與全息面共形，如圖1（a）所示，其聲壓可由式（1）表示：

[pr=SEq（rE）G（r，rE）dSE（rE）]. （1）

其中：[r]為場點位置矢量；[SE]為與全息面共形的矩形虛擬邊界；[q（rE）]為虛擬源的源強；[rE]為虛擬邊界上虛擬源位置矢量；[G（r，rE）]為自由場格林函數(shù)。若將虛擬邊界離散成[N]個大小形狀均相同的矩形常數(shù)單元，如圖1（b）所示，則式（1）可寫為：

[pr=i=1Nqi0Lx0LyG（r，rE）dSEi（rE） ，i=1，2，…，N]. （2）

式（2）即為傳統(tǒng)波疊加法的聲壓計算公式，其中，[qi]為第[i]個單元的虛擬源強；[Lx]和[Ly]分別為離散單元的邊長；[SEi]為第[i]個離散單元。

由式（2）可以看出，計算任意場點的聲壓時，需要對每個單元進(jìn)行積分計算，因此，使用傳統(tǒng)波疊加法來計算聲場的效率很低。為了解決傳統(tǒng)波疊加法計算效率的問題，通常采用虛擬點源代替單元積分的方法來計算聲場[14-15]。如圖1（c）所示，若在虛擬邊界上布置[N]個單極子虛擬點源，則任意場點的聲壓可由該[N]個單極子點源所產(chǎn)生的聲場疊加表示[16]：

[pr=i=1NqiG（r，rEi）][， ][i=1，2，…，N]. （3）

式（3）即為單極子等效源法，其中，[rEi]為第[i]個等效源的位置。利用聲源近場區(qū)域的全息面測得聲壓數(shù)據(jù)，并代入式（2）或式（3），求出虛擬源強[17]：

[q=G-1HEpH]. （4）

其中：[q]為虛擬源強列向量；[G-1HE]為全息測點與虛擬源間的傳遞矩陣逆矩陣；[pH]為全息測點聲壓數(shù)據(jù)列向量。求出虛擬源的源強之后，再代入式（2）或式（3）即可獲得重建面上場點聲壓，這便是基于波疊加法的近場聲全息原理。

通過上述對基于波疊加法近場聲全息原理的分析可以看出，傳統(tǒng)波疊加法從波疊加方程導(dǎo)出，理論上計算精度更高，但是存在計算效率低的問題。單極子等效源法是對傳統(tǒng)波疊加法的簡化算法，該方法不需要單元插值處理和積分運算，計算效率高，在位置和數(shù)量合適的情況下可以達(dá)到較高的精度[18]。然而從傳統(tǒng)波疊加法到單極子等效源法的近似過程中存在著一定的積分誤差，這些誤差一定程度上影響著單極子等效源法的計算精度。為兼顧波疊加法的計算精度和單極子等效源法的計算效率，本文提出一種波函數(shù)方法。

2 波函數(shù)的構(gòu)造

對于圖1的離散模型，每一個離散的矩形單元大小形狀均相同，因而其聲場必然都是相同的。對此，任取一個單元，并以單元的形心[O]為原點，建立坐標(biāo)系[Oxyz]，其[x]軸和[y]軸分別平行于單元的長和寬，如圖2所示。

由式（2）可知，相對于該單元中心的場點[r]在單元下的積分形式為：

[p（r）=-Lx2Lx2-Ly2Ly2G（r，r）dxdy. ] （5）

式中，[r]是單元上的點的坐標(biāo)矢量。若能將式中的積分形式表示成波函數(shù)[W（r，θ，φ）]的形式：

[p（r）=W（r，θ，φ）]. （6）

則可以在計算任意場點的聲場時，只需把其相對坐標(biāo)[r、θ、φ]代入波函數(shù)[W（r，θ，φ）]即可，從而避免積分。由球面波輻射理論[2]可知，在聲源的外域，聲壓滿足球面坐標(biāo)系下Helmholtz方程的解，從而式（6）可用球面波展開：

[pr=n=0∞ m=-nnPmnrYmnθ，φ]. （7）

式中：[Pmn（r）]是半徑為[r]的球面上的第[（n，m）]階聲壓球面波譜；[Ymnθ，φ]為歸一化球諧函數(shù)。

[Ymnθ，φ=2n+14πn-m！n+m！Pmncosθeimφ]. （8）

其中：[Pmncosθ]為連帶勒讓德函數(shù)，i=[-1]。由式（7）可知，只要求出球面波譜[Pmn（r）]即可獲得[p（r）]的解。

由球面波譜理論可知，球面波譜[Pmn（r）]可通過已知半徑為[R]的球面波譜[Pmn（R）]來求解：

[Pmn（r）= h（2）n（kr）h（2）n（kR）Pmn（R）]. （9）

式中，[h（2）n]（*）表示第二類球[Hankel]函數(shù)，其中，球面波譜[Pmn（R）]可利用球諧函數(shù)的完備正交性進(jìn)行求解[19]：

[Pmn（R）= 02π0πpRYmnθ，φsinθdθdφ]. （10）

式中：[pR]為該單元在半徑[R]的球面邊界上的聲壓分布；[Ymnθ，φ]是[Ymnθ，φ]的復(fù)共軛。因此，只要確定了[pR]，即可由式（10）求得球面波譜[Pmn（R）]。對于矩形常數(shù)單元，一般有2種確定[pR]的方法，一種方法是采用式（5）計算：

[p（R）=-Lx2Lx2-Ly2Ly2G（R，r）dxdy.] （11）

該方法適用于計算任意半徑球面邊界的聲壓，但需要計算數(shù)值積分。將式（11）獲得的[pR]代入式（10）可得到任意半徑[R]球面的球面波譜[Pmn（R）]：

[Pmn（R）=][-Lx2Lx2-Ly2Ly202π0πG（R，r）Ymn（θ，φ）sinθdθdφdxdy.]

（12）

另一種方法是直接采用格林函數(shù)的遠(yuǎn)場近似表達(dá)式獲得 [20]：

[pRF= LxLyexp-ikRF4πRFsinckxLx2sinckyLy2].

（13）

式中：[kx= ksinθcosφ]，[ky= ksinθsinφ]，i=[-1]。該方法適用于計算半徑[RF]遠(yuǎn)大于單元尺寸的遠(yuǎn)場球面邊界的聲壓。將獲得的[pRF]代入式（10）可得到遠(yuǎn)場球面波譜[Pmn（RF）]：

[Pmn（RF）=02π0π LxLyexp（-ikRF）4πRFsinckxLx2]·

[sinckyLy2Ymnθ，φsinθdθdφ]. （14）

式中，[sinc?=sin?/?]?？梢钥闯觯嬎闳我馇蛎孢吔绲那蛎娌ㄗV[Pmn（R）]需要四重積分，而計算遠(yuǎn)場球面邊界的球面波譜[Pmn（RF）]僅需要二重積分，求解效率更高。因此，本文均采用遠(yuǎn)場球面波譜[Pmn（RF）]來構(gòu)造波函數(shù)。將式（7）、式（9）、式（14）聯(lián)立得單元的波函數(shù)：

[Wr，θ，φ=n=0∞m=-nnh（2）nkrh（2）nkRFPmn（RF）Ymnθ，φ]. （15）

于是，任意場點的聲壓可由[N]個單元的波函數(shù)疊加求得：

[p（r）=i=1NqiWri，θi，φi ， i=1，2，…，N]. （16）

其中：[qi]為第[i]個波函數(shù)的源強；[ri、θi、φi]分別為場點相對于第[i]單元的相對坐標(biāo)值。

3 數(shù)值仿真

在采用波函數(shù)進(jìn)行計算時，必須先確定截斷項數(shù)。理論上來說，選取的階數(shù)越多，波函數(shù)越逼近精確的積分聲場，但項數(shù)過多會導(dǎo)致計算效率降低。經(jīng)大量試算發(fā)現(xiàn)，僅需取[n=2]，就可以保證波函數(shù)具有較好的計算效果。而且由于矩形聲源的對稱特性，只有當(dāng)[n]和[m]均為偶數(shù)時，球面波譜為非零解，即僅需求解[4]個球面波譜：[P0，0RF]、[P0，2RF]、[P-2，2RF]、[P2，2RF]。下文的算例均采用該參數(shù)進(jìn)行計算。

在實際工程中，連續(xù)分布的結(jié)構(gòu)振動聲源較為常見，如板、殼等結(jié)構(gòu)振動產(chǎn)生的聲場。因此，本算例采用受間歇激勵作用的四周無限大矩形簡支板聲源來比較波函數(shù)法、單極子等效源法和傳統(tǒng)波疊加法的聲場重建效果，仿真模型如圖3所示。全息面、重點面、虛擬面及聲源的相對位置如圖3（a）所示，矩形簡支板為鋁板，其長寬為[0.5 m×0.5 m]，厚度為[0.002 m]，楊氏模量[E=7.1×1010 Pa]，泊松比[υ=0.33]，密度為[2.7×103 kg/m3]。簡支板位于[Z=0]處的[XY]平面上，對簡支板施加幅值為[1 N]的中心激勵，激勵點位置為[（0.25，0.25，0）]；全息面大小與簡支板相同，位于簡支板上方[0.05 m]，長寬方向均按[0.025 m]間隔劃分共計[21×21]個測點；重建面位于簡支板上方[0.03 m]處，由于該矩形簡支板為中心激勵的正方形簡支板，其產(chǎn)生的聲場有中心對稱的特點，因此在重建面上分別選取A（0.25，0.25，0.03）、B（0.35，0.35，0.03）、C（0.45，0.45，0.03）、D（0.25，0.45，0.03）、E（0.35，0.25，0.03）等5個點作為重建點，如圖3（b）所示。虛擬面單元劃分方式與全息面相同，共計[20×20]個大小相同的正方形虛擬單元，并將其設(shè)置在重建面下方[1]倍虛擬單元距離處，等效源點取為各單元的中心點。全息面和重建點的解析復(fù)聲壓和振速值由瑞利積分獲得 [21]。分別用3種方法重建5個點處的聲壓與振速，重建頻率范圍設(shè)置為[f∈［100，4 500］] Hz，重建頻率間隔設(shè)置為[10 Hz]，重建誤差（[δ]）由式（17）計算：

[δ= E-ErE×100% ]. （17）

其中：[E]為解析值，[Er]為重建值。

圖4顯示了3種方法重建[5]個場點處的聲壓和振速幅值的相對誤差。可以看出，對于不同位置的重建點，波函數(shù)法和傳統(tǒng)波疊加法的總體重建效果（尤其是振速）均優(yōu)于單極子等效源法。其中，當(dāng)重建點位于重建面中心時，波函數(shù)法和傳統(tǒng)波疊加法的重建誤差略低于單極子等效源法；當(dāng)重建點逐漸遠(yuǎn)離中心時，波函數(shù)法和傳統(tǒng)波疊加法的重建誤差明顯低于單極子等效源法；在低頻段，波函數(shù)法的聲壓重建精度與傳統(tǒng)波疊加法相當(dāng)，振速重建精度略低于傳統(tǒng)波疊加法；但在中高頻段，波函數(shù)法的聲壓和振速重建效果甚至優(yōu)于傳統(tǒng)波疊加法。這是由于本文在采用波疊加法時，僅采用了[4×4]點的高斯勒讓德積分，對于低頻，這些積分點是足夠的，所以其精度很高，但隨著頻率的升高，積分?jǐn)?shù)量的不足導(dǎo)致波疊加法的精度有所下降。但本文方法由顯示解析表達(dá)式構(gòu)成，并不受限于單元積分點的數(shù)量，因此在高頻具有一定的優(yōu)勢。

表1中記錄了該算例中3種方法所消耗的時間?？梢园l(fā)現(xiàn)，形式最為簡單的等效源法在所有方法中的耗時最短，其計算時間僅為波疊加法的1/16。波函數(shù)法的計算時間雖長于等效源法，但由于無需數(shù)值積分，因而其計算時間僅為波疊加法的1/6。綜合考慮計算精度與計算效率，波函數(shù)方法相較于等效源法和波疊加法均具有一定優(yōu)勢。

4 結(jié)論

針對傳統(tǒng)波疊加法數(shù)值積分計算效率低和單極子等效源法過度簡化的問題，提出了一種波函數(shù)代替積分的方法，該方法兼顧了傳統(tǒng)波疊加法的計算精度和單極子等效源法的計算效率。文中對所提方法進(jìn)行了詳細(xì)推導(dǎo)和闡述，并通過數(shù)值仿真進(jìn)行了驗證。計算結(jié)果表明：本文方法在全頻段的計算精度，尤其是振速的計算精度均高于單極子等效源法，在高頻段的計算精度高于波疊加法，計算效率遠(yuǎn)高于波疊加法。目前本文方法僅應(yīng)用于平面近場聲全息中，對于其在任意形狀近場聲全息中的應(yīng)用，有待進(jìn)一步深入研究。

參考文獻(xiàn)

[1] WILLIAMS E G，MAYNARD J D. Holographic imaging without the wavelength resolution limit[J]. Physical Review Letters，1980，45（7）：554-557.

[2] WILLIAMS E G. Fourier acoustics：sound radiation and nearfield acoustical holography[M]. Pittsburgh：Academic Press，1999.

[3] 陳心昭，畢傳興. 近場聲全息技術(shù)及其應(yīng)用[M]. 北京：科學(xué)出版社，2013.

[4] VERONESI W A，MAYNARD J D. Digital holographic reconstruction of sources with arbitrarily shaped surfaces[J]. Journal of the Acoustical Society of America，1989，85（2）：588-598.

[5] BAI M R. Application of BEM （boundary element method） based acoustic holography to radiation analysis of sound sources with arbitrarily shaped geometries[J]. Journal of the Acoustical Society of America，1992，92（1）：533.

[6] 吳飛，余神志，劉蘇行，等. 車用制動空氣壓縮機結(jié)構(gòu)輻射噪聲特性研究[J]. 聲學(xué)技術(shù)，2021，40（2）：228-233.

[7] 張建民，安俊英. 深海波導(dǎo)中目標(biāo)低頻聲散射特性研究[J]. 聲學(xué)技術(shù)，2022，41（3）：419-425.

[8] KOOPMANN G H，SONG L，F(xiàn)AHNLINE J B. A method for computing acoustic fields based on the principle of wave superposition[J]. Journal of the Acoustical Society of America，1989，86（6）：2433-2438.

[9] 張永斌，畢傳興，陳劍，等. 基于等效源法的平面近場聲全息及其實驗研究[J]. 聲學(xué)學(xué)報，2007，32（6）：489-496.

[10] 陳漢濤，郭文勇，韓江桂，等. 船舶機艙內(nèi)高頻弱聲源近場聲全息方法[J]. 艦船科學(xué)技術(shù)，2019，41（21）：138-143，147.

[11] CHELLIAH K，RAMAN G，MUEHLEISEN R T. An experimental comparison of various methods of nearfield acoustic holography[J]. Journal of Sound and Vibration，2017，403：21-37.

[12] VALDIVIA N P，WILLIAMS E G. Study of the comparison of the methods of equivalent sources and boundary element methods for near-field acoustic holography[J]. Journal of the Acoustical Society of America，2006，120（6）：3694-3705.

[13] VALDIVIA N P. Advanced equivalent source methodologies for near-field acoustic holography[J]. Journal of Sound and Vibration，2019，438：66-82.

[14] 石梓玉，向宇，陸靜，等. 一種提高聲場重構(gòu)穩(wěn)定性的射線等效源法[J]. 廣西科技大學(xué)學(xué)報，2019，30（3）：1-7，21.

[15] 陳巖豪，向宇，石梓玉. 基于波疊加近場聲全息射線波函數(shù)的信息補償[J]. 廣西科技大學(xué)學(xué)報，2021，32（3）：6-12.

[16] WU S W，XIANG Y. Location optimization of monopole equivalent sources in wave superposition method[J]. International Journal of Acoustics and Vibrations，2018，23（2）：254-263.

[17] 張陽，向宇，石梓玉. 二維聲全息的δ函數(shù)約束型射線波疊加法[J]. 廣西科技大學(xué)學(xué)報，2020，31（4）：19-26，48.

[18] BAI M R，CHEN C C，LIN J H. On optimal retreat distance for the equivalent source method-based nearfield acoustical holography[J]. Journal of the Acoustical Society of America，2011，129（3）：1407-1416.

[19] 丁丹丹. 基于球形陣的封閉空間噪聲源定位[D]. 哈爾濱：哈爾濱工程大學(xué)，2017.

[20] 王斌，湯渭霖，范軍. 一種輻射聲場近似計算方法：單元輻射疊加法[J]. 聲學(xué)學(xué)報，2008，33（3）：226-230.

[21] WILLIAMS E G，MAYNARD J D. Numerical evaluation of the Rayleigh integral for planar radiators using the FFT[J]. Journal of the Acoustical Society of America，1982，72（6）：2020-2030.

Planar near-field acoustic holography based on rectangular element wave function superposition method

FAN Shaojie1，2， XIANG Yu*2， SHI Ziyu2

（1. School of Mechanical and Automotive Engineering， Guangxi University of Science and Technology， Liuzhou 545616， China; 2. Guangxi Key Laboratory of Automobile Component and Vehicle Technology

（Guangxi University of Science and Technology）， Liuzhou 545616， China）

Abstract： Both wave superposition and equivalent source methods are commonly used in planar near-field acoustic holography. The conventional wave superposition method needs to calculate numerical integration for all elements， which is less efficient; the equivalent source method directly simplifies the element integration to a monopole point source， and the computational efficiency is improved， but its over-simplification leads to the loss of accuracy. To address the above problems， a wave function superposition method based on the wave superposition method for planar near-field acoustic holography is proposed as an alternative to rectangular element integration. The method represents the integrated sound field of the rectangular element as the wave function form of the spherical wave spectrum of Helmholtz equation， and combines the nature of the spherical wave spectrum to interpolate the far-field spherical wave spectrum of the cell to the near field， thus obtaining the analytical expression of the explicit wave function of the external sound field equivalent to the integrated sound field of the rectangular element， effectively avoiding the numerical integration of the conventional wave superposition method and the oversimplification of the traditional equivalent source method. The reconstruction effects of the proposed method， the conventional wave superposition method and the conventional equivalent source method are compared using the near-field acoustic holographic algorithm of a rectangular simply supported plate. The results show that the reconstruction accuracy of the proposed method is higher than that of the conventional equivalent source method in all frequency bands， and the reconstruction accuracy in the high frequency bands is higher than that of the conventional wave superposition method， and its computational efficiency is significantly higher than that of the conventional wave superposition method.

Key words： near-field acoustic holography; wave superposition method; equivalent source method; wave function

（責(zé)任編輯：黎 婭）