基于稀疏貝葉斯學習的高分辨率Patch近場聲全息

扈 宇， 胡定玉， 方 宇， 肖 悅

(1. 上海工程技術大學 城市軌道交通學院，上海 201620；2. 南昌工程學院 江西省精密驅動與控制重點實驗室，南昌 330099)

基于空間二維Fourier變換的近場聲全息[1-2](Near-field Acoustic Holography, NAH)技術能夠充分利用包含高空間頻率的倏逝波信息實現高分辨率的聲場重建，是一種重要的聲場可視化與噪聲源識別工具。但NAH的極限分辨率取決于全息面上的測點間隔，因此要獲得較高的分辨率，就需要極高的測量成本，這極大影響了該技術的實際應用價值。針對這種問題，Harris等[3]提出了利用聲壓梯度為指導的全息聲壓的三次Hermite插值法提高聲場重建的分辨率，但該方法需要測量質點振速信息，增大了測量工作量。Xu等[4]提出了利用帶限信號恢復算法的NAH分辨率增強方法，但該方法需要經過多次迭代計算，增加了計算時間。隨后，徐亮等[5]提出了一種基于正交球面波插值的方法，但是對于形狀較為復雜的聲源插值效果并不理想。張小正等[6]通過波疊加法原理構建虛源面對全息數據插值，一定程度緩解了聲源形狀對重建結果的影響。最近，毛榮富等[7]利用支持向量回歸算法實現了全息數據插值，有效地減少了測量工作量。孫超等[8]提出了一種基于極限學習機的全息數據插值方法，通過神經網絡訓練數據提高全息分辨率并對全息數據外推。

基于空間二維Fourier變換的NAH技術在應用中遇到的另一個問題是測量孔徑對NAH技術重建精度的影響[9-10]。由于NAH的理論是建立在測量面為無窮大的基礎上，因此當測量孔徑較小時，窗函數的截斷效應會導致波數譜泄露等一系列影響重建精度的問題。Saijyou等[11]在研究大尺寸聲源重建問題時通過全息面補零和波數域反復迭代濾波有效抑制了測量孔徑的影響。Williams等[12]進一步提出了基于空間二維Fourier變換的Patch NAH，以解決小測量孔徑下的全息重建問題。徐亮等[13-14]在波數域進行全息數據外推并通過極小化加權范數優化了外推方法，Scholte等[15-16]對重建過程的低通濾波方法進行優化并利用線性邊界預測法實現了空間域上的全息數據外推，有效提高了小孔徑全息數據的重建精度。在全息數據外推方法不斷發展的同時，基于波疊加法[17]、分布源邊界點法[18]、統計最優法[19]與邊界元法[20]的Patch NAH技術相繼被提出，擴大了Patch NAH的適用范圍。

本文針對基于空間二維Fourier變換的NAH重建精度受到測量孔徑影響較大的問題，提出一種基于稀疏貝葉斯學習(Sparse Bayesian Learning, SBL)算法的高分辨率Patch NAH方法。該方法利用測點的空間坐標構建高斯核函數，通過SBL算法構建全息面聲壓擬合模型實現全息面聲壓數據的插值和外推，并進一步將插值和外推的數據用于NAH重建，有效抑制了測量孔徑與測點間隔對NAH的限制。本文首先給出基于SBL的全息面聲壓插值與外推理論模型，然后通過數值仿真對該方法用于Patch NAH重建的有效性、抗噪性與分辨率提升等性能進行分析，最后通過實驗進一步驗證該方法的有效性。

1 基于空間二維Fourier變換的NAH

基于空間二維Fourier變換的NAH的基本原理是利用均勻分布在靠近聲源表面的全息面上有限個測點的聲場信息，通過全息面和重建面聲壓在波數域的傳遞關系重建聲源輻射區域任意平面的聲壓及其他聲學物理量。如圖1所示，波數域聲壓從全息面到重建面的傳遞關系為

P(kx,ky,zr)=P(kx,ky,zh)e-ikz(zh-zr)

(1)

(2)

由于在實際計算中采用了快速Fourier變換算法，且測量孔徑與測量點數有限，因此重建分辨率受到測量點數的影響。另外將式(2)分解可得到

(3)

在重建過程中通常會利用窗函數抑制噪聲的放大，其中最常用的經典窗函數是Maynard等[21]提出的指數窗，其形式為

(4)

圖1 聲源面、重建面與全息面分布圖Fig.1 A diagram of NAH

2 基于SBL算法的全息面聲壓插值與外推

2.1 高斯核函數

高斯核函數是徑向基函數的一種，二維高斯核函數的數學表達式為

(5)

2.2 基于SBL算法的全息面聲壓插值

高斯核函數可以將低維數據映射到高維并對其進行線性求解，在升維過程中保證了以高斯核函數為基的解的稀疏性，理論上可以通過較少測點擬合聲壓在全息面上的分布，從而實現超高分辨率[24]。由于在全息重建中使用的是復聲壓，且復聲壓實部與虛部的誤差雖都服從零均值高斯分布，但方差可能會有不同，因此可以分別對復聲壓實部與虛部插值并將兩者合并獲取復聲壓。下面僅以實部插值為例介紹全息面聲壓插值過程。

通過高斯核函數映射后，全息面上任意一點ri處的聲壓實部理論值可由高斯核函數表示為

(6)

pre=Φw+ε

(7)

(8)

(9)

式中：α=[α1,α2,…,αN]T，αj為權重系數wj的超參數。

在α與β2皆為未知的情況下， 通過測量的聲壓求解權重系數w的問題即可轉化為最大化聯合后驗概率Prob(w,α,β2pre)的問題， Prob(w,α,β2pre)最大值對應的w即為最優的權重系數。 但由于Prob(w,α,β2pre)很難通過積分直接獲取， 因此將Prob(w,α,β2pre)分解為

Prob(w,α,β2pre)=
Prob(wpre,α,β2)Prob(α,β2pre)

(10)

式中：w的后驗概率分布Prob(wpre,α,β2)可以根據貝葉斯定理表示為

(11)

式中： 分母Prob(preα,β2)可由似然函數式(8)對權重系數w進行邊緣積分得到，即

(12)

式中：A=diag(α1,α2,α3,…,αN)，C=β2I+ΦA-1ΦT。

將式(8)、式(9)和式(12)代入式(11)可得

Prob(wpre,α,β2)=N (μ,Σ)

(13)

式中：Σ=(β-2ΦTΦ+A)-1，μ=β-2ΣΦTpre， Prob(wpre,α,β2)滿足均值為μ、方差為Σ的N維高斯分布。

得到w的后驗分布Prob(wpre,α,β2)后， 稀疏貝葉斯學習過程即被轉化為超參數后驗分布Prob(α,β2pre)∝Prob(preα,β2)Prob(α)Prob(β2)關于α與β2的最大值問題。 在貝葉斯模式中， Prob(preα,β2)被稱為邊緣似然函數。 對邊緣似然函數進行第Ⅱ類型最大似然參數估計可求出α與β2的估計值， 即利用式(12)分別對α與β2求偏導并使導數等于零， 可求出α與β2的更新公式

(14)

(15)

2.3 基于SBL算法的全息面聲壓外推

全息面聲壓外推過程和插值過程類似，都是利用高斯核函數實現，區別在于外推點的位置和插值點位置的不同。由于外推點均位于測量孔徑外部，而全息面聲壓是采用以測點位置為中心峰值位置的高斯核函數進行擬合，因此根據高斯核函數的衰減特性，測量孔徑外的外推點聲壓平滑衰減，最終可得到邊緣為零或無限近似于零的全息面聲壓數據，因此利用外推后的數據進行重建，可以極大避免窗函數截斷效應的影響。需要說明的是，在應用過程中可以同時進行全息面聲壓的插值和外推，可以在幾乎不增加計算時間的同時提高重建結果的分辨率。

3 數值仿真

(16)

式中：Pcal為插值或重建的聲壓，Pth為對應的理論值。

3.1 全息面聲壓插值

為便于描述，定義插值倍數為插值前后分辨率的比值。對全息面聲壓進行5倍插值，即插值后全息面聲壓的分辨率為1 cm。圖2顯示的是頻率為900 Hz時全息面上的測量聲壓(圖2(a))、插值后的聲壓(圖2(a))及其理論值(圖2(c))。從圖中可以看出，插值結果與理論值幾乎完全一致，利用式(16)計算出插值后全息數據的插值誤差僅為1.88%，證明了本文插值模型的有效性。

圖2 900 Hz時全息面聲壓插值效果對比圖Fig.2 The pressure on the hologram at 900 Hz

為研究插值倍數與插值精度之間的關系，對全息面聲壓進行不同倍數的插值。頻率為900 Hz時插值誤差與插值倍數的關系如圖3所示。從圖3中可以看出，雖然插值點間距縮小使得全息數據點數成倍增長，但插值誤差的變化很小，說明插值點數對于插值結果的影響很小。圖3同時也顯示了全息面聲壓的測量數據與理論值之間的相對誤差。比較測量誤差和插值誤差可以看出，基于SBL算法的全息數據插值方法對測量數據中的噪聲有一定的抑制能力。圖4給出了頻率為500 Hz和900 Hz時，在不同的信噪比(SNR)的情況下插值誤差與測量誤差的對比。從圖中可以看出，在測量條件較差的情況下，使用插值模型可以較好地抑制噪聲，證明了插值模型具有較高的工程實用性。

圖3 插值誤差與插值倍數的關系Fig.3 The relative errors of the measurement and theinterpolation versus the ratios of the spatial resolutionsbefore and after interpolation

圖4 插值誤差與測量誤差隨信噪比變化曲線Fig.4 The relative errors of the measurement and theinterpolation versus the SNR

3.2 高分辨率聲場重建

對全息面聲壓進行5倍插值并同時進行外推，并利用插值和外推后的全息面聲壓重建距板表面0.02 m處的聲壓。圖5顯示的是頻率為900 Hz時重建結果與理論值的比較。從圖中可以看出，重建的聲壓與理論值吻合地很好，通過式(16)計算兩者之間的誤差僅為3.32%，由此可見，對全息面聲壓進行插值和外推處理可以有效地避免有效孔徑效應對重建精度的影響，同時可以獲取高分辨率的重建結果。

圖6顯示的是在100～2 000 Hz的頻率范圍內聲壓的重建誤差隨頻率的變化。從圖中可以看出，聲壓的重建誤差均在10 %以下，證明了本文方法的有效性。

圖5 900 Hz時重建結果與重建面處實際聲壓幅值分布對比圖Fig.5 The pressure on the reconstruction plane at 900 Hz

為進一步驗證高分辨率Patch HAH對小測量孔徑的適應性，在全息面上取三個不同大小的測量孔徑，其位置如圖7所示。分別對三個孔徑的全息數據進行插值和外推，并進一步利用NAH重建距板表面2 cm處的聲壓。圖8顯示的是頻率為900 Hz時利用三個測量孔徑內聲壓重建的結果與其理論值的比較。從圖中可以看出，三個測量孔徑的邊緣附近均存在較大峰值，測量孔徑產生的截斷效應會對基于二維空間Fourier變換的NAH的重建精度產生很大影響。而利用本文方法對全息面聲壓進行插值和外推后，利用三個測量孔徑的全息數據都獲得了較好的重建結果。其中，三個孔徑的重建誤差分別為孔徑1∶10.95%，孔徑2∶15.53%和孔徑3∶16.35%。證明了本文插值外推方法對較小的測量孔徑仍具有較好的穩定性。

圖7 全息面測量孔徑分布位置圖Fig.7 Positions of the measurement apertures

圖8 900 Hz時各孔徑重建結果對比圖Fig.8 The pressure on the three reconstruction plane at 900 Hz

4 實驗驗證

實驗在半消聲室中進行，半消聲室本底噪聲小于16 dB。實驗裝置實物照片如圖9所示。聲源為80×66 cm2的固支鋼板，鋼板厚度為2 mm。對鋼板施加簡諧信號產生聲場。全息面位于鋼板上方5 cm處，測點間隔為5 cm，共15×11個測點。重建面位于鋼板上方2 cm處。為檢驗重建結果，同時以2.5 cm的測量間隔測量重建面上的聲壓作為其真實值，共31×21個測點。在15×11個測點得到全息面中，選取2個大小不同、測點數不同的測量孔徑并利用孔徑內的數據進行重建，孔徑在全息面上的位置分布如圖10所示。算法迭代過程中判斷收斂條件與仿真相同。

圖9 實驗裝置實物圖Fig.9 A photo of the experimental setup

圖11顯示的是頻率為700 Hz、900 Hz與1 100 Hz時，對兩個測量孔徑的全息數據進行外推并插值后重建結果與真實值的比較。從圖中可以看出，利用本文方法對不同尺寸測量孔徑進行重建皆獲得了較好的重建效果。不同激勵下各測量孔徑重建誤差見表1。重建結果表明，本文方法可以準確地重建出聲場分布，并且有效地提升了重建結果分辨率，證明了利用SBL算法對全息數據插值與外推后可以實現高分辨Patch NAH。

圖10 全息面測量孔徑位置Fig.10 Positions of the measuring apertures

激勵頻率/Hz7009001 100孔徑1重建誤差/%13.8810.3218.69孔徑2重建誤差/%7.7414.1511.94

圖11 不同激勵頻率下各孔徑重建結果對比圖Fig.11 Comparison of reconstruction results with the corresponding true values

5 結 論

本文建立了一種基于SBL的高分辨率Patch NAH方法，解決了測量孔徑與測點數對基于空間二維Fourier變換的NAH重建效果影響較大的問題.研究結果表明：利用高斯核函數和SBL算法對全息面數據進行插值和外推，可以有效地避免測量孔徑對NAH重建結果的影響，同時獲取高分辨率的重建結果。另外，利用插值模型還可以在一定程度上減少測量噪聲，具有較高的工程應用價值。