基于InSAR地表形變約束的瑪多MS7.4地震孕育發(fā)生機(jī)理數(shù)值模擬研究

祝愛玉, 王永哲*, 李永華, 張東寧

1 中國(guó)地震局地球物理研究所, 北京 100081 2 北京白家疃地球科學(xué)國(guó)家野外科學(xué)觀測(cè)研究站, 北京 100095 3 中國(guó)地震局震源物理重點(diǎn)實(shí)驗(yàn)室, 北京 100095

0 引言

根據(jù)中國(guó)地震臺(tái)網(wǎng)中心測(cè)定,2021年5月22日2時(shí)4分在青海省果洛州瑪多縣發(fā)生7.4級(jí)地震,震中位于北緯34.59°,東經(jīng)98.34°,震源深度17 km.震后的地質(zhì)科考結(jié)果表明,瑪多MS7.4地震發(fā)震構(gòu)造西段發(fā)現(xiàn)大型地表破裂帶.該破裂帶總體走向NW230°—270°,由兩條NWW向的剪切破裂左階羽列組成.破裂帶由一系列的擠壓鼓包和張裂隙相間排列,張裂隙從十幾厘米到2～3 m不等,裂縫通常呈右階雁列狀排列,并顯示明顯垂直位錯(cuò),垂直位錯(cuò)可達(dá)約30 cm,并可清晰看到水平擦痕,顯示為具有左旋走滑兼正斷位移的同震斷層.

大地測(cè)量是獲取地震同震地表形變的重要手段,并為地震發(fā)生機(jī)理的研究提供重要基礎(chǔ)數(shù)據(jù)(袁霜等, 2020; 方進(jìn)等, 2019; 申文豪等, 2019; 陳威等, 2018; 季靈運(yùn)等, 2017; Xu, 2017; Song et al., 2019; Jónsson et al., 2002; Funning et al., 2005).同震地表形變觀測(cè)中,以全球定位系統(tǒng)(GPS)技術(shù)和合成孔徑雷達(dá)干涉測(cè)量(InSAR)技術(shù)應(yīng)用最廣泛.GPS觀測(cè)雖然得到的形變數(shù)據(jù)精度較高,但其空間分辨率取決于觀測(cè)站的分布和數(shù)量,瑪多地震所在區(qū)域的自然環(huán)境惡劣,用于進(jìn)行地表形變觀測(cè)的GPS觀測(cè)站非常少(朱亞戈等, 2021),難以準(zhǔn)確反映本次地震的地表詳細(xì)形變特征.InSAR技術(shù)由于可以獲得米級(jí)乃至分米級(jí)空間分辨率且覆蓋幾十公里范圍內(nèi)的形變信息,使得該技術(shù)在1993年首次應(yīng)用于強(qiáng)震地表形變觀測(cè)(Massonnet et al., 1993)以來(lái)得到廣泛應(yīng)用(Wang et al., 2019; Wen et al., 2021; Yang et al., 2021; Xu et al., 2016).該技術(shù)還具有全天時(shí)、全天侯、受天氣影響小等特點(diǎn),可以在地表植被覆蓋稀疏的環(huán)境條件下保持高相干性,獲取雷達(dá)視線方向1～2 cm精度的形變 (孫建寶等,2007),盡管觀測(cè)精度不如GPS技術(shù),但對(duì)于強(qiáng)震引起的地表形變觀測(cè)已經(jīng)足夠.此次地震區(qū)域地表植被稀疏,非常適合利用InSAR技術(shù)獲取此次地震的同震地表形變場(chǎng).采用InSAR觀測(cè)瑪多地震的同震形變場(chǎng),不僅能深入的研究地震造成的地表影響,同時(shí)能為數(shù)值模擬提供有力的約束.

瑪多地震是2008年汶川地震之后中國(guó)發(fā)生的震級(jí)最高的一次強(qiáng)震,該地震發(fā)生在巴顏喀拉地塊內(nèi)部的昆侖山口—江錯(cuò)斷裂上(王未來(lái)等, 2021).其北部和南部的東昆侖斷裂帶和甘孜—玉樹斷裂帶,為青藏高原東北緣的巨型走滑斷裂帶(鄧起東等, 2002; 劉光勛, 1996; Fu and Awata, 2005; Kirby et al., 2007; Karplus et al., 2011; 李陳俠等, 2011; 張培震等, 2008).根據(jù)歷史地震記錄,1900年以來(lái),瑪多地震震中周邊300 km范圍內(nèi)共發(fā)生7級(jí)以上地震5次,依次為1937年托索湖7.5級(jí)地震(李陳俠等, 2006)、1947年達(dá)日7.7級(jí)地震(劉雷等, 2021; 梁明劍等, 2020)、1963年阿蘭湖7.0級(jí)地震、1990年共和7.0級(jí)地震(李旭和陳運(yùn)泰, 1996)和2010年玉樹7.3級(jí)地震(張勇等, 2010; 孫鑫喆, 2012).強(qiáng)震導(dǎo)致的庫(kù)侖應(yīng)力變化和后續(xù)地震的發(fā)生在時(shí)間、空間上具有一定的相關(guān)性,斷層帶上庫(kù)侖應(yīng)力的增加對(duì)后續(xù)中強(qiáng)地震的發(fā)生具有顯著促進(jìn)作用(Toda et al., 2011a; Xiong et al., 2010; 石耀霖和曹建玲, 2010; Harris, 1998).為此開展歷史強(qiáng)震對(duì)瑪多地震的影響研究,對(duì)于分析瑪多地震的孕育機(jī)理具有重要的意義.

地震精定位時(shí)空序列顯示瑪多地震呈現(xiàn)出典型的主震-余震型特征,余震序列總長(zhǎng)度約170 km.主震位于余震區(qū)中間,呈雙側(cè)破裂特征 (王未來(lái)等, 2021).強(qiáng)震的同震位錯(cuò)會(huì)導(dǎo)致鄰近區(qū)域應(yīng)力狀態(tài)的改變(Toda et al., 2008; 單斌等, 2009),庫(kù)侖應(yīng)力變化大于0.01 MPa就會(huì)觸發(fā)地震(Toda et al., 2011b; Ma et al., 2005), 為此,通過(guò)研究瑪多地震序列導(dǎo)致的應(yīng)力場(chǎng)及庫(kù)侖應(yīng)力場(chǎng)變化,可為確定本次地震對(duì)其周邊斷裂帶的影響,進(jìn)而預(yù)判余震趨勢(shì)具有重要的意義(Li et al., 2021; Liu et al., 2014; Stein and Okal, 2011; Luo and Liu, 2010; 嚴(yán)珍珍等, 2009; 朱守彪等, 2008; Rundle et al., 2006)(圖1).

圖1 瑪多地震周邊構(gòu)造Fig.1 Geological structure around Madoi earthquake

本文采用合成孔徑雷達(dá)干涉測(cè)量(InSAR)技術(shù)得到地表形變場(chǎng)及地震破裂面分段結(jié)構(gòu),并建立能反映瑪多地震分段破裂面的黏彈性有限元模型,以InSAR觀測(cè)的地表形變、該地區(qū)深部結(jié)構(gòu)、巖石圈流變特性、GPS觀測(cè)給出的塊體運(yùn)動(dòng)為約束,模擬計(jì)算了瑪多地震引起的三維同震位移及同震應(yīng)力變化,討論了瑪多地震的發(fā)生對(duì)區(qū)域應(yīng)力場(chǎng)及周邊斷層的影響; 以瑪多破裂面為接收斷層模擬了周邊5次強(qiáng)震的同震庫(kù)侖應(yīng)力變化及震后松弛,討論了周邊強(qiáng)震對(duì)瑪多地震的觸發(fā)作用.

1 InSAR同震地表形變獲取

地震發(fā)生后,為全面了解地震的地表形變特征,我們收集了歐空局Sentinel-1A/1B兩顆衛(wèi)星升、降軌的合成孔徑雷達(dá)影像,組成兩個(gè)干涉對(duì).這兩幅干涉對(duì)的時(shí)間跨度均為6 d,即震前觀測(cè)時(shí)間為2021年5月20日和震后2021年5月26日.美國(guó)Scripps海洋研究所Sandwell教授團(tuán)隊(duì)開發(fā)的開源GMTSAR軟件(Sandwell et al., 2011 )作為本文干涉處理的軟件,采用的處理策略為二軌差分方法,利用水平分辨率3弧秒的SRTM地表高程數(shù)據(jù)(Farr et al., 2007)作為差分基礎(chǔ),從而去除原始干涉圖中的地形相位貢獻(xiàn).為保證獲取的地表形變的可靠性,在使用SNAPHU算法進(jìn)行差分相位解纏時(shí),我們只對(duì)相干性高于0.1的區(qū)域進(jìn)行了解纏,而忽略了低于0.1的區(qū)域形變信息.最終,將解纏后的雷達(dá)坐標(biāo)系下的結(jié)果地理編碼到WGS-84坐標(biāo)系下,得到此次地震的升、降軌雷達(dá)視線方向同震形變場(chǎng).結(jié)果顯示,此次地震引起的升軌方向最大抬升量和沉降量分別為: 1.398 m和1.482 m,而降軌則為: 1.291 m和1.506 m.為提高建模時(shí)的計(jì)算效率,我們對(duì)升、降軌同震形變場(chǎng)進(jìn)行了重采樣.采樣的方法是每3 km取一個(gè)觀測(cè)數(shù)據(jù)點(diǎn),同時(shí),不考慮相干性小于0.3的觀測(cè)數(shù)據(jù)點(diǎn).用于有限元模型的地表形變約束數(shù)據(jù),如圖2所示.

圖2 瑪多地震地表InSAR同震形變觀測(cè)(a) 升軌方向; (b) 降軌方向.Fig.2 The observed surface deformation of Madoi earthquake by InSAR(a) Ascending direction; (b) Descending direction.

為確定斷層所在位置以及對(duì)地震發(fā)生機(jī)制進(jìn)行定性分析(楊亞夫等, 2016),本文利用升、降軌雷達(dá)視線方向形變測(cè)量結(jié)果,根據(jù)SAR圖像的成像幾何獲得了東西方向和垂直方向的地表形變結(jié)果.具體的,利用式(1)所表征的地表三維形變與InSAR雷達(dá)視線向形變的關(guān)系,考慮南北方向形變對(duì)于雷達(dá)視線方向貢獻(xiàn)最小,InSAR技術(shù)對(duì)南北向形變的識(shí)別也最不敏感(王永哲等, 2012).因此,在忽略雷達(dá)視線方向中的南北方向形變分量情況下,利用升、降軌視線方向形變解算出東西方向和垂直方向的同震形變場(chǎng)(圖3).

圖3 瑪多地震的地表東西向(a)和垂直向(b)的InSAR同震形變場(chǎng)Fig.3 The coseismic surface deformation in the (a) East-West and (b) Vertical directions of the Madoi earthquake observed by InSAR

從圖中可以看出,地表的向西運(yùn)動(dòng)和向東運(yùn)動(dòng)存在明顯的分界線,可推斷此次地震破裂出地表,形成了較長(zhǎng)的破裂跡線; 破裂跡線兩側(cè)存在明顯的東西反向運(yùn)動(dòng),這與左旋走滑的機(jī)制非常吻合.另外,垂直向運(yùn)動(dòng)相對(duì)于東西向運(yùn)動(dòng)來(lái)講量級(jí)較小,可以推斷此次地震伴有少量的兩側(cè)相對(duì)升降運(yùn)動(dòng).根據(jù)東西向地表的錯(cuò)動(dòng),勾畫了斷層在地表的破裂跡線,將整個(gè)斷層共分成6段用于數(shù)值模擬的分析.考慮到東西向和垂直向形變的誤差,本文在建模時(shí),未將東西向和垂直向形變作為約束進(jìn)行擬合.

(1)

其中,dlos為雷達(dá)視線方向形變,de、du分別為東西向和垂直向形變,θ、α分別為雷達(dá)入射角和飛行方向角.

2 有限元數(shù)值模擬

2.1 數(shù)值模型

InSAR形變觀測(cè)給出了瑪多地震的破裂面分為6段(圖3),本文在建立幾何模型時(shí)充分考慮瑪多地震破裂面分段,建立了能反映分段破裂面及周邊斷層(鄧起東等, 2002)的三維黏彈性有限元模型,模型范圍為95°E—102.5°E,31.5°N—37°N,深度為100 km(圖4).研究區(qū)包含2條大型的走滑斷裂帶,分別為東昆侖斷裂帶及甘孜—玉樹斷裂帶，瑪多地震的發(fā)震構(gòu)造即昆侖山口—江錯(cuò)斷裂,以及一些較小的構(gòu)造,例如: 瑪多—甘德斷裂、達(dá)日斷裂、鄂拉山斷裂、日月山斷裂、共和盆地的南緣和北緣斷裂等.模型為XYZ坐標(biāo)系,X軸朝東、Y軸朝北、Z軸朝上.

圖4 三維有限元模型(a) 幾何模型; (b) 彈性模量; (c) 黏滯系數(shù).Fig.4 Three dimensional finite element model(a) Geometric model; (b) Elastic modulus; (c) Viscosity coefficient.

以Wang和Shen(2020)給出的中國(guó)大陸及其鄰區(qū)地殼GPS速度場(chǎng)(相對(duì)于歐亞大陸)作為模型邊界條件和檢驗(yàn)?zāi)M結(jié)果的依據(jù).采用GPS的觀測(cè)速度值插值到模型邊界給出模型四周邊界的水平速度值.該速度值乘以計(jì)算時(shí)間步長(zhǎng)作為有限元模型的水平位移約束條件,且假定模型四周的水平約束不隨深度變化,即假設(shè)從地表到100 km深度保持一致,垂直方向位移可以保持自由.三維模型底面所有節(jié)點(diǎn)垂直方向(z方向)上位移約束為0,水平向可以自由運(yùn)動(dòng).模型上表面為自由表面,即法向應(yīng)力和剪應(yīng)力均為零.

采用多物理場(chǎng)耦合有限元軟件COMSOL進(jìn)行模擬計(jì)算,時(shí)間步長(zhǎng)500年(依據(jù)模型給定彈性和流變材料參數(shù)計(jì)算出的數(shù)百年至數(shù)千年的黏彈性松弛時(shí)間),根據(jù)時(shí)間步長(zhǎng)逐步施加邊界位移約束.經(jīng)過(guò)了40個(gè)時(shí)間步,即20000年的加載作用,以最后得到的研究區(qū)穩(wěn)定地殼構(gòu)造應(yīng)力場(chǎng)作為研究地區(qū)的背景構(gòu)造應(yīng)力場(chǎng).通過(guò)模擬計(jì)算的地表形變和現(xiàn)今GPS觀測(cè)形變對(duì)比分析(圖5),可以看出計(jì)算值和測(cè)量值從方向和大小來(lái)看都比較相近,說(shuō)明了計(jì)算模型的可靠性.

圖5 地殼形變速度的計(jì)算值和測(cè)量值的對(duì)比(圖中紅色箭頭為模擬結(jié)果; 藍(lán)色箭頭為Wang和Shen(2020)發(fā)表的相對(duì)歐亞板塊的GPS觀測(cè)速度結(jié)果)Fig.5 Comparison of calculated and observed horizontal deformation rate (The red arrows are the simulated results;The blue arrows are GPS observations (Wang and Shen, 2020))

2.2 同震位移模擬結(jié)果

基于黏彈性三維有限元模型,構(gòu)建符合實(shí)際的背景應(yīng)力場(chǎng).通過(guò)降低瑪多地震分段破裂區(qū)域的彈性模量,采用試錯(cuò)法確定破裂段1,2,3,4的彈性模量降低為7.8 GPa,破裂段5,6的彈性模量降低為5.8 GPa,以地表的InSAR形變觀測(cè)為約束,從而確定瑪多地震引起的位移場(chǎng)和應(yīng)力場(chǎng)的變化.由于本文是將瑪多地震的主震和震后大約8 d的時(shí)間內(nèi)余震作為一個(gè)整體進(jìn)行求解,因此本文的應(yīng)力場(chǎng)變化能反映主震和大部分余震形成的應(yīng)力釋放.

為了將模擬的地表形變與Sentinel-1升、降軌InSAR同震地表形變比較,采用公式(1)將三維地表形變投影到雷達(dá)視線方向(升、降軌)上,并與實(shí)際觀測(cè)的形變(圖2a,6b)相減,得到觀測(cè)值與模擬值之間的殘差如圖6a(升軌)和6c(降軌)所示.根據(jù)殘差數(shù)據(jù)分布的直方圖(圖6b,6d)可以看出升軌的殘差超過(guò)80%的值小于0.1 m,95%的殘差不超過(guò)0.2 m;降軌的殘差分布85%的值不超過(guò)0.2 m.在進(jìn)行模型構(gòu)建時(shí),本文將整個(gè)斷裂簡(jiǎn)單地分成了6段.然而,此次瑪多地震實(shí)際的地表破裂比較復(fù)雜,并非簡(jiǎn)單的直線分段,這會(huì)造成本來(lái)位于上盤的觀測(cè)點(diǎn)在模擬時(shí)按照下盤來(lái)模擬,下盤的觀測(cè)點(diǎn)按照上盤來(lái)模擬,這必然會(huì)造成模擬值與觀測(cè)值存在較大的殘差.從殘差圖中可以看出,盡管有少數(shù)幾個(gè)觀測(cè)點(diǎn)的殘差值較大,但這些殘差較大的點(diǎn)對(duì)于此次地震的整個(gè)形變場(chǎng)模式影響不大,模擬的升、降軌形變場(chǎng)與實(shí)際觀測(cè)值具有較好地一致性.根據(jù)以上觀測(cè)值和模擬值的對(duì)比分析,可以看出模擬得到的地表形變分布特征和InSAR觀測(cè)的地表形變基本吻合.

圖6 模擬值和InSAR觀測(cè)的地表形變殘差(a) 升軌方向觀測(cè)值與模擬值的殘差; (b) 升軌方向殘差直方圖; (c) 降軌方向觀測(cè)值與模擬值的殘差; (d) 降軌方向殘差直方圖.Fig.6 Surface deformation residuals of simulated results and InSAR observations(a) The residuals of the observed and simulated values in the ascending direction; (b) The histogram of the residual error in the ascending direction; (c) The residual error between the observed value and the simulated value in the descending direction; (d) Histogram of residual error in the descending direction.

本文模擬得出的三維同震位移(圖7)結(jié)果顯示,東西向分量最大,其最大和最小值分別為1.1 m和-1.38 m;其次是南北向分量,其最大和最小值分別為0.7 m和-0.74 m; 垂向分量相對(duì)較小,最大為0.3 m.瑪多地震總的同震位移(圖7d)呈現(xiàn)出“四象限”分布,是走滑斷層錯(cuò)動(dòng)的典型圖像.同震位移結(jié)果還顯示,此次地震發(fā)震斷層的北側(cè),位移方向?yàn)楸蔽飨? 在發(fā)震斷層的南側(cè),位移方向?yàn)闁|南向; 發(fā)震斷層的南側(cè)位移要略大于北側(cè)位移,比較大的錯(cuò)動(dòng)主要集中在東側(cè)和西側(cè),中間部位的錯(cuò)動(dòng)相對(duì)較小.

其次，觀察非對(duì)角線的子群關(guān)系：子群1到2，1到3， 3到5，8到7等子群聯(lián)系緊密。結(jié)合上文子群內(nèi)部密度及往返聯(lián)系分析，從子群2與子群1密度達(dá)到0.640，表明在子群1與2的游客行為中，游客經(jīng)成都、重慶都市旅游景區(qū)向著西部更遠(yuǎn)景區(qū)頻繁流動(dòng)，如至武隆、至九寨-黃龍，若悲觀估計(jì)，三峽旅游可能正逐步發(fā)展成為成渝西部地區(qū)旅游的跳板，而非旅途終點(diǎn)站，游客的最終旅游目的地愈加朝著西部、多種形態(tài)的內(nèi)陸旅游目的地延伸。

圖7 三維同震位移模擬結(jié)果(a) 東西向分量; (b) 南北向分量; (c) 垂直向分量; (d) 總位移.Fig.7 The simulation results of 3D coseismic displacement (a) East-West component; (b) North-South component; (c) Vertical component; (d) Total displacement.

2.3 同震應(yīng)力場(chǎng)變化模擬結(jié)果

以地表同震位移為約束,本文給出了10 km深度的同震應(yīng)力變化和庫(kù)侖應(yīng)力變化(圖8).庫(kù)侖應(yīng)力變化ΔCFS為

ΔCFS=Δτ+μΔσn,

(2)

其中Δτ為斷層面上的剪切力變化; Δσn為斷層面上的正應(yīng)力變化,壓應(yīng)力為負(fù)數(shù);μ為斷層面上的有效摩擦系數(shù),一般取值范圍為0.2～0.8,本文取值為0.4(Stein et al., 1992; King et al., 1994; Freed, 2005; 石耀霖和曹建玲, 2010).有效摩擦系數(shù)的取值會(huì)在一定程度上改變庫(kù)侖應(yīng)力的大小,但不會(huì)影響庫(kù)侖應(yīng)力正負(fù)的改變.

圖8a顯示在整個(gè)破裂面上是以壓應(yīng)力為主,在破裂面的西南向和北東向有拉張分量,并呈現(xiàn)不規(guī)則的四象限分布.圖8b顯示破裂面的北側(cè)和南側(cè)剪應(yīng)力增強(qiáng),而東側(cè)和西側(cè)剪應(yīng)力減低,因此地震呈現(xiàn)左旋走滑趨勢(shì).以上正應(yīng)力和剪應(yīng)力的分布特征都揭示了瑪多地震為典型的走滑型地震的特征.由于本文考慮了區(qū)域彈性模量等參數(shù)的非均勻性,以及區(qū)域斷層分布的影響,而應(yīng)力場(chǎng)分布受區(qū)域構(gòu)造和深部結(jié)構(gòu)的影響,因此本文結(jié)果更能反映區(qū)域應(yīng)力場(chǎng)分布特征.

根據(jù)瑪多地震的震源機(jī)制解(郭祥云,http: ∥www.cea-igp.ac.cn/kydt/278249.Html; 萬(wàn)永革等,https:∥mp.weixin.qq.com/s/60hFBoaBXa-c_Q9BVCNYHw),其破裂面的整體走向?yàn)?01°,傾向?yàn)?7°,滑動(dòng)角為-7°,采用公式(2)計(jì)算得到10 km處的同震庫(kù)侖應(yīng)力變化(圖8c,8d,8e).圖8d給出了震后3個(gè)小時(shí)的地震精定位結(jié)果,可以看出余震主要分布在西部和東部,而中間部位有一段地震空區(qū).圖8e給出了震后8天的余震精定位,隨著時(shí)間的推移,余震數(shù)量越來(lái)越多,可以看出余震向西部和東部繼續(xù)擴(kuò)展,東部的余震逐漸轉(zhuǎn)向?yàn)楸睎|向,中間部位的地震空區(qū)也在逐漸變小.對(duì)比庫(kù)侖應(yīng)力變化可以看出,余震都分布在庫(kù)侖應(yīng)力影區(qū),說(shuō)明瑪多地震已經(jīng)充分釋放了應(yīng)變，其和庫(kù)侖應(yīng)力影區(qū)有很好的一致性,同時(shí)中間的地震空區(qū)和庫(kù)侖應(yīng)力增強(qiáng)區(qū)域也有很好的對(duì)應(yīng)關(guān)系.

圖8 同震應(yīng)力變化(a) 正應(yīng)力變化; (b) 剪應(yīng)力變化; (c) 庫(kù)侖應(yīng)力變化; (d) 庫(kù)侖應(yīng)力變化和震后3個(gè)小時(shí)的余震分布; (e) 庫(kù)侖應(yīng)力變化和震后8天內(nèi)的余震分布.Fig.8 Coseismic stress variation(a) Normal stress variation; (b) Shear stress variation; (c) Coulomb stress variation; (d) Coulomb stress variation and aftershock distribution in 3 hours after the Madoi earthquake; (e) Coulomb stress variation and aftershock distribution within 8 days after the Madoi earthquake.

從庫(kù)侖應(yīng)力變化(圖8c)看出,瑪多地震發(fā)震斷層兩側(cè)的震后庫(kù)侖應(yīng)力變化為負(fù)數(shù).說(shuō)明瑪多地震的發(fā)生,使得發(fā)震斷層南北兩側(cè)區(qū)域應(yīng)力都得到了充分的釋放.而庫(kù)侖應(yīng)力的增加區(qū)域是發(fā)震斷層的西南和東北端,推測(cè)未來(lái)的余震的發(fā)生主要向西南和東北方向發(fā)展.

3 討論

3.1 區(qū)域構(gòu)造應(yīng)力對(duì)瑪多地震的影響

根據(jù)三維應(yīng)力和應(yīng)變模擬計(jì)算結(jié)果給出了未考慮歷史強(qiáng)烈地震活動(dòng)時(shí)研究區(qū)最大剪應(yīng)力分布模擬結(jié)果(圖9),其中最大剪應(yīng)力為(σ1-σ3)/2,σ1為最大主壓應(yīng)力,σ3為最小主壓應(yīng).圖9中的黑色五角星表示自1900年以來(lái)發(fā)生在研究區(qū)域M7級(jí)以上地震(表1).首先從不同深度的最大剪應(yīng)力分布可以看出,10 km處的剪應(yīng)力最大.根據(jù)地震精定位的研究結(jié)果(王未來(lái)等,2021),瑪多地震的主震深度為10 km左右,這與我們的計(jì)算結(jié)果一致.其次,區(qū)域剪應(yīng)力高值區(qū)的分布,主要集中在東昆侖斷裂帶西部、昆侖山口—江錯(cuò)斷裂、甘孜—玉樹斷裂的西部,這些剪應(yīng)力高值區(qū)和歷史的強(qiáng)震有很好的對(duì)應(yīng)關(guān)系,尤其是1937年托索湖地震、1963年阿蘭湖地震、2010年玉樹地震和2021年瑪多地震等.同時(shí),高剪應(yīng)力區(qū)域和彈性模量低值區(qū)域相對(duì)應(yīng)(圖4b),即在周邊塊體的擠壓作用下,較低的彈性模量更容易累積應(yīng)變,因而使得斷層更加容易滑動(dòng).此次瑪多7.4級(jí)地震明顯處于應(yīng)力增強(qiáng)區(qū)域,因此地震的發(fā)生和周邊塊體的擠壓作用及地處低彈性模量區(qū)域密不可分.

3.2 周邊歷史強(qiáng)震對(duì)2021年瑪多MS7.4地震的觸發(fā)作用

本文根據(jù)歷史地震目錄統(tǒng)計(jì)了瑪多地震震中周邊300 km范圍內(nèi),1900年以來(lái)發(fā)生7級(jí)以上地震5次(表1),其中1990年共和地震是逆沖型地震,其他4次強(qiáng)震發(fā)生在巴顏喀拉地塊的邊界或者內(nèi)部,都是走滑型地震,而1937年的托索湖地震破裂尺度達(dá)到240 km,影響范圍極廣.

表1 區(qū)域1900年以來(lái)研究區(qū)7級(jí)以上地震目錄和破裂方向數(shù)據(jù)Table 1 Catalogue and fracture direction of earthquakes with magnitude 7 and above in the study area since 1900

圖9 不同深度的最大剪應(yīng)力(a) 5 km處; (b) 10 km處;(c) 15 km處; (d) 20 km處.Fig.9 Maximum shear stress at different depths(a) 5 km; (b) 10 km; (c) 15 km; (d) 20 km.

本文計(jì)算了以上5次7級(jí)以上地震事件投影到瑪多地震破裂面的同震庫(kù)侖應(yīng)力變化(圖10).根據(jù)同震庫(kù)侖應(yīng)力的影響可以看出(圖10a),1937年托索湖7.5級(jí)地震、1963年阿蘭湖7.0級(jí)地震和1990年共和7.0級(jí)地震在瑪多地震主破裂面上的投影都為負(fù)數(shù),說(shuō)明其對(duì)瑪多地震破裂面有卸載的作用,使得瑪多地震的發(fā)生推遲.2010年玉樹7.3級(jí)地震的發(fā)生在瑪多地震主破裂面的庫(kù)侖應(yīng)力投影為正數(shù),說(shuō)明玉樹地震對(duì)瑪多地震具有促進(jìn)作用.

圖10 1900年以來(lái)區(qū)域5次強(qiáng)震投影在瑪多地震破裂面上的庫(kù)侖應(yīng)力變化(a) 同震庫(kù)侖應(yīng)力變化; (b) 考慮構(gòu)造應(yīng)力的庫(kù)侖應(yīng)力變化.Fig.10 Coulomb stress variation of 5 strong earthquakes projected on the fracture surface of Madoi earthquake(a) Coseismic Coulomb stress variation; (b) Coulomb stress variation considering tectonic stress.

其次,我們考慮了震后的黏彈松弛效應(yīng),計(jì)算了以上5次強(qiáng)震在瑪多地震主破裂面的庫(kù)侖應(yīng)力變化(如圖10b),可以看出隨著震后的松弛效應(yīng),瑪多主破裂面的庫(kù)侖應(yīng)力變化在塊體邊界的作用力下,應(yīng)變?cè)诓粩嗟睦鄯e,隨著區(qū)域的地震的發(fā)生,例如1963年的阿蘭湖地震、1990年的共和地震使得瑪多地震發(fā)震斷層的應(yīng)力得到了部分釋放,從而推遲了瑪多地震的發(fā)生,而2010年的玉樹地震使得瑪多地震發(fā)震斷層的應(yīng)力突然增加,從而加速了瑪多地震的發(fā)生.

3.3 瑪多地震序列對(duì)周邊斷層的影響

由于以上的同震庫(kù)侖應(yīng)力的計(jì)算是以瑪多地震發(fā)震斷層為接收斷層,因此只能反映和發(fā)震斷層相同走向和傾向的周邊斷層的庫(kù)侖應(yīng)力影響,而不能反映其他周邊斷層的影響.為了進(jìn)一步研究瑪多地震對(duì)周邊主要斷層的庫(kù)侖應(yīng)力影響,本文給出了以周邊主要斷層為接收斷層的庫(kù)侖應(yīng)力變化(圖11所示),可以看出瑪多地震序列的發(fā)生,使得南北向的應(yīng)力得到了充分的釋放,導(dǎo)致東昆侖斷裂帶的西部、瑪多—甘得斷裂帶的中部偏西部位、達(dá)日斷裂帶的西部以及甘孜—玉樹斷裂帶的東部區(qū)域的庫(kù)侖應(yīng)力變化為負(fù); 庫(kù)侖應(yīng)力變化明顯增強(qiáng)的區(qū)域?yàn)槔錾娇凇e(cuò)斷裂的北部、瑪多—甘德斷裂的西部及東部部分區(qū)域、東昆侖斷裂帶中部、共和盆地南緣斷裂和北緣斷裂,以上區(qū)域的庫(kù)侖應(yīng)力變化明顯大于0.01 MPa,說(shuō)明瑪多地震序列的發(fā)生對(duì)以上斷裂帶有很大影響,大大的增加了上述斷裂帶滑動(dòng)的可能; 甘孜—玉樹斷裂帶西部及日月山斷裂帶的庫(kù)侖應(yīng)力變化為正,說(shuō)明瑪多地震會(huì)加快以上斷裂帶滑動(dòng).

圖11 周邊主要斷裂帶的同震庫(kù)侖應(yīng)力變化Fig.11 Variation of coseismic Coulomb stress of main faults around the Madoi earthquake

4 結(jié)論

在瑪多地震發(fā)生的區(qū)域,Sentinel-1影像數(shù)據(jù)具有較好的覆蓋,本文利用較短時(shí)間基線的四景雷達(dá)影像數(shù)據(jù),通過(guò)干涉處理獲得了瑪多地震的升、降軌InSAR同震形變場(chǎng),得到了兩個(gè)雷達(dá)視線方向的同震形變場(chǎng).通過(guò)忽略南北向形變,解算了東西方向和垂直方向的同震形變場(chǎng).根據(jù)東西向地表的錯(cuò)動(dòng),勾畫了斷層在地表的跡線,將破裂面共分成6段用于數(shù)值模擬的分析.東西方向的地表形變表明,此次地震發(fā)生了較大的東西向運(yùn)動(dòng),破裂跡線南側(cè)向東、北側(cè)向西運(yùn)動(dòng),且兩側(cè)的相對(duì)運(yùn)動(dòng)明顯,垂直向形變量級(jí)較小,沿地表破裂跡線無(wú)明顯地表錯(cuò)動(dòng).由此可知,此次地震為左旋走滑為主的地震,這與地震學(xué)所得結(jié)果一致.

本文模擬計(jì)算得到了瑪多地震周邊區(qū)域三維的地表形變,通過(guò)與InSAR觀測(cè)結(jié)果進(jìn)行比較,證實(shí)了模擬結(jié)果的可靠性.在此基礎(chǔ)上,給出了三維地表同震形變場(chǎng),可以看出瑪多地震引起了較大的地表變形,其中東邊部分區(qū)域的地表形變較大,錯(cuò)動(dòng)量達(dá)到2.5 m,中間有部分區(qū)域形變相對(duì)較小,同時(shí)也可以看出瑪多地震為典型的左旋走滑型地震,模擬結(jié)果與震源機(jī)制及地質(zhì)考察的結(jié)果一致.

本文以Sentinel-1升、降軌InSAR同震地表形變等多個(gè)結(jié)果作為約束,給出了瑪多地震引起的應(yīng)力場(chǎng)變化,結(jié)果顯示瑪多地震的余震都分布在庫(kù)侖應(yīng)力為負(fù)的區(qū)域,說(shuō)明瑪多地震序列釋放了發(fā)震斷層上積累的應(yīng)力; 而發(fā)震斷層的西南和東北端庫(kù)侖應(yīng)力有所增加,推測(cè)未來(lái)余震的發(fā)生可能會(huì)朝西南和東北方向發(fā)展.

為了探討歷史地震對(duì)瑪多地震的影響,本文模擬了周邊5次7級(jí)以上強(qiáng)震對(duì)瑪多地震斷層面的影響.結(jié)果顯示周邊地塊的加載使得瑪多地震斷層面的應(yīng)變不斷的累積,2010年玉樹地震對(duì)瑪多地震的發(fā)生有很大的促進(jìn)作用.本文也給出了瑪多地震的發(fā)生對(duì)周邊斷層的影響,模擬結(jié)果顯示,瑪多地震造成了昆侖山口—江錯(cuò)斷裂的西部、東昆侖斷裂帶中部、共和盆地南緣斷裂、北緣斷裂等多個(gè)應(yīng)力強(qiáng)加載區(qū)(庫(kù)侖應(yīng)力變化明顯大于0.01 MPa),推測(cè)瑪多地震的發(fā)生將大大增加上述斷裂帶滑動(dòng)的可能; 瑪多地震同時(shí)造成了甘孜—玉樹斷裂帶西部、日月山段帶等多個(gè)弱加載區(qū)(庫(kù)侖應(yīng)力>0),由此推測(cè)瑪多地震可能會(huì)對(duì)上述斷裂帶有一定的促進(jìn)作用.

同時(shí),本文也存在一定的不足之處.本文采用軟弱夾層反映斷層的屬性,利用減低地震破裂區(qū)域的彈性模量的方式模擬地震引起的應(yīng)力場(chǎng)的變化.因此本文的模擬結(jié)果主要是在InSAR觀測(cè)形變的約束下,反映瑪多地震的孕育發(fā)震原因及瑪多地震的發(fā)生可能會(huì)對(duì)周邊斷層的影響,對(duì)于精細(xì)的模擬瑪多地震的滑動(dòng)特征,接下來(lái)我們需要考慮具體的地震破裂模型.

致謝本文采用的最新的波速和密度數(shù)據(jù)是中國(guó)地震局地球物理研究所石磊博士提供,黏滯系數(shù)參數(shù)是國(guó)土資源部孫玉軍博士提供.歐洲空間局(ESA)免費(fèi)發(fā)布Sentinel-1A/1B SAR數(shù)據(jù).在此表示感謝.同時(shí)感謝匿名審稿專家的建議.