三峽庫區仙女山斷裂周緣地震類型特征識別

徐玉聰 朱建 董建輝 周魯

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

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

中圖法分類號：P315 文獻標志碼：A DOI：10.15974/j.cnki.slsdkb.2022.02.002

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

0 引 言

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

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

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

本文利用基于集合經驗模態分解算法（EEMD）的HHT算法，對仙女山斷裂周緣的構造、礦震、塌陷地震的波形特征進行了系統的分析，從時頻分布和邊際譜分布兩方面提出3種地震的判別依據。

1 理論和方法

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

1.1 EEMD和希爾伯特變換

在進行Hilbert變換之前首先對信號進行EEMD分解，其中集合經驗模態分解算法（EEMD）的分解流程如圖2所示，具體步驟包括：

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

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

（3） 重復步驟（1）和（2），每次加入新的正態分布白噪聲序列;

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

各個單分量IMF信號cj（t）由經過經驗模態分解后不同時間尺度的信號組成，且最終分解結果還包含一個信號趨勢項r（t）。具體實現步驟如下。

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

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

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

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

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

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

1.2 最小二乘算法

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

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

2 數據選取

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

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

依據地震事件的發震時刻，所有事件選取三分向中垂直向波形數據，截取P波到時前1～2 s作為波形分析的開始，數值均采用未經過電壓轉換的counts值，利用python腳本程序對原始的波形數據進行處理，去除儀器響應部分，并進行去均值和去傾斜校正，然后選取30 s左右的時長進行分析，選取的數據充分包含了地震波形的震蕩和衰減信息。

3 頻譜特征分析

3.1 時頻分析

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

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

3種地震類型的EEMD分解結果均按照原始數據的尺度，根據高頻到低頻的順序分解出了IMF分量。同種地震類型的EEMD分解結果曲線形態是一致的，這是因為在分解過程中加入了白噪聲，同時這也避免了計算Hilbert譜過程中的模態混疊現象。

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

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

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

3.2 邊際譜分析

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

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

4 結 論

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

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

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

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

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

參考文獻：

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

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

[3] 毛世榕，管振德，閻春恒. 基于小波包分形和神經網絡的地震與巖溶塌陷識別[J].地震學報，2018，40（2）： 195-204.

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

[5] 黃漢明，邊銀菊，盧世軍，等.天然地震與人工爆破的波形小波特征研究[J].地震學報，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]. 地球物理學進展， 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]. 系統工程與電子技術，2004，26（9）：1324-1326.

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

（編輯：高小雲）