聯合去噪算法在輸氣管道泄漏聲信號中的應用

師曉開，孫首群

(上海理工大學 機械工程學院，上海 200093)

近年來，低碳環保的觀念深入人心，天然氣作為一種環保、便捷的能源越來越受到人們的青睞以及國家的重視，天然氣需求也日益增長。隨著北方農村開始普遍采用天然氣取暖政策的實施，天然氣在中國能源結構中的占比不斷增加。管道輸送因其運輸量大、連續實時、經濟安全、占地小等優點，成為了天然氣的最主要輸送方式。中國地貌廣闊，東西距離長，管道距離長，隨著使用年限的增長，發生泄漏事故是不可避免的，又因天然氣易燃易爆的特性，所以當發生天然氣泄漏時非常容易發生重大事故，并且也是溫室效應的一個主要來源。為保證天然氣管道安全運行，將危險降到最低，需要對其進行實時監測，及時排除故障，避免造成更大的損失[1-2]。為此筆者提出了一種改進的變分模態分解VMD(variational mode decomposition)算法識別天然氣管道的泄漏故障。

針對天然氣管道的泄漏預警研究，前人已經做出很多成果。曹雪偉[3]等通過小波分析對管道的泄漏信號去噪進行了研究。經驗模態分解EMD(empirical mode decomposition)算法根據數據本身的時間尺度特征的信號分解方法，在處理非線性非平穩的信號數據上具有明顯的優勢，目前也有許多學者采用EMD算法處理管道泄漏信號[4-5]。文獻[6]等采用改進的EMD算法對管道泄漏聲信號進行了增強去噪。文獻[7]提出了一種VMD算法和改進的巴氏距離的去噪方法，信號經VMD算法分解為若干個模態，然后利用改進的BD(bhattacharyya distance)優選模態，將有效模態進行重構得到去噪后的信號。文獻[8]提出了基于VMD-WT的聯合去噪方法，信號經過VMD算法分解為若干個模態，然后用小波閾值去噪(WT)算法對每個本征模態函數(IMF)進行二次去噪，使各IMF的周期特征更加明顯。文獻[9]提出了基于SG-VMD-SVD的去噪算法研究，信號通過SG平滑濾波處理后濾除了部分高頻噪聲，然后進行VMD算法分解后根據曼哈頓距離選取信號分量，并對噪聲分量進行了奇異值去噪。

1 變分模態分解

1.1 VMD算法介紹

KONSTANTIN基于經驗模態分解與變分原理，提出了一種可變尺度的信號分解方法。VMD算法是通過將輸入信號f分解為離散的子信號(模態)uk，即IMF。這些子信號在重現輸入時具有特定的稀疏性。VMD算法將本征模態函數uk(t)重新定義成一個調幅-調頻信號，如式(1)所示：

uk(t)=Ak(t)cos(φk(t))

(1)

式中：Ak(t)——uk(t)的瞬時幅值，uk(t)的瞬時頻率為ωk(t)，ωk(t)=φ′k(t)，且Ak(t)和φk(t)的導數均為非負數。

每個模態的稀疏度優先選擇為其頻譜域中的帶寬。換句話說，假設每個模態uk(t)在中心頻率ωk(t)附近大部分都是緊湊的。為了評估每一個模態的帶寬，提出以下方案：

1)對于每個IMF，經過希爾伯特變換來計算相關的解析信號，獲得其對應的邊際譜。

2)對于每個IMF，通過與調諧到各自中心頻率的估值與其指數混合，將模態的頻譜移至“基帶”。

3)通過解調信號的H高斯平滑度，即梯度的L2范數的平方來估計帶寬。產生如式(2)所示的約束變分問題：

(2)

式中：uk——各本征模態函數IMF，uk={u1，u2， …，uK}；ωk——各階IMF的中心頻率，ωk={ω1，ω2， …，ωK}； j2=-1； ?t——導數算子；δ(t)——狄拉克函數，也被稱作脈沖函數。

重構約束可以采用不同的方式解決。通常在存在噪聲的情況下，即高斯噪聲，二次懲罰是重構保真度的最常用方法。另一方面，拉格朗日乘數是嚴格執行約束的常用方法。這兩種方法的結合保證了在有限權重下二次懲罰較好的收斂性質。因此，如果同時使用二次懲罰項和拉格朗日乘數λ，這樣就將約束變分問題簡化成非約束變分問題，得到引入增廣拉格朗日矩陣，如式(3)所示：

(3)

在交替方向乘子方法(ADMM)的迭代次優化序列中，發現了原始最小化問題的解決方案是獲取到增廣拉格朗日算法L的“鞍點”。可結合Parseval/Plancherel傅里葉等距法，在譜域上解決上述問題。計算公式如下所示：

(4)

(5)

(6)

拉格朗日乘數的作用是強制約束，而二次懲罰α則改善了收斂性。如果不是以精確的重構作為目標，特別是存在不應該包括在分解中的(強烈)噪聲的情況下，則僅在降低拉格朗日乘數的同時使用二次懲罰是合適的選擇。實際上，二次懲罰本身代表與加性高斯噪聲相關的最小二乘數據保真度。拉格朗日乘數可通過將其值保持為0來有效關閉，最簡單的方法是選擇其更新參數τ=0。

基于以上，VMD算法的完整過程如下：

Step2： 令n=n+1執行整個循環；

Step3： 根據式(4)執行第一個內循環，更新uk；

Step4： 根據式(5)執行下一個內循環，更新ωk；

Step5： 根據式(6)更新λ；

Step6： 重復步驟2～5的操作，直到滿足式(7)所示條件：

(7)

式中：ε——收斂準則的閾值，ε=1×10-7。

1.2 一種參數自適應VMD方法

VMD與EMD及其改進算法的優點是基于數學變分學，其理論推導要扎實得多，對數據采樣和噪聲具有較強的魯棒性。但是，VMD分解算法中有兩個關鍵參數需要確定，即分解模式數K和二次懲罰α。建議α的中等值為2 000[6]。為了避免過分解或模態混疊現象的產生，因此必須選擇合適的K值。本文通過對比信號經快速傅里葉變換(FFT)后進行VMD分解后各模態分量的中心頻率，根據中心頻率相近會發生模態混疊原理，來確定模態數K值。

1.3 多尺度排列熵

排列熵PE(permutation entropy)是Bandt和Pompe[10]提出的，它可以表示時間序列的復雜性。這是一種分析系統復雜性的非線性動力學方法，該算法具有快速、易于實現的特點。但PE只是將相位空間矢量中的分量進行簡單的排列，導致了時間序列中振幅信息的丟失。Aziz和Arifm[11]在多尺度算法的基礎上提出了多尺度排列熵MPE(multiscale permutation entropy)概念，可以在PE的基礎上考慮時間序列的幅值信息，大幅提高了PE的魯棒性。

(8)

(9)

(10)

(11)

Hp=Hp(m)/ln(m!)

(12)

由式(12)可以明顯看出，0≤Hp≤1，且Hp的值越小的話，表示時間序列越規則。Hp的值可以反應并放大時間序列的變化。

1.4 小波變換

小波變換[12]是一種基于短時傅里葉變換形式而發展出來的信號分析方法，主要用來分析非平穩信號。小波去噪其本質是基于小波分解系數來進行的。信號的各不同頻帶上的信號和噪聲，由于其對應的小波分解系數的強度不同，篩選、保留原信號的小波分解系數，去除噪聲分布在個頻帶上的小波分解系數，然后將處理后的信號進行小波重構，即可得到去噪信號。

目前，小波去噪最主要有模極大值小波去噪、相關性去噪(屏蔽去噪法)、小波閾值去噪，平移不變量法去噪四種方法。根據曹雪偉等人[3]的研究，小波軟閾值去噪法，去噪的信噪比(SNR)值大、均方誤差(MSE)值小，且去噪時間短等優點。本文選用小波軟閾值去噪方法對VMD分解后的含噪信號進行去噪研究。

小波閾值去噪的主要步驟如下：

1)選取合適的小波基函數以及分解層，對原始輸入信號進行小波變換，這樣便得到一組小波分解系數。

2)選擇合適的閾值對小波分解系數處理。

3)重構小波系數得到去噪信號。

2 VMD-MPE-WT聯合去噪算法

首先對管道泄漏含噪信號進行分解，得到一系列包含有用分量和噪聲分量的IMF分量。因此，不能簡單地將這些IMF組件作為噪聲丟棄，或作為純信號保留。因此，這里引入MPE來區分這些IMF成分。MPE的Hp取值范圍為0～1，定義0≤Hp≤0.4為純信號，將其保留；0.4

圖1 VMD-MPE-WT聯合去噪算法流程示意

1)將原始信號通過FFT變換，對比中心頻率不混頻原則來確定模態數K，進行自適應VMD去噪。同時，通過VMD將噪聲信號分解為一系列的BLIMF，其頻率逐漸增大。由于噪聲頻率較高，主要集中在高頻區域；有用信號頻率較低，其主要集中在前幾個BLIMF區域[13]。

2)計算各BLIMF的MPE。如果MPE值小于0.4，則認為相對應的BLIMF被認為是信號主導分量。反之，相對應的BLIMF將被視為噪聲主導分量。

3)將噪聲含量較高的IMF分量通過小波閾值去噪后，再將前幾個有用信號占比高的IMF分量和經過小波去噪后的噪聲占比高的IMF分量重構，得到經VMD-WT去噪后的純凈信號。

3 仿真信號分析

為驗證筆者提出算法的實用性，仿真研究了輸氣管道的泄漏聲信號，該類信號具有典型的非平穩與隨機性。當管道發生泄漏時，往往會在泄漏發生的一瞬間產生突變，并且會在一定的時間范圍內恢復原樣。

因此，制定了在壓力為2.5 MPa的條件下進行仿真模擬泄漏研究，假定泄漏位置隨機選在第一萬個采樣點處，原信號的表達式如式(13)所示：

(13)

對原信號y(t)加入一個η的高斯白噪聲，SNR為5 dB，采樣點數為3×104點。未加入噪聲的原始泄漏信號的波形如圖2所示，在該波形上加入高斯白噪聲后的信號波形如圖3所示。

圖2 未加入高斯白噪聲的純凈信號示意

圖3 加入高斯白噪聲的含噪信號示意

根據中心頻率不混頻原則，取不同K值，得到不同分解本征模態函數的中心頻率，見表1所列。由表1可以清晰地看出，當模態個數K等于6的時候，IMF5和IMF6的中心頻率開始接近，非常容易發生混頻現象。因此，模態個數應該選取為5。

表1 不同K值對應的本征模態函數的中心頻率 Hz

對模擬泄漏信號進行VMD分解后得到5個本征模態函數時域圖，如圖4所示，以及各IMF分量的頻率幅值如圖5所示。由圖4，圖5可以看出，IMF1分量為低頻分量且與原信號相似度最高，但圖中無法直觀清楚地分辨出IMF2，IMF3，IMF4與IMF5的頻率大小。通過計算對比各分量與原信號的相關系數與方差貢獻率如圖6所示，IMF1和IMF2之間的相關系數和方差貢獻率急劇下降，是噪聲分量與信號分量的轉折點，因此選取IMF1為有效分量，其他分量定義為噪聲分量。

圖4 各本征模態函數時域示意

圖5 各本征模態函數頻率幅值示意

圖6 各分量與原信號的相關系數與方差貢獻率示意

計算各分量的多尺度排列熵如圖7所示，多尺度排列熵(MPE)Hp的取值范圍為0～1，定義0≤Hp≤0.4為純信號，將其保留，0.40.9時，其有效分量幾乎為零，可直接刪去。

圖7 各IMF多尺度排列熵示意

由圖6和圖7相關系數、方差貢獻率以及各IMF多尺度排列熵可知，IMF1，IMF2和IMF4為含噪信號，需要對其進行下一步去噪處理，而IMF3，IMF5及余項的MPE均大于0.9，證明其有用成分幾乎為零，直接刪去。

選擇不同的小波基函數可以不同程度地對信號進行去噪，通過各小波簇系的基本函數特性和小波函數圖的變化趨勢選取最佳的小波基函數。通過對管道泄漏信號消噪后SNR和MSE對比，可以選出最佳的小波去噪參數。

SNR是指系統中有用信號與噪聲的比例，定義如式(14)所示：

(14)

式中：PS——信號的有效功率；PN——噪聲的有效功率。

若SSNR=0時，則表示信號與噪聲的平均功率相等，此時信號的可靠性極低；若SSNR為負則表示噪聲信號超過信號功率的10%，有用信號已經淹沒在噪聲中，無意義。

MSE是衡量“平均誤差”的最簡單最有效的方法，反映了數據的變化程度，定義如式(15)所示：

(15)

均方誤差的數值越小，則預測模型和真實模型的誤差越小，預測精度越高。

選擇不同的小波基函數及軟閾值法以及不同的分解尺度，去噪效果會有所不同，對比db1～db10小波基，以及每個小波基對應的分解尺度1～9，對泄漏音波的仿真信號進行了去噪分析[14]。根據SNR和MSE不同結果得到，選用db2小波基函數，在其分解尺度為9的情況下，去噪效果最好。

綜上所述，小波去噪，VMD去噪，VMD-WT去噪和本文提出的VMD-MPE-WT去噪的SNR值和MSE值對比見表2所列。

表2 不同去噪方法對比

根據表2，以上三種去噪方法，根據信號去噪質量與SNR大小為正相關，與MSE大小為負相關。可以看出本文提出的VMD-MPE-WT聯合去噪方法去噪效果最好。VMD-MPE-WT聯合去噪結果與未加噪聲的仿真信號對比如圖8所示。

圖8 VMD-MPE-WT聯合去噪結果與未加噪聲的仿真信號對比示意

由圖8分析可以得出，由VMD分解信號，結合MPE進行篩選后，采用db2小波在第9層分解尺度下采用軟閾值法對泄漏信號去噪濾波效果最好，其與原始信號相似度最高，原始信號的特征信息可以得到較好的保留。

4 結束語

筆者提出的VMD-MPE-WT去噪方法，針對VMD分解算法，根據K值選取原則，在實際應用領域無可靠的理論根據的情況下，本文采用中心頻率不混頻原則對變分模態分解算法的K值進行選取，結合MPE對分解信號進行進一步篩選，去除噪聲含量過高的分量，提高了SNR；再通過結合小波軟閾值去噪算法對小波基函數和分解尺度進行了尋優對比，對含噪分量進行進一步去噪、重構，可以盡可能地減少高頻率分段的有用信號的丟失，也盡最大程度地去除泄漏信號中混入的噪聲。通過仿真信號的試驗，與其他方法對比后，根據SNR，MSE去噪性能指標驗證了筆者提出的去噪算法在油氣管道泄漏聲信號去噪領域有優越性。