基于重力异常时变的重力基准图时变修正方法及系统与流程

本发明涉及属于控制理论领域，主要针对现有重力基准图导航匹配精度随时间推移会降低的问题，而提出了一种基于重力异常时变的重力基准图时变修正方法及系统。

背景技术：

直观认为，除去地震等因素的影响外，地球重力场是静止不变的，并不随时间发生变化，这也是重力值可以作为一种无源辅助导航匹配要素的关键所在和基本出发点。

随着重力测量技术的进步以及对地球重力场的深入研究，地球物质分布变化引起重力场随时间发生变化的现象被发现。同时，随着重力卫星的相继发射，研究时变重力场成为现实。图1所示为同一区域2002年和2006年的局部重力异常示意图，可以明显的看到该区域的重力异常值发生了变化。

目前，现有技术可以精确测量全球范围的重力数据，精度较高，但是不具有时变性，只有通过多次测量才可以发现重力在变化。那么，通过测量得到的重力数据去绘制的重力基准图，也不具有时变性，在导航时且在初始导航效果良好的情况下，由于时效性缺乏，制备重力基准图后，若较长时间不对其进行时变修正，当再次进行导航匹配，就会出现较大的误差，定位精度明显降低。

技术实现要素：

本发明的目的是提供一种基于重力异常时变的重力基准图时变修正方法及系统，通过造成重力异常变化的影响因素分类抽象得到重力异常拟合公式，并且结合消除南北条带误差的重力异常变化公式，得到重力异常时变公式，并将其用于重力基准图的时变修正，降低重力基准图的误匹配率。

为实现上述目的，本发明提供了如下方案：

一种基于重力异常时变的重力基准图时变修正方法，所述重力基准图时变修正方法包括：

根据grace卫星数据解算获得重力数据，并根据所述重力数据，得到月重力场模型；

选取任意一个月重力场模型为背景重力场模型，并将所述月重力场模型的表达式和所述背景重力场模型的表达式进行差值计算，得到所述背景重力场下的重力异常变化公式；

对所述重力异常变化公式进行空间平滑处理，确定平滑后的重力异常变化公式；

对重力场时变影响因素进行分类拟合，得到重力异常拟合公式；

根据所述平滑后的重力异常变化公式和所述重力异常拟合公式，确定重力异常长期变化率；

根据所述重力异常长期变化率，确定当前时刻的重力异常时变公式；

根据所述重力异常时变公式，计算当前时刻的重力异常值；

根据所述重力异常值，对制备重力基准图时所采用的原始数据进行时变修正，得到时变修正重力基准图。

可选的，所述对所述重力异常变化公式进行空间平滑处理，确定平滑后的重力异常变化公式，具体包括：

将各向同性的高斯平滑滤波函数代入所述重力异常变化公式中，得到平滑后的重力异常变化公式，以消除利用grace卫星数据解算重力数据时所存在的南北条带误差干扰。

可选的，在执行将各向同性的高斯平滑滤波函数代入所述重力异常变化公式之前，还包括：

确定重力场模型的截断阶数；

根据所述截断阶数，确定高斯平滑半径。

可选的，所述平滑后的重力异常变化公式为：

其中，f为万有引力常数，m为地球质量，r是观测点到地球质心的距离，ae为地球椭球长半径，λ、分别对应于观测点的经度和纬度，l，m为时变重力场模型球谐展开的l阶、m次，wl为高斯平滑滤波函数，为月重力场模型与背景重力场模型的球谐系数之差，为完全规格化的缔合勒让德多项式。

可选的，所述重力异常拟合公式为

其中，δt为十进制月重力场模型的时间与背景重力场模型的时间的差值，a为常值，近似于直线拟合的常值，b为重力异常长期变化率，i＝1代表年周期变化项，i＝2代表半年周期变化项，则年变化的振幅为i＝3表示与s2半日波相关的161天周期项，ωi为具体拟合时的角变化率，正弦函数和余弦函数拟合一条曲线时的变化率，ξ表示噪声和误差影响，属于不确定项。

可选的，所述重力异常时变公式为

其中，t0＝2003.0为基准时间，为基准时间的重力异常值，为t时刻的重力异常值。

可选的，所述根据所述重力异常值，对制备重力基准图时所采用的原始数据进行时变修正，得到时变修正重力基准图，具体包括：

根据所述重力异常值，对制备重力基准图所用的原始数据进行时变修正，得到时变修正数据；

根据所述时变修正数据，利用插值算法进行网格化制图，得到时变修正重力基准图。

本发明还提供了一种基于重力异常时变的重力基准图时变修正系统，所述重力基准图时变修正系统包括：

月重力场模型得到模块，用于根据grace卫星数据解算获得重力数据，并根据所述重力数据，得到月重力场模型；

重力异常变化公式得到模块，用于选取任意一个月重力场模型为背景重力场模型，并将所述月重力场模型的表达式和所述背景重力场模型的表达式进行差值计算，得到所述背景重力场下的重力异常变化公式；

平滑后的重力异常变化公式确定模块，用于对所述重力异常变化公式进行空间平滑处理，确定平滑后的重力异常变化公式；

重力异常拟合公式得到模块，用于对重力场时变影响因素进行分类拟合，得到重力异常拟合公式；

重力异常长期变化率确定模块，用于根据所述平滑后的重力异常变化公式和所述重力异常拟合公式，确定重力异常长期变化率；

重力异常时变公式确定模块，用于根据所述重力异常长期变化率，确定当前时刻的重力异常时变公式；

当前时刻重力异常值计算模块，用于根据所述重力异常时变公式，计算当前时刻的重力异常值；

时变修正重力基准得到模块，用于根据所述重力异常值，对制备重力基准图时所采用的原始数据进行时变修正，得到时变修正重力基准图。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明提供了一种基于重力异常时变的重力基准图时变修正方法及系统，包括：首先由grace卫星数据解算获得重力数据计算得到月重力场模型，选取任意一个月重力场模型为背景重力场模型，并将月重力场模型的表达式和背景重力场模型的表达式进行差值计算得到背景重力场下的重力异常变化公式，再对重力异常变化公式进行空间平滑处理，以消除利用grace卫星数据解算重力数据时所存在的南北条带误差干扰；然后对重力场时变影响因素进行分类拟合得到重力异常拟合公式，并根据平滑后的重力异常变化公式和重力异常拟合公式确定重力异常长期变化率，实现对重力异常变化的预测；最后根据重力异常长期变化率确定当前时刻的重力异常时变公式，计算当前时刻的重力异常值，对制备重力基准图时所采用的原始数据进行时变修正，得到时变修正重力基准图，降低重力基准图的误匹配率。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明2002年和2006年的局部重力异常图；

图2为本发明实施例重力基准图时变修正方法的流程示意图；

图3为本发明受南北条带误差影响的地球重力场变化图；

图4为本发明重力异常阶方差与截断阶数的关系图；

图5为本发明重力非潮汐变化示意图；

图6为本发明实施例重力基准图时变修正系统的结构示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的目的是提供一种基于重力异常时变的重力基准图时变修正方法及系统，以用于重力基准图的时变修正，降低重力基准图的误匹配率。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

本发明创新部分包括，第一，在重力异常拟合公式的求解中，对重力异常产生变化的影响因素实现数学拟合是关键。造成重力变化的因素繁多，且机理各不相同，因此对其实现拟合，要按影响的原理不同进行分类拟合。首先对每一类的机理进行合理分析，抽象出数学拟合的表达式，其次利用已有的重力数据资料，求解得到未知参数，进而得到拟合表达式，最后将各项求和，即得到重力异常拟合公式。

第二，得到重力异常时变公式，不仅需要重力异常拟合公式和重力异常时变公式进行组合求解，还需要对卫星数据求解的重力异常时变公式进行滤波，已消除南北条带误差的影响，进而得到重力异常长期变化率。利用重力异常长期变化率，结合背景重力场模型的重力异常值，即可推导得到重力异常时变公式。

图2为本发明实施例重力基准图时变修正方法的流程示意图，如图2所示，本发明实施例提供的重力基准图时变修正方法包括以下几个步骤：

步骤101：根据grace卫星数据解算获得重力数据，并根据所述重力数据，得到月重力场模型。

步骤102：选取任意一个月重力场模型为背景重力场模型，并将所述月重力场模型的表达式和所述背景重力场模型的表达式进行差值计算，得到所述背景重力场下的重力异常变化公式。

步骤103：对所述重力异常变化公式进行空间平滑处理，确定平滑后的重力异常变化公式。

步骤104：对重力场时变影响因素进行分类拟合，得到重力异常拟合公式。

步骤105：根据所述平滑后的重力异常变化公式和所述重力异常拟合公式，确定重力异常长期变化率。

步骤106：根据所述重力异常长期变化率，确定当前时刻的重力异常时变公式。

步骤107：根据所述重力异常时变公式，计算当前时刻的重力异常值。

步骤108：根据所述重力异常值，对制备重力基准图时所采用的原始数据进行时变修正，得到时变修正重力基准图。

步骤101具体包括：由grace卫星数据解算获得的重力数据计算可得到月重力场模型。

某点处的重力异常δg可由公式(1)表达：δg＝g-g'(1)。

在公式(1)中，g为实际重力值，g'为正常重力值。

进而由布隆斯公式可得:

在公式(2)中，t为扰动位，n为大地水准面高度。

由于高程h的计算是沿法线方向，其正方向即是外法线的方向，由此可得到物理大地测量学的基本微分方程，即：

由基本微分方程可以看出，要对重力异常进行计算，首先要对扰动位进行研究。为简化计算，假设正常椭球的旋转轴、中心分别与实际地球的旋转轴、质心重合，且两者具有相同的旋转角速度，在此条件下对重力扰动位的计算进行讨论。

由假设条件可得，对于正常椭球，其离心力位等于实际地球的离心力位，则可得其扰动位为：t＝u-vh(4)。

在公式(4)中，u为地球引力位，vh为正常引力位，两者均是调和函数，可用球谐函数级数表示为：

在公式(5)和公式(6)中，f为万有引力常数，m为地球质量，ae为地球椭球长半径，为l阶m次完全规格化的球谐函数系数，为完全规格化的缔合勒让德多项式，r是观测点到地球质心的距离，λ、分别对应于观测点的经度和纬度，为正常引力位系数，常值，与正常椭球参数有关，为2l阶0次完全规格化缔合勒让德多项式。

其中，正常引力位的二阶带谐系数表达式为：

在公式(7)中，e为椭球第一偏心率，e'为椭球第二偏心率，为地球的自转角速度。由正常椭球参数，取即可求得椭球带谐系数及其引力位。

将公式(5)和公式(6)式代入公式(4)，扰动位球谐函数的级数展开式为：

式中，

同时，由球近似下的边界条件为：

将公式(8)代入公式(2)中，可得到重力异常的球谐表达式，即月重力场模型的表达式，为：

在公式(10)中，f为万有引力常数，m为地球质量，r是观测点到地球质心的距离，ae为地球椭球长半径，l，m为时变重力场模型球谐展开的l阶、m次，和为完全规格化的球谐函数系数，为l阶m次完全规格化的球谐函数系数，为完全规格化的缔合勒让德多项式，λ、分别对应于观测点的经度和纬度。

步骤102具体包括：选取某月重力场模型为背景重力场模型，并将月重力场模型的表达式与背景重力场模型的表达式进行差值计算，即可得到背景重力场下的重力异常变化公式，为：

其中，为月重力场模型与背景重力场模型的球谐系数之差。

步骤103具体包括：

研究表明，利用grace卫星数据进行计算时，会存在严重的南北条带误差干扰，如图3所示的受南北条带误差影响时的地球重力场变化图。因此，为提取有效的重力数据，需要分析其变化特征，进行空间平滑滤波处理。

为降低高阶项误差的影响，准确的分析地球重力场时变情况，本发明采用各向同性的高斯平滑滤波函数作为空间平滑处理的手段。将高斯平滑滤波函数代入公式(11)中，可得平滑后的重力异常变化公式：

其中，为月重力场模型与背景重力场模型的球谐系数之差，wl为高斯平滑滤波函数。各向同性的高斯平滑滤波函数为：

在公式(13)中，γ为地表点与之间的角距，满足关系式：

其递推公式为：

在将高斯平滑滤波函数代入式(11)之前，还需要：

s1：确定重力场模型的截断阶数。

利用grace卫星数据恢复重力场模型的阶次是有限制的，当采用的球谐系数超过90阶时，重力场模型的误差会迅速增大。

因此，要进行阶方差的计算，在误差允许的范围内选取合适的截断阶数。阶方差的计算公式为

对2002年04月到2006年12月间共53个月的重力数据进行计算，得到了不同阶数时重力变化的分布情况，清晰的反映了重力异常阶方差与截断阶数的关系，为计算时选取截断阶数提供了参考，如图4所示。计算所用重力场模型为csr提供的grace月重力场模型。

s2：高斯平滑半径的确定。

选取合理的高斯平滑半径，可以有效减小干扰，对时变重力场的计算有着重要影响，所以应该选取合适的高斯平滑半径b。又因为当高斯平滑半径b过大时，重力场可变信号损失严重。因此，对于高斯平滑半径的选取应该在截断阶数确定后确定。

研究表明，可以根据重力场模型的阶数与高斯平滑半径b的关系式，计算高斯平滑半径b；重力场模型的阶数与高斯平滑半径b的关系式为n为时变重力场模型的阶数。

步骤104具体包括：

研究表明，造成重力随时间变化的因素较多，主要是地球质量的分布变化引起的，大致可分为潮汐变化和非潮汐变化，影响因素包括大气质量变化、海水潮汐、水文状况变化以及地震、火山喷发等自然因素。重力非潮汐变化示意图如图5所示。

以上影响因素，除去可进行周期性观测的外，仍存在一些不可控因素，无法实现统一的数学表达。因此对于重力场时变影响因素进行分析量化时，不可利用简单的线性实现拟合，要结合其影响因素的具体性质按类进行划分拟合。本发明把造成重力场时变影响因素分为3类进行拟合，具体为：

1)线性拟合项，主要反映的是重力场随时间的长期缓慢变化，具体拟合为常值，表示重力随时间的长期变化率。

2)周期项，主要包括造成重力场变化的可进行周期性观测的影响因素，利用正余弦函数进行组合，实现拟合。

3)随机项，包括误差、噪声等无法进行预测和表达的影响因素，由ξ表示。

进而得到的重力异常拟合公式为：

在公式(16)中，δt为十进制月重力场模型的时间与背景重力场模型的时间的差值，a为常值，近似于直线拟合的常值，b为重力异常长期变化率，i＝1代表年周期变化项，i＝2代表半年周期变化项，则年变化的振幅为i＝3表示与s2半日波相关的161天周期项，ωi为具体拟合时的角变化率，正弦函数和余弦函数拟合一条曲线时的变化率，ξ表示噪声和误差影响，属于不确定项。

步骤105具体包括：结合由grace卫星数据解算获得月重力数据，选取局部所需区域的月重力数据，并按公式(12)计算后，可以得到局部区域一系列的某个时间段内的重力值的变化量，即(δg,δt)值，这些值反映了某月重力异常变化与背景月重力异常及时间的关系。

结合公式(12)和公式(16)，可以得到重力异常长期变化率b。利用该重力异常长期变化率b实现对重力异常变化的预测。

步骤106具体包括：

基于重力异常长期变化率b的t时的重力异常时变公式为：

公式(17)中，t0＝2003.0为基准时间，为基准时间的重力异常值。

步骤108具体包括：

利用求得的t时的重力异常值对制备重力基准图时所采用的原始数据进行时变修正，得到时变修正数据，得到时变修正重力基准图。

具体为：输入时间t，计算十进制时间t'，利用公式(17)计算得到对应的重力异常值根据重力异常值利用插值算法进行网格化制图，得到时变修正重力基准图。

为实现上述目的，本发明还提供了一种基于重力异常时变的重力基准图时变修正系统。

图6为本发明实施例重力基准图时变修正系统的结构示意图，如图6所示，本发明实施例提供的所述重力基准图时变修正系统包括：

月重力场模型得到模块1，用于根据grace卫星数据解算获得重力数据，并根据所述重力数据，得到月重力场模型。

重力异常变化公式得到模块2，用于选取任意一个月重力场模型为背景重力场模型，并将所述月重力场模型的表达式和所述背景重力场模型的表达式进行差值计算，得到所述背景重力场下的重力异常变化公式。

平滑后的重力异常变化公式确定模块3，用于对所述重力异常变化公式进行空间平滑处理，确定平滑后的重力异常变化公式。

重力异常拟合公式得到模块4，用于对重力场时变影响因素进行分类拟合，得到重力异常拟合公式。

重力异常长期变化率确定模块5，用于根据所述平滑后的重力异常变化公式和所述重力异常拟合公式，确定重力异常长期变化率。

重力异常时变公式确定模块6，用于根据所述重力异常长期变化率，确定当前时刻的重力异常时变公式。

当前时刻重力异常值计算模块7，用于根据所述重力异常时变公式，计算当前时刻的重力异常值。

时变修正重力基准得到模块8，用于根据所述重力异常值，对制备重力基准图时所采用的原始数据进行时变修正，得到时变修正重力基准图。

本发明结合卫星时变重力场的研究成果，对重力场的时变现象进行定量研究和近似拟合，推到得到重力异常时变公式，利用重力异常时变公式对所用数据进行拟合后得到某一时刻的最新重力数据，从而实现对重力基准图的时变修正，进而制得导航精度更高的重力基准图。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。