稀疏臺網(wǎng)的中小地震震源機制解反演研究——以2017年1月4日西藏仲巴4.7級地震為例

梁姍姍 徐志國 劉杰 黃志斌 鄒立曄 陳楠

1)中國地震臺網(wǎng)中心，北京 100045 2)中國科學院大學計算地球動力學重點實驗室，北京 100049 3)國家海洋環(huán)境預報中心，北京 100081

0 引言

測定中小地震的震源機制解是地震學應用和研究領(lǐng)域中一項重要的工作。由于中小地震(M<5.0)發(fā)生頻率較高，其震源機制為區(qū)域構(gòu)造和構(gòu)造應力場等研究提供重要信息。震源機制解反映了震源區(qū)附近的應力狀態(tài)，對大地震前后發(fā)生的中小地震震源機制的研究，有助于認識大地震的發(fā)震構(gòu)造特征，探索大地震的孕育過程，并為研究大地震的動態(tài)破裂模型及強余震觸發(fā)模型提供重要的約束條件(王勤彩等，2009)。利用中小地震矩張量反演得到的地震矩、震源機制解、震源深度等參數(shù)可以解釋斷層在地震前后具體的運動情況，并為地震災害評估、震后地震應急服務(郭祥云等，2014)。

中小地震震源機制解求解方法主要有P波初動法、振幅比法以及矩張量反演法(俞春泉等，2009；崔效鋒等，2011)。對于中小地震來說，由于震級小、能量弱，通常難以利用遠場波形反演地震矩張量解(趙翠萍等，2008；唐蘭蘭等，2012)，在觀測資料不足的情況下，近場波形反演結(jié)果也存在著較大的不確定性。而對于稀疏臺網(wǎng)區(qū)域，只有少量臺站記錄到清晰資料，僅僅依靠P波初動很難得到可靠的震源機制解，其在很大程度上依賴于P波初動資料的數(shù)量及其分布狀況(胡幸平等，2008)。因此，反演稀疏臺網(wǎng)記錄的中小地震的震源機制解，是一項非常具有挑戰(zhàn)性的工作。

本文以2017年1月4日西藏區(qū)域臺網(wǎng)記錄的仲巴縣M4.7地震為例，首先利用西藏區(qū)域臺網(wǎng)近震波形資料，采用ISOLA(ISO lated asperities)方法對該地震進行震源機制反演，并探討了結(jié)果的穩(wěn)定性和可靠性，確定其最佳震源機制解；然后利用清晰的P波初動資料，采用CSPS(cyclic scanning of the polarity solution)初動掃描法，利用P波初動極性資料和少量臺站近震波形約束，獲取最佳震源機制結(jié)果，并與ISOLA近震全波形反演結(jié)果相比較，驗證聯(lián)合約束反演方法的可行性。本文的工作將有助于中小地震震源機制的研究，尤其是稀疏臺網(wǎng)記錄的網(wǎng)緣或網(wǎng)外地震，為深入了解稀疏臺網(wǎng)地區(qū)的區(qū)域構(gòu)造和構(gòu)造應力場等研究工作提供有力的數(shù)據(jù)支撐。

圖1 西藏仲巴地震的震中位置和本研究使用的地震臺站分布五角星表示2017年1月4日西藏仲巴M4.7地震；紅色三角形表示ISOLA近震全波形反演所用臺站；灰色三角形和圓圈內(nèi)紅色三角形分別表示利用P波初動和近震波形聯(lián)合約束反演所用臺站；震源機制解為本文近震全波形反演結(jié)果

1 數(shù)據(jù)和方法

1.1 資料選取與處理

本文使用了西藏自治區(qū)數(shù)字地震臺網(wǎng)提供的臺網(wǎng)觀測報告給出的定位結(jié)果和區(qū)域地震波形資料，選取了12個清晰的地震臺站的波形記錄參與震源機制計算，圖1 給出了本研究使用的臺站分布，震中距為107～760km，在500km范圍內(nèi)僅有5個臺站，其余臺站主要集中分布在震中以東地區(qū)，臺站分布稀疏且方位角覆蓋較差。

計算過程中，首先由波形觀測資料讀取P波初動極性，選取初動清晰的臺站資料用于波形約束反演；然后將波形記錄去除儀器響應、去傾斜以及去平均值等預處理，應用帶通濾波器檢查波形低頻干擾成分(Zahradnik，2005；Zahradnik et al，2010)，去除信號異常的臺站記錄，主要是因為低頻干擾成分會對反演結(jié)果造成影響，降低反演結(jié)果的準確性。在實際反演過程中，為了突出有效信號，減小地殼精細結(jié)構(gòu)對反演結(jié)果造成的影響，我們根據(jù)震中距范圍選取合適的濾波頻帶。

當采用ISOLA雙力偶約束反演模式反演震源機制時，反演過程采用點源模型。選擇震中距在500km范圍內(nèi)的5個臺站的寬頻帶波形記錄(最大震中距約為419km)，對反演波形進行帶通濾波，臺站分布如圖1 中紅色三角形所示。為了減少地殼速度模型不確定性對研究結(jié)果的影響，濾波頻帶的下限頻率應盡可能低，其主要受臺站背景噪聲影響，上限頻率主要受震中距、地震大小和地殼速度模型的準確程度影響。一般來說，我們能夠反演的波形記錄的震中距不大于10倍MSW(MSW為最小剪切波長度)(Fojtíkovet al，2014)。假設(shè)剪切波速度vs大約為 3.5km/s，反演的最高頻率為0.08Hz，其對應的MSW為12.5 s×3.5km/s=43.75km，所以，參加反演臺站的最大震中距應小于10MSW(即437.5km)。本研究中通過對每個臺站進行波譜分析，確定了每個臺站的濾波頻帶范圍為0.05～0.08Hz，采用二階Butterworth濾波器進行濾波。我們采用圖2 中的速度模型和離散波數(shù)法(Bouchon，1981)計算格林函數(shù)，采樣間隔為0.2Hz。

圖2 地殼速度模型

當應用CSPS初動掃描法，利用P波初動資料和近震波形聯(lián)合約束反演此次地震的震源機制時，我們選用了圖1 中12個臺站的P波初動極性資料，并聯(lián)合距離震中最近的2個臺站近震全波形臺站資料(圖1 中圓圈內(nèi)的2個紅色三角形)。在震源機制反演中，若選取的速度模型嚴重偏離真實地殼模型，其反演的結(jié)果具有較低的可信度，本文使用的一維速度模型參考了Crust2.0全球速度模型(Bassin et al，2000)。

1.2 研究方法

本研究主要采用Sokos等(2008、2013)的ISOLA(ISO lated asperities)全波形矩張量反演方法，該方法采用單一點源或多點源模型，通過網(wǎng)格搜索和最小二乘法反演得到地震的最佳矩張量解。該方法具有純雙力偶(DC)、偏量矩張量(雙力偶分量DC+補償線性矢量偶極分量CLVD)和全矩張量(DC+CLVD+ISO分量)和固定震源機制4種求解方式，可以在點、線、面3種不同維度空間范圍搜索最佳解，其中固定震源機制解求解方式僅反演得到矩心時間、深度和地震矩。同時可以根據(jù)臺站分布選取單一或多重的地殼速度模型，降低震源-臺站復雜射線路徑(李飛等，2017)對反演結(jié)果的影響；反演過程中，該方法給出了多種量化參數(shù)來評估反演解的穩(wěn)定性和可靠性，如VR、CN等參數(shù)。VR(Variance Reduction)為方差減少量，定義為

VR=1-∑(o-s)2/∑o2

(1)

式中，o和s分別表示觀測波形和理論波形；VR為表征觀測波形和理論波形的擬合度的重要參數(shù)。CN(Condition Number)為條件數(shù)，其與參與計算的臺站分布、濾波頻段和速度模型有關(guān)，而與實際參與反演波形無關(guān)。反演過程中，需要對矩陣方程進行求解計算，即

u=Gm

(2)

為了反演稀疏臺網(wǎng)記錄的中小地震震源機制解，特別是在僅有少量近震臺站波形記錄，但較大震中距范圍內(nèi)臺站記錄P波明顯可用的情況下，F(xiàn)ojtíkov等(2014)等提出了CSPS方法，該方法利用P波初動極性資料和少量臺站近震波形約束求解震源機制解的方法，其思路為：采用P波初動資料在空間范圍內(nèi)搜索震源機制結(jié)果，反演結(jié)果作為輸入初始結(jié)果，聯(lián)合近震波形數(shù)據(jù)，采用固定震源機制求解方式，計算波形擬合方差，選擇波形擬合程度最高對應的震源機制結(jié)果為最優(yōu)解。該方法最大的優(yōu)勢就是結(jié)合了P波初動極性解信息，約束了由于少量臺站參與波形反演引起解的不穩(wěn)定性情形。

為了驗證P波初動和少量近震全波形約束反演中小地震震源機制解的有效性和可用性，我們采用ISOLA中的初動極性循環(huán)掃描法CSPS(Fojtíkovet al，2014)對西藏仲巴M4.7地震求解震源機制解，采用ISOLA中純雙力偶約束的固定震源機制解求解方式，應用圖1中12個臺站的P波初動極性資料，聯(lián)合少量的近震全波形臺站資料(圖1中圓圈內(nèi)的紅色三角形)，計算波形擬合程度VR值，應用VR值量化固定震源機制反演結(jié)果。本文中應用FOCMEC(Snoke，2003)計算得到的多組走向角、傾角和滑動角值，作為CSPS初動掃描法的輸入值。FOCMEC方法利用P波、SV波和SH波的初動方向以及振幅比聯(lián)合測定震源機制解。通過比較理論計算和實際觀測所得的P波、SV波、SH波的初動方向和振幅比矛盾數(shù)最小的方式得到震源機制解。在本研究中未采用S波初動和振幅比進行約束，僅采用P波初動極性資料計算初始震源機制解，通過近場全波形對輸入初始解進行約束。

2 結(jié)果分析

2.1 近震全波形反演

在反演過程中，首先參考西藏區(qū)域臺網(wǎng)觀測報告內(nèi)容給出的初始地震定位參數(shù)(30.590°N，83.784°E；深度7km)。固定震中位置，在深度方向進行搜索，搜索范圍為0～20km，搜索步長為2km；矩心時間搜索范圍設(shè)為發(fā)震時刻前后4s，搜索步長為0.1s，反演不同深度的震源機制解，以波形最大擬合方差相應的震源深度為最佳矩心深度；其后，固定最佳矩心深度，以初始震中和最佳矩心深度為中心點，搜索范圍為30km×30km，網(wǎng)格步長為5km×5km，共49個網(wǎng)格點，在平面網(wǎng)格內(nèi)搜索最佳震源機制解。

由于參加反演的臺站數(shù)量、臺站分布方位等因素影響計算震源機制結(jié)果的可靠性，為了得到最佳的反演結(jié)果，我們根據(jù)震中距范圍由近及遠分別采用不同的臺站數(shù)量來反演此次地震的震源機制，反演得到矩心深度、矩心時間、斷層節(jié)面解參數(shù)、矩震級MW以及相關(guān)解的質(zhì)量評估參數(shù)等，詳細信息如表1所示。圖3為由不同臺站反演得到的仲巴地震的震源機制，由圖3可以看出，臺站數(shù)量少于2個的反演結(jié)果(圖3(a)、3(b))，與哈佛大學GCMT結(jié)果差別很大，震源機制解性質(zhì)完全不同；臺站數(shù)量大于2的反演結(jié)果(圖3(c)、3(d)、3(e))與哈佛大學結(jié)果相差不多。由此表明，少量臺站波形記錄反演的震源機制解結(jié)果可靠性低，并不能真實地反映地震的斷層面解。

表1不同臺站反演結(jié)果及質(zhì)量評估參數(shù)

結(jié)果來源數(shù)據(jù)類型|臺站數(shù)量|臺站名稱節(jié)面Ⅰ節(jié)面Ⅱ矩心深度/km矩震級VRCN走向/(°)傾角/(°)滑動角/(°)走向/(°)傾角/(°)滑動角/(°)GCMT體波|14|—面波|85|—10876163202731518.74.8——ISOLA全波形 |1|ZBA2012119161808444.60.9025.4全波形|2|ZBA,GZE9248154200714544.60.876.5全波形|3|ZBA,GZE,PLA1571212777016044.50.742.8全波形|4|ZBA,GZE,PLA,NMA10568175197852244.50.672.7全波形|5|ZBA,GZE,PLA,NMA,SQHE10985-1771988-464.60.713.0

由表1可知，少量臺站參與矩張量反演結(jié)果對應的臺站觀測波形和理論波形擬合程度非常好，方差VR>0.8，但條件數(shù)CN>6較大，反演矩陣的穩(wěn)定性差，其對應的結(jié)果并不能表示真實的地震斷層面參數(shù)。這也表明，基于波形的矩張量反演方法并不能僅僅通過波形擬合程度來評估反演結(jié)果的可靠性，要盡可能多的臺站參與反演計算，并聯(lián)合臺站初動方向來檢驗反演結(jié)果的可靠性。

圖3 由不同臺站反演得到的仲巴地震的震源機制圖(a)、(b)、(c)給出了不同臺站數(shù)量(1～5個)時ISOLA反演震源機制結(jié)果和哈佛大學反演結(jié)果；黑色線為ISOLA反演結(jié)果對應的節(jié)線；黑色圓為P 軸；白色圓為T 軸；陰影區(qū)對應的節(jié)面為哈佛大學反演結(jié)果

圖 4 不同深度上仲巴地震的震源機制解及其反演方差

圖5 固定震中位置垂直深度方向搜索西藏仲巴地震震源機制解對應的三分量理論(紅色)和實際(黑色)波形擬合圖波形下方數(shù)字表示波形擬合方差量；左側(cè)大寫字母表示臺站名

我們利用5個臺站波形反演不同深度上的震源機制解，以波形擬合方差最大的震源深度和震源機制解作為最佳結(jié)果，圖4給出了不同震源深度下波形互相關(guān)之間的關(guān)系。從圖4可以看出，震源機制解在不同深度上變化不大，均顯示為走滑型地震，震源深度6km處，波形互相關(guān)系數(shù)最大。我們得到此次地震最佳矩心深度為6km，與臺網(wǎng)初始定位深度相差1km。圖5為固定震中位置垂直深度方向搜索西藏仲巴地震震源機制解對應的波形擬合圖，可以看出方差減少量VR為0.71，其中，參與反演的5個臺站波形的15個分量中，波形擬合方差減少量達0.5以上的有8個，占總分量的0.53%，可見臺站波形擬合較好。

為了減少初始震中定位偏差對反演結(jié)果的影響，我們以初始震中和最佳矩心深度為中心點(25號網(wǎng)格點)，在矩心深度平面網(wǎng)格區(qū)域內(nèi)搜索空間最佳震源機制解，圖6顯示了矩心深度平面上49個不同網(wǎng)格點上震源機制解及波形擬合相關(guān)系數(shù)值。由圖6可知，25號網(wǎng)格點對應的波形互相關(guān)系數(shù)最大，其對應的結(jié)果為最佳震源機制解，以最大震源球表示。矩心深度平面網(wǎng)格搜索的最佳雙力偶機制解為，節(jié)面Ⅰ的走向/傾角/滑動角為109°/85°/-177°，節(jié)面Ⅱ的走向/傾角/滑動角為19°/88°/-4°，最佳矩心位置為30.590°N和83.784°E。在矩心深度平面網(wǎng)格區(qū)域內(nèi)搜索最佳震源機制解對應的波形擬合圖，方差減少量VR為0.71，與固定震中位置搜索結(jié)果一致，同時也表明，對于中小地震大多數(shù)初始破裂點與破裂的矩心位置是基本一致的(趙翠萍等，2008)。

圖 6 不同網(wǎng)格點上震源機制解及波形擬合相關(guān)系數(shù)五角星為震源機制搜索網(wǎng)格中心點在地表上的投影點；彩色條紋表示波形波形擬合相關(guān)系數(shù)；最大的震源球位置表示平面網(wǎng)格空間搜索的最佳震源機制解

為了估算震源機制反演斷層面解參數(shù)的不確定度，我們采用刀切法(Jackknifing method)對反演的參數(shù)進行不確定度分析(Sokos et al，2013；Boyd et al，2015)。考慮到不同方位和震中距地震記錄對反演結(jié)果的影響，我們對參加反演的5個臺站，每次減少1個臺站共5次；每次減少1個分量共15次，合計進行20次反演，并把所有的反演結(jié)果投影在同一震源球上(圖7)，由圖7可以看出，斷層面節(jié)面線分布相對集中，說明多次反演結(jié)果比較穩(wěn)定。

圖 7 斷層面參數(shù)不確定度的估計黑色圓為P 軸；白色圓為T 軸

圖 8 利用P波初動資料單獨反演(a)和P波初動和波形約束反演(b)得到的仲巴地震的震源機制解陰影區(qū)為ISOLA反演結(jié)果對應的壓縮區(qū)；實心方形表示P 波初動向上；空心方形表示P 波初動向下；黑色圓為P 軸；白色圓為T 軸

2.2 P波初動極性循環(huán)掃描結(jié)果

首先，我們采用FOCMEC方法，通過網(wǎng)格搜索，計算Crust2.0地殼速度結(jié)構(gòu)模型下的震源機制解。反演過程中，設(shè)置的P波初動的矛盾數(shù)為0，得到500組可接受解(圖8(a))。由圖8所示，由于可用的P波初動資料少，且大部分臺站分布集中，且初動向上，反演結(jié)果節(jié)面分布非常分散，很難判斷哪組解為最佳解。

為了得到最佳震源機制解，我們應用FOCMEC計算得到的多組走向角、傾角和滑動角值作為初動極性循環(huán)掃描法CSPS的輸入值，采用ISOLA中純雙力偶約束的固定震源機制解求解方式，聯(lián)合ZBA和GZE兩個臺站近震波形資料，計算觀測波形與理論波形擬合程度VR值。通過選取0.9*VRopt～VRopt范圍的解作為P波初動極性掃描計算得到的一組可行解，VRopt為最佳波形擬合震源機制解對應的方差減少量。本研究采用P波初動極性循環(huán)掃描法得到的最佳震源機制解為，節(jié)面Ⅰ的走向/傾角/滑動角為110°/67°/-172°，節(jié)面Ⅱ的走向/傾角/滑動角為17°/83°/-23°。

P波初動掃描計算的VRopt值為0.62，與全波形反演計算的VR值相比較小，波形擬合程度較低，但是通過P波初動循環(huán)掃描方法得到的可行解與ISOLA全波形反演結(jié)果相差不多，解的節(jié)面分布相對集中，P波初動方向與震源球體的壓縮和膨脹特性相對應(圖8(b))，能夠較好地反映小地震的震源特性。

總體來說，這2種方法給出的震源機制解具有很好的一致性。這在一定程度上驗證了利用P波初動聯(lián)合近震全波形求解震源機制解的可行性和可靠性。因此，我們認為應用極性掃描法方法約束反演震源機制解時，在臺站數(shù)量較少時也能得到穩(wěn)定可靠的結(jié)果，該方法對中小地震比較適用，有助于反演稀疏臺網(wǎng)記錄的中小地震震源機制解。

3 討論與結(jié)論

本研究利用ISOLA全波形反演方法，以西藏仲巴M4.7地震為例，計算稀疏臺網(wǎng)記錄的中小地震的最佳雙力偶震源機制解，節(jié)面Ⅰ：走向109°/傾角85°/滑動角-177°，節(jié)面Ⅱ：走向19°/傾角88°/滑動角-4°，最佳矩心位置為30.590°N、83.784°E，最佳質(zhì)心深度為6km，矩震級MW4.6。

圖9 震源區(qū)附近歷史地震震源機制解黑色圓點代表研究區(qū)域歷史地震(M≥5.0)；圓形為本次仲巴地震初始震中位置；五角星形為矩心位置；紅色沙灘球為本研究結(jié)果；黃色沙灘球為美國哈佛大學GCMT結(jié)果；黑線代表研究區(qū)內(nèi)的主要斷裂分布；歷史震源機制來自于哈佛大學GCMT結(jié)果

西藏仲巴M4.7地震發(fā)生在喜馬拉雅山西段，該區(qū)域地質(zhì)構(gòu)造較為復雜，發(fā)育著多條大型斷裂(圖9)。從地質(zhì)力學的角度來看，這個地方是構(gòu)造應力易集中的地區(qū)，也是地震的易發(fā)區(qū)，震中區(qū)及周邊發(fā)生的地震均為淺源地震，震源機制類型以正斷層性質(zhì)和走滑型性質(zhì)地震為主。此次震中位于帕龍錯-倉木錯和嘉黎斷裂交匯處。近年來，沿著這條斷裂帶曾發(fā)生過2004年7月6.7級地震、2005年4月6.5級地震(圖9)。震源機制結(jié)果表明，2017年1月4日西藏仲巴4.6級地震為走滑性質(zhì)，這一結(jié)果與震源區(qū)域附近歷史地震震源機制解具有相同性質(zhì)。

為了比較近震全波形反演結(jié)果的穩(wěn)定性和可靠性，我們進行多次反演，并對反演結(jié)果的不確定度進行定量分析，研究表明在少量臺站記錄波形參加計算的震源機制解存在較大的不確定性，反演過程中，不能單單依靠實際波形和理論波形擬合程度來判斷結(jié)果的好壞，我們需要結(jié)合臺站初動資料或更多臺站參加反演來提高結(jié)果的準確性。

本文還采用P波初動和近震波形聯(lián)合約束反演本次地震的震源機制，給出了最佳震源機制解，并與ISOLA近震全波形反演給出的結(jié)果具有較好的一致性。結(jié)果表明，當?shù)卣鹫鸺壭?，特別是臺站較稀疏，初動資料少且方位角覆蓋差，僅僅依賴P波初動資料很難確定震源機制解，在少量近臺觀測資料的情況下，聯(lián)合少量臺站的三分量波形記錄進行約束反演，也能得到可靠的震源機制解。

ISOLA方法和現(xiàn)在較為廣泛使用的CAP方法相同，均屬于全波形反演方法，所不同的是反演時CAP方法是分別對體波和面波分別進行反演(張廣偉等，2014；梁姍姍等，2017)，而ISLOA方法是整個波形參與反演，并可將震源設(shè)置為單點源，也可以設(shè)置為多點源，且是在四維空間中搜索最優(yōu)解。

綜上，對于稀疏臺網(wǎng)中小地震震源機制求解時，利用P波初動資料和近震波形聯(lián)合約束反演方法是一種簡潔明確的震源機制解量化方法，該方法可以穩(wěn)定、可靠地計算中小地震的震源機制解。

致謝：本研究使用了由希臘帕特雷大學提供的ISOLA近震全波形矩張量反演軟件，大部分圖件使用GMT軟件包繪制，在此一并表示感謝。