向下延拓是航空重力的关键问题^[1]，是将空中重力异常值归算到大地水准面，再求解边值问题得到大地水准面高等其他重力场元^[2-3]。经过大地测量学家多年理论研究^[4]，发现当地球外部扰动位满足边值条件后^[5]，可将重力异常“向上延拓”至球外，此过程可用泊松积分方程表示，而多年来受限于泊松积分逆算子的研究进展，多采用迭代算法^[2]，文献[6]便使用逆泊松积分迭代求解法，并给出了向下延拓的“良态条件”(向下延拓高度和格网间隔的比值)，但迭代使观测噪声影响累积剧增，因此被直接求逆的代替。此后有了在求逆过程中加入正则化参数来抑制噪声的方法^{[1, 7-8]}，正则化方法研究已较为成熟，具体算法有Tikhonov正则化算法^[9-10]、基于Landweber迭代法的正则化算法和截断奇异值分解法等^[11]，文献[7-11]中正则化算法在模拟条件下精度可到达毫伽级(mGal=10^-5 m/s²)。而当超越“良态条件”时，即使正则化算法也无法进行有效计算，需使用广义岭估计算法求解病态问题^[12]，在病态条件下，文献[11]利用广义岭估计提高了向下延拓精度。此后，国内外专家为避开直接求逆，先后研究使用最小二乘配置法^{[4, 13-15]}、虚拟点质量法^[2-3]等间接求逆方法，但在过程中仍涉及矩阵求逆过程，仍存在不稳定问题^[16]。另外为减少重力场“中长波”向下延拓误差影响，使用“移去恢复法”^{[1, 7-8, 10-12, 17]}。此外基于差分思想提出联合超高阶模型和地形信息求取向下延拓改正数，并在文献[16, 18-19]的试验结果中得到较高精度。

向下延拓的目标是得到高精度的地面重力数据，确定高精度大地水准面，上述研究方法多种多样，但是距此仍有距离。查阅航空重力文献资料^{[13, 17, 20-21]}，发现实际向下延拓过程中可能受到偶然误差和系统误差影响，但尚未有资料研究处理向下延拓过程中的系统误差和偶然误差^{[1, 7-8, 12]}。就处理系统误差而言^[22]，常用方法有交叉点平差、重复测线比较和相邻测线比较，首先利用外部重力异常确定交叉测线、重复测线和相邻测线的精度，然后比较它们与正常测线之间的差值改正系统误差，这些方法本质是需要“外部数据”，但极区或者高山地区进行航空重力测量时，并没有高精度的“外部数据”，系统误差改正甚至发现成为难题；此外，偶然误差是由测量过程中多种原因引起，其在向下延拓过程中亦会放大^{[17, 21]}。因此，本文将研究分析两类误差在向下延拓过程中的变化特性，利用误差传播规律估计不同误差影响的数值，分析两类误差影响的变化规律并提出处理方法。采用半参数模型处理系统误差^[23-28]，采用正则化算法处理偶然误差^[8-10]，通过试验验证各方法的有效性，提出两步法处理向下延拓过程中的系统误差和偶然误差。

1 逆泊松积分和移去恢复法试验

使用移去恢复法可减少重力场“中长波”向下延拓误差，但有其“局限性”，以逆泊松积分模型为例说明^[6]。

1.1 逆泊松积分

泊松积分方程为

式中，Δg_h(r, φ, λ)是空中重力异常，其中r是R(地球平均半径)和飞行高度h的和, φ、λ为球坐标下的纬度和经度；Δg_g(R, φ′, λ′)是地面重力异常值，φ′、λ′分别是、地面点的纬度和经度；l为球距离。

逆泊松积分模型可表示为

式中，B代表设计矩阵; i=1, 2, …, M; M为观测点(空中重力异常点)个数，N是待求点个数。

非对角元素为

式中，Δσ_j为待求点所处网格的面积；N_C是处在积分计算内区的流动点的个数；ψ₀为1°，通过最小二乘法平差得到向下延拓值。

1.2 试验数据说明

本文试验的飞行高度、格网间隔和测量范围的设置参照美国大地测量局CS01航空重力项目。具体来说，飞行高度为6.3 km，格网间隔为6′，纬度范围27°N-31°N，经度范围为87°W-88°W，向下延拓高度差为6.3 km。采用重力场模型模拟空中重力异常：首先使用EGM2008模型(最高2190阶)计算飞行高度处无误差的空中重力异常Δg_{h_ggm}，可理解为“无观测误差”；然后根据实际情况定量加入偶然误差和系统误差形成含误差的空中重力异常值Δg_h，具体为Δg_{h_ggm}加入系统误差s_simu和偶然误差(标准差2 mGal(1 Gal=10^-2 m/s²))得到Δg_h，计算s_simu参考文献[1]中航空重力系统误差类型和计算公式，其中s_simu类型包括恒值系统误差和变值系统误差(线性变化系统误差)，s_simu标准偏差约为1.3 mGal(见表 1)，线性变化系数约为0.002。EGM2008模型在此试验区域的精度较高，因此将其生成的地面重力异常作为本文此区域试验的检核值。需要说明的是，在实际测量中，建议采用实测地面重力数据作为检核值。

1.3 移去恢复法试验

对空中重力异常进行向下延拓时，常常用到移去恢复法。为了说明向下延拓结果随移去恢复阶次增加的变化情况，使用逆泊松积分分别对Δg_{h_ggm}和Δg_h进行试验。

表 2为移去恢复法试验的差值统计情况。其中Δg_h的结果中(见表 2、图 1)，移去恢复360阶相比0阶(未移去恢复)有显著改善，而使用中高阶重力场模型(如720阶或1080阶)进行向下延拓移去恢复时，精度提高相对较小。当选用Δg_{h_ggm}(无观测误差，见表 2、图 2)，结果精度较高。究其原因，应与偶然误差和系统误差的影响有关，含系统误差影响和偶然误差影响的结果精度较差，说明实际计算中单一使用移去恢复法有局限性，需加入对系统误差和偶然误差的处理策略。

需要补充的是，并不是所有地区均可使用中高阶重力场模型移去恢复。从根本上来说，其原因与重力场模型建模时加入该地区地面重力数据量有关，加入数据多的地区中高阶精度高，未加入地面数据的地区精度差。因此，地球重力场模型在不同地区精度相差较大，同时移去恢复法可使模型阶次上限不同。例如EGM08模型建模时加入大量北美地区地面重力数据^[29-30]，可使用中高阶次模型进行向下延拓移去恢复，而在中国特别是青藏高原区域加入地面重力数据稀少，只可使用低阶次。

2 系统误差影响

航空重力系统误差包括多种类型^[31]。可分为恒值系统误差(偏差)和变值系统误差(漂移)。本节利用误差传播规律，推导公式并模拟计算系统误差影响。采用半参数模型估计系统误差。

2.1 系统误差影响分析

当观测值中既含有偶然误差D，又含有系统误差s时，称观测值含有综合误差W，可表述为

观测值含有综合误差时，综合误差的方差简称综合方差。根据系统误差与偶然误差的联合传播律，当观测值中含有综合误差时，逆泊松积分向下延拓值的综合方差为

式中，J=(B^TPB)^-1B^TP，其中P为权，Q为协因数(单位阵)；s₀²单位权方差。是逆泊松积分参数估值；s是系统误差。

式(6) 右边第2项即系统误差对综合方差的影响，因此系统误差影响可表示为

式中，t为必要观测数。上式中系统误差s确定时，变量为设计矩阵B。通过分析设计矩阵公式，发现其与格网间隔、向下延拓高度相关。

为分析系统误差影响与格网间隔(图 3中为grid)、延拓高度(图 3中为h_gap)之间的关系，设计试验，将系统误差分为3类：恒值系统误差(图中为cse)、变值系统误差(图中为vse)和总系统误差(图中为hse)，三者的数值、设定方法与1.2节一致。图 3(a)为系统误差影响与格网间隔的关系，它的延拓高度6.3 km。格网间隔为3′时，3 mGal恒值系统误差的影响约4.3 mGal，从中看出，随格网间隔增大，恒值和变值系统误差影响减弱(见表 3)。图 3(b)为系统误差影响与延拓高度的关系，格网间隔为6′，4 km延拓高度(总系统误差影响)约4.1 mGal，6.3 km时(文中试验高度)约6 mGal。而系统误差影响随着延拓高度增加而变大，说明减小向下延拓高度有利于减弱系统误差影响。

2.2 半参数模型应用于向下延拓

半参数模型形式为

式中，L为n维观测向量；X为t维参数向量，t是必要观测数；Δ为n维偶然误差向量；非参数向量S，可用自然样条函数表示^{[24-25, 28]}。选择它的原因是：航空重力系统误差包括偏差、漂移和周期性误差，表现为一条光滑的曲线，需选择曲线函数去拟合，而自然样条函数是拟合最好的曲线插值函数。

计算采用补偿最小二乘原理

式中左边第1项残差平方和，α_S>0为光滑参数，第2项自然样条函数补偿项。自然样条函数补偿项和本文光滑参数α_S的求解采参见文献[24]。

计算时采用条件极值的拉格朗日乘数法构造函数。参数和系统误差的估值为

式中，M=(P+α_SK)^-1P。

半参数模型应用于向下延拓(简称半参数模型)：当区域内有“外部重力数据”时，半参数模型用空中重力异常减去“外部数据”得到剩余重力异常，然后把它作为观测值估计系统误差，估计系统误差和向下延拓两步完成；当区域内无“外部重力”时，半参数模型将空中重力异常值即为观测值，利用自然样条函数可表示系统误差^[24-26]，最重要的是，它利用补偿最小二乘法和光滑参数的优化算法，实现“不用外部数据”得到系统误差和地面重力异常，估计系统误差和向下延拓一步完成。另查阅文献发现半参数模型在其他应用中^[32-35]，亦可同时求解系统误差和参数，即在“求解过程中”估计系统误差，一步完成。

2.3 半参数模型试验

为检验半参数模型(“有外部数据”和“无外部数据”两种情况下)估计系统误差的有效性，设计对比试验。利用半参数模型和逆泊松积分两种方法，空中重力异常依然采用Δg_h，延拓高度和数据范围等与1.2节中描述一致，移去恢复模型选前360阶EGM2008模型。

有“外部数据”时，传统方法可估计偏差，也就是恒值系统误差，而无“外部数据”时，无法估计系统误差(见图 4和表 4)。但是，半参数模型在两种情况下，均可估计系统误差，有“外部数据”时，半参数模型较精确地估计了系统误差，而无“外部数据”时，系统误差估计值精度虽然下降，但相比传统方法，依然优势显著。更重要的是，“无外部数据时”，半参数模型“在向下延拓过程中”不但可估计系统误差，亦可提高向下延拓精度(见表 5)。

3 偶然误差影响 3.1 偶然误差影响分析

根据公式(12)，偶然误差影响可定义为

当s₀²确定时，偶然误差影响与逆泊松积分设计矩阵有关，而设计矩阵与格网间隔、延拓高度有关。为分析偶然误差影响与格网间隔(图中为grid)、延拓高度(图中为h_gap)的关系，设计试验，3组中偶然误差量级(标准差)分别为2 mGal、4 mGal和6 mGal。图 5左图(延拓高度为6.3 km)中，随格网间隔增大3组偶然误差影响减小。右图(格网间隔6′)中，随延拓高度升高3组偶然误差影响增加。因此，偶然误差影响变化规律与系统误差影响一致。

3.2 正则化算法试验

在参数估计中加入正则化参数可以减小偶然误差影响^{[1, 7]}，本文选用Tikhonov正则化算法，Tikhonov正则化解为

式中，I为单位阵；α_R为大于零的正则化参数。正则化参数使用广义交互确认法选择^[11-12]。

试验采用数据Δg_{h_ggm}，但每组加入偶然误差量级不同，分别为2 mGal、4 mGal和6 mGal。移去恢复模型选前360阶EGM2008模型。

表 6分别为正则化解和逆泊松积分解的差值统计。随偶然误差量级的增加，逆泊松积分解的差值均方根从约6.3 mGal增加到约8.9 mGal，这是由于向下延拓导致的偶然误差放大，而正则化解的差值均方根误差变化并不显著，向下延拓正则化解受偶然误差影响比逆泊松积分小。因此，在逆泊松积分模型中加入正则化参数可减弱偶然误差影响。

当空中重力异常改用Δg_h(含有系统误差和偶然误差)。由于受到系统误差影响，正则化算法的差值均值约为2~3 mGal(见表 7)，最小值与最大值绝对值相差较大，说明正则化算法受到了系统误差影响。考虑半参数模型可估计系统误差，应使用基于半参数模型和正则化算法的两步法，首先使用半参数模型估计系统误差，再使用正则化算法，这样可先后有效减弱系统误差影响和偶然误差影响，显著提高向下延拓结果精度，此试验中精度最高可达约2.3 mGal。

上述结果说明两步法较单一方法更为有效，为验证此结论在其他地区的适用性，另随机选取某验证区域实施对比试验，此区域位于海洋区域，向下延拓高度11 km，格网间隔6′，所加误差等信息之前给出信息一致。模拟观测值仍采用EGM2008模型计算的重力异常加模拟观测误差的策略，模拟观测误差和上述试验一致，检核值仍旧使用EGM2008计算的地面重力异常值。利用EGM2008模型前360阶次进行移去恢复，结果(表 8)中两步法较单一方法的结果优势显著。具体来看，逆泊松积分结果受误差影响较大，半参数模型的结果受系统误差影响较小，正则化算法可显著减弱误差的影响，但系统误差难以剔除，两步法可同时减弱系统误差和偶然误差的影响，数据精度显著提高。

4 结语

本文分析向下延拓过程中系统误差和偶然误差的变化特性，并给出处理方法，具体如下：

(1) 通过试验说明移去恢复法的有效性和局限性，说明使用中高阶重力场模型进行移去恢复法时，需处理系统误差和偶然误差。

(2) 推导公式并计算系统误差影响数值，分析其变化特性及与格网间隔、延拓高度之间的关系，增大格网间隔和降低飞行高度可减弱系统误差影响，应用半参数模型处理系统误差影响，设计试验说明其有效估计系统误差数值。

(3) 推导偶然误差影响公式，试验分析其与格网间隔、延拓高度之间的关系，说明偶然误差变化特性与系统误差一致，试验说明正则化算法可减弱偶然误差影响。

(4) 综合分析各类误差影响后，提出基于半参数模型和正则化算法的两步法，试验说明新方法有效作用于各类误差影响，较单一方法精度更优。本文试验中，使用EGM2008模型计算重力异常数据，加入标准差为2 mGal的偶然误差和模拟系统误差，在向下延拓高度6.3 km和格网间隔6′的条件下，两步法的模拟试验精度可达2.3 mGal。

考虑实际航空重力测量精度为mGal级，为减小误差对向下延拓结果的影响，应采用两步法：先用半参数模型估计系统误差，然后加入正则化参数抑制偶然误差影响。本文方法分析处理了航空重力向下延拓过程中的系统误差和偶然误差，为得到高精度地面重力异常提供思路，且为研究航空重力系统误差和偶然误差提供方法。