石梓玉 向宇 陸靜 王玉江

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

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

引 言

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

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

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

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

1 理論部分

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

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

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

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

式中 qΦ（rEj）為第j個(gè)離散等效源的源強(qiáng)；Φ（r，rEj）為該等效源對(duì)應(yīng)的波函數(shù)。

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

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

求解方程（3），得源強(qiáng)向量QG為：

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

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

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

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

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

（1）波函數(shù)計(jì)算效率方面的缺陷。文獻(xiàn)［22?23］均是利用格林函數(shù)的方向?qū)?shù)作為射線波函數(shù)，該類型波函數(shù)的高階導(dǎo)數(shù)（超過(guò)6階后）解析表達(dá)式非常復(fù)雜，導(dǎo)致計(jì)算效率較低甚至無(wú)法計(jì)算。

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

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

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

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

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

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

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

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

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

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

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

利用式（6），可得rEj處單位強(qiáng)度等效源在場(chǎng)點(diǎn)r處輻射的聲壓為：

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

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

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

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

求解出等效源強(qiáng)后，再利用下式即可重建場(chǎng)點(diǎn)r處的聲壓：

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

2 方法驗(yàn)證

2.1 仿真模型說(shuō)明

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

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

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

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

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

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

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

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

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

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

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

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

2.4 仿真3：矩形簡(jiǎn)支板輻射聲壓的重建

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

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

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

3 結(jié) 論

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

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

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

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

參考文獻(xiàn)

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劉延善， 曾向陽(yáng)， 王海濤. 封閉空間聲場(chǎng)重構(gòu)的多層等效源法［J］. 聲學(xué)學(xué)報(bào)， 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張永斌. 基于等效源法和質(zhì)點(diǎn)振速測(cè)量的近場(chǎng)聲全息技術(shù)［D］. 合肥： 合肥工業(yè)大學(xué)， 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向宇， 石梓玉， 陸靜， 等. 基于波疊加法的非共形近場(chǎng)聲全息波函數(shù)的構(gòu)造與選擇［J］. 振動(dòng)與沖擊， 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石梓玉， 向宇， 陸靜， 等. 一種提高聲場(chǎng)重構(gòu)穩(wěn)定性的射線等效源法［J］. 廣西科技大學(xué)學(xué)報(bào)， 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畢傳興， 陳心昭， 陳劍， 等. 正交球面波源邊界點(diǎn)法及其在聲全息中的應(yīng)用［J］. 科學(xué)通報(bào)， 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畢傳興， 胡定玉， 徐亮， 等. 采用球面波疊加法還原自由聲場(chǎng)的方法［J］. 聲學(xué)學(xué)報(bào)， 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