经济高效的航空重力测量能够覆盖地面重力测量难以到达地区，可提供精确的地球重力场中、短波信息，已经成为大规模重力数据采集的主要手段。大地测量和地球物理勘探等领域需要高精度的航空重力测量数据，对航空重力空中测线数据质量进行客观评价是保证航空重力数据应用水平的重要环节。

航空重力数据的内符合精度一般采用重复测线对比和交叉点不符值分析来评价^[1-5]，外符合精度则需要利用其他独立观测值进行评定^[6-8]。一种评价航空重力数据外符合精度的方案是将航空重力观测值向下延拓至地面或大地水准面上，与地面重力观测进行对比检核^[9-16]。该方案的问题在于向下延拓是一个数据噪声得到放大的不适定过程，无论采用何种正则化解法或平滑方法都只是在一定程度上削弱或控制向下延拓误差，与地面重力观测的对比结果必然包含向下延拓误差，其量级与飞行高度和观测噪声有关，难以真实反映航空重力空中测线观测值的精度水平。

另外一种方法是将地面重力数据向上延拓至空中，与航空重力数据进行对比。不同于向下延拓，向上延拓是非常稳定的，在没有积分边界效应的情况下，即数据范围相对于计算范围有足够的冗余，向上延拓误差主要来源于地面重力格网化插值误差和代表性误差。此外，由于重力场谱组合理论和谱权确定方法取得的最新进展^[17-18]，直接利用空中重力数据，不经过显式的向下延拓，已能够稳定解算区域重力场和大地水准面。因此，如何对航空重力空中测线数据外符合精度进行客观评价值得研究。

本文介绍3种基于地面重力向上延拓的航空重力测线数据外符合精度评估方法：基于梯度、基于泊松积分和基于快速傅里叶变换(FFT)的地面重力向上延拓。基于这3种方法，提出交叉检验方法估计地面重力数据误差及其空中误差传播，对毛乌素测区的航空重力空中测线数据进行外符合精度评价。

1 基于梯度的地面重力向上延拓

若在航空重力测区布设有密集的地面重力点，由于地面重力观测精度高，可利用地面重力数据对航空重力空中测线观测值的精度进行验证。一种方法是将地面重力扰动向上延拓至空中，与实测空中测线重力扰动进行比较。航空重力测量中采用差分GNSS动态定位方式获取测线采样点处的大地高，因此本文中的测线重力观测值均为重力扰动，而不是重力异常。

地面重力扰动向上延拓可采用梯度法，计算地面重力扰动沿椭球法线方向上的梯度，按照点对点延拓方式将地面重力扰动向上延拓至空中。空中与地面重力扰动满足如下Taylor级数展开式

式中，δg^Ter为地面重力扰动观测值；h为地面点的大地高；δg^TerU为地面重力扰动向上延拓值；h′为空中点的大地高；为地面重力扰动梯度，可采用下式由地面重力扰动计算^[19]

式中，l₀为计算点到流动点的空间距离。式(2)属于积分计算，应先将地面重力扰动格网化后再代入计算。

地面重力点数量一般远小于航空重力测线点数量，因此需要将航空重力测线观测值内插到地面重力的空中延拓点处，再对二者做差

2 基于泊松积分的地面重力向上延拓

除了基于梯度的点对点向上延拓外，还可利用泊松积分将地面重力数据向上延拓至空中测线处^[19]。理论上泊松积分需要输入全球均匀分布的数据，但实际中地面重力数据仅覆盖局部区域，因此一般需采用移去恢复法进行计算，选定某个地球重力场模型作为参考重力场，地面重力数据覆盖范围以外由重力场模型贡献，高频地形影响采用残余地形模型(RTM)方法计算^[20]。移去恢复法的另外一个优势在于，移去参考重力场和高频地形影响后的残余地面重力扰动量级很小，在格网化插值、泊松积分计算时引入的误差也很小。

移去参考重力场和高频地形影响的残余地面重力扰动δg_Res^Ter写为

式中，δg_EGM^Ter为由参考重力场模型计算的地面重力扰动；δg_RTM^Ter为地面点的高频地形影响，由下式计算^[20]

式中，ρ为平均地形密度；(x_P, y_P, z_P)和(x, y, z)分别为计算点与流动点的局部直角坐标；H_P为计算点的真实地形高度；H_Ref为计算点的平均地形高度。

图 1描述了基于泊松积分的地面重力数据向上延拓原理。地面重力数据分布在具有地形起伏的地球表面，需将地面重力扰动解析延拓归算到当地平均地形高度面上。首先利用式(2)由残余地面重力扰动计算得到梯度，再根据地面重力点与平均地形高度面之间的高度差，类似于式(1)计算解析延拓改正量，将残余地面重力扰动统一归算到当地平均地形高度面。

在球近似条件下，将平均地形高度面上的残余地面重力扰动δg_Res^TerM向上延拓得到空中残余重力扰动δg_Res^TerU的泊松积分表示为

式中，R_E为地球平均半径；h_M为地面点的平均大地高；h′为空中点的大地高；σ表示单位球积分区域；dσ=cos ϕ′dϕ′dλ′为单位球积分面元；l为球面σ上的流动点到空中计算点的距离。

残余地面重力扰动向上延拓到空中后，还需要在空中点位处恢复参考重力场模型和高频地形影响的贡献，得到最终的地面重力扰动向上延拓值δg^TerU

式中，δg_EGM^air是由参考重力场模型计算的空中点位处的重力扰动；δg_RTM^air为空中点位处的高频地形影响。

3 基于FFT的地面重力向上延拓

在平面直角坐标系中，泊松积分公式写为

式中, (x, y)表示平面坐标；Δh为向上延拓高度，即平均地形高度面到平均飞行高度面的距离。式(10)可写成空域内卷积的形式^[21]

式中，*为卷积号；w(x, y)为泊松核

利用快速傅里叶变换(FFT)将空域内的卷积化为频域内的乘积

式中，F₂表示二维傅里叶变换；W(u, v)为泊松核在频域内的表达式

式中，u、v分别表示x、y方向上单位距离内的周期数；表示总频率；i²=－1。

因此在频域内重力扰动向上延拓表达为

式中，F₂^－1表示逆二维傅里叶变换。FFT方法本质上是泊松积分的平面近似，其优点在于计算效率高，对于大规模重力数据处理更为有利。

4 地面重力误差估计及其空中传播

在飞行高度处，航空重力扰动与地面重力扰动向上延拓值之间的差值由航空重力数据误差和地面重力向上延拓值误差两部分构成。地面重力向上延拓值误差主要来源于：①单点地面重力测量误差；②地面重力解析延拓误差，即将地球表面的重力观测值归算到平均地形高度面引起的误差；③地面重力格网插值误差，即对地面重力数据进行格网化产生的插值误差；④地面重力代表性误差，即用格网平均值代替离散点值引起的代表性误差^[22]；⑤向上延拓方法本身引起的误差。此外，若数据范围有限，靠近数据范围边界的计算点还会受到边界效应的影响。

单点地面重力测量误差为10~20 μGal，可以忽略不计。由于采用了移去恢复法，残余地面重力扰动量级很小，地面重力解析延拓误差也可忽略^[23]。边缘计算点的边界效应可采用扩大数据点范围或内缩计算点范围的方式来消除。因此，为了估算航空重力数据误差，需要分离的主要是地面重力格网插值误差和代表性误差在飞行高度处的影响。

地面重力格网插值误差和代表性误差可采用交叉检验法进行估计。具体方法是：若欲估计某个地面重力点的误差，先利用所有离散地面重力点格网化得到重力格网，再由重力格网内插得到该观测点处的内插重力值，原始重力值与内插重力值之差τ为

式中，δg_Ori为原始重力扰动值；δg_Int为内插重力扰动值；差值τ近似由地面重力格网插值误差和代表性误差组成。

利用交叉检验法得到地面重力格网插值误差和代表性误差后，还需要利用向上延拓方法得到飞行高度处的地面重力误差影响。根据误差传播定律，扣除地面重力误差影响后，空中航空重力扰动数据的中误差可表示为

式中，σ_Δ为航空重力与地面重力向上延拓值互差的标准差；σ_τ为地面重力误差向上延拓传播到空中点后的标准差。图 2为航空重力数据误差估计的流程。

5 毛乌素测区航空重力空中测线数据外符合精度评价 5.1 测区和数据情况

2015年4月17日至6月11日，国家测绘地理信息局第一大地测量队在毛乌素地区开展了航空重力测量试验性项目。毛乌素测区位于陕西与内蒙古交界处，地貌以沙漠丘陵为主，地形起伏如图 3所示，绘图采用SRTM数字高程模型数据^[24]，该区域平均地形高度1248 m，最小高度730 m，最大高度2420 m。飞机选用塞斯纳208机型，机上搭载GT-2A航空重力测量系统和1台Ashtech双频GPS接收机，地面架设2个GPS基准站。选取鄂尔多斯机场作为飞机起降保障机场，机场距毛乌素测区的飞行距离约为80 km。

航空重力空中测线覆盖范围为纬线(38.77°N—39.32°N)和经线(109.23°E—110.52°E)围成的区域，布设有32条测线，其中东西方向11条，南北方向21条，平均飞行高度为2087 m，平均飞行速度约260 km/h，数据采样率为2 Hz，测线间距为5 km，共包含71 773个数据点，测线数据采用100 s的卡尔曼滤波消除高频噪声。在地面上按照5 km平均间距布设了459个地面重力点，向航空重力数据范围四周各延伸10 km，覆盖范围为纬线(38.67°N—39.41°N)和经线(109.12°E—110.63°E)围成的区域。在每个地面重力点上实施了约30 min的GPS静态定位观测以测定其大地高，GPS数据处理采用WayPoint软件，与附近的IGS站点组网解算得到地面点坐标，大地高精度约为10 cm，按照索米里安公式^[19]严密计算观测点处的正常重力值，参考椭球选取CGCS2000椭球。航空重力测线与地面重力点的空间分布如图 4所示，其中绿线为航空重力空中测线，红点为地面重力点。表 1列出了航空和地面重力扰动观测值的统计信息。

5.2 地面重力误差估计及其空中传播

利用交叉检验法估算地面重力数据的格网插值误差和代表性误差，格网化方法选取反距离加权插值，选择这一方法的原因是毛乌素测区的地面重力点分布比较合理、密集，简单的反距离加权插值即可达到较好的格网化效果，另外该方法对数据的处理较为简单，对数据的干预少，更能保持数据本身的特性。对459个地面点进行反距离加权插值格网化，得到2′×2′重力扰动格网，再内插到离散观测点处，将内插重力扰动与原始重力扰动做差得到地面重力误差估值。分别采用梯度法、泊松积分法和FFT法，计算得到地面重力误差对空中点的影响，其中泊松积分法和FFT法需要先将地面重力离散点的误差估值格网化得到2′×2′误差格网。表 2列出了地面重力扰动误差估计和3种方法空中误差传播的统计信息，地面重力数据误差的标准差为1.69 mGal，向上延拓传播到空中点的误差影响分别为0.66、0.78和0.92 mGal，不能忽略。图 5—图 7分别绘制了3种向上延拓方法所得地面重力扰动误差对空中点的影响。

5.3 航空重力数据外符合精度

利用3种向上延拓方案，对毛乌素航空重力空中测线重力扰动观测值的外符合精度进行评价。第1种方案是利用梯度法将地面重力扰动向上延拓至空中相应点位，将航空重力测线重力扰动内插到空中延拓点处，内插值与地面重力扰动向上延拓值进行比较，式(1)中的梯度直接采用2′×2′地面重力扰动按照式(2)计算，积分半径选取0.5°。第2种方案是利用泊松积分法将地面重力扰动向上延拓至空中测线处，再将航空重力测线重力扰动与地面重力扰动向上延拓值进行比较，泊松积分计算采用移去恢复法，球冠积分半径选取0.5°，参考重力场模型选取EIGEN-6C4^[25]，地形数据采用3″分辨率的SRTM数字高程模型^[26]。第3种方案是利用FFT方法将地面重力扰动向上延拓至空中，采用与方案2相同的移去恢复处理。基于这3种向上延拓方案，航空重力扰动与地面重力扰动空中延拓值的差值统计信息见表 3。图 8—图 10分别绘制了3种向上延拓方案所得航空与地面重力扰动的差值。如前所述，表 3中的差值主要包括航空重力数据误差和地面重力误差影响，采用式(17)扣除地面重力误差影响，3种向上延拓方案得到航空重力数据的中误差分别为4.7、2.94和2.83 mGal。可以看出，扣除了地面重力误差影响后，在航空与地面重力扰动差值的基础上，3种方案估计得到的航空重力扰动误差均有所减小，梯度法所得误差最大，泊松积分法与FFT法非常接近。

毛乌素测区数据范围在纬线方向跨越0.74°、经线方向横跨1.51°，泊松积分公式的球冠积分半径选取为0.5°，因此上述对比结果不可避免地包含了因数据范围有限导致的边界效应，见图 8—10。为对航空重力测线观测值精度进行准确评估，需消除因地面重力向上延拓带来的边界效应，但因为毛乌素测区地面重力数据范围有限，本文试验中无法完全消除边界效应，只能削弱其影响。在数据点范围基础上，将参与对比的计算点范围内缩，可削弱边界效应。图 8—图 10中的矩形区域为内缩后的对比计算点范围。表 4分别列出了3种方案内缩计算范围后航空重力扰动与地面重力向上延拓值的差值统计信息。可以看出，3种方案的对比结果均受到边界效应的影响，梯度法的边界效应更为显著。将计算点范围在纬向内缩0.1°、经向内缩0.25°后，航空与地面重力的差值显著减小。在内缩后的区域，地面重力误差向上延拓传播到空中点后的标准差分别为0.31、0.59和0.71 mGal。结合表 4，利用式(17)扣除地面重力误差影响，3种向上延拓方案所得航空重力数据的中误差分别为1.62、1.5和1.42 mGal。

需要指出的是，这一内缩处理无法完全消除边界效应，只是在一定程度上削弱了边界效应，估计得到的航空重力数据误差中还包含有残余边界效应。此外，从表 4中的平均差值可以看出，毛乌素测区的航空重力与地面重力数据之间基本不存在系统偏差，这对于GT-2A航空重力数据的后续应用，如大地水准面建模是非常有益的。

6 结论

利用毛乌素测区实测航空和地面重力数据，分别对基于梯度、基于泊松积分和基于FFT的地面重力向上延拓这3种航空重力空中测线数据外符合精度评价方法进行了对比分析，结果表明：

(1) 航空重力空中测线数据外符合精度评定必须扣除地面重力格网插值误差和代表性误差对空中点的影响，交叉检验法可有效估计这两项主要误差影响。毛乌素测区地面重力数据误差的标准差为1.69 mGal，采用3种向上延拓方案传播到空中点的误差影响分别为0.66 mGal、0.78 mGal和0.92 mGal，不能忽略。

(2) 基于泊松积分和基于FFT的地面重力向上延拓方案均能够客观评价航空重力空中测线观测值的外符合精度，二者表现相当，梯度法表现稍差。在包含残余边界效应的情况下，扣除地面重力误差影响后，毛乌素测区GT-2A航空重力空中测线重力扰动的外符合精度优于1.42 mGal。

(3) GT-2A航空重力观测值与地面重力观测值之间的系统偏差很小，这对大地水准面建模等航空重力数据的后期应用是非常有益的。