周成， 馬張烽， 蔣彌

1 河海大學(xué)地球科學(xué)與工程學(xué)院， 南京 211100 2 中山大學(xué)測繪科學(xué)與技術(shù)學(xué)院， 廣東 珠海 519082

0 引言

2017年4月5日6時9分(協(xié)調(diào)世界時，UTC)，伊朗東北部桑塞菲德(Sang sefid)地區(qū)發(fā)生了MW6.1地震，震中位于35.776°N 60.436°E，震源深度13 km(https:∥earthquak-e.usgs.gov/earthquakes).據(jù)當(dāng)?shù)孛襟w報道，此次地震造成兩人死亡，數(shù)十人受傷，靠近震中西北部的是伊朗第二大人口城市馬什哈德(Mashhad).此次地震附近沒有發(fā)現(xiàn)有明顯位移的已知斷層，分析構(gòu)造背景認(rèn)為此次地震屬于一次板塊間俯沖帶地塊的擠壓造成的盲沖事件(Aflaki et al.，2019)，俯沖帶的擠壓對Binalud斷層帶的東南部和Kopeh Dagh斷層帶的東部(Hollingsworth et al.，2008)造成了一定的沖擊；從大范圍的構(gòu)造機制來看，地震的發(fā)震原因可以歸結(jié)為：阿拉伯板塊相對于歐亞板塊以31 mm·a-1北移(DeMets et al.，2010)，與南移的歐亞板塊擠壓碰撞造成板塊間分布式縮短與剪切作用，同時阿拉伯板塊的北移導(dǎo)致伊朗東部的Lut地塊隨之北向推進，進而與桑塞菲德地區(qū)堅硬巖石層擠壓碰撞而造成桑塞菲德內(nèi)部剪切應(yīng)力作用發(fā)震(Jackson，1992).

伊朗桑塞菲德地區(qū)被一系列WNW-ESE走向逆沖走滑斷層包圍(圖1)，且由于阿拉伯板塊的北移俯沖運動，造成了該區(qū)域存在強烈的地殼運動、構(gòu)造變形與地震事件.其中桑塞菲德北部是伊朗東北部的主要構(gòu)造變形區(qū)，主要活動斷裂以一條WNW-ESE走向線型山脈Kopeh Dagh為主，這條由阿拉伯板塊北向運動孕育的長約600 km的逆沖系統(tǒng)隔絕了伊朗與土庫曼斯坦地區(qū)；桑塞菲德西部為著名的Binalud斷層帶和Mashhad斷層帶與伊朗東北部的一系列WNW-ESE走向逆沖走滑斷層共同吸收了絕大部分的阿拉伯板塊匯聚速率，這條NW走向范圍約130 km的Binalud山脈是伊朗地區(qū)構(gòu)造變形最嚴(yán)重的斷層帶之一(Shabanian et al.，2012)，它也是穩(wěn)定的伊朗中部板塊與活躍的東北部的邊界，形成了伊朗東北部活躍的構(gòu)造帶(Shabanian et al.，2009，2010)；南部是以Lut地塊為主的N-S走向的剪切帶(Foroutan et al.，2014)，它在阿拉伯板塊的推動下逐漸北移與桑塞菲德地塊擠壓扭轉(zhuǎn)造成Lut地塊東西側(cè)逆向運動，形成了Lut地塊東西兩側(cè)N-S走向的右旋走滑斷層以及北側(cè)Khaf和Jangal逆沖系統(tǒng)，這也被認(rèn)為是伊朗—阿富汗東部邊界地震帶的主要構(gòu)造(Walker and Jackson，2004).然而由于桑塞菲德地區(qū)的沙漠性氣候且地處偏遠，又未發(fā)生過大型事件導(dǎo)致了該地區(qū)構(gòu)造資料的不足，也阻礙了學(xué)者們對該地區(qū)盲斷層的識別.

伊朗桑塞菲德地震發(fā)生后，一些國內(nèi)外學(xué)者和專業(yè)機構(gòu)對此次地震的構(gòu)造背景和發(fā)震斷層的傾向進行了研究(表1)對于發(fā)震斷層的傾向目前仍存在爭議.Ghayournajarkar等(2020)使用哨兵一號和ALOS-2的InSAR數(shù)據(jù)獲取了精確的同震形變場，并利用PSO(Particle Swarm Optimization)非線性優(yōu)化反演方法推斷該地震的盲斷層為向東北傾斜的右旋逆沖斷層.Xu 等(2018)使用哨兵一號數(shù)據(jù)用于同震反演，通過分析恢復(fù)的二維形變場并研究了主震庫倫應(yīng)力的分布，推斷該地震的發(fā)震斷層為向西南傾斜的逆沖斷層.Su 等(2019)基于均勻斷層模型采用拉普拉斯平滑約束的方式反演推斷斷層滑動分布主要為右旋剪切位移，進而判斷該地震的發(fā)震斷層為東北傾向的逆沖兼具右旋走滑斷層.Aflaki 等(2019)綜合分析余震分布與震區(qū)現(xiàn)場地質(zhì)特征考察結(jié)果，推斷該地震發(fā)震斷層為西北-東南走向且向東北傾斜的逆沖斷層.綜上，不同學(xué)者對于桑塞菲德地震的斷層構(gòu)造機制的認(rèn)識充滿了爭議，其中的主要原因在于基本都是使用傳統(tǒng)的二步法反演斷層滑動分布特征，這也導(dǎo)致獲取的滑動分布特征的單一，無法分析不同反演方法之間的各向異性結(jié)果，以得到更全面的印證.因此，本文采用相同的斷層模型和地表觀測數(shù)據(jù)作為約束，運用不同的反演方法更辯證地分析滑動分布的合理性，進而為該盲斷層的構(gòu)造機制提供更為全面的認(rèn)識.

本文首先根據(jù)震源機制解(USGS，GCMT，IRSC)，結(jié)合獲取的InSAR升降軌同震形變場的空間變化形態(tài)及余震分布，假設(shè)了發(fā)震斷層?xùn)|北傾和南傾的兩種傾向可能，并基于彈性半空間構(gòu)建了單一平面斷層模型，使用二維InSAR形變數(shù)據(jù)為約束，采用基于約束最小二乘的動力學(xué)反演方法(Wang et al.，2008)以及貝葉斯反演方法(Amey et al.，2018)結(jié)合解算得到潛在發(fā)震斷層非線性斷層參數(shù)和斷層面上的滑動分布特征.最后，本文具體討論了反演得到的斷層參數(shù)和滑動分布特征之間的各向異性，并結(jié)合余震空間分布特征以及地質(zhì)構(gòu)造背景進行綜合分析，對桑塞菲德盲沖斷層的構(gòu)造機制進行了解釋.

1 InSAR數(shù)據(jù)及盲斷層傾向確定

1.1 InSAR數(shù)據(jù)

本研究采用兩對Sentinel-1A寬幅數(shù)據(jù)和兩對ALOS-2影像分別獲取了桑塞菲德地震的同震形變場.其中，選取了哨兵一號T13升軌和T93降軌，ALOS-2 T171升軌和T64降軌數(shù)據(jù).圖1為伊朗東北部構(gòu)造背景以及桑塞菲德主震位置及余震分布，表1為本文使用SAR影像詳細參數(shù).

圖1 伊朗桑塞菲德地震構(gòu)造背景(a) 伊朗地區(qū)相對于歐亞大陸參考框架二維GPS速度場; (b) 本文研究區(qū)域. 其中震源球為不同機構(gòu)和學(xué)者提供的震源機制解，深黃色點為IRSC(Iranian Seismological Center)記錄的桑塞菲德主震后兩個月內(nèi)的余震(http:∥irsc.ut.ac.ir)，紅色線條來自USGS記錄的伊朗已知斷層(https:∥catalog.data.gov/dataset/major-faults-in-iran-flt2cg)，KDF為Kopeh Dagh斷層帶，MF為Mashhad斷層帶，BF為Binalud斷層帶，KF為Khaf斷層帶，JF為Jangal斷層帶，深藍色為伊朗地區(qū)2006年至2015年相對于歐亞大陸參考框架二維GPS速度場(Khorrami et al.，2019).Fig.1 Tectonic background of Sang sefid earthquake in IRAN(a) GPS velocities relative to a Eurasian reference frame in Iran; (b) The Area Of Interest of this study. The focal mechanisms indicate different results from different institutions and scholars，The dark yellow dots represent the aftershocks two months after the mainshock recorded by IRSC (Iranian Seismological Center；http:∥irsc.ut.ac.ir)，The red lines represent faults from the USGS World Energy Project (https:∥catalog.data.gov/dataset/major-faults-in-iran-flt2cg) ，Abbreviations for faults are KDF：Kopeh Dagh fault，MF：Mashhad fault，BF：Binalud fault，KF：Khaf fault，JF：Jangal fault，The dark blue vectors are the GPS velocity field relative to a Eurasian reference frame for the period 2006—2015 in Iran.

表1 本文使用的Sentinel-1A數(shù)據(jù)參數(shù)Table 1 Parameters of Sentinel-1A data used in this study

首先利用精密軌道以及SRTM 1 sec數(shù)字高程模型采用幾何配準(zhǔn)的方式進行哨兵數(shù)據(jù)參考影像的配準(zhǔn)，然后在進行差分干涉處理之前對配準(zhǔn)好的堆棧進行增強譜分集(Enhanced Spectral Diversity，ESD)處理，以校正子帶中每個burst重疊部分的相位不連續(xù)(Ma et al.，2019)，從而使得burst重疊區(qū)域的相位差小于3°；在進行精配準(zhǔn)之后，利用分辨率為30 m的數(shù)字地形高程模型(SRTM 1 sec)進行去地形處理，然后采用傳統(tǒng)的Goldstein濾波器(Goldstein and Werner，1998)對干涉圖進行濾波處理，以抑制相位噪聲，提高干涉圖信噪比，使干涉條紋更光滑，從而減小相位解纏的難度.最后使用SNAPHU方法(Chen and Zebker，2000)對干涉相位進行二維解纏，ALOS-2數(shù)據(jù)采用的L波段受電離層與軌道誤差影響較大，在解纏后采用二階曲面去除電離層與軌道誤差代表的長波信號，并將干涉相位地理編碼為WGS84地理坐標(biāo)，最終得到清晰的同震干涉條紋圖和同震形變圖(圖2—3).

圖2為進行D-InSAR處理后的Sentinel-1A衛(wèi)星和ALOS-2 SAR傳感器獲取的升軌和降軌同震形變干涉圖，它清晰地反映出了桑塞菲德地震引起的地表形變，整個升降軌干涉條紋圖空間形態(tài)特征明顯，大致呈西北-東南走向橢球牛眼形分布，同震干涉圖分布均勻，覆蓋范圍達到30 km×20 km.干涉條紋沿包裹中心向外擴散，較為稀疏，形變梯度較小，但也足夠表明地震引起了明顯的地表位移；其次，升降軌干涉圖都為單葉干涉條紋，說明桑塞菲德地震引起的地表形變以沿視線向單一隆升或遠離，表明此次地震發(fā)震斷層以逆沖為主的斷層性質(zhì).

1.2 盲斷層模型確定

當(dāng)?shù)卣饦?gòu)造背景已知或者發(fā)震斷層確定時，利用大地測量觀測數(shù)據(jù)(InSAR，GNSS)進行單斷層或多斷層參數(shù)反演的研究已經(jīng)趨于成熟(Kuang et al.，2019；Feng et al.，2020；屈春燕等，2017)，但當(dāng)遇到中型地震(6≤MW≤6.5)且地表破裂情況未知，發(fā)震構(gòu)造背景模糊時，確定發(fā)震斷層參數(shù)尤其是傾向一直是一個難點，這也阻礙了盲斷裂地帶的地震危險性分析(Biggs et al.，2006；Materna et al.，2019).針對如何在多數(shù)地震情況未知時確定盲斷層的傾向和位置這個問題，本文提出了相應(yīng)的解決方案.眾所周知，InSAR數(shù)據(jù)具有很高的空間分辨率，這也為確定斷層的位置提供了十分寶貴的信息，通過InSAR觀測數(shù)據(jù)和重新定位的余震信息可以假設(shè)出潛在的盲斷層位置，然后由InSAR觀測數(shù)據(jù)反演約束出最佳擬合斷層面，進而獲得最佳斷層模型.

表2為國內(nèi)外學(xué)者和科研機構(gòu)反演獲得的桑塞菲德地震的震源機制解.其中，USGS，GCMT和IRSC認(rèn)為該地震的發(fā)震斷層具有兩種潛在的可能斷層，主要為東北傾和南傾.不同學(xué)者對該發(fā)震斷層的傾向也有不同的認(rèn)識，Xu等(2018)認(rèn)為該盲斷層向西南傾的可能性更大，斷層為西北-東南走向，方位角和傾角分別為120°和40°；具有爭議的是其他大多數(shù)學(xué)者都認(rèn)為該盲斷層為東北傾，且反演結(jié)果的滑動角差異較小，均反映了該斷層逆沖為主且具有右旋分量的特性；但不同學(xué)者給出的方位角和傾角均有一定的差異，方位角最大差9.2°，傾角最大差9.9°，這也造成了該發(fā)震斷層具體傾向的模糊.綜合分析不同學(xué)者給出的結(jié)果，可以認(rèn)為桑塞菲德地震的發(fā)震斷層主要有東北傾和南傾兩種潛在可能，關(guān)于本文的結(jié)果將做進一步分析和解釋.

通常，發(fā)震斷層的地表跡線投影位于干涉條紋的邊緣附近.因此，根據(jù)桑塞菲德地震獲取的同震形變干涉圖，本文假設(shè)了WNW-ESE走向和W-E走向兩種潛在斷層(圖2).同時，根據(jù)USGS和IRSC給出的震源機制解，設(shè)置了斷層初始走向、傾角和滑動角搜索范圍，并假設(shè)地球是一個均勻的彈性半空間，泊松比為0.25，基于Okada85彈性位錯模型(Okada，1985)構(gòu)建了東北傾和南傾兩種單一平面斷層模型，斷層模型深度為15 km，面積為30 km×20 km，以覆蓋整個滑動區(qū)域.

圖2 基于InSAR數(shù)據(jù)獲取的同震形變干涉圖(a)，(b) 哨兵一號的升(左)降(右)軌干涉圖； (c)，(d) ALOS-2的升(左)降(右)軌干涉圖.其中黑色和紅色分別代表東北傾斷層和南傾斷層地表跡線及震源機制解.Fig.2 The wrapped interferogram of the co-seismic deformation based on InSAR data(a)，(b) The wrapped interferogram of Sentinel-1A ascending and descending track; (c)，(d) The wrapped interferogram of ALOS-2 ascending and descending track. The black and red represent the fault plane and focal mechanisms of NE-dipping and South-dipping fault respectively.

表2 不同學(xué)者和機構(gòu)反演的斷層參數(shù)Table 2 Fault parameters retrieved by different scholars and institutions

2 二維形變場分解與滑動分布反演

2.1 二維形變場分解

通過目視判別桑塞菲德地震同震形變場(圖3)，可以觀察到升降軌的同震形場在空間變化形態(tài)具有一定的相似性.從分布范圍來看，升降軌同震形變場都呈西北-東南向分布，這與該區(qū)域北部一系列剪切斷層(Binalud斷層和Kopeh Dagh斷層)平行，且形變場的分布區(qū)域也大致相似；從分布特征來看，升降軌同震形變場都呈橢球牛眼形狀，且都展現(xiàn)了沿衛(wèi)星視線方向(line of sight，LOS)方向隆升位移的特征；其中，升軌和降軌形變圖各自展現(xiàn)出沿視線向靠近衛(wèi)星的最大位移量達10.2 cm和13.4 cm；對于視線向位移，降軌的東向位移會對視線向位移產(chǎn)生正向作用，升軌的東向位移會對視線向位移產(chǎn)生負(fù)作用(Wang et al.，2017)，而降軌的視線向隆升要稍大于升軌，因此可知桑塞菲德地震引起的地表形變以東向為主.從地震活動物理機制來看，該地震升降軌同震形變場沿視線方向一致性隆起的特征，也更加印證了此次地震發(fā)震斷層的逆沖特性.

圖3 地震同震形變場LOS向解纏干涉位移(a)，(b) 哨兵一號的升(左)降(右)軌; (c)，(d) ALOS-2的升(左)降(右)軌. 其中黑色和紅色分別代表東北傾斷層和南傾斷層地表跡線及震源機制解.Fig.3 Unwrapped Interferometric displacement in LOS direction of earthquake event(a)，(b) The unwrapped interferometric displacement of Sentinel-1A ascending and descending track; (c)，(d) The unwrapped interferometric displacement of ALOS-2 ascending and descending track. The black and red represent the fault plane and focal mechanisms of NE-dipping and South-dipping fault respectively.

因為升降軌視線向位移不具備極性特征，僅使用視線向位移無法確定桑塞菲德地震的發(fā)震斷層，本文采用結(jié)合應(yīng)力應(yīng)變模型的方差分量估計方法(Liu et al.，2018)聯(lián)合哨兵和ALOS-2數(shù)據(jù)恢復(fù)二維形變.根據(jù)衛(wèi)星的成像模式可知，雷達視線向位移由東西向、南北向和垂直向位移構(gòu)成，而由于哨兵一號和ALOS-2的近極軌道飛行(與經(jīng)線差10°左右)，導(dǎo)致雷達觀測對南北向位移極不敏感且計算矩陣奇異，因此為保證東西向和垂直向的分解精度，將南北向位移做忽略不計處理，此時將雷達視線分解如式(1)：

(1)

式中DLOS為視線向位移，UEU為分解矩陣，λ為雷達入射角，θ為衛(wèi)星升降軌衛(wèi)星方位角，vE和vU為東西和垂直位移.通過以上的關(guān)系式，并結(jié)合應(yīng)力應(yīng)變模型(袁霜等，2020；He et al.，2019)即可解算出東西和垂直位移.由于桑塞菲德地震區(qū)大多以荒漠巖石層為主，相干性較高，因此在解算過程中，為提高解算精度，提取相干性高于0.4的像素點進行解算，得到二維形變場(圖4).

分析二維形變場可以發(fā)現(xiàn)，桑塞菲德地震引起了明顯的東向位移和隆升位移，其中隆升位移主要位于震源中心附近(圖4b)，最大值達到13.4 cm，并且在隆升區(qū)域西南部和東北部展現(xiàn)出輕微的下沉；東向位移主要位于形變場的東北部(圖4a)，最大值達到6.8 cm，在形變場的西南部也顯示了部分西向位移，最大值達到4.1 cm；其次，分析本文假設(shè)的兩種潛在斷層，水平形變場主導(dǎo)的東向位移即為兩種潛在斷層上盤的東向運動，可解釋為東北傾盲斷層的右旋運動或者南傾盲斷層的左旋運動；而水平形變場西南部的西向運動量與范圍相對于東北部的東向運動較小，分析局部構(gòu)造認(rèn)為這可能是由于震區(qū)南部Lut地塊向北推擠造成的塊體擠壓縮短而導(dǎo)致的東西形變場的極性特征(圖1)，這與假設(shè)東北傾盲斷層為發(fā)震斷層時的右旋運動并不矛盾.兩種潛在斷層都可以很好的解釋二維形變場，因此，僅從二維形變場仍然無法確定盲斷層的傾向，還需進行反演分析.

2.2 滑動分布反演

2.2.1 約束最小二乘動力學(xué)反演

本文采用解算升降軌一維形變獲取的二維形變場數(shù)據(jù)為約束，并基于彈性半空間構(gòu)建的東北傾和南傾兩種矩形位錯模型，沿長度和寬度將斷層劃分成30×20個矩形塊，矩形位錯大小設(shè)置為1 km×1 km，使用SDM(Steepest Descent Method)反演獲得了與地表測量數(shù)據(jù)具有最佳吻合度的非線性斷層參數(shù)和發(fā)震斷層面上的最優(yōu)滑動分布特征.InSAR的高空間分辨率數(shù)據(jù)雖然有利于描述斷層細節(jié)，但也加大了反演難度.因此，在反演之前，為提高計算效率，減小因數(shù)據(jù)冗雜而產(chǎn)生的收斂問題，對獲取的二維形變場數(shù)據(jù)進行四叉樹降采樣(Jónsson et al.，2002)處理，采用GBIS軟件(Bagnardi and Hooper，2018)提供的降采樣包，設(shè)置形變梯度閾值為0.005 m，并根據(jù)相干性加權(quán)以提高反演精度，充分利用質(zhì)量較好的數(shù)據(jù)點，最終降采樣得到2645個數(shù)據(jù)點參與反演.

圖4 二維形變場(a) 東西方向形變; (b) 垂直方向形變, 其中黑色和紅色分別代表東北傾斷層和南傾斷層地表跡線及震源機制解.Fig.4 2-dimensional displacement field(a) East and West displacement; (b) Vertical displacement. The black and red represent the fault plane and focal mechanisms of NE-dipping and South-dipping fault respectively.

經(jīng)過近10000次迭代達到收斂，得到兩種潛在斷層的非線性斷層參數(shù)(表3)和滑動分布特征(圖5a—b).分析滑動分布特性可以發(fā)現(xiàn)，東北傾和南傾斷層破裂都未到達地表，這與野外地質(zhì)考察結(jié)果一致.東北傾斷層滑動較為集中，沿走向6～22 km處，主要盲沖作用在7 km深度(圖5a)，最大滑動量達到1.09 m；且反演得到的平均滑動角為128.80°(表3)，顯示該斷層是以逆沖為主的兼具右旋走滑的逆斷層，這與二維形變場分析的結(jié)果保持一致.南傾斷層盲沖作用主要分布在3～11 km深度，沿走向8～20 km處，最大滑動量出現(xiàn)在7 km深處達到0.69 m(圖5b)；平均滑動角為77.50°(表3)，表明該斷層是以逆沖運動為主的逆斷層，這也與二維形變場分析的結(jié)果一致.另外，反演得到震源參數(shù)與各機構(gòu)(USGS，GCMT，IRSC)基本一致，其中各幾何參數(shù)的誤差主要是因為各機構(gòu)給出的矩震級是利用主震發(fā)生時地震波反演而得，它描述的是主震發(fā)生時的瞬時強度，而桑塞菲德地震發(fā)生后的很短時間內(nèi)在主震附近發(fā)生了至少三次MW5.0以上的余震(Aflaki et al.，2019)，這會導(dǎo)致主震的發(fā)震斷層進一步破裂，而InSAR數(shù)據(jù)又是在震前和震后獲取的靜態(tài)數(shù)據(jù)，其中ALOS-2的升軌數(shù)據(jù)的時間基線達到154天，這必然也會涵蓋了由于余震導(dǎo)致的震后位移以及震前的部分非構(gòu)造信號；其次，ALOS-2數(shù)據(jù)在研究區(qū)右上角的不足也會導(dǎo)致InSAR空間分辨率的降低而造成一定反演結(jié)果的差異.

圖6和圖7為反演獲取的二維形變場的模擬值和殘差值，從殘差圖可以看出東北傾斷層和南傾斷層模型與二維同震形變場都有著較好的適配，兩種潛在斷層模型的擬合殘差都較小，其中，東北傾和南傾斷層模型反演得到的二維形變值和模擬值的平均RMS(Root Mean Square)殘差分別為東西向0.42 cm和0.50 cm，垂直向0.31 cm和0.49 cm，觀測數(shù)據(jù)與模型的相關(guān)系數(shù)分別為0.997和0.992.可見SDM反演的結(jié)果更支持東北傾斷層為可能的發(fā)震斷層，但由于兩種斷層模型擬合性都較好，本文將基于此斷層模型進行貝葉斯反演滑動分布，獲取最優(yōu)滑動分布特性的后驗概率密度，并做進一步驗證分析.

圖5 SDM反演和貝葉斯反演的滑動分布(a) SDM反演的東北傾斷層滑動分布; (b) SDM反演的南傾斷層滑動分布; (c) 貝葉斯反演的東北傾斷層滑動分布; (d) 貝葉斯反演的南傾斷層滑動分布. 其中藍色實心點為USGS-NEIC記錄的主震后三個月內(nèi)(2017年4月5日至2017年7月5日)MW≥3.5的余震數(shù)據(jù)，黑色五角星為主震震中.Fig.5 Slip distribution inverted by SDM and a Bayesian inversion method(a) Slip distribution of NE-dipping fault inverted by SDM; (b) Slip distribution of South-dipping fault inverted by SDM; (c) Slip distribution of NE-dipping fault inverted by Bayesian inversion method; (d) Slip distribution of South-dipping fault inverted by Bayesian inversion method. The blue dots represent the aftershocks three months (5 April 2017 to 5 July 2017) after the mainshock recorded by USGS-NEIC，the black star represents the epicenter of mainshock.

表3 SDM反演獲取的斷層參數(shù)Table 3 Fault parameters inverted by SDM

圖6 反演獲得的東北傾斷層模型的二維形變觀測值、模擬值和殘差Fig.6 2D observed，modeled and residual displacement for NE-dipping fault

圖7 反演獲得的南傾斷層模型的二維形變觀測值、模擬值和殘差Fig.7 2D observed，modeled and residual displacement for South-dipping fault

2.2.2 貝葉斯反演

利用貝葉斯方法反演大地測量數(shù)據(jù)集可以快速地獲取源斷層模型參數(shù)的后驗概率密度，進而獲取斷層模型上具有最大后驗概率的滑動分布特征(Amey，2018).根據(jù)貝葉斯定理，利用似然函數(shù)和先驗概率構(gòu)建后驗概率函數(shù)如式(2)：

(2)

式中，m為包含滑動量和滑動角的斷層模型參數(shù)，d為大地測量數(shù)據(jù)集(InSAR數(shù)據(jù))，p(m|d)為斷層模型參數(shù)的后驗概率，p(d|m)是描述大地測量數(shù)據(jù)與斷層模型參數(shù)匹配度的似然值，p(m)為斷層模型參數(shù)的先驗概率，p(d)為不考慮模型參數(shù)的標(biāo)準(zhǔn)化常量.反演之前，通過設(shè)置斷層模型參數(shù)的可能范圍來計算模型參數(shù)的先驗概率，其中東北傾斷層的滑動角的搜索范圍設(shè)置在90°～180°之間，南傾斷層的滑動角的搜索范圍設(shè)置在0°～90°之間，然后使用蒙特卡羅方法進行先驗樣本采樣，并通過對使用模型參數(shù)預(yù)測的觀測值與實際觀測值的擬合度加權(quán)計算似然函數(shù)如式(3)：

(3)

式中，N為InSAR數(shù)據(jù)點的總個數(shù)，d為InSAR觀測值，Σd為InSAR數(shù)據(jù)的方差-協(xié)方差，它描述的是InSAR觀測值之間的相關(guān)性，用來對InSAR數(shù)據(jù)加權(quán)，s為滑動量，G為基于彈性半空間內(nèi)由傾角和滑動角計算的設(shè)計矩陣(Bagnardi and Hooper，2018)，模擬形變值由G與s計算而得.反演通過迭代的方式進行，在每次迭代后設(shè)計矩陣G會隨著重新計算的滑動角而更新，進而會不斷更新似然函數(shù)以計算斷層參數(shù)的后驗概率，以獲取與觀測值擬合度最高的斷層滑動參數(shù).

采用SDM反演得到的具有最佳擬合的斷層模型作為貝葉斯反演的先驗斷層模型，通過迭代反演獲取了具有最大后驗概率的斷層參數(shù)(表4)和滑動分布特征(圖5c—d).其中，東北傾和南傾斷層的平均滑動角分別為129.07°和58.94°，最大滑動量分別為1.02 m和0.66 m(表4).滑動分布結(jié)果表明東北傾斷層逆沖兼具右旋走滑以及南傾斷層逆沖兼具左旋走滑的運動特性，東北傾斷層的滑動主要集中在5～9 km深度范圍，最大滑動量出現(xiàn)在7 km深處(圖5c)，與SDM反演結(jié)果基本一致；而南傾斷層的滑動集中在3～9 km深度范圍，最大滑動量分布在6 km深處(圖5d)，與SDM反演結(jié)果略有差異.其次，對比分析USGS-NEIC提供的震后的數(shù)次余震與滑動分布(圖5)，余震主要分布在東北傾斷層主震滑動分布范圍四周(圖5a，c)，與主震的滑動分布形成良好的互補，而在南傾斷層上卻顯示出與主震滑動范圍更多的重合(圖5b，d)，這與地震學(xué)余震機理不符，因此可以認(rèn)為USGS給出的余震結(jié)果與東北傾斷層更符合.綜上所述，由于采用相同的斷層模型與大地測量觀測數(shù)據(jù)，貝葉斯反演的結(jié)果與SDM反演具有一定的相似性，但仍然可以看到差異性，關(guān)于差異以及斷層傾向本文將做進一步分析和討論.

表4 貝葉斯反演獲取的斷層參數(shù)Table 4 Fault parameters inverted by a Bayesian inversion method

3 討論與分析

3.1 發(fā)震斷層傾向判斷

為了更好地討論桑塞菲德地震發(fā)震斷層的傾向，本文將約束最小二乘動力學(xué)反演方法與貝葉斯反演方法獲取的結(jié)果進行對比與討論.兩種反演方法得到的東北傾斷層參數(shù)基本一致(表3—4)，其中滑動角僅差0.27°，滑動量相差0.07 m；而南傾斷層參數(shù)的差異較大，滑動角相差18.56°，滑動量相差0.03 m.其次，兩種反演方法獲得的東北傾斷層的滑動分布特征也具有較高的一致性(圖5a，c)，滑動都集中在7 km深度上下，且滑動角的方向都顯示出了發(fā)震斷層逆沖兼具右旋走滑的運動性質(zhì).對于南傾斷層，利用SDM方法反演的滑動分布范圍更大(圖5b)，較圖5d破裂范圍要大2 km，且展現(xiàn)的幾乎是純逆沖的運動性質(zhì)，而貝葉斯反演滑動分布結(jié)果卻顯示出了更多的左旋分量(圖5d)；此外在11 km至13 km處顯示出部分余滑，與SDM結(jié)果有一定差異.利用地震學(xué)方法得出的震源機制解是研究地震性質(zhì)，剖析斷層模型的主要手段，大地測量數(shù)據(jù)與斷層模型的擬合度也是作為確定盲斷層模型決策的重要信息，反演殘差結(jié)果也顯示南傾斷層反演獲得東西和垂直向位移的殘差0.50 cm和0.49 cm都要大于東北傾斷層的殘差0.42 cm和0.31 cm，并且從殘差圖也可以清晰的看到南傾斷層在主要形變區(qū)較東北傾斷層顯示出了更多的系統(tǒng)性殘差(圖6—7)，因此認(rèn)為東北傾斷層為此次地震的發(fā)震斷層.

利用升降軌聯(lián)合解算的桑塞菲德地震二維形變場可以觀察到更多地表位移信息，東西和垂直形變場都顯示出WNW-ESE走向的形變區(qū)域(圖4a)，并且東西向形變場中主要東向位移和部分西向位移也表明了該發(fā)震斷層走滑特征，這與兩種反演方法獲取的東北傾斷層的滑動分布特征一致，而SDM反演的南傾斷層卻顯示出幾乎純逆沖的運動特性(圖5b)，與二維形變場特征不符；其次，東北傾向斷層的WNW-ESE走向與伊朗東北部斷層帶已知的逆沖活動斷層(Binalud斷層帶，Kopeh Dagh斷層帶)走向基本一致(圖1).

圖8 桑塞菲德地震余震分布黑色表示東北傾斷層，紅色表示南傾斷層.Fig.8 The aftershocks of Sang sefid earthquakeBlack and red planes represent NE-dipping and South-dipping respectively.

通過構(gòu)建桑塞菲德主震后進行重新定位的59次余震分布(Aflaki et al.，2019)與兩種潛在斷層相對位置模型(圖8)，可以看到余震在反演獲得的東北傾斷層模型震源中心附近分布較為密集，且余震的震源中心有沿著WNW-ESE走向發(fā)展的趨勢(圖8)，這與東北傾斷層的走向基本一致.從兩種潛在斷層與余震的相對深度位置來看(圖8)，余震的深度有隨著WNW-ESE走向斷層傾向而逐漸增加的趨勢，且余震的深度分布與東北傾向斷層模型更吻合.另外，發(fā)現(xiàn)重新定位的余震深度都系統(tǒng)性地大于本文反演的主震深度，這可能是由于該地區(qū)地殼介質(zhì)不均勻，巖石新生成且?guī)r體破碎導(dǎo)致反演算法假設(shè)的均勻介質(zhì)的地殼速度要大于實際地殼速度，從而使反演的主震深度偏高.

綜合以上原因，本文認(rèn)為伊朗桑塞菲德地震的發(fā)震斷層為WNW-ESE走向，向東北傾斜的逆沖兼具右旋走滑分量的盲斷層.同時也能說明對于構(gòu)造背景模糊又缺乏大地震事件記錄，現(xiàn)場考察地質(zhì)特征不明顯的盲斷層事件，僅利用高空間分辨率的InSAR觀測數(shù)據(jù)來判斷盲斷層的傾向仍然非常具有挑戰(zhàn)性，未來隨著InSAR技術(shù)觀測精度的進一步提高，或者更豐富的觀測數(shù)據(jù)如GNSS數(shù)據(jù)參與進行聯(lián)合約束，以及分布密度更高的余震數(shù)據(jù)輔助判斷，盲斷層傾向的判斷的可靠性將進一步提升.

3.2 盲斷層地震構(gòu)造意義

伊朗桑塞菲德地震發(fā)生在伊朗東北部，震區(qū)周圍分布著已知的活躍走滑斷層和逆沖斷層(圖1a)，但震區(qū)附近缺乏活動斷層以及重大地震事件的記錄造成了該震區(qū)構(gòu)造背景的模糊.眾所周知，斷層位置的確定是進行地震危險性分析的首要任務(wù)，因此桑塞菲德地震的發(fā)生以及該盲斷層的識別對于補充該地區(qū)地震構(gòu)造背景具有重要意義.

此次地震發(fā)生后短時間內(nèi)發(fā)生了數(shù)次MW5.0以上的余震，也充分表明了該盲斷層的活躍性，并且本文利用InSAR觀測數(shù)據(jù)反演的滑動分布結(jié)果也顯示該斷層在深度7 km處出現(xiàn)了大于1 m的明顯滑動，但都未使該斷層破裂至地表，這可能是由于桑塞菲德地區(qū)存在覆蓋層較厚的沉積物以及堅硬的巖石層(Shabanian et al.，2010)，這也可能導(dǎo)致本文反演獲得的震源深度存在高估以及地震震級存在低估，進而也解釋了圖8重新定位的余震深度分布(10～20 km)要比反演獲得的震源深度(7 km)更深的結(jié)果.其次余震整體空間分布呈廣泛分散形態(tài)，符合盲斷層活動斷裂帶余震分布的特征(Stein and Ekstr?m，1992)，因此該地區(qū)可能普遍存在盲斷層或次生斷層而形成活動盲斷裂帶，這樣一個未被識別的斷裂帶也對該地區(qū)的地震風(fēng)險造成了隱患，將來應(yīng)對構(gòu)造背景模糊的桑塞菲德地區(qū)做進一步關(guān)注與研究，以豐富該地區(qū)地震活動機制并掌握此活動斷裂的更多細節(jié).

3.3 桑塞菲德地震破裂模式

伊朗桑塞菲德地震被認(rèn)為是阿拉伯板塊北移與歐亞大陸板塊碰撞擠壓而造成的逆沖事件(Su et al.，2019)，局部來看認(rèn)為是伊朗東部Lut地塊北移與桑塞菲德地區(qū)堅硬的巖石層碰撞而孕育的剪切斜沖事件，因此該斷層破裂模式可以認(rèn)為是阿拉伯板塊與歐亞大陸板塊間俯沖帶逆沖斷層破裂(Hyndman and Wang，1993).板塊間俯沖帶逆沖事件通常有三個特征：同震時的斷層上盤前緣隆升、后端下沉和下盤后端下沉，這是由于地震破裂期間逆斷層上下盤相向運動，上盤前緣受到下盤的俯沖推覆而縮短隆起，后端受到拉張而下沉，下盤后緣受到上盤的仰沖擠壓而下沉，這種破裂模式造成的地表形變與本文獲取的二維形變場中垂直位移空間分布形態(tài)一致，主要表現(xiàn)為震區(qū)的主要隆升與西南部和東北部的輕微下沉.

震前的大地測量學(xué)以及地質(zhì)學(xué)研究表明伊朗東北部的構(gòu)造變形主要集中在塊體邊界處，包括Lut地塊邊界以及Kopeh Dagh構(gòu)造帶附近(Mousavi et al.，2013)，GPS速度場的結(jié)果也表明Lut地塊較歐亞大陸北移的速率已達到5.7 mm·a-1，這一北向作用幾乎都被伊朗東北部的一系列WNW-ESE的走滑逆沖系統(tǒng)所吸收(Walpersdorf et al.，2014)，這也是Kopeh Dagh與Binalud構(gòu)造帶西北走向褶皺地貌與逆沖構(gòu)造的主要發(fā)育成因(Lyberis and Manby，1999)，與此相似的是，桑塞菲德地區(qū)也以西北走向褶皺地貌為主且逐漸展現(xiàn)出背斜上升的趨勢(Aflaki et al.，2019)，此外桑塞菲德地區(qū)地處Lut地塊北部，受其北移俯沖的直接推擠碰撞作用效果，因此可以認(rèn)為桑塞菲德地區(qū)也正逐漸形成與Kopeh Dagh與Binalud構(gòu)造帶相似的構(gòu)造帶，即WNW-ESE方向走滑逆沖盲斷層為主的構(gòu)造系統(tǒng)，這也進一步印證了本論文認(rèn)為的桑塞菲德地震WNW-ESE盲斷層為發(fā)震斷層的判斷.由于發(fā)震斷層未破裂至地表，本文只能通過地表地貌特征以及構(gòu)造形變來分析桑塞菲德地震的破裂模式，將來隨著該發(fā)震斷層的構(gòu)造背景逐漸明晰，桑塞菲德斷層的破裂模式也將得到進一步解釋.

4 結(jié)論

本文利用哨兵一號和ALOS-2 SAR影像數(shù)據(jù)獲取了伊朗桑塞菲德MW6.1地震高空間分辨率的同震形變場，并結(jié)合四幅干涉圖和應(yīng)力應(yīng)變模型的方差分量估計方法恢復(fù)了二維形變場.通過構(gòu)建東北傾和南傾兩種潛在單一斷層模型，基于彈性半空間以二維形變數(shù)據(jù)為約束利用約束最小二乘與貝葉斯反演算法獲取了兩種潛在斷層的滑動分布特征，并結(jié)合余震數(shù)據(jù)和已知活動斷層跡線判斷了發(fā)震斷層的傾向，得到以下幾點認(rèn)識：

(1) 獲取的哨兵和ALOS-2升降軌同震干涉條紋都呈WNW-ESE走向橢球牛眼形分布，由于地震未破裂至地表，形變場相干性較好，條紋清晰，且升降軌都為單葉條紋表明此次地震斷層主要逆沖的特性.解纏后的同震形變場都以LOS向隆升為主.利用多源數(shù)據(jù)融合恢復(fù)的二維形變場主要顯示出震源附近的主要東向位移和隆升位移，以及西南部的少量西向位移，形變主要發(fā)生在斷層上盤，解釋為東北傾斷層的逆沖右旋走滑運動以及南傾斷層的逆沖左旋走滑運動，且水平方向的西向位移與震區(qū)南部Lut地塊的北移俯沖推擠構(gòu)造機理一致，二維形變場的空間分布特征也與俯沖帶逆沖斷層破裂模式吻合.

(2) 使用約束最小二乘動力學(xué)反演和貝葉斯反演的東北傾斷層的滑動分布具有較高的一致性，結(jié)果顯示滑動主要分布在地下5～9 km處并未破裂至地表，在地下7 km處孕育了超過1 m的最大滑動量，平均滑動角約128.8°，矩張量達1.427×1018N·m，相當(dāng)于震級MW6.07的地震，顯示出主要逆沖兼具右旋走滑分量的運動性質(zhì)；而兩種方法反演的南傾斷層的滑動分布結(jié)果在破裂范圍和滑動角具有一定的差異性；東北傾斷層WNW-ESE走向也與伊朗東北部已知的活動斷層平行(圖1a)，通過將余震在斷層深度上的投影也發(fā)現(xiàn)余震的分布與東北傾斷層滑動分布具有更好的符合，并且余震深度的擴展方向也與東北傾斷層的傾向保持一致，反演殘差也顯示東北傾斷層具有更好的擬合結(jié)果.綜合分析東北傾斷層與已知活動斷層構(gòu)造關(guān)系，地震學(xué)反演獲取的震源機制解的重要性，余震分布和滑動分布合理性結(jié)果認(rèn)為東北傾斷層為桑塞菲德地震的發(fā)震斷裂.

(3) 通過分析余震的空間分布形態(tài)，震前伊朗東北部地塊構(gòu)造位移以及桑塞菲德地區(qū)褶皺背斜隆升的地貌認(rèn)為該地區(qū)可能普遍存在盲斷層而形成盲斷裂帶，桑塞菲德地震的發(fā)生也讓該活動斷裂帶受到關(guān)注，未來還應(yīng)繼續(xù)觀察該地區(qū)的長期震后松弛效應(yīng)以獲取更細節(jié)性的構(gòu)造特征，為獲取該地區(qū)置信度更高的地震危險性評估為服務(wù).本文著眼于對盲斷層的識別，發(fā)現(xiàn)僅通過地表觀測數(shù)據(jù)來識別盲斷層依舊充滿挑戰(zhàn)性，而采用不同反演方法進行比較分析，再通過余震分布和已知活動斷裂資料加以輔助有助于對盲斷層構(gòu)造機制的研究，也為進一步對盲斷層的識別研究提供了新思路.

致謝本文所使用的Sentinel-1數(shù)據(jù)由ESA/Copernicus提供.GFZ德國地學(xué)中心汪榮江老師提供了SDM約束最小二乘反演軟件.Dr.Amey提供的SlipBERI貝葉斯滑動分布反演工具包.Dr.Bagnardi提供的GBIS貝葉斯反演軟件工具包.中山大學(xué)地球科學(xué)與工程學(xué)院周宇教授和日本東北大學(xué)理學(xué)院地球物理系Dr.Ghayournajarkar對本文提出的寶貴建議.中南大學(xué)地球科學(xué)與信息物理學(xué)院劉計洪博士提出的關(guān)于應(yīng)力應(yīng)變模型的建議.本文中的全部圖形由GMT 6.1.0繪制.在此一并感謝.