基于CEEMDAN+二次時延估計的管道泄漏定位技術

楊岳鵬 徐晨 姚世聰 岳景寅 李志強 王海林

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

油氣管道泄漏往往造成巨大的經濟損失和環境污染，嚴重影響企業的社會形象，因此對油氣管道進行泄漏定位檢測顯得尤為重要[1]。目前，管道泄漏檢測的方法主要有流量平衡法、瞬態模型法、神經網絡法、負壓波法和音波法[2]，其中負壓波法和音波法的應用最為廣泛。

1 油氣管道泄漏檢測方法及篩選

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

針對以上不足，以次聲波泄漏檢測原理為基礎，采用改進的自適應白噪聲完全集合經驗模態分解（CEEMDAN）對原始信號進行降噪分解，利用相關系數選取與原始信號相關度較大的固有模態分量（IMF），去除噪聲和虛假信息，運用二次時延估計實現時差計算[9]，對泄漏點進行定位。通過現場試驗對算法進行驗證，以期為管道泄漏定位提供理論依據和實際參考。

2 泄漏定位原理

2.1 次聲波泄漏定位原理

當輸氣管道發生泄漏時，氣體從泄漏處噴出，管內氣體的流動參數發生改變，泄漏處氣體密度減小，壓力降低，泄漏點相鄰兩側的氣體在壓差的作用下向泄漏點處補充，進而使泄漏點相鄰兩側的氣體密度進一步減小，壓力進一步降低，依次類推向更遠處傳播。輸氣管道泄漏產生的次聲波本質上屬于湍流脈動引起的四極子聲源和偶極子聲源，聲波衰減遵循指數規律，安裝在管道上、下游兩端的聲波傳感器監聽并捕捉傳來的聲波信號，根據式（1）對泄漏點進行定位。

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

一般a遠大于V，故式（1）可簡化為公式（2）：

可見，泄漏定位的關鍵在于Δt的計算精度，而信號中含有的大量噪聲會對時延估計造成較大影響，故需對信號進行降噪處理[10-11]。

2.2 CEEMDAN信號降噪原理

CEEMDAN 是在EMD、EEMD（集合平均經驗模態分解）算法基礎上的改進，其重構誤差幾乎為0，解決了不同信號添加噪聲后出現不同模態的問題。算法步驟如下：

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

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

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

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

（3）重復以上步驟，當殘差信號為單調函數時，停止計算，得到n+1 個固有模態分量和殘差，如公式（5）、（6）所示：

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

原始信號分解后得到噪聲分量、冗余分量和有效分量，通過相關函數分析頻域內兩個信號的相關程度，得到原始信號中的有效分量，相關系數公式如公式（7）所示：

式中：R(i)為第i個固有模態分量與原始信號x(t)的相關系數；IMFi(t)為固有模態分量的平均值；為原始信號的平均值。

目前，關于相關系數的閾值尚沒有統一認定，但通常認為R(i)的絕對值小于0.2 為微相關，0.2～0.5 為中相關，0.5～0.8 為高相關，0.8～1.0 為顯著相關。

2.3 二次時延估計

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

式中：τ即為兩個信號的時間差。

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

式中：φ12(w)為廣義互加權函數，當取1時即為基本窗函數，此時R12(τ)表示基本互相關傅里葉變換；X1(w)為信號x1(t)的傅里葉變換；X2*(w)為信號x2(t)的傅里葉變換共軛。

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

3 實例分析

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

圖1 次聲波信號波形圖Fig.1 Waveform diagram of infrasonic signal

3.1 泄漏信號分解

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

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

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

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

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

表1 原始信號降噪效果對比Tab.1 Comparison of noise reduction effects of original signals

3.2 泄漏定位

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

圖5 互相關時域分析Fig.5 Cross-correlation time domain analysis

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

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

3.3 對比驗證

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

表2 不同定位算法結果對比Tab.2 Comparison of results of different positioning methods

4 結論

（1）采用CEEMDAN算法對次聲波信號進行了降噪處理，處理后SNR最大，且RMSE最小，說明該算法在改善模態混疊上效果顯著，信號重構后的誤差更小。

（2）采用二次時延互相關算法進行定位計算，并與其余算法進行了對比，CEEMDAN+二次時延估計算法在微小泄漏定位檢測方面具有良好的適應性，最大相對誤差0.84%。