地球重力场及其时间变化综合反映了地球表面和内部物质的分布状态、地球的形变以及质量迁移、各圈层之间物质交换和耦合。反映在不同时间尺度上，包括潮汐运动引起的重力短期变化，地表流体(包括大气、海洋和地下水等)引起的重力季节性变化以及构造运动、冰川均衡调整和地震活动等引起的重力长期变化^[1-3]。

地球重力场的时间变化可以通过地面连续重力观测、地面重复重力观测(包括绝对重力观测和相对重力观测)和卫星重力观测获得，区域重力网的时间变化主要通过地面重复重力观测获得。文献[4]对“中国地壳运动观测网络”项目1998—2008年共5期重复重力观测资料分别进行了静态平差解算，并利用各期静态平差结果分析了重力场的动态变化特征及其与地壳构造活动和大震活动的关系。目前，重力网平差方法主要采用静态平差方法，即假设观测期间各测点的重力值固定不变^[5-7]，未考虑观测过程中测点重力值随时间的变化。静态平差对观测区域小、观测时间跨度短、构造活动较弱的局部重力网是适用的，此时，单期观测中测点重力随时间的变化可以忽略。而在大区域特别是构造环境复杂，构造活动较强的区域进行重力网平差时，由于单期观测时间跨度大，观测期间测点重力值已发生变化，若仍采用传统静态平差方法，假设观测过程中重力未发生变化，必然会引入较大的系统误差。为了满足重力场时变研究的需要，亟须研究基于动态平差理论的数据处理方法。文献[8]利用重力随时间变化的线性速率模型对加拿大重力标准网超过40年的观测数据进行了动态平差，并为该地区冰川均衡调整模型的改进提供约束。文献[9]利用相同的平差模型对滇西重力网1986—2009年共52期流动重力观测数据进行了整体平差，得到了该地区重力场长期变化趋势。

我国的地面流动重力测量起步于20世纪60年代，经过近半个世纪的发展，流动重力测网已经覆盖了全国21个地震重点监视防御区中的18个，几乎实现了对全国主要地震监视防御区的全面监测^[10]。重力仪的观测精度从20世纪80年代中期的30×10^-8 m/s²提高到目前的10×10^-8 m/s²，最新的CG-5型重力仪观测精度可以达到6×10^-8 m/s^2[11]。随着“中国地壳运动观测网络”以及“中国大陆构造环境监测网络”项目的实施，覆盖整个中国大陆区域的重力网已经基本布设完成，其观测周期为2~3年，绝对重力测量精度不低于5×10^-8 m/s²，相对重力测量精度优于20×10^-8 m/s^2[12-13]，采用适当的数据处理方法可以得到高精度中国大陆重力场变化特征。然而，由于我国的全国重力网(“中国大陆构造环境监测网络”)覆盖整个中国大陆，具有覆盖面广、地震多发以及陆内构造活动复杂(不同地块构造活动差异较大)等特点^[14]，且该网络刚刚布设完成，观测资料较少，各期之间重力变化较大且呈现明显的非线性变化^[4]。因此，在全国重力网数据处理中应用静态平差模型和线性速率动态平差模型均会引入较大的系统误差。针对形变监测网动态平差存在的类似问题，文献[15]提出了形变网的分段线性动态平差法，文献[16]在分段线性模型中加入高程变化速率参数，提出了水准网分段速率整体平差模型。

考虑到我国全国重力网布设和观测的现状，为获得精确的中国大陆的重力场时变信息，本文研究了基于分段线性模型的数据处理方法，并应用于全国重力网的动态平差，为验证模型的可靠性，采用实测和模拟数据进行验证，该模型有利于获得真实的重力场变化信息。

1 动态平差模型

流动重力测量数据处理分为数据预处理和平差处理两部分。其中预处理包括仪器格值系数改正、潮汐改正、气压改正、仪器高改正和零漂改正等。

目前野外测量应用最多的是LCR-G型重力仪和CG-5重力仪。传统的平差模型^[6-8]将仪器在观测周期内的零漂率看作一个常数，并作为未知参数加入平差模型中。然而，由于LCR-G型重力仪零漂率是非线性的，且正负交替变化，不能将其看作常数；而CG-5重力仪零漂率在使用初期呈近似线性的持续增大态势，变化量约为20×10^-8 m·s^-2/h，之后会逐步下降、稳定^[17]，其零漂率也不宜在整个测量过程中看作一个常数。基于此，在平差模型中未将仪器零漂率作为未知参数，而是在数据预处理阶段按闭合观测时间段分段求取仪器零漂率，将经过预处理各项改正后的段差值作为平差的观测值。

1.1 流动重力数据预处理

由于受仪器结构以及野外观测所处环境噪声等的影响，流动重力观测数据中除偶然误差外还有系统误差，在平差前需要对数据进行预处理以消除或削弱这些系统误差的影响。根据CG-5重力仪的特殊结构以及各测站观测值需要的各项改正，测站上的重力观测值可用式(1)表示^[18]

式中，右边第1项为仪器格值系数转换，即由仪器读数转换为重力值；第2项为潮汐改正；第3项为气压改正；第4项为仪器高改正。其中，g_i为测站i重力仪的原始读数(对于LCR-G型重力仪需要先由仪器读数转换为相应的重力值)；E_k为仪器的格值系数；δ_th为重力潮汐因子；δf_c为永久性潮汐对重力的影响；G(t)为理论重力潮汐值；p为测点实测气压值；p_h为测点标准气压值；θ为重力垂直梯度；h为重力仪面板高度。

经过以上模型计算的各测点重力值，再按闭合观测时间段分段求取各测段的零漂改正即可得到观测段差值。零漂改正公式为^[18]

式中，；δg_d为仪器零漂改正值；Δt为观测段差终点与起点的观测时间差；d为仪器零漂率；g、g′、t、t′分别为测线起始点的往返观测值和观测时刻；g_i、g′_i、t_i、t′_i分别为需要计算静态漂移的观测点到达和离开的观测值和相应的观测时刻。

1.2 分段线性动态平差原理与模型

1.2.1 分段线性模型

由于重力场的时变特性，可以将重力网中各测点的重力值看作时间t的函数g(t)，由于重力场的时间变化反映了地壳运动、地震、火山喷发、大气海洋流动、冰雪融化、地下水变化以及地球内部构造变化等复杂的地球物理和地球动力学现象，因此，仅由有限次观测无法精确地求取函数g(t)的表达式。该模型假设重力随时间的变化率λ是分段线性的，即相邻两期观测之间重力变化率λ是线性的，则g(t)=g₀+λt。假设有m期观测，则可划分m-1个时段(图 1)，其中T_k(1≤ k≤m)为第k期观测的中心时刻，T₀为m期观测的中心时刻，λ^(k)(1≤k≤m-1)为第k期观测到第k+1期观测时间段的重力变化率，即可求取m-1个时段的重力变化率。

如图 1所示，假设T₃<T₀<T₄,则可以列出各期在观测点上的重力值计算公式为

1.2.2 分段线性模型的最小二乘平差

综合以上各式，可得到分段线性模型的误差方程

式中，V_ij为测点i、j之间的段差改正数； 和 分别为测点i、j在观测中心时刻T₀的平差重力值；g_i和g_j为测点i、j的实测重力值，其观测时间分别为t_i和t_j； 和 表示为如下向量形式

式中，(k=1,2,…,m-1)分别为测点i、j在T_k至T_k+1的重力变化率平差值；b_ik(k=1,2,…,m-1)、b_jk(k=1,2,…,m-1)分别为观测时刻t_i和t_j至T₀所经各观测时段的时间长度，未经过的时段元素为0，其取值与T₀的选择及测点观测时刻t_i和t_j有关。如图 1所示，当t_i<T₂、t_j<T₂、T₃<T₀<T₄，则

将测网中有绝对重力观测值的测点作为基准点，则根据每期的绝对重力观测值可列误差方程为

式中，为已知点i在中心时刻重力值的平差值；G_i为已知点i当期实测绝对重力值；V_i为平差改正数； 同式(4)。

综合式(4)、式(5)，误差方程式为

式中，V为改正数向量；A为系数矩阵； 为未知数向量；L为观测值向量，计算公式为

由于各仪器的观测精度不同，假设第k台仪器的观测中误差为σ_k，且其所有观测段差的精度相同，根据协方差传播律有

若采用n台仪器观测，且各台仪器的观测值独立，其协方差矩阵为

式中，D(g_A)为基准点绝对重力观测值的协方差阵。

设单位权中误差为σ₀，观测值协因数阵Q_LL= ，则观测值向量L的权阵为

根据最小二乘原理，改正数向量需满足条件

得到最小二乘解为

1.2.3 粗差剔除及平差模型的总体检验

由于仪器自身原因、测点的外界环境以及人为影响，在平差中所采用的多期重复重力观测资料不可避免会存在粗差。为了有效地消除或减弱粗差对参数估值的不良影响，本文采用抗差估计的选权迭代法——选取适当的权函数，通过不断的迭代，使粗差观测值的权为零(或接近于零)，从而达到抗差的目的。

由抗差最小二乘估计原理^[19-20]，误差方程式(6)的抗差解为

式中，为等价权矩阵，其矩阵内元素 定义为

式中，p_ij是权矩阵P的元素；w_i和w_j是以标准化残差 和 为自变量代入等价权函数所得结果，称为权因子。标准化残差为

式中，v_i、v_j是观测值残差；m_vi、m_vj是残差中误差。权因子的计算利用文献[21]提出的IGGⅢ模型，即

式中，v为标准化残差；k₀、k₁分别为保权临界值和零权临界值，其取值视具体情况而定。

相关抗差估计在实际计算时以迭代形式表示为

式中，是用第k次迭代的残差计算的等价权。

当 时，停止迭代，即为参数的抗差解，其中ε是收敛条件。

消除粗差的影响以后，由于平差过程所采用的数学模型及随机模型不正确，也会导致错误的平差结果。为验证平差模型的正确性，本文采用后验方差检验方法对采用的平差模型进行总体检验。

根据下列检验条件^[6]，若成立，则表示所采用的平差模型正确

式中，分别为先验中误差和后验中误差；m=n-t为平差的自由度；χ²(1-α;m)为当置信水平为1-α，自由度为m时的χ²分布临界值。

2 平差模型的模拟试验

分段线性动态平差模型可以在观测期数较少的情况下得到更真实的重力变化情况，既解决了传统静态平差模型未考虑测点重力随时间的变化的问题，又克服了线性动态平差模型参数过于理想化的问题。为验证该模型的可靠性，根据青藏高原与四川盆地的相互运动情况模拟了一组数据进行平差计算。由于受到印度板块的碰撞挤压，青藏高原东缘的巴颜喀拉块体向东南运动，受到四川盆地阻挡，造成块体内部在垂直方向的运动。据此，模拟了相邻两块体的运动情况，如图 2所示。假设block A地势较高，平均高程4000 m左右，重力值较小，并受到block C的碰撞而向block B运动，受到block B的阻挡导致地壳缓慢上升；block B地势较低，平均高程100 m，是一个稳定的克拉通陆块，阻挡了block A东向的继续运动，并且在整个过程中构造保持稳定，重力值较大，地壳缓慢下沉。block A和block B之间是一宽100 km左右的断裂带。

为了监测block A和block B的地壳垂直运动情况，在两块体一共布设了18个重力测点，每隔2~3年观测一次。具体点位及测线情况如图 3所示。

假设block A受到block C和block B的碰撞挤压而缓慢隆升，由于受到block C的碰撞，由西向东隆升速率逐渐减小，block A由剖面a、b分为3个区域(1)、(2)、(3)，假设各区域内地壳垂直运动情况及重力变化速率相近；block B中心区域沉降速率较大，周边区域相对较小，块体内部各测点垂直运动速率及重力变化速率基本为线性。参考Ekman等利用芬诺斯堪的亚地区1966—1984年流动重力和水准观测结果得出的重力-高程变化的比值-2.2×10^-8 m·s^-2/cm^[22]。模拟数据如表 1所示。

结合图 3中的测线情况，模拟了测网从1998—2008年共5期观测时间，见表 2。利用各期段差观测时间以及表 1中各测点的初始重力值和重力变化率，即可以得到各测段理论的观测段差数据。

为了模拟实际含有误差的数据，在理论段差值中加入随机误差。目前CG-5的实际精度可以达到6×10^-8 m/s^2[11]，即测点精度为6×10^-8 m/s²，依据误差传播定律，观测段差精度为8.5×10^-8 m/s²。野外实际测量时一般采用3台仪器同时观测，假设模拟数据采用3台仪器的平均值作为观测数据，根据误差传播定律，则模拟的段差精度为4.9×10^-8 m/s²。因此在理论段差数据中加入5×10^-8 m/s²的随机误差。

对测网采用上述模型进行平差时，选取重力段差最大的两个点作为已知点，分别是block A的01号点和block B的18号点，这样能够对测网中其他点做更好的约束。假设已知点的初始值和各期重力变化率均为已知，模型假设block A各区域内测点重力变化规律相近，block B内部各测点重力变化速率接近线性。在实际计算时，模型中相同的参数可以作为一个未知参数，这样以适当减少未知参数，获得更高的平差精度。

对上述试验进行整体动态平差解算，其平差结果的单位权中误差为2.9×10^-8 m/s²，重力变化率精度的平均值为1.3×10^-8 m·s^-2/a。为验证平差模型的可靠性，计算了该模型在95%的置信度的检验统计量——该试验模型自由度为m=169，置信度1-α=0.95，求得 =42.25、χ²(1-α;m)=200.34、χ_(m)²<χ²(1-α;m)，因此该模型符合统计检验模型，该动态平差模型有效。

为比较分段线性动态平差模型与传统静态平差模型的优劣，对5期模拟数据分别进行静态平差解算，各期平差结果均符合上述统计检验模型，可以进一步对两种平差结果进行分析比较。计算了两种平差方法得到的重力变化率与理论值的差值，其分布直方图分别见图 4和图 5。可以看出，动态平差结果中所得差值的分布符合正态分布，且80%以上的平差值与理论重力变化率的差值在1×10^-8 m·s^-2/a以内，只有两个差值超过2×10^-8 m·s^-2/a。而静态平差结果所得差值分布较散，平差值与理论值差值在1×10^-8 m·s^-2/a以内的占44.4%，差值超过2×10^-8 m·s^-2/a的占21%，相较动态平差结果误差较大。结果表明，该平差模型与传统的静态平差模型相比更能够反映实际重力场的变化情况，适用于时变重力场的研究。

3 全国重力网平差

我国全国重力网资料具有观测期数较少、覆盖区域构造环境复杂以及观测周期跨度较长的特点，与模拟数据相似，分段线性动态平差模型适用于全国网实测资料的平差处理。本文采用“中国地壳运动观测网络”项目1998—2005年共4期观测资料，并分别进行了传统静态平差和分段线性动态平差处理。

中国地壳运动观测网络(crustal movement observation network of China,CMONOC)是我国“九五”国家重大科学工程建设项目之一，地面重力观测包括23个绝对重力观测点以及400多个相对重力观测点。该项目绝对测量工作由中国科学院测量与地球物理研究所采用FG-5绝对重力仪观测，每个基准站上绝对重力测量精度优于5 μGal(1 Gal=0.01 m/s²)。相对重力联测工作由中国地震局、总参测绘局和国家测绘局联合完成，均由LCR-G型重力仪进行施测。所用4期观测资料详情见表 3。

利用静态平差方法对以上4期观测资料分别进行平差计算，不考虑测点重力值随时间的变化，利用各期绝对以及相对重力观测资料分别进行平差，其中第3期由于无绝对重力观测，其起算基准采用前一期观测时间较晚的7个绝对重力观测值，分别为长春、海拉尔、绥阳、广州、厦门、昆明、泸州，其观测时间为2001年10月8日—2002年1月29日，这些点与第3期相对重力观测时间比较接近。各期平差结果详见表 4。

利用分段线性动态平差方法对4期观测资料进行整体平差，首先利用表 3中各期相对重力观测时间中值T_i对重力变化率进行分段处理，结合总观测时间中值T₀及T₀时刻测点重力值即可列立误差方程，平差结果见表 4。

表 4中的点值精度代表了平差所求测点重力值的平均精度，对于静态平差而言，点值精度表示各期测点重力值的平均精度，对于整体平差而言，点值精度表示在总观测时间中值时刻测点重力值的平均精度。另外，整体平差得到的两期之间重力变化率的平均精度为10.1×10^-8 m·s^-2/a。结果表明，动态平差结果比静态平差精度更高。为进一步比较两种平差方法的优劣，计算了两种方法所得重力变化率的差值，其分布如图 6所示。

由图 6可以看出，两种平差结果最大差值为50×10^-8 m·s^-2/a，平均差值为13.4×10^-8 m·s^-2/a，其中大部分差值在20×10^-8 m·s^-2/a以内。造成这种差异的原因可能有以下几个方面：①单期观测时间跨度大，动态平差考虑了单期观测期间内测点重力随时间的变化，而静态平差未考虑此参数；②观测数据包括4期相对重力观测数据和3期绝对重力观测数据，其中2002年只进行了相对重力联测，没有进行绝对重力观测，所以静态平差时，第3期的基准点重力值只能采用观测时间相近的上一期少量绝对重力观测数据，此种情况下，采用静态平差对结果影响较大。然而由于测区构造环境复杂，很难获得真实的重力场变化特征，目前，仅由实测数据无法准确比较两种平差方法的优劣。

4 结语及讨论

本文介绍了重力网的分段线性动态平差模型，将测点各时段的重力变化率作为未知参数引入平差模型，并对重力网多期观测数据进行整体动态平差。该模型与传统的静态平差模型相比作了几点改进：①考虑仪器零漂率对平差结果的影响，在动态平差模型中不再将仪器零漂率看作常数，而是作为未知参数，在预处理阶段按闭合观测时段分段计算仪器零漂率，并将预处理后的段差值作为平差元素；②考虑了观测过程中测点重力随时间的变化，将测点的重力变化率作为未知参数加入平差模型中，解决了由于观测周期较长导致的测点重力变化对平差结果的影响；③将测点的重力变化率看作分段线性，有效减弱了由于观测时段内测点重力剧烈变化对平差结果的影响，动态平差模型更符合实际。同时，该平差模型也存在不足，由于引入重力变化率作为未知参数并联合多期观测资料整体平差，观测方程与平差参数均显著增多，法方程系数相比单期数据平差时更接近病态，平差计算困难，计算时间更长。

平差得到测点各时段内的重力变化率，为了满足地震监测、地壳垂直运动、地球内部物质迁移及冰川均衡调整等地球动力学研究的不同需要，可以对平差结果做不同的处理：①对各时段重力变化求和得到重力变化的累积结果；②分析不同时段重力变化率的变化情况；③对各时段重力变化率进行拟合求取测点长期的平均重力变化率。

为验证平差模型的有效性，本文分别利用分段线性动态平差模型和传统的静态平差模型对我国全国重力网数据和模拟数据进行了处理。结果表明，该平差模型相较传统静态平差模型可以更好地反映测点的实际重力变化情况，有利于真实重力场变化信息的获取，在目前流动重力观测仪器精度越来越高以及观测手段越来越多的情况下，采用适当的数据处理方法，可以进一步满足地球动力学研究的需要。