0 引言

一直以来，在使用自由落体式绝对重力仪进行绝对重力测量时，需要利用高精度相对重力仪测定重力垂直梯度.无法获取测站点重力垂直梯度时，只能采用重力归算中的空间改正系数-308.6 μGal·m^-1代替(李建国，2012；陈钊，2014；Xing et al., 2009；肖凡等，2013).实际上，绝对重力仪测观测数据，即自由落体运动的下落时间-距离序列，包含重力垂直梯度的信息.如果能够直接利用绝对重力仪观测数据确定重力垂直梯度，一方面可以减少绝对重力测量中的工作量和所必需的仪器设备，提高观测数据的使用率，降低测量成本；另一方面重力垂直梯度也可应用于反演地质构造形态、寻找油气田等领域，并且相比于重力异常具有更高的分辨力(Marson and Klingele, 1993；郇恒飞等，2012).

Niebauer等(1995)和Hanada等(1996)提出了FG5测量中可能存在的误差源，分析了观测数据中可能存在的有色噪声.Charles和Hipkin(1995)对实测数据进行分析，探测到了FG5观测数据中的有色噪声.随后众多学者研究了FG5观测数据中有色噪声的探测方法，包括Lomb-Scargle谱分析法(Orlob and Braun, 2012；张俊华等，2006)和高斯钟求和法(Orlob and Braun, 2013).同时，通过多台FG5绝对重力仪的比对观测(张为民等，2008)及与超导重力仪的观测结果进行比较(Van Camp et al., 2005)，验证了不同绝对重力仪之间存在系统偏差.Hipkin(1999)首次利用自由落体式绝对重力仪14400次下落测量数据解算了测站的重力垂直梯度，精度达到了3 μGal·m^-1，并且与相对重力仪的重力垂直梯度观测结果一致.然而，Hipkin(1999)在解算重力垂直梯度时没有考虑观测数据中的有色噪声，解算方法的适用性不强.所以至今利用绝对重力测量数据解算重力垂直梯度仍然处于实验阶段，而未能形成稳定有效的数据处理方法.

本文首先建立利用FG5实测数据求解重力垂直梯度的数据处理模型与算法；然后分析FG5绝对重力仪观测噪声(有色噪声和高斯噪声)和下落时段的选取对解算重力垂直梯度的影响，在此基础上对有色噪声进行探测和建模，并量化下落时段的选取依据；最后分别对FG5-214绝对重力仪在两个测站上的观测数据进行处理和分析，验证本文方法的正确性和有效性.

1 数据处理模型与算法

自由落体式绝对重力仪的观测数据是在质体下落过程中采集得到的下落距离-时间序列(x_i-t_i).由此解算重力加速度和重力垂直梯度的基本观测方程为(Microg-LaCoste，2012)：

(1)

(2)

其中，x₀和v₀分别是t=0时刻下落质体的位置和速度，g₀是x=0处的重力加速度，γ为测站处的重力垂直梯度.是时间观测值经过光传播延迟改正后(公式(2))得到的下落时间.在传统的绝对重力测量中，γ由相对重力仪测定.本文试图将γ作为未知参数与x₀、v₀和g₀同时求解.

FG5绝对重力仪使用了基于相位调制的激光干涉测距仪，采集的下落距离观测量中包含了频率为激光调制频率Ω的正弦项(或余弦项)，因此，观测方程(1)改写为(Orlob and Braun, 2012)：

(3)

当重力垂直梯度作为未知参数时，观测方程不再是未知参数的线性方程，需要采用迭代计算方法，具体的迭代计算步骤可以参考Hipkin(1999).本文依据该方法处理FG5-214的观测数据，但所求解的重力垂直梯度与相对重力仪的测量结果相差超过50 μGal·m^-1，未能达到预期精度.为此，本文提出一种改进的解算方法，算法流程如图 1所示.

对于第m次自由落体实验，由公式(4)求出下落距离拟合残差序列{v_i}_m.由n次下落实验的距离拟合残差序列叠加平均得到残差均值序列{v_i}(公式(5)).对残差均值序列进行分析，建立有色噪声模型，如公式(6)所示，采用有色噪声模型对残差均值序列进行拟合逼近求出模型参数，得到有色噪声序列{v_{i, s}}.由残差均值序列减去有色噪声序列得到剩余残差序列{v_{i, r}}(公式(7)).依据{v_{i, r}}选取可靠的下落时段内的观测数据，减去有色噪声{v_{i, s}}，再采用Hipkin(1999)的方法重新解算.统计多次下落解算结果的均值和标准差，得到重力垂直梯度最终结果和精度.

(4)

(5)

(6)

(7)

公式(4)中、、、、和分别为采用Hipkin(1999)方法解算的下落初始位置、速度、加速度、重力垂直梯度和周期项振幅参数的估值.采用该方法同时解算重力和重力垂直梯度时，依据期望达到的精度，重力和重力垂直梯度的迭代逼近阈值分别选为0.1 μGal和1 μGal·m^-1.

依据对实测数据的分析，认为有色噪声模型包括指数衰减项和周期项.采用公式(6)拟合逼近残差均值序列{v_i}求解有色噪声模型参数时，先将{v_i}分为数段(7~10段)，选取每一段中的残差幅值最大值对指数衰减项进行拟合，得到公式(6)中的参数k.{v_i}扣除衰减项的影响后，主要包含周期项误差，通过对其插值形成等时间间隔残差序列，采用傅里叶变换探测周期项频率f_j，最后通过线性拟合求解周期项参数A_j和D_j(j=1, 2, …, N).

一般认为下落开始和结束阶段的观测数据中噪声较大，前人在解算时凭借经验剔除下落开始和结束阶段的观测数据确定下落时段.本文依据剩余残差{v_{i, r}}选取下落时段，其准则为：选取的下落时段要尽可能长，以保留足够多的观测数据；同时选取的下落时段内剩余残差的幅值不得超过0.1 nm(具体依据见2.1节).

2 实验结果与分析 2.1 观测噪声的影响

采用模拟实验研究有色噪声和高斯噪声对解算结果的影响.首先给定一套参数，设计出一组不含噪声的观测数据，然后分别加入高斯噪声和余弦噪声计算重力垂直梯度.对于高斯噪声，研究重力垂直梯度随着噪声标准差的变化，结果如图 2a所示；对于余弦噪声，研究重力垂直梯度随着余弦项振幅和频率的变化，结果如图 2b所示.

图 2a和图 2b中的直线表示289±200 μGal·m^-1，是重力垂直梯度值的参考范围.由图 2a可知，只有当高斯噪声的标准差在0.1 nm以内时，其对结果的影响才能忽略不计；标准差超过1 nm时，结果的大小就有可能超出合理范围.由图 2b可知，余弦噪声对结果的影响随着振幅的减小而降低，随着频率的增加而减小.然而，即便是振幅为0.1 nm、频率为150 Hz的余弦噪声对结果的影响也不能忽略.因此，剩余残差超过0.1 nm的下落时段对应的观测数据应当舍弃，同时有必要对观测数据中的有色噪声进行探测和建模.

2.2 下落时段选取的影响

利用FG-214绝对重力仪在A测站上采集了2400次下落实验的观测数据，每次下落实验可以获得700对距离-时间观测数据.设计每次实验的时段起点在1~30号观测数据中选取，时段终点在601~700号观测数据中选取，总共3000种下落时段选取方案.每一种方案可以解算得到2400个重力垂直梯度，剔除明显受到粗差影响的结果后取其均值，3000个均值结果如图 3所示.

图 3中的编号表示不同的下落时段选取方案，直线表示相对重力仪测定的重力垂直梯度(289.0 μGal·m^-1).选取n号数据为起算点，m号数据为下落终点，对应的结果编号则为(n－1)×100+m.图 3显示，当选取不同下落时段时，解算得到的重力垂直梯度在180~310 μGal·m^-1之间变化.由此可知，下落时段的选取对重力垂直梯度的解算结果影响很大，有必要对下落时段的选取给出合理的量化依据.

2.3 实测数据处理结果

由FG5-214绝对重力仪在A、B两个测站上采集的两组观测数据作为实验数据，依据本文提出的数据处理方法解算重力垂直梯度，分别与使用相对重力仪在两个测站上测定的重力垂直梯度进行比较.

以处理A测站上的观测数据为例，对A测站上2400次下落数据独立解算，得到2400个距离拟合残差序列，通过叠加求均值获取残差均值序列，如图 4所示.

由图 4可知，2400次下落的距离拟合残差均值序列的幅值超过了0.2 nm，相比于单次下落的距离拟合残差(幅值不超过3 nm)只减小了一个量级，显示距离拟合残差中含有有色噪声.

使用残差均值序列对公式(6)进行拟合，得到有色噪声模型：

(8)

有色噪声建模过程及结果如图 5所示.图 5a中上下两条曲线为衰减项的拟合结果，基本符合残差均值序列的衰减趋势.图 5b为残差均值序列剔除衰减项的影响后所得残差.图 5b中的残差序列经过插值得到等时间间隔序列，然后进行傅里叶变换得到频谱图，如图 5c所示.获取周期项频率后，通过线性拟合计算出周期项系数，最终构建了有色噪声模型(公式(8))，结果如图 5d所示.

从残差均值序列中剔除有色噪声得到剩余残差序列，如图 6所示.图 6中，两条横线为±0.1 nm的阈值线，阈值线与残差序列有多个交点，下落时段的起点和终点将在这些交点中选择.

剩余残差超出阈值线的下落时段对应的观测数据含有粗差，连续的粗差观测数据将引入低频噪声，对结果影响较大；而少量离散的粗差观测数据的存在对解算结果基本没有影响.因此，依据本文提出的时段选取准则，选取图 6中两条竖线之间的观测数据参与解算.

采用选取的下落时段内的观测数据，扣除有色噪声，解算得到2400个重力垂直梯度.统计重力垂直梯度结果的均值和中误差，并以3倍中误差作为阈值，进一步筛选结果，最终认为可靠的结果如图 7所示.

采用同样的方法处理FG5-214在B测站的观测数据.整理A和B两个测站的数据处理结果，如表 1所示.

表 1中的参考值是利用相对重力仪测定的测站重力垂直梯度，精度优于3 μGal·m^-1.表 1中的原有方法是未考虑观测数据中有色噪声并依据经验选取下落时段的解算方法，其解算结果与参考值的差值超过了50 μGal·m^-1.采用本文提出的改进方法处理A、B测站的观测数据，解算的重力垂直梯度与参考值分别相差7.82 μGal·m^-1和2.48 μGal·m^-1.相比于原有方法的处理结果，本文解算结果的精度没有明显变化，但是均值更接近于参考值.此外，测站B的解算结果明显优于测站A的结果，这主要是因为前者的观测条件和数据质量均好于后者所致.

在A、B两个测站上，统计可得单次下落数据解算的重力垂直梯度的精度分别为376.00 μGal·m^-1和216.84 μGal·m^-1.Hipkin(1999)统计的单次下落重力垂直梯度的标准差为354.9 μGal·m^-1，与本文结果基本一致.Hipkin(1999)采用了14400次下落实验的数据进行解算，统计得到重力垂直梯度的均值精度为3.0 μGal·m^-1.如果采用同样数据量的观测数据依据本文的方法进行解算，所得重力垂直梯度均值的统计精度应该可以达到3.0 μGal·m^-1.

3 结论与讨论

本文利用绝对重力测量数据解算重力垂直梯度，提出了一套新的数据处理方法，主要优势在于探测出了观测数据中的有色噪声并成功建模，同时对参与解算的下落时段的选取初步给出了量化的依据.以相对重力仪的观测结果为参考标准对解算的重力垂直梯度进行评价，结果表明本文方法的处理结果相比于原有方法的处理结果更接近于参考值，说明本文提出的数据处理模型和方法相比于原有方法有较大的改进，具有一定的参考价值.同时，利用本文的方法处理了FG5-232绝对重力仪观测数据，所得重力垂直梯度结果与相对重力仪观测结果的偏差达到了40 μGal·m^-1，表明本文方法不适用于处理FG5-232绝对重力仪观测数据，原因有待查明.另外，由本文的研究结果还可以总结出以下结论：

(1) 自由落体式绝对重力仪在下落开始和结束阶段的观测数据不稳定且可能连续出现粗差，因此下落时段的选取对解算高精度重力垂直梯度影响较大.本文通过对多次自由落体实验的下落距离拟合残差叠加求均值，不仅可以探测观测数据中的粗差和有色噪声，而且为下落时段的选取提供了量化的依据.

(2) 要保证重力垂直梯度的均值结果中不存在较大偏差，观测噪声中有色噪声必须予以剔除或削弱，同时随机噪声的幅值应该控制在0.1 nm以内.由于重力仪之间的硬件差异可能导致观测噪声的特性不同，所以在处理其他重力仪观测数据时应当具体分析，本文的数据处理方法可供参考.

(3) 绝对重力测量中重力垂直梯度的测量精度应当优于3 μGal·m^-1.由本文的实验结果可知，单次下落实验的重力垂直梯度解算结果的标准差大约为300 μGal·m^-1.要达到3 μGal·m^-1的统计精度，需要在一个测站上进行10000次下落实验.然而，绝对重力测量一般测定25组有效数据，每组观测100次下落，实施10000次以上的下落实验会导致测量效率大大降低，所以本文确定重力垂直梯度的方法还无法满足工程化应用.

致谢

感谢国家测绘地理信息局第一大地测量队提供了FG5绝对重力仪观测数据.感谢许闯博士在数据处理方面提供了宝贵意见.