船舶機(jī)艙內(nèi)高頻弱聲源近場(chǎng)聲全息方法

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

(海軍工程大學(xué) 動(dòng)力工程學(xué)院，湖北 武漢 430033)

0 引 言

在船舶機(jī)艙環(huán)境下，高壓流體泄漏產(chǎn)生的噪聲與機(jī)艙內(nèi)其他動(dòng)力裝置、艙體振動(dòng)等產(chǎn)生的噪聲相比，其頻率高、聲壓級(jí)低、指向性強(qiáng)、隨傳播距離衰減快，符合高頻弱聲源的特征。船員往往無(wú)法及時(shí)感知此種高頻弱聲源的出現(xiàn)，而通過(guò)近場(chǎng)聲全息（Nearfield Acoustical Holography，NAH）[1]的聲成像技術(shù)可以獲知一定空間內(nèi)聲場(chǎng)分布情況并直觀地定位出噪聲源，還可對(duì)相應(yīng)設(shè)備、系統(tǒng)運(yùn)行情況進(jìn)行監(jiān)測(cè)。

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

本文結(jié)合高頻弱聲源的特征與聲源在平面上的分布特點(diǎn)，提出在船舶機(jī)艙內(nèi)的低信噪比條件下，結(jié)合經(jīng)驗(yàn)?zāi)B(tài)分解（Empirical Mode Decomposition，EMD）法與壓縮感知的高頻弱聲源近場(chǎng)聲全息方法。通過(guò)仿真實(shí)驗(yàn)，與Fourier變換方法和BEM的重構(gòu)聲壓效果相比較，所得結(jié)果表明該方法可作為船舶機(jī)艙內(nèi)低信噪比條件下的高頻弱聲源近場(chǎng)聲全息方法的一種有效補(bǔ)充，具有較大的應(yīng)用價(jià)值。

1 理論模型及算法實(shí)現(xiàn)

1.1 平面等效源近場(chǎng)聲輻射計(jì)算模型

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

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

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

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

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

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

虛擬源強(qiáng)的反演公式為：

進(jìn)而將實(shí)際聲源面的聲壓表示為：

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

已知L和Ph，便可重構(gòu)出近場(chǎng)平面上的聲壓分布。

1.2 陣元信號(hào)權(quán)值計(jì)算

實(shí)際情況中，原信號(hào)和噪聲信號(hào)往往為非平穩(wěn)信號(hào)，能量或聲壓有效值不知，無(wú)法得到信號(hào)的信噪比。現(xiàn)對(duì)各陣元采樣所得信號(hào)運(yùn)用EMD法進(jìn)行分析，該方法能使復(fù)雜信號(hào)分解為有限個(gè)本征模函數(shù)（Intrinsic Mode Function，IMF），所得各IMF分量包含了原信號(hào)不同時(shí)間尺度的局部特征信號(hào)。

圖2 信號(hào)時(shí)頻圖Fig.2 Time-frequency diagrams of each signal

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

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

可知，該混合信號(hào)的高頻部分（即f=5 kHz）分布在IMF1中。頻域圖中，縱坐標(biāo)數(shù)值越大，代表對(duì)應(yīng)此頻率的信號(hào)部分在該分量中的貢獻(xiàn)越大。將每個(gè)IMF頻域圖上各頻率點(diǎn)對(duì)應(yīng)的幅值累加，共得3個(gè)值，再判斷每個(gè)分量主要成分對(duì)應(yīng)頻率是否大于4.9 kHz。此信號(hào)中，IMF1的主要成分對(duì)應(yīng)頻率大于4.9 kHz。將所得3個(gè)值累加作為分母，IMF1對(duì)應(yīng)的值累加作為分子，即得到該混合信號(hào)所需要乘的權(quán)值α=0.317，進(jìn)而得到此信號(hào)中高頻部分的聲壓值。α值與實(shí)際高頻部分聲壓貢獻(xiàn)值占比0.316十分接近，表明該方法可行。

1.3 聲壓陣信號(hào)的稀疏表示模型

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

圖4為傳聲器陣列信號(hào)采集模型，在全息平面Sh上，有一個(gè)按規(guī)則排布的M元傳聲器陣列。聲源點(diǎn)通常數(shù)目有限，故聲源點(diǎn)在聲源面Ss上為稀疏分布。可對(duì)聲源面Ss進(jìn)行均勻劃分格點(diǎn)，使聲源點(diǎn)位置在其上能夠稀疏表示[9]。

如圖5所示，聲源面Ss被均分為R個(gè)區(qū)域，實(shí)心區(qū)域表示該處存在聲源，空心區(qū)域表示該處不存在聲源。可將規(guī)則分布的M元傳聲器陣列接受到的聲壓信號(hào)稀疏表示為：

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

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

1.4 結(jié)合壓縮感技術(shù)的平面等效源近場(chǎng)聲全息方法實(shí)現(xiàn)

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

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

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

2 仿真分析

2.1 圓盤(pán)振動(dòng)模型

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

圖6 聲場(chǎng)三維模型Fig.6 Sound field three-dimensional model

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

2.2 聲場(chǎng)全息仿真成像分析

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

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

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

本文聲陣列采樣在時(shí)域內(nèi)進(jìn)行，選取聲陣列中某一傳聲器，其在時(shí)域內(nèi)采樣所得聲壓變化如圖10所示。可知，在0.4 s后該點(diǎn)聲壓變化趨于穩(wěn)定。故均選用0.4 s后的數(shù)據(jù)進(jìn)行運(yùn)算。

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

圖10 陣元信號(hào)時(shí)域圖Fig.10 Time-domain map of a single element signal

采用聲陣列對(duì)圖8（b）平面上聲壓進(jìn)行采樣，各傳聲器采樣所得時(shí)域信號(hào)時(shí)長(zhǎng)均為0.005 s。計(jì)算得到121個(gè)點(diǎn)的聲壓值，各采樣點(diǎn)平均聲壓幅值為0.038 Pa。

將虛擬聲源面設(shè)置在實(shí)際聲源面上，并均勻劃分為16×16個(gè)網(wǎng)格。將聲壓重構(gòu)面設(shè)在距離圓盤(pán)0.1 m的平面上，均勻劃分為21×21個(gè)網(wǎng)格。聲陣列在距圓盤(pán)0.2 m的全息面上，采樣得到的聲壓分布如圖11（a）所示。

在實(shí)際船舶機(jī)艙內(nèi)，各干擾噪聲源距離全息面相對(duì)較遠(yuǎn)，可視為遠(yuǎn)場(chǎng)聲源。當(dāng)干擾噪聲到達(dá)全息面時(shí)可近似視為平面波，故聲陣列平面上的噪聲可視為均勻分布在各采樣點(diǎn)上。此處，在每個(gè)采樣點(diǎn)的數(shù)據(jù)中增加時(shí)長(zhǎng)同為0.005 s的干擾噪聲模擬噪聲干擾。

現(xiàn)分為3組進(jìn)行實(shí)驗(yàn)，每組添加不同的噪聲信號(hào)，具體情況為：

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

2）添加時(shí)長(zhǎng)為0.005 s，表達(dá)式為

3）添加時(shí)長(zhǎng)為0.005 s，表達(dá)式為

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

當(dāng)添加噪聲信號(hào)表達(dá)式為y=0.0 1 9×?xí)r，運(yùn)用Fourier變換法、BEM和權(quán)值法進(jìn)行0.1 m平面聲壓重構(gòu)，所得結(jié)果分別如圖12所示。

圖12 聲壓重構(gòu)圖Fig.12 Sound pressure reconstruction maps

當(dāng)添加噪聲信號(hào)表達(dá)式為y=0.0 3 8×?xí)r，運(yùn)用Fourier變換法、BEM和權(quán)值法進(jìn)行0.1 m平面聲壓重構(gòu)，所得結(jié)果分別如圖13所示。

當(dāng)添加噪聲信號(hào)表達(dá)式為y=0.1 9×?xí)r，運(yùn)用Fourier變換法、BEM和權(quán)值法進(jìn)行0.1 m平面聲壓重構(gòu)，所得結(jié)果分別如圖14所示。

2.3 結(jié)果與討論

仿真實(shí)驗(yàn)中，各組采樣所得聲信號(hào)平均信噪比分別為6.02 dB，0 dB和-13.98 dB。當(dāng)信噪比為零或負(fù)值時(shí)，表明此時(shí)信號(hào)功率小于噪聲功率。

圖13 聲壓重構(gòu)圖Fig.13 Sound pressure reconstruction maps

圖14 聲壓重構(gòu)圖Fig.14 Sound pressure reconstruction maps

添加噪聲信號(hào)表達(dá)式為y=0.0 1 9×?xí)r，對(duì)圖12進(jìn)行分析。僅采用權(quán)值法所得重構(gòu)聲壓分布情況與圖8（a）一致性較好，能清晰分辨出聲源數(shù)目與分布位置。而采用Fourier變換法和BEM所得結(jié)果不夠理想，只能大概判斷聲源點(diǎn)所在區(qū)域，聲源數(shù)目不能夠清晰分辨。

添加噪聲信號(hào)表達(dá)式為y=0.0 3 8×?xí)r，對(duì)圖14進(jìn)行分析。各方法所得結(jié)果與添加噪聲信號(hào)表達(dá)式為時(shí)基本一致。

添加噪聲信號(hào)表達(dá)式為y=0.1 9×?xí)r，對(duì)圖14進(jìn)行分析。此時(shí)僅采用權(quán)值法所得的重構(gòu)聲壓分布情況與圖8（a）一致性較好，能清晰分辨出聲源數(shù)目與分布位置。而采用Fourier變換法和BEM所得結(jié)果此時(shí)受噪聲影響巨大，不再具有參考價(jià)值。

得到上述結(jié)果的原因在于，F(xiàn)ourier變換法和BEM均未對(duì)含強(qiáng)干擾噪聲的信號(hào)進(jìn)行處理，使得進(jìn)行運(yùn)算的數(shù)據(jù)包含較強(qiáng)噪聲。而權(quán)值法針對(duì)各采樣信號(hào)高頻部分，計(jì)算所需乘上的權(quán)值，故在相同條件下權(quán)值法對(duì)高頻弱聲源的定位與聲壓分布重構(gòu)效果優(yōu)于Fourier變換法和BEM。現(xiàn)對(duì)添加噪聲情況下，聲陣列上各采樣點(diǎn)是否乘上權(quán)值與理想值進(jìn)行量化對(duì)比，具體公式如下：

其中：Ph（i）為距離聲源0.2 m的全息面上各點(diǎn)實(shí)際測(cè)得聲壓值；Pf（i）為對(duì)應(yīng)平面上不添加干擾噪聲時(shí)的理論聲壓；Ph（i）和Pf（i）各有121個(gè)值。現(xiàn)采用式（14），對(duì)是否乘權(quán)值的差異進(jìn)行量化表示，計(jì)算所得結(jié)果如表1所示。

表1 采樣點(diǎn)平均差值比較Tab.1 Difference comparison

可知，乘上權(quán)值后各采樣點(diǎn)數(shù)據(jù)明顯更接近沒(méi)有干擾噪聲時(shí)的理想聲壓值。但由重構(gòu)圖12（c），13（c），14（c）也可看出，權(quán)值法重構(gòu)聲源分布接近實(shí)際情況，但其聲壓分布值與實(shí)際情況存在一定的差距。出現(xiàn)此情形原因在于EMD法本身存在的模態(tài)混疊現(xiàn)象，該現(xiàn)象會(huì)導(dǎo)致無(wú)法根據(jù)特征尺度有效地分離出不同的模態(tài)成分，出現(xiàn)現(xiàn)有的IMF包含不同時(shí)間尺度組分的情況。故通過(guò)EMD法計(jì)算得到各點(diǎn)的權(quán)值不夠精確。

3 結(jié) 語(yǔ)

為解決船舶機(jī)艙內(nèi)嘈雜環(huán)境下高壓流體泄漏不易被察覺(jué)的問(wèn)題，本文提出了一種高頻弱聲源平面近場(chǎng)聲全息方法。并利用三維聲場(chǎng)模型與Fourier變換法和BEM進(jìn)行了比較分析，得到相關(guān)結(jié)論如下：

1）高頻弱聲源平面近場(chǎng)聲全息方法，虛擬面不需要與實(shí)際聲源面共形，計(jì)算過(guò)程更方便且通用性更強(qiáng)。其網(wǎng)格劃分可以為均勻或非均勻，適用性較靈活。

2）在船舶機(jī)艙等低信噪比條件下，該方法對(duì)采樣信號(hào)進(jìn)行乘權(quán)值處理，能較好的對(duì)高頻弱聲源進(jìn)行定位與實(shí)現(xiàn)近場(chǎng)平面聲壓重構(gòu)，且重構(gòu)效果較好。

本文提出的高頻弱聲源近場(chǎng)聲全息方法，從理論到仿真實(shí)驗(yàn)結(jié)果表明該方法具有較強(qiáng)的抗噪能力，可以作為船舶機(jī)艙內(nèi)信噪比較低條件下高頻弱聲源近場(chǎng)聲全息方法的有效補(bǔ)充。