基于壓縮感知和快速波疊加譜方法的近場聲全息

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

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

0 引 言

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

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

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

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

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

其中，ε表示一個和噪聲相關的誤差項。

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

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

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

圖2 場點、源點位置示意圖Fig.2 Schematic diagram of field point and source point

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

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

3 數值算例

3.1 脈動球源聲場重建驗證

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

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

式中：ρ為介質密度，c為介質中的聲速。

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

圖3 結點劃分圖Fig.3 Node partition diagram

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

圖4 f=1 000 Hz時，不同采樣點數下本文方法和壓縮等效源法的重建誤差曲線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范圍內的聲場重建效果進行了比較。算例中取 14個全息采樣點，同樣在測量聲壓時添加信噪比為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可以看出，隨著頻率的增加，兩種方法的重建誤差都隨之增大，但在整個計算頻率內，本文方法的重建誤差均低于壓縮等效源法，并能夠保證準確地求解結果，體現了本文方法的優越性。

3.2 長條形聲源的聲場重建

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

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

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

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

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

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

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

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

圖9 f=1 500、2 500 Hz，采樣點數為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 結 論

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

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

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

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

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