韋 進(jìn) 郝洪濤 韓宇飛 胡敏章 江 穎 劉子維

1)中國地震局地震研究所，武漢 430071 2)中國地震局地震大地測量重點(diǎn)實(shí)驗(yàn)室，武漢 430071 3)中國地震臺網(wǎng)中心，北京 100045 4)湖北省地震局，武漢 430071

0 引言

地面連續(xù)重力臺能記錄到環(huán)境變化、 地表形變和內(nèi)部物質(zhì)運(yùn)移產(chǎn)生的重力變化。 該重力變化通常由2部分組成： 重力臺下方的物質(zhì)再分布產(chǎn)生的密度變化和地表垂直位移變化(Xingetal.，2021)。 大地震會導(dǎo)致震中所在斷層沿線發(fā)生破裂，而破裂的斷層兩側(cè)會出現(xiàn)大幅形變，斷層所在地塊的地下物質(zhì)分布也會發(fā)生變化(Vanetal.，2017)。 這些變化會在地表形成永久的形變(同震形變)和重力變化(同震重力變化)。 其中同震重力變化的振幅主要受地震的震級和震中距影響(Prasadetal.，2017)。 該物理現(xiàn)象已多次被觀測精度優(yōu)于1×10-8m·s-2的遠(yuǎn)場超導(dǎo)重力儀捕捉到(Imanishietal.，2004，2009)。 此外，近年來一些彈簧重力儀也捕捉到同震重力信號(Weietal.，2013； Prasadetal.，2017)。

1 數(shù)據(jù)與處理方法

gH(ti-Tg)+hH(ti-Th)ti+kexp[-(ti-Tk)/τ]H(ti-Tk)ti+vi

(1)

式(1)中，gobs(ti)為重力固體潮分鐘值觀測數(shù)據(jù)，ti為間隔為分鐘的歷元時刻，i取為1，…，N；d0和d1為重力殘差值的常數(shù)項(xiàng)和零漂項(xiàng)；f為觀測站平均潮汐因子，gtheory為基于天頂距的重力固體潮計(jì)算值； Δp為同址氣壓觀測數(shù)據(jù)變化值，p為氣壓導(dǎo)納值；aj和bj為重力殘差中由于彈簧重力儀非線性零漂產(chǎn)生的周日和半日周期的零漂項(xiàng)，由于只考慮周日和半日波信號的非線性項(xiàng)，故j最大為2。Tj分別為1440和2880；H為階躍函數(shù)，Tg為階躍歷元，Th為速率改變歷元，Tk為開始弛豫變化歷元，g為階躍值，h為重力速率改變后的速率，k為弛豫變化的振幅；τ為弛豫時間；vi為誤差項(xiàng)。

本文所分析數(shù)據(jù)的采樣率均為分鐘/次，1min內(nèi)P波可傳播的距離為330～420km。 震中距約800km范圍內(nèi)的重力臺在地震發(fā)生后1～2min內(nèi)可記錄到地震波，因此設(shè)定階躍歷元Tg，速率改變歷元Th和開始弛豫變化歷元Tk均為瑪多MS7.4地震的發(fā)震時刻歷元。 式(1)中的階躍值g就是重力臺觀測的同震重力變化值。

2 觀測結(jié)果和數(shù)據(jù)分析

2.1 重力臺觀測的同震重力變化

本文利用式(1)計(jì)算了圖 1 距離震中800km內(nèi)的5個重力臺的同震重力變化值。 圖 2 為震中距最近的瑪沁重力臺潮汐改正和同震重力變化結(jié)果。 利用弛豫分析方法(簡稱指數(shù)函數(shù)法)計(jì)算出瑪沁臺同震重力變化為(2.95±0.70)×10-8m·s-2(圖2c中的藍(lán)色線條)。

圖2 瑪沁重力臺重力改正和同震重力變化Fig. 2 The gravity correction and the co-seismic gravity change in Maqin gravity station.

圖 3 為利用式(1)計(jì)算的其余4個遠(yuǎn)場重力臺的同震重力變化結(jié)果。 為比較地震后重力臺是否包含弛豫變化過程，本文基于式(1)也給出了一組不考慮指數(shù)項(xiàng)只考慮階躍項(xiàng)(簡稱階躍函數(shù)法)的結(jié)果(圖2c和圖3a—d中的紅色線條)。 結(jié)果表明，兩者的同震重力變化量級差異不大，僅為(0.1～1)×10-8m·s-2。 5個臺的最優(yōu)弛豫時間也僅11min，和蘆山地震郫縣臺GS15重力儀記錄的約2000min的結(jié)果(Weietal.，2013)相比，前者的地震弛豫過程非常短，且林芝臺階躍函數(shù)法估計(jì)的結(jié)果和模擬的一致性更好(表2)。 由此可見，和指數(shù)函數(shù)法相比，階躍函數(shù)法的估計(jì)結(jié)果更符合本次地震對gPhone重力儀的震后影響特征。

圖3 4個遠(yuǎn)場重力臺的同震重力變化結(jié)果Fig. 3 The co-seismic gravity change observed by the 4 far-field gravity stations.

2.2 同震重力變化模擬與比較

USGS在協(xié)調(diào)世界時2021年5月26日16時40分23秒給出了利用24個遠(yuǎn)場寬頻帶P波波形、 28個寬頻帶SH波波形以及79個長周期面波波形反演的本地震的斷層滑動模型。 該地震的破裂長度約為170km，寬25km，子斷層的最大破裂長3.5km(2)https:∥earthquake.usgs.gov/earthquakes/eventpage/us7000e54r/finite-fault。。 本次地震發(fā)生在昆侖山口-江錯斷裂的延伸斷層附近(華俊等，2021)。 為檢驗(yàn)和研究觀測和理論模擬結(jié)果的正確性和相關(guān)性，本文利用Okubo的無限半空間模型(Okubo，1992)和孫文科的地球分層模型(Sunetal.，2006，2009)(簡稱Sun)分別正演了研究區(qū)域的同震重力變化空間分布(圖 4)，并比較了5個重力臺位置的觀測和模型值，結(jié)果見表1。

圖4 基于Okubo和Sun同震位錯模型的模擬結(jié)果比較Fig. 4 The simulation of co-seismic gravity change by the dislocation model of Okubo and Sun.

表1 gPhone重力儀位置的同震重力變化和模擬結(jié)果的比較Table1 Comparison of the observed and simulated co-seismic gravity change

基于Okubo和Sun位錯模型模擬的研究區(qū)同震重力變化空間分布特征(圖 4)，在地震發(fā)震斷層兩側(cè)的近場附近有(-15～20)×10-8m·s-2的北正南負(fù)、 相互平行的同震重力變化異常區(qū)域，在遠(yuǎn)場變化較小。 2個模型的模擬結(jié)果都表現(xiàn)為四象限分布特征。 和5個重力臺觀測的同震重力變化分布特征比較(表1)認(rèn)為，遠(yuǎn)場結(jié)果和Sun位錯模型的模擬結(jié)果一致性更好。

張喆等(2021)研究認(rèn)為瑪多MS7.4地震是一次走向?yàn)?81°且斷層面近直立的走滑型地震。 高臺和林芝臺位于發(fā)震斷層S、 N兩側(cè)，瑪沁和松潘位于發(fā)震斷層走向的E側(cè)，中甸臺位于震中的SE側(cè)(圖 4)。

在稍近區(qū)域，發(fā)震斷層走向E側(cè)的2個重力臺(瑪沁、 松潘)不僅和震中同處巴顏喀拉塊體，而且其同震重力變化(≥2.9×10-8m·s-2)遠(yuǎn)大于垂直走向的重力臺結(jié)果。 其中，震中距僅175km的瑪沁臺的觀測值((2.95～4.00)×10-8m·s-2)和模擬結(jié)果((0.21～0.26)×10-8m·s-2)雖相差約1個量級，但方向一致。

走滑型地震形成的斷層破裂在理論上會沿破裂的走向延伸。 該破裂過程不僅會在近場區(qū)域產(chǎn)生形變導(dǎo)致大幅同震重力變化，還會在沿斷層面走向的遠(yuǎn)場區(qū)域產(chǎn)生比垂直走向方向更大的同震重力變化。 瑪多MS7.4地震的發(fā)震斷層面近直立的走滑特征會使得上述規(guī)律更加明顯。 在如圖4a、 b所示的 2個模型的模擬結(jié)果中，地震遠(yuǎn)場區(qū)域瑪沁臺SW側(cè)有明顯的重力異常增大的現(xiàn)象。 從觀測結(jié)果看，位于斷層面走向方向的重力臺的同震重力變化結(jié)果比垂直于走向方向的重力臺結(jié)果高出1個數(shù)量級。 同時，無論是觀測值還是模型模擬都符合本次左旋走滑地震的運(yùn)動特征。

3 結(jié)果和討論

3.1 重力潮汐觀測和同震重力變化精度討論

表2 氣壓導(dǎo)納值和平均潮汐因子的統(tǒng)計(jì)比較Table2 Statistical comparison of barometric admittance and mean tidal factor

3.2 瑪沁和松潘臺同址GNSS和重力同震觀測結(jié)果的比較

圖5 瑪沁臺GNSS和重力同震觀測結(jié)果的比較Fig. 5 Comparison of the co-seismic vertical displacement and gravity change in Maqin station.

瑪多MS7.4地震，瑪沁臺觀測的同震重力變化為(2～4)×10-8m·s-2，松潘臺則為9.13×10-8m·s-2。 雖然瑪沁臺的GNSS和重力觀測的同震重力變化量級相當(dāng)，但松潘臺的同震重力變化比GNSS觀測以及模型模擬的結(jié)果都大，且松潘和瑪沁臺同震重力變化分布特征也不符合震中距越大同震重力變化越小的規(guī)律。

3.3 九寨溝MS7.0和瑪多MS7.4地震的同震重力變化模擬和比較

對于遠(yuǎn)場而言，尼泊爾MW7.8地震時麗江(7.3×10-8m·s-2)和武漢臺的結(jié)果(8.0×10-8m·s-2)(Xingetal.，2021)與松潘臺的同震重力變化結(jié)果(9.13×10-8m·s-2)相似。 而就震級和震中位置而言，瑪多MS7.4地震和2017年8月8日的九寨溝MS7.0地震更相似。 為研究松潘臺同震重力變化特征是否還受到其他因素影響，本文收集了王衛(wèi)民(3)http:∥www.itpcas.ac.cn/xwzx/zhxw/201708/t20170809_4840737.html。和鄭緒君等(2017)給出的九寨溝MS7.0地震的滑動斷層模型結(jié)果，分別利用Okubo和Sun位錯模型模擬九寨溝MS7.0地震周邊區(qū)域的同震重力變化空間分布情況，并與重力臺的觀測結(jié)果進(jìn)行了比較(圖 6)。

圖6 九寨溝MS7.0地震的同震重力模擬和觀測結(jié)果Fig. 6 The simulated and observed co-seismic gravity change in the Jiuzhaigou MS7.3 earthquake.

3.4 漾濞MS6.4地震對中甸臺瑪多MS7.4地震同震重力變化估計(jì)結(jié)果的影響討論

Xing等(2021)認(rèn)為在2015年尼泊爾MW7.8地震中，超導(dǎo)重力儀記錄到的同震信號中有約0.1×10-8m·s-2的地下水位波動產(chǎn)生的影響，這表明同震重力變化會受到環(huán)境因素的影響。 而在瑪多MS7.4地震發(fā)生前5h出現(xiàn)的最顯著的環(huán)境變化是漾濞MS6.4地震。 在5個臺的重力時序數(shù)據(jù)中，只有中甸臺在漾濞MS6.4地震發(fā)生時產(chǎn)生了明顯的同震重力信號(圖3c)。 為消除漾濞MS6.4地震對瑪多MS7.4同震重力變化估計(jì)的影響，本文在進(jìn)行同震重力變化計(jì)算時還加入了漾濞地震參數(shù)，用2種方法再次估計(jì)中甸臺2次地震同震重力變化，結(jié)果如圖 7 所示。

圖7 中甸重力臺記錄的漾濞MS6.4和瑪多MS7.4地震的同震重力變化結(jié)果Fig. 7 The co-seismic gravity changes by Yangbi MS6.4 and Maduo MS7.4 in Zhongdian station.

在利用指數(shù)函數(shù)估計(jì)同震重力變化時，漾濞MS6.4地震的弛豫時間為100min，瑪多MS7.4地震的弛豫時間為20min。 利用該方法得到的漾濞地震的估計(jì)結(jié)果為(1.31±0.98)×10-8m·s-2，瑪多地震為(1.09±0.97)×10-8m·s-2； 而利用階躍函數(shù)法估計(jì)的2個地震的同震重力變化結(jié)果分別為(0.53±0.77)×10-8m·s-2和(1.63±0.76)×10-8m·s-2。 階躍函數(shù)計(jì)算的漾濞地震結(jié)果誤差大于觀測值，而指數(shù)函數(shù)法則未出現(xiàn)該現(xiàn)象，同時指數(shù)函數(shù)計(jì)算的結(jié)果更接近Sun模型的模擬值。 因此中甸臺記錄的瑪多MS7.4地震的同震重力變化應(yīng)為(1.09±0.97)×10-8m·s-2，和未改正前的(1.78±0.24)×10-8m·s-2(表1)相比更接近理論值(0.08×10-8m·s-2)。

3.5 松潘臺的同震重力變化分析及和理論模型結(jié)果差異的討論

本文在分析松潘臺的同震重力變化時發(fā)現(xiàn)，無論是瑪多MS7.4地震(9.13×10-8m·s-2)還是九寨溝MS7.0地震(11.28×10-8m·s-2)都產(chǎn)生了大于其他遠(yuǎn)場重力臺的結(jié)果。 從精度分析來看，該臺地震前后的觀測值都符合潮汐的理論模型。 從圖 1 可知，該臺處于巴顏喀拉塊體的東邊界、 岷江斷裂附近，屬于塊體邊界區(qū)域，而其他重力臺多處于塊體內(nèi)部。 歷史上該塊體邊界及其鄰區(qū)發(fā)生過多次MS>7.0地震(松潘-平武2次MS7.2地震和汶川M8.0地震)，表明該邊界區(qū)域?qū)儆谖镔|(zhì)積累和應(yīng)力釋放活躍區(qū)域。 重力觀測結(jié)果也顯示巴顏喀拉塊體東邊界的E、 W兩側(cè)一直是重力差異變化較為明顯的區(qū)域(祝意青等，2017)。 發(fā)生在巴顏喀拉塊體內(nèi)部的瑪多MS7.4強(qiáng)震不僅加快了塊體內(nèi)部近場區(qū)域的物質(zhì)積累和調(diào)整(圖 4)，同時也會觸發(fā)本就構(gòu)造活動背景極強(qiáng)的巴顏喀拉塊體東邊界區(qū)域的物質(zhì)運(yùn)移和調(diào)整，進(jìn)而出現(xiàn)塊體邊界重力臺的同震重力變化大于震中距較小的塊體內(nèi)部重力臺(瑪沁臺)觀測結(jié)果的現(xiàn)象。 這一推測和蘆山MS7.0地震中姑咱臺的結(jié)果情況相似(鄧明莉等，2014)。

此外，GNSS僅能測量地表垂直位移產(chǎn)生的重力變化，而重力觀測對物質(zhì)積累和密度變化更加敏感。 瑪多MS7.4地震除對松潘臺產(chǎn)生了垂直位移的影響外，應(yīng)該還包含該塊體邊界處不斷積累的荷載影響。 這一假設(shè)和Koyna和Warna關(guān)于水庫的負(fù)荷加載效應(yīng)會產(chǎn)生大于同震模擬結(jié)果的情況相似(Prasadetal.，2017)。

對于同震重力觀測變化大于位錯模型模擬結(jié)果的情況，一方面認(rèn)為模擬結(jié)果僅反映了地表形變垂直方向變化產(chǎn)生的重力變化，不能完全反映物質(zhì)運(yùn)移和密度變化帶來的影響，另一方面認(rèn)為位錯模型的模擬是基于一個有限單元的正演結(jié)果，反映該有限單元的平均特征，而地面重力臺觀測結(jié)果更突出某個斷層或斷裂的局部信息。 此外，同震位錯模型在進(jìn)行遠(yuǎn)場信號模擬時并不包含遠(yuǎn)場區(qū)塊邊界的地震活動信息和構(gòu)造背景信息，因此重力觀測和位錯理論模擬結(jié)果會存在差異。

4 結(jié)論

本文對瑪多MS7.4地震震中800km范圍內(nèi)地震前后共15d正常運(yùn)行的5個連續(xù)重力臺同震重力變化進(jìn)行分析，震中距175km的瑪沁臺記錄到了瑪多地震產(chǎn)生的(2.9～4.0)×10-8m·s-2同震重力變化信號。 通過與GNSS觀測及九寨溝地震同震重力變化空間分布進(jìn)行比較后認(rèn)為，松潘臺在瑪多MS7.4地震時刻產(chǎn)生的約9.13×10-8m·s-2的信號中除同震重力信號外，還包含有其他因素的影響。 中甸臺在消除漾濞MS6.4地震同震影響后，同震重力變化信號約為1.09×10-8m·s-2。 垂直于發(fā)震斷層走向的林芝臺和高臺臺的觀測結(jié)果和誤差相當(dāng)(0.5×10-8m·s-2)，但信號方向和Sun模型模擬結(jié)果的一致性更高。

本文的重力臺觀測結(jié)果均為遠(yuǎn)場數(shù)據(jù)。 觀測和模擬結(jié)果的方向一致，表明遠(yuǎn)場重力臺觀測到的同震重力變化信號中應(yīng)該包含地震對遠(yuǎn)場地表形變引起的重力變化成分。 而觀測和模擬量級巨大的差異除與觀測系統(tǒng)自身的精度有關(guān)外(高臺臺、 林芝臺)，可能還與重力臺所處區(qū)域的其他因素有關(guān)(如松潘臺等)，可能包括背景構(gòu)造特征。 從本次瑪多MS7.4地震的數(shù)據(jù)來看，震中距175～800km范圍內(nèi)的gPhone重力儀能夠觀測到約(0.5～4)×10-8m·s-2的同震重力變化信號。

致謝國家重力臺網(wǎng)中心提供了瑪多地震前后gPhone連續(xù)重力臺的觀測數(shù)據(jù)； 湖北省地震局GNSS數(shù)據(jù)共享平臺提供了瑪沁和松潘臺的日解算結(jié)果； 王衛(wèi)民和張勇教授提供了九寨溝MS7.0地震的斷層滑動模型數(shù)據(jù)； 審稿專家為本文提出了寶貴意見和建議。 在此一并表示感謝!