基于最大切应变的震前GRACE重力异常信息提取方法

宋冬梅 王 慧 单新建 王 斌 崔建勇

1)中国石油大学(华东)，海洋与空间信息学院，青岛 266580 2)中国石油大学(华东)，研究生院，青岛 266580 3)中国地震局地质研究所，北京 100029

0 引言

破坏性地震会造成国家经济建设和人民生命财产安全的巨大损失(Allenetal.，2019； Fanetal.，2019)。因此，地震监测对于抗震设防具有重要的现实意义(Tanetal.，2019)。为此，世界上各个地震频发的国家或地区均建立了多个地震监测台站(申俊等，2018)。然而，利用地面台网传统的地震监测手段只能获取有限区域范围内的地震前兆信息，难以从宏观上反映孕震过程(Paveletal.，2017； Ziegeretal.，2018)。20世纪后期，遥感技术快速发展，使得大面积观测和获取实时信息得以实现，这有效地弥补了定点观测等传统地震监测手段的不足(范一大等，2016； Zhaoetal.，2021)，遥感技术已经成为地震监测的重要手段(彭令等，2017； Piscinietal.，2017； Jiaoetal.，2018； Baiketal.，2019)。

Massonnet等(1994)采用卫星雷达干涉测量技术获得了1992年美国加利福尼亚州兰德斯地震震后的地表形变信息。单新建等(2014)利用InSAR和GPS测量技术获取了汶川地震垂直连续形变场，结果显示发震断层两侧的垂直形变衰减较快，横跨断裂带的形变量>30cm，宽度≤50km，且沿断层垂直形变高值区主要集中分布在发震断裂的南段、 中段和北端。以上提到的SAR技术和GPS测量虽然能够很好地监测地表的地壳运动，但是却难以探测地壳深部的物质运移和质量变化(王武星等，2014)。而地球重力场可反映地球的内部结构和组成，被视作是固体地球动力过程的历史再现(Sunetal.，2017)。GRACE卫星作为重力卫星技术的代表，其获得的重力场数据包含了地球的密度结构信息，可用于震前板块间和板块内部作用机制等方面的研究(Peidouetal.，2019； Velicognaetal.，2020)。Wahr等(1998)在GRACE卫星发射之前便已给出了时变重力场的基本理论，并进行了数值模拟计算。Sun等(2004)基于位错理论就GRACE卫星能否探测到同震重力变化开展了研究，得出结论： GRACE重力卫星能够探测出7.5级以上地震的同震形变信号。Han等(2006)基于GRACE RL04 Level-1数据对2004年苏门答腊地震进行了研究，获得了世界上第一个由重力卫星数据探测得到的同震重力变化。

虽然GRACE卫星能够探测到大地震的重力变化，但由于重力场模型中的位系数含有较大噪声，使得反演结果会出现严重的南北条带状波纹现象，因此许多学者提出如高斯滤波等多种滤波方法对GRACE卫星数据进行去噪(Guetal.，2017； 丁一航等，2018； 郭飞霄等，2018； Chaoetal.，2019； 王陈燕等，2019； 刘潇，2021)。然而，进行滤波处理后发现，尽管噪声可被去除，但同时也造成有用信息的丢失，使得地震重力场异常信息的反演精度下降(Peidouetal.，2020)。针对目前基于GRACE卫星数据的重力异常信息提取技术的不足，本研究提出了一种基于最大切应变的震前重力异常信息提取方法，该方法在抑制噪声的同时，减少了有用信息的损失。本文以2008年汶川地震和2015年尼泊尔地震为例，基于GRACE卫星数据反演的重力最大切应变进行了断裂带的震前构造活动分析。

1 研究区域及相关数据介绍

1.1 研究区域与实验震例

青藏高原位于亚洲大陆中部(26°00′～39°47′N，73°19′～104°47′E)，是世界上海拔最高、 面积最大且地形最复杂的高原(张培震等，2003)。新生代以来，印度板块与欧亚板块通过碰撞使得青藏高原大规模地壳隆起，并由此产生了厚约80km的地壳(Cuietal.，2001)。由于地壳内的升温导致其内部介质呈较强的黏塑性，且高原持续受到印度板块向N的推挤作用力，使得下地壳不断驱动着上地壳脆性层的运动与变形，由此进一步形成了青藏高原地区复杂的孕震环境(李吉均等，1998)。

青藏高原内部存在较多断裂带(图 1)，是全球陆地强震爆发的主体区。21世纪以来，该区域的地质构造运动十分活跃，且强震频发(张培震等，2003； 徐锡伟等，2017)。据统计，在2000—2021年期间，该区域共发生了43个6级及以上强震，其中包括2次7.5级以上的大地震，即MW7.9 汶川地震和MW7.8 尼泊尔地震，各自对应的发震断裂带为龙门山断裂带和南喜马拉雅造山带，如图 1 所示(1)http： ∥earthquake.usgs.gov/earthquakes/search/。。其中，龙门山断裂带处于青藏高原东部，是一条呈NE-SW向延伸的构造带，长度>500km，宽30～50km(Dirksetal.，1994； Xuetal.，2009； Zhangetal.，2009； 孙闯，2017)。喜马拉雅主碰撞造山带是一条向S突出的EW向弧形构造带，长约2500km，宽300～500km，属于青藏高原的南边界，并位于印度板块和欧亚板块的交界处(Burgetal.，1997)。需要说明的是，汶川地震和尼泊尔地震爆发后的3个月内，其发震断裂带上还发生了多个6级以上地震，各地震的详细信息见表1。

图 1 龙门山断裂带(黄色矩形框①)和南喜马拉雅造山带(黄色矩形框②)在青藏高原上的位置Fig. 1 The location of the Longmenshan fault zone(yellow rectangle ①) and the southern Himalayan orogenic belt(yellow rectangle ②)in Tibetan plateau.F1喜马拉雅主冲带； F2喀喇昆仑-嘉里断裂带； F3玛尼-玉树-鲜水河断裂带； F4昆仑-玛沁断裂带； F5阿尔金-海原断裂带； F6金沙江-红河断裂带。上述断裂带将青藏高原划分为6个不同的活动块体(Taylor et al.，2009)，分别为拉萨块体 (B1)、 羌塘块体(B2)、 巴颜喀拉块体(B3)、 柴达木块体(B4)、 祁连块体(B5)和川滇块体(B6)

1.2 数据

1.2.1 GRACE数据

GRACE卫星由德国空间飞行中心DLR和美国宇航局NASA联合设计研制，并于2002年3月成功发射(Inceetal.，2019)。GRACE由2颗低轨卫星组成，工作模式采用卫星跟踪卫星模式，包括低低模式(SST-LL)和高低模式(SST-HL)。本研究所用的重力数据是由美国得克萨斯大学空间研究中心(CSR)于2018年4月公开发布的GRACE RL06版本Level-2产品(2)http： ∥www2.csr.utexas.edu/grace/RL06.html。，该产品包含了2002年4月—2017年7月共计163个月的重力场数据，重力场模型的最大阶次为96，其时间分辨率为1个月，空间分辨率约为300km。需要说明的是，该数据已经消除了由潮汐及非潮汐引起的大气和海洋质量变化方面的影响(Swensonetal.，2008)。由于GRACE为极轨卫星，卫星轨道倾角较大，使得GRACE重力场模型的2阶项精度较低，故本文采用CSR提供的卫星激光测距(SLR)的测量值对2阶项进行了替换(Chengetal.，2004)。

1.2.2 GLDAS数据

全球陆地数据同化系统(GLDAS)由美国宇航局(NASA)、 美国国家环境预报中心(NCEP)和国家海洋大气局(NOAA)联合研制(Spennemannetal.，2015)。GLDAS采用数据同化技术将卫星观测数据和地表观测数据整合到统一的模型中。模型以网格数据的形式集成了多种陆地表面场信息，如全球的降水量(降雨和降雪)、 蒸散(土壤水分蒸发和植被蒸腾)、 地表径流、 地下径流、 土壤湿度、 地表温度和地表热流等。本研究采用的GLDAS-2.1模型数据(3)https： ∥ldas.gsfc.nasa.gov/gldas/。包括地下0～2m土壤含水量和雪水，其时间分辨率为1个月，空间分辨率为1°。

2 方法

为尽可能地凸显出构造活动信息，本文提出了一种基于GRACE卫星数据的震前重力异常信息提取方法，通过最大切应变对重力异常信息进行表征。该方法在抑制GRACE重力卫星数据中存在的南北条带误差的同时，还避免了与构造活动相关的重力信号的丢失。该方法的详细流程如图 2 所示。

图 2 本文方法的流程图Fig. 2 Flow chart of the proposed method.

本方法具体包括3个步骤：

(1)去除水文信息季节性变化的影响。GRACE卫星观测到的重力场变化主要由陆地水储量变化、 构造变形和地下物质流动等因素引起，其中陆地水储量变化引起的重力变化相对较大(Deggimetal.，2021)。因此，首先需要去除由水文的季节性变化引起的重力场变化，得到去除水文干扰的扰动引力位增量。具体过程见2.1节。

(2)计算最大切应变。本研究通过引入重力最大切应变来描述与构造活动相关的信息。在步骤(1)的基础上，计算了扰动引力位增量的2阶梯度，并进一步根据重力应变张量获得了重力最大切应变，以凸显出与构造活动相关的重力变化。其中，扰动引力位的2阶梯度对南北条带误差有明显的压制作用。具体过程见2.2节。

(3)基于偏移指数K进行震前异常信息时空分析。为进一步探究震前重力异常的时空演变规律，需要将最大切应变在时间域上进行标准化，得到偏移指数K，并在空间域上对K值的分布特征进行分析。具体过程见2.3节。

2.1 去除水文信息

GRACE时变重力场可以通过球谐系数恢复。球谐系数实际上是扰动位常微分方程的解，满足以下拉普拉斯方程(Eshaghetal.，2012)：

(1)

然而，GRACE重力卫星观测到的地壳表层质量变化信息是由地质构造活动以及水文的季节性变化共同作用所引起的。同时，陆地水储量的季节性变化(如四川盆地在夏季降雨量偏多)引起的重力变化比较显著(王武星等，2014)。因此，本研究利用GLDAS水文模型提供的土壤水和雪水数据可实现部分水文信号的分离。为此，首先将GLDAS数据进行球谐系数展开并截止到96阶(与GRACE的阶数相同)，获得水文球谐系数Cgldas和Sgldas，并与GRACE数据的球谐系数Cgrace、Sgrace分别作差，得到去除土壤水和雪水影响后的球谐系数Ccorrect和Scorrect。然后，对上述得到的相邻月份的球谐系数进行差分处理，以进一步滤除水文信号的季节性变化的影响：

(2)

设ΔCti、 ΔSti和ΔCti+1、 ΔSti+1为经过差分处理后的相邻月差分球谐系数，将它们分别代入式(1)中即可得到各自对应的扰动引力位增量(即重力变化量)ΔTi和ΔTi+1。

2.2 重力最大切应变

去除水文干扰以后的重力变化包含地壳运动、 变形和物质运移活动信息。为最大限度地凸显出与构造活动相关的重力异常信息，本研究引入重力最大切应变。该应变通过重力应变张量导出，而重力应变张量是基于拉格朗日方法由扰动引力位增量的2阶梯度计算得到的，其对南北条带误差有明显的压制作用(Lietal.，2015； Daietal.，2016)。

(3)

(4)

上述重力应变张量E的最大特征值与最小特征值之差即为最大切应变MSH：

MSH=λmax-λmin

(5)

MSH在地质学中的物理意义为地壳某点处在应力作用下产生的切应变的最大值(许才军等，2006)，因此本研究将重力最大切应变作为描述构造活动强弱的表征量。

2.3 重力异常的刻画

为进一步提取出震前异常变化，更好地认识断裂带的孕震过程，本文借助偏移指数K对最大切应变时间序列的数据进行刻画。通过设定偏移指数K的阈值，实现各像素点位置处不同时间序列的异常探测。偏移指数K的计算公式为

(6)

式中，(x，y)表示像素位置，t为序列数据的时间维，MSH为重力最大切应变值，μMSH为MSH数据序列的均值，σMSH为MSH数据序列的标准差。

按照统计学的3倍标准差原则，当一组数据中的某个数值与均值的差值大于3倍标准差时，则这个数据被判定为异常数据(李光强，2009)。因此，根据式(6)可知，>3的K值即为最大切应变的异常数据。需要注意的是，偏移指数值的计算针对图像中各像素点在其时间序列上的变化，与空间域上的其他像素点无关。此外，在对各像素点的时间序列进行异常探测的基础上，还对每个像素点的K值在空间域上的分布特点进行观察，从而得到对异常时空演变规律的认识。

3 实验结果与分析

本研究以汶川地震和尼泊尔地震为例，根据GRACE扰动引力位增量的2阶梯度计算得到了重力应变张量，进而获得各像素在震前2a的重力最大切应变时间序列，并进一步通过偏移指数K探测时间序列中存在的异常。基于此，本研究在断裂带观测尺度上完成了震前重力异常时空变化特征的分析，并以此探讨震前断裂带可能存在的构造活动。此外，为验证重力异常为震前特有的现象，本研究还开展了非震年份的同期对比实验。其中，非震年份的偏移指数值是利用GRACE卫星月数据无缺失值的2005—2010年的同期数据计算得到的偏移指数K值的平均值，详细计算过程为： 基于重力最大切应变分别计算断裂带区域中每个0.5°×0.5°网格上在2005—2006年、 2007—2008年和2009—2010年的K值，然后可得到这3个K值的平均值K0，最后选取与震期相同月份的K0表征非震年份的同期状态。

3.1 汶川震例

考虑到2008年5月12日、 5月25日及8月5日的MW6.1 、MW6.1 和MW6.0 3次余震和汶川MW7.9 主震皆发生在龙门山断裂带上，且发震时间相隔在3个月以内，因此将以上4个地震视作一个地震序列，详细信息见表1。图 3 显示了震前2a龙门山断裂带区域的重力异常值(>3的K值)总和的变化曲线。从中发现，在汶川地震前3个月即地震序列震前半年，断裂带上出现了重力异常值迅速上升的现象，并达到了2年内的最大值，由此推测在2008年2—4月期间，断裂带区域的构造活动尤为活跃。为了更加清晰直观地获取震前重力异常演变规律，本研究还显示了重力异常震前半年的时空分布情况，如图 4 所示。图中最大的K值高达4.6，如图4b 中的紫色框所示。

表 1 本研究中各震例的详细信息(数据来源： USGS官网)Table1 Details of each earthquake case involved in this study(data source： USGS)

图 3 龙门山断裂带上震前2a的重力异常值总和Fig. 3 The sum of gravity anomalies on the Longmenshan fault zone two years before the Wenchuan earthquake.红色虚线表示 MW7.9 汶川地震的发震时间

图 4 龙门山断裂带的偏移指数K在震前的时空变化结果Fig. 4 Pre-earthquake spatio-temporal variation of the offset index K on the Longmenshan fault zone.该图像的像素是对去除水文影响后的扰动引力位进行插值所得的0.5°×0.5°的结果。图幅中的最大值位于紫色矩形框所在位置。红色标注日期为发震期，红色圆点表示发震位置

图 4b 表明，在汶川MW7.9 地震及其2个余震(汶川MW6.1 地震和青川MW6.1 地震)发震前的1～3个月，也即为青川MW6.0 地震震前的4～6个月间，龙门山断裂带的中段和北段出现了大面积异常区，且异常区的面积为该地区上个时间段(即图4a)异常面积的3倍。由此可知，此时龙门山断裂带的中段和北段出现了显著的重力异常现象，结合本研究计算得到的该时段的GRACE重力变化值为正值这一事实，可推测出该时间段里龙门山断裂带已进入应力累积过程，积累了大量的应变能，失稳性已达到较高水平。而在此之前，仅在龙门山断裂带的南部区域出现少量异常值(图4a)。值得注意的是，地震序列的震中位置皆在上述显著异常区内，其中MW6.1 和MW6.0 青川地震的震中位置距研究区内最大切应变的最大K值所在的位置均在50km范围内。

图 4d 为汶川MW7.9 地震震后1个月内的重力异常值的空间分布，从图中可知龙门山断裂带上再次出现了明显的重力异常现象，主要表现为高异常值沿断裂带走向分布，并多出现在断裂带北部。由图4d 再结合图4e 中的重力异常值的空间分布可推测，此时龙门山断裂带的主要构造活动向N迁移，而1个月后在龙门山断裂带的东北端发生的青川MW6.0 地震验证了我们的推测。这也表明相比于单个地震而言，针对地震序列的时空分析可以更加全面地认识断裂带的构造活动特征。图 5 显示，在非震年份的同期中，龙门山断裂带上出现的最大K值为2.8，并未出现重力异常。

图 5 龙门山断裂带在非震年份的同期偏移指数K的时空变化结果Fig. 5 Temporal and spatial variation results of the simultaneous offset index K of the Longmenshan fault zone in non-earthquake years.图中圆点的含义与图 4 一致

图 6 南喜马拉雅造山带上震前2a的重力异常值总和Fig. 6 The sum of gravity anomalies in the southern Himalayan orogenic belt two years before the MW7.8 Nepal earthquake.红色虚线表示 MW7.8 尼泊尔地震的发震时间

3.2 尼泊尔震例

尼泊尔地震序列包含了6个强震，这6个地震的发震时间间隔不到1个月，详细信息见表1。本研究计算了震前2年南喜马拉雅造山带区域的所有重力异常值总和，计算结果见图 6。从图中可发现，在震前6个月，区域上出现了重力异常值显著升高的现象，并达到了2013—2015年间的最大值。图 7 显示了震前6个月南喜马拉雅造山带上重力异常值的空间分布。

图 7 南喜马拉雅造山带上的偏移指数K的震前时空变化结果Fig. 7 Pre-earthquake spatio-temporal variation of the offset index Kin the southern Himalayan orogenic belt.标注*的时间段内存在GRACE原始数据缺失的现象，即无2014年12月的GRACE数据

如图7b 所示，在MW7.8 尼泊尔地震的震前6个月，南喜马拉雅造山带区域出现了大面积异常区，占所在图幅面积的三分之一。该异常区全部围绕着6个震中位置展布，并且出现最大K值的位置距离MW7.8 尼泊尔地震的震中约50km。上述异常现象表明，这6个震中附近区域皆存在较高的应力作用。由此推测，在2014年10月—2015年1月期间，印度板块对青藏高原的挤压作用明显增强，断裂带两侧的构造活动异常活跃。如图7c 所示，沿着南喜马拉雅造山带出现了多个异常值，且尼泊尔地震序列里的6个震中皆位于高异常值出现的位置，结合该时段GRACE数据中的重力变化值为正这一事实，可以推断出此时的断裂带上已完成了应力积累。而在震时断裂带及其以南区域的GRACE重力变化为负值，由此推测断裂带所受的应力达到临界值后产生了破裂，前期积累的应力得到释放，并最终导致发生了尼泊尔地区的6个地震。与地震期不同的是，在非震年份的同期中，南喜马拉雅造山带上存在的最大K值为2.6，这表明该区域并未出现重力异常现象(图 8)。

图 8 南喜马拉雅造山带在非震年份的同期偏移指数K的时空变化结果Fig. 8 Temporal and spatial variation results of the simultaneous offset index K of the southern Himalayan orogenic belt in non-earthquake years.图中圆点的含义与图 7 一致

综上所述，基于地震序列能够更好地揭示断裂带在震前的构造活动过程。研究发现，震前往往伴随着强烈的断裂带构造活动，应力集中现象明显，具体表现为研究区中异常的最大值出现在距离震中约50km范围内，且该异常发生的时间往往在震前半年，这可能是一种地震前兆信号。特别值得注意的是，不同于以往研究中将GRACE RL06数据60阶以后的球谐系数直接舍弃或者对其进行简单滤波处理，本研究所采用的重力最大切应变的计算方法较好地利用了该卫星数据的中高阶信息，与单纯借助重力值本身相比更能灵敏地探测到地球内部的密度变化信息，从而更好地凸显出了由构造活动所引起的震前重力异常变化。

4 讨论

以往对震前重力异常的研究方法通常直接基于重力值本身(即扰动位增量)的变化进行分析，但这种方法需要首先利用不同半径的高斯滤波对南北条带噪声进行消除，容易造成数据的失真与信息的丢失(Hasegawaetal.，2011； Zhaoetal.，2015)。而本研究所使用的重力最大切应变是在扰动位增量的基础上进一步通过计算其2阶梯度所得到的，该梯度可以在不损失重力信息的同时又能够达到抑制噪声的效果(Lietal.，2015； Daietal.，2016)。为比较以上2种方法对构造活动探测能力的差别，选取龙门山断裂带为研究区，对基于最大切应变的异常提取方法和传统的基于扰动位增量的异常提取方法所得的汶川地震震前2a时间序列里的总K值进行了比对分析。该处的总K值为研究区内每个像素的K值的总和。实验结果如图 9 所示。

图 9a 显示了基于最大切应变时间序列的研究区总K值的变化情况，可以发现>3的总K值在震前3个月内迅速增大，结合此时K值的分布与断裂带的空间展布基本一致这一事实(图4b)，上述现象极有可能是断裂带强烈的构造活动所引起的。图9b、 9c和9d显示了基于扰动引力位增量时间序列的总K值变化，其各自对应的高斯滤波平滑半径分别为250km、 300km和500km。在这3张图中可以发现，除了以250km为平滑半径的高斯滤波所得的异常提取结果在震前2a有>3的K值出现外，其余半径的滤波结果均未出现K>3的情况。对>1和>2的K值变化曲线进行分析可以发现，经过平滑半径为250km和300km的高斯滤波处理后，扰动位增量的总K值在震前半年出现了较为明显的局部增长现象(图9b，c)，但随着滤波半径的进一步增大，该现象随之消失，如图9d 所示，这说明随着滤波半径的增大，在去噪强度增加的同时也造成了信号的丢失。

图 9 汶川地震前2a龙门山断裂带上最大切应变和扰动位的K值时间序列Fig. 9 Total K value of maximum shear strain and disturbance potential in the Longmenshan fault zone two years before the Wenchuan earthquake.a 基于最大切应变震前2a的时间序列计算的龙门山断裂带区域里所有像素的K值(K>1、 2、 3)总和。b—d 基于重力值本身(相邻月份扰动位增量)在震前2a的时间序列计算的龙门山断裂带区域里所有像素的K值(K>1、 2、 3)总和，其各自对应的高斯滤波平滑半径分别为250km、 300km和500km。图中红色箭头表示汶川地震的发震时间。由于K值是由GRACE卫星的3个相邻月重力场数据计算得到，为了便于表示，在图中的时间坐标轴上仅标注这3个月中最后1个月的时间

表 2 震前2a最大切应变和扰动位时间序列中各自的最大K值Table2 Maximum K value in the time series of maximum shear strain and disturbance potential two years before the Wenchuan earthquake

此外，单个像素的重力最大切应变时间序列里的最大K值为4.6，出现在震前3个月内。而扰动位增量时间序列的最大K值出现的时间在发震期的1a以前，且由滤波半径为300km和500km的高斯滤波处理后的扰动位增量时间序列里的最大K值皆<3，详细信息见表2。此外，由本文方法计算得到的最大K值位于龙门山断裂带的北端，相比于经典方法得到的最大K值，该值所在的位置最接近青川地震的震中。需要说明的是，表2 中最大K值所在像素的颜色取决于其值的大小，颜色标尺参照图 5。

综上所述，相比于常用的扰动位方法，本研究基于最大切应变的方法对重力场的异常信息更为敏感，因此对与构造活动相关的重力异常信息具有较强的探测能力。

5 结论

GRACE卫星能够探测到地震重力场的变化，因此被广泛应用于地震监测研究中。然而，以往学者们在提取地震的重力变化信息时，均采用不同方法对GRACE卫星数据进行滤波处理，造成了重力场中高阶信息的损失，因此很难提取到8级及以下地震的重力异常。对此，本研究提出了一种基于最大切应变的震前重力异常信息提取方法，该方法既能够避免重力信息的损失，又能够达到抑制噪声的效果，以凸显与构造活动相关的重力异常信息。本文以汶川地震和尼泊尔地震为例，分析了龙门山发震断裂带的震前重力异常的时空分布特征。主要研究结论如下：

(1)本文所提出的基于最大切应变的方法比传统方法对重力场的异常信息提取能力更强，因此能够提取到汶川地震和尼泊尔地震这2个8级以下地震的震前异常信息，这为重力数据应用于地震监测研究提供了一个新方法。

(2)在震前，重力异常值呈现出沿断裂带分布的现象，而且大范围地存在于震中附近，最大异常值出现在距离震中约50km的范围内，且上述显著异常发生的时间为震前半年内。

(3)通过地震期与非震期的对比分析可以发现，断裂带在震前出现了大面积的重力异常，而非震期则表现得较为平静，未出现明显的重力异常，表明断裂带在震前的构造活动比非震期更为活跃。