尾波干涉技術在2017年重慶武隆MS 5.0地震后震源區地殼介質變化的應用

汪 建 向 涯 楊亞運 李翠平

1）中國重慶 400020 重慶市地震局

2）中國武漢 430071 中國地震局地震研究所

3）中國武漢 430071 中國地震局地震大地測量重點實驗室

0 引言

地震波在地球介質中經多次散射形成尾波。尾波在地球介質中傳播時間長，由于重復采樣、疊加和放大作用，其對介質的微小變化更加敏感，能識別直達波不能識別的介質的細微變化（張金川等，2014）。基于該特性，可利用尾波干涉方法識別震源、散射體位置、地震波速度與介質的細微變化，并廣泛應用于斷層活動性、工程材料檢測、冰川活動、火山監測等領域（肖卓等，2015）。近年來，諸多地震學者利用尾波干涉方法開展了一系列研究，如：宋麗莉等（2012）通過室內試驗，利用尾波干涉方法觀測到直達波難以分辨的巖石物性變化；張金川（2014）利用尾波干涉技術研究廣西龍灘水庫重復地震，發現散射體運移可能與庫水向地下裂隙中的運移有關；肖卓等（2016）運用尾波干涉技術分析2014 年云南盈江雙震期間地殼介質狀態的變化；李樂等（2017）利用重復地震尾波研究汶川強震破裂區震后波速的微變化；鮑子文（2019）利用重復地震尾波研究安徽霍山地區介質變化，發現地震波速度存在季節性變化，與大氣降水有較好的相關性。

2017 年11 月23 日17 時43 分，重慶武隆區發生MS5.0 地震。該地震事件是有地震儀器記錄以來重慶渝東南地區發生的最大地震。文中利用尾波干涉原理，研究武隆MS5.0地震發生后，余震期間震源區地下介質的變化，有助于了解強震發生時地殼介質性質和應力狀態的變化。

1 資料選取

武隆區地處重慶市東南部，位于烏江下游峽谷區，巖石成分以碳酸鹽、砂巖和頁巖為主，巖溶地貌發育，形成大量溶洞群、天坑群、地縫等自然景觀。該地區附近分布方斗山基底斷裂、七曜山—金佛山基底斷裂及彭水基底斷裂，中小斷裂有三會沖斷裂、芙蓉江斷裂、文復斷裂、馬武斷裂、火石埡斷裂、郁山斷裂、和尚巖斷裂（戴應洪等，2020）（圖1）。李翠平等（2019）分析認為，震中附近流體的深部滲透作用導致局部斷層活化，誘發此次武隆MS5.0 地震，發震斷層與文復斷裂關系較為密切。

圖1 研究區地震臺站、震中及區域構造Fig.1 Seismic stations,epicenters,and regional structure in the study area

武隆MS5.0 地震發生后，在震中附近架設3 套流動臺（L5541、L5550、L5551），均配備短周期地震計（采樣率100 Hz），進行地震監測。各流動臺臺基噪聲測試參數見表1，可見環境地噪聲水平均在Ⅱ級或以上。截至12 月23 日，重慶數字地震臺網共記錄到武隆及鄰區（29°00′—29°42′N，107°12′—108°24′E）發生ML≥1.0 地震108 次（含MS5.0 主震），最大余震為2017 年11 月29 日10 時34 分武隆ML3.1 地震。選取108 次ML≥1.0 地震波形數據和觀測報告，分析地殼介質變化。

表1 流動臺站臺基噪聲參數Table 1 Noise parameters of three mobile stations

2 重復地震挑選

重復地震概念于20 世紀60 年代被提出，目前尚無統一定義（李宇彤等，2011），Nadeau 等（1995）將發生在同一斷層位置，震級和復發間隔相近，震源機制與波形高度相似的一組地震稱為重復地震。Schaff等（2004）則將重復地震定義為被至少一個臺站記錄到，且波形相關系數不小于0.8的地震對事件。文中通過計算波形互相關系數（cross-correlation簡稱cc）和約束震源位置的方法，來挑選武隆MS5.0 地震發生后的重復地震。

2.1 地震序列重定位

汪建等（2020）利用武隆地區2010—2017 年重復地震和觀測報告資料，發現該區地震定位和震相拾取精度較高。為得到更精確的震源位置，利用雙差定位法，采用武隆地區精細地殼速度結構模型（表2）（王小龍等，2013），對所選108 次地震序列進行重定位。雙差定位結果顯示，MS5.0 主震發生在三會沖、芙蓉江、文復斷裂與馬武斷裂交匯部位，震中位于烏江上游彭水電站與下游銀盤電站水位段之間，與烏江水系最近距離約2 km，震源深度約10 km，余震序列呈條帶分布，與附近斷裂走向平行，地震序列破裂面呈NE—SW向單側破裂（圖1）。

表2 武隆地區速度結構模型（據王小龍等，2013）Table 2 The seismic velocity model used in the location (According to Wang Xiaolong et al,2013)

2.2 重復地震挑選

將地震波形轉換為SAC 格式，選取垂直分量，進行去傾斜、去線性趨勢及0.5—10 Hz 帶通濾波預處理（汪建等，2018），挑選重復地震，步驟如下：①計算同一臺站中任意2個地震波形的互相關系數，計算起點為P 波到時前1 s，計算長度選擇4 倍Sg、Pg 震相走時差；②找出同時被3 個臺站（L5541、L5550、L5551）記錄且各臺互相關系數（cc）≥0.9的地震對事件作為相似地震；③結合雙差定位結果，尋找間距（不考慮震源深度的影響）小于1 km 的相似地震。

經挑選，得到3 組滿足條件的地震對，均包含2 個地震事件。重復地震以D1、D2、D3編號，主震設為D0。重復地震與流動臺站分布見圖2，可見D1、D2、D3均分布在主震D0南側，由北向南依次擴展，結合發震時刻可知，重復地震對均發生在余震活動相對減少階段。

圖2 臺站與重復地震分布Fig.2 Distribution of seismic stations and repeating earthquakes

對武隆MS5.0 地震序列重復地震對參數信息進行統計，結果見表3，可知：D1、D2、D3的震級均小于ML2.0，震級差值最大ML0.5，最小ML0.3；震源深度差值最大2 km，最小1 km；地震對重復周期具有明顯變化，時間間隔從十幾小時至幾天不等，最大間隔近4 天，D1、D2時間間隔相近。

表3 武隆MS 5.0 地震及重復地震對的參數統計Table 3 Parameter statistics of the Wulong MS 5.0 earthquake and repeating earthquakes

3 尾波干涉分析

3.1 尾波干涉原理

Snieder（2002）基于前人研究提出尾波干涉方法，并對其原理進行了系統闡述。尾波干涉法基于路徑疊加理論，認為波場是沿可能散射路徑傳播的子波的疊加，重復地震與臺站之間的傳播路徑大致相同，波形差異由介質變化所引起，重復地震尾波的變化反映了以臺站和重復地震為焦點的橢球體上介質的變化，可以利用重復地震尾波研究地殼介質的波速變化，進而分析引起波速變化的原因（肖卓等，2016）。

波形數據原始采樣率為100 Hz，即兩列波形走時延遲最小為0.01 s。為獲得更小的走時延遲，在尾波干涉測量之前，采用插值法將原始信號的采樣率提升至10 000 Hz（Peng et al，2006），即兩列波形的走時延遲最小為0.000 1 s。當波速在介質中發生變化時，重復地震的尾波穿過該介質時也將發生走時擾動，在接收點表現為波形的走時偏移。文中采用移動窗口互相關法（Niu et al，2003；肖卓等，2016）計算重復地震之間的走時偏移，即將發震時刻較早的重復地震當作參考地震，從P 波初至到時起算，對重復地震進行尾波干涉分析。通過移動時間窗，計算每一個時間窗內波形的走時延遲，得到走時延遲及不相關系數隨流逝時間的分布。結合前人經驗（張金川，2014；肖卓等，2016），選取P 波初至到時至S 波到時的2 倍加4 s 作為尾波干涉測量窗口，移動窗口長度為1 s，移動步長為0.05 s，即每1 s 內有20 個干涉測量點。

假設介質中波速變化是均勻的，若波速下降（增大），則兩列波形的走時偏移將隨著流逝時間線性上升（減?。⊿nieder，2006）。設原始波速為v，介質整體發生大小為δv的波速變化，則走時偏移τ與波速變化δv的對應關系（Snieder，2006；Payan et al，2011）為

式中，走時偏移τ與流逝時間t經線性擬合得到的直線斜率即為介質波速的相對變化。在計算每小時時間窗j內波形的走時偏移τj時，其測量誤差估計（Liu et al，2010）為

式中，tj為小時間窗j的中心位置，M為時間窗的數量，ητ為擬合殘差，表示為

測量誤差下限στ需滿足Cramer-Rao Lower Bound 法則（Quazi，1981；Carter，1987），公式如下

式中，f0為信號主頻，B為信號的頻寬與主頻之比，T為窗口長度，ρ為波形相關系數，SNR 為信噪比。f0和B由波形數據控制，測量誤差主要取決于互相關系數與信噪比，選取互相關系數及信噪比較高的尾波段進行尾波干涉測量較為合理（肖卓等，2016）。

3.2 尾波干涉結果

利用上述原理，采用尾波干涉技術，分別對重復地震對D1、D2、D3在3 個流動臺的波形記錄進行尾波干涉處理。

3.2.1 L5541 臺站記錄尾波干涉處理。分別對重復地震對D1、D2、D3在L5541 臺的記錄進行尾波干涉處理，結果見圖3、圖4、圖5（0 s 設為P 波初至到時）。

圖3 L5541 臺站D1 波形記錄及其不相關系數和走時延遲隨流逝時間的變化Fig.3 The seismic waveforms,uncorrelation coefficient and delay time of D1 repeating earthquakes with time at the L5541 station

圖4 L5541 臺站D2 波形記錄及其不相關系數和走時延遲隨流逝時間的變化Fig.4 The seismic waveforms,uncorrelation coefficient and delay time of D2 repeating earthquakes with time at the L5541 station

圖5 L5541 臺站D3 波形記錄及其不相關系數和走時延遲隨流逝時間的變化Fig.5 The seismic waveforms,uncorrelation coefficient and delay time of D3 repeating earthquake with time at the L5541 station

（1）重復地震對D1。D1的2 次地震分別發生在主震D0后第7 和第8 天，分析L5541 臺站的2 次地震預處理波形及其不相關系數和走時延遲隨流逝時間（從P 波初至到時起算）的變化，結果見圖3。重復地震對D1的2 次地震事件震級相差不大（0.5），但其相關系數卻達0.937。由圖3 可見：①在P 波后續部分（1.00—1.95 s），走時延遲基本保持不變，后迅速下降至零值附近波動；②在S 波振幅較大部分（1.95—3.20 s），走時延遲緩慢下降，S 波后續部分（3.20—4.70 s）的走時延遲呈先增加后減小的變化趨勢；③隨著流逝時間的增加，S 波早期尾波段（5.25—6.75 s）走時延遲明顯增加，意味著波速的減小，S 波尾波后續部分（6.75—7.75 s）的走時延遲明顯減小，最大走時延遲約65×10-3s；④結合移動窗口內波形的信噪比、不相關系數和走時延遲的線性變化規律，選取S 波早期尾波段（5.25—6.75 s），由式（1）計算得到，其相對波速變化約為2.1‰。

（2）重復地震對D2。D2的2 次地震分別發生在主震D0后第8 和第9 天，圖4 給出重復地震對D2在L5541 臺站的尾波干涉結果。由圖4 可見：①在P 波后續部分（1.00—1.55 s），走時延遲迅速下降至零值附近波動；②在S 波振幅最大部分（1.55—2.40 s），走時延遲基本保持不變，隨著流逝時間的增加，S 波后續部分（2.50—3.75 s）的走時延遲有緩慢增加的趨勢；③在S 波早期尾波段（3.80—4.90 s），走時延遲明顯增加，變化較為線性，S 波尾波后續部分（4.90—6.35 s）的走時延遲明顯減小，最大走時延遲約19.1×10-3s；④選取S 波早期尾波段（3.80—4.90 s）計算其相對波速變化，結果約為13.4‰，波速變化較大。

重復地震對D2距L5541 臺站最近，傳播路徑不穿過烏江流域，可能受到斷層擠壓作用，巖石裂隙增多，導致地下介質劇烈變化，波速下降值較大。若選取S 波后續部分（2.50—3.75 s）計算其相對波速變化，結果約1.4‰，約為S 波早期尾波段計算結果的1/10，正是由于尾波的多次散射、疊加、放大，以及比S 波更敏感的特點，尾波識別的波速變化比S 波放大近10 倍。

（3）重復地震對D3。D3的2 次地震分別發生在主震D0后第9 和第14 天，圖5 給出重復地震對D3在L5541 臺站的尾波干涉結果。由圖5 可見：①S 波衰減較快，在S 波部分（2.35—3.25 s）走時延遲緩慢增加，在S波早期尾波段（3.25—4.60 s）走時延遲增加較為明顯，S 波后續尾波部分（4.60—7.35 s）走時延遲變化復雜，整體呈下降趨勢。選取S 波早期尾波段（3.25—4.60 s），計算得到其相對波速變化約為2.2‰；②在P波尾波部分（1.25—1.45 s）發現走時延遲和不相關系數均存在1 個短周期“尖峰”（spike）變化，且變化比較尖銳，相比于P 波尾波，S 波及其尾波變化更為平滑。

Niu 等（2003）在美國圣安德列斯斷層帕克菲爾德區域，對重復微震進行數值模擬，發現若只是單一散射體的位置發生變化，走時擾動與不相關系數均將出現一個短周期的“尖峰”變化，若震源位置或介質波速發生變化，隨著時間的流逝，走時擾動與不相關系數將呈現長周期變化趨勢，根據實驗結果，P 波尾波的“尖峰”變化可能與地下介質中散射體的運移有關。由圖2 可知，重復地震對D3到L5541 臺站的傳播路徑穿過烏江流域。該區巖溶發育，遍布溶洞和地下暗河，因S 波不能在液體中傳播，P 波尾波出現的“尖峰”變化可能與震后庫水在地下裂隙中的運移有關。

3.2.2 L5550 臺站記錄尾波干涉處理。分別對重復地震對D1、D2、D3在L5541 臺的記錄進行尾波干涉處理，可知D1在該臺站尾波干涉結果未見明顯變化，D2、D3的尾波干涉結果見圖6、圖7。

圖6 L5550 臺站D2 波形記錄及其不相關系數和走時延遲隨流逝時間的變化Fig.6 The seismic waveforms,uncorrelation coefficient and delay time of D2 repeating earthquakes with time at the L5550 station

圖7 L5550 臺站D3 波形記錄及其不相關系數和走時延遲隨流逝時間的變化Fig.7 The seismic waveforms,uncorrelation coefficient and delay time of D3 repeating earthquakes with time at the L5550 station

（1）重復地震對D2。圖6 給出重復地震對D2在L5550 臺站的尾波干涉結果。由圖6可見：①在P 波尾波部分（1.00—2.55 s）走時延遲無明顯變化；②走時延遲在S 波早期部分（2.55—4.00 s）減小，后續部分（4.00—6.05 s）增加，在S 波早期尾波段（6.05—7.25 s）明顯減小，意味著波速的增加，在S 波后續尾波部分（7.25—9.40 s）保持在零值附近波動，最大走時延遲約為6.2×10-3s；③選取S波早期尾波段（6.05—7.25 s）計算其相對波速變化，結果約為6.7‰，可能是由于此段波形相關系數不高，由式（2）—（3）計算，測量誤差與波速變化在同一個數量級。

（2）重復地震對D3。圖7給出重復地震對D3在L5550臺站的尾波干涉結果。由圖7可見：①P 波尾波部分（1.00—2.70 s）的相關系數較低，走時延遲先增加后減??；②在S 波部分（2.70—5.65 s），走時延遲變化復雜，呈減小—增加—減小的變化趨勢，S 波早期尾波段（5.75—6.85 s）走時延遲明顯增加，變化較為線性，S 波后續尾波部分（6.85—8.25 s）走時延遲下降至零值保持不變；③在不同波段（P、P 尾波、S、S 尾波），走時延遲的變化趨勢不同且復雜，但在S 波早期尾波段（5.65—6.85 s），走時延遲有明顯增加，且較為線性，最大走時延遲約為9.5×10-3s，其相對波速變化約為6.1‰。

3.2.3 L5551 臺站記錄尾波干涉處理。分別對重復地震對D1、D2、D3在L5551 臺的記錄進行尾波干涉處理，可知D3在該臺尾波干涉結果未見明顯變化，D1、D2的尾波干涉結果見圖8、圖9。

（1）重復地震對D1。圖8 給出重復地震對D1在L5551 臺站的尾波干涉結果。由圖8可見：①在P 波尾波部分（1.35—1.75 s），發現走時延遲和不相關系數均存在2 個短周期“尖峰”變化，S 波及其尾波的走時延遲也有數個“尖峰”變化，但其最大變化幅度約為P 波的1/3，根據Niu 等（2003）的實驗結果，此現象可能與震源位置或介質波速變化有關，若震源位置發生變化，直達波與尾波變化的量級應基本一致（Grêt et al，2006），假設由于介質波速變化導致P 波尾波、S 波及其尾波出現數個“尖峰”變化，則可能與地下介質的不均勻性或復雜性有關；②在S 波及其尾波段（2.10—7.40 s），走時延遲整體呈上升趨勢；③選取S 波早期尾波段（4.00—6.05 s）計算其相對波速變化，結果約為2.1‰。

圖8 L5551 臺站D1 波形記錄及其不相關系數和走時延遲隨流逝時間的變化Fig.8 The seismic waveforms,uncorrelation coefficient,and delay time of D1 repeating earthquakes with time at the L5551 station

（2）重復地震對D2。圖9 給出重復地震對D2在L5551 臺站的尾波干涉結果。由圖9 可見：①P 波尾波部分（1.00—2.10 s）相關系數不高，走時延遲變化較為復雜；②S 波振幅較大部分（2.10—4.65 s）走時延遲基本保持不變，在S 波早期尾波段（4.65—5.95 s）走時延遲明顯減小，意味著波速的增加，在S 波后續尾波部分（5.95—7.50 s）走時延遲保持在零值附近波動；③選取S 波早期尾波段（4.65—5.95 s）計算其相對波速變化，結果約為4.5‰。

圖9 L5551 臺站D2 波形記錄及其不相關系數和走時延遲隨流逝時間的變化Fig.9 The seismic waveforms,uncorrelation coefficient,and delay time of D2 repeating earthquakes with time at the L5551 station

4 結論與討論

綜上，絕大多數的測量誤差比波速變化測量值小1 個數量級，所得波速變化結果較為可靠。圖2 可見，在L5541 臺站（距震源區最近），3 組重復地震記錄的S 波尾波部分走時延遲均有明顯增加現象，可能意味著地震波速的持續減小，其中D1和D3的波速相對變化小于3‰，D2的波速相對變化最大，約為13.4‰，意味著震源區附近巖石遭受破壞，結構弱化，應力得到較大釋放；在D3的P 波尾波部分發現了走時延遲和不相關系數均存在1個短周期“尖峰”變化，可能與庫水在地下裂隙中的運移有關。在L5550 臺站，D2的S 波尾波部分走時延遲明顯減小，波速上升約6.7‰；D3的S 波尾波部分走時延遲明顯增加，波速下降約6.1‰，二者得到的波速變化值接近，可能表明震后局部地區地震波速度存在上升—下降的恢復過程，該階段應力狀態可能處于由累積到釋放的快速調整過程。在L5551 臺站，D1的S 波尾波部分走時延遲明顯增加，波速下降約2.1‰；D2的S 波尾波部分走時延遲明顯減小，波速上升約4.5‰。

選取布設在武隆MS5.0 地震震中附近3 個流動臺的地震事件波形和觀測報告，利用波形互相關技術，結合雙差定位結果，選出3 組符合條件的重復地震。利用尾波干涉原理，通過移動窗口互相關方法，計算了每組重復地震之間的走時延遲變化。尾波干涉結果顯示，在P 波的尾波部分，走時延遲變化較為復雜，可能與P 波和S 波不同的變化規律和PS 轉換波復雜的震相有關。在S 波振幅較大部分，相關系數較高，走時延遲變化不大，可能是由于S 波在重復地震識別過程中參與互相關計算的權重較高。在S 波早期的尾波部分，走時延遲變化最為明顯且較為線性，可能是由地殼介質的地震波速變化引起的。尾波干涉結果反映的是整個傳播路徑上介質變化的平均結果，距離震源區較遠的XNS、LAX、XIT 等臺站，走時延遲變化并不明顯，說明震源區外的地殼介質變化微弱。

選擇流動臺站的數據進行尾波干涉分析，優勢在于流動臺站距震中位置足夠近，布局更加合理，且識別出的重復地震數量優于固定臺站，更能捕捉到地下介質的細微變化。因武隆MS5.0 地震前近450 天研究區內無地震記錄，不能通過尾波干涉方法得到震前地下介質變化信息。

感謝江西省地震局查小惠工程師為本研究提供波形互相關和尾波干涉程序。