基于球谐模型与多传感器融合的高精度重力扰动补偿方法

引用格式：刘宇鑫，王新龙，王勋，等.基于球谐模型与多传感器融合的高精度重力扰动补偿方法［J］.航空兵器，2023，30（1）：104-113.

LiuYuxin，WangXinlong，WangXun，etal.AnAccurateGravityDisturbanceCompensationMethodBasedonSphericalHarmonicModelandMultiSensorFusion［J］.AeroWeaponry，2023，30（1）：104-113.（inChinese）

摘要：随着对高精度长航时惯导系统性能要求的提升，重力扰动已成为惯导系统的主要误差源，能否对其进行有效补偿是提升导航精度的关键因素。传统基于球谐模型的补偿方法无法反映地球重力场的细节信息，对中短波重力扰动分量补偿效果不佳；传统状态估计法中马尔科夫状态模型对长波重力扰动分量描述精度较差，因此对长波分量的补偿效果不佳。以上方法均无法对波段较宽的实际重力扰动进行精确补偿，针对这一问题，本文设计了一种综合利用球谐模型与多传感器信息融合的高精度重力扰动补偿方法。该方法一方面利用低阶球谐模型对长波重力扰动分量计算精度较高的特点，对长波分量进行补偿；另一方面利用捷联惯导/激光多普勒测速/气压高度计构成完全自主的组合导航系统，将残余的中短波重力扰动分量建立为高精度的马尔科夫模型，从而实现对中短波分量的状态估计与补偿。仿真结果表明，所提出的重力扰动补偿方法可有效改善重力扰动估计效果，进而实现高精度的重力扰动补偿，具有较高的导航精度。

关键词：捷联惯导系统；重力扰动补偿；重力场球谐模型；状态估计；组合导航

中图分类号：TJ765;V249.32+2

文献标识码：A

文章编号：1673-5048（2023）01-0104-10

DOI：10.12132/ISSN.1673-5048.2022.0038

0引言

惯性导航系统具有高度的自主性、隐蔽性以及信息完备等特点，目前广泛应用与军事、工程、科学研究以及民用领域中［1］。在惯导解算过程中，通常采用理论重力模型计算理论重力并从比力测量值中对其进行补偿［2］。但理论重力与实际重力并不相等，该差异即为重力扰动，其导致的补偿残差会造成惯导解算误差。多数地区重力扰动的大小在15～80mGal之间［3］，是高精度惯导系统的主要误差源之一［4-5］，对导航精度有较大的影响。

目前重力补偿方法可根据是否已知载体运动区域的重力网格数据而分为两大类。在已知区域重力网格数据库时，通过对该数据库进行插值获取当前位置的重力扰动矢量并对其进行补偿［6-8］。这种方法在局部区域内的补偿精度较高，且具有较好的实时性，但无法在该数据库未覆盖的区域适用。在工程中高精度的重力扰动测量数据往往较难获取，因此这种补偿方法的应用区域受限。在缺乏区域重力网格数据库时，重力扰动补偿方法又可细分为两种：基于重力场模型的补偿方法以及基于状态估计的补偿方法。

在基于重力场模型的补偿方法中，通常利用EGM2008球谐模型计算重力扰动并进行补偿。这种方法适用于全球任意位置的重力扰动补偿，但其计算结果无法反映地球重力场的细节信息［9-10］，因此导致该方法在山地等重力扰动变化较为显著区域的补偿精度较低，以中国南部高原区为例，其计算误差的标准差达到20.74mGal［11］，难以满足高精度的补偿需求。此外，高阶球谐模型的计算负担较大，模型参数占用存储空间大，计算耗时长［12］，不适用于传统导航系统的硬件配置环境，并且难以满足惯导系统解算实时性的要求。

状态估计法将捷联惯导系统与其他导航传感器结合成组合导航系统，通过引入重力扰动状态量，建立其状态空间模型，利用最优估计算法实现对重力扰动的估计与补偿［13-14］。这类方法在补偿重力扰动的同时，能有效抑制惯导系统误差的发散，因此具有较好的补偿效果。但在状态估计过程中，重力扰动与惯导系统误差耦合，需要借助高精度的重力扰动状态模型才能对其进行有效估计，即重力扰动的建模精度是制约该方法重力补偿精度的重要因素［15］。实际重力扰动可视为一系列不同幅值和波长的重力扰动分量信号的叠加，其所覆盖的波段很广，包含波长约1～10000km的分量信号［16］。然而常用的马尔科夫重力扰动模型所能覆盖的波段有限，当相关距离取几十千米时，能对100km以下的重力扰动分量进行较为精确的描述［15，17］，但对长波分量的描述效果不佳，因此状态估计法通常只能对中短波重力扰动分量进行较为有效的估计，对全波段的重力扰动估计精度有限。

针对以上问题，本文设计了一种球谐模型补偿与状态估计补偿相互结合的重力扰动补偿方法，在缺乏区域重力网格数据库时能完全自主地实现重力扰动的估计与补偿，并对该方法的性能进行了验证。

1重力扰动对惯导解算影响分析

1.1重力扰动成因

在地球物理学中通常人为地选择一个形状规则、密度均匀的旋转椭球体（即参考椭球）对地球做近似，由该参考椭球产生的重力为理论重力，记为γ0。以常用的1980年国际大地测量参考椭球为例，理论重力加速度计算公式为［18］

γ0=9.780325×［1+a1sin2L+a2sin2（2L）］－a3h（1）

式中：a1=0.00530240；a2=-0.00000582；a3=0.00003086；L为地理纬度；h为载体高度。

不过实际上地球并非是理想的旋转椭球体，其形状与参考椭球并不完全相同，而且地球内部物质分布结构不规则，导致局部地区的密度不均，从而使得某一点处的实际重力加速度g与理论重力加速度γ0存在一定差异，由此可将这一差异定义为重力扰动δg，即

δg=g－γ0（2）

为了便于研究，通常将某点处实际重力方向与理论重力方向之间的夹角称为垂线偏差，将其在卯酉圈上的分量记为η，在子午圈方向上的分量记为ξ；将重力扰动矢量在理论重力矢量方向上的分量称为重力异常，记为δgU，如图1所示。

由几何关系可得重力扰动矢量在地理坐标系中的投影：

δg=［－γ0η－γ0ξ－δgU］T（3）

由此可见，重力扰动会引起实际重力与理论重力在当地地理坐标系下三个方向的误差。

1.2重力扰动对惯导影响机理

捷联惯导系统利用固联于载体上的加速度计和陀螺仪分别测量载体相对惯性空间的比力f～b和角速度ω～bib，并以此为基础进行导航解算，其原理如图2所示。

在解算姿态时，利用陀螺测得角速度结合载体的位置和速度计算得到姿态角速度，进而通过姿态微分方程实时更新载体本体系b系与导航坐标系n系之间的坐标转换矩阵C^nb，s（为与后续其他系统进行区分，利用下标s表示SINS相关计算参数），从而解算得到姿态角。在解算位置时，对加速度测量所得比力中的重力等有害加速度进行补偿得到载体运动加速度，经过两次积分运算后可分别解算得到当前的速度和位置。

考虑惯性器件测量误差以及初始误差的影响，可得惯导误差方程：

δv·ns=f^n×φns－（2ω^nie，s+ω^nen，s）×δvns－（2δωnie+δωnen）×v^ns－δgn0+C^nb，sδfb（4）

δL·s=δvN，sRM+h^s－v^N，s（RM+h^s）2δhs

δλ·s=secL^sRN+h^sδvE，s+v^E，ssecL^stanL^sRN+h^sδLs－v^E，ssecL^s（RN+h^s）2δhsδh·s=δvU，s（5）

φ·ns=δωnie，s－ω^nin，s×φns+δωnen，s－Cnb，sδωbib，s（6）

式中：RM为地球的子午圈半径；RN为地球的卯酉圈半径；ω^nie，s为地球自转角速度；δωnie，s为地球自转角速度的计算误差；ω^nen，s为导航系相对地球系的转动角速度；δωnie，s为导航系相对地球系的转动角速度的计算误差；δvns=［δvE，sδvN，sδvU，s］为速度误差；δpns=［δLsδλsδhs］T为位置误差；φs为姿态失准角；δg0为理论重力加速度的计算误差；δfb为比力测量误差；δωbib为陀螺仪的测量误差。

然而由于重力扰动的影响，惯导工作时用于补偿比力测量值中有害加速度的理论重力加速度与实际重力加速度并不相等，实际重力加速度为理论重力加速度与重力扰动之和，即

gn=gn0+δgn0（7）

仅采用理论重力加速度公式对比力测量值进行补偿，会导致计算结果中存在一定的重力扰动残差，进而引起额外的加速度解算误差，此时加速度解算误差δv～·ns为

δv～·ns=δv·ns－δgn（8）

将式（4）代入式（8），则可得考虑重力扰动时加速度解算误差的具体形式：

δv～·ns=f^n×φns－（2ω^nie，s+ω^nen，s）×δvns－δgn－

C^nb，s

Δb－（2δωnie，s+δωnen，s）×vns－δgn0（9）

式中：Δb为加速度计测量误差。

受重力扰动影响的加速度误差经过两次积分后分别引起速度误差δv^ns与位置误差δp^ns。速度和位置误差又进一步引起角速度计算误差δω^nie和δω^nen：

δω^nie=0－ωiesinL^s·δL^s

ωiecosL^s·δL^s（10）

δω^nen=－δv^N，sRM+h^s+v^N，sδh^s（RM+h^s）2δv^E，sRN+h^s－v^E，sδh^s（RN+h^s）2

δv^E，stanL^sRN+h^s+v^E，sδL^ssec2L^sRN+h^s－v^E，sδh^stanL^s（RN+h^s）2（11）

由式（8）～（11）可知，重力扰动所引起额外的加速度解算误差会导致额外的角速度计算误差，进而通过式（6）引起额外的姿态解算误差。姿态误差又进一步通过式（4）～（5）造成位置误差和速度误差，形成闭环回路，持续对惯导系统的位置、速度和姿态解算产生影响。重力扰动对捷联惯导系统的影响机理与加速度计零偏类似，其所造成的惯导解算误差随时间积累。由以上分析可知，重力扰动是影响高精度长航时惯导系统性能的关键因素之一。

2球谐模型与状态估计相结合的补偿方案设计

众所周知，捷联惯导系统通常采用理论重力加速度模型对比力测量值中的实际重力加速度进行粗略补偿，因此其导航解算结果中包含由理论与实际重力差异（即重力扰动）所引起的导航误差。

为了实现重力扰动的估计与补偿，可借助其他导航传感器提供不受重力扰动影响的导航参数，利用状态估计算法从导航解算误差中分离得到重力扰动估计值，并对其进行补偿。与此同时，针对状态估计法中马尔科夫状态模型对长波重力扰动描述精度有限的问题，考虑到导航系统的硬件配置限制，截取低阶EGM2008球谐模型对长波的重力扰动进行补偿，从而提升整体的补偿效果。

2.1球谐模型补偿

球谐模型补偿部分根据EGM2008模型计算当前解算所得位置处的长波重力扰动分量，并在补偿有害加速度时直接对加速度解算值进行校正，其原理如图3所示。

通常用空间分辨率λres表示球谐模型对重力扰动描述的精细程度，其计算公式如下：

λres=πaN（12）

式中：a为地球的半长轴长度；N为球谐函数的最高阶数。球谐模型的最高阶数越大，相应计算结果的分辨率越高，但受导航平台硬件存储空间以及计算实时性的限制，EGM2008模型阶数小于12阶时才能满足导航系统计算资源要求［5］，因此截取前12阶的参数计算重力扰动的长波分量，其具体计算公式可参见文献［19］。