地球重力场的空间异常特征和动态时间变化可为地球内部构造和动力学研究提供重要的基础信息。区域重力场与地壳运动、环境变化和深部物质运移等因素密切相关^[1-3]。随着重力观测仪精度的提高以及观测技术的不断发展，利用重力测量获得区域重力场随时间的非潮汐变化，成为研究地震孕育、发生和调整以及地球内部物质运移过程的重要方式^[4-6]。

以监测地壳运动及研究地震活动相关的重力场变化为目的，定期开展固定点位高精度重力重复测量的工作，称为地震重力测量或流动重力测量^[5]。重力测量的数据需经过平差获得各测点的重力值。其中，经典平差和拟稳平差是常用的两种地震重力网平差方法。经典平差以测网内已知重力值为起算基准；拟稳平差将测网中的重力点分为相对稳定和不稳定两部分，以相对稳定点的重心作为起算基准^[6]。对于有同步观测的绝对重力点分布的测网，经典平差方法最简单实用。然而，由于历史资料及部分重力网的绝对重力点分布较少，测网的绝对重力控制得不到保证，因此拟稳平差在地震重力网平差中同样被广泛使用 ^[1，5]。

重力网平差对地震重力数据处理的可靠性起着至关重要的作用^[7]。不同的平差方法获得的结果并不一致，甚至可能超过重力值随时间的变化量，成为研究重力场非潮汐变化的重要干扰之一^[8-9]。然而，目前对地震重力网经典平差和拟稳平差的深入研究较少，因选取与重力网不匹配的平差方法而导致所得重力场变化与真实结果不符的情况时有发生。因此，本文以南北地震带南段重力网的测量数据为基础，通过定义重力点值变化的指标量，定性与定量地分析经典平差和拟稳平差点值精度、段差残差和点值变化的差异，并探讨经典平差的起算基准及拟稳平差的拟稳点对区域重力场变化的影响，为地震重力网平差提供参考依据。

1 平差原理及数据处理

1.1   重力网经典平差和拟稳平差原理

地震重力测量是沿一条测线的同址往返测量，相邻两测点间的重力段差是测量的基本单元和重力网平差的最基本元素。根据图 1所示的测量过程，可得出误差方程和权矩阵为^[5]：

图  1  地震重力测量示意图

Figure  1.  Schematic Diagram of Mobile Gravity Measurement

下载: 全尺寸图片 幻灯片

式中，$\mathit{V}=(V_{\mathrm{2.1}}{V}_{\mathrm{3.2}}\cdots V_{\mathrm{2}i-\mathrm{1.2}i-\mathrm{2}}{)}^{\mathrm{T}}$为误差向量；$\mathit{A}$为系数矩阵；$\mathit{X}=(g_{\mathrm{1}}^{\text{'}}{g}_{\mathrm{2}}^{\text{'}}\cdots g_i^{\text{'}}{)}^{\mathrm{T}}$为未知数向量；$\mathit{L}=(∆g_{\mathrm{2.1}}∆g_{\mathrm{3.2}}\cdots ∆g_{\mathrm{2}i-\mathrm{1.2}i-\mathrm{2}}{)}^{\mathrm{T}}$为观测值向量；$\mathit{D}(\mathit{L})$为观测方程协方差矩阵；$\mathit{P}$为权矩阵；${\mu }^{\mathrm{2}}$为单位权方差。

若测网存在起算基准，即附加了约束条件，根据${\mathit{V}}^{\mathrm{T}}\mathit{P}\mathit{V}=\mathrm{m}\mathrm{i}\mathrm{n}$原则，式（1）中的系数矩阵$\mathit{A}$列满秩，则未知数$\mathit{X}$有唯一解：

式中，$\mathit{N}={\mathit{A}}^{\mathrm{T}}\mathit{P}\mathit{A}$；$\mathit{W}={\mathit{A}}^{\mathrm{T}}\mathit{P}\mathit{L}$；n、t分别为观测方程个数和未知数个数。

若测网无起算基准，则式（1）系数矩阵$\mathit{A}$秩亏，解不唯一，此时可采用拟稳平差。设不稳定点的未知数为${\mathit{X}}_{\mathrm{1}}$，拟稳点的未知数为${\mathit{X}}_{\mathrm{2}}$，则式（1）变为：

当附加最小范数条件${\mathit{X}}_{\mathrm{2}}^{\mathrm{T}}{\mathit{X}}_{\mathrm{2}}=\mathrm{m}\mathrm{i}\mathrm{n}$时，由式（4）的法方程，可求得拟稳点的未知数${\mathit{X}}_{\mathrm{2}}$以及非稳点的未知数${\mathit{X}}_{\mathrm{1}}$^[9]：

式中，${\mathit{a}}^{\mathrm{T}}={\mathit{A}}_{\mathrm{2}}^{\mathrm{T}}-{\mathit{N}}_{\mathrm{21}}{\mathit{N}}_{\mathrm{11}}^{-\mathrm{1}}{\mathit{A}}_{\mathrm{1}}^{\mathrm{T}}$；${\mathit{M}}_m^{-}$为最小范数逆矩阵；${\mathit{N}}_{\mathrm{11}}={\mathit{A}}_{\mathrm{1}}^{\mathrm{T}}\mathit{P}{\mathit{A}}_{\mathrm{1}}$；${\mathit{N}}_{\mathrm{12}}={\mathit{A}}_{\mathrm{1}}^{\mathrm{T}}\mathit{P}{\mathit{A}}_{\mathrm{2}}$。

1.2   数据预处理及一次项系数评定

本文选用南北地震带南段重力网2015—2016年数据进行经典平差和拟稳平差研究，相对重力和绝对重力测量分别由CG5型相对重力仪和FG5型绝对重力仪完成。重力网分布、各重力点观测次数及单程观测时间如图 2所示。

图  2  研究区地震重力网及测量概况

Figure  2.  Distribution of Mobile Gravity Network and Measurements in Study Area

下载: 全尺寸图片 幻灯片

地震重力测量受环境变化等因素影响，测量数据中存在偶然误差和系统误差。由于两期测量采用的仪器一致、观测人员相同，可有效地控制不同观测人员观测技术及习惯等引起的偶然误差^[10]。为消除或削弱系统误差，本文对相对重力数据进行固体潮、一次项系数和零漂等校正^[1]，绝对重力数据进行固体潮、气压、极移和垂直梯度等校正。相对重力仪一次项系数测前利用段差202 mGal的贵阳、水城绝对点进行标定（系数Ⅰ），然而测网重力点值跨度达到782 mGal，标定的重力段差并未完全覆盖测网读数段^[5]。因此，本文在数据预处理时以测网8个绝对点为基准对两台仪器的一次项系数进行重新标定（系数Ⅱ），结果和精度如表 1所示。由表 1可知，2015年两台仪器系数Ⅱ相比系数Ⅰ增大1.01×10^-4~1.36×10^-4，变化较为一致；而2016年系数Ⅱ相比系数Ⅰ，CG5 #1229仪器减小3.98×10^-4，CG5 #1235仪器增大1.16×10^-4，变化相差较大。

表  1  相对重力仪及一次项系数标定结果

Table  1.  Relative Gravimeter with Its One Degree Term of Chromatic Polynomial

测量时间	测量时长/d	仪器型号	一次项系数
			系数Ⅰ	系数Ⅱ
2015-06-27—2015-08-15	50	CG5 #1229	0.999 674±0.000 047	0.999 775±0.000 038
		CG5 #1235	0.999 973±0.000 047	1.000 109±0.000 035
2016-08-27—2016-10-02	37	CG5 #1229	0.999 998±0.000 015	0.999 600±0.000 023
		CG5 #1235	1.000 066±0.000 015	1.000 182±0.000 023

下载: 导出CSV 

| 显示表格

本文采用系数Ⅰ、Ⅱ对2016年的重力数据进行平差，平均点值精度分别为11.6 μGal和10.8 μGal。基于系数Ⅰ、Ⅱ计算的重力段差和互差相关性的结果如图 3所示。由图 3可知，系数Ⅱ拟合的线性相关性明显小于系数Ⅰ，说明系数Ⅱ的两台仪器观测值互差较小，系统性误差更低。

图  3  重力段差与互差相关性

Figure  3.  Correlation Between the Section Difference and the Mutual Difference

下载: 全尺寸图片 幻灯片

此外，本文利用测网中1个绝对点为起算基准进行平差，统计计算剩余7个绝对点的平差值与实测绝对值差异，结果如图 4所示。

图  4  起算基准的平差值与实测绝对值差异

Figure  4.  Difference Between the Adjustment Values and the Real Gravity Values at Absolute Gravity Stations

下载: 全尺寸图片 幻灯片

由图 4可以看出，系数Ⅱ相比系数Ⅰ，绝对点的平差值在整体上更接近于实测的绝对值。因此，本文采用整网平差标定的系数Ⅱ进行经典平差和拟稳平差研究，其中，经典平差采用精度为±5 μGal的弱基准。

1.3   重力点值定量参数化指标量

地震重力测量的目的是利用平差后的重力点值提取监测区域多期重力点值的非潮汐变化信息，并识别其中的异常点，为地震趋势预报、地壳形变研究提供依据。在重力场变化分析时，经常采用会引入虚假异常信号的插值与滤波方法获得监测区域整体的时变重力场，忽略了点值误差的影响。因此，本文以重力点值变化叠加点值误差的可视化形式对重力场的变化情况进行客观分析和讨论，并在重力段差指标量^[11-13]的研究基础上，基于重力点值变化及误差，定义两种定量参数化指标量P值和R值，为分析区域重力场变化的显著性程度提供定量依据。

P值的计算如下：

式中，N为有效点值个数；$\mathrm{\delta }{g}_i$为第i个点值的变化；$e_i$为第i个点值变化的误差。指数P反映的是测网重力点值变化的整体显著性程度，当0 < P≤1时，认为测网的区域重力场变化不明显，当P > 1时，则认为测网整体存在显著的重力场变化。

R值的计算如下：

式中，i=1，2…N；N_i为大于i倍误差的有效测点个数；R_i为大于i倍误差变化的测点比例。

2 平差结果

2.1   精度评定

对2016年重力数据以8个绝对点为基准进行经典平差，点值精度如图 5所示。由图 5可知，绝大部分测点的平差精度优于12 μGal，精度的差异性分布与测网网形结构密切相关，支线测点平差精度低于闭合环测点平差精度。对于一个由支线和闭合环组成的测网，根据误差传播定律，支线上第i个测点的平差精度可以表示为${\sigma }_i=\sqrt[]{i}{\sigma }_{\mathrm{0}}$（${\sigma }_{\mathrm{0}}$为每个测段的观测精度）^[13]。对于支线测点来说，平差精度将随着支线测段数量的增加逐渐降低，闭合环则将误差按权重分摊到各个测段，从而避免了误差的积累。因此，为提升重力网平差的精度，在重力测量和数据处理时应尽量将过长的支线连成闭合环。

图  5  重力网经典平差点值精度

Figure  5.  Accuracy of Gravity Stations in Classic Adjustment of Gravity Network

下载: 全尺寸图片 幻灯片

为选取相对稳定重力点作为拟稳平差的拟稳点，本文首先开展测网绝对重力点的稳定性分析。图 6为绝对点重力值的时序变化热力图，由图 6可知，河池、水城、思茅、蒙自、泸州和贵阳的点值变化较为平稳，与其相对稳定的构造特征相吻合。耿马和勐海的点值时序变化较为剧烈，可能反映了其下方较为活跃的构造运动^[14]。因此，确定相对稳定的河池、水城、思茅、蒙自、泸州和贵阳测点为拟稳点。

图  6  绝对重力值时序变化热力图

Figure  6.  Heatmap Chart for the Sequential Variation of Absolute Gravity Values

下载: 全尺寸图片 幻灯片

此外，综合测网近年来重力场的动态变化结果^[15]和各测点的构造背景以及环境变化等信息，进一步选取砚山、关岭、西双版纳和安龙共10个重力点作为拟稳平差的拟稳点。拟稳平差的点值精度如图 7所示，可以看出，测网东部重力点的平差精度明显优于西部测点，这是因为大部分拟稳点及拟稳点的重心分布在测网东部，而拟稳平差以拟稳点的重心为起算基准。根据误差传播定律，远离起算基准的测网西部测点平差精度不及起算基准附近的测点。因此，为平衡拟稳平差点值精度的区域性差异，拟稳点应尽可能相对均匀地分布在重力网中。

图  7  重力网拟稳平差点值精度

Figure  7.  Accuracy of Gravity Stations in Quasi-Stable Adjustment of Gravity Network

下载: 全尺寸图片 幻灯片

段差残差可以直观地反映平差后段差与原始段差的差异，是衡量平差方法质量和有效性的重要指标^[16]。图 8为经典平差和拟稳平差的段差残差和统计直方图，由图 8可知，两种平差方法获得的段差残差波动较小，残差值符合正态分布，并且绝大部分残差分布在±5 μGal以内。两种平差方法计算的重力段差均能较好地还原真实段差值，其中，经典平差的残差标准差为2.99 μGal，拟稳平差的残差标准差为2.06 μGal。经典平差由于受基准点的控制，能更好地控制误差传递，减少重力段差的误差，而拟稳平差以拟稳点重心为基准，缺少控制，因此更与原始测量重力段差值接近。可见，经典平差比拟稳平差能更好地进行段差值修正，更真实地反映重力段差值的客观变化。

图  8  重力段差残差及统计直方图

Figure  8.  Residuals and Its Histograms of the Gravity Section

下载: 全尺寸图片 幻灯片

2.2   重力点值变化

利用绝对点（见图 5）和拟稳点（见图 7）对2015—2016年的重力数据进行经典平差和拟稳平差，获得两种平差方法一年尺度的重力网点值变化和点值变化的差异情况，结果如图 9所示。其中，黑色和灰色重力点由于点位破坏和环境变迁，导致重力点值变化失真，因此数据分析时将不再讨论。

图  9  2015—2016年经典平差和拟稳平差的重力网点值变化

Figure  9.  Gravity Values Variation in Gravity Network of Classical and Quasi-Stable Adjustment Methods from 2015 to 2016

下载: 全尺寸图片 幻灯片

由图 9可知，经典平差和拟稳平差的点值在测网北部和西南部均存在突出的重力值正变化，而勐海周缘测点则以重力值负变化为主；两种平差方法在整体上能较为一致地呈现测网重力场变化的形态和趋势，但在某些局部区域仍然存在一定的差异，最大差异可达30 μGal。区域A经典平差重力值的正变化比拟稳平差剧烈，甚至存在大于3倍误差的重力值正变化；区域B拟稳平差的河池点重力变化值为负值，且附近存在突出的负变化，而经典平差的河池重力值为正变化。结合图 6可知，2016年第2期的耿马、河池点的绝对重力值相比2015年同期大，增量分别为36 μGal、11 μGal。经典平差以起算点值作为其他各点值的计算基准，起算点值的变化直接影响其他各点值的变化；拟稳平差以全部拟稳点的重心为起算基准，某一拟稳点的点值变化对其他测点的点值影响较小，造成了两种平差方法在起算点附近点值变化的显著差异。表 2给出了两种平差点值变化的定量参数化指标量P值和R值，其中，J1和J2分别表示时变和不变基准的经典平差的指标量结果，N1和N2分别表示稳定和不稳定拟稳点的拟稳平差的指标量结果。

表  2  经典平差和拟稳平差的指标量P值和R值

Table  2.  Indexes P and R of Classical and Quasi-Stable Adjustment Method

统计项	平差方法	点/（值）	P	R
				R1	R2	R3
J1	经典平差	基准值（时变）	1.70	0.41	0.09	0.02
N1	拟稳平差	拟稳点（稳定点）	0.55	0.16	0	0
J2	经典平差	基准值（不变）	0.96	0.32	0.06	0.01
N2	拟稳平差	拟稳点（非稳定点）	0.41	0.11	0	0

下载: 导出CSV 

| 显示表格

由表 2可知，经典平差J1的P值和R值均明显大于拟稳平差N1；经典平差J1的P=1.70 > 1，表明测网整体存在较为突出的重力点值变化，而拟稳平差N1的P=0.55 < 1，反映测网整体的重力点值变化不显著。

尽管两种平差方法在重力场点值变化的显著性程度以及局部区域的点值变化上存在差异，但对于绝对点较少或者无绝对点分布的测网，选择合理恰当的拟稳点，拟稳平差的重力场变化趋势整体与经典平差差别不大。这也为处理无绝对点分布的历史重力网资料及绝对点分布较少的重力网资料提供了一种参考。

2.3   起算点/拟稳点对重力点值变化影响

地震重力测量中绝对重力与相对重力测量时间同步，为平差计算提供随时间变化的时变基准是较为理想的情况。实际上由于国内高精度绝对重力仪数量偏少、仪器故障等原因，绝对重力与相对重力的测量时间并不完全同步^[10]。采用经典平差可能会将某期的绝对重力值用于相邻的多期资料中，使得平差的起算基准恒定不变。图 10为采用相同起算基准进行经典平差计算的点值变化，以及与时变基准平差的点值变化的差异。

图  10  经典平差不变基准的点值变化及其与时变基准结果的差值

Figure  10.  Gravity Values Variation in Gravity Network Calculated by Classical Adjustment Method with Fixed Benchmark and Difference Between It and that with Changing Benchmark

下载: 全尺寸图片 幻灯片

由图 10可知，与采用时变基准的重力点值变化相比，测网整体的重力点值变化形态较为一致。但在耿马和泸州等绝对点值变化较大的区域差异性显著，差异大小与绝对值的变化幅度基本一致。此外，结合定量参数指标量（见表 2）可知，相同点值精度情况下，不变基准J2的指标量远小于时变基准J1的结果。J2的指标量P=0.96，表明测网重力点值变化的异常程度不高，与J1的指标量P=1.70所反映的测网整体存在较为突出的重力点值变化的结果截然相反。可见，经典平差的起算基准对重力网点值变化的显著性程度密切相关，不变的起算基准大幅度减小了定量参数化指标量的P值和R值，降低了其反映的重力网点值变化显著性程度，对分析重力场变化与地壳运动及地震趋势预报的关系^[17]造成一定的干扰。因此，为了真实客观地反映重力场变化的显著性程度，重力网优化改造时应继续完善与扩大绝对重力测量的时空监测能力，为重力网平差计算提供高精度的时变基准。

为分析非稳定的拟稳点对拟稳平差点值变化的影响，本文选取了测网中地质构造复杂、断裂活动较为活跃^[15]的10个非稳定点作为拟稳点，获得2015—2016年拟稳平差的重力网点值变化结果。非稳定点及断裂分布如表 3所示。

表  3  拟稳点及断裂分布^[15]

Table  3.  Distribution Quasi-Stable Points and Fractures

拟稳点	断裂分布		拟稳点	断裂分布
耿马	南宁河断裂（F1）		石屏	建水-石屏断裂（F6）
澜沧	龙陵-澜沧断裂（F2）		百色	百色-合浦断裂（F7）
勐海	打洛-景洪断裂（F3）		永川	华荣山断裂（F8）
普洱	镇远-普洱断裂（F4）		盘江	曲靖断裂（F9）
扬武	红河断裂（F5）		贵定	贵定-富泉断裂（F10）

下载: 导出CSV 

| 显示表格

图 11（a）为非稳定点拟稳平差的重力网点值变化结果，结合图9（a）、9（b）可知，利用非稳定点作为拟稳点进行拟稳平差的点值变化结果，与经典平差以及稳定点作为拟稳点的拟稳平差的点值变化结果均存在显著差异。测网西南部区域C的点值变化存在严重的扭曲和畸变，未能客观真实地反映区域重力场的变化趋势。图 11（b）给出了两种拟稳点下拟稳平差的重力点值变化差异情况，发现拟稳点选取不当将可能导致重力点值变化存在整体系统性偏差。

图  11  非稳定点拟稳平差的重力网点值变化及其与稳定点结果的差值

Figure  11.  Gravity Values Variation Calculated by Quasi-Stable Adjustment Method with Unstable Points and Difference Between It and that with Stable Points

下载: 全尺寸图片 幻灯片

图 12为利用非稳定点作为拟稳点进行拟稳平差的段差残差结果，由图 12可知，非稳定点拟稳平差的段差残差相比稳定点拟稳平差（见图 8（c））上下波动幅度大，对重力段差的还原不及稳定点的平差结果。此外，表 2中非稳定点N2指标量P值和R值明显小于稳定点N1，表明非稳定的拟稳点削弱了重力场点值变化的显著性异常程度。可见，拟稳点的选择直接关系到拟稳平差结果，对重力场变化的形态至关重要。拟稳点的选取要结合测网以往的重力点值变化资料和测点的区域构造背景，避开构造活动、断裂运动等地质背景较为复杂的测点，以及地下水开采及工业活动等易造成观测环境和重力值发生较大变化的区域^[18-19]。

图  12  非稳定点拟稳平差的段差残差

Figure  12.  Gravity Segment Difference Residuals Calculated by Quasi-Stable Adjustment Method of Unstable Points

下载: 全尺寸图片 幻灯片

3 结语

重力场的动态演化特征是地震预测和地壳运动等研究的重要信息来源，根据区域重力网的实际情况，选用合适的平差方法才能获取接近真实的重力场变化信息。本文基于南北地震带南段重力网2015—2016年的地震重力测量数据，有效地削弱偶然误差和系统误差的影响，对经典平差和拟稳平差开展深入的研究，得到以下结论：（1）提出了利用平差点值精度、重力段差与互差相关性、绝对点平差值与绝对值的差异结果来综合评定一次项系数优劣的方法。经典平差和拟稳平差的段差残差波动较小，对重力段差的恢复情况较好，标准差分析表明经典平差比拟稳平差能更真实地反映重力段差变化。（2）基于重力点值变化及其误差定义了两种定量参数化指标量P值和R值，为分析区域重力场点值变化的显著性程度提供定量依据。经典平差和拟稳平差的2015—2016年重力网点值整体变化具有较为近似的重力场变化形态和趋势，拟稳平差的指标量P=0.55（$<$1），远小于经典平差的P=1.70（> 1）。（3）经典平差起算基准对重力网的点值变化息息相关，并直接影响定量参数化指标量P值和R值。在点值变化误差相同的前提下，不变基准的指标量P=0.96（< 1）小于时变基准的指标量P=1.70（> 1），降低了其反映的重力场变化显著性程度。拟稳点决定了拟稳平差的重力场变化形态与真实重力场变化的近似程度，非稳定的拟稳点将导致重力点值变化发生扭曲和畸变。