方 超, 楊兆發(fā), 楊 威

(安徽理工大學(xué) 電氣與信息工程學(xué)院, 安徽 淮南 232001)

由于年限較長(zhǎng)、地質(zhì)沉降、化學(xué)腐蝕等原因，埋地輸氣管道經(jīng)常發(fā)生穿孔泄漏而造成爆炸等重大安全事故，造成嚴(yán)重的生命財(cái)產(chǎn)損失。因此，對(duì)輸氣管道泄漏檢測(cè)和精準(zhǔn)定位技術(shù)的研究具有非常重要的意義。對(duì)于大多數(shù)埋入地下的輸氣管道，基于泄漏聲波的泄漏點(diǎn)地表檢測(cè)技術(shù)是當(dāng)前研究的熱點(diǎn)。埋地輸氣管道泄漏聲波信號(hào)沿土壤介質(zhì)傳播衰減較大，信號(hào)幅值、頻率均較低[1]，采集時(shí)極易受周圍復(fù)雜環(huán)境噪聲影響。特別是埋于城市周圍的管道，通常存在交通噪聲、工業(yè)噪聲等強(qiáng)噪聲干擾，造成采集的聲波信號(hào)失真，影響后續(xù)信號(hào)處理，有必要對(duì)泄漏信號(hào)進(jìn)行降噪處理。

埋地輸氣管道泄漏聲波信號(hào)具有典型的非平穩(wěn)性。小波分析[2]具有時(shí)頻同時(shí)局部化能力，可以較好地刻畫(huà)信號(hào)的非平穩(wěn)性，常用于信號(hào)降噪。小波閾值(Wavelet Threshold, WT)是一種典型的降噪方法，但分解層數(shù)和小波基函數(shù)的選擇對(duì)其降噪效果影響很大。2014年，DRAGOMIRETKIY等[3]提出變分模態(tài)分解算法(Variational Mode Decomposition, VMD)，通過(guò)迭代搜尋變分模型最優(yōu)解來(lái)確定信號(hào)每個(gè)窄帶本征模態(tài)分量(BIMF)的中心頻率和帶寬，可實(shí)現(xiàn)信號(hào)頻域帶的有效分離。該方法解決了經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition, EMD)方法的模態(tài)混疊、端點(diǎn)效應(yīng)以及缺乏相關(guān)數(shù)學(xué)理論依據(jù)等缺點(diǎn)，同時(shí)具有更好的噪聲魯棒性和維納濾波特性。在較強(qiáng)隨機(jī)噪聲干擾下，VMD去除高頻噪聲效果良好，但仍然會(huì)殘留中低頻段噪聲。奇異譜分析(Singular Spectrum Analysis, SSA)降噪算法[4]穩(wěn)定性較高、降噪效果較為理想。相較于受閾值函數(shù)和小波基影響的小波閾值降噪，SSA法通過(guò)選取合適的奇異值進(jìn)行信號(hào)重構(gòu)，能最大限度地消除噪聲并保留有用信息。綜上所述，考慮到實(shí)際埋地輸氣管道泄漏聲波采集過(guò)程中存在較強(qiáng)隨機(jī)噪聲干擾，提出了一種基于變分模態(tài)分解和奇異譜分析(VMD-SSA)的聯(lián)合降噪方法。

1 原理方法

1.1 VMD分解原理

VMD分解過(guò)程本質(zhì)上是變分問(wèn)題的求解過(guò)程[5]，目的是在所有模態(tài)之和等于輸入信號(hào)f(t)的前提下使所有模態(tài)帶寬之和最小。通過(guò)估計(jì)各模態(tài)帶寬引出約束性變分模型：

(1)

式中：uk={u1...uK}和ωk={ω1...ωK}分別表示VMD分解的K個(gè)模態(tài)和其相應(yīng)的中心頻率。

引入Lagrange乘法算子λ(t)和懲罰因子α將上述約束性變分問(wèn)題轉(zhuǎn)換為非約束性變分問(wèn)題得到推廣的Lagrange表達(dá)式。利用乘法算子交替方向法(Alternate Direction Method of Multipliers, ADMM)迭代更新各模態(tài)分量及其頻率中心，最終求得變分問(wèn)題的最優(yōu)解。

1.2 SSA法降噪原理

含噪信號(hào)的時(shí)間序列l(wèi)=(l1,l2,…,ln)，選擇合適的窗口長(zhǎng)度W(2

L=USVT

(2)

式中：S=diag(λ1,λ2,…,λW)，U、V為正交矩陣，λ1,λ2,…,λW為奇異值。

由SSA理論可知信號(hào)中的有用分量對(duì)應(yīng)前F個(gè)奇異值，噪聲分量對(duì)應(yīng)其后數(shù)值較小的奇異值?？衫们癋個(gè)奇異值進(jìn)行矩陣重構(gòu)實(shí)現(xiàn)信號(hào)降噪，選擇的F值過(guò)小會(huì)導(dǎo)致重構(gòu)信號(hào)中有用信號(hào)成分丟失，過(guò)大則會(huì)導(dǎo)致噪聲成分增多。

常用的選取F值的方法是奇異值差分譜法[6]，其定義為：

βi=λi-λi+1,i=1,2,…,W-1

(3)

根據(jù)差分譜的定義，有用信號(hào)和噪聲的相關(guān)性不同導(dǎo)致差分譜中產(chǎn)生譜峰，差分譜譜峰蘊(yùn)含著信噪分離界限。選取最大差分譜譜峰對(duì)應(yīng)的第F個(gè)奇異值，將F之后的奇異值置零進(jìn)行SVD重構(gòu)即可實(shí)現(xiàn)信號(hào)降噪。

2 VMD-SSA聯(lián)合降噪算法

2.1 BIMF分量篩選

為了提取含噪泄漏聲波信號(hào)中的有用信息，需要篩選出VMD分解的有效BIMF分量進(jìn)行重構(gòu)。常用的互相關(guān)系數(shù)篩選原則不適用于原始信號(hào)未知的情況，具有一定局限性。本文引入能量熵[7]的變化來(lái)判別分解所得BIMF分量中信號(hào)與噪聲的分界點(diǎn)。確定臨界噪聲分量的方法是：首先計(jì)算從高頻到低頻各頻帶模態(tài)分量的能量熵；然后根據(jù)能量熵分布圖找到第一個(gè)局部最小值對(duì)應(yīng)的模態(tài)即臨界噪聲分量；最后以該臨界分量為界濾除高頻噪聲分量實(shí)現(xiàn)一次降噪。能量熵的計(jì)算公式

Hk=-PklgPk

(4)

式中：Pk表示第k個(gè)模態(tài)能量在總能量中的比值。

2.2 有效奇異值數(shù)目確定

通過(guò)奇異值差分譜理論選取F值重構(gòu)降噪信號(hào)，存在有用信息丟失現(xiàn)象。改進(jìn)算法結(jié)合主成分分析[8]的思想在奇異差分譜理論上增加新的F值選取約束條件，具體步驟如下：

(1)VMD分解所得從高頻到低頻依次排列的各BIMF分量中，噪聲對(duì)每個(gè)模態(tài)的支配力逐漸減弱。計(jì)算出臨界噪聲分量的能量E0以及濾除高頻噪聲分量后剩余BIMF分量的總能量E1，并假設(shè)每個(gè)BIMF分量中噪聲能量為E0，進(jìn)而估計(jì)出剩余模態(tài)中的信號(hào)能量貢獻(xiàn)

R=(E1-d×E0)/E1

(5)

式中：d為剩余BIMF分量數(shù)目。

(2)對(duì)余下模態(tài)重構(gòu)后進(jìn)行SSA分析并畫(huà)出奇異值差分譜，找到每個(gè)譜峰對(duì)應(yīng)的奇異值序列fi，計(jì)算對(duì)應(yīng)SVD分量的能量貢獻(xiàn)Pi(i為奇異值差分譜從左到右的譜峰序號(hào)):

(6)

(3)當(dāng)首個(gè)波峰f1對(duì)應(yīng)的能量貢獻(xiàn)小于R時(shí)，計(jì)算下一個(gè)波峰對(duì)應(yīng)的能量貢獻(xiàn)，當(dāng)Pi

2.3 聯(lián)合算法降噪步驟

埋地輸氣管道泄漏聲波信號(hào)VMD-SSA降噪步驟如下：

(1)對(duì)含噪泄漏聲波信號(hào)進(jìn)行VMD分解得到K個(gè)BIMF分量，通過(guò)不同的K值對(duì)應(yīng)的各個(gè)模態(tài)中心頻率來(lái)確定最優(yōu)分解次數(shù)；

(2)計(jì)算各模態(tài)能量熵，根據(jù)其分布圖篩選出臨界噪聲分量，以該臨界分量為界去除高頻噪聲模態(tài)分量；

(3)對(duì)余下的模態(tài)分量重構(gòu)之后進(jìn)行SSA分析，通過(guò)奇異值差分譜和能量貢獻(xiàn)確定有效奇異值個(gè)數(shù)F；

(4)將F之后的奇異值置零并進(jìn)行SVD重構(gòu)得到最終的降噪信號(hào)。

3 仿真驗(yàn)證

通過(guò)仿真驗(yàn)證改進(jìn)F值選取方法的可行性。

構(gòu)造信號(hào)：

X=2sin(2π×t)+sin(2π×1 500×t)+0.5×sin(2π×2 500×t)

(7)

信號(hào)采樣點(diǎn)1 024，采樣頻率1 024 Hz。加入2 dB的高斯白噪聲，VMD分解層數(shù)為6，懲罰因子為1 000。各模態(tài)能量熵分布如圖1所示，極小值對(duì)應(yīng)的分量為BIMF4。即BIMF4、BIMF5可視為噪聲分量加以去除，前3個(gè)含有用信號(hào)的BIMF分量重構(gòu)后進(jìn)行SSA降噪。根據(jù)式(5)計(jì)算出余下BIMF分量的信號(hào)能量貢獻(xiàn)估計(jì)值R=0.836。畫(huà)出奇異值差分譜以及進(jìn)行SSA降噪前信號(hào)頻譜如圖2所示，計(jì)算出前2個(gè)譜峰對(duì)應(yīng)的SVD分量的能量貢獻(xiàn)為0.648和0.830(均小于0.836)，第3個(gè)譜峰能量貢獻(xiàn)為0.912(大于0.836)。

圖1 各模態(tài)能量熵分布圖

(a)奇異值差分譜 (b)信號(hào)頻譜

不同F(xiàn)值重構(gòu)得到的降噪信號(hào)頻譜如圖3所示，依據(jù)奇異值差分譜最大峰值點(diǎn)選取的F值為2，結(jié)合能量貢獻(xiàn)的奇異值差分譜法搜索的F值為6，另外選取F值為5、7做對(duì)比。與圖2中SSA降噪前信號(hào)頻譜對(duì)比可知：圖3(a)中明顯丟失了300 Hz、2 500 Hz的頻率分量，圖3(b)中2 500 Hz頻率峰幅值下降了近一半，即F值為2、5重構(gòu)的降噪信號(hào)存在有用信息丟失現(xiàn)象；圖3(c)和圖3(d)中有用頻率成分尖峰都完整保留，圖3(c)頻譜更加平滑、毛躁減少，計(jì)算兩者與原始信號(hào)的相關(guān)度分別為0.961 8、0.954 4，說(shuō)明本文改進(jìn)的奇異值差分譜法能搜索最優(yōu)的F值，能較好的保留有用信號(hào)并實(shí)現(xiàn)降噪。

4 實(shí)驗(yàn)結(jié)果與分析

4.1 VMD-SSA原理驗(yàn)證

埋地輸氣管道泄漏聲波檢測(cè)實(shí)驗(yàn)使用長(zhǎng)1.3 m的鍍鋅鋼管模擬管道泄漏部分，管壁預(yù)設(shè)1.5 mm圓形泄漏孔。模擬輸氣管道埋土環(huán)境，設(shè)置埋土深度為80 cm，土質(zhì)濕度為6%。采集信號(hào)的傳感器為INV9822通用型壓電加速度傳感器，放置在土層表面。傳感器靈敏度為500 mV/g，頻率響應(yīng)范圍為0.2～2.5 kHz，采樣率設(shè)置為10 kHz?？諌簷C(jī)按照1 MPa標(biāo)準(zhǔn)壓力為儲(chǔ)氣罐充氣，充滿氣體后關(guān)閉空壓機(jī)，同時(shí)調(diào)節(jié)減壓閥至實(shí)驗(yàn)?zāi)繕?biāo)壓力，直至形成穩(wěn)定氣體泄漏。信號(hào)采集裝置采集原始泄漏聲波信號(hào)，將采集數(shù)據(jù)送至計(jì)算機(jī)處理。同一實(shí)驗(yàn)條件下，開(kāi)啟另一臺(tái)空壓機(jī)產(chǎn)生持續(xù)的強(qiáng)噪聲，用來(lái)模擬管道周圍復(fù)雜的環(huán)境噪聲。重復(fù)上述步驟，進(jìn)而獲得含噪泄漏信號(hào)。為了研究不同位置泄漏聲波降噪效果的區(qū)別，在土層表面設(shè)置3個(gè)采集點(diǎn),采集點(diǎn)與泄漏點(diǎn)的空間距離關(guān)系為：采集點(diǎn)1<采集點(diǎn)2<采集點(diǎn)3。對(duì)各采集點(diǎn)信號(hào)截取5 000點(diǎn)(0.5 s)作為降噪算法的輸入數(shù)據(jù)。

按照聯(lián)合降噪步驟，對(duì)采集點(diǎn)1含噪信號(hào)進(jìn)行VMD分解，懲罰因子α設(shè)為1 800。不同分解層數(shù)的中心頻率如表1所示，模態(tài)數(shù)K為5時(shí)各模態(tài)中心頻率相差較大，發(fā)生了欠分解現(xiàn)象；而模態(tài)數(shù)為7時(shí)出現(xiàn)了多個(gè)中心頻率相近的模態(tài)，出現(xiàn)了過(guò)分解現(xiàn)象，因此最終選取最優(yōu)分解次數(shù)為6。

表1 不同分解層數(shù)的中心頻率

根據(jù)式(4)計(jì)算每個(gè)BIMF分量的能量熵，并根據(jù)分布圖確定BIMF4、BIMF5、BIMF6為噪聲分量加以去除。對(duì)余下的BIMF分量重構(gòu)之后進(jìn)行SSA分析，通過(guò)奇異值差分譜結(jié)合能量貢獻(xiàn)搜索的F值為4，將其后的奇異值置零后進(jìn)行SVD重構(gòu)得到最終降噪信號(hào)。為了表現(xiàn)改進(jìn)方法的優(yōu)越性，使用小波閾值(WT)、變分模態(tài)分解(VMD)兩種方法進(jìn)行采集點(diǎn)1的泄漏信號(hào)降噪處理做對(duì)比。3種方法的降噪信號(hào)的時(shí)頻分析圖如圖4所示，從上至下依次是WT、VMD、VMD-SSA方法降噪后的波形圖與頻譜圖。

由圖4可知，相較于WT方法，改進(jìn)算法(VMD-SSA)降噪信號(hào)的時(shí)域波形幅值更加擬合原始信號(hào)；WT方法降噪信號(hào)在0～500 Hz頻段內(nèi)明顯失真，1 500 Hz附近的譜峰被濾去，喪失了部分有用信號(hào)成分，且在1 000～2 000 Hz頻段仍存在大量噪聲；VMD方法降噪信號(hào)中高頻噪聲基本去除，但在1 000～2 000 HZ頻段內(nèi)，VMD方法降噪后存在大量噪點(diǎn)；改進(jìn)算法可以進(jìn)一步濾除噪聲；相較于WT、VMD方法，采用改進(jìn)算法的降噪信號(hào)和原始信號(hào)契合度更高，雙譜峰特征明顯，較好實(shí)現(xiàn)對(duì)中低頻段噪聲的二次濾除。

(a)WT算法降噪信號(hào)波形圖 (b)WT算法降噪信號(hào)頻譜圖

綜上，對(duì)于噪聲干擾較大的泄漏聲波信號(hào)，VMD-SSA方法濾除噪聲的能力優(yōu)于其他2種方法。為了研究VMD-SSA方法對(duì)不同位置的泄漏聲波信號(hào)的降噪效果，對(duì)采集點(diǎn)2、采集點(diǎn)3的含噪信號(hào)進(jìn)行VMD-SSA方法降噪，畫(huà)出降噪前后信號(hào)的時(shí)頻幅值圖如圖5所示。由圖5可知兩個(gè)采集點(diǎn)的降噪信號(hào)噪點(diǎn)幾乎全部去除，時(shí)頻幅值圖尖峰明顯，較好的保留了信號(hào)特征。

圖5 采集點(diǎn)2、采集點(diǎn)3泄漏聲波信號(hào)降噪前(左)后(右)時(shí)頻幅值圖

4.2 降噪性能指標(biāo)

評(píng)價(jià)降噪效果的指標(biāo)一般有信噪比(SNR)和均方根誤差(RMSE)。降噪信號(hào)的SNR越高、RMSE越低表示降噪效果越好。SNR和RMSE的計(jì)算公式如下：

(8)

(9)

式中：f1(n)表示原始信號(hào)，f2(n)表示去噪后的信號(hào)，N表示信號(hào)長(zhǎng)度。

定量分析幾種方法的降噪性能，分別計(jì)算出3個(gè)采集點(diǎn)的含噪泄漏聲波信號(hào)通過(guò)小波閾值(WT)、變分模態(tài)分解(VMD)以及改進(jìn)算法(VMD-SSA)降噪前后的SNR和RMSE值,如表2所示。由表2可知對(duì)于不同采集點(diǎn)的泄漏信號(hào)，VMD-SSA方法相較于WT、VMD方法降噪信號(hào)的SNR有明顯提高、RMSE有明顯下降。計(jì)算出3種方法降噪信號(hào)的SNR平均提高量和RMSE平均下降量(相較于降噪前)，相較于WT、VMD方法，改進(jìn)算法降噪后的SNR增量分別提升32.3%、21.9%；RMSE減量分別達(dá)到17.6%、10.1%。

表2 不同方法降噪性能指標(biāo)對(duì)比

綜上所述：對(duì)于不同位置噪聲干擾下的泄漏信號(hào)，WT、VMD方法都有一定的降噪能力，但降噪效果相近，不能表現(xiàn)出明顯區(qū)別和提高。改進(jìn)算法(VMD-SSA)的降噪效果明顯提高，降噪優(yōu)越性得到證明。

5 結(jié) 語(yǔ)

針對(duì)埋地輸氣管道泄漏聲波信號(hào)沿土壤介質(zhì)傳播過(guò)程中易受強(qiáng)噪聲干擾這一問(wèn)題，提出了變分模態(tài)分解結(jié)合奇異譜分析(VMD-SSA)的降噪方法。實(shí)驗(yàn)分析表明該方法降噪效果良好，主要結(jié)論如下：

(1)通過(guò)能量熵值確定含噪BIMF分量并根據(jù)優(yōu)化F值選取的奇異值差分譜實(shí)現(xiàn)進(jìn)一步降噪，有用信息得到保留的同時(shí)，可以有效去除噪聲；

(2)對(duì)于不同位置采集的泄漏聲波信號(hào)，VMD-SSA法降噪能力較之WT、VMD方法，無(wú)論是直觀降噪效果還是定量降噪指標(biāo)都有一定提升；

(3)改進(jìn)方法可以從較強(qiáng)噪聲中提取出原始泄漏聲波信號(hào)，且不會(huì)造成信號(hào)失真。實(shí)驗(yàn)過(guò)程契合實(shí)際管道泄漏聲波面臨噪聲干擾這一環(huán)境，具有一定的實(shí)際應(yīng)用價(jià)值。