邢志祥，朱一龍，郝永梅，蔣軍成，盛 璘，楊 健，楊 克

(1.常州大學(xué)環(huán)境與安全工程學(xué)院，江蘇 常州 213164；2.常州港華燃?xì)庥邢薰荆K 常州 213161)

城市管道在長期使用過程中，因受內(nèi)外腐蝕等多種因素的影響而容易引發(fā)泄漏[1-2]，尤其是早期發(fā)生的微小泄漏，由于其泄漏孔徑小、泄漏流量小等特點(diǎn)，導(dǎo)致微小泄漏難以通過視覺、聽覺和管道檢測系統(tǒng)捕捉。實(shí)際工況中，一段被測管道往往存在多個泄漏源，而多點(diǎn)泄漏具有信號干擾大、信號包含內(nèi)容多等特點(diǎn)，導(dǎo)致多點(diǎn)泄漏信號難以提取、定位誤差大，給實(shí)際檢測帶來諸多問題[3]。

近年來，國內(nèi)外學(xué)者針對管道泄漏檢測和定位做了大量的研究，如章沖[4]利用光纖檢測法對管道多點(diǎn)泄漏進(jìn)行了定位檢測，但該法造價(jià)昂貴、定位精度較低；李玉星等[5]基于希爾伯特黃變換與盲源分離技術(shù)對管道泄漏聲波進(jìn)行提取、分離、定位，但其泄漏孔徑為0.6 mm，并不適用于微小泄漏的研究；楊凱[6]利用改進(jìn)的CEEMD(互補(bǔ)集合經(jīng)驗(yàn)?zāi)B(tài)分解)方法對地震含噪信號進(jìn)行去噪，解決了模態(tài)混疊問題并提高了信噪比，但未對分解產(chǎn)生的無效冗余分量進(jìn)行剔除，其精度有待提高；Ozevin等[7]基于聲發(fā)射和幾何連通法，由交叉相關(guān)函數(shù)計(jì)算出多維管網(wǎng)泄漏的位置，但該研究僅在理論上得到證實(shí)，并未出現(xiàn)實(shí)施案例；Scholey等[8]利用傳感器間距和波速依賴性對管道泄漏信號進(jìn)行時(shí)間差計(jì)算，從而產(chǎn)生定位結(jié)果，但該法需要大量彈性波訓(xùn)練點(diǎn)；Mostafapour等[9]利用小波變換和交叉時(shí)間頻譜理論，對管道聲發(fā)射信號進(jìn)行了精確定位，但無法做到多點(diǎn)泄漏定位。

以上研究雖然不同程度地解決了聲發(fā)射信號的處理和定位問題，但對于城市管道的微小泄漏，尤其是多點(diǎn)微小泄漏，尚未形成較為可靠、完善的定位技術(shù)[10]。因此，本文在管道泄漏定位現(xiàn)有研究成果的基礎(chǔ)上，基于管道泄漏聲波理論，提出了一種基于改進(jìn)CEEMD的城市管道多點(diǎn)泄漏定位分析新方法，以解決多點(diǎn)泄漏源信號的提取，并進(jìn)一步精確定位。

1 改進(jìn)的CEEMD方法

1.1 CEEMD方法原理

CEEMD是在經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)基礎(chǔ)上發(fā)展起來的一種信號分解方法，然而EMD不穩(wěn)定且存在模態(tài)混疊問題[11]，于是Wu等[12]提出了集合經(jīng)驗(yàn)?zāi)B(tài)分解(Ensemble Empirical Mode Decomposition，EEMD)方法，利用高斯白噪聲頻譜均勻分布統(tǒng)計(jì)特性加至信號中，使得信號在不同尺度上具有連續(xù)性。EEMD方法雖解決了模態(tài)混疊問題，但卻不能完全中和人為添加的白噪聲，文獻(xiàn)[13]在EEMD方法的基礎(chǔ)上，將互補(bǔ)集合經(jīng)驗(yàn)?zāi)B(tài)分解(Complementary Ensemble Empirical Mode Decomposition，CEEMD)方法應(yīng)用信號分析中，減小了EEMD人為添加白噪聲引起的重構(gòu)誤差。CEEMD方法的具體算法步驟如下[14]：

定義運(yùn)算符號Ej(·)為EMD分解的第j個固有模態(tài)函數(shù)，ωi(t)(i=1,2,…,n)是單位方差為零均值的高斯白噪聲，εk是每個階段設(shè)定的信噪比系數(shù)。

(1) 首先對原始信號x(t)加入高斯白噪聲ωi(t)，對目標(biāo)信號x(t)+ε0ωi(t)進(jìn)行n次EMD分解，獲得第一個固有模態(tài)函數(shù)IMF1(t):

(1)

(2) 然后計(jì)算得到一階殘差r1(t):

r1(t)=x(t)-IMF1(t)

(2)

(3) 在一階殘差r1(t)中添加白噪聲信號，構(gòu)造信號r1(t)+ε1E1[ωi(t)]并再次進(jìn)行分解，獲得第二個固有模態(tài)函數(shù)IMF2(t):

(3)

(4) 以此類推，計(jì)算第k階殘差rk(t)(k=2,3,…,K):

rk(t)=r(k-1)(t)-IMFk(t)

(4)

(5) 再繼續(xù)分解rk(t)+εkEk[ωi(t)]，得到第k+1個固有模態(tài)函數(shù)分量為

(5)

(6) 繼續(xù)分解直至滿足結(jié)束條件，即殘差極值最多不超過兩個，否則返回至步驟(4)～(5)中進(jìn)行計(jì)算，最終目標(biāo)信號x(t)可表示為

(6)

1.2 改進(jìn)的CEEMD方法

CEEMD方法不僅解決了EMD引起的模態(tài)混疊等問題，而且在保證EEMD分解效果的前提下，還減小了分解過程中人為添加白噪聲不能完全中合的缺陷。但CEEMD方法仍然存在兩個問題:其一是分解后的信號模態(tài)中仍殘留隨機(jī)噪聲；二是存冗余分量。

通過引用信號局部均值算子M(·)改進(jìn)CEEMD方法分解過程，對局部均值直接進(jìn)行估算，既提高了運(yùn)算效率，又可分解得到冗余分量。冗余分量的計(jì)算過程可表示為

(7)

式中：r(t)表示冗余分量，k=1,2，…，n;M(·)表示產(chǎn)生信號局部均值的算子；E(·)表示產(chǎn)生第k個模態(tài)的算子；ω是單位方差為零均值的高斯白噪聲；ε為每個階段設(shè)定的信噪比系數(shù)。

同時(shí)，通過在CEEMD過程中編譯帶通濾波器代碼，能夠把噪聲集中在較少的IMF分量上，由此，通過改進(jìn)的CEEMD方法計(jì)算整體局部均值，最終可使信號分解結(jié)果為噪聲分量、有效信號分量和冗余分量，即：

(8)

1.3 互相關(guān)峭度理論

互相關(guān)函數(shù)是描述隨機(jī)信號m(t)、n(t)在任意兩個不同時(shí)刻t1、t2的取值之間相關(guān)程度的函數(shù)，互相關(guān)函數(shù)給出了頻域內(nèi)兩個信號是否相關(guān)的判斷指標(biāo)，通過聯(lián)系兩個信號的互譜與自譜，可精確判定輸出信號來自于輸入信號[15]的相關(guān)程度，其相關(guān)系數(shù)的計(jì)算公式為

(9)

目前在統(tǒng)計(jì)學(xué)界，相關(guān)系數(shù)大小所表示的意義尚未形成統(tǒng)一觀點(diǎn)，但通常按表1確定[16]。

表1 相關(guān)系數(shù)-相關(guān)程度表示[16]Table 1 Correlation coefficient-correlation degree representation

根據(jù)互相關(guān)峭度理論，通過計(jì)算各階IMF分量與原始粗定位信號的相關(guān)程度，進(jìn)一步驗(yàn)證改進(jìn)的CEEMD方法的分解效果。

2 基于改進(jìn)CEEMD的管道多點(diǎn)泄漏定位

2.1 基于改進(jìn)CEEMD的管道多點(diǎn)泄漏定位模型建立

將改進(jìn)的CEEMD方法結(jié)合互相關(guān)峭度理論和盲源分離技術(shù)，建立了管道多點(diǎn)泄漏定位模型，見圖1。

圖1 管道多點(diǎn)泄漏定位模型Fig.1 Pipeline multi-point leakage location model

通過改進(jìn)的CEEMD方法可將管道泄漏原始信號直接分解為包含噪聲分量、有效信號分量和冗余分量三種成分的IMF分量；經(jīng)互相關(guān)峭度理論對提取的有效信號分量進(jìn)行驗(yàn)證，確保改進(jìn)的CEEMD方法的分解效果；再結(jié)合小波重構(gòu)和盲源分離技術(shù)，獲得獨(dú)立泄漏源信號；最后通過管道聲發(fā)射時(shí)差定位公式，獲得管道泄漏的精確定位。

2.2 基于改進(jìn)CEEMD的管道多點(diǎn)泄漏定位計(jì)算步驟

對于管道多點(diǎn)泄漏，上、下游傳感器接收到的信號為混合粗定位信號，基于改進(jìn)的CEEMD方法對管道多點(diǎn)泄漏進(jìn)行信號提取和定位的計(jì)算步驟如下：

(1) 通過管道多點(diǎn)泄漏模擬試驗(yàn)獲取原始粗定位信號x(t)。

(2) 采用改進(jìn)的CEEMD方法分解原始信號x(t)，獲得IMF有效信號分量。

(3) 根據(jù)互相關(guān)峭度理論，驗(yàn)證有效信號來自于真實(shí)信號的相關(guān)程度。

(4) 對驗(yàn)證的有效信號分量進(jìn)行小波軟閾值去噪，獲得觀測信號s(t)。

(5) 利用Fast ICA算法對混合的觀測信號s(t)進(jìn)行盲源分離，獲得獨(dú)立泄漏源信號h(t)，并計(jì)算獨(dú)立泄漏源的信號定位時(shí)差Δt。

(6) 計(jì)算管道泄漏點(diǎn)的位置x。

3 實(shí)例應(yīng)用與分析

3.1 金屬管道多點(diǎn)泄漏模擬試驗(yàn)

圖2 管道泄漏模擬試驗(yàn)系統(tǒng)圖Fig.2 Pipeline leakage simulation system

試驗(yàn)管道長度為50 m，管道材質(zhì)為鋼材，管道規(guī)格為DN150，運(yùn)輸介質(zhì)為壓縮空氣，見圖2。1號上游傳感器位于零點(diǎn)處，2號下游傳感器位于40 m處，距離1號傳感器16 m和30 m處設(shè)有2.0 mm的泄漏孔，見圖3。本試驗(yàn)分別在0.3 MPa、0.2 MPa、0.1 MPa壓力下進(jìn)行，共做3組試驗(yàn)，每組試驗(yàn)進(jìn)行4次，以獲得泄漏信號的原始數(shù)據(jù)。

圖3 管道泄漏模擬試驗(yàn)示意圖Fig.3 Schematic diagram of pipeline leakage simulation test

根據(jù)管道聲發(fā)射檢測原理獲取試驗(yàn)數(shù)據(jù)，隨機(jī)提取12個峰值相對集中且幅值較高的泄漏信號原始數(shù)據(jù)作為后續(xù)處理原始粗定位信號x(t)，見表2。

表2 管道兩點(diǎn)泄漏的原始數(shù)據(jù)Table 2 Pipeline two-point leakage original data

3.2 管道多點(diǎn)泄漏分析定位處理

由表2可知，管道泄漏粗定位存在較大誤差，需進(jìn)一步進(jìn)行分析處理。本文隨機(jī)選取表2中序號為1的數(shù)據(jù)進(jìn)行具體計(jì)算。

3.2.1 改進(jìn)的CEEMD方法對原始信號進(jìn)行分解

將表2中序號為1的數(shù)據(jù)上游通道的原始粗定位信號作為CEEMD的輸入量進(jìn)行分解，其分解結(jié)果見圖4。

圖4 改進(jìn)的CEEMD對上游原始信號的處理結(jié)果圖Fig.4 Processing result graph of upstream original signal by improved CEEMD

由圖4可見，上游原始信號被自適應(yīng)地分解為11階IMF分量。其中，噪聲分量主要集中在前3階IMF分量；4～7階IMF分量保持泄漏信號的幅值特性，屬于有效泄漏信號分量；8～11階IMF分量的信號特性逐漸消失，屬于冗余分量。因此，本文選擇4～7階IMF有效信號分量進(jìn)行提取。

3.2.2 互相關(guān)峭度驗(yàn)證

根據(jù)前述的互相關(guān)峭度理論，可對改進(jìn)的CEEMD方法的分解結(jié)果進(jìn)行驗(yàn)證：對所有分解得到的IMF分量，求其本身與原始粗定位信號x(t)的互相關(guān)程度，將相關(guān)性為微相關(guān)的信號分量認(rèn)定為噪聲分量和冗余分量，通過對比微相關(guān)信號分量與改進(jìn)的CEEMD方法分解得到的噪聲分量和冗余分量是否一致，以進(jìn)一步驗(yàn)證改進(jìn)的CEEMD方法的有效性。本文選取第4階清晰度較高的IMF分量進(jìn)行互相關(guān)系數(shù)計(jì)算，其計(jì)算結(jié)果見圖5。上游共11階IMF分量的互相關(guān)系數(shù)的計(jì)算結(jié)果匯總于圖6。

圖5 第4階IMF互相關(guān)系數(shù)的計(jì)算結(jié)果Fig.5 Calculation result coefficient of the forth-order IMF cross-correlation

圖6 上游各階IMF分量互相關(guān)系數(shù)的計(jì)算結(jié)果Fig.6 Calculation results of cross-correlation coefficients of upstream component IMF '

由圖6可見，前3階IMF分量的互相關(guān)系數(shù)在0.3以下，互相關(guān)系數(shù)較低；4～6階IMF分量的互相關(guān)系數(shù)較高，均高于0.8，尤其第4階IMF分量的互相關(guān)系數(shù)幾乎接近于1；第7階IMF分量的互相關(guān)系數(shù)驟降至0.3左右；8～11階IMF分量均與原始泄漏信號無相關(guān)性，其互相關(guān)系數(shù)幾乎為0。該驗(yàn)證結(jié)果與改進(jìn)的CEEMD方法的計(jì)算結(jié)果基本一致，從而證明了改進(jìn)的CEEMD方法的有效性。

3.2.3 小波軟閾值去噪

改進(jìn)的CEEMD方法分解出來的噪聲主要是來自于環(huán)境、介質(zhì)等的噪聲，而有效信號分量中還含有人為添加的白噪聲，需要借助小波軟閾值去噪。針對管道聲發(fā)射信號的特性，采用可實(shí)現(xiàn)離散小波變換的小波基函數(shù)[17]，即sym8函數(shù)；分解尺度由泄漏聲波中含有的信號種類數(shù)決定，即分解尺度N取3。對4～7階IMF真實(shí)信號分量進(jìn)行去噪處理，并選取第4階清晰度較高的IMF分量進(jìn)行小波軟閾值去噪計(jì)算結(jié)果展示，見圖7。

圖7 第4階IMF分量小波軟閾值去噪的計(jì)算結(jié)果Fig.7 Calculation results of wavelet soft threshold denoising of the forth-order IMF component

影響小波去噪效果的因素很多，選擇不同小波基、閾值、分解尺度，其去噪效果都不盡相同。因此，本文采用信噪比(SNR)和均方誤差(RMSE)指標(biāo)對小波去噪的效果進(jìn)行衡量，小波去噪前后指標(biāo)的對比結(jié)果見，表3。

表3 小波去噪前后指標(biāo)的對比Table 3 Comparison before and after wavelet denoising

由表3可知，去噪前后的SNR和RMSE指標(biāo)均發(fā)生了變化，具體來說，去噪后信號的SNR指標(biāo)變大，RMSE指標(biāo)變小，表示小波軟閾值去噪有一定的效果。另外，由圖7可見，去噪后信號的平滑度明顯增加，噪聲干擾明顯減少。對去噪后的4～7階IMF真實(shí)信號分量做重構(gòu)，得到上游觀測信號s1(t)，即：s1(t)=IMF4+…+IMF7。同理，可得到下游觀測信號s2(t)。

3.2.4 Fast ICA算法的信號盲源分離

盲源分離是指在信號與其他無關(guān)信號相互混合，且混合方式未知的情況下，將信號從混合的復(fù)雜信號中獨(dú)立地分離出來的處理方法。獨(dú)立分量分析(Fast ICA)，又稱為快速獨(dú)立成分分析算法，是迭代速度最快、效率最高的一種盲源分離算法，利用Fast ICA算法提取泄漏信號的特征，可實(shí)現(xiàn)快速收斂，無需考慮不必要的參數(shù)。

利用Fast ICA算法對經(jīng)CEEMD處理后得到的觀測信號進(jìn)行時(shí)域分析，從而可得到盲源分離后的上、下游泄漏信號，見圖8和圖9。

圖8 經(jīng)Fast ICA算法盲源分離后上游泄漏信號的 計(jì)算結(jié)果圖Fig.8 Calculation result map of upstream leakage signals after Fast ICA algorithm blind source separation

圖9 經(jīng)Fast ICA算法盲源分離后下游泄漏信號的 計(jì)算結(jié)果圖Fig.9 Calculation result graph of downstream leakage signals after Fast ICA algorithm blind source separation

3.2.5 時(shí)差計(jì)算及泄漏定位

通過計(jì)算上、下游泄漏信號采樣點(diǎn)的差值，即可求得泄漏信號傳播到兩個傳感器的時(shí)差Δt。由圖8和圖9中的泄漏信號采樣點(diǎn)可計(jì)算表2中序號為1數(shù)據(jù)兩點(diǎn)混合信號中的其中一個獨(dú)立泄漏源的信號定位時(shí)差Δt，有Δt=(957-264)÷105=0.006 93 s。

已知管道泄漏點(diǎn)的位置x由泄漏信號傳播到兩個傳感器的時(shí)差Δt和管道泄漏聲波速度v所決定，因此只要求得Δt和v即可對管道泄漏點(diǎn)進(jìn)行定位。管道泄漏聲波速度不僅受管道材料的影響，還受到不同介質(zhì)、工況的影響，且使用不同方法檢測到的管道泄漏聲波速度都有所差別。目前，國內(nèi)外尚未形成統(tǒng)一的管道泄漏聲波速度計(jì)算方法，不同研究人員檢測計(jì)算得到的管道泄漏聲波速度值也不一致。根據(jù)文獻(xiàn)[18]，壓縮空氣介質(zhì)所產(chǎn)生的聲發(fā)射波在鋼制管道中的傳播波速為850～1 050 m/s，并通過大量試驗(yàn)對典型管道泄漏聲波信號速度進(jìn)行了測量，其中對65 m長的φ159×60 mm鋼制管道的測量恰好與本文工況一致，因此最終確定管道泄漏聲波速度平均值為950 m/s。于是，根據(jù)時(shí)差定位公式計(jì)算得到定位結(jié)果，即管道泄漏點(diǎn)的位置x=(40-0.006 93×950)÷2=16.71 m。

3.2.6 未剔除冗余分量的定位計(jì)算對比

楊凱[6]采用CEEMD方法對含噪地震信號進(jìn)行了分解處理，但在其算法過程中并未對冗余分量做進(jìn)一步處理。本文運(yùn)用文獻(xiàn)[6]的信號處理方法對表2中12個原始泄漏數(shù)據(jù)進(jìn)行了定位計(jì)算處理，得到未剔除冗余分量的定位結(jié)果。將原始粗定位結(jié)果、未剔除冗余分量的定位結(jié)果與本文的精確定位結(jié)果進(jìn)行了對比，三種定位結(jié)果的相對誤差三維柱狀圖，見圖10。

圖10 三種處理方法定位的相對誤差三維柱狀圖Fig.10 Three-dimensional histogram of the relative error of three processing methods

由圖10可見，三種處理方法定位的相對誤差有著明顯的大小區(qū)分，具體來說，原始粗定位方法的相對誤差較高，其平均相對誤差為8.99%，未剔除冗余分量的定位方法的平均相對誤差為6.36%，其定位精度與原始粗定位方法相比有了一定的提升；本文精確定位方法的平均相對誤差為2.53%，與其他兩種處理方法相比，其定位精度明顯更高。

4 結(jié) 論

(1) 針對管道多點(diǎn)泄漏信號混雜、易受噪聲干擾，而造成真實(shí)泄漏信號提取困難、定位不準(zhǔn)等問題，本文提出了基于改進(jìn)的CEEMD方法的管道多點(diǎn)泄漏定位分析方法。

(2) 通過對CEEMD分解方法進(jìn)行改進(jìn)，將各階IMF分量顯示為有效信號分量、噪聲分量和冗余分量三種分量，可直接獲取有效信號，從而提高了管道泄漏檢測的效率，實(shí)現(xiàn)了管道泄漏檢測的精確定位。

(3) 將本文精確定位方法與未剔除冗余分量的定位方法和原始粗定位方法進(jìn)行了對比，結(jié)果表明本文精確定位方法和未剔除冗余分量的定位方法都能在一定程度上降低定位誤差，但本文的定位方法能將相對誤差控制在更小的范圍內(nèi)，同時(shí)適用于單點(diǎn)和多點(diǎn)泄漏，更符合泄漏管道的實(shí)際情況。

(4) 本文以圓孔為泄漏模擬點(diǎn)，未考慮泄漏孔不同形狀、大小對管道泄漏定位的影響，這可作為進(jìn)一步研究的方向。