韓文翔， 章蘭珠

(華東理工大學(xué) 機(jī)械與動力工程學(xué)院，上海 200237)

管道作為運(yùn)輸介質(zhì)的載體，具有經(jīng)濟(jì)、安全、不間斷性、運(yùn)輸量大等多重優(yōu)勢。在提供便利的同時，由于腐蝕、老化和人為破壞等因素的影響，一旦發(fā)生泄漏事故，后果頗為嚴(yán)重，損壞國民經(jīng)濟(jì)甚至威脅生命安全[1]。針對管道泄漏的檢測已經(jīng)有示蹤法、負(fù)壓波法、聲學(xué)檢測法等較成熟的檢測方法。其中聲發(fā)射(acoustic emission，AE)技術(shù)作為一種具有在線動態(tài)監(jiān)測特點(diǎn)的新型檢測技術(shù)，可以及時發(fā)現(xiàn)管道的早期缺陷，提高檢測效率，具有良好的工程實(shí)際應(yīng)用價值[2]。

目前在管道泄漏領(lǐng)域的研究較多集中在單點(diǎn)泄漏場景，利用時差定位法、神經(jīng)網(wǎng)絡(luò)等方法可有效解決管道泄漏的定位問題。而真實(shí)工況下往往會發(fā)生多點(diǎn)泄漏的現(xiàn)象，多點(diǎn)泄漏由于泄漏信號之間的相互干擾給研究帶來了一定的難度，目前還沒有形成成熟的管道多點(diǎn)泄漏定位方法。高華[3-4]等提出一種改進(jìn)差分進(jìn)化算法的管道泄漏定位方法，將采集到的信號經(jīng)過小波消噪后并根據(jù)信號幅值建立改進(jìn)差分進(jìn)化算法的目標(biāo)優(yōu)化函數(shù)，通過對其求解獲得泄漏源位置，但是該方法未考慮彎管、變徑等工況對信號幅值的干擾。而聲發(fā)射信號可以避免復(fù)雜工況對信號特性的影響，管道在發(fā)生泄漏時由于內(nèi)外壓差導(dǎo)致內(nèi)部介質(zhì)與外部環(huán)境之間的摩擦從而形成高頻應(yīng)力波[5]，這種機(jī)械波組成的聲發(fā)射信號與語音信號同為非平穩(wěn)信號，可借鑒語音信號處理中的盲辨識法應(yīng)用于管道多點(diǎn)泄漏的檢測。盲辨識問題實(shí)質(zhì)上是盲源分離中混合矩陣的估計問題，根據(jù)信號源數(shù)量與傳感器數(shù)量之間的關(guān)系可分為超定、適定和欠定3種類型[6]。目前該領(lǐng)域的研究熱點(diǎn)主要在欠定盲源分離，而超定和適定兩種情況下的盲源分離技術(shù)相對成熟，本文的研究重點(diǎn)在于利用信號處理技術(shù)提取有效的泄漏特征信號并結(jié)合盲源分離算法實(shí)現(xiàn)管道的多點(diǎn)泄漏源定位，因此選用基于快速獨(dú)立成分分析(fast independent component analysis，F(xiàn)ast ICA)的適定盲源分離算法。

在盲源分離的應(yīng)用中，往往將噪聲分量當(dāng)做獨(dú)立成分分離，考慮到去噪效果和盲源分離后排序未知等問題的限制，本文先對采集信號進(jìn)行去噪預(yù)處理。通過對變分模態(tài)分解(variational mode decomposition，VMD)分解后有效模態(tài)選取依據(jù)的研究提出一種結(jié)合互譜的信號去噪算法，并通過試驗(yàn)驗(yàn)證了該去噪方法的有效性和多點(diǎn)泄漏應(yīng)用的可靠性。

1 基于聲發(fā)射技術(shù)的時差定位原理

聲發(fā)射技術(shù)是近年來得到快速發(fā)展的一種動態(tài)無損檢測新技術(shù)，該技術(shù)可以在不破壞被檢結(jié)構(gòu)的情況下對構(gòu)件進(jìn)行全天在線監(jiān)控，并能對結(jié)構(gòu)中出現(xiàn)的疲勞及失效做出早期預(yù)警，具有實(shí)時監(jiān)測、動態(tài)監(jiān)測、覆蓋面廣等優(yōu)點(diǎn)[7]。管道發(fā)生泄漏時，釋放的氣體與管道裂口通過摩擦所產(chǎn)生的能量以應(yīng)力波的形式釋放，此過程稱為聲發(fā)射現(xiàn)象，通過傳感器對聲發(fā)射信號的采集和分析可實(shí)現(xiàn)管道泄漏的檢測。但是由于聲發(fā)射信號的本身特性及環(huán)境因素的影響，需要對采集信號進(jìn)行適當(dāng)?shù)姆治龊吞幚聿拍苓M(jìn)一步精確確定泄漏點(diǎn)的位置。

在泄漏源定位方面，常采用時差定位法[8]，其原理如下

(1)

式中：x為被檢測管道泄漏點(diǎn)的位置，即泄漏點(diǎn)到上游傳感器距離，m；Δt為泄漏信號到達(dá)兩個傳感器的時間差，s；v為泄漏信號在管道中的傳播速度，m/s。

由于兩傳感器之間的距離是已知的，因此只需知道Δt便可得出泄漏源的位置。互相關(guān)函數(shù)是計算兩個信號相關(guān)性的常用方法，根據(jù)互相關(guān)時域分析圖中最大相關(guān)峰值進(jìn)而計算得到上，下游有效泄漏信號到達(dá)兩個傳感器的時差Δt。

2 泄漏特征信號提取

管道運(yùn)行過程中環(huán)境噪聲的干擾會對定位結(jié)果造成顯著影響，因此有必要先對采集到的信號進(jìn)行去噪處理，目前信號去噪領(lǐng)域已具很多成熟的算法，例如小波分解、小波包分解等信號處理算法。雖然小波分解在工業(yè)中的應(yīng)用已經(jīng)有了長遠(yuǎn)的發(fā)展和進(jìn)步，但是小波基函數(shù)的選擇和分解層數(shù)的確定始終是擺在分析人員面前的問題[9]。為此Huang等[10]提出一種完全基于數(shù)據(jù)驅(qū)動的自適應(yīng)分解方法——經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition，EMD)，它將信號分解成有限個從高頻到低頻的固有模態(tài)函數(shù)(intrinsic mode function，IMF) 和殘余分量之和, 能夠很好地分析處理非線性、非平穩(wěn)信號。由于信號中存在間斷信號和脈沖干擾等異常事件的干擾，使得EMD 存在模態(tài)混疊的問題，這會降低EMD 的降噪效果?；诖?，研究人員又陸續(xù)發(fā)明了集合經(jīng)驗(yàn)?zāi)B(tài)分解(ensemble empirical mode decomposition，EEMD)，VMD等信號處理算法[11]。

2.1 互譜分析原理

互譜可由互相關(guān)函數(shù)經(jīng)傅里葉變換求得從而在頻域內(nèi)描述兩信號s1(k)和s2(k)之間的相關(guān)性[12]。對s1(k)和s2(k)的相關(guān)函數(shù)作傅里葉變換得互譜公式如下

Ss1s2(f)=αSss(f)ej2πfτ=αSss(f)ejθs1s2(f)

(2)

式中，Sss(f)為采集信號的自譜。由式(2)可得時間延遲為τ，則相位差為

θs1s2(f)=2πfτ

(3)

在采集信號的頻帶內(nèi)，θs1s2(f)與頻率成線性關(guān)系，斜率為2πτ?；プV的計算公式如下

(4)

式中：S1(f)和S2(f)分別為s1(k)和s2(k)的離散傅里葉變換系數(shù)；*為復(fù)共軛。

互譜估計常采用Welch平均周期圖法表示為

Gs1s2(f)=Gs1s2(f)+jQs1s2(f)

(5)

式中，Gs1s2和Qs1s2(f)分別為共譜密度和正交譜，可得相位譜如下

(6)

將s1(k)和s2(k)的互譜相位譜作差可得互譜相位差譜，差譜值越接近水平表明頻率相關(guān)性越高。

利用Welch平均周期圖法可得S1(f)和S2(f)的自譜Gx1x2(f)和Gx2x2(f)，則s1(k)和s2(k)的幅值平方相干函數(shù)為

(7)

相干函數(shù)可從頻域的角度反應(yīng)兩段信號的相關(guān)程度，頻帶的相干函數(shù)值越大表明該頻段的信噪比越高，因此可根據(jù)互譜相位差譜和相干函數(shù)估計泄漏特征頻帶。

試驗(yàn)采用北京軟島科技研發(fā)的聲發(fā)射儀器和實(shí)驗(yàn)室搭建的管道裝置，在泄漏源上下游兩端各放置一個聲發(fā)射傳感器，選用2.5 M的采樣率，計算兩段信號的互譜相位差譜和相關(guān)函數(shù)如圖1所示。

圖1 s1(t)和s2(t)的相干函數(shù)及互譜相位差譜Fig.1 Coherence function and cross-spectral phase difference spectrum of s1(t) and s2(t)

由圖1(a)可知，相關(guān)性較高的頻段主要分布在0～500 kHz；由圖1(b)可0～500 kHz內(nèi)未發(fā)現(xiàn)明顯水平差譜頻段，分析其原因可能是泄漏特征信號與噪聲信號相互重疊從而影響了互譜效果。為此提出一種互譜結(jié)合變分模態(tài)分解的去噪算法。

2.2 結(jié)合互譜的變分模態(tài)分解去噪算法

VMD通過迭代搜尋變分模型最優(yōu)解來確定每個分量的頻率中心及帶寬，相比于EMD，VMD具有堅實(shí)的理論基礎(chǔ)，有效地彌補(bǔ)了 EMD和EEMD等不足[13]。采集信號經(jīng)過VMD后得到若干個IMF(FIM)分量和殘余分量，公式如下

(8)

式中:k為VMD分解后的IMF分量個數(shù)；γ(t)為分解后的殘余分量。

2.2.1 分解層數(shù)的確定

相比于EMD自適應(yīng)分解信號的特點(diǎn)，VMD的分解層數(shù)K由人為輸入，如果K過小，低頻段將被劃分為高頻段，給模式函數(shù)的求解帶來了困難。但是如果K值過大，會導(dǎo)致過度分解，計算時間過長。因此，分解層數(shù)的確定是VMD在應(yīng)用過程中必須考慮的問題[14]。

由于采樣率選用2.5 MHz，選用經(jīng)驗(yàn)值K=10將采集信號進(jìn)行VMD分解后肘頻圖，如圖2所示。

圖2 VMD分解結(jié)果時頻圖Fig.2 Time-frequency diagram of VMD decomposition result

由圖2(b)可知，采集信號從低頻到高頻依次分解，并且只要分解層數(shù)選取適當(dāng)可有效避免模態(tài)混疊等問題。由于高頻部分的分解毫無意義，不僅消耗算力并且容易造成過度分解。因此可以通過互譜大致定位有效頻段并結(jié)合濾波實(shí)現(xiàn)信號的去噪預(yù)處理，通過巴特沃斯濾波器濾出其中的高相似度信號，再對濾波后的信號進(jìn)行VMD分解。

將信號進(jìn)行互譜分析并設(shè)計低通濾波器濾出0～500 kHz的頻域圖，如圖3所示。

圖3 濾波后頻域圖Fig.3 Frequency domain graph after filtering

可通過是否發(fā)生過度分解為界限確定分解層數(shù)，對濾波后的信號進(jìn)行4層分解和5層分解，得出的頻譜圖如圖4所示。

圖4 VMD分解4層和5層頻域圖Fig.4 Result diagram of VMD decomposition layer 4 and layer 5

由圖4可知，當(dāng)分解至5層時第4個IMF和第5個IMF出現(xiàn)了過度分解現(xiàn)象，因此將分解層數(shù)定位4層。

2.2.2 有效模態(tài)選取

基于時差定位的管道泄漏定位試驗(yàn)中，主要的泄漏特征信號往往包含在VMD分解后的一個或多個模態(tài)中，因此針對有效模態(tài)的選取進(jìn)行了深入研究。

目前已有的研究中，往往通過相關(guān)系數(shù)、能量值和信息熵[15]等熵值判斷有效模態(tài)的選取，而以上方法都有一定的局限性。由于高頻噪聲的干擾，通過相關(guān)系數(shù)和能量值選取有效模態(tài)往往導(dǎo)致結(jié)果誤差較大。而后者依據(jù)的理論是聲發(fā)射作為一種動態(tài)檢測方法，泄漏特征成分越多，對應(yīng)信號的隨機(jī)性越強(qiáng)，從而信息熵、樣本熵等熵值較大，該方法具有一定的適用性，但是在泄漏率低的情況下仍然無法辨識出有效模態(tài)。

將采集信號經(jīng)過VMD分解成10個模態(tài)后，分別計算每個模態(tài)的能量值，相關(guān)系數(shù)和信息熵結(jié)果如圖5所示。

由圖5可知，能量值最高的模態(tài)為IMF1，信息熵最大的模態(tài)是IMF6，相關(guān)系數(shù)最大的模態(tài)是IMF4?？梢娀谛畔㈧睾拖嚓P(guān)系數(shù)選擇的的IMF頻段并不是泄漏特征信號所在的主要頻段。為保證該試驗(yàn)的科學(xué)性，選取0.1 L/min(小泄漏)，0.5 L/min(中泄漏)，1 L/min(大泄漏)3種泄漏率下各10組試驗(yàn)，按經(jīng)驗(yàn)值將信號分解至10層后，分別計算各個模態(tài)的能量值，相關(guān)系數(shù)和信息熵，記錄每個模態(tài)得到最大值的次數(shù)，得到的結(jié)果如圖6～圖8所示

圖6 0.1 L/min 泄漏率下各指標(biāo)對應(yīng)模態(tài)選取系數(shù)Fig.6 Modal selection coefficient of each index under 0.1 L/min leakage rate

圖7 0.5 L/min 泄漏率下各指標(biāo)對應(yīng)模態(tài)選取系數(shù)Fig.7 Modal selection coefficient of each index under 0.5 L/min leakage rate

圖8 1.0 L/min 泄漏率下各指標(biāo)對應(yīng)模態(tài)選取次數(shù)Fig.8 Number of modal selection corresponding to each index under 1.0 L/min leakage rate

由圖6～圖8可知，基于單一指標(biāo)的情況下容易出現(xiàn)高頻模態(tài)被選用的情況從而導(dǎo)致定位結(jié)果存在較大的誤差。為此提出一種結(jié)合相關(guān)系數(shù)、信息熵和能量值的綜合指標(biāo)，其中相關(guān)系數(shù)的選用為避免噪聲的干擾采用基于互譜濾波后的信號與分解后的IMF進(jìn)行互相關(guān)計算。將3個指標(biāo)分別歸一化后相加，取值最大的IMF為包含最多泄漏特征信號的IMF。

設(shè)能量值為a，信息熵為b，相關(guān)系數(shù)為c，將每個指標(biāo)對應(yīng)的數(shù)值進(jìn)行歸一化后再將不同指標(biāo)的數(shù)值相加即得到綜合指標(biāo)xzonghe如下

(9)

信號經(jīng)過VMD并通過綜合指標(biāo)選取IMF的結(jié)果分布如圖9所示。

圖9 不同泄漏率下綜合指標(biāo)選取次數(shù)Fig.9 Selection times of comprehensive indicators under different leakage rates

由圖9可知，采用綜合指標(biāo)后選取的模態(tài)基本上集中在低頻模態(tài)，結(jié)果與理論泄漏特征信號所在頻段相符合，所以采用綜合指標(biāo)選取有效的IMF比單一指標(biāo)具有更高的可靠性。

通過以上信號處理方式保證了去噪后的信號包含了盡可能多的泄漏特征信號，進(jìn)而在此基礎(chǔ)上通過FastICA算法將不同泄漏源產(chǎn)生的泄漏信號分離出來。

3 基于FastICA的盲源分離應(yīng)用

FastICA其目標(biāo)就是計算矩陣A的逆矩陣B，使其作用在觀測信號x(t)上，從而可以分離解耦出具有非高斯性的接近原始信號s(t)的源信號[16]，基本原理如圖10所示。

圖10 ICA模型Fig.10 ICA model

FastICA算法有基于峭度、基于似然最大、基于負(fù)熵最大等形式，這里采用基于負(fù)熵最大的FastICA算法，該算法以負(fù)熵最大作為一個搜尋方向，可以實(shí)現(xiàn)順序地提取獨(dú)立源，可以有效地將不動點(diǎn)迭代所帶來的優(yōu)良算法特性與負(fù)熵所帶來的更好統(tǒng)計特性結(jié)合起來，使得收斂更加快速、穩(wěn)健。

綜合以上分析，本文提出一種互譜結(jié)合VMD的信號處理方法，并選用綜合指標(biāo)最大值作為選取有效IMF的依據(jù)，在管道同一側(cè)對去噪后的信號進(jìn)行盲源分離計算得到兩個包含不同泄漏源特征的信號，再根據(jù)相關(guān)系數(shù)匹配出管道兩側(cè)由同一泄漏源產(chǎn)生的信號，最后根據(jù)時差定位得到定位結(jié)果，完整的泄漏定位方法如圖11所示。

圖11 多點(diǎn)泄漏定位算法Fig.11 Multi-point leak location algorithm

4 多點(diǎn)泄漏定位試驗(yàn)

為證明所提出方法的有效性，研究以兩泄漏點(diǎn)為例，利用北京軟島科技公司研發(fā)的聲發(fā)射采集設(shè)備對實(shí)驗(yàn)室搭建的管道泄漏模型展開研究，在管道上、下游分別放置兩個聲發(fā)射傳感器，通過轉(zhuǎn)動泄漏閥模擬泄漏源，設(shè)置兩端傳感器距離為1 300 mm，以左端傳感器為基準(zhǔn)，兩泄漏源位置分別為300 mm和1 000 mm，泄漏率大小分別為1 L/min和0.5 L/min,設(shè)置管道內(nèi)壓為0. 3 MPa。試驗(yàn)方案如圖12所示。

圖12 多點(diǎn)泄漏定位試驗(yàn)方案(mm)Fig.12 Experimental scheme of multi-point leak location(mm)

試驗(yàn)裝置如圖13所示。

圖13 管道泄漏模型Fig.13 Pipeline leakage model

以上游的兩個傳感器為例，兩傳感器采集的混合信號s1和s2頻域圖如下圖14、圖15所示。

圖14 混合信號s1頻域圖Fig.14 Frequency domain diagram of mixed signal s1

圖15 混合信號s2頻域圖Fig.15 Frequency domain diagram of mixed signal s2

將信號進(jìn)行互譜結(jié)合VMD去噪并根據(jù)綜合指標(biāo)選取的IMFs1和IMFs2頻域圖如圖16、圖17所示。

圖16 去噪信號IMFs1頻域圖Fig.16 IMFs1 frequency domain diagram of denoising signal

圖17 去噪信號IMFs2頻域圖Fig.17 IMFs2 frequency domain diagram of denoising signal

對IMF1和IMF2進(jìn)行FastICA計算后得到兩個獨(dú)立的信號源q1和q2如圖18、圖19所示。

圖18 信號源q1頻域圖Fig.18 Frequency domain diagram of signal source q1

圖19 信號源q2頻域圖Fig.19 Frequency domain diagram of signal source q2

同理計算管道下游采集的兩路信號s3和s4，經(jīng)過去噪處理和FastICA計算后得到信號源q3和q4，對q1，q2和q3，q4兩兩計算相關(guān)系數(shù)值如表1所示。

表1 信號源相關(guān)系數(shù)計算

由表1可知，(q1，q3)和(q2，q4)的相關(guān)系數(shù)較高，分別計算其相關(guān)延遲時間并得到最終的定位結(jié)果。在該工況下進(jìn)行10組試驗(yàn)得到的試驗(yàn)結(jié)果如圖20、圖21所示。

圖20 泄漏源位置為300 mm的定位結(jié)果Fig.20 Locate the leakage source at 300 mm

圖21 泄漏源位置為700 mm的定位結(jié)果Fig.21 Locate the leakage source at 700 mm

為衡量該方法的可靠性，對試驗(yàn)結(jié)果進(jìn)行數(shù)據(jù)分析，箱形圖可以用來觀察數(shù)據(jù)整體的分布情況，利用中位數(shù)，25%分位數(shù)，75%分位數(shù)，上邊界，下邊界等統(tǒng)計量來描述數(shù)據(jù)的整體分布情況。基于試驗(yàn)結(jié)果繪制箱型如圖22、圖23所示。

圖22 泄漏源位置為300 mm的定位統(tǒng)計結(jié)果Fig.22 Locate the leakage source at 300 mm

圖23 泄漏源位置為700 mm的定位統(tǒng)計結(jié)果Fig.23 Location statistics for leak source 700 mm

由圖22、圖23可知，泄漏源為300 mm的定位結(jié)果較為密集得分布在300 mm左右，誤差率在2%～15%，泄漏源為700 mm的定位結(jié)果存在一個偏離分布的異常值，可將其舍棄，總體誤差率在1%～10%?；谝陨辖Y(jié)果分析可得，該多點(diǎn)泄漏定位方法具有較高的可靠性和有效性。

5 結(jié) 論

本文通過對信號進(jìn)行變分模態(tài)分解后模態(tài)選取方法的研究，提出一種結(jié)合多指標(biāo)的模態(tài)選取方法，并根據(jù)互譜和變分模態(tài)分解各自的局限性將兩種算法結(jié)合應(yīng)用，利用該算法對采集信號進(jìn)行去噪預(yù)處理，將去噪后的信號通過FastICA實(shí)現(xiàn)源信號的分離，進(jìn)而通過時差定位實(shí)現(xiàn)管道的多點(diǎn)泄漏定位。

最后，通過實(shí)驗(yàn)室的管道泄漏模型開展了多點(diǎn)泄漏定位試驗(yàn)，通過試驗(yàn)結(jié)果的分析驗(yàn)證了該方法的有效性，為工業(yè)實(shí)況提供了一種管道多點(diǎn)泄漏的理論思路。