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

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

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

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

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

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

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

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

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

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

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

2 重力信號提取

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

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

(1)

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

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

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

(2)

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

將地震前后的等效水高按式(2)轉(zhuǎn)化為球諧系數(shù)，為了達(dá)到和GRACE一致的空間分辨率，將展開后的球諧系數(shù)截?cái)嘀?0階，并加以300 km的扇形濾波，獲得地震前后兩年1°×1°的等效水高之差，如圖1(b)所示。空間平滑后，部分陸地信號泄露到海洋地區(qū)，空間分辨率降低，等效水高振幅減小，由原來的-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)空間對應(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)基本對應(yīng)，空間平滑后，9 mm的等效水高變化大致與0.4 μgal的重力變化相對應(yīng)。總體上看，土壤水含量和雪水的變化引起的重力變化約為-0.34～0.09 μgal。

2.3 同震和震后重力變化

按照前文所述方法，在不考慮陸地水儲量效應(yīng)時(shí)獲得的日本Mw9.0地震震中周圍地區(qū)的同震重力變化，如圖2(a)所示，黑色五角星代表震中，黑色方框?yàn)閃ei給出的斷層模型邊界[15]。由提取的重力信號可知，在斷層破裂面兩側(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)對日本Mw9.0地震同震重力的影響不明顯，但在汶川Ms8.0地震[16]、尼泊爾M8.1地震[17]等同震重力變化較小的情況下，提取同震重力變化信號需要考慮水文效應(yīng)的影響。

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

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

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

為了反映震后重力的詳細(xì)變化過程，在斷層上下盤各取一個(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ù)均值為背景場，分別計(jì)算2009-03—2016-03各月重力場相對于背景場的差值，并加以300 km的扇形濾波，對于缺失月份的重力場數(shù)據(jù)，取其前后月重力場模型的均值來代替，獲取A、B兩點(diǎn)相對于背景場的月重力變化時(shí)間序列并采用最小二乘擬合的方法扣除時(shí)序重力變化中周期性水文及S2潮汐波的影響[18]

(3)

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

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

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

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

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

為驗(yàn)證GRACE衛(wèi)星測得同震信號的可靠性，采用分層黏彈性半空間位錯(cuò)模型PSGRN/PSCMP[19]及Wei發(fā)布的斷層模型[15]計(jì)算2011年日本Mw9.0地震的同震重力變化，計(jì)算過程中采用5層分布的地球模型，見表2。PSGRN/PSCMP位錯(cuò)模型計(jì)算的是變形地表的重力變化，與地面重力測量結(jié)果對應(yīng)，為了與GRACE衛(wèi)星的空間觀測結(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重力場模型計(jì)算的同震重力變化在量級和空間分布上具有較好的一致性。斷層上盤的負(fù)重力異常，理論結(jié)果為-6.6 μgal，GRACE衛(wèi)星觀測結(jié)果為-5.3 μgal；斷層下盤的正重力異常，理論結(jié)果為3.7 μgal，GRACE衛(wèi)星觀測結(jié)果為4.2 μgal。理論重力變化與GRACE觀測結(jié)果存在的差別可能是因?yàn)榀B積法計(jì)算的同震重力變化中包含震后重力變化，同震和震后信號未分離。

4 結(jié) 語

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

表2日本Mw9.0地震5層分布的地球模型(各層參數(shù)來自于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)增加趨勢，斷層上盤所在的日本海域重力增加的最大值為2.6 μgal，下盤所在的太平洋區(qū)域重力增加的最大值為4.5 μgal，比同震重力大0.3 μgal，可能與斷層下盤所在地層的黏滯性相對較低有關(guān)。利用最小二乘擬合的方法計(jì)算了兩特征點(diǎn)在地震前后85個(gè)月的重力年變率，震后兩個(gè)特征點(diǎn)的重力年變率分別為0.19與0.77 μgal/a，重力變化均呈現(xiàn)增加趨勢。

在利用疊積法計(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檢測出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].測繪學(xué)報(bào)，2012，41(2):171-176.

[10] 張克亮，甘衛(wèi)軍，周新.GRACE衛(wèi)星重力場同震變化的經(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í)變重力梯度探測大地震的同震變化[D].武漢：武漢大學(xué)，2011.

[17] 鄒正波，李輝，吳云龍,等.尼泊爾M8.1地震震前衛(wèi)星重力場時(shí)變特征[J].大地測量與地球動力學(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.