2019年巴基斯坦新米爾普爾MW6.0地震的同震形變場與斷層滑動分布反演

賈 蕊 張國宏 解朝娣 單新建 張迎峰 李成龍 黃自成

1)中國地震局地質研究所，地震動力學重點實驗室，北京 100029 2)云南大學地球科學學院，地球物理系，昆明 650500

0 引言

2019年9月24日11時01分(Universal Time Coordinated，UTC)，巴基斯坦北部距新米爾普爾(New Mirpur)7km處發生MW6.0地震。美國地質調查局(United States Geological Survey，USGS)利用遠場地震波數據進行測定，給出的地震震中位于(33.078°N，73.794°E)，震源深度為11.5km(表1)。

表1 不同機構給出的2019年12月24日巴基斯坦新米爾普爾MW6.0地震震源機制解Table1 Source parameters of the MW6.0 earthquake in New Mirpur，Pakistan，24 December，2019 from different agencies

本次巴基斯坦MW6.0地震很可能是由于印度板塊和歐亞板塊碰撞系統內的逆沖斷裂造成的，巴基斯坦北部、印度和阿富汗鄰近地區的地震活動性和活斷層都是印度板塊與歐亞板塊碰撞系統的直接結果。印度板塊持續以45～52mm/a的速率向N移動，從而在該碰撞區引起地殼縮短和變形，這使得包括巴基斯坦西部和北部在內的喜馬拉雅地區成為世界上地震活動性最活躍的地區之一(鄧起東等，2014)。喜馬拉雅造山帶作為印度板塊與歐亞板塊碰撞形成的產物，主要由4條斷裂帶組成，從北向南依次為藏南拆離系(Southern Tibetan detachment system，STDS)、主中央逆沖帶(main central thrust，MCT)、主邊界逆沖帶(main boundary thrust，MBT)以及主前緣逆沖帶(main frontal thrust，MFT)(楊曉平等，2016)。本次地震的震中靠近喜馬拉雅主前緣逆沖次級斷裂(圖1)。MFT是地球上最大的活動碰撞造山帶中最長的活動收縮構造帶，大體呈EW走向。自全新世以來，MFT作為板塊邊界斷裂帶，吸收了兩大板塊間至少 1/3 的碰撞縮短量。如此大規模的高速率會聚活動使得MFT斷裂帶上的地震活動頻度高、強度大，歷史上曾發生過多次7級以上大地震，如1897年印度西隆8.7級地震、2005年巴基斯坦克什米爾7.8級地震以及2015年尼泊爾8.1級地震等，這些地震分別發生在MFT的不同段落上(吳中海等，2016)。陸地上的活動斷層造成了多次大型破壞性地震，產生了地表斷裂，促進了構造地貌的發育。同震地表破裂的詳細觀測不僅為了解地震本身和估計地震重現周期提供了基本信息，而且為解釋其他地區構造地貌特征提供了幫助。然而，逆斷層地表斷裂比走滑斷層或正斷層更少見。2005年克什米爾地震是喜馬拉雅地區歷史上第1次地表斷裂事件，也是該地區已知的最大地震。在2019年巴基斯坦6.0地震后，至今還沒有學者使用InSAR技術觀測研究此次地震。研究此次地震為評估活動逆斷層的地震危險性以及了解巴基斯坦北部的地震活動性提供了難得的機會。

圖1 2019年12月24日巴基斯坦新米爾普爾MW6.0地震周邊構造背景圖Fig. 1 Tectonic background of the MW6.0 earthquake in New Mirpur，Pakistan，24 December，2019.a 區域構造背景圖； 紅色實線代表活動斷層，紅色沙灘球代表本研究所得的震源機制解，黑色沙灘球代表GCMT確定的本次地震的震源機制解，黃色沙灘球代表USGS確定的本次地震的震源機制解，黃色圓圈代表1919年至今4級及以上地震震中，黑色箭頭代表板塊的運動方向，黑色小方框表示SAR數據的覆蓋范圍。b 紅色五角星指示本次地震震中(USGS)，白色 方框內為軌道T100和T107的Sentinel-1A影像范圍

如圖1b所示，該地區在歷史上發生過多次中強地震。在全球范圍內，平均每年將發生約10次7級以上破壞性地震及超過100萬次破壞性有限甚至沒有地表破裂但仍能被地震儀器記錄到的小地震(震級范圍為2.5～4.5級)，而在這兩者之間，更有成百上千次中強地震(震級范圍為5.5～6.5級)發生，破壞性顯著。中強地震往往容易被忽視，只有極少數由于震源深度很淺而造成人員傷亡及結構破壞，得以被人們關注。構造應力如何沿著斷層累積并轉化為大地震引發能量的突然釋放，這仍然是一個懸而未決的問題。對大量中強地震的震源破裂以及破裂過程開展研究，可以更詳細地了解構造載荷如何從靜止期發展至活躍期，對分析地震活動性以及大地震的突然爆發之間的關系也有一定參考價值。

利用區域或遠場地震數據可以反演出中強地震的震源機制，用于分析目標地震的斷層特征。震源機制解提供2個節面，然而由于中強地震往往不會產生地表破裂，故很難從2個節面中區分出發震斷層。盡管地震臺站的儀器時間采樣精度非常高，但其空間覆蓋率往往不足以區分出真實的震源機制解。綜上所述，針對中強地震震源斷層的進一步研究具有實際意義，需要利用衛星及大地測量等傳統地震學以外的更獨立的研究方法。作為應用最廣泛的衛星大地測量技術之一，干涉SAR(InSAR)衛星影像是獲取與地震有關的同震變形的常用手段。InSAR雖然時間采樣率低，視衛星重訪時間長短而定，但具有非常高的空間采樣率。InSAR的同震變形精度可達2～3cm，能夠得到清晰的中強地震引起的地表變形結果。InSAR形變數據具有上升和下降等多種觀測模式，可協助識別地表破裂及震源斷層面。近年來，使用InSAR技術反演盲斷層地震的發震斷層參數和滑動分布已成為研究盲斷層地震的重要手段。此外，使用InSAR重建的同震形變場能夠獲得高精度的同震觀測結果，該方法已成功應用于2015年MW6.4皮山地震(Fengetal.，2016)、2016年門源地震(Lietal.，2016)以及2017年MW6.5九寨溝地震(單新建等，2017)等地震的研究中。

從表1 可看出，2個機構利用遠場地震波數據給出的巴基斯坦MW6.0地震的震源機制解差異較大，并且無法得知此次地震真實的發震斷層。因此，使用InSAR技術重建巴基斯坦新米爾普爾MW6.0地震的同震形變場，并反演這次地震的發震斷層參數與滑動分布、深入了解其破裂特征與發震構造，對了解該區域內的逆沖斷裂帶活動性與發震構造機制有重要意義。此外，在新米爾普爾北部的吉拉姆河(Jhelum)上建有一座曼格拉大壩(Mangla)，該水庫位于上新世—更新世西瓦利克(Siwalik)系砂巖和黏土層交替形成的低矮山丘上。西瓦利克堆積物是由河流沉積的淡水沉積物，在中新世時期開始的一次大的抬升使其恢復了活力。西瓦利克系下的默里(Murree)地層為同一地槽中新世時期的咸淡水和淡水沉積。較老的結晶巖形成了海拔較高的內喜馬拉雅山脈，而較年輕的默里和西瓦利克系形成了外喜馬拉雅山脈或山麓丘陵。該水庫距離巴基斯坦新米爾普爾MW6.0地震的震中不到10km，因此該水庫對該地震是否有觸發作用是值得注意的一點。在本研究中，我們將結合水庫周圍的地質構造以及地震活動性對該問題加以討論。

本研究中，我們利用歐洲航空局(European Space Agency，ESA)獲取的升降軌Sentinel-1A SAR數據重建了此次地震的InSAR同震形變場，并在升、降軌形變場的聯合約束下反演獲得了發震斷層參數和同震滑動分布，從而對這次地震造成的形變場特征、斷層運動性質和發震構造機制進行分析與探討，為進一步研究巴基斯坦北部的地震構造活動性提供一定依據。

1 InSAR同震形變觀測

我們獲取了2019年9月24日巴基斯坦MW6.0地震的Sentinel-1A衛星升、降軌SAR影像(TOPS模式)。震前Sentinel-1A升軌SAR影像的觀測日期為2019年9月16日，震后升軌影像的觀測日期為2019年9月28日； 震前Sentinel-1A降軌SAR影像的觀測日期為2019年9月5日，震后降軌影像的觀測日期為2019年10月11日。SAR影像干涉對的垂直基線分別為-44m(升軌)和-9m(降軌)(表2)。

表2 Sentinel-1A衛星影像參數信息Table2 Details of Sentinel-1A interferograms used in this study

基于D-InSAR技術，我們使用了GAMMA軟件平臺進行干涉處理。采用美國宇航局發布的SRTM 3sec數字高程模型(SRTM DEM 90m)消除地形相位誤差。處理時對干涉圖進行了10×2(距離向×方位向)多視處理以降低干涉圖的噪聲。采用最小費用流算法進行相位解纏(Werneretal.，2000； Gallandetal.，2016)，先對高質量的區域進行解纏，得到可靠的解纏相位值，然后以這些可靠相位值構建相位模型，并以此為基礎對其他低相干區域進行解纏(Goldsteinetal.，1998； Farretal.，2007)。最后，將解纏后的干涉影像地理編碼至WGS84坐標系下，得到升、降軌同震InSAR干涉條紋和LOS向形變場，如圖2 所示。

圖2 2019年12月24日巴基斯坦新米爾普爾MW6.0地震升、降軌同震形變場Fig. 2 The coseismic deformation field of the MW6.0 earthquake in New Mirpur，Pakistan，24 December 2019 observed by Sentinel-1A.紅色沙灘球為USGS給出的震源機制解，紅色五角星為此次地震的震中，天藍色為水庫水體的范圍

根據升、降軌同震形變場的長軸方向可確定發震斷層的走向為NW-SE，其結構呈非對稱性，分布于喜馬拉雅前緣斷裂的次級逆沖斷裂帶北緣。斷層SW盤的干涉條紋明顯密集于NE盤(圖2a，c)，且SW盤的形變值為正、NE盤為負(圖2b，d)。同一地區的升、降軌形變量相同，圖2 顯示出該發震斷層的運動性質以逆沖運動為主，觀測獲得的衛星視線向(LOS)最大隆升形變量約為0.1m(升軌)和0.05m(降軌)，最大沉降形變量約為0.03m(升軌)和0.05m(降軌)。由于形變中心北部有一座水庫，震中附近失相干嚴重，導致形變場條紋不連續，對獲取對此次地震的震源性質產生了一定影響。

2 斷層滑動模型的反演

同震觀測數據過多會增加計算量且降低數據的信噪比，使反演結果難以收斂。對于形變梯度較大的區域，采用均勻采樣法能有效降低其對整體結果的影響。因此，在反演之前，為避免數據量過大導致計算效率降低，我們采用均勻采樣方法對升降軌InSAR形變場進行降采樣處理。最終，分別得到4i195(升軌)和3i311(降軌)個形變數據點用以反演發震斷層面上的滑動分布。

巴基斯坦新米爾普爾MW6.0地震屬于未破裂到地表的盲斷層事件，發震區域的地質與地震活動性研究較少，難以從現有資料中精準地確定發震斷層的位置。分析InSAR形變場可知(圖2)，形變區域大致沿NW-SW向延伸，再結合震中區域的地質構造特征猜想本次地震的發震斷層有2種可能，即傾向SW的斷裂帶或傾向NE的斷裂帶(MFT的次級斷裂)。因此，我們建立了2種初始斷層模型，并對初始斷層模型的幾何參數進行多次測試，搜索斷層參數的最優解以反演斷層面上的精細滑動分布。然后，根據觀測和模擬的擬合情況以及不同發震斷層面的滑動情況來推測這次地震的可能發震斷層。

2.1 模型1(傾向SW)斷層滑動模型及滑動分布反演

我們將升、降軌形變場的觀測數據作為約束條件，根據模型1獲取斷層面上的精細滑動分布。首先，固定模型1的斷層跡線位置，再測試不同的斷層傾角對數據的擬合情況，得出當傾角為90°時擬合效果相對較好，擬合度達89%。升、降軌數據的均方根誤差約為0.01m，擬合殘差≤2cm。圖3 為模型1的觀測-模擬-殘差圖，由圖3c、f可看出擬合度較好。模型1反演得到的最佳斷層參數詳見表3。得到的最佳傾角為90°，但這種情況在自然界中實際并不存在，逆沖型地震的傾角一般為30°～50°，現已知發生過的逆沖型地震的最大傾角為65°，該模型不符合逆沖型地震的客觀規律，因此排除了該模型。

圖3 模型1(傾向SW)的升、降軌觀測值、模擬值和殘差Fig. 3 Observation，simulation and residual from ascending(a—c)and descending(d—f)InSAR data and inversion by model 1.黑色虛線為模型1的斷層跡線，紅色五角星為此次地震震中(USGS)

表3 均勻滑動反演發震斷層的幾何和運動學參數(模型1)Table3 The source parameters for the Pakistan earthquake(model 1)

2.2 模型2(傾向NE)斷層模型及滑動分布反演

將InSAR同震形變場作為約束條件，根據發震斷層的幾何參數，采用模擬退火法反演發震斷層的同震滑動分布。模擬退火法是一種非線性直接反演的優化算法，能夠避免結果為局部極大值而非全局最大值的情況(齊浪等，2011)。同樣以升、降軌形變場的觀測數據作為約束，獲得模型2斷層面上的精細滑動分布。采用LOG函數確定最佳傾角和平滑因子(馮萬鵬等，2010)。首先，根據USGS給出的震源機制解固定NE傾的斷層傾角，并不斷平移斷層的位置，將擬合殘差最小的斷層位置作為最佳發震斷層的位置，然后再固定該發震斷層，在不斷調整斷層傾角的同時改變平滑因子，并計算每次反演結果的均方根誤差(root-mean-square error，RMSE)(式(1))。設定的傾角范圍為2°～60°，以2°為步長，設定平滑因子的范圍為0～0.3，以0.05為步長，所得結果如圖4 所示。

(1)

其中，zobs為觀測值，zmod為模擬值(馮萬鵬等，2010)。

根據反演結果來看，模擬退火法在基于InSAR數據計算發震斷層參數方面的效果較好，所得結果基本滿足了反演需要。

圖4 傾角與RMS的關系圖Fig. 4 Diagram of inclination and RMS.a 傾角與平滑因子的曲線關系； b 傾角與平滑因子的函數關系

圖5 模型2(傾向NE)的升、降軌觀測值、模擬值和殘差Fig. 5 Observation，simulation and residual from ascending(a—c)and descending(d—f)InSAR data and inversion by model 2.黑色虛線為模型2的斷層跡線，紅色五角星為此次地震震中(USGS)

圖6 模型2(傾向NE)斷層面上的滑動分布Fig. 6 Slip distribution of the northeast-dipping fault model 2.

表4 均勻滑動反演發震斷層的幾何和運動學參數(模型2)Table4 The variable source parameters for the Pakistan earthquake(model 2)

由圖4 可看出，當傾角約為10°～20°平滑因子約為0.05時擬合效果最好，擬合度達89.37%。圖5 為升、降軌擬合效果圖，無論是升軌還是降軌，擬合得到的形變場都有明顯的一升一降盤，擬合誤差≤4cm，符合一般逆沖型地震的特點。根據所確定的斷層最佳參數，畫出模型2的斷層滑動分布圖，如圖6 所示，滑動整體呈橢圓形集中分布。滑動分布沿走向集中在10～14km范圍內，沿發震斷層面傾向主要集中在10～20km內，滑動分布表現為逆沖性質。該模型反演得到的最佳斷層參數詳見表4。分析認為，該斷層模型與此次地震的發震斷層特征更相符。但反演得到的矩震級與USGS給出的震級差距較大，這可能是由于此次地震的震級及形變區域的形變量太小，未能被InSAR測量捕捉。大氣延遲、軌道、DEM等誤差也對反演結果產生了影響。

3 討論

基于此次新米爾普爾MW6.0地震的InSAR同震形變場數據，反演計算得到的震中位置為(33.09°N，73.75°E)，與USGS確定的結果相差4.48km，與GCMT所得結果相差68.65km。經分析認為，USGS和GCMT給出的發震地點和震源機制解是基于遠場地震波資料反演得出的，但由于臺站資料具有不確定性，導致其給出的震中位置差異較大。USGS等機構在確定震源深度時主要利用了地震波的走時，但這類方法對地震臺的間距要求很高，只有當最小震中距與震源深度達到一定比例時，所得的震源深度才能足夠準確。對于6級以上地震，遠震體波能夠較準確地測定出震源深度，但對于5級以下的中小地震而言，遠場臺站數據信噪比較低，不利于準確測定震源深度。

總體而言，傳統的地震學方法在確定震源參數時對臺站的分布密度及觀測數據的質量等要求很高。對于中小地震，地震學方法費時費力，且結果誤差較大，不適合進行震源參數研究。InSAR是一種新興的大地觀測技術，具有分辨率高、形變敏感度高且覆蓋范圍廣等特點，以InSAR數據為基礎可利用彈性位錯模型獲取震源參數，所反演的斷層參數和滑動分布結果更具可靠性。中小地震的形變量較小，而InSAR測量技術的精度很高，能識別出mm級的形變。InSAR空間分辨率較高，能對發震位置起到很好的約束作用。此外，近年來發展后的時序InSAR技術能很好地處理震前、震后時間序列的形變場，在很大程度上提升了中小地震地面形變的觀測精度，為研究中小地震斷層面三維幾何形態和力學屬性提供更加可靠的地面位移數據，為研究發震區域的活斷層的運動性質和地震周期提供更加清晰的認識。

圖7 同震滑動矢量與地形(a)、重力(b)的分布關系圖Fig. 7 Diagram of co-seismic sliding vector(a)and gravity(b)distribution.黃色沙灘球是本次地震震源機制解，藍色和黑色箭頭均為滑動矢量

新米爾普爾MW6.0地震也為研究青藏高原南緣的大規模構造提供了機會。震中區域處于歐亞板塊與印度板塊碰撞區，印度板塊持續向N推擠，使得該碰撞區發生變形。由圖7 可知，歐亞板塊與印度板塊之間的相對運動方位角雖然變化不大，但逆沖滑動方位角確實發生了變化。因此，我們推測除了驅動板塊邊界運動的力之外，重力也控制著逆沖運動的方向(Copley，2012)。當山地與低地之間存在重力勢能差則意味著山地會變形以降低勢能，從而導致橫向擴展至剛性下邊界之上，這種響應重力驅動力的橫向運動被稱為重力流(Copleyetal.，2015)，在重力驅動力的作用下地殼向垂直于地形等高線的下坡方向移動，且導致地殼沿青藏高原南緣的弧形山脈徑向逆沖。因此，此次地震的滑動矢量結果表明，重力驅動力對于控制青藏高原南緣的弧形邊緣而言十分重要(Copleyetal.，2007)。

新米爾普爾MW6.0地震北部約10km的吉拉姆河上建有一座中型水庫及曼格拉大壩。該水庫區域的地質構造為基底巖石上覆蓋著厚約3km的西瓦利克系易碎砂巖和黏土，這些巖石對施加的荷載幾乎沒有抵抗作用，并且會逐漸變形以適應荷載，這樣的地質構造特征導致該水庫沒有出現明顯的地震效應，在水庫發生過的幾次大地震都與其東北部的斷層有關。另外，在此次地震的發生地新米爾普爾建有地震觀測站，對蓄水前、后記錄到的地震事件的頻率進行分析，發現蓄水導致地震活動性出現明顯增加，然而經進一步分析發現，地震的增多實則是因為觀測儀器的靈敏度提高了，而非水庫蓄水觸發了地震(Brown，1974)。Brown(1974)認為未來地震活動性可能由于構造力的累積再次增加，而非由于水庫水位波動引起的荷載變化。綜合水庫區域的地質構造分析以及古地震與歷史地震的研究成果，我們認為此次地震并非由水庫蓄水觸發的。

4 結論

本文以InSAR升、降軌數據為約束，重建了2019年9月24日巴基斯坦MW6.0地震同震形變場，通過對不同斷層模型的反演試驗，給出了發震斷層可靠的幾何參數和滑動分布特征，并得到了以下結論：

基于2019年巴基斯坦地震的升、降軌SAR數據，得到了該地震的同震形變場，其形態呈現明顯的一升一降兩盤分布，顯示此次地震為逆沖型，兼具少量右旋走滑，LOS向最大隆升量為0.1m，最大沉降量為0.05m。巴基斯坦MW6.0地震沒有出現地表破裂位置的失相干條帶，這可能是由于此次地震發生在低傾角的印度板塊和歐亞板塊俯沖帶上，而板塊俯沖帶上逆斷層的運動和破裂方式都與板內的逆斷層不同。斷層面的傾角較緩，斷層破裂時上盤向S仰沖，北盤向N俯沖，由于擠壓推覆構造造成了上盤的前端隆升，而其后端受到拉張成為沉降區(Savage，1983； Hyndmanetal.，1993； 屈春燕等，2017)。僅根據InSAR數據無法準確得到發震斷層的傾向，因此我們建立了傾向SW和傾向NE 2種斷層模型，并分別進行多次反演實驗，再結合USGS等機構給出的震源機制解及巴基斯坦北部區域的地質構造背景信息，最終認為傾向NE的斷層模型與此次地震的發震斷層特征更相符。該模型表明，2019年巴基斯坦地震的斷層模型為主前緣逆沖斷裂(MFT)的次級斷裂。斷層走向為304°，傾角為15°，最大滑動量為0.64m，平均滑動角為124.85°，震源深度為2.77km，反演得到的矩震級為MW6.0。

此次地震產生的破裂并未出露地表，是一次盲斷層型的低角度逆沖活動。根據反演結果，分析認為該地震是由喜馬拉雅主前緣逆沖斷裂的次級斷裂引起，再結合區域地貌特征推測此次事件是印度板塊向歐亞板塊持續低角度俯沖所造成的。東部西太平洋-菲律賓板塊向W的俯沖作用、西部印度板塊向N的俯沖碰撞作用使得中國大陸地殼強烈變形。汪素云等(1996)發現造成這種地殼強烈變形的主要動力來源于印度板塊，印度板塊與歐亞板塊的碰撞作用使中國大陸地殼呈現出 “西強東弱”的形變特點，大多數的地震活動集中在青藏高原及其鄰區(吳中海等，2016)。位于青藏高原南部的喜馬拉雅造山帶是印度板塊與歐亞板塊碰撞的產物，其規模巨大，且正持續高速運動，由此導致高頻率、大強度地震的不斷發生。因此，了解喜馬拉雅造山帶上的地震活動特征有利于研究青藏高原的地殼形變特征，同時也有助于探索中國未來的大地震潛在的危險性以及活動趨勢，助力中國的防震減災事業。

致謝歐洲航空局(ESA)為本研究提供了Sentinel-1A雷達衛星影像； 德國GMZ的汪榮江老師提供了斷層滑動分布反演SDM程序包； 本文部分圖件使用GMT軟件繪制。在此一并表示感謝!