利用聲壓幅值相關性識別聲源特性研究

潘邦雄，羅瀚波，趙秀亮，趙曉丹

（江蘇大學汽車與交通工程學院，江蘇鎮江 212013）

聲源識別可以從聲學系統中得到聲源位置，能量分布，頻率特性等信息，是噪聲控制工程的重要組成部分，在機械設備故障檢測中有廣泛應用[1]。波束形成是一種目前主要使用的聲源識別技術[2]，延時求和是一種經典的波束形成算法[3-4]，它可以有效地識別平面波入射方向。隨著發展，波束形成方法逐步拓展到球面波聲源識別[5-6]、表面噪聲源識別與排序[7]、風扇空氣動力噪聲源識別等[8]，并由固定點聲源識別進化到運動點聲源識別[9-10]。為了提高波束形成的識別效果，國內外學者在延時求和算法的基礎上做出改進，提出了MUSⅠC 聲源定位[11-12]、互功率譜波束形成[13-14]和除自譜的互功率譜波束形成等算法[15-16]，以上這些算法雖然在操作形式上與延時求和有所不同，但都仍然以相位信息為核心進行聲源識別，識別過程中傳聲器陣列需要進行嚴格的相位匹配，對測量儀器要求很高。

本文將延時求和算法的輸出結果轉化成內積運算的形式，然后提出利用聲壓幅值信息通過內積運算來對聲源進行識別。以球面波為例，分析推導用聲壓幅值信息進行聲源識別的方法，繼而建立具體的聲源識別步驟，對輻射簡諧波的聲源進行識別并與波束形成法的識別特性相比較，最后結合仿真計算與消聲室實驗對該方法進行驗證。

1 延時求和波束形成算法內積形式推導

延時求和算法通過對各傳聲器接收到的聲壓信號進行相位匹配與累加求和來識別聲源入射方向，圖1 為平面波模型下的延時求和算法示意圖，根據波束形成理論[17]，該方法的輸出為：

圖1 波束形成方法示意圖

其中：B（,ω）是陣列聚焦方向為k?，聲源角頻率為ω時，波束形成方法的總輸出；M表示陣列中的傳聲器個數；Pm（ω）表示第m個傳聲器所接收到的聲壓信號；為波數；c為聲速表示陣列聚焦方向表示基準傳聲器到第m個傳聲器的距離；是向量在向量方向上的投影，即聲波到基準傳聲器與到第m個傳聲器的聲程差。

將每個傳聲器接收到的信號相對參考傳聲器進行時延，然后求和，當陣列聚焦到聲波入射方向時，各個傳聲器接收到的聲壓信號相位相同，該方向上入射的信號因為同相而增強，總輸出B（,ω）在方向上取到最大值，從而判定聲源入射角度。

對輸出結果B（,ω）做形式上的轉換，各傳聲器接收的聲壓信號Pm（ω）可以表示為：

其中：Q表示基準傳聲器接收到的聲壓信號值；表示入射方向為的聲波到基準傳聲器與到第m個傳聲器的聲程差。將各個傳聲器接收到的聲信號Pm（ω）整合成M維聲壓向量：

將聚焦方向對應的延時量整合成M維時延向量b：

輸出結果B（,ω）可以表示為：

其中：a·b表示向量a與向量b的內積。

由以上轉換可知，延時求和波束形成算法本質上是一種通過內積運算進行的相關診斷，當陣列聚焦方向與聲波入射方向一致時，聲壓向量a與時延向量b達到同一相位，兩者構成相關關系，兩個向量之間的內積取到最大值，由此判斷出聲波入射方向。

2 算法理論分析與識別步驟

前面分析，延時求和波束形成算法是通過求聲信號相位信息的最大相關性來進行聲源識別的。實際上，傳聲器陣列接收的聲壓信號中，除了相位信息以外還包括其他信息，比如聲壓幅值信息。課題組前期針對聲壓有效值進行白噪聲聲源識別進行了探究，本文進一步嘗試利用聲壓幅值信息做相關運算來識別簡諧波聲源。

球面波聲源具有典型性，工程實踐中常常將復雜聲源簡化成一個或多個球面波聲源[18]，以其為例進行分析。假設在陣列面的上方空間內存在一個球面波聲源S，示意圖見圖2，球面波聲源S的聲壓輻射公式p（t）是：

圖2 利用聲壓幅值進行聲源識別示意圖

其中：A反映聲源強度；r為接收點到聲源的距離；φ為聲源初相位。

假設使用一組M×N維的傳聲器陣列進行測量，位于陣列第i行第j列的傳聲器Mij（i=1,…,M;j=1,…,N）測得的聲壓信號幅值pij為：

其中：（xs,ys,zs）表示聲源的空間位置；（xi,yj,0 ）表示傳聲器Mij在XOY平面上的坐標。

將每個傳聲器測得的聲壓幅值按其在陣列中的位置整合成M×N維的聲壓幅值向量：

在聲源可能存在的空間區域內設置虛擬聲源S'，強度為1Pa·m，利用公式（7）計算虛擬聲源輻射到傳聲器Mij處的聲壓幅值p'ij：

其中：（x′S,y′S,z′S）表示虛擬聲源的空間位置，同樣將p′ij整合成虛擬聲壓幅值向量：

向量pv進行歸一化運算得到歸一化向量e：

將向量p與向量e做內積運算并取模值，根據柯西-施瓦茲不等式定理[19]：向量p與向量e的內積模值小于等于這兩個向量各自模的乘積：

當且僅當兩個向量滿足如下關系時：

其中：q為常數，向量e與向量p成線性相關關系，式（13）中的等號成立，內積模取到最大值由公式（14）可知，此時虛擬聲源的位置與目標聲源重合，即x′S=x′S,y′S=ys,z′S=zs。以上理論分析表明，通過搜索的最大值，使用聲壓幅值信息進行相關運算即能識別聲源。識別出聲源位置時，向量p與向量pv有以下關系：

作向量p與向量e的內積：

可得聲源強度：

基于以上理論分析進行聲源識別，建立識別步驟如下：

步驟（1）使用傳聲器陣列進行測量，得到陣列中各傳聲器處的聲壓幅值pij，按公式（9）組成聲壓幅值向量p；

步驟（2）在聲源可能存在的區域內先劃分粗略的網格，在網格的節點上設置聲源強度為1 Pa·m的虛擬球面聲源S′；

步驟（3）用公式（10）計算虛擬聲源輻射到各傳聲器的聲壓信號幅值p′ij，按公式（11）組成虛擬聲壓幅值向量pv；

步驟（4）向量pv進行歸一化運算得到向量e，計算向量p與向量e的內積模值；

步驟（5）虛擬聲源在網格上遍歷，重復步驟（3）～步驟（4）得到各網格點的內積模值；

步驟（6）搜索內積模值的最大值，找到最大值對應的網格點位置；

步驟（7）采用逐級細化在最大值位置附近做更密集的網格劃分，重復步驟（2）～步驟（6）做進一步的細化搜索；

步驟（8）判斷精度是否滿足要求，未達到精度要求時返回步驟（7），達到精度要求時停止搜索，得到聲源位置。

步驟（9）根據公式（17）計算出聲源強度。

3 簡諧波聲源識別結果對比

3.1 識別平面上的點聲源

仿真算例1：設聲源面距離測量面0.5 m，聲源理論位置為（0.451 3 m，0.322 2 m），聲源強度0.707 Pa·m，聲源頻率200 Hz。傳聲器陣列維數為6×6，傳聲器間距為0.3 m×0.3 m，要求聲源識別精度為0.001 m。分別使用延時求和波束形成方法與本文方法對該簡諧波聲源進行識別。兩種方法的識別結果如表1所示，輸出分布如圖3所示。

表1 平面上點聲源識別結果

從表1和圖3可以看出，兩種方法的輸出結果在聲源面上的分布都存在唯一一個最大值與聲源位置對應，且兩種方法識別出的聲源強度與理論值也基本符合，證明兩種方法均能準確識別出聲源位置。但是，波束形成法計算聲源強度的過程需要考慮各傳聲器接收到聲壓信號幅值的差異性，引入互譜函數算法，計算過程相對復雜[20]。而本文方法利用公式（17）進行運算可以得到準確的聲源強度，計算過程更簡便。

圖3 兩種方法的輸出分布

3.2 識別空間中的點聲源

仿真算例2：設聲源位于1 m×1 m×1 m的空間區域內（如圖2），傳聲器陣列與識別區域下表面之間的距離是0.2 m，聲源理論位置為0.535 6 m，0.468 5 m，0.735 4 m，聲源強度0.707 Pa·m，頻率200 Hz，要求聲源識別精度為0.001 m。

同時，為了考察本文方法的抗噪能力，在仿真算例2 的聲壓信號中加入高斯白噪聲信號，觀察在信噪比分別為20 dB，15 dB，10 dB下的識別效果，結果如表2所示。

表2 空間上點聲源識別結果

由表2可知，在不同強度的隨機噪聲干擾下，本文方法對聲源位置的識別結果基本穩定，當信噪比在15 dB 以上時，識別值與理論值的偏差低于5%，在工程測量允許的誤差范圍以內。這是由于本文方法使用的內積運算具有一定的濾波作用，能夠較好地減弱隨機噪聲的干擾。而延時求和等傳統波束形成算法可以識別二維平面上的聲源，但在識別聲源深度方面計算過程復雜，搜索效率低下[21—24]。本文方法對三維空間中聲源的識別只比二維平面多一維搜索，無需調整傳聲器陣列形狀，計算簡便，識別速度更快。

3.3 聲源頻率對識別效果的影響

仿真算例3：使聲源輻射200 Hz 和4 000 Hz 的簡諧波，其余參數的設定與仿真算例1 中相同。分別使用兩種算法對聲源進行識別，識別結果如表3所示，兩種方法在聲源面上的輸出分布如圖4所示。

表3 不同頻率下兩種方法的識別結果

從圖4中可以看出，識別低頻聲源時，以上兩種方法的輸出都只存在唯一的峰值與聲源位置對應；識別高頻聲源時，波束形成方法在聲源位置外也出現了輸出峰值，這是由于傳聲器間距大于聲源所輻射聲波波長而導致的“鬼影”現象，而本文方法仍然只存在唯一一個峰值與之對應。這說明本文方法識別高頻聲源時，能有效避免“鬼影”現象的產生，從而提高了聲源位置的識別精度。

圖4 不同頻率下兩種方法的輸出分布

3.4 傳聲器陣列安裝誤差對識別效果的影響

仿真算例4：假設簡諧波輻射的點聲源強度0.65 Pa·m，頻率1 000 Hz，聲源位置為（0.256 4 m，0.856 7 m）。傳聲器陣列設為6×6的矩形陣列，間隔為0.2 m×0.2 m，聲源深度為1 m。陣列中第i行第j列傳聲器的實際位置為（xij±εx,yij±εy,zij±εz），其中εx,εy,εz表示傳聲器在XYZ方向上的位置偏差，均滿足μ=0.01m,σ=0.005 m的正態分布，分別使用本文方法與波束形成方法識別，結果如表4所示。

從表4 可以看出，在傳聲器實際位置與標量之間存在一定誤差的情況下，本文方法識別出的聲源位置與強度信息誤差均小于1%，準確識別出了聲源，而波束形成方法在考慮傳聲器位置偏差時的識別誤差大于5%，這是由于使用以相位信息為核心的聲源識別方法進行識別時，傳聲器位置的偏差對相位的影響較大（尤其是在識別高頻聲源時），而本文方法使用聲壓信號有效值為信息基礎進行聲源的識別，傳聲器位置偏差帶來的聲源強度誤差較小，在一定傳聲器位置誤差下仍然能保持較好的識別效果。

表4 傳聲器位置誤差影響的結果

4 實驗驗證

在半消聲室內進行實驗，圖5 展示了實驗裝置的現場圖，球聲源用鉸鏈懸吊在空中，球聲源球心與陣列之間的距離為1 m，球聲源在聲源面的坐標是0.15 m，0.32 m。使用YG201 型傳聲器組成的4×8通道陣列進行測量，傳聲器間距為0.1 m×0.1 m?？刂芆S003 型無指向型球聲源分別產生300 Hz 和3 000 Hz 的簡諧聲波，聲源強度為0.127 3 Pa·m，使用UT3408FR3-ⅠCP 型信號采集系統采集各測量點處聲壓信號，采樣頻率為12 800 Hz，采樣時間2.58 s，分別使用波束形成方法與本文方法進行識別。觀察距離陣列1 m平面上的波束形成輸出與內積模值分布情況如圖6 所示，聲源識別實驗結果如表5所示。

圖5 實驗現場照片

從圖6 的聲學成像可看出，在聲源存在的區域兩種方法的輸出分布均存在峰值與聲源位置對應，聲源頻率的改變對內積模值的分布幾乎沒有影響，但對波束形成的輸出分布影響較大，當聲源頻率達到3 000 Hz時，在沒有聲源存在的位置也出現了“鬼影”。由表5的識別結果可看出，兩種方法對不同頻率聲源的識別結果穩定，識別誤差均小于0.05 m，在工程測量可接受的范圍之內。證明本文方法能有效識別空間中的簡諧波聲源，并有效避免了傳統波束形成方法識別高頻聲源時出現的“鬼影”現象。

圖6 不同頻率下兩種方法輸出分布

表5 不同頻率下兩種方法識別結果

5 結語

（1）本文提出一種利用聲壓幅值信息進行聲源識別的方法，該方法通過測量聲壓信號的幅值，組成聲壓幅值向量，將聲壓幅值向量與虛擬聲壓幅值向量作內積，進行相關運算來識別聲源。

（2）利用該方法對簡諧波聲源進行識別并分析其抗噪能力，將識別結果與延時求和波束形成的識別結果作比較。結果表明，該方法能準確識別出聲源位置。相比于傳統波束形成算法，該方法無需對傳聲器陣列進行嚴格的相位匹配；同時，消除了傳統波束形成算法識別高頻聲源時產生的“鬼影”現象，提高了識別精度。

（3）本文方法在傳聲器陣列存在安裝誤差時，比波束形成法擁有更準確的識別效果，具有一定的工程實用性。