华南陆地时变重力观测数据质量评估

杨锦玲，陈 石，王林海，卢红艳，李红蕾，张 贝

1. 中国地震局地球物理研究所，北京 100081； 2. 福建省地震局，福建 福州 350003； 3. 北京白家疃国家地球科学野外观测研究站，北京 100095

地震、地下水、火山活动及各种构造运动和地壳垂直变形等地球动力学过程，都会引起一定程度的地球重力场时空微变化。高精度时变微重力监测手段常用于观测火山活动期间岩浆的运动[1-2]、地下水储量变化[3]、油气库含水层[4]及地震危险性分析[5-9]。与卫星重力[10]和海洋重力测量[11]相比，高精度的地表重力重复观测具有距离地壳内部场源近、观测位置可重复性强、观测仪器精度高等特点，适用于监测与地壳内部场源直接相关的微重力信号。

21世纪，依赖绝对重力观测技术逐渐成熟，为研究重力场时变信号提供了更加可靠的时空基准[12-16]。例如：2004年至2009年的绝对和相对重力联合在地表重复观测到的重力场变化与菲律宾板块俯冲引起的慢地震相关[5]；2008年汶川8级地震和2015年尼泊尔7.8级地震从绝对重力重复测量结果提取到震前重力变化信号[13,15]。但是，由于陆地重力观测特别是时变重力场信号的获取，仪器观测能力限制和场源非唯一性仍是当前制约重力方法研究深部介质变化的核心难题。重力信号对质量变化敏感，微伽级重力数据常受到观测方法、仪器精度以及潮汐、环境、地下水和地表沉降等多种不确定性因素的影响，这也使得陆地重力观测得到的重力场变化存在较多不确定性。

陆地重力观测仪器的不确定性主要来源于弹簧重力仪的零漂特性和格值系数[17]。目前，国内外陆地重力观测的主流设备以LaCoste&Romberg型、CG-5型和Burris型弹簧重力仪为主，但是弹簧型重力仪不具备绝对观测能力，且CG-5相对重力仪每天零漂达50～200×10-8m/s2[18]。对于较大空间尺度和较长野外作业时间的重力观测，非线性漂移信号的估计是提高重力信号解算精度的关键。文献[19]提出了一种专用于陆地高精度微重力观测的贝叶斯平差方法，该方法可以对多仪器联合观测自动定权、自适应量化重力仪非线性漂移观测误差，为提高重力平差计算的精度提供了一种新的解决方案。

华南地区历史上中强地震活跃，其近岸海域有6～7级地震的强震构造背景。在华南地区开展高精度的时变重力观测，对于识别潜在强震风险源和研究强震孕育全周期过程的微重力时空变化特征具有重要的科学意义。此外，全面对华南相对重力仪器非线性漂移和格值系数等不确定性因素进行量化评价，在此前还基本属于空白，但该工作是进一步应用地震重力资料分析区域地震危险性的基础，因此本文将采用最新的贝叶斯平差方法，对华南地区近年的时变重力观测资料进行评价，为后续资料的分析和应用提供科学依据。

本文以华南地区2015年—2018年的时变重力观测资料为基础，引入贝叶斯平差技术，给出华南地区的重力观测数据精度分析，提出通过重力段差残差和绝对重力交叉验证方法来评定重力观测资料的精度指标方案。第2节详细介绍了华南构造背景和时变重力测网概况；第3节给出了研究方法原理；第4节给出了华南地区的重力数据质量评价；最后给出了研究结论。本文研究结果对评价陆地时变重力监测数据质量，更合理和客观地解释重力场变化与场源介质属性之间的关系，具有较好的参考意义。

1 华南构造背景与时变重力观测网络

华南位居太平洋西缘，其南西侧与东南亚块体接触，南东侧为西太平洋构造区[20]。现今华南地区由扬子和华夏地块经过多期复合演变形成现今的基本形态[21]，其大地构造与演化受欧亚板块、太平洋板块和菲律宾海板块共同作用的影响，动力学背景复杂，研究区内发育多条活动断裂。华南地区历史地震活动强烈，1604年曾发生泉州海外8.0级地震。2010年以来该区域以中强地震为主，先后发生了2013年广东东源5.1级地震、2013年福建仙游5.2级地震和2018年台湾海峡6.2级地震等多次中强地震。图1(a)显示了2010年以来的地震主要集中在长乐-诏安断裂带、邵武-河源断裂带以及滨海断裂带。

时变重力观测网有助于获取较高时空分辨率的重力场信息,并有效监测地震重点危险区地球物理场的动态变化特征，以识别区域潜在强震风险源，此前已在青藏高原中强地震[7-9]的中长期危险区预测和地球科学研究中发挥了重要的作用。中国地震局自20世纪70至80年代开始在华南地区开展时变重力观测。经过多次升级和改造，至2019年华南地区已建成了测点分布较均匀、有绝对基准控制的高质量陆面时变重力观测系统，并基本覆盖福建和广东陆地主要的构造和断裂带，如图1(b)所示。由图1(b)可知，测网内有厦门、武夷山、平潭、广州和韶关绝对重力基准点，平均两年重复观测一次，精度优于5×10-8m/s2。相对重力重复测量点有550个，测段640段，测点间距在30～50 km范围内，平均6个月重复观测一次，采用4台CG-5型相对重力仪，两两一组进行复测。

2 贝叶斯重力平差模型

陆地重力数据处理的核心是平差计算[22]。目前常用的经典重力平差方法采用最小二乘方法获得观测点的最佳估值，并假定参与观测的每台重力仪漂移率是线性的[23-24]。经典平差方程可表示为

(1)

(2)

式(1)和式(2)也可合并简化为

(3)

式中，Normal(·)表示正态分布，均值为0，方差为W-1；A、D、G分别为观测顺序矩阵、段差观测时间矩阵和绝对点矩阵；x、v代表点值和漂移率，为待求参数；y为观测段差。

SX=Y

(4)

式(4)为绝对重力约束下的观测方程，Δy为同一台仪器在相邻两测点的重力读数差，T为理论固体潮，α为潮汐因子；P为大气压力负荷引起的重力变化，β为气压导纳。

通过求上述平差方程的最大似然估计，可以得到x、v的最优估计。首先，可以由式(1)—式(3)，给出多台相对重力仪和绝对重力观测的联合概率密度分布

(5)

区别于经典平差对漂移的线性假设以及分段零漂的处理方式[25]，本文采用的贝叶斯重力平差方法[19]提供了非线性漂移的计算模型。该方法假定相对重力仪的漂移率随时间变化光滑，将其作为先验条件，可表示为

(6)

(7)

式中，l为仪器待求格值系数。

基于贝叶斯公式，后验概率似然函数可表示为

(8)

(9)

式中，L(y|v)可由方程式(5)得到；π(v)为漂移率的先验分布，并采用ABIC方法[26]进行贝叶斯估计模型的拟合优良性分析。

对于相对重力仪的格值优化问题，测网中如果同时存在两个或多个独立观测的绝对重力点，通过上述贝叶斯优化计算，即可以实现格值系数的优化计算[27]。

不同区域的数据应根据实际情况选择不同平差方法进行数据处理。对于小范围的重力测量，按经典平差的线性漂移模型解算即可。但对于双程闭合时间长，空间测点地理跨度大的华南区域流动重力测量，则需重点考虑仪器的非线性漂移特性，以提高数据精度。因此本文重点对比经典与贝叶斯平差方法在华南区域的差异和方法的适用性。

3 华南时变流动重力观测质量评价

在陆地时变重力数据处理中，相对重力仪的非线性漂移特性和格值系数变化是两个主要的误差来源。本文分别采用经典平差和贝叶斯重力平差方法处理了华南测网2015—2018年的重力观测数据，基于段差残差、段差互差和绝对点值差参数，量化相对重力仪非线性漂移特征和格值系数的不确定性，并基于绝对重力基准观测的交叉检验评价华南时变重力数据质量，分析区域的重力场时空特征。

3.1 仪器不确定性评价

3.1.1 漂移率变化

选取2018 C1期华南测网的观测数据进行分析，该测期绝对重力观测有厦门、武夷山、平潭、广州和韶关5个测点，标准差优于5×10-8m/s2，与流动重力观测时间准同步。分别采用CG5-1316、CG5-814(福建测网)和CG5-232、CG5-369(广东测网)4台相对重力仪进行同步观测，并对广福、蕉岭等公共测点进行联测。为便于讨论，下文以FJ-1316、FJ-814、GD-232和GD-369命名4台仪器，将3月和9月的测期分别命名为C1期和C2期。

图2为4台相对重力仪的线性和非线性漂移率变化。图2可见4台相对重力仪存在非线性漂移，且彼此间差异性明显。仪器的漂移特性随测量进程呈不同的变化趋势。FJ-1316的漂移率达0.8×10-5(m/s2)1 d(图2(a))，观测周期内漂移率最小和最大值相差约0.15×10-5m/s2，其每日漂移率表现为逐渐增大后趋于稳定的过程，而FJ-814的漂移率变化较稳定。随着仪器观测时长增加，GD-232在测量后半程的非线性漂移率明显大于测量前期(图2(b))。

图1 华南地区地震地质构造背景与时变重力观测网络Fig.1 The geologic structure seismic activity and time-varying gravity observation network in South China

图2 华南测网弹簧重力仪的漂移率变化Fig.2 Estimated drift rates of four gravimeters at Fujian and Guangdong survey campaign

图3为采用经典和贝叶斯重力平差方法计算的段差残差变化。平差后的段差残差基本分布在±30×10-8m/s2范围。贝叶斯重力平差方法的残差(图3(b))基本呈随机分布状态，幅值略小于经典平差结果(图3(a))。4台仪器的段差残差统计直方图(图3(c))显示出贝叶斯平差方法中小于±5×10-8m/s2的段差残差数量多于经典平差法。两种方法的段差残差分布都符合高斯分布特征，经典平差的段差残差标准差为6.8×10-8m/s2，贝叶斯平差的段差残差标准差为6.17×10-8m/s2。从段差残差结果和统计参数特征看，在华南测网贝叶斯重力平差结果略好于经典方法。

3.1.2 格值系数变化检测

相对重力仪格值是影响重力测量误差的重要因素。目前地震重力观测中相对重力仪的格值系数通常采用基线场标定结果[28]，但基线场标定间隔时间长[27]，其周期长达3～5 a，不能满足重力观测半年周期的需求。本文中所采用的4台重力仪，FJ-814和FJ-1316初始格值由仪器厂家给定。GD-232查无基线场标定记录，GD-369在2010年进行过基线场标定，标定格值为1.000 456。由于福建测网由厂家给定的初始格值是在实验室观测条件下获得，广东测网的基线场标定周期过长，因此每次开展重力观测前仪器都进行了实地测区的格值标定。但实际测区的格值标定仍存在着由于段差过小、观测人员操作差异等不确定性因素，为此本文应用贝叶斯重力平差优化模型，通过高精度的绝对重力基准点，在平差方程中采用贝叶斯优化方法进行格值系数估计，依据ABIC准则寻找最优超参数，进而得到绝对重力约束的格值系数。

FJ-1316和FJ-814在研究时段格值系数不变(图4)，FJ-1316和FJ-814的初始格值分别为1.000 655和0.999 934，GD-232的格值系数变化在2016 C2至2018 C1的格值系数为0.998 610，GD-369在2016 C2至2018 C1期间为0.999 902。由于研究时段内仅有两期绝对重力观测(2015年12月和2018年1月)，因此2016 C1至2017 C2这4期选用2015年12月的绝对重力观测为基准，2018 C1以2018年1月的绝对重力观测作为基准。2016 C1和2018 C1两期绝对和相对观测时间准同步。

基于各期的绝对重力点基准，采用贝叶斯重力平差方法计算2015年—2018年4台仪器优化的格值系数。表1列出了各期的优化格值系数和格值偏差(即初始格值与优化格值之差)统计。图4为华南测网相对重力仪初始和优化格值系数随时间的变化，图中用紫色虚线标示出绝对重力观测的时间。华南测区贝叶斯格值优化的格值系数与格值偏差统计见表1。如表1和图4所示，研究时段内福建测网的格值系数偏差明显大于广东测网。FJ-1316的优化格值与初始格值的平均偏差约6.28×10-4，最大偏差约7.02×10-4。广东测网内GD-369和GD-232的优化格值与初始格值的平均偏差为4.6～6.23×10-5，最大偏差为2017 C2期GD-369的2.26×10-4。从图4亦可见，在2016 C1—2016 C2期和2018 C1期，当绝对重力和相对重力观测时间近于同步时，GD-232和GD-369的优化格值系数和初始格值系数更为接近，一致性较好(二者的差值小)。

表1 华南测区贝叶斯格值优化的格值系数与格值偏差统计

福建测网格值偏差较大的主要原因主要在于仪器性能及标定周期过长。从图5段差互差的相关性统计可见，福建重力测网内主要段差分布在±80×10-5m/s2范围，而仪器格值系数出现5.89×10-4的偏差，将可能引起47×10-8m/s2的误差，相比仪器非线性漂移变化，格值系数的误差可能严重影响整网平差结果的准确性和可靠性。

虽然FJ-1316和FJ-814仪器的优化格值与初始格值偏差较大，但二者格值偏差量级接近(表1)。以2018 C1期为例，格值优化前后福建与广东段差的互差范围都分布在±0.025×10-5m/s2范围，仍符合地震重力测量规范的要求。从优化前后段差互差相关性(图5)看，格值优化后段差互差的变化不明显。因此，当两台仪器格值同步发生变化，仅依据仪器之间的互差不足以评价结果可靠性，需要进一步通过测网内的绝对重力测量结果开展交叉验证。

图3 段差残差图Fig.3 Gravity difference residuals

图4 华南测网相对重力仪初始和优化格值系数的变化(图中虚线为绝对重力观测时间)Fig.4 Changes of calibration scale factor of relative gravimeters in Huanan network

图5 分别采用初始和优化格值系数时测点段差与仪器互差变化Fig.5 Gravity difference of two adjacent stations and gravity difference of the simultaneous observation instruments in the case of optimized and calibrated scale factor

在既有绝对观测也有相对观测的混合重力测网中，高精度、独立的绝对重力观测可用来验证相对重力的观测结果，因此进一步采用绝对重力交叉验证的方法对格值优化结果进行评估。交叉验证(cross-validation)基本思想是将原始数据进行分组，一部分作为训练集来训练模型，另一部分作为验证集测试评价训练模型。首先将绝对重力点分组，选择通过测网中部分绝对点作为平差起算点参与平差，解算得到相对重力推算的绝对重力点值，并与绝对重力实测的绝对点值相减得到绝对点值差。因此，绝对点值差可以作为相对观测与绝对观测一致性的检验指标，反映重力仪格值的准确程度。

平潭和韶关作为起算点，厦门、武夷山和广州作为检验点(对应图6(b))的绝对点值差见表2。由表2可知，优化格值的绝对点值差相比初始格值的绝对点值差更小。其中，厦门和武夷山变化较为显著，优化格值的绝对点值差减小20～60×10-8m/s2，即采用优化格值推算的重力值更符合用于检验的绝对重力实际测量值。因此，与仪器间段差互差的相关性方法相比，采用绝对重力交叉验证格值准确性的方法更有效。

表2 华南测区贝叶斯格值优化前后的绝对重力抽样检验

3.2 绝对重力交叉验证

在上文基础上，对测网内5个绝对点依次进行交叉验证，以全面检验贝叶斯重力平差方法的有效性和可靠性。以2018 C1测期准同步的绝对重力观测，设计4组方案(表3)对相对重力平差结果进行交叉验证。

表3 绝对重力交叉验证

图6是分别基于表3方案采用经典和贝叶斯重力平差方法得到的4组绝对点值差。由图6可知，贝叶斯重力平差方法(采用优化格值)的结果明显优于经典平差(采用初始格值的)的结果，这与仪器格值误差得到有效抑制有关。而随着用于平差的绝对点增加，两种方法的点值差较差逐步变小。在图6(a)中，采用平潭1个点估算时，两种平差方法计算的点值差较差为20～60×10-8m/s2。而当采用平潭、韶关两个点估算时(图6(b))，位于测网中部的广州和厦门贝叶斯重力平差相比于经典平差的点值差都减小约20×10-8m/s2，测网边缘的武夷山改善更显著，其贝叶斯重力平差相比于经典平差的点值差减小57×10-8m/s2。从图6(a)与6(b)经典平差的武夷山点值差的变化亦可见，受FJ-814和FJ-1316两台CG-5弹簧重力仪的格值系数不准确的影响，相对重力平差的结果与绝对重力观测差异较大。在图6(c)和6(d)中，两种方法的平潭绝对点值差较差约10×10-8m/s2。综合以上，测网内绝对重力的交叉验证结果表明贝叶斯重力平差方法得到的优化格值系数比初始格值系数更准确。

为进一步研究不同平差方法和绝对点约束下的点值结果空间差异性，采用表3中方案b—方案d分别对华南测网中不同绝对点起算下贝叶斯重力平差相比于经典平差的点值差空间分布进行了可视化，如图7所示，可知福建测网的因平差方法产生的影响明显大于广东测网，笔者认为这主要与测量仪器的格值误差有关(表1)。值得注意的是，随着用于平差的绝对点增加，平差方法引起的结果差异逐渐减小，且在测网边缘处的点值精度改善最明显。因此，重力观测网内绝对点的布置应尽量均匀，且为保证点值精度，至少有两个以上的绝对点进行平差。当选择平潭和韶关两个绝对点为起算点(图7(a))时，测网边缘的武夷山点值差都明显大于其他点。219个测点的点值差为10～20×10-8m/s2，75%测点点值差分布在±30×10-8m/s2范围，均值为20×10-8m/s2。当测网内有韶关、武夷山、厦门和广州4个绝对点作为起算(图7(c))时，284个测点的点值差在±10×10-8m/s2。而且点值差的分布更接近正态分布，均值为18.8×10-8m/s2。由图6—图7可知，华南测网观测资料由于相对重力仪非线性漂移和格值系数不准确共可引起约20×10-8m/s2的不确定性，且通过贝叶斯方法可以有效抑制。

图6 不同起算点下经典和贝叶斯平差的绝对点点值估算误差Fig.6 Differences between estimated gravity values at each absolute gravity station by classical adjustment and bayesian adjustment approach when different absolute gravity stations are using in the base observations

图7 不同绝对点起算下贝叶斯重力平差相比于经典平差的点值差空间分布Fig.7 Differences between estimated gravity values by classical adjustment and Bayesian adjustment approach when different absolute gravity stations are using in the base observations

4 结 论

陆地时变微重力观测是研究地壳内部物质变化的重要手段。华南重力网具有测点间距大、覆盖范围广、观测时间长和对重力资料解算精度要求高等特点。本文采用华南地区2015年—2018年的实测流动重力观测资料，引入贝叶斯平差方法，量化分析了由于相对重力仪的非线性漂移、一次项系数不准确的影响，得到以下结论。

(1) 华南重力测网使用的4台相对重力仪存在非线性漂移，且彼此间差异性明显，其中福建测区使用的FJ-1316仪器漂移率最高达0.8×10-8m/s2/d。从段差残差结果和统计参数特征看，贝叶斯平差结果略好于经典方法。

(2) 对仪器格值系数误差的评估结果发现，福建测网FJ-1316和FJ-814的初始与优化格值偏差平均值为5.89～6.28×10-4，广东测网GD-232和GD-369初始与优化的格值系数偏差平均值为4.2～6.23×10-5。应用绝对重力交叉检验法对格值优化结果进行评价，优化后的检验绝对重力点值差可减小约20～60×10-8m/s2重力变化，采用贝叶斯方法优化后格值的平差结果更符合用于检验的绝对重力测量值。值得注意的是，用于福建测网中的两台相对重力仪，由于两台仪器格值系数同步变化且相互间差异较小，导致仅从仪器间互差来验证结果可靠性很难发现存在的格值系数误差。

(3) 通过测网中的绝对重力基准开展交叉检验，随着用于平差的绝对点增加，经典和贝叶斯平差方法的差异逐渐减小，且在测网边缘处的点值精度改善最明显，分析认为这主要是由于仪器格值误差得到有效抑制。华南测网观测资料由于相对重力仪非线性漂移和格值系数不准确共可引起约20×10-8m/s2的不确定性，且通过贝叶斯方法可以有效抑制。

本文基于最新的贝叶斯平差方法，着重从相对重力仪的不确定参数优化(漂移非线性特征、格值系数)和绝对重力交叉检验技术对华南地区的时变重力数据进行了分析。本文研究方法和思路同样适合于中国大陆其他地区的重力数据质量评价，也有助于更合理和客观地开展重力场变化与断层运动、地壳形变、强震孕育和物质变迁等方面的研究。