引言

重力场作为最基本的地球物理场，其静态异常特征和动态时间变化均可为深地壳结构及其性质的研究提供重要的科学依据. 21世纪以来，随着重力观测技术的迅速发展，高精度绝对重力仪器的商业化和便携化，使得高精度的微重力观测成为可能. “微重力”或“微重力测量”一词最早出现于20世纪60年代初，是指重力测量的精度和测量探查的对象引起的重力效应以微伽的数量级来度量.以20世纪80年代美国的Lacoste型重力仪为代表的相对重力仪的出现，使研究重力场的时间变化成为可能(王谦身，1995).

区域性的重力场微变化与地表形变、近地表物质运移和深部介质变形等因素密切相关，可表示为测量时间和测点位置的函数，具有明确的地球物理意义，因此可用于研究地球内部的物质变化过程(Li，Fu，1983；Crossley et al，2013).重力场的变化与地壳深处的地震孕育发生过程是紧密相联的，以研究地壳深处的地震孕育发生过程为目的的地震重力测量，通常使用多台相对重力仪按设定的路线和周期对固定测点进行重复观测，以获取可能与地壳变形和震源介质物性变化相关的重力信号(许厚泽，2003；孙和平，2004).重力测量具有速度快、耗资少、覆盖面大、灵活机动等特点，能有效地监测地壳运动和地球内部的构造活动(贾民育，1996).通过对流动重力测量数据进行绘图得到的重力场动态变化图像是地震监测研究的基本信息来源之一，利用该图像可研究监测区域重力场动态演化特征，为地震孕育、发展、发生过程的研究提供科学依据(陈运泰等，1980；李辉等, 2009, 2010; 陈石等，2015).此外，通过重力测量还可以获得构造活动区的重力场随时间的非潮汐变化，研究地震孕育、发生和调整过程中重力场的时、空、强变化，而且还可以利用地球重力场的空间分布特征和非潮汐变化研究地球内部构造和地球动力学过程(国家地震局，1979).通常以陆地运输方式，定期开展固定点位重力测量的监测工作可以称为地震重力或流动重力，是一种典型的时变微重力信号获取途径.

以往分析重力场时变信号时多采用等值线方法，但由于地表实际重力观测受限于环境、交通、仪器、成本等因素，重力测点空间分布不均匀，时间采样高度离散化.另外，由于重力场随时间变化的量级很小，与测量误差之间的区别通常不显著.常规方法中，对高度离散的测点进行插值后，再进行等值线绘制，往往会引入一定的虚假异常信号.针对这些问题，本文拟提出一种新的可视化方法以研究时变重力场的变化特征，并定义相关指标量辅以定量分析.

基于本文提出的方法，利用首都圈重力测网近几年的多期平差结果，尝试性地对首都圈地区的重力场动态变化进行分析.由于近年来该区地震活动性较弱(张培震等，2013)，自1998年张北M_S6.2地震后，M6.0以上地震多年平静，仅有M4.0左右的较小地震发生，本文选取2014年9月6日河北省张家口市涿鹿县M_S4.3地震(40.3°N，115.4°E)前后的重力测量结果，并辅以相关分析，以期精确地刻画出发生变化的重力测段的位置、大小及其显著性程度，为研究重力场随时间的微变化信号与地震活动性之间的关系提供有力的技术支持.

1.   流动重力测量与数据

流动重力测量数据具有空间高度离散、时间变化量级小等特点.在分析流动重力观测所获得的重力场变化信号时，一般需要对多期重力测量结果进行平差计算先得到重力点值，再对同一重力测点不同测量时间的结果进行差分，这样获得的结果常被用于研究不同时空尺度的重力场变化特征.在重力点值分析过程中，由于受到测网分形特征(贾民育，1996；徐如刚等，2007；李辉等，2010)、平差计算方法(董运洪等，2006)、仪器特性以及野外观测所处环境条件等因素的影响，重力场动态变化图像并不能完全如实地反映区域重力异常特征.其图示方法也常采用等值线形式，这种基于重力点值的等值线表示方法存在明显的缺陷：在重力等值线的绘制过程中，需将不规则分布的重力点值进行插值后才能完成，其等值线形态不仅受限于插值算法，而且还与测点的空间分布有关，因此可能会引入虚假的重力变化, 而重力点值随时间变化的量级很小，通常变化范围仅为0—100×10^-8 m/s²，虚假信号严重地干扰了微重力变化分析.另外，等值线方法也不能在可视化重力场变化的同时表示观测误差，重力点值变化与测量误差的区别并不显著，测量误差本身可能会引起重力变化的不确定性.

因此，为了解决上述影响因素对重力场变化分析所造成的影响，本文重新定义了一种基于流动重力段差变化的可视化方法.流动重力测量主要使用相对重力仪在不同测点测量获取测点间的重力段差，因此从数据上来讲，相比平差后的重力点值，段差值是直接测量结果，更具客观性，能直接反映所测量的物理场信息，同时辅助传统的等值线表示方法，综合判断时变微重力信号的可靠性程度.

2.   可视化方法及指标量定义

2.1   流动重力段差变化表示法

流动重力测量是在固定测点间进行往返测量，相邻两个测点之间的重力段差作为测量的最基本单元，具有一定的独立性.重力段差的测量误差除了与扣除仪器漂移、固体潮等因素后的往返测量闭合差有关外，还与两个测点之间的段差值大小相关.而对一个测段进行多次测量获取的重力段差变化相比段差值本身仅为1/10³—1/10⁴.针对这些特点，本文研究重新定义基于流动重力段差变化的可视化方法.

如图 1所示，图 1a中A点和B点表示一个测段的相邻两个测点，O点表示测段AB连线的中心点.以O点为椭圆中心，椭圆的两个轴所在直线分别用p轴和v轴表示，p轴平行于测段AB，v轴垂直于测段AB.箭头所在直线与p轴重合，箭头的尾部起始于椭圆中心O点，箭头指向段差变化增大的方向.对重力段差值(单位为10^-5 m/s²)、重力段差的时间变化值(单位为10^-8 m/s²)和测量误差(单位为10^-8 m/s²)进行归一化处理，使得段差、段差变化和测量误差这些物理量以相同比例变化.归一化处理时，令

图  1  流动重力变化段差表示

(a)图示所含基本要素; (b)箭头位于椭圆内部; (c)椭圆沿p轴伸长且箭头位于其外部; (d)椭圆沿v轴伸长且箭头伸出其外部; (e)段差变化指向一点; (f)段差连续变化

Figure  1.  Illustration of the campaigned gravity difference between adjacent stations

(a) Elements contained in the diagram; (b) The arrow is inside the ellipse; (c) The ellipse extends along the p axis and the arrow is out of the ellipse; (d) The ellipse extends along the v axis and the arrow protruding ellipse; (e) Segment difference points to a point; (f) Continuous gravity differences between adjacent stations

下载: 全尺寸图片 幻灯片

式中，y_i为相应段差归一化处理后的数值, x_i为段差的测量值(同样方法处理段差变化和测量误差).图 1a中，椭圆v轴的轴长表示相邻两测点间的重力段差大小，椭圆p轴的半轴长度表示测量误差大小，箭头的长度表示相邻两期测段AB的重力段差变化大小.

如图 1c和1d所示，当箭头的长度超出椭圆的边界，说明段差变化大于测量误差，这意味着发生了显著性的重力变化，并且箭头长度大于椭圆p轴半轴长度越多，说明测量结果越可靠；反之，如图 1b所示，当箭头在椭圆内部，说明测量误差大于段差变化，意味着该测段的测量结果在分析重力场变化时所反映的信号误差较大，结果不可靠.如图 1d所示，椭圆在v轴方向越长，在p轴方向越扁，说明段差值虽然很大，但是观测误差反而较小，这样的数据说明了观测时精度较高，数据可靠.因此椭圆在p轴方向的扁平程度可以用来衡量测网某一期测量结果的整体误差大小.

另外，在实际的重力段差变化图像中，常出现以下两种类型的异常：如图 1e所示，周围所有测点的段差变化箭头方向都指向位于中间的测点E，这类“孤立型”变化可能说明该测点周围的环境发生了变化，测量结果存在不确定性；如图 1f所示，几个测点均连续地指向测点E，说明该处的重力场变化具有趋势性和连续性，反映了重力异常的区域性变化趋势.这两种典型特征有助于判别区域重力异常特征，且对重力异常可靠性的判断十分重要.

2.2   指标量定义

以往对重力场动态变化特征分析时，通常采用定性分析方法来描述重力场变化图像.本文结合段差变化分析方法，结合测量的重力段差变化及其误差来定义两种定量参数化指标，即G值和C值，用于分析区域重力场变化的显著性程度，可作为分析重力场变化与强震之间关系的定量依据.

G值定义为

式中，N为有效测段个数，dg_i为第i个测段的重力变化大小，e_i为第i个测段的测量误差.指数G反映的是测网整体产生重力异常的程度，可以用来判断测网整体是否发生了可观测的区域重力场变化.当0 < G < 1时，认为测网没有发生能够观测到的区域重力场变化；当G>1时，则认为测网整体存在重力场变化信号.

C值定义为

式中N_i为大于i倍误差变化的测段数. C_i反映了大于i倍误差变化的测段所占的比例.

2.3   多期流动重力数据的可视化

针对多期重力场变化的时空分布特征分析多个测段的时空特征时，借鉴了信息科学中常使用的热力图(heatmap).热力图是常用的一种数据可视化方法，可直观地反映研究对象某一特征的空间分布态势，即用热力图中某一点的位置来表示不同空间位置的物理量分布情况，用每个点的颜色变化来反映该点某个属性的强度变化大小(张俊涛等，2012).通过热力图进行可视化的最大优势是可以简单地聚合大量数据，并可以使用一种渐变色标表示，直观地展现空间数据的疏密程度或频率高低.

随着测量期次的增多，流动重力测量数据量也随之增大，密集程度提高，为了更加直观地反映多期重力数据随空间位置的分布情况，通过热力图以一种直观的形式来呈现空间位置分布不同的测段的重力变化情况；进而由此得到测网整体随时间的变化，明确哪些区域的测点重力随时间变化较大、较明显，以及哪些区域的测点随时间变化的幅值较小、较稳定.

图 2给出了多期流动重力测段数据的热力图可视化结果示意图, 可以直观地展示不同测量时间的流动重力测量结果，图中对应色标的颜色表示对应测段的重力测量结果.从热力图示意图中，可以直观地展示区域内时变重力场的变化特征，包括重力正异常、负异常出现的时间和空间信息，同时有助于分析重力场的时变变化特征.此外，在实际测量中因各种因素导致的数据缺失，也可在热力图中展示.

图  2  流动重力测量结果的热力图

Figure  2.  Heatmap for the campaigned gravity measurement results

下载: 全尺寸图片 幻灯片

3.   实例分析

本文结合以往常用的表现重力场时空变化特征的等值线方法(殷秀华等，1982；贾民育，詹洁晖，2000；祝意青等，2012)对比说明重力点值分析方法与重力段差分析方法的异同，并应用本文提出的重力段差可视化方法、多期重力段差数据热力图、G和C指标量随时间的变化来说明重力场变化的显著性程度.

本文使用首都圈地区流动重力测网结果，地理范围为(39°N—41.5°N，113°E—117.5°E)，包括北京、河北、天津、山西等地.首都圈测网位于中国大陆北部的华北构造区，区内以NNE向坳陷和隆起带为主，其北部以呈E-W向分布的燕山—阴山构造带为界，南部以秦岭—大别山构造带北缘为界，西部边界包含银川盆地，东部以渤海为界(顾功叙等，1997).首都圈地区由于受断层形变影响，较大程度上反映了张家口—渤海构造活动带的活动情况(张鹤等，2010)，研究区内郯庐断裂带沿冀北和华北断块东部边缘伸展，并且在较老的大规模隆升带和盆地上发育有小规模的地堑地垒构造(顾功叙等，1997).首都圈的断裂分布及流动重力联测路线如图 3所示，图中绿色方框所圈出的部分为数据量较多且测点较为密集的区域，下文所绘流动重力段差分析图即针对这一区域.

图  3  首都圈地区重力联测路线及断裂分布

Figure  3.  Distribution of gravity survey routes and fault structures in the capital region

下载: 全尺寸图片 幻灯片

首都圈地区每年进行两期测量，测量时间分别集中在春季和秋季.由于秋季降水量较春季小，受雨水影响较小，所以本文选用首都圈重力测网2010—2015年每年第二期的测量数据，复测时间为每年8—11月份，共6期测量平差结果和段差变化.测量过程中，2010—2012年使用的是LCR-G型重力仪，2013年开始使用CG-5仪器，观测精度均优于5×10^-8 m/s².数据处理使用地震局广泛推广使用的高精度重力测量资料处理系统LGADJ软件(李辉等，2009；祝意青等，2012), 并将测网内相对稳定的测点作为拟稳点进行拟稳平差，统一起算基准(祝意青等, 2010, 2012)；对各期观测数据计算结果整体分析，了解仪器观测精度后，合理确定仪器的先验方差，从而得到最佳结算结果.

3.1   重力场变化的段差表示

将相邻两期平差结果进行差分得到2010—2015年期间一年尺度的重力场时空变化，并用重力段差表示法得到一年尺度的段差变化大小和方向，图 4和图 5为截取数据较多的区域绘制而成.现以2010—2011期次的测量平差结果和段差变化为例，对比流动重力段差分析法与重力点值分析法，如图 4所示.可以看出：测网重力变化较小，区域重力场变化较为平缓，重力场自西北向东南呈现正、负、正的有序性变化，其中延庆南部测点的箭头都向外发散，这种情况可能是由于单个测点的环境变化所导致的四周测点产生较大的段差变化.重力等值线由重力点值分析方法插值得到，从图 4可以看出，测网北部观测数据很少，但是由于插值却出现了许多条重力等值线，产生了虚假的重力异常值.在流动重力段差图中，区域的重力场变化通常呈连续的趋势性变化，但图 4中宣化东部的两个测点的箭头长度明显大于观测误差，这是由于此处的两个测点周围的测点非常少，无法反映测网的重力场趋势性变化，因此这两个测点所产生的重力异常具有很大的不确定性.通过对比流动段差分析法与重力点值分析法，可以清楚地看到流动段差分析法有助于分析观测数据中不可靠的部分，便于筛选、剔除数据，尤其是对重力测量数据较为稀疏区域的作用更为突出.

图  4  2010—2011年首都圈地区重力场差分变化

图中绘图数据为2011.9与2010.9两期的差分结果；所显示地震为1970年1月至2014年10月期间发生的地震

Figure  4.  Differential variation of gravity field in capital area during 2010 to 2011

The data of the figure is from the difference of gravity measurement between the issues 2011.9 and 2010.9, and the earthquakes of capital region from 1970 to 2014 are also shown

下载: 全尺寸图片 幻灯片

图  5  首都圈地区流动重力变化段差结果结果

图中所显示的地震为1970年1月至2014年10月期间发生的地震
(a) 2011—2012年; (b) 2012—2013年; (c) 2013—2014年; (d) 2014—2015年

Figure  5.  The gravity difference results of mobile gravity measurement in capital region

The earthquakes of capital region from 1970 to 2014 are also shown
(a) 2011—2012; (b) 2012—2013; (c) 2013—2014; (d) 2014—2015

下载: 全尺寸图片 幻灯片

根据上述流动重力段差分析方法，对首都圈地区2011—2015年的多期流动重力段差变化数据进行分析，结果如图 5所示.由图 5a，b中箭头与椭圆的相对大小可知，2014年涿鹿地震前，整个测网重力场变化较为平缓，段差变化较小.由于涿鹿M_S4.3地震发生在2014年9月6日，而2014年第二期大部分测点是在9月6日之后测量的，因此，2014—2013期的流动重力段差分析图像主要反映的是涿鹿地震震后的情况，而延庆附近区域的测点是在涿鹿地震前测量的.图 5c中，延庆周围测点箭头长度所反映的段差变化异常程度相比于2012—2013年显著提高，由于这些测点是在涿鹿地震前测量的，因此这些显著性异常程度的提高可以反映涿鹿地震前物质不断向震源区发生迁移的过程；而涿鹿附近的测点是在2014年9月6日涿鹿地震之后测量的，该区域段差增加方向不同于2012—2013年，箭头指向较为发散，地下物质向四周各方向运移.涿鹿M_S4.3地震后(图 5d)，2015—2014期次的重力场变化由东向西呈现出由正到负的变化，整个测网的段差变化异常程度较2014—2013期次明显降低，箭头也较为发散，反映了地震后物质不断调整的过程.

3.2   指标量分析

由上述重力场变化的两个评价显著性异常的指标G和C值的定义，将首都圈测网2010—2015年每年第二期(8—11月份)的测量结果进行计算得到测网整体的重力场变化的显著性异常程度，如图 6所示.由G值和C值的变化曲线可以看出，G值和C值随时间的变化趋势基本一致，2014年前均呈现出较为平缓的趋势，在2014年达到最大值，之后则大幅快速降低，回落至与2013年前的数值持平或更低. 2014年之前，即M_S4.3涿鹿地震发生前，G值和C值比较稳定或表现为较小幅度地均匀增加；在涿鹿M_S4.3地震(红色虚线)发震前、后的测量结果表明，反映测网显著性变化测段比例的C₁(大于1倍误差的测段比例)、C₂(大于2倍误差的测段比例)和异常指数G值突然升高，相比之前增幅很大，达到最大值，C₁在2013年突然升高一直持续到2014年；2015年曲线迅速回落，G值和C值均大幅降至地震发生前的水平或稍低于之前的水平.

图  6  2010—2015年G值和C值变化曲线

Figure  6.  The change curves of the indices G and C during 2010 to 2015

下载: 全尺寸图片 幻灯片

通过以上分析可以得出，重力场变化的评价指标量G值和C值变化曲线可以较好地反映地震前后的重力场变化特征.震前，两个指标量均表现为大幅上升的变化趋势，可以观测到明显的异常；震后两个指标量均呈急剧下降.

3.3   多期数据的可视化分析

上述对重力场动态变化图像的分析将每两期差分的结果分别绘制成图，本节则在此基础上，将首都圈测网2010—2015年得到的多期流动重力数据结果进行可视化分析，从而获得测网整体随时间的变化情况，明确测网中哪些测点随时间更易发生重力场变化，以及哪些区域的测点随时间的变化比较稳定.

现将首都圈测网2010—2015年的流动重力数据结果以热力图的形式表示.由于测段较多，为图示说明方便，将首都圈测网人为地分成几个测线分别进行热力图绘图，并选取其中的测线A和B分析说明.测线A和B的分布如图 3所示，热力图结果如图 7所示，其中横轴所表示的测量期次①为2011—2010数据的差分结果，以此类推.

图  7  测线A (a)和B (b) 2010—2015年重力段差变化的时序表示法

Figure  7.  Time series representation of gravity difference of the measuring lines A (a) and B (b) during 2010 to 2015

下载: 全尺寸图片 幻灯片

可以看出：测线A整体变化较为平稳，在2010—2015年期间的段差变化值均不大，大部分测段为负值变化；测段1，3，4，6随时间呈负正交替变化的规律，并且在2014—2013期次④，即涿鹿地震发生期间，大多数测段的变化达到最小，在零值附近.测线B整体变化非常剧烈，测段随时间以正值变化为主.从空间上看，变化较为剧烈的测段主要出现在测段2，5，6，7，即靠近涿鹿地震震源区的相关测段；从时间上看，变化较大的时期出现在2013—2012③和2014—2013④期次.在2015—2014⑤期次，即涿鹿地震后，测段又呈现出很小的负值变化.

综上可得，2010—2015年首都圈测网在靠近涿鹿M_S4.3地震震源区，即涿鹿地区，测线所得到的重力场变化比较显著.从时间上看，许多测段由2013—2012③期次的负值转变为2014—2013④期次的正值.但由于变化量级较小，且无法排除由于降雨或地下水变化产生的重力场变化因素，在2015—2014⑤期次测段基本呈很小的负值变化，反映了等效的物质亏损状态.因此，热力图可以将多期流动重力数据整体进行绘图表示，并且其分析结果与流动重力段差变化图所反映的结果具有很好的一致性.

4.   讨论与结论

近年来，以地震行业需求为目的的地震流动重力测量工作发展迅速，全国性的重力测量网络已经在中国大陆初步建成.这种测点空间分布不规则、疏密程度不一致的观测系统产出的测量成果如果仍简单地套用等值线方法来表示区域重力场变化，将会出现很多问题.本文在分析基于重力点值的重力等值线图示方法的基础上，考虑到重力测点空间分布不均匀、地表观测的重力场变化信号复杂、重力场时间变化量与测量误差之间的区别不显著等问题，提出了基于重力段差值的重力场变化可视化方法，并且定义了两个指标量以衡量重力场的变化程度.基于实际数据测试和对比分析，主要结论如下：

1) 本文提出的基于流动重力段差的可视化方法，适合于表示时空分散性较强的重力场动态变化信息，能够避免由于空间插值或测量误差引起的虚假信号.

2) 本文定义的G值和C值两个指标量，可用于评价区域重力场变化的显著性程度.根据对涿鹿M_S4.3地震前后指标变化的研究，该地震发生之前的很长一段时期内G值和C值随时间的变化平稳；震前，G值和C值出现突增；震后，G和C值又迅速回落.

3) 热力图的可视化方法适用于表示多期流动重力数据，其结果与流动重力段差表示方法给出的绘图结果有很好的一致性.

综上所述，本文提出的基于重力段差的可视化方虽然无法用于预测地震的孕育发生过程，但是段差分析方法可以作为重力等值线分析方法的一个较好的补充；尤其是对于重力测点分布稀疏、数据较少区域的重力等值线，该方法有助于排除与插值方法和测量误差等相关的虚假异常信息，从而辅助判断重力变化等值线的可靠性.本文定义的G值和C值两个指标量可以定量化地反映区域性重力场变化的显著性程度，用于研究重力场变化现象与多种地球动力学因素之间的因果关系.由于重力测量是一种物理测量，除了与地表变形和质量变化过程相关外，还可以反映来自地壳深部的质量变化信号，是多种信号源的叠加.虽然本文提出的可视化方法可以很直观地给出区域重力场的变化特征，但是影响重力段差变化的信号来源很多，如何分离不同场源的重力信号本文并未讨论，但对于解释重力场变化问题十分重要，需要相关的深入研究.