孟思晨,孟秋,陳啟志,胡才博

(中國(guó)科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院 中國(guó)科學(xué)院計(jì)算地球動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室, 北京 100049) (2021年11月14日收稿; 2022年4月26日收修改稿)

2021年5月22日2時(shí)4分,青海果洛州瑪多縣(34.59°N,98.34°E)發(fā)生MS7.4地震(圖1),震源深度17 km(中國(guó)地震臺(tái)網(wǎng)正式測(cè)定)[1]。該地震極震區(qū)的烈度最大值為Ⅹ度,Ⅵ度及以上區(qū)域面積約為53 704 km2,宏觀震中和微觀震中分別為瑪多縣瑪查理鎮(zhèn)和黃河鄉(xiāng)[2]。從發(fā)震構(gòu)造上看,青?，敹嗟卣鸢l(fā)生在青藏高原中部的Ⅱ級(jí)地塊巴顏喀拉塊體北部,距其北邊界東昆侖斷裂帶約70 km。近年來,巴顏喀拉塊體周緣發(fā)生一系列強(qiáng)震,如2008年汶川8.0級(jí)地震、2011年玉樹7.1級(jí)地震、2013年蘆山7.0級(jí)地震、2017年九寨溝7.0級(jí)地震,而瑪多地震卻是唯一發(fā)生在巴顏喀拉塊體內(nèi)部的強(qiáng)震[3],因此受到科研人員的廣泛關(guān)注。人們從震源機(jī)制、震源破裂過程、發(fā)震斷裂、余震序列、同震滑動(dòng)分布、震區(qū)構(gòu)造應(yīng)力場(chǎng)、庫(kù)侖應(yīng)力變化和地震危險(xiǎn)性等方面對(duì)瑪多地震進(jìn)行了研究。

F1：達(dá)日斷裂帶;F2：瑪多—甘德斷裂帶;F3：東昆侖斷裂帶;F4：昆中斷裂帶;F5：鄂拉山斷裂帶,斷層數(shù)據(jù)引自詹艷等[4]。圖1 瑪多地震的構(gòu)造背景和地震活動(dòng)性Fig.1 The tectonic setting of the Maduo earthquake and the historical seismicity in research area

中國(guó)地震臺(tái)網(wǎng)、美國(guó)地質(zhì)調(diào)查局等機(jī)構(gòu)對(duì)此次地震的震源機(jī)制和震源破裂過程反演后認(rèn)為該地震為一次左旋走滑型地震事件,兼有正斷分量[4],為雙側(cè)破裂方式,震區(qū)構(gòu)造應(yīng)力場(chǎng)最大主應(yīng)力軸方向?yàn)镹E-SW,最小主應(yīng)力軸方向?yàn)镹NW-SSE[1]。李志才等[5]對(duì)GNSS數(shù)據(jù)進(jìn)行反演確定的地震矩、斷裂的幾何參數(shù)及同震滑動(dòng)分布,也揭示了發(fā)震斷裂左旋走滑特性,同時(shí)得到同震滑動(dòng)量的數(shù)值及分布。余震雙差重定位研究發(fā)現(xiàn),余震帶與昆侖山口—江錯(cuò)斷裂距離最近且兩者走向有一定的重合度,說明瑪多地震的發(fā)震斷裂為昆侖山口—江錯(cuò)斷裂[3,6],與野外地震考察結(jié)果相一致[7-10]。前人對(duì)于瑪多地震同震及震后變形研究可以分為大地測(cè)量數(shù)據(jù)反演及同震、震后變形及相關(guān)應(yīng)力變化模擬多個(gè)方面。研究者基于InSAR、GPS等大地測(cè)量數(shù)據(jù)與圖像,揭示出最大滑動(dòng)量為6.0 m[11-12],并確定了同震變形、靜態(tài)滑動(dòng)分布以及斷層破裂過程[13],量化了地震間應(yīng)變率、斷層幾何以及同震滑動(dòng)分布[14]。基于InSAR數(shù)據(jù)反演的斷層滑動(dòng)模型和PREM地球模型,研究者計(jì)算了此次地震的理論同震位移、大地水準(zhǔn)面、重力及應(yīng)變的變化,得到地殼變形的空間分布,為研究地震斷層及解釋大地測(cè)量數(shù)據(jù)提供了參考[15]。在同震、震后變形模擬及庫(kù)侖應(yīng)力計(jì)算方面,楊光遠(yuǎn)等[16]使用Coulomb3.3程序[17-18]計(jì)算瑪多地震引起的地震破裂面及周邊主要斷層上的庫(kù)侖應(yīng)力變化;談洪波等[19]和Li等[20]使用均勻半空間解析解PSGRN/PSCRM程序[21]對(duì)瑪多地震同震、震后變形、重力變化及發(fā)震斷層及周圍區(qū)域的庫(kù)侖應(yīng)力變化進(jìn)行了計(jì)算模擬。

前人使用解析法或半解析法對(duì)同震、震后及庫(kù)侖應(yīng)力變化進(jìn)行理論計(jì)算,模型比較簡(jiǎn)單,計(jì)算速度快,因此,在大地震發(fā)生不久,通常會(huì)有相關(guān)的同震、震后及庫(kù)侖應(yīng)力變化成果發(fā)布,對(duì)大地震周圍區(qū)域的地震危險(xiǎn)性做出第一時(shí)間的評(píng)估和判斷。發(fā)震斷層周圍介質(zhì)的非均勻性、幾何的復(fù)雜性,對(duì)大地震的同震變形影響很大,簡(jiǎn)單的解析解和半解析解很難考慮這些復(fù)雜因素,而數(shù)值模擬可以發(fā)揮獨(dú)特的優(yōu)勢(shì)。在數(shù)值模擬中,有限元方法可以根據(jù)實(shí)際地質(zhì)體的幾何、材料及邊界條件靈活地進(jìn)行網(wǎng)格劃分,同時(shí)進(jìn)行并行計(jì)算,實(shí)現(xiàn)對(duì)復(fù)雜地質(zhì)體高分辨率模擬。分裂節(jié)點(diǎn)法通過將位錯(cuò)加入到載荷向量而不改變剛度矩陣,可以簡(jiǎn)潔有效地將斷層位錯(cuò)引入有限元計(jì)算。為此,本文采用C語言編寫了分裂節(jié)點(diǎn)法的三維并行彈性有限元程序,用該程序計(jì)算了青?，敹嗟卣鸬耐鹱冃魏蛻?yīng)力變化,并對(duì)研究區(qū)域主要斷層地震危險(xiǎn)性進(jìn)行探討。

1 有限元方法

1.1 使用分裂節(jié)點(diǎn)法的有限元模型構(gòu)建

用有限元方法模擬地震變形時(shí),首先面臨的問題是如何將斷層的同震位錯(cuò)引入有限元數(shù)值計(jì)算中。分裂節(jié)點(diǎn)法可以簡(jiǎn)單有效地解決該問題[22]。如圖2所示,單元1、2 分別位于斷層的左右兩側(cè),U為位移,上標(biāo)表示單元號(hào),下標(biāo)表示單元內(nèi)的節(jié)點(diǎn)號(hào);ΔU表示位錯(cuò)量[22]。

圖2 分裂節(jié)點(diǎn)法的示意圖Fig.2 Schematic diagram of split node technique

則有

(1)

(2)

(3)

(4)

單元1的單元?jiǎng)偠染仃嚪匠虨?/p>

(5)

同樣,單元2的單元?jiǎng)偠染仃嚪匠虨?/p>

(6)

總體剛度矩陣方程為

(7)

由此可以看出,在分裂節(jié)點(diǎn)法中,斷層位錯(cuò)的引入體現(xiàn)在單元、總體的載荷向量上,而對(duì)于單元、總體的剛度矩陣是沒有影響的,大大簡(jiǎn)化了剛度矩陣和載荷向量的組裝,且對(duì)于計(jì)算的準(zhǔn)確性和穩(wěn)定性沒有影響[19],體現(xiàn)了該方法的簡(jiǎn)潔性和有效性。

1.2 有限元與理論解結(jié)果對(duì)比

為驗(yàn)證我們自主研發(fā)的采用分裂節(jié)點(diǎn)技術(shù)的三維并行彈性有限元程序的有效性,我們建立了一個(gè)三維走滑斷層地震模型(模型1)和一個(gè)三維逆斷層地震模型(模型2),分別利用有限元程序和基于半無限空間位錯(cuò)理論[23]的Coulomb 3.3[17-18]及edgrn/edcmp[21]計(jì)算2個(gè)斷層模型的地表同震位移場(chǎng),并對(duì)比兩者的結(jié)果。2個(gè)斷層模型的參數(shù)如表1所示。

表1 走滑和逆沖斷層地震模型的材料參數(shù)與斷層參數(shù)Table 1 The material parameters and fault parameters of strike-slip and reverse fault earthquakes

圖3分別展示了理論解和有限元計(jì)算得到的走滑斷層和逆斷層的同震地表位移的等值線圖和剖面結(jié)果圖。圖3(a)～3(f)左側(cè)為Coulomb3.3[17-18]解析解結(jié)果,右側(cè)為有限元模擬結(jié)果,表明有限元得到的走滑斷層地震和逆斷層地震引起的地表同震位移的分布形態(tài)與理論解有很好的一致性。圖3(g)～3(i)分別給出了有限元和理論解的過斷層中心點(diǎn)垂直斷層走向的剖面上走滑斷層地震的南北向同震位移、逆斷層地震的東西向位移和垂直向位移的對(duì)比圖,FEM表示有限元模擬結(jié)果,Theory1表示Coulomb3.3[17-18]計(jì)算結(jié)果,Theory2表示edgrn/edcmp[21]計(jì)算結(jié)果,表明有限元與理論解有很好的一致性。

圖3 走滑斷層和逆斷層地震引起的同震位移的有限元模擬結(jié)果與解析解對(duì)比Fig.3 Comparison of the analytical solutions with finite element models simulation results for the coseismic displacement components of the strike-slip fault and thrust fault

2 青?，敹嗟卣鸬耐鹱冃魏蛻?yīng)力變化數(shù)值模擬與比較

2.1 有限元模型的設(shè)置

采用前人反演得到的非均勻位錯(cuò)模型以及我們基于分裂節(jié)點(diǎn)技術(shù)[22]的有限元方法,本文計(jì)算了瑪多地震的同震位移和同震應(yīng)力場(chǎng)。有限元模型東西長(zhǎng)1 100 km,南北長(zhǎng)1 000 km,深度100 km,共有64個(gè)并行計(jì)算分區(qū),有限元網(wǎng)格數(shù)量為2 171 715,可以實(shí)現(xiàn)較高分辨率的有限元并行計(jì)算。本文采用中國(guó)地震局地質(zhì)研究所的非均勻位錯(cuò)有限斷層模型[15],如圖4所示。

圖4 瑪多地震主震斷層同震位錯(cuò)分布[15]Fig.4 Coseismic dislocation distribution of Maduo earthquake[15]

本文建立了2個(gè)有限元模型：均勻模型和分層模型,其材料參數(shù)如表2所示,楊氏模量和泊松比可由表中的密度、P波速度和S波速度計(jì)算得到。由于有限元模型的尺寸大大超過主震斷層,本文假定除模型上表面(地表)外,其余邊界上的法向位移為0。

表2 瑪多地震同震變形三維有限元分層模型的材料參數(shù)(自USGS)Table 2 Material parameters of the 3D finite element layered model of Maduo earthquake(USGS)

2.2 瑪多地震地表同震變形和同震應(yīng)力變化

根據(jù)上述三維彈性有限元模型,本文模擬得到的瑪多地震地表同震位移場(chǎng)見圖5。東西向位移u大于零的部分位于斷層南側(cè),小于零的部分主要集中于斷層北側(cè),說明斷層南側(cè)以向東運(yùn)動(dòng)為主,斷層北側(cè)則以向西運(yùn)動(dòng)為主(圖5(a));南北向位移v大于零的部分位于斷層北側(cè),小于零的部分主要集中于斷層?xùn)|南部,說明斷層南側(cè)以向南運(yùn)動(dòng)為主,斷層北側(cè)以向北運(yùn)動(dòng)為主,呈現(xiàn)拉張的運(yùn)動(dòng)方式(圖5(b))。以上同震位移場(chǎng)分布特征與瑪多地震具有明顯的左旋走滑機(jī)制一致;垂向位移w大于零的部分位于斷層北部,小于零的部分位于斷層南部,斷層北部為上升盤,南部為下降盤,上升盤為斷層下盤,下降盤為上盤(圖5(c)),表明發(fā)震斷層具有正斷分量。

圖5 瑪多地震的地表同震位移分量等值線圖Fig.5 Contour maps of surface coseismic displacement components of Maduo earthquake

本文計(jì)算結(jié)果中瑪多地震引起的研究區(qū)域內(nèi)主要斷層10 km深度的應(yīng)力變化,如圖6所示。圖中正應(yīng)力變化大于零的區(qū)域表示拉應(yīng)力增加,主要位于斷層南西側(cè)和北東側(cè),與拉應(yīng)力減小的區(qū)域交錯(cuò)分布(圖6(a));剪應(yīng)力變化在發(fā)震斷層兩側(cè)是小于零的,在斷層尖端向外是大于零的(圖6(b));庫(kù)侖應(yīng)力變化在斷層兩側(cè)基本上是小于零的,庫(kù)侖應(yīng)力變化大于零的區(qū)域以斷層兩端為中心(圖6(c)和6(d)中的白色曲線分別是庫(kù)侖應(yīng)力變化的0和0.001 MPa的等值線),在斷層兩側(cè)呈蝴蝶狀對(duì)稱分布(圖6(c)、6(d))。9～11 km深度范圍內(nèi)除主震斷層帶兩側(cè)大于1.0級(jí)的余震基本位于10 km深度水平面上平行于主震的庫(kù)侖應(yīng)力變化大于零的區(qū)域(圖6(d))。

黑色空心圓圈表示9～11 km深度的M≥1.0的地震。圖6 瑪多地震引起的深度為10 km水平面與主震震源機(jī)制一致的微元面上的同震應(yīng)力變化Fig.6 The coseismic stress changes on the plane consistent with the focal mechanism of the main earthquake of Maduo earthquake at 10 km depth

本文同時(shí)計(jì)算了研究區(qū)域內(nèi)與瑪多地震發(fā)震斷層相鄰的5條斷層上的同震庫(kù)侖應(yīng)力變化,如圖7,計(jì)算使用的有效摩擦系數(shù)為0.4[24]。結(jié)果表明,瑪多地震引起鄂拉山斷裂南段、昆中斷裂東段、東昆侖斷裂帶的大部分(除中部靠西的小部分外)、瑪多—甘德斷裂北西段、南東段及主震震中附近以及達(dá)日斷裂帶北西段庫(kù)侖應(yīng)力增加,地震危險(xiǎn)性增加,這些地方是需要密切關(guān)注、加強(qiáng)地震監(jiān)測(cè)的區(qū)域。

F1：達(dá)日斷裂帶;F2：瑪多—甘德斷裂帶;F3：東昆侖斷裂帶;F4：昆中斷裂帶;F5：鄂拉山斷裂帶,斷層數(shù)據(jù)引自詹艷等[4]。五角星為2021年5月22日青海瑪多MS7.4地震的震中。圖7 瑪多地震震中周圍主要斷裂上的庫(kù)侖應(yīng)力變化Fig.7 Coulomb stress change on the major faults around the epicenter of Maduo earthquake

3 討論

3.1 瑪多地震地表同震變形與大地測(cè)量觀測(cè)數(shù)據(jù)及與前人計(jì)算結(jié)果的對(duì)比

為驗(yàn)證本文采用分裂節(jié)點(diǎn)有限元方法模擬實(shí)際大地震同震變形的有效性,本文設(shè)計(jì)了瑪多地震的2種不同介質(zhì)模型：均勻介質(zhì)模型和分層介質(zhì)模型,2種模型材料參數(shù)如表2。本文分別計(jì)算了2種介質(zhì)模型的同震變形場(chǎng),并將其處理為衛(wèi)星視線方向LOS位移量,與GPS觀測(cè)的同震水平位移、InSAR升降軌LOS位移量觀測(cè)結(jié)果進(jìn)行對(duì)比,結(jié)果如圖8所示。圖8(a)、8(b)中的黑線表示瑪多地震發(fā)震斷層,藍(lán)線表示有限元模型模擬結(jié)果,紅線表示GPS觀測(cè)結(jié)果(低頻),粉線表示GPS觀測(cè)結(jié)果(高頻),左側(cè)是均勻有限元模型模擬結(jié)果,右側(cè)是分層有限元模型模擬結(jié)果。圖8(c)、8(e)、8(g)分別是InSAR升軌的觀測(cè)數(shù)據(jù)、均勻和分層有限元模型模擬結(jié)果。圖8(d)、8(f)、8(h)分別是InSAR降軌的觀測(cè)數(shù)據(jù)、均勻和分層有限元模型模擬結(jié)果,表明均勻有限元和分層有限元模型得到的升軌和降軌的衛(wèi)星視線方向LOS位移量與InSAR觀測(cè)得到的衛(wèi)星視線方向LOS位移量有較好的一致性。

本文采用統(tǒng)計(jì)均方根(root mean square, RMS)方法對(duì)有限元模擬結(jié)果和觀測(cè)變形量之間的差異進(jìn)行定量說明。RMS[15]定義如下

(8)

表3 瑪多地震同震位移的有限元模擬結(jié)果與大地測(cè)量觀測(cè)位移量之間的RMSTable 3 RMS comparison of finite element simulation results and geodetic observed coseismic displacement of Maduo earthquake

從表3可以看出,有限元模擬同震位移的結(jié)果與大地測(cè)量觀測(cè)結(jié)果之間的RMS在合理范圍之內(nèi),與楊君妍等[15]的計(jì)算結(jié)果一致,說明本文所采用的有限元方法對(duì)大地震同震變形的模擬是有效的。

本文通過分裂節(jié)點(diǎn)法將位錯(cuò)數(shù)據(jù)加入有限元模型計(jì)算的同震位移結(jié)果表明,發(fā)震斷層南側(cè)以向東、向南及向下運(yùn)動(dòng)為主,斷層北側(cè)則以向西、向北及向上運(yùn)動(dòng)為主,總體表現(xiàn)為有正斷分量的左旋走滑斷層。該計(jì)算結(jié)果與楊君妍等[15]使用球形分層地球模型位錯(cuò)理論、談洪波等[19]使用均勻半空間位錯(cuò)理論進(jìn)行同震位錯(cuò)計(jì)算的解析法或半解析法所得結(jié)果基本一致。

3.2 瑪多地震對(duì)巴顏喀拉地塊及周邊地震危險(xiǎn)性的影響

本文瑪多地震同震庫(kù)侖應(yīng)力計(jì)算結(jié)果表明,發(fā)震斷層的兩側(cè)沿走向方向上庫(kù)侖應(yīng)力下降,而在斷層兩端及以兩端為中心呈蝴蝶狀發(fā)散的區(qū)域內(nèi),庫(kù)侖應(yīng)力增加。這個(gè)現(xiàn)象表明瑪多地震對(duì)發(fā)震斷層及斷層周邊地殼產(chǎn)生了顯著的應(yīng)力調(diào)整：瑪多地震部分釋放了發(fā)震斷層所積累的應(yīng)力(發(fā)震斷層兩側(cè)的庫(kù)侖應(yīng)力變化大部分都是顯著降低),也造成斷層兩端及周邊部分區(qū)域的地殼應(yīng)力積累(圖6(c))。該結(jié)果與楊光遠(yuǎn)等[16]利用Coulomb 3.3程序在同等條件(10 km深度處和主震震源機(jī)制一致的微元面)下的計(jì)算結(jié)果對(duì)比,同震庫(kù)侖應(yīng)力分布基本一致。我們的模擬結(jié)果還表明,9～11 km深度范圍內(nèi)除主震斷層帶兩側(cè)大于1.0級(jí)的余震基本位于10 km深度水平面上平行于主震的庫(kù)侖應(yīng)力變化大于零的區(qū)域(圖6(d))。研究區(qū)域內(nèi)鄂拉山斷裂、昆中斷裂、東昆侖斷裂、瑪多—甘德斷裂以及達(dá)日斷裂這5條主要斷裂的同震庫(kù)侖應(yīng)力結(jié)果表明：瑪多地震引起鄂拉山斷裂南段,昆中斷裂東段,東昆侖斷裂帶的大部分、瑪多—甘德斷裂西北段和東南段以及達(dá)日斷裂帶北西段庫(kù)侖應(yīng)力增加(圖7)。鄂拉山斷裂南段與昆中斷裂東段位于巴顏喀拉塊體北部的東昆侖—柴達(dá)木塊體,歷史上該地區(qū)曾發(fā)生過M5～M6及M6～M7地震,此次瑪多地震引起該地區(qū)庫(kù)侖應(yīng)力增加,具有較高的地震危險(xiǎn)性,應(yīng)對(duì)其密切關(guān)注。東昆侖斷裂帶歷史上曾發(fā)生過如1937年MS7.5地震(阿拉克湖—托索湖段和下大武段)、1963年MS7.0地震(秀溝—阿拉克湖段)以及2001年MS7.1地震(庫(kù)賽湖段),是一條活躍斷裂帶[25]。瑪多地震使瑪沁—瑪曲段東南段和托索湖—下大武段東南段庫(kù)侖應(yīng)力增加,這與Li等[20]的結(jié)果一致。這2個(gè)斷裂段在瑪多地震時(shí)產(chǎn)生應(yīng)力積累,加之瑪沁—瑪曲段近年來未發(fā)生大地震,導(dǎo)致該地區(qū)地震危險(xiǎn)性大大增加,應(yīng)當(dāng)加強(qiáng)地震監(jiān)測(cè)、做好防震防災(zāi)工作。瑪多—甘德斷裂中段庫(kù)侖應(yīng)力變化分布較為復(fù)雜,庫(kù)侖應(yīng)力顯著增加與顯著降低的區(qū)段交錯(cuò)分布,可能與其復(fù)雜的斷層結(jié)構(gòu)有關(guān)[4]?，敹唷实聰嗔巡糠謪^(qū)段庫(kù)侖應(yīng)力顯著降低可能與瑪多地震釋放了其部分應(yīng)力有關(guān),但其兩端的庫(kù)侖應(yīng)力顯著增加,地震危險(xiǎn)性有所提高。位于瑪多地震發(fā)震斷層南部的達(dá)日斷裂帶北段及中段庫(kù)侖應(yīng)力增加,而其余部分庫(kù)侖應(yīng)力水平下降,考慮到1947年達(dá)日M7.7地震所帶來的達(dá)日斷裂帶西北段庫(kù)侖應(yīng)力增加及其較高的應(yīng)力積累速率[26],其地震危險(xiǎn)性仍需進(jìn)一步研究與評(píng)估。

4 結(jié)論

本文采用分裂節(jié)點(diǎn)技術(shù)開發(fā)三維并行彈性有限元程序,并通過地震位錯(cuò)模型的理論解驗(yàn)證了有限元程序的正確性。該有限元程序能夠考慮復(fù)雜幾何、非均勻材料和位錯(cuò)非均勻分布,可以模擬計(jì)算大地震同震變形。利用該程序模擬計(jì)算了瑪多MS7.4地震的同震位移場(chǎng)、應(yīng)力場(chǎng)及庫(kù)侖應(yīng)力變化,計(jì)算得到的地表同震位移結(jié)果與GPS和InSAR大地測(cè)量數(shù)據(jù)一致。模擬結(jié)果確認(rèn)了瑪多MS7.4地震的發(fā)震斷層為一條左旋走滑兼有正斷分量的斷層。計(jì)算得到的庫(kù)侖應(yīng)力增加的區(qū)域與瑪多地震的余震分布有較好的對(duì)應(yīng)性。庫(kù)侖應(yīng)力結(jié)果表明,瑪多地震引起了鄂拉山斷裂南段,昆中斷裂東段,東昆侖斷裂帶的大部分區(qū)段、瑪多—甘德斷裂西北段和東南段以及達(dá)日斷裂帶北西段庫(kù)侖應(yīng)力增加,這些區(qū)域的地震危險(xiǎn)性增加,是未來需重點(diǎn)關(guān)注的地區(qū)。

石耀霖院士和周永發(fā)博士對(duì)本研究自主研發(fā)的三維有限元程序編制給予了很多建設(shè)性的指導(dǎo)和幫助,周元澤教授對(duì)論文初稿提出了很多建設(shè)性的意見,在此一并表示感謝。