張 超, 侯 男, 路敬祎, 王 闖

(東北石油大學(xué) a. 人工智能能源研究院; b. 黑龍江省網(wǎng)絡(luò)化與智能控制重點實驗室, 黑龍江 大慶 163318)

0 引 言

管道運輸作為能源存儲與使用的常用手段, 具有運輸成本低、 適用性強等優(yōu)勢, 但管道存在設(shè)備老化以及人為破壞等多種外部因素, 管道泄漏時有發(fā)生, 甚至對人們的生命財產(chǎn)安全造成傷害[1]。為了避免或降低管道泄漏事故的發(fā)生, 對管道泄漏檢測技術(shù)的研究非常必要。目前常用的去噪方法主要有傅里葉變換、 小波變換等, 但這些方法不能依據(jù)信號自身的特點自行確定對應(yīng)的分解條件, 在先驗經(jīng)驗的影響下會產(chǎn)生誤差而導(dǎo)致精度下降。Huang等[2]提出了一種稱為經(jīng)驗?zāi)B(tài)分解(EMD: Empirical Mode Decomposition)的自適應(yīng)信號處理技術(shù), 該方法能較好地處理非線性和非平穩(wěn)信號。但是, EMD存在端點效應(yīng)和模態(tài)混疊的弊端, 使信號分解效果不理想。Smith[3]在經(jīng)驗?zāi)B(tài)分解方法的基礎(chǔ)上提出了一種新型自適應(yīng)的時頻分析方法----局域均值分解(LMD: Local Mean Decomposition)。LMD雖然可以對信號進行濾波, 但是解調(diào)過程會導(dǎo)致信號出現(xiàn)突變。

Dragomiretskiy等[4]提出VMD(Variational Mode Decomposition)算法, 是一種自適應(yīng)信號處理方法, 可有效克服EMD算法的端點效應(yīng)和模態(tài)混疊問題。預(yù)設(shè)尺度k和懲罰因子α對分解效果有很大的影響, 而人為經(jīng)驗設(shè)定的k過大或過小通常會導(dǎo)致過分解或欠分解現(xiàn)象;α越大模態(tài)分量帶寬越小; 反之,α越小模態(tài)分量的帶寬就越大。由此可知VMD參數(shù)的設(shè)置決定是否能準(zhǔn)確提取出信號的有用信息。為減少人為設(shè)置參數(shù)可能導(dǎo)致的不良影響, 蔣麗英等[5]采用粒子群優(yōu)化算法對VMD算法的最佳影響參數(shù)組合[k,α]進行搜索, 利用參數(shù)優(yōu)化后的VMD算法分解采集的齒輪數(shù)據(jù), 并應(yīng)用到齒輪故障特征參數(shù)提取。

筆者采用粒子群優(yōu)化算法對VMD參數(shù)組進行尋優(yōu)處理, 而粒子群算法在后期容易收斂緩慢, 陷入局部最優(yōu), 無法得到較為準(zhǔn)確的結(jié)果。針對上述問題, 筆者提出了一種利用混沌算法、 Sigmoid函數(shù)改進PSO(Particle Swarm Optimization)優(yōu)化VMD的改進算法, 能得到k和α的最優(yōu)解, 然后利用VMD對仿真信號進行分解, 分別計算仿真信號的概率密度函數(shù)(PDF: Probability Density Function)和每個模態(tài)分量的PDF之間的歐氏距離, 根據(jù)相似性大小選擇主要模態(tài)分量進行重構(gòu), 從而消除噪聲的干擾。通過算法仿真和管道數(shù)據(jù)實驗分析, 結(jié)果表明改進的PSO-VMD算法具有較好的去噪性能, 并驗證了該方法在管道泄漏檢測中應(yīng)用的可行性。

1 變分模態(tài)分解基本原理

VMD算法將輸入的信號進行模態(tài)分離處理后, 每個IMF(Intrinsic Mode Function)分量的帶寬和頻率中心由連續(xù)迭代確定, 由k個IMF分量組成集合uk={u1,u2,…,uk},k=1,2,…,k, 獲得的每個IMF分量具體步驟如下：

1) 利用希爾伯特變換, 計算每個uk的解析函數(shù), 得到與其對應(yīng)的邊際譜;

2) 通過指數(shù)混合調(diào)制到估計的中心頻率, 將每個模態(tài)uk的頻譜轉(zhuǎn)移到對應(yīng)的基帶;

3) 由解調(diào)信號的高斯光滑度和梯度平方準(zhǔn)則方法估計帶寬大小。

上述步驟獲得的約束變分問題為

(1)

(2)

其中*為卷積;f為分解的主信號;ωk={ω1,ω2,…,ωk}為其各自對應(yīng)的中心頻率; ?t為對函數(shù)求時間的導(dǎo)數(shù);δ(t)為單位脈沖函數(shù)。

在求解約束問題時, 需要考慮懲罰因子α和用λ(T)表示的拉格朗日乘子這兩項。將引入的兩項結(jié)合得到擴展的拉格朗日表達式如下

通常采用乘法算子交替方向解決上述問題, 對uk+1,ωk+1,λk+1不斷迭代優(yōu)化求取擴展拉格朗日表達式的“鞍點”。

對信號進行VMD分解的具體步驟如下:

2) 執(zhí)行循環(huán), 令n=n+1, 更新uk和ωk

(4)

(5)

3) 更新乘法算子λ

(6)

2 粒子群優(yōu)化算法

粒子群優(yōu)化算法是在人類對鳥類捕食行為研究的基礎(chǔ)上得到的一種進化計算方法[6]。粒子群優(yōu)化算法作為一種啟發(fā)式智能算法, 整個群體由d維搜索空間中移動的N個粒子組成[7], 第i個粒子在空間中的位置表示為Xid=(Xi1,Xi2,…,Xid), 第i個粒子歷史最優(yōu)位置表示為Pid=(Pi1,Pi2,…,Pid), 全部粒子群歷史最優(yōu)位置表示為Gid=(Gi1,Gi2,…,Gid), 第i個粒子的歷史速度表示為Vid=(Vi1,Vi2,…,Vid)。

在每次迭代過程中, 粒子通過個體極值和群體極值更新自身的速度和位置, 更新公式如下

Vid(t+1)=ωVid(t)+c1r1(Pid(t)-Xid(t))+c2r2(Gid(t)-Xid(t))

(7)

Xid(t+1)=Xid(t)+Vid(t+1)

(8)

其中ω為慣性系數(shù),c1和c2為加速系數(shù), 代表認知系數(shù)和社會系數(shù)。r1和r2為[0,1]之間服從均勻分布的兩個獨立隨機數(shù)。由式(7)和式(8)可以看到, 參數(shù)ω、c1、c2的取值對粒子的更新速度和位置有很大的影響。文獻[8]表明當(dāng)慣性權(quán)值較大時, 粒子搜索能力增強, 局部細調(diào)能力減弱; 反之, 當(dāng)慣性權(quán)值較小時,粒子的局部開發(fā)能力增強, 探測能力減弱。適當(dāng)調(diào)整加速因子不僅便于粒子實現(xiàn)快速搜索, 提高全局搜索能力, 同時利于局部的精細化, 獲得更好的全局最優(yōu)解[9]。由此可見, 選取恰當(dāng)?shù)膮?shù)是研究和改進粒子群優(yōu)化算法性能的關(guān)鍵。

3 改進的PSO算法

3.1 混沌算法

混沌是一種看似比較凌亂, 但實際卻是非常有規(guī)則的運動[10], 具有遍歷性、 隨機性和規(guī)律性等特性?；煦缢惴ǖ闹饕枷胧窍葘⒎N群進行混沌優(yōu)化, 再將其映射到原解所在的空間范圍內(nèi)?；煦鐑?yōu)化后使結(jié)果跳出局部最優(yōu), 并且收斂速度加快。筆者主要使用混沌算法對粒子群的加速因子進行初始優(yōu)化, 使其跳出局部最優(yōu), 加快尋優(yōu)?；煦缬泻芏鄤恿W(xué)模型, 筆者采用一維非線性映射模型Logistic映射[11]

Y(t+1)=μY(t)(1-Y(t))

(9)

其中Y(t)是自變量,μ是控制變量, 當(dāng)μ=4時, 系統(tǒng)處于完全混沌狀態(tài), 式(9)產(chǎn)生的序列為混沌變量。

3.2 Sigmoid函數(shù)

Sigmoid函數(shù)常被用于神經(jīng)網(wǎng)絡(luò)的激活函數(shù), 將變量映射到(0,1)之間

(10)

Sigmoid函數(shù)擁有平滑的上下邊界域, 當(dāng)x<-9.903 438時,f(x)=0; 當(dāng)x≥9.903 438時,f(x)=1。

3.3 優(yōu)化PSO參數(shù)

對粒子群中的加速因子進行混沌優(yōu)化, 映射到原解空間為

c1s=(0.5-0.9)(i-j)/i+0.9

(11)

其中i=20為最大迭代次數(shù),j為當(dāng)前迭代次數(shù),c1s為加速因子終值。

式(11)是在區(qū)間[0.5,0.9]上變化的單調(diào)遞增函數(shù), 代入到混沌映射模型式(9)中變成一個在[1,0.36]上的遞減函數(shù), 經(jīng)過反復(fù)仿真實驗可知, 式(11)代入混沌映射模型后乘上常數(shù)θ=2得到

c1=4.0θc1s(1-c1s)

(12)

其中常數(shù)θ=2, 使c1數(shù)值由2遞減至0.72時, 算法可以獲得最優(yōu)的適應(yīng)度值

c2s=(0.9-0.5)(i-j)/i+0.5

(13)

c2=4.0θc2s(1-c2s)

(14)

其中c2s為加速因子終值,c2原理同c1。這樣的改進可使粒子在迭代初期,c1值較大,c2值較小, 有利于粒子增強全局搜索能力; 隨著迭代次數(shù)的增加,c1值變得較小,c2值變得較大, 這樣有利于局部的細節(jié)化, 使粒子局部搜索能力增強。因為在收斂后期個體極值影響變小, 全局極值影響加大, 所以這兩個函數(shù)能更加細化這一影響。

采用Sigmoid函數(shù)優(yōu)化慣性權(quán)重, 即

(15)

由式(15)可見, 在迭代初期,ω值為1, 隨著迭代次數(shù)的不斷增加,ω逐漸減小, 最后減小到0.5。當(dāng)慣性權(quán)值較大時, 粒子全局搜索能力增強; 迭代后期, 慣性權(quán)值較小時, 粒子的速度不斷減小, 局部開發(fā)能力增強, 有利于粒子局部搜索的精細化。

4 改進的PSO-VMD-ED算法

VMD對信號進行分解前, 需要人為經(jīng)驗設(shè)定預(yù)設(shè)尺度k和懲罰因子α, 而參數(shù)組合[k,α]對信號分解效果有很大的影響。通過對上述算法的研究, 筆者采用混沌算法和Sigmoid函數(shù)改進PSO算法, 再利用改進的PSO優(yōu)化VMD參數(shù), 而PSO在對VMD參數(shù)進行尋優(yōu)搜索時需要確定一個適應(yīng)度函數(shù), 粒子通過與新粒子適應(yīng)度值作比較以更新粒子位置, 筆者選取包絡(luò)熵作為適應(yīng)度函數(shù)。零均值信號f(j)(j=1,2,…,N)的包絡(luò)熵EP如下

(16)

其中Pj為a(j)的歸一化形式;a(j)為信號f(j)經(jīng)Hilbert解調(diào)后得到的包絡(luò)信號。包絡(luò)熵反應(yīng)了信號的稀疏程度, 如果得到的IMF分量所含噪聲越多, 則信號的稀疏程度越小, 包絡(luò)熵值就越大; 如果IMF分量的噪聲越少, 規(guī)律性越強, 則信號的稀疏程度越大, 包絡(luò)熵值也就越小。包絡(luò)熵值最小的一個分量為含有故障特征信息的最佳分量。

當(dāng)?shù)趇個粒子處于某一位置Xid時, 計算此位置條件下VMD處理得到的所有IMF分量的包絡(luò)熵值, 將最小包絡(luò)熵值稱為局部極小熵值minEP。改進的PSO-VMD算法流程圖如圖1所示。

圖1 改進PSO優(yōu)化VMD算法流程圖Fig.1 Flow chart of improved PSO optimization VMD algorithm

改進的PSO-VMD算法基本步驟如下：

1) 初始化粒子群;

2) 令個體極值和全局極值為粒子初始位置;

3) 更新慣性權(quán)重、 加速因子以及速度模型;

4) 在不同粒子位置下對信號作VMD處理, 計算每個粒子的適應(yīng)度值minEP, 比較適應(yīng)度值大小, 求出個體極值Pbestd, 若y(Pid)

5) 計算全局極值Gbestd, 若y(Pbestd)

6) 判斷是否達到了設(shè)定的最大迭代次數(shù), 若滿足則輸出全局最優(yōu)位置及適應(yīng)度值, 算法結(jié)束； 否則, 跳轉(zhuǎn)到步驟3)。

信號濾波即降噪, 主要利用有用信號和噪聲信號在時頻域分布的不同特性進行篩選并剔除背景噪聲信號, 達到增強有用信號的目的。振蕩幅度是信號的一個重要特性, 而概率密度函數(shù)包含了信號分解模態(tài)的完整特征信息。因此, 需計算原始信號與每個模態(tài)的PDF, 通過相似性測量方法評估分解模態(tài)與原始信號的PDF之間的相似性, 故選用歐氏距離的幾何測量方法。

歐氏距離是易于理解的一種計算距離的方法, 指在M維空間中兩個點之間的真實距離, 或向量的自然長度。在二維和三維空間中的歐氏距離就是兩點之間的實際距離。M維空間中A點與B點之間的歐氏距離可表示為

(17)

其中A點坐標(biāo)為(a1,a2,…,an),B點坐標(biāo)為(b1,b2,…,bn)。

利用各模態(tài)分量與原始信號的概率密度函數(shù)之間的歐氏距離表示其與原始信號的相似程度L, 如下

L(i)=D[P(x(t)),P(Hi(t))]

(18)

其中P為概率密度函數(shù)(PDF),H為本征模態(tài)函數(shù)(BLIMF: Band-Limited Intrinsic Mode Function)。

當(dāng)坡度顯著增加時, 表明在該BLIMF之后發(fā)生的相似性迅速降低。將θ定義為相鄰兩個BLIMF之間的最大斜率

θ=max|L(i)-L(i+1)|,i=1,2,…,(N-1)

(19)

通過比較仿真信號與VMD各模態(tài)分量的PDF之間的歐氏距離, 以ED增量最大的兩個相鄰模態(tài)分量作為有效分量選取的轉(zhuǎn)折點, 將ED突變最大及其之后的模態(tài)分量都視為噪聲信號。這樣的選取方法可以篩選出真正的有效模態(tài), 達到去噪的效果。

5 實驗與分析

5.1 改進PSO-VMD-ED算法的仿真實驗分析

為了驗證改進算法的有效性, 選取含各分量為5 Hz,50 Hz,200 Hz 3個頻段的余弦信號作為原始信號并加入噪聲強度為20 dB的高斯白噪聲進行仿真分析, 仿真信號為

f(x)=1.2cos(10πt)+cos(100πt)+0.6cos(400πt)+λ

(20)

原始仿真信號如圖2a所示, 其中, 由改進的PSO-VMD算法得到的最優(yōu)解帶寬參數(shù)α=4 000, 預(yù)設(shè)尺度k=5, 信號采樣頻率為1 000 Hz, 對該信號進行分解, 經(jīng)VMD分解后得到的5個模態(tài)分量結(jié)果如圖2b～圖2f所示。VMD分解后得到的5個IMF分量的頻譜曲線如圖3所示。

圖2 原始復(fù)合信號和VMD分解得到的5個模態(tài)Fig.2 Five modes obtained from the original composite signal and VMD decomposition

由圖3可知, VMD算法具有很強的中心頻率捕獲能力, 從式(20)可知, 該復(fù)合信號的中心頻率分別為5·2π Hz, 50·2π Hz, 200·2π Hz, 對應(yīng)圖3中前3個模態(tài)分量的頻譜曲線圖的峰值： 5·2π、 50·2π、 200·2π。

圖3 VMD分解得到的5個IMF分量的頻譜曲線Fig.3 Spectrum curves of five IMF components obtained by VMD decomposition

加噪后的仿真信號經(jīng)VMD分解, 得到了5個模態(tài)分量, 分別計算加噪后的仿真信號和5個模態(tài)分量的幅值概率密度函數(shù), 如圖4所示。

由圖4可以看出, 不同模態(tài)分量的概率密度函數(shù)差別很大, 為了進一步比較各個成分間的差別, 分別計算VMD各個模態(tài)分量的PDF和加噪后仿真信號的PDF之間的歐氏距離, 如圖5所示。

圖4 加噪仿真信號和VMD模態(tài)分量的概率密度函數(shù)Fig.4 Probability density functions of the noisy simulation signal and VMD modal components

由圖5可以看出, 前3個模態(tài)分量的歐氏距離明顯小于其他分量, 這與仿真信號的構(gòu)成相符合。根據(jù)式(17)～式(19)選取增量最大的兩個相鄰模態(tài)分量作為有效分量選取的轉(zhuǎn)折點, 因此, 選擇前3個有效模態(tài)對信號進行重構(gòu)。類似地, 采用文獻[12]中提出的基于VMD與相關(guān)系數(shù)的去噪方法(VMD-CORR： Variational Mode Decomposition-Correlation Coeffificient), 將EMD自適應(yīng)分解后得到的k作為VMD的參數(shù), 輸入人為經(jīng)驗設(shè)定α值, 設(shè)置閾值為0.3, 選取相關(guān)系數(shù)大于0.3的有效分量進行重構(gòu), 對加噪后的信號進行去噪。采用Komaty等[13]提出的基于EMD與歐氏距離的去噪方法(EMD-ED), 對仿真信號先進行EMD分解, 再與歐氏距離結(jié)合, 選取有效模態(tài)進行重構(gòu)。為全面地評估以上方法的性能, 筆者分別計算以上方法重構(gòu)信號的信噪比(SNR： Signal Noise Ratio)和均方誤差(MSE： Mean Square Error), 結(jié)果對比如圖6、 圖7所示。

圖5 VMD各模態(tài)分量的PDF的歐氏距離Fig.5 Euclidean distance of PDF for each modal component of VMD

信噪比又稱訊噪比, 指系統(tǒng)中信號與噪聲的比例, 定義為

(21)

其中PS和PN分別為信號和噪聲的有效功率。當(dāng)SSNR=0, 此時信號平均功率與噪聲平均功率相等, 信號可靠性極低; 當(dāng)SSNR<0, 此時噪聲超過信號功率的1/10, 即信號已淹沒在噪聲中。

(22)

從圖6可以看出, VMD-CORR算法和EMD-ED算法對信號去噪后的信噪比分別為20.213 6 dB、 13.893 6 dB, 而筆者提出的改進VMD算法可以提高信號的信噪比, 達到23.199 6 dB。從圖7可以看出, VMD-CORR算法和EMD-ED算法對信號去噪后的均方誤差分別為0.114 2 dB、0.432 6 dB, 筆者改進的VMD算法的均方誤差則最小, 為0.082 0 dB。故改進后的算法能有效去除信號中的噪聲, 去噪效果優(yōu)于未改進的VMD-CORR算法以及EMD-ED算法, 得到了較為理想的去噪效果。

圖6 SSNR對比圖 Fig.6 SSNR contrast diagram

5.2 改進算法在管道泄漏檢測信號中的應(yīng)用

為了驗證筆者算法的有效性, 將此算法應(yīng)用到實際管道上, 并與VMD-CORR算法以及EMD-ED算法做對比實驗。筆者的實驗數(shù)據(jù)采用天然氣管道泄漏檢測模擬實驗平臺, 該實驗平臺的管道長160 m, 管徑為DN50, 管壁厚1 cm, 管道上設(shè)置了多個泄漏點, 各泄漏點間距為10 m, 由4分球閥連接。其中氣體壓力為0.5 MPa, 流量60 m3/h。檢測實驗平臺如圖8所示。筆者采用Labview編程環(huán)境對數(shù)據(jù)采集, 用于收集管道數(shù)據(jù)的采集卡是NI-9215采集卡, 采樣頻率為3 kHz。通過迅速切換球閥開關(guān)模擬管道泄漏, 圖9為實驗室采集到的泄漏數(shù)據(jù)。

圖8 天然氣管道泄漏檢測模擬實驗平臺Fig.8 Natural gas pipeline leak detection simulation test platform

管道原始信號一般為復(fù)合信號, 且存在大量的背景噪音, 筆者將采集到的管道原始信號作為改進后的PSO-VMD的輸入信號, 得到最優(yōu)解k=5,α=1 563, 再作為VMD-ED的輸入?yún)?shù), 分解后各個模態(tài)分量的歐氏距離如表1所示。

表1 管道原始信號模態(tài)分量的歐氏距離

根據(jù)各模態(tài)分量概率密度函數(shù)的歐氏距離, 由式(17)～式(19)可計算出第3個分量與第4個分量之間的增量最大, 選取增量最大的兩個相鄰模態(tài)分量作為有效分量選取的轉(zhuǎn)折點, 因此選取第1、2、3個模態(tài)分量, 對管道信號進行重構(gòu), 并計算重構(gòu)信號的信噪比和均方誤差作為性能評價指標(biāo)。計算改進算法與VMD-CORR算法以及EMD-ED算法的管道重構(gòu)信號的信噪比和均方誤差, 其對比圖如圖10和圖11所示。

圖10 SSNR對比圖Fig.11 MMSE contrast diagram

由圖10可以看出, VMD-CORR算法和EMD-ED算法的重構(gòu)信號的信噪比分別為13.511 0 dB、 13.353 6 dB, 而改進的VMD-ED算法重構(gòu)信號的信噪比顯著提高, 為20.441 8 dB。由圖11可以看出, VMD-CORR算法和EMD-ED算法的重構(gòu)信號的均方誤差分別為0.060 6 dB、 0.061 7 dB, 均高于改進算法的0.027 3 dB。圖10、 圖11的SNR和MSE的結(jié)果對比, 驗證了筆者提出的改進算法可以有效提高信號的信噪比, 降低均方誤差, 能較好地處理復(fù)雜的管道信號, 驗證了其在管道泄漏檢測中的有效性。

6 結(jié) 語

針對VMD算法分解后有效模態(tài)分量選擇困難以及去噪效果不理想等問題, 筆者采用一種改進PSO算法的VMD參數(shù)優(yōu)化方法, 有效獲取了參數(shù)組[k,α], 并與歐氏距離結(jié)合, 篩選有效的模態(tài)分量, 進行信號重構(gòu)。通過對仿真信號和實驗室管道采集信號進行試驗, 結(jié)果表明該方法相較于其他算法, 有效提高了信號的信噪比, 且均方誤差較小, 去噪效果更為顯著, 驗證了該方法在長輸管道泄漏檢測中應(yīng)用的可行性。但此算法對VMD參數(shù)組優(yōu)化所需時間較長, 可以根據(jù)該問題進行下一步研究。