山西原平ML4.7地震矩張量反演

李 斌，李自紅，宋美琴，關鵬虎

(1.太原理工大學 a.礦業(yè)工程學院，b.地震與地質(zhì)災害防治研究所，太原 030024；2.山西省地震局，太原 030021)

北京時間2016年4月7日04時49分，在山西省原平市蘇龍口鎮(zhèn)(38.87°N，112.91°E)發(fā)生了ML4.7地震，震源深度11 km，山西原平、代縣、定襄、太原、晉中、陽泉以及朔州等地震感明顯。這是自2016年1月1日以來，繼2016年3月12日運城M4.4地震、3月15日運城M3.0地震和3月27日運城M4.0地震后在山西省境內(nèi)發(fā)生的第4次顯著地震事件，也是2009年3月28日原平ML4.5地震打破山西地震帶近5年的M4地震平靜后在該地區(qū)發(fā)生的又一次顯著地震。以往震例研究表明，山西地震帶中強地震“南北遷移”活動特征明顯[1-2]。此次原平ML4.7地震的發(fā)生，是山西地震帶顯著地震活動由南向北遷移的表現(xiàn)，對山西地震帶中強地震具有重要的預警作用，預示著山西地震帶可能進入了新一輪的4級地震活躍時段[3]。因此，深入理解此次顯著地震事件的震源機制，有助于理解地震發(fā)生時震源區(qū)的應力狀態(tài)與發(fā)震構造，以及山西地震帶中強地震的孕育規(guī)律與活動特征，也是跟蹤震情發(fā)展趨勢的重要依據(jù)。

目前顯著地震震源機制的求解方法較多，如P波初動法[4]、振幅比法[5]、P波初動聯(lián)合振幅比法[4,6]、CAP全波形反演法[7-8]與矩張量反演法[9-11]等，且各種方法的研究結果都對世界上很多地震作出過較為合理的解釋[12-16]。地震矩張量相比于其他方法，對于顯著地震的震源機制解能有較為準確的描述[11]，且具有一定的優(yōu)勢。相對于P波初動與振幅比法，矩張量反演對臺站布局要求低，且只需數(shù)字波形資料，無需進行震相分析，避免了人工量取震相誤差帶來的不確定性；反演結果反映的是整個破裂過程的信息，而非僅僅是初始破裂信息[17]。相比于CAP方法，矩張量反演方法的計算結果中含補償線性矢量偶極 (CLVD)和各向同性分量(ISO)，因此地震矩張量可以分解為表示斷層面剪切錯動的雙力偶分量(DC)、震源體膨脹或收縮的各向同性分量(ISO)、震源體優(yōu)勢方向的張裂或擠壓變形(CLVD)等[9-10]。通常構造地震是由斷層兩盤的剪切位錯觸發(fā)，因此矩張量能更好地反映斷層破裂信息[18]。

本文利用山西區(qū)域地震臺網(wǎng)寬頻帶數(shù)字波形資料，在采用Hypo2000對地震重新定位的基礎上，采用TDMT_INV時間域矩張量反演方法對2016年4月7日山西原平ML4.7地震進行矩張量反演計算，通過逐個增減參與反演的臺站數(shù)量進行多次反演，結合震源機制隨反演深度的變化特征進行參數(shù)的穩(wěn)定性評價。將獲得的震源機制結果與不同臺站垂向分量的P波極性在震源球上的分布位置，以及與利用Snoke方法和CAP方法得到的震源機制解進行對比，以驗證參數(shù)的合理性與可靠性。研究結果將有助于理解山西地震帶中強地震的孕育規(guī)律與活動特征。

1 方法與原理

國內(nèi)外學者對地震矩張量的求解方法進行了大量的研究。PATTON et al[19]首先將遠震矩張量反演算法做了修改，進行了區(qū)域地震矩張量反演。DREGER et al[9]提出了在時間域利用區(qū)域地震波形(Pnl)進行地震矩張量反演的方法(time domain moment tensor inversion，TDMT_INV).本文采用的是DREGER于2003年改進理論波形計算方法后的TDMT_INV方法[10]，即在計算理論格林函數(shù)時加入了SAIKIA[20]改進的離散波數(shù)積分法。目前該方法在國際上已經(jīng)得到廣泛應用[21-24]。隨著寬頻帶數(shù)字臺網(wǎng)的建立，國內(nèi)許多研究人員也利用該方法反演震源機制解：趙翠萍等[25]、唐蘭蘭等[26]、屠泓為等[27]將此方法應用于伽師震源區(qū)；王勤彩等[28]反演了汶川大地震余震序列88個地震的矩張量解。

TDMT_INV方法基于雙力偶點源模型，在自由表面距離震源x處t時刻的實際觀測記錄dn(x,t)可表示為[29-31]：

(1)

(2)

式中：oi和si分別為觀測數(shù)據(jù)與基于格林函數(shù)計算的時間域理論數(shù)據(jù)。VR的計算不僅考慮了波形的相似性，還考慮了絕對振幅的大小；VR值越大，表明波形擬合程度越高，求解出來的地震矩張量越可靠。在不同深度上搜索VR的最大值，以獲得最佳震源深度與震源機制的最優(yōu)解。

2 資料選取與數(shù)據(jù)處理

2.1 地震波形數(shù)據(jù)

目前，山西區(qū)域測震臺網(wǎng)由72個實時記錄與傳輸?shù)臄?shù)字化寬頻帶地震臺站組網(wǎng)而成(山西省內(nèi)臺站57個，鄰省臺站15個)，基本上均勻分布于山西省及周邊地區(qū)(圖1(a)).現(xiàn)有臺網(wǎng)條件下，ML≈2.0級及以上的地震，通常能同時被周圍幾十個臺站記錄到[32]。對于此次原平ML4.7地震，山西區(qū)域測震臺網(wǎng)72個地震臺站全部記錄到此次地震，38個臺站的垂向分量有清晰的P波初動。本文共選取了震中距120 km范圍內(nèi)8個臺站的波形數(shù)據(jù)參與反演，臺站方位及分布如圖1(b)所示。選取數(shù)據(jù)的基本標準為：記錄波形質(zhì)量好，信噪比高，二階Butterworth濾波器濾波后波形清晰，滿足矩張量反演數(shù)據(jù)可靠性的要求。

HAVSKOV et al[33]研究表明，對于質(zhì)量較好的波形數(shù)據(jù)，利用速度波形數(shù)據(jù)或積分至位移數(shù)據(jù)進行矩張量反演,其結果具有較好的一致性。本文在數(shù)據(jù)預處理過程中對比發(fā)現(xiàn)，速度波形數(shù)據(jù)在三個分量旋轉(zhuǎn)及濾波后波形表現(xiàn)得更加穩(wěn)定。鑒于此，本文選用了速度波形數(shù)據(jù)進行后續(xù)的反演計算。數(shù)據(jù)處理與計算過程主要包括：

1) 對原始波形數(shù)據(jù)進行去均值、去傾斜、去儀器響應的預處理。

2) 將3個分量分別旋轉(zhuǎn)至垂向、切向和徑向分量。

3) 選取二階Butterworth濾波器進行濾波。頻帶范圍的選取與震級有關；對于本次ML4.7地震，選取的頻帶范圍為0.02～0.05 Hz.

4) 設定計算起始深度1.1 km(避免計算深度位于速度結構分層界面上)及步長1.0 km，基于區(qū)域速度結構等參數(shù)計算各臺站不同深度范圍的理論格林函數(shù)。

5) 進行矩張量反演及反演結果置信度評估。

圖1 山西地震臺網(wǎng)臺站(a)及參與計算的8個地震臺站(b)分布圖Fig.1 Distribution of stations of Shanxi Seismic Network (a) and the 8 selected stations(b)

2.2 地殼速度結構

以往研究表明，無論哪種震源機制求解方法，地殼速度結構都是影響震源機制解的重要因素之一[17,24]。為了減少地殼速度結構對矩張量反演結果的影響，本文在Hpyo2000重新定位與矩張量反演過程中均采用了山西地震帶最新的1D速度結構模型[34]，地殼分層及數(shù)據(jù)詳見表1.

表1 本文計算采用的山西地震帶1D速度結構模型Table 1 1D velocity model of Shanxi seismic belt applied in this study

3 計算結果與分析

3.1 重新定位結果

矩張量反演需要首先計算各參與反演臺站的理論格林函數(shù)。除地殼速度結構外，震源位置特別是震源深度誤差會對理論格林函數(shù)的計算結果產(chǎn)生影響，進而可能影響矩張量反演的最終結果，因此準確的定位結果是可靠矩張量反演的基礎數(shù)據(jù)之一[6,17]。對于本次原平ML4.7地震，山西地震臺網(wǎng)與中國地震臺網(wǎng)均給出了定位結果，見表2.兩個結果的震中位置基本一致，但山西臺網(wǎng)給出的震源深度為11 km，中國臺網(wǎng)為16 km，偏差較大。鑒于此，本文基于山西地震帶最新的1D地殼速度結構，采用Hpyo2000[35]定位方法對本次地震進行了重新定位，得到主震震中位置為：38.86°N，112.92°E，震源深度14.7 km. Hpyo2000定位方法由Klein開發(fā)，采用傳統(tǒng)的Geiger法基本思路，更適合于觀測臺站分布均勻的網(wǎng)內(nèi)近震及地方震。重新定位結果將為后續(xù)矩張量反演提供保障。

表2 2016年4月7日山西原平ML4.7地震定位結果Table 2 Location of the Yuanping ML4.7 earthquake on April 7, 2016

3.2 震源機制結果

本文共選取了震中距120 km范圍內(nèi)8個地震臺站記錄的三分量速度波形數(shù)據(jù)進行地震矩張量反演。采用的濾波頻帶為0.02～0.05 Hz，起算深度1.1 km、步長1.0 km，共反演了1.1～25.1 km范圍內(nèi)不同深度的矩張量結果。根據(jù)反映觀測數(shù)據(jù)與理論波形擬合情況的約化方差VR和Var/Pdc值隨深度變化的情況，搜索出擬合程度最好、雙力偶成分相對較高的解，即為最佳震源機制解。圖2為矩張量反演過程中觀測數(shù)據(jù)與理論波形最佳擬合情況下的反演結果。參與反演的8個地震臺站中，除SHC臺垂直分量和SHY臺南北分量因波形有一定畸變在反演過程中被刪除外，其余各分量的觀測波形與理論波形的擬合程度尚好，平均VR=56.3%，結果符合信度要求[6,24,33]。

圖2 山西原平ML4.7地震震源機制的波形擬合情況及矩張量反演結果Fig.2 Waveform fitting and results of moment tensor inversion of Yuanping ML4.7 earthquake

反演得到的最佳雙力偶解節(jié)面I參數(shù)為走向143°，傾向63°，滑動角-66°；節(jié)面II走向278°，傾向36°，滑動角-128°；震源機制類型屬于正斷兼有右旋走滑分量。雙力偶分量Percent DC=93，線性補償分量Percent CLVD=7，震源機制符合雙力偶點源的假設條件。主壓應力P軸方位角為94°，傾角為64°；主張應力T軸方位角為216°，傾角為15°.標量地震矩M0= 2.492 64×1015N·m，矩震級結果為MW=4.2，與MS震級基本一致[3]。

3.3 震源深度的確定

圖3是約化方差VR與震源機制解隨不同反演深度的變化圖。可以看出，隨著反演深度的不斷增加，震源機制結果逐漸從正斷向正斷兼有走滑分量、再向逆斷趨勢變化；但在5～16 km深度范圍內(nèi)震源機制性質(zhì)及結果基本穩(wěn)定，表現(xiàn)為正斷兼有右旋走滑分量，且在11～14 km范圍內(nèi)VR相對穩(wěn)定且取得最大值，為擬合程度最好、雙力偶成分相對較高的解。因此，判斷本次地震的最佳震源深度在11～14 km；這與重新定位結果的震源深度14.7 km具有一定可比性。

圖3 約化方差VR隨反演深度的變化Fig.3 Variance of variance reduction versus inversion depth

3.4 參數(shù)穩(wěn)定性分析

影響地震矩張量反演結果穩(wěn)定性的因素較多，如區(qū)域地殼速度結構、地震定位誤差(特別是震源深度誤差)、參與反演臺站的個數(shù)、臺站與震中的方位角誤差，等等[17]。為進一步分析與驗證反演參數(shù)的穩(wěn)定性，在反演過程中我們首先對參與反演的8個臺站，按每次增減一個臺站進行反演；同時，基于Hypo2000重新定位結果，利用本區(qū)域最新的地殼速度模型進行了多次反演。具體反演結果詳見圖4.圖4顯示，在5、6、7、8個臺站參與反演的情況下，最佳震源機制解節(jié)面I的走向最大偏差15°，傾角偏差8°，滑動角偏差10°.在參與反演臺站數(shù)量減少的情況下，約化方差VR有所增高，但其余反演參數(shù)并未出現(xiàn)太大的變化。結合圖3中5～16 km深度范圍內(nèi)震源機制性質(zhì)基本穩(wěn)定，及11～14 km范圍內(nèi)VR值穩(wěn)定且取得最大值的結果，認為本文基于8個地震臺站反演得到的參數(shù)基本穩(wěn)定。

圖4 不同數(shù)量臺站參與矩張量反演時最佳震源機制解及其他參數(shù)的取值Fig.4 Focal mechanisms and other parameters from moment tensor inversion based on different stations

3.5 可靠性與誤差分析

矩張量反演結果的合理性與可靠性通常可利用各記錄臺站垂直分量上的P波初動極性在震源球上的分布位置來分析判定。一般的做法是：首先利用區(qū)域地殼模型、震中距及震源深度計算記錄臺站的離源角；然后根據(jù)臺站方位與離源角將其投影到震源球上，位于壓縮區(qū)的P波初動向上，反之向下；通過考察P波極性與反演得到的震源機制解的一致性來驗證結果的合理性與可靠性[18]。

Strike=113°,Dip=53°,Rake=-74°,Source=Snoke圖5 原平ML4.7地震的矩張量反演結果(a)以及基于P波初動的Snoke方法求解的結果(b)Fig.5 Focal mechanism of Yuanping ML4.7 earthquake by moment tensor inversion (a) and the focal mechanism result by Snoke method (b)

對于本次原平ML4.7地震，山西測震臺網(wǎng)中共有38個地震臺站垂直分量記錄到清晰的P波初動，各個臺站的P波初動極性及在震源球上的分布位置情況如圖5(a)所示。除了少數(shù)幾個臺站外，其余32個臺站的P波極性與基于8個臺站波形數(shù)據(jù)反演得到的最佳震源機制解基本一致。

此外，為進一步分析結果的可靠性與可能誤差，基于38個地震臺站垂直分量記錄的清晰的P波初動，利用Snoke方法求解了本次地震的震源機制斷面解(圖5(b))，并將矩張量反演得到的震源機制解與Snoke方法結果，及文獻[3]利用CAP方法得到的解進行比較，具體震源機制解參數(shù)詳見表3.可以看到：3種方法得到的震源機制解存在一定的差異，但震源機制性質(zhì)一致，節(jié)面參數(shù)基本趨勢一致，總體上具有可比性；存在差異的主要原因在于，3種方法的理論原理不同，邏輯算法及利用的數(shù)據(jù)類型不同，等等[33]。

表3 原平ML4.7地震震源機制參數(shù)Table 3 Parameters of Yuanping ML4.7 earthquake focal mechanism solution

4 討論與結論

本文利用山西地震臺網(wǎng)寬頻帶數(shù)字波形資料，在Hypo2000重新定位的基礎上，采用TDMT_INV時間域矩張量反演方法，對2016年4月7日山西原平ML4.7地震進行了矩張量反演計算與參數(shù)的穩(wěn)定性評價，并將獲得的最佳雙力偶解與利用Snoke方法及CAP方法得到的結果進行對比，驗證了結果的合理性與可靠性。初步得到以下認識和結論：

1) 基于震中距120 km范圍內(nèi)的8個觀測臺站的波形資料反演得到的地震矩張量結果顯示：最佳雙力偶分量Percent DC=93，線性補償分量Percent CLVD=7，震源機制符合雙力偶點源的假設條件；節(jié)面Ⅰ參數(shù)為走向143°，傾向63°，滑動角-66°；節(jié)面Ⅱ走向278°，傾向36°，滑動角-128°；最佳震源深度11～14 km.主壓應力P軸方位角94°，傾角64°；主張應力T軸方位角216°，傾角15°；標量地震矩M0=2.492 64×1015N·m，矩震級結果為MW=4.2，與MS震級基本一致。

2) 影響地震矩張量反演結果穩(wěn)定性的因素較多，如區(qū)域地殼速度結構、地震定位誤差(特別是震源深度誤差)、參與反演臺站的個數(shù)、臺站與震中的方位角誤差，等等。本文基于研究區(qū)最新的地殼速度結構，在Hypo2000重新定位的基礎上，按逐次增減一個臺站的方式進行多次反演，分析與驗證了結果的穩(wěn)定性，認為本文基于8個地震臺站反演得到的參數(shù)基本穩(wěn)定。此外，各臺站垂向分量上清晰的P波極性在震源球上的分布與反演得到的最佳震源機制解基本一致。

3) 本文矩張量反演得到的最佳震源機制解與基于Snoke方法和CAP方法的結果存在一定的差異，但總體具有可比性。結果間存在誤差的原因在于，不同方法的理論基礎不同、利用的數(shù)據(jù)類型不同等。基于P波初動的斷層面解反映的是初始破裂情況，而矩張量解能更好地詮釋整個地震破裂過程。

致謝感謝山西地震臺網(wǎng)中心為本研究提供地震波形數(shù)據(jù)。