利用航空重力测量技术可获得高精度的中高频地球重力场信息，研究表明，载体运动加速度的测量精度是制约航空重力测量技术发展的主要障碍之一^[1-3]。当前航空重力测量系统普遍采用动态差分GPS技术（differential GPS，DGPS）获得载体位置序列，对位置序列分别进行一次和二次微分得到载体速度和运动加速度^{[1, 4]}。然而该方法对位置的精度要求较高（cm级）^[1-2]，在动态长基线的情况下模糊度往往难以较好地固定，进而影响定位精度，因此需要在测区布设一定数量的基准站，以确保测线上所有点位的解算精度。文献[2]利用双差相位观测值构造了多谱勒观测值及其变化率，并用该观测值直接求解载体加速度，该方法消除了模糊度参数，因而大大降低了对位置的精度要求，但由于实验仅采用了4颗GPS卫星参与解算，因此最终效果并不理想。之后，一些学者对该方法做了进一步完善，文献[3]在计算观测值间协方差模型时考虑了电离层和对流层的影响，而文献[5-6]为提高加速度解算的鲁棒性和精度，采用了多基站的数据处理模式。此外，文献[1]进一步比较了位置差分法与相位差分法求解的载体速度与加速度，并指出当载体处于低动态时，位置差分法较相位差分法解算的载体运动加速度精度更高。

上述方法均采用了DGPS技术，这不仅增加了人力、物力的开销，而且在高原、山区、荒漠等困难测区布设基站也存在较大难度。为此，文献[7]尝试通过对精密单点定位(precise point positioning，PPP)技术获得的载体位置二次微分求解载体加速度，并联合惯性测量单元解算重力扰动。其实验结果表明，去除异常值后，重力扰动交叉点（共24个）不符值的均方根(root mean square，RMS)为3.3 mGal (1 mGal=1×10^-5 m/s²)，并指出利用PPP技术解算载体位置需要较长的收敛时间。文献[8]同样采用PPP浮点解解算的加速度获得了约1.0 mGal的网平差后重力扰动估计精度。近年来，文献[9-10]提出了利用单台GPS接收机通过历元间相位观测值之差求解高精度载体速度的方法VADASE，该方法在静态环境下可获得cm/s甚至mm/s级的测速精度。文献[11]从数学模型上证明了VADASE方法与PPP具有等价性。然而相较于PPP，该方法消除模糊度的同时可直接估计载体速度，通过一次差分便可获得加速度信息。考虑到单站标量航空重力测量具有很高的实用性和市场价值，目前国际上还鲜有研究，本文首先将VADASE方法扩展至航空重力测量应用中，着重分析卫星轨道、卫星钟差产品对加速度估计的影响，然后联合GT-2A航空重力测量系统的实测数据求解测线重力扰动，最后分别采用交叉点不符值以及地面实测重力值评定其内外符合精度。

1 单站标量航空重力测量基本原理 1.1 标量航空重力测量数学模型

以标量航空重力测量为例，基本观测方程为^[12]：

${\rm{ \mathsf{ δ} }}g = \dot v - \left( {f - {f_0} + {g_b}} \right) + {\rm{\Delta }}{\alpha _E} - \gamma $

(1)

式中，δg为重力扰动；为GPS测得的重力仪参考中心的垂向运动加速度；f为垂直加速度计测得的比力；f₀为飞机静止时在停机坪处垂直加速度计测得的比力；g_b为停机坪处的绝对重力值；γ为正常重力；Δα_E为各项误差改正，包括厄特弗斯改正、水平加速度改正、自由空间改正等，具体计算可参考文献[13-14]。

为使重力扰动测量精度优于1.0 mGal，载体定位精度需优于3.0 m，测速精度为cm/s级。载体速度可通过VADASE获得，其精度可达cm/s甚至mm/s级。通过对速度积分获得高精度的坐标位移，加上已知的坐标初值，便可获得载体动态位置信息，其定位误差为0.1 m/h^[15]。因此，由式（1）可见，影响重力扰动解算精度的一项主要误差来自GPS获得的载体运动加速度。

1.2 VADASE求解加速度

基于无电离层组合的历元间相位观测值之差可表示为^[10]：

${\rm{\Delta }}\varphi _r^s = {\rm{\Delta }}\rho _r^s + {\rm{\Delta }}{t_r} - {\rm{\Delta }}{t^s} + {\rm{\Delta }}T_r^s + {\rm{\Delta }}p_r^s + {\rm{\Delta }}E_{r, j}^s$

(2)

式中，Δ为历元间单差算子；φ_r^s为无电离层组合相位观测值；ρ_r^s为站星间几何距离；t^s和t_r分别为卫星和接收机钟差；E_{r, j}^s为卫星轨道钟差误差、多路径误差以及观测噪声之和。其中几何距离变化量Δρ_r^s由3部分构成：

${\rm{\Delta }}\rho _r^s = \left[ {{\rm{\Delta }}\rho _r^s{]_{{\rm{OR}}}} + } \right[{\rm{\Delta }}\rho _r^s{]_{{\rm{Etol}}}} + e_r^s{\rm{\Delta }}x$

(3)

式中，[Δρ_r^s]_OR表示由卫星轨道以及地球自转引起的变化量；[Δρ_r^s]_Etol为地球潮汐影响，包括固体潮和海潮；Δx为载体历元间位移变化量，即速度；e_r^s为站星间方向余弦。

由于航空重力测量不受潮汐影响，因此式（3）可简化为：

$ {\rm{\Delta }}\rho _r^s = {[{\rm{\Delta }}\rho _r^s]_{{\rm{OR}}}} + e_r^s{\rm{\Delta }}x $

(4)

将式（4）代入式（2）可得：

$ {\rm{\Delta }}\varphi _r^s = \left( {e_r^s{\rm{\Delta }}x + {\rm{\Delta }}{t_r}} \right) + ({[{\rm{\Delta }}\rho _r^s]_{{\rm{OR}}}} - {\rm{\Delta }}{t^s} + {\rm{\Delta }}T_r^s + {\rm{\Delta }}p_r^s) + {\rm{\Delta }}E_{r, j}^s $

(5)

式中，([Δρ_r^s]_OR-Δt^s+ΔT_r^s+Δp_r^s)可用已知模型以及卫星轨道和钟差产品计算获得。则待估参数Δx和Δt_r可用最小二乘法计算获得。

利用高精度的速度，可将载体位置表示为：

$ {P_{i + 1}} = {P_i} + {\rm{\Delta }}x $

(6)

其中，i+1时刻的位置P_i+1可由i时刻的位置P_i加上历元间位移变化量表示。载体初始坐标可由起飞前一段时间的静态精密单点定位获得。

对速度作一次差分，便可获得加速度：

$ {\rm{\Delta }}{{\dot x}_i} = \sum\limits_{j = - M}^M {{h_M}\left( j \right){\rm{\Delta }}{x_{i - j}}} $

(7)

此处选用一阶泰勒多项式进行差分，因此差分阶数M设为1，差分系数为h=[0.5 0-0.5]^{T[1, 15]}。

2 GT-2A航空重力测量系统飞行实验

本文实验选用2015年4月25日至4月26日中国陕西秦岭测区两个架次的航空重力测量数据。第一架次由东西方向共8条检测线构成，编号001~008，每条测线长约110 km，测线间距约为5 km。第二架次由南北方向共11条主测线构成，编号009~019，每条测线长约60 km，测线间距约为3 km。为减少气流影响，实测时间均为当地时间04:00:00至10:00:00左右，平均飞行高度约为4 000 m，飞行速度约为250 km/h。同时在机场布设一个基站，距测区约100 km，用于解算DGPS的结果作为参考值。此外，测区内共有21个地面重力观测点，用以评估重力扰动的外符合精度。

采用GT-2A航空重力测量系统采集重力数据^[16]时，在飞机顶部中轴线上安置航空天线A，与其相连的为一台ProFlex800 Ashtech接收机，如图 1所示。GNSS接收机可接收GPS L₁/L₂频率的观测数据，数据采样率为2 Hz。采用FIR低通滤波器对GNSS加速度滤波，本文采用滤波长度T为100 s^[17]，对应空间分辨率R为3.5 km（R=TV/2，V为载体速度），滤波截止频率设为0.01 Hz。之后联合重力数据、姿态数据解算重力扰动。

3 误差分析

对流层和电离层二次差误差在短时（0.5 s）内认为可忽略不计，重力信号主要集中在低频段0~0.01 Hz内^{[1, 18]}，因此本节采用傅里叶变换，重点分析卫星轨道和卫星钟差对加速度估计的影响。

3.1 卫星轨道影响

目前，国际GPS服务(international GPS service，IGS)最终精密轨道产品精度可达2~3 cm^[19-22]，本文将该产品计算的加速度作为参考值，分析广播轨道的广播星历(broadcast ephemeris，BRDC)和超快速轨道的快速预报星历(IGS ultra rapid products，IGU)（观测）对加速度估计的影响。为保证结果不受卫星钟差结果影响，实验中统一采用IGS 5 min采样间隔的精密钟差产品。采用机场基准站8 h的静态观测数据通过静态模拟动态处理，以评估卫星轨道对加速度的影响。广播轨道和超快速轨道计算结果与参考值之差如图 2所示。

由图 2可见，3种卫星轨道产品之间的差异在所示频段内对加速度估计影响相对较小，其振幅基本小于0.1 mGal。此外，IGU和BRDC频谱基本重合，差异并不显著。主要原因是，卫星轨道变化连续光滑，即便是广播轨道产品，其定轨精度约为1 m，该项误差也可通过短时历元间差分予以消除。将BRDC和IGU计算结果互差，其互差结果的均值为0 mGal，RMS为0.1 mGal。再次说明，卫星轨道对加速度估计的影响可以忽略。

3.2 卫星钟差影响

相较于卫星轨道，卫星钟差的稳定性要差。因此，本文设计了5组实验方案，如表 1所示，分析钟差采样率以及不同机构提供的钟差产品对加速度估计的影响。IGS目前只发布5 min和30 s的卫星钟差产品，欧洲定轨中心(Center for Orbit Determination in European，CODE)可提供5 s钟差产品^[23]。为保证产品的一致性，每组方案的精密产品均由同一机构提供。由于测站静止，因此认为加速度“真值”为0。

图 3给出了5种方案加速度解算结果，并用不同颜色区分。由图 3可见，方案5计算的加速度在所示频段内振幅最小；而方案1和方案2计算结果最差，即使在低频段0~0.01 Hz内，卫星钟差误差也有明显的振幅；方案3和方案4的加速度频谱基本相当，在0~0.01 Hz内振幅仅次于方案5，但在0.01 Hz之后，振幅快速增大。

图 4中加速度均已经过有限冲激响应滤波器(finite impulse response，FIR)的低通滤波，滤波截止频率为0.01 Hz。时域上同样可见，方案1和方案2结果最差，两者的RMS为9.53 mGal；方案3和方案4结果次之，其RMS分别为6.02 mGal和6.17 mGal；而方案5结果最优，RMS为0.92 mGal；5种方案加速度解算结果的均值均为0。这说明卫星钟差精度与采样率对加速度估计具有很大的影响。其主要原因为高采样率的卫星钟差能捕获到卫星钟差的高频信号。

4 载体加速度计算实例

利用本文所提方法计算载体加速度，并结合重力数据与姿态数据解算测线上的重力扰动，以GT-2A航空重力测量系统自带的商业软件解算的空间重力扰动作为参考值（标称精度为0.6 mGal）^[24]，结合交叉点重力扰动不符值评定内符合精度^[25]。最后，采用泊松积分法将本文方法以及DGPS方法计算的空间重力扰动向下延拓至地面，其中积分半径为0.5°，远区贡献采用EGM2008重力场模型计算，将延拓值与地面实测重力值对比以评定外符合精度。

分别采用第3节中的5种实验方案解算重力扰动，并与参考值互差。图 5为第一条测线的互差结果，并用不同颜色的线条区分。其中BRDC和IGS 5 min解算结果最差，与参考值互差基本在5 ~ 10 mGal左右，IGS 30 s和CODE 30 s结果基本一致，而采用CODE 5 s的精密钟差结果最优。

图 6给出了不同实验方案下解算的重力扰动与参考值互差的统计结果。统计结果显示，所有互差结果具有0均值的特性。从RMS结果来看，方案1和方案2的RMS要明显大于其他方案，基本在5~10 mGal之间；方案3和方案4的RMS次之；而方案5的结果最优，其RMS基本小于1 mGal。此外，本文还采用无电离层组合观测值对方案5重新解算，其与L₁解算结果的互差均值为-0.06 mGal，RMS约为0.45 mGal，说明电离层的影响相对较小。

对100 s滤波长度下的测线重力扰动进行网平差后结果如图 7所示。图 7中横、纵坐标分别为相对经、纬度，不同颜色和大小的圆圈表示不符值绝对值的大小，共88个。可见，绝大部分交叉点不符值的绝对值处于0~2 mGal之间，平差后解得交叉点不符值的RMS为1.13 mGal。

采用傅里叶变换法，将本文方法以及DGPS方法计算的空间重力扰动延拓至地面^[25]，其与地面重力实测值之差如图 8所示。由图 8可见，延拓值与地面重力实测值之差90%在±5 mGal以内，个别接近6 mGal。统计差值的RMS，对本文方法和DGPS法而言，分别为3.10 mGal和2.99 mGal。产生与地面实测重力值差异的主要原因有：一方面来自参数估计和仪器的量测误差；另一方面，向下延拓会放大噪声。上述实验结果说明，采用本文方法估计航空重力载体加速度是可行的。

5 结语

针对传统航空重力测量作业方式存在的弊端，本文将VADASE方法扩展至航空重力测量中，用于单站求解载体加速度。首先介绍了VADASE方法的基本原理，接着重点分析了卫星轨道和钟差对加速度估计的影响。结果表明，卫星轨道对加速度估计的影响很小，可忽略不计；相同采样率不同机构提供的钟差产品对加速度估计的影响差异不明显，而高采样率的卫星钟差可减少钟差误差的影响。其主要原因是，卫星钟差较轨道而言，稳定性较差，高采样率能捕获到钟差的高频信号。然而，本文方法受精密钟差产品时延的制约，相较于DGPS，仍无法及时解算载体加速度。

采用CODE 5 s精密钟差解算载体加速度，联合重力数据和姿态数据解算测线重力扰动，并与DGPS解算的参考值对比，两者差值优于1 mGal。网平差后重力扰动交叉点不符值的RMS为1.13 mGal。将本文方法和DGPS法解算的重力扰动向下延拓至地面，并与地面实测重力参考值对比，统计差值的RMS，分别为3.10 mGal和2.99 mGal。上述实验结果表明，单站标量航空重力测量的方法是可行的。