混響背景下基于波束成形的聲成像技術

龐凱元，劉桂峰，陳思余，夏 菁

(海軍工程大學 動力工程學院, 湖北 武漢 430000)

0 引 言

混響指聲波在傳播過程中遇到墻壁等障礙物的反射和吸收作用所形成的物理現象。波束成形技術是一種基于傳聲器陣列測量的噪聲源識別技術[1]，可以對采集到的聲場信息進行聲源面的聚焦重構，通過加強目標聲源的輸出能量實現聲源定位[2-3]。當在機艙等密閉的環境中進行聲成像時，傳感器除了接收到聲源產生的聲波信號外，還會接收到混響場中聲波經過多次反射后的反射聲波，并且隨著反射的次數增加，接收到的信號能量也逐漸衰減，最終影響成像效果。

因此，對混響背景下的聲學成像技術進行研究具有重要意義。

Bees[4]提出將去混響問題分兩步實現，最小均方準則設計去混響逆濾波器。Sehr 等[5]提出了一種基于語音特征參數域去混響的方法。Jiang[6]針對混響背景中聲源定位不準的問題，提出了一種將深度融合和卷積神經網絡相結合的新型算法。王小軍[7]提出了一種基于混響時間感知的特征提取方法。孟維鑫等[8]提出了一種新的多通道聯合去混響波束形成器。張宏宇[9]針對船艙中的混響環境，提出用房間脈沖響應函數代替自由空間格林函數，有效地抑制了混響的干擾。趙中華等[10]為了提高混響條件下時延算法的精度，提出了基于自適應特征值傳播的時間延時算法。魏元澤[11]為了有效濾除噪聲和混響抑制，提出了一種基于數學形態學和EMD 法的聲源定位方法。韓哲等[12]提出了一種新的分布式自適應去混響算法。

本文主要針對常規波束成形算法（CBF）、基于正交匹配追蹤的反卷積波束成形算法（OMP-DAMAS）以及基于脈沖響應函數的自適應波束成形算法（RIRSC-DAMAS）3 種不同的聲成像技術，模擬混響環境，通過對比定位精度、算法運行速度等，確定混響背景中進行聲源定位的最佳方法。

1 波束成形

波束成形是一種空間聲場的可視化技術，在聲源識別領域得到了廣泛的研究和應用[13]。波束成形的核心思想是將聲源面離散化形成聚焦網格點，利用傳聲器陣列采集聲信號反向聚焦各網格點，并按聚焦點位置進行相位補償后加和輸出，真實聲源所在聚焦點的輸出量被加強，其他聚焦點的輸出量被衰減，從而有效識別聲源。

1.1 常規波束成形

常規波束成形算法是利用多個陣元接收到的信號具有的相位差，通過在接收端進行延時求和、加權輸出得到輸出功率，確認波達方向。

在進行信號處理時，通常將信號從時域轉到頻域上進行處理，從而提高運行效率。陣列接收到的信號表達式為：

式中：X(t)為在t 時刻的陣列接收信號，A為聲源信號與陣列間的傳遞矩陣，S(t)和N(t)分別為聲源信號和噪聲信號。

對式（1）進行加權后得到輸出信號為：

通過傅里葉變換，將輸出信號轉換至頻域：

式中，ω=[ω1,ω2...ωM]T為權矢量，Y(k)和X(k)分別為Y(t)和X(t)在頻域上的表達式。

分別計算時域、頻域的陣列輸出平均功率：

式中：L為快拍數即采樣點數：Rt和Rk分別為信號在時域和頻域上的協方差矩陣，進行最大似然估計得到陣列接收信號的協方差矩陣為：

最終得到聲源面的聲壓值為：

1.2 反卷積波束成形

反卷積波束成形算法主要以DAMAS[14-15]算法為代表，利用點擴散函數（PSF）以及高斯賽德爾迭代法逆向計算目標聲源面的聲壓分布，從而識別聲源。該算法較常規波束成形算法提高了空間分辨率以及定位精度。

首先將目標聲源面劃分為若干網格聚焦點，測量面上測得的聲壓可表示為：

式中：G為目標聲源面與測量面之間的傳遞矩陣；q為聲源源強。

聲壓互譜矩陣的表達式為：

當目標聲源為不相干聲源時，qqH可簡化表示為。

各聚焦網格點的輸出功率為：

式中：ω為導向向量；gn為傳遞矩陣G對應的列向量。

1.3 基于脈沖響應函數的自適應波束成形

在混響背景下，聲波在傳播時受到障礙物的反射和吸收，導致陣列接收到的信息復雜多樣。而自適應波束成形算法可以依據自適應壓縮采樣重構算法(SAMP-COSAMP)，能夠自動搜索聲源個數，將測得的聲源信息進行加權處理，達到重構聲源的目的。該算法較反卷積波束成形克服了稀疏度的限制，無需提前預知聲源個數，但缺點是運行時間有所延長。

1.3.1 基于鏡像源法的脈沖響應函數

聲波在傳播過程中往往會受到壁面的反射和吸收作用，這就形成了一個復雜的混響場，當利用陣列接收傳播信號時也會因為反射聲波的干擾對最終目標聲源的定位效果產生較大的影響，而為了在混響背景下仍然能夠精準定位聲源，首先要計算聲波在混響場中的傳遞函數，即室內脈沖響應函數。

鏡像源法是利用鏡像的原理，計算聲波的傳播路程，再結合壁面對聲波的吸收和反射作用，計算再傳播過程中的能量衰減，得到室內脈沖響應函數：

通過傅里葉變換，將其轉換至頻域：

式中：βx,1，βy,1和βz,1為靠近坐標原點的壁面反射系數；βx,2，βy,2和βz,2表示遠離坐標原點的壁面反射系數；Rp為實際聲源或鏡像聲源到各陣元的距離；Rr表示多階反射對應的虛擬空間尺寸。

根據聲波在混響場中的傳遞函數h(t)，計算陣列接收到的信號：

通過傅里葉變換，將輸出信號轉換至頻域：

1.3.2 自適應壓縮采樣重構算法

陣元接收到信號之后，利用自適應壓縮采樣重構算法對信號進行重構，該算法基本步驟為：

步驟1令初始稀疏度K0=1，估計步長L=s，支撐集F=φ；

步驟2根據得到相關系數，并提取出K0個最大值對應的索引值存入索引集F；

步驟3若則K=K0+L，轉步驟2；

步驟4計算初始余量r=y-；

步驟5令初始估計信號x=0，初始化階段stage=1，迭代次數k=1，候選集C=φ，索引值集合S=φ;

步驟6根據步驟2計算相關系數，并按uj≥0.5max|u|的標準將索引值存入S中；

步驟 7令C=F∪S，根據j=1,2,...N}計算索引集C中原子與余量的相關系數，并提取出K0個最大值對應的索引值存入Fnew中，根據得到估計信號xnew，再根據rnew=y-ΦFx式得到余量；

步驟8當‖xnew-x‖2≤ε，則停止迭代，否則轉步驟9；

步驟9當‖rnew‖2≥‖r‖2，則stage=stage+1，L=stage×s，轉步驟6；否則，F=Fnew，r=rnew，k=k+1，轉步驟6。

1.3.3 混響場中基于脈沖響應函數的自適應波束成形算法

在混響場中，陣列接收信號也滿足式（1），其中的傳遞函數A即為室內脈沖響應函數，結合自適應壓縮采樣重構算法，便形成了可在混響場中進行聲學成像的波束成形算法，該算法主要流程如圖1 所示。

圖1 算法流程圖Fig.1 Flow diagram of the algorithm

2 仿真分析

常見機械故障特征頻率分布在1 800～4200 Hz 頻帶內，故選取1000 Hz，2000 Hz，4000 Hz 分別為低、中、高頻進行仿真研究分析。利用Comsol Multiphysics 對聲場信息進行建模和仿真實驗，其主要思路為：選擇射線聲學模塊進行建模仿真，首先構建出混響場模型，選擇射線從點釋放作為聲源，然后在距離聲源面一定位置處添加測量陣列，用于捕捉聲場信息，最后用上文3 種波束成形算法對捕捉到的聲壓數據進行處理，完成聲源定位工作，并與真實的聲源位置進行對比，確定混響背景下最佳的聲成像方法。

2.1 混響場模型

首先在環境中構建一個3 m×4 m×3 m 的長方體混響區域，將整個區域的介質設置為相對濕度為50%的空氣，并將流體模型設置為大氣衰減。然后在（1.5 m,3 m, 1.5 m）的位置構建一點聲源，設置總源功率為1 W，波矢空間射線數為100 000。將陣列測量面設置在距離聲源面1 m 處，陣元數為11×11，陣元間距為0.1 m。構建的聲場模型及陣列如圖2 和圖3 所示。由于是混響環境，需要對各個壁面的吸聲系數進行設置，整個長方體混響場的正面及左側面材料為玻璃，右側面及頂部為墻壁，地面材料為瓷磚，后側面為上窗下柜，不同壁面材料在不同頻率下的吸聲系數如表1 所示。

圖2 混響場模型Fig.2 Reverb-field model

圖3 陣列位置示意圖Fig.3 Schematic representation of the array position

表1 各壁面吸聲系數Tab.1 Sound absorption coefficient of each wall surface

根據表1 對各壁面進行設置，在進行Comsol 研究時，以0.01 s 為步長，停止時間為0.1 s 進行參數化掃描，得到各頻率下的聲源射線軌跡點樣式以及測量面的聲壓級分布。圖4 和圖5 分別為頻率為1000 Hz 時仿真得到的射線軌跡和聲壓級分布圖。

圖4 射線軌跡圖Fig.4 Ray track diagram

圖5 測量面聲壓級分布圖Fig.5 Sound pressure level distribution diagram of the measurement surface

2.2 仿真實驗與結果分析

利用Comsol 仿真軟件導出頻率分別為1000 Hz，2000 Hz 以及4000 Hz 下的聲壓數據，將其分別利用常規波束成形算法（CBF）、基于OMP 的反卷積波束成形算法（OMP-DAMAS）以及基于壓縮感知的自適應波束成形算法（RIR-SC-DAMAS）進行處理，得到相應的聲成像如圖6～圖8 所示。3 種算法在不同頻率下的聲源定位結果對比如表2所示。

圖6 1000 Hz 下各算法聲學成像圖Fig.6 Acoustic imaging diagram of each algorithm at 1000 Hz

圖8 4000 Hz 下各算法聲學成像圖Fig.8 Acoustic imaging diagram of each algorithm at 4000 Hz

圖7 2000 Hz 下各算法聲學成像圖Fig.7 Acoustic imaging diagram of each algorithm at 2000 Hz

通過圖6～圖8 及表2 可知，在混響場中相同頻率下CBF 算法的定位結果旁瓣大，成像多，定位效果較差；OMP-DAMAS 算法旁瓣小，但是只在中低頻有較好的定位效果，高頻下會影響定位精度；RIR-SC-DAMAS 算法無論是在中低頻還是高頻下都有良好的定位效果。

表2 不同頻率下成像結果對比表Tab.2 Comparison table of the imaging results at different frequencies

各算法分別在不同頻率下運行5 次，5 次運行時間的平均值如表3 所示。

表3 各算法運行時間表Tab.3 Operation schedule of each algorithm

可以看出，每種算法在不同頻率下的運算時間相差不大，但不同算法在同一頻率下的運算時間有很大的差別。其中CBF 算法在運行時僅需普通的數據計算，所以其運行時間最短；OMP-DAMAS 算法較CBF 算法需要利用反卷積以及壓縮感知對聲源進行重構，大大增加了運行時間；RIR-SC-DAMAS 算法在前兩者的基礎上還要模擬混響環境，重新構造室內脈沖響應函數，故運行時間較前兩者有明顯的增加。

為進一步比較3 種算法在復雜環境下的成像精度，在混響場中加入snr 為0 的高斯白噪聲，分別利用CBF，OMP-DAMAS 以及RIR-SC-DAMAS 算法對陣列接收的聲場信息進行處理，得到不同頻率下對應的聲成像如圖9～圖11 所示。

圖9 1000 Hz 下各算法聲學成像圖(SNR=0)Fig.9 Acoustic imaging diagram of each algorithm at 1 000 Hz (SNR=0)

圖10 2000 Hz 下各算法聲學成像圖(SNR=0)Fig.10 Acoustic imaging diagram of each algorithm at 2000 Hz (SNR=0)

在加入SNR=0 的高斯白噪聲條件下各算法在不同頻率下的聲源定位結果對比如表4 所示。

通過圖9～圖11 及表4 可知，在加入干擾噪聲后，CBF 算法的成像精度大大降低，并且在聲成像圖上出現大量干擾聲源；OMP-DAMAS 算法在中低頻下受影響不大，但在高頻下的聲成像效果遠遠偏離了實際聲源位置；RIR-SC-DAMAS 算法在加入干擾噪聲后受到的影響較小，仍能保持較高的成像精度。

表4 不同頻率下成像結果對比表(SNR=0)Tab.4 Comparison table of imaging results at different frequencies (SNR=0)

圖11 4000 Hz 下各算法聲學成像圖(SNR=0)Fig.11 Acoustic imaging of each algorithm at 4000 Hz (SNR=0)

計算各算法在加入干擾噪聲后不同頻率下的運行時間，取5 次運行時間的平均值如表5 所示。

表5 與表3 所得結論基本相同，在同一頻率下，CBF 算法用時最短，OMP-DAMAS 算法次之，RIR-SCDAMAS 算法用時最長。

表5 各算法運行時間表(SNR=0)Tab.5 Operation schedule of each algorithm (SNR=0)

綜上所述，在混響場中，CBF 算法以及OMP-DAMAS 算法只在中低頻下具有一定的定位精度，而RIRSC-DAMAS 算法較前兩者的定位精度有極大的精準度，特別是當干擾噪聲存在時，該算法受到的影響較小，而另外2 種算法基本無法準確定位.但從運行時間來看，RIR-SC-DAMAS 算法的運行時間遠遠超過其他2 種算法。

3 結 語

本文首先介紹了3 種基于波束成形算法的聲學成像技術，即常規波束成形法（CBF）、基于正交匹配追蹤的反卷積波束成形算法（OMP-DAMAS）以及基于脈沖響應函數的自適應波束成形算法（RIR-SC-DAMAS），然后構建混響場模型，分別利用3 種算法對不同頻率下的目標聲源進行聲學成像，最后對仿真結果進行對比分析并得出結論，即在有干擾噪聲存在的混響背景中，基于脈沖響應函數的自適應波束成形算法具有更好的定位精度和效果。