引言

地震通过能量释放及地球内部与表面的永久形变和渐变形变，能够引起重力场数千千米的空间变化以及数秒至几十年的时间变化(Han et al，2013)，美德合作的GRACE(gravity recovery and climate experiment)重力卫星可以观测到巨大地震的重力场时间和空间变化. 该卫星通过十余年高低跟踪与低低跟踪技术组合观测，为科学界提供了长期近连续的高精度时变地球重力场信息，为研究大型板块运动及其对地震孕育发生的影响提供了重要的科学观测依据，为研究地震过程中重力场的时间和空间演化奠定了基础 .

目前相关研究已表明，GRACE卫星能够探测到大地震前重力场的时空变化，已观测到2004年苏门达腊M_W9.3和2011年日本M_W9.0等地震的同震重力变化(Han et al，2011；Zhou et al，2012；邹正波等，2013；张克亮等，2014). 2011年日本M_W9.0地震发生在日本海沟底部，与太平洋板块、北美板块和欧亚板块构造活动密切相关；而且该震震级大，处于中国大陆边缘构造带上，对中国大陆地震活动的影响深远.

经验正交函数方法早期主要应用于气象气候和海洋学，后逐渐扩展至地磁、重力等领域. De Viron等(2008)首次将该方法应用于卫星重力地震研究领域，给出了GRACE卫星在现有精度下对地震的探测能力；Schmeer等(2012)利用该方法从GRACE月重力场中定量地分离出海洋、大气及陆地水对地球质量变化的贡献；周江存等(2013)利用该方法对日本M_W9.0地震进行了研究，但是由于当时仅有震后两个月的资料(2008年1月—2011年5月)，未能提取出同震跳变以及与模型相对应的沿断层两侧的正负重力变化；张克亮等(2014)基于2003年1月—2013年6月的GRACE数据采用最小二乘法和经验正交函数方法，提取了日本M_W9.0地震震区及其周边的同震变化特征，并对其进行了重点分析，但未将陆地水储量变化这一对非构造重力变化影响较大的因素考虑在内.

本文将在前人研究的基础上，选取更长时段(2002年4月—2015年4月)的卫星资料，在充分考虑陆地水储量变化影响后，对日本M_W9.0地震震中及其周边地区的时变重力场进行计算分析，给出区域长期累计重力变化和差分重力变化结果，并利用经验正交函数对区域重力场进行时空分离，以综合分析该地震震前孕震、同震释放和震后调整对该区域的影响以及该地震对重力场变化的贡献，从而探查日本地震重力场演变的时间和空间特征.

1.   数据

1.1   GRACE卫星时变重力场模型数据

GRACE卫星自2002年3月17日发射以来，不断对外发布产品，其中应用最广泛的是月尺度球谐系数组成的重力场模型，即月重力场模型. 本文采用美国德克萨斯大学空间研究中心发布的RL05月重力场模型，该模型在处理过程中已扣除固体潮、极潮、大气、海洋等因素的影响. 选取2002年4月—2015年4月GRACE重力场模型对重力时空信号进行分析，其中2003年6月、2011年1月和6月、2012年5月和10月、2013年3月和8—9月因缺乏关键有效载荷数据而无模型，另外地震发生时段2011年3月的数据也一并剔除. 根据产品说明将模型中的C₂₀替换为美国德克萨斯大学空间研究中心发布的由卫星激光测距数据计算所得到的C₂₀.

1.2   全球陆面数据同化系统陆地水数据

全球陆面数据同化系统(global l and data assimilation system，简写为GLDAS)是由美国宇航局利用近实时地空观测系统数据约束陆地表面的状态，采用数据同化技术，提取陆地表面状态及流量变量(如等效雪水、土壤水湿度、总植冠水、降雨、降雪、蒸发、径流等)的全球水文模型之一，其对外发布的陆地水模型最大时间分辨率为3小时，空间分辨率为0.25°×0.25°(Rodell et al，2004). GRACE卫星信号中的陆地水储量变化对总质量变化的影响显著，利用GLDAS的1°×1°月尺度产品可有效提取陆地水储量变化(包括土壤水湿度、雪水储量和植冠表面水储量)的影响，有助于分析GRACE结果中构造因素对重力场的影响.

2.   重力时变信号提取

2.1   重力变化参数选取

重力变化可用重力扰动δg和重力异常Δg来表示. 前者为任意点P的实测重力值g_P与该点正常重力值γ_P之差，即δg=g_P-γ_P；后者为P点的实测重力g_P与该点对应的正常椭球面上Q点的正常重力值γ_Q之差，即Δg=g_P-γ_Q. 两者之间存在如下关系(Hofmann-Wellenhof，Moritz，2005)：

其球谐表达式分别为

式中，G为万有引力常数，M为地球质量，a为地球平均半径，n和m分别为地球重力场模型的阶和次，r为该点到地球质心的距离，N_PQ为P点与Q点之间的距离，λ为经度，θ为余纬，P_nm为规格化勒让德函数，ΔC_nm和ΔS_nm为真实重力场与正常重力场的完全规格化球谐系数之差.

考虑到重力扰动表征的是测点自身的重力观测值与该点的正常重力之差，易于与其它观测资料如断层模型模拟的重力变化比较. 由式(1)—(3)可以看出，重力扰动与重力异常相比存在一个(n+1)/(n-1)的比例因子，且该值越大，越容易从卫星微小信号中提取异常变化，故本文选取重力扰动这一参数来表征重力时变特征.

2.2   GRACE重力时变信号提取

为有效减小重力变化中明显的非地球物理信号性南北向异常条带的影响，避免去相关处理所造成的误差，本文采用扇形滤波技术(Zhang et al，2009)来处理异常条带，并选取统一基准作为背景场，计算GRACE月重力扰动δg_Gm为

式中，W为高斯滤波，r_fan为扇形滤波的平滑半径.

本文选取GRACE平均月重力场模型作为背景场，利用式(4)计算其重力变化δg_Gb，由此得到便于解释的GRACE重力变化dg_GRACE为

2.3   GLDAS陆地水重力时变信号提取

陆地水储量变化是GRACE信号中最重要的非构造因素之一，可利用模型定量地研究其与地震相关的重力变化. 本文利用GLDAS的1°×1°产品，综合土壤水湿度、雪水储量和植冠表面水储量等提取陆地水储量. 由于任意三维函数均可用球函数表示(郭俊义，1994)，因此由陆地水储量变化可得到与之对应的球谐系数为(Wahr et al，1998)

其中，ρ_w为水密度，k_n为复合勒夫数，Δσ(θ，λ)为网格点的陆地水储量，C^GLDAS_nm和S_nm^GLDAS均为基于GLDAS模型得到的球谐系数.

将解算得到的球谐系数代入式(4)—(5)即可得到相应的由陆地水储量变化引起的重力变化. 与1°×1°水文模型对应的重力场模型球谐系数最大阶次为180，而GRACE重力场模型对应的最大阶次为60，其分辨率为3°×3°. 为便于与GRACE结果比较，本文选取60阶次的球谐系数进行计算和后续讨论.

以2011年3月的GLDAS陆地水储量变化向重力变化转化为例，给出利用GLDAS提取重力变化信号的过程. 首先，对原始信号(图 1a)利用式(6)进行球谐处理，结合式(2)得到等空间分辨率的重力信号(图 1b); 然后，对球谐系数60阶截断后，利用式(2)取得与GRACE结果空间分辨率相当的重力信号(图 1c)；最后，在此基础上利用300 km扇形滤波[式(4)]处理后得到可以与GRACE结果比较的重力信号(图 1d). 比较图 1a-c与图 1d可知，GLDAS数据处理虽然使原始信号的振幅减小，空间分辨率降低；但仍能有效地捕捉到高频信号的影响，为更好地分析GRACE重力变化的原因提供了有力的支持.

图  1  由GLDAS陆地水储量变化到重力变化(dg_GLDAS)的转换过程处理

(a) GLDAS 1°×1°陆地水储量; (b) 截断至180阶后的重力变化; (c) 截断至60阶后的重力变化; (d) 截断至60阶，300 km扇形滤波后的重力变化

Figure  1.  Process of conversion from variation of terrestrial water storage to gravity change (dg_GLDAS) based on the GLDAS

(a) Terrestrial water storage with 1°×1° resolution; (b) Gravity changes with truncation at 180; (c) Gravity changes with truncation at 60; (d) Gravity changes with truncation at 60 and 300 km fan filter

下载: 全尺寸图片 幻灯片

2.4   扣除陆地水储量影响的重力时变信号提取

利用2.3节所述方法提取与GRACE观测时间匹配的陆地水储量引起的重力变化dg_GLDAS，结合GRACE卫星重力变化dg_GRACE，可得到研究区域内非陆地水储量因素引起的重力变化dg=dg_GRACE-dg_GLDAS. 下面以2011年3月地震发生时的数据为例进行计算，以分析采用GLDAS数据对GRACE重力变化结果的影响.

图 2给出了2011年3月由GRACE观测到的重力变化，由陆地水储量引起的重力变化，以及扣除陆地水储量变化这一非构造性因素影响得到的重力变化. 可以看出：GRACE卫星观测到震区两侧呈明显的正负异常，日本海附近重力减少，其变化幅值最大为-5×10^-8 m/s²；震中东南局部重力增加，增幅小于4×10^-8 m/s²(图 2a)；而陆地水信号并未出现相关响应，仅在中国东北及俄罗斯边界区域出现了重力增加(图 2b)；扣除陆地水影响后，在日本岛两侧出现了以“负—正”为主要特征的重力变化结果(图 2c)，更有利于研究地震等构造因素对区域重力场的影响.

图  2  2011年3月区域重力场变化

(a) GRACE观测到的重力变化; (b) GLDAS陆地水储量引起的重力变化; (c) GRACE扣除GLDAS影响后的重力变化

Figure  2.  The regional gravity changes in March of 2011

(a) From GRACE; (b) From GLDAS; (c) From GRACE minus GLDAS

下载: 全尺寸图片 幻灯片

2.5   经验正交函数方法提取区域重力时空变化

经验正交函数(empirical orthogonal function，简写为EOF)方法，也称为主分量分析方法，是利用数学手段提取物理场主要特征的一种重要方法. 该方法主要利用数据集的方差贡献来识别和强化信号的相似性，通过将随时间变化的地球物理场分解为不随时间变化的空间部分和不随空间变化的时间部分，从而有效实现时空信号的分离.

经验正交函数方法是通过对观测矩阵Y进行分解实现的. Y为扣除平均背景场后的重力场时空变化，大小为n×p，其中n为观测值的个数(即格网点数)，p为观测值的时间长度. 对时空相关观测矩阵Y进行经验正交函数分解，使之满足

式中：E是由Y的n×n信号协方差矩阵C=YY^T/p的特征向量组成的矩阵；C=EΛE′，其中Λ为包含C矩阵n个特征值的对角矩阵，且特征值按从大到小的顺序排列；D=E^TY与表征时间的主分量(时间系数)对应.

这样，通过经验正交函数分解可将观测场Y分离为不随时间变化的空间特征矩阵E和不随空间变化的时间特征矩阵D. 因此，经验正交函数分解可以理解为对矩阵Y按照新的正交基(E)进行重新投影，使之表征为空间函数与时间函数的线性组合.

3.   累计重力变化和差分重力变化结果

通过GRACE卫星长期对地重复观测，高效获取了全球地球重力场的时空变化，为本文研究区域重力场模型的重力变化提供了可靠的资料. 然而GRACE观测到的重力变化是包括地震在内的多种地球物理过程的综合效应，因此分析日本M_W9.0地震震中及其周边区域重力场时空变化特征时，需认真分析各因素的影响. Matsuo和Heki(2011)认为研究区域陆地水的季节影响和年季水变化(主要是冬季雪)影响是除地震外的最重要因素，因此本文利用GLDAS数据扣除上述两种因素的影响，便于后续分析区域重力场与地震的关系.

图 3和图 4分别给出了日本震区周边(30°N—50°N，125°E—160°E)2003—2014年扣除背景性周期性信号影响后的累积重力变化和差分重力变化. 图中年度重力变化均基于当年3月至第二年2月的数据，其中2011年重力变化扣除了地震发生月份(3月)的数据. 可以看出，日本地震震区及其周边区域重力场变化的时空演变具有以下特征：

图  3  2003—2014年日本地震震区及其周边区域地震前(a-h)、后(i-l)的累积重力变化

图中红色方框为断层边界，黑色圆点为日本M_W9.0地震震中

Figure  3.  Accumulated gravity changes of the seismic zone and its surrounding areas from 2003 to 2014 before (a-h) and after (i-l) the M_W9.0 Janpan earthuqake. Each panel presents the evolution of annual regional gravity changes. The red rectangle denotes the border of fault, and the dot denotes the epicenter of the M_W9.0 Japan earthquake

下载: 全尺寸图片 幻灯片

图  4  2003—2014年日本地震震区及其周边区域的差分重力变化

(a)-(g) M_W9.0地震震前差分重力变化；(h) 同震重力变化; (i)-(k) 震后差分重力变化

Figure  4.  Differential gravity change of the seismic zone of M_W9.0 Japan earthquake and its surrounding area from 2003 to 2014

(a)-(g) Pre-seismic differential gravity change; (h) Coseismic gravity change; (i)-(k) Post-seismic differential gravity change

下载: 全尺寸图片 幻灯片

1)研究区域范围内的重力基本在-1×10^-8—1×10^-8 m/s²内波动，尤其是2003—2008年重力变化相对较小，自2009年起重力变化出现高于1×10^-8 m/s²的异常波动，并持续增加. 从累计重力变化(图 3)可以看出，2011年后重力变化基本被地震同震信号类似的特征控制，而差分重力变化(图 4)也在2011年后呈现截然不同的特征. 引起重力变化的原因一般与地球物理事件相关，而GRACE重力变化主要源于陆地水储量和地下水的变化以及冰后回弹、板块运动、地震等；本文中陆地水储量变化已经从结果中扣除，日本及其周边区域2008—2011年并未发生除M_W9.0地震以外的重大事件(Matsuo，Heki，2011)，由此可推测区域重力场变化与板块运动及地震具有相关性.

2)日本地震震源区及其周缘震前重力呈现异常变化，在累积重力变化结果中尤其明显. 从2009年起中国东北、朝鲜半岛以东及日本海以西区域出现大范围的重力增加，2010年重力增幅有所减弱，并向日本岛迁移，直至2011年发生日本M_W9.0地震.

3)日本M_W9.0地震的同震效应对区域重力场影响显著. 由2003—2014年累计变化结果(图 3)可以看出，该地震发生前后重力变化的空间分布特征存在显著差异，相对于震前，震后累积重力(图 3i-l)的空间分布与同震重力变化(图 4h)特征基本一致，说明区域重力场变化特征以M_W9.0同震变化为主. 从差分重力变化(图 4i-k)可以看出，震后震区重力开始增加，但增加幅度及位置有所差异.

4)震后重力场出现明显调整. 由2012年相对于2011年的震后重力差分结果(图 4i)及2012年累积重力变化(图 3j)可以看出，该地震所造成的断层以东区域重力负变化已经开始有所补偿，且持续至2014年(图 3l，4k)，而在日本岛以西震中所在区域的重力仍在增加.

上述重力场变化特征均为研究日本M_W9.0地震震后区域重力场调整及其动态演变过程提供了宝贵的观测证据，可用于后续黏滞系数等地球参数的反演.

4.   经验正交函数方法提取的重力时间变化结果

针对2002年4月—2015年4月GRACE卫星重力时变结果，选择不同的时间长度对卫星重力资料进行经验正交函数分解，形成5个方案. 通过分析各方案的时变重力场空间分布特征，给出各主分量所表征的物理事件对区域重力场的贡献，以揭示M_W9.0地震对区域重力场的影响. 虽然经验正交函数方法在提取数据集主要信号方面具有优势，但其各主分量未必具有明确的物理含义，且对结果的物理解释需要相应的先验条件. 本文研究区域的先验信息物理解释依据主要来源于：① 位错模型模拟的理论同震重力场变化(Zhou et al，2012；邹正波等，2013；张克亮等，2014)；② GRACE观测到的区域时变重力场(图 3，4). 此外，重力变化在时间分布上是否存在同震跳变以及在空间分布上是否存在与同震变化相似的特征，也是为经验正交函数结果提供合理物理解释的重要依据.

表 1列出了不同方案的前7个主分量方差的贡献. 可以看出：第一主分量的方差贡献百分比最大，基本为27%—36%；第二主分量有所减少，为16%—18%；第三主分量则为10%—12%，因此前3个主分量的累计方差贡献达到55%—72%，前两个主分量的累积方差贡献为45%—52%. 比较5套方案可知，5种时间方案的方差贡献差异并不大，相对而言最长时段的方案1信号相对不集中. 虽然各方案选取的时间长度不同，但经验正交函数分解后的时空变化特征显示各主分量表征的物理含义基本相同，故本文选取时间相对较长的方案2即地震前后4年的数据进行分析.

表  1  不同方案的前7个主分量方差贡献

Table  1.  Variance contribution of the first seven principle components for five schemes

方案	时段	各主分量方差贡献
		第一主分量	第二主分量	第三主分量	第四主分量	第五主分量	第六主分量	第七主分量
1	2002-04—2015-04	27%	18%	10%	9%	6%	5%	3%
2	2007-03—2015-03	34%	17%	12%	6%	5%	4%	3%
3	2008-03—2014-03	36%	17%	9%	6%	6%	4%	3%
4	2009-03—2013-03	34%	17%	11%	7%	6%	5%	3%
5	2010-03—2012-03	36%	16%	12%	7%	6%	5%	3%

下载: 导出CSV 

| 显示表格

4.1   重力场时间变化特征

通过分析各主分量(即时间系数)能够获取重力场随时间变化的规律. 图 5给出了方案2(2007年3月—2015年3月)的重力场时间变化，可以看出：第一时间系数振幅最大，为-80—100，存在显著的同震跳变，地震前后变化率差异显著；第二时间系数存在明显的周期性变化特征，地震前后一年出现剧烈变化，震前增加，震后减小，发震时刻变化量达到最大；第三时间系数的幅值则较小，为-40—60，仅2015年后期数据出现扰动. 总体看来，日本地震前(2008年下半年开始)时间系数出现明显的趋势性增加，震前一年(2010年)该值减小，这与唐山地震前的重力变化情况相一致(李瑞浩，1997).

图  5  采用方案2所得到的研究区域前3个时间系数的时间变化

Figure  5.  The first three time coefficient variation of EOF for the scheme 2 in the studied area

下载: 全尺寸图片 幻灯片

4.2   重力场空间变化特征

图 6给出了利用经验正交函数方法所得到的前3个主分量的空间特征场，可以看出: 第一空间特征场(图 6a)区域几乎以发震断层为界分成东西两部分，日本海出现范围较小的重力场负变化，正变化区则相对较大；第二空间特征场(图 6b)显示研究区域重力场基本为负，说明重力场变化特征一致，即重力均增加或减小，结合图 3和图 4，可知其变化与震前(日本岛以东)、震后(集中在本州岛及周边)的重力异常变化密不可分；第三空间特征场(图 6c)的空间分布呈正—负—正—负的特点，震区最东侧正变化依然显著，说明该特征场仍然可能包含同震及震前信号，且出现了一定的南北向条带异常，证实了经验正交方法也可以扣除GRACE的噪声，即将南北向条带异常留在相对贡献较小的第三分量.

图  6  研究区域经验正交函数的前3个空间特征向量场

(a) 第一空间；(b) 第二空间; (c) 第三空间

Figure  6.  The first three spatial eigenvector fields of EOF in the studied area

(a) The first; (b) The second; (c) The third

下载: 全尺寸图片 幻灯片

5.   讨论与结论

基于2002年4月—2015年4月GRACE卫星的长期重力变化资料，扣除全球陆面数据同化系统陆地水储量变化的影响，得到了日本M_W9.0地震震区及其周边的累积重力变化和差分重力变化. 利用经验正交函数方法从海量时空数据中有效地分离出时间空间信号的物理意义不明确的特征，以位错模型模拟所得的同震信号及区域累积重力变化和差分重力变化作为先验信息，获取了区域重力场的空间和时间演变特征.

基于经验正交函数方法的区域重力场结果显示，重力时空变化的第一主分量与日本M_W9.0地震的同震重力变化相关性强，即时间系数第一主分量变化幅度大且出现同震跳变，空间场第一主分量出现与理论模拟结果类似的近垂直断层方向的由负到正的重力变化. 这说明日本地震对该区域重力变化信号的贡献显著，其在GRACE卫星2002—2015年长达13年的观测结果中的贡献为27%，在地震前后1—4年的重力场信号中其贡献更大，可达34%—36%，充分体现了地震对区域重力场所造成的显著影响.

本研究进一步证实了日本地震在造成区域重力场变化的诸多因素(气象、水文及构造等)中起主导作用. 但是，本文利用经验正交函数方法分析所得结果的前3个主分量仅能解释约50%的区域重力场信号，下一步将考虑利用旋转经验正交函数等方法使地震相关重力信号更加集中，以最少主分量提取最佳地震重力场时空变化信息.