利用重力衛星監測尼泊爾Ms8.1地震前后重力變化

崔立魯， 周 甜， 張 誠， 鄒正波， 李 盼， 宋 哲

(1.成都大學建筑與土木工程學院， 成都 610106； 2.武漢大學測繪學院， 武漢 430079；3.中國地震局地震研究所， 中國地震局地震大地測量重點實驗室， 武漢 430071)

2015年4月25日14時11分尼泊爾發生了8.1級地震，震源深度為15 km，震中位于尼泊爾中東部郎塘地區(Lamjung)。這次地震造成人員傷亡共計31 089人，經濟損失可能超過50億美元。強震的發生會引起大范圍內的地表物質遷移，從而發生劇烈的重力變化，因此通過監測地震地區及其周邊重力變化探測地震孕震機制的可能性是存在的[1]。Han等[2]首次利用重力恢復與氣候實驗(gravity recovery and climate explorer, GRACE)時變重力場數據提取了2004年蘇門答臘Mw8.1地震的同震重力變化，其沿斷裂帶呈現出“正-負”分布。隨后大量的研究成果表明GRACE衛星能夠探測到8.0級及以上的地震所引起的重力變化[3-5]。

Matsuo等[6]利用GRACE衛星數據成功探測到2011日本大地震的同震重力變化；Tanaka等[7]則基于GRACE時變重力場數據提取了震后重力變化；Dai等[8]、Zhou等[9]通過GRACE衛星數據反演了2011年日本大地震震源參數；鄒正波等[10]利用GRACE數據提取了尼泊爾地震震前重力變化，與位錯模型得到的同震重力變化保持一致；瞿偉等[11]通過GRACE衛星重力資料計算得到尼泊爾地震震中及周邊重力和地表密度變化，為地震孕震機制的研究提供了有益的參考；楊九元等[12]利用地殼分層模型對尼泊爾地震同震及震后形變場及重力變化進行了模擬計算，結果表明與實測數據一致。

本文利用2012年1月—2017年6月的GRACE RL06 Level-2時變重力場數據對2015年尼泊爾地震前后及同震重力變化信號進行提取，詳細分析地震前后震中及周邊重力變化的時空變化特征，并對震中和若干特征點的重力變化時間序列進行多項式擬合分析。根據數據分析結果探討大地震地震孕震機制及其構造形變。

1 理論方法

1.1 GRACE數據處理

利用地球重力位球諧系數計算地表重力變化的表達式為[13-14]

(1)

為了削弱衛星軌道誤差、雙星K波段測距誤差、加速度計誤差等影響以及重力場模型高階次項相關性誤差，本文采用300 km 高斯濾波加P3M8多項式濾波的組合濾波算法對GRACE時變重力場球諧系數進行平滑處理[15-18]。同時由于GRACE衛星軌道的原因，得到的C20項精度較低，一般采用由衛星激光測距(satellite laser ranging，SLR)測量得到的C20項進行替代[19]，并利用Swenson等[20]的成果對一階項進行了地球質心變化改正，在扣除冰后回彈的影響[21]。

在提取地震前后重力變化信號時，本文利用2004—2010年的平均重力場作為背景場，然后將地震前后每年的平均重力場減去背景場。其中，因尼泊爾地震發生在2015年4月，在計算平均重力場時需要扣除相應年份4月的數據。

1.2 扣除陸地水重力時變信號

陸地水儲量變化是構成GRACE時變重力信號的主要部分之一，因此在提取由地震所引起的時變重力信號必須從GRACE信號中扣除由水文部分造成的影響。首先利用全球陸地資料同化系統(global land data assimilation system, GLDAS)水文模型的全球格網數據計算得到相應的球諧系數，其表達式為[22]

Δh(θ,λ)sinθdθdλ

(2)

選取與GRACE衛星數據相同月份的GLDAS水文模型數據，對該數據按照式(2)進行球諧展開；然后從GRACE時變重力場球諧系數中扣除GLDAS水文模型對應的球諧系數；最后將扣除水文影響后的GRACE時變重力場球諧系數按照1.1節中的數據處理流程，計算得到重力變化。

2 數據說明

2.1 GRACE時變重力場數據

本文采用由美國得克薩斯大學空間研究中心(Center for Space Research of Texas University in Austin, CSR)發布的RL06版本Level-2的月重力場模型，截斷階數為60。該模型可從該機構網站上免費下載，同時GRACE月重力場模型保持著實時更新。該重力場模型已扣除了由潮汐及非潮汐引起的大氣和海洋質量變化影響，在此基礎上減去地球長期靜態重力場的影響，可以得到地球表面質量變化的結果，其主要包括陸地水儲量變化、冰川質量變化和固體地球質量變化等[23]。

2.2 GLDAS水文模型

GLDAS水文模型是由美國航空航天局(National Aeronautics and Space Administration, NASA)戈達德空間飛行中心(Goddard Space Flight Center, GSFC)與美國海洋和大氣局(National Oceanic and Atmospheric Administration, NOAA)國家環境預報中心(National Centers for Environmental Prediction, NCEP)聯合發布的。

3 實驗結果

3.1 濾波處理結果

首先利用GRACE RL06數據按照式(1)的計算原理得到全球重力變化分布，如圖1所示。從圖1(a)中可以看出，未經過濾波處理前存在著明顯的南北條帶誤差影響，無法提取正常信號，而經過濾波以后可以清楚地看到全球重力變化分布情況。

圖1 2015年3月全球重力變化分布圖Fig.1 The map of gravity changes in global gravity on March 2015

3.2 水文信號的影響

圖2為利用GRACE RL06數據計算得到的2015年3月尼泊爾地震震中及其周邊重力變化。圖中紅色五角星表示震中位置。從圖2可知，在震中附近出現了明顯的正負重力變化對稱分布，在震中的西南和東南方向出現了明顯的正重力變化，而西北和正東方向出現了較大的負重力變化。當扣除水文信號影響以后，西南和東南方向的正重力變化在數值上明顯地較小了，結合地理知識發現，這兩處正重力變化出現的地方剛好是印度的恒河流域和孟加拉國的恒河入海口，這里存在著大量的地表徑流，因此這兩處發生的正重力信號變化正是由水文信號所引起的，而圖2(b)的結果正好驗證了這一結論。

圖2 2015年3月尼泊爾地震及周邊重力變化Fig.2 The gravity changes of Nepal earthquake and it’s surrounding on March 2015

3.3 地震前后重力變化

3.3.1 年際變化

圖3為2011—2017年尼泊爾震中及周邊地區地震前后重力年際變化情況。由圖3(a)～圖3(e)可知，在尼泊爾地震震中及周邊地區出現了明顯地正負重力信號的時空差異分布，其中正負重力信號以震中為中心呈現出明顯的四象限分布特征。其中負重力信號出現在喜馬拉雅地區和緬甸板塊北部地區，兩者的連線基本上與印度板塊和青藏高原的分界線相重合，其值呈逐年遞減的趨勢。這是由于印度板塊不斷地向歐亞板塊下方俯沖，引起印度板塊物質向青藏高原底部聚集造成的；而在青藏高原內部、印度板塊北部和華南板塊東側為正重力信號，前兩者的連線基本垂直與印度和青藏高原的分界線，其重力值呈逐年遞增趨勢。其中青藏高原內部的正重力信號在2013年達到了最大值，這與青藏高原底部出現物質聚集相關，而華南板塊東側的正重力信號增加，可能是由于青藏高原物質東移時，受到華南板塊阻擋所引起的，這些與尹鵬的研究結果相符[24]。

以2015年為分界線，圖3(a)～圖3(d)為地震發生前，而圖3(e)～圖3(g)為地震發生后。無論地震發生前后，正負重力信號的空間分布基本上是沒有變化的，都是以震中為中心呈現出四象限分布。在地震以后，正負重力信號的變化呈現出對稱趨勢，即負重力信號呈現出先增加后減少的趨勢，而正重力信號則表現為先減少后增加的趨勢。

從圖3中正負重力信號和震中位置的相對關系可知，尼泊爾地震發生在重力變化的零值線附近，這與尹鵬[24]和鄒正波等[10]的研究成果一致。這表明利用GRACE衛星監測地震發生前重力變化的時空分布有助于為地震發生位置的確定提供一定的幫助。尼泊爾地震的孕震背景為印度板塊沿北東方向以45 mm/a的速率不斷地擠壓歐亞板塊，在地震發生的斷層南部呈現出物質壓縮、質量增加、重力增加的趨勢，而北部地區則表現為物質受到拉伸，重力減少[10]。

3.3.2 季節變化

本文通過分別將2004—2010年第一至四季度的平均重力場作為背景場，用2011—2017年第一至四季度的平均場減去相應的背景場，從而得到尼泊爾地震震中及周邊地區震前5年和震后2年重力變化，結果分別如圖4和圖5所示。相對于年際變化，季度變化以更高的時間分辨率觀測和探索地震發生前后震中及周邊地區重力變化的時空分布特征，進一步了解地震的孕震機制。

圖3 尼泊爾震中及周邊地區地震前后重力年際變化Fig.3 The gravity interannual changes before and after earthquake of Nepal and it’s surrounding area

從圖4中可以看到很明顯的重力變化時空特征：①2011年第一到四季度該地區主要呈現以負重力變化為主的信號，在第四季度時印度板塊和緬甸板塊北部呈現出明顯的正重力變化，其中印度板塊從負重力變化轉變為正重力變化。②從2012—2015年第一季度可知，沿著青藏板塊和印度板塊分界線，以震中為對稱中心出現了明顯的負重力變化，其重力值逐漸減少，在地震前一個季度兩塊負重力變化區域面積和形狀較為相似。同時在印度板塊北部和緬甸板塊東部出現了正重力信號，并逐漸增強。到了地震發生前一個季度，信號開始減弱。③結合2012—2015年第二季度的重力變化情況可以看出，印度板塊的重力變化較為顯著，其中位于印度板塊和緬甸板塊之間的負重力變化在不斷地擴大，其值也在不斷減少，其中整體呈現向東移動的趨勢。④從2011—2014年第三季度可知，位于喜馬拉雅地區和緬甸北部的負重力信號先增強再減小，覆蓋范圍也是先擴大后縮小。⑤根據2011—2014年第四季度的情況可知，印度北部和緬甸中部出現了明顯的正重力信號，并表現為逐漸增大的趨勢，到了2013年第四季度達到峰值，而2014年第四季度則開始變小。⑥在2011—2015年間，青藏高原內部一直表現為正重力信號的特征，且大小和位置都較為穩定。

圖4 尼泊爾震中及周邊地區震前重力季節變化Fig.4 The gravity seasonal changes before earthquake of Nepal and it’s surrounding area

圖5為2015—2017年尼泊爾地震震后重力季度變化。整體上，喜馬拉雅地區和緬甸板塊北部的負重力信號在震后一直存在，在震后第一個季度負重力信號最弱，之后一直呈增大趨勢，且覆蓋范圍也不斷擴大。而之前出現的正重力信號一度減弱，又在2016年第三季度開始增強。這說明在震后的一段時間內震中及周邊板塊內部物質在不斷地調整。

3.4 特征點重力變化

為了進一步研究震中及周邊區域震前震后重力變化情況，本文選取了5個重力變化較為顯著的點，加上震中位置(28.2°N，84.7°E)，共計6個點，即A(22.0°N，78.0°E)、B(31.0°N，77.0°E)、C(21°N，100.0°E)、D(29.0°N，94.0°E)和E(35.0°N，87.0°E)，其分布如圖4(l)所示。A位于印度板塊，B點位于印度板塊和歐亞板塊之間，C點位于緬甸板塊，D、E位于歐亞板塊的青藏地區。同樣的，以2004—2010年1—12月的平均重力場作為背景場，將2011年1月—2017年6月每個月的重力場減去背景場，獲得以上6個點2011—2017的年重力變化的時間序列，結果如圖6(a)所示。

圖5 尼泊爾震中及周邊地區震后重力季節變化Fig.5 The gravity seasonal changes after earthquake of Nepal and it’s surrounding area

圖6 點A～E及震中點重力變化時間序列Fig.6 Time series of gravity changes at A～E and epicenter

由圖6可知，點A～E及震中點的重力變化時間序列具備顯著的周期性，這是因為GRACE衛星所探測到重力變化除了由地震引起的之外，還有陸地水循環和地球內部物質移動等所引起的周期性變化。為了從該時間序列中提取由地震引起的異常重力變化，本文扣除了上述周期性變化，得到新的重力變化時間序列，如圖6(b)所示。

由圖6可知，除了點E之外，其余點在震后都發生重力值的躍遷，這與所在區域的重力累計變化情況保持了一致。其中A點從正值變為了負值，C點從3.70×10-8m·s-2急劇下降到1.13×10-8m·s-2，震中點從1.36×10-8m·s-2降到了-2.76×10-8m·s-2。

①位于印度板塊A點從2011年的-0.73×10-8m·s-2開始緩慢下降到2012年的-1.54×10-8m·s-2，然后開始緩慢增加，到2013年9月達到峰值，其值為4.63×10-8m·s-2，接著開始下降，直至0.4×10-8m·s-2，呈現出減少、增加、減小的變化特征。其震后一年左右的時間處于負重力信號狀態，這與震前長期處于正信號狀態形成對比。②位于印度板塊和青藏板塊之間的B點，整體上呈現出減小趨勢，但是在震前一年的時間內出現了先減小后增加的特征，這與點A的變化趨勢保持一致。震后繼續呈較小趨勢。③位于緬甸板塊的C點與A點重力變化趨勢基本相同，其最高峰值出現在2014年1月，震后重力變化較為平穩，但是一年內波動較大。④D點位于青藏板塊的南部邊緣，整體呈現出減小趨勢，這與B點基本一致，在震后一個月內出現了異常的負重力信號，且震后第二年出現較大的波動。⑤而位于青藏板塊內部的E點呈現出先增大再減小的趨勢，在震后表現為一直減小的趨勢，與D點相同的是出現了異常負重力變化信號。⑥而震中點在震前重力變化的變化趨勢與點A～C保持一致，即呈現出先減少，再增加，其后減少的變化特征，震后出現了異常負重力變化信號，直到兩年后才恢復正重力變化信號，這說明震后兩年震區內部的物質調整完成，恢復震前狀態。

4 結論

本文利用GRACE時變重力場數據反演尼泊爾地震震中及周邊地區震前震后連續6年多的重力變化，包括年際變化和季度變化，并對變化較為顯著的特征點時間序列進行了數值分析。實驗結果表明：長期重力變化具有明顯的時空分布特征，正負重力變化信號在印度板塊和青藏板塊的弧形分界線附近呈四象限對稱分布，而震中位置位于重力變化的零值線附近。各特征點在震前和震后表現出明顯的不同趨勢，其中在震后均出現了不同大小的躍遷和較大振幅的波動，這些均說明震區內部物質在進行均衡調整。著充分說明利用GRACE計算得到的重力變化可以有效地反映青藏高原地區大尺度的構造動力學特征，對于地震孕震機制的解釋說明也具有一定的科學價值。