意位置的卫星反演图（任意β值），要满足估计是无偏的（即上式恒等于0)，必须保证：
[0115] FTc-fx = 0 (6)
[0116] 由公式（3)、（5)可得，插值误差的方差为：
tong] 要使得估计是最优的，必须使得φχ取值最小，结合公式（6)的约束条件，采用拉 格朗日乘数法，建立拉格朗日函数：
[0120] Lc, λ = σ 21+cTRc-2cTr-λ TFTc-f
[0121] 式中，c和λ为待求解量。对c求导，并令其等于〇,得到：
[0122] Lc' c，λ = 2 σ 2Rc-r_F λ = 〇 (7)
[0123] 结合公式（6)和（7)，可以解得：
[0124] λ I = FTR-lF-lFTR-lr-f c = R-Ir-F λ 1
[0125] 其中，λ I = - λ2σ 2
[0126] 将解得的λ和c代入公式（5)，可以获得待插值位置重力场的解析表达式：
[0127] yx = r-F λ1TR-1Y = rTR-lY-FTR-lr-fTFTR-lF-lFTR-lY (8)
[0128] 对于公式（8)，按照公式（I)的结构提取合并，可以得到：β i，I
[0129] = FTR-1F-1FTR-1Y (9)
[0130] 将β i，1代入公式（8)，整理后可以得到：
[0131] yx = rTR-lY-FTR-lr-fT β i, I = f (x) T β i, 1+r (x) T γ (10)
[0132] 式中，γ = R-I(Y-F) β i，I (11)
[0133] 观察公式（10)和公式（1)，f(x)和r(x)分别对应于事先拟合出的回归模型和 相关模型。一旦确定了背景图的经炜度范围，权重矩阵β i，1及γ也可计算得到（公式 (9)、公式（11))。因此，将待插值位置X代入公式（10)，即可得到待插值位置的重力异常 值。仔细分析公式（11)，在 γ = R-I (Y-F) β i, 1 中，Ri, j = Ezxizxj = O 2R Θ，xi，Xj 描 述了重力场在任意两个采样点xi和xj之间的空间分布特点。而相关矩阵rx = Ezxzxj = σ 2R Θ，X，Xj则表征了待插值位置X与采样点之间的重力场空间分布特点。因此，该算法 在满足最优无偏的基础上，综合考虑了重力场空间分布特点以及采样点之间的相互作用影 响，完整地描述了空间重力异常分布。
[0134] 下面给出一个仿真验证，对本发明的效果进行验证。
[0135] 利用实测重力异常数据中的1#、2#、3#、5#、6#、7#测线，构造重力背景图；分别用 三次样条插值法和本发明提出的算法，在构造出的背景图上插值计算4#测线重力异常，并 与4#测线实测数据相比对，如图5所示。图中曲线分别代表实测重力异常、三次样条插值 结果及本专利提出插值算法的插值结果。
[0136] 无论哪种插值算法，其插值结果均存在误差：如果插值算法专门针对重力场的相 关性进行了建模描述，则插值误差只包含常值信号和随机信号，这两种信号的变异函数均 为常数；如果插值算法在建模过程中忽略了重力场的相关性，插值误差就包含随机游走信 号，其变异函数取值就不是常数。因此，可以通过计算插值误差变异函数的方法，来评判插 值算法的性能：如果误差的变异函数取值为常数，则可以判定所得到的模型准确反映了重 力场的相关性，该算法的性能较好。
[0137] 两种插值算法的误差，其变异函数分别如图6及图7所示。由变异函数可以看出， 本专利算法的插值误差变异函数为常值，所建立起来的模型准确的描述了重力场的相关 性。而三次样条插值误差中存在明显的随机游走误差，效果较差。
[0138] 两种插值算法的精度为：
[0140] 由上表看出，本专利算法的插值结果满足无偏最优的要求。
[0141] 本发明不仅局限于上述的【具体实施方式】，只要插值算法特别针对重力场的相关性 进行了建模描述，对其确定性部分和相关部分分别采用不同的方式建模，并应用最优化方 法确定两部分的权重系数，均落入本专利的保护范围内。
【主权项】
1. 一种基于泛克里金法的重力场插值方法，其特征在于包括以下步骤： 步骤1、针对空间中的重力异常模型建立回归模型； 步骤2、针对空间中的重力异常模型建立相关模型； 步骤3、计算回归模型和相关模型的权重系数，代入重力异常模型，便可获得待插值位 置上的重力异常空间分布表达式，再将待插值位置代入重力异常模型，便可计算出待插值 位置的重力值。2. 根据权利要求1所述的基于泛克里金法的重力场插值方法，其特征在于：所述重力 异常模型为： yx = f lx*** fpx^ 0 i, 1+zx 式中，X为空间位置，是一个1*2维位置矩阵，第一列保存经度，第二列保存炜度；yx为 模型解算出的重力值。3. 根据权利要求1或2所述的基于泛克里金法的重力场插值方法，其特征在于：所述 步骤1建立回归模型的方法包括以下步骤： (1) 建立由全部因子xl, x2,…，xn组成的自变量集合，设自变量集合包含n个因子， 并设定因子被选入回归模型的显著性水平为a 1，被剔除出回归模型的显著性水平为a 2， 0<a I ^ a 2<1 ； (2) 在背景图上抽取m组经炜度与重力值数据，作为回归模型的样本； (3) 从全部因子中，抽取出任意两个因子xi，xj，依据m对样本值，利用最小二乘法计算 其参数向量T i, jl，得到用这两个因子描述的空间重力场ykl = xi, xj* T i, jl，并计算这 个模型的剩余平方和Si, jl = k = lmyk_ykl2 ; (4) 重复执行步骤（3)，遍历全部因子，计算任意两个因子所组成的模型的剩余平方和 Si, jl，组成离差矩阵Sl ; (5) 计算各个因子在离差矩阵中的贡献Vjl = (Sj，yl)2Sj，jl，设贡献最大的因子编号 为kl，对其做显著性水平为a 1的显著性F检验： F = VjlSEl(M-I-I) 其中，SEl = k = lmyk-y2_Vjl，y = k = Imykm ; 当F>Fa ll，m-1-l时，该因子选入回归模型； (6) 在kl因子选入回归模型后，调整离差矩阵S2 : 按照步骤（5)的方法重新计算剩余因子在新离差矩阵中的贡献，并对贡献最大的因子 k2做显著性水平为a 1，自由度为m-2-l的显著性F检验，作为其是否被引入回归模型的判 据； (7) 在选入两个因子后，首先按照步骤（5)的方法，将离差矩阵S2调整为S3,在选入新 因子之前，先对已经选入的两个因子，按照被剔除出回归模型的显著性水平a 2做显著性F 检验，以判断是否需要将其从回归模型中剔除，F>Fa 21，m-3-l时保留，F〈Fa 21，m-3-l时 剔除，之后再引入新的因子； (8) 按照步骤（7)筛选的方法，对自变量集合中的因子进行遍历，直到新的因子不再被 引入回归模型，同时回归模型内的因子也不再被剔除时，即可获得重力场的回归模型。4. 根据权利要求1或2所述的基于泛克里金法的重力场插值方法，其特征在于：所述 步骤2建立相关模型的方法包括以下步骤： (1) 对重力矩阵Y进行平差处理，得到Y' ； (2) 计算Y'的相关矩阵，对于空间中任意两个采样点位置xi与xj，Y'的相关矩阵为： EY' xi*Y' xj = Ri, j = 〇 2R 9，xi, xj 上式中，〇 2为背景图重力值的方差； (3) 使用尚斯相关t旲型对Z进彳丁描述，尚斯相关_旲型的相关函数形式如下： CO = 〇 2 = Cl Ch = Cle-h 0 2 Chh 0 (4) 计算相关尺度0的最优解： min 0 = c 0 = Rlm 〇 2 其中，R为矩阵Ri，j所对应的行列式值，〇 2为背景图重力值的方差；通过调整相关尺 度9，结合高斯模型，来逼近Z的相关矩阵Ri，j，构造和描述相关模型。
【专利摘要】本发明涉及一种基于泛克里金法的重力场插值方法，其主要技术特点是包括以下步骤：针对空间中的重力异常模型建立回归模型；针对空间中的重力异常模型建立相关模型；计算回归模型和相关模型的权重系数，代入重力异常模型，便可获得待插值位置上的重力异常空间分布表达式，再将待插值位置代入重力异常模型，便可计算出待插值位置的重力值。本发明采用逐步回归法对重力场的非平稳性及相关性分别建模描述，使用多元逐步回归法，结合拉格朗日乘数法，完整的描述了重力场的非平稳性，保证了插值结果的无偏性，具有较高的插值精度。
【IPC分类】G01C21/20
【公开号】CN105043390
【申请号】CN201510371479
【发明人】张辰, 周贤高, 于浩, 赵晶晴, 欧阳永忠
【申请人】中国船舶重工集团公司第七0七研究所
【公开日】2015年11月11日
【申请日】2015年6月29日