基于區域寬頻帶和近震波形數據的多點源反演
——以2016 MW7.1日本熊本地震為例

施其斌， 韋生吉*

1 南洋理工大學亞洲環境學院， 新加坡 639798 2 南洋理工大學新加坡地球觀測中心， 新加坡 639798

0 引言

地震震源的運動學過程是地震學研究的一個基本內容.隨著地震波觀測數據的質量的提高——包括臺站密度的增加和觀測頻段范圍的拓展等——震源運動學的發展經歷了從簡單點源(例如地震的位置、大小、震級和矩張量解(Aki and Richards, 2002))到更接近真實地震的復雜的有限破裂模型(姚振興和紀晨, 1997; Ji et al., 2002; Xu and Chen, 2002; 許力生和陳運泰, 2004; 張勇等, 2008; Zhang et al., 2009; Yue et al., 2012; 岳漢等, 2019)的進步.點源對于中強度(MW6-7)地震的描述往往過于簡單，然而有限破裂模型又常常用過多的參數來非唯一地表示地震的破裂過程.對后者而言，研究者一般很難估計斷層幾何等參數的誤差.地震的多點源解剛好填補了點源和有限破裂模型之間的空缺，彌補了上述兩者的不足之處.通常情況下，多點源解所適用的反演頻段也介于單點源和有限破裂模型的頻段之間.在此頻段，多點源解通常只需要用幾個(例如3～8個)雙力偶位錯點源來描述地震的破裂過程，要遠比單矩張量解更能代表地震的復雜性.同時，多點源解的參數空間遠小于有限破裂模型，因而不會過度擬合觀測數據，更易于對震源參數進行誤差估計.近年來，多點源反演已經被多次應用于大地震研究中(Chen et al., 2014; Duputel and Rivera, 2017; Yue and Lay, 2020; Jia et al., 2020)，此類研究通常使用遠震體波和長周期的面波.對于那些密集強震記錄到的中強度地震，多點源反演也獲得了很好的效果 (Shi et al., 2018; Shi and Wei, 2020).然而，對中強度地震的多點源研究還沒有充分利用區域距離(<1000 km)的臺站記錄.實際上，除了少數地震頻發的發達地區(例如日本和美國西海岸)，絕大部分中強度地震頻發的地區并沒有密集的強震臺網.相比之下，區域距離范圍內臺站數量對全球絕大部分地區都更加可觀.因此，是否能利用此類數據對中強度地震的破裂過程進行多點源反演非常值得探究.為此，本文通過對2016MW7.1熊本地震的研究來測試利用區域寬頻帶地震數據進行多點源地震反演.我們將利用前/余震來對區域臺站的不同分量進行路徑校正，以得到適用于主震反演的頻段.然后對比加入了更為密集的強震和高頻GPS數據得到的多點源解和現有的有限破裂模型，來進一步驗證多點源解.

我們改進了Shi等(2018)所提出的基于馬爾可夫鏈蒙特卡洛(Markov Chain Monte Carlo，MCMC)的反演方法.該方法借鑒了前人所提出的“剪切黏貼”(Cut-And-Paste，CAP)理念 (Zhao and Helmberger, 1994; Zhu and Helmberger, 1996)，將每個臺站的三分量地震波數據截成Pnl波和面波部分，在使用理論地震波擬合數據時允許對不同分量使用不同的時移進行擬合，然后計算時移后的波形擬合誤差.Shi等(2018)采用該方法，利用強震數據研究了2016熊本地震序列的MW6.2前震，結果揭示了發震區斷層幾何形態的復雜性.在本文中，我們將該多點源反演方法應用到包含高頻GPS數據的密集近震觀測，并拓展到對區域地震數據的分析.在此基礎上，提出了使用正則化的Yoffe震源時間函數來代替之前的三角形函數，并且還提出了基于臺站和波形分量來選擇濾波參數，從而最大化地利用地震波形所包含的信息.為了對每個分量選擇合適的濾波頻段，我們利用比研究地震震級更小的地震來進行路徑校正.這些地震的震級應當小到在主震反演的濾波頻段可以被合理地簡化成一個點源，但同時也要大到可以被足夠多的臺站記錄到，不管是強震儀還是寬頻帶地震儀.在后文將給出更多相關的技術細節.

1 2016年熊本地震序列

2016年的MW7.1熊本地震序列發生在日本九州島的Beppu-Shimabara地塹地區，起始于Futagawa斷層和Hinagu斷層相交處.該斷層系統的形成主要是由菲律賓板塊向西北方向略微斜向俯沖所導致，板塊俯沖的速度5～6 cm·a-1，斜沖投影在Futagawa斷層和Hinagu斷層上的滑動速率為7～8 mm·a-1(Nishimura and Hashimoto, 2006)(圖1a).

序列始于4月14日發生的MW6.2前震，約2.5 h后，發生了一個MW6.0的余震，MW7.1主震發生在當地時間4月16日的凌晨，震中位于前震的東北方向(Kato et al., 2016).該地震序列的余震極其活躍，產生了至少30000個Mj(JMA(Japan Meteorological Agency)震級，即日本氣象局測定的地震震級)大于2.5的余震 (Kato et al., 2016).地震序列還造成了復雜的地表破裂，以右旋走滑為主但有明顯的正斷層滑動 (Shirahama et al., 2016)，尤其在兩個斷層相交處和破裂的東北段(圖1b).

該地震序列有著高質量的密集強震觀測、衛星影像觀測以及實地考察記錄，為震源運動學研究提供了寶貴的數據.前人利用地震學數據或者聯合測地學數據對MW6.2 前震(Asano and Iwata, 2016)和MW7.1主震(Asano and Iwata, 2016; Yagi et al., 2016; Hao et al., 2017; Kobayashi, 2017; Yue et al., 2017)進行了有限破裂模型的反演.但是不同研究者所得到的主震破裂模型之間仍存在不小的差異.例如，Asano和Iwata(2016)的結果顯示主震的破裂主要集中在8～10 km 深度，位于震中東北約20 km 處；而Yue等(2017)和Kobayashi(2017)的結果則顯示主震主要的位錯比較淺，已經破裂到了地表.這些有限破裂模型使用2～4個斷層面來近似主震較為復雜的斷層幾何.相對的，前人對MW6.2前震進行的有限破裂模型反演采用了相對簡單的斷層幾何形態，僅為一個矩形的斷層面(Asano and Iwata, 2016).不同于有限破裂模型，Shi等(2018)提出了一種基于MCMC的地震多點源反演方法，用于在相對高頻反演MW6.2前震所產生強震波形.他們使用了一個MW5.4 地震的波形記錄來對所使用的波形分量進行篩選并確定最佳濾波頻段，在此基礎上的多點源反演結果揭示了MW6.2前震較為復雜的斷層幾何，即發現該地震發生在沿相反方向傾斜的多個斷層上.

2 方法和數據

如前所述，Shi等(2018)所提出的多點源反演方法已經在MW6.2前震中得到了成功的應用.該方法盡管對強震數據反演非常有效，但還沒有在高頻GPS數據和區域寬頻帶地震波的反演中嘗試過.為此，本文將該方法拓展到高頻GPS和區域寬頻帶地震記錄來對MW7.1 熊本主震的破裂過程進行研究，并將結果和強震+高頻GPS記錄反演的多點源解以及已經發表的有限破裂模型進行對比，旨在了解該方法在中強度地震研究中的優勢和局限性.我

們還對該方法做了進一步的改進.Shi等(2018)的方法中對所有觀測臺站相同的波形分量采用相同的濾波頻段(例如用0.02～0.2 Hz對Pnl統一進行濾波)，但是由于每個臺站所對應的速度結構不相同，其波形分量對應的最佳擬合頻段也有所不同.這里我們采用了臺站相關的濾波參數，以最大限度地從波形數據中提取地震破裂過程的信息.另外，還改進了反演所用的震源時間函數，不同于之前的等腰三角形函數，這里使用正則化的Yoffe函數(Tinti et al., 2005)來作為震源時間函數(圖2)，以測試所用的數據是否能進一步分辨地震位錯發生的時間過程.實際上Yoffe 函數更符合滑動弱化震源動力學破裂模型中斷層上應力的演化過程(Yoffe, 1951; Tinti et al., 2005)，這對于理解震源破裂物理過程具有重要的意義.新的震源時間函數用α和β兩個參數來定義函數的形狀(圖2)，數學表達如下：

(1)

圖2 正則化Yoffe震源時間函數的數學定義該震源時間函數表示為黑色，由三個對稱軸相同但周期不同的cosine/sine函數組成，分別顯示為紅色、藍色和綠色.Fig.2 The mathematical composition of the reguralized Yoffe-type source time functionThe Yoffe function is shown in black, composed of the three cosine/sine functions with the same symmetric axis but different periods and shown in red, blue and green, respectively.

其中f1(t),f2(t)和f3(t)分別用α和β來表示：

(2)

(3)

(4)

表達式(1)具有更多的靈活性，既可以表示一個對稱的函數(α=0，β=1)，也可以表示非對稱的“長尾”的函數(例如α=0.5，β=0.2).

所使用的波形數據包括強震、高頻GPS和寬頻帶地震數據(圖3).三分量強震數據是從KiK-Net和K-net網站上(https:∥www.kyoshin.bosai.go.jp)直接下載，其原始波形為加速度記錄，將其積分至速度和位移波形用于反演.寬頻帶波形數據是從NIED(https:∥www.fnet.bosai.go.jp/top.php)上下載，去儀器響應后轉化成位移和速度波形兩種記錄.我們還從日本的Nippon GPS Data Service Corporation購買了1 Hz采樣的主震高頻GPS數據.除了主震(MW7.1)，還從這些地震臺網上下載了兩個校正地震事件的波形數據.這兩個校正事件分別發生在Futagawa斷層(MW5.4，32.776°N，130.815°E，2016-04-14-13∶07∶35.290)和Hinagu斷層上(MW6.0，32.7007°N，130.778°E，2016-04-14-15∶03∶46.450).這兩個事件均為MW6.2前震的余震，其破裂過程相對簡單，在相對較高(例如0.3 Hz)的頻段可以看作點源.利用這兩個校正事件，我們確定了所使用的臺站和波形分量，以及臺站相關的濾波參數(例如以80%互相關系數來選擇濾波頻段和分量).最終選取了110個強震臺，包含276個Pnl波形分量和251個面波波形分量，36個寬頻帶臺站，包含90個Pnl波形分量和83個面波波形分量，以及26個高頻GPS臺站.圖3給出了這三個地震以及MW6.2前震在幾個代表性臺站上的0.2 Hz低通濾波濾波波形.可以看到兩個校正事件(MW5.4，MW6.0)的波形相對簡單，而MW6.2前震則要復雜得多，說明其破裂過程更加復雜.值得指出的是，在KMM005臺站上，MW6.2前震的波形甚至比MW7.1主震還要復雜.這一方面是因為主震是向東北方向的單側破裂，而KMM005臺站(方位角38°)正好位于破裂的傳播方向上，因此多普勒效應將波形壓縮顯得更為簡單；另一方面是因為主震的震源時間函數更長、更平滑，導致了主震的波形的頻率成分更低一些，看起來沒那么復雜.而背離主震破裂方向的KMM008臺站其波形振幅要低很多，但是持續時間更長，波形也更為復雜.這些都符合主震單側破裂的特征.同時寬頻帶臺站的波形數據則顯示MW7.1主震的波形要比其他地震復雜得多(圖3c，d)，因此將這些數據用于多點源反演可以對其破裂過程進行很好的約束.相對于強震臺站，我們需要選擇更低頻段的濾波才能對區域臺站的校正事件波形實現較高互相關系數的擬合.區域臺站的路徑校正主要是通過MW6.0前震實現的.

圖3 近Futagawa斷層的熊本主震及前震的強震和寬頻數據對比(a)和(b)為強震儀KMM006和KMM005的觀測數據低通濾波后的波形； (c)和(d)為寬頻地震儀TKO和SIB記錄低通濾波后的波形.7.1級主震的波形由藍色特別顯示.由上至下分別為豎直、南北和東西分量的原始地震波形.所有波形均通過了角頻率為0.2 Hz的低通濾波.Fig.3 Strong-motion and broadband observations of the Kumamoto mainshock and foreshocks near the Futagawa fault(a) and (b) show the low-pass filtered waveforms of the strong-motion stations KMM006 and KMM005, respectively. (c) and (d) show the low-pass filtered waveforms of the broadband stations TKO and SIB, respectively. The mainshock waveforms are highlighted in blue. The vertical, north-south and east-west waveforms are arranged from the top to the bottom. All waveforms are low-pass filtered with the corner frequency of 0.2 Hz.

3 反演結果

我們使用經過路徑校正得到的波形分量和濾波頻段，來對主震的多點源解進行反演.從單個點源反演開始逐漸增加點源的數量，直到反演得到的波形互相關系數的統計分布非常接近校正事件的波形擬合效果，這樣可以避免對數據的過度擬合.圖4給出了聯合使用強震數據和區域寬頻帶數據，和只使用區域寬頻帶數據得到的多點源反演結果，包括用速度波形和位移波形得到的結果.包含強震數據的反演結果顯示我們需要至少4個子事件來足夠準確的表示主震的運動學破裂過程.我們也將三點源反演的結果和四點源解進行比對(圖5)，發現四點源解包含更多細節.如圖6所示，當只用3個子事件時，得到的互相關系數的統計直方圖和用校正事件得到的直方圖還有明顯的差距，而當用到4個子事件時統計的結果和校正事件更為接近.我們注意到3個和4個子事件的波形擬合結果最明顯的差別集中在離地震最近的幾個強震臺站上，把幾個代表性的臺站在圖7中標注了出來.對這幾個臺站，3個子事件明顯對波形和振幅的擬合都差于4個子事件的解.總體而言，這兩種數據組合的反演結果非常穩定和可靠，這點從反演的誤差分布可以看出來(圖8).包含強震數據的反演所對應的誤差分析和統計顯示，各個子事件的走向、傾角和滑動角的誤差都在5°以內，震級的誤差小于0.1，深度和水平位置的誤差則控制在2 km以內(圖8).而區域寬頻帶數據所對應的誤差水平要相對更大一些.同時我們還看到用區域寬頻帶速度波形得到的反演結果和強震數據得到的結果更一致，也更穩定.這樣的差別是合理的，因為首先區域數據所對應的濾波頻段要更低，而且臺站距離地震更遠，對地震的時空分辨率不可避免的要更低；其次在相同的濾波頻段速度記錄所對應的波形要比位移波形更復雜，因此對各子事件約束更為可靠.總之，區域寬頻帶數據的結果和強震數據結果非常一致，所對應的誤差范圍也足以分辨斷層幾何的變化和破裂過程的時空演化，驗證了可以將該方法擴展到區域臺站的反演中的有效性.

圖4 矩震級7.1級的主震破裂過程以及早期余震分布(a)強震和寬頻數據聯合反演所得的四點源解將主震分為四個由西南向東北傳播的子事件，其水平位置用較大圓圈表示，相應的震源機制解繪制成紅色沙灘球形狀.余震則由較小圓圈表示，并且震源深度用不同顏色顯示.黑色和紅色線段分別表示地震前所測繪的活動斷層和此次地震所產生的地表破裂 (Shirahama et al., 2016).近地震臺站用黑色三角形表示.(b)該四點源解(同(a))與測地學資料導出的有限斷層模型的對比.此次地震序列的主要事件由紅色五角星標記.右下角地震矩張量投影依次為四點源解的合成解、Global Centroid-Moment-Tensor (CMT)反演解和美國地質調查局的公布結果.(c)聯合高采樣GPS數據，強震地震數據和寬頻地震數據的反演結果.左上角為該四點源解的子事件震源時間函數(黑色)以及累加震源時間函數(灰色陰影).地震儀和GPS臺站分別由黑色三角形和白色方塊表示.(d)只利用寬頻地震儀反演的四點源解.其中黑色和灰色分別表示位移波形和速度波形的反演結果.左上角為位移波形反演得到的震源時間函數，而右下角為速度波形得到的震源時間函數.Fig.4 Rupture process of the MW7.1 mainshock and early aftershocks(a) The four-point-source solution of the joint inversion of the strong-motion and broadband data divides the mainshock into four subevents propagating from south-west to north-east. Four large circles denote their locations and four red beachballs indicate their focal mechanisms. Aftershocks are indicated by small circles that are color-coded by event depth. Black and red lines represent the previously mapped active faults and the surface ruptures (Shirahama et al., 2016) of the MW7.1 earthquake, respectively. Near-fault strong-motion stations are denoted by black triangles. (b) The comparison between the four-point-source solution (same as (a)) and the finite-fault model derived from the geodetic data. The red stars show the epicenters of the major earthquakes that occurred in this sequence. The lower-right inset shows the total moment tensor of the four-point-source solution in comparison with the Global Centroid-Moment-Tensor (CMT) and United States Geological Survey (USGS) solutions. (c) The four-point-source solution obtained from the inversion incorporating the high-rate GPS, strong-motion and broadband seismic data. Top-left inset shows the source time functions of subevents (black) and the total sum (gray-shaded area). The seismic and GPS stations are denoted by black triangles and white squares, respectively. (d) Four-point-source solutions obtained from the inversions with only F-net broadband data. Black and gray color indicates inversion results obtained from the displacement and velocity waveforms, respectively. The top-left shows the source time functions corresponding to the solution obtained from the displacement waveform inversion. The lower-right is corresponding to the solution obtained from the velocity waveform inversion.

圖5 強震和寬頻數據聯合反演所得的四點源解和三點源解的對比(a)四點源解的震源參數和地圖展示.空心圓圈表示子事件的水平位置，與其相對應的震源機制解相連.空心和實心的三角形分別表示朝向和遠離破裂方向的地震臺站.底部小圖顯示四個子事件震源時間函數(黑色)以及累加震源時間函數(灰色區域).(b)與(a)相似，但為三點源解的震源參數和地圖表示.底部小圖中綠色虛線為(a)中的累加震源時間函數.Fig.5 Comparison between the four-point-source and three-point-source solutions obtained from the joint inversions of the strong-motion and broadband data(a) Source parameters and the map view of the four-point-source solution. The subevent locations are denoted by open circles, connected with the beachballs that show the corresponding focal mechanisms. The open and solid triangles indicate the strong-motion stations that are toward and away from the rupture directivity. The bottom inset shows the source time functions of four subevents (black) and the sum of them (gray-shaded area). (b) Similar to (a), but for the source parameters and the map view of the three-point-source solution. The green dashed line in the bottom inset shows the cummulative source time function in (a).

圖6 路徑校正事件和主震三點源以及四點源反演的波形擬合效果比較此處只對比了利用強震和寬頻數據聯合反演的情況.四幅直方圖分別統計不同波形分量的互相關系數.Fig.6 Comparison of the quality of the waveform fits between the path calibration event, three-point-source inversion and four-point-source inversionHere the strong-motion and broadband waveforms are jointly used. The four panels of histograms show the statistics of the cross-correlation coefficients of different waveform components.

圖7 部分波形擬合效果(a)和(b)分別為圖5中四點源解和三點源解的波形擬合.帶通濾波后的觀測數據和合成數據分別顯示為黑色和紅色，其擬合所需時移以及表示擬合效果的互相關系數分別標記在波形下方.臺站名、臺站震中距、臺站方位角標記在波形左側，波形分量名標記在頂部.濾波頻段因臺站不同而有所區別，Pnl波的頻段總體在0.01～0.30 Hz，面波的頻段總體在0.01～0.22 Hz.Fig.7 Selected waveform fits(a) and (b) are waveform fits of the four-point-source and three-point-source shown in Fig.5. The bandpass filtered data and synthetic waveforms are shown in black and red, respectively. Numbers below each set of waveforms are the time shift used for alignment and cross-correlation coefficient of the two waveforms. Station names are shown on the left with the distance and azimuth indicated above and below respectively. The names of component are shown on the top. The filtering frequency is station-dependent, within 0.01～0.30 Hz for Pnl wave and 0.01～0.22 Hz for surface waves.

圖8 多組數據反演所得四點源解的不確定性分析(a)表示強震儀波形和寬頻地震儀波形數據的聯合反演不確定性.(b)表示GPS，強震儀波形和寬頻地震儀波形數據的聯合反演不確定性.(c)和(d)各代表只使用寬頻地震儀位移波形和只用寬頻地震以速度波形的反演的不確定性分析.以上四幅圖通過相似方式呈現MCMC采樣過程來反映其不確定性.四個子事件按照發震順序分別由紅、藍、紫、綠四色顯示.地圖上主要繪制了四個子事件水平位置采樣和反演得到的最佳震源機制.左上角繪制四個子事件斷層解抽樣.左側中間小圖為子事件震源時間函數的抽樣，其中粗黑線為反演最終震源時間函數，并且四點源累加震源時間函數由灰色填充區域顯示.右下角小圖為PP′豎直剖面上投影的子事件的深度抽樣，最佳深度顯示為黑色圓圈.Fig.8 Uncertainty analysis of the four-point-source inversions using multiple datasets(a) shows the uncertainty of the joint inversion of the strong-motion and broadband waveforms. (b) shows the uncertainty of the joint inversions incorporating the GPS, strong-motion and broadband data. (c) and (d) are for the inversions using only the broadband displacement and only the velocity waveforms, respectively. The uncertainty is represented by MCMC samples of the four subevents. The four subevents are shown in red, blue, purple and green in the order of occurrence. The sampled horizontal locations and the optimal focal mechanisms is shown in map view. The top-left inset shows the sampled fault-plane solutions of the subevents. The middle-left inset shows the sampled source time functions of the subevents with the final solution denoted by the thick black lines. The gray shaded area is the total source time function of the four-point-source solution. The bottom-right inset shows sampled depths projected to the PP′ cross-section with the optimal depths denoted by black circles.

此外，我們注意到當把高頻GPS數據包括在反演中時，得到的結果和強震數據得到的結果幾乎完全相同(圖9)，這說明數據之間的一致性非常好，也進一步驗證多點源反演的穩定性和可靠性.

圖9 GPS、強震儀和寬頻地震數據的聯合反演所得到的五點源解(a)五點源解的地圖表示.最靠近主斷層的GPS臺站由白色方塊表示.右下角為GPS臺0701附近區域的放大視圖.左上角為五個子事件對應的震源時間函數以及累積震源時間函數.(b)部分GPS臺站位移波形擬合，臺站位置如(a)所示.帶通濾波后的觀測數據和合成數據分別顯示為黑色和紅色，其擬合所需時移以及表示擬合效果的互相關系數分別標記在波形下方.臺站名、臺站震中距、臺站方位角標記在波形左側，波形分量名標記在頂部.濾波頻段因臺站不同而有所區別，總體在0.001～0.1 Hz.(c)0701觀測波形(黑色)與五點源解的合成波形(紅色)以及圖4a中四點源解的合成波形(淺藍)的對比.Fig.9 Five-point-source solution obtained from the joint inversion of GPS, strong-motion and broadband data(a) The map view of the five-point-source solution. The nearest GPS stations are denoted by the white squares. The lower-right inset is a enlarged view near station 0701. The top-left inset shows the source time functions of the five subevents and the cumulative source time function. (b) Selected waveform fits of GPS stations shown in (a). The bandpass filtered data and synthetic waveforms are shown in black and red, respectively. The numbers below each set of waveforms are the time shift used for alignment and cross-correlation coefficient of the two waveforms. Station names are shown on the left with the distance and azimuth indicated above and below respectively. The names of component are shown on the top. The filtering frequency is station dependent and is within 0.001～0.1 Hz. (c) The comparison between the observation (black) of 0701, five-point-source synthetics (red), and four-point-source synthetics (cyan) that was calculated with the solution in Fig.4a.

近震(強震+高頻GPS)和區域數據聯合反演的結果顯示需要至少4個子事件來表示主震.第一個子事件(MW6.75)位于12 km 深度，非常靠近其震中(大約在Hinagu斷層和Futagawa斷層的交界處)，具有一個高傾角(78°)的純走滑震源機制解，其走向(208°)和Hinagu斷層的走向非常一致.第二個(MW6.67)和第三個(MW6.62)子事件的深度幾乎相同(7 km)，分別位于震中東北側12 km和20 km處，它們的震源機制解(圖5a)也非常相似，值得指出的是該機制解的斷層走向(230°±3°)和Futagawa斷層的走向(232°)更為接近.從位置和斷層幾何來看，這兩個子事件應當發生在Futagawa斷層上，顯示少量的正斷層分量.最后一個子事件(MW6.54)最淺(2 km)，其震源機制解顯示斜向的滑動角，其走滑和正斷層分量幾乎相當.整個地震的震源持續時間大約20 s，最大的地震矩釋放發生在～8 s.我們得到的反演結果和Yue等(2017)以及Hao等(2017)的主震有限破裂模型具有很好的一致性，尤其是和Hao等(2017).我們的多點源解所揭示的斷層幾何形態，包括最北段的低傾角性質也與Yue等(2017)通過余震分布約束得到的幾何狀態吻合.我們將多點源解相加得到的等效矩張量解和其他機構給出的單個點源矩張量解也具有高度一致性(圖4b)，該一致性是對多點源解的有效驗證(例如Wei et al.,2011; Avouac et al.,2014；Shi et al.,2018).

4 討論

反演的結果顯示，隨著地震破裂的發展，子事件(1～4)質心深度不斷變淺并且越往后正斷滑移的分量越大，其水平位置處的余震的深度也隨之變淺.我們的子事件深度和Hao等(2017)和Yue等(2017)的有限破裂模型結果非常一致.Hao等(2017)將其解釋為斷層發震層深度(或者溫度)受附近的Aso火山影響所致，越靠近火山3～400 ℃等溫線的深度越淺，Lythgoe等(2021)對2018年Lombok 地震序列的解釋也類似.第二和第三個子事件所在的位置(圖4a)正好對應了余震的空區，余震則分布在更深和更淺處，這和其他地震(例如Wei et al., 2011)的觀測一致，可以很好地用應力變化來解釋；Yue等(2017)也指出了這一點.然而，第四個子事件附近幾乎沒有余震，而且在沿著斷層走向的方向約15 km處余震才重新出現.Yue等(2017)提出了地形應力加載對斷層上正應力和剪應力的影響，其模型計算發現，盡管設置了傾角為65°的正斷層，地形產生的正應力卻顯著升高反而不利于正斷層滑動.除溫度和地形應力的影響之外，不均勻分布的孔隙流體壓力的潛在影響也不能忽視.Aizawa等(2017)對于熊本地震區域電阻率研究結果顯示，Futagawa斷層淺部近Aso火山部分存在明顯高電導率，與高流體含量的特征相符.Matsumoto等(2018)利用Asano和Iwata(2016)的有限斷層模型，結合歷史地震活動，計算了Futagawa斷層上的理論孔隙流體壓力，其結果表明Aso火山附近的淺部孔隙流體壓力超過了該斷層的理論強度.這對斷層黏滑特性的改變是比較可觀的，同時也促進了主震破裂過程中該區域彈性勢能的充分釋放.值得注意的是，盡管本次地震序列中Aso火山口西北緣缺乏余震，但從2002年至2016年歷史地震活動記錄中，Uchide等(2016)發現這一區域實際上發生了大量地震.因此，這一區域的地震活動性有待未來進一步研究.

我們通過多點源反演所得到的四個子事件顯示出Futagawa斷層相比于Hinagu斷層具有更多的正斷層成分，尤其是其北段靠近Aso火山口附近的最后一個子事件擁有幾乎同等大小的正斷層和走滑分量.Yue等(2017)則指出該部分主震滑移包含較大比例的正斷層成分.實際上，在相對遠離Aso火山口的Futagawa斷層部分，第二個和第三個子事件由走滑主導.如此差異可能是由于火山附近地形增加了正斷層性質的相對剪切應力，而未同時增加走滑方向的剪切應力.精細的震后實地勘測(Shirahama et al., 2016)表明，Futagawa斷層地表破裂主要為走滑性質，而其東南側大約2 km處的Idenokuchi斷層的地表破裂為正斷層性質.結合野外測量和衛星遙感數據，Toda等(2016)提出了傾斜滑移分解模型，即Futagawa斷層主要釋放了長期彈性積累的走滑分量，而Idenokuchi斷層主要釋放了正斷分量，從而形成小規模地塹.而靠近火山口的Futagawa斷層可能缺乏這種走滑和正斷層的分解機制，從而傾滑過程(正斷+走滑)主導了同震位錯.長期的地殼尺度的地表形變觀測也可以為我們分析斷層性質成因提供參考.通過對2000—2010年分布在九州島的GNSS數據的反演，Mochizuki和Mitsui(2016)發現Aso火山附近尤其是其西緣的GPS臺站(32.87°N，131.00°E)的下移分量會被均一結構模型遠遠低估；而在Aso下方添加一塊與火山巖漿活動相關的席狀形變源則可以充分模擬出該臺站的下移.該結果表明，由于特殊的深部地質構造，Aso附近的Futagawa斷層可能在長期地殼形變中積累比遠離Aso區域更多的豎直方向的彈性應變，從而提高其地震活動釋放的正斷層分量.盡管從衛星干涉影像中并未獲得最后一個子事件附近的正斷層所產生的明顯地表形變，但我們認為最后一個子事件所揭示的淺部傾滑的特性是可靠的，并且可能并未在地表產生容易被勘測到的正斷層成分.

我們注意到位于斷層上或者離斷層很近的臺站(GPS臺0701和強震儀KMMH16)的波形擬合不是很理想(圖9c，淺藍色波形).我們認為這可能是多點源解對有限破裂過程近似所造成的影響.因為斷層上的臺站更容易受到最近破裂的影響，而這些局部破裂振幅雖然不比其他破裂更大，但由于非常靠近臺站，其近場項(隨距離的立方衰減)波場振幅在該臺站的記錄更大，因此沒有被反演.作為對比，我們參考了其他學者對于熊本地震震源過程的研究，發現大部分有限斷層模型都會明顯低估0701臺的總水平位移；而Zhang等(2018)的模型則是包含了靠近0701的近地表高滑移區，從而可以充分擬合0701的靜態位移，盡管其他GPS臺站位移量擬合不是很好.為了更好的擬合在斷層上的波形數據，必須要考慮到有限破裂過程的影響，使用有限破裂模型反演對這些近臺數據的擬合將會有進一步的提高.為了測試GPS臺站對于此次地震破裂究竟對破裂的矩心位置還是對最近的破裂更加敏感，我們增加了三組模擬測試(圖10).對于每組測試，先用六個點源構建一個輸入震源模型.該輸入震源模型包含兩個滑移集中區，其中一個位于西北部，靠近熊本主震的震中，由一個MW6.7的點源表示；另一個位于東北部，靠近GPS觀測站0701，由線性排列的五個MW6.3的點源組成(如圖10中紅色圓圈所示).我們用構建的輸入震源模型合成理論波形，并且加入相對強度為20%的高斯噪聲，用于反演測試.三組模擬測試分別針對了第二個滑移集中區的矩心位置相對臨近(圖10a)和遠離GPS臺0701(圖10b和10c)的情況，其中圖10b和10c分別對應該矩心位置位于0701的西側和東側的輸入震源模型.在接下來的反演測試中，我們只用兩個點源去恢復輸入震源模型，并且使用了和實際數據一樣的臺站選擇(包括近發震斷層的GPS臺站)和濾波頻段.三組反演結果(圖10)很好地恢復了輸入震源模型的矩心位置，(近似)震源時間函數以及震源機制解.因此，盡管兩個點源簡化了震源模型，加入接近斷層的波形數據依然能有效恢復輸入模型.仍應注意，在一些特殊情況下，即使使用了有限斷層模型，擬合非常靠近斷層的地震觀測也存在挑戰.比如在對斷層結構復雜性缺乏充分認知的情況，或是臺站所在位置的地質條件對波場傳播產生巨大影響的情況.例如，Lee等(2019)對2018年MW6.4中國臺灣花蓮地震的有限斷層反演中發現，盡管設置了三個斷層來模擬近震數據，但是該模型遠遠低估緊靠斷層的尤其是其中兩個斷層之間的臺站的速度振幅.他們將這解釋為該臺站所具有的較強的場地效應以及厚沉積層的影響.對于MW7.1熊本地震，我們試圖添加一個額外的子事件來模擬0701的波形.如圖9所示，通過聯合反演包含0701在內的GPS數據、強震儀波形和寬頻地震儀波形，獲得五點源解(圖9a).和四點源解類似，五點源解依舊顯示了自西北向東南方向的連續破裂過程，并且主要破裂區域的斷面幾何、震級大小和震源深度等都非常相似.但主要區別在于，五點源解將四點源解的最后一個子事件分為兩個位于Aso火山口旁的更小的正斷層為主的子事件(MW6.2-6.3)，并且這兩個子事件代表不同的斷層走向(圖9a及其右下角小圖).考慮這兩個子事件的位置誤差在內，它們反映出Futagawa斷層在Aso火山口附近形成了兩條分支斷層，其中一條保持主斷層的走向(大致232°，由第五個子事件表示)，另一條更加近東西走向(大致250°，由第四個子事件表示).此處的斷層分支很可能對地震破裂的停止起到了關鍵性作用.該斷層分支與Aso火山活動以及火山特殊地質結構的關系值得進一步探究，包括對長期歷史地震數據、火山監測數據、Aso區域地殼結構以及巖石學特征等可以提供較為全面的研究角度.上述五點源解可以非常好地擬合臺站0701的波形，并且合成圖和觀測波形的互相關系數達到100%(圖9c).值得注意的是，盡管先前所得到的四點源解無法理想地擬合0701的觀測波形，但是可以與五點源解幾乎同等程度地擬合其余GPS臺站(圖9b).總之，這印證了多點源反演方法的優勢：在波形數據指導下，通過提高點源數量來測試觀測數據對地震破裂過程的敏感性.比如在不使用緊靠Futagawa斷層的0701GPS臺站，四點源解可以充分擬合其余近震GPS、強震儀和區域寬頻地震儀的波形數據；而更加反映局部斷層破裂過程的五點源解則需要利用非常靠近該局部斷層的0701臺站的波形.這一經驗突出了緊靠潛在的發震斷層布置觀測臺站的重要性，這對普遍的地震破裂過程的研究非常關鍵(比如2018年MW6.4花蓮地震).

圖10 有限破裂效應的模擬測試(a) 輸入震源模型由六個點源組成(紅色)，用來合成模擬波形，用于后續兩點源反演.反演結果(藍色)有效地恢復了輸入震源模型中的兩個滑移集中區的矩心位置、震源機制解以及震源時間函數(見右下小圖).輸入震源模型和反演結果所對應的GPS臺站0701 的波形三分量分別由左上小圖中的紅線和藍線表示.濾波頻段為0.001～0.1 Hz. 圖(b)和(c)與(a)的繪制方法類似，對應另外兩組模擬測試.(b)中第二個滑移集中區相對于(a)朝西南方向平移了大約2 km.相反，(c)中第二個滑移集中區相對于(a)朝東北方向平移了大約2 km.Fig.10 Synthetic tests for the finite rupture effect(a) The input six-point-source model (red) are used to compute the synthetic waveforms for a two-point-source inversion. The inversion result (blue) well recovered the centroid locations, focal mechanisms and source-time functions (lower-right inset) of the two asperities of the input source model. The three-component waveforms computed for the input and inverted source models are shown in red and blue, respectively, for GPS station 0701 and filtered to 0.001～0.1 Hz. (b) and (c)Similar to (a), but for another two synthetic tests. The centroid locations of the second asperity in (b) is shifted to southwest by about 2 km relative to (a). The centroid locations of the second asperity in (c) is shifted to northeast by about 2 km relative to (a).

此外，0701臺的強振幅波形有可能包含了其他同震地質事件的貢獻，比如該地震序列引發的最大的山體滑坡(Hung et al., 2018; Song et al., 2019)剛好發生在0701西北方約1 km處.Hung等(2018)認為該滑坡附近的強震儀明顯記錄到了滑坡造成的位移波形，并由此確定滑坡對應的波形起始時間為主震發震時刻之后18.45～21 s.這一結果支持該最大滑坡是由主震動態觸發.但是是否可以從0701臺站的波形數據中分離出明顯的對應滑坡的成分，這需要對滑坡的運動學過程進行進一步模型設置，這超出了本文的探究范圍.我們用剪切位錯(雙力偶)源對0701波形的擬合已接近完美，并且得到的斷層走向和地表破裂極其吻合，我們認為滑坡對0701波形的貢獻應當比較有限.

熊本地震的寬頻帶臺站數據質量相當高，并且主要的傳播路徑是在大陸上，如果地質條件更為復雜，例如地震波傳播要經過大陸和海洋，尤其是經過海溝，或者是發生在海溝附近的地震，其對應的地震波記錄會更加復雜.在此情況下，一維格林函數就不再適用，如果仍舊使用一維格林函數則必須要使用更長周期來對數據進行濾波，那么對震源破裂過程的分辨率就會下降.而Qian等(2019)和Wu等(2018)所提出的方法和應用表明，三維格林函數對源區結構復雜的波場傳播在10～20 s的周期是不可忽略的，這些研究也表明在沒有很厚沉積層的情況下考慮地形和水層就可以對波形擬合大為改進.因此使用三維格林函數進行多點源反演是將來一項重要的工作.我們注意到Duputel和Rivera(2017)在地震多點源反演中使用了全球SEM(Spectral Element Method，譜元法)程序(Komatitsch和Tromp, 2002)計算的三維格林函數，但是全球版SEM程序對海水使用了應力近似邊界條件，并不能模擬水波及其多次波.可能也是由于這個原因，Duputel和Rivera(2017)盡管使用了三維格林函數，依然用了50 s以及更長周期的地震波來對2016 年發生的MW7.8 Kaikouro地震進行反演.

我們所使用的正則化的Yoffe函數在多點源地震反演更符合震源破裂動力學(Tinti et al., 2005).在地震破裂過程的實際研究中，由于震源時間函數和其他震源參數之間存在取舍關系，并且高頻地震波的快速衰減不容忽視，因此破裂過程的震源時間函數往往極難分辨.顯然，在波形擬合可靠(經過校正)的情況下，反演所用的頻段越高就越能分辨破裂過程的震源時間函數.我們的結果顯示前三個子事件的震源時間函數更接近于高斯函數，在速度和位移波形反演中均是如此.最后一個子事件在速度波形反演結果中更接近于Yoffe函數.對比于離該子事件最近的幾個臺站波形(圖7)，我們認為這是一個可靠的結果.這說明淺部的地震破裂很可能受到前面子事件造成的動態應力傳播的影響.在反演所用的頻段，我們不需要具有強不對稱性的Yoffe函數來擬合數據，這有兩種解釋：(1)熊本地震的主震的震源時間函數確實更接近于高斯函數；(2)需要用更高頻段的地震波來進一步分析.我們提出的正則化的Yoffe函數震源時間函數可以在其他地震的破裂過程研究中使用，其靈活性將有助于我們進一步了解破裂形成和發展的物理機理.

5 結論

我們將一個基于MCMC的地震多點源波形反演方法擴展到了區域距離的地震波形數據，并對2016年MW7.1熊本地震進行了應用.結果顯示區域地震數據可以較好地分辨此類中強度(M6-7)地震的時空破裂過程.利用離斷層很近的0701高頻GPS數據可以分辨出地震破裂最后階段更為精細的破裂過程.我們可以將該方法進一步推廣到全球中強度地震的運動學過程研究中，來更準確地確定基本的震源物理參數，以理解震源物理過程.

致謝謹此祝賀陳颙先生從事地球物理教學科研工作60周年.感謝兩位匿名審稿人寶貴的建設性意見，感謝日本Kik-net，K-net，F-net為本文研究所提供的波形數據.