石梓玉 向宇 陸靜 王玉江

摘要 在基于傳統(tǒng)波疊加法的近場聲全息技術中，多采用輻射球面波的單極子作為等效源，易導致傳遞矩陣病態(tài)，利用射線波函數(shù)替換球面波函數(shù)可有效改善傳遞矩陣病態(tài)性。然而，以往的射線波函數(shù)法采用格林函數(shù)的方向導數(shù)作為波函數(shù)，其解析表達式復雜，計算效率低。此外，以往方法的波函數(shù)指向設置對節(jié)點分布方式要求較高，限制了其應用的靈活性。針對上述問題，采用（n， 0）階的球面波源重新構造了一系列射線波函數(shù)，該射線波函數(shù)可利用球Hankel函數(shù)和Legendre多項式的遞推形式方便地計算出其任意階的表達式，大幅提高了效率。通過改進射線波函數(shù)的主指向設置，使其在實際使用中更加靈活，提出一種基于正交球面波源的射線波函數(shù)波疊加法。利用正四面體輻射體、兩端帶球帽的圓柱輻射體和簡支矩形鋼板聲源3個數(shù)值仿真，對比驗證了傳統(tǒng)方法和所提方法在聲場重建中的效果。仿真結果表明：在聲場細節(jié)信息較為豐富的高頻下，即便采用正則化方法求解，傳統(tǒng)波疊加法由于傳遞矩陣病態(tài)嚴重，在3個仿真中均難以保證重建精度，誤差為20%左右；而建立的射線波函數(shù)法則有效降低了傳遞矩陣的條件數(shù)，改善了系統(tǒng)病態(tài)性，獲得了更高的重建精度，誤差在5%～10%之間，說明了所提方法更具優(yōu)越性。

關鍵詞 近場聲全息; 波疊加法; 正交球面波源; 射線波函數(shù)

引 言

近場聲全息技術（NAH）是一種有效的噪聲源定位、識別和聲場重建技術，經(jīng)過近幾十年的研究，在算法方面已相繼提出了基于空間Fourier變換的NAH方法［1］、基于Kirchhoff?Helmholtz邊界積分方程的邊界元方法［2?4］、源強模擬方法［5?6］、波疊加法［7］等。其中，波疊加法作為一種無奇異性、精度高且適用于任意形狀結構的聲場計算方法，自1989年提出以來就已被廣泛應用于各種聲學問題的計算中［8?11］。其原理是在聲源面內(nèi)縮的一個封閉虛擬曲面上布置連續(xù)分布的等效源來表示聲源向外輻射的聲場，不僅克服了基于空間Fourier變換算法只能計算規(guī)則形狀聲源的缺點，又避免了邊界元法所帶來的復雜插值運算和奇異積分處理。但將波疊加法應用于聲全息計算時，由于等效源面與全息測量面之間距離的影響，其離散后形成的傳遞矩陣通常是一個大條件數(shù)的病態(tài)矩陣，導致源強求解穩(wěn)定性較差［12］。

為了提高源強求解的穩(wěn)定性，通常需采用正則化方法，目前最常用的正則化方法有截斷奇異值方法（TSVD）、Tikhonov正則化方法等［13］。但這些正則化方法本質上都是將傳遞矩陣中對測量誤差非常敏感的小奇異值項進行截斷或濾除，該過程必然會損失一部分聲場細節(jié)信息。如果傳遞矩陣病態(tài)嚴重，那么正則化時就必須選取較大的截斷點或正則化參數(shù)以過濾更多的小奇異值項，這會加劇聲場細節(jié)信息的丟失。對于復雜程度不高的低頻聲場，小奇異值項對聲場的貢獻相對較小，在正則化后一般均可保證聲場重建的分辨率。但對于聲場細節(jié)信息豐富的高頻聲場，小奇異值項的濾除會對聲場重建精度造成很大影響。因此，為保證波疊加法的重建精度，即便采用正則化方法，也應盡量改善傳遞矩陣的病態(tài)性。

為改善傳遞矩陣的病態(tài)性，以往的研究大多側重于優(yōu)化等效源布置面的位置［12，14?18］。但由于聲源形狀和性質、全息測量面形狀、測點分布方式等的復雜性，等效源面最佳分布和位置的選擇是一個非常復雜的問題，且至今尚無一個成熟有效的方法［12，19］。文獻［20?23］在對傳統(tǒng)單極子波疊加法傳遞矩陣的病態(tài)性進行分析后發(fā)現(xiàn)，由于單極子向外輻射的波函數(shù)是以球面形式衰減的自由場格林函數(shù)，因此當?shù)刃г袋c之間或全息測量點之間的位置僅有微小改變時，格林函數(shù)的大小變化過于平緩，導致傳遞矩陣因不同行或列近似相等而病態(tài)。進而，文獻［22?23］從優(yōu)化波函數(shù)波陣面衰減速度的角度，提出了一種利用射線波函數(shù)替代傳統(tǒng)球面波函數(shù)以改善傳遞矩陣病態(tài)性的射線波函數(shù)法，在改善傳遞矩陣病態(tài)性方面取得了一定效果。然而，該方法采用格林函數(shù)的方向導數(shù)作為射線波函數(shù)，這種類型射線波函數(shù)的解析表達式較為復雜，不僅計算效率低且難以計算到高階導數(shù)。此外，該方法還將射線波函數(shù)的主指向設置為各等效源到其對應測點的方向，即要求等效源和測點無論是在數(shù)量還是分布方式上均必須一一對應，這無疑限制了射線波函數(shù)法在實際工程中的應用。

為了提高射線波函數(shù)法的計算效率和靈活性，本文采用（n， 0）階正交球面波源構造了一種新型的射線波函數(shù)，該波函數(shù)無需求導運算，僅利用球Hankel函數(shù)和Legrend多項式的遞推形式即可方便地計算出其任意階的表達式，不僅計算效率得到大幅提高，而且不會增加計算難度和復雜性。與此同時，本文通過分析射線波函數(shù)對傳遞矩陣病態(tài)性的改善機理，提出一種射線波函數(shù)的主指向設置方法，使等效源與測點的布置不再受到數(shù)量和分布均須一一對應的限制，大大提高了射線波函數(shù)法在應用中的靈活性。最后，通過3個不同聲源的數(shù)值仿真算例對比驗證了本文方法在聲場重建中的有效性和準確性。

1 理論部分

1.1 傳統(tǒng)單極子波疊加法的病態(tài)性及射線波函數(shù)法

波疊加法的基本思想是：聲源向外輻射的聲場可由連續(xù)分布于其內(nèi)部的等效源所輻射的聲場疊加代替，該思想可用下式描述［7］：

式中 p（r）為聲源在空間r處輻射的聲壓；σΦ（rE）為位于rE處的等效源源強，下標Φ表示其波函數(shù)為Φ（r，rE），根據(jù)等效源的類型不同，可為單極子等效源、偶極子等效源、單?偶極子組合型等效源等；Ω為等效源在聲源內(nèi)部的分布區(qū)域，SE為該區(qū)域的邊界，如圖1（a）所示。

為了便于計算，通常將等效源布置在虛擬邊界SE上，并在積分離散時將每個單元內(nèi)的等效源源強和波函數(shù)均視為常數(shù)，且配置在每個單元的中點，如圖1（b）所示。式（1）離散為如下形式［5?6］：

式中 qΦ（rEj）為第j個離散等效源的源強；Φ（r，rEj）為該等效源對應的波函數(shù)。

由式（2）可知，只要確定了N個等效源的源強qΦ（rE1），qΦ（rE2），…，qΦ（rEN），即可計算出聲源在任意場點r處輻射的聲壓。在基于波疊加法的近場聲全息技術中，求解源強所需的方程組可以通過測量聲源近場區(qū)域的聲壓或振速來建立。假設全息測量信息為聲壓，等效源為單極子，利用式（2）可得如下矩陣方程：

式中 pH=[p（rH1）p（rH2）…p（rHM）]T為M×1的測量聲壓列向量；[GH]ij=G（rHi，rEj）為M×N的全息測量數(shù)據(jù)與等效源強間的傳遞矩陣，G（rHi，rEj）為自由場格林函數(shù)，即單極子波函數(shù)；QG=[qG（rE1）qG（rE2）…qG（rEN）]T為N×1的單極子等效源強列向量。

求解方程（3），得源強向量QG為：

式中 G+H表示矩陣GH的廣義逆。

為獲得式（4）的最小二乘解，要求測點數(shù)M不小于等效源數(shù)目N，即M≥N。此外，在傳統(tǒng)波疊加法中，一般采用單極子作為等效源，其波函數(shù)為自由場格林函數(shù)G（rHi，rEj）。文獻［20?23］在對由該函數(shù)所構成的傳遞矩陣GH的病態(tài)性進行分析后指出，由于格林函數(shù)G（rHi，rEj）是一個只與兩點距離有關且以球面衰減的波函數(shù)，因此當不同測點或等效源點之間的位置僅有微小改變時，矩陣的不同行或列的元素將近似相等，導致傳遞矩陣GH因向量間的強線性相關性而病態(tài)。對于該病態(tài)問題，一般是借助正則化方法抑制測量誤差的放大并以此穩(wěn)定求解過程。常用的正則化方法包括截斷奇異值方法（TSVD）、Tikhonov正則化方法［13］等。但正則化方法本質上是通過類似濾波的方法將傳遞矩陣中對誤差敏感的小奇異值項濾除。如果傳遞矩陣病態(tài)嚴重，在正則化時就需選取較大的截斷點或正則化參數(shù)等以濾除更多的奇異值項，這將會加劇聲場細節(jié)信息丟失，導致重建復雜聲源或高頻聲場時精度下降。因而，為保證重建精度，應盡量改善波疊加法傳遞矩陣的病態(tài)性。

為了改善上述問題，文獻［22?23］中提出了一種利用強指向性波函數(shù)替換球面形式波函數(shù)以提高聲場重建穩(wěn)定性的方法，并稱之為射線波函數(shù)法。其基本原理是將傳統(tǒng)波疊加法中單極子等效源輻射的球面波函數(shù)替換為滿足Helmholtz方程和Sommerfield輻射條件且主值指向等效源對應測點的射線波函數(shù)，如圖2所示。這樣一來，等效源輻射的聲波將僅在其對應測點處具有較大的聲波激勵，并生成較大的主對角元素，而在非對應測點處的聲波激勵則快速衰減，即生成較小的非對角元素，進而得到一個主對角元素占優(yōu)的良態(tài)傳遞矩陣，以此提高聲場重建的穩(wěn)定性。

1.2 以往射線波函數(shù)法的缺陷及改進

在對文獻［22?23］中提出的射線波函數(shù)法進行深入研究后發(fā)現(xiàn)，該方法雖然在改善傳遞矩陣病態(tài)性方面具有一定優(yōu)勢，但仍有以下兩個缺陷：

（1）波函數(shù)計算效率方面的缺陷。文獻［22?23］均是利用格林函數(shù)的方向導數(shù)作為射線波函數(shù)，該類型波函數(shù)的高階導數(shù)（超過6階后）解析表達式非常復雜，導致計算效率較低甚至無法計算。

（2）射線波函數(shù)主指向設置方面的缺陷。文獻［22?23］中均是將射線波函數(shù)的主指向設為等效源的對應測點方向，這樣設置雖然可以形成主對角占優(yōu)形態(tài)良好的傳遞矩陣，但要求等效源數(shù)量與測點數(shù)量相同，且在運算過程中須保證這兩組節(jié)點的編號始終一一對應。而在實際應用中，為獲得更高的計算精度，等效源最好均勻布置在與聲源面共形的虛擬面上，同時為了便于制造和降低成本，全息測量面則通常為規(guī)則形狀且測點規(guī)則分布，這必然難以保證全息測點與等效源點之間的一一對應。

本文的主要內(nèi)容則是針對以上兩點提出如下改進辦法：

（1）對波函數(shù)計算效率方面的改進——重新構造射線波函數(shù)。源模擬技術的研究表明，只要是滿足Helmholtz方程和Sommerfeld輻射條件的解析函數(shù)Φ（r，rE）均可作為等效源波函數(shù)［5?6］。因而射線波函數(shù)的選取并不局限于格林函數(shù)的導數(shù)，本文將利用Helmholtz方程在球坐標系下的基本解，即正交球面波源重新構造射線波函數(shù)。

（2）對射線波函數(shù)主指向設置方面的改進。實際上，由1.1節(jié)中圖2的分析可知，射線波函數(shù)能夠降低傳遞矩陣線性相關性的主要原因在于它在非主指向的衰減速度遠快于格林函數(shù)的球面波。因此，即便不采取等效源與主測點一一對應的設置方式，理論上仍應能顯著降低傳遞矩陣因線性相關性過強導致的病態(tài)。在應用中，為了保證等效源輻射的聲波在空間中分布均勻，只要各射線波函數(shù)的主指向均勻向外分散即可。例如，對于封閉的等效源面和全息面，可將射線波函數(shù)的主指向設置為從等效源面所包含空間的幾何中心到各等效源連線的方向，如圖3（a）中的l1，l2，l3，…，lN所示；對于非封閉的等效源面和全息面，則可以設置為等效源面的外法向方向，如圖3（b）中的l1，l2，l3，…，lN所示。在后文的方法驗證部分將通過仿真算例進行驗證。

1.3 基于正交球面波源的射線波函數(shù)

正交球面波源是由不同階次的球諧函數(shù)和球Hankel函數(shù)構成，其具體表達式為［24］：

式中 hn（?）為第一類或第二類n階球Hankel函數(shù)，本文采用第二類n階球Hankel函數(shù)，即h（2）n（?）；Ymn（?）為歸一化的球諧函數(shù)；k為波數(shù)；（r，θ，?）表示場點在坐標系中的位置，如圖4所示。

由式（5）可知，球面波源pnm（r，θ，?）的指向性取決于球諧函數(shù)項Ymn（θ，?）。又由Ymn（θ，?）的性質，當m=0時，其表達式中的角度變量?將被消除，指向性僅由θ決定，此時該函數(shù)的指向形態(tài)必然關于z軸回轉對稱且主值指向z軸。圖5給出了m=0，n分別為0，1，3，7，10，15時球諧函數(shù)的指向形態(tài)圖。由圖5可見，當m=0時，球諧函數(shù)在z軸方向具有強指向性，且n越大，其指向性越強。因而，可取m=0時的球面波源pnm（r，θ，?）作為射線波函數(shù)。令式（5）中m=0，略去常系數(shù)后將其記為Dn（r，θ）：

上式即為帶有參數(shù)n的球面波源型射線波函數(shù)，其中Pn（cosθ）為n階Legendre多項式。由于式（6）可利用球Hankel函數(shù)和Legendre多項式的遞推形式方便地計算到任意階，因此相較于格林函數(shù)導數(shù)型的射線波函數(shù)，不僅計算效率得到大幅提高，而且不會增加計算難度和復雜性。

1.4 射線波函數(shù)在波疊加法中的應用

由1.3節(jié)的分析可知，將射線波函數(shù)應用于波疊加法時，需設置各等效源所輻射的射線波函數(shù)指向其對應主指向l1，l2，l3，…，lN，如圖3中所示。下面假設第j個等效源在全局坐標系Oxyz中的位置為rEj，其對應的主指向為lj，如圖6所示。由于射線波函數(shù)關于自身坐標的z軸回轉對稱，因此，欲使其主瓣指向lj，可將位置rEj作為坐標原點、主指向lj為z軸作一局部坐標系Ox'y'z'。此時，只要將射線波函數(shù)Dn（r，θ）中的變量θ和r替換為局部坐標系Ox'y'z'中的變量θ'和r'，即可使其指向z'軸，即lj方向。

利用式（6），可得rEj處單位強度等效源在場點r處輻射的聲壓為：

為方便起見，將上式中的局部坐標變量r'和θ'用全局坐標變量rEj，r和lj表示如下：

將上式中的場點位置r替換為全息面測點位置rHi，即可得到等效源與全息面間的傳遞矩陣DH，其元素為：

進而式（3）和（4）的矩陣方程可分別改寫為：

式中 QD=[qD（rE1）qD（rE2）…qD（rEN）]T為射線波函數(shù)所對應的等效源強向量。

求解出等效源強后，再利用下式即可重建場點r處的聲壓：

值得一提的是，在n=0和n=1時，式（6）分別對應單極子等效源和偶極子等效源［25］。此時，式（12）將變?yōu)閭鹘y(tǒng)的單層勢和雙層勢波疊加法。

2 方法驗證

2.1 仿真模型說明

為驗證本文所提方法在聲場重建中的效果，該部分設計了3種不同類型聲源的仿真模型，分別為正四面體輻射體、兩端帶球帽的圓柱輻射體和四邊簡支矩形鋼板聲源。需要指出的是，由于球面波源通常在單點多極法中被用于計算球形或近似球形的聲源［5?6，24，26?27］，而本文則是將其應用于波疊加法中，因此為了驗證由球面波源構造的射線波函數(shù)在波疊加法中依然適用于任意形狀聲源的優(yōu)勢，上述3種聲源均為非球形聲源。其中，正四面體輻射體主要考察等效源面和全息面均為封閉且共形，等效源數(shù)量與測點數(shù)量相同情況下的重建效果；兩端帶球帽的圓柱輻射體主要考察等效源面和全息面均為封閉，但僅近似共形，且等效源與測點并非一一對應情況下的重建效果；簡支板聲源則考察等效源面和全息面均為非封閉情況下對空間連續(xù)型結構聲源的重建效果。3個仿真中均添加信噪比為25 dB的高斯白噪聲，并使用Tikhonov正則化方法求解，正則化參數(shù)采用L曲線進行選擇。同時，由前文所述，傳遞矩陣的病態(tài)性會導致正則化對聲場細節(jié)信息的濾除增加，從而在重建高頻復雜聲場時精度下降，因此3個仿真均在1000 Hz以上的較高頻下進行。

2.2 仿真1：正四面體輻射體的聲壓重建

如圖7（a）所示為一幾何中心位于坐標原點O，底面的一邊與y軸平行且邊長為a=1 m的正四面體輻射體。假設該輻射體表面的聲壓由置于其中心O處的點聲源S1和z軸正方向上與原點O距離為d=0.3 m的點聲源S2疊加產(chǎn)生。為了增加聲場的復雜度，將兩個點聲源設置為不同的多極子聲源，其具體的表達式為：

式中 坐標變量rS1，θS1，?S1，rS2，θS2，?S2均為以點源所在位置為坐標原點、坐標軸方向與全局坐標系平行的局部坐標系中的位置標量。

仿真中采用三角形單元均勻劃分四面體表面，共計340個節(jié)點。測量面設置為包裹輻射體且所有邊長均為aH=1.2 m的共形四面體面，測點的數(shù)量和分布方式與聲源表面節(jié)點相同，如圖7（b）所示。等效源布置在聲源表面向內(nèi)以0.5比例縮進的虛擬面上，等效源的數(shù)量和分布方式與聲源表面節(jié)點相同。由于該模型中等效源與全息測點一一對應，因此設置射線波函數(shù)指向各等效源的對應測點。重建波數(shù)設置為k=35（1910.6 Hz）。

圖8為正四面體聲源的表面聲壓幅值分布。由圖8可見，即便已經(jīng)采用了正則化方法進行求解，除了n=4階波函數(shù)的重建結果與解析聲壓吻合較好以外，其余階波函數(shù)的重建聲壓均產(chǎn)生了明顯偏差。為了量化重建結果，圖9給出了各階波函數(shù)對應的重建聲壓相對誤差曲線及傳遞矩陣條件數(shù)曲線，其中，相對誤差由下式計算：

式中 pn表示第n階波函數(shù)的重建聲壓向量；p表示解析聲壓向量。

由圖9（a）可以發(fā)現(xiàn)，當波函數(shù)階數(shù)n=0和n=1時（即傳統(tǒng)波疊加法），重建誤差較大，超過了20%。而且當波函數(shù)的階數(shù)在0～4之間時，重建誤差基本呈下降趨勢，并在4階時最小，誤差小于10%。通過對比圖9（b）的條件數(shù)曲線可知，重建誤差下降是由于傳遞矩陣條件數(shù)逐漸減少，因此正則化后損失的聲場信息隨之減少，重建結果更為精確和穩(wěn)定。而當波函數(shù)的階數(shù)大于4階時，傳遞矩陣的條件數(shù)雖然仍保持下降趨勢，但重建誤差卻反而逐漸上升。這是因為過高的階數(shù)會使得射線波函數(shù)的指向性太強，并由此導致各等效源輻射的聲波過于集中在其對應主指向的射線束內(nèi)，而在非主指向上將因波函數(shù)衰減太快而產(chǎn)生較大誤差，因此無法正確描述真實聲場。由此可見，射線波函數(shù)的指向性并不是越強越好，而是應在合理的范圍內(nèi)選取。一個簡單的方法是利用“輔助面法［24］”進行選取，即在聲源面與全息面間設置若干輔助測量點，然后計算各階射線波函數(shù)（通常只需要計算到前10～15階）在輔助測點處的重建聲壓，最后將重建聲壓與該點的測量聲壓進行對比，取兩者誤差為最小時所對應的階數(shù)即可。由于篇幅所限，本文不展開討論波函數(shù)階數(shù)的選取。實際上，由球面波源構造而成的射線波函數(shù)作為Helmholtz方程在球坐標系下的基本解，其階數(shù)的選擇需要綜合考慮聲源特性、NAH模型配置及噪聲性質等［28］。

2.3 仿真2：兩端帶球帽的圓柱輻射體的聲壓重建

以坐標原點O為幾何中心設置一兩端帶球帽的圓柱殼長條輻射體，其中圓柱部分長度為2a=1 m，圓柱和球帽的半徑均為a=0.5 m，如圖10（a）所示。假設該輻射體表面聲壓由分別置于兩個球帽球心處的點聲源S1和S2疊加產(chǎn)生。與仿真1類似，將兩個點聲源設置為多極子，其中：

式中的坐標變量說明見仿真1。仿真中仍采用三角形單元劃分聲源表面，共計629個節(jié)點，如圖10（b）所示。測量面設置為包裹聲源的長方體表面，其長寬高尺寸為2.1 m×1.1 m×1.1 m，在長方向等間隔分布17個測點，寬和高方向等間隔分布9個測點，測點數(shù)共計642個，如圖10（a）所示。等效源布置在聲源表面向內(nèi)以0.5比例縮進的虛擬面上，等效源的數(shù)量和分布方式與聲源表面節(jié)點相同。注意，該模型中等效源與全息測點的數(shù)量和分布方式均不同，因此不存在一一對應關系，射線波函數(shù)的指向按照圖3（a）所示進行設置。重建波數(shù)設置為k=25（1364.8 Hz）。

與仿真1類似，該仿真中分別給出了各階波函數(shù)的聲源表面重建聲壓幅值分布云圖、聲壓相對誤差曲線及傳遞矩陣條件數(shù)曲線，如圖11和12所示。由圖12（a）可以發(fā)現(xiàn)，與仿真1不同，該重建模型中最小誤差對應的波函數(shù)階數(shù)為n=5。由圖12（b）的條件數(shù)曲線可以發(fā)現(xiàn)，雖然射線波函數(shù)并未指向等效源的對應測點，但傳遞矩陣的條件數(shù)仍然呈現(xiàn)良好的下降趨勢，這表明本文在1.3節(jié)提出的射線波函數(shù)指向設置方法是有效的。

2.4 仿真3：矩形簡支板輻射聲壓的重建

在波疊加法求解聲場外問題的理論中，通常要求給定邊界條件的全息面為一封閉曲面，但在實際工程中，重建由板、殼等結構振動產(chǎn)生的聲場是一種常見的應用場景，此時一般采用非封閉的平面全息面測量聲場信息。這一節(jié)的仿真將利用一個被無限大剛性障板圍繞的簡支板聲源驗證本文方法在全息面為非封閉平面時的聲場重建效果。

如圖13所示為一左下角位于坐標原點、長方向和寬方向分別與y軸和x軸重合的長方形鋼板。長寬厚尺寸設置為1 m×0.5 m×0.003 m，楊氏模量為E=2.1×1011 Pa，泊松比為μ=0.23，密度為ρ1=7.8×103 kg/m3。邊界條件設置為四邊簡支，并在板的（0.25 m，0.75 m）位置處施加一幅值為1 N，波數(shù)為k=63（3439.2 Hz）的簡諧激勵力。全息測量面的大小與簡支板相同并布置在其正上方0.1 m處，測點在長方向和寬方向的掃描間隔均為0.02 m，共1326個測量點。等效源面位于簡支板正下方0.1 m處，其分布方式與測點相同。射線波函數(shù)的指向按圖3（b）中的平面情形進行設置。仿真中重建簡支板上方0.05 m處的聲壓，并與解析聲壓對比。簡支板的解析聲壓是利用Rayleigh積分計算得到，詳細的推導及表達式見文獻［25］。

該仿真中除了給出與仿真1和仿真2相同的重建聲壓幅值分布云圖、聲壓相對誤差曲線及傳遞矩陣條件數(shù)曲線外，為了突出各階波函數(shù)對簡支板各區(qū)域聲場的重建效果，還給出了重建面上的相對誤差分布云圖，上述結果如圖14～16所示。從圖14的聲壓幅值分布云圖上看，各階波函數(shù)的重建聲壓均與解析聲壓吻合得較好。但通過對比圖15中的相對誤差云圖則可以發(fā)現(xiàn)，當波函數(shù)階數(shù)為n=0和n=1時（即傳統(tǒng)波疊加法），重建誤差云圖中有大面積誤差≥15%的深紅色。而當波函數(shù)為n=2～8時，云圖中深紅色的大誤差區(qū)域顯著減少，并主要分布在解析聲壓接近于0處和簡支板邊緣處，這是由于計算相對誤差時分母接近0對誤差產(chǎn)生的放大作用和聲場信息泄露所導致的，其余大部分區(qū)域均為誤差0%～5%的深藍色和淺藍色。由圖16（a）可以發(fā)現(xiàn)，該重建模型的最小誤差對應的波函數(shù)階數(shù)為n=4，最小誤差在5%左右，且階數(shù)在n=2～8的范圍內(nèi)重建誤差均相差不大，都在5%～10%之間。從圖16（b）中可以看到，傳遞矩陣的條件數(shù)曲線仍然呈下降趨勢。該仿真表明，對于非封閉的重建模型，本文方法仍具有一定優(yōu)勢。

3 結 論

針對波疊加法傳遞矩陣的病態(tài)問題，在射線波函數(shù)法的基礎上，提出了一種由（n， 0）階正交球面波源構造而成的射線波函數(shù)，并改進了射線波函數(shù)的主指向設置。文中對所提方法進行了詳細的推導和闡述，并通過數(shù)值仿真進行了驗證，主要結論有：

（1）由（n， 0）階正交球面波源構成的射線波函數(shù)，其指向性隨著階數(shù)n的增大逐漸增強，并可利用Hankel函數(shù)和Legendre多項式的遞推形式方便地獲得任意階表達式，相較于以往格林函數(shù)導數(shù)型的射線波函數(shù)，其計算難度和復雜程度大大降低；

（2）在射線波函數(shù)法中，射線波函數(shù)的主指向無需嚴格指向等效源的對應測點。只需各射線波函數(shù)的主指向在空間中均勻向外分散，即可有效地減少傳遞矩陣條件數(shù)，并保證聲場重建精度，大大提高了射線波函數(shù)法在應用中的靈活性；

（3）在1000 Hz以上的含噪聲高頻聲場中，由于聲場細節(jié)信息豐富，即便采用Tikhonov正則化方法求解，傳統(tǒng)波疊加法仍無法保證重建精度，誤差為20%左右，而射線波函數(shù)法則能獲得更高的重建精度，誤差在5%～10%之間。

參考文獻

1Williams E G， Maynard J D， Skudrzyk E. Sound source reconstructions using a microphone array［J］. The Journal of the Acoustical Society of America， 1980， 68（1）： 340?344.

2Veronesi W A， Maynard J D. Digital holographic reconstruction of sources with arbitrarily shaped surfaces ［J］. The Journal of the Acoustical Society of America， 1989， 85（2）： 588-598.

3BAI M. Application of BEM?based acoustic holography to radiation analysis of sound sources with arbitrarily shaped geometries［J］. The Journal of the Acoustical Society of America， 1992， 92（1）： 533-549.

4Cunefare K A， Koopmann G， Brod K. A boundary element method for acoustic radiation valid for all wavenumbers ［J］. The Journal of the Acoustical Society of America， 1989， 85（1）： 39?48.

5Ochmann M. The source simulation technique for acoustic radiation problems［J］. Acta Acustica United with Acustica， 1995， 81（6）： 512?527.

6Ochmann M. The full?field equations for acoustic radiation and scattering［J］. The Journal of the Acoustical Society of America， 1999， 105（5）： 2574?2584.

7Koopmann G H， Song L M， Fahnline J B. A method for computing acoustic fields based on the principle of wave superposition ［J］. The Journal of the Acoustical Society of America， 1989， 86（6）： 2433?2438.

8Du B K， Zeng X Y， Vorl?nder M. Multizone sound field reproduction based on equivalent source method［J］. Acoustics Australia， 2021， 49（2）： 317-329.

9Chen J， Liu C， Zhang X Z， et al. An approach for indoor prediction of the pass?by noise of a vehicle based on the time?domain equivalent source method［J］. Mechanical Systems and Signal Processing， 2021， 146：107037.

10劉延善， 曾向陽， 王海濤. 封閉空間聲場重構的多層等效源法［J］. 聲學學報， 2020， 45（3）： 367?376.

Liu Yanshan， Zeng Xiangyang， Wang Haitao. 3D sound field reconstruction for the enclosed cavity using the multilayer equivalent sources method［J］. Acta Acustica， 2020， 45（3）： 367?376.

11Valdivia N P. Advanced equivalent source methodologies for near?field acoustic holography［J］. Journal of Sound and Vibration， 2019， 438： 66?82.

12BAI M R， CHEN C C， Lin J H. On optimal retreat distance of equivalent source method?based near?field acoustical holography［J］. The Journal of the Acoustical Society of America， 2011， 129（3）： 1407?1416.

13Williams E G. Regularization methods for near?field acoustical holography［J］. The Journal of the Acoustical Society of America， 2001， 110（4）： 1976?1988.

14Valdivia 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］. The Journal of the Acoustical Society of America， 2006， 120（6）： 3694?3705.

15Gounot Y J R， Musafir R E. On appropriate equivalent monopole sets for rigid body scattering problems［J］. The Journal of the Acoustical Society of America， 2007， 122（6）： 3195?3205.

16Gounot Y J R， Musafir R E. Genetic algorithms： a global search tool to find optimal equivalent source sets ［J］. Journal of Sound and Vibration， 2009， 322（1?2）： 282?298.

17Gounot Y J R， Musafir R E. Simulation of scattered fields： some guidelines for the equivalent source method ［J］. Journal of Sound and Vibration， 2011， 330（15）： 3698?3709.

18WU S W， XIANG Y. Location optimization of monopole equivalent sources in wave superposition method ［J］. International Journal of Acoustics and Vibration， 2018， 23（2）： 254?263.

19LEE S. Review： the use of equivalent source method in computational acoustics［J］. Journal of Computational Acoustics， 2017， 25（1）： 1630001.

20張永斌. 基于等效源法和質點振速測量的近場聲全息技術［D］. 合肥： 合肥工業(yè)大學， 2010.

Zhang Yongbin. Research on nearfield acoustic holography based on the equivalent source method and the measurement of particle velocity［D］. Hefei： Hefei University of Technology， 2010.

21Fernandez?Grande E， Xenaki A， Gerstoft P. A sparse equivalent source method for near?field acoustic holography［J］. The Journal of the Acoustical Society of America， 2017， 141（1）：532?542.

22向宇， 石梓玉， 陸靜， 等. 基于波疊加法的非共形近場聲全息波函數(shù)的構造與選擇［J］. 振動與沖擊， 2020， 39（15）： 183?192.

Xiang Yu， Shi Ziyu， Lu Jing， et al. Construction and selection of nonconformal near?field acoustic holography wave function based on wave superposition method［J］. Journal of Vibration and Shock， 2020， 39（15）： 183?192.

23石梓玉， 向宇， 陸靜， 等. 一種提高聲場重構穩(wěn)定性的射線等效源法［J］. 廣西科技大學學報， 2019， 30（3）： 1?7.

Shi Ziyu， Xiang Yu， Lu Jing， et al. A ray equivalent source method for improving the stability of acoustic field reconstruction［J］. Journal of Guangxi University of Science and Technology， 2019， 30（3）： 1?7.

24畢傳興， 陳心昭， 陳劍， 等. 正交球面波源邊界點法及其在聲全息中的應用［J］. 科學通報， 2004， 49（13）： 1322-1331.

Bi Chuanxing， Chen Xinzhao， Chen Jian， et al. Orthogonal spherical wave source boundary point method and its application to acoustic holography［J］. Chinese Science Bulletin， 2004， 49（13）：1322-1331.

25Williams E G. Fourier Acoustics： Sound Radiation and Nearfield Acoustical Holography［M］. London： Cambridge University Press， 1999.

26Gomes J， Hald J， Juhl P， et al. On the applicability of the spherical wave expansion with a single origin for near?field acoustical holography［J］. The Journal of the Acoustical Society of America， 2009， 125（3）： 1529?1537.

27畢傳興， 胡定玉， 徐亮， 等. 采用球面波疊加法還原自由聲場的方法［J］. 聲學學報， 2014， 39（3）： 339?346.

Bi Chuanxing， Hu Dingyu， Xu Liang， et al. Recovery of the free field in a noisy environment by using the spherical wave superposition method［J］. Acta Acustica， 2014， 39（3）： 339?346.

28Wu H J， Jiang W K， Zhang H B. A mapping relationship based near?field acoustic holography with spherical fundamental solutions for Helmholtz equation［J］. Journal of Sound and Vibration， 2016， 373： 66?88.

Wave superposition method with ray wave functions for near?field acoustic holography based on orthogonal spherical wave sources

SHI Zi?yu 1，2 ?XIANG Yu 1，2 ?LU Jing 1，2WANG Yu?jiang 1，2

1. Guangxi Key Laboratory of Automobile Components and Vehicle Technology， Guangxi University of Science and Technology， Liuzhou 545006， China；

2. School of Mechanical and Automotive Engineering， Guangxi University of Science and Technology， Liuzhou 545006， China

Abstract In the near-field acoustic holography based on the traditional wave superposition method， the monopole radiating the spherical wave is often used as the equivalent source， which is easy to cause the ill conditioned transfer matrix. Replacing the spherical wave function with the ray wave function can effectively improve the ill conditioned transfer matrix. However， the previous ray wave function method uses the directional derivative of Green's function as the wave function， which has complex analytical expression and low computational efficiency. In addition， the wave function direction setting of the previous methods has high requirements for the node distribution， which limits the flexibility of its application. In order to solve the above problems， a series of ray wave functions are reconstructed by using （n， 0） order spherical wave source， this ray wave function can be easily calculated to any order by using the recursive form of spherical Hankel function and Legendre polynomial， and the efficiency is greatly improved. By improving the main direction setting of ray wave function to make it more flexible in practical use， a ray wave function wave superposition method based on orthogonal spherical wave source is proposed. Three numerical simulations， i.e.， radiator of regular tetrahedral， radiator of cylinder with two spherical caps and simply supported rectangular steel plate acoustic source， are performed to verify the sound field reconstruction effect in both traditional method and proposed method. The simulation results show that in the high frequency with abundant sound field details， due to the serious ill conditioned transfer matrix， the traditional wave superposition method cannot guarantee the reconstruction accuracy in the three simulations even with the help of regularization method， and the error is about 20%. The established ray wave function method effectively reduces the condition number of the transfer matrix and improves the ill condition of the system. Therefore， higher reconstruction accuracy is obtained， and the error is between 5% and 10%， which shows that the method in this paper is more superior.

Keywords near?field acoustic holography; wave superposition method; orthogonal spherical wave source; ray wave function