利用GRACE數(shù)據(jù)檢測(cè)日本Mw9.0地震同震和震后重力變化

尹 鵬，張永志，焦佳爽，南康康，王慧玲

(長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院，陜西 西安 710054)

2011年3月11日在日本東北海域發(fā)生Mw9.0級(jí)地震，此次地震發(fā)生在太平洋板塊和歐亞板塊的邊界上。根據(jù)美國(guó)地質(zhì)勘查局(USGS)發(fā)布的地震監(jiān)測(cè)信息，該地震的震中位于東經(jīng)142.4°、北緯38.3°的日本海溝底部。

強(qiáng)震的孕育和發(fā)生必然引起震源區(qū)和外圍地區(qū)一定范圍內(nèi)地球重力場(chǎng)的變化。Sun & Okubo[1]提出大于Mw9.0的剪切型或大于Mw7.5的張裂型地震所產(chǎn)生的重力變化可以被GRACE觀測(cè)到，并得到證實(shí)[2-7]。2011年日本Mw9.0地震發(fā)生后，Matsuo & Heki[8]利用地震后兩個(gè)月的GRACE數(shù)據(jù)檢測(cè)到了-7 μgal重力負(fù)變化，同時(shí)發(fā)現(xiàn)GRACE檢測(cè)到的大地震重力負(fù)變化幅度與地震矩的大小存在近似線性關(guān)系；鄒正波等[9]采用不同的時(shí)間尺度，詳細(xì)分析研究了地震發(fā)生前后的重力動(dòng)態(tài)演化過(guò)程，反映出該區(qū)域的物質(zhì)運(yùn)動(dòng)和質(zhì)量遷移；張克亮等[10]采用最小二乘法和經(jīng)驗(yàn)正交函數(shù)方法，提取了日本Mw9.0地震震區(qū)及其周邊的同震變化特征，但未考慮陸地水儲(chǔ)量變化的影響，并且缺乏斷層模型正演結(jié)果的檢驗(yàn)。

前人已利用GRACE數(shù)據(jù)對(duì)日本地震同震和震后重力變化進(jìn)行研究，但震后數(shù)據(jù)的時(shí)段不長(zhǎng)，考慮到數(shù)據(jù)長(zhǎng)度及陸地水儲(chǔ)量的變化在同震和震后研究的重要性，本文利用GRACE數(shù)據(jù)，扣除陸地水儲(chǔ)量的影響，結(jié)合位錯(cuò)模型正演的同震重力變化，計(jì)算日本Mw9.0地震震中及其周邊地區(qū)的同震和震后5 a的重力場(chǎng)演化過(guò)程。

1 數(shù)據(jù)選取

1.1 GRACE衛(wèi)星月重力場(chǎng)模型數(shù)據(jù)

德國(guó)地學(xué)中心(GFZ)發(fā)布的RL05月重力場(chǎng)模型給出完全正則化的球諧系數(shù)，在解算過(guò)程中已經(jīng)扣除固體潮、極潮、大氣、海洋等因素的影響。該模型的最高階次為90階，大于60階的數(shù)據(jù)會(huì)引入較大的誤差，故將球諧系數(shù)截?cái)嘀?0階[11]，選取2009年3月至2016年3月的重力場(chǎng)模型進(jìn)行分析，由于GRACE衛(wèi)星進(jìn)入生命末期，存在數(shù)據(jù)缺失情況，同時(shí)地震發(fā)生時(shí)段(2011年3月)的數(shù)據(jù)也一并舍去不用。由于GRACE衛(wèi)星對(duì)于低階項(xiàng)不敏感，采用人工激光測(cè)距得到的二階項(xiàng)代替球諧系數(shù)中的C20項(xiàng)。

1.2 GLDAS陸地水?dāng)?shù)據(jù)

GLDAS(global land data assimilation system)是美國(guó)國(guó)家航空航天局(NASA)戈達(dá)德空間飛行中心(GSFC)發(fā)布的陸地水儲(chǔ)量模型。GLDAS陸地水儲(chǔ)量包括土壤水分和雪水變化，最大時(shí)間分辨率為3 h，最大空間分辨率為0.25°×0.25°。GRACE觀測(cè)的重力變化包含多種信號(hào)共同作用，土壤水分和雪水變化是除地震外最重要的影響因素[8]，本文選取1°×1°的月水儲(chǔ)量模型，提取陸地水含量引起的重力變化，并將其從GRACE月重力場(chǎng)模型中扣除，去除這一非構(gòu)造因素對(duì)同震和震后重力變化的影響。

2 重力信號(hào)提取

2.1 GRACE數(shù)據(jù)處理

由于GRACE衛(wèi)星觀測(cè)本身的限制及解算中采用的大氣和海洋等模型誤差因素的影響，GRACE時(shí)變重力場(chǎng)在圖像上表現(xiàn)出明顯的“南北條帶”誤差，本文采用300 km的扇形濾波處理?xiàng)l帶誤差，加入扇形濾波后橢球面上的GRACE時(shí)變重力變化可表達(dá)為[12]

(1)

提取同震重力變化信號(hào)時(shí)，為了削弱時(shí)變重力場(chǎng)的季節(jié)性影響，突出同震重力變化，選取日本Mw9.0地震前2009-11—2010-02、2010-11—2011-02的平均重力場(chǎng)為背景場(chǎng)，利用震后兩年相同月份的平均重力場(chǎng)與背景場(chǎng)差分，其中2011-01數(shù)據(jù)缺失，為了保證數(shù)據(jù)的對(duì)稱性，舍去2010年1月的數(shù)據(jù)。提取震后重力演化信號(hào)時(shí)，利用震后每年相同月份的平均重力場(chǎng)與背景重力場(chǎng)進(jìn)行差分，其中2014-02、2014-12、2015-11數(shù)據(jù)缺失，分別用2014—2016年1月的數(shù)據(jù)代替。

2.2 GLDAS數(shù)據(jù)處理

球面上的任意函數(shù)f(θ,λ)可以展開為球諧級(jí)數(shù)[13]，由此可將土壤水分和雪水含量轉(zhuǎn)化為與之對(duì)應(yīng)的球諧系數(shù)[14]

(2)

選取與GRACE數(shù)據(jù)相同月份的GLDAS數(shù)據(jù)計(jì)算陸地水儲(chǔ)量對(duì)同震重力的影響，如圖1所示，圖中黑色五角星表示震中。圖1(a)為日本地震震中周圍地區(qū)震后2年與震前2年1°×1°的平均等效水高之差，主要為負(fù)值。在黑龍江東部的中俄交界處、俄羅斯濱海邊疆區(qū)南部及日本本州東北部等地有比較明顯的水文效應(yīng)，俄羅斯濱海邊疆區(qū)南部的等效水高之差正負(fù)值均達(dá)到最大，分別為-115.2和59.8 mm；在日本的北海道、本州西南部、四國(guó)及九州地區(qū)的水文效應(yīng)為正值，日本本州東北部和北海道的最大等效水高分別為-44.3、23.2 mm。

將地震前后的等效水高按式(2)轉(zhuǎn)化為球諧系數(shù)，為了達(dá)到和GRACE一致的空間分辨率，將展開后的球諧系數(shù)截?cái)嘀?0階，并加以300 km的扇形濾波，獲得地震前后兩年1°×1°的等效水高之差，如圖1(b)所示?？臻g平滑后，部分陸地信號(hào)泄露到海洋地區(qū)，空間分辨率降低，等效水高振幅減小，由原來(lái)的-115.2 mm～59.8 mm，降至-8.5 mm～1.6 mm，在俄羅斯濱海邊疆區(qū)南部、日本本州東北部的負(fù)水文效應(yīng)比較顯著，分別為-8.5和-4.2 mm，與圖1(a)的負(fù)水文效應(yīng)空間對(duì)應(yīng)關(guān)系較好；本州東南呈現(xiàn)正的水文效應(yīng)，最大值為1.6 mm。

圖1 日本Mw9.0地震前后2年土壤水分和雪水變化及其引起的重力變化

將式(2)得到的球諧系數(shù)截?cái)酁?0階后代入式(1)，即可求得土壤水分和雪水變化引起的重力變化，如圖1(c)所示，在空間分布上與圖1(b)基本對(duì)應(yīng)，空間平滑后，9 mm的等效水高變化大致與0.4 μgal的重力變化相對(duì)應(yīng)?？傮w上看，土壤水含量和雪水的變化引起的重力變化約為-0.34～0.09 μgal。

2.3 同震和震后重力變化

按照前文所述方法，在不考慮陸地水儲(chǔ)量效應(yīng)時(shí)獲得的日本Mw9.0地震震中周圍地區(qū)的同震重力變化，如圖2(a)所示，黑色五角星代表震中，黑色方框?yàn)閃ei給出的斷層模型邊界[15]。由提取的重力信號(hào)可知，在斷層破裂面兩側(cè)重力呈正-負(fù)變化分布，量值范圍在-5.4～4.1 μgal，震中北方向出現(xiàn)了以(140.8°E，45.4°N)為中心，量值最高達(dá)3.7 μgal的正重力異常極值區(qū)。

圖2(b)所示為扣除土壤水分和雪水影響之后重力變化的空間分布，重力變化幅度為-5.3～4.2 μgal，北海道北部(140.8°E，45.4°N)的正重力異常極大值仍然有3.9 μgal，由于采用疊積法計(jì)算的同震重力變化包含震后重力變化信息，此處出現(xiàn)的重力異常可能與震后殼幔物質(zhì)調(diào)整有關(guān)[6]。圖2(b)在空間分布上與圖2(a)基本一致，但有細(xì)微差別，北海道附近的重力極大值區(qū)擴(kuò)張，日本本州西部日本海域的負(fù)極值區(qū)收縮，雖然水文效應(yīng)對(duì)日本Mw9.0地震同震重力的影響不明顯，但在汶川Ms8.0地震[16]、尼泊爾M8.1地震[17]等同震重力變化較小的情況下，提取同震重力變化信號(hào)需要考慮水文效應(yīng)的影響。

圖2 日本Mw9.0地震的同震重力變化

圖3所示為日本Mw9.0地震震后重力時(shí)空演化過(guò)程，均扣除了土壤水分和雪水的影響。(a)—(e)分別為震后1～5 a相對(duì)于背景重力場(chǎng)的重力變化量，日本Mw9.0地震震后震中及其周邊地區(qū)的重力變化基本與同震重力變化特征基本一致，震后5 a重力整體上呈增加趨勢(shì)，斷層上盤所在的日本海域重力增加的最大值為2.6 μgal，斷層下盤所在的太平洋區(qū)域重力增加較明顯，重力增加的峰值為4.5 μgal，比同震重力變化大0.3 μgal。

2.4 特征點(diǎn)時(shí)序重力變化

為了反映震后重力的詳細(xì)變化過(guò)程，在斷層上下盤各取一個(gè)重力變化極值點(diǎn)A(136.2°E，38.9°N)、B(143.4°E，35.3°N)，如圖2(b)中所示的白色三角形。以2009-03—2011-02共23個(gè)月(2011-01數(shù)據(jù)缺失)的GRACE數(shù)據(jù)均值為背景場(chǎng)，分別計(jì)算2009-03—2016-03各月重力場(chǎng)相對(duì)于背景場(chǎng)的差值，并加以300 km的扇形濾波，對(duì)于缺失月份的重力場(chǎng)數(shù)據(jù)，取其前后月重力場(chǎng)模型的均值來(lái)代替，獲取A、B兩點(diǎn)相對(duì)于背景場(chǎng)的月重力變化時(shí)間序列并采用最小二乘擬合的方法扣除時(shí)序重力變化中周期性水文及S2潮汐波的影響[18]

(3)

式中，A為常數(shù)項(xiàng)；B為年變率；Δt為相對(duì)參考?xì)v元的時(shí)間差；wi為振幅Ci和Di的周期項(xiàng)角頻率；i=1,2,3分別表示周年、半周年項(xiàng)及與S2潮汐波相關(guān)的161 d混頻項(xiàng)；ε為殘差。

圖3 日本Mw9.0地震震后重力演化過(guò)程

圖4(a)、(b)分別表示A、B點(diǎn)的月時(shí)序重力變化，方塊曲線表示僅進(jìn)行扇形濾波的時(shí)間序列，圓圈曲線表示最小二乘擬合的時(shí)間序列，黑色虛線表示地震發(fā)生的2011年3月，黑色五角星表示地震發(fā)生當(dāng)月的重力場(chǎng)與背景重力場(chǎng)的差值，黑色實(shí)線表示地震前后扣除周期性信號(hào)之后重力的線性變化趨勢(shì)。由圖4可以看出，A、B兩點(diǎn)都有明顯的同震階躍，A、B兩點(diǎn)的年變率見表1，震前A點(diǎn)重力增加，B點(diǎn)基本保持不變；震后均呈增加趨勢(shì)，B點(diǎn)的重力年變率為0.77 μgal/a，重力增加較A快，可能是因?yàn)锽點(diǎn)所在的太平洋地層黏滯性相對(duì)較低。

表1 日本Mw9.0地震前后A、B兩點(diǎn)重力的年變率

3 日本Mw9.0地震理論同震重力變化

為驗(yàn)證GRACE衛(wèi)星測(cè)得同震信號(hào)的可靠性，采用分層黏彈性半空間位錯(cuò)模型PSGRN/PSCMP[19]及Wei發(fā)布的斷層模型[15]計(jì)算2011年日本Mw9.0地震的同震重力變化，計(jì)算過(guò)程中采用5層分布的地球模型，見表2。PSGRN/PSCMP位錯(cuò)模型計(jì)算的是變形地表的重力變化，與地面重力測(cè)量結(jié)果對(duì)應(yīng)，為了與GRACE衛(wèi)星的空間觀測(cè)結(jié)果保持一致，需要在位錯(cuò)模型計(jì)算結(jié)果的基礎(chǔ)上加入地表垂直形變引起的空間改正[20]，由于此次地震發(fā)生在海域，還需加入海水質(zhì)量改正[21],將改正后的結(jié)果加以300 km的扇形濾波，以達(dá)到和GRACE一致的空間分辨率。理論同震重力變化空間分布如圖5所示。

圖4 日本Mw9.0地震前后A、B兩點(diǎn)的時(shí)序重力變化

圖5 日本Mw9.0地震理論同震重力變化

由圖5可知，理論同震重力變化與利用GRACE重力場(chǎng)模型計(jì)算的同震重力變化在量級(jí)和空間分布上具有較好的一致性。斷層上盤的負(fù)重力異常，理論結(jié)果為-6.6 μgal，GRACE衛(wèi)星觀測(cè)結(jié)果為-5.3 μgal；斷層下盤的正重力異常，理論結(jié)果為3.7 μgal，GRACE衛(wèi)星觀測(cè)結(jié)果為4.2 μgal。理論重力變化與GRACE觀測(cè)結(jié)果存在的差別可能是因?yàn)榀B積法計(jì)算的同震重力變化中包含震后重力變化，同震和震后信號(hào)未分離。

4 結(jié) 語(yǔ)

日本Mw9.0地震前后兩年，震中及周邊地區(qū)土壤水分和雪水的平均等效水高變化為-8.5～1.6 mm，引起的重力變化為-0.34～0.09 μgal，9 mm的等效水高變化與0.4 μgal的重力變化相對(duì)應(yīng)，雖然水文效引起的重力變化較小，但扣除水文效應(yīng)能更準(zhǔn)確地反映同震重力變化?？鄢男?yīng)的影響之后，同震重力變化的幅度為-5.3～4.2 μgal，PSGRN/PSCMP模型計(jì)算的同震重力變化為-6.6～3.7 μgal，二者在空間分布上的對(duì)應(yīng)關(guān)系比較好，在數(shù)值上存在細(xì)微差別的原因是沒(méi)有將震后重力變化與同震重力變化分離。

表2日本Mw9.0地震5層分布的地球模型(各層參數(shù)來(lái)自于CRUST1.0)

深度/km密度/(kg·m-3)P波速度/(km·s-1)S波速度/(km·s-1)0～2．0518201．740．332．05～6．0427206．003．506．04～12．0328606．603．8012．30～22．0030507．103．9022．00～∞32607．904．00

震后重力變化與同震重力變化特征基本保持一致，震后5 a震中附近區(qū)域的重力整體上呈現(xiàn)增加趨勢(shì)，斷層上盤所在的日本海域重力增加的最大值為2.6 μgal，下盤所在的太平洋區(qū)域重力增加的最大值為4.5 μgal，比同震重力大0.3 μgal，可能與斷層下盤所在地層的黏滯性相對(duì)較低有關(guān)。利用最小二乘擬合的方法計(jì)算了兩特征點(diǎn)在地震前后85個(gè)月的重力年變率，震后兩個(gè)特征點(diǎn)的重力年變率分別為0.19與0.77 μgal/a，重力變化均呈現(xiàn)增加趨勢(shì)。

在利用疊積法計(jì)算日本Mw9.0地震同震和震后重力變化時(shí)，采用了300 km的扇形濾波，較好地消除了南北條帶誤差，但仍然存在不太明顯的條帶，可以嘗試去相關(guān)濾波和扇形濾波相結(jié)合的方法。

[1] SUN W，OKUBO S.Coseismic Deformations Detectable by Satellite Gravity Missions：A Case Study of Alaska(1964，2002)and Hokkaido(2003)Earthquakes in the Spectral Domain[J].Journal of Geophysical Research Solid Earth，2004，109(B4):327-341.

[2] HAN S C，SHUM C K，BEVIS M，et al.Crustal Dilatation Observed by GRACE after the 2004 Sumatra-Andaman Earthquake[J].Science，2006，313(5787):658.

[3] CHEN J L，WILSON C R，TAPLEY B D，et al.GRACE Detects Coseismic and Postseismic Deformation from the Sumatra-Andaman Earthquake[J].Geophysical Research Letters，2007，34(13):173-180.

[4] RYOKO O，KOSUKE H.Slow Postseismic Recovery of Geoid Depression Formed by the 2004 Sumatra-Andaman Earthquake by Mantle Water Diffusion[J].Geophysical Research Letters，2007，34(6):160-166.

[5] ISABELLE P，VALENTIN M，MICHEL D，et al.Coseismic and Post-Seismic Signatures of the Sumatra 2004 December and 2005 March Earthquakes in GRACE Satellite Gravity[J].Geophysical Journal of the Royal Astronomical Society，2007，171(1):177-190.

[6] DE LINAGE C，RIVERA L，HINDERER J，et al.Separation of Coseismic and Postseismic Gravity Changes for the 2004 Sumatra-Andaman Earthquake from 4.6 yr of GRACE Observations and Modelling of the Coseismic Change by Normal-Modes Summation[J].Translated World Seismology，2010，176(3):695-714.

[7] 周新，孫文科，付廣裕.重力衛(wèi)星GRACE檢測(cè)出2010年智利Mw8.8地震的同震重力變化[J].地球物理學(xué)報(bào)，2011，54(7):1745-1749.

[8] KOJI M，KOSUKE H.Coseismic Gravity Changes of the 2011 Tohoku-Oki Earthquake from Satellite Gravimetry[J].Geophysical Research Letters，2011，38(7):113-120.

[9] 鄒正波，羅志才，吳海波,等.日本MW9.0地震前GRACE衛(wèi)星重力變化[J].測(cè)繪學(xué)報(bào)，2012，41(2):171-176.

[10] 張克亮，甘衛(wèi)軍，周新.GRACE衛(wèi)星重力場(chǎng)同震變化的經(jīng)驗(yàn)正交函數(shù)分解:以日本MW9.0地震為例[J].地震地質(zhì)，2014，36(3):763-774.

[11] CHEN J L，RODELL M，WILSON C R，et al.Low Degree Spherical Harmonic Influences on Gravity Recovery and Climate Experiment (GRACE) Water Storage Estimates[J].Geophysical Research Letters，2005，32(14):57-76.

[12] ZHANG Z，CHAO B F，YANG L，et al.An Effective Filtering for GRACE Time-Variable Gravity：Fan Filter[J].Geophysical Research Letters，2009，36(17):1397-1413.

[13] HEISKANEN W A，MORITZ H.Physical Geodesy[J].Bulletin Géodésique (1946-1975)，1967，86(1):491-492.

[14] WAHR J，MOLENAAR M，BRYAN F.Time Variability of the Earth’s Gravity Field：Hydrological and Oceanic Effects and Their Possible Detection Using GRACE[J].Journal of Geophysical Research Solid Earth，1998，103(B12):30205-30230.

[15] WEI S，SLADEN A，ARIA group. Updated Result[EB/OL].2011-11-03. http:∥www.tectonics.caltech.edu/slip_history/2011_taiheiyo-oki/.

[16] 李進(jìn).基于GRACE時(shí)變重力梯度探測(cè)大地震的同震變化[D].武漢：武漢大學(xué)，2011.

[17] 鄒正波，李輝，吳云龍,等.尼泊爾M8.1地震震前衛(wèi)星重力場(chǎng)時(shí)變特征[J].大地測(cè)量與地球動(dòng)力學(xué)，2015，35(4):547-551.

[18] STEFFEN H，GITLEIN O，DENKER H，et al.Present Rate of Uplift in Fennoscandia from GRACE and Absolute Gravimetry[J].Tectonophysics，2009，474(1-2):69-77.

[19] WANG R，MARTN F，ROTH F.PSGRN/PSCMP—A New Code for Calculating Co- and Post-seismic Deformation，Geoid and Gravity Changes Based on the Viscoelastic-Gravitational Dislocation Theory[J].Computers & Geosciences，2006，32(4):527-541.

[20] SUN W，OKUBO S，F(xiàn)U G，et al.General Formulations of Global Co-seismic Deformations Caused by an Arbitrary Dislocation in a Spherically Symmetric Earth Model-Applicable to Deformed Earth Surface and Space-fixed Point[J].Geophysical Journal International，2009，177(3):817-833.

[21] HEKI K，MATSUO K.Coseismic Gravity Changes of the 2010 Earthquake in Central Chile from Satellite Gravimetry[J].Geophysical Research Letters，2010，37(24):701-719.