鐘慧婷, 靳志同,2, 萬永革,2

(1. 防災(zāi)科技學(xué)院, 河北 三河 065201;2. 河北省地震動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室, 河北 三河 065201)

0 引言

確定斷層的破裂方式對(duì)研究區(qū)的整個(gè)地震孕育過程研究及其震后地震趨勢(shì)的發(fā)展均起到非常重要的作用,而重力場(chǎng)是反映地球介質(zhì)密度變化和各種環(huán)境下動(dòng)力學(xué)特征的最基本、最直接的物理量,重力-地震的聯(lián)合反演具有顯著的互補(bǔ)效果和更完備額的地質(zhì)學(xué)解釋[1],所以重力資料在反演震源參數(shù)和震源破裂中起著不可或缺的作用。如賈可可等[2]利用同震重力變化反演出日本“3·11”地震的震源破裂參數(shù)。鄭增記等[3]基于GRACE記錄的重力資料反演得到2012年蘇門答臘地震斷層參數(shù)。楊君妍等[4]通過GRACE重力衛(wèi)星獲取的數(shù)據(jù)反演出了地震位錯(cuò)Love數(shù),并且得到了更準(zhǔn)確的格林函數(shù),反演出孕震區(qū)域局部構(gòu)造信息。

除此之外,在地震預(yù)報(bào)領(lǐng)域,重力資料也發(fā)揮了一定的作用。祝意青等指出強(qiáng)震易發(fā)生在與構(gòu)造活動(dòng)有關(guān)聯(lián)的重力變化正、負(fù)異常區(qū)過渡的高梯度帶上,重力變化等值線的拐彎部位或四象限中心附近[5]。2001年昆侖山口西8.1級(jí)地震發(fā)生在重力變化高梯度帶上,重力負(fù)值變化低于-40 μgal[6];2021年云南漾濞6.4級(jí)地震震中處于重力變化的高梯度帶區(qū)域,判斷漾濞地區(qū)地下活動(dòng)的能量已經(jīng)完全釋放完,震后短期內(nèi)無更大地震發(fā)生的可能[7]。這些研究表明重力資料可以為地球動(dòng)力學(xué)過程、地震危險(xiǎn)性判斷提供重要信息。由于重力監(jiān)測(cè)資料中包含了地震斷層錯(cuò)動(dòng)產(chǎn)生的重力變化值,而地震產(chǎn)生的重力變化會(huì)對(duì)揭示地震孕育等長(zhǎng)期動(dòng)力學(xué)行為造成干擾,因此將地震導(dǎo)致的重力變化扣除更有利于揭示地震孕育等地球動(dòng)力學(xué)過程。

2022年1月8日,青海門源發(fā)生了6.9級(jí)地震,震中位于101.28°E、37.82°N,震源深度為10 km,造成了房屋建筑受損、蘭新高鐵受損停運(yùn)等問題。此次地震位于青藏高原東北部邊緣的祁連山地震帶,由一系列的深大活動(dòng)斷裂組成。作為地震活動(dòng)主體區(qū)域的祁連山次級(jí)地塊,四周被深大走滑活動(dòng)斷裂包圍,形成一個(gè)相對(duì)獨(dú)立的較為活動(dòng)的次級(jí)地塊[8]。此次6.9級(jí)地震的發(fā)震構(gòu)造既不是冷龍嶺斷裂,也不是托萊山斷裂,而是發(fā)育在它們之間階區(qū)中的道溝斷裂[9-10],主要是由淺層的左旋走滑斷層造成的。

(紅色沙灘球表示本次地震主震的震源機(jī)制解,藍(lán)色三角號(hào)表示研究區(qū)域的前兆臺(tái),紅色線代表研究區(qū)域內(nèi)不同的斷層,左下角子圖表示研究地區(qū)所在的位置)圖1 研究區(qū)地形圖Fig.1 Topographic map in the study region

近100年來,該地區(qū)已發(fā)生高于6.0級(jí)的地震20余次,其中發(fā)生的距離此次地震最近的一次大地震是2016年門源縣6.4級(jí)地震。地震斷層錯(cuò)動(dòng)會(huì)導(dǎo)致地面重力變化,而中國(guó)地震前兆臺(tái)網(wǎng)的重力儀可以觀測(cè)到這種變化。這種變化信息對(duì)確定斷層的錯(cuò)動(dòng)方式、大小等都起到非常重要的作用。地震前后重力變化的原因一般有兩方面:一是震前地震孕育過程中的質(zhì)量遷移和地表形變引起的重力變化,二是同震位錯(cuò)引起周圍物理密度的變化;前一種原因往往是重力變化的主要因素,后一種原因?qū)е碌闹亓ψ兓话阍诰嚯x震中幾十千米以外處只有約10 μgal量級(jí)且隨著距離的增加迅速衰減[11]。

2022年青海門源發(fā)生6.9級(jí)地震,該地震發(fā)生在周圍地震前兆臺(tái)布設(shè)比較密集的地區(qū),為了甄別重力前兆數(shù)據(jù)中的該地震錯(cuò)動(dòng)引起的重力變化,從而為提取有意義的重力地震前兆提供基礎(chǔ)數(shù)據(jù),并了解該地震對(duì)周圍物質(zhì)運(yùn)移和密度變化的影響,本文采用Okubo地震位錯(cuò)模型,分析了各種類型的地震錯(cuò)動(dòng)引起的地表重力變化的特征,并計(jì)算了2022年青海門源縣6.9級(jí)地震在周圍地區(qū)引起的重力變化。

1 Okubo理論及計(jì)算公式

國(guó)內(nèi)外許多學(xué)者對(duì)位錯(cuò)和斷層引起的重力變化進(jìn)行過深入研究。Stekettee[12]最早將位錯(cuò)理論引入地震學(xué),Maruyama[13]給出了垂直和水平伸展斷層引起的地表位移場(chǎng)的完整解析解,陳運(yùn)泰等[14-15]討論了結(jié)合半空間無限空間位錯(cuò)理論和地表形變資料進(jìn)行震源反演的一般方法,Okada[16]總結(jié)了前人的工作,給出了一套完整簡(jiǎn)潔實(shí)用的同震變形計(jì)算公式。Okubo[11,17-18]研究了半無限空間介質(zhì)內(nèi)剪切和張裂斷層的重力變化問題,導(dǎo)出了點(diǎn)源和有限斷層的同震重力位和重力變化的解析表達(dá)式,Sun[19]解決了球形地球分層模型中地震位錯(cuò)的同震重力變化計(jì)算問題。李振洪等利用InSAR技術(shù)獲取2022年青海門源地震的同震形變場(chǎng),可以發(fā)現(xiàn)本次地震形變僅限局部地區(qū),同時(shí)利用了彈性半空間的位錯(cuò)模型確定了本次地震事件的斷層幾何參數(shù)[20],鑒于2022年青海門源地震產(chǎn)生的重力變化僅限于局部地區(qū),本文介紹Okubo的半空間線彈性介質(zhì)中地震位錯(cuò)產(chǎn)生的重力變化計(jì)算方法。

Okubo理論主要計(jì)算了有限矩形斷層的重力變化,根據(jù)最后的推導(dǎo)結(jié)果,用雙豎約定符號(hào)標(biāo)示為簡(jiǎn)潔模式[21]:

f(ξ,η)‖≡f(x1,p)-f(x1,p-W)-
f(x1-L,p)+f(x1-L,p-W)

(1)

重力變化可以表示為:

Δg(x1,x2)={ρG[U1Sg(ξ,η)+U2Dg(ξ,η)+
U3·Tg(ξ,η)]+ΔρGU6g(ξ,η)}‖-
βΔh(x1,x2)

(2)

式中:G為萬有引力常數(shù);ρ為介質(zhì)密度;Δρ=ρ′-ρ,ρ′為引張破裂后的填充物質(zhì)的密度。S、D、T、C則表示走滑、傾滑、引張和引張破裂填充物的貢獻(xiàn),(Sg,Dg,Tg,Cg)是對(duì)(S,D,T,C)進(jìn)行微分得到的,β≈0.309 mgal/m為自由空氣重力梯度,Δh為高程變化,它們的具體表達(dá)式為:

(3)

(4)

(5)

Cg(ξ,η)=2I2cosδ-sinδ·lg(R+ξ)

(6)

q=x2sinδ-(d-x3)cosδ

(7)

(8)

同樣,高程變化可以表示為:

(9)

具體表達(dá)式為:

(10)

(11)

(12)

其中:

(13)

(14)

(15)

I5(ξ,η)=2(1-2v)I1secδ

(16)

當(dāng)cosδ=0時(shí):

(17)

(18)

式中:δ為傾角;v為泊松系數(shù);d為有限矩形斷層的深度。

2 典型錯(cuò)動(dòng)斷層在其周圍產(chǎn)生的重力變化規(guī)律

為了解典型斷層引起的重力變化規(guī)律,本小節(jié)計(jì)算了單個(gè)左旋走滑斷層、正斷層、逆斷層的重力變化情況,并且進(jìn)行了分析。走向上單個(gè)斷層長(zhǎng)度為10 km,傾角方向上的斷層寬度為10 km,在滑動(dòng)方向上的位錯(cuò)為5 m,取介質(zhì)密度為2 670 kg/m3,計(jì)算結(jié)果分析如下:

2.1 走滑斷層重力變化特征

圖2為單個(gè)走滑斷層引起的重力變化分布圖,走向?yàn)檎龞|,傾角為90°,滑動(dòng)角為0°,其中AA′為計(jì)算圖3的剖面位置。圖3為破裂中心位于不同深度的地表AA′剖面的重力變化曲線。

圖2 走滑斷層重力變化(AA′為圖3的剖面位置)Fig.2 Gravity changes induced by strike-slip fault (AA′is the position of profile in Fig.3)

圖3 同一剖面破裂中心位于不同深度的走滑斷層產(chǎn)生的重力變化(黑色水平線代表零重力值)Fig.3 Gravity changes along profile induced by strike-slip fault with fracture center at different depths (The black horizontal line represents the zero gravity value)

由圖2可以看出,垂直走滑斷層產(chǎn)生的重力變化總體呈現(xiàn)出四象限對(duì)稱分布;在近處,以斷層跡線逆時(shí)針?biāo)闫?一三象限重力變化為正值,二四象限為負(fù)值;在遠(yuǎn)處,一三象限重力變化為負(fù)值,二四象限為正值。圖3顯示,同一剖面重力變化曲線關(guān)于斷層所在位置對(duì)稱,重力變化大小隨著遠(yuǎn)離斷層而逐漸縮小;對(duì)于相同破裂尺度和滑動(dòng)量的垂直走滑斷層,破裂中心深度越深,震中附近的重力變化幅度越小,而自斷層位置向遠(yuǎn)處的重力變化幅度衰減越慢。反之,震源破裂中心越淺,斷層產(chǎn)生的重力變化局限在斷層附近,峰值較大,但很快衰減。也就是說,相同破裂尺度和破裂量的淺表震源引起的重力變化值局限在震中的局部地區(qū),而相同破裂尺度和破裂量的深部震源引起的重力變化會(huì)產(chǎn)生更大范圍的重力變化,但重力變化值較小。

2.2 傾滑斷層重力變化特征

圖4與圖7為單個(gè)正斷層(滑動(dòng)角為-90°)和逆斷層(滑動(dòng)角為90°)引起的重力變化分布圖,走向?yàn)檎龞|,傾角為45°。以正斷層為例,圖5為相同破裂中心位于不同深度的AA′剖面上的重力變化曲線,圖6為破裂中心位于6 km處,不同傾角的斷層破裂在AA′剖面產(chǎn)生的重力變化。逆斷層同理,圖像為7～9。

圖4 正斷層重力變化(AA′為圖5的剖面位置)Fig.4 Gravity variation induced by normal fault (AA′ is the position of profile in Fig.5)

圖5 同一剖面破裂中心位于不同深度的正斷層產(chǎn)生的重力變化Fig.5 Gravity changes along profile induced by normal fault with fracture center at different depths

圖6 同一剖面破裂中心位于6 km處不同傾角正斷層重力變化Fig.6 Gravity changes along profile induced by normal fault with fracture center at a depth of 6 km and different dip angles

圖7 逆斷層重力變化(AA′為計(jì)算8的剖面位置)Fig.7 Gravity variation induced by reverse fault (AA′ is the position of profile in Fig.8)

由于傾滑斷層在震源近處引起的重力變化值過大,為了更好地描述遠(yuǎn)處的重力變化特征,本文將重力變化大于5 μgal的值不顯示,如圖4和7?？梢园l(fā)現(xiàn),在遠(yuǎn)處,正斷層在震源南北方向產(chǎn)生的重力變化為正值,在東西方向產(chǎn)生的重力變化為負(fù)值,逆斷層反之。圖5和圖8顯示,上盤重力變化的數(shù)值總體上比下盤重力變化的數(shù)值小;對(duì)于相同破裂尺度和滑動(dòng)量的傾滑斷層,破裂中心深度越深,遠(yuǎn)處的重力變化最大值越小。

圖8 同一剖面破裂中心位于不同深度的逆斷層的重力變化Fig.8 Gravity changes along profile induced by reverse fault with fracture center at different depths

從圖6和圖9可看出,在傾角小于45°時(shí),斷層上盤產(chǎn)生的重力變化峰值較大,且隨著斷層的傾角的增大,在斷層上盤重力變化波峰越大;在傾角大于45°時(shí),斷層下盤產(chǎn)生的重力變化峰值較大,且隨著傾角的增大,斷層下盤重力變化的峰值變小。當(dāng)正斷層傾角較小時(shí),上盤重力變化主要以負(fù)值為主,下盤重力變化主要以正值為主;當(dāng)傾角較大時(shí),上盤重力變化主要為正值,下盤為負(fù)值;逆斷層反之。

圖9 同一剖面同一深度不同傾角逆斷層重力變化Fig.9 Gravity variation along profile induced by reverse fault with fracture center at a depth of 6 km and different dip angles

3 2022年門源地震產(chǎn)生的重力場(chǎng)

本節(jié)采取Okubo理論,震中周圍介質(zhì)密度選取為2 670 kg/m3,地震破裂模型采用李振洪等[20]根據(jù)InSAR資料聯(lián)合了升、降軌地表形變場(chǎng)反演的地震破裂模型;反演模型包含走向104°和109°兩個(gè)斷層段,走向104°的斷層段長(zhǎng)度為10 km,寬度為16 km,傾角為80°,滑動(dòng)角為0°,最大滑動(dòng)量為2.5 m;走向?yàn)?09°的斷層段長(zhǎng)度為20 km,寬度為16 km,傾角為80°,滑動(dòng)角為5°,最大滑動(dòng)量為3 m[20]。將發(fā)震斷層離散成子斷層反演得到:該斷層主要受左旋滑動(dòng)控制,屬于高傾角走滑型地震;子斷層最大滑動(dòng)量為3.5 m,出現(xiàn)在地下約4 km處,滑動(dòng)主要集中在地下2～7 km的區(qū)域[20]。本文利用多個(gè)子斷層疊加的方式,計(jì)算出重力變化情況,繪制出圖10。

(黑點(diǎn)表示城市,藍(lán)色三角號(hào)表示未能監(jiān)測(cè)到重力變化的臺(tái)站,紅色三角號(hào)表示能夠檢測(cè)到重力變化的臺(tái)站,黑色花瓣輪廓線是重力變化為1μgal的分界線,沙灘球是本次地震的震源機(jī)制解)圖10 2022年1月8日門源地震后重力變化Fig.10 Gravity changes generated by Menyuan earthquake on January 8, 2022

圖10顯示,本次地震重力增加的區(qū)域范圍要比重力減小的區(qū)域范圍小。在數(shù)值上,重力變化負(fù)值區(qū)最大為-8.82 μgal,正值區(qū)最大變化為28.28 μgal。震中東北和西南方向的重力減小,震中西北和東南方向的重力增加;即在距離震中遠(yuǎn)處斷層錯(cuò)動(dòng)前方重力變化為正,后方重力變化為負(fù),與前文敘述單個(gè)走滑斷層引起的重力變化特征基本一致。

重力變化一般在距離震中幾十千米以外處只有約10 μgal的量級(jí),并隨著距離的增加迅速減小[11]。本次門源地震斷層錯(cuò)動(dòng)導(dǎo)致的重力變化從震源處向四周衰減得很快,在數(shù)百公里之外幾乎不發(fā)生重力變化。當(dāng)前CG-5高精度重力儀在重力觀測(cè)上分辨率可以達(dá)到1 μgal[22],圖10中,黑色花瓣輪廓線是重力變化為1 μgal的分界線,花瓣內(nèi)是能夠采用重力儀監(jiān)測(cè)到地震錯(cuò)動(dòng)產(chǎn)生重力變化的地區(qū),震中東北和西南方向、西北和東南方向的花瓣輪廓最遠(yuǎn)大約在100 km處,超出該花瓣范圍的區(qū)域?qū)㈦y以監(jiān)測(cè)到重力變化。

圖10中有63個(gè)用三角號(hào)表示的前兆臺(tái),紅色三角(門源臺(tái))標(biāo)志表示能夠監(jiān)測(cè)到此次門源地震重力變化的臺(tái)站,藍(lán)色為不能夠監(jiān)測(cè)到重力變化的臺(tái)站。需要指出的是,在已知臺(tái)站中,有一個(gè)S9(門源臺(tái))前兆臺(tái)站能夠監(jiān)測(cè)到本次門源地震的重力變化,其余臺(tái)站處重力變化低于1 μgal,超過了重力儀監(jiān)測(cè)的范圍,所以未能監(jiān)測(cè)到重力變化。即便在本次門源地震中只有一個(gè)前兆臺(tái)能監(jiān)測(cè)到重力變化,在監(jiān)測(cè)和分析該地區(qū)長(zhǎng)期重力變化從而揭示構(gòu)造活動(dòng)時(shí),也應(yīng)該扣除該地震造成的重力變化。

重力場(chǎng)的變化,能較好地反映地殼厚度的差異、地殼密度的變化和深部物質(zhì)遷移等構(gòu)造活動(dòng)信息,重力場(chǎng)隨時(shí)間變化與地震的形成和發(fā)展有著內(nèi)在聯(lián)系[23],所以重力資料在反演震源參數(shù)中也起著重要的作用。當(dāng)然,單靠重力資料也是遠(yuǎn)遠(yuǎn)不夠的,需要聯(lián)合其他多種資料共同約束,從而得到更加準(zhǔn)確的震源參數(shù)。

4 討論與結(jié)論

為分析地震錯(cuò)動(dòng)在地表導(dǎo)致的重力變化,本研究根據(jù)Okubo位錯(cuò)理論公式,計(jì)算分析了典型斷層(走滑、傾滑)在地表的重力變化特征:(1)垂直走滑斷層產(chǎn)生的重力變化總體呈現(xiàn)出四象限對(duì)稱分布;對(duì)于同一剖面,重力變化大小隨著遠(yuǎn)離斷層而逐漸縮小;對(duì)于相同破裂尺度和滑動(dòng)量的垂直走滑斷層,深度越深,震中附近的重力變化幅度越小,而自斷層位置向遠(yuǎn)處的重力變化幅度的衰減越慢。反之,震源越淺,斷層產(chǎn)生的重力變化局限。(2)傾滑斷層上盤重力變化的數(shù)值總體上比下盤重力變化的數(shù)值小;對(duì)于相同破裂尺度和滑動(dòng)量的傾滑斷層,破裂中心深度越深,遠(yuǎn)處的重力變化最大值越小。同時(shí)本文對(duì)2022年1月8日門源6.9級(jí)地震的破裂分布進(jìn)行重力變化的計(jì)算和分析,得到:本次地震重力增加的區(qū)域范圍要比重力減小的區(qū)域范圍小;在數(shù)值上,重力變化負(fù)值區(qū)最大為-8.82 μgal,正值區(qū)最大變化為28.28 μgal,震中東北和西南方向重力降低,而震中西北和東南方向的重力增加。位于震中附近的門源臺(tái)的重力增加3.40 μgal,能夠被重力儀器監(jiān)測(cè)到。

本文計(jì)算采用是李振洪等[20]的破裂模型,隨著資料的公開和共享,更加精確的震源模型會(huì)改變本文的計(jì)算結(jié)果,但本文基于InSAR的數(shù)據(jù)反演模型抓住了主要的特征,即使后續(xù)有更為精確的破裂模型,對(duì)我們的計(jì)算結(jié)果不會(huì)較大影響。

余震也會(huì)對(duì)本文的計(jì)算結(jié)果也會(huì)有一定的影響。萬永革等在計(jì)算蘭德斯地震余震產(chǎn)生的位移場(chǎng)時(shí)發(fā)現(xiàn),地震斷層面及其附近余震產(chǎn)生了厘米量級(jí)的位移變化[24],對(duì)計(jì)算地表重力變化的數(shù)值影響較小;并且目前沒有更小地震的震源機(jī)制和破裂模型數(shù)據(jù),所以本文中沒有計(jì)算更小地震對(duì)地表重力變化的影響。

除了上文所說的震源破裂模型、余震會(huì)對(duì)計(jì)算結(jié)果有影響,本文采取震源周圍介質(zhì)為均勻彈性半空間,密度統(tǒng)一為2 670 kg/m3,這也會(huì)對(duì)計(jì)算結(jié)果造成一定的影響。在計(jì)算門源地震后的重力變化時(shí),同樣還采用了密度為2 630 kg/m3、2 650 kg/m3、2 690 kg/m3進(jìn)行計(jì)算,但是最后的計(jì)算結(jié)果影響都比較小,重力變化特征也基本一致。在不同密度的計(jì)算結(jié)果中,也是只有一個(gè)S9(門源臺(tái))前兆臺(tái)站能夠監(jiān)測(cè)到本次門源地震的重力變化,其余臺(tái)站處重力變化比較小。同樣震后余滑或介質(zhì)黏彈性松弛等產(chǎn)生的影響也較小,所以可以忽略不計(jì)。

雖然上述種種問題還有待于今后進(jìn)一步解決,但本文在一級(jí)近似情形下給出了不同類型斷層錯(cuò)動(dòng)在地表產(chǎn)生的重力變化模式,并根據(jù)2022年門源地震的破裂模型計(jì)算得到了該地震產(chǎn)生的重力變化,這對(duì)采用重力資料反演該地震的破裂模型和研究長(zhǎng)期重力變化從而揭示其他動(dòng)力學(xué)過程是有益的。