盲逆斷層震后變形耦合機制研究
——以2017年伊朗MW7.3地震為例

王洵,王衛民,何建坤,周云

1 中國地震局地震預測研究所,北京 100036 2 中國科學院青藏高原研究所大陸碰撞與高原隆升重點實驗室,北京 100101 3 中國地震局地球物理研究所,北京 100081

0 引言

高精度長周期地表變形監測的最終目標是從動力學角度揭示變形產生的機理,從而為地質災害預警提供現實指導和理論依據(Wang et al.,2009; Gualandi et al.,2017).強震發生后,余震和斷層面的震后蠕滑表明地震發生后巖石圈還在延續應力狀態調整,甚至蠕滑釋放的能量能達到主震釋放的能量20%以上,遠遠超過余震釋放的能量的總和(Savage and Svarc,1997; Freed,2007),震后蠕滑在地震周期中扮演著重要的角色,被認為影響局部甚至更大范圍內的地震觸發效應(Lohman and McGuire,2007; Feng et al.,2020).研究表明產生震后蠕滑的物理機制不同,所引起的應力變化也不同,而且對走滑型地震和傾滑型地震,同一物理機制如黏彈松弛效應所產生的應力變化極性與同震應力變化存在差異(Freed,2005),即對于傾滑型地震而言,某些落在同震應力“陰影區”的區域可能在較長時間的下地殼黏彈松弛后成為應力加載區域(Nostro et al.,2001).因此,探討震后變形機制,可以為人們評估重復地震地震危險性、“地震空區”等問題提供理論依據.

當前震后形變研究更多聚焦于大型走滑斷裂帶和俯沖帶強震周期的研究,這些區域發生破壞性強震的概率大,構造清晰、形變信號較強、便于空間對地觀測技術的測量(Wen et al.,2019; Huang et al.,2020; Zhao et al.,2021; Fukuda and Johnson,2021).對于盲逆斷層強震而言,不便于確定發震斷層幾何形態,而且多處于大型塊體之間,構造復雜,自然環境也不便于人工勘測.近年來,隨著全球定位系統(Global Positioning System,GPS)及合成孔徑雷達干涉測量(Interferometric Synthetic Aperture Radar,InSAR)等觀測能力的提升,對陸內逆沖型強震如2008年汶川地震及2015年尼泊爾地震等也逐步開展震后變形機制、巖石圈流變結構的研究(Shao et al.,2011; Diao et al.,2021).2017年11月12日18時18分17秒(UTC)在伊朗伊拉克邊境扎格羅斯造山帶山前褶皺內盲逆斷層內發生了一次MW7.3地震(隨后統稱伊朗地震),震中經緯度為34.911°N、45.959°E,震源深度約19 km(美國地質調查局,USGS),阿拉伯板塊的持續北東向推擠(Nilforoushan et al.,2003),在扎格羅斯斷褶帶內發育有一系列與兩個板塊碰撞方向垂直、呈NW-SE走向的逆沖推覆斷裂(Allen et al.,2011).推測此次地震的成因是持續的板塊運動碰撞,使該區域的應力不斷增加,最終導致深部基底滑脫(Barnhart et al.,2018).扎格羅斯山前斷褶帶西北段地震活動性較低,有完整地震目錄以來,在2017年MW7.3地震之前震源區約200 km內沒有MW6.0級以上地震發生,而震后約一年內,在扎格羅斯山前斷裂帶發生了兩次MW6.0以上的走滑型地震,分別是2018年8月25日MW6.0地震和2018年11月25日MW6.3地震(圖1),預示著需要加強關注該區域的地震危險性.

圖1 2017年MW7.3伊朗地震構造背景

伊朗地震InSAR觀測資料豐富,前人對主震的斷層滑動分布、和早期震后變形特征及產生機制進行了較深入的研究(Barnhart et al.,2018; Ding et al.,2018; Feng et al.,2018; He et al.,2019; Wang and Bürgmann,2020).這些研究的震源特征是基本一致的——逆沖略帶右旋分量,但滑動分布模型存在一定的差異,如滑動量的大小、滑動分布的深度、分布形態是集中還是分為兩個較大區域等、對震后變形機制解釋從一個或兩個方面進行探討,研究的時間尺度和震后余滑特征也存在爭議.為了從應力的角度系統地研究震后形變中各物理機制的作用,獲取一個詳細、高度自洽的震源模型及高精度震后變形是必要條件(Barbot et al.,2009; Hong and Liu,2021; Zhao et al.,2021).伊朗地震發生在盲逆沖斷層,為了更好地約束發震斷層的空間位置,本文使用這次地震的多視角InSAR同震形變資料,結合遠場地震波形基于分層地殼模型對地震震源過程開展研究,解算同震三分量(3D)位移場,采用InSAR時序分析方法(Chaussard et al.,2014; Tong and Schmidt,2016)獲取震源區震后形變特征,基于更精細的震源破裂模型探討震后余滑及巖石圈黏彈松弛對震后形變的貢獻,認識扎格羅斯山前逆沖斷裂擴展方式以及巖石圈形變模式.

1 同震震源模型

基于分層地殼模型,采用有限斷層方法反演研究地震的破裂過程(Hartzell and Heaton,1983; 王衛民等,2008),引入不同觀測資料約束震源過程模型,聯合反演獲取更準確的震源破裂的空間和時間信息(王洵等,2019).從美國地震學研究聯合會數據中心獲取伊朗地震的遠場臺站(震中距30°～90°)波形資料,按照臺站方位角均勻分布,挑選高信噪比的26個P波,24個SH波數據;從歐空局獲得Sentinel-1A覆蓋震源區的三個軌道多景雷達影像應用于同震位移場解算,其中升軌072影像應用于震后變形時間序列研究(表1).

表1 2017年伊朗地震Sentinel-1影像參數

1.1 D-InSAR形變場

星載InSAR具備全時全天候的優勢,可以獲取連續完整的地表變形信息,廣泛應用于地震、火山等研究(Massonnet et al.,1993; 單新建等,2002).目前歐空局運行的Sentinel-1系列雷達重軌拍攝間隔為12天或6天,減少時間失相干的影響,有助于及時獲得地震產生的高精度地表形變測量數據 (Elliott et al.,2015b; Wang et al.,2022).本文使用GMTsar程序(Sandwell et al.,2011a,b)對震前和震后影像對進行差分干涉處理,利用歐空局在衛星拍攝20天后提供的精密定軌星歷數據進行軌道校正,引入外部地形數據——分辨率為3″(弧度秒)的SRTM3 DEM去除地形效應,采用Goldstein濾波處理降低噪聲影響,得到伊朗地震的同震形變場干涉圖,相位解纏使用Snaphu模塊(Chen and Zebker,2002),最后地理編碼得到WGS-84大地坐標系下的地表形變場.

伊朗地震斷層近SE-NW走向,北東方向是地形復雜的扎格羅斯山脈,西南方向是阿拉伯平原,地勢較為平緩.在地震震中區域即扎格羅斯山前,獲得了完整清晰的InSAR干涉條紋圖(圖2),顯示伊朗地震造成的地表形變范圍約450 km2,在東北和西南區域產生了兩個形變符號相反的扇區,分界特征與扎格羅斯山前斷裂高度吻合.干涉圖條紋連續、平滑,相位解纏得到的震中區域地表形變質量較高,由于Sentinel-1是C波段雷達,一條完整干涉條紋對應約2.8 cm的LOS向形變量,此次伊朗地震震級較大,干涉條紋十分密集,而且在局部區域出現條紋紊亂的現象,因此,需要剔除相干系數不高的區域,對震源區近場位移資料重采樣處理,減少數據量的同時還保證了準確性.震源機制解(表2)表明,伊朗地震為逆沖型地震,由圖2a升軌干涉場可知,地震上盤雷達視線向(LOS)抬升量約93 cm,表明該區域地表運動垂直向上或向西運動;在上盤北東區域雷達視線向沉降約17 cm,表示該區域有垂向向下或東向位移.降軌雷達干涉場(圖2c)顯示有三個形變中心,LOS向抬升區域與較大的沉降區范圍及量級相當,結合升降軌干涉場特征可以初步判定同震位移具有較大的西向運動分量.

表2 伊朗地震震源機制解

圖2 伊朗地震InSAR同震形變干涉場

1.2 3D同震位移場解算

D-InSAR是一種側視測量技術,存在視線向(LOS)模糊問題,對不同方向的形變敏感性也存在差異(查顯杰等,2006; He et al.,2019).伊朗地震震源區雷達影像資料豐富,Sentinel-1雷達平臺提供了多組不同視角的測量影像,可以直接解算獲取地震產生的真實位移值.以垂向向上(U)、東向(E)、北向(N)三個形變分量來表達真實位移場,由成像(雷達右視)的幾何關系,并約定地物目標靠近雷達時雷達視線向位移dLOS為正,遠離雷達時dLOS為負,三分量分別用UU、UE、UN表示:

升軌模式:

dlos=UUcosθ-(UNsinφ+UEcosφ)sinθ,

(1)

降軌模式:

dlos=UUcosθ-(UNsinφ-UEcosφ)sinθ,

(2)

其中,θ為雷達入射角,φ為雷達飛行方向與南北軸夾角.應用D-InSAR技術獲取軌道號T072、T174、T006同震干涉場如圖2所示,分別設其同震LOS向位移場為LOS1、LOS2、LOS3,則由公式(1)、(2)得

(3)

式中,θi為雷達入射角,φi為雷達飛行方向與南北軸夾角,等式右邊第一個矩陣為地表三分量投影到LOS的方向矢量.因此由三個不同視角的LOS向位移值,即可由公式(3)直接解算得到同震三分量地表真實位移.然而,無論是升軌還是降軌雷達,其飛行方向與南北向夾角φi約為13°,相較垂向和東西向位移,南北向位移對LOS向的貢獻值要小～1個數量級,方向矢量矩陣求逆之后,南北向系數比垂向、東西向大～2個數量級,顯著放大了誤差.加之InSAR測量時本身對方位向(近南北向)的位移敏感性較差,分辨率較低,直接解算時南北向位移會出現較大誤差(劉計洪等,2022),因此不能直接用于斷層滑動分布反演,本文將聯立三對InSAR 測量LOS向形變場進行震源模型反演,進而與真實3D位移場對比分析.詳細3D位移場特征見討論部分.

1.3 震源破裂過程聯合反演

獲取覆蓋震源區的震前、震后的多視角SAR影像,聯合遠場波形數據對震源過程模型進行約束,反演獲取更準確的伊朗地震震源破裂的空間和時間信息.伊朗地震發生在伊朗—伊拉克邊界的扎格羅斯山褶皺沖斷帶內,由于缺乏近場地震臺站及GNSS觀測,本文結合中國地震臺網中心發布震中信息和InSAR同震位移場選定震中位置為45.856°E,34.750°N,由Crust2.0提取伊朗地震地殼分層信息(Bassin et al.,2000),震中位置海拔約580 m,地殼厚度41 km,其中軟沉積層厚度1 km,硬沉積層厚度1 km,P波平均速度為6.24 km·s-1,S波平均速度為3.48 km·s-1,平均密度為2.87 g·cm-3.

研究中采用分層介質模型(Crust2.0)計算靜態位移,廣義射線計算遠場體波,利用快速模擬退火技術實現震源模型的非線性反演(王衛民等,2005,2015).目標函數的數學表達式為

(4)

首先基于點源模型,利用遠場P波反演獲得伊朗地震震源機制解(表2),結合InSAR干涉圖(圖2),構建有限斷層模型模擬伊朗地震的發震斷層,為方便后續震后變形機制研究,擬構建較大尺寸的斷層模型:長115 km、寬85 km、走向335°、傾向10°劃分為5 km×5 km的391個子斷層,起始破裂點位于(13,11),深度為18 km.利用InSAR資料、和遠場波形資料聯合反演震源過程,由形變干涉場圖2可知,震源區的InSAR干涉圖質量較高,覆蓋完全,相干系數高,反演中可相應加大權重;較InSAR近場資料而言,遠場波形資料對斷層面的破裂滑動分布約束較低,反演時參數的變化對其影響較大,但可以很好地約束震源的時間過程,而且最大振幅能更好地約束斷層面最大滑動量.基于上述考慮,波形資料的與InSAR資料的權重范圍分別為(0.1～0.4),(0.9～0.6),試算步長為0.05,三個視角InSAR測量值反演中比重一致,當InSAR同震形變場地震遠場波形資料這兩種數據在反演中所占的權重比分別為0.75、0.25,光滑因子為0.15時,InSAR資料和波形記錄的觀測值和理論值擬合情況均較好.

圖3為斷層面的破裂滑動分布地表投影及三維模型,由圖示滑動矢量沿著斷層面往西南方向滑動,表現為低傾角逆沖略帶右旋分量,滑動分布較為集中.圖4為遠場波形擬合結果,P波和S波波形相關性均較好,整體擬合程度較高.圖5分別給出了三個不同軌道的InSAR同震形變位移場觀測值、有限斷層模型合成LOS向理論值以及二者的殘差情況,可以看出升軌及降軌影像資料整體擬合情況均較優,二者大部分區域殘差值約2 cm以內.

圖3 破裂滑動分布投影(a)和三維示意圖(b)

圖5 有限斷層模型LOS向同震形變場擬合

圖6給出了地震破裂的時間過程,顯示斷層破裂由震中往四周擴展,較大滑動分布區域位于震中南側,破裂沒有達到地表,大部分滑動在15 s內完成,整個破裂滑動持續約25 s.結合區域地質背景及相關研究,通過InSAR與遠場波形資料的聯合反演獲取2017 年MW7.3級伊朗地震震源模型,表明此次地震的發震斷層是扎格羅斯山前盲逆沖斷層,斷層東傾,滑動分布主要集中地下深度14～20 km處,寬度約40 km,最大滑動量～7 m,位于初始破裂點以南,計算得到地震累計釋放地震矩能量為MW7.30,是一次逆沖為主兼具右旋走滑分量的地震事件.地震上盤即扎格羅斯山前緣與阿拉伯板塊相互擠壓致使地殼縮短,是此次地震的主要成因,伊朗地震促進了扎格羅斯山逆沖褶皺帶的生長和擴展.

圖6 伊朗地震破裂時間過程

2 震后變形研究

2.1 震后變形時序分析

震后變形是地震活動周期中一個重要階段,對活動斷層不同階段變形的觀測和研究有助于我們認識地震的演化過程.為了獲取2017年MW7.3地震造成的震后地表形變,本文采用小基線集法(SBAS)進行InSAR時序分析.獲取軌道號T072震后雷達影像22景(2017-11-023—2019-04-11),圖7為以20180110號影像為主影像時其他影像的相對基線信息和時間信息,多數影像的相對基線小于50 m,所有影像對的相對基線長度均在200 m以內,有利于降低基線長度過長帶來的干涉誤差,考慮影像對的時間相干性和相對基線長度,為避免時間去相干性導致影像干涉質量較差帶來整體誤差,設定時間基線為75天,相對基線長度為200 m,配對成60個影像對(圖7),對干涉圖質量刪選后處理得到形變速率.

圖7 干涉影像對時空基線信息

圖8顯示地理坐標下伊朗MW7.3地震震后震源區的地表LOS向形變位移速率,位移速率較大區域雷達傾角約42.3°,扎格羅斯山前即地震上盤出現明顯的震后抬升現象,其整體抬升趨勢與T072軌道雷達所測量的同震變形十分相似,形變區域略往西南偏移,如圖8右下插圖所示AA′剖線的速率最大達約8 cm·a-1.表明伊朗地震震后應力調整顯著,主震沒有完全釋放該區域所累積的能量,震后形變對褶皺的隆升起到促進作用.

2.2 運動學余滑模型

在當前研究中,震后余滑和黏彈性應力松弛對震后地表短期形變的貢獻還無法有效區分.為深入分析震后變形機制,首先構建一個純運動學余滑模型來擬合地表震后約半年(20180510)變形特征,余滑斷層幾何模型與同震斷層模型一致,反演方法除了去掉波形反演項,地殼模型及其他參數均與同震反演一致.該模型假設余滑是產生地表變形的唯一因素,如果變形是由黏彈松弛或孔隙彈性回彈導致的,那么可能無法獲取較好解釋觀測的運動學余滑模型(Rousset et al.,2012).圖9為基于震后序列20180510形變場反演運動學余滑模型,斷層面余滑分布較為集中、平滑,基本上都在同震滑動區域周邊,往地表方向及西南方呈現互補的趨勢,而且理論震后形變整體特征與測量值吻合較好,二者殘差也相對較小;斷層面最大余滑位移約67 cm,反演得到矩震級約為MW6.7.推測主震半年內,主要由斷層余滑現象控制著震后地表形變.

圖9 運動學余滑模型反演

2.3 應力驅動余滑及黏彈松弛耦合模型

合理的同震破裂模型是構建有效的應力驅動余滑或黏彈松弛模型的基礎.Barbot等(2009)通過構建斷層余滑數值模型研究2004年MW6.0 Parkfield地震震后GPS形變特征,發現余滑釋放的地震矩能量超過了主震,而且斷層面余震分布受控于無震蠕變區域.為進一步分析應力驅動余滑機制及殼下物質黏彈效應對震后變形的貢獻,同時利用空間對地測量技術獲取的震后地表形變,構建合理的模型反演該區域下地殼和上地幔的流變學結構(邵志剛等,2008).本文主要嘗試利用應力驅動震后余滑與黏彈性松弛耦合模型分析震后形變的驅動機制,并初步探討該區域巖石圈黏滯系數的范圍.使用半解析程序Relax進行震后余滑及黏彈性松弛模擬計算,利用發震斷層滑動分布作驅動余滑序列的初始條件來估算速率強化參數,在模擬計算時假設已經發生較大同震滑動的子斷層是速率弱化區不參與震后蠕滑,將同震滑動量為厘米級的子斷層元設置為速率強化-弱化的過渡區(即滑動量設為0),以確保這些區域也可以發生蠕滑現象(Barbot et al.,2009).由于研究區地震剖面研究資料缺乏,沒有很好的速度結構對巖石圈模型進行約束,參考Crust2.0分層地殼模型、發震斷層深度及莫霍面深度,不考慮地表沉積層的影響,構建三層黏彈性模型:彈性上地殼厚25 km、黏性下地殼厚20 km,黏性下地幔厚20 km,采用上述聯合反演獲取的同震斷層滑動分布作為輸入模型,為了充分模擬斷層面下部的余滑,將同震斷層模型往下擴展至125 km,如圖10所示.前人研究表明震后較短時間內(<5年),各類黏滯模型計算的結果差異較小,一般也較多采用Maxwell體作為黏彈介質模型(邵志剛等,2007),本文利用Maxwell體計算伊朗地震的震后變形序列,與InSAR獲取的震后變形時間序列擬合分析,設下地殼及上地幔黏滯系數(η1/η2)范圍為5×1017～1×1020Pa·s,試算步長0.25×1017Pa·s,模擬計算中發現當η值小于5×1018時,模擬形變場無法重現InSAR觀測形變場特征.圖10所示為當摩擦參數(a-b)σ=1.3 MPa、V0=7 m·a-1,黏滯系數組合η1/η2為1×1019Pa·s/1.25×1019Pa·s時,震源區主震后0.5 a、1.0 a、1.42 a耦合模型中斷層面余滑分布及模擬的震后變形與同時期T072軌道InSAR時間序列對比情況,同位置AA′剖面模擬形變值較為接近同期InSAR觀測值,理論計算的整體形變特征與觀測值較為吻合;斷層面余滑分布顯示震后余滑主要分布在同震滑動區域的上方淺部,下部余滑相對較弱,震后1.42 a余滑最大值約1.1 m,通過分析震后約一年半內耦合模型中各物理機制的貢獻量,發現余滑對地表的震后形變貢獻超過85%,結合純運動學余滑模型可以推斷伊朗地震震后一年半時間內,震后余滑機制是震后變形的主要因素.

圖10 應力余滑及黏彈松弛耦合模型模擬斷層面余滑分布及地表震后變形

3 討論

扎格羅斯山前斷褶帶西北段地震活動性較低,有完整地震目錄以來,在2017年地震之前震源區約200 km內沒有MW6.0級以上地震發生,此次地震提供了一個新的視角來認識該區域的斷層特征.本文利用Sentinel-1雷達資料獲取了震源區完整清晰的同震形變場,聯合遠場波形數據反演得到此次地震的破裂過程,震源模型顯示斷層較大滑動分布深度14～20 km處,沿走向長度約40 km,有兩塊滑動分布峰值,與Chen等(2018)、Feng等(2018)斷層模型總體趨勢較一致,均為沿近南北向優勢擴展,本文與Chen等(2018)破裂時間過程較類似,本文滑動分布更深更集中,最大滑動量與初始破裂點處于同一深度,相對較大且深,發震斷層傾角在10°左右,與前人研究結果相差不大,走向略有差異;表明該褶皺帶內的盲逆斷層傾角十分平緩,與喜馬拉雅造山帶發生的尼泊爾地震傾角較接近(王衛民等,2015; 張勇等,2015),而扎格羅斯山脈褶皺帶其他部位的逆沖地震的震源機制解顯示,發震斷層傾角均在30～60°(Elliott et al.,2015a),扎格羅斯造山帶是阿拉伯板塊與歐亞板塊的碰撞造山帶,也是地球上最年輕的陸-陸碰撞造山帶,至今構造作用活躍(McQuarrie,2004),不同區域逆沖斷層傾角存在較大的差異可能預示著扎格羅斯造山帶分段處于不同的演化階段.

通過InSAR技術獲得視線向形變值,分解得到三分量變形值,分析垂直地表變化和水平運動信息,可以進一步剖析地震的發震特征及機理(Hooper and Zebker,2007; Xu et al.,2020).D-InSAR只能獲取視線向的一維地表形變量,僅由一幅含視線向位移信息的差分干涉圖無法定量分析地表的真實位移,因此需要利用多視角差分干涉圖解算得到地表形變三分量信息(Elliott et al.,2016).由于雷達平臺的限制,大多數震例研究均基于同一衛星平臺不同飛行姿態(升軌、降軌)的測量資料,再聯合方位向位移或理論位移模擬值解算得到3D位移場(孫建寶等,2006).伊朗地震InSAR觀測資料豐富,本文利用Sentinel-1平臺提供的三對不同視角干涉影像,將三個LOS向位移值代入對公式(3),等式兩邊分別左乘方向矢量的逆矩陣,即可得到三分量地表真實位移場,如圖11所示.由圖11可知震源區垂直方向上的位移基本上分為兩個中心區域,東北區域向下沉降,最大值約50 cm,西南區域向上隆升,最大值約65 cm,隆升區域范圍較大約25 km×30 km,與逆沖斷層運動模式相吻合;東西方向:主要表現為西向運動,分為東北與西南方向兩個中心區,最大西向位移45～55 cm,西北和東南區域靠近斷層處有小范圍的東向位移,最大值約5 cm.利用上述聯合反演中基于有限斷層震源模型獲得的三分量理論值與觀測值擬合對比發現(圖11),垂向及東向位移場特征一致,位移極值也較接近,整體殘差較小,表明聯合反演獲得的斷層滑動分布模型高度自洽,同時表明采用合適的計算策略聯立多組InSAR觀測獲取地表真實三分量位移場,直接進行滑動分布反演,可以減少反演所需的計算資源.其中東向、垂向形變范圍及數值大小與前人研究(Ding et al.,2018; He et al.,2019)基本一致,另外,解算的南北向位移誤差較大,但整體以向南運動為主,同垂向和東向的規律性形變不一致,在震源區沒有獲得形變中心,這與理論值差異較大.利用上述聯合反演獲取的震源模型,計算得到同震南北向位移最大值約0.65 m(整體向南運動),與He等(2019)利用ALOS-2及Sentinel-1五對干涉影像,采用多種方法聯立解算的南向位移整體形態較為一致,小于其最大值(0.76 m),原因可能是InSAR測量中含有部分震后形變信息的影響,及與震源模型、地形信息、沉積層等有關.由解算的真實三分量位移值及理論值可知,伊朗地震以逆沖為主兼略帶右旋走滑分量,地震上盤總體往西南運動,最大垂直隆升約65 cm左右,對扎格羅斯山前褶皺帶起到生長和伸展作用.

圖11 伊朗地震同震解算位移場及理論值

對于大陸淺源強震而言,其震前、同震和震后的直接表現形式就是地表的變化,對斷裂帶及周邊區域地殼形變進行長期監測是地震預警的一項重要工作(張培震等,2003; He et al.,2013).當前研究中震后變形的解釋機制主要有震后余滑(Freed,2007),震后孔隙彈性回跳(He and Peltzer,2010)和殼下介質黏彈性松弛(Hu et al.,2016),這些不同的活動機制表現出的特征反映了地球深部物質的不同屬性(Lasserre et al.,2005).探討巖石圈流變學對震后蠕變隨時間變化的影響也是研究熱點之一(Montési,2004).一般而言,震后較長的周期內,余滑和黏彈性松弛均是必要的物理過程,但在當前的研究中,二者對地表形變貢獻的比例是很難精確區分的(溫揚茂等,2014).在震后較短的時間內,這些顯著不同的物理過程可產生相似的地表變形特征,因此構建一個與震后形變觀測場整體特征自洽的耦合模型在物理意義上是更合理的(Rousset et al.,2012).

盲逆斷層地震較難確定發震斷層的位置及幾何特征,利用GNSS等測量手段不易獲取全局震后變形特征,Rousset等(2012)利用GPS數據研究1999年臺灣集集地震震后變形中發現,觀測數據的空間位置會顯著影響對震后變形機制的約束.在盲逆型地震研究中,獲得大面積測量反映遠-近場整體變形特征的InSAR觀測資料尤其重要.Hong和Liu(2021)在2015年尼泊爾地震(盲逆斷層型)震后黏彈性模型研究中,發現擬合遠場GPS觀測的黏彈模型無法擬合近場變形特征:理論震后變形與近場InSAR及GPS觀測值相反,通過對比分析余滑機制可以更好解釋尼泊爾地震的震后變形特征,并表明余滑主要分布在同震斷層深部,在斷層淺部仍然存在再次發生強震的風險.本文的純運動學余滑模型及應力驅動余滑—黏彈耦合模型均表明伊朗地震震后余滑主要分布在同震滑動區域的上方斷層淺部,推測斷層上部再次發生較大地震的可能性較小.本次地震發生在結晶基地內部,為探討是否對上部沉積蓋層存在明顯的應力觸發作用,計算得到地下10 km處(沉積蓋層與結晶基底分隔界面)(Barnhart et al.,2018)的同震及震后產生的庫侖應力變化,接收斷層參數設置中傾角設為5°,其他與發震斷層一致,如圖12所示.同震庫侖應力主要以應力下降為主,庫侖應力范圍為(-0.26,0.06)MPa,在震后余滑及黏彈松弛的作用下,震后一年及一年半的庫侖應力改變范圍分別為(-0.062,0.024)MPa和(-0.064,0.026)MPa,雖然同震及震后庫侖應力增加的最大值都超過了地震觸發閾值0.1 bar(0.01 MPa),但量級都較小,而且地震一年之后庫侖應力的改變量十分微弱,結合震后余滑主要發生在發震斷層淺部,會釋放部分應力,綜合推測此次地震對上部沉積蓋層的地震觸發影響不顯著.

圖12 地下10 km處同震及震后庫侖應力變化

前人研究發現黏彈性松弛在伊朗地震前期震后變形中不占主要貢獻(Liu and Xu,2019),本文在震后耦合模型試算中,同樣進行了單獨物理機制的模型擬合,發現如果不引入余滑模型,僅考慮分層地殼黏彈模型的松弛效應,當黏滯系數降至～1×1018Pa·s時LOS向抬升才能達到同期觀測值～10 cm的量級,但模擬計算的整體變形范圍遠較觀測值大,而且在北東方向出現顯著的LOS向沉降區域,與InSAR實際觀測不符合.另外,一般內陸板塊如果沒有特殊地質構造如火山或熱地幔上涌等,其黏滯系數不會這么低.本文所觀測的震后形變特征主要體現在震源區的LOS向抬升,可以推斷伊朗主震一年半內同震應力驅動余滑機制在震后變形中起主導作用,單獨用下地殼、上地幔黏彈性松弛無法完全解釋震后變形現象.半解析模型推測震源區巖石圈黏滯系數的量級下限約為1019Pa·s,略大于Zhao等(2023)獲得的下地殼黏度下限 3×1018Pa·s,其原因可能與同震破裂模型及黏彈地殼模型有關,對余滑是引起前期震后變形主要因素的觀點是一致的.此外,前人研究中均指出在發震斷層淺部余滑現象較為顯著,但也存在一定的差異(Barnhart et al.,2018; Feng et al.,2018; Liu and Xu,2019; Zhao et al.,2023),包括發生余滑斷層幾何形態,余滑主要分布區域與同震滑動的位置關系.為探討斷層幾何模型對余滑影響,設置同震斷層模型在淺部結晶基地處傾角為0°,下部傾角保持不變,運動學余滑模型和應力驅動-黏彈耦合模型的反演和計算方法同上,發現總體特征與原模型較為一致.余滑主要發生在同震破裂區域的上方,主要表現為逆沖,水平斷層面和斷層深部余滑現象較弱.

4 結論

伊朗地震發生在扎格羅斯山前褶皺帶盲逆沖斷層內,破裂模型顯示發震斷層東傾、走向335°,表現為低傾角(10°)逆沖略帶右旋滑動,斷層面滑動分布主要集中深度14～20 km處,最大滑動量約7 m,較大滑動區域位于破裂點南側,震源動態過程顯示由初始破裂點往兩側擴展,大部分滑動約15 s完成,破裂沒有到達地表處.3D同震位移場顯示地震上盤垂向隆升最大值約65 cm,上盤東北區域沉降值最大值約50 cm,位移場整體西向運動,兩個中心區最大值約55 cm,理論計算位移場整體南向運動最大值約65 cm.地震上盤即扎格羅斯山前緣與阿拉伯板塊相互擠壓致使地殼縮短,深部基底滑脫是發生地震的主要成因,該地震促進了扎格羅斯山逆沖褶皺帶的生長和西南向擴展.

震后約一年半內的InSAR時序分析顯示震源區震后調整顯著,地震上盤出現明顯的LOS向抬升,形變場特征與同震相似,剖線AA′最大LOS向抬升速率最大值約8 cm·a-1,震后變形促進了褶皺的隆升,運動學余滑模型顯示震后半年釋放能量相當于MW6.7地震.通過半解析方法構建應力驅動震后余滑及黏彈松弛耦合模型,表明震后一年半內斷層面同震應力驅動余滑主要分布在淺部,在震后變形中起到主導作用,扎格羅斯褶皺帶西北段巖石圈下地殼的黏滯系數量級下限約為1019Pa·s,考慮到震源區位于阿拉伯板塊與歐亞板塊的匯聚部位,如要更準確分析巖石圈的流變特征,需要利用震后更長時間的InSAR形變時間序列,構建三維橫向不均一黏彈性模型進一步分析.

致謝感謝審稿專家提出的寶貴建議.本文遠場波形資料由IRIS提供,SENTINEL-1影像資料由歐空局(ESA)提供,本文圖件由GMT軟件生成,在此一并致謝.