高 涵,袁希平,甘 淑,張 明

1. 昆明理工大學(xué)國土資源工程學(xué)院，云南 昆明 650093； 2. 云南省地震局信息中心，云南 昆明 650225； 3. 云南省測繪產(chǎn)品檢測站，云南 昆明 650034

由于具備全天候、高精度、三維定位的優(yōu)勢，GNSS連續(xù)觀測技術(shù)的應(yīng)用已逐漸由傳統(tǒng)空間大地測量學(xué)領(lǐng)域向地球物理學(xué)學(xué)科過渡，被地震學(xué)家廣泛應(yīng)用于地殼運動監(jiān)測當(dāng)中，逐漸成為地震地殼監(jiān)測不可或缺的補充手段[1-9]。伴隨現(xiàn)代信息技術(shù)的發(fā)展，人們不再將地震活動當(dāng)成瞬時狀態(tài)改變，而更加傾向于通過地殼活動的歷史規(guī)律及動態(tài)運動過程來理解地震孕育活動。

國內(nèi)外眾多學(xué)者通過地表監(jiān)測時間序列的方法對地震活動進行了深入研究。文獻[10]通過電磁觀測網(wǎng)在云南漾濞地震前夕從電磁時間序列中得到了地震前兆信息。文獻[11]通過Swarm衛(wèi)星建立了一種從衛(wèi)星地磁時間序列測量中檢測異常并進行地震預(yù)報的分析方法，通過研究得出高質(zhì)量的長期觀測數(shù)據(jù)是確認(rèn)與地震或其他自然災(zāi)害有關(guān)的異常的必要條件。文獻[12]通過綜合檢校法考慮同震階躍和RMS減小率，判斷可疑站點是否受地震影響,識別和估計了時間序列中的同震和震后形變。文獻[13]通過ARMA方法利用GPS時間序列對尼泊爾地殼形變和應(yīng)變進行了測量和建模分析，研究發(fā)現(xiàn)從尼泊爾北部到南部邊界的8～10 mm/a的主推力是該區(qū)域構(gòu)造應(yīng)變變化的主要原因。文獻[14]通過相對強度指數(shù)自動探測了單站GPS時間序列中的瞬間變形，驗證了該方法是一種能夠從地球物理位移時間序列中自動檢測異常的穩(wěn)定客觀的方法。然而，相對于地表位移時間序列，較少有學(xué)者專注于應(yīng)變時序歷史信息研究，通過云南省地震局的相關(guān)研究發(fā)現(xiàn)[15]，GNSS應(yīng)變反映的變化過程可能與強震的孕育-發(fā)生之間存在聯(lián)系，應(yīng)變及其異常變化的研究對未來強震危險區(qū)的判識具有重要意義。然而傳統(tǒng)的時頻分析方法中傅里葉變換是一種全局的變換，不在時域即在頻域，無法表述信號的時頻局部性質(zhì)，而時頻局部性質(zhì)恰好是非平穩(wěn)信號最基本和最關(guān)鍵的性質(zhì)。傳統(tǒng)的小波變換雖然具有自動改變窗口長短的功能，可以較好地把信號在時間上和頻域上局部化，從而呈現(xiàn)了信號的局部奇異性，但選擇小波基困難，小波基不能對信號進行自適應(yīng)。

綜上所述，為了深入挖掘應(yīng)變時序中所攜帶的信號特征，本文針對應(yīng)變時序選取專門適用于非線性非平穩(wěn)信號處理，又克服了傳統(tǒng)傅里葉變換和小波變換缺點的時頻分析方法，即整體經(jīng)驗?zāi)B(tài)分解的希爾伯特-黃變換分析方法(Hilbert-Huang transform-ensemble empirical mode decomposition，HHT-EEMD)來探索云南區(qū)域中強地震前GNSS應(yīng)變時序中的時間-頻率-能量聯(lián)合分布特征，嘗試挖掘應(yīng)變時頻信號中所攜帶的孕震信息，為未來云南區(qū)域強震的判定提供一定的參考。

1 HHT-EEMD數(shù)學(xué)模型

希爾伯特-黃變換(HHT)[16-18]主要包含2個步驟：經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition,EMD)與Hilbert譜分析(Hilbert spectral analysis,HSA)。首先，EMD可以自適應(yīng)地將原始信號分解為頻率由高至低的固有模態(tài)分量(intrinsic mode functions,IMF)和一個殘差趨勢項，分解得到的IMF分量具有一定的物理意義；其次，再對各個IMF分量進行Hilbert變換，獲得其對應(yīng)的瞬時振幅和瞬時頻率，進而把瞬時振幅表示在時頻表面上，就得到了Hilbert譜，這一過程稱為希爾伯特-黃變換，即HHT。

同時需要指出的是，EMD方法容易產(chǎn)生模式混疊現(xiàn)象，造成邊界效應(yīng)，而整體經(jīng)驗?zāi)B(tài)分解(ensemble empirical mode decomposition,EEMD)可利用隨機噪聲削弱這種現(xiàn)象[19-20]。因此，本文采用EEMD方法來將原始信號分解成IMF。

1.1 EEMD

EMD的原理是首先搜索原始信號s(t)所有的極大值點與極小值點，然后將所有極大值點和所有極小值點分別用一條曲線擬合得到s(t)上、下包絡(luò)線；若上、下包絡(luò)線的均值記作m(t)，用s(t)減去m(t)得到的新序列記作ci(t)，若ci(t)滿足一定的條件時，定義ci(t)為IMF的一個分量，否則視ci(t)為新的s(t)繼續(xù)重復(fù)以上操作，直至ci(t)滿足要求為止。由此可分解得到n個IMF分量和一個殘差趨勢項，即

(1)

式中，ci(t)為第i個IMF分量；rn(t)為殘余趨勢項。

而EEMD就是在進行每一次EMD分解之前，在原始信號s(t)中加入一定量的高斯白噪聲wk(t)，第k次待分解的信號表示為

sk(t)=s(t)+wk(t)

(2)

由于高斯白噪聲是一個均值為零的正態(tài)隨機序列，因此當(dāng)達到一定重復(fù)分解次數(shù)時，即可消除加入的白噪聲影響，使得分解后的IMF均來自s(t)本身。

利用EEMD可以實現(xiàn)原始時間序列中噪聲的分離，EEMD是將組成原始信號的各尺度分量不斷從高頻到低頻進行提取，則分解得到的IMF順序是按頻率由高到低進行排列的，即首先得到最高頻的分量，然后是次高頻的，最終得到一個頻率接近為0的低頻殘余分量。在地震監(jiān)測時間序列分析中，高頻信息頻率較高代表噪聲、誤差等短周期抖動；低頻IMF分量相對高頻分量較為平緩，主要代表區(qū)域地殼中長趨勢運動特征。

1.2 Hilbert譜分析

對給定的連續(xù)曲線X(t)，其Hilbert變換定義為

(3)

式中，PV表示奇異積分的主值；τ為平行移動窗口參數(shù)。通過Hilbert變換，解析信號定義為

Z(t)=X(t)+iY(t)=a(t)ei θ(t)

(4)

(5)

(6)

式中，a(t)是瞬時振幅；θ是相位函數(shù)；瞬時頻率ω(t)可表示為

(7)

在得到IMF分量后，很容易將Hilbert變換應(yīng)用于每個IMF分量，然后根據(jù)式(3)—式(7)計算瞬時頻率。對每個IMF分量X(t)執(zhí)行Hilbert變換后，H(ω,t)可以表示為

(8)

式中，Re[?]表示“?”的實部；aj(t)和ωj(t)分別是第j個IMF的瞬時振幅和頻率。振幅的頻率-時間分布稱為“Hilbert振幅譜”H(ω,t)，或者簡稱為“Hilbert譜”。定義Hilbert譜后，進而可以將邊際譜H(ω)定義為

(9)

2 數(shù)據(jù)處理及分析思路

本文選取了包括陸態(tài)網(wǎng)絡(luò)及云南省地震局自建站在內(nèi)的云南區(qū)域59個GNSS連續(xù)跟蹤站2013—2019年所觀測到的30 s采樣間隔的觀測數(shù)據(jù)，首先利用GAMIT軟件進行基線解算，獲得測站的區(qū)域單日松弛解；然后通過GLOBK軟件將區(qū)域單日松弛解和SOPAC提供的IGS全球單天無約束解聯(lián)合平差，得到ITRF2014框架下的測站4日解坐標(biāo)時間序列，進而獲取得到區(qū)域GNSS位移場；最后，將研究區(qū)域(97°E—107°E，21°N—29°N)進行1°×1°格網(wǎng)化，用網(wǎng)格化后的格網(wǎng)最終代表該格網(wǎng)經(jīng)緯度范圍內(nèi)的應(yīng)變情況，并按照從S至N，W至E的順序進行編號，如圖1所示。其中，震源球表示研究區(qū)域內(nèi)2013—2019年Ms≥5級地震的震源機制解。

圖1 網(wǎng)格劃分及站點分布[21]Fig.1 Grid division and site distribution[21]

劃分格網(wǎng)后，參考文獻[22—23]中的平面應(yīng)變計算公式以計算不同格網(wǎng)的應(yīng)變分量。首先需要通過Kriging插值方法對區(qū)域內(nèi)分布不均勻的站點位移進行插值，得到不同格網(wǎng)的位移數(shù)據(jù)，進而通過應(yīng)變計算得到80個格網(wǎng)的應(yīng)變參數(shù)時間序列，部分應(yīng)變時間序列如圖2所示。其中，黑色豎線標(biāo)注為格網(wǎng)時間點對應(yīng)的地震。由于地震活動是大氣-陸地相互作用、海洋-陸地相互作用、人類-地面相互作用、地球內(nèi)部之間相互作用、火山活動等因素共同作用的結(jié)果，這必然導(dǎo)致地震觀測中存在各種雜亂無章的噪聲信號。如圖2所示，由于應(yīng)變時間序列中包含了各類噪聲干擾，有可能掩蓋地震前的異常信號。因此，在分別獲取得到80個格網(wǎng)的應(yīng)變時間序列后，還需利用適用于非線性非平穩(wěn)信號處理的HHT-EEMD分析方法對結(jié)果進行分析去噪，進一步嘗試挖掘孕震信息。

圖2 格網(wǎng)應(yīng)變時間序列Fig.2 Strain time series of grid

HHT-EEMD分析思路如下。

(1) 對每個格網(wǎng)應(yīng)變時間序列進行EEMD分解，將應(yīng)變時序分解為不同頻率的IMF分量和殘差趨勢項。

(2) 殘差趨勢項分析：殘差趨勢項主要代表格網(wǎng)長周期運動趨勢，對部分典型格網(wǎng)及周邊格網(wǎng)進行聯(lián)合分析，可判斷格網(wǎng)與周邊格網(wǎng)應(yīng)力狀態(tài)、張壓特性等，進而初步從宏觀了解格網(wǎng)區(qū)域內(nèi)應(yīng)變積累態(tài)勢。

(3) IMF分量異常曲線識別：通過顯著性檢驗從不同頻率的IMF分量中選取出包含有效信息的IMF分量，再劃定2倍標(biāo)準(zhǔn)差曲線作為選取的IMF分量的異常指標(biāo)閾值，超過即視為異常。

(4) 進一步對選取的IMF分量進行Hilbert變換，獲取該分量的時間-頻率-振幅聯(lián)合分布特征(Hilbert譜)，以及瞬時頻率、瞬時振幅等局部特征(Hilbert瞬時能量譜)。

(5) 聯(lián)合云南區(qū)域的歷史震例，結(jié)合上述HHT—EEMD分析結(jié)果進行震例總結(jié)，從應(yīng)力、構(gòu)造活動等方面對格網(wǎng)內(nèi)歷史地震進行分析，嘗試挖掘應(yīng)變時頻信號中所攜帶的潛在信息。分析流程如圖3所示。

圖3 處理及分析思路Fig.3 Processing and analysis ideas

3 結(jié)果分析與討論

3.1 殘差趨勢項分析

如前所述，首先對獲取得到的80個格網(wǎng)應(yīng)變時間序列進行EEMD分解，得到與之對應(yīng)的IMF分量和殘差趨勢項，殘差趨勢項如圖4所示。圖中每個格網(wǎng)的x方向代表時間，y方向代表格網(wǎng)應(yīng)變時序?qū)?yīng)進行EEMD分解后得到的長周期殘差趨勢項。由圖4可以發(fā)現(xiàn)，殘差趨勢項消除了高頻IMF分量噪聲的影響，體現(xiàn)了原應(yīng)變時間序列的總體運動趨勢。對格網(wǎng)及周邊格網(wǎng)進行殘差趨勢項聯(lián)合分析，可以判斷格網(wǎng)相對于周邊格網(wǎng)應(yīng)力狀態(tài)、張壓特性等特征隨時間的相對變化趨勢，進而對比發(fā)現(xiàn)并找出格網(wǎng)間的差異性運動，從宏觀角度了解區(qū)域應(yīng)變相對積累態(tài)勢。不同格網(wǎng)應(yīng)變-時間分析方法與傳統(tǒng)的靜態(tài)應(yīng)變相比，具有豐富時間維度信息，突出應(yīng)變變化細(xì)節(jié)的優(yōu)勢。

圖4 2013—2019年網(wǎng)格應(yīng)變時序經(jīng)EEMD分解后的殘差趨勢項Fig.4 Residual trend term of each grid stain time series after decomposition by EEMD during 2013—2019

由圖4最大剪應(yīng)變殘差趨勢項可以發(fā)現(xiàn)，除靠近東部區(qū)域華南板塊的格網(wǎng)以外，云南區(qū)域剪應(yīng)變運動整體活躍，小江斷裂帶(F5)、蓮峰斷裂(F6)相關(guān)格網(wǎng)為整個云南的剪應(yīng)變相對高值區(qū)，與云南區(qū)域其他格網(wǎng)相比，該斷裂帶區(qū)域剪切運動更為突出，這與歷史上這一斷裂帶周邊頻發(fā)走滑性質(zhì)強震的現(xiàn)象相符。另外，瀾滄江斷裂(F1)北段呈顯著面應(yīng)變高值，而南段剪應(yīng)變呈現(xiàn)出相對高值，說明這個區(qū)域的剪切運動較為活躍，這與瀾滄江斷裂帶在新生代強烈活動,且以走滑逆沖運動為主相符。這可能是在太平洋板塊向西推擠作用下的揚子華南板塊與印度板塊發(fā)生強烈碰撞擠壓、向北運移，在北部強烈擠壓，在南部強烈運移共同作用的產(chǎn)物。同時，結(jié)合圖1和表1整理收集的云南區(qū)域2013—2019年Ms≥5級的歷史震例震源機制解信息來看，該區(qū)域的剪應(yīng)變活動高值異常與統(tǒng)計年間區(qū)域內(nèi)地震活動活躍較為相符。

由圖4面應(yīng)變殘差趨勢項可以看到，以紅河斷裂(F2)為界，北部面膨脹活動較為活躍，川滇菱形塊體西北部拉張明顯，由西向東逐漸由張轉(zhuǎn)壓，整個云南區(qū)域主要呈現(xiàn)張壓交替的時空演化特征，這種現(xiàn)象表明紅河斷裂在云南區(qū)域的面膨脹構(gòu)造運動中起重要作用。鶴慶-洱源斷裂(F3)、小金河斷裂(F4)附近，面應(yīng)變伴隨時間的推移而展現(xiàn)出局部相對高值的運動趨勢，由表1可知，云南區(qū)域內(nèi)地震活動主要以走滑斷層破裂為主，而該區(qū)域統(tǒng)計期間卻于2013年發(fā)生了正斷型斷層破裂的Ms5.9德欽地震，這種發(fā)震斷層活動性質(zhì)與區(qū)域格網(wǎng)運動相符的現(xiàn)象表明地震孕育活動與區(qū)域內(nèi)剪應(yīng)力和張壓應(yīng)力伴隨時間發(fā)展的相對運動趨勢高度相關(guān)。

表1 云南區(qū)域Ms≥5.0級地震統(tǒng)計表(2013—2019年)Tab.1 Earthquake statistical table above Ms5.0 in Yunnan region(2013—2019)

3.2 震例總結(jié)

前文雖然基于EEMD分解的應(yīng)變時序殘差趨勢項從宏觀角度分析了云南區(qū)域長趨勢應(yīng)變積累態(tài)勢，但無法獲取應(yīng)變信號中所包含的瞬時活動特性，因此下文將進一步結(jié)合HHT分析方法，通過顯著性檢驗選取出較為合適的IMF分量進行Hilbert變換，獲得其時間-頻率-振幅(能量)的聯(lián)合分布特征(Hilbert譜)，以及瞬時頻率和瞬時振幅等局部特征(Hilbert瞬時能量譜)，并結(jié)合區(qū)域構(gòu)造信息對部分格網(wǎng)進行震例總結(jié)，嘗試挖掘應(yīng)變時序時頻信號中所攜帶的孕震信息。

3.2.1 震例1: 23號格網(wǎng)對應(yīng)地震

23號格網(wǎng)內(nèi)，分別發(fā)生了2014年10月7日景谷Ms6.6級，以及12月7日景谷Ms5.8、Ms5.9級地震，震中均位于云南省普洱市景谷傣族彝族自治縣。另外，2018年9月8日的墨江Ms5.9級地震震中在云南普洱市墨江縣，位于23號格網(wǎng)右側(cè)相鄰的24號格網(wǎng)。由圖4的面應(yīng)變殘差趨勢圖可以看到，該區(qū)域的面應(yīng)變活動稍有壓性但整體較為平緩，而剪應(yīng)變活動較為劇烈，甚至較少有衰減現(xiàn)象，這說明該區(qū)域內(nèi)的走滑斷裂中長期處于活躍狀態(tài)，其構(gòu)造應(yīng)力主要受剪應(yīng)力支配。區(qū)域內(nèi)剪應(yīng)變變化活躍、面應(yīng)變變化平緩的應(yīng)變特性與景谷地震發(fā)震于右旋走滑的隱伏斷層和墨江地震發(fā)震于右旋走滑的阿墨江斷裂帶西支斷裂的現(xiàn)象相符，這從斷裂活動特性與斷層調(diào)查性質(zhì)一致的角度驗證了本文中應(yīng)變殘差趨勢分析的適用性。

圖5(a)和圖6(a)分別是23號格網(wǎng)最大剪應(yīng)變時序和面應(yīng)變時序進行EEMD分解后的8個IMF分量結(jié)果?？梢钥吹剑琁MF的各個分量按照從高頻到低頻的順序依次排列，高頻IMF頻率較高代表噪聲、誤差等短周期抖動，殘差趨勢項代表區(qū)域地殼中長趨勢背景運動特征，可見EEMD具有分頻剖面的類似特征，可以利用信號與噪聲在不同固有模態(tài)函數(shù)剖面上的差異來去噪，分離高頻噪聲和低頻趨勢，進而突出異常信號特征。

獲得IMF分量之后，再對其進行顯著性檢驗選出包含有效信息較多的IMF分量，如圖7所示。由文獻[24]可知，IMF分量于置信曲線上且越遠離置信區(qū)間曲線，說明IMF中所攜帶的有用信息越多。從最大剪應(yīng)變顯著性檢驗圖中可以看出，IMF5和IMF6分量相對其他分量距離回歸曲線較遠；而從面應(yīng)變顯著性檢驗圖中可以看出，IMF5分量相對其他分量距離回歸曲線較遠，因此，筆者認(rèn)為這些中低頻的IMF分量中可能存在有用信息，下面將嘗試從上述IMF分量中提取異常信息。震例2的IMF分量選取方法同上，后續(xù)不再贅述。

圖7 23號格網(wǎng)最大剪應(yīng)變和面應(yīng)變的IMF分量顯著性假設(shè)檢驗Fig.7 The significance of dilation strain and maximum shear strain IMFs for No. 23 grid cell

選定IMF分量后，首先對其進行異常曲線識別。由于地塊并非“完美”剛性塊體，受外部作用力后，塊體的運動狀態(tài)將以應(yīng)力強弱的形式表現(xiàn)出來，如果通過塊體的歷史應(yīng)變狀態(tài)建立線性回歸模型，則當(dāng)塊體本身的變形或狀態(tài)異常時，必定會與原有的回歸模型產(chǎn)生較大的偏差，統(tǒng)計后發(fā)現(xiàn)去趨勢后的應(yīng)變值服從正態(tài)分布，95%以上的數(shù)據(jù)均在2倍標(biāo)準(zhǔn)差之內(nèi)[22-23,25]，因此，本文以2倍中誤差為限設(shè)置異常指標(biāo)線，超出2倍中誤差則認(rèn)為塊體運動異常。而在計算2倍中誤差過程中，常規(guī)的統(tǒng)計方法不考慮數(shù)據(jù)的預(yù)測時間節(jié)點，是一種依賴時間的未來函數(shù)，本文采用對時間逐漸累積，數(shù)據(jù)逐漸參與計算的方法，來進行異常指標(biāo)線獲取，即對于歷史信息T日只用T-1日前的數(shù)據(jù)參與計算，劃定2倍標(biāo)準(zhǔn)差作為異常指標(biāo)閾值，最終獲得不依賴于T+1日后的數(shù)據(jù)得到的異常指標(biāo)曲線。選定的IMF分量及其異常指標(biāo)曲線如圖5(b)、(c)和圖6(b)所示。

圖5 23號格網(wǎng)最大剪應(yīng)變時間序列HHT-EEMD結(jié)果Fig.5 The maximum shear strain HHT-EEMD results for the No. 23 grid cell

然后，對選定的IMF分量進行Hilbert變換，以獲得瞬時頻率、瞬時振幅等局部時頻信息，圖5(d)和圖6(c)是最大剪應(yīng)變時序和面應(yīng)變時序分別對應(yīng)的總Hilbert譜。圖5(e)—圖5(h)和圖6(d)、圖6(e)是選定的IMF分量經(jīng)Hilbert變換獲取的Hilbert譜、瞬時振幅與瞬時頻率。從整體來看，剔除了高頻噪聲和低頻趨勢后的中頻應(yīng)變IMF曲線，無論是最大剪應(yīng)變還是面應(yīng)變均在景谷兩次地震和墨江地震前夕出現(xiàn)異常波動，均超過了圖中的2倍標(biāo)準(zhǔn)差曲線，且均呈現(xiàn)出瞬時振幅突然增大的現(xiàn)象，這可能與當(dāng)前IMF分量頻率對應(yīng)的地殼復(fù)雜運動和地殼震前地震能量的早期釋放有關(guān)。同時，也表明HHT-EEMD在地震信號分析中具有較好的描述信號局部時-頻特性的能力，能夠一定程度描述應(yīng)變時序的時變特性，具有較高的時間和頻率分辨率。假設(shè)選定IMF中地震前夕瞬時能量的突然增加是有效地震信號的前提下，由圖5(f)、(h)和圖6(e)可知，包含信號的IMF會長時間處于正常波動狀態(tài)，并且只會在地震發(fā)生前不久產(chǎn)生異常。上述現(xiàn)象從時頻分析的角度解釋了為什么只能找到地震所屬板塊的運動趨勢，而較早進行地震預(yù)測非常困難。在地震發(fā)生前的較短時間內(nèi)雖然可能會有異常擾動信號但如果要在地震前有效地挖掘到孕震信息，則相關(guān)信號時間序列的計算周期必須小于地震信號出現(xiàn)的時間。另外，格網(wǎng)最大剪應(yīng)變的異常比面應(yīng)變突出，這也與景谷地震以及墨江地震發(fā)震于右旋走滑斷裂的性質(zhì)相符。

圖6 23號格網(wǎng)面應(yīng)變時間序列HHT-EEMD結(jié)果Fig.6 The dilation HHT-EEMD results for No. 23 grid cell

3.2.2 震例2： 42號格網(wǎng)對應(yīng)地震

42號格網(wǎng)內(nèi)，發(fā)生了2013年3月3日、4月17日的洱源Ms5.5、Ms5.0級地震，2015年10月30日的昌寧Ms5.1級地震，以及2017年3月27日漾濞Ms5.1級地震。洱源地震發(fā)震于右旋走滑兼正斷性質(zhì)的維西-喬后斷裂；通過震源機制解評估昌寧地震可能發(fā)震于42號格網(wǎng)內(nèi)的正斷性質(zhì)隱伏斷裂；漾濞地震發(fā)震于右旋走滑性質(zhì)的維西-喬后斷裂中南段。

圖8和圖9分別是42號格網(wǎng)最大剪應(yīng)變時序和面應(yīng)變時序HHT-EEMD的結(jié)果?？梢钥闯觯m然最大剪應(yīng)變的IMF5、面應(yīng)變的IMF4也出現(xiàn)了在地震前瞬時振幅突然增大的現(xiàn)象，但與震例1中不同的是：在42號格網(wǎng)內(nèi)敏感IMF分量并未出現(xiàn)超過2倍標(biāo)準(zhǔn)差現(xiàn)象，但敏感IMF在經(jīng)過Hilbert變換后，卻在地震前呈現(xiàn)出一定的異常信號，這表明Hilbert變換能夠有效地將信號分解在時頻維度，在IMF分量異常曲線識別無效的情況下，仍然能夠更深層次地通過局部的瞬時頻率、瞬時振幅等方式凸顯異常，對于應(yīng)變時序的異常識別具有一定的適用性。

圖8 42號格網(wǎng)最大剪應(yīng)變時間序列HHT-EEMD結(jié)果Fig.8 The maximum shear strain HHT-EEMD results for the No. 42 grid cell

從震例結(jié)果來看，經(jīng)過EEMD分解后的每個IMF分量代表了原始信號相互獨立的各頻率分量，其能將原始應(yīng)變時間序列中的不同頻率信號清晰地剝離，揭示信號里蘊含的高頻、中頻和低頻等多尺度振蕩特征。EEMD分解得到的IMF分量中：高頻分量可能主要由各種噪聲組成，因此離散度高。中頻分量可能包含一定的前兆信息，但其中常伴有季節(jié)性、年變等周期擾動。當(dāng)判斷中頻IMF分量顯著性時，可通過顯著性檢驗方法選取包含較多信息的IMF分量后，再進一步對其進行異常曲線識別和判斷，以及Hilbert變換，可獲得信號的瞬時頻率、瞬時振幅等局部特征，進而更深層次地凸顯信號中所含的異常信息，為中強地震的判定提供一定的參考輔助。低頻IMF分量則多代表信號的長周期趨勢性變化。EEMD算法雖然是一種完備性好，對平穩(wěn)信號、非平穩(wěn)信號、非線性信號都適用的自適應(yīng)的信號處理方法，其與EMD相比，能在一定程度上克服模態(tài)混疊現(xiàn)象，但該方法添加的白噪聲殘留可能給信號帶來噪聲干擾，如圖8(e)中震后峰值、圖9(e)中的第二高峰值。雖然白噪聲有助于削弱模態(tài)混疊現(xiàn)象，但在信號中引起了噪聲殘留，因此當(dāng)異常信號出現(xiàn)時，需要多種異常判據(jù)進行相互佐證，綜合判斷異常指標(biāo)曲線、瞬時頻率和瞬時振幅的共振程度，是否有關(guān)的異常判據(jù)均處于極值、拐點或其他態(tài)勢，并通過不斷積累數(shù)據(jù)建立信號回歸經(jīng)驗方程以不斷提高孕震分析的可靠性。此外，分解得到的IMF分量偶爾會存在上、下包絡(luò)在數(shù)據(jù)序列兩端發(fā)散的現(xiàn)象，且這種發(fā)散會隨著運算的進行而逐漸向內(nèi)，從而使得整個數(shù)據(jù)序列受到影響，造成信號的端點效應(yīng)。端點效應(yīng)是影響EEMD分解精度的主要因素，如圖5(b)MJ5.9后面的峰值信息，在靠近結(jié)束邊緣區(qū)域出現(xiàn)了非常大的抖動。當(dāng)遇到端點異常信號時，不能完全確定端點處的異常信號是有效信號，需要進一步通過鏡像法、極值延拓法、神經(jīng)網(wǎng)絡(luò)預(yù)測、多項式外延方法、平行延拓法、邊界局部特征尺度延拓法等算法來削弱抑制端點效應(yīng)，綜合分析區(qū)域內(nèi)其他地震地質(zhì)相關(guān)資料，利用除此IMF分量外的其他異常因素來推測其發(fā)震的可能性。

圖9 42號格網(wǎng)面應(yīng)變時間序列HHT-EEMD結(jié)果Fig.9 The dilation HHT-EEMD results for No.42 grid cell

4 結(jié) 語

本文基于云南區(qū)域2013—2019年GNSS格網(wǎng)應(yīng)變時間序列，利用HHT-EEMD分析方法，結(jié)合震例探索云南區(qū)域中強地震前GNSS應(yīng)變時序的時-頻-能量分布特征，研究其作為孕震信息的有效性。結(jié)果顯示：

(1) EEMD可依據(jù)數(shù)據(jù)的時間特征尺度進行分解，將不同頻率的信號依次剝離，能有效分離高頻噪聲和低頻趨勢，較好地剖析地震信號在不同頻率尺度上的變化特征。

(2) 通過繪制、對比不同格網(wǎng)的殘差趨勢項，能夠發(fā)現(xiàn)不同格網(wǎng)間隨時間變化的相對應(yīng)變差異運動趨勢，與傳統(tǒng)的靜態(tài)應(yīng)變圖像相比，殘差趨勢項圖豐富了時間維度信息。

(3) 通過顯著性檢驗選取包含信息較多的IMF分量，并對其進行Hilbert變換，可獲得瞬時頻率、瞬時振幅等局部時頻信息。在地震前，選取的IMF分量瞬時振幅突然增大的現(xiàn)象可能與當(dāng)前IMF分量頻率對應(yīng)的地殼復(fù)雜運動和地殼震前地震能量的早期釋放有關(guān)。

(4) Hilbert變換能夠通過信號的瞬時局部特性描述信號隨時間的細(xì)微變化，在異常曲線識別無效的情況下仍能更深層次凸顯異常，對于時序的異常識別具有一定的適用性。

(5) 通過EEMD、殘差趨勢項分析、IMF分量異常識別和Hilbert變換的方法綜合動態(tài)分析應(yīng)變時間序列，能夠在部分地震前夕發(fā)現(xiàn)一些潛在異常信息，為未來云南區(qū)域強震危險地點的判定提供一定的參考，對于未來更加科學(xué)地挖掘應(yīng)變時序中蘊含的孕震信息有一定探索意義。

致謝：感謝云南省地震局信息中心提供的數(shù)據(jù)支撐；感謝云南地震臺洪敏老師給予的支持與幫助。