基于EMD和SST算法的閃電電場信號去噪研究

(1.西北師范大學(xué) 物理與電子工程學(xué)院，甘肅 蘭州 730070； 2.甘肅省智能信息技術(shù)與應(yīng)用工程研究中心，甘肅 蘭州 730070； 3.國網(wǎng)蘭州供電公司，甘肅 蘭州 730070)

在閃電電場信號的采集、傳輸過程中不可避免地要受到噪聲的干擾[1]，如環(huán)境噪聲、儀器噪聲、信道噪聲等，這些噪聲的存在會直接影響到對一些微小閃電放電過程的認(rèn)識，也會給后期的數(shù)據(jù)處理和分析帶來不便，所以對閃電電場信號的去噪處理成為了亟待解決的問題。

目前，閃電信號的去噪主要依賴于硬件濾波電路，但這種技術(shù)存在局限性，去噪效果并不理想，因此對閃電信號去噪算法的研究越來越受到人們的重視。李鵬等人應(yīng)用傳統(tǒng)數(shù)字濾波和小波閾值法進(jìn)行了閃電信號去噪的比較研究[2]。高太長等人針對回?fù)暨^程中電場變化情況，提出了一種基于多小波變換的閃電信號去噪方法[3]。劉志田等人提出了基于提升小波變換的閃電電場變化信號去噪方法[4]。火元蓮等人提出了一種基于雙密度雙樹小波變換的閃電信號去噪方法[5]。近年來，由Huang[6]等人提出的經(jīng)驗?zāi)B(tài)分解算(Empirical Mode Decomposion,EMD)和由Daubechies[7]等人提出的同步壓縮小波變換算法(Synchrosqueezing Wavelet Transform,SST)在非平穩(wěn)信號的處理中得到了蓬勃發(fā)展。EMD主要依據(jù)原始時間序列自身的時間尺度特征分析序列，而不必預(yù)先設(shè)定任何基函數(shù)[8]。其主要思想是它能將復(fù)雜信號分解為若干個本征模態(tài)函數(shù)(Intrinsic Mode Function,IMF)，由于分解是基于信號序列時間尺度的局部特性[9]，因此具有自適應(yīng)性，相比于短時傅里葉變換、小波分解等方法其在處理非平穩(wěn)信號上有明顯優(yōu)勢，所以一經(jīng)提出就在不同工程領(lǐng)域得到了廣泛的應(yīng)用[10-12]。SST是基于小波變換提出的一種新的時頻分析方法，可將一維信號變換到二維空間，即可將時域信號轉(zhuǎn)化成高分辨率的時頻譜，再結(jié)合時頻譜重排的思想，通過壓縮任一中心頻率附近區(qū)間值，從而得到同步壓縮小波變換量值[13]。本文將EMD算法和SST算法結(jié)合起來，利用EMD算法能夠自適應(yīng)分解信號和SST算法可將噪聲壓縮為點狀噪聲或顆粒狀噪聲并集中分布的特點[14]，從而選用中值濾波可達(dá)到對噪聲的抑制[15-16]，實現(xiàn)對閃電電場信號的去噪。

1 經(jīng)驗?zāi)B(tài)分解和同步壓縮小波變換

1.1 經(jīng)驗?zāi)B(tài)分解

經(jīng)驗?zāi)B(tài)分解可將待分析信號自適應(yīng)地分解為若干個不同時間尺度的IMF分量，對信號x(t)進(jìn)行EMD分解的過程簡介如下[6]。

首先，分別采用三次樣條插值函數(shù)來擬合原始信號x(t)的上下包絡(luò)線xmax(t)、xmin(t)，通過對得到的上下包絡(luò)線求均值得到原始信號的平均包絡(luò)線m1(t)

m1(t)=(xmax(t)+xmin(t))/2

(1)

用x(t)減去m1(t)得到h1(t)

h1(t)=x(t)-m1(t)

(2)

通過判斷得到的h1(t) 分量是否滿足構(gòu)成IMF分量的2個條件，來決定其是否為第一階IMF分量。一般情況下，第一次分解得到的h1(t)信號分量很難滿足判定條件，因此需要將h1(t)信號分量作為新的原始信號，重復(fù)上述步驟，即可以得到

h11(t)=h1(t)-m11(t)

(3)

式中，m11(t)為h1(t)的上下包絡(luò)線的均值。繼續(xù)對h11(t)進(jìn)行循環(huán)判斷和篩選，重復(fù)以上的步驟。若上述過程循環(huán)了k次，記為h1k(t)

h1k(t)=h1(k-1)(t)-m1k(t)

(4)

當(dāng)h1k(t)滿足結(jié)束條件時則終止篩選。

1.2 同步壓縮小波變換

同步壓縮小波變換[17-18]是在連續(xù)小波變換(CWT)的基礎(chǔ)上提出的一種新的時頻分析方法[19]。通過對CWT結(jié)果進(jìn)行壓縮，降低了噪聲的能量，一定程度上壓制了噪聲，且SST也是一種可逆變換，即通過反變換可完全恢復(fù)原始信號，其過程如下。

首先對時域信號f(t)進(jìn)行連續(xù)小波變換得到小波系數(shù)wf(a,b)

(5)

式中，a為尺度因子;b為時間因子;φ*為母小波的共軛。

令諧波信號f(t)=Acos(wt)，根據(jù)Plancherel定理，式(5)的頻率域變換為

(6)

將f(t)的傅里葉變換代入式(6)得

(7)

通過對小波系數(shù)求偏導(dǎo)可估計瞬時頻率

(8)

通過式(8)，將小波系數(shù)wf(a,b)從時間-尺度平面映射到時間-頻率平面，即wf?wf(a,b),b」。

(9)

式中，ak為離散尺度，k為尺度個數(shù)。

同步壓縮小波變換的反變換為

(10)

式中，c-1φ表示取實部。

2 仿真信號的去噪分析

2.1 組合去噪方法流程

所提出的組合算法利用了EMD算法能夠自適應(yīng)分解信號和SST算法可將噪聲壓縮為點狀噪聲或顆粒狀噪聲并集中分布的特點，從而利用中值濾波法可以實現(xiàn)對閃電電場信號的去噪處理。其具體步驟如下：

① 對含噪信號進(jìn)行EMD自適應(yīng)分解，得到各階IMF分量；

② 利用原始信號與各階IMF分量的相關(guān)系數(shù)大小，挑選出IMF的優(yōu)勢分量；

③ 利用IMF優(yōu)勢分量重構(gòu)出信號；

④ 將重構(gòu)后的信號進(jìn)行SST變換；

⑤ 對SST變換時頻圖進(jìn)行中值濾波；

⑥ 用SST逆變換重構(gòu)出去噪后的信號。

算法流程圖如圖1所示。

從圖1可知，本算法的關(guān)鍵流程是對信號進(jìn)行EMD分解和SST變換。

圖1 去噪流程圖

2.2 仿真信號實驗結(jié)果

根據(jù)M.W.Wik等人的研究，可以把閃電電磁脈沖(LEMP)歸結(jié)為雙指數(shù)衰減型脈沖波形。因此，標(biāo)準(zhǔn)閃電波的時域波形可表示為[20]

E(t)=E0(e-αt-e-βt)

(11)

式中，E0為脈沖波形的幅值系數(shù)；α、β為波前衰減系數(shù)和波尾衰減系數(shù)。在仿真過程中E0=30 V/m，α=2.0×107s-1，β=2.0×106s-1，采樣頻率fs=60 MHz。在雙指數(shù)衰減脈沖(如圖2(a))上疊加一信噪比(SNR)為30 dB 的高斯白噪聲(如圖2(b))。分別利用小波閾值法(小波函數(shù)取db5小波，分解層數(shù)取為6)、EMD算法、SST算法和所提出的組合方法對含噪信號進(jìn)行消噪處理。在組合算法中，先用EMD算法自適應(yīng)地將含噪信號分解為7個IMF分量，計算出各階IMF分量與原始信號的相關(guān)系數(shù)大小，選擇相關(guān)系數(shù)大于0.3的各階IMF分量重構(gòu)出信號(EMD算法去噪結(jié)果)，然后將重構(gòu)信號進(jìn)行SST變換，將變換后的時頻圖進(jìn)行中值濾波，之后進(jìn)行SST逆變換，逆變換后的結(jié)果即為本文去噪方法的去噪結(jié)果。去噪結(jié)果如圖3所示，圖中依次為小波閾值去噪、EMD去噪、SST去噪、組合方法去噪的結(jié)果。

圖2 原始信號和含噪信號

圖3 4種方法去噪效果圖

為了定量地比較組合方法和另外3種方法的去噪效果，分別計算了信噪比(SNR)、相關(guān)系數(shù)(CC)和均方誤差(MSE)。這3種參數(shù)均能從不同角度來度量去噪效果，其中信噪比和相關(guān)系數(shù)越大，均方誤差越小，說明去噪效果越好。計算結(jié)果如表1所示，結(jié)果表明基于EMD和SST算法的信噪比和相關(guān)系數(shù)比其余3種方法要大，均方誤差較另外3種方法要小，說明組合去噪方法去噪效果優(yōu)于其他3種算法的去噪效果。

表1 4種去噪方法的濾波結(jié)果比較

3 實測閃電電場信號的去噪分析

3.1 數(shù)據(jù)來源及實驗結(jié)果

本文選用的實驗數(shù)據(jù)是2009年在青海大通地區(qū)用快電場變化測量儀記錄的觀測資料。大通地區(qū)海拔較高，屬于雷電多發(fā)區(qū)，每個測站安裝有用于閃電輻射脈沖三維定位的閃電VHF輻射源到達(dá)時間差(TOA)定位系統(tǒng)和GPS同步的高精度時鐘(±25 ns) ，快電場變化探測儀帶寬為100 Hz～5 MHz，時間常數(shù)為1 ms、采樣速率為2.5 MHz。

在對閃電電場信號進(jìn)行去噪之前，先進(jìn)行去均值和歸一化處理。去均值即對閃電電場信號進(jìn)行零均值化處理。由于采集到的閃電電場信號的距離大小不同，閃電放電強度不同，為了減少數(shù)據(jù)的分散性，對閃電電場信號進(jìn)行歸一化，以得到幅值范圍統(tǒng)一的信號。考慮到EMD算法需要在整個信號長度范圍內(nèi)作樣條插值，當(dāng)信號采樣點數(shù)較多時，特別是極值點數(shù)目多的情況下，采用三次樣條擬合法計算量特別大，而采集到的一維閃電電場信號記錄時間超過800 ms，包含了2096000個采樣點，因此，直接對閃電信號進(jìn)行EMD分解是很難實現(xiàn)的。為了提高EMD算法的計算速率，首先對閃電電場信號以1∶1000進(jìn)行了重采樣，對重采樣后的信號分別用小波閾值法、EMD算法、SST算法和組合算法對地閃個例信號(如圖4所示)進(jìn)行了去噪處理。

圖4 地閃輻射場信號

在所提組合算法中，首先利用EMD算法將經(jīng)過去均值和歸一化后的地閃輻射場信號(如圖4)自適應(yīng)地分解為13個IMF分量，計算出原地閃輻射場信號與各階IMF分量相關(guān)系數(shù)的大小，挑選出相關(guān)系數(shù)大于0.2的各階IMF分量，利用這些優(yōu)勢分量重構(gòu)出地閃信號，然后對重構(gòu)后的信號進(jìn)行SST變換，將變換后的SST時頻面進(jìn)行中值濾波，之后再進(jìn)行SST逆變換，逆變換后的結(jié)果即為所提算法的去噪結(jié)果。4種方法去噪效果如圖5所示。

圖5 4種方法去噪效果圖

從圖5中可以看出，單獨用EMD算法和單獨用SST算法去噪效果并不明顯，信號還存在著大量噪聲，這會影響到對閃電信號的時頻分析和特征提取，進(jìn)而影響對閃電電場信號的自動化識別。小波閾值去噪的結(jié)果要比另外3種方法去噪效果平滑，這是因為小波閾值法在去除噪聲的同時，也將信號中的細(xì)節(jié)部分平滑了，由于閃電信號的突變和尖峰處包含了大量信息，該方法在去噪的同時也將信號中的有用成分濾除掉了，因此這種方法對后期的數(shù)據(jù)處理是很不利的，而基于EMD和SST算法在有效去噪的同時也保留了信號的細(xì)節(jié)部分，便于對閃電微小放電過程如回?fù)舻鹊恼J(rèn)識，也便于后期數(shù)據(jù)處理。因此，組合算法在對閃電電場信號的去噪處理方面具有優(yōu)越性。

3.2 實驗結(jié)果定量分析

分別選取5例云閃輻射場信號(IC)和5例地閃輻射場信號(CG)，計算在4種去噪方法下的相關(guān)系數(shù)、信噪比、均方誤差。其中，相關(guān)系數(shù)計算結(jié)果如圖6所示，信噪比計算結(jié)果如表2所示，均方誤差計算結(jié)果如圖7所示。相關(guān)系數(shù)和信噪比越大，均方誤差越小說明去噪效果越好。

圖6 4種去噪方法下的相關(guān)系數(shù)

閃電信號小波閾值法EMD算法SST算法本文算法GC115.123.420.226.8GC217.320.419.223.5GC314.716.915.421.3GC419.521.620.126.9GC520.123.521.228.2IC119.123.421.626.3IC220.423.522.928.4IC313.615.914.820.3IC49.413.511.218.9IC52.93.13.08.4

圖7 4種去噪方法下的均方誤差

觀察這3幅圖表可發(fā)現(xiàn)，最大的信噪比和相關(guān)系數(shù)、最小的均方誤差均是本文所提方法計算出的。因此，本文算法的去噪效果要優(yōu)于文中另外3種方法的去噪效果，為閃電電場信號的去噪提供新思路。

4 結(jié)束語

本文著眼于處理非平穩(wěn)信號的EMD算法和SST算法，將其結(jié)合起來發(fā)揮其各自的優(yōu)勢，并運用于閃電電場信號的去噪研究，實驗結(jié)果表明，相比于傳統(tǒng)小波閾值去噪和單獨使用EMD算法和單獨使用SST算法去噪，本文所提算法的去噪效果更好，這對閃電電場信號的后期處理，例如特征提取和分類識別具有重要的意義。