朱一龍，邢志祥，郝永梅，吳 潔，嚴(yán)欣明，岳云飛

(1.常州大學(xué) 環(huán)境與安全工程學(xué)院，江蘇 常州 213164；2.江蘇特種設(shè)備安全監(jiān)督檢驗(yàn)研究院常州分院，江蘇 常州 213161)

0 引言

城市管道發(fā)生泄漏往往會(huì)造成較大的經(jīng)濟(jì)損失及環(huán)境污染，甚至可能危及人身安全。因此，管道泄漏檢測(cè)技術(shù)受到廣泛關(guān)注，并得到快速發(fā)展。目前，管道泄漏檢測(cè)方法主要有聲發(fā)射檢測(cè)法、壓力梯度法、負(fù)壓波法等，其中聲發(fā)射檢測(cè)法憑借其靈敏度高、適應(yīng)性強(qiáng)等特點(diǎn)得到廣泛應(yīng)用[1]。但聲發(fā)射信號(hào)中混雜有大量噪聲，定位誤差大，需對(duì)其進(jìn)行信號(hào)處理，處理方法有小波閾值去噪、經(jīng)驗(yàn)?zāi)B(tài)分解、奇異值分解等。然而以上方法都存在的缺陷有：小波閾值去噪會(huì)在重構(gòu)信號(hào)時(shí)產(chǎn)生振蕩，容易丟失有用信號(hào)[2]；經(jīng)驗(yàn)?zāi)B(tài)分解存在模態(tài)混疊、端點(diǎn)效應(yīng)等缺點(diǎn)[3]；奇異值分解通過(guò)逆變換重構(gòu)信號(hào)，雖然有較好地去噪效果，但計(jì)算效率過(guò)低[4]。

針對(duì)以上不足，提出1種改進(jìn)的集合經(jīng)驗(yàn)?zāi)B(tài)分解(Improved Ensemble Empirical Mode Decomposition，IEEMD)樣本熵分析的管道泄漏定位法。運(yùn)用改進(jìn)的集合經(jīng)驗(yàn)?zāi)B(tài)分解處理管道原始泄漏信號(hào)，得到各階信號(hào)分量；然后采用樣本熵(Sample Entropy，SampEn)分析法對(duì)各階信號(hào)分量?jī)?yōu)選重構(gòu)，獲得有效泄漏信號(hào)；最后結(jié)合互相關(guān)時(shí)延計(jì)算對(duì)泄漏點(diǎn)進(jìn)行精確定位。

1 泄漏信號(hào)提取

1.1 IEEMD信號(hào)分解

EEMD(集合經(jīng)驗(yàn)?zāi)B(tài)分解)是1種針對(duì)非線性、非平穩(wěn)信號(hào)的分析方法，此分析方法將信號(hào)分解為一系列振蕩函數(shù)，強(qiáng)大的自適應(yīng)特征使其廣泛應(yīng)用于信號(hào)處理領(lǐng)域[5]。但EEMD存在2個(gè)問(wèn)題：1)分解后的信號(hào)模態(tài)中仍殘留噪聲分量和冗余分量；2)運(yùn)行時(shí)間過(guò)長(zhǎng)[6]。

針對(duì)以上問(wèn)題，在EEMD中添加自相關(guān)函數(shù)計(jì)算和EMD過(guò)程，提出IEEMD算法。已知在EEMD過(guò)程中，高頻噪聲信號(hào)分量會(huì)被首先分離出來(lái)，且噪聲信號(hào)的自相關(guān)函數(shù)具有快速衰減到零的特點(diǎn)。因此，通過(guò)計(jì)算各階信號(hào)分量的自相關(guān)函數(shù)識(shí)別噪聲信號(hào)分量(判斷噪聲分離與非噪聲分離的分界)。當(dāng)這些會(huì)引起模態(tài)混疊噪聲信號(hào)分量被識(shí)別、分離后，剩余信號(hào)無(wú)需再進(jìn)行完整的EEMD分解，于是改用EMD對(duì)剩余信號(hào)分解，旨在節(jié)省時(shí)間、提高信號(hào)分解效率。IEEMD算法的具體步驟為：

1)在管道原始泄漏信號(hào)x(t)中加入1對(duì)符號(hào)相反的白噪聲，添加白噪聲的對(duì)數(shù)為N，得到2N個(gè)添加白噪聲后的信號(hào)：

(1)

式中：mi(t)為第i個(gè)添加白噪聲后的信號(hào)；x(t)為原始泄漏信號(hào)；ni(t)為添加的第i個(gè)白噪聲信號(hào)；“+”和“-”表示符號(hào)相反。

2)利用EEMD篩分原理分解mi+(t),mi-(t)，得到2N個(gè)一階信號(hào)分量ui,+1和ui,-1。原始泄漏信號(hào)x(t)第1個(gè)信號(hào)分量u1為2N個(gè)一階信號(hào)分量的集成平均值，即：

(2)

3)計(jì)算u1的自相關(guān)函數(shù)，判定其是否為噪聲信號(hào)分量。若是噪聲信號(hào)分量，繼續(xù)進(jìn)行下一步迭代；若不是噪聲信號(hào)分量，停止迭代。

4)對(duì)剩余函數(shù)ri+(t)及ri-(t)進(jìn)行更新，并作為新的分解對(duì)象重復(fù)以上步驟，得到第2個(gè)信號(hào)分量u2：

(3)

5)當(dāng)?shù)趐個(gè)信號(hào)分量up不是噪聲信號(hào)分量時(shí)，停止迭代。從原始泄漏信號(hào)x(t)中減去步驟1)～ 4)中得到的p-1個(gè)信號(hào)分量后，改用EMD分解剩余信號(hào)，最終可得到：

(4)

1.2 基于樣本熵的有效泄漏信號(hào)選取

Shannon第1次將熵的概念引入到信號(hào)處理中，證明熵與信息內(nèi)容的不確定程度有等價(jià)關(guān)系，并把此種不確定程度稱為信息熵[7]。Koimogorov在此基礎(chǔ)上，提出K熵[8]。但K熵計(jì)算條件要求苛刻，對(duì)噪聲過(guò)于敏感?；诖?，Pinus提出了抗噪能力強(qiáng)的近似熵，并取得良好效果[9]。但近似熵存在自匹配問(wèn)題，導(dǎo)致其計(jì)算結(jié)果依賴于數(shù)據(jù)長(zhǎng)度。為了克服這一缺點(diǎn)，Richman提出樣本熵[10]，其算法比近似熵更簡(jiǎn)單，計(jì)算量更小，所需時(shí)間更短。同時(shí)，由于克服了自匹配問(wèn)題，樣本熵對(duì)數(shù)據(jù)長(zhǎng)度的依賴更小。

根據(jù)樣本熵上述特點(diǎn)，將樣本熵作為管道泄漏故障信息提取的特征參數(shù)。任意x(t)={x(1),x(2),…,x(N)}的樣本熵計(jì)算過(guò)程為：

1)由信號(hào)序列構(gòu)成1個(gè)m維向量：

(5)

2)定義xm(i)和xm(j)之間的距離：

(6)

3)給定容限r(nóng)，統(tǒng)計(jì)每個(gè)xm(i)對(duì)應(yīng)的d[xm(i),xm(j)]≤r的數(shù)目，記為Ai。把Ai與N-m+1的比值記為：

(7)

4)求出Bim(r)的平均值Bm(r)：

(8)

5)用同樣的方法求出Bm+1(r)，則該時(shí)間序列的樣本熵定義為：

(9)

經(jīng)過(guò)IEEMD分解后，真實(shí)信號(hào)分量中還夾雜有冗余分量，冗余分量會(huì)干擾真實(shí)信號(hào)提取的完備性，必須將其剔除[11]。通過(guò)計(jì)算各階信號(hào)分量樣本熵實(shí)現(xiàn)信號(hào)分量的優(yōu)選重構(gòu)：樣本熵值越大，對(duì)應(yīng)其非穩(wěn)定時(shí)間序列的復(fù)雜程度越高，即片段中包含的有效泄漏信息越多；冗余分量表現(xiàn)為周期震蕩函數(shù)，其復(fù)雜程度不高，即樣本熵值較低。因此，選擇樣本熵值較大的信號(hào)分量重構(gòu)，即：

h(t)=ua+ub+…+un

(10)

式中：h(t)為有效泄漏信號(hào)；ua,ub,un為優(yōu)選的信號(hào)分量。

2 泄漏定位方法

2.1 互相關(guān)時(shí)延計(jì)算

互相關(guān)時(shí)延法是計(jì)算2個(gè)信號(hào)相關(guān)性的常用方法[12]，采用互相關(guān)時(shí)延法分析時(shí)間延遲與頻率之間的關(guān)系，進(jìn)而計(jì)算得到上、下游有效泄漏信號(hào)到達(dá)2個(gè)傳感器的時(shí)差Δt。

互相關(guān)時(shí)延計(jì)算管道泄漏信號(hào)時(shí)差的具體步驟為：

1)利用IEEMD分解法，分別得到上、下游有效泄漏信號(hào)h1(t)和h2(t)；然后計(jì)算2個(gè)有效泄漏信號(hào)的互相關(guān)函數(shù)：

(11)

式中：τ代表2個(gè)信號(hào)的時(shí)間延遲。

2)根據(jù)互相關(guān)時(shí)域分析圖中最大相關(guān)峰值計(jì)算得到2路有效泄漏信號(hào)時(shí)差Δt：

Δt=M×fs

(12)

式中：M為得到的互相關(guān)函數(shù)中幅值最大相關(guān)峰所對(duì)應(yīng)的采樣點(diǎn)數(shù)；fs為信號(hào)采集儀器設(shè)置的采樣頻率。

2.2 聲發(fā)射時(shí)差定位法

管道泄漏定位采用聲發(fā)射時(shí)差定位原理進(jìn)行計(jì)算[13]，見(jiàn)式(13)：

(13)

式中：x為被檢測(cè)管道泄漏點(diǎn)的位置，即泄漏點(diǎn)到上游傳感器距離，m；L為上下游傳感器之間距離，m；Δt為泄漏信號(hào)到達(dá)2個(gè)傳感器的時(shí)間差，s；v為泄漏信號(hào)在管道中的傳播速度，m/s。

3 管道多點(diǎn)泄漏試驗(yàn)

3.1 試驗(yàn)簡(jiǎn)介

試驗(yàn)現(xiàn)場(chǎng)有3根100 m長(zhǎng)的鋼制管道，其規(guī)格由上至下依次為DN65,DN90和DN150。本文選取管徑為DN65的鋼管作為研究對(duì)象，并通過(guò)空氣壓縮機(jī)向管內(nèi)提供0.3 MPa的恒定低壓氣流，以模擬管道正常運(yùn)行，如圖1所示。

圖1 管道試驗(yàn)系統(tǒng)Fig.1 Pipeline test system

3.2 有效泄漏信號(hào)的提取

本文通過(guò)聲發(fā)射檢測(cè)系統(tǒng)采集管道泄漏原始數(shù)據(jù)，主要設(shè)備包括：1臺(tái)美國(guó)PCA32通道聲發(fā)射采集儀、2條光纜和2個(gè)PXR04高靈敏度聲發(fā)射感器。其中，1號(hào)傳感器位于管道首端0 m處，2號(hào)傳感器位于管道末端100 m處，在距管道首端49 m處的位置設(shè)有1個(gè)直徑為2.0 mm的泄漏孔，如圖2所示。

圖2 模擬管道試驗(yàn)系統(tǒng)Fig.2 Simulated pipeline test system

通過(guò)聲發(fā)射傳感器接收管道泄漏信號(hào)并由光纜傳遞至聲發(fā)射采集儀，通過(guò)聲發(fā)射采集儀控制管道泄漏原始數(shù)據(jù)的采集、顯示、分析和保存。將聲發(fā)射采集儀的采樣精度設(shè)置為8 bit，采樣點(diǎn)數(shù)設(shè)置為2 048，采樣頻率設(shè)置為100 kHz，待管內(nèi)氣體平穩(wěn)流動(dòng)時(shí)采集管道泄漏原始數(shù)據(jù)。隨機(jī)提取1組數(shù)據(jù)進(jìn)行處理分析，其上游原始泄漏信號(hào)x1(t)如圖3所示。

圖3 上游原始泄漏信號(hào)Fig.3 Upstream original leakage signals

通過(guò)IEEMD分解，將上游原始泄漏信號(hào)x1(t)自適應(yīng)地由高頻到低頻分解為8階信號(hào)分量，如圖4所示。計(jì)算1～8階信號(hào)分量的樣本熵，結(jié)果如圖5所示。

圖4 上游原始泄漏信號(hào)IEEMD分解Fig.4 IEEMD decomposition of upstream original leakage signals

圖5 樣本熵計(jì)算結(jié)果Fig.5 Calculation results of sample entropy

由圖5可知，前3階信號(hào)分量的樣本熵值較小，均在0.3以下；4～6階的樣本熵值較大，幾乎接近于1；7～8階的樣本熵值非常小，可忽略不計(jì)。因此，根據(jù)基于樣本熵的有效泄漏信號(hào)選取法對(duì)各階信號(hào)分量?jī)?yōu)選，選擇4～6階信號(hào)分量重構(gòu)獲得上游有效泄漏信號(hào)h1(t)，如圖6所示。

圖6 上游有效泄漏信號(hào)Fig.6 Upstream effective leakage signals

同理，對(duì)下游原始泄漏信號(hào)同樣進(jìn)行如上所述的處理，得到下游有效泄漏信號(hào)h2(t)。

3.3 泄漏定位計(jì)算

根據(jù)互相關(guān)時(shí)延與聲發(fā)射時(shí)差定位法，通過(guò)互相關(guān)時(shí)域分析計(jì)算上、下游有效泄漏信號(hào)h1(t),h2(t)的時(shí)差，結(jié)果如圖7所示。

圖7 互相關(guān)時(shí)域分析Fig.7 Cross-correlation time domain analysis

由圖7可知，在0～400 kHz頻率范圍內(nèi)，采樣點(diǎn)數(shù)為963時(shí)幅值相關(guān)峰取得最大值，即互相關(guān)時(shí)延計(jì)算所需的采樣點(diǎn)數(shù)M=963。

將采樣點(diǎn)數(shù)M=963代入式(12)求得上、下游有效泄漏信號(hào)時(shí)差Δt=0.009 63 s，再將時(shí)差Δt和管道泄漏聲速v代入式(13)即可精確計(jì)算泄漏點(diǎn)位置。

管道泄漏聲速v不僅受管道材質(zhì)的影響，還受到不同輸送介質(zhì)、不同運(yùn)行工況的影響，且使用不同方法檢測(cè)到的泄漏聲速都有所差別。目前，國(guó)內(nèi)外尚未形成統(tǒng)一的泄漏聲速計(jì)算方法，各研究人員檢測(cè)計(jì)算得到的數(shù)值也不一致。根據(jù)沈功田《聲發(fā)射檢測(cè)技術(shù)及應(yīng)用》[14]可知，壓縮空氣介質(zhì)所產(chǎn)生的聲發(fā)射波在鋼制管道中的傳播速度為850～1 050 m/s，并通過(guò)大量實(shí)驗(yàn)對(duì)典型管道泄漏聲速進(jìn)行了測(cè)量，其中對(duì)1根95 m長(zhǎng)的DN65鋼制管道的測(cè)量與本文工況基本一致，最終確定其聲速平均值為980 m/s。于是，根據(jù)式(13)計(jì)算得到泄漏點(diǎn)精確位置，即x=(100-0.009 63×980)/2=45.28 m。

3.4 結(jié)果對(duì)比與誤差分析

本文采用基于IEEMD樣本熵分析的方法提取有效泄漏信號(hào)，影響有效泄漏信號(hào)提取結(jié)果的因素較多，選擇不同分解尺度，設(shè)定不同閾值，提取到的有效泄漏信號(hào)均不同[15]。

為驗(yàn)證所述方法的信號(hào)提取效果，針對(duì)章節(jié)3.2中的上游原始泄漏信號(hào)x1(t)，采用3種方法進(jìn)行有效泄漏信號(hào)提取，并采用信噪比(SNR)和均方誤差(RMSE)指標(biāo)對(duì)提取效果進(jìn)行衡量，3種方法及有效泄漏信號(hào)提取效果對(duì)比結(jié)果見(jiàn)表1。其中，SNR為衡量真實(shí)信號(hào)與噪聲之間所占比例大小的指標(biāo)，RMSE為反映原始泄漏信號(hào)與有效泄漏信號(hào)之間差異程度的度量。

表1 有效泄漏信號(hào)提取效果對(duì)比Table 1 Comparison on extraction effect of effective leakage signals

由表1可知，3種信號(hào)提取方法的SNR和RMSE指標(biāo)均不同。IEEMD+樣本熵分析法(本文方法)的SNR指標(biāo)最大且RMSE指標(biāo)最小。這說(shuō)明經(jīng)本文方法處理后，信號(hào)中真實(shí)信號(hào)所占比例最大，且該信號(hào)與原始泄漏信號(hào)的相關(guān)程度最大，即本文方法有效泄漏信號(hào)提取效果最好。

為進(jìn)一步驗(yàn)證所述方法的泄漏定位效果，隨機(jī)選取30組聲發(fā)射原始泄漏信號(hào)，仍采用表1的3種方法進(jìn)行泄漏定位并計(jì)算其相對(duì)定位誤差，結(jié)果如圖8所示。

圖8 相對(duì)定位誤差統(tǒng)計(jì)匯總Fig.8 Statistical summary of relative localization error

由圖8可知，3種定位方法相對(duì)誤差大小有明顯區(qū)分。本文提出方法(IEEMD+樣本熵分析法)的平均相對(duì)定位誤差4.06%；傳統(tǒng)EEMD分解與樣本熵分析相結(jié)合的方法平均相對(duì)定位誤差15.23%；IEEMD分解與傳統(tǒng)近似熵分析相結(jié)合的方法平均相對(duì)定位誤差20.17%。對(duì)比結(jié)果表明，本文提出的基于IEEMD樣本熵分析的管道泄漏定位方法與其他方法相比精度明顯更高。

4 結(jié)論

1)針對(duì)管道泄漏信號(hào)的混雜性、易受噪聲干擾性而造成泄漏定位精度低的問(wèn)題，提出基于IEEMD樣本熵分析的管道泄漏定位方法。

2)采用IEEMD與樣本熵分析結(jié)合的方法提取有效泄漏信號(hào)，有效提高了信號(hào)的信噪比，降低了信號(hào)的均方誤差。

3)基于IEEMD樣本熵分析的泄漏定位法可將定位誤差降低至4.06%，實(shí)現(xiàn)了管道泄漏精確定位。

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