CN105043390A - 基于泛克里金法的重力场插值方法 - Google Patents

Abstract

本发明涉及一种基于泛克里金法的重力场插值方法，其主要技术特点是包括以下步骤：针对空间中的重力异常模型建立回归模型；针对空间中的重力异常模型建立相关模型；计算回归模型和相关模型的权重系数，代入重力异常模型，便可获得待插值位置上的重力异常空间分布表达式，再将待插值位置代入重力异常模型，便可计算出待插值位置的重力值。本发明采用逐步回归法对重力场的非平稳性及相关性分别建模描述，使用多元逐步回归法，结合拉格朗日乘数法，完整的描述了重力场的非平稳性，保证了插值结果的无偏性，具有较高的插值精度。

Description

基于泛克里金法的重力场插值方法

技术领域

本发明属于惯性/重力组合导航系统技术领域，具体涉及一种基于泛克里金法的重力场插值方法。

背景技术

导航用背景图是重力匹配导航的基础，背景图的精度及分辨率直接影响匹配精度。目前大部分海域的重力异常信息，都采用卫星反演的手段获得。但该方法的分辨率较低，最高只有1’*1’。在匹配过程中，为了保证重力异常的测量精度，载体必须以低速航行，从实用性的角度考虑，匹配时间不宜过长，在长度有限的航迹上，卫星反演背景图包含的信息量太少，无法满足匹配定位需求，因此必须对其进行插值处理。

目前常用的背景图插值算法包括距离倒数插值法及三次样条插值法。距离倒数插值法首先搜索最靠近待插值点的四个背景图采样点(东北、西北、东南、西南方向)，之后计算待插值点到四个采样点的距离，以距离的倒数作为权重，取四个采样点重力异常值的加权和，作为待插值点的重力异常。三次样条插值法对采样点之间的空间区域，采用连续、二阶可导，且导数连续的三次函数拟合重力异常分布，其插值结果更加光滑。这两种算法计算量较小，但均没有针对重力场的特点进行针对性的建模描述。

由于地质结构千变万化，重力异常随之剧烈起伏，这会导致在不同区域内，重力异常的均值和方差并非常数，也即重力场是一个非平稳场，上述两种算法人为设定非平稳的重力场在空间上按照特定的模型分布，因此其插值结果不可避免的存在偏差。此外，地质结构在空间上是存在连续性的，这种连续性表现在重力异常值上，就是重力异常在空间上存在相关性，而这种相关性在不同空间范围内，又是各异的，上述两种算法均没有针对重力场的相关性进行建模描述，因此无法准确表现出地质结构的连续性。最后，上述两种插值算法，其插值的权重系数仅由间隔距离计算得到，亦不是最优的。综上所述，上述两种插值算法，既无法达到无偏最优的插值目标，又没有描述重力场的相关性，因此插值精度较低。

发明内容

本发明的目的在于克服现有技术的不足，提供一种设计合理且插值精度高的基于泛克里金法的重力场插值方法。

本发明解决其技术问题是采取以下技术方案实现的：

一种基于泛克里金法的重力场插值方法，于包括以下步骤：

步骤1、针对空间中的重力异常模型建立回归模型；

步骤2、针对空间中的重力异常模型建立相关模型；

步骤3、计算回归模型和相关模型的权重系数，代入重力异常模型，便可获得待插值位置上的重力异常空间分布表达式，再将待插值位置代入重力异常模型，便可计算出待插值位置的重力值。

而且，所述重力异常模型为：

yx＝f1x…fpx*βi,1+zx

式中，x为空间位置，是一个1*2维位置矩阵，第一列保存经度，第二列保存纬度；yx为模型解算出的重力值。

而且，所述步骤1建立回归模型的方法包括以下步骤：

(1)建立由全部因子x1,x2,…,xn组成的自变量集合，设自变量集合包含n个因子，并设定因子被选入回归模型的显著性水平为α1，被剔除出回归模型的显著性水平为α2，0<α1≤α2<1；

(2)在背景图上抽取m组经纬度与重力值数据，作为回归模型的样本；

(3)从全部因子中，抽取出任意两个因子xi,xj，依据m对样本值，利用最小二乘法计算其参数向量τi,j1，得到用这两个因子描述的空间重力场yk1＝xi,xj*τi,j1，并计算这个模型的剩余平方和Si,j1＝k＝1myk-yk12；

(4)重复执行步骤(3)，遍历全部因子，计算任意两个因子所组成的模型的剩余平方和Si,j1，组成离差矩阵S1；

(5)计算各个因子在离差矩阵中的贡献Vj1＝(Sj,y1)2Sj,j1，设贡献最大的因子编号为k1，对其做显著性水平为α1的显著性F检验：

F＝Vj1SE1(M-1-1)

其中，SE1＝k＝1myk-y2-Vj1,y＝k＝1mykm；

当F>Fα11,m-1-1时，该因子选入回归模型；

(6)在k1因子选入回归模型后，调整离差矩阵S2：

按照步骤(5)的方法重新计算剩余因子在新离差矩阵中的贡献，并对贡献最大的因子k2做显著性水平为α1，自由度为m-2-1的显著性F检验，作为其是否被引入回归模型的判据；

(7)在选入两个因子后，首先按照步骤(5)的方法，将离差矩阵S2调整为S3，在选入新因子之前，先对已经选入的两个因子，按照被剔除出回归模型的显著性水平α2做显著性F检验，以判断是否需要将其从回归模型中剔除，F>Fα21,m-3-1时保留，F<Fα21,m-3-1时剔除，之后再引入新的因子；

(8)按照步骤(7)筛选的方法，对自变量集合中的因子进行遍历，直到新的因子不再被引入回归模型，同时回归模型内的因子也不再被剔除时，即可获得重力场的回归模型。

而且，所述步骤2建立相关模型的方法包括以下步骤：

(1)对重力矩阵Y进行平差处理，得到Y'；

(2)计算Y'的相关矩阵，对于空间中任意两个采样点位置xi与xj，Y'的相关矩阵为：

EY'xi*Y'xj＝Ri,j＝σ2Rθ,xi,xj

上式中，σ2为背景图重力值的方差；

(3)使用高斯相关模型对Z进行描述，高斯相关模型的相关函数形式如下：

C0＝σ2＝C1

Ch＝C1e-hθ2

Chh→∞0

(4)计算相关尺度θ的最优解：

minθ＝cθ≡R1mσ2

其中，R为矩阵Ri,j所对应的行列式值，σ2为背景图重力值的方差；通过调整相关尺度θ，结合高斯模型，来逼近Z的相关矩阵Ri,j，构造和描述相关模型。

本发明的优点和积极效果是：

1、本发明采用逐步回归法对重力场的非平稳性及相关性分别建模描述，使用多元逐步回归法，结合拉格朗日乘数法，完整的描述了重力场的非平稳性，保证了插值结果的无偏性，具有较高的插值精度。

2、本发明拉格朗日乘数法，保证了插值结果的最优性。

3、本发明采用回归模型的引入方式，使得构建的重力异常模型完整表述重力场的相关性。

附图说明

图1为逐步回归法的计算步骤；

图2为某重力异常序列的相关部分；

图3为重力异常相关部分的相关函数；

图4为高斯模型协方差函数；

图5为4#测线插值结果；

图6为4#测线采用三次样条插值误差变异函数；

图7为4#测线采用本专利算法插值误差变异函数；

具体实施方式

以下结合附图对本发明实施例做进一步详述：

一种基于泛克里金法的重力场插值方法，是基于泛克里金法来设计的。泛克里金法是一种针对非平稳场的插值算法，该算法是建立在变异函数理论分析基础上，对有限区域内变量取值进行无偏最优估计的一种方法。本发明的算法假定空间中的重力异常符合如下模型：

yx＝f1x…fpx*βi,1+zx(1)

式(1)中，x为空间位置，是一个1*2维位置矩阵，第一列保存经度，第二列保存纬度；yx为模型解算出的重力值。

将m个背景图采样点位置x1…xm及与之相对应的m个重力异常值yx1…yxm代入式(1)，可以得到：

yx1＝f1x1…fpx1*βi,1+zx1

…(2)

yxm＝f1xm…fpm*βi,1+zxm

并且设定：X＝x1…xmT，Y＝yx1…yxmT,Z＝zx1…zxmT，β＝βi,1，

f1x1…fpx1

F＝…

f1xm…fpxm

于是，式(2)可简写为：Y＝F*β+Z(3)

一种基于泛克里金法的重力场插值方法，包括以下步骤：

步骤1、针对空间中的重力异常模型建立回归模型

由于重力场是非平稳场，在不同的区域内，重力异常的均值并非常数。针对这个特点，本算法采用以经度、纬度为自变量，重力异常值为因变量的多项式，来拟合回归模型，并利用逐步回归法，确定组成回归模型的多项式因子(下面简称为因子)及其阶数。

本步骤是对公式(1)中的f1x…fpx建模，它由p个多项式组成，多项式的自变量为经度、纬度，因变量为重力异常，每个多项式被称为一个因子，这部分模型被称为回归模型，主要用来描述重力场的非平稳性。

为了保证建模的精度，在建模过程中使用了多元逐步回归法，其本质思想就是将采样点的空间位置X及其所对应的重力异常值Y，视为回归模型的样本，根据该样本，从全部因子(也即多项式)集合中，筛选出对重力异常有显著影响的因子(也即多项式)，构成回归模型。多元逐步回归法的引入，保证了被选入回归模型的多项式，其阶数、结构形式完全由重力场空间分布决定，结合步骤三中采用的拉格朗日乘数法计算出的权重系数β，即可保证回归模型对重力场估计的无偏性，进而保证公式(1)中剩余的重力异常信号z＝Y-Fβ均值为零，使得用z来描述重力场自相关性成为可能。多元逐步回归法与拉格朗日乘数法相结合，有效地将“非平稳场建模问题”转化为“广义平稳随机过程建模问题”。

建立回归模型的具体计算步骤如下：

(1)建立由全部因子x1,x2,…,xn组成的自变量集合，设自变量集合包含n个因子(n>p)；并设定因子被选入回归模型的显著性水平为α1，被剔除出回归模型的显著性水平为α2，0<α1≤α2<1。

(2)在背景图上抽取m组数据(经纬度与重力值，如公式(3)所示)，作为回归模型的样本。

(3)从全部n个因子中，抽取出任意两个因子xi,xj(1≤i≤n,1≤j≤n)，假定空间重力场可以用选取出的两个因子描述，依据m对样本值，利用最小二乘法计算其参数向量τi,j1，得到用这两个因子描述的空间重力场yk1＝xi,xj*τi,j1，并计算这个模型的剩余平方和Si,j1＝k＝1myk-yk12。

(4)按照步骤(3)的方法，遍历全部因子，计算任意两个因子所组成的模型的剩余平方和Si,j1(1≤i≤n,1≤j≤n)，组成离差矩阵S1：

S1＝S1,11…S1,n1S1,y1…………Sn,11…Sn,n1Sn,y1

离差矩阵的最后一列表征了以单一因子描述重力场时的绝对误差，误差越小，对应的因子对重力场空间分布的影响越显著。

(5)为了便于比对不同因子对重力场影响的显著程度，计算各个因子在离差矩阵中的贡献Vj1＝(Sj,y1)2Sj,j1，设贡献最大的因子编号为k1，对其做显著性水平为α1的显著性F检验：

F＝Vj1SE1(M-1-1)

其中，SE1＝k＝1myk-y2-Vj1,y＝k＝1mykm。

当F>Fα11,m-1-1时，该因子选入回归模型，步骤(3)-(5)为第一轮筛选。

(6)在k1因子选入回归模型后，需要调整离差矩阵：

S2＝S1,12…S1,n2S1,y2…………Sn,12…Sn,n2Sn,y2

其中，Si,j2＝Sk1,j1Sk1,k11当i＝k1,j≠k1Si,j1-(Si,k11*Sk1,j1)Sk1,k11

当i≠k1,j≠k11Sk1,k11当i＝j＝k1(-Si,k11)Sk1,k11当i≠k1,j＝k1

按照步骤(5)的方法，重新计算剩余因子在新离差矩阵中的贡献，并对贡献最大的因子k2做显著性水平为α1，自由度为m-2-1的显著性F检验，作为其是否被引入回归模型的判据，此为第二轮筛选。

(7)在选入两个因子后，首先按照步骤(5)的方法，将离差矩阵S2调整为S3，进行第三轮筛选。在选入新因子之前，先对已经选入的两个因子，按照被剔除出回归模型的显著性水平α2做显著性F检验，以判断是否需要将其从回归模型中剔除(F>Fα21,m-3-1时保留，F<Fα21,m-3-1时剔除)，之后再引入新的因子。

(8)按照第三轮筛选的方法，对自变量集合中的因子进行遍历，直到新的因子不再被引入回归模型，同时回归模型内的因子也不再被剔除时，即可获得重力场的多项式回归模型。逐步回归法的计算步骤如图1所示。

步骤2、针对空间中的重力异常模型建立相关模型

本步骤是对公式(1)中的z(x)建模，这部分模型被称为相关模型，主要用来描述重力场的相关性。

由于在计算插值系数之前，无法获得z(x)的准确数值，因此需要借助重力矩阵Y，对重力场的自相关性进行建模描述。为了消除均值变化对相关模型精度的影响，建模时需要首先对重力矩阵Y规范化，进行平差处理，使得插值范围内的重力异常均值为0。规范化后的重力矩阵Y'，可以用广义平稳随机过程来描述，可以代替z(x)，来表征重力场的相关性。

相关模型建模的本质思想，就是将Y'的相关函数视为相关模型的样本，通过选取合适的相关模型Ch,θ，来满足下列约束条件：

Ch,θ＝CovY'x+h,Y'x＝E{Y'x+h-EY'x+h*Y'x-EY'x}＝EY'x+h*Y'x

可以看出，采用相关模型表征重力场的相关性，不仅描述了采样点与待插值点之间的相对位置及相关关系，还考虑了采样点的空间分布结构特点，完整地描述了重力场的自相关性。为了行文方便，下文中，Z与Y'依推导需求出现，二者可以相互替代。

相关模型的建立步骤如下所示：

(1)对重力矩阵Y进行平差处理，得到Y'。

(2)计算Y'的相关矩阵，它是相关模型的一个样本估计。对于空间中任意两个采样点位置xi与xj，Y'的相关矩阵如下式所示：

EY'xi*Y'xj＝Ri,j＝σ2Rθ,xi,xj(4)

上式中，σ2为背景图重力值的方差。

(3)选取合适的相关模型。

在获得Y'的相关矩阵(公式(4))后，需要寻找合适的相关模型对其进行拟合。一般而言，对于一个强度在空间上连续可导的矢量场，其相关模型近似于指数形式。相反地，其相关模型近似于线性形式。当载体沿东西方向航行时，某航迹上重力异常的相关部分，如图2所示。经过计算，其相关函数，如图3所示。从图3可以看出，由于重力场在空间上连续可导，因此Z的相关矩阵接近于指数相关模型，本算法中，使用高斯相关模型对Z进行描述。高斯相关模型的相关函数形式如下：

C0＝σ2＝C1

Ch＝C1e-hθ2

Chh→∞0

其相关函数如图4所示(图中σ2＝1，θ＝1.7308h)。

从图3可以看出，随着间隔距离增大，Y'的相关性都会减小。在高斯相关模型中，包含一个参数θ，称为相关尺度：θ越大，Y'的相关矩阵取值随距离增大减小得越快，Y'越接近随机信号；相反，θ越小，Y'越接近常值。

(4)计算相关尺度θ。

结合高斯相关模型的样本估计(公式(4))，依据最优原则，利用最小二乘法，θ的最优解可以表示为下式：

minθ＝cθ≡R1mσ2

其中，R为矩阵Ri,j所对应的行列式值，σ2为背景图重力值的方差。通过调整相关尺度θ，结合高斯模型，来逼近Z的相关矩阵Ri,j，构造和描述相关部分。

步骤3、计算回归模型和相关模型的权重系数，代入公式(1)，即可获得待插值位置上的重力异常空间分布表达式，再将待插值位置代入公式(1)，即可计算出待插值位置的重力值。

本步骤在建立了重力场的空间分布模型后，该算法基于最优化原理，利用拉格朗日乘数法计算确定性部分的权重βi,1及相关部分的权重γ，求解权重矩阵。设得到的权重矩阵为c，为一个1*m维的权重矩阵，则插值算法在待测点处的估计值为：

yx＝cTY(5)

插值算法的目的，就是求解出矩阵c的取值，进而计算出回归模型的权重βi,1及相关模型的权重γ，以满足估计的无偏和最优原则。

由公式(3)、(5)可得，插值结果的误差为：

yx-yx＝cTY-yx＝cTFβ+Z-fxTβ+z＝cTZ-z+FTc-fxTβ

由于Z为广义平稳随机过程，因此Z和z均值为相等的常数，且对于位于地球上的任意位置的卫星反演图(任意β值)，要满足估计是无偏的(即上式恒等于0)，必须保证：

FTc-fx＝0(6)

由公式(3)、(5)可得，插值误差的方差为：

式中rx＝Ezxzxj＝σ2Rθ,x,xj

要使得估计是最优的，必须使得取值最小，结合公式(6)的约束条件，采用拉格朗日乘数法，建立拉格朗日函数：

Lc,λ＝σ21+cTRc-2cTr-λTFTc-f

式中，c和λ为待求解量。对c求导，并令其等于0，得到：

Lc'c,λ＝2σ2Rc-r-Fλ＝0(7)

结合公式(6)和(7)，可以解得：

λ1＝FTR-1F-1FTR-1r-fc＝R-1r-Fλ1

其中，λ1＝-λ2σ2

将解得的λ和c代入公式(5)，可以获得待插值位置重力场的解析表达式：

yx＝r-Fλ1TR-1Y＝rTR-1Y-FTR-1r-fTFTR-1F-1FTR-1Y(8)

对于公式(8)，按照公式(1)的结构提取合并，可以得到：βi,1

＝FTR-1F-1FTR-1Y(9)

将βi,1代入公式(8)，整理后可以得到：

yx＝rTR-1Y-FTR-1r-fTβi,1＝f(x)Tβi,1+r(x)Tγ(10)

式中，γ＝R-1(Y-F)βi,1(11)

观察公式(10)和公式(1)，f(x)和r(x)分别对应于事先拟合出的回归模型和相关模型。一旦确定了背景图的经纬度范围，权重矩阵βi,1及γ也可计算得到(公式(9)、公式(11))。因此，将待插值位置x代入公式(10)，即可得到待插值位置的重力异常值。仔细分析公式(11)，在γ＝R-1(Y-F)βi,1中，Ri,j＝Ezxizxj＝σ2Rθ,xi,xj描述了重力场在任意两个采样点xi和xj之间的空间分布特点。而相关矩阵rx＝Ezxzxj＝σ2Rθ,x,xj则表征了待插值位置x与采样点之间的重力场空间分布特点。因此，该算法在满足最优无偏的基础上，综合考虑了重力场空间分布特点以及采样点之间的相互作用影响，完整地描述了空间重力异常分布。

下面给出一个仿真验证，对本发明的效果进行验证。

利用实测重力异常数据中的1#、2#、3#、5#、6#、7#测线，构造重力背景图；分别用三次样条插值法和本发明提出的算法，在构造出的背景图上插值计算4#测线重力异常，并与4#测线实测数据相比对，如图5所示。图中曲线分别代表实测重力异常、三次样条插值结果及本专利提出插值算法的插值结果。

无论哪种插值算法，其插值结果均存在误差：如果插值算法专门针对重力场的相关性进行了建模描述，则插值误差只包含常值信号和随机信号，这两种信号的变异函数均为常数；如果插值算法在建模过程中忽略了重力场的相关性，插值误差就包含随机游走信号，其变异函数取值就不是常数。因此，可以通过计算插值误差变异函数的方法，来评判插值算法的性能：如果误差的变异函数取值为常数，则可以判定所得到的模型准确反映了重力场的相关性，该算法的性能较好。

两种插值算法的误差，其变异函数分别如图6及图7所示。由变异函数可以看出，本专利算法的插值误差变异函数为常值，所建立起来的模型准确的描述了重力场的相关性。而三次样条插值误差中存在明显的随机游走误差，效果较差。

两种插值算法的精度为：

	本文提出的插值算法(mgal)	三次样条插值(mgal)
			误差均值	-0.0043	0.5316
误差峰峰值	3.0045	3.3046
			误差标准差	0.564	0.6298

由上表看出，本专利算法的插值结果满足无偏最优的要求。

本发明不仅局限于上述的具体实施方式，只要插值算法特别针对重力场的相关性进行了建模描述，对其确定性部分和相关部分分别采用不同的方式建模，并应用最优化方法确定两部分的权重系数，均落入本专利的保护范围内。

Claims (4)

1.一种基于泛克里金法的重力场插值方法，其特征在于包括以下步骤：

步骤1、针对空间中的重力异常模型建立回归模型；

步骤2、针对空间中的重力异常模型建立相关模型；

步骤3、计算回归模型和相关模型的权重系数，代入重力异常模型，便可获得待插值位置上的重力异常空间分布表达式，再将待插值位置代入重力异常模型，便可计算出待插值位置的重力值。

2.根据权利要求1所述的基于泛克里金法的重力场插值方法，其特征在于：所述重力异常模型为：

yx＝f1x…fpx*βi,1+zx

式中，x为空间位置，是一个1*2维位置矩阵，第一列保存经度，第二列保存纬度；yx为模型解算出的重力值。

3.根据权利要求1或2所述的基于泛克里金法的重力场插值方法，其特征在于：所述步骤1建立回归模型的方法包括以下步骤：

(1)建立由全部因子x1,x2,…,xn组成的自变量集合，设自变量集合包含n个因子，并设定因子被选入回归模型的显著性水平为α1，被剔除出回归模型的显著性水平为α2，0<α1≤α2<1；

(2)在背景图上抽取m组经纬度与重力值数据，作为回归模型的样本；

(3)从全部因子中，抽取出任意两个因子xi,xj，依据m对样本值，利用最小二乘法计算其参数向量τi,j1，得到用这两个因子描述的空间重力场yk1＝xi,xj*τi,j1，并计算这个模型的剩余平方和Si,j1＝k＝1myk-yk12；

(4)重复执行步骤(3)，遍历全部因子，计算任意两个因子所组成的模型的剩余平方和Si,j1，组成离差矩阵S1；

(5)计算各个因子在离差矩阵中的贡献Vj1＝(Sj,y1)2Sj,j1，设贡献最大的因子编号为k1，对其做显著性水平为α1的显著性F检验：

F＝Vj1SE1(M-1-1)

其中，SE1＝k＝1myk-y2-Vj1,y＝k＝1mykm；

当F>Fα11,m-1-1时，该因子选入回归模型；

(6)在k1因子选入回归模型后，调整离差矩阵S2：

按照步骤(5)的方法重新计算剩余因子在新离差矩阵中的贡献，并对贡献最大的因子k2做显著性水平为α1，自由度为m-2-1的显著性F检验，作为其是否被引入回归模型的判据；

(7)在选入两个因子后，首先按照步骤(5)的方法，将离差矩阵S2调整为S3，在选入新因子之前，先对已经选入的两个因子，按照被剔除出回归模型的显著性水平α2做显著性F检验，以判断是否需要将其从回归模型中剔除，F>Fα21,m-3-1时保留，F<Fα21,m-3-1时剔除，之后再引入新的因子；

(8)按照步骤(7)筛选的方法，对自变量集合中的因子进行遍历，直到新的因子不再被引入回归模型，同时回归模型内的因子也不再被剔除时，即可获得重力场的回归模型。

4.根据权利要求1或2所述的基于泛克里金法的重力场插值方法，其特征在于：所述步骤2建立相关模型的方法包括以下步骤：

(1)对重力矩阵Y进行平差处理，得到Y'；

(2)计算Y'的相关矩阵，对于空间中任意两个采样点位置xi与xj，Y'的相关矩阵为：

EY'xi*Y'xj＝Ri,j＝σ2Rθ,xi,xj

上式中，σ2为背景图重力值的方差；

(3)使用高斯相关模型对Z进行描述，高斯相关模型的相关函数形式如下：

C0＝σ2＝C1

Ch＝C1e-hθ2

Chh→∞0

(4)计算相关尺度θ的最优解：

minθ＝cθ≡R1mσ2

其中，R为矩阵Ri,j所对应的行列式值，σ2为背景图重力值的方差；通过调整相关尺度θ，结合高斯模型，来逼近Z的相关矩阵Ri,j，构造和描述相关模型。