VMD-SCT-GMF濾波算法

王冬梅何 彬路敬祎肖建利

(東北石油大學(xué)a.電氣信息工程學(xué)院;b.黑龍江省網(wǎng)絡(luò)化與智能控制重點(diǎn)實(shí)驗(yàn)室,黑龍江大慶 163318)

0 引 言

長輸管道作為天然氣的主要運(yùn)輸手段,管道的安全性一直是關(guān)注的重點(diǎn)。由于管道受到外界自然環(huán)境和人為破壞的影響容易發(fā)生泄漏,管道泄漏會(huì)導(dǎo)致資源浪費(fèi)和環(huán)境污染等問題,嚴(yán)重情況下還會(huì)發(fā)生爆炸事故,監(jiān)測管道是否發(fā)生泄漏已成為當(dāng)前的重點(diǎn)問題[1]。然而室外管道發(fā)生泄漏時(shí)其信號(hào)往往混雜大量噪聲信號(hào),有效濾除噪聲成為管道泄漏檢測的重要環(huán)節(jié)。

長期以來,小波分析(WT:Wavelet Transform)、經(jīng)驗(yàn)?zāi)B(tài)分解[2](EMD:Empirical Mode Decomposition)和總體平均經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD:Ensemble Empirical Mode Decomposition)廣泛應(yīng)用于處理信號(hào)降噪分析。Wang等[3]提出利用小波閾值去噪,再將去噪后的信號(hào)進(jìn)行EMD分解,并對(duì)分解后的信號(hào)進(jìn)一步地頻譜濾波后重構(gòu)得到去噪信號(hào)。但由于小波分析、EMD和EEMD處理信號(hào)時(shí)都存在一些問題,小波分析在選取小波基和閾值上有不確定性,EMD在分析信號(hào)時(shí)存在模態(tài)混疊等缺陷,EEMD雖然有效解決了EMD的模態(tài)混疊問題[4],但其分解的效率仍然不高。Dragomiretskiy等[5]提出的VMD(Variational Mode Decomposition)是一種自適應(yīng)信號(hào)分解方法,有效提高了信號(hào)的處理效率和抑制信號(hào)分解過程中產(chǎn)生的模態(tài)混疊現(xiàn)象。

王秀芳等[6]提出利用VMD處理信號(hào)得到若干個(gè)模態(tài)分量,分別計(jì)算每個(gè)模態(tài)分量的概率密度函數(shù)的能量值,根據(jù)相鄰兩個(gè)能量值之間變化選取有效模態(tài)進(jìn)行信號(hào)重構(gòu)。尹淑欣等[7]提出一種通過計(jì)算VMD分解后各個(gè)模態(tài)分量與原信號(hào)之間的互信息值,選擇其中互信息值較大的模態(tài)分量為有效模態(tài)分量重構(gòu)得到去噪信號(hào)。為了更精準(zhǔn)提取有用信號(hào),提高對(duì)信號(hào)的去噪性能,找到有用信號(hào)和噪聲信號(hào)的準(zhǔn)確分界點(diǎn)以及使用不同結(jié)構(gòu)元素組成的形態(tài)濾波器是解決問題的關(guān)鍵。筆者提出一種根據(jù)變點(diǎn)理論尋找噪聲模態(tài)分量和有效模態(tài)分量分界點(diǎn),準(zhǔn)確選取有效模態(tài)分量重構(gòu)信號(hào),然后通過采用兩種不同的結(jié)構(gòu)元素構(gòu)成廣義形態(tài)濾波器對(duì)信號(hào)進(jìn)行濾波得到去噪后的信號(hào)。

1 VMD算法原理

VMD本質(zhì)是非遞歸性和自適應(yīng)的多分辨率變分模態(tài)分解方法,具有分解精度高,分解層數(shù)少,抗模態(tài)混疊等優(yōu)點(diǎn),VMD將一個(gè)實(shí)際信號(hào)同時(shí)分解為多個(gè)有限帶寬的模態(tài),同時(shí)估計(jì)出對(duì)應(yīng)的各個(gè)模態(tài)分量的中心頻率[8],VMD算法過程由構(gòu)造變分模型和求解變分模型兩部分組成。

變分模型是在原輸入信號(hào)f與分解后所有模態(tài)函數(shù)之和的約束條件下,尋找k個(gè)有限帶寬模態(tài)函數(shù),使每個(gè)模態(tài)函數(shù)的估計(jì)帶寬之和最小[9-10]。該約束最優(yōu)化變分問題可表達(dá)為

為解決上述的約束最優(yōu)化問題,利用二次懲罰項(xiàng)和拉格朗日乘子法的優(yōu)勢,引入了增廣Lagrangian函數(shù),如下所示

其中α為罰參數(shù),λ為Lagrange乘子。

利用交替方向乘子算法進(jìn)行一系列的迭代尋找最優(yōu)解,從而將原始信號(hào)f分解為k個(gè)IMF(Intninsic Mode Function)分量,具實(shí)現(xiàn)步驟如下。

Step 2n=n+1,執(zhí)行整個(gè)循環(huán)。

Step 3k=k+1,直至k=K。對(duì)所有ω≥0,更新泛函

更新泛函ωk

Step 4 對(duì)所有ω≥0,進(jìn)行雙重提升

其中γ為噪聲容限,當(dāng)信號(hào)含有強(qiáng)噪聲時(shí),可設(shè)定γ=0達(dá)到更好的去噪效果。

Step 5 重復(fù)Step2~Step4,直至滿足如下的迭代約束條件

2 VMD-SCT-GMF濾波算法研究

2.1 歸一化自相關(guān)函數(shù)絕對(duì)值均值變點(diǎn)

管道泄漏信號(hào)經(jīng)過VMD分解后得到若干模態(tài)分量,采用歸一化的自相關(guān)函數(shù)計(jì)算噪聲模態(tài)和有效模態(tài)分量的自相關(guān)函數(shù)值,并引入統(tǒng)計(jì)學(xué)中的變點(diǎn)理論區(qū)分有效模態(tài)和噪聲模態(tài)。隨機(jī)信號(hào)的自相關(guān)函數(shù)反映了信號(hào)與其自身在不同時(shí)刻的相似度,假設(shè)Z(t)為隨機(jī)信號(hào),則其歸一化自相關(guān)函數(shù)表示為

信號(hào)經(jīng)過VMD分解后得到的噪聲模態(tài)分量,由于其在各個(gè)時(shí)刻的相關(guān)度低,使其歸一化自相關(guān)函數(shù)值在0點(diǎn)處最大,而在其他時(shí)刻數(shù)值迅速減小至0。而信號(hào)分量則具有在0點(diǎn)處最大,在其他時(shí)刻點(diǎn)具有周期性緩慢減小到0的特點(diǎn)[11]。根據(jù)噪聲與有效模態(tài)分量歸一化自相關(guān)函數(shù)值的特點(diǎn),可以通過繪制二者歸一化自相關(guān)函數(shù)圖的方法進(jìn)行有效模態(tài)與噪聲模態(tài)分界點(diǎn)的主觀判定。為客觀合理地對(duì)有效模態(tài)和噪聲模態(tài)的分界點(diǎn)進(jìn)行判定,筆者通過引入統(tǒng)計(jì)學(xué)中的變點(diǎn)理論構(gòu)建模態(tài)分界點(diǎn)檢測模型。采用CUSUM(Cumulative Sum)型估計(jì)量對(duì)數(shù)據(jù)進(jìn)行估計(jì)檢測變點(diǎn)值,假設(shè)輸入一個(gè)序列X(t),該序列的變點(diǎn)可以根據(jù)

估計(jì)求得[12]。

其中m1為變點(diǎn)值,Yτ為該序列,k為變點(diǎn)時(shí)刻,m為序列變化前的長度,T為序列的長度。

2.2 廣義數(shù)學(xué)形態(tài)濾波算法

2.2.1 形態(tài)學(xué)基本原理

設(shè)輸入f(n)為定義在F=(0,1,…,N-1)上的離散函數(shù),結(jié)構(gòu)元素g(n)為定義在G=(0,1,…,M-1)上的離散函數(shù),其中N≥M。則f(n)對(duì)g(n)的腐蝕,膨脹表達(dá)式如下[13]

其中n=0,1,…,N-1;m=0,1,…,M-1。

f(n)對(duì)g(n)開閉運(yùn)算分別為

其中?,⊕,?,·分別為腐蝕,膨脹,開和閉運(yùn)算。

2.2.2 廣義數(shù)學(xué)形態(tài)濾波器

傳統(tǒng)形態(tài)濾波器的降噪效果取決于選取的單一結(jié)構(gòu)元素的尺寸和形狀,開-閉濾波器和閉-開濾波器級(jí)聯(lián)采用相同結(jié)構(gòu)元素會(huì)使結(jié)果出現(xiàn)輸出偏倚的問題[14]。針對(duì)上述問題,筆者在開-閉濾波器和閉-開濾波器級(jí)聯(lián)選擇不同尺寸的結(jié)構(gòu)元素(三角形和余弦型兩種結(jié)構(gòu)元素),不但解決了輸出偏倚問題,而且有更佳的降噪效果。

假設(shè)g1(n)和g2(n)分別為兩種不同類型的結(jié)構(gòu)元素,則廣義形態(tài)開-閉和閉-開濾波器表示如下

廣義形態(tài)濾波器中包含傳統(tǒng)形態(tài)濾波器中的開,閉兩種運(yùn)算,所以任選式(14)、式(15)中其一濾波器依然存在輸出偏倚的問題,式(14)使輸出偏小,式(15)使輸出偏大,為抑制輸出偏倚采取加權(quán)平均的方法,如下

2.2.3 基于變點(diǎn)理論的VMD結(jié)合廣義形態(tài)濾波降噪算法

廣義形態(tài)濾波器降噪效果與結(jié)構(gòu)元素的選取有直接關(guān)系。根據(jù)管道泄漏時(shí)信號(hào)特點(diǎn),筆者采用三角形和余弦結(jié)構(gòu)元素組合形成聯(lián)合濾波器的方法進(jìn)行濾波。經(jīng)上述分析,為提高對(duì)實(shí)際工況中管道泄漏信號(hào)檢測的準(zhǔn)確率,筆者首先對(duì)管道泄漏信號(hào)進(jìn)行VMD分解,得到一系列模態(tài)分量,然后采用變點(diǎn)理論判定區(qū)分模態(tài)分量中的有效模態(tài)和噪聲模態(tài),將噪聲模態(tài)剔除后,留下有效模態(tài)重構(gòu)信號(hào),最后將有效模態(tài)分量重構(gòu)的信號(hào)通過廣義形態(tài)濾波器進(jìn)行降噪,得到去噪后的信號(hào)。筆者方法流程如圖1所示。

圖1 基于變點(diǎn)理論的VMD結(jié)合廣義形態(tài)濾波算法流程Fig.1 VMD combined with generalized morphological filtering algorithm based on variable point theory

2.2.4 去噪評(píng)價(jià)指標(biāo)

筆者選擇信噪比(SNR:Signal to Noise Ratio,RSNR),均方誤差(MSE:Mean Square Error,EMSE),均方根誤差(RMSE:Root Mean Squared Error,ERMSE)作為評(píng)價(jià)指標(biāo)。

其中PS為信號(hào)功率,PN為噪聲功率。

其中y(t)為去噪信號(hào),Y為原始信號(hào)。

3 管道泄漏信號(hào)去噪分析

筆者實(shí)驗(yàn)數(shù)據(jù)均來源于東北石油大學(xué)天然氣管道泄漏檢測模擬實(shí)驗(yàn)平臺(tái)。該實(shí)驗(yàn)臺(tái)借助HD-II型管道泄漏檢測系統(tǒng)模擬實(shí)際現(xiàn)場管道,是一種高壓氣體管道泄漏檢測裝置。實(shí)驗(yàn)環(huán)境溫度為24.3℃,管道壓力為0.5 MPa,流量為60 m3/h。實(shí)驗(yàn)軟件采用Labview編程環(huán)境,數(shù)據(jù)采集使用NI-9215采集板卡,采樣頻率為3 kHz。數(shù)據(jù)采集時(shí),利用壓電式聲波傳感器將采集的數(shù)據(jù)上傳到計(jì)算機(jī)。通過實(shí)驗(yàn)?zāi)M實(shí)際現(xiàn)場管道中發(fā)生泄漏時(shí)的工況,并把采集的數(shù)據(jù)作為本次實(shí)驗(yàn)研究的原始信號(hào),筆者根據(jù)EMD自適應(yīng)分解原始信號(hào)得到8個(gè)模態(tài)分量,于是將VMD分解參數(shù)K設(shè)為8,α設(shè)為1 000[15]。原始信號(hào)經(jīng)過VMD分解得到8個(gè)模態(tài),如圖2所示。

圖2 泄漏信號(hào)VMD分解的各模態(tài)分量Fig.2 The modal components of leakage signal VMD decomposition

由于無法直觀地從圖2中區(qū)分有效模態(tài)分量和噪聲模態(tài)分量,筆者方法首先計(jì)算出VMD分解后各模態(tài)的自相關(guān)函數(shù)絕對(duì)值均值如表1所示,結(jié)合式(8)計(jì)算出變點(diǎn)值為1.163 1,對(duì)應(yīng)表1中IMF5分量,即選取IMF1~IMF4為有效模態(tài)。圖3為各個(gè)模態(tài)分量與原始信號(hào)的互相關(guān)系數(shù),通過設(shè)定閾值判定相關(guān)系數(shù)大于閾值的模態(tài)分量作為有效模態(tài)。若閾值設(shè)定為0.3,則選取IMF1、IMF2為有效模態(tài)。閾值設(shè)定為0.14時(shí),選取IMF1~IMF3為有效模態(tài)。因此存在設(shè)定的閾值不同則選取的有效模態(tài)不同的情況。

圖3 泄漏信號(hào)互相關(guān)系數(shù)圖Fig.3 Leakage signal correlation number

表1 模態(tài)分量歸一化自相關(guān)函數(shù)絕對(duì)值均值Tab.1 The mean value of normalized autocorrelation function of modal components

圖4為各個(gè)模態(tài)分量的歸一化自相關(guān)函數(shù)值,根據(jù)噪聲信號(hào)和有用信號(hào)隨著時(shí)間變化衰減程度存在明顯差異的特點(diǎn)可見,圖4中只有IMF1的隨著時(shí)間變化緩慢衰減至0,IMF2~IMF8表現(xiàn)為迅速衰減至0,由此可判定IMF1為有效模態(tài)。

圖4 模態(tài)分量的歸一化自相關(guān)函數(shù)值Fig.4 Normalized autocorrelation function values of modal components

圖5為各模態(tài)分量對(duì)應(yīng)的頻譜圖,根據(jù)實(shí)驗(yàn)發(fā)現(xiàn)采集的泄漏信號(hào)能量通常集中在100 Hz以下[16],經(jīng)上述分析得知互相關(guān)系數(shù)法和歸一化自相關(guān)函數(shù)選取的有效模態(tài)分量至多是3個(gè)模態(tài)分量,由圖5可知未能充分選出包含有用信息的模態(tài)分量。筆者方法選擇IMF1~IMF4為有效模態(tài),能充分提取有用信息。

圖5 模態(tài)分量頻譜圖Fig.5 Modal component spectrum

綜上可知,采用互相關(guān)系數(shù)和歸一化自相關(guān)函數(shù)選擇有效模態(tài)存在一定的主觀性和不準(zhǔn)確性,筆者提出的自相關(guān)絕對(duì)值均值準(zhǔn)則選取有效模態(tài)的方法具有定量和客觀性的優(yōu)點(diǎn)。為驗(yàn)證筆者提出的去噪方法對(duì)實(shí)際泄漏信號(hào)去噪的有效性,將采集的泄漏信號(hào)利用VMD-SCT-GMF進(jìn)行去噪得到去噪后信號(hào),如圖6所示。并且表2中列出筆者方法與VMD結(jié)合豪斯多夫距離(HD:Hausdorff Distance)選擇有效模態(tài)重構(gòu)信號(hào)去噪方法(VMD-HD)、將VMD與互相關(guān)(CC:Cross-Correlation)結(jié)合選取有效模態(tài),再用小波去噪的方法(VMD-CC-WT),以及采用互信息(MI:Mutual Information)選取有效模態(tài)重構(gòu)信號(hào),通過小波降噪的方法(VMD-MI-WT)或VMD與互信息結(jié)合降噪方法(VMD-MI)的實(shí)驗(yàn)結(jié)果對(duì)比。

圖6 去噪后的信號(hào)Fig.6 Signal after denoising

表2 實(shí)驗(yàn)結(jié)果對(duì)比Tab.2 Comparison of experimental results

由表2可知,去噪性能指標(biāo)中以信噪比為例,筆者算法VMD-SCT-GMF比VMD-HD算法相比去噪性能提高113%,與VMD-CC-WT算法相比去噪性能提高3%。在均方誤差方面VMD-HD、VMD-MI、VMDMI-WT、VMD-CC-WT分別為0.033 9、0.033 6、0.012 7、0.012 0,筆者所提方法VMD-SCT-GMF的均方誤差為0.011 4,其值均小于其他方法。筆者方法均方根誤差為0.633 9也是最小值。通過上述分析可知VMD-SCT-GMF算法的去噪效果更佳,驗(yàn)證了筆者提出的去噪方法對(duì)實(shí)際信號(hào)處理的有效性。

4 結(jié) 語

筆者通過引入統(tǒng)計(jì)學(xué)變點(diǎn)理論結(jié)合自相關(guān)函數(shù)的方法,有效地解決了當(dāng)VMD分解后無法客觀地、定量區(qū)分其中有效模態(tài)和噪聲模態(tài)的問題。對(duì)比分析采用互相關(guān)函數(shù)和自相關(guān)函數(shù)選取有效模態(tài)時(shí)存在判定差異的情況,筆者方法采用變點(diǎn)理論能定量地判定有效模態(tài)和噪聲模態(tài)的分界點(diǎn),為選取有效模態(tài)提供了一種可行方法。在提高泄漏信號(hào)檢測精度,提升實(shí)際泄漏信號(hào)的去噪效果方面,采用定量選取有效模態(tài)結(jié)合廣義形態(tài)濾波器對(duì)重構(gòu)信號(hào)進(jìn)行降噪的方法。經(jīng)過實(shí)驗(yàn)對(duì)比分析,驗(yàn)證了筆者方法的有效性。