重力场恢复与气候试验卫星(GRACE)发射于2002年，设计运行轨道高度约500 km，可提供覆盖全球大部分地区的时变重力场数据.GRACE卫星重力场的时变特性、良好的海陆覆盖等特点显示出在大范围质量变化监测中的优势.利用GRACE卫星的时变重力场数据，可对水文等多种地表质量迁移现象进行定量化研究，如青藏高原的冰川质量变化(Yi and Sun, 2014；刘杰等, 2015)、亚马逊流域的降水(Tapley et al., 2004)、极地冰盖的质量变化(Velicogna and Wahr, 2005；鄂栋臣等, 2009；鞠晓蕾等, 2013)等.同时，GRACE卫星也可应用于对固体地球内部的质量迁移的研究，特别是地震引起的重力变化(Han et al., 2006；Chen et al., 2007；王武星等, 2010；Matsuo and Heki, 2011；邢乐林等, 2012；邹正波等, 2012).研究表明GRACE卫星对M7.5以上的地震具有理论上的探测能力(Sun and Okubo, 2004).从实际应用看，GRACE可探测MW8.0以上的地震(Tanaka et al., 2015).许多研究对GRACE的地震学应用进行了探索，如GRACE可被应用于同震重力变化(Matsuo and Heki, 2011)、震后重力变化(Tanaka and Heki, 2014)、震源参数反演(Dai et al., 2014；Zhou et al., 2018)等研究中，并得出了许多的重要成果.其中，在对2004年苏门答腊—安达曼地震(以下简称苏门答腊地震)的研究中，GRACE卫星重力观测到超过10μGal的同震重力变化(Han et al., 2006).对同震、震后变化的监测将有利于对地震活动机制的研究，同时也可作为基于大地测量的俯冲带地震循环研究(Wang et al., 2012)的重要补充部分.

2010年MW8.8智利Maule地震(以下简称智利地震)是GRACE运行期间震级最高的三次地震之一(另外两次分别是2004年MW9.0苏门答腊地震和2011年MW9.0日本东北地震).地震引发了震中附近区域的大规模地表形变.震中附近的GPS台站记录到了地震引起的广泛的同震位移以及震后长期位移(Heki, 2011).在Auraco半岛观测到1.8~2.5 m的最大隆升(Farías et al., 2010)，断层的运动将海岸线向海洋方向推动了0.5km.而在中央山谷(Central Valley)则出现了大范围沉降，GPS观测到沉降值30~70 cm(Vigny et al., 2011).震后12天的震后余滑大约产生了相当于同震4%的地震矩(Vigny et al., 2011).地震甚至引发了南安第斯火山地区超过15 cm的沉降(Pritchard et al., 2013).地震亦给智利中部地区造成了巨大破坏.大规模的形变伴随着急剧的重力变化，此前对智利地震的重力变化已有相关研究.一部分研究者对同震重力变化进行了观测和建模解释，例如Heki和Matsuo(2010)对智利地震的同震变化进行了计算，并使用位错模型进行了解释.周新等(2011)也进行了类似的工作，观测到最大~5 μGal的同震重力变化.Han等(2010)对智利地震的同震重力变化进行了分析，对于智利地震同震重力变化分布的不对称性利用膨胀模型进行了解释.另一部分研究者则对震后重力变化的规律进行了探索，例如Tanaka和Heki(2014)对智利地震的单点震后重力变化进行了研究，并与其他俯冲型大地震进行了比较.

本文基于最新的GRACE RL05 Level-2时变重力数据对2010年智利地震的同震及震后长期重力变化进行了研究.从时间上，对震后每年的重力变化结果进行计算；从空间上，对智利地震的单点和大范围的长期震后重力变化进行多项式拟合分析，并基于位错理论对同震及震后重力变化进行了解释.

本文使用美国德克萨斯大学空间研究中心(CSR)发布的从2002年到2016年共153个月的GRACE RL05 Level-2时变重力场数据.数据由60阶球谐系数组成，对应空间分辨率~330 km.RL05数据相较上一代RL04数据，对于噪声的处理更加优化，同时信号也保持了与RL04相当的强度(Bettadpur, 2012).

为去除噪声和提取地震信号，本文对GRACE RL05 Level-2重力场数据进行了多项处理.首先，使用卫星激光测距(SLR)得到的C20项(Cheng and Tapley, 2004)替代了受制于GRACE极地轨道而精度较低的原GRACE RL05 Level-2数据中的C20项.其次，使用330 km高斯滤波(Colombo, 1981)对重力场结果进行了滤波，去除高频噪声的影响.在此基础上，使用P3M6多项式拟合的方法(Swenson and Wahr, 2006)进行了去相关滤波，去除GRACE信号中存在的南北条带噪声.最后，对非地震因素(如水文等)等周期信号引起的影响，在同震或震后重力变化的计算中使用年平均做差或者多项式拟合的办法予以去除.

同震和震后重力变化的观测量为地表处的重力异常的变化，其计算公式(Hofmann-Wellenhof and Moritz, 2006)如下所示：

其中，g为地表处重力加速度，G为地球的引力常数，r为地球平均半径(本文取r = 6371 km)，P()为归一化缔合勒让德函数，θ和λ分别为余纬和经度，Cnm和Snm为GRACE平均月重力场完全规格化系数.本文使用的截断阶数为Nmax=60.

本文利用上述数据和方法计算智利地震的同震和震后重力变化.同震变化的研究可用于监测大地震在短期内对震中周围地区造成的地表变形、优化震源模型、了解地球内部在地震发生后的密度变化.震后变化的研究则将有助于探索地震发生区域的地下结构的黏弹性性质.

同震变化的一般计算方法为前后两年的平均重力场结果做差，平均做差计算方法的公式见式(2)—(5).这一方法能够有效地去除周年和周年以内的周期性噪声信号.本文除一般的前后两年平均做差方法外，对智利地震前后一年和三年的重力场也进行了平均做差.式中gi为某月的GRACE RL05 Level-2月重力场数据，n为实际使用的月份数.在计算中为减小误差，降低非地震的周期信号(如周年水文信号)等因素的影响，使用的数据不包含地震发生的当月的重力场数据、缺失的数据以及缺失的数据对应的他年同月数据.MW8.8智利地震发生于GMT 2010年2月27日，按照上述平均做差的方法对2008—2011年四年的数据进行处理，得到最终的同震变化结果.本文研究中所使用的GRACE数据的月份分布情况如图 1所示.

本文对智利地震的年总震后重力变化(2011—2016年)以及单点的重力变化序列进行了计算.

(1) 总震后重力变化

年总震后重力变化的计算使用震后某年的年平均重力场与2008—2009年(震前)的平均重力场做差得到(如式(6)).本文将2008—2009年的平均重力作为震前的平均重力场：

(2) 单点重力变化序列

对于单点的重力变化，计算出153个月的重力值后，为凸显震后变化趋势，对其进行多项式拟合，可提取出包含的线性项、周期项以及震后衰减项.其中线性项指重力随时间的线性变化，周期项指重力信号中包含的周年、半年以及其他周期信号，震后衰减项则是用指数衰减函数模拟的震后重力变化.

本文的拟合公式基于Tanaka和Heki(2014)提出的双参数震后拟合公式，并添加一项162天S2周期项作为补充，具体如下：

其中，A、B为震前的线性拟合系数，C、D、E分别为水文信号周年、半年、S2潮汐项强度，H(Δt)是Heaviside阶梯函数(Δt为发震后的历时)，F为同震变化强度，G1、G2表示震后变化强度，指数部分是震后的衰减项，τ为常数.式中t为时间，Δt为距震后的时间，ω、ωS2为对应的周年以及S2潮汐水文信号的变化速度，θ1、θ2、θ3为对应的信号初始相位.在对智利地震的同震和震后重力变化研究中，本文同时对有两个衰减项(双参数)和一个衰减项(单参数)的两种假设进行了探究.

智利地震震中位于36.290°S 73.239°W (根据：https://earthquake.usgs.gov/earthquakes/eventpage/official20100227063411530_30)，震中附近的同震重力变化现象较为明显.同震重力变化(如图 2所示)关于震中对称，在俯冲带大陆侧形成负重力变化，在俯冲带海洋侧形成正重力变化.负重力变化呈椭圆状，极值位于中心，约-6 μGal.正重力变化的极值同样位于其中心，约为2 μGal.震中附近的重力变化约为0.正负重力变化的面积差异较大，负重力变化面积显著大于正重力变化.同时，在正重力变化以西的太平洋上存在一个较弱的负重力变化信号，大小在2 μGal以内，呈条带状.在地震同震变化以东和以北的地区可以观测到在亚马逊河以及拉普拉塔河区域可能由于河流和降水造成的水文信号.本文基于GLDAS水文数据模型对该地区的水文重力信号进行了验证.计算中使用了该区域2009—2012年的年平均土壤湿度数据进行了做差(结果如图 4所示)，并模拟出水文变化引起的重力变化信号.模拟方法参考John Wahr的相关论文(Wahr et al., 1998).通过水文数据的模拟，发现在图幅东北方向的拉普拉塔河流域存在一个1~2 μGal的水文重力信号.除去此水文重力信号后，原图的东北方向的正重力变化信号消失，说明此信号主要由拉普拉塔河流域的土壤水质量变化导致.

利用从智利地震发生到目前已累积较多的GRACE时变重力场数据，本文对智利地震震前震后一年、三年的数据进行了平均做差(如图 3a、3c所示).不同年数累积计算的同震重力变化的结果数值相近但有所差异.对比一年、两年和三年平均做差结果，正重力变化的范围波动较大，前后三年平均的结果大于一年平均.这表明同震正重力变化的计算可能受到了震后信号、水文信号等其他信号的影响.另一个显著特征是在三年平均结果中，在一年、两年的结果的图幅东北方向均出现拉普拉塔河水文信号消失.三年平均做差的结果中数值波动较少，与一年、两年结果相比，同震信号更为清晰.

使用前文所述的方法，本文计算了2010年智利地震震后2011—2016年每年的总震后变化(含有同震变化)结果(如图 5所示).总震后重力变化具有两个明显特征：1)震后重力变化随时间整体呈减小的趋势，在2013—2016年尤为明显，重力正负变化区面积迅速缩小；2)震后重力在震后的2—3年内有增大的趋势，主要表现为2012—2013年与2011年相比较重力变化超过-4 μGal的区域面积有所增大.这一重力变化增强的现象可能与震后余滑(afterslip)、黏弹性松弛等多种原因有关(Tanaka and Heki, 2014).2011年是地震后的第一年，此时的总震后重力变化与地震的同震重力变化几乎相同，正负重力变化的位置和强度都相似.主要的差异在于正重力变化稍强，达到2 μGal以上的面积更大.在正重力变化以西的海区仍可以看到一个较弱的负重力变化条带，约-0.5 μGal.负重力变化以东的区域与同震结果相比有微小差异，负重力变化的面积更大.2012年，正重力变化的面积和位置几乎无变化.海洋上的弱重力变化减弱，面积减小.负重力变化的面积增大，与周围的梯度也增大了.重力值在-4 μGal以下的面积增大.同时，注意到负重力变化以东的区域存在其他因素引起的正重力变化，达到1~2 μGal.2013年，正负重力变化的面积均增大，整体图形更加平缓.在负重力变化东北处出现了一个约-2 μGal的新的负重力变化.2014年，正负重力变化均明显减弱.正重力变化的极值减小为~2 μGal，负重力变化的极值仍然在~-4 μGal左右.这一年相对于震前，正负重力变化周围信号无明显变化.2015年，正重力变化消失.负重力变化面积进一步变小，极值也在变小.拉普拉塔河流域水文变化相关的正重力变化进一步凸显.2016年相较于2015年而言，正重力变化重新出现，且达到4 μGal，远大于原来2 μGal左右的正重力变化.这个新的正重力变化与负重力变化的大小相当，图形相似，关于震中对称.震中对称的负重力变化也重新出现，与正重力变化相似，达到-2 μGal，面积略小一些.为研究2016年重新出现的正负重力变化，将研究区域放大到整个南美洲南部时发现，在结果中主要包括的信号有：亚马逊流域和拉普拉塔流域的水文信号、巴塔哥尼亚地区的一个较强的冰川信号和其他噪声信号.此时出现的正负重力变化与周围的重力变化达到~4 μGal，从整体看是比较明显的信号.

MW8.8智利地震的同震和震后变化中存在正负对称变化的特征，本文选择具有代表性的正负极值点提取重力变化时间序列，分别为正变化点B(77°W 37°S)和负变化点A(69°W 35°S)(如图 2).将A、B点从2002年首个GRACE重力场月结果到2016年第153个月重力场结果的重力值提出后，使用多项式拟合(Chen et al., 2007)的方法去除了年周期、半年周期和162天S2固体潮项的影响.最终结果如图 6所示，纵轴是重力加速度，横轴是年份，灰色直线为地震发生的时间.拟合基于震后变化分为短期和长期的假设，故设两个衰减时间τ1和τ2.根据拟合中的最小残差值得到拟合最优的τ组合(其结果见3.2节)(图 7).

点A和点B的重力变化时间序列均反映了明显的同震和震后重力变化.正变化点B处的同震变化相对较小，在~2μGal左右.震后的点B处曲线呈缓慢上升的趋势.负变化点A处的同震变化比较剧烈，曲线在2010年地震时急剧下降，达~10μGal.震后点A处亦出现与点B处相似的增大趋势，A点的增大速度更快.A点在震前还具有线性减小的趋势，推测可能受到全球尺度气候变化/质量迁移等的影响.两点处的同震变化均与做差得到的同震变化值相符合.A、B点震后的变化时间都比较长，增大的趋势一直延续到2016年，这与前文的震后总变化结果相符合.

本文的观测结果基于GRACE RL05 Level-2数据，整体上与其他研究机构利用GRACE RL04 Level-2数据(以下简称RL04)得到的智利地震同震变化结果相近.2010年智利地震发生后，许多研究小组发表了基于GRACE Level-2数据的同震观测结果.Heki和Matsuo(2010)使用CSR RL04数据，基于扇形滤波(Zhang et al., 2009)后计算出负重力变化的极值在-4 μGal左右，而未发现明显的正重力变化.周新等(2011)基于RL04数据和高斯300 km滤波，得到的智利地震同震变化的极值达7 μGal.在正重力变化区选择的特征点变化为2 μGal，在负重力变化区为-5 μGal.在类似研究中，可以看到和本文类似的在负重力变化西南方向较弱的正重力变化.本文计算结果与以上结果在极值方面相当(2 μGal，-4 μGal)，而在正负重力变化的具体形态上还有所差异(表 1).

智利地震的GRACE同震重力变化观测结果也反映出大型俯冲型地震的共同特点.本文利用GRACE RL05 Level-2数据，同样使用年平均做差法对2004年苏门答腊地震、2011年日本东北地震进行了同震变化的计算(如图 8).三次地震的计算结果表现出相似性：(1)地震的同震变化的正重力变化和负重力变化均关于震中呈对称分布；(2)负重力变化的面积均大于正重力变化的面积；(3)负重力变化呈椭圆状，正重力变化呈条带状；(4)负重力变化均位于俯冲带大陆侧，而正重力变化位于俯冲带大洋侧.

智利地震的GRACE同震重力变化有其自身的特点.由于震级的关系(MW8.8)，智利地震的同震变化幅度相对于其他两次地震(苏门答腊、日本)较小.苏门答腊的同震正负极值为+5 μGal、-10 μGal左右，而智利地震的极值仅在+2 μGal、-5 μGal左右，约为前者的一半.智利地震观测到的正重力变化较小，而负重力变化较大，两者对称性较弱，未出现其他地震所具有的明显对称分布(见表 2)，这一特点与海水的补偿作用有关(Heki and Matsuo, 2010).

本文基于单参数和双参数两种模式对智利地震的震后重力变化进行了拟合.在单参数拟合下，计算得到观测区域内各个点的最小残差参数如图 9a所示.除图幅东部和南部区域外，其余大部分地区的震后变化特征时间约为1.1年左右.同时求得的单参数震后重力变化强度G1分布如图 9b所示，除东北部水文信号之外，在震中附近区域出现了~1 μGal的震后重力变化，与其他区域相比具有明显差异.这一方法亦得到了最佳拟合的同震的重力变化结果(如图 9c)，与平均做差得到的同震重力变化相似.

在双参数拟合下，在整个图幅区域得到的长期重力变化特征时间均为τ1=1.1年，而短期重力变化则呈条带状分布(如图 10)，在震中附近区域主要为τ2=0.6年.这一方法也得到了拟合的同震重力变化(如图 11a)，从同震变化看，正重力变化几乎消失，而在长期和短期重力变化上未表现出很大的变化值(图 11b、11c).与单参数的拟合结果对比，主要结果反映出智利地震未出现与其他地震(苏门答腊、日本东北)相似的明显的两段震后重力变化过程.智利地震的震后变化可概括为特征时间约1.1年的长期衰减过程，主要发生于震中附近，大小约为1 μGal.

为研究2010年MW8.8智利地震的同震和震后重力变化内在机制，本文基于位错理论对地震的重力变化过程进行数值模拟.本文使用的位错理论基于黏弹性地球模型，此前已广泛应用于地震同震和震后重力变化的模拟中(Pollitz et al., 2006；De Linage et al., 2009；Li et al., 2016).计算使用的程序为德国地学中心(GFZ)汪荣江研究员开发的黏弹性位错计算程序PSGRN/PSCMP(Wang et al., 2006).计算使用的地球模型参数如表 3所示.模拟得到地表重力变化经过自由空气改正和海洋区域的布格重力改正后，按球谐函数展开至60阶(与本文使用的GRACE数据相同)后，进行330 km高斯滤波得到最终结果.本文模拟中使用的断层模型采用Tong等(2010)基于地震波数据的计算结果.

黏弹性位错理论计算得到的2010年智利地震总震后重力变化结果如图 12所示.与重力卫星GRACE观测到的结果相比较，黏弹性模型模拟的重力变化在震中附近符合的较好.震后第一年(2011年)，震中附近区域的正重力变化极值约2 μGal，负重力变化极值约3 μGal，与GRACE的观测结果相近.同时，正重力变化与负重力变化等值线形态上有所差异，负重力变化强于正重力变化，与GRACE卫星重力的结果基本相符.从模拟结果中还可以看出：2010年MW8.8智利地震的震后重力变化整体呈衰减趋势；震后第2年(2012年)年正重力变化有所增强，与GRACE观测到震后第2年正重力变化增强基本相符.

据此前的相关研究，造成震后重力变化的机制主要有：震后余滑、地球的黏弹性以及孔隙流体回弹，其中前两者是主要的震后重力变化机制.从位错理论计算与本文的重力卫星GRACE观测结果的对比来看，黏弹性模型较好地拟合了震后的重力变化过程.此前有研究指出，余滑造成的重力变化与震后的短期重力变化相联系(Tanaka and Heki, 2014)，而本次智利地震中未出现明显的短期震后重力变化，故可认为地球的黏弹性性质是造成智利地震震后重力变化的主要机制.

(1) 正负重力变化的变化不同步

从整体上看，智利地震的总震后变化呈现出重力逐渐减小的恢复趋势.但是注意到正重力变化和负重力变化的恢复速度并不完全相同.2013年以前正负重力变化均无明显的变化.2014年，负重力变化大大减小，原来的区域在-4 μGal以上的变为-2 μGal，而此时正重力变化并未随之减小.2015年则是正重力变化消失，负重力变化与2014年相当.

(2) 正负重力变化的对称性

从震后总重力变化图(图 5)来看，正负重力变化的对称性随时间逐渐变好.正重力变化的面积(>2 μGal的区域)逐年增大，而负重力变化对应区域(< -2 μGal)则未明显增大.与同震变化比较可以推知，正重力变化区域在经历了同震变化以后，震后逐年仍在以一定速度增长，且短期内增长速率快于负重力变化区域.这一推测与单点序列的重力变化(图 6)相吻合.

(3) 2016年的重力变化信号

经过2011—2015年的变化，重力变化的强度已相对减弱，而2016年再次出现较强的重力变化信号，并且该信号呈现出关于震中的良好对称性.信号的正值约4 μGal，负值约-2 μGal.由于GRACE运行到2016年精度已有所下降，部分月份缺失，该信号的真实来源有待研究.

2010年MW8.8智利地震引发了震中区域的同震和震后长期重力变化.本文使用2002—2016年共153个月的GRACE RL05 Level-2月平均时变重力场数据，进行330 km高斯滤波、P3M6去相关滤波后，采用平均做差、最小二乘拟合的方法，获得了2010年智利地震的同震和震后长期重力变化结果.在重力变化的计算中，为减小水文信号的干扰而采用3年平均的方法可以获得好的效果.从重力变化结果中可得出以下结论：

(1) 同震重力变化特征

与同为俯冲型大地震的2004年苏门答腊地震、2011年日本东北地震相比，2010年智利地震的同震变化呈现出一正一负类似的对称分布.在俯冲带陆地一侧，重力变化呈负变化的趋势，在海洋一侧则呈现出正变化.2010年智利地震与同类型俯冲带地震的显著差异是其同震变化的绝对值较小.苏门答腊地震的同震变化极值在约-11 μGal，而智利地震的同震变化的极值仅约-5 μGal.智利地震的重力变化分布对称性不显著，正重力变化的分布区域较小.

(2) 震后长期重力变化特征

对2011到2016年的GRACE年重力变化分析发现，智利地震的震后重力变化呈现先增长后衰减的趋势.此前有研究指出俯冲型大地震存在震后重力变化的短期和长期两个阶段，而本文中计算的智利地震的震后长期重力变化未出现明显的短期、长期阶段的变化特征.对震后重力变化的最小二乘拟合结果表明，智利地震震后仅存在单个时长约1.1年的重力变化过程.这一过程呈指数衰减，重力变化值约为1 μGal.

(3) 震后重力变化的解释

在基于GRACE卫星对2010年MW8.8智利地震的研究中未发现明显的短期重力变化过程，同时基于黏弹性地球模型的位错理论模拟与GRACE观测的震后重力变化较为相符，故认为2010年MW8.8智利地震的震后重力变化主要由地球的黏弹性性质导致.