基于時頻分析和改進(jìn)差分進(jìn)化的多點泄漏定位方法

郎憲明， 朱永強(qiáng)， 李金娜， 曹江濤， 李 平， 宋華東

(1.遼寧石油化工大學(xué) 信息與控制工程學(xué)院，遼寧 撫順 113001； 2.東北大學(xué) 信息科學(xué)學(xué)院，沈陽 110819； 3.沈陽儀表科學(xué)研究院有限公司，沈陽 110043)

管道是輸送石油、天然氣等流體的重要工具[1]。目前，有多種管道泄漏檢測方法已被應(yīng)用于管道泄漏監(jiān)測。然而，現(xiàn)有的泄漏定位方法一般都是針對某一個泄漏點進(jìn)行的，當(dāng)流體管道局部管段出現(xiàn)壓力過高、腐蝕、磨損、振動等，有時會出現(xiàn)2個或2個以上泄漏點。因此，針對單點泄漏的檢測方法對管道進(jìn)行多點泄漏定位往往失效[2]。

目前，基于聲波的管泄漏檢測方法主要是基于單點泄漏聲波信號的智能分析[3-4]。利用Hilbert-Huang變換對管道首末站的聲發(fā)射信號進(jìn)行去噪處理，精確獲得聲波信號突變點時間，提高了泄漏定位的精度[5]。針對壓力信號下降緩慢的缺點，提出了一種信號壓縮方法，將壓力信號轉(zhuǎn)換為超聲波聲速信號，突出聲波拐點的信息[6]。采用低通濾波與小波濾波相結(jié)合的方法，對管道首末站的壓力信號進(jìn)行實時高精度濾波，保證泄漏信息不被過濾掉[7]。由于小波去噪需要根據(jù)經(jīng)驗選擇合適的小波基函數(shù)，而經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition，EMD)不用憑經(jīng)驗選擇基函數(shù)的優(yōu)點，因此采用EMD對管道泄漏信號進(jìn)行處理，提取管道泄漏信號的突變點信息，以減小外界環(huán)境噪聲對泄漏點定位的影響[8]。變分模態(tài)分解(variational mode decomposition，VMD)是一種新的信號分解方法[9]，它可以降低EMD模態(tài)混疊對信號消噪的影響。

由于多數(shù)研究是關(guān)于管道的單點泄漏情況[10-15]，當(dāng)管道同時發(fā)生多點泄漏，產(chǎn)生的瞬態(tài)泄漏聲波信號混疊在一起，由于泄漏聲波信號相互影響，以及外部噪聲(如施工、管道閥門噪聲等)疊加在管道泄漏聲波信號上，從而影響瞬態(tài)聲波信號在信號傳播中的衰減規(guī)律，因此在混雜信號的情況下，多點泄漏信號不能有效地去噪。文獻(xiàn)[16]采用擴(kuò)展卡爾曼濾波方法對多點連續(xù)泄漏進(jìn)行檢測。文獻(xiàn)[17]利用超聲波聲速的變化來進(jìn)行多點連續(xù)泄漏進(jìn)行定位。

當(dāng)多點泄漏聲波向管道首末站傳播時，很難區(qū)分各個泄漏聲波信號到達(dá)管道首末站的突變點時間。因此，基于單點泄漏的泄漏定位方法不適用于多點泄漏。然而，時頻分析(time-frequency analysis，TFA)可以從時域和頻域分析來分析非平穩(wěn)的多點泄漏聲波信號[18]，通過時頻分析來確定多點泄漏的數(shù)量。

差分進(jìn)化(differential evolution，DE)算法是一種新的群體智能優(yōu)化算法[19]。它具有簡單、易于實現(xiàn)、無梯度信息、參數(shù)少等優(yōu)點，被廣泛應(yīng)用于工業(yè)過程參數(shù)優(yōu)化中。然而，進(jìn)化算法在進(jìn)化后期容易陷入局部最優(yōu)，收斂速度較慢。針對此問題，文獻(xiàn)[20]引入了混合自適應(yīng)局部搜索算法來提高DE算法的收斂速度。為了加快離散進(jìn)化算法的收斂速度，文獻(xiàn)[21]采用正交交替的方法來提高收斂速度。為了提高差分算法的收斂速度和精度，本文采用粒子群優(yōu)化算法對差分進(jìn)化算法進(jìn)行改進(jìn)，根據(jù)建立的多點泄漏定位目標(biāo)函數(shù)，采用改進(jìn)差分進(jìn)化算法估計多點泄漏位置。

1 多點泄漏定位原理

1.1 改進(jìn)變分模態(tài)分解方法

當(dāng)VMD用于信號分解時，根據(jù)經(jīng)驗選擇適當(dāng)?shù)哪B(tài)分解個數(shù)m[22]。由于m決定了頻率分辨率，如果選擇m過小，則會出現(xiàn)信號欠分解，對泄漏聲波信號變化不敏感，如果選擇m太大，會產(chǎn)生偽分量，增加計算量。

將采集到的聲波信號y(t)采用VMD分解，計算如下

(1)

式中，γ(t)是信號的重構(gòu)誤差。

能量函數(shù)的計算如下

Eγ=γ(t)2

(2)

Eγ能量越小，兩個信號的相似度越高，如果誤差能量為0，則說明兩個信號相同。

由于VMD分解模態(tài)的數(shù)量對信號重構(gòu)有影響，因此采用能量函數(shù)對VMD進(jìn)行改進(jìn)。改進(jìn)變分模態(tài)分解(improved variational mode decomposing, IVMD)算法的計算步驟如下：

(1) VMD算法中的分解模態(tài)數(shù)量設(shè)定為m。

(2) 經(jīng)VMD信號分解后，獲得了一定數(shù)量的本征模態(tài)函數(shù)(intrinsic mode functions, IMFs)。

(3) 計算不同數(shù)量分解本征模態(tài)函數(shù)的誤差能量函數(shù)。

(4) 繪制誤差能量圖，首先計算出各誤差能量函數(shù)的平均值閾值。一般情況下，如果誤差能量值遠(yuǎn)小于閾值，稱為有效本征模態(tài)函數(shù)。有效本征模態(tài)函數(shù)包含泄漏聲波信號的突變點時刻、信號幅值變化等信息。

因此，泄漏聲波信號經(jīng)IVMD選擇有效本征模態(tài)函數(shù)對信號進(jìn)行消噪。

1.2 多點泄漏定位算法

當(dāng)管道出現(xiàn)多點泄漏同時發(fā)生時，多點泄漏引起的聲波幅值變化混合在一起[23-24]，由于多點泄漏信號之間的相互作用，在管道首末站獲得多點泄漏聲波的疊加信號，多點泄漏聲波傳播和泄漏衰減的振幅描述為

(3)

(4)

式中：xi是i個泄漏點產(chǎn)生的泄漏聲波信號原始振幅;a是衰減因子;xa和xb分別是管道首末站獲取的多點泄漏聲波信號;li是第i個泄漏點;L是管道的長度。

當(dāng)出現(xiàn)單點泄漏時，根據(jù)式(3)和式(4)，按照首末站獲得的泄漏聲波信號幅值變化可計算出單點泄漏位置。

根據(jù)多點泄漏聲波信號幅值變化的特點，假設(shè)每個泄漏點產(chǎn)生的聲波信號原始振幅大小一樣，根據(jù)式(3)和式(4)，則多點泄漏的定位計算如下

(5)

根據(jù)管道首末站的獲得的多點泄漏聲波信號振幅變化，根據(jù)式(5)可以計算出多點泄漏位置。然而，同時發(fā)生的多點泄漏引起的振幅變化是重疊的，無法區(qū)分多點泄漏引起的振幅變化。因此，根據(jù)已知多點泄漏數(shù)量的情況下，根據(jù)式(5)建立多點泄漏定位函數(shù)，然后通過最小化管道首末站的多點泄漏聲波信號振幅誤差來估計多點泄漏位置。

2 多點泄漏定位方法

2.1 時頻分析

采用Wigner分布對管道多點泄漏聲波信號振幅的變化進(jìn)行分析[25-27]。

Sx(t,ω)=ASPW(x)

(6)

式中：x是多點泄漏信號的振幅變化;Sx(t,ω)為多點泄漏聲波振幅變化的時頻分布。

Wigner的分布描述如下

ASPW(t,ω)=

(7)

式中，g(t)和h(t)分別是多點泄漏聲波信號時域和頻域的光滑函數(shù)。

多點泄漏信號的振幅變化時頻分布具有不同的分布曲線，可以根據(jù)分布曲線確定泄漏點的數(shù)量。

2.2 差分進(jìn)化算法

DE算法是一種隨機(jī)并行的搜索算法[28]。它是從隨機(jī)的初始個體開始，按照一定的迭代規(guī)則。根據(jù)每個個體的適應(yīng)度值，保留較好的個體，引導(dǎo)搜索過程達(dá)到最優(yōu)解。

對于具有n個變量的全局優(yōu)化問題。全局優(yōu)化問題可轉(zhuǎn)化為最小值問題，求解如下函數(shù)

(8)

式中：D是問題空間解的維數(shù)；aj和bj分別是χj的上限和下限。

DE算法主要包括變異操作、交叉操作和選擇操作。

(1) 突變操作

(9)

式中：t是進(jìn)化次數(shù)；χr1,χr2和χr3是隨機(jī)選取的三個單獨(dú)個體；ui是修正因子；F是突變因子。

(2) 交叉操作

為了增加種群多樣性，采用了交叉。

(10)

式中：j=1,2,…,D，D是維度；RC是交叉率，RC∈[0,1]。

(3) 選擇操作

由式(10)從試驗向量和原始向量中選出較好的個體。

(11)

式中，f(χ)是χ的目標(biāo)函數(shù)。

2.3 改進(jìn)差分進(jìn)化(improved differential evolution，IDE)算法

在優(yōu)化過程中，DE存在收斂速度慢、早熟等缺點，嚴(yán)重影響了算法的性能。因此，采用了一種改進(jìn)的差分進(jìn)化算法。在DE算法中引入PSO原理，用PSO搜索整個種群，使其快速跳出局部最優(yōu)，避免早熟現(xiàn)象[29]。

粒子群優(yōu)化的個體變異操作如下

(12)

其中，χbest是優(yōu)勢個體。

IDE的流程圖如圖1所示。

圖1 IDE流程圖Fig.1 Flow chart of IDE

2.4 多點泄漏定位方法

多點泄漏定位方法的流程圖如圖2所示。從圖2中，可以使用IVMD對泄漏聲波信號進(jìn)行去噪。傳統(tǒng)的泄漏定位方法很難區(qū)分多點泄漏個數(shù)，多點泄漏聲波信號經(jīng)IVMD去噪后，然后采用TFA算法來估計泄漏點數(shù)量。然后，根據(jù)式(5)建立的多點泄漏定位目標(biāo)函數(shù)，采用IDE估計多點泄漏的位置。

圖2 多點泄漏定位的流程圖Fig.2 Flow chart of multiple leaks localization method

3 現(xiàn)場試驗

在中國石油某個成品油管線進(jìn)行了多次多點泄漏試驗。管道長度為20 000 m，管道內(nèi)徑為100 mm。管道上游壓力為2 MPa，管道下游壓力為0.5 MPa，流量為200 L/min，密度是830 kg/m3，環(huán)境溫度是25 ℃，通過聲傳感器(PCB 106B)獲取管道首末站多點泄漏的聲波信號。聲波數(shù)據(jù)的采集試驗裝置原理圖如圖3所示。管道試驗系統(tǒng)由成品油管道和多點泄漏定位系統(tǒng)組成。多點泄漏定位系統(tǒng)主要由PCB106b、工控機(jī)和數(shù)據(jù)采集器NI CDAQ-9181組成。泄漏定位方法的關(guān)鍵是采集管道首末站多點泄漏聲波信號的時間同步，本試驗時間同步采用GPS衛(wèi)星時鐘信號校時，管道首末站的泄漏聲波信號采樣頻率是1 000 Hz。為了模擬多點泄漏，3個泄漏位置分別距離管道首站14 200 m、15 500 m和15 800 m處，根據(jù)式(5)，模擬3點泄漏的泄漏閥門孔徑分別為15 mm、15 mm和15 mm。實際管道運(yùn)行過程中，超過3點同時泄漏會給管道運(yùn)行造成非常嚴(yán)重的事故。另外，3點同時泄漏是目前研究熱點問題，本文只針對管道3點同時泄漏情況進(jìn)行相關(guān)試驗。

圖3 試驗裝置示意圖Fig.3 Schematic diagram of test apparatus

在距離首站14 200 m位置，模擬泄漏孔徑為15 mm，泄漏持續(xù)時間為2 s，管道首末站分別獲得的單點泄漏的聲波信號變化如圖4(a)所示。在距離首站14 200 m和15 500 m的2個泄漏位置，模擬泄漏孔徑分別是15 mm和15 mm，泄漏持續(xù)時間為2 s，管道首末站分別獲得的2點泄漏的聲波信號變化如圖4(a)所示。另外，在距離首站14 200 m、15 500 m和15 800 m的3個泄漏位置，模擬泄漏孔徑分別是15 mm、15 mm和15 mm，泄漏持續(xù)時間為2 s，管道首末站分別獲得的3點泄漏的聲波信號變化如圖4(c)所示。

(a) 單點泄漏

(b) 2點泄漏

(c) 3點泄漏圖4 獲取的多點泄漏聲波信號Fig.4 Acquired acoustic wave signal under multiple leaks

從圖4可以看出，隨著泄漏點數(shù)量的增加，聲壓信號幅值逐漸增大。管道末站獲得單點泄漏的聲壓信號幅值波動范圍為0～2.2，2點泄漏時的聲壓信號幅值波動范圍為0～6.2，3點泄漏的聲壓信號幅值波動范圍為0～7.8。

根據(jù)PCB106b傳感器獲得的泄漏聲波信號，經(jīng)VMD進(jìn)行5層本征模態(tài)分解，5個IMFs可以較好地消除被測信號的噪聲。根據(jù)聲波信號的5個IMFs，計算誤差能量函數(shù)并繪制誤差能量圖，本征模態(tài)函數(shù)的誤差能量圖如圖5所示。

從圖5可以看出，當(dāng)本征模態(tài)的數(shù)目為2個時，誤差能量函數(shù)值為0.92，這些本征模式包含更有效的泄漏聲波信號突變點、聲波信號幅值變化等信息。因此，我們可以設(shè)定閾值為0.93，并且在本試驗中選擇了2個本征模態(tài)函數(shù)。泄漏聲波信號經(jīng)VMD的5層本征模態(tài)分解和經(jīng)IVMD的2個本征模態(tài)函數(shù)的聲波信號消噪對比曲線如圖6所示。

圖5 本征模態(tài)函數(shù)的誤差能量函數(shù)值Fig.5 Error energy diagram of intrinsic mode

圖6 經(jīng)VMD和IVMD消噪后的聲波信號對比曲線Fig.6 Comparison curve of acoustic signals de-noised with VMD and IVMD

從圖6中可以看出，IVMD的2個IMFs可以較好地恢復(fù)泄漏聲波信號的變化，并且相比VMD的計算量要小。

根據(jù)誤差能量函數(shù)計算，IVMD選擇2個IMFs對多點泄漏聲波信號進(jìn)行信號消噪處理。通過IVMD對單點泄漏、2點泄漏、3點泄漏三種情況下的泄漏聲波信號進(jìn)行處理。根據(jù)管道首末站傳感器PCB106b獲得的泄漏聲波信號進(jìn)行IVMD信號消噪處理，經(jīng)IVMD和EMD對單點泄漏、2點泄漏、3點泄漏聲波信號的去噪對比曲線如圖7所示。

(a) 單點泄漏

(b) 2點泄漏

(c) 3點泄漏圖7 經(jīng)IVMD和EMD消噪后的多點泄漏聲波信號Fig.7 Multiple leaks signals de-noised with IVMD and EMD

通過對單點泄漏、2點泄漏、3點泄漏信號消噪后進(jìn)行TFA分析，其時頻分析曲線如圖8所示。

從圖8可以看出，單點泄漏、2點泄漏、3點泄漏泄漏引起的頻率分量變化明顯，可以通過TFA有效地區(qū)分。因此，可以通過對多點泄漏聲波信號的TFA分析，估計出泄漏點的數(shù)量。

通過TFA方法估計泄漏點的數(shù)量，根據(jù)式(5)建立多點泄漏定位目標(biāo)函數(shù)，采用IDE對多點泄漏定位目標(biāo)函數(shù)進(jìn)行求解，估計出多點泄漏的位置。在DE算法中，突變因子F=0.5、迭代次數(shù)t=100和維數(shù)D，其中維數(shù)D由TFA估計出泄漏點的數(shù)量決定，交叉率RC=0.9和種群數(shù)量為100。另外，將IDE的參數(shù)與DE的數(shù)值一樣，粒子群數(shù)量是40，維數(shù)是20，學(xué)習(xí)因子是2，慣性權(quán)重是1.2。

當(dāng)管道同時發(fā)生3點泄漏時，IDE和DE的適應(yīng)度函數(shù)收斂性能曲線如圖9所示。

圖9 IDE和DE的適應(yīng)度函數(shù)收斂曲線Fig.9 Fitness function convergence performance of IDE and DE

從圖9可以看出，當(dāng)IDE算法的迭代次數(shù)為75時，適應(yīng)度函數(shù)最終收斂到0.35左右，而DE算法的適應(yīng)度函數(shù)為72，但適應(yīng)度函數(shù)最終收斂到6，因此，IDE算法在估計泄漏點位置求解速度和求解精度方面相比DE算法更具優(yōu)勢。由于IDE算法通過引入PSO算法對其進(jìn)行了改進(jìn)，使得IDE算法在求解多點泄漏定位函數(shù)時跳出局部極值，收斂到全局最優(yōu)解。因此，IDE算法求解多點泄漏定位函數(shù)。

管道首末站采用PCB106b傳感器獲得多點泄漏聲波信號數(shù)據(jù)，信號采樣率為1 000 Hz。建立管道單點泄漏、2點泄漏、3點泄漏共3種工況數(shù)據(jù)，每種工況進(jìn)行30組試驗。采用IVMD方法對單點泄漏、2點泄漏、3點泄漏的泄漏聲波信號進(jìn)行降噪處理。利用2個本征模態(tài)函數(shù)分別對單點泄漏、2點泄漏、3點泄漏的泄漏聲波信號進(jìn)行重構(gòu)。同時，EMD算法選取6～12個IMF分別對單點泄漏、2點泄漏、3點泄漏的泄漏聲波信號進(jìn)行消噪。然后，利用TFA和IDE對采集到的多次泄漏聲數(shù)據(jù)進(jìn)行處理。不同方法的泄漏試驗定位結(jié)果如表1所示。由于多點泄漏是一個小概率事件，并且超過3點同時泄漏會引起流體管道出現(xiàn)重大事故發(fā)生。目前，3點泄漏是研究的熱點問題，在現(xiàn)場試驗中，不考慮管道出現(xiàn)超過3點同時泄漏的情況。

從表1可以看出，多點泄漏的互相關(guān)分析中只有一個峰值。因此，無法計算2點或3點泄漏的泄漏位置。TFA可以用來確定多個泄漏點的數(shù)量，然后使用IDE來估計多點泄漏的位置。在距離管道首站14 200 m，15 500 m和15 800 m位置同時發(fā)生3點泄漏，經(jīng)過IVMD對泄漏聲波信號進(jìn)行消噪，然后將消噪后的信號進(jìn)行TFA計算，估計出泄漏點數(shù)量，并根據(jù)式(5)，建立多點泄漏定位目標(biāo)函數(shù)，采用IDE算法求解此目標(biāo)函數(shù)。經(jīng)TFA和IDE算法估計的3點泄漏位置如圖10所示。

表1 不同模型的多點泄漏定位結(jié)果對比

圖10 經(jīng)TFA和IDE估計的3點泄漏位置Fig.10 Three leaks localization by TFA and IDE

從圖10可以看出，用IDE算法迭代次數(shù)為14，估計的第1點泄漏定位值為14 230 m，定位誤差為30 m；第2點泄漏的定位估計值為15 530 m，IDE算法迭代次數(shù)為10次，定位誤差為30 m；第3點漏定位估計值為15 818 m，IDE算法迭代次數(shù)為16次，定位誤差為18 m。因此，該方法可以估計多點泄漏的位置。由于管道實際運(yùn)行中，管道首末站的次聲波傳感器位置及模擬多點泄漏位置不能隨意調(diào)整。因此，在接下來的研究中會根據(jù)實際工況調(diào)整傳感器及泄漏點的位置，體現(xiàn)多點泄漏定位方法的實用性。

4 結(jié) 論

本文提出一種基于時頻分析和改進(jìn)的差分進(jìn)化算法的多點泄漏定位方法。通過IVMD對泄漏聲波信號進(jìn)行消噪處理，將消噪信號進(jìn)行TFA分析，估計多點泄漏的泄漏點數(shù)量，根據(jù)估計的泄漏點數(shù)量及泄漏定位目標(biāo)函數(shù)，采用IDE算法估計多點泄漏位置。在某成品油管線進(jìn)行了多點泄漏定位試驗，泄漏點定位誤差最小是18 m。試驗結(jié)果表明，該方法可以在3點泄漏同時發(fā)生的情況下，有效地估計3點泄漏的位置。在后續(xù)的研究中，課題組將對管道3點以上同時泄漏點的定位進(jìn)行研究。