張 陽，向 宇，石梓玉，陸 靜，王玉江

(1. 廣西科技大學(xué)廣西汽車零部件與整車技術(shù)重點(diǎn)實驗室，廣西柳州545006；2. 廣西科技大學(xué)機(jī)械與交通工程學(xué)院，廣西柳州545006)

0 引 言

近場聲全息技術(shù)(Near-Field Acoustic Holography, NAH)是一種噪聲源識別、定位及聲場可視化的有效技術(shù)[1]。經(jīng)過幾十年的發(fā)展，各國學(xué)者相繼提出了空間傅里葉變換算法(Spatial Fourier Transform, SFT)[2]、邊界元法(Boundary Element Method,BEM)[3]、等效源法(Equivalent Source Method,ESM)[4]等。其中，基于等效源法的近場聲全息因原理簡單、適應(yīng)性強(qiáng)等優(yōu)點(diǎn)得到了廣泛應(yīng)用[5-7]。為了保證該方法的重建精度，必須要有足夠數(shù)量的測量點(diǎn)和源強(qiáng)點(diǎn)。但在NAH技術(shù)中，通常情況下是通過測量外域聲場的信息來重建振動體的近場聲學(xué)量，受限于測量成本和測量條件，一般很難獲得測量點(diǎn)數(shù)大于源強(qiáng)點(diǎn)數(shù)的超定方程組。此時若采用基于l2范數(shù)的最小二乘方法求解欠定方程組，難以得到理想的結(jié)果。因此，如何在減少測點(diǎn)的情況下，仍能穩(wěn)定獲得較高的重建精度，一直是NAH技術(shù)重點(diǎn)關(guān)注的問題。

Patch近場聲全息方法[8-9]是解決上述問題的一種有效途徑，而壓縮感知技術(shù)(Compressive Sensing,CS)[10-11]則提供了另外一種手段。該技術(shù)利用信號的稀疏性來解決欠定問題，若信號本身稀疏，或在某個變換域稀疏，即可通過一個與變換基無關(guān)的觀測矩陣，在低于奈奎斯特采樣定理的采樣頻率下高概率地重構(gòu)出原始信號。其實質(zhì)上是一種尋找欠定方程組稀疏解的技術(shù)，由于利用了基于l0或l1范數(shù)的重構(gòu)算法替代傳統(tǒng)基于l2范數(shù)的最小二乘法，因此在減少測點(diǎn)的情況下也能保證較高的重建精度。2012年，Chardon等[12]首次將壓縮感知應(yīng)用到NAH中，并通過實驗證明了該方法在減少測點(diǎn)數(shù)時仍具有較高的計算精度。Fernandez-Grande等[13]將壓縮感知理論與等效源法相結(jié)合，提出了一種壓縮等效源法(Compressive Equivalent Source Method,C-ESM)。Hald等[14]在壓縮等效源法的基礎(chǔ)上，結(jié)合最速下降法提高了計算效率，并對其在高頻段下的聲場重建效果進(jìn)行了研究。畢傳興等[15]提出了一種基于功率阻矩陣特征分解的壓縮模態(tài)等效源法，并通過數(shù)值仿真與實驗證明了該方法在減少測點(diǎn)數(shù)方面優(yōu)于壓縮等效源法。胡定玉等[16]基于等效源法與奇異值分解法，提出了一種用于空間擴(kuò)展聲源的稀疏采樣聲場重建方法。然而，上述方法都是將虛擬源強(qiáng)進(jìn)行單元數(shù)值離散，這將對最終的求解精度造成一定影響，尤其是在高頻聲場的重建時，其對精度的影響更加嚴(yán)重。

作者在前期研究中，將等效源法中的虛擬源強(qiáng)和積分核函數(shù)在軸對稱虛擬面上沿周向和子午線方向利用傅里葉級數(shù)進(jìn)行展開，提出了一種求解軸旋轉(zhuǎn)空穴三維諾伊曼(Neumann)聲輻射問題的快速波疊加譜方法[17-18]。由于該方法是一種傅里葉級數(shù)展開形式的半解析、半數(shù)值方法，因而其計算精度高于傳統(tǒng)的離散源強(qiáng)等效源法，同時作者在進(jìn)一步的研究中發(fā)現(xiàn)，待求的傅里葉展開系數(shù)向量總是一個僅含少量非零值的稀疏向量，該特性正好滿足了壓縮感知理論對稀疏性的要求。因此，本文將文獻(xiàn)[18-19]的快速波疊加譜方法與壓縮感知技術(shù)進(jìn)行有機(jī)結(jié)合，推廣用于NAH技術(shù)中，提出了一種基于壓縮感知和快速波疊加譜方法的近場聲全息方法。文中首先推導(dǎo)等效源強(qiáng)和積分核函數(shù)在虛擬面上的傅里葉級數(shù)展開形式，并分析待求展開系數(shù)向量的稀疏性，然后在此基礎(chǔ)上結(jié)合壓縮感知算法求解，最后通過數(shù)值仿真驗證所提方法在聲場重建中的有效性與優(yōu)越性。

1 基于壓縮感知的等效源法

如圖1所示，考慮一任意形狀的三維振動體，其中S為振動體邊界，E+為振動體外部輻射域，振動體內(nèi)部布置一虛擬曲面SE，其內(nèi)部區(qū)域記為E。

其中，ε表示一個和噪聲相關(guān)的誤差項。

2 結(jié)合壓縮感知的快速波疊加譜方法

文獻(xiàn)[17-18]中，將源強(qiáng)和積分核函數(shù)在軸對稱的虛擬面上采用傅里葉級數(shù)進(jìn)行展開，提出了一種求解軸旋轉(zhuǎn)空穴三維 Neumann聲輻射問題的快速波疊加譜方法，并取得了較好的聲場計算結(jié)果。本文借鑒文獻(xiàn)[17-18]的方法與思路并將其應(yīng)用于近場聲全息的聲場計算中。

假設(shè)等效源點(diǎn)rE和場點(diǎn)r在柱坐標(biāo)系下的位置坐標(biāo)分別為，如圖 2(a)所示。它們在各自子午面內(nèi)的二維坐標(biāo)分別為，如圖2(b)所示。

圖2 場點(diǎn)、源點(diǎn)位置示意圖Fig.2 Schematic diagram of field point and source point

聲波一般是由結(jié)構(gòu)振動激發(fā)聲介質(zhì)產(chǎn)生的，根據(jù)結(jié)構(gòu)動力學(xué)理論，結(jié)構(gòu)在任意載荷激勵下產(chǎn)生的動力響應(yīng)主要由低階模態(tài)疊加而成，因而結(jié)構(gòu)振動輻射的聲波也主要集中在低階模態(tài)。此外，由文獻(xiàn)[20]可知，聲音在傳播的過程中，低階模態(tài)衰減慢、輻射能力強(qiáng)，高階模態(tài)衰減快、輻射能力弱。換言之，隨著傳播距離的增大，聲源對聲場的貢獻(xiàn)主要集中在低階模態(tài)，高階模態(tài)類似于基于空間傅里葉變換的近場聲全息中的倏逝波模態(tài)，階數(shù)越高衰減越快，也就是高階的Cmn越來越小。因此，只要式(21)的求和截斷項數(shù)mf和nf取得足夠大，由該式表示的聲場即包含了輻射能力強(qiáng)的低階模態(tài)和足夠多的倏逝波模態(tài)，因而傅里葉展開系數(shù)向量Cmn必然是一個具有一定稀疏性的稀疏向量，且截斷項數(shù)mf和nf取得越大，其稀疏性越強(qiáng)。

由于Cmn是一個稀疏向量，因此可直接利用壓縮感知技術(shù)進(jìn)行求解，即在最小l1范數(shù)下求解如下問題：

3 數(shù)值算例

3.1 脈動球源聲場重建驗證

在添加噪聲的情況下，利用脈動球源對比不同采樣點(diǎn)數(shù)時本文方法與壓縮等效源法的聲場重建效果。其中，壓縮等效源法與本文方法皆采用SPGL1工具箱[21]中的基追蹤降噪算法(Basis Pursuit De-Noising, BPDN)進(jìn)行求解。

已知脈動球源半徑為rS，周邊產(chǎn)生均勻徑向振速幅值為v0的脈動球源在距離球心 r處所輻射聲壓的解析解[22]

式中：ρ為介質(zhì)密度，c為介質(zhì)中的聲速。

仿真中設(shè)置脈動球的半徑為 1 m，球心在坐標(biāo)原點(diǎn)處。等效源點(diǎn)、全息采樣點(diǎn)分別布置在半徑為0.3、1.1 m且與聲源面同心的球面上。在壓縮等效源法中設(shè)置虛擬等效源離散間隔與全息面一致；在本文方法中取虛擬面周向、子午方向抽樣數(shù)為N=M= 28，周向、子午線方向求和截斷項數(shù)分別為nf=30、 mf= 5 0。全息面、等效源面結(jié)點(diǎn)劃分方式見圖 3。在進(jìn)行計算時，兩種重建方法均添加信噪比為20 dB的高斯白噪聲。重建脈動球表面的114個離散結(jié)點(diǎn)聲壓，并計算其相對誤差。相對誤差 e定義為

圖3 結(jié)點(diǎn)劃分圖Fig.3 Node partition diagram

圖4顯示了頻率為1 000 Hz、采樣點(diǎn)數(shù)范圍為14～314時本文方法和壓縮等效源法的重建誤差曲線曲線上的數(shù)值為對應(yīng)的采樣點(diǎn)數(shù)?？梢钥闯觯?dāng)采樣點(diǎn)數(shù)增加到 62時，壓縮等效源法的重建誤差趨向于 2%。但本文方法的誤差曲線始終低于壓縮等效源法，即使是在 14個采樣點(diǎn)的情況下，其重建誤差仍一直低于 2%。同時，本文方法的重建精度也不依賴于采樣點(diǎn)數(shù)，采樣點(diǎn)數(shù)的增加并未明顯改善本文方法的重建精度，這是因為對于脈動球源，利用本文方法得到的傅里葉展開系數(shù)具有良好的稀疏性，因此 14個采樣點(diǎn)就能夠保證本文方法對于脈動球表面聲壓的重建精確。

圖4 f=1 000 Hz時，不同采樣點(diǎn)數(shù)下本文方法和壓縮等效源法的重建誤差曲線Fig.4 Reconstruction error curves of the proposed method and C-ESM under different sampling points at f=1 000 Hz

如圖 5對本文方法與壓縮等效源法在頻率f=100～3 000 Hz范圍內(nèi)的聲場重建效果進(jìn)行了比較。算例中取 14個全息采樣點(diǎn)，同樣在測量聲壓時添加信噪比為20 dB的高斯白噪聲。

圖5 f=100～3 000 Hz時，本文方法和壓縮等效源法的重建誤差曲線Fig.5 Reconstruction error curve of the proposed method and C-ESM at f=100～3 000 Hz

由圖5可以看出，隨著頻率的增加，兩種方法的重建誤差都隨之增大，但在整個計算頻率內(nèi)，本文方法的重建誤差均低于壓縮等效源法，并能夠保證準(zhǔn)確地求解結(jié)果，體現(xiàn)了本文方法的優(yōu)越性。

3.2 長條形聲源的聲場重建

為進(jìn)一步驗證本文方法的優(yōu)越性，通過重建長條形聲源的表面聲壓，將本文方法與壓縮等效源法進(jìn)行了比較。

仿真中，長條形聲源的長寬高尺寸設(shè)為0.4 m×0.4 m×1.2 m，其表面聲壓由置于內(nèi)部的 16個不同強(qiáng)度的點(diǎn)聲源通過“替代法”生成。這些點(diǎn)聲源布置在尺寸為 0.025 m×0.025 m×0.1 m 的長方體框架上，如圖6所示。全息面設(shè)置為半徑rH=0.8 m的球面，共布置366個測點(diǎn)。當(dāng)使用壓縮等效源法計算時，設(shè)置等效源面與聲源面共形，其尺寸較聲源面縮小為 0.16 m×0.16 m×0.48 m(縮小比例為1:2.5)，并在x和y方向均勻布置6個等效源、z方向均勻布置16個等效源，共352個等效源，如圖6(a)所示。當(dāng)使用本文方法計算時，為使等效源面盡可能地與聲源面共形，將其設(shè)置為與全息面同心的橢球面，橢球面z向的長半軸a取0.3 m，x和y向的短半軸b、c均取0.075 m，如圖6(b)所示。其中，抽樣數(shù)M、N及求和截斷項數(shù)mf、nf均與第3.1節(jié)一致。在仿真計算中，兩種重建方法均添加信噪比為30 dB的高斯白噪聲。

圖6 長條形聲源和兩種方法的等效源布置示意Fig.6 Layout of long strip source and equivalent source of two methods

圖 7給出了兩種方法在頻率 f=100～3 000 Hz范圍內(nèi)重建長條聲源表面聲壓的誤差曲線，其中重建誤差可根據(jù)式(25)計算得到。從圖7中可以看出，當(dāng)采樣點(diǎn)數(shù)足夠時，兩種方法在整個分析頻段內(nèi)均能夠保證精確的重建結(jié)果，且誤差都低于4%。

圖7 采樣點(diǎn)數(shù)為366時，兩種方法的重建誤差隨頻率變化的曲線Fig. 7 The curves of reconstruction error versus frequency for two methods with 366 sampling points

分析在采樣點(diǎn)數(shù)減少的情況下兩種方法的聲場重建效果。圖8給出了在114個采樣點(diǎn)的情況下，兩種方法在頻率 f=100～3 000 Hz范圍內(nèi)重建長條聲源表面聲壓的誤差曲線。由圖8可以看出，在100～2 000Hz頻段內(nèi)，兩種方法均能得到準(zhǔn)確的重建結(jié)果，且重建誤差都低于 20%。隨著頻率的增加、當(dāng)f=2 500～3 000 Hz時，壓縮等效源法的誤差逐漸增大，而本文方法仍保持著相對穩(wěn)定的重建結(jié)果。這是因為本文方法是一種傅里葉級數(shù)展開形式的解析計算方法，只要求和截斷項數(shù)nf、mf取得足夠大，由式(21)所表達(dá)的聲場可無限近似真實聲場，同時也使得該式包含了輻射能力強(qiáng)的低階模態(tài)和足夠多的倏逝波模態(tài)，獲得了足夠的聲場細(xì)節(jié)信息；待求的傅里葉展開系數(shù)向量還具有一定的稀疏性，滿足了壓縮感知理論中的稀疏性要求。因此在相同采樣點(diǎn)數(shù)的聲場重建中(尤其是在高頻時)，本文方法的重建精度優(yōu)于壓縮等效源法的重建精度。

圖8 采樣點(diǎn)數(shù)114時，兩種方法的重建誤差隨頻率變化的曲線Fig. 8 The curves of reconstruction error versus frequency for two methods with 114 sampling points

圖9給出了采樣點(diǎn)數(shù)為 114、頻率分別為1 500 Hz和2 500 Hz時上述兩種方法的重建誤差隨信噪比變化的曲線。從圖9中可以看出，當(dāng)頻率為1 500 Hz時，兩種方法均可以保證較好的重建精度。但當(dāng)頻率增大時(f=2 500 Hz)，即使增大信噪比，壓縮等效源法也得不到很好的重建結(jié)果，但本文方法的誤差曲線在整個信噪比分析范圍內(nèi)(10～50 dB)始終低于壓縮等效源法，并且重建誤差始終保持在20%以下，表明本文方法不僅在相同采樣點(diǎn)數(shù)的高頻聲場計算中獲得了更高的重建精度，而且在不同的信噪比時均具有良好的抗噪性能。

圖9 f=1 500、2 500 Hz，采樣點(diǎn)數(shù)為114時，兩種方法的重建誤差隨信噪比變化的曲線Fig.9 The curves of reconstruction error versus signal to noise ratio for two methods with 114 sampling points at f =1 500 and 2 500 Hz

4 結(jié) 論

本文在等效源法的基礎(chǔ)上，將等效源強(qiáng)和積分核函數(shù)在軸對稱的虛擬面上進(jìn)行雙向傅里葉級數(shù)展開，所建立的聲壓表達(dá)式是一種半解析半數(shù)值的形式，因而相對于其他數(shù)值離散方法具有更高的求解精度。另外，由于待求的傅里葉展開系數(shù)向量還具有一定稀疏性，因此在結(jié)合壓縮感知的求解方法后，即使在少測點(diǎn)的欠定條件下也獲得了滿意的重建結(jié)果。文中利用數(shù)值仿真對比了所提方法與傳統(tǒng)壓縮等效源法的聲場重建效果。結(jié)果表明：

(1) 對于重建脈動球聲源，兩種方法在整個計算頻段下均有較好的重建效果，但所提方法的重建效果略優(yōu)于壓縮等效源法；

(2) 對于重建長條形聲源，兩種方法在 100～1 500 Hz的頻段內(nèi)重建效果相當(dāng)，但在 2 500～3 000 Hz的高頻段時，本文方法的重建精度均優(yōu)于壓縮等效源法；

(3) 在1 500 Hz和2 500 Hz頻率時，本文方法的抗噪性能優(yōu)于壓縮等效源法。

雖然本文僅在軸對稱虛擬面上將源強(qiáng)進(jìn)行展開，但若進(jìn)一步采用雙向正交曲線坐標(biāo)下的傅里葉展開，即可將本文方法推廣到任意形狀封閉聲源的近場聲全息計算中。