孔祥瑞,翟麗娜,李夢瑩,張志宏,黃明威
(遼寧省地震局,遼寧 沈陽110034)
形變觀測是通過精密儀器,對地震前發生的與地震孕育和發生過程有關的構造運動、震前變形、同震形變、震后恢復等相關現象的觀測。因此,形變觀測產生的數據,極有可能蘊含著地震前的前兆信息。一般情況下,地震前兆觀測數據中含有不同的周期成分有如下幾種:(1)趨勢變化,即不同變化的長期變化,時間尺度為幾年,甚至幾十年不等;(2)年周期變化,表現為各種因素的綜合影響引起的,年變一般每年存在極大值、極小值;(3)固體潮變化,前兆觀測資料含有的固體潮汐資料,整點值數據可以刻畫固體潮變化;(4)各種隨機因素引起的觀測數值的隨機變化,一般為短期[1-2]。
由于形變觀測儀器在獲取震前形變信號變化的同時,也會記錄到人類生產、活動的干擾,以及觀測自身系統的噪聲影響所產生的干擾信號,這些都會使得觀測數據具有非平穩性、非線性特征。隨著數字化形變觀測儀器大量投入觀測,數據觀測頻率由小時值數據提高到分鐘值數據,采樣率大幅度提高。但與此同時,各種高頻干擾信息如風擾、同震變化,使得觀測儀器會記錄到更多的干擾因素,為我們分析提取形變異常信息造成了一定的困難,因此排除這類干擾因素越發重要。傳統的處理方法如FFT 或者小波變換的前提條件,都要求所處理的信號必須具有線性和平穩性的特征,因此客觀上講,傳統的處理方法處理前兆資料具有一定的局限性[3]。
本文采用近年出現的經驗模態分解模型(EmpiricalModeDecomposition,EMD)[4-6],將形變觀測信號經過EMD分解為一組不同尺度內在的、客觀的固有模態函數(IntrinsicMode Function,IMF)。基于EMD的優良性能,本文首先利用EMD 模型對形變觀測數據進行一定分解,在獲得IMF 之后再把某些IMF 進行適當的組合,便可以構成特定意義下的低通、帶通與高通自適應濾波器,然后進行數值仿真試驗,最后處理形變觀測資料,探索消除形變觀測資料中的干擾因素。
經驗模態分解方法(EMD)是將原始信號分解成若干固有模態函數(IMF),在這個過程中原始信號的各種頻率的成分以不同固有模態函數的形式從原始時問信號中分解出來。其中固有模態函數滿足:(1)整個數據中,零點數與極點數相等或至多相差1;(2)信號上任意一點,由局部極大值點確定的包絡線和由局部極小值點確定的包絡線的均值均為0。
IMF 求取的步驟如下:對一原始信號X(t),先求出這個序列的極大值點構成的上包絡線和極小值點構成的下包絡線的均值即為m,將原始信號X(t)減去包絡均值m,得到一個新的數據序列h,即:

重復以上操作,直至h滿足IMF 條件時,可將h視為第一個IMF,即:

將r 視為新的X(t),重復以上步驟,即可依次獲得IMF1、IMF2…最終獲得分解結果:

其中,X(t)為原始信號,IMF 代表固有模態函數,r 為殘余函數,即剩余的趨勢項部分。
基于以上分解獲得的IMF 序列,我們可以完全進行重構,將IMF 序列按頻率由高到低進行排列,將干擾所影響的優勢頻率進行剔除,從而進行高頻干擾去除。
通過對已知模擬信號進行EMD分解,是了解其效果最為簡單的方法,為此設計包含隨機噪聲的已知信號,其中圖1(a)基礎函數為x=cos(2×π×0.5×t-π/2);添加隨機噪聲參數為0.3倍符合標準正態分布的一維矩陣的函數圖1(b),最終可以獲得隨機噪聲的時間序列圖1(c)。

圖1 隨機噪聲信號時序曲線Fig.1Sequentialcurvesof random noisesignal

圖2 隨機噪聲信號時序分解及結果曲線Fig.2 Sequencedecompositionandresult curvesof random noise signal
將含有隨機噪聲信號時序經由EMD 模型分解,可以獲得由6重IMF 組成的重組濾波(圖2),而我們可以將分解的6重信號進行疊加,可以獲得原始噪聲序列曲線(圖1(c))。各個分量所對應的時頻譜圖像如圖3所示,主頻率從高到低依次排列,其中高頻信息主要集中在前三個分量(a)、(b)、(c)上,所以干擾噪聲也就集中在此部分,因此我們可以排除這些信號進行重組,這樣我們既可以保持原曲線的形態,又可以剔除高頻干擾信息。

圖3 各個IMF 對應時頻譜圖像Fig.3The variousIMF imagesof theseasonal spectrum
為了驗證獲得數據的有效性,我們可以將圖2(a)、(b)、(c)進行疊加,與原噪聲曲線所添加的噪聲圖1(b)進行相關性分析,分析所提取出的高頻信息與添加的噪聲的相關性,經過對比,其相關性系數達到0.930,基本滿足提取高頻信息的需求,使用此方法我們可以提取震前形變高頻異常,或者在保存數據趨勢變化的基礎下,去除例如風擾等外界環境和儀器本身的高頻信息,為分析預報提供幫助。
同樣道理我們將圖2(d)、2(e)、2(f)進行疊加,與未添加噪聲的曲線圖1(a)進行相關性分析,其相關性系數為0.999。此方法可以在保證趨勢變化的基礎上,對數據進行濾波,濾波的尺度我們可以自行掌控。由此可以說明信號經過合理的分解與重構,基本可以獲得高頻和低頻信息,為進一步的數據分析提供幫助。
隨著數字化觀測儀器的普及,分鐘值形變觀測數據中高頻干擾是經常出現的現象,例如風擾、降雨、同震等影響,這些干擾并不會影響數據的曲線形態,而是對數據產生高頻干擾,這些干擾或是瞬時的,如同震;或是持續的,如風擾。我們對這類信號主要是去除高頻信息,按上述的例子,只需將原始信號減去先剔除的高頻信號即可,以丹東水管EW 分量為例,丹東水管傾斜自觀測以來,數據一直良好,年變形態清晰,但該測項分量受風擾干擾尤為明顯,具體體現為高頻波動,尤其在每年的3—6月,風擾更加明顯,本文選取2016年3—4 月份數據,在3月21日存在一次同震響應(圖4)。

圖4 丹東水管EW 分量觀測曲線Fig.4The EW component observationcurveof Dandong water pipe
按照上述方法進行分解,可以獲得不同頻率的多個IMF 函數,本例子獲得13階IMF 函數,通過對分解出的各個階IMF 函數分析,可發現大風擾動和同震的高頻干擾信息主要集中在前六層IMF 函數,故將其前六層重新組合可獲得高頻干擾信息,因此我將前六層信息進行組合,獲取干擾信息(圖5),將其他剩余的IMF函數進行重新組合,便可獲取去除高頻干擾的觀測曲線,為進一步的分析提供幫助(圖6)。

圖5 丹東水管EW 分量高頻干擾曲線Fig.5High frequency interferencecurve of EW component of Dandong water pipe

通過對丹東水管EW 分量的風擾高頻干擾和同震響應的提取和去除,我們可以基本去除其影響信息,并且去除干擾后的函數也并沒有因為濾波的因素而改變和影響本序列原有的信息,因此該方法盡可能的在保留固體潮等數據特征的基礎上去除干擾,為數據進一步分析提供支持。
本文首先通過對EMD方法及原理進行簡單介紹,然后對理想下的已知模擬信號進行人為加入噪聲信息,用該方法進行實驗,最后把EMD方法技術應用于形變信號處理領域之中。通過以上綜合分析可得到如下結論:
(1)EMD分解中首先被提取出來的IMF 頻率最高,第一個被EMD方法提取出來的IMF函數表示原信號的最高頻率成分,再此之后每一層的信號頻率逐漸降低,直到分解到最低頻停止分解,每一個IMF 序列均穩定在一個頻率之間,也正因此,每一個IMF 序列都具有一定的物理意義。并且該方法并沒有對引號進行過多干預,保證信號的原本特性的真實性。
(2)通過對存在的典型大風干擾的丹東臺水管觀測時間序列進行EMD分解和IMF組合濾波,進一步表明基于EMD分解在形變信號處理時的方法與思路,并且通過對數據處理之后的形態來看,該方法最大幅度的保存了信號的原始信息,如固體潮變化。因此認為,該方法對形變信號的提取高頻信息具有有效性和可靠性。