2015, Vol. 58
2. 中国科学院计算地球动力学重点实验室, 北京 100049;
3. 大地测量与地球动力学国家重点实验室, 武汉 430077;
4. 东京大学地震研究所, 东京 113-0032
2. Key Laboratory of Computational Geodynamics, Graduate University of Chinese Academy of Sciences, Beijing 100049, China;
3. State Key Laboratory of Geodesy and Earth's Dynamics, Wuhan 430077, China;
4. Earthquake Research Institute, University of Tokyo, Tokyo 113-0032, Japan
We adopt the RL05 data provided by UT/CSR (Center for Space Research, University of Texas), and the data spans from January 2003 to December 2014. We replaced the Earth's oblateness values (C20) with those from Satellite Laser Ranging because of their poor accuracy. After the de-correlation filter using polynomials of degree 3 for coefficients with orders 6 or higher to alleviate longitudinal stripes and the Gaussian smoothing with averaging radius of 350 km to reduce short wavelength noises, we obtain the post-seismic gravity changes due to the 2004 Sumatra earthquake (Mw9.3) based on the difference method and time series of two points respectively located on the hanging side and the heading side. We adopt the visco-elastic spherical dislocation theory to validate the GRACE-derived result. We first use the integration of the finite fault models to solve the problem of large errors in the near-field computation using the visco-elastic spherical dislocation theory, and comparing the modified results with the results of elastic spherical dislocation theory.
After that, we retrieve post-seismic gravity changes induced by the earthquake. The result shows that the peak of the gravity changes is from -8.3 to 4.6 μGal in the period of 2004—2005 based on the GRACE, which is consistent with the result derived by the elastic spherical dislocation theory in the spatial distribution and the magnitude. Based on the spatial distribution of the gravity changes and the analysis of time series, we find that the rate of the post-seismic gravity change at the heading side is 0.55 μGal/yr and the rate is 0.16 μGal/yr on the other side, which demonstrates that the viscosity of the two sides of the fault is different.
The differences of coseismic vertical displacement and gravity changes based on the modified visco-elastic theory are respectively 0.06% and 0.1%, compared with those derived from the elastic spherical dislocation theory. We also analyze the Green functions of the vertical displacement and gravity changes at the source depth of 32 km in two different periods of 0~0.7 years and 0~7 years, which demonstrate that we have successfully solved the problem of large errors in the near-field computation using the visco-elastic spherical dislocation theory. Based on the corrected method, we simulate the post-seismic gravity changes due to the 2004 Sumatra earthquake (Mw9.3), finding that only when the viscosity of the heading side is 1×1018Pa·s and the viscosity of the hanging side is 8×1018Pa·s, the simulating results are consistent with the GRACE-derived results in the spatial distribution and the magnitude. Based on the theory results, we find that the rate of the post-seismic gravity changes at the heading side is 0.52 μGal/yr and the rate of the post-seismic gravity changes at the hanging side is -0.12 μGal/yr. Comparing the observing results with the simulating results, it demonstrates that the viscosity of the mantle at the Sumatra region is laterally inhomogeneous. The viscosity at the left side (the heading side) of the fault is smaller than the one at the right side (the hanging side) of the fault, and the viscosity is respectively 1×1018Pa·s and 8×1018Pa·s.
1 引言
众所周知，一次大型地震除了会引起地壳的同震永久变形外，还会因断层余滑或/和黏弹性地幔产生显著的震后形变，这些变形均可被大地测量手段 如GPS和InSAR观测到(Deng et al., 1998; Sheu and Shieh， 2004; Jónsson et al., 2003).自Steketee(1958)将弹性位错理论引入地震学以来，经过几十年的发展，一维径向分层地球模型的弹性位错理论已经较为成熟，该理论被频繁用于解释地表的观测数据以及利用大地测量资料研究地震的震源机制.例如：Okada(1985)总结了前人对均匀半无限空间的位错变形研究工作后，得到了一套关于点源和矩形有限断层的同震位移和应变的解析表达式，并被后人广泛应用；陈运泰等(1975，1979)讨论了结合半无限空间位错理论和地表形变资料进行震源反演的一般方法，并实际反演了1966邢台地震和1976唐山地震的震源破裂过程；Rundle(1980)利用贝塞尔函数得到了忽略自重的层状半无限空间的同震变形；Pollitz(1996)给出了无自重层状球形地球的同震变形；Sun和Okubo(1993)和Sun 等(1996)分别用龙格—库塔方法得到了自重、可压缩的层状球形地球的同震重力变化和位移.Wang等(2003)将Rundle和Okada的理论扩展得到了分层半无限空间的位错理论，并采用传递矩阵的方法得到地表的同震变形；郝金来和姚振兴(2012)研究了均匀分层介质模型的同震位移、应变和应力的计算，使用了梯形积分和Filon积分相结合的方法加快同震位移的计算速度.
将弹性位错理论略微扩展，便可研究震后黏弹性变形，该理论可用于解释黏弹性地幔引起的震后地壳变形和地球的质量迁移.具体地，对黏弹性的本构关系经过Laplace变换后得到与弹性相同的形式，再按照弹性位错理论的解法得到Laplace域的变形，然后对其做Laplace反变换便可得到空间域的黏弹性变形.例如：Rundle(1982)研究了层状地球内矩形逆断层的黏弹性变形；Pollitz(1992)解决了在一个黏弹性无自重的地球模型下的位错源产生的位移和应变；Ma和Kusznir(1994)导出了3层弹性介质模型的位移计算公式，并将其应用于计算张裂断层引起的同震和震后位移；Piersanti等(1995)、Sabadini等(1995)研究了自重不可压缩层状地球模型内位错源在地表产生的位移；Pollitz(2003)从理论上研究了无自重的横向不均匀地球下的黏弹性变形问题；Wang(1999)将Sun和Okubo(1993)的弹性位错理论推广到可以计算震后黏弹性位移和重力变化，但地球模型只能分有限层数.Wang等(2006)给出了自重可压缩分层半无限空间黏弹地球下的变形，并开发了一套可以计算 同震和震后位移、重力、应变和应力的软件PSGRN/PSCMP，目前被广泛应用.由于黏弹性简振模的Laplace反变换是一个积分问题，一般可以化为留数定理寻找奇点来解决，可是该方法的缺陷是随着地球黏弹性层数的增加，奇点个数也随之增加，难以找到全部的奇点，因此只能将地球分为有限层.Tanaka等(2006，2007)利用在复平面路径积分的方法解决该难点，但是需要保证数值积分的精度.Spada和Boschi(2006)利用Post-Widder公式来克服数值积分的困难.该方法提供了一个估计Laplace逆变换的简便方法，避免了Bromwich路径积分和留数的计算，但需要高性能计算来实现.为了避免计算的困难，本文利用Tanaka等(2006)的理论模型来计算震后重力变化和垂直位移.
随着大地测量技术的发展，除了GPS和InSAR能够观测到地震引起的地壳变形之外，重力卫星 GRACE(Gravity Recovery and Climate Experiment)能够观测到特大的俯冲型地震引起的~300 km空间分辨率的长波长重力变化，如2004年苏门答腊地震(Mw9.3)(Han et al., 2006; Panet et al., 2007; de Linage et al., 2009; Cambiotti et al., 2011; Broerse et al., 2011)，2010年智利地震(Mw8.8)(周新等，2011)，2011年日本东北大地震(Mw9.0)(Zhou et al., 2012; Wang et al., 2012; Cambiotti and Sabadini， 2013). 这些研究表明地壳垂直变形引起的海水质量重新分布与固体地球的重力变化为同量级.因此，用理论模型模拟重力卫星GRACE的重力变化数据时，需要考虑海水质量重新调整的影响.除了同震重力变化外，GRACE还能够观测到上地幔黏弹性松弛引起的重力场变化(Panet et al., 2010).其中研究得较多的是2004年苏门答腊地震(Mw9.3).Ogawa和Heki(2007)利用350 km的各向同性高斯滤波器得到了同震和震后大地水准面变化，并用流体扩散模型解释震后大地水准面变化，得到松弛时间约为0.6a，但他们并没有考虑2005年 Nias地震(Mw8.6)的影响.Panet等(2007)利用多尺度小 波分析了CNES(Centre National d′Etudes Spatials)的50阶重力场解，发现在不同的空间尺度下，震后松弛时间也不尽相同.Han等(2008)利用Slepian函数对GRACE观测数据进行处理，得到震后松弛时间为150天.De Viron等(2007)利用经验正交函数对CNES的GRACE观测数据进行处理，结果表明由于苏门答腊地震产生的空间重力信号在Nias地震附近的区域太强而无法检测出2005 Nias地震.De Linage等(2009)分别从CNES和CSR(Center for Space Research University of Texas)的GRACE重力场系数解提取了关于苏门答腊地震的同震和震后重力变化信号，并利用自重、弹性、非旋转的球体地球模型模拟了同震观测信号.Cannelli等(2008)假设一个黏度为1×1018 Pa·s 的Maxwell体软流圈模型模拟了震后的低阶重力场，并与激光测距数据进行了对比.Panet等(2010)结合GRACE重力和远场GPS观测资料，对一个Burgers体模型的上地幔双黏度系数进行约束，得到在220~660 km 深度区域的地幔黏度为8×1018 Pa·s.Hoechner等(2011)认为用一个双黏度Burgers体可更好地拟合GPS观测时间序列，且一个单一黏度Maxwell体模型加上震后余滑与GRACE观测无法吻合，因为断层余滑的影响只集中在非常局部的区域.
本文利用GRACE卫星RL05数据提取了2004年苏门答腊地震(Mw9.3)的震后重力变化，分析断层两侧地幔黏滞性差异，并利用Tanaka等(2006)的黏弹性球体位错理论结合GRACE观测数据模拟该地区地幔黏滞性系数，分析该地区地幔横向黏滞性构造，为该地区地球内部构造研究和地震孕育机理研究等提供支持.
2 GRACE卫星数据处理本文采用CSR发布的RL05数据，此数据为完全正则化的最高达60阶的球谐系数，精度比RL04数据更高，因此RL05数据能提取更高精度的重力变化信号(鞠晓蕾等，2013).为了研究2004年苏门答腊地震(Mw9.3)引起的震中周围地区重力场变化，本文采用2003年1月—2011年12月的数据(其中，2003年6月，2011年1月，2011年6月，同时，由于地震发生在2004年12月，故此月数据舍去不用)进行分析.数据处理过程中，C20系数用人卫激光测距得到的C20项替代.
由于GRACE卫星轨道高度在500 km左右，重力信号随轨道高度增加而衰减，所以GRACE卫星获得的重力变化信号仅仅是中长波分量，其高频信号很弱，以至于小于观测误差.由于高阶噪声的影响，GRACE时变重力场在图像上表现为沿经度方向的条带，为了获取有用的重力信号，需要对球谐系数的高阶项进行平滑处理，提高信噪比.同时，重力场高阶球谐系数解表现出奇偶阶相关性(Swenson and Whar， 2006)，需要去除.本文采用去相关P3M6方法，即对m=6以上球谐系数的奇偶阶分别用三次多项式拟合，然后将拟合多项式从原始球谐系数中扣除.另外，本文采用各向同性高斯滤波器，平滑半径根据d=20000/n确定(周新，2012)，所用数据最高阶为n=60，空间分辨率为330 km.我们采取平滑半径为350 km的平滑因子进行高斯滤波.经平滑后，高阶噪声明显降低，信号也相应变弱.
为了消除季节性非地震因素的影响，我们采用差分法(周新等，2011；Chen et al., 2007)提取震后7年的重力变化.所用数据为2003年至2011年的9—11月数据，背景场取2003年和2004年9—11月的均值.例如，用2009年9—11月的重力场均值减去背景重力场即为震后五年的重力变化.图 1给出了GRACE监测到的2004年苏门答腊Mw9.3级地震的震后重力时空分布演化过程.
 |
图 1 GRACE监测到的2004年苏门答腊地震(Mw9.3)引起的震后重力变化Fig. 1 Post-seismic gravity changes induced by the 2004 Sumatra earthquake(Mw9.3)detected by GRACE |
震后断层左侧地区(下盘)重力呈上升趋势，变化平均速率为0.55 μGal/yr，但是在断层右侧(上盘)重力没有明显变化，平均年变化速率只有0.16 μGal/yr，在GRACE观测误差之内.断层两侧震后重力变化速率的差异性说明断层两侧存在较大的横向黏滞性差异.断层上盘所在区域地幔黏滞性较低，断层下盘所在区域地幔黏滞性较高.
为了更清晰地分析该地震导致的震后重力变化时间特征，本文取震中两侧重力变化最大和最小的两点A(95.0°E，1.5°N)和B(97.0°E，7.0°N)进行分析，计算其2003年1月至2011年12月期间的月重力相对于2003—2004年平均重力场变化的时间序列.这两个点由地震导致的重力变化信号最大，较 其他点有更高的信噪比.此外，由于该结果经过350 km 的高斯平滑，能反应该区域震后重力变化趋势.本文采用Chen等(2007)方法，即利用最小二乘拟合剔除年、半年周期的水文影响以及161天周期的S2潮汐波的影响.图 2给出了最小二乘拟合得到的地震前后A、B两点重力变化的时间序列.图 2a和2b分别表示A、B点重力变化时间序列，红色表示对GRACE数据进行去相关滤波和高斯平滑后的时间序列，蓝色表示最小二乘拟合后的结果，从图中可以 看出，A、B两点在地震前后均有明显的重力阶跃，但是A、B两点所在区域震后变化存在较大差异，震后A点所在区域重力有明显增加，而B点重力变化速率较A点要小.该结果与Chen等(2007)得到的震后重力变化趋势一致，但与王武星等(2011)的地震前(91.5°E，4.5°N)点(与A点接近)下降，震后一段时间上升，2007年9—10月份后又下降的结论不同，其原因是由于王武星等(2011)未对GRACE观测时间序列剔除季节变化的影响.
 |
图 2 GRACE监测的2004年苏门答腊地震(Mw9.3)前后A(95.0°E，1.5°N)、B(97.0°E，7.0°N)两点重力变化时间序列. 红色曲线为GRACE数据经过去相关滤波和高斯平滑后得到的时间序列，蓝色为在红色曲线的基础上，剔除季节性变化后的重力变化时间序列Fig. 2 The gravity change series of points A(95.0°E，1.5°N) and B(97.0°E，7.0°N)before and after of the 2004 Sumatra earthquake(Mw9.3). The red curves represent the gravity changes achieved through de-correlation filtering and Gaussian smoothing. The blue curves represent the gravity changes excluding the seasonal signals |
结果表明，A、B两点所在区域震后重力变化趋势相同，均有所增加，但是，地震后A点重力变化幅度较大，达到5.6 μGal，而B点相对变化较小，只有2.3 μGal.下盘区域的重力上升比上盘区域更显著，表明下盘区域黏滞性较高，地层相对柔软，其地层由于地震的作用持续变形.而上盘所在的地层接近陆地地壳特性，地层相对脆硬，震后受地震的影响小而导致的重力变化较小.
3 苏门答腊地震(Mw9.3)震后重力变化模拟Tanaka等(2006，2007)提出的黏弹球体位错理论是对弹性球体位错理论(Sun and Okubo， 1993; Sun et al., 2009)的扩展，即对黏弹性的本构关系经过拉普拉斯变化得到与弹性变形格林函数相同的形式，按照弹性位错理论的求解方法得到拉普拉斯域的变形格林函数，然后对其进行反拉普拉斯变换便可得到空间域的黏弹性变形.此外，Tanaka等(2006)利用复平面路径积分的方法解决了地球只能有限分层的缺陷，并且利用Poster-Widder公式克服了复平面路径积分数值积分精度难以保证的困难.该黏弹性球体位错理论可以同时考虑地球的可压缩性和无限分层的影响，是迄今为止较为完备的计算震后变形的理论之一，用其计算远场变形可保证较高精度(Tanaka et al., 2006)，但是计算近场变形时存在较大误差，图 3为分别利用黏弹性位错程序和弹性位错程序计算的2004年苏门答腊地震(Mw9.3)的同震重力变化空间分布，图 1a是依据Tanaka的黏弹位错程序模拟得到同震重力变化空间分布，其结果沿断层方向有很强的正负相间变化.图 1b是由弹性位错程序计算得到的同震重力变化空间分布，该分布在断层两侧呈现单一变化.图 1a与图 1b差异很大，且与周新(2012)利用球体位错程序计算的逆冲型地震导致的空间重力变化分布不一致.
 |
图 3 2004年苏门答腊地震(Mw9.3)引起的近场同震重力变化理论结果.黑色五角星为地震震中 (a)基于Tanaka等(2006)的黏弹位错程序的计算结果；(b)基于弹性位错程序(付广裕和孙文科，2012)的计算结果Fig. 3 The oretical near-field gravity changes caused by the 2004 Sumatra Earthquake(Mw9.3). Black pentagram is the epicenter (a)The result calculated by the viscoelastic dislocation program(Tanaka et al，2006); (b)The result calculated by the elastic dislocation program(Fu and Sun， 2012) |
黏弹位错程序近场计算误差较大的原因是其在计算地震导致的变形时把子断层视为点震源，这对于远场(子断层尺寸相对于震中距可以忽略)变形计算能保证其精度，但计算近场变形会导致较大误差甚至错误(Sun and Okubo， 1998; 付广裕和孙文科，2012).为了克服黏弹位错程序计算近场变形误差较大的缺陷，本文采用Sun和Okubo(1998)提出的有限断层细分方法对黏弹位错程序进行改造，很好地解决了黏弹位错程序近场计算误差较大的问题.图 4为分别利用弹性和改正后的黏弹位错程序计算得到的同震重力变化、垂直位移和对应变形之间的差异，用弹性位错程序和改造后的黏弹位错程序分别计算的同震重力变化和垂直位移平均差异仅为0.06%和0.1%，说明利用Sun和Okubo(1998)的有限断层细分方法可以很好地解决Tanaka等(2006)的黏弹性球体位错程序计算同震近场变形误差较大的问题.
 |
图 4 基于弹性位错程序和黏弹位错程序计算的2004年苏门答腊地震(Mw9.3)引起的同震重力和垂直位移变化及两种位错理论结果的差异 (a)利用弹性位错程序计算获得的同震重力变化；(b)利用黏弹位错程序计算获得的同震重力变化； (c)为(a)和(b)的差异；(d)、(e)和(f)是两种算法的同震垂直位移及其差异. Fig. 4 Co-seismic gravity changes and vertical displacement induced by the 2004 Sumatra earthquake(Mw9.3)calculated by the elastic spherical dislocation theory and the viscoelastic spherical dislocation theory， and differences between two theories in gravity changes and vertical displacement (a)Co-seismic gravity changes calculated by the elastic spherical dislocation theory;(b)Co-seismic gravity changes calculated by viscoelastic spherical dislocation theory;(c)The difference between(a) and (b);(d)、(e) and (f)are the vertical displacement calculated by two theories and their difference. |
为了验证改进的黏弹性位错震后格林函数的正确性，我们计算了四种独立震源在震源深度为32 km，地幔黏度为1×1018Pa·s的0~0.7年和0~7年的重力变化格林函数.图 5为0~7年的重力变化格林函数，当时间尺度较长时，几年甚至几十年，形变速率随着时间变化逐渐减小，将会在某个时间以后趋于稳定，即不再变化.
 |
图 5 震源深度为32 km时四种独立震源的重力变化格林函数在0~7年的变化趋势，
图中黑色圆点为弹性球体位错理论格林函数 (a)垂直倾滑；(b)水平走滑；(c)水平引张；(d)垂直引张.Fig. 5 Green Functions of gravity changes for four independent sources(32 km depth)from 0yr to 7yr. The black dot represents the Green Functions of gravity changes based on the elastic dislocation theory (a)Vertical slip;(b)Strike slip;(c)Horizontal tensile;(d)Vertical tensile. |
图 6为利用0~0.7年的重力变化格林函数，从图中可以看出，在一年内震后变形随着时间越靠近同震时刻，格林函数也逐渐逼近同震格林函数.
 |
图 6 震源深度为32 km时四种独立震源的重力变化格林函数在0~0.7年的变化趋势，
图中黑色圆点为弹性球体位错理论格林函数 (a)垂直倾滑；(b)水平走滑；(c)水平引张；(d)垂直引张.Fig. 6 Green Functions of gravity changes for four independent sources(32 km depth)from 0yr to 0.7yr. The black dot represents the Green Functions of gravity changes based on the elastic dislocation theory (a)Vertical Slip;(b)Strike Slip;(c)Horizontal Tensile;(d)Vertical Tensile. |
本文利用改正后的黏弹位错程序计算了2004年苏门答腊地震(Mw9.3)引起的震后重力变化，计算网格密度为0.5°×0.5°，采用的地球分层模型为Tanaka等(2006)的模型(表 1)；断层数据为Hoechner等(2008)利用GPS数据反演得到的(图 7).根据前人研究结果，该地区地幔黏滞系数介于1×1018~1×1019 Pa·s之间.王武星等(2011)利用GRACE数据得到苏门答腊地区地幔黏滞系数在1×1018 Pa·s左右；Tanaka等(2009)利 用谱元法模拟认为该地区黏滞系数为1×1019 Pa·s左右. 并且由于温度、岩性等的差异，海洋地区较陆地地区地幔黏滞性系数偏小(Cadek和Fleitout，2003)；苏门答腊地区海洋一侧黏滞系数要小于陆地一侧才能使得模拟震后重力变化与GRACE观测结果一致(王武星等，2011 ).因此，本文在利用黏弹球体位错程序进行模拟时断层两侧分别取不同的黏滞系数进行计算，发现当断层左侧(下盘)区域地幔黏滞系数取1×1018 Pa·s，断层右侧(上盘)地幔黏滞性系数取8×1018 Pa·s时，模拟结果与GRACE观测结果拟合度较高.需要指出的是：本文采用的断层两侧分别计算格林函数的方法是一个近似方法，其基本思想是按照断层一侧的地幔黏滞性构造构建层状地球模型，然后依此层状地球模型计算格林函数，这样的地球模型与实际地球模型存在差异，因此该计算方法并不严密.但是，由于格林函数表示的是地震震源在地表产生的变形，在断层模型与地表观测站之间的整个变形传播路径之上，地球模型并无误差，因此，我们的数据处理方法计算获得的格林函数，误差应该不大.关于地球横向不均匀构造对地震形变的影响的严格计算，需参照更为严格的基于三维不均匀地球的地震位错理论(Fu and Sun， 2010).
|
表 1 苏门答腊地区地球分层模型参数Table 1 Parameters of the layered Earth model at Sumatra region |
 |
图 7 2004年苏门答腊地震(Mw9.3)同震断层滑动模型(Hoechner et al., 2008)Fig. 7 The co-seismic slip model of the 2004 Sumatra Earthquake(Mw9.3)based on Hoechner et al.,(2008) |
此外，由于GRACE观测到的重力变化是质量迁移的综合响应，包含由于海底垂直变形导致海水重新分布引起的重力变化，因此要与GRACE观测结果对比，需要对模拟结果进行海水改正(周新，2012)，即将海底垂直变形导致的海水扰动看成不可压缩的布格层，计算重力扰动并对模拟值加以改正.海水负荷变形相对于海水扰动可以忽略不计.同时，由于我们对GRACE数据进行了P3M6去相关技术和350 km的高斯平滑处理，为了与其比较，我们 也对模拟结果进行相应的滤波处理.震后七年重力 变化如图 8.
 |
图 8 黏弹位错程序计算的2004苏门答腊地震(Mw9.3)引起的震后重力演化过程Fig. 8 Sketch The co-seismic slip model of the 2004 Sumatra Earthquake(Mw9.3)based on Hoechner et al.,(2008) |
模拟结果和GRACE观测结果比较，模拟值要比GRACE观测结果偏小1~2 μGal，这是由于导致震后变形的机制一般包含黏滞性、震后余滑和空隙流体回弹等(王武星等，2011)，观测值同时包含了这几种因素共同的影响，几种因素耦合在一起，且不一定都是正影响，而模拟结果只包含黏滞性的影响.理论重力变化峰值与GRACE卫星观测结果在空间分布和变化趋势上有很好的可比性，即断层下盘震后重力持续增加，断层上盘震后重力基本保持稳定，说明在震后较长时间尺度内黏滞性松弛是震后变形的主要因素.模拟结果表明，断层下盘震后重力变化平均速率约为0.52 μGal/yr，断层上盘变化平均速率约为-0.12 μGal/yr.断层上盘震后重力平均变化速率观测值与理论值基本一致，而断层下盘观测到的平均重力变化为0.16 μGal/yr，这是该区域地层接近陆地地壳特性，地层相对脆硬，震后受地震的影响而导致的变化较小，且GRACE存在观测误差导致的.从理论和观测上同时验证了断层两侧黏滞性有很大差异，断层上盘区域地幔黏滞系数大约为8×1018 Pa·s，断层下盘区域地幔黏滞系数大约为1×1018 Pa·s.
4 结论与讨论本文利用GRACE数据提取了苏门答腊地震(Mw9.3)导致的震后重力变化.2004至2005年的重力变化峰值约为-8.3~4.6 μGal(见图 1)，该结果与利用分层弹性球体位错计算程序计算得到的同震重力变化(张国庆等，2015)在空间分布和变化量级上基本一致，震后重力变化时间序列与Chen等(2007)利用GRACE数据提取该地震的重力变化(等效水高)趋势一致.震后重力变化空间分布观测结果表明，该地区断层两侧震后重力变化速率存在明显差异，断层左侧(下盘)区域变化速率约为 0.55 μGal/yr，而断层右侧(上盘)区域变化速率约为0.16 μGal/yr.说明该地区地幔黏滞性存在显著的差异，断层左侧(下盘)区域地幔黏滞系数小于断层右侧(上盘)区域地幔黏滞系数.
Tanaka等(2006)的黏弹球体位错理论由Sun和Okubo(1993)的弹性球体位错理论的基础上扩展得到，其配套计算程序在计算地震导致的近场震后变形时存在较大误差，因此我们采用Sun和Okubo(1998)的有限断层细分方法对黏弹球体位错计算程序进行改进，并利用改正后的黏弹格林函数与弹性格林函数(付广裕和孙文科，2012)分别模拟了同震重力变化和同震垂直位移.结果显示，两种计算程序得到的重力变化和垂直位移之间的差异分别为0.06%和0.1%，证明了我们对Tanaka等(2006)的黏弹性球体位错计算程序进行改造的正确性和必要性.同时，为了验证震后格林函数改正的正确性，我们计算了震源深度为32 km的四种独立震源引起的震后0~0.7年和0~7年的重力变化格林函数，结果表明，震后时间越靠近同震时刻，相应的格林函数也逐渐逼近同震格林函数；当震后时间尺度较长时，变化速率随着时间增加而逐渐减小，并将在 某个时间点以后趋于稳定，不再变化，这与黏滞性松弛理论在变化趋势上是一致的.
本文利用改正后的黏弹性位错程序计算了苏门答腊地震(Mw9.3)引起的震后重力变化，经过反复拟合发现，当下盘所在区域地幔黏滞系数取1×1018 Pa·s，上盘所在区域地幔黏滞系数取8×1018 Pa·s时，拟合结果与GRACE观测结果在峰值空间分布和震后变化趋势上基本一致，说明在较长时间尺度内地幔黏滞性是导致震后变形的主要因素.此外，模拟结果整体比观测结果小1~2 μGal.震后重力变化观测结果是孔隙流体回弹、余滑等因素对震后变形的影响与黏滞性效应耦合在一起的综合效应，但理论计算只考虑了地幔黏滞效应的影响，未考虑震后余滑、地幔孔隙流体回弹等因素的影响，这是模拟结果整体比观测结果小的原因所在.
致谢 感谢两位匿名评审对于本文的建设性意见.本文采用的GRACE卫星RL05数据是由德国地球 科学研究中心网站提供(http://icgem.gfz-potsdam.de/ICGEM/TimeSeries.html)| [1] | Broerse D B T, Vermeersen L L A, Riva R E M, et al. 2011. Ocean contribution to coseismic crustal deformation and geoid anomalies: Application to the 2004 December 26 Sumatra—Andaman earthquake. Earth Planet. Sci. Lett., 305(3-4): 341-349,doi: 10.1016/j.epsl.2011.03.011. |
| [2] | Cadek O, Fleitout L. 2003. Effect of lateral viscosity variations in the top 300 km on the geoid and dynamic topography. ,Geophys. J. Int. 152: 566-580. |
| [3] | Cambiotti G, Bordoni A, Sabadini R, et al. 2011. GRACE gravity data help constraining seismic models of the 2004 Sumatran earthquake. J. Geophys. Res., 116: B10403, doi: 10.1029/2010JB007848. |
| [4] | Cambiotti G, Sabadini R. 2013. Gravitational seismology retrieving Centroid-Moment-Tensor solution of the 2011 Tohoku earthquake. J. Geophys. Res., 118, 183-194, doi: 10.1029/ 2012JB009555. |
| [5] | Cannelli V, Melini D, Piersanti A, et al. 2008. Postseismic signature of the 2004 Sumatra earthquake on low-degree gravity harmonics. J. Geophys. Res., 113, B12414, doi: 10.1029/2007JB005296. |
| [6] | Chen J L, Wilson C R, Tapley B D, et al. 2007. GRACE detects coseismic and postseismic deformation from the Sumatra-Andaman earthquake. Geophys. Res. Lett., 34(13):L13302, doi:10.1029/2007GL030356. |
| [7] | Chen Y T, Huang L R, Lin B H, et al. 1979. A dislocation model of the Tangshan earthquake of 1976 from the inversion of geodetic data. Chinese J. Geophys. (in Chinese), 22(3): 201-215. |
| [8] | Chen Y T, Lin B H, Lin Z Y, et al. 1975. The focal mechanism of 1966 Hsingtai earthquake as inferred from the ground deformation observations. Chinese J. Geophys. (in Chinese), 18(3): 164-181. |
| [9] | Deng J S, Gurnis M, Kanamori H, et al. 1998.Viscoelastic low in the lower crust after the 1992 Landers, California, earthquake. Science, 282(5394): 1689-1692. |
| [10] | De Linage C, Rivera L, Hinderer J, et al. 2009. Separation of coseismic and postseismic gravity changes for the 2004 Sumatran earthquake from 4.6 yr of GRACE observations and modelling of the coseismic change by normal mode summation. Geophys. J. Int., 176(3): 695-714. |
| [11] | De Viron O, Panet I, Mikhailov V, et al. 2007. Retrieving earthquake signature in GRACE gravity solutions. Geophys. J. Int., 174, doi: 10.1111/j.1365-246X.2008. 03807.x. |
| [12] | Fu G Y, Sun W K, Fukuda Y, et al. 2010. Coseismic displacements caused by point dislocations in a three dimensional heterogeneous, spherically earth model. Geophys. J. Int., 183(2): 706-726, doi: 10.1111/j.1365-246X.2010.04757.x. |
| [13] | Fu G Y, Sun W K. 2012. Overall design and specific structures of the computing codes for coseismic deformations on a layered spherical earth. Earthquake (in Chinese), 32(2): 73-87. |
| [14] | Han S C, Shum C K, Bevis M, et al. 2006. Crustal dilatation observed by GRACE after the 2004 Sumatra-Andaman earthquake. Science, 313(5787): 658-662. |
| [15] | Han S C, Sauber J, Luthcke S B, et al. 2008. Implications of postseismic gravity change following the great 2004 Sumatra—Andaman earthquake from the regional harmonic analysis of GRACE intersatellite tracking data. J. Geophys. Res., 113: B11413, doi:10.1029/2008JB005 705. |
| [16] | Hao J L, Yao Z X. 2012. The coseismic displacement, strain and stress in the layered elastic model. Chinese J. Geophys. (in Chinese), 5(5): 1682-1694, doi: 10.6038/j.issn.0001-5733.2012.05.0 25. |
| [17] | Hoechner A, Babeyko A Y, Sobolev S V. 2008. Enhanced GPS inversion technique applied to the 2004 Sumatra earthquake and tsunami. Geophys. Res. Lett., 35: L08310, doi: 10.1029/2007G L033133. |
| [18] | Hoechner A, Sobolev S V, Einarsson I, et al. 2011. Investigation on afterslip and steady state and transient rheology based on postseismic deformation and geoid change caused by the Sumatra 2004 earthquake. Geochem. Geophys. Geosyst., 12: Q07010, doi: 10.1029/2010GC 003450. |
| [19] | Ju X L, Shen Y Z, Zhang Z Z. 2013. Antarctic ice mass change analysis based on GEACE RL05 data. Chinese J. Geophys. (in Chinese), 56(9): 2918-2927. |
| [20] | Jónsson S, Segall P, Pedersen R, et al. 2003. Post-earthquake ground movements correlated to pore-pressure transients. Nature, 424(6945): 179-183. |
| [21] | Ma X Q, Kusznir N J. 1994. Effects of rigidity layering, gravity and stress relaxation on 3-D subsurface fault displacement fields. Geophys. J. Int., 118(1): 201-220. |
| [22] | Ogawa R, Heki K. 2007. Slow postseismic recovery of geoid depression formed by the 2004 Sumatra-Andaman Earthquake by mantle water diffusion. Geophys. Res. Lett., 34: L06313, doi: 10.1029/2007GL029340. |
| [23] | Okada Y. 1985. Surface deformation caused by shear and tensile faults in a half-space. Bull. Seism. Soc. Am., 75(4): 1135-1154. |
| [24] | Panet I, Mikhailov V, Diament M, et al. 2007. Coseismic and post-seismic signatures of the Sumatra 2004 December and 2005 March earthquakes in GRACE satellite gravity. Geophy. J. Int., 171(1): 177-190. |
| [25] | Panet I, Pollitz F, Mikhailov V, et al. 2010. Upper mantle rheology from GRACE and GPS postseismic deformation after the 2004 Sumatra—Andaman earthquake. Geochem. Geophys. Geosyst., 11(6): Q06008, doi: 10.1029/2009GC002905. |
| [26] | Piersanti A, Spada G, Sabadini R, et al. 1995. Global post-seismic deformation. Geophys. J. Int., 120(3): 544-566. |
| [27] | Pollitz F F. 1996. Coseismic deformation from earthquake faulting in a layered spherical Earth. Geophys. J. Int., 125: 1-14 |
| [28] | Pollitz F F. 1992. Postseismic relaxation theory on the spherical earth. Bulletin of the Seismological Society of America, 82(1): 422-453. |
| [29] | Pollitz F F. 2003. Post-seismic relaxation theory on a laterally heterogeneous viscoelastic model. , Geophys. J. Int.155(1): 57-78. |
| [30] | Rundle J B. 1980. Static elastic-gravitational deformation of a layered half space by point couple sources. J. Geophys. Res.: Solid Earth (1978—2012), 85(B10): 5355-5363. |
| [31] | Rundle J B. 1982. Viscoelastic-gravitational deformation by a rectangular thrust fault in a layered Earth. J. Geophys. Res., 87(B9): 7787-7796. |
| [32] | Sabadini R, Piersanti A, Spada G. 1995. Toroidal/poloidal partitioning of global post-seismic deformation. Geophys. Res. Lett., 22(8): 985-988. |
| [33] | Sheu S Y, Shieh C F. 2004. Viscoelastic-after slip concurrence: A possible mechanism in the early post-seismic deformation of the Mw7.6, Chi-Chi (Taiwan) earthquake. Geophys. J. Int., 159(3): 1112-1124. |
| [34] | Spada G, Boschi L. 2006. Using the Post-Widder formula to compute the Earth's viscoelastic Love numbers. Geophys. J. Int., 166(1): 309-321. |
| [35] | Steketee J A. 1958. On Volterra's dislocations in a semi-infinite elastic medium. Canadian Journal of Physics, 36(2): 192-205. |
| [36] | Sun W K, Okubo S. 1993. Surface potential and gravity changes due to internal dislocations in a spherical Earth—I. Theory for a point dislocation. Geophys. J. Int., 114(3): 569-592. |
| [37] | Sun W K, Okubo S, Vaníek P. 1996. Global displacements caused by point dislocations in a realistic Earth model. J. Geophys. Res., 101(B4): 8561—8577, doi:10.1029/95JB03536. |
| [38] | Sun W K, Okubo S. 1998. Surface potential and gravity changes due to internal dislocations in a spherical Earth—II. Application to a finite fault. , Geophys. J. Int.132(1): 79-88. |
| [39] | Sun W K, Okubo S, Fu G Y, et al. 2009. 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. Geophys. J. Int., 177(3): 817-833, doi: 10.1111/j.1365-246X.2009.04113.x. |
| [40] | Swenson S, Wahr J. 2006. Post-processing removal of correlated errors in GRACE data. Geophys. Res. Lett., 33: L08402, doi: 10.1029/2005GL025285. |
| [41] | Tanaka T, Okuno J, Okubo S. 2006. A new method for the computation of global viscoelastic post-seismic deformation in a realistic earth model (I)—vertical displacement and gravity variation. Geophys. J. Int., 164(2): 273-289, doi: 10.1111/j.1365-246X.2005.02821.x. |
| [42] | Tanaka T, Okuno J, Okubo S. 2007. A new method for the computation of global viscoelastic post-seismic deformation in a realistic earth model (II)—horizontal displacement. Geophys. J. Int., 170(3): 1031-1052, doi: 10.1111/j.1365-246X.2007.03486.x. |
| [43] | Tanaka Y, Klemann V, Fleming K, et al. 2009. Spectral finite element approach to postseismic deformation in a viscoelastic self-gravitating spherical Earth. Geophys. J. Int., 176:715-739. |
| [44] | Wang H S. 1999. Surface vertical displacements, potential perturbations and gravity changes of a viscoelastic earth model induced by internal point dislocations. Geophys. J. Int., 137(2): 429-440. |
| [45] | Wang R J, Lorenzo Martín F, Roth F. 2003. Computation of deformation induced by earthquakes in a multi-layered elastic crust-FORTRAN programs EDGRN/EDCMP. Computers and Geosciences, 29(2): 195-207. |
| [46] | Wang R J, Lorenzo-Martín F, Roth F. 2006. PSGRN/PSCMP—A new code for calculating co-and post-seismic deformation, geoid and gravity changes based on the viscoelastic-gravitational dislocation theory. Computers and Geosciences, 32(4): 527-541. |
| [47] | Wang W X, Shi Y L, Sun W K, et al. 2011. Viscous lithospheric structure beneath Sumatra inferred from post-seismic gravity changes detected by GRACE. Sci. China Earth Sci., 54(8): 1257-1267. |
| [48] | Wang L, Shum C K, Simons F J, et al. 2012. Coseismic and postseismic deformation of the 2011 Tohoku-Oki earthquake constrained by GRACE gravimetry. Geophys. Res. Lett., 39: L07301, doi: 10.1029/2012GL051104. |
| [49] | Zhang G Q, Fu G Y, Zhou X. 2015. The evolution process of gravitational field after the Sumatra Mw9.3 earthquake from GRACE RL05 data. Journal of Geodesy and Geodymanics (in Chinese), 35(2):303-308. |
| [50] | Zhou X, Sun W K, Fu G Y. 2011. Gravity satellite GRACE detects coseismic gravity changes caused by 2010 Chile Mw8.8 earthquake. Chinese J. Geophys., 54(7): 1745-1749. |
| [51] | Zhou X. 2012. Coseismic deformation and slip characteristics from modern geodetic data for the 2011 Tohoku-Oki (Mw9.0) megathrust earthquake (in Chinese). Beijing: University of Chinese Academy of Sciences. |
| [52] | Zhou X, Sun W K, Zhao B, et al. 2012. Geodetic observations detecting coseismic displacements and gravity changes caused by the Mw=9.0 Tohoku-Oki earthquake. J. Geophys. Res., 117: B05408, doi: 10.1029/2011JB008849. |
| [1] | 陈运泰, 黄立人, 林邦慧等. 1979. 用大地测量资料反演的1976年唐山地震的位错模式. 地球物理学报, 22(3): 201-215. |
| [2] | 陈运泰, 林邦慧, 林中洋等. 1975. 根据地面形变的观测研究1966年邢台地震的震源过程. 地球物理学报, 18(3): 164-181. |
| [3] | 付广裕, 孙文科. 2012. 球体位错理论计算程序的总体设计与具体实现. 地震, 32(2): 73-87. |
| [4] | 郝金来, 姚振兴. 2012. 均匀弹性分层介质模型中的同震位移, 应变以及应力. 地球物理学报, 55(5): 1682-1694, doi: 10.6038/j.issn.0001-5733.2012.05.025. |
| [5] | 鞠晓蕾, 沈云中, 张子占. 2013. 基于 GRACE 卫星 RL05 数据的南极冰盖质量变化分析. 地球物理学报, 56(9): 2918-2927. |
| [6] | 张国庆, 付广裕, 周新. 2015. 利用GRACE卫星RL05数据提取苏门答腊Mw9.3级地震同震和震后重力变化. 大地测量与地球动力学,35(2): 303-308. |
| [7] | 周新, 孙文科, 付广裕. 2011. 重力卫星 GRACE 检测出 2010 年智利Mw8.8地震的同震重力变化. 地球物理学报, 54(7): 1745-1749. |
| [8] | 周新. 2012. 利用现代大地测量数据研究2011日本东北大地震(Mw9.0)的同震变形和断层滑动分布［博士论文］. 北京: 中国科学院大学. |