一种适用于高精度惯导系统的重力场建模方法技术方案

本发明专利技术公开了一种适用于高精度惯导系统的重力场建模方法，解决了传统惯导系统中正常重力模型由于椭球近似造成的精度不足问题，以及高阶重力场球谐模型在惯导系统中的不适用问题。本方法利用球心投影的一一映射关系，通过将球谐模型计算得到的高精度格网重力加速度信息映射至立方体平面的方式，使复杂球面模型递推计算转化为二维平面插值计算，实现精确重力场模型的快速计算，具有适用于高精度惯导系统实时解算的优势。

【技术实现步骤摘要】
一种适用于高精度惯导系统的重力场建模方法
本发时涉及的是一种长航时载体中惯性导航装置适用的重力模型，特别涉及的是一种高精度惯性导航系统中三维重力矢量快速、精确的计算方法。
技术介绍
惯性导航系统作为一种自主式导航设备，已被广泛应用于航天、航空、航海以及地面载体等领域中。随着惯性器件(陀螺、加速度计)精度的不断提高，高精度纯惯性长航时的导航能力成为了新一代惯导系统的发展方向。特别地，随着冷原子技术的发展，惯性器件引起的导航位置误差每小时可降至米级，这就使得系统级误差成为了制约导航精度提升的主要因素。重力加速度是加速度计有害误差的重要补偿项，其模型计算精度对导航定位结果的影响也日渐突出。目前，惯性导航系统广泛采用的是参考椭球表面的正常重力模型，以WGS84椭球为例，模型表达式如下：其中，为大地纬度。重力加速度方向沿参考椭球面的法线方向，水平分量为零。而实际重力矢量垂直于大地水准面，二者之间的垂直分量偏差一般可达60～70mGal，水平垂线偏差在全球范围内最大可达100"，这一误差水平已经大大超过了现阶段加速度计零偏大小(10μg，相当于2"的垂线偏差)。另一种对于地球重力场的精确描述方法，是广泛应用于测绘领域的重力场球谐级数表达式，表示的重力场引力位公式为：式中，f为引力常数；M表示地球质量；a表示地球赤道半径；ρ表示点的向径；n为调和项阶数；m为调和项次数；λ表示经度；θ表示余纬，即Pnm表示缔合勒让德函数；Cnm与Snm为调和系数，采用精密的地球重力卫星测量数据以及地面测量数据计算得到。引力位VE加上离心力位Φ即可表示重力位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上相应格网节点位置上的重力矢量值，每个面上共3N2个数据；步骤3：在平面S1上取四个正方形截面1,2,3,4；在平面S2上取两个正方形截面5,6。分别对应地球表面的固定区域，其映射关系如下表1所示：表1映射关系步骤4：采用二维快速傅里叶变换将区间S1、S2内重力矢量各分量值换算为新模型节点系数sj,k(j,k＝1,2,…,N)；步骤5：利用B样条插值对格网区间重力场进行平滑计算，样条阶数记为κ，建立重力模型计算式：其中，βκ(*)表示B样条插值基函数。步骤6：根据预定轨迹范围在导航计算机内事先存储一定数量的节点系数，获取载体运行位置信息，利用式(2)表示的模型实时计算当前位置下的重力矢量值。所述的步骤2中，格网节点位置上的重力矢量值，其计算方式可以有以下两种：方式一：利用地球引力位球谐级数展开式：式中，f为引力常数；M表示地球质量；a表示地球赤道半径；ρ表示点的向径；n为调和项阶数；m为调和项次数；θ表示余纬，即Pnm表示缔合勒让德函数；Cnm与Snm为调和系数。取无穷项截短到nmax次，再加上离心力位Φ(Φ＝Ω2(x2+y2)/2，其中，Ω为地球自转角速率，坐标(x,y,z)定义在地心坐标系下)即可得到重力近似势能W。利用如下关系式转换至地理坐标系：式中，RN表示卯酉圈曲率半径；a表示椭球长半轴；b表示椭球短半轴；h表示相对参考椭球面的高度。重力矢量即在地理标系中的分量为：式中，gE、gN、gU分别表示重力的东向、北向、天向分量。方式二：根据扰动重力定义，得到扰动重力位的球谐表达式：式中，Cnm表示经椭球参数修正后的调和系数。利用地面重力异常、垂线偏差与扰动位间的关系：其中，Δg为地面重力异常值；ξ、η分别为地面垂线偏差子午面分量和卯酉面分量；γ为正常重力。得到地理坐标系下的重力矢量表示为：通过上述方式给定节点位置坐标即可得到相应的重力矢量参数。所述的步骤4中，节点系数计算方式如下：首先给定B样条曲线表示的平面S1、S2上重力矢量计算式为：其中，sj,k为待求的节点系数；βκ(*)为B样条插值基函数，且当样条阶数κ确定时表达式已知。根据节点处重力矢量值，在平面S1、S2上的重力矢量还可以表示成如下拉格朗日多项式的形式：式中，系数Aj,k即表示节点处的重力矢量值，为已知量；L(*)表示狄拉克函数，当n＝0时，L(n)＝1；当n≠0时，L(n)＝0；n为整数，表示节点。由于在节点处L(n)可以表示为插值基函数值β(n)的线性组合，在一维情况下二者的傅里叶变换有如下关系：式中，显然，a(ω)的计算可以通过β2κ+1(n)的离散傅里叶变换得到。因此，在频域空间利用二维N点FFT作用于系数Aj,k，除以样条关系因子a(ω)(在离散傅里叶变换后得到的仍为频域下的N个离散点，其中一点可以表示为a(j/N))，再经过傅立叶逆变换即可得到sj,k。本专利技术一种适用于高精度惯导系统的重力场建模方法的优点在于：(1)利用球面至正方体表面一一对应关系，将球面复杂连续函数转化为平面格网数据，进而使得重力模型的复杂度得到化简；(2)利用二维样条插值算法，建立简化后的重力模型，在保证模型计算精度的同时，使得模型适用于惯导系统实时计算。附图说明图1为本专利技术方法的实施步骤流程图；图2为本专利技术平面S1和S2上地面重力场元计算值及立方体表面映射关系图；图3为本专利技术格网点上B样条模型系数计算结果图；图4为本专利技术B样条模型与求谐函数计算用时本文档来自技高网...

【技术保护点】
一种适用于高精度惯导系统的重力场建模方法，包括以下几个步骤：步骤1：将地球表面进行平面延展，得到经度、纬度范围均为[0,2π)的二维平面S1，再将地球整体旋转90度，使极区转至赤道位置，经延展得到经度、纬度范围均为[0,2π)的二维平面S2，对纬度范围进行一周期延拓，对应以实际纬度范围[‑π/2,π/2]计算的重力场信息扩展结果：式中，为以大地坐标系表示的地球引力位，表示大地坐标系下的纬度和经度，ρ表示点的向径，平面坐标系(p,q)与球面经纬度坐标映射关系：步骤2：根据导航精度要求确定重力格网分辨率N，利用球谐模型位场与场元关系式计算平面S1与S2上相应格网节点位置上的重力矢量值，每个面上共3N2个数据；步骤3：在平面S1上取四个正方形截面，分别为截面1、截面2、截面3、截面4；在平面S2上取两个正方形截面，分别为截面5、截面6；分别对应地球表面的固定区域，其映射关系表1所示：表1 映射关系步骤4：采用二维快速傅里叶变换将平面S1、S2内重力矢量各分量值换算为新模型节点系数sj,k，j,k＝1,2,…,N；步骤5：利用B样条插值对格网区间重力场进行平滑计算，样条阶数记为κ，建立重力模型计算式：S(p,q)=Σj=0N-1Σk=0N-1sj,kβκ(NP-j)βκ(Nq-k)---(2)]]>其中，βκ(*)表示B样条插值基函数；步骤6：根据预定轨迹范围在导航计算机内预先存储节点系数，获取载体运行位置信息，利用式(2)的模型实时计算当前位置下的重力矢量值。...

【技术特征摘要】
1.一种适用于高精度惯导系统的重力场建模方法，包括以下几个步骤：步骤1：将地球表面进行平面延展，得到经度、纬度范围均为[0,2π)的二维平面S1，再将地球整体旋转90度，使极区转至赤道位置，经延展得到经度、纬度范围均为[0,2π)的二维平面S2，对纬度范围进行一周期延拓，对应以实际纬度范围[-π/2,π/2]计算的重力场信息扩展结果：式中，为以大地坐标系表示的地球引力位，表示大地坐标系下的纬度和经度，ρ表示点的向径，平面坐标系(p,q)与球面经纬度坐标映射关系：步骤2：根据导航精度要求确定重力格网分辨率N，利用球谐模型位场与场元关系式计算平面S1与S2上相应格网节点位置上的重力矢量值，每个面上共3N2个数据；步骤3：在平面S1上取四个正方形截面，分别为截面1、截面2、截面3、截面4；在平面S2上取两个正方形截面，分别为截面5、截面6；分别对应地球表面的固定区域，其映射关系表1所示：表1映射关系步骤4：采用二维快速傅里叶变换将平面S1、S2内重力矢量各分量值换算为新模型节点系数sj,k，j,k＝1,2,…,N；步骤5：利用B样条插值对格网区间重力场进行平滑计算，样条阶数记为κ，建立重力模型计算式：其中，βκ(*)表示B样条插值基函数；步骤6：根据预定轨迹范围在导航计算机内预先存储节点系数，获取载体运行位置信息，利用式(2)的模型实时计算当前位置下的重力矢量值。2.根据权利要求1所述的一种适用于高精度惯导系统的重力场建模方法，所述的步骤2中，格网节点位置上的重力矢量值具体为：利用地球引力位球谐级数展开式：式中，f为引力常数，M表示地球质量，a表示地球赤道半径，ρ表示点的向径，n为调和项阶数，m为调和项次数，θ表示余纬，即Pnm表示缔合勒让德函数，Cnm与Snm为调和系数；取无穷项截短到nmax次，再加上离心力位Φ，Φ＝Ω2(x2+y2)/2，其中，Ω为地球自转角速率，坐标(x,y,z)定义在地心坐标系下，最终得到重力近似势能W，利用如下关系式转换至地理坐标系：式中，RN表示卯酉圈曲率半径，a表示椭球长半轴，b表示椭球短半轴，h表示相对参考椭球面的高度；重力矢量即在地理标系中的分量为：式中，gE、gN、gU分别表示重力的东向、北向、天向分量。3.根据权利要求1所述的一种适用于高精度惯导系统的重力场建模方法，所述的步骤2中，格网节点位置上的重力矢量值具体为：根据扰动重力定义，得到扰动重力位的球谐表达式：

【专利技术属性】
技术研发人员：杨功流，王晶，李晶，刘元元，周潇，
申请(专利权)人：北京航空航天大学，
类型：发明
国别省市：北京;11