徐玉聰 朱建 董建輝 周魯

摘要：為識別仙女山斷裂周緣發(fā)生的構(gòu)造地震、礦震、塌陷地震特征，基于集合經(jīng)驗(yàn)?zāi)B(tài)分解時頻分析方法（EEMD）對地震波形信號進(jìn)行了分解，通過信號的IMF分量的時頻譜和邊際譜分析了其差異。結(jié)果表明：① 構(gòu)造地震的時頻譜頻率在0～40 Hz均有分布，頻率成分復(fù)雜，邊際譜集中在0～20 Hz;② 礦震時頻譜在0～30 Hz均有分布，頻率成分相對簡單，但在地震波幅值變化最大的區(qū)域存在20～30 Hz的頻率成分空白區(qū)，邊際譜存在一個0～5 Hz的主頻帶和2 Hz左右的主頻;③ 塌陷地震時頻譜頻率成分分布在0～30 Hz，在Pg波到達(dá)前存在0～5 Hz左右低頻區(qū)域，邊際譜的主頻帶和主頻集中在0～3 Hz和1 Hz。該研究成果可以作為識別構(gòu)造、礦震、塌陷地震事件的判別依據(jù)。

關(guān)鍵詞：EEMD; 仙女山斷裂; 時頻分析; 邊際譜; 三峽水庫

中圖法分類號：P315 文獻(xiàn)標(biāo)志碼：A DOI：10.15974/j.cnki.slsdkb.2022.02.002

文章編號：1006 - 0081（2022）02 - 0009 - 08

0 引 言

仙女山斷裂是三峽庫區(qū)的主要斷裂帶之一，位置如圖1所示，位于黃陵地塊西南側(cè)，北起秭歸荒口以北的風(fēng)吹埡，南止于五峰漁洋關(guān)，全長80余千米，距三峽壩址最近處約19 km。自從三峽水庫蓄水后，仙女山斷裂周緣發(fā)生了大量的微小地震活動，大致可以分為構(gòu)造地震、礦震、塌陷地震。有效識別這些微小地震類型對于仙女山斷裂附近的地震活動特征和趨勢研究具有重要意義。

國內(nèi)外研究人員針對不同地震類型地震波的時頻譜特征識別開展了大量的研究工作，其中主要集中在短時傅里葉變換（STFT）[1-2]、小波變換（WT）[3-5]、戈勃（Gabor）展開、維格納-威爾分布（Wigner-Ville）[6]等。雖然上述時頻分析方法在應(yīng)用中取得了一定成果，但是這些方法本身也存在一些局限性，主要由于以下幾種原因造成：這些方法都基于線性理論、需選取不同時窗、需引入交叉項(xiàng)或者依賴于前驗(yàn)的基函數(shù)等，而地震波作為一種非線性、非平穩(wěn)的隨機(jī)信號，需要依賴于信號自身非線性、非平穩(wěn)的特點(diǎn)進(jìn)行分析計(jì)算，因此，合適的時頻譜分析方法對于分析研究結(jié)果也會產(chǎn)生較大影響。

希爾伯特-黃變換（Hilbert-Huang，HHT）作為一種非線性、非平穩(wěn)信號的處理方法，通過經(jīng)驗(yàn)?zāi)B(tài)分解（Empirical Mode Decomposition，EMD） [7]，可基于信號自身特點(diǎn)，將復(fù)雜地震信號自適應(yīng)地分解為一系列固有模態(tài)函數(shù)（Intrinsic Mode Function，IMF），通過對每個IMF分量進(jìn)行Hilbert變換，可以準(zhǔn)確獲取反映信號時頻特征的瞬時頻率和振幅、邊際譜等結(jié)果，目前在地震分析、機(jī)械診斷等信號處理中得到了廣泛應(yīng)用[8-11]。

本文利用基于集合經(jīng)驗(yàn)?zāi)B(tài)分解算法（EEMD）的HHT算法，對仙女山斷裂周緣的構(gòu)造、礦震、塌陷地震的波形特征進(jìn)行了系統(tǒng)的分析，從時頻分布和邊際譜分布兩方面提出3種地震的判別依據(jù)。

1 理論和方法

HHT時頻分析方法由經(jīng)驗(yàn)?zāi)B(tài)分解和希爾伯特變換組成，具有能夠穩(wěn)定、自適應(yīng)地處理非線性和非平穩(wěn)信號的特點(diǎn)。其中本文采用的集合經(jīng)驗(yàn)?zāi)B(tài)分解算法（Ensemble Empirical Mode Decomposition，EEMD）是對Huang提出的EMD方法的優(yōu)化改進(jìn)[12]。該方法主要原理是在原始數(shù)據(jù)中加入白噪聲，以此補(bǔ)充EMD分解過程中尺度不足的問題;不同時間尺度的信號可以自動分離到與其相適應(yīng)的參考尺度上去，從而有效改進(jìn)EMD分解過程中的模態(tài)混疊現(xiàn)象，得到一系列的IMF和趨勢項(xiàng)。此外，該方法同時保留了EMD方法簡單高效處理非線性、非平穩(wěn)信號的特點(diǎn);每個IMF分量經(jīng)過Hilbert變換，可以獲取信號的瞬時頻率和振幅、邊際譜等數(shù)據(jù)。

1.1 EEMD和希爾伯特變換

在進(jìn)行Hilbert變換之前首先對信號進(jìn)行EEMD分解，其中集合經(jīng)驗(yàn)?zāi)B(tài)分解算法（EEMD）的分解流程如圖2所示，具體步驟包括：

（1） 將正態(tài)分布的白噪聲加到原始信號x（t）;

（2） 將加入白噪聲εμi（t）的信號作為一個整體，進(jìn)行EMD分解，得到各IMF分量cij（t）;

（3） 重復(fù)步驟（1）和（2），每次加入新的正態(tài)分布白噪聲序列;

（4） 將每次得到的IMF分量cij（t）做集成平均處理后作為最終結(jié)果。

各個單分量IMF信號cj（t）由經(jīng)過經(jīng)驗(yàn)?zāi)B(tài)分解后不同時間尺度的信號組成，且最終分解結(jié)果還包含一個信號趨勢項(xiàng)r（t）。具體實(shí)現(xiàn)步驟如下。

式中：x（t）代表原始信號;εμi（t）代表加入的白噪聲信號，其中ε代表加入白噪聲的噪聲水平;xi（t）是加入噪聲后的信號，i=1，2，…，n。

對加入白噪聲后的信號xi（t）進(jìn)行EMD分解，得到m個IMF分量cij（t）和趨勢項(xiàng)ri（t），其中j=1，2，…，m，如式（2）所示：

最后通過集合平均，分別得到原始信號x（t）的分解結(jié)果cj（t）和趨勢項(xiàng)r（t），如式（3）～（5）所示。

通過EEMD分解得到一系列單分量IMF信號之后，對各個IMF分量進(jìn)行信號的Hilbert變換，得到H（t），如式（6）所示，其中αj為第j個IMF的單分量信號的幅值，ω為角頻率。

對H（t）信號進(jìn)行實(shí)部展開，即得到Hilbert譜H（ω，t），Hilbert譜的頻率是瞬時頻率，是隨時間變化的函數(shù)，其表示為

對Hilbert譜H（ω，t）在時間尺度t上進(jìn)行積分，如式（8）所示，可以計(jì)算得到Hilbert邊際譜h（ω），其幅值可以真實(shí)地反應(yīng)頻率在信號中是否存在[13]，并且可以精細(xì)地表示信號幅值在整個頻率段上隨頻率的變化情況。

1.2 最小二乘算法

最小二乘法的基本理論是：給定數(shù)據(jù)集{（yk， zk）}（k=0，1，2，…，p），求取合適的函數(shù)f（y），使誤差的平方和E2最小，其中E2由下面公式（9）得到，其主要數(shù)學(xué)意義即把擬合值和觀察值之間的最小化距離平方和作為目標(biāo)來優(yōu)化。

通過最小二乘法得到擬合函數(shù)f（y），可以有效表征原始數(shù)據(jù)集的變化趨勢和變化特征。

2 數(shù)據(jù)選取

為了最大限度減小地震波在傳播過程中受到介質(zhì)、幾何擴(kuò)散、頻散和儀器等因素的影響[14]，盡量選取同一地震臺記錄到的地震事件，同時考慮到礦震的震級大部分較小，所以針對構(gòu)造、礦震、塌陷地震均選取ML 2.0～3.0級左右地震進(jìn)行分析。

基于上述原則，選取發(fā)生于2007～2016年之間的仙女山斷裂南麓周坪地震臺的數(shù)字化觀測波形事件作為研究對象，其中構(gòu)造、礦震、塌陷地震事件各10次，相關(guān)地震資料已經(jīng)過地震波資料分析和現(xiàn)場實(shí)地考察確認(rèn)，震中位置和周坪地震臺距離Δ≤15 km，震級分布于ML 2.0～3.0級之間，采樣頻率100 Hz，頻帶范圍0～40 Hz。

依據(jù)地震事件的發(fā)震時刻，所有事件選取三分向中垂直向波形數(shù)據(jù)，截取P波到時前1～2 s作為波形分析的開始，數(shù)值均采用未經(jīng)過電壓轉(zhuǎn)換的counts值，利用python腳本程序?qū)υ嫉牟ㄐ螖?shù)據(jù)進(jìn)行處理，去除儀器響應(yīng)部分，并進(jìn)行去均值和去傾斜校正，然后選取30 s左右的時長進(jìn)行分析，選取的數(shù)據(jù)充分包含了地震波形的震蕩和衰減信息。

3 頻譜特征分析

3.1 時頻分析

本文對10次構(gòu)造地震、礦震、塌陷地震事件的EEMD分解和Hilbert時頻圖譜進(jìn)行了計(jì)算分析，時頻圖譜主要反映頻率成分隨時間的變化特征，這里僅選取3種地震各2次事件做詳細(xì)說明。其中構(gòu)造地震事件1（Tectonic_EQ1），震級為ML 2.1級，EEMD分解結(jié)果和時頻圖如圖3（a）～（b）所示，構(gòu)造地震事件2（Tectonic_EQ2），震級為ML 2.0級，EEMD分解結(jié)果和時頻圖如圖3（c）～（d）所示;礦震事件1（Mine_EQ1），震級為ML 2.2級，EEMD分解結(jié)果和時頻圖如圖4（a）～（b）所示，礦震事件2（Mine_EQ2），震級為ML 2.1級，EEMD分解結(jié)果和時頻圖如圖4（c）～（d）所示;塌陷地震事件1（Collapse_EQ1），震級為ML 2.6級，EEMD分解結(jié)果和時頻圖如圖5（a）～（b）所示，塌陷地震事件2（Collapse_EQ2），震級為ML 2.8級，EEMD分解結(jié)果和時頻圖如圖5（c）～（d）所示。上圖中，Imf1～5代表1～5階IMF分量，Res代表殘差。

通過對比分析3種類型地震事件的原始波曲線可知，波形曲線差異比較大，對于不同地震事件類型，Pg波到達(dá)后均有明顯的不同跳變出現(xiàn)，這和不同的發(fā)震因素有關(guān);同種類型的地震事件， Pg波初至之后其波形的形態(tài)基本一致。不同地震事件分解后各IMF分量的特征如圖6所示，其中對于相同類型地震事件，其IMF分量的最大振幅曲線走勢基本一致;最大振幅對應(yīng)時間和各個分量占原始信號的百分比因發(fā)震因素不同而規(guī)律性不明顯，但是不同地震事件Imf1最大振幅對應(yīng)的頻率差別明顯：構(gòu)造地震事件頻率分布在15～25 Hz之間，礦震事件頻率分布在5～15 Hz之間，塌陷事件頻率分布在5 Hz左右。該分析結(jié)果也說明了不同地震事件類型在頻率特征上存在明顯差異。

3種地震類型的EEMD分解結(jié)果均按照原始數(shù)據(jù)的尺度，根據(jù)高頻到低頻的順序分解出了IMF分量。同種地震類型的EEMD分解結(jié)果曲線形態(tài)是一致的，這是因?yàn)樵诜纸膺^程中加入了白噪聲，同時這也避免了計(jì)算Hilbert譜過程中的模態(tài)混疊現(xiàn)象。

分析以上3種類型地震的時頻分布圖，Pg波到達(dá)之前的時頻譜具有如下特征：① 構(gòu)造地震在Pg波到達(dá)之前隨著時間變化出現(xiàn)了能量較低但頻帶較寬、成分復(fù)雜的頻率，并且一致延伸到Pg波到達(dá);② 礦震在Pg波到達(dá)之前，隨著時間推移，出現(xiàn)了能量較低但頻帶較窄、成分簡單的頻率，并且在Pg波到達(dá)之前有一個不規(guī)則的低頻帶區(qū)域;③塌陷地震在Pg波到達(dá)之前，只存在極低頻率成分的頻帶區(qū)。

進(jìn)一步分析Pg波到達(dá)之后的時頻譜，構(gòu)造地震在地震波增大和衰減的過程中，0～40 Hz的頻率成分均有分布;礦震在地震波增大和衰減的過程中，0～30 Hz的頻率成分也均有分布，但是在地震波幅值變換最大的區(qū)域，20～30 Hz頻率成分出現(xiàn)了空白區(qū);塌陷地震在地震波增大和衰減的過程中，頻率成分主要集中在0～30 Hz。

通過分析仙女山斷裂周緣3種地震的時頻圖分得出，在Pg波到達(dá)前后，不同地震類型的頻譜存在明顯差異，主要是因?yàn)榘l(fā)震原因的不同：①構(gòu)造地震由地殼介質(zhì)中巖石破裂引起，是巖層中應(yīng)力從積累到釋放的過程，其震源機(jī)制較為復(fù)雜，因而地震波頻率成分復(fù)雜，能量釋放緩慢，震動持續(xù)時間較長，這與分析的構(gòu)造地震頻帶寬、頻率成分分布均勻是一致的;②礦震是人工震源，震源較淺，震動幅值大，能量衰減快，傳播介質(zhì)疏松，高頻成分易被吸收，所以時頻譜圖在最大震動幅值前后出現(xiàn)了高頻成分的空白區(qū);③塌陷地震是由于塌陷體墜落撞擊底部巖體而導(dǎo)致的瞬時壓縮源，其地震波傳播路徑簡單，從圖5（b）， （d）的時頻譜圖分析可知，與構(gòu)造地震和礦震相比，塌陷地震的頻帶范圍窄、高頻成分少。

3.2 邊際譜分析

采用公式（8）中的計(jì)算方法，對3.1節(jié)中詳細(xì)分析的構(gòu)造地震、礦震、塌陷地震各10次事件進(jìn)行計(jì)算得到原始地震信號的邊際譜，并采用公式（9）對邊際譜的計(jì)算結(jié)果做最小二乘擬合。選取與3.1節(jié)中進(jìn)行時頻分析的6次事件相對應(yīng)的計(jì)算結(jié)果進(jìn)行如下分析說明。

圖7是兩次構(gòu)造地震的邊際譜及擬合曲線，從圖中可以直觀看到邊際譜集中在0～20 Hz的范圍，擬合曲線呈現(xiàn)出逐漸下降的趨勢，這與構(gòu)造地震能量緩慢釋放的原因是一致的。礦震事件的邊際譜及擬合曲線如圖8所示，礦震有一個明顯的主頻帶，集中在0～5 Hz，并且存在一個2 Hz左右的主頻。分析圖9所示塌陷地震的邊際譜可知，塌陷地震也存在一個主頻帶，集中在0～3 Hz，主頻略低于礦震，在1 Hz左右。上述3種地震的邊際譜，其頻帶分布和主頻均存在明顯差異，主要由于不同地震類型的發(fā)震因素和傳播路徑不同引起，而邊際譜分析結(jié)果準(zhǔn)確地反映出了不同地震類型的差異，所以上述邊際譜分布特點(diǎn)可以作為識別仙女山斷裂周緣3種地震類型的重要依據(jù)。

4 結(jié) 論

通過基于EEMD的HHT時頻譜和邊際譜算法，對三峽庫區(qū)周坪地震臺記錄到的仙女山斷裂周緣大量的構(gòu)造、礦震、塌陷地震波形的時頻譜和邊際譜進(jìn)行分析，發(fā)現(xiàn)3種地震類型的時頻分布和邊際譜特征存在明顯差異，并且主要集中以下幾個方面。

（1）在Pg波到達(dá)之前，從時頻譜分析，構(gòu)造地震的頻帶較寬，頻率在0～40 Hz均有分布，并且頻率成分復(fù)雜;礦震的頻帶較窄，基本分布在頻率0～30 Hz之間，頻率成分相對簡單;而塌陷地震頻帶極窄，分布在0～5 Hz間，頻率成分集中在低頻區(qū)域。

（2）從時頻譜分析Pg波到達(dá)后至地震波衰減的時頻分布可知，構(gòu)造地震和礦震在0～40 Hz均有頻率成分分布，但是礦震在地震波幅值變化最大的區(qū)域出現(xiàn)了20～30 Hz頻率成分空白區(qū)，這是構(gòu)造地震和礦震時頻分布的主要差異;而塌陷地震的頻率成分僅分布在0～30 Hz。

（3）從邊際譜擬合結(jié)果分析發(fā)現(xiàn)，構(gòu)造地震邊際譜集中在0～20 Hz，邊際譜曲線在空間上呈現(xiàn)逐漸下降的趨勢;礦震的邊際譜存在一個0～5 Hz的主頻帶和2 Hz左右的主頻;而塌陷地震的主頻帶和主頻低于礦震，集中在0～3 Hz和1 Hz;三者的邊際譜擬合結(jié)果存在明顯差異。

本文利用EEMD方法對三峽庫區(qū)仙女山斷裂周緣的構(gòu)造、礦震、塌陷地震事件波形信號進(jìn)行分解，得到不同時間尺度信號的時頻特征和邊際譜，并基于上述計(jì)算結(jié)果，系統(tǒng)地分析和總結(jié)了該區(qū)域構(gòu)造、礦震、塌陷地震事件的規(guī)律。該研究方法可以作為識別其他研究區(qū)域構(gòu)造、礦震、塌陷地震事件的判別依據(jù)，幫助地震分析人員對實(shí)時記錄的地震波形進(jìn)行快速識別、分類，進(jìn)而提高數(shù)據(jù)分析處理的精度和效率，為抗震救災(zāi)提供技術(shù)支撐;同時，該方法也可以應(yīng)用于堰塞湖、滑坡、泥石流等災(zāi)害過程所產(chǎn)生地震動信號的處理和分析，能夠與其他研究手段相結(jié)合，最終實(shí)現(xiàn)災(zāi)害過程的重構(gòu)，進(jìn)而有效指導(dǎo)地質(zhì)災(zāi)害的監(jiān)測預(yù)警和快速評估工作。

參考文獻(xiàn)：

[1] 崔鑫，許力生，許忠雄，等. 小地震與人工爆破記錄的時頻分析[J]. 地震工程學(xué)報(bào)，2016，38（1）：72-77.

[2] 靳玉貞，林木金，范曉瑜，等. 山西地區(qū)爆破、塌陷（礦震）特殊地震動特征識別[J]. 地震地磁觀測與研究，2015，36（3）： 64-66.

[3] 毛世榕，管振德，閻春恒. 基于小波包分形和神經(jīng)網(wǎng)絡(luò)的地震與巖溶塌陷識別[J].地震學(xué)報(bào)，2018，40（2）： 195-204.

[4] 曾憲偉，趙為明，盛菊琴，等.應(yīng)用小波包識別寧夏及鄰區(qū)的地震和爆破[J].地震研究，2008（2）： 142-148，198.

[5] 黃漢明，邊銀菊，盧世軍，等.天然地震與人工爆破的波形小波特征研究[J].地震學(xué)報(bào)，2010，32（3）：270-276.

[6] SANJIT K， DASH G， SASIBHUSAN R， Arrhythmia detection using Wigner-Ville distribution based neural network[J].? IEEE Transaction on Information， 2016， 85（5）： 806-811.

[7] HUANG N E， SHEN Z， LONG S R， et al. A confidence limit for the empirical mode decomposition and Hilbert spectral analysis[J].? ?Proceedings of the Royal Society of London， 2003， 4510（2037）： 2317-2345.

[8] CHEN Y， ZHANG G， GAN S， et al. Enhancing seismic reflections using empirical mode decomposition in the flattened domain[J].Journal of Applied Geophysics， 2015， 56（71）： 99-105.

[9] LI H， LI Z， MO W. A time varying filter Empirical Mode Decomposition[J]. Signal Processing，2017，138（12）： 146-158.

[10] JAN N Y， YING L，? HUANG N E， et al. System identification of linear structures based on Hilbert-Huang spectral analysis. Part 1：normal modes[J].? Earthquake Engineering & Structural Dynamics，2016，87（32）： 1443-1467.

[11] 楊培杰，印興耀，張廣智. 希爾伯特-黃變換地震信號時頻分析與屬性提取[J]. 地球物理學(xué)進(jìn)展， 2007， 22（5）： 1585-1590.

[12] WANG Y H， YEH C H， Young H W V， etc. On the computational complexity of the empirical mode decomposition algorithm[J]. Physica A，2014（400）：159-167.

[13] 鐘佑明，秦樹人，湯寶平. 希爾伯特-黃變換中邊際譜研究[J]. 系統(tǒng)工程與電子技術(shù)，2004，26（9）：1324-1326.

[14] 張帆，朱新運(yùn)，熊丹，等.基于非線性時頻分析的地震和爆破識別[J]. 華南地震，2014，34（2）：56-63.

（編輯：高小雲(yún)）