基于改進最小方差無失真響應角度譜算法的氣體泄漏定位方法?

王天琪 楊祥國 晏 園 劉成龍 張金磊

(武漢理工大學船海與能源動力工程學院 武漢 430063)

0 引言

氣體的儲存及運輸在工業領域應用非常普遍，一旦發生泄漏，就會帶來巨大的浪費和能量損失，甚至威脅到人的生命安全[1-2]。因此，對氣體泄漏源進行快速準確的定位顯得尤為重要。氣體泄漏定位是一個廣泛的領域，在過去的幾十年中得到了深入的研究[3]?；诼晫W的泄漏源定位方法是通過采集泄漏源釋放的超聲信號，進而利用聲源定位算法定位泄漏源的位置，其中包括時延估計(Time delay estimation,TDE)[4]、機器學習[5]、波束形成[6]和多信號分類算法(Multiple signal classification,MUSIC)[7]等。韋娟等[8]提出的基于模擬信號時域壓縮的氣體泄漏源波達方向估計算法，通過構造形式相對簡單的等效測量矩陣，有效降低了計算復雜度，實現了較高的估計精度。Eret等[9]通過傳聲器陣列檢測泄漏所產生的可聽聲頻段(低于20 kHz)的泄漏聲波，利用波束形成方法得到泄漏聲源的波達方向(Direction of arrival,DOA)估計，然而在該頻段內噪聲會將泄漏聲源掩蓋，導致算法失效。Bian等[10]提出基于聲學信號的實時氣體泄漏定位方法，該方法使用超聲傳感器中不同陣元所接收的信號之間的時空域相關性，解決了連續泄漏超聲信號的實時定向問題。Yan 等[11]提出一種將聲發射(Acoustic emission,AE)傳感器結合MUSIC算法和小波數據包分析的新方法，對大型壓力容器上多個泄漏源的定位進行了實驗研究，具有良好的定位多個泄漏源的潛力?；诮嵌茸V[12-13]的聲源測向方法適用于噪聲環境下的聲源DOA 估計。Blandin 等[13]提出了一種基于最小方差無失真響應(Minimum variance distortionless response,MVDR)的角度譜定位算法，實現了聲源的準確定位，但該方法在進行高頻聲源定位時獲得的聚合角度譜中因空間混疊及其他實際因素影響，容易導致較多的漏警和虛警。

針對上述問題，本文利用泄漏信號的超聲特性，通過引入信噪比(Signal noise ratio,SNR)追蹤加權和基于時頻稀疏性的分頻帶處理的方式，提出基于改進MVDR角度譜算法的氣體泄漏定位方法，以期在復雜的噪聲環境中提升氣體泄漏定位的精度，增強工程的實用性。

1 MVDR角度譜方法原理

假設空間中有N個泄漏源信號，由M個陣元所組成的超聲傳聲器陣列所接收到的泄漏聲學信號可表示為

式(1)中，xm(t,f,τ) (m=1,2,···,M)表示第m個傳聲器陣元對應的觀察信號，sn(t,f) (n=1,2,···,N)表示第n個目標泄漏源的聲學信號；n(t,f)表示接收信號中包含混響和噪聲部分的加性復合噪聲；A(f,τ)=[a1(f,τ1),a2(f,τ2),···,an(f,τn)]是所有聲源信號到傳聲器陣列之間的傳播模型矩陣，an(f,τn)表示與第n個源對應的導向矢量：

式(2)中，τn為所求聲源到傳聲器兩個陣元間的到達時間差。

MVDR 波束形成是一種自適應波束形成算法，旨在使用轉向矢量對空間內所有方位角進行“掃描”，通過調整權值矢量使得在能量一定的前提下的輸出信號最大[14]。加權后的輸出響應可表示為

式(3)中，WH=[w1,w2,···,wm(t,f)]T為權值向量。

假設信號與噪聲干擾不相關，將式(1)代入式(3)：

根據式(4)可以得到聲源信號與噪聲信號的功率，如下所示：

式(6)中，w(·)是局域二維窗函數，本文采用的是海明窗。

根據式(6)可以得出作MVDR 波束形成后，各個方向的期望信號能量：

聚合角度譜的時頻域角度譜函數：

所有網格點內的函數?MVDR(t,f,τ)在時頻域上積累較大的值：

2 改進算法的構建

2.1 SNR追蹤加權

由于噪聲不可避免，聲源數較多的混響環境下，會導致MVDR 波束形成器在抑制干擾信號的同時削減期望信號，降低定位精度。本文提出一種基于Softplus 函數的SNR 追蹤加權方式，首先根據每個時間幀某一陣元接收信號的頻域協方差矩陣估計SNR，通過SNR追蹤提取受噪聲影響較小且單個聲源能量占優的時頻支撐域，進而通過Softplus 激活函數自適應的調整不同SNR 下氣體泄漏信號的頻域分量對期望信號能量P(t,f,τ)的貢獻，在充滿低頻噪聲的環境中增加氣體泄漏信號超聲頻段的頻譜權重，達到抑制加性復合噪聲的目的。

超聲傳聲器陣列的SNR輸出可表示為

假設噪聲穩定不變且前L段為純噪聲，則有

結合泄漏源所在環境給出經驗閾值μs，若該時頻支撐域的χSNR高于給定閾值，則被認為包含一個主導聲源，從而得到滿足局部SNR追蹤的時頻支撐域集合為

提取受環境噪聲影響較小的時頻域集合?SNR(t,f,τ)，進而獲得滿足泄漏能量占優的時頻域參與角度譜函數的累加。

本文采取Softplus 激活函數，Softplus 函數其理論表達式為ln(1+ex)，利用此函數對x求微分，Softplus 函數的導函數是1/(1+e-x)，可以看出該函數是一種非線性連續可微函數，其變化平緩且光滑，收斂速度相對較快，不會出現強制的稀疏處理而屏蔽一部分的有效特征。利用該函數將SNRχSNR映射為每個頻點的加權值，即

將加權值G(t,f,τ)用于角度譜函數的計算，以調整不同SNR的頻域分量對角度譜函數的貢獻，加權后的角度譜函數中泄漏信號SNR 低的頻率權重小，而SNR高的頻率權重大。

2.2 基于時頻稀疏性的分頻帶處理

聲學信號的時頻稀疏性，即W-DO 特性，一般指將時域觀察聲學信號經由短時傅里葉變換(Short time Fourier transform,STFT)到時頻域信號后，各個聲源的能量只分布在占比很小的時頻支撐域中，不同的聲源間的聲學信號時頻成分幾乎是完全錯開的，每個時頻點上僅有一個聲源占優[15-16]。利用聲學信號的時頻稀疏性，可將整個頻域內的多聲源定位問題轉化為局部頻帶內的單聲源定位問題。

將式(1)聲學信號經過STFT 后可表示為如下數學模型：

其中，Am,n(f,τ)包含直達波與混響成分，反射波的強度均勻分布于所有可能的傳播方向上，此時傳遞函數Am,n(f,τ)中的混響成分可以看作空間散布的噪聲，故可分解為直達波成分(f,τ)及混響成分(f,τ)兩部分的求和，即：

當滿足時頻稀疏特性時，且僅考慮直達波成分：

其中，αm,n和im,n/fs分別表示第n個聲源與第m個傳聲器陣元之間傳遞函數中與直達波成分相對應的傳播衰減系數和傳播時間。

設泄漏源定位時的頻率范圍設置在[fL,fH]之間，由空域奈奎斯特采樣定理，可得最大不模糊頻率：

為了避免高頻混疊現象，將頻帶劃分為[fL,fc]與[fc,fH]兩個子區間。當fH>2fc時，[fc,fH]又可分為[fc,2fc]和[2fc,fH]兩個部分，那么根據

式(18)中，U(·)為向上取整運算符，將整個頻帶劃分為[fL,fc],[fc,2fc],···,[(P-1)fc,fH]。此時結合式(8)中MVDR 角度譜函數的定義式，可得單個子頻帶內的局部角度譜函數：

基于聲學泄漏信號的時頻稀疏特性，利用頻率區間劃分的方式計算局部角度譜函數，通過時頻域累積產生角度譜函數時，抑制由低頻段能量的積累造成的噪聲源過估計，而且可以有效避免高頻混疊現象。

2.3 改進后的MVDR角度譜算法

式(19)中僅包含單個子頻帶內單個聲源的位置信息，為了準確獲得整個頻帶內所有聲源的位置信息，將2.1 節中所提到的SNR 追蹤加權運用到經分頻帶處理的局部角度譜算法中，增加高SNR的子頻帶角度譜函數權重，即增加泄漏聲源的有效權重。通過加權值G(t,f,τ)與局部角度譜函數?MVDR聯乘，并對所有子頻帶的局部角度譜函數進行累加，即可得到整個時頻域的角度譜函數：

式(20)中，NpH與NpL分別表示最高頻率與最低頻率在子頻帶中所對應的序號。

綜上所述，本文算法的相應處理步驟歸納為如下：

步驟1 將傳聲器陣列收集到的聲學信號做STFT，得到離散時頻域的表達式為xm；

步驟2 由式(10)通過SNR 追蹤求得各個時頻域的SNRχSNR，通過式(12) 提取受噪聲影響較小的時頻域集合?SNR；

步驟3 根據式(18)將時頻域集合?SNR劃分為不存在高頻混疊現象的P個子頻帶，并由式(19)求得各個離散頻帶的局部角度譜函數?MVDR；

步驟4 將每個子頻帶的加權值G(t,f,τ)與局部角度譜函數?MVDR聯乘，并將所有子頻帶的局部角度譜函數累加求和，得到ΦMVDR，再通過式(9)求得即為估計的泄漏源位置。

3 數值仿真及分析

為了驗證算法性能的優越性，利用仿真實驗平臺，參考文獻[17]，將氣體泄漏聲源的檢測頻段設為0～48 kHz。設置傳聲器陣列為直徑0.1 m 的均勻圓陣，陣元個數為16，方位角與俯仰角的角度搜索范圍均為(0?～180?)。在計算機仿真環境中，分析本文算法、傳統MVDR 角度譜算法和MUSIC算法的有效性，并對比在不同SNR和不同采樣點數下的定位性能。

3.1 算法有效性分析

設置泄漏源數N=1，采樣點數1024 個，泄漏源SNR 設為10 dB，泄漏源方向角和俯仰角分別為φ=75.00?、θ=20.00?。其仿真結果如圖1所示。

圖1 泄漏源定位仿真結果Fig.1 Leak source location simulation results

仿真定位結果顯示3 種方法均能對泄漏源進行定位。通過3 種算法的結果圖對比，可以看出傳統MVDR角度譜算法的主瓣寬度較寬，且抑制旁瓣效果較差(圖1(a))；MUSIC 算法相較傳統MVDR 角度譜算法略有優化(圖1(b))；本文算法成像結果的主瓣寬度變窄，旁瓣的能量幅值降低(圖1(c))，仿真結果顯示本文方法對噪聲的魯棒性更強，削弱了由混響造成的旁瓣，抑制了低頻噪聲信號的能量積累，增加了泄漏聲源的時頻域權重，表明定位算法的抗干擾能力得到進一步提升。

不同方位的DOA估計性能如表1 所示。由表1可以看出，在條件相同的環境下，本文算法、傳統MVDR 角度譜算法和MUSIC 算法都能進行有效的DOA 估計，但本文算法與預設結果偏差較小，改進后的MVDR 角度譜算法的DOA 估計精度得到了很大的提高。方向角的平均定位誤差約為0.60?，俯仰角的平均定位誤差約為0.62?，相較于其他兩種算法，本文算法在相同條件下的穩定性明顯提高。

表1 不同方位的泄漏源定位仿真性能Table 1 Simulation performance of leak source location in different directions

3.2 算法定位性能分析

為進一步驗證本文算法性能優勢，保證其他參數不變，比較在不同SNR 和采樣點數情況下3 種算法的定位性能。為評價算法性能，本文選擇均方根誤差(Root mean square error,RMSE)作為評價標準。

3.2.1 SNR對算法性能影響

將聲學泄漏信號的采樣點數設置為1024次，使SNR 從0～20 dB 以間隔5 dB 變化，對3 種算法分別進行多次實驗。圖2 給出了不同算法的RMSE 隨SNR變化的曲線圖。

圖2 算法RMSE 與SNR 關系Fig.2 Relation between algorithm RMSE and SNR

根據圖2 結果可知，在不同SNR 環境下，傳統MVDR 角度譜算法、MUSIC 算法和本文算法均能進行泄漏源定位，但本文方法的RMSE 均優于傳統MVDR 角度譜算法與MUSIC 算法，尤其在中低SNR 環境中，本文算法的定位RMSE 相較于其他算法提升明顯，說明該算法在低SNR 環境中的適應性較強，定位效果得到明顯改善，符合算法設計初衷。

3.2.2 采樣數對算法性能影響

將聲學泄漏信號的SNR 設置為10 dB 不變，使采樣點數從100～1000 次以間隔100 變化，對3 種算法分別進行多次實驗。圖3 給出了不同算法的RMSE隨采樣點數變化的曲線圖。

圖3 算法RMSE 與采樣點數關系Fig.3 Relation between algorithm RMSE and sampling points

由圖3 仿真結果可知：當采樣點數不同時，3 種算法都能進行有效定位，本文算法在低采樣點數下可以更好地識別出真實泄漏源位置。在相同采樣點數下，本文算法RMSE 均小于傳統MVDR 角度譜算法和最大可控響應功率波束形成法算法的RMSE，整體優于其他兩種算法，由此可見本文算法在搜索目標泄漏源時具有成功率高、穩定性強等優勢，體現出了很好的定位性能。

綜上所述，在采樣數較少和SNR 較低的情況下，引入SNR追蹤加權和基于時頻稀疏性分頻帶處理的MVDR 角度譜算法可以獲得較為滿意的定位精度，進而說明其DOA估計性能的優越性。

4 實驗驗證

4.1 實驗設置

為了保證算法在工程應用中的可行性，驗證本文算法的有效性和正確性，建立氣體泄漏源定位系統，如圖4所示。定位系統由空氣壓縮機、計算機、超聲傳聲器陣列及攝像頭組成。

圖4 氣體泄漏源定位系統Fig.4 Gas leak source location system

文中研究對象為壓力氣體泄漏源，為了模擬出實際的泄漏源，選用了型號為QTX1500X2 的空氣壓縮機，它可以提供7 bar 的泄漏壓力，裝配一個0.5 mm的噴嘴作為泄漏孔，用于模擬氣體泄漏產生聲波信號。為保證最佳響應，實驗中選用16個型號為SPH0641LU4H-1 壓電式超聲傳聲器組成傳聲器圓形陣列，用以采集泄漏源產生的超聲信號，該傳感器的采樣頻率為96 kHz，靈敏度為-26 dBFS，尺寸為8 mm×7.5 mm，具有優秀的寬帶聲頻性能和射頻抗干擾性能。同時利用聲場可視化技術[18]，通過攝像頭標定，建立聲源空間坐標與圖像像素點之間的映射關系，將空間位置的聲功率圖與攝像頭拍攝的實景圖像融合后進行顯示，形成氣體泄漏信號可視化圖，用以觀測算法估計位置與實際泄漏點的一致性。

實驗是在一個7 m×4 m×4 m 的混響室中進行的，噪聲為空間白噪聲及信號發生器模擬噪聲。在實際的氣體泄漏定位實驗中，首先采用上述的傳聲器陣列采集由模擬泄漏源產生的泄漏超聲信號，然后再通過本文定位算法對泄漏源進行方位估計。實驗前保持空壓機噴嘴和超聲傳聲器陣列在同一高度上，實驗過程中通過改變兩者之間的相對位置來模擬不同角度下的泄漏源。

4.2 實驗結果

利用本文算法分別對不同位置、不同SNR以及不同采樣點數下的實際氣體泄漏進行定位。本節設計的實驗工況及定位結果如表2 所示。以不同位置泄漏中心坐標為基準，可得本文算法的不同工況下平均方位角定位誤差分別為0.95?，平均俯仰角定位誤差0.87?，總誤差始終小于3.5?，說明該算法穩定性高且估計效果較好。

表2 實驗工況及定位結果Table 2 Experimental conditions and positioning results

取工況1、工況6、工況7、工況9 為例具體分析，由超聲傳聲器陣列第一個陣元所接收到原始信號的時域圖及頻譜圖依次如圖5 所示。從圖5(a)中可以看出，陣元與泄漏源的位置及SNR變化導致時域信號的幅值存在差異。從圖5(b)中可以看出，空壓機所模擬出的氣體泄漏信號為寬頻信號，包含了多個頻率成分，但采集到的泄漏信號中能量主要分布在超聲頻段內，證實了分頻帶處理抑制低頻噪聲的有效性。

圖5 泄漏信號的時域與頻譜Fig.5 Time domain diagram and spectrum diagram of leakage signal

對于工況1、工況6、工況7、工況9，通過改進MVDR 角度譜算法得到氣體泄漏源的定位結果如圖6 所示，搜索由方位角與俯仰角組成的整個二維平面，其中角度譜函數輸出功率的譜峰位置即為泄漏源的位置，由圖6(a)、圖6(c)可知，在SNR 較大且采樣點數較多時，算法受噪聲影響較小，能夠精準定位泄漏源位置。由圖6(b)、圖6(d)可知，在SNR較小且采樣點數較少的情況下，本文算法均能進行有效的定位，具有較強的魯棒性。氣體泄漏可視化圖如圖7 所示，根據可視化圖可以觀察到泄漏源空間位置與算法定位位置保持一致，超聲傳聲器陣列采集的實驗數據經過本文算法計算可以較為準確估計氣體泄漏位置，證實了本文算法在實際工程應用中的可行性。

圖6 氣體泄漏定位結果Fig.6 Gas leak location results

圖7 氣體泄漏可視化圖Fig.7 Gas leakage visualization

5 結論

本文針對在現有算法在噪聲環境中難以準確定位微小氣體泄漏的問題，根據氣體泄漏聲源的時頻域特性，提出了一種基于改進MVDR角度譜算法的氣體泄漏定位方法，得到了以下結論：

(1) 本文在傳統MVDR 角度譜算法的基礎上引入SNR 追蹤加權和基于時頻稀疏性的分頻帶處理，提取受噪聲影響較小且單個聲源能量占優的時頻支撐域，自適應地調整不同頻率分量對角度譜函數的貢獻，并且可以抑制低頻段能量的積累而造成的低頻噪聲源過估計，避免高頻混疊現象的產生。

(2) 仿真實驗驗證了在不同方位、不同SNR 以及不同采樣點數的情況下，對比傳統MVDR角度譜算法和MUSIC 算法，本文算法均能保證更精確穩健的定位性能。

(3) 實驗環境下對不同工況的氣體泄漏源進行定位均可以得到較為準確的結果，且總定位誤差始終小于3.5?。證明了本文算法在實際工程應用中的準確性與有效性。