它由p个多项式组成，多项式的自变量 为经度、炜度，因变量为重力异常，每个多项式被称为一个因子，这部分模型被称为回归模 型，主要用来描述重力场的非平稳性。
[0067] 为了保证建模的精度，在建模过程中使用了多元逐步回归法，其本质思想就是将 采样点的空间位置X及其所对应的重力异常值Y，视为回归模型的样本，根据该样本，从全 部因子（也即多项式）集合中，筛选出对重力异常有显著影响的因子（也即多项式），构 成回归模型。多元逐步回归法的引入，保证了被选入回归模型的多项式，其阶数、结构形式 完全由重力场空间分布决定，结合步骤三中采用的拉格朗日乘数法计算出的权重系数β， 即可保证回归模型对重力场估计的无偏性，进而保证公式（1)中剩余的重力异常信号Z = Y-FP均值为零，使得用ζ来描述重力场自相关性成为可能。多元逐步回归法与拉格朗日乘 数法相结合，有效地将"非平稳场建模问题"转化为"广义平稳随机过程建模问题"。
[0068] 建立回归模型的具体计算步骤如下：
[0069] (1)建立由全部因子xl, x2,…，χη组成的自变量集合，设自变量集合包含η个因 子（η>ρ);并设定因子被选入回归模型的显著性水平为α 1，被剔除出回归模型的显著性水 平为 α 2,0〈 α 1 < α 2〈1。
[0070] (2)在背景图上抽取m组数据（经炜度与重力值，如公式（3)所示），作为回归模 型的样本。
[0071] (3)从全部η个因子中，抽取出任意两个因子xi, xj (1彡i彡n, 1彡j彡η)，假定 空间重力场可以用选取出的两个因子描述，依据m对样本值，利用最小二乘法计算其参数 向量τ i，jl，得到用这两个因子描述的空间重力场ykl = xi, xj* τ i, jl，并计算这个模型 的剩余平方和Si, jl = k = lmyk-ykl2。
[0072] ⑷按照步骤⑶的方法，遍历全部因子，计算任意两个因子所组成的模型的剩余 平方和Si, jl (1彡i彡n，1彡j彡η)，组成离差矩阵Sl :
[0073] SI = Sl, 11 …SI, nl SI, yl............Sn, 11... Sn, nl Sn, yl
[0074] 离差矩阵的最后一列表征了以单一因子描述重力场时的绝对误差，误差越小，对 应的因子对重力场空间分布的影响越显著。
[0075] (5)为了便于比对不同因子对重力场影响的显著程度，计算各个因子在离差矩阵 中的贡献Vjl = (Sj，yl)2Sj，jl，设贡献最大的因子编号为kl，对其做显著性水平为a 1的 显著性F检验：
[0076] F = VjlSEl(M-I-I)
[0077] 其中，SEl = k = lmyk-y2_Vjl，y = k = lmykm。
[0078] 当F>Fa 时，该因子选入回归模型，步骤（3)-(5)为第一轮筛选。
[0079] (6)在kl因子选入回归模型后，需要调整离差矩阵：
[0080] S2 = Sl，12-Sl，n2Sl，y2............Sn, 12... Sn, n2 Sn, y2
[0081] 其中，Si, j2 = Ski, jlSkl，kll 当 i = kl，j 乒 klSi，jl-(Si，kll*Skl，jl)Skl，kll
[0082] 当 i 辛 kl，j 辛 kllSkl，kll 当 i = j = kl (-Si, kll) Ski, kll 当 i 辛 kl，j = kl
[0083] 按照步骤（5)的方法，重新计算剩余因子在新离差矩阵中的贡献，并对贡献最大 的因子k2做显著性水平为α 1，自由度为m-2-l的显著性F检验，作为其是否被引入回归模 型的判据，此为第二轮筛选。
[0084] (7)在选入两个因子后，首先按照步骤（5)的方法，将离差矩阵S2调整为S3,进行 第三轮筛选。在选入新因子之前，先对已经选入的两个因子，按照被剔除出回归模型的显著 性水平α 2做显著性F检验，以判断是否需要将其从回归模型中剔除（F>Fa 21，m-3-l时保 留，F〈F α 21，m-3-l时剔除），之后再引入新的因子。
[0085] (8)按照第三轮筛选的方法，对自变量集合中的因子进行遍历，直到新的因子不再 被引入回归模型，同时回归模型内的因子也不再被剔除时，即可获得重力场的多项式回归 模型。逐步回归法的计算步骤如图1所示。
[0086] 步骤2、针对空间中的重力异常模型建立相关模型
[0087] 本步骤是对公式（1)中的Z(X)建模，这部分模型被称为相关模型，主要用来描述 重力场的相关性。
[0088]由于在计算插值系数之前，无法获得ζ (X)的准确数值，因此需要借助重力矩阵Υ， 对重力场的自相关性进行建模描述。为了消除均值变化对相关模型精度的影响，建模时需 要首先对重力矩阵Y规范化，进行平差处理，使得插值范围内的重力异常均值为0。规范化 后的重力矩阵Υ'，可以用广义平稳随机过程来描述，可以代替Z (X)，来表征重力场的相关 性。
[0089] 相关模型建模的本质思想，就是将Υ'的相关函数视为相关模型的样本，通过选取 合适的相关模型Ch，Θ，来满足下列约束条件：
[0090] Ch, Θ = CovY，x+h，Y，X = E {Y，x+h-EY，x+h*Y，x-EY，x} = EY，x+h*Y，X
[0091] 可以看出，采用相关模型表征重力场的相关性，不仅描述了采样点与待插值点之 间的相对位置及相关关系，还考虑了采样点的空间分布结构特点，完整地描述了重力场的 自相关性。为了行文方便，下文中，Z与Y'依推导需求出现，二者可以相互替代。
[0092] 相关模型的建立步骤如下所示：
[0093] (1)对重力矩阵Y进行平差处理，得到Y'。
[0094] (2)计算Y'的相关矩阵，它是相关模型的一个样本估计。对于空间中任意两个采 样点位置Xi与xj，Y'的相关矩阵如下式所示：
[0095] EY' xi*Y' xj = Ri，j = 〇 2R Θ，xi，xj (4)
[0096] 上式中，σ 2为背景图重力值的方差。
[0097] (3)选取合适的相关模型。
[0098] 在获得Υ'的相关矩阵（公式（4))后，需要寻找合适的相关模型对其进行拟合。一 般而言，对于一个强度在空间上连续可导的矢量场，其相关模型近似于指数形式。相反地， 其相关模型近似于线性形式。当载体沿东西方向航行时，某航迹上重力异常的相关部分，如 图2所示。经过计算，其相关函数，如图3所示。从图3可以看出，由于重力场在空间上连 续可导，因此Z的相关矩阵接近于指数相关模型，本算法中，使用高斯相关模型对Z进行描 述。高斯相关模型的相关函数形式如下：
[0099] CO = σ 2 = Cl
[0100] Ch = Cle_h θ 2
[0101] Chh-① O
[0102] 其相关函数如图4所示（图中〇2 = 1，θ = 1. 7308h)。
[0103] 从图3可以看出，随着间隔距离增大，Y'的相关性都会减小。在高斯相关模型中， 包含一个参数Θ，称为相关尺度：Θ越大，Y'的相关矩阵取值随距离增大减小得越快，Y'越 接近随机信号；相反，Θ越小，Y'越接近常值。
[0104] ⑷计算相关尺度Θ。
[0105] 结合高斯相关模型的样本估计（公式（4))，依据最优原则，利用最小二乘法，Θ的 最优解可以表示为下式：
[0106] min Θ = c Θ E Rlm 〇 2
[0107] 其中，R为矩阵Ri，j所对应的行列式值，σ 2为背景图重力值的方差。通过调整相 关尺度Θ，结合高斯模型，来逼近Z的相关矩阵Ri，j，构造和描述相关部分。
[0108] 步骤3、计算回归模型和相关模型的权重系数，代入公式（1)，即可获得待插值位 置上的重力异常空间分布表达式，再将待插值位置代入公式（1)，即可计算出待插值位置的 重力值。
[0109] 本步骤在建立了重力场的空间分布模型后，该算法基于最优化原理，利用拉格朗 日乘数法计算确定性部分的权重β i，1及相关部分的权重γ，求解权重矩阵。设得到的权 重矩阵为c，为一个l*m维的权重矩阵，则插值算法在待测点处的估计值为：
[0110] yx = cTY (5)
[0111] 插值算法的目的，就是求解出矩阵C的取值，进而计算出回归模型的权重β i，1及 相关模型的权重γ，以满足估计的无偏和最优原则。
[0112] 由公式⑶、（5)可得，插值结果的误差为：
[0113] yx-yx = cTY-yx = cTF β +Z-fxT β +z = cTZ-z+FTc-fxT β
[0114] 由于Z为广义平稳随机过程，因此Z和z均值为相等的常数，且对于位于地球上的 任