楊岳鵬 徐晨 姚世聰 岳景寅 李志強(qiáng) 王海林

渤海鉆探第四鉆井工程分公司

油氣管道泄漏往往造成巨大的經(jīng)濟(jì)損失和環(huán)境污染，嚴(yán)重影響企業(yè)的社會(huì)形象，因此對(duì)油氣管道進(jìn)行泄漏定位檢測(cè)顯得尤為重要[1]。目前，管道泄漏檢測(cè)的方法主要有流量平衡法、瞬態(tài)模型法、神經(jīng)網(wǎng)絡(luò)法、負(fù)壓波法和音波法[2]，其中負(fù)壓波法和音波法的應(yīng)用最為廣泛。

1 油氣管道泄漏檢測(cè)方法及篩選

負(fù)壓波法主要依靠壓力傳感器對(duì)泄漏信號(hào)進(jìn)行捕捉，頻率范圍在60～100 Hz，泄漏中出現(xiàn)的大幅壓降在液體管道上具有較高的靈敏性，但對(duì)于可壓縮的氣體管道適應(yīng)性較差；音波法按照音波頻率不同分為次聲波（頻率范圍0.000 1～20 Hz）和超聲波（頻率范圍16 000～1012Hz）法等，BRODETSKY 等[3]通過(guò)研究音波的衰減情況，證明了頻率越大，音波信號(hào)的衰減越大，故次聲波較其余音波在介質(zhì)中傳播距離遠(yuǎn)且能量損失低。次聲波法檢測(cè)的關(guān)鍵在于信號(hào)奇異點(diǎn)的獲取和次聲波到達(dá)上、下游傳感器的時(shí)延估計(jì)，但管道在泄漏的過(guò)程中，系統(tǒng)本體噪聲、管道運(yùn)行噪聲和外部隨機(jī)噪聲等均會(huì)對(duì)次聲波信號(hào)造成影響，使信號(hào)湮滅在噪聲中，尤其對(duì)緩慢泄漏和微小泄漏的檢測(cè)成功率不高。近年來(lái)，國(guó)內(nèi)外學(xué)者針對(duì)信號(hào)降噪和時(shí)延估計(jì)進(jìn)行了大量研究，主要采用小波分析[4]、經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）[5-7]和局部均值分解（LMD）[8]等方法。但上述方法在信號(hào)處理中仍存在一定缺陷，如小波分析會(huì)造成重構(gòu)信號(hào)振動(dòng)失真，丟失有效信息；EMD 存在模態(tài)混疊和殘余噪聲問(wèn)題；LMD 存在端點(diǎn)效應(yīng)、平滑次數(shù)較多的現(xiàn)象，且平滑時(shí)步長(zhǎng)無(wú)法最優(yōu)化處理，計(jì)算速度較慢。

針對(duì)以上不足，以次聲波泄漏檢測(cè)原理為基礎(chǔ)，采用改進(jìn)的自適應(yīng)白噪聲完全集合經(jīng)驗(yàn)?zāi)B(tài)分解（CEEMDAN）對(duì)原始信號(hào)進(jìn)行降噪分解，利用相關(guān)系數(shù)選取與原始信號(hào)相關(guān)度較大的固有模態(tài)分量（IMF），去除噪聲和虛假信息，運(yùn)用二次時(shí)延估計(jì)實(shí)現(xiàn)時(shí)差計(jì)算[9]，對(duì)泄漏點(diǎn)進(jìn)行定位。通過(guò)現(xiàn)場(chǎng)試驗(yàn)對(duì)算法進(jìn)行驗(yàn)證，以期為管道泄漏定位提供理論依據(jù)和實(shí)際參考。

2 泄漏定位原理

2.1 次聲波泄漏定位原理

當(dāng)輸氣管道發(fā)生泄漏時(shí)，氣體從泄漏處噴出，管內(nèi)氣體的流動(dòng)參數(shù)發(fā)生改變，泄漏處氣體密度減小，壓力降低，泄漏點(diǎn)相鄰兩側(cè)的氣體在壓差的作用下向泄漏點(diǎn)處補(bǔ)充，進(jìn)而使泄漏點(diǎn)相鄰兩側(cè)的氣體密度進(jìn)一步減小，壓力進(jìn)一步降低，依次類推向更遠(yuǎn)處傳播。輸氣管道泄漏產(chǎn)生的次聲波本質(zhì)上屬于湍流脈動(dòng)引起的四極子聲源和偶極子聲源，聲波衰減遵循指數(shù)規(guī)律，安裝在管道上、下游兩端的聲波傳感器監(jiān)聽(tīng)并捕捉傳來(lái)的聲波信號(hào)，根據(jù)式（1）對(duì)泄漏點(diǎn)進(jìn)行定位。

式中：L1為泄漏點(diǎn)距離上游傳感器的距離，m；a為次聲波聲速，m/s；V為流體速度，m/s；L2為上、下游傳感器的距離；Δt為上、下游接收到次聲波的時(shí)間差。

一般a遠(yuǎn)大于V，故式（1）可簡(jiǎn)化為公式（2）：

可見(jiàn)，泄漏定位的關(guān)鍵在于Δt的計(jì)算精度，而信號(hào)中含有的大量噪聲會(huì)對(duì)時(shí)延估計(jì)造成較大影響，故需對(duì)信號(hào)進(jìn)行降噪處理[10-11]。

2.2 CEEMDAN信號(hào)降噪原理

CEEMDAN 是在EMD、EEMD（集合平均經(jīng)驗(yàn)?zāi)B(tài)分解）算法基礎(chǔ)上的改進(jìn)，其重構(gòu)誤差幾乎為0，解決了不同信號(hào)添加噪聲后出現(xiàn)不同模態(tài)的問(wèn)題。算法步驟如下：

（1）設(shè)采用EMD算法分解后得到的第i個(gè)固有模態(tài)分量為imfi()，在原始信號(hào)x(t)中添加了標(biāo)準(zhǔn)差為kε、均值為0 的高斯白噪聲n(t)，得到混合信號(hào)x(t)+kεn(t)，其中ε為原始信號(hào)標(biāo)準(zhǔn)差，k為白噪聲幅度。對(duì)混合信號(hào)進(jìn)行分解，按時(shí)間尺度由大到小得到M次分解時(shí)的不同模態(tài)分量，將M個(gè)模態(tài)分量求總體平均后，得到CEEMDAN算法的第一個(gè)固有模態(tài)分量IMF1(t)，如公式（3）所示：

其中，殘差r1(t)=x(t)-IMF1(t)。

（2）將殘差與白噪聲經(jīng)EMD 算法分解后的模態(tài)分量繼續(xù)組成混合信號(hào)r1(t)+imf1[kεnj(t)]，按照公式（3）分解后，得到CEEMDAN 算法的第一個(gè)固有模態(tài)分量IMF2(t)，如公式（4）所示：

其中，殘差r2(t)=r1(t)-IMF2(t)。

（3）重復(fù)以上步驟，當(dāng)殘差信號(hào)為單調(diào)函數(shù)時(shí)，停止計(jì)算，得到n+1 個(gè)固有模態(tài)分量和殘差，如公式（5）、（6）所示：

CEEMDAN算法加入了由EMD算法分解后的噪聲信號(hào)，而不是將噪聲直接添加到原始信號(hào)中，解決了白噪聲從高頻到低頻轉(zhuǎn)移傳遞的問(wèn)題。但仍存在白噪聲無(wú)法完全中和的現(xiàn)象，故采用以正負(fù)成對(duì)的形式將白噪聲添加到信號(hào)中，則每次分解得到的模態(tài)分量中，一部分殘差分量含有正噪聲，另一部分殘差分量含有負(fù)噪聲，兩者平均后正負(fù)消除，使算法更準(zhǔn)確。

原始信號(hào)分解后得到噪聲分量、冗余分量和有效分量，通過(guò)相關(guān)函數(shù)分析頻域內(nèi)兩個(gè)信號(hào)的相關(guān)程度，得到原始信號(hào)中的有效分量，相關(guān)系數(shù)公式如公式（7）所示：

式中：R(i)為第i個(gè)固有模態(tài)分量與原始信號(hào)x(t)的相關(guān)系數(shù)；IMFi(t)為固有模態(tài)分量的平均值；為原始信號(hào)的平均值。

目前，關(guān)于相關(guān)系數(shù)的閾值尚沒(méi)有統(tǒng)一認(rèn)定，但通常認(rèn)為R(i)的絕對(duì)值小于0.2 為微相關(guān)，0.2～0.5 為中相關(guān)，0.5～0.8 為高相關(guān)，0.8～1.0 為顯著相關(guān)。

2.3 二次時(shí)延估計(jì)

時(shí)延估計(jì)算法常采用廣義互相關(guān)算法，將經(jīng)相關(guān)分析后的IMF分量疊加后進(jìn)行重構(gòu)，得到降噪后的信號(hào)，將上、下游信號(hào)進(jìn)行互相關(guān)函數(shù)計(jì)算，通過(guò)相關(guān)峰值的位置得到從泄漏點(diǎn)傳播到兩端的時(shí)間差。設(shè)上、下游降噪后的信號(hào)分別為x1(t)、x2(t)，計(jì)算其互相關(guān)函數(shù)，如公式（8）所示：

式中：τ即為兩個(gè)信號(hào)的時(shí)間差。

求τ即為求R12(τ)的最大值，如公式（9）所示。定義信號(hào)x1(t)、x2(t)的互功率譜為G12，則R12(τ)可用G12表示，如公式（10）所示：

式中：φ12(w)為廣義互加權(quán)函數(shù)，當(dāng)取1時(shí)即為基本窗函數(shù)，此時(shí)R12(τ)表示基本互相關(guān)傅里葉變換；X1(w)為信號(hào)x1(t)的傅里葉變換；X2*(w)為信號(hào)x2(t)的傅里葉變換共軛。

廣義互相關(guān)分析中假設(shè)聲源信號(hào)與噪聲信號(hào)互不相關(guān)，但在實(shí)際中，噪聲信號(hào)并不完全獨(dú)立，干擾項(xiàng)依然存在，故在此采用二次時(shí)延估計(jì)相關(guān)算法。將信號(hào)x1(t)進(jìn)行自相關(guān)后得到R11(τ)，將R11(τ)與R12(τ)互相關(guān)后得到R(τ)，信號(hào)經(jīng)過(guò)一次自相關(guān)后可提高信噪比，信號(hào)噪聲的多次相關(guān)也抑制了噪聲對(duì)時(shí)差結(jié)果的影響。

3 實(shí)例分析

以兩個(gè)集氣站間管道為例進(jìn)行實(shí)驗(yàn)，驗(yàn)證算法準(zhǔn)確性。管輸介質(zhì)為凈化天然氣，管徑Φ273 mm×6 mm，全長(zhǎng)10.154 km，運(yùn)行壓力5.5～6.5 MPa，在首末站設(shè)置聲波傳感器，利用中間5.231 km處的閥室放空閥作為泄漏點(diǎn)，采用實(shí)際放空的方式模擬泄漏進(jìn)行測(cè)試，閥門開度10%。一個(gè)完整的泄漏信號(hào)波形圖中共采集104個(gè)數(shù)據(jù)點(diǎn)，信號(hào)時(shí)長(zhǎng)為20 s，次聲波信號(hào)波形圖見(jiàn)圖1，縱坐標(biāo)是歸一化后的信號(hào)幅值。當(dāng)管道處于正常工況時(shí)，聲波傳感器捕捉到的信號(hào)為背景噪聲，當(dāng)管道一旦發(fā)生泄漏，產(chǎn)生的泄漏信號(hào)和正常工況下的背景噪聲會(huì)一同被傳感器捕捉。圖1中未泄漏時(shí)信號(hào)相對(duì)平穩(wěn)，泄漏后信號(hào)幅值有較大波動(dòng)，先下降后上升，符合放空閥門啟閉的特征，同時(shí)下游信號(hào)的幅值明顯高于上游信號(hào)，這是由于泄漏點(diǎn)距離下游傳感器更近，信號(hào)衰減小，與實(shí)際情況相符。但信號(hào)波形中均含有不同程度的噪聲，需對(duì)信號(hào)降噪處理。

圖1 次聲波信號(hào)波形圖Fig.1 Waveform diagram of infrasonic signal

3.1 泄漏信號(hào)分解

按照2.2 段的流程對(duì)泄漏信號(hào)進(jìn)行自適應(yīng)，以下游次聲波信號(hào)為例，得到10 組IMF 分量和殘差分量（圖2）。其中，前3階IMF分量呈現(xiàn)寬、高頻特性，是噪聲分量的主要集中區(qū)域；4～7 階IMF 分量的幅值與原始信號(hào)保持一致，屬于有效泄漏信號(hào)分量；8～10階IMF分量的信號(hào)特征逐漸消失，呈單調(diào)函數(shù)，屬于冗余分量。為驗(yàn)證分解結(jié)果的有效性，計(jì)算各階IMF分量與原始信號(hào)的相關(guān)系數(shù)，用于提取有效信號(hào)（圖3）。前3階的相關(guān)系數(shù)均在0.2以下，屬于微相關(guān)；4-7階的相關(guān)系數(shù)較大，均在0.5 以上，屬于高相關(guān)，其中IMF6 的相關(guān)系數(shù)最大，為0.712；第8 階的相關(guān)系數(shù)突然降低至0.2 左右，9～10階的IMF分量與原始信號(hào)的相關(guān)系數(shù)幾乎為0，說(shuō)明9～10階的IMF分量與原始信號(hào)沒(méi)有相關(guān)性。綜上所述，選取4～7 階的IMF 分量和殘差分量，將其疊加進(jìn)行重構(gòu)形成降噪后的信號(hào)（圖4）。

圖2 泄漏信號(hào)分解的各階IMF分量和殘差分量（下游信號(hào)）Fig.2 IMF components and residual components of leakage signal decomposition(downstream signals)

圖3 各階IMF分量與原始信號(hào)的相關(guān)系數(shù)（下游信號(hào)）Fig.3 Correlation coefficients between IMF components of various orders and original signals(downstream signals)

圖4 原始信號(hào)與降噪信號(hào)波形圖Fig.4 Waveform diagram of original signal and denoising signal

信號(hào)降噪后，波形曲線更加光滑，保留了原始信號(hào)的波形趨勢(shì)，泄漏時(shí)刻的波形更加明顯，且其余未泄漏時(shí)刻的幅值波動(dòng)明顯減小。為評(píng)價(jià)CEEMDAN 算法的降噪效果，對(duì)比了EMD、EEMD算法分解重構(gòu)后的信號(hào)，采用信噪比(SNR)和均方根誤差(RMSE)衡量去噪效果，SNR反映有效分量與噪聲分量之間的比值，RMSE反映有效分量與原始信號(hào)之間的差異程度。因此在濾波后，應(yīng)預(yù)期提高信噪比，減小均方根誤差，即SNR越大、RMSE越小，降噪效果越好（表1）。其中，CEEMDAN算法的SNR最大且RMSE最小，說(shuō)明該算法在改善模態(tài)混疊上效果顯著，自適應(yīng)地區(qū)分了噪聲和有效信號(hào)，保證了分解的完整性，信號(hào)重構(gòu)后的誤差更小。

表1 原始信號(hào)降噪效果對(duì)比Tab.1 Comparison of noise reduction effects of original signals

3.2 泄漏定位

分別采用廣義互相關(guān)算法和二次時(shí)延互相關(guān)算法進(jìn)行定位計(jì)算，互相關(guān)時(shí)域分析見(jiàn)圖5。廣義互相關(guān)算法得到的延時(shí)采樣點(diǎn)數(shù)為281，則泄漏時(shí)差為0.562 s；二次時(shí)延互相關(guān)算法得到的延時(shí)采樣點(diǎn)數(shù)為326，則泄漏時(shí)差為0.652 s。兩種算法的泄漏時(shí)差均為正值，這是由于泄漏點(diǎn)距離上游更遠(yuǎn)，因此傳播至上游傳感器處花費(fèi)的時(shí)間更長(zhǎng)，時(shí)差量為正。次聲波波速a與管輸介質(zhì)的密度、溫度、壓力和質(zhì)點(diǎn)速度等參數(shù)相關(guān)，采用公式（11）計(jì)算如下：

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

式中：k為介質(zhì)壓縮系數(shù)，無(wú)量綱，根據(jù)管輸壓力和天然氣壓縮因子計(jì)算；ρ為氣體密度，kg/m3；D為管道直徑，m；E為管材彈性模量，Pa；e為管道壁厚，m；C為與管道約束條件相關(guān)的修正系數(shù)，按照第二種約束條件計(jì)算，C=1-μ2，μ為泊松系數(shù)取0.3。

根據(jù)不同工況多次計(jì)算取平均值，得到a=405.2 m/s。最終，按照公式（4）計(jì)算廣義互相關(guān)算法的定位結(jié)果為5.190 km，定位誤差41 m，相對(duì)誤差0.78%；二次時(shí)延互相關(guān)的定位結(jié)果為5.209 km，定位誤差22 m，相對(duì)誤差0.42%。

3.3 對(duì)比驗(yàn)證

為驗(yàn)證CEEMDAN+二次時(shí)延估計(jì)算法的可靠性，將其與文獻(xiàn)[8]中的LMD+小波分析、文獻(xiàn)[4]中的小波變換+領(lǐng)域差值兩種算法進(jìn)行對(duì)比，并選取閥門開度10%、5%和1%進(jìn)行測(cè)試（表2）。隨著泄漏點(diǎn)孔徑的減小，互相關(guān)分析的峰值點(diǎn)與真實(shí)峰值點(diǎn)之間的差距逐漸增大，定位誤差也逐漸增大，其中當(dāng)閥門開度為1%時(shí)，只有CEEMDAN+二次時(shí)延估計(jì)算法對(duì)泄漏進(jìn)行了定位，而其余兩種算法由于信號(hào)下降的拐點(diǎn)不易識(shí)別，故未識(shí)別出報(bào)警事件，說(shuō)明CEEMDAN+二次時(shí)延估計(jì)算法未出現(xiàn)誤報(bào)或漏報(bào)，對(duì)微小泄漏的檢測(cè)具有良好的適應(yīng)性。當(dāng)閥門開度為10%時(shí)，CEEMDAN+二次時(shí)延估計(jì)算法的定位相對(duì)誤差較其余兩種算法分別降低了71.03%、60%；當(dāng)閥門開度為5%時(shí)，定位相對(duì)誤差較其余兩種算法分別降低了84.21%、73.38%，對(duì)10 km 左右的管線總體定位誤差可控制在50 m以內(nèi)。

表2 不同定位算法結(jié)果對(duì)比Tab.2 Comparison of results of different positioning methods

4 結(jié)論

（1）采用CEEMDAN算法對(duì)次聲波信號(hào)進(jìn)行了降噪處理，處理后SNR最大，且RMSE最小，說(shuō)明該算法在改善模態(tài)混疊上效果顯著，信號(hào)重構(gòu)后的誤差更小。

（2）采用二次時(shí)延互相關(guān)算法進(jìn)行定位計(jì)算，并與其余算法進(jìn)行了對(duì)比，CEEMDAN+二次時(shí)延估計(jì)算法在微小泄漏定位檢測(cè)方面具有良好的適應(yīng)性，最大相對(duì)誤差0.84%。