江海,莫躍,何其愚
(1.中國(guó)南方電網(wǎng)超高壓輸電公司天生橋局,貴州 興義 562400;(2.云南新天地人工環(huán)境工程有限公司,云南 昆明 650000)
近年來(lái),水系統(tǒng)管道因受自然環(huán)境因素以及人為因素影響導(dǎo)致泄漏事故頻發(fā),其不僅造成水資源的浪費(fèi)及經(jīng)濟(jì)損失,還帶來(lái)極大的消防安全隱患[1-4]。因此,管道泄漏監(jiān)測(cè)逐漸成為研究熱門。至今為止,管道泄漏監(jiān)測(cè)方法,主要有負(fù)壓波法、次聲波法、壓力法等[5-8]。管道泄露監(jiān)測(cè)過(guò)程中,采集的信號(hào)中不僅包含著泄露聲信號(hào),而且存在著很多的環(huán)境噪聲,導(dǎo)致泄漏信號(hào)難以被提取識(shí)別從而造成管道泄漏監(jiān)測(cè)系統(tǒng)誤報(bào)警、漏報(bào)警,影響系統(tǒng)監(jiān)測(cè)以及定位性能[9-11],所以為了確保泄露檢測(cè)與定位的可靠性與準(zhǔn)確性,必須對(duì)采集的信號(hào)進(jìn)行去噪處理。因此,對(duì)采集信號(hào)進(jìn)行去噪處理的研究具有重要意義。
傳統(tǒng)的信號(hào)去噪方法通過(guò)設(shè)計(jì)濾波器將某一特定頻帶寬度的噪聲信號(hào)濾除以達(dá)到去噪目的,而當(dāng)噪聲信號(hào)為瞬態(tài)信號(hào)或者非平穩(wěn)信號(hào)時(shí),其頻譜通常較寬,該方法難以處理[12]。小波閾值去噪方法可以將采集信號(hào)轉(zhuǎn)換到小波域,通過(guò)設(shè)定相應(yīng)的閾值將泄露信號(hào)以及噪聲信號(hào)分離,但小波基函數(shù)以及小波閾值的選取通常較為困難[13-14]。部分學(xué)者基于信號(hào)局部的特征時(shí)間尺度構(gòu)建一種自適應(yīng)分解方法-經(jīng)驗(yàn)?zāi)B(tài)分解法(Empirical Mode Decomposition,EMD),該方法無(wú)須設(shè)置基函數(shù)以及閾值即可對(duì)信號(hào)按照不同頻率進(jìn)行分解[15-16]。2014年,Dragomiretskiy等學(xué)者提出變分模態(tài)分解方法(Variational Mode Decomposition,VMD),該方法有效解決了EMD方法分解導(dǎo)致的端點(diǎn)效應(yīng)以及頻譜混疊問(wèn)題[17-19]。目前,該方法被廣泛應(yīng)用于多個(gè)研究領(lǐng)域,如醫(yī)學(xué)圖像處理、機(jī)械故障診斷等[20-25]。
但是,VMD分解只是將原始信號(hào)分解為多個(gè)具有不同中心頻率的固有模態(tài)分量和殘余分量,而如何從VMD分解所得多個(gè)固有模態(tài)分量中挑選有效分量重構(gòu)原始信號(hào)達(dá)到去噪目的仍是一個(gè)難題[26]。此外,VMD分解效果受人為設(shè)定分解層數(shù)K值影響較大。單獨(dú)將VMD算法用作降噪,難度較大,且降噪結(jié)果受人為因素影響,需對(duì)VMD算法進(jìn)行改進(jìn)。
因此,本文提出一種基于反饋機(jī)制的VMD與巴氏距離(Bhattacharyya Distance,BD)的聯(lián)合去噪方法(VMDF-BD),該方法無(wú)須預(yù)先設(shè)定K值即可完成對(duì)有效分量的篩選。
為驗(yàn)證本文所提方法有效性,采用VMDF-BD方法對(duì)仿真信號(hào)進(jìn)行了去噪處理,并分別與多種聯(lián)合去噪算法進(jìn)行比較。仿真結(jié)果表明,該方法去噪效果優(yōu)于其他四種聯(lián)合去噪算法,即VMD-SC、VMD-HD、EMD-SC、EMD-HD。此外,利用該方法對(duì)實(shí)際采集管道信號(hào)進(jìn)行去噪處理,實(shí)驗(yàn)結(jié)果表明,該方法能夠有效去除信號(hào)中存在的噪聲。
VMD是一種自適應(yīng)、完全非遞歸的模態(tài)變分和信號(hào)處理的方法,是一種新的可變尺度的信號(hào)分解方法,其能將實(shí)值信號(hào)分解為K個(gè)具有特定稀疏特性的模態(tài)分量。該方法具有可以確定模態(tài)分解個(gè)數(shù)的優(yōu)點(diǎn),其自適應(yīng)性表現(xiàn)在根據(jù)實(shí)際情況確定所給序列的模態(tài)分解個(gè)數(shù),隨后的搜索和求解過(guò)程中可以自適應(yīng)地匹配每種模態(tài)的最佳中心頻率和有限帶寬,并且可以實(shí)現(xiàn)固有模態(tài)分量(IMF)的有效分離、信號(hào)的頻域劃分、進(jìn)而得到給定信號(hào)的有效分解成分,最終獲得變分問(wèn)題的最優(yōu)解。其本質(zhì)是通過(guò)迭代求解所構(gòu)建約束變分問(wèn)題的最優(yōu)解從而得到估計(jì)帶寬之和最小的K個(gè)模態(tài)分量,核心思想是構(gòu)建和求解變分問(wèn)題。其具體的分解步驟如下:
首先,通過(guò)希爾伯特變換求解各個(gè)模態(tài)分量的單邊頻譜,如式(1)
其中,*為卷積運(yùn)算,δ(t)為單位脈沖函數(shù),uk(t)為第k個(gè)模態(tài)分量函數(shù),k=1,2,…,K。之后通過(guò)引入各模態(tài)的估計(jì)中心頻率指數(shù)項(xiàng),將各模態(tài)頻譜分別調(diào)制到相應(yīng)的基頻帶上,如式(2)
其中,·為乘法運(yùn)算,ωk為第k個(gè)模態(tài)分量的中心頻率。然后,通過(guò)求式(2)的梯度平方L2范數(shù),即求得各個(gè)模態(tài)分量的估計(jì)帶寬,如式(3)
其中,?t為對(duì)函數(shù)求時(shí)間t的導(dǎo)數(shù)。最后,加入約束條件,則該約束變分問(wèn)題最終表示為
其中,s為原始信號(hào)。為較好求解該約束變分問(wèn)題,引入二次懲罰因子α以及拉格朗日乘法算子λ(t)將其轉(zhuǎn)換為無(wú)約束變分問(wèn)題,則式(4)可表示為
最終,通過(guò)交替方向乘子算法(Alternate Direction Method of Multipliers,ADMM)求解式(5),得到各參量迭代更新公式分別如式(6)(7)和(8)所示
其中,m為迭代次數(shù),τ為更新參數(shù),其控制拉格朗日乘法算子的收斂速度,迭代終止條件為,ε為收斂精度。
它克服了EMD方法存在端點(diǎn)效應(yīng)和模態(tài)分量混疊的問(wèn)題,并且具有更堅(jiān)實(shí)的數(shù)學(xué)理論基礎(chǔ),可以降低復(fù)雜度高和非線性強(qiáng)的時(shí)間序列非平穩(wěn)性,分解獲得包含多個(gè)不同頻率尺度且相對(duì)平穩(wěn)的子序列,適用于非平穩(wěn)性的序列。
巴氏距離主要用于統(tǒng)計(jì)兩概率分布的相似性,其結(jié)果與巴氏系數(shù)密切相關(guān)。巴氏距離的具體計(jì)算方法如下:假設(shè)在同一定義域x內(nèi),存在兩個(gè)離散概率分布分別為p和q,則巴氏距離為
其中,BD(p,q)為Bhattacharyya系數(shù),其主要用于度量?jī)山y(tǒng)計(jì)樣本之間的重疊量,表示為
本文中,由于采集信號(hào)及其VMD分解所得各模態(tài)分量的概率分布未知,因此在計(jì)算各分量的巴氏距離之前,需要對(duì)各個(gè)分量進(jìn)行密度估計(jì),根據(jù)密度估計(jì)的結(jié)果再求巴氏距離。
本文基于反饋機(jī)制的VMD與巴氏距離構(gòu)建的聯(lián)合去噪方法,首先將采集信號(hào)作為輸入信號(hào)進(jìn)行2層VMD得到兩個(gè)模態(tài)分量,之后分別計(jì)算兩模態(tài)分量與采集信號(hào)之間的距離,將距離較小的作為較純凈的模態(tài)分量反饋回輸入端從輸入信號(hào)中減去作為新的輸入信號(hào),如此迭代直至信號(hào)完全分解,最后,利用歷次迭代所得較純凈模態(tài)分量相加重構(gòu)原始采集信號(hào),達(dá)到去噪目的。其主要步驟如下:
第一步,初始化分解層數(shù),也即VMD分解所得模態(tài)數(shù)K=2,初始化迭代次數(shù)kk=1;初始化輸入信號(hào)為采集信號(hào),Sinkk=S。
第二步,根據(jù)式(6)、式(7)、式(8)對(duì)輸入信號(hào)Sinkk進(jìn)行VMD,得到兩模態(tài)分量ukk1,ukk2。
第三步,根據(jù)式(9)、式(10)分別計(jì)算兩模態(tài)分量與采集信號(hào)之間的巴氏距離,得到BDkk1,BDkk2,取二者中最小距離對(duì)應(yīng)模態(tài)分量作為較純凈模態(tài)分量,即upurekk=ukk1,i=argmin(BDkk1,BDkk2),其中argmin()表示取最小值的索引。
第四步,判斷當(dāng)前采集信號(hào)是否完全分解,其判別條件為:當(dāng)當(dāng)前分解所得兩模態(tài)分量對(duì)應(yīng)巴氏距離的最小值仍大于前次分解所得兩模態(tài)分量對(duì)應(yīng)巴氏距離的最大值時(shí),則該信號(hào)完全分解,即min(BDkk1,BDkk2)>max(BDkk-11,BDkk-12)。若信號(hào)完全分解,則停止迭代,并利用歷次迭代所得較純凈模態(tài)分量重構(gòu)原始信號(hào),即,;否則,轉(zhuǎn)第五步;
第五步,將純凈模態(tài)分量返回輸入端并從原輸入信號(hào)中減去,得新的輸入信號(hào),即,Sinkk+1=Sinkkupurekk,并重復(fù)第2步~第4步。
為驗(yàn)證本文所提方法的去噪效果以及有效性,本文將對(duì)VMDF-BD算法進(jìn)行驗(yàn)證。采用用仿真信號(hào)中添加白噪聲,用降噪方法對(duì)該疊加信號(hào)進(jìn)行降噪,并從多角度對(duì)比其降噪結(jié)果。從而達(dá)到驗(yàn)證VMDF-BD算法降噪效果是否可靠的目的。
本文在正弦信號(hào)中加入已知的高斯白噪聲模擬管道采集信號(hào)。仿真實(shí)驗(yàn)運(yùn)行環(huán)境為MATLAB R2020b,仿真實(shí)驗(yàn)所用正弦信號(hào)為x(t)=sin(2π×13×t)+cos(2π×96×t)+cos(2π×145×t),該仿真信號(hào)中所疊加的頻率分別為13Hz、96Hz、145Hz,這三種頻率也是判斷降噪效果的標(biāo)準(zhǔn)。采樣頻率為1 kHz,采樣點(diǎn)數(shù)為1200點(diǎn)。
針對(duì)上述的仿真信號(hào)進(jìn)行VMD分解,根據(jù)分解不同分解層數(shù)的中心頻率來(lái)確定最優(yōu)分解層數(shù)。
仿真實(shí)驗(yàn)正弦信號(hào)由上述公式仿真而成,分別由頻率為13Hz、96Hz、145Hz的正弦信號(hào)疊加。本實(shí)驗(yàn)采用手動(dòng)確定k值的大小來(lái)找出最佳得分解層數(shù)。采用不同K值對(duì)仿真信號(hào)進(jìn)行VMD,然后對(duì)分解所得的模態(tài)分量求其中心頻率,所得各模態(tài)分量對(duì)應(yīng)中心頻率如下表所示。
表1 不同K值VMD所得各模態(tài)分量中心頻率Tab.1 Center frequencies of modal components obtained by VMD with different K values
由表1可知,當(dāng)K=2,即分解層數(shù)小于仿真信號(hào)所含不同頻率成分?jǐn)?shù)時(shí),可以看到其只分解出中心頻率為13Hz附近的分量,高頻分量并不存在,屬于欠分解;當(dāng)k=3,可以看到其中心頻率包括12.99999,95.9918,145.0022,根據(jù)上述仿真可知,該正玄信號(hào)的組成頻率為13 Hz、96Hz、145Hz,因此可只k=3的分解尺度合適;當(dāng)K=4以及K=5,即分解層數(shù)大于仿真信號(hào)所含不同頻率成分?jǐn)?shù),可以看到不僅分解得到了與仿真信號(hào)頻率成分相近的模態(tài)分量,而且還存在中心頻率為149.4952Hz、66.4038Hz、149.7513Hz的成分,屬于過(guò)分解。由此可知,VMD分解性能受K值影響較大。因此對(duì)于如何確定分解層數(shù)成為影響去噪效果的重要因素。
據(jù)此,本文提出VMDF-BD方法,該方法是基于反饋機(jī)制的VMD與巴氏距離構(gòu)建的聯(lián)合去噪方法,該方法的反饋機(jī)制能夠較好地迭代終止判斷從而能夠較好地達(dá)到去噪的目的。
迭代何時(shí)終止是去噪是否有效果的重要影響因素,因此為驗(yàn)證本文所提VMDF-BD方法所設(shè)定迭代終止條件成立,對(duì)上述仿真實(shí)驗(yàn)正弦信號(hào)應(yīng)用VMDF-BD方法進(jìn)行分解,根據(jù)其分解后的結(jié)果求其各次迭代后所得分量中心頻率及巴氏距離如表2所示。
表2 歷次迭代各模態(tài)分量中心頻率以及巴氏距離Tab.2 Central frequency and pasteurization distance of each modal component in previous iterations
由表2可知,當(dāng)?shù)恋谌螘r(shí),本次分解所得兩模態(tài)分量對(duì)應(yīng)巴氏距離分別為0.0824、1.4973,其最小值依然比上次分解所得距離最大值0.0817大,因此根據(jù)本文所設(shè)迭代終止條件,當(dāng)前信號(hào)已經(jīng)完全分解。而從中心頻率的變化上看,第三次迭代產(chǎn)生的新分量的中心頻率與正弦信號(hào)中所含三頻率成分均無(wú)關(guān),根據(jù)1.3中的分解完成判別方法,說(shuō)明當(dāng)前該信號(hào)已經(jīng)完全分解,其分解所得各分量為歷次迭代分解所得較純凈模態(tài)分量。因此,VMDF-BD方法所設(shè)計(jì)迭代終止條件成立。
對(duì)于信號(hào)處理來(lái)說(shuō),去噪效果往往會(huì)在很大程度上影響結(jié)果的判斷,本實(shí)驗(yàn)的最終目的也是找到較合適的去噪方法。因此本文在上述仿真實(shí)驗(yàn)正弦信號(hào)中加入snr=5dB的高斯白噪聲,分別采用VMDF-BD,VMD-SC,EMD-SC,VMD-HD,EMDHD對(duì)其進(jìn)行去噪處理。為保證實(shí)驗(yàn)的科學(xué)性,VMD-SC以及VMD-HD中VMD分解層數(shù)K值設(shè)定與EMD-HD以及EMD-SC自適應(yīng)分解層數(shù)一致,均為9層。將加入白噪聲的信號(hào)采用上述幾種方法去噪之后重構(gòu),各方法所得重構(gòu)信號(hào)與原始信號(hào)時(shí)域?qū)Ρ葓D如圖1所示,頻域?qū)Ρ葓D如圖2所示。
圖1 時(shí)域?qū)Ρ葓DFig.1 Time domain comparison diagram
從圖1可知,VMD-HD以及EMD-HD所得去噪結(jié)果在時(shí)域方向即與原始信號(hào)相去甚遠(yuǎn),在去掉噪聲信號(hào)的同時(shí)也去掉了較多的有效信號(hào),雖然信號(hào)較為平滑,但同時(shí)也導(dǎo)致信號(hào)嚴(yán)重失真。VMDSC以及EMD-SC與原始信號(hào)波形基本一致,但其內(nèi)仍包含有較多噪聲信號(hào)。而本文提出的VMDFBD所得去噪信號(hào)與原始信號(hào)相比,波形基本一致,且更平滑,去噪效果顯著。為了更清楚的表征各種方法的去噪效果,本實(shí)驗(yàn)也在頻域方向做了對(duì)比。
從圖2中可知,與時(shí)域?qū)Ρ葓D相似,VMD-HD以及EMD-HD所得去噪后信號(hào)頻譜與原始信號(hào)頻譜相比,均缺少了兩個(gè)頻率較高的分量,即這兩種去噪方法不可靠。VMD-SC所得去噪后信號(hào)頻譜與原始信號(hào)頻譜相比,其去除了部分分布在低頻的噪聲,高頻噪聲并不能很好的去除,而EMD-SC所得去噪后信號(hào)頻譜與原始信號(hào)頻譜相比,基本一致,去噪效果不明顯。本文所提VMDF-BD所得去噪后信號(hào)頻譜與原始信號(hào)頻譜相比,在原有頻率分量處基本保持一致,在其余處頻譜幅值趨于0,很好地達(dá)到了去噪的要求,從頻域上觀察,其去噪效果顯著。
圖2 頻域?qū)Ρ葓DFig.2 frequency domain comparison diagram
去噪效果的好壞一般采用均方根誤差RMSE和信噪比SNR兩個(gè)指標(biāo),RMSE是預(yù)測(cè)值與真實(shí)值偏差的平方與觀測(cè)次數(shù)n比值的平方根,一般情況其值越小越好,SNR指有用信號(hào)的功率和噪聲信號(hào)功率的比值,其比值越大,說(shuō)明去噪效果越好。因此本實(shí)驗(yàn)為進(jìn)一步說(shuō)明去噪性能,引入去噪性能指標(biāo)對(duì)各去噪算法進(jìn)行比對(duì),即RMSE和SNR,分別定義為
其中,s(t)為原始信號(hào),s′(t)為去噪后信號(hào),T為信號(hào)總長(zhǎng)度。則各去噪方法SNR以及RMSE計(jì)算結(jié)果,如表3所示。
表3 各去噪方法SNR以及RMSETab.3 SNR and RMSE of each denoising method
一般來(lái)說(shuō),信噪比越大,均方根誤差越小,則去噪效果越好。從表3可知看出:RMSE指標(biāo)的排序?yàn)?VMDF-BD<VMD-SC<EMD-SC<VMD-HD<EMD-HD;SNR指標(biāo)的排序?yàn)?VMDF-BD>VMD-SC>EMD-SC>VMD-HD>EMD-HD。本文所提方法VMDF-BD的SNR最大,RMSE最小。VMD-SC與EMD-SC的結(jié)果次之,去噪效果相似,均優(yōu)于VMDHD與EMD-HD的去噪效果。與上述時(shí)域和頻域的判斷結(jié)果一致,在上述五種去噪算法中,本文所提VMDF-BD方法去噪效果最好。
將上述方法用于實(shí)際采集的信號(hào)分析。本次實(shí)驗(yàn)裝置如圖3所示,主要包括三個(gè)部分:次聲波傳感器,數(shù)據(jù)采集及傳輸模塊,PC端。
傳感器安裝于管道之上采集信號(hào)并通過(guò)同軸電纜傳輸給數(shù)據(jù)采集與傳輸模塊,并經(jīng)其轉(zhuǎn)換為數(shù)字信號(hào)后經(jīng)雙絞線發(fā)送至PC機(jī),在PC機(jī)上可進(jìn)行信號(hào)分析,從而進(jìn)行管道泄露判斷及泄露點(diǎn)定位。
圖3 實(shí)驗(yàn)裝置示意圖Fig.3 Schematic diagram of experimental device
本實(shí)驗(yàn)的實(shí)際實(shí)驗(yàn)參數(shù)為:管道為約200米長(zhǎng)型號(hào)為DN20的PVC直水管道,管內(nèi)水壓為0.1Mpa。數(shù)據(jù)采集與傳輸模塊內(nèi)采樣頻率為500Hz,采樣精度為24 bit。
通過(guò)閥門開(kāi)度控制管道有無(wú)泄漏。為加快信號(hào)處理的速度,從采集信號(hào)中截取1200點(diǎn)進(jìn)行分析。對(duì)該信號(hào)的分析包括時(shí)域和頻域兩個(gè)方向進(jìn)行去噪前后對(duì)比。本實(shí)驗(yàn)所采用的去噪方法即上述的VMDF-BD法。原始信號(hào)的去噪前后時(shí)域?qū)Ρ葓D如圖4所示。
圖4 去噪前后信號(hào)時(shí)域?qū)Ρ葓DFig.4 time domain comparison of signals before and after denoising
從圖4可知,去噪后信號(hào)與原始信號(hào)相比,波形基本一致,保留了較多的信號(hào)細(xì)節(jié),但更平滑,可知該去噪方法效果顯著。為了清楚地展現(xiàn)去噪效果,本實(shí)驗(yàn)將原始信號(hào)和去噪后的信號(hào)均進(jìn)行了時(shí)頻轉(zhuǎn)換,其頻域?qū)Ρ葓D如圖5所示。
圖5 去噪前后頻域信號(hào)對(duì)比圖Fig.5 Comparison of frequency domain signals before and after denoising
據(jù)資料顯示,管道泄露所產(chǎn)生的頻率處較低頻的范圍,因此高頻部分為噪聲信號(hào)。從圖5可看出,去噪后的信號(hào)的幅值均低于去噪前的信號(hào),尤其在高頻部分,本文提出VMDF-BD法用于管道泄露聲波信號(hào)去,能夠達(dá)到去噪目的。因此,本文提出的方法可以有效去除實(shí)際管道采集信號(hào)中引入的寬帶噪聲。
針對(duì)管道采集信號(hào)中存在的噪聲干擾問(wèn)題,提出了一種基于變分模態(tài)分解與巴氏距離的聯(lián)合去噪算法。通過(guò)對(duì)仿真信號(hào)降噪的方法,分析得出了該算法具有較好降噪效果的結(jié)論。并將該算法實(shí)際應(yīng)用在管道采集的信號(hào)中進(jìn)行降噪處理,其降噪效果較好。
該算法取得了較好的去噪效果。其主要優(yōu)點(diǎn)如下:
(1)采用巴氏距離來(lái)度量各分解所得模態(tài)分量與原始信號(hào)之間的相似性,使去噪信號(hào)減少失真。
(2)引入反饋機(jī)制,無(wú)須人為設(shè)定VMD分解層數(shù),確保采集信號(hào)完全分解,以盡可能去除噪聲。
(3)仿真實(shí)驗(yàn)中,本文所提方法去噪效果明顯優(yōu)于其余四種去噪算法。實(shí)際測(cè)試中,本文所提方法能夠有效去除實(shí)測(cè)數(shù)據(jù)的噪聲,并且保留原始信號(hào)的基本特性。
因此,本文所提去噪方法具備一定的研究?jī)r(jià)值,并且能夠?yàn)橄乱徊降男盘?hào)識(shí)別、分類提供幫助。