侯 森，胡長青

(1.中國科學(xué)院聲學(xué)研究所東海研究站，上海 201815；2.中國科學(xué)院大學(xué)，北京 100049)

0 引 言

氣體泄漏廣泛存在于海洋環(huán)境之中，其中海底甲烷水合物自然泄漏及相關(guān)開采運(yùn)輸造成的泄漏，往往都具有長時(shí)間持續(xù)泄漏的特征，不僅會(huì)造成經(jīng)濟(jì)損失，更容易引發(fā)安全和環(huán)境問題[1-3]。依據(jù)泄漏氣泡產(chǎn)生的聲學(xué)信號作為檢測識別方法，具有易于實(shí)現(xiàn)、適用性強(qiáng)等特點(diǎn)。

為了實(shí)現(xiàn)在實(shí)際海洋環(huán)境中對氣體泄漏的聲學(xué)檢測，很多人對氣泡聲信號進(jìn)行了研究。王鑫等[4]對氣泡聲發(fā)射信號進(jìn)行了分類，并指出當(dāng)泄漏量較小時(shí)，氣泡聲發(fā)射信號以振蕩聲為主。侯森等[5]基于水合物氣泡泄漏時(shí)的物理特性，對泄漏氣泡振蕩特性進(jìn)行了修正，計(jì)算了不同環(huán)境下泄漏氣泡半徑與能量和半徑與頻率的關(guān)系。杜非[6]指出海洋中持續(xù)性氣泡泄漏具有較強(qiáng)的周期性和重復(fù)性，并研究了微量氣泡泄漏時(shí)的聲學(xué)特性。吳連軍等[7]研究了微弱信號的提取方法，指出頻域積累法可以有效提高窄帶信號的信噪比。此外，還有很多人對微弱聲信號的提取和檢測做出了研究[8-10]。

為了解決氣泡信號容易被環(huán)境噪聲遮蔽的問題，本文從氣泡振蕩特性的研究出發(fā)，采用了頻域積累的方法進(jìn)行降噪處理，并針對連續(xù)泄漏氣泡具有強(qiáng)周期性的特點(diǎn)，通過一種基于時(shí)頻分析的自適應(yīng)檢測算法，判別有無氣泡聲信號并對氣泡振蕩主頻率進(jìn)行檢測。通過仿真和實(shí)驗(yàn)，驗(yàn)證了本文方法可以在較強(qiáng)干擾的環(huán)境中實(shí)現(xiàn)對持續(xù)泄漏氣泡的檢測、判別。

1 泄漏氣泡聲學(xué)特性

1.1 氣泡動(dòng)力學(xué)模型

當(dāng)泄漏量較小時(shí)，持續(xù)泄漏氣泡的聲信號主要由氣泡振蕩產(chǎn)生[4]。在構(gòu)建氣泡振蕩模型時(shí)，做了以下假設(shè)：氣泡為理想球形、液體不可壓縮、忽略海流等因素引起的壓力場變化并忽略泡內(nèi)蒸氣壓。同時(shí)在模型中考慮了液體黏性以及氣泡內(nèi)氣體非理想時(shí)對氣泡振蕩的影響。氣泡內(nèi)壓強(qiáng)為

其中：P0為環(huán)境壓強(qiáng)，R0為氣泡平衡半徑；Pi0為氣泡平衡時(shí)的內(nèi)部壓強(qiáng)；σ為張力系數(shù)?？紤]到實(shí)際氣泡內(nèi)為非理想氣體，引入描述氣泡可壓縮性的系數(shù)β進(jìn)行修正[5]，則任意時(shí)刻下的氣體狀態(tài)方程如式(2)所示：

其中：Pi為任意時(shí)刻的氣體內(nèi)部壓強(qiáng)，R為任意時(shí)刻下的氣泡半徑大小，κ為多方指數(shù)。引入粘滯系數(shù)u，則任意時(shí)刻下氣泡外部壓強(qiáng)Pout為

假設(shè)液體不可壓縮，根據(jù)能量守恒定理，液體動(dòng)能的增量等于氣泡對液體所做的功，如式(4)所示：

由于式(5)難以求取解析解，本文通過四階龍格庫塔(Runge-kutta)算法進(jìn)行了特定初值下的數(shù)值求解。設(shè)環(huán)境參數(shù)如下：P0=1.01×105Pa；β=0.1；ρ=1 000 kg·m-3；κ=1.4 ；σ=72×10-3N·m-1；u= 8×10-4Pa·s-1。氣泡初始半徑為平衡半徑的 1.2倍，即Ri/R0=1.2。則不同半徑氣泡的振動(dòng)如圖 1所示，氣泡在泄漏時(shí)獲得初始能量，并圍繞平衡半徑做振蕩衰減。從圖 1中可以看出，氣泡振動(dòng)聲信號單頻特性明顯，且能量衰減很快，環(huán)境一定時(shí)，氣泡振蕩頻率隨氣泡半徑增大而降低；氣泡振幅則隨氣泡半徑增大而增大。

圖1 不同半徑氣泡振蕩曲線Fig.1 Vibrations of bubbles with different radii

海洋環(huán)境中發(fā)生的泄漏通常具有持續(xù)性，并且有位置固定，壓差、漏孔等參數(shù)在一定時(shí)間內(nèi)穩(wěn)定的特點(diǎn)[6]。故同一泄漏點(diǎn)產(chǎn)生的氣泡往往具有氣泡半徑相近、泄漏氣泡間隔相對恒定的特點(diǎn)。為模擬海洋環(huán)境中持續(xù)泄漏氣泡的聲信號，除采用上文的環(huán)境參數(shù)外，設(shè)生成氣泡的平衡半徑R0為1.3～1.8 mm間的隨機(jī)值，氣泡生成間隔為0.1 s。圖 2所示是理想無噪聲的連續(xù)氣泡聲信號，由于氣泡能量衰減極快，在時(shí)域上表現(xiàn)為振動(dòng)頻率相近的一連串脈沖信號。而當(dāng)環(huán)境噪聲能量較強(qiáng)時(shí)，氣泡聲信號易被淹沒。在加入信噪比為?10 dB的加性高斯噪聲后，噪聲背景下的氣泡聲信號如圖 3所示，氣泡聲信號被環(huán)境噪聲淹沒。在強(qiáng)噪聲干擾下，很難在時(shí)域上觀察出氣泡振蕩引起的脈沖波形。

圖2 無噪聲干擾時(shí)持續(xù)泄漏氣泡聲信號Fig.2 Bubble sound signal without noise

圖3 噪聲背景下氣泡聲信號Fig.3 Bubble sound signal under noise background

1.2 頻域積累

針對氣泡聲信號易被環(huán)境噪聲淹沒的問題，由于氣泡聲和環(huán)境噪聲頻段重疊，傳統(tǒng)的濾波方法難以有效濾除噪聲干擾。當(dāng)氣泡持續(xù)泄漏時(shí)，泄漏產(chǎn)生的聲信號主要集中在一定頻段范圍內(nèi)，采用頻域積累的方法，可以使氣泡聲能量得到明顯加強(qiáng)；而環(huán)境噪聲的相關(guān)性差，經(jīng)過頻域累加后，其能量不容易被積累，故可以有效提高信噪比。頻域積累法的主要步驟為：將待處理的信號x(n)分割成M段等長度的信號片段，依次將各信號片段做傅里葉變換并在頻域進(jìn)行累加，基于離散傅里葉變換的的頻域累加公式為[12]

其中：xm(l)為第m段信號，L為離散傅里葉變換的點(diǎn)數(shù)，表示旋轉(zhuǎn)因子，X(L)為傅里葉變換后的頻域信號。將經(jīng)過頻域積累后的X(L)進(jìn)行反傅里葉變換即可得重構(gòu)的時(shí)域信號波形。對圖3進(jìn)行頻域積累降噪處理，重構(gòu)后的時(shí)域信號如圖4所示。

由圖4和圖2的對比可以看出，經(jīng)過頻域積累處理，重構(gòu)后的聲信號可以較好地還原出原信號的時(shí)域波形，脈沖峰的間隔及振蕩幅值等信息得到了較好保留。由圖4和圖3的對比可以看出，圖4中信號的信噪比有明顯的提升，且降噪效果隨頻域積累次數(shù)的增加而提升。

圖4 基于頻域積累法降噪處理Fig.4 Noise reduction based on frequency domain accumulation

1.3 時(shí)頻分析

對頻域積累前后的聲信號做時(shí)頻分析，如圖 5所示。圖5(a)為圖3中信號的時(shí)頻圖，從圖中可以看出，當(dāng)加入寬帶強(qiáng)噪聲干擾后，在時(shí)頻圖中難以分辨氣泡聲信號。圖5(b)為圖4(c)中信號的時(shí)頻分析，從圖中可以看出，進(jìn)行 50次頻域積累處理后，信噪比明顯提高，同時(shí)氣泡的振蕩頻率以及泄漏時(shí)間間隔等信息保留較好。

圖5 頻域積累前后時(shí)頻分析Fig.5 Time-frequency analyses before and after frequency domain accumulation

從圖5(b)中可以看出，氣泡聲信號疊加噪聲的模型中，能量隨時(shí)間的變化在不同的頻段處有明顯的差異，為此對圖5(b)進(jìn)行了展開討論，圖6和圖7分析了不同頻段處的能量變化規(guī)律。圖6(a)為氣泡振蕩主頻段處的時(shí)間-能量變化圖，由于持續(xù)泄漏氣泡具有很強(qiáng)的重復(fù)性，所以在能量變化上表現(xiàn)出了強(qiáng)周期性。為了描述能量的周期性變化，對時(shí)間-能量曲線做了頻譜分析，結(jié)果如圖6(b)所示，時(shí)間-能量曲線的頻譜具有明顯的基頻峰和諧波峰，基頻和氣泡泄漏時(shí)間間隔吻合。而非氣泡振蕩頻段的能量變化如圖 7所示，可以發(fā)現(xiàn)能量隨時(shí)間的變化具有隨機(jī)性，時(shí)間-能量的頻譜中也沒有基頻峰和諧波峰的對應(yīng)關(guān)系。所以檢測特定頻段下的能量變化規(guī)律可以作為氣泡檢測的依據(jù)，具體方法將在第2節(jié)中展開討論。

圖6 氣泡聲信號主頻段處聲信號能量隨時(shí)間和頻率的變化Fig.6 Variation of sound signal energy at bubble resonance frequencies with time and frequency

圖7 低頻段聲信號能量隨時(shí)間和頻率的變化Fig.7 Variation of sound signal energy at low frequencies with time and frequency

2 實(shí)驗(yàn)數(shù)據(jù)處理

2020年項(xiàng)目組在某海域進(jìn)行了海上實(shí)驗(yàn)。實(shí)驗(yàn)設(shè)計(jì)圖如圖8(a)所示，利用高壓氣泵向單孔導(dǎo)氣管供氣來模擬海洋中持續(xù)泄漏氣泡，在出氣口下方 1 m處布置水聽器和溫深儀(Temperature Depth,TD)，并選取了 4個(gè)深度測量氣泡聲信號數(shù)據(jù)。實(shí)驗(yàn)深度變化如圖8(b)所示，分別在5.5、16.0、24.0和 36.5 m深度處進(jìn)行了氣泡聲信號測量。自容式水聽器和溫深儀如圖8(c)示。實(shí)驗(yàn)平臺及供氣氣泵如圖8(d)所示，實(shí)驗(yàn)過程中通過改變供氣壓強(qiáng)，保持實(shí)驗(yàn)過程壓差基本恒定在0.45 MPa。

圖8 實(shí)驗(yàn)場景簡介Fig.8 Introduction to the experimental scene

以5.5 m處的聲信號數(shù)據(jù)處理為例，經(jīng)過濾波等預(yù)處理后，接收到的信號如圖 9所示。由圖 9可見，在實(shí)際海洋環(huán)境下，氣泡信號被環(huán)境噪聲完全淹沒，很難觀察到氣泡信號特征。運(yùn)用頻域積累的方法進(jìn)行 30次處理后，結(jié)果如圖 10所示，從時(shí)域上看，信噪比有所提高，可以在部分時(shí)段見到氣泡形成的脈沖峰，但噪聲干擾依然嚴(yán)重。從圖 11中的時(shí)頻圖上看，除了一些可能由氣泡信號引起的亮點(diǎn)外，低頻段的噪聲干擾明顯，較難區(qū)分出氣泡信號的能量頻段。由此可見，在非理想噪聲環(huán)境下，頻域積累方法結(jié)合傳統(tǒng)的時(shí)頻分析不能有效檢測氣泡聲。

圖9 預(yù)處理后海洋環(huán)境中聲信號圖Fig.9 Sound signal in marine environment after preprocessing

圖10 頻域積累后氣泡聲信號圖Fig.10 The bubble sound signal after frequency domain accumulation

圖11 頻域積累后氣泡時(shí)頻圖Fig.11 Time-frequency analysis of the bubble sound signal after frequency domain accumulation

為了在較強(qiáng)干擾下對氣泡振蕩頻段進(jìn)行檢測，本文基于1.3節(jié)中的內(nèi)容采用一種自適應(yīng)檢測算法對氣泡聲信號進(jìn)行搜索判別。為了方便區(qū)分，下述算法中F表示聲信號的頻率，f表示為時(shí)間-能量頻譜圖中的頻率。

(1) 假設(shè)Fi為待檢測的聲信號頻率，計(jì)算聲信號在Fi處的時(shí)間-能量變化關(guān)系并做頻譜分析。對時(shí)間-能量曲線的頻譜做歸一化處理，并設(shè)[fmin,fmax]為搜索區(qū)間，進(jìn)行峰值檢測。

(2) 從低頻峰開始遍歷搜索區(qū)間內(nèi)的峰值，將第一個(gè)峰值f0設(shè)為基頻，并進(jìn)行諧波檢測。搜索[nf0??f,nf0+?f]區(qū)間內(nèi)是否存在對應(yīng)的第n次諧波峰(n=1,2,??,5)。若5次檢測中存在3次以上諧波峰，則認(rèn)為可能存在氣泡并記錄下來；若諧波峰少于3次，則記結(jié)果為空，更換基頻f0重復(fù)此步驟，直到遍歷搜索區(qū)間內(nèi)所有峰值。若最后結(jié)果仍為空，則認(rèn)為該頻率Fi下無氣泡。

(3) 對輸出結(jié)果為可能存在氣泡的頻率值Fi進(jìn)行結(jié)果檢驗(yàn)。如果其鄰近頻段輸出結(jié)果同樣為可能存在氣泡，則認(rèn)為該頻段為氣泡振蕩頻率；如果鄰近頻段結(jié)果是無氣泡，則認(rèn)為Fi處的輸出結(jié)果是虛警予以剔除。

圖12 氣泡聲信號自適應(yīng)檢測流程圖Fig.12 Flowchart of adaptive detection of bubble sound signal

上述檢測算法流程如圖12所示。其中考慮到步驟(2)中，基頻f0的誤差會(huì)給高次諧波檢測帶來較大影響，在每次檢測到諧波后，對基頻f0進(jìn)行修正，修正公式為

其中：fn0為第n次諧波檢測后對基頻的修正值；fn為檢測到的第n次諧波的頻率。

根據(jù)聲信號主要能量的頻段范圍，對頻率F在0.1～1.5 kHz范圍內(nèi)的聲信號進(jìn)行了檢測。設(shè)置參數(shù)如下：fmin=2 Hz，fmax=100 Hz，?f=1 Hz。代入上述算法，對不同頻率下有無氣泡進(jìn)行了檢測。

不同頻段處的聲信號時(shí)間-能量及頻譜變化如圖 13、14所示。通過圖 13可以發(fā)現(xiàn)，氣泡振蕩頻段附近的時(shí)間-能量曲線，波形有較明顯的周期性，時(shí)間-能量的頻譜具有較為明顯的基頻峰和諧波峰，如圖中箭頭所指；而在噪聲干擾較強(qiáng)的頻段，如圖 14所示，能量隨時(shí)間變化沒有明顯的周期性，且其頻譜無諧波特性。通過檢測算法可以明顯區(qū)分。

圖13 存在氣泡頻段的聲信號能量隨時(shí)間和頻率的變化Fig.13 Variation of sound signal energy at bubble resonance frequencies with time and frequency

圖14 低頻段處聲信號的能量隨時(shí)間和頻率的變化Fig.14 Variation of sound signal energy at low frequencies with time and frequency

經(jīng)過本文算法在不同深度下檢測得到的氣泡聲頻率如表1所示。表1中第二列為檢測結(jié)果存在氣泡的頻段范圍；第三列為第二列數(shù)據(jù)的中值，近似認(rèn)為是氣泡振蕩的主頻率；第四列數(shù)據(jù)為不同深度氣泡聲頻率與5.5 m深處氣泡頻率的比值。Minneart[13]指出氣泡頻率是氣泡半徑和壓強(qiáng)的函數(shù)。由于漏孔大小及內(nèi)外壓差恒定，不同深度的氣泡半徑僅由環(huán)境壓強(qiáng)單一決定。故在本實(shí)驗(yàn)環(huán)境下，可以由不同深度下的環(huán)境壓強(qiáng)經(jīng) Minneart頻率公式計(jì)算得到頻率比值，如表 1中第五列所示。通過表 1中數(shù)據(jù)可以看出，本文方法檢測得到的氣泡聲頻率隨深度變化的趨勢與理論值較為吻合，表明了本文方法在實(shí)際海洋環(huán)境中檢測提取氣泡聲信號的有效性。

表1 不同深度下氣泡聲頻率Table 1 Bubble sound frequencies at different depths

3 結(jié) 論

實(shí)際海洋環(huán)境中的泄漏氣泡聲信號，聲能量較弱，頻譜與其它環(huán)境噪聲頻譜重疊，傳統(tǒng)方法較難將氣泡聲信號分離檢測出來。本文首先通過頻域積累的方法，進(jìn)行了降噪處理，提高了信噪比。并在此基礎(chǔ)上，基于持續(xù)泄漏氣泡在氣泡振蕩主頻段附近能量存在強(qiáng)周期性的特點(diǎn)，采用一種基于時(shí)頻分析的自適應(yīng)檢測算法，對不同頻率的聲信號進(jìn)行有無氣泡聲的檢測，最后對所有頻率檢測結(jié)果進(jìn)行匯總，得出氣泡聲的頻段范圍。對海上實(shí)驗(yàn)數(shù)據(jù)的處理表明，本文方法可以較好地實(shí)現(xiàn)對氣泡聲信號的檢測和提取。