張子鑫， 肖友洪， 魏富康

(哈爾濱工程大學 動力與能源工程學院，黑龍江 哈爾濱 150001)

近場聲全息技術(shù)(near-field acoustic holography, NAH)是一種有效的聲場可視化技術(shù)，可以通過聲場中的二維平面上測量得到的聲壓獲得三維空間中的聲學量?；贔ourier變換[1](fast Fourier transform, FFT)、邊界元法[2](boundary element method, BEM)、Helmholtz最小二乘法[3](Helmholtz equation least squares, HELS)等的NAH已經(jīng)得到廣泛研究。然而當全息面孔徑較小時，重建精度受到影響，原有的方法受到限制，局部近場聲全息技術(shù)逐漸發(fā)展起來，如基于Fourier變換[4]和邊界元法[5]的局部近場聲全息技術(shù)。統(tǒng)計最優(yōu)近場聲全息[6](statistical optimal near field acoustic holography，SONAH)同樣也屬于局部近場聲全息技術(shù)，通過平面波與倏逝波的線性疊加重建聲場，突破了傳統(tǒng)近場聲全息技術(shù)的孔徑限制。盡管類似SONAH的局部近場聲全息技術(shù)不再依賴全息面的孔徑，但算法仍受限于Nyquist采樣定理，導致在中高頻時即使采樣面積較小，也需要大量采樣點。

壓縮感知技術(shù)[7]為上述問題提供了解決辦法，將壓縮感知與近場聲全息技術(shù)相結(jié)合進行聲場重建，可實現(xiàn)較少的采樣點得到較好的重建效果[8-9]。本文在壓縮感知和邊界元法結(jié)合的基礎上，運用SONAH進行全息面上聲壓的數(shù)據(jù)外推，并將其與未外推前的重建誤差進行了對比分析。通過仿真與實驗研究，驗證了該組合近場聲全息技術(shù)的有效性，并對影響因素進行分析，確定了獲得較好重建效果時的外推面積大小。

1 組合近場聲全息算法

1.1 邊界元法與壓縮感知

基于邊界元法的近場聲全息技術(shù)的理論基礎為Helmholtz積分方程，若已知振動體表面的聲壓和法向振速，即可得到聲場中任意點處的聲壓。將場點移到聲源表面S，建立聲源表面聲壓pS與法向振速vn之間的關系[10]：

HSpS=GSvn

(1)

利用Helmholtz積分方程與式(1)建立場點聲壓pf與振動體表面法向振速vn之間的關系：

pf=Gvn

(2)

vn=(GHG+λLHL)-1GHpf

(3)

式中λ為正則化參數(shù)，通常采用L曲線法[11]，廣義交叉驗證法[12]等方法進行選取。由全息面聲壓得到振動體表面振速后，利用式(2)即可得到聲場中任意點處的聲壓。

壓縮感知理論應用的前提是信號是稀疏的，或者在某一個基函數(shù)下的分解系數(shù)是稀疏的，可以通過一個與該基函數(shù)不相干的觀測矩陣以遠低于Nyquist采樣定理的采樣頻率獲取信號的全部信息。在邊界元法中表面振速是需要獲取的信號，選用聲輻射模態(tài)作為聲源表面振速的基函數(shù)[13-14]，其最早被用來分析輻射聲功率。聲源輻射聲功率表示為：

(4)

(5)

(6)

式中矩陣U=[u1u2…uN]即可作為表面振速的基函數(shù)，此時，表面振速表示為：

(7)

式中：N為聲源表面節(jié)點數(shù)，此時表面振速的求解過程可轉(zhuǎn)化為先對分解系數(shù)α進行求解，因為聲源表面振速在基函數(shù)U下的分解系數(shù)α是K稀疏的，且K?N，所以可通過稀疏正則化進行求逆，將其轉(zhuǎn)化為最小l0范數(shù)進行求解：

arg min‖α‖0， s.t.pf=GUα

(8)

同樣，可以采用更加簡單的最小l1范數(shù)代替最小l0范數(shù)，直接利用CVX工具箱進行求解。求得α后即可通過式(7)得到聲源表面振速，再利用式(2)得到聲場中任意點處的聲壓。

1.2 SONAH外推

SONAH的前提是將聲場中任意點聲壓看作無數(shù)個空間波數(shù)域的平面?zhèn)鞑ゲê唾渴挪ǖ寞B加[15]，其數(shù)學表示為：

(9)

式中：M表示離散后的波數(shù)矢量總數(shù)；km=(kxm,kym,kzm)為第m階的波數(shù)矢量；Φ(km,r)=exp[i(kxmx+kymy+kzmz)]為第m階單元平面波。

SONAH外推的模型示意圖如圖1所示，假設全息面h上有n個測量點，外推后的全息面H上有N個測量點，根據(jù)波疊加原理，H上任意點Φ(km,rE)處的單元平面波Φ(km,rE)可由h上同階的單元平面波線性疊加得到：

圖1 模型布置Fig.1 Diagram of model layout

(10)

式中：ci(rE)為全息面h上第i個測量點處各階單元平面波的權(quán)重系數(shù)。離散后的波數(shù)矢量總數(shù)為M，由此可得到M個線性方程，組成矩陣形式為：

b=Ac(rE)

(11)

通過正則化方法，可得到系數(shù)向量的解為：

c(rE)=(AHA+λLHL)-1AHb

(12)

聯(lián)合式(9)、(10)和式(12)可以將全息面H上一點rE用全息面h上各點聲壓表示為：

p(rE)=phc(rE)=ph(AHA+λLHL)-1AHb

(13)

本文提出的組合近場聲全息方法，即基于外推后的全息面H上各點的聲壓，并利用邊界元法結(jié)合壓縮感知，實現(xiàn)了較少場點情況下的聲場重建。

2 板聲場重建仿真

仿真中，選用0.8 m×0.8 m×0.003 m的鋁板作為聲源，邊界施加固支約束，以鋁板中心為坐標原點建立模型，分別在(-0.1,-0.18,-0.001 5) m和(0.03,0.165,-0.001 5) m處施加大小為10 N的激勵力。全息面h選在距原點高0.04 m處，大小為0.64 m×0.64 m，離散后的點數(shù)為81個，陣元間距為0.08 m。重建面S選在距離原點高0.02 m處，大小為0.8 m×0.8 m，重建點個數(shù)為121個。

為了定量評價該方法的結(jié)果，誤差被定義為：

e=‖pr-pt‖2/‖pt‖2

(14)

式中：pr表示算法得到重建面的聲壓；pt表示仿真得到的重建面聲壓理論值。

2.1 所提方法的驗證

全息面h上場點數(shù)為9×9，外推后全息面H大小為0.8 m×0.8 m，場點數(shù)為11×11，利用外推后全息面H上的場點聲壓，通過邊界元法進行聲場重建，分別比較了邊界元法與傳統(tǒng)正則化方法和壓縮感知方法相結(jié)合的重建誤差，并將結(jié)果與未外推前的重建誤差進行比較。結(jié)果如圖2所示，圖中TR與CS分別為傳統(tǒng)正則化方法與壓縮感知方法的縮寫。

圖2 采用組合近場聲全息前后重建誤差的比較Fig.2 Comparison of reconstruction errors before and after using combined near-field acoustic holography

場點選擇81個時，使用壓縮感知與邊界元法相結(jié)合的重建誤差整體小于傳統(tǒng)正則化方法，說明前者對于較少場點的重建效果較好。但是隨著頻率的增加，4種情況下的重建誤差都隨之增加，這是因為頻率的增加導致波長變短，倏逝波衰減迅速，聲場信息減少，從而降低了重建精度。比較利用外推前后的聲壓重建聲場的誤差，可以發(fā)現(xiàn)經(jīng)過外推，重建誤差明顯降低，且在整個分析頻段內(nèi)重建誤差均小于30%，驗證了外推方法與邊界元法結(jié)合后的算法的有效性。但是從圖中可以看出，利用外推后的場點聲壓進行聲場重建時，傳統(tǒng)正則化方法與壓縮感知方法的重建誤差十分接近，而且在個別頻率傳統(tǒng)正則化方法的重建誤差略小于壓縮感知方法的重建誤差，說明外推后的聲壓產(chǎn)生的誤差對傳統(tǒng)正則化方法的影響更小。

圖3分別展示了1 050 Hz時的重建面聲壓的理論云圖與利用組合近場聲全息技術(shù)前后重建面的聲壓云圖。從云圖中可以看出，測量場點較少時，利用邊界元法結(jié)合傳統(tǒng)正則化方法并不能對聲場進行有效重建，并且重建的聲壓幅值也大于理論值，而雖然利用邊界元法結(jié)合壓縮感知改善了重建結(jié)果，但在中心位置重建的幅值較小。利用組合聲全息方法，先對全息面聲壓進行外推，可以對場點聲壓進行有效地補充，可以看出重建結(jié)果與理論值吻合較好，較外推前利用傳統(tǒng)正則化方法的結(jié)果有很大的改善，相較于少場點時的壓縮感知方法，外推后的結(jié)果在中心位置的幅值也有所增加，驗證了該方法的有效性。

圖3 1 050 Hz外推前后重建面聲壓云圖比較Fig.3 Sound pressure nephogram of reconstructed plane before and after extrapolation at 1 050 Hz

2.2 外推點數(shù)的影響

這一節(jié)討論陣元間距相同，外推點數(shù)不同對重建效果的影響。仿真模型與上一節(jié)相同，外推后全息面H的尺寸分別選為0.8 m×0.8 m、0.96 m×0.96 m、1.12 m×1.12 m，對應的外推后的場點分數(shù)別為11×11、13×13、15×15。首先計算了3種情況下外推得到的聲壓與仿真得到的聲壓之間的誤差，其結(jié)果如圖4所示。

圖4 不同外推點數(shù)的估計誤差比較Fig.4 Estimation error of different extrapolation points

圖4中可以看出隨著外推點數(shù)的增加，外推的整體誤差也在增加，且在外推點數(shù)為15×15時，400 Hz以上的頻率外推誤差已經(jīng)超過40%，此時外推方法已經(jīng)失去應有的作用，這是外推方法的局限性。外推方法本質(zhì)上是通過已知的場點聲壓來對未知的場點聲壓進行估計，外推的距離越大，預測場點離已知場點越遠，聲壓的變化就越大，外推也就失去效果。而外推點數(shù)為13×13時，低頻時的外推誤差較小，而在中高頻時外推誤差則逐漸增加，這是因為隨著頻率的增加，板的振動情況逐漸復雜，外推點的聲壓與已知點的聲壓變化較大，導致外推誤差增加。

將外推后的聲壓作為輸入計算重建面的重建誤差，因為外推點數(shù)為15×15的外推誤差在中高頻時已經(jīng)超過50%，再用其進行重建也沒有意義。因此在此只比較外推點數(shù)為11×11和13×13的重建誤差，重建算法選擇壓縮感知與邊界元法的結(jié)合，重建誤差的結(jié)果如圖5所示。在200 Hz以下時外推點數(shù)為13×13的重建誤差更小，隨著頻率增加，各點聲壓的估計誤差也增加，此時引入后續(xù)重建的估計誤差的比重超過了其對場點信息的補充，導致重建誤差也隨之增加。因此對于外推方法來說，將外推點數(shù)限制在沿x，y方向分別外推兩個陣元間距時，可得到較好的重建結(jié)果，有效降低重建誤差。

圖5 不同外推點數(shù)的重建誤差比較Fig.5 Reconstruction error of different extrapolation points

2.3 外推-內(nèi)插與僅外推方法的比較

前2節(jié)主要討論了外推前后陣元間距不變時，組合近場聲全息技術(shù)的優(yōu)勢以及影響因素。本節(jié)將外推后陣元間距變?yōu)?.04 m，則全息面H上的場點數(shù)為21×21，大小為0.8 m×0.8 m，此時是在外推的基礎上同時進行了內(nèi)插，二者運用的方法相同，仿真模型布置圖如圖6所示。

圖6 外推-內(nèi)插仿真模型布置Fig.6 Diagram of extrapolation-interpolation

首先比較外推-內(nèi)插與僅外推時，全息面H上的聲壓的估計誤差，結(jié)果如圖7所示，在800 Hz以下時，外推-內(nèi)插的誤差小于只進行外推的誤差，而在頻率大于1 200 Hz時，外推-內(nèi)插的誤差卻急劇增加。這是因為在低頻時板的振動情況較為簡單，對全息面h上的數(shù)據(jù)進行外推-內(nèi)插，雖然需要估計的場點增多，但是各點估計值較為準確，而在高頻時板的振動情況變得復雜，各點估計誤差累積導致整體誤差增大。1 300 Hz時外推-內(nèi)插的誤差達到了63.8%。

圖7 外推-內(nèi)插與僅外推方法的估計誤差比較Fig.7 Estimation error comparison between extrapolation-interpolation and extrapolation methods

圖8中分別展示了1 300 Hz時只進行外推的全息面H的理論和外推聲壓云圖，以及進行外推-內(nèi)插時的全息面H的理論與外推聲壓云圖。相較于121個場點，441個場點時聲壓云圖的分辨率更高。

圖8 1 300 Hz時全息面理論與外推聲壓云圖Fig.8 Holographic plane theoretical and extrapolated sound pressure nephogram at 1 300 Hz

從81個場點經(jīng)過外推-內(nèi)插得到441個場點的聲壓，雖然云圖內(nèi)幅值較大位置沒有發(fā)生偏移，但與理論值相比略顯粗糙，因此估計誤差遠大于僅外推時的誤差。

利用拓展后的全息面聲壓分別通過邊界元法結(jié)合傳統(tǒng)正則化方法與壓縮感知方法進行聲場重建，重建的誤差如圖9所示。盡管在1 150 Hz以后，相較于僅外推時，外推-內(nèi)插對全息面聲壓的估計誤差急劇增加，但是其重建誤差卻與僅外推時的誤差基本相同。出現(xiàn)上述現(xiàn)象的原因在于與理論值相比，外推-內(nèi)插得到的全息面聲壓幅值較大的位置能夠較為準確地顯示出來，且重建面的點數(shù)為121個，其的分辨率小于全息面的分辨率。而在1 000 Hz以下的中低頻，利用外推-內(nèi)插得到的聲壓通過邊界元法結(jié)合壓縮感知方法所得的重建誤差在4種情況下最低。但外推-內(nèi)插在豐富了場點信息的同時，也為后續(xù)傳遞矩陣的計算增加了時間成本。

圖9 外推-內(nèi)插與僅外推時重建誤差的比較Fig.9 Reconstruction error comparison between extrapolation-interpolation and extrapolation methods

3 板聲場重建實驗

為驗證算法的性能，針對板聲源展開實驗研究，實驗布置如圖10所示。將0.8 m×0.8 m×0.003 m的板固定在箱體上，箱體內(nèi)懸空設置一個激振器對板施加激勵，使其振動作為目標聲源，激振器的信號頻率設置為100～600 Hz，間隔為100 Hz，分析不同頻率下算法的重建結(jié)果。全息面與重建面的大小均為0.66 m×0.66 m，位置分別距離板聲源0.04 m和0.02 m。全息面與重建面聲壓均由6×6的麥克風陣列分4次測得，點數(shù)為144個，另設參考麥克風保證測量聲壓的相位，相鄰麥克風之間的間距為0.06 m。測量結(jié)束后，提取全息面中心64個點作為外推前的場點，組合成全息面h，大小為0.42 m×0.42 m。選取外推后全息面H的大小為0.54 m×0.54 m，場點數(shù)為100個。保持全息面H的大小不變，減小陣元間距，則外推-內(nèi)插的場點數(shù)為361個。

圖10 板實驗布置Fig.10 Experimental setup of plate

圖11為600 Hz時采用組合近場聲全息技術(shù)前后的重建面的聲壓分布云圖，并將其與實驗測得的云圖進行比較。從圖中可以看出，選取64個場點時，全息面h面積較小，雖然在板中心的聲壓分布情況可以重建，但是對于邊緣處則由于缺少必要的場點信息而無法得到較好的重建效果。而采用組合近場聲全息技術(shù)以后，分別對全息面h上的場點聲壓進行了外推以及外推-內(nèi)插，從圖中可以看出，在重建面邊緣位置可以得到較好的重建結(jié)果，且相較于僅外推時，采用外推-內(nèi)插后重建結(jié)果的幅值與測量值更為接近，重建效果有所改善。

圖11 600 Hz重建面聲壓云圖Fig.11 Sound pressure nephogram of reconstructed surface at 600 Hz

圖12給出了各頻率下采用組合近場聲全息技術(shù)前后的重建誤差，可以看出除了100 Hz，其他頻率處經(jīng)過外推后重建誤差都有所降低，且隨著頻率增加，誤差降低的幅度也隨之增加。和前文分析相一致，低頻時較小的全息面即可進行聲場重建，在中高頻時波長變短，場點包含的信息變少，且板的振動情況也隨之變得復雜。此時，采用外推方法可對全息面聲壓進行有效補充，重建誤差也隨之降低。采用外推-內(nèi)插則是在外推的基礎上更加豐富了場點信息，在中低頻時可改善測量結(jié)果。

圖12 采用組合近場聲全息前后重建誤差比較Fig.12 Comparison of reconstruction errors before and after using combined near-field acoustic holography

4 結(jié)論

1) 組合近場聲全息技術(shù)針對傳統(tǒng)近場聲全息技術(shù)少場點測量重建精度不高的問題，通過SONAH進行聲場外推，將邊界元法與壓縮感知相結(jié)合，有效降低了聲場重建誤差。

2) 采用組合近場聲全息技術(shù)，分析外推點數(shù)與外推面積的選定，通常沿x，y方向分別外推2個陣元間距，即可得到較好的重建效果。

3) 相較于對全息面僅外推時，外推-內(nèi)插的重建結(jié)果在中低頻時有一定的改善，但同時增加了計算的時間。