一种适用于高精度惯导系统的重力场建模方法与流程

本发时涉及的是一种长航时载体中惯性导航装置适用的重力模型，特别涉及的是一种高精度惯性导航系统中三维重力矢量快速、精确的计算方法。

背景技术：

惯性导航系统作为一种自主式导航设备，已被广泛应用于航天、航空、航海以及地面载体等领域中。随着惯性器件(陀螺、加速度计)精度的不断提高，高精度纯惯性长航时的导航能力成为了新一代惯导系统的发展方向。特别地，随着冷原子技术的发展，惯性器件引起的导航位置误差每小时可降至米级，这就使得系统级误差成为了制约导航精度提升的主要因素。重力加速度是加速度计有害误差的重要补偿项，其模型计算精度对导航定位结果的影响也日渐突出。

目前，惯性导航系统广泛采用的是参考椭球表面的正常重力模型，以WGS84椭球为例，模型表达式如下：

其中，为大地纬度。重力加速度方向沿参考椭球面的法线方向，水平分量为零。而实际重力矢量垂直于大地水准面，二者之间的垂直分量偏差一般可达60～70mGal，水平垂线偏差在全球范围内最大可达100"，这一误差水平已经大大超过了现阶段加速度计零偏大小(10μg，相当于2"的垂线偏差)。

另一种对于地球重力场的精确描述方法，是广泛应用于测绘领域的重力场球谐级数表达式，表示的重力场引力位公式为：

式中，f为引力常数；M表示地球质量；a表示地球赤道半径；ρ表示点的向径；n为调和项阶数；m为调和项次数；λ表示经度；θ表示余纬，即P_nm表示缔合勒让德函数；C_nm与S_nm为调和系数，采用精密的地球重力卫星测量数据以及地面测量数据计算得到。引力位V_E加上离心力位Φ即可表示重力位U。现阶段已有大量以球谐级数式表示的高阶重力场模型问世，如EGM系列、ICGEM系列模型等，模型最高阶次数达到2160阶、分辨率约为9.25km，这些模型可以足够精确的描述地球重力场。然而，将球谐级数重力场模型应用在高频率实时解算的惯导系统中，还需要解决如下问题：

1.模型复杂度问题

重力场球谐级数式精度随阶数增加而提高，对于阶数为n的模型，其系数总数为(n-1)(n+3)个，即，如果采用2160阶的模型，需要计算机提供额外大于200MB的存储空间，并完成4669917个数据的实时调用。同时，在每个阶数下都会涉及到完全正则化勒让德函数的递归计算。这使得保证模型精度条件下，复杂度也相对提高，进而不利于惯导系统的实时解算。因此，直接将高阶重力模型替代传统模型无法适用于惯导系统。

2.计算时效性问题

常规惯导系统导航结果输出频率一般在100-200Hz，重力加速度矢量计算作为导航解算过程中的一个运算步骤要求其具有较高的实效性。而球谐模型的计算速度与模型阶数呈二次方递增的关系，即模型阶数每提高10倍，计算量就会增大100倍，无法满足导航计算机更新频率的要求。

技术实现要素：

本发明的目的是为了解决上述问题，提出一种适用于高精度惯导系统的重力场建模方法，采用投影的技术手段，通过将球谐模型计算得到的格网重力加速度信息映射至立方体平面内的方式，使复杂球面模型转化为二维插值计算，实现精确重力场模型的快速计算，具有适用于高精度惯导系统实时解算的优势。

本发明一种适用于高精度惯导系统的重力场建模方法，包括以下几个步骤：

步骤1：将地球表面进行平面延展，得到经度、纬度范围均为[0,2π)的二维平面S1；再将地球整体旋转90度，使极区转至赤道位置，经延展得到经度、纬度范围均为[0,2π)的二维平面S2；这里对纬度范围进行了一周期延拓，对应以实际纬度范围[-π/2,π/2]计算的重力场信息扩展结果：

式中，为以大地坐标系表示的地球引力位；表示大地坐标系下的纬度和经度；ρ表示点的向径。平面坐标系(p,q)与球面经纬度坐标映射关系：

步骤2：根据导航精度要求确定重力格网分辨率N，利用球谐模型位场与场元关系式计算平面S1与S2上相应格网节点位置上的重力矢量值，每个面上共3N ²个数据；

步骤3：在平面S1上取四个正方形截面1,2,3,4；在平面S2上取两个正方形截面5,6。分别对应地球表面的固定区域，其映射关系如下表1所示：

表1映射关系

步骤4：采用二维快速傅里叶变换将区间S1、S2内重力矢量各分量值换算为新模型节点系数s_j,k(j,k＝1,2,…,N)；

步骤5：利用B样条插值对格网区间重力场进行平滑计算，样条阶数记为κ，建立重力模型计算式：

其中，β_κ(*)表示B样条插值基函数。

步骤6：根据预定轨迹范围在导航计算机内事先存储一定数量的节点系数，获取载体运行位置信息，利用式(2)表示的模型实时计算当前位置下的重力矢量值。

所述的步骤2中，格网节点位置上的重力矢量值，其计算方式可以有以下两种：

方式一：

利用地球引力位球谐级数展开式：

式中，f为引力常数；M表示地球质量；a表示地球赤道半径；ρ表示点的向径；n为调和项阶数；m为调和项次数；θ表示余纬，即P_nm表示缔合勒让德函数；C_nm与S_nm为调和系数。取无穷项截短到n_max次，再加上离心力位Φ(Φ＝Ω²(x²+y²)/2，其中，Ω为地球自转角速率，坐标(x,y,z)定义在地心坐标系下)即可得到重力近似势能W。利用如下关系式转换至地理坐标系：

式中，R_N表示卯酉圈曲率半径；a表示椭球长半轴；b表示椭球短半轴；h表示相对参考椭球面的高度。

重力矢量即在地理标系中的分量为：

式中，g_E、g_N、g_U分别表示重力的东向、北向、天向分量。

方式二：

根据扰动重力定义，得到扰动重力位的球谐表达式：

式中，C_nm表示经椭球参数修正后的调和系数。利用地面重力异常、垂线偏差与扰动位间的关系：

其中，Δg为地面重力异常值；ξ、η分别为地面垂线偏差子午面分量和卯酉面分量；γ为正常重力。得到地理坐标系下的重力矢量表示为：

通过上述方式给定节点位置坐标即可得到相应的重力矢量参数。

所述的步骤4中，节点系数计算方式如下：

首先给定B样条曲线表示的平面S1、S2上重力矢量计算式为：

其中，s_j,k为待求的节点系数；β_κ(*)为B样条插值基函数，且当样条阶数κ确定时表达式已知。根据节点处重力矢量值，在平面S1、S2上的重力矢量还可以表示成如下拉格朗日多项式的形式：

式中，系数A_j,k即表示节点处的重力矢量值，为已知量；L(*)表示狄拉克函数，当n＝0时，L(n)＝1；当n≠0时，L(n)＝0；n为整数，表示节点。由于在节点处L(n)可以表示为插值基函数值β(n)的线性组合，在一维情况下二者的傅里叶变换有如下关系：

式中，

显然，a(ω)的计算可以通过β_2κ+1(n)的离散傅里叶变换得到。因此，在频域空间利用二维N点FFT作用于系数A_j,k，除以样条关系因子a(ω)(在离散傅里叶变换后得到的仍为频域下的N个离散点，其中一点可以表示为a(j/N))，再经过傅立叶逆变换即可得到s_j,k。

本发明一种适用于高精度惯导系统的重力场建模方法的优点在于：

(1)利用球面至正方体表面一一对应关系，将球面复杂连续函数转化为平面格网数据，进而使得重力模型的复杂度得到化简；

(2)利用二维样条插值算法，建立简化后的重力模型，在保证模型计算精度的同时，使得模型适用于惯导系统实时计算。

附图说明

图1为本发明方法的实施步骤流程图；

图2为本发明平面S1和S2上地面重力场元计算值及立方体表面映射关系图；

图3为本发明格网点上B样条模型系数计算结果图；

图4为本发明B样条模型与求谐函数计算用时对比曲线图；

具体实施方式

本发明在于提供一种适用于高精度惯导系统的重力场建模方法，采用投影的技术手段，通过将球谐模型计算得到的格网重力加速度信息映射至立方体平面内的方式，使复杂球面模型转化为二维插值计算，实现精确重力场模型的快速计算，具有适用于高精度惯导系统实时解算的优势。

如图1所示，本发明包括以下几个步骤：

步骤1：参照载体预定的运行轨迹(经度范围纬度范围[λ₀,λ₁])，将地球表面进行平面延展，得到经度、纬度范围均为[0,2π)的二维平面S1；再将地球整体旋转90度，使极区转至赤道位置，经延展得到经度、纬度范围均为[0,2π)的二维平面S2。这里对纬度范围进行了一周期延拓，对应以实际纬度范围[-π/2,π/2]计算的重力场信息扩展结果：

式中，为以大地坐标系表示的地球引力位；表示大地坐标系下的纬度和经度；ρ表示点的向径。所得两个平面如附图2所示(图中编号1-6所示区域对应实际地球表面)。平面坐标系(p,q)与球面经纬度坐标映射关系为：

步骤2：根据导航精度要求确定重力格网分辨率N＝8192(等同于球面分辨率约为2.6')，利用如下方式计算平面S1与S2上相应格网节点位置上的重力矢量值：

根据扰动重力定义，得到扰动重力位的球谐表达式：

式中，f为引力常数；M表示地球质量；a表示地球赤道半径；ρ表示点的向径；n为调和项阶数；m为调和项次数；θ表示余纬，即P_nm表示缔合勒让德函数；C_nm与S_nm为已知调和系数，C_nm表示经椭球参数修正后的调和系数，以WGS84椭球为参考椭球，C_nm的低阶项偶阶带谐系数更新为：

式中，J₂，J₄，J₆，J₈，J₁₀表示参考椭球的偶阶带谐系数，为已知量。取无穷项截短到n_max次，利用地面重力异常、垂线偏差与扰动位间的关系：

其中，Δg为地面重力异常值；ξ、η分别为地面垂线偏差子午面分量和卯酉面分量；γ为正常重力。得到地理坐标系下的重力矢量表示为：

通过上述方式给定节点位置坐标即可得到相应的重力矢量参数。得到的两个平面上重力异常如附图2(a)、(b)所示。

步骤3：在平面S1上取四个正方形截面1,2,3,4；在平面S2上取两个正方形截面5,6。分别对应地球表面的固定区域，如附图2(a)、(b)内编号1-6的正方形截面，其映射关系如下表所示：

步骤4：采用二维快速傅里叶变换将区间S1、S2内重力矢量各分量值换算为新模型节点系数s_j,k(j,k＝1,2,…,N)。节点系数计算过程如下：

(1)利用如下递推关系得到κ阶(令κ＝11)和2κ+1阶B样条节点n上的函数值：

其中，

(2)对β_2κ+1(n)按式(4)做N点FFT变换，得a(ω)(在离散傅里叶变换后得到的仍为频域下的N个离散点，其中一点可以表示为a(j/N))：

(3)建立以拉格朗日多项式表示的平面S1、S2上的重力矢量计算式：

其中，A_j,k表示节点处的重力矢量值；L(*)表示狄拉克函数。

(4)对式(5)表示的节点系数A_j,k做二维N点FFT，得到

(5)利用上述结果得到B样条重力矢量计算式节点系数的傅里叶变换值。

(6)对式(6)中进行傅立叶逆变换得到s_j,k，数据结果如附图3所示。

步骤5：用B样条插值对格网区间重力场进行平滑计算，样条阶数记为κ，建立重力矢量计算式：

其中，β_κ(*)表示B样条插值基函数。

步骤6：根据预定轨迹范围在导航计算机内事先存储一定数量的节点系数，获取载体运行位置信息，利用步骤5中式(7)表示的模型实时计算当前位置下的重力矢量值。

本发明通过以上步骤得到的有益效果如附图4所示。利用本发明所述的方法，重力模型计算用时不随模型阶数而改变，相对计算时间较球谐模型有了较大提高。