安連鎖,馮 強,沈國清,姜根山,張世平,王 鵬,周 鑫
(1.華北電力大學電站設備狀態(tài)監(jiān)測與控制教育部重點實驗室,北京102206;2.北京電力經(jīng)濟技術研究院,北京100055)
鍋爐壓力管道泄漏直接影響熱力發(fā)電站的安全經(jīng)濟運行,對壓力管道初期泄漏進行及時檢測能夠降低運行事故概率,并帶來一定的經(jīng)濟性.采用聲學法檢測泄漏信號雖然具有安裝和維護方便、非接觸式測量、實時性以及不需要外加能量等優(yōu)點,被廣泛應用于電站鍋爐、石油和天然氣檢查等方面,但是對于電站鍋爐爐膛復雜的聲環(huán)境,泄漏信號衰減較大,通過聲壓和頻譜分析很難檢測到初期的泄漏信號,且僅通過聲壓和頻譜分析檢測泄漏信號很容易造成誤判和漏判.若在早期發(fā)現(xiàn)電站鍋爐壓力管道泄漏,可及時上報,給電網(wǎng)調(diào)控預留出足夠時間,合理安排停爐.
經(jīng)驗模態(tài)分解(Empirical Mode Decomposition,EMD)[1-2]可對信號進行平穩(wěn)化處理,并具有自適應的優(yōu)點,能夠?qū)⑿盘柗纸鉃槿舾杀菊髂B(tài)函數(shù)(Intrinsic Mode Function,IMF).對包含特征信號的IMF進行研究能夠達到一定的降噪效果.在微弱信號檢測方面,由于Duffing振子對噪聲具有較顯著的免疫力,并且其對有用信號敏感的特性被廣泛應用于各領域.針對爐膛壓力管道初始微弱信號檢測問題,采用EMD 聯(lián)合自相關差分Duffing振子(CD-EMD)方法進行了研究.
美國NASA 的Nordon E Huang 工作組于1998年提出EMD 方法[3].該方法能夠?qū)⑷我庑盘柗纸鉃槿舾蒊MF和一個殘余項,因為依據(jù)自身時間尺度特征進行信號分解,所以不必像小波處理一樣預先設置基函數(shù).在對信號進行平穩(wěn)化處理的過程中,根據(jù)待檢測信號的特征,選取適當?shù)腎MF 能夠達到降噪的效果.從本質(zhì)上說,EMD 方法將時域信號按照頻率尺度逐級分解,產(chǎn)生若干IMF.
采用EMD 方法將信號分解為IMF 時需要滿足2個條件:一是對于整個數(shù)據(jù)序列,存在相等或者至多差1的極值點個數(shù)和過零點個數(shù);二是對于任意一點,信號局部極大值和局部極小值所定義的上包絡線和下包絡線的均值為0.
在滿足上述條件的前提下,EMD 方法可將信號x(t)篩選為若干IMF和一個殘余項,具體步驟如下:
(1)確定數(shù)據(jù)序列的所有極大值點,利用三次樣條插值擬合出上包絡線u(t),采取同樣方法對所有極小值點擬合出下包絡線d(t),并記上包絡線和下包絡線的均值為m(t),即m(t)=[d(t)+u(t)]/2.
(2)定義h1(t)=x(t)-m(t),對于非線性信號和非平穩(wěn)信號,通常h1(t)不滿足IMF 所要求的2個條件,對于不滿足條件的h1(t)需要重復執(zhí)行步驟(1)和步驟(2),直至h1(t)符合IMF 所要求的條件為止,將滿足條件的h1(t)記作h1imf(t).
(3)將原數(shù)據(jù)序列減去h1imf(t),得到剩余數(shù)據(jù)序列r1(t):
(4)將剩余數(shù)據(jù)序列r1(t)作為原始數(shù)據(jù)序列,重復步驟(1)~步驟(3),依次分解得到:
篩選過程的停止準則[4]為:限制兩連續(xù)處理結(jié)果的標準差Sd,通常Sd的取值為0.2~0.3.Sd的表達式如下:
式中:T為數(shù)據(jù)序列時間尺度;hi(k-1)(t)和hik(t)為IMF篩選過程中兩連續(xù)處理結(jié)果的數(shù)據(jù)序列;k為第i個IMF篩選過程中的次數(shù)序列.
信號經(jīng)過上述步驟被分解為n個IMFhiimf(t)和一個殘余函數(shù)rn(t),即
其中,各IMF分量包含了不同時間尺度的特征信號,殘余函數(shù)表示了原數(shù)據(jù)趨勢量信息,如圖1所示.
圖1 電站鍋爐實測噪聲信號的EMD方法分解示意圖Fig.1 EMD method diagram of background noise in power plant boiler
Holmes型Duffing系統(tǒng)方程為
式中:k為阻尼系數(shù);fcos(ωt)為周期項;-x+x3為非線性項,非線性項部分決定了Duffing振子具有非線性動力學特性.
先前的學者在研究Duffing振子時通常采用小頻率參數(shù)來檢測低頻信號.當檢測高頻信號時需要增大阻尼系數(shù)k,并且同時調(diào)整周期項系數(shù)f,然而實際運用中k應取較小值,因為當k較大時,含噪聲待檢測信號輸入Duffing振子不足以使系統(tǒng)狀態(tài)躍變,即不存在混沌狀態(tài)到大尺度周期狀態(tài)的相變[5].對于不同的待檢測信號,需要根據(jù)大量實驗數(shù)據(jù)調(diào)整臨界參數(shù)k和f,其過程較為繁瑣.針對上述問題,采用賴志慧等[6]提出的變尺度檢測方法,該方法對任意頻率的信號在時間尺度上進行變換,然后將變換后的信號輸入到Duffing振子,其核心在于將任意頻率的信號轉(zhuǎn)換為角頻率為1rad/s的信號.
構造Duffing振子檢測系統(tǒng):
式中:Ⅰ(t)為輸入信號,包含了待檢測周期信號Acos(ωt)和噪聲n(t),即Ⅰ(t)=Acos(ωt)+n(t).
由于Duffing振子對噪聲具有免疫特性[7],因此對于不含待檢測信號的噪聲,系統(tǒng)不會發(fā)生狀態(tài)躍變,仍處于混沌狀態(tài),但當輸入信號含有微小幅值的待檢測信號時,系統(tǒng)會發(fā)生狀態(tài)躍變,進入大尺度周期狀態(tài),如圖2所示.對于電站鍋爐壓力管道的微弱泄漏信號,正是應用Duffing振子對特定頻率信號的敏感性來檢測的.
圖2 Duffing振子檢測系統(tǒng)相變圖Fig.2 Phase diagram of Duffing oscillator detection system
在上述分析中,將待檢測信號初始相位均假設為0,然而實際工程應用中,初始相位為0的情況幾乎不存在.由文獻[6]可知,由于待檢測信號初始相位的不同而存在一定的檢測誤區(qū),使得Duffing振子微弱泄漏信號檢測在工程應用中受到限制.
差分法采用的是雙振子模式,構造方程如下:
式中:Ⅰ(t)=Acos(ωt+φ)+n(t),φ∈[0,2π];ξ通常為1~1.5內(nèi)的數(shù)值,若ξ=1,則式(7)等同于一維Duffing振子.
ξ決定了差分后共模成分和混沌成分的凸顯,當其值趨于1時,兩振子差分后共模噪聲得到抑制,但經(jīng)驗閾值越小,影響混沌狀態(tài)判別的準確度;而當其值趨于1.5時,經(jīng)驗閾值增大,Duffing振子輸出相位差增大,分析難度增加,在本文中ξ取1.004.
對Ⅰ(t)進行自相關運算,則
式中:v為自相關時延值;R為自相關輸出值;下標xx表示實際接受信號,ss表示純信號,nn表示噪聲信號.
對于爐膛背景噪聲n(t),若v較大,Rxx(v)表現(xiàn)為Rss(v)特征.那么對于待檢測信號Acos(ωt+φ),則有
通過式(9)可知,輸入信號經(jīng)過自相關運算后,可消除待檢測信號相位差時干擾.
小波理論是在傅里葉變換基礎上發(fā)展而來的,由于其對時域及頻域都具有局部特征表征能力,因此廣泛應用于各個行業(yè)[8-9].進行小波變換時首先要選取小波函數(shù),小波函數(shù)必須滿足:
式中:φ(ω)為φ(t)的傅里葉變換,且φ(ω)∈L2(R),φ(ω)ω=0=0;L2(R)表示二次可積函數(shù)組成的空間,φ(t)為小波基函數(shù).
對小波基函數(shù)進行伸縮、平移,可以得到小波基函數(shù)φa,τ(t):
式中:a為尺度因子,a∈R,a>0;τ為平移因子,τ∈R.
實際工程應用中采集到的均為離散信號,需要進行尺度因子和平移因子的二進制變換,取a=2i,τ=2ij,j∈Z,帶入式(11),對應的離散小波基函數(shù)為
對于任意L2(R)空間,在小波基函數(shù)φa,τ(t)下展開函數(shù)f(t),稱為函數(shù)f(t)的連續(xù)小波變換,其表達式為
式中:φ*((t-τ)/a)為φ((t-τ)/a)的共軛復數(shù).
函數(shù)f(t)經(jīng)過連續(xù)小波變換后,可由下式重構恢復原信號,即進行小波逆變換:
式中:Cφ為小波基函數(shù)φa,τ(t)的可容許條件.
電站鍋爐壓力管道微弱泄漏信號檢測系統(tǒng)由聲感知設備、聲波導管、信號調(diào)理器和功率放大器等設備組成,其中聲波導管起到隔離燃燒和隔熱的作用,對聲感知設備起到保護作用.聲感知設備的布置方式需要結(jié)合爐膛復雜的實際環(huán)境進行調(diào)整,在此不討論,僅給出一般性示例,如圖3所示.
圖3 電站鍋爐立體圖Fig.3 Power plant boiler stereogram
聲感知設備設定為每隔一段時間采集一次,并將采集數(shù)據(jù)與先前一次的采集數(shù)據(jù)進行分析對比.對采集到的數(shù)據(jù)先進行EMD 分解,對包含泄漏信號特征量的IMF進行頻譜分析,同時將對應的IMF經(jīng)自相關處理后輸入至差分Duffing 振子檢測系統(tǒng),通過經(jīng)驗閾值進行微弱泄漏信號判別.
為評估EMD、小波、EMD 聯(lián)合Duffing 振子(D-EMD)和CD-EMD 方法的性能,在待檢測信號中加入現(xiàn)場采集的爐膛背景噪聲,爐膛背景噪聲采樣頻率設定為102 400Hz,采樣個數(shù)為65 536.采用不同信噪比對各方法分別進行信號檢測仿真實驗.信噪比計算公式采用,其中為待檢測信號的均方值,為噪聲信號的均方值.
為了從時域及頻域圖上顯示各方法檢測效果,待檢測聲源選取0.6 MPa間歇射流氣動聲源(以下簡稱0.6 MPa氣動聲源).由于高頻信號在爐膛實際高溫粉塵環(huán)境中衰減較嚴重,因此選取10kHz以下的4個峰值頻率(2 583Hz、3 588Hz、6 309Hz和9 152Hz)為信號特征頻率,如圖4所示.
圖4 0.6 MPa氣動聲源頻譜圖Fig.4 Frequency spectrum of the 0.6 MPa pneumatic sound
將0.6 MPa氣動聲源信號加入爐膛背景噪聲,當信噪比為-5dB 時,從時域圖上EMD 方法和小波分解均能較好地辨別氣動聲源信號,但實際泄漏信號為非間歇射流模型,因此在時域圖上進行判別存在技術難題.在頻域方面,EMD 方法可以較好地檢測到信號特征頻率,如圖5所示(其中FFT 表示傅里葉變換).小波分解采用db5小波基進行5層分解,能夠檢測到2 583Hz、3 588Hz和6 309 Hz頻率,但出現(xiàn)主峰值7 645 Hz頻率干擾,且無法有效檢測到9 152Hz頻率,如圖6所示.
圖5 信噪比為-5dB時EMD方法分解高頻重構時域圖及對應的FFT 頻譜圖Fig.5 High-frequency reconstruction time domain and FFT frequency domain by EMD method at-5dB signal to noise ratio
圖6 信噪比為-5dB時小波分解高頻重構時域圖及對應的FFT 頻譜圖Fig.6 high-frequency reconstruction time domain and FFT frequency domain by wavelet decomposition at-5dB signal to noise ratio
當信噪比為-14.5dB 時,在時域方面,EMD方法及小波分解均能較好地辨別氣動聲源信號;在頻域方面,EMD 方法分解后的峰值頻率仍為特征頻率,但出現(xiàn)旁瓣干擾,如圖7所示.小波分解出現(xiàn)嚴重的主頻干擾,且2 583Hz畸變?yōu)? 547Hz,如圖8所示,因此在-14.5dB信噪比環(huán)境下,小波分解不適用于爐膛壓力管道微弱泄漏信號檢測.
圖7 信噪比為-14.5dB時EMD方法分解高頻重構時域圖及對應的FFT 頻域圖Fig.7 High-frequency reconstruction time domain and FFT frequency domain by EMD method at-14.5dB signal to noise ratio
當信噪比為-24dB時,EMD 方法分解已無法正常識別特征頻率,如圖9所示,因此引入D-EMD方法,并改進為CD-EMD 方法,將2種方法進行對比分析.在無初始相位干擾情況下,信噪比為-24 dB時,D-EMD 和CD-EMD方法的輸出相圖為大尺度周期狀態(tài),如圖10和圖11所示.
圖8 信噪比為-14.5dB時小波分解含噪聲待檢測信號時域圖及FFT 頻譜圖Fig.8 Time domain and frequency domain of noisy signals by wavelet decomposition at-14.5dB signal to noise ratio
圖9 信噪比為-24dB時EMD方法分解含噪聲待檢測信號時域圖及FFT 頻譜圖Fig.9 Time domain and frequency domain of noisy signals by EMD method at-24dB signal to noise ratio
圖10 D-EMD方法輸出相圖Fig.10 D-EMD method output phase diagram
圖11 CD-EMD方法輸出相圖Fig.11 CD-EMD method output phase diagram
5.2.1 初始相位的干擾
根據(jù)文獻[8]的結(jié)論,當φ∈[0,π/3]∪[5π/3,2π]時,Duffing振子可以對微弱周期信號進行有效識別,因此選取不同待檢測信號初始相位(π/6、5π/6、7π/6 和11π/6),分別采用D-EMD 和CD-EMD方法進行仿真實驗.
圖12和圖13分別給出了D-EMD 和CD-EMD方法不同初始相位的輸出相圖.表1給出了2種方法不同初始相位狀態(tài)的對比.由圖12、圖13和表1可以看出,D-EMD 方法受待檢測信號初始相位的干擾嚴重,容易出現(xiàn)誤判;而CD-EMD方法通過引入自相關方法,能夠?qū)⑷我獬跏枷辔淮龣z測信號轉(zhuǎn)換為0相位,避免了初始相位的干擾,消除了由初始相位導致的特征頻率誤判.
圖12 D-EMD方法不同初始相位的輸出相圖Fig.12 Output phase diagram of D-EMD method with different initial phases
圖13 CD-EMD方法不同初始相位的輸出相圖Fig.13 Output phase diagram of CD-EMD method with different initial phases
表1 初始相位的影響對比Tab.1 Comparison of influence among different initial phases
5.2.2 CD-EMD方法特征頻率的判別
CD-EMD 方法對于混沌狀態(tài)判別的優(yōu)點在于可以通過設定經(jīng)驗閾值進行判別,而且有效地避開了lyapunov等數(shù)值計算方法,提高了運算效率,增強了系統(tǒng)運行的實時性.當輸入信號不包含特征頻率時,差分波形圖如圖14(a)所示,幅值在2.5Pa以上;當輸入信號包含特征頻率時,差分波形圖如圖14(b)所示,幅值在0.1Pa以下.在實際工程應用時可以通過設定經(jīng)驗閾值的方式,當差分幅值高于經(jīng)驗閾值時即為混沌狀態(tài),輸入信號不包含特征頻率;當差分幅值低于經(jīng)驗閾值時即為大尺度周期狀態(tài),輸入信號包含特征頻率,若多個CD-EMD陣列同時檢測到特征頻率,即可判定爐膛壓力管道出現(xiàn)微弱泄漏信號.
圖14 差分波形圖Fig.14 Difference waveform
5.2.3 實際運行中特征頻率的選取
本研究針對泄漏初期微弱泄漏信號,高頻泄漏信號在原始未經(jīng)處理的頻譜圖上難以識別,而經(jīng)EMD 方法分解后,爐膛背景噪聲也存在少量高頻信號,但是總體上與泄漏信號存在差別.
對1mm、2mm、3mm、5mm 和9mm 口徑的噴嘴在8個大氣壓下使用空氣壓縮機模擬泄漏信號,采用前文提到的控制信噪比的方式,使信噪比為-24dB,將不同口徑的泄漏信號加入爐膛背景噪聲,EMD 方法分解結(jié)果如圖15所示.
對比各口徑噴嘴EMD 方法分解輸出的2imf圖形可知,爐膛背景噪聲輸出圖形在主峰值右側(cè)頻率,頻率分布驟降,總體分布呈三角分布形式;含噪聲待檢測信號輸出圖形在主峰值右側(cè)頻率,頻率有一段平坦部分,總體分布為梯形分布.判別泄漏信號特征頻率從主峰值右側(cè)頻率分布入手,從主峰值右側(cè)尋幾組次峰值(以下簡稱右次峰值),列出差分Duffing振子陣列,爐膛背景噪聲右次峰值不足以驅(qū)動差分Duffing振子,而含噪聲待檢測信號能夠驅(qū)動差分Duffing振子發(fā)生相變.
圖15 不同口徑噴嘴的EMD方法分解圖Fig.15 EMD method decomposition of differently sized nozzles
(1)通過控制信噪比的方式,選取0.6 MPa氣動聲源加入爐膛背景噪聲,當信噪比為-5dB 時,EMD 方法可以較好地檢測到信號特征頻率,小波分解出現(xiàn)主峰值7 645 Hz頻率干擾,且9 152 Hz頻率無法得到有效檢測.當信噪比為-14.5dB 時,EMD 方法分解后峰值頻率仍為特征頻率,但特征頻率的最大幅值與爐膛背景噪聲EMD 方法分解后的最大幅值無明顯差別,小波分解則出現(xiàn)嚴重干擾.EMD 方法在信噪比大于-14.5dB 時的檢測效果明顯優(yōu)于小波分解.
(2)在D-EMD方法的基礎上加以改進后提出CD-EMD 方法,將信噪比檢測下限降低至-24dB,并且通過各象限取不同初始相位進行對比研究,DEMD 方法出現(xiàn)嚴重的初始相位干擾,不適用于實際工程應用,而CD-EMD方法能夠有效地克服初始相位引起的誤判.
(3)CD-EMD 方法可以通過經(jīng)驗閾值的方式判別有無特征頻率信號,而且有效地避開了lyapunov等數(shù)值計算方法,提高了系統(tǒng)實際運行的實時性及判別準確度.
[1]邵忍平,曹精明,李永龍.基于EMD 小波閾值去噪和時頻分析的齒輪故障模式識別與診斷[J].振動與沖擊,2012,31(8):96-101.
SHAO Renping,CAO Jingming,LI Yonglong.Gear fault pattern identification and diagnosis using timefrequency analysis and wavelet threshold de-noising based on EMD[J].Journal of Vibration and Shock,2012,31(8):96-101.
[2]王文波,張曉東,汪祥莉.脈沖星信號的經(jīng)驗模態(tài)分解模態(tài)單元比例萎縮消噪算法[J].物理學報,2013,62(6):534-542.
WANG Wenbo,ZHANG Xiaodong,WANG Xiangli.Pulsar signal denoising method based on empirical mode decomposition mode cell proportion shrinking[J].Acta Physica Sinica,2013,62(6):534-542.
[3]HUANG N E,SHEN Z,LONG S R,etal.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J].Proceedings:Mathematical,Physical and Engineering Sciences,1998,454:903-995.
[4]RILLING G,F(xiàn)LANDRIN P,GON?ALVèS P.On empirical mode decomposition and its algorithms[C]//In Proceedings of the 6th IEEE/EURASIP Workshop on Nonlinear Signal and Image Processing.Italy:NSIP,2003.
[5]楊紅英,葛傳虎,葉昊,等.基于Duffing振子的天然氣管道泄漏檢測方法[J].高技術通訊,2010,20(6):628-631.
YANG Hongying,GE Chuanhu,YE Hao,etal.Leak detection based on Duffing oscillators for gas pipelines[J].High Technology Letters,2010,20(6):628-631.
[6]賴志慧,冷永剛,孫建橋,等.基于Duffing振子的變尺度微弱特征信號檢測方法研究[J].物理學報,2012,61(5):60-68.
LAI Zhihui,LENG Yonggang,SUN Jianqiao,etal.Weak characteristic signal detection based on scale transformation of Duffing oscillator[J].Acta Physica Sinica,2012,61(5):60-68.
[7]許雪梅,戴鵬,楊兵初,等.光聲池中微弱光聲信號檢測[J].物理學報,2013,62(20):1-9.
XU Xuemei,DAI Peng,YANG Bingchu,etal.Weak photoacoustic signal detection in photoacoustic cell[J].Acta Physica Sinica,2013,62(20):1-9.
[8]胡三高,安宏文,馬志勇,等.基于小波奇異值分析的汽輪機碰磨特征提?。跩].動力工程學報,2013,33(3):184-188.
HU Sangao,AN Hongwen,MA Zhiyong,etal.Feature extraction of rubbing fault for steam turbines based on wavelet singularity analysis[J].Journal of Chinese Society of Power Engineering,2013,33(3):184-188.
[9]許小剛,王松嶺,劉錦廉.基于小波包能量分析及改進支持向量機的風機機械故障診斷[J].動力工程學報,2013,33(8):606-612.
XU Xiaogang,WANG Songling,LIU Jinlian.Mechanical fault diagnosis of fan based on wavelet packet energy analysis and improved support vector machine[J].Journal of Chinese Society of Power Engineering,2013,33(8):606-612.