Wang等(2012)提出的重力梯度的计算公式如下:

式中:V为地震前后重力场变化; ΔC_nm和ΔS_nm是n阶m次的球谐系数变化值; P_nm为勒让德函数; G为引力常数; M为地球质量; a_e为赤道半径; r、θ、λ为地心距离、余纬、经度。

在地方指北坐标系下,重力梯度进一步可以推导为:

式中:V_r、V_θ、V_λ为V关于r、θ、λ的一阶导数; V_θθ、V_λθ、V_rr、V_rλ、V_θλ为V的二阶导数。

本文选择差分法消除平均引力及周期性信号影响获取同震重力变化后,通过式(2)计算获得重力梯度,再选用高斯滤波对条带噪声进行平滑处理。差分法是在地震前后各选取一段时间GRACE月模型数据求取其平均值并做差,以求出地震前后因地震引起的同震阶跃值。周期变化信号是GRACE卫星所得的时变重力场信号的重要组成部分,其伴随重力场变化长期存在,如年际信号等,所以采用差分法对各时段的平均场重力进行差分处理,从而较好地减弱周期信号的影响。

Steketee(1958)最早在地震学中引入位错模型,并迅速建立起基于断层研究的地球位错理论。Okada(1992)推导得出基于半无限空间均匀介质地球模型的位错公式。Wang等(2006)对该公式进行改进,提出考虑到层状黏弹性介质的矩形位错理论,独立研发了PSGRN/PSCMP计算源程序。Sun等(2010)提出球体位错模型理论,能够高精度解译200 km外的远场同震重力变化。

本文采用Wang等(2006)提出的黏弹性半空间层状位错模型及其相应的配套计算程序PSGRN/PSCMP,对2007年明古鲁M_W8.4大地震的同震重力变化作了正演计算。对通过上述程序计算得到重力异常后进行球谐展开,将得到的系数带入式(2)即可计算得到重力梯度。为使位错模型计算结果与GRACE卫星数据计算结果具有一致的空间分辨率,将位错正演结果进行球谐截断至GRACE卫星数据相同阶次,并采用相同滤波半径的高斯滤波进行处理。

图2给出了使用GRACE卫星数据计算得到明古鲁M_W8.4地震前后12个月不同重力梯度变量的重力梯度均值变化。从图中可以看到,在震中两侧重力梯度变量均具有非常显著的正、负的变化。各梯度分量重力变化的极大值、极小值也均出现在震中两侧附近。从图2d～f中,还能观测到较明显的高斯滤波未消除的南北异常条带噪声。若消除异常条带噪声,这将极大地减弱地球的物理信号,从而使空间分辨率下降,干扰滤波分析,因此本文后续未作进一步的滤波处理。

为了与 GRACE卫星数据计算结果进行比较分析,本文使用黏弹性半空间层状位错模型进行正演计算,得到了明古鲁M_W8.4地震的同震重力梯度变化分布(图3)。由图3可知,在断层的上盘隆起带中,各重力梯度变量存在正向变化,观测到的重力梯度有显著升高的迹象,而在断层的上盘俯冲带存在重力负极值区。

图2 用GRACE卫星数据计算得到的重力梯度变化空间分布
Fig.2 Spatial distribution of the data signals of gradient variables calculated by GRACE satellite data

图3 用位错模型正演得到的重力梯度变化空间分布
Fig.3 Spatial distribution of the data signals of the gradient variables by the forward modeling of the dislocation model

图4给出了基于GRACE卫星数据和位错模型的重力梯度变化差值分布。根据图4可知,6个梯度变量的重力梯度差值都在同一量级且数值都较小。总体看,正、负重力梯度差异的空间分布范围类似,这表明 GRACE卫星能够探测到地震重力梯度信号,但在具体数值上,用GRACE卫星探测得到的重力梯度变化与用位错模型的正演计算结果有一定差别。 需要说明的是,高斯滤波半径的选择,会影响梯度变量南北条带噪声的削弱情况,如果适当增大滤波半径,可以进一步消除异常条带,但会严重削弱地球真实物理信号,从而导致空间分辨率降低,不利于滤波分析。本文选用的滤波半径为300 km,导致V_yy、V_zy、V_zz梯度变量差值存在南北条带噪声。

图4 基于GRACE卫星数据和位错模型的重力梯度变化差值空间分布
Fig.4 Spatial distribution of the difference of coseismic gravity gradient variations between GRACE satellite data and the dislocation model

表2 基于GRACE卫星数据和位错模型的重力梯度最大值、最小值、均值、标准差及差值
Tab.2 Maximum value,minimum value,average value,standard deviation and difference of forward results of gradient variables based on the GRACE satellite data and the dislocation model

表2为使用GRACE卫星数据和位错模型所得的各梯度变量的重力梯度最大值、最小值、标准差及差值。由表可知,各梯度变量最大值和最小值的数量级基本一致,差值在同一数量级或更小。进一步对这两种方法的处理结果进行相关性分析,得到GRACE计算结果与位错模型正演结果的相关系数矩阵及相关性见表3,相关性随区域的变化如图5所示。

表3 基于GRACE卫星数据和位错模型的各梯度变量相关性
Tab.3 Correlation between results of gradient variables by the GRACE satellite data and the dislocation model

图5 基于GRACE卫星数据和位错模型的重力梯度各梯度变量相关性空间分布
Fig.5 Spatial distribution of the correlation between the gradient components based on GRACE satellite and dislocation model

根据皮尔逊设计的统计指标(谢海花等,2021; 陈天异,2021)可知,相关系数0.00～±0.30是微相关,±0.30～±0.50是实相关,±0.50～±0.80是显著相关,±0.80～±1.00是高度相关。因此,由于GRACE计算结果和位错模型正演结果的6个梯度变量相关系数为0.38～0.51(表3),可以确定两种方法处理地震数据的结果具有相关性,且相关关系的密切程度较高。为衡量相同区域使用这两种方法得出的重力梯度之间的相关性,计算各梯度变量之间的皮尔逊相关系数,图5给出了两种方法各梯度变量相关性的区域分布,在(95°～105°E,0°～10°S)范围内相关性较高,震中西南方向的海水域相关性效果比较好。为进一步探究距震中远近与相关性的关系,选取(100°～102°E,2.4°～4.4°S)为初始范围,以1°为步长向四周扩散,对得到的13个形状相同大小不一的矩形进行相关性分析,依次得到关于GRACE卫星数据和位错模型正演结果的各梯度变量的相关系数,各梯度变量的相关系数随震中距的变化情况如图6所示。

图6 基于GRACE卫星数据与位错模型的各梯度变量的相关系数随震中距的变化
Fig.6 Variation of the correlation coefficient between the gradient by the GRACE satellite data and the gradient by the dislocation model along with epicentral distance

以V_xx分量为例,在包括了震中的初始区域范围中,两种方法得到的各梯度变量相关性最大,为0.98。随着区域逐渐扩大,相关系数值越来越小,最终整个区域的相关系数约为0.50。6个梯度变量相关系数变化幅度较为接近,除V_zx外,其余分量均表现出距离震中越近的区域相关性越高的特点,震中区域的相关性均在0.8之上; V_zx分量在震中区域相关性稍低,但也高于0.6。上述结果表明:GRACE卫星数据探测结果与位错模型正演结果在震中区域具有强相关性,能够用于地震的同震重力梯度变化研究。