小波多尺度分析尼泊尔MW 7.8 地震前卫星重力变化

刘 潇

（中铁第四勘察设计院集团有限公司， 武汉 430063）

0 引言

2015年 4 月 25 日 14 时 11 分， 尼泊尔境内发生了MW7.8 地震， 美国地质调查局 （United States Geological Survey， USGS）测定的此次地震的震中位于尼泊尔博卡拉（28.15°N， 84.71°E）， 即位于首都加德满都西北77 km 处。 此次地震的震级较大，并引发多处滑坡、 冰崩雪崩； 造成了大量人员伤亡以及建筑物倒塌和损坏。

能够用重力手段监测地震是由于在大地震孕育和发生过程中， 震源区附近会产生物质迁移和质量再分布现象， 重力场也随之变化[1]； 2002 年3月发射成功的GRACE 重力卫星， 具有高精度、 高时间空间分辨率， 能够监测出重力变化， 其时变重力场能够提高人类认识地球内部物质迁移。 重力场的变化量揭示了地球 系统内部的物质运动、分布及变化， 它主要是由地球流体圈层的质量迁移引起的， 也可能与冰川均衡调整和地震现象造成的质量运动有关。

孙文科[2]尝试将GRACE 运用在日本大地震的研究中， 取得了不错的结果； 邹正波[3]等研究GRACE 重力变化的计算和分析， 对球谐系数选取以及高斯平滑半径进行了讨论； 王武星[4]等利用GRACE 研究了汶川地震发震前的重力异常； 已有研究成果表面GRACE 卫星观测到的重力变化包含了浅部到深部的各种质量改变成分， 可以通过位场分离的方法去识别与提取。 小波多尺度分析方法将重力异常分解到不同的尺度空间中， 尺度大小决定了重力异常所反映的地质体埋深[5-6]。

本文利用GRACE 的时变重力特征和小波多尺度分解的优点， 对尼泊尔地震前的重力变化进行研究， 并试图分析和解释重力的变化与地震的孕育及地球内部物理的流动之间的关系。

1 数据处理与方法

1.1 GRACE 的数据处理

目前提供GRACE 产品的机构有法国的GRGS（Space Geodesy Research Group）、 美国的 CSR、 德国的 GFZ （GeoForschungsZentrum）和美国的 JPL（Jet Propulsion Laboratory）， 其 中 除 了 GRGS 的RL02 产品为10 日解， 其他产品均为月解。

根据Wahr[6]的理论， 定义重力扰动计算公式：

式（1）中： R 地球的半径， ΔC1m和 ΔS1m都是球谐系数的变化值， θ 和 φ 余纬（距离 z 轴的角度， 称作余纬度或顶角， 角度从 0 到 180°）和经度正则化的勒让德函数

1.2 二维小波分解

小波多尺度分析方法作为一种位场分离途径，可以将重力异常分解到不同的尺度空间中， 尺度大小决定了重力异常所反映的地质体埋深。 在实际的GRACE 重力数据处理当中， 常用的是二维小波变化， 所以有必要将小波变化从一维推广到二维， 这里直接给出二维离散型小波变换公式。 对于函数（信号）f（t）∈L2（R2）， 是一个二维信号， x，y 分别表示其纵横坐标， ψ（x，y）表示二维基本小波， 满足条件[5-7]：

二维小波变换具有旋转的能力， 不但具有放大的功能， 还会产生极化的现象。 在实际应用中，二维小波变换通过不同阶数的分解提取出来各个方向的分量。

2 尼泊尔地震GRACE 重力场的变化

选取的空间尺度为 15°～45°N， 70°～100°E，分布于震中四周， 南北、 东西均大于3000 km， 以更大的空间尺度显示该震前的重力变化信息， 时间尺度为震前 5 年， 即 2010-01 至 2015-04， 以反映震前中短期的重力变化信息[8-10]。

本文使用 UTCSR 提供的 RL05 （Level-2 Release-05）月重力场产品GSM（GRACE Satellite only Model）， 最大阶数 为 60 阶； 数据 处理 中采 用Swenson[11]的一阶项结果代替； 采用高斯平滑滤波方法， 选取平滑半径为400 km； GRACE 轨道存在系统误差， 使用 SLR 计算的 C20 替换12]； 利用Duan[13]方法做去相关性处理。 获得年度月重力累计变化（图 1）。

图1 尼泊尔及邻区年度重力变化Fig. 1 Annual gravity variation in Nepal and its adjacent areas

从图1 年度累计重力变化连续图发现：

此次大地震前， 尼泊尔南北出现了比较明显的重力异常变化， 分别是北部的正异常变化， 峰值幅度达到25×10-8m/s2， 以及南部地区的负重力异常变化， 峰值幅度达到-30×10-8m/s2。

为了定量的比较震前的重力值变化[14]， 利用减去2006—2014 年的GRACE 数据的平均系数值，计算了尼泊尔地震中心的卫星重力异常时间序列（图2）。 GRACE 反演的结果中有明显的地表周期性水文影响， 为消除这些变化的影响， 通过最小二乘法扣除年、 半年、 季节尺度重力场变化， 以突显地震相关重力信息。

从图2 中可以看出， 在2012 年5 月之后重力异常值处于累积的（正值）， 出现这样情况的原因是， 大地震发生前在震源区可能有物质迁移， 使得该地区出现了明显的重力积累效应；

图2 尼泊尔地震中心（28.15°N，84.71°E）的卫星重力变化时间序列Fig. 2 Satellite gravity center of the Nepal earthquake center（28.15 ° N，84.71 ° E）

3 尼泊尔地震GRACE 重力场小波多尺度分解及结果分析

利用GRACE 卫星观测到的重力场主要是由水储量的变化、 构造形变和地下物质流动等因素引起的， 水储量变化、 冰盖融化等引起的地表密度变化相对较大且变化较快。

由于水储量、 构造形变以及深部物质流动在空间尺度上具有差异， 由浅入深可能表现为不同波长成分。

一般情况下， 小波分解最大尺度选取的越大，那么被滤除的细节部分就越多， 分离出的成分越多； 但是在进行信号重构的时候系数置零的部分就越多， 重构后的信号就会丢失很多有用信号，造成信号的失真。 一般， 对于信噪比较大的信号，观测信号比有用信号的量级要大， 这时候较小的分解阶数就能很好的解决问题， 相反的情况就需要分解的阶数较高。对于重力数据理论模型， 不同尺度分解出来的重力数据， 对应不同深度、 不同密度的不均匀分布重力异常多尺度分析的最理想结果， 是通过功率谱分析， 计算各阶细节对应的场源深度。 同时与小波多尺度分解出的细节部分对比， 从地球物理的角度， 确定分解阶次的合理性。 刁博在2007 的文章中发现5 阶细节反映的场源深度为17～34 km，而6 阶细节只是对5 阶细节的伸缩， 是布格重力异常在小波分析纯数学意义上所产生的假象[15]。 因此，利用小波多尺度分析方法对2011-2014 年GRACE卫星观测到的重力场进行5 阶分解， 结果如图3。小波分解的2～4 阶细节反应的是 7～20 km 左右的浅表扰动或者上地壳密度变化， 5 阶细 节功率谱场源似深度为17～34 km， 主要反映的是中、 下地壳密度变化； 5 阶近似主要反映的是34～100 km 范围下地壳及深部地幔物质运移引起的密度变化， 基本上剔除了中、 浅部质量迁移引起的重力变化[16]。从5 阶近似可以看到从2010 年到 2014 年， 尼泊尔震区呈现大范围的重力上升变化， 幅值为10μGal， 这可能是深部高原物质移动产生的结果。

图3 2011—2014 年GRACE 重力变化5 阶小波多尺度分解细节及近似部分Fig.3 Details and approximate part of GRACE gravity change 5th order wavelet multi-scale decomposition from 2011 to 2014

4 结语

利用GRACE 重力卫星RL05 月重力场数据，获取尼泊尔MW7.8 地震前后震源区周缘2010-2015年每月差分重力变化， 以及震中点位重力时间序列变化。 并利用小波多尺度分析方法对卫星重力场进行分解， 得到了反映不同深度的重力场细节和近似部分。 研究结果表明： 尼泊尔地震前5 年（2010-2014 年）内在震源区周边出现了比较明显的卫星重力异常正负交替和迁移现象， 2014 至2015年间， 震区周边形成了明显正负异常区， 正重力异常区重力增加现象明显。 震中的重力时间序列分布指出从2013 年开始， 重力变化处于稳定状态， 并都处在一个较高位； 同时通过功率谱分析，计算多尺度分解各阶细节对应的场源深度， 2 阶、3 阶和4 阶小波分析后的重力异常细节， 反映了浅部地质体的位置和异常情形； 5 阶小波细节更加侧重于刻画地区构造形变和深部物质流动引起的重力变化。 反映出大地震前震源区周边地下物质运动、 质量迁移和能量积累等问题。