基于改進CEEMDAN?熵方法的管道泄漏工況識別

李傳憲，逯雯雯，石亞男，杜世聰，鄭琬郁，李鵬宇

管道運輸具有輸量大、成本低、受環(huán)境影響小的優(yōu)點，在石油天然氣行業(yè)有著越來越重要的作用。管道腐蝕、施工破壞，以及人為破壞等因素使管道頻繁發(fā)生泄漏，造成了嚴(yán)重的人員傷亡、巨大的財產(chǎn)損失和環(huán)境破壞[1]，因此提高輸油管道泄漏檢測的準(zhǔn)確性十分必要。目前常用的輸油管道泄漏檢測方法主要有光纖檢測法、實時模型法及負(fù)壓波法[2]，其中負(fù)壓波法具有精度高、投資少、原理簡單的優(yōu)點，應(yīng)用最為廣泛。在管道泄漏檢測中，現(xiàn)場往往存在各種不同工況操作引起的干擾以及傳感器本身的噪聲干擾等，都會導(dǎo)致負(fù)壓波的信噪比較低，進而影響管道泄漏檢測的準(zhǔn)確性。

針對管道泄漏引起的負(fù)壓波信號的去噪研究，小波分析[3]和EMD分解[4]都是常用的方法，然而小波分析不具有自適應(yīng)性，需要選擇合適的小波基函數(shù)和閾值函數(shù)[5]，EMD分解方法雖然具有自適應(yīng)性，但是分解過程會造成模態(tài)混疊和端點效應(yīng)[6]，針對這一問題，本文提出了改進的添加成對白噪聲的完全集合經(jīng)驗?zāi)B(tài)分解算法（改進的CEEMDAN）對負(fù)壓波信號進行預(yù)處理，該方法不僅具有自適應(yīng)性，克服了模態(tài)混疊，而且消除了輔助白噪聲，提高了運算效率，使噪聲信號和泄漏信號可以更好的分離開，提高了負(fù)壓波信號的去噪效果。針對管道不同工況下信號的特征提取，能量特征[7?9]和峭度特征[10]常用于正常運行、泄漏以及調(diào)泵等差別較大工況的特征描述，而不適用于差別較小工況的特征描述，為了全面表征管道不同工況的壓力信號特征，本文提出了基于熵的特征向量提取方法，計算負(fù)壓波信號的能量熵、峭度熵以及排列熵，從能量、峭度、時間序列復(fù)雜性三個角度全面提取不同工況下的信號特征，并第一次將排列熵應(yīng)用到管道不同工況的識別中。

1 EMD、EEMD與 CEEMDAN方法介紹及仿真模擬

1.1 EMD及EEMD方法

經(jīng)驗?zāi)B(tài)分解（EMD）是希爾伯特黃變換(HHT)方法的關(guān)鍵部分，它依據(jù)信號自身的時間尺度特征進行分解，而不依賴于任何基函數(shù)，具有自適應(yīng)性[11]。EMD的實質(zhì)是按照時間尺度由小到大的順序把不同頻率的分量逐步從原始信號中分離出來的過程，得到不同頻率的分量稱為固有模態(tài)分量(IMF)[7]。

原始信號x(t)經(jīng)EMD分解后可表示為:

式中，IMFi(t)為包含不同時間尺度的固有模態(tài)分量，r(t)為殘余分量，N為固有模態(tài)分量的個數(shù)。

實際檢測到的信號具有復(fù)雜性，當(dāng)信號中出現(xiàn)間斷或跳躍性變化時，傳統(tǒng)EMD分解無法將具有不同特征時間尺度的成分完全分離開來，會使不同時間尺度的信號出現(xiàn)在同一分量中，產(chǎn)生模態(tài)混疊現(xiàn)象，進而導(dǎo)致各IMF失去原來的物理意義。集合經(jīng)驗?zāi)B(tài)分解（EEMD）利用白噪聲功率譜密度均勻分布的特點，引入輔助白噪聲，在一定程度上削弱了模態(tài)混疊現(xiàn)象[12]。

EEMD的步驟為：

（1）在原始信號x(t)中添加標(biāo)準(zhǔn)差為kε的零均值白噪聲 n(t)，得到混合信號 x（t）+kεn（t），其中，ε為原始信號的標(biāo)準(zhǔn)差，k為白噪聲幅度系數(shù)。一般k取0.01～0.50較合適，經(jīng)過大量實驗，本文取其為0.20。

（2）利用EMD對上述混合信號進行分解。

（3）每次加入不同的n(t)，重復(fù)步驟（1）、（2）。

（4）最后計算同一特征時間尺度的所有分量的均值- -- --------IMFi(t)，并將其作為最終的IMFi(t)。

1.2 CEEMDAN及其改進算法

EEMD通過多次加入白噪聲求均值的方法，在一定程度上減少了模態(tài)混疊現(xiàn)象的產(chǎn)生，但是加入的白噪聲所帶來的誤差不會因為求均值的過程完全消除，并且隨著重復(fù)加入白噪聲次數(shù)的增多，導(dǎo)致運算量大大增加。針對這個問題，引入添加自適應(yīng)白噪聲的完全集合經(jīng)驗?zāi)B(tài)分解（CEEMDAN）[13?14]。

利用EMD進行分解得到的第i個模態(tài)分量用imfi()表示，利用CEEMDAN方法分解所得到的第i個模態(tài)分量記為- -- -- ------IMFi(t)，那么 CEEMDAN方法步驟如下：

（1）CEEMDAN方法第一個固有模態(tài)分量的獲得與EEMD算法相同，即對混合信號x（t）+kεnj（t），j=1，2，…，M進行分解，并保留M次分解的第一個固有模態(tài)分量，將所得到的分量取平均值，作為CEEMDAN方法的第一個固有模態(tài)分量。

（2）白噪聲經(jīng)EMD分解得到的模態(tài)分量表示為imfi(kεnj(t))，將其與上一步計算得到的殘余信號組合得到混合信號為 ri(t)+imfi(kεnj(t))，i=1。利用EEMD對得到的混合信號進行分解，保留每次分解的第一個固有模態(tài)分量，取平均值可以得到CEEMDAN方法的第二個固有模態(tài)分量。

（3）繼續(xù)利用步驟（2）的算法對殘余分量進行分解，可得第n+1個固有模態(tài)分量及其殘余分量為：

（4）當(dāng)殘余分量滿足終止條件時，算法終止。

在CEEMDAN的算法步驟中可以發(fā)現(xiàn)，多次利用了EEMD方法計算每階模態(tài)分量。在求得每階固有模態(tài)分量時仍然需要重復(fù)大量運算以消除添加白噪聲的影響，并且添加的白噪聲很難被完全消除。針對這種現(xiàn)象，提出添加成對白噪聲的方法，也就是將所加的每個噪聲都采用正負(fù)成對的形式添加到信號中，如式（6）所示：

式中，S是原始信號；N是添加的白噪聲；M1是原始信號和正噪聲形成的混合信號；M2是原始信號和負(fù)噪聲形成的混合信號。

通過添加成對白噪聲進行輔助的方法，每添加一次白噪聲所得到的固有模態(tài)分量中，其中一部分包含了正噪聲的殘余量，另一部分包含了負(fù)噪聲的殘余量，在利用均值求最終的固有模態(tài)分量時，大部分噪聲會正負(fù)抵消，這樣可以在運算量較少的情況下，較精確的去除添加輔助噪聲的影響。

1.3 仿真模擬

以模擬信號為例對以上幾種方法進行驗證，對余弦信號x(t)=2cos(20πt)加入間歇性高頻振蕩信號，采樣點數(shù)為2 000，得到如圖1所示波形。

圖1 合成模擬信號Fig.1 The synthetic analog signals

分別利用EMD、EEMD、CEEMDAN方法和改進的CEEMDAN方法對上述合成模擬信號進行分解，其中EEMD、CEEMDAN和改進的CEEMDAN三種方法都加入幅度系數(shù)為0.20的白噪聲運算100次，得到如圖2所示的4種方法分解結(jié)果。從圖2（a）中可以發(fā)現(xiàn)，經(jīng)典EMD方法無法完全把高頻分量和低頻分量分離開，出現(xiàn)了嚴(yán)重的模態(tài)混疊現(xiàn)象。而圖 2（b）、（c）、（d）基本上都把高頻分量和低頻分量分離開了，說明EEMD、CEEMDAN和改進的CEEMDAN三種方法基本上消除了模態(tài)混疊現(xiàn)象。分別將EEMD、CEEMDAN和改進的CEEMDAN分解后的各個IMF分量重構(gòu)，并計算原始信號能量與重構(gòu)后的信號能量，Eo=4.040 7×103， EEEMD=4.043 8×103， ECEEMD=4.040 7 × 103，E改CEEMD=4.040 7× 103。發(fā)現(xiàn)EEMD的重構(gòu)信號引入了很多噪聲，重構(gòu)信號的能量與原始信號的能量差別較大，影響后續(xù)處理，而CEEMDAN和改進的CEEMDAN幾乎不引入噪聲，重構(gòu)信號的能量與原始信號的能量一樣，保證了算法的完備性。

圖2 四種方法對圖1復(fù)合信號的分解結(jié)果Fig.2 The decomposing results of Fig.1 composite signal by four methods

對比圖2（c）、（d）可以發(fā)現(xiàn)，CEEMDAN分解得到的前4個IMF分量以噪聲為主，無法識別出其中的高頻分量，而經(jīng)過改進的CEEMDAN分解后，前4個IMF分量包含的噪聲比較少，從第4個IMF分量開始已經(jīng)明顯包含高頻振蕩分量，并且IMF 5、IMF6分量的光滑性更好，說明加入成對噪聲的方法，既可以克服模態(tài)混疊，又能消除大量噪聲。對比CEEMDAN及改進的CEEMDAN方法計算每階IMF分量的迭代次數(shù)，如圖3所示，圖中藍色盒子里的紅色短線代表迭代次數(shù)的中位數(shù)，而紅色的“+”表示偏離中位數(shù)較大的點，其位置越靠上代表迭代次數(shù)越多，運算時間越長，可以發(fā)現(xiàn)改進的CEEMDAN方法的迭代次數(shù)明顯降低，運算效率大幅提高。綜合以上分析，改進的CEEMDAN方法，在克服模態(tài)混疊，消除輔助白噪聲和提高運算效率三個方面明顯好于其他三種方法，因此本文選擇改進的CEEMDAN方法進行后續(xù)處理。

圖3 CEEMDAN方法改進前后迭代次數(shù)Fig.3 Iteration times of CEEMDAN method and impr oved CEEM DAN method

2 現(xiàn)場數(shù)據(jù)中的應(yīng)用

實驗以西部某輸油管線為研究對象，輸油管線長52.6 km，管徑557 mm，設(shè)計壓力8 MPa。壓力傳感器安裝在管線兩端。在實驗中，采取人工下油的方式模擬泄漏，將管道首末兩端采集到的負(fù)壓波信號作為檢測對象。為了區(qū)別泄漏與其他正常工況操作引起的負(fù)壓波，分別對切罐、切泵等工況引起的負(fù)壓波進行采集，進行分析處理。對不同工況下的數(shù)據(jù)進行統(tǒng)一，下文所有用到的負(fù)壓波信號數(shù)據(jù)均采用了式（7）進行標(biāo)準(zhǔn)化處理：

式中，x(t)代表原始信號；x′(t)代表標(biāo)準(zhǔn)化以后的信號；min(x(t))代表原始信號的最小值；max(x(t))代表原始信號的最大值。不同工況的時域波形如圖4所示。

圖4 泄漏、切罐、切泵三種不同工況下的時域波形Fig.4 Waveforms of leakage,cutting tanks and cutting pumps

2.1 信號預(yù)處理

除了不同的工況操作會對泄漏的識別產(chǎn)生干擾以外，現(xiàn)場環(huán)境噪聲以及獲取信號的軟硬件設(shè)備的噪聲都會導(dǎo)致采集到的負(fù)壓波信噪比較低，影響泄漏工況的識別精度。因此需要對采集到的負(fù)壓波信號進行去噪預(yù)處理，以提高工況識別準(zhǔn)確性。

改進的CEEMDAN可以有效抑制模態(tài)混疊現(xiàn)象，并且可以最大程度地降低外加白噪聲對信號分解的影響，但是實際管道泄漏信號中包含了大量的噪聲，只有提取主要包含泄漏信號的分量即有效IMF分量，去除噪聲影響，才能最大程度得到有意義的特征值。管道一旦發(fā)生泄漏，其產(chǎn)生的負(fù)壓波會同時向管道上下游傳播，本文提出利用管道上下游的兩路信號進行相關(guān)分析的方法來提取有效IMF分量。在管道上下游采集到的負(fù)壓波信號中，既包含泄漏信號又包含噪聲信號，但是只有泄漏信號具有相關(guān)性，而噪聲信號互不相關(guān)。利用這一規(guī)律，將一路信號分解得到的IMF分量分別與另一路信號進行相關(guān)分析就可以提取含有泄漏信號的有效分量[15]。

以一次實驗獲得的泄漏負(fù)壓波信號為例進行分析，圖5所示分別是管道上下游采集到的壓力信號x1（n）和x2(n)。圖6所示的IMF分量圖是利用改進的CEEMDAN對上游負(fù)壓波信號x1(n)進行分解得到的，可以發(fā)現(xiàn)前5個IMF分量以噪聲為主，第6個IMF分量及以后包含的噪聲明顯減少，但是哪些是有效IMF分量不能進行明確判斷，因此需要進一步分析。

分別計算各個IMF分量與x1(n)和x2(n)的相關(guān)系數(shù)，得到表1，表中r表示原始信號經(jīng)過分解得到的最終殘余分量。從表1可以發(fā)現(xiàn)，IMF 1和IMF2與x1(n)的相關(guān)系數(shù)明顯偏低，而IMF3-IMF9與x1(n)的相關(guān)系數(shù)差別不大，這可能是因為IMF3-IMF9中包含了與x1(n)相關(guān)性較大的噪聲部分，因此無法辨別哪些是相關(guān)分量，而IMF1-IMF5與x2(n)的相關(guān)系數(shù)較其他明顯偏低，并且與圖7所示的信息相吻合，因此前5個IMF分量以噪聲為主，是應(yīng)當(dāng)去除的不相關(guān)分量，而IMF6-IMF11是應(yīng)當(dāng)與殘余分量一起參與重構(gòu)的有效分量。

圖5 管道上下游采集到的負(fù)壓波信號Fig.5 Negative pressure wave signals collected fr om upstream and downstream pipes

圖6 上游負(fù)壓波信號的IMF分量圖Fig.6 IMF components’diagrams of upstream negative pr essur e wave signal

2.2 特征向量提取

準(zhǔn)確的特征向量提取可以大大提高輸油管道泄漏檢測的準(zhǔn)確性。本文提出了基于熵的特征提取，Shannon提出的信息熵[16]是衡量系統(tǒng)隨機性的指標(biāo)，信息熵越大代表系統(tǒng)的隨機性越大。系統(tǒng)X包括多個事件 X={x1，x2，…，xn}，并且事件發(fā)生的概率為 P={p1，p2，…，pn}。那么事件所攜帶的信息量：I(xi)=-pilg pi，全部事件的信息量累加為系

表1 各個IMF分量分別與x 1(n)和x2(n)的相關(guān)系數(shù)Table 1 Correlation coefficients between IMF components and x1(n)/x2(n)

當(dāng)管道發(fā)生泄漏時，其負(fù)壓波信號所攜帶的能量與其他工況相比有很大的差異性，并且在不同的頻率段中能量分布情況會發(fā)生變化。因此，本文利用能量熵對管道泄漏工況進行識別。經(jīng)過2.1信號預(yù)處理，已經(jīng)獲得了有效IMF分量，然后對其提取能量熵，其步驟如下：

（1）分別計算各個有效IMF分量的能量：

式中，N為信號的采樣長度；xi為信號的幅值。

（2）計算有效IMF分量的能量總和：

式中，m為有效IMF的個數(shù)。

（3）用單個IMF的能量占總能量的比例作為其概率，按照Shannon信息熵的形式計算每個樣本的能量熵：

圖7所示為泄漏、切泵、切罐的22個樣本的能量熵。

圖7 泄漏、切泵和切罐的22個樣本的能量熵Fig.7 Energy entropy of 22 samples of leakage,cutting pumps and cutting tanks

由圖7可以發(fā)現(xiàn)，切泵產(chǎn)生的負(fù)壓波的能量熵小于泄漏和切罐產(chǎn)生的負(fù)壓波的能量熵，但是泄漏和切罐的能量熵沒有明顯區(qū)別。因此，能量熵可以一定程度上描述三種工況的負(fù)壓波的特征，但是不夠全面。

管道發(fā)生泄漏時，產(chǎn)生具有顯著沖擊特點的負(fù)壓波。峭度可以描述信號的沖擊特性，但是峭度值一般只能用來判斷有沒有泄漏發(fā)生，而不能用來區(qū)分其他工況。因此本文仿照能量熵提出了峭度熵的概念，既可以突出峭度對沖擊信號靈敏度高的特點，又能考察不同頻率下的沖擊特性。峭度熵的求解過程如下：

（1）分別計算各個有效IMF分量的峭度：

式中，xi為信號的幅值，為信號的均值，σ為信號的標(biāo)準(zhǔn)差，N為信號的采樣長度。

（2）計算有效IMF分量的峭度總和：

式中，m為有效IMF的個數(shù)。

（3）將單個IMF的峭度占總峭度的比例作為其概率，按照Shannon信息熵的形式計算每個樣本的峭度熵：

圖8為泄漏、切泵、切罐的22個樣本的峭度熵。

圖8 泄漏、切泵和切罐的22個樣本的峭度熵Fig.8 Kurtosis entropy of 22 samples of leakage,cutting pumps and cutting tanks

由圖8可以發(fā)現(xiàn)，切泵產(chǎn)生的負(fù)壓波的峭度熵最小，泄漏產(chǎn)生的負(fù)壓波的峭度熵最大，而切罐產(chǎn)生的負(fù)壓波的峭度熵處于中間大小，因此峭度熵可以在一定程度上用于三種工況的模式識別。

排列熵[17?18]（PE）不同于上文所述的能量熵和峭度熵，它是時間序列的復(fù)雜度的評價指標(biāo)，其值越大，時間序列的復(fù)雜度越高，而且它對時間序列的突變現(xiàn)象具有很高的敏感性。輸油管道發(fā)生泄漏或者進行其他工況操作時，壓力波往往會發(fā)生突變現(xiàn)象，因此可以用排列熵表征不同工況的特征。其求解過程如下：

（1）將2.1中得到的有效IMF分量與殘余分量一起重構(gòu)得到新序列{s(n)，n=1，2，…，N}。

（2）對新序列進行相空間重構(gòu)時，采用互信息法確定最優(yōu)時延τ為17，采用偽近鄰法[19]確定最優(yōu)嵌入維數(shù)m為2，得到重構(gòu)序列Xj=[s(j)，s(j+τ)，…，s(j+(m-1)τ)]， 其 中 下 標(biāo) j=1，2，…，N-(m-1)τ。

（3）以重構(gòu)序列Xj為行向量構(gòu)造矩陣X，并對矩陣X的每一行向量進行升序排列，得到s[j+(j1-1)τ]≤ s[j+(j2-1)τ]≤ … ≤ s[j+(jm-1)τ]。其中 j1，j2，…，jm表示重新排列后各個元素所在列的下標(biāo)，是一組符號序列。

（4）矩陣X的每一行向量，經(jīng)過重新排列后都可以得到一組符號序列 S（g）=（j1，j2，…，jm）。其中g(shù)=1，2，…，l，l≤ m!。統(tǒng)計每一組符號序列出現(xiàn)的次數(shù)，進而可以計算出每一種符號序列出現(xiàn)的概

（5）仿照Shannon熵的形式，定義時間序列{s(n)，n=1，2，…，N}的 PE為：

圖9為泄漏、切泵、切罐的22個樣本的排列熵，可以發(fā)現(xiàn)，泄漏產(chǎn)生的負(fù)壓波的排列熵明顯小于切泵和切罐產(chǎn)生的負(fù)壓波的排列熵，而切泵和切罐產(chǎn)生的負(fù)壓波的排列熵幾乎重合。這說明，排列熵可以用于泄漏工況的識別，而對于其他工況的識別效果不好。

2.3 工況識別

2.2 中提出基于熵的特征提取，分別對能量熵、峭度熵、排列熵進行了研究，將能量熵、峭度熵、排列熵作為特征向量，可以全面刻畫不同工況產(chǎn)生的負(fù)壓波的特性，利于后續(xù)的模式識別，分別以三種工況下的一組信號為例展示特征向量如表2所示。由于支持向量機（SVM）適用于小樣本數(shù)據(jù)的分類[20?21]，而現(xiàn)場實際輸油管道能提供的泄漏樣本很少，因此選用支持向量機具有較好的分類精度。構(gòu)造兩個支持向量機分類模型，分別用于區(qū)分泄漏和切罐，泄漏和切泵，將其命名為SVM 1和SVM 2。SVM 1輸出值為1表示泄漏，-1表示切罐；SVM 2輸出值為1表示泄漏，-1表示切泵。其中每種工況的22個樣本用于訓(xùn)練模型，另外22個樣本用于測試模型的分類效果。

圖9 泄漏、切泵、切罐的22個樣本的排列熵Fig.9 Permutation entropy of 22 samples for leakage,cutting pumps and cutting tanks

表2 改進的CEEMDDAN分解三種工況的特征向量Table 2 The eigenvectors of three working conditions decomposed by improved CEEMDDAN method

利用改進的CEEMDAN和EEMD方法分別對原始負(fù)壓波信號進行分解，然后求其有效IMF分量的能量熵、峭度熵以及排列熵，利用支持向量機對三種工況進行分類，得到表3。

表3 改進的CEEMDAN?SVM（C?S）和EEMD?SVM（E?S）診斷結(jié)果Table 3 Diagnostic results of improved CEEMDAN?SVM(C?S)and EEMD?SVM(E?S)

由表3可知，改進的CEEMDAN得到的分類結(jié)果均為100%，而EEMD對于泄漏和切泵的分類準(zhǔn)確率較高，但是對泄漏和切罐的分類準(zhǔn)確率較低，說明本文改進的CEEMDAN方法效果較好。

為了驗證所選特征向量的有效性，在利用改進的CEEMDAN分解后，分別單獨利用能量熵和峭度熵，能量熵和排列熵，峭度熵和排列熵對三種工況進行分類得到表4。由表4可知，能量熵和峭度熵對泄漏和切泵的識別效果較好，對泄漏和切罐的識別效果較差，說明其對不同工況的適應(yīng)性較差；能量熵和排列熵組合的工況識別結(jié)果均較差，峭度熵和排列熵組合的工況識別效果較好，但仍然不如能量熵、峭度熵和排列熵組合的診斷效果好，說明本文提出的特征向量可以全面描述不同工況的負(fù)壓波特性，具有較強的適應(yīng)性。

3 結(jié) 論

分別對EMD、EEMD、CEEMDAN和改進的CEEMDAN方法進行仿真模擬，并對現(xiàn)場實際輸油管道不同工況下的負(fù)壓波信號進行預(yù)處理，選用不同的特征向量進行分類。得到以下結(jié)論：

（1）通過對4種方法的仿真模擬研究發(fā)現(xiàn)，改進的CEEMDAN方法能有效抑制模態(tài)混疊，消除噪聲，信號所包含的信息可以比較真實的通過IMF的特征來反映。將改進的CEEMDAN應(yīng)用于輸油管道泄漏檢測中，發(fā)現(xiàn)基于改進的CEEMDAN的管道泄漏識別方法比基于EEMD的泄漏識別方法準(zhǔn)確率更高。

（2）提出基于熵的特征向量，通過對不同特征向量的模式識別效果進行對比，發(fā)現(xiàn)能量熵、峭度熵、排列熵作為不同工況的特征向量可以全面地反映不同工況的特征信息，比任何單獨兩種熵組成的特征向量效果更好，對不同工況具有更強的適應(yīng)性。