根据中国地震台网测定，2017-08-08四川省西北部九寨沟县发生M_S7.0地震，震中位于33.2°N、103.82°E，震源深度约20 km。中国地震局地球物理研究所、美国地质调查局和相关文献计算结果均表明，本次九寨沟地震为走滑型地震^[1-2]。九寨沟地震周边地形和构造复杂(图 1)，近10 a来，其南部和北部先后发生2008年汶川M_S8.0地震、2013年芦山M_S7.0地震、2013年岷县-漳县M_S6.6地震和2016年门源M_S6.4地震等破坏性强震。本次地震的发震断层为树正断裂^{[1, 3]}，其位于塔藏断裂、岷江断裂和虎牙断裂西北延长线的交会区域，也是东昆仑断裂带由走滑到虎牙断裂逆冲构造的过渡区域^[4]。

地震孕震过程中的质量迁移必然伴随着重力场的改变，震前重力场变化可以由GRACE重力卫星和地表重力观测得到^[5-6]。通过绝对重力仪和相对重力仪组合的方式可以进行流动重力观测，不同期观测得到的重力变化特征可以作为强震发生前的可靠前兆^[7]。目前利用流动重力资料判断地震发生的可能地点和时间主要依据重力场时空动态演化特征和相关经验信息。识别重力场动态演化规律时，目标区域经常会出现多个“正负重力变化过渡带”等震前重力场变化特征，但受限于目前对地震孕育机制的认识和观测点位分布等因素，多数情况下难以对其进行取舍。本文尝试通过将重力变化数据进行球谐展开并滤波，压制局部重力场变化特征，突出主要演化规律，以便于目标特征的识别。

1 重力变化球谐展开方法

空间任意一点(θ, λ)的引力位球谐展开表达式^[8]为:

$ \begin{gathered} V=\frac{G M}{r} \sum\limits_{n=0}^{\infty}\left(\frac{a}{r}\right)^{n} \sum\limits_{m=0}^{n}\left[\bar{C}_{nm} \cos m \lambda+\right. \\ \left.\bar{D}_{nm} \sin m \lambda\right] \bar{P}_{n}^{m}(\cos \theta) \end{gathered} $

(1)

式中，λ为经度，θ为余纬，G为万有引力常数，M为地球质量，r为点(θ, λ)至球心的距离，a为椭球长半轴长，${{{\bar C}_{nm}}}$和${{{\bar D}_{nm}}}$为规格化球谐系数，$\bar P_n^m$(cosθ)为完全正规化勒让德函数。引力位变化ΔV可以用球谐系数变化Δ${{{\bar C}_{nm}}}$和Δ${{{\bar D}_{nm}}}$来表示：

$ \begin{array}{r} \Delta V=\frac{G M}{r} \sum\limits_{n=0}^{\infty}\left(\frac{a}{r}\right)^{n} \sum\limits_{m=0}^{n}\left[\Delta \bar{C}_{n m} \cos m \lambda+\right. \\ \left.\Delta \bar{D}_{n m} \sin m \lambda\right] \bar{P}_{n}^{m}(\cos \theta) \end{array} $

(2)

重力仪在地球表面观测到的重力变化δg为:

$ \delta g = \Delta g + \beta {u_r} $

(3)

式中，Δg为空间固定点的重力变化，β为布格梯度，u_r为径向位移。采用球模型近似，空间固定点重力变化Δg和引力位变化ΔV有如下关系：

$ \begin{gathered} \Delta g=-\frac{\partial \Delta V}{\partial r}=\frac{G M}{a^{2}} \sum\limits_{n=0}^{\infty}(n+1) \cdot \\ \sum\limits_{m=0}^{n}\left[\Delta \bar{C}_{n m} \cos m \lambda+\Delta \bar{D}_{n m} \sin m \lambda\right] \bar{P}_{n}^{m}(\cos \theta) \end{gathered} $

(4)

当Δg已知时，可以直接进行球谐展开：

$ \begin{aligned} \Delta g=& \frac{G M}{a^{2}} \sum\limits_{n=0}^{\infty} \sum\limits_{m=0}^{n}\left[\Delta \bar{A}_{n m} \cos m \lambda+\right.\\ &\left.\Delta \bar{B}_{n m} \sin m \lambda\right] \bar{P}_{n}^{m}(\cos \theta) \end{aligned} $

(5)

式(4)和式(5)的球谐系数满足：

$ \left[\begin{array}{c} \Delta \bar{C} \\ \Delta \bar{D} \end{array}\right]=\frac{1}{n+1}\left[\begin{array}{l} \Delta \bar{A} \\ \Delta \bar{B} \end{array}\right] $

(6)

2 重力变化计算结果

流动重力观测结果可以很好地反映地震孕育过程中产生的重力变化^[9-11]。本文重力变化数据来自中国地震局重力学科管理部，时间变化范围为2015-04~2017-04。每期重力数据处理过程是采用绝对重力结果作为起算基准，然后将绝对重力和相对重力进行联合平差，计算过程中采用绝对重力基准对相对重力仪格值进行解算以削弱格值误差。不同期数据之间作差即可得到相应时间段内的原始重力变化。

为了获取研究区域重力场动态变化图像，需对重力变化数据进行后处理，流程如下：1)剔除极少数误差较大的数据，最终使用的数据点位平均精度优于15 μGal。2)由于提供的重力变化数据是地表值，而根据上节方法计算处理时需要空间固定点重力变化值，为此进行布格梯度归算。布格梯度改正包含垂直形变效应和布格层效应。垂直梯度可取-3.086 μGal/cm，地壳密度采用平均值2.67 g/cm³，由布格改正公式得到布格改正为1.1 μGal/cm。布格梯度是垂直改正和布格改正之和，大小为-2.0 μGal/cm。研究区域每年隆升速率变化范围在1~6 mm/a^[12]，以最大隆升速率6 mm/a计算，每年因为垂直形变引起的重力变化绝对值不超过1.2 μGal，该数值小于相对重力仪观测到的重力变化，因此本文不区分重力仪观测的地表重力变化和空间固定点重力变化。3)实际数据处理时，选择整个南北地震带重力观测数据(重力观测点位分布见文献[13])。将区域内重力变化数据进行格网化，区域外重力变化值设为0，根据式(5)进行球谐分析，得到60阶球谐系数。利用球谐系数进行求和并作3°范围平滑处理，得到的重力变化见图 2。

图 2显示，九寨沟地震前2 a青藏高原东缘重力场动态变化剧烈。下面具体分析重力场变化特征：

1) 2015-04~2015-10期间，全区域绝大部分地区为正重力变化。在九寨沟、黑水、汶川、丹巴和雅安一带是重力增加峰值区域，最大正重力变化达到20 μGal。除庄浪东北地区重力变化接近0外，其余地区重力变化幅度沿着峰值区域向西北和东南两侧递减。峰值区域走向与龙门山断裂带走向基本一致，说明重力场变化与构造运动密切相关。

2) 2015-10~2016-04期间，大部分地区重力变化趋势与2015-04~2015-10相反，表现为负重力变化。巴颜喀拉块体东部和柴达木块体的贵南至泽库一带是主要的重力减小地区，周边的川滇块体东北部和华南地块西北部重力变化平缓，但在漳县和西和东北地区出现正重力变化。区域整体重力场变化特征表现为，巴颜喀拉块体东缘沿东北方向至庄浪附近，负重力变化逐渐过渡为正重力变化，形成负-正重力变化过渡带。

3) 2016-04~2016-10期间，重力场变化最明显的特点是巴颜喀拉块体负重力变化在2015-10~2016-04期间变化基础上继续减小，并且负重力变化分布向西南延伸至川滇块体东北部。阿坝、马尔康至新龙一带是负重力变化的峰值区域，最大负重力变化达到-35 μGal。在柴达木块体和华南块体内部，重力变化分布特征分别是向北侧和东侧逐渐由负重力变化向正重力变化过渡，在贵南附近、资阳和遂宁东南区域出现重力正变化特征。区域整体重力场变化特征为川滇块体东北部负重力变化向东部、东南部和北部逐渐过渡为正重力变化，形成重力变化过渡带。

4) 2016-10~2017-04期间，重力场变化在前2期基础上出现反转回调，具体表现为2个特征：一是负重力变化，在阿坝、若尔盖、岷县和漳县西北地区，重力变化由0逐渐减小，在贵南附近数值达到-20 μGal；二是正重力变化，在汶川、绵阳、平武和广元一带出现正重力变化峰值区域，最大重力变化为18 μGal，并且正重力峰值区域走向与龙门山断裂带走向一致。在重力变化峰值区域两侧，即西北方向至阿坝、若尔盖、岷县和漳县一带，东南方向至遂宁和资阳一带，正重力变化幅度逐渐减小。区域重力场变化特征是贵南附近的负重力变化沿东南方向逐渐过渡至正重力变化，形成重力变化过渡带。

3 重力变化与九寨沟地震

地震孕育过程中伴随着质量迁移和能量交换，可通过分析震前重力场动态演化特征来捕捉地震前兆信息。分析图 2可知，九寨沟地震前，震源周边重力场演化可以分为3个阶段。第1阶段为重力增大阶段，2015-04~2015-10巴颜喀拉块体东部边界重力场均表现为正变化，沿着边界向四周正变化幅度逐渐减小。第2阶段为重力减小阶段，2015-10~2016-04重力场开始出现负变化，2016-04~2016-10重力场继续减小，并且变化幅度增大。该阶段内，九寨沟地震震源均位于正负重力变化过渡带上。第3阶段为重力场恢复调整阶段，巴颜喀拉块体东部表现出与第2阶段截然不同的重力变化，即重力增加。并且在该阶段，震源仍然位于正负重力变化过渡带上。仔细对比第2阶段和第3阶段的正负过渡带可以发现，前者由负向正过渡的方向为东北向，后者为东南向，两者方向发生约90°旋转。正负重力变化过渡带旋转，表明岩石圈内部物质迁移剧烈。祝意青等^[13]指出，九寨沟地震之前震源周边重力场变化原因有2个：1)2013年四川芦山7.0级地震和甘肃岷县-漳县6.6级地震使青藏高原东部岩石圈物质重新分布；2)川西高原及邻区深部壳幔边界与上地幔物质和能量强烈交换。九寨沟地震前重力场变化体征显示，区域正负重力变化过渡带空间走向发生旋转后，紧接着发生了九寨沟地震。因此，捕捉正负重力变化过渡带空间走向旋转变化信息对于寻找地震可能发生的时间具有重要前兆意义。

本文分析重力场动态演化规律时，将经重力网平差后的数据进行球谐展开至60阶，并进行3°范围平滑滤波，这样处理后可以在图 2(b)、2(c)、和2(d)所示的3个时间段重力场变化图像中清晰显示重力正负变化过渡带。3期正负重力变化过渡带的交会区域基本集中在四川西北部，九寨沟地震震源即位于该区域内部。破坏性强震孕育会引起重力场在较大范围内出现改变，本文中所用的球谐分析方法可以较好地突出大尺度重力场变化特征，抑制局部细节特征，这种特点有利于确定地震可能发生的大概位置。已有研究表明，强震易发生在重力变化正负梯度过渡带的零值线附近^{[9, 11, 14-15]}，本文给出的结果显示，九寨沟M_S7.0地震震源位于正负过渡带内，但与零值线存在一定距离，这与本文中采用的球谐分析技术突出主要重力场变化特征、观测误差、点位分布等因素有关。实际采用流动重力观测资料进行地震前兆分析时，建议采用重力网平差后直接得到的结果和球谐分析方法得到的结果同时进行分析。

4 结语

本文采用球谐分析方法研究九寨沟地震前重力场动态变化特征，取得的主要认识有：

1) 2015-04~2017-04期间，震源周边区域重力场变化经历增大、减小和恢复调整3个阶段，后2个阶段在震源附近出现正负重力变化过渡带，且重力变化过渡带空间走向出现旋转现象。重力场演变特征较好地反映了九寨沟地震的前兆现象。

2) 本文使用的球谐分析方法，可以有效突出区域重力场演化的主要特征，便于地震前兆特征识别。但同时也会抑制局部细节变化，实际分析重力场变化数据时，建议与传统方法同时使用。