船舶機艙內高頻弱聲源近場聲全息方法

陳漢濤，郭文勇，韓江桂，張宏宇

(海軍工程大學 動力工程學院，湖北 武漢 430033)

0 引 言

在船舶機艙環境下，高壓流體泄漏產生的噪聲與機艙內其他動力裝置、艙體振動等產生的噪聲相比，其頻率高、聲壓級低、指向性強、隨傳播距離衰減快，符合高頻弱聲源的特征。船員往往無法及時感知此種高頻弱聲源的出現，而通過近場聲全息（Nearfield Acoustical Holography，NAH）[1]的聲成像技術可以獲知一定空間內聲場分布情況并直觀地定位出噪聲源，還可對相應設備、系統運行情況進行監測。

Maynard等[2]利用空間二維傅里葉變換實現聲場重構，使聲全息理論走向成熟。Koopmann等[3]于1989年提出用多個等效聲源以波疊加方式來代替振動物體產生的聲場。Jeon等[4]于2005年進一步提出等效源方法（Equivalent Source Method，ESM）來重構振動聲場，并與邊界元方法（Boundary Element Method，BEM）進行了比較，達到同等聲輻射效果時等效源方法所需的測點數量更少，在工程上更易實現。Candes[5]和Donoho等[6]于2006年正式提出“壓縮感知（Compressed Sensing，CS）”的概念。Gilles Chardon[7]于2012年將壓縮感知理論首先引入平面近場聲全息中，驗證了該方法的可行性。Matteo Kirchner[8]將壓縮感知應用于圓柱體聲源的近場聲全息中，相較于傳統方法其在減少測量時使用傳聲器數目的同時，仍能得到較好的成像效果。然而，高頻弱聲源是噪聲源識別領域的難題，直接使用傳聲器陣列的傳統近場聲全息方法測量無法獲得有效的識別結果。

本文結合高頻弱聲源的特征與聲源在平面上的分布特點，提出在船舶機艙內的低信噪比條件下，結合經驗模態分解（Empirical Mode Decomposition，EMD）法與壓縮感知的高頻弱聲源近場聲全息方法。通過仿真實驗，與Fourier變換方法和BEM的重構聲壓效果相比較，所得結果表明該方法可作為船舶機艙內低信噪比條件下的高頻弱聲源近場聲全息方法的一種有效補充，具有較大的應用價值。

1 理論模型及算法實現

1.1 平面等效源近場聲輻射計算模型

將虛擬聲源點設置在一個平面上，具有較好的便捷性與通用性。如圖1所示，Sv為虛擬聲源面，Ss為聲源平面，Sf為聲場中某一平面，Sh為全息面，n為聲源面外法線方向，Ph為全息面上的聲壓分布，Ps為聲源面上聲壓分布，Qv為虛擬聲源面上聲強。

圖1 平面等效源法近場聲全息示意圖Fig.1 A schematic diagram of plane ESM.

全息面上聲壓分布可由積分方程表示如下：

式中：Ps為聲源表面聲壓列向量；為等效源序列的源強列向量；為等效源序列與聲源表面間的聲壓匹配矩陣：

同理，全息面上聲壓進行離散化處理后可表示為：

式中：Ph為全息面上聲壓列向量；Qv為等效源序列的源強列向量；Ghv為等效源序列與全息面間聲壓匹配矩陣：

虛擬源強的反演公式為：

進而將實際聲源面的聲壓表示為：

將聲源面到全息面的聲壓傳遞矩陣表示為：

已知L和Ph，便可重構出近場平面上的聲壓分布。

1.2 陣元信號權值計算

實際情況中，原信號和噪聲信號往往為非平穩信號，能量或聲壓有效值不知，無法得到信號的信噪比。現對各陣元采樣所得信號運用EMD法進行分析，該方法能使復雜信號分解為有限個本征模函數（Intrinsic Mode Function，IMF），所得各IMF分量包含了原信號不同時間尺度的局部特征信號。

圖2 信號時頻圖Fig.2 Time-frequency diagrams of each signal

該信號混有5 kHz，450 Hz和1450 Hz的3種頻率信號，若將5 kHz對應的高頻信號作為原信號，顯然該混合信號信噪比較低。故可對該信號進行預處理，乘以一個權值α，α取值為0～1之間，對應高頻部分在此信號中的占比。此過程可提高陣列分辨能力，提高對高頻弱聲源的定位精度。對此混合信號運用EMD法，得到有限個IMF分量，結果如圖3所示。

圖3 IMF分量時頻圖Fig.3 The time-frequency diagrams of IMF components

可知，該混合信號的高頻部分（即f=5 kHz）分布在IMF1中。頻域圖中，縱坐標數值越大，代表對應此頻率的信號部分在該分量中的貢獻越大。將每個IMF頻域圖上各頻率點對應的幅值累加，共得3個值，再判斷每個分量主要成分對應頻率是否大于4.9 kHz。此信號中，IMF1的主要成分對應頻率大于4.9 kHz。將所得3個值累加作為分母，IMF1對應的值累加作為分子，即得到該混合信號所需要乘的權值α=0.317，進而得到此信號中高頻部分的聲壓值。α值與實際高頻部分聲壓貢獻值占比0.316十分接近，表明該方法可行。

1.3 聲壓陣信號的稀疏表示模型

圖4 傳聲器陣列信號采樣模型Fig.4 Microphone array signal sampling model

圖4為傳聲器陣列信號采集模型，在全息平面Sh上，有一個按規則排布的M元傳聲器陣列。聲源點通常數目有限，故聲源點在聲源面Ss上為稀疏分布。可對聲源面Ss進行均勻劃分格點，使聲源點位置在其上能夠稀疏表示[9]。

如圖5所示，聲源面Ss被均分為R個區域，實心區域表示該處存在聲源，空心區域表示該處不存在聲源。可將規則分布的M元傳聲器陣列接受到的聲壓信號稀疏表示為：

其中：y為M×1維矩陣，表示全息面上傳聲器陣列接收到的信號；A為M×R維矩陣，表示傳感矩陣；x為R×1維矩陣，表示包含位置信息的聲源信號，其中只有N（N＜＜R）個非零數值；為接收到的噪聲信號。

圖5 聲源面上聲源稀疏表示Fig.5 Sound source surface sparse representation

1.4 結合壓縮感技術的平面等效源近場聲全息方法實現

實際聲源數量與分布相對于待測空間平面往往是稀疏的，選取適當的觀測矩陣，可在少測量值情況下實現聲源定位。平面等效源近場聲輻射模型中，對虛擬聲源面Sv進行了離散化處理。虛擬聲源點在虛擬聲源平面Sv上稀疏分布，將平面等效源法中的傳遞矩陣L直接作為傳感矩陣A，此矩陣符合約束等距性條件（Restricted Isometry Property，RIP）[10]。先將虛擬聲源面上的虛擬聲源點位置通過壓縮感知理論的正交匹配追蹤（Orthogonal Matching Pursuit，OMP）算法求出，再結合已知的傳遞矩陣重構出近場中平面聲壓分布。由于非均勻網格的小面元的面積計算不方便，故一般情況下取均勻網格。

虛擬聲源平面到全息面的聲壓傳遞矩陣L已由式（8）給出，此即為傳感矩陣，將其代入式（10）中得：

求出帶有稀疏信息的聲源面信息x’后，已知傳遞矩陣L，可以由式（13）求得相距為d的任意近場平面的聲壓分布圖。

2 仿真分析

2.1 圓盤振動模型

在Comsol中設置2個受到高頻點激勵的金屬圓盤，圓盤半徑均為0.05 m，厚度同為0.005 m，材料均為合金結構鋼。圓盤中軸線與y軸平行，在x-z平面上分布的圓心坐標分別為（0.2，0.2）和（-0.2，-0.2），具體分布情況如圖6所示。在金屬圓盤四周施加固定約束，在兩金屬圓盤中心同時施加N的激勵力，頻率為3 kHz。

圖6 聲場三維模型Fig.6 Sound field three-dimensional model

建立0.6 m×1 m×1 m的長方體自由聲場空間，并在四周添加厚度為0.2 m的完美匹配層。圖7為圓盤振動產生的噪聲，在時間t=0.2 s時在空間中形成的聲壓等值面分布圖。

2.2 聲場全息仿真成像分析

設置2個距離金屬圓盤分別為0.1 m和0.2 m的平面，此2個平面理想聲壓分布情況分別如圖8所示。采用11×11均勻排布的傳聲器陣列對圖8（b）所示聲壓分布進行采樣，陣列具體形式如圖9所示。

圖7 聲壓等值面分布圖（t=0.2 s）Fig.7 Sound pressure isosurface distri-bution（t=0.2 s）

圖8 平面聲壓分布圖Fig.8 Sound pressure distribution maps

本文聲陣列采樣在時域內進行，選取聲陣列中某一傳聲器，其在時域內采樣所得聲壓變化如圖10所示。可知，在0.4 s后該點聲壓變化趨于穩定。故均選用0.4 s后的數據進行運算。

圖9 傳聲器排布示意圖Fig.9 Microphones arrangement diagram

圖10 陣元信號時域圖Fig.10 Time-domain map of a single element signal

采用聲陣列對圖8（b）平面上聲壓進行采樣，各傳聲器采樣所得時域信號時長均為0.005 s。計算得到121個點的聲壓值，各采樣點平均聲壓幅值為0.038 Pa。

將虛擬聲源面設置在實際聲源面上，并均勻劃分為16×16個網格。將聲壓重構面設在距離圓盤0.1 m的平面上，均勻劃分為21×21個網格。聲陣列在距圓盤0.2 m的全息面上，采樣得到的聲壓分布如圖11（a）所示。

在實際船舶機艙內，各干擾噪聲源距離全息面相對較遠，可視為遠場聲源。當干擾噪聲到達全息面時可近似視為平面波，故聲陣列平面上的噪聲可視為均勻分布在各采樣點上。此處，在每個采樣點的數據中增加時長同為0.005 s的干擾噪聲模擬噪聲干擾。

現分為3組進行實驗，每組添加不同的噪聲信號，具體情況為：

圖11 不同噪聲條件下全息面聲壓采樣圖Fig.11 Sound pressure sampling data with different noise in the holographic plane

2）添加時長為0.005 s，表達式為

3）添加時長為0.005 s，表達式為

在距離圓盤0.2 m的全息面上，分別得到如圖11（b）～圖11（d）所示的采樣聲壓分布圖。直觀比較各圖差異不大，主要區別在于聲壓值范圍不同。為方便表示，后文將本文所提方法簡稱為權值法。

當添加噪聲信號表達式為y=0.0 1 9×時，運用Fourier變換法、BEM和權值法進行0.1 m平面聲壓重構，所得結果分別如圖12所示。

圖12 聲壓重構圖Fig.12 Sound pressure reconstruction maps

當添加噪聲信號表達式為y=0.0 3 8×時，運用Fourier變換法、BEM和權值法進行0.1 m平面聲壓重構，所得結果分別如圖13所示。

當添加噪聲信號表達式為y=0.1 9×時，運用Fourier變換法、BEM和權值法進行0.1 m平面聲壓重構，所得結果分別如圖14所示。

2.3 結果與討論

仿真實驗中，各組采樣所得聲信號平均信噪比分別為6.02 dB，0 dB和-13.98 dB。當信噪比為零或負值時，表明此時信號功率小于噪聲功率。

圖13 聲壓重構圖Fig.13 Sound pressure reconstruction maps

圖14 聲壓重構圖Fig.14 Sound pressure reconstruction maps

添加噪聲信號表達式為y=0.0 1 9×時，對圖12進行分析。僅采用權值法所得重構聲壓分布情況與圖8（a）一致性較好，能清晰分辨出聲源數目與分布位置。而采用Fourier變換法和BEM所得結果不夠理想，只能大概判斷聲源點所在區域，聲源數目不能夠清晰分辨。

添加噪聲信號表達式為y=0.0 3 8×時，對圖14進行分析。各方法所得結果與添加噪聲信號表達式為時基本一致。

添加噪聲信號表達式為y=0.1 9×時，對圖14進行分析。此時僅采用權值法所得的重構聲壓分布情況與圖8（a）一致性較好，能清晰分辨出聲源數目與分布位置。而采用Fourier變換法和BEM所得結果此時受噪聲影響巨大，不再具有參考價值。

得到上述結果的原因在于，Fourier變換法和BEM均未對含強干擾噪聲的信號進行處理，使得進行運算的數據包含較強噪聲。而權值法針對各采樣信號高頻部分，計算所需乘上的權值，故在相同條件下權值法對高頻弱聲源的定位與聲壓分布重構效果優于Fourier變換法和BEM。現對添加噪聲情況下，聲陣列上各采樣點是否乘上權值與理想值進行量化對比，具體公式如下：

其中：Ph（i）為距離聲源0.2 m的全息面上各點實際測得聲壓值；Pf（i）為對應平面上不添加干擾噪聲時的理論聲壓；Ph（i）和Pf（i）各有121個值。現采用式（14），對是否乘權值的差異進行量化表示，計算所得結果如表1所示。

表1 采樣點平均差值比較Tab.1 Difference comparison

可知，乘上權值后各采樣點數據明顯更接近沒有干擾噪聲時的理想聲壓值。但由重構圖12（c），13（c），14（c）也可看出，權值法重構聲源分布接近實際情況，但其聲壓分布值與實際情況存在一定的差距。出現此情形原因在于EMD法本身存在的模態混疊現象，該現象會導致無法根據特征尺度有效地分離出不同的模態成分，出現現有的IMF包含不同時間尺度組分的情況。故通過EMD法計算得到各點的權值不夠精確。

3 結 語

為解決船舶機艙內嘈雜環境下高壓流體泄漏不易被察覺的問題，本文提出了一種高頻弱聲源平面近場聲全息方法。并利用三維聲場模型與Fourier變換法和BEM進行了比較分析，得到相關結論如下：

1）高頻弱聲源平面近場聲全息方法，虛擬面不需要與實際聲源面共形，計算過程更方便且通用性更強。其網格劃分可以為均勻或非均勻，適用性較靈活。

2）在船舶機艙等低信噪比條件下，該方法對采樣信號進行乘權值處理，能較好的對高頻弱聲源進行定位與實現近場平面聲壓重構，且重構效果較好。

本文提出的高頻弱聲源近場聲全息方法，從理論到仿真實驗結果表明該方法具有較強的抗噪能力，可以作為船舶機艙內信噪比較低條件下高頻弱聲源近場聲全息方法的有效補充。