利用重力卫星监测尼泊尔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计算得到的重力变化可以有效地反映青藏高原地区大尺度的构造动力学特征，对于地震孕震机制的解释说明也具有一定的科学价值。