基于hilmann變換的非平穩地震動信號分析

在強地震記錄的地面運動中，地震波是不穩定的，包括地表土層和其他介質的非線性信息。在傳統的不規則地震記錄處理方法中，frr可以在頻域中獲得非常高的分辨率，但fower譜不能反映信號的瞬態信息。因此，時域失去了區分能力。短期frr可以在一定程度上描述信號的瞬時頻率含量，但不受不確定性原則的限制，因此無法同時獲得高的時間和頻率范圍。小波變換是一種多分辨率信號治理方法，可以同時在時間和頻率范圍內獲得較高的分辨率，但其分辨率仍然受到一定的限制。基于Hilbert變換的Hilbert譜分析,在處理某些特殊的非平穩信號時,能夠獲得非常好的結果;但是,Hilbert譜分析所面臨的問題是:對于大部分非平穩的信號,Hilbert譜分析失去了本來的物理意義.為此,1998年Huang提出了一種新的信號處理技術,其基本的思路是:為了將Hilbert變換應用到任意非平穩信號,首先利用Huang提出的經驗模態分解(empiricalmodedecomposition,EMD)方法,將給定信號分解成若干個本征模函數(intrinsicmodefunction,IMF);然后,再對每個IMF進行Hilbert譜分析(Hilbertspectralanalysis,HSA),得到每個IMF的Hilbert譜;最后,匯總所有IMF的Hilbert譜,得到原始非平穩信號的Hilbert譜.這種新的信號處理方法被稱為Hilbert-Huangtransform,簡稱HHT.在EMD方法中,信號兩端的邊界效應所帶來的誤差會向內傳播,進而“污染”整個數據序列,使得最后的結果失去意義,尤其對于低頻的IMF分量來說,這種邊界效應所引起的誤差更加嚴重.Huang針對這個問題,提出用“特征波”對原始信號進行延拓的方法,但并未公開具體的處理方法,并且已將該方法在美國申請了專利.此外,Huang指出,EMD方法所面臨的邊界延拓問題還沒有完全解決.因此,解決邊界延拓問題對于HHT理論、以致對其他信號處理,都具有理論和實際意義.鄧擁軍等提出了使用神經網絡方法對原始信號進行延拓來解決EMD中的邊界問題.本文針對同一問題,應用傳統的自回歸(auto-regressive,AR)模型對原始數據進行延拓,提出了“邊篩分,邊延拓”的邊界處理方法,計算結果表明,該方法對于EMD的準確分解非常有效,對于低頻IMF分量也得到了很好的結果.1hhd的信號處理在EMD方法中,Huang將具有如下性質的信號定義為本征模函數:(1)該信號的極值點的數目與零交點的數目相等或至多相差一個;(2)該信號的極大值點與極小值點關于零軸對稱.Huang認為Hilbert變換作用到具有上述性質的信號上時,便能給出物理意義明確的結果來,即能將IMF中所蘊涵的波內調解機制(intra-wavemodulation)提取出來,所得到的Hilbert譜與小波譜以及其他基于Fourier變換時變頻譜相比,具有非常高的時頻分辨率.EMD方法通過一種被Huang成為“篩分”(sifting)的過程,對數據逐步進行分解,最后得到一系列IMF分量.具體處理方法是:給定實信號x(t),找出x(t)所有的極大值點并將其用三次樣條函數擬合成原始信號的上包絡線;找出x(t)所有的極小值點并將其用三次樣條函數擬合成原始信號的下包絡線;上下包絡線的均值為原始信號的平均包絡線m1(t);將原始信號x(t)減去m1(t)后即可得到一個新的信號h1(t):這個過程稱為篩分,原始信號x(t)經過一次篩分后變為h1(t).一般說來,h1(t)仍然不是一個IMF,為此需要對它重復上述篩分處理.重復處理k次后,若最后所得到的信號滿足IMF條件,就得到了原始信號的第1個IMF分量c1(t):c1(t)代表原始信號中最高頻的IMF分量.將原始信號x(t)減去c1(t)就可以得到去除高頻成分的殘余信號r1(t);對r1(t)進行上述篩分處理后可以得到第2個IMF分量c2(t),然后將r1(t)減去c2(t)后可以得到r3(t);如此重復下去直到最后一個殘余信號rn(t)不可再分解為止.最后,原始信號x(t)可以表示成通常,原始信號x(t)經過上述EMD方法分解后所得到的IMF分量的數目是很少的,這說明EMD方法的效率是很高的.原始信號x(t)經過EMD分解后,分別對每個IMF分量進行Hilbert譜分析(HSA)得到每個IMF分量的Hilbert譜,匯總所有IMF分量的Hilbert譜就得到了原始信號x(t)的Hilbert譜.即給定任意一個IMF,c(t),其Hilbert變換定義為其中P表示Cauchy主值.c(t)與y(t)可以合成解析信號z(t):從而,可以定義時變的幅值a(t)和相位θ(t)如下:進一步可以定義瞬時頻率ω(t):時變幅值a(t)的時頻分布就定義為分量c(t)的Hilbert譜:最后,匯總所有分量的Hilbert譜,就得到原始信號的Hilbert譜:綜上所述,HHT對信號的處理可用下面的流程圖來表示:式中HS(t,ω)表示時變Hilbert譜(Hilbertspectrum).上述基于Hilbert變換及經驗模態分解的HHT方法在處理非平穩信號時,得到的Hilbert譜由于不受不確定性原理的制約,它能同時在時域和頻域內獲得很高的分辨率.但是,困擾EMD方法的一個大問題就是數字信號的邊界問題.下面舉一個簡單的例子說明.通常我們所處理的信號都局限在有限區間.為此,考慮下式給出的模擬信號:圖1給出了限定在指定區間的原始模擬信號.可以看出,在區間[0.3,1.0]內,信號x(t)存在3個極大值點、3個極小值點;利用三次樣條插值分別連接極大值點和極小值點,所得到的樣條曲線自然延拓到信號的邊界處就得到了如圖所示的上下包絡線.圖2中的真實包絡線是按照如下方式得出的:首先,根據(13)式求出原始信號在區間[0.3,1.0]之外的部分;然后,利用三次樣條插值分別連接所有極大值點和極小值點;最后取所得樣條曲線在區間[0.3,1.0]內的部分.從圖2中可以看出,如果事先不知道信號在所給定區間之外的部分,只是利用給定區間內的極大值和極小值通過三次樣條插值來擬合信號的上下包絡線,所得的三次樣條包絡線與真實樣條包絡線相比會出現嚴重失真.在本例中,上下包絡線都出現了失真,尤其是上包絡線,在信號兩端出現了嚴重的失真.這種包絡線失真引起的誤差最初只會影響信號兩端,但隨著篩分過程的進行,邊界處的誤差會向內傳播,進而“污染”到內部的數據,使得最后的結果失去意義.而我們實際處理的信號都是有限長度的離散信號x(n),(n=1,…N),而且,我們事先無法知道給定數據之外的信號.因此,在進行EMD處理時,我們就需要根據已知數據,按照一定規則加以延拓,進而在信號兩端分別獲得附加的極大值和極小值;然后根據信號自身的所有極大值點和兩個附加的極大值點構造三次樣條上包絡線;根據信號自身的所有極小值點和兩個附加的極小值點構造三次樣條下包絡線;這兩條包絡線如果與真實樣條包絡線接近,則我們就能有效地控制數據邊界所帶來的誤差.本文提出“邊篩分、邊延拓”的邊界處理方法,利用自回歸模型(AR模型),通過線性預測對數據進行延拓.由于自回歸模型最后所得到的線性方程組的系數矩陣是對稱的Toeplitz矩陣,采用Levinson-Durbin遞推算法能夠快速地求得AR模型的系數,因此,在篩分的過程中,采用此算法對數據加以延拓,即“邊篩分、邊延拓”,能夠有效地對數據進行篩分,而且也能獲得較高的計算速度.2種邊端誤差的消除設x(n)之前的p個數據{x(n-p),x(n-p+1),…,x(n-1)}已知,我們希望利用這p個數據的線性組合來預測x(n)的值.記(n)是對真實值x(n)的預測,則有根據線性預測方法,通過使得預測值(n)與真實值x(n)之間總的預測誤差功率最小,系數{αk}可以求出.系數{αk}求出后,根據(14)式求出x(n)的預測值(n),然后根據新的數據序列{x(n-p+1),…,x(n-1),(n)}來預測n+1時刻x(n+1)的值(n+1),依次類推,可以求出x(n-1)以后任意時刻離散信號的預測值.當然,距離x(n-1)越遠的預測值與真實值之間的誤差也就越大,但是我們所要預測的是已知數據外最近的極大值和極小值,因此需要預測的數據點不會離已知數據的邊界太遠,從而預測的樣條包絡線能夠具有較好的精度.利用同樣的方法可以對同一數據序列進行反向預測.圖3給出了利用線性預測得出的(13)式信號的樣條上下包絡線及真實的包絡線,從中可以看出,兩者之間誤差非常小.圖4給出了相應的預測平均包絡線與真實平均包絡線,同樣可以看出,兩者非常接近.圖4所給出的平均包絡線參加最后的篩分運算,就可以得到物理意義明確的結果.在對一個給定信號(它可以為原始信號x(t),也可以為前一次篩分處理所得的殘余信號hij(t),見(2)式)進行篩分處理之前,首先利用AR模型對其進行線性預測,得出上下三次樣條包絡,然后再對其進行篩分,這種計算方法即為“邊篩分,邊延拓”的算法.按照此算法,對(13)式所給出的模擬信號進行經驗模態分解,所得結果與原始的兩個簡諧波信號如圖5所示.從中可以看出,對于高頻分量來說,按照本文所建議的算法,利用經驗模態分解所得的結果與真實分量相差不大,邊界效應僅在信號兩端引起微小的誤差,對中間大部分數據沒有影響;對于低頻分量來說,盡管數據邊端所引起的誤差傳播到了數據內部,但是,按照本文算法所得的低頻分量與真實分量之間非常接近,而且,此分量只有半個波存在,這種分解效果非常令人滿意.通過上面模擬信號的結果可以看出:數據邊端所引起的誤差對低頻分量的影響較大,對高頻分量影響較小;按照“邊篩分、邊延拓”的算法,利用AR模型,通過線性預測對數據進行延拓,可以有效地抑制數據邊端誤差對EMD分解結果中低頻分量的影響.下面,利用上述方法對一實際地震動記錄ElCentro波的分解為例加以說明.圖6給出了原始ElCentro波的加速度記錄.在此,首先使用Huang的EMD程序,對此地震波加以分解,得到8個IMF分量;然后使用本文提出的“邊篩分,邊延拓”的計算方法對同一地震波加以分解,同樣可以得到8個IMF分量;最后將兩種結果加以比較,以驗證本文方法的有效性.比較結果示于圖7中,由于數據邊端所引起的誤差對高頻分量影響不大,主要是影響低頻分量,所以圖中只給出了最后4個IMF分量的比較結果,對于前4個較高頻IMF分量來說,兩者幾乎完全一致,在此不再給出.圖7中,實線為本文所得結果;虛線為Huang的EMD軟件計算所得結果.從圖中可以看出,對于較高頻分量,兩者符合程度非常好;對于較低頻分量,盡管兩者存在一定誤差,但相差不是很大;雖然最后一個分量在邊界處相差較大,但是與高頻分量相比,它的尺度非常小,它的誤差可以忽略.3基于自回歸線性模型的原始數據建模本文提出一個EMD方法中邊界處理的“邊篩分、邊延拓”方法,利用AR模型,通過線性預測對信號兩端加以延拓,分別獲得附加的極大值點和極小值點,然后利用三次樣條插值將附加的極值點與信號本身的極值點連接起來,從而擬合出原始信號的上下包絡線