談洪波 王嘉沛 楊光亮 陳正松吳桂桔 申重陽 黃金水

1)中國科學(xué)技術(shù)大學(xué)，地球和空間科學(xué)學(xué)院，合肥 230026 2)中國地震局地震研究所，武漢 430071

0 引言

地震發(fā)生后，中國地震局迅速組織實(shí)施了地震科學(xué)考察，開展了多手段聯(lián)合觀測(包括重力和GNSS等)，為震后地震趨勢判定、 地震應(yīng)急救援等工作提供資料。 震源理論的發(fā)展及觀測數(shù)據(jù)的不斷豐富，使我們可在地震發(fā)生后快速確定震源的基本參數(shù)，也使得將地震位錯(cuò)理論應(yīng)用于震后趨勢分析成為可能。 Zhang等(2014)發(fā)展了一套反演地震破裂過程的迭代反褶積方法，并聯(lián)合高頻GPS和強(qiáng)震數(shù)據(jù)對(duì)日本MW9.0地震的破裂過程進(jìn)行了反演。 劉剛等(2015)利用遠(yuǎn)場地震波和GPS數(shù)據(jù)給出了2015年尼泊爾MS8.1地震的破裂過程。 Feng等(2017)利用InSAR和GPS數(shù)據(jù)反演了2015年尼泊爾MW7.8地震及其MW7.2余震的斷層滑動(dòng)分布。 談洪波等(2013)模擬了玉樹MS7.1地震的同震重力變化并利用絕對(duì)重力觀測結(jié)果進(jìn)行了檢驗(yàn)。 Zhang等(2016)對(duì)日本MW9.0地震的同震模擬結(jié)果及觀測形變、 重力變化進(jìn)行了比較，發(fā)現(xiàn)兩者符合較好。 Riguzzi等(2019)模擬了均勻彈性半空間中2009—2016年意大利中部4次MW6地震的同震重力效應(yīng)，并就模擬結(jié)果進(jìn)行了合理解釋。 但實(shí)際地球要復(fù)雜得多，地球曲率、 層狀構(gòu)造、 黏滯性、 自重效應(yīng)、 三維構(gòu)造等均會(huì)對(duì)計(jì)算結(jié)果產(chǎn)生不同程度的影響(Pollitz，1996； Sunetal.，2002； 孫文科，2008，談洪波等，2010)。 同震只是介質(zhì)在完全彈性情況下的瞬時(shí)效應(yīng)，對(duì)于長時(shí)間尺度的震后效應(yīng)而言，其形變和重力變化有可能完全改變同震效應(yīng)的形態(tài)，利用震后數(shù)據(jù)獲取斷層和地球內(nèi)部信息也是非常有意義的研究內(nèi)容(張晁軍等，2008； 賀鵬超等，2018； 趙斌等，2020)。 因此，本文將采用更為接近實(shí)際情況的彈性—黏彈分層介質(zhì)模型，并考慮自重的影響，從形變和重力的角度對(duì)該地震的同震及其震后效應(yīng)進(jìn)行理論模擬，為利用實(shí)際形變和重力資料解釋此次地震孕震過程(包括同震和震后)的相關(guān)研究提供理論依據(jù)。

1 理論基礎(chǔ)

均勻彈性半無限空間地球模型的位錯(cuò)理論由Okada(1985)和Okubo(1992)完善和發(fā)展。 假設(shè)地震破裂發(fā)生在一個(gè)平面上，斷層如圖 2 所示，長寬分別為L和W，傾角為δ； 斷層底部深度為d； 斷層面上的走滑、 傾滑和張裂方向錯(cuò)動(dòng)量分別為U1、U2、U3，在地表(x1，x2，0)處引起的形變?chǔ)和重力變化Δg可分別表示為(Okada，1985； Okubo，1992)

圖1 研究區(qū)主要構(gòu)造分布Fig. 1 Distribution of the main structures in the study area.黑點(diǎn)為截至2021年5月30日的ML>2余震

圖 2 矩形斷層位錯(cuò)模型示意圖Fig. 2 Dislocation model of a rectangle fault.

圖 3 SLS流變模型Fig. 3 The SLS rheological model.

(1)

Δg(x1，x2)={ρG[U1Sg(ξ，η)+U2Dg(ξ，η)+U3Tg(ξ，η)]+ΔρGU3Cg(ξ，η)}‖-βΔh(x1，x2)

(2)

其中，G為萬有引力常數(shù)；ρ、 Δρ分別為介質(zhì)密度、 張裂紋內(nèi)密度與介質(zhì)密度之差；S(ξ，η)、D(ξ，η)、T(ξ，η)和Sg(ξ，η)、Dg(ξ，η)、Tg(ξ，η)、Cg(ξ，η)為系數(shù)； 自由空氣重力梯度β=0.309×10-5m/s2，Δh為地表高程變化，Chinnery記號(hào)‖表示置換(Okada，1985)。

考慮地球自重影響時(shí)，Rundle(1980)首次給出了基于無限小靜態(tài)形變理論動(dòng)力平衡方程：

·Γ+ρψ+u·g-g(ρu)=0

(3)

(4)

Γ=(λ·u)I+μ[u+(u)T]

(5)

式中，λ為拉梅常數(shù)；μ為切變模量；I為單位張量。 根據(jù)Adams-Williamson條件(Longman，1963)，方程(3)可變換為(Wangetal.，2006)

Γ+ρguz-(ρ

(6)

對(duì)于黏彈介質(zhì)，一般基于彈性-黏彈線性對(duì)應(yīng)原理(Fung，1965)： 黏彈體在頻率域的應(yīng)力-應(yīng)變關(guān)系與彈性體在時(shí)間域的一致，故在頻率域中黏彈問題可采用彈性問題的處理方式計(jì)算后再通過逆傅里葉變換回到時(shí)間域(Matsuuraetal.，1981)。 Wang等(2006)采用SLS(Standard Linear Solid)流變體將彈性體和黏彈性體用統(tǒng)一的模型表示(圖 2)，這時(shí)復(fù)切變模量和復(fù)拉梅常數(shù)可表示為

(7)

(8)

聯(lián)合方程(5)—(8)，根據(jù)正則化Haskell(1953)傳播算法，可將邊值問題轉(zhuǎn)變成代數(shù)問題，根據(jù)給定的地球參數(shù)求出頻率域中對(duì)應(yīng)的格林函數(shù)(即單位點(diǎn)脈沖位錯(cuò)響應(yīng))，再采用抗假頻技術(shù)(Kindetal.，1982)對(duì)頻率域中的格林函數(shù)進(jìn)行傅里葉逆變換回到時(shí)間域，最后通過式(1)和(2)的線性疊加求出地表形變和重力變化( Wang，1999； Wangetal.，2006)。

2 介質(zhì)模型和發(fā)震斷層模型

前期的研究結(jié)果表明地殼分層對(duì)地表同震效應(yīng)的影響可達(dá)20%以上(Sunetal.，2002； 孫文科，2008； 談洪波等，2010)，因此需要考慮垂向分層結(jié)構(gòu)。 本文參考Crust1.0(4)https:∥igppweb.ucsd.edu/～gabi/crust1.html。地殼模型，將研究區(qū)地球分為上、 中、 下地殼和上地幔，具體參數(shù)見表1。 此外，也考慮了下地殼和上地幔介質(zhì)的黏滯性，黏滯系數(shù)η參考張晁軍等(2008)、 賀鵬超等(2018)和趙斌等(2020)對(duì)青藏高原的研究結(jié)果。

表1 地殼—上地幔介質(zhì)分層結(jié)構(gòu)模型Table1 The layered structure of the crust-upper-mantle based on Crust 1.0

圖4 斷層滑動(dòng)分布(據(jù)USGS結(jié)果繪制)Fig. 4 The slip distribution on the fault plane(according to USGS results).

3 同震效應(yīng)

圖5 地表同震效應(yīng)Fig. 5 Surface co-seismic effects.a 經(jīng)度方向位移； b 緯度方向位移； c 垂直位移； d 重力變化

4 震后效應(yīng)

同時(shí)我們也對(duì)瑪多、 稱多和達(dá)日縣城處的震后形變和重力變化趨勢進(jìn)行了模擬(表2，圖 7)，這3個(gè)點(diǎn)的震中距分別為36km、 172km和156km。 圖7b為經(jīng)度方向的位移，瑪多縣城同震時(shí)向W運(yùn)動(dòng)了126mm，400a后位移量將達(dá)197mm，增加了71mm(同震效應(yīng)的56%)； 同震時(shí)稱多和達(dá)日縣城向E分別運(yùn)動(dòng)了8mm和19mm，400a后位移量將達(dá)39mm和43mm，分別增加了31mm(同震效應(yīng)的387%)和24mm(同震效應(yīng)的126%)； 整體而言，瑪多縣城震后EW向位移量在300a后基本穩(wěn)定不變，而稱多和達(dá)日縣城在50a后基本穩(wěn)定不變。 圖7c為緯度方向的位移，瑪多和稱多縣城同震時(shí)分別向N運(yùn)動(dòng)了151mm和6mm，400a后位移量將達(dá)198mm和18mm，分別增加了47mm(同震效應(yīng)的31%)和12mm(同震效應(yīng)的200%)； 同震時(shí)達(dá)日縣城向S運(yùn)動(dòng)了22mm，400a后位移量將達(dá)67mm，增加了45mm(同震效應(yīng)的204%)； 瑪多縣城震后SN向位移量在100a后基本穩(wěn)定不變，而稱多和達(dá)日縣城的位移量基本穩(wěn)定所需時(shí)間分別為20a和100a。

表2 部分縣城位置的震后效應(yīng)Table2 The post-seismic effects on some of the cities

圖 7d 為垂向位移，同震時(shí)瑪多縣城下沉3mm，5a后累計(jì)下沉量達(dá)5.6mm，隨后地表開始持續(xù)抬升，400a后抬升量將達(dá)46.5mm，400a間地表抬升量達(dá)49.5mm(同震效應(yīng)的16.5倍)； 稱多縣城在同震時(shí)地表抬升了3mm，震后5a間有微弱的下沉，隨后持續(xù)抬升，400a后抬升量達(dá)31.7mm，抬升量增加了28.7mm(同震效應(yīng)的9.5倍)； 而達(dá)日在同震時(shí)下沉4.3mm，震后5a間有微弱的抬升，隨后持續(xù)下沉，400a后下沉量達(dá)65.9mm，下沉量增加了61.6mm(同震效應(yīng)的14.3倍)； 3個(gè)縣城震后的垂直位移量在300a后基本穩(wěn)定不變。 圖7e為重力變化，與垂直位移震后變化趨勢類似并呈現(xiàn)一定的負(fù)相關(guān)性： 同震時(shí)瑪多縣城的重力值增加5.0μGal，5a后累計(jì)增加到5.4μGal，隨后開始持續(xù)減小，400a后重力值減小為4.4μGal，400a間重力值減小了9.4μGal； 稱多縣城在同震時(shí)地表重力值減小了1μGal，震后5a間有微弱的增加(0.08μGal)，隨后持續(xù)減小，400a后重力值減小為6.5μGal，400a間重力值減小了5.5μGal； 而達(dá)日在同震時(shí)的重力值增加2.0μGal，震后5a間有微弱的減小(0.19μGal)，隨后持續(xù)增大，400a后重力值增大為13.5μGal，400a間重力值增大了11.5μGal； 3個(gè)縣城震后地表重力變化在300a后基本穩(wěn)定不變。

圖7 震后長期時(shí)變效應(yīng)Fig. 7 Long-term time-varying effects after the earthquake.a 所選點(diǎn)與斷層的相對(duì)位置； b 經(jīng)度方向位移； c 緯度方向位移； d 垂向位移； e 重力變化

5 討論與結(jié)論

文中基于彈性-黏彈分層半空間中的矩形位錯(cuò)理論，采用USGS發(fā)布的由遠(yuǎn)場地震波反演的瑪多MS7.4地震發(fā)震斷層模型，從理論上模擬了此次地震引起的地表同震及震后形變和重力變化，通過分析可得到如下主要結(jié)論：

圖8 震后400a間的黏彈性松弛效應(yīng)與同震效應(yīng)比值分布Fig. 8 The distribution of the ratios between the viscoelastic relaxation effects during 400 years after the Maduo earthquake and the coseismic effects.a 經(jīng)度方向位移； b 緯度方向位移； c 垂向位移； d 重力變化

圖9 震后5a間的黏彈性松弛效應(yīng)與同震效應(yīng)比值分布Fig. 9 The distribution of the ratios between the viscoelastic relaxation effects during 5 years after the Maduo earthquake and the coseismic effects.a 經(jīng)度方向位移； b 緯度方向位移； c 垂向位移； d 重力變化

圖10 同震模擬結(jié)果與GNSS觀測結(jié)果的比較Fig. 10 Comparison between the coseismic simulation and GNSS observation results.藍(lán)色和紅色箭頭為GNSS觀測結(jié)果，紫色和綠色為模擬結(jié)果

致謝USGS提供了發(fā)震斷層模型，中國地震局地震研究所趙斌博士提供了同震GNSS數(shù)據(jù)，德國波茨坦地球科學(xué)中心汪榮江教授提供了PSGRN/PSCMP程序； 審稿專家為本文提供了寶貴建議。 在此一并表示感謝!