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

刘宇鑫 王新龙 王勋 高文宁 胡晓东

引用格式：刘宇鑫，王新龙，王勋，等.基于球谐模型与多传感器融合的高精度重力扰动补偿方法［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］。

2.2状态估计补偿

通过将SINS与其他不受重力扰动影响的导航传感器结合形成组合导航系统，对惯导系统误差进行估计，同时通过最优估计方法实现中短波重力扰动分量的估计与补偿。

2.2.1陀螺/激光多普勒测速系统

在其他导航传感器的选择上，激光多普勒测速仪（LDV）具有测量精度高、动态响应快、测量范围大、测量线性度好等优点，能提供高精度的速率测量信息［20］，通过与陀螺仪组合形成的导航系统可提供完备的导航信息，其原理如图4所示。

后文中用下标D表示陀螺/多普勒测速仪系统的相关计算参数，以与SINS相关计算参数进行区分。利用激光多普勒测速仪测量载体相对地面的运动速率，利用陀螺仪测量所得角速度实时计算载体的姿态矩阵C^nb，D，再根据已知的安装矩阵Cbm，将多普勒测速仪在其安装坐标系下所测的速度v～mD投影到地理坐标系下得到

v^nD=v^E，Dv^N，Dv^U，DT=C^nb，DCbmv～mD（13）

由此可根据位置微分方程解算得到位置信息：

p^nD=L^Dλ^Dh^DT（14）

式中：L^·D=v^N，DRM+h^D；λ^·D=secL^DRN+h^Dv^E，D；h^·D=v^U，D。

与捷联惯导系统类似，利用陀螺仪测量所得角速度ω～bib以及由位置和速度计算所得的ω^nie，D和ω^nen，D，可对姿态矩阵C^nb，D进行更新，即

C^·nb，D=C^nb，D（Ω～bib）－（Ω^nie，D+Ω^nen，D）C^nb，D（15）

式中：Ω～bib，Ω^nie，D和Ω^nen，D分别为ω～bib，ω^nie，D和ω^nen，D构成的反对称矩阵。

由式（13）～（15）可知，陀螺/激光多普勒测速导航系统直接利用激光多普勒测速仪测得的速度进行解算，不涉及重力的补偿，因此其导航结果不受重力扰动的影响，可在状态估计法中与惯导解算结果进行匹配，用于重力扰动的估计。

2.2.2状态估计补偿系统

利用SINS与陀螺/激光多普勒测速系统构成组合导航系统，与此同时，由于两导航系统均依赖陀螺的测量值，误差会随时间累积，因此用气压高度计对系统进行阻尼，其原理如图5所示。

将SINS和陀螺/LDV系统解算所得导航参数之差以及气压高度计提供的高度信息作为量测值，同时建立中短波重力扰动分量的状态模型，进而利用最优估计算法对重力扰动以及导航误差进行状态估计与补偿。

2.3重力扰动补偿方案设计

综合球谐模型以及状态估计补偿法，可得如图6所示的重力扰动补偿方案。

所设计方案同时采取球谐模型以及状态估计两种途径对重力扰动进行补偿。由解算所得位置，利用球谐模型计算长波重力扰动分量并进行补偿；基于组合导航系统对导航参数误差、传感器测量误差以及残余的中、短波重力扰动分量进行估计与补偿，进一步降低重力扰动对惯导解算的影响。

3重力扰动补偿系统模型的建立

3.1状态模型

在所设计重力扰动补偿系统中，滤波器的状态模型由4部分组成，包括惯性器件测量误差模型、重力扰动状态模型、捷联惯导系统误差模型以及陀螺/多普勒测速仪导航系统误差模型。

（1）惯性器件误差模型

对于高精度的惯性器件，可将其测量误差建模为常值误差与白噪声之和，即

ε=εb+wg（16）

Δ=Δb+wa（17）

式中：ε为陀螺仪元件的测量误差；εb和wg分别为陀螺仪测量的常值漂移和白噪声；Δ为加速度计的测量误差；Δb和wa分别为加速度计测量的常值偏置和白噪声。

（2）重力扰动状态模型

在利用EGM2008球谐模型对长波重力扰动分量进行补偿后，重力扰动残差在长波段的能量较小，因此可以利用带阻尼的二阶高斯-马尔可夫模型来增强模型在长波段的衰减，从而与实际的重力扰动补偿残差更加贴合。该模型可表示为［15］

δg¨i（t）=－2ξω0，iδg·i（t）－ω20，iδgi（t）+qi（t）（18）

式中：i分别表示E，N，U；ξ为阻尼系数，通过调节阻尼值可以得到具有不同功率谱特性的重力扰动模型；ω0，i为马尔科夫模型的中心频率，可表示为2πv/di，di为重力扰动相关距离，v为载体的运动速度；qi（t）为高斯白噪声，其方差Qi=4ξ（2πv/di）3σ2δg，i，σ2δg，i为各向重力扰动分量的方差。本文中取ξ=0.7，dE=245.4km，dN=213.3km，dU=155.5km，三向σ2δg，i由球谐模型近似计算为1069.29mGal2。

（3）捷联惯导系统误差模型

考虑重力扰动影响的惯导速度误差方程如式（9）所示，位置误差方程和姿态误差方程如式（5）～（6）所示。

（4）陀螺/激光多普勒测速系统误差模型

a.速度誤差方程

激光多普勒测速仪测量的是载体在其安装坐标系下相对地面的运动速度，其解算所得载体在地理系中的速度受到姿态误差角φD、三向安装方位误差角α以及标度因数误差δk的影响，速度误差方程为

δvnD=v^nD×φD+Cnmαv^mD+δkCnmv^mD（19）

式中：δvnD=［δvE，DδvN，DδvU，D］T为陀螺/多普勒测速仪导航系统解算的速度误差。

b.位置误差方程

对式（14）的位置更新方程等式两边同时求微分，可得位置误差方程：

δL·D=δvN，DRM+h^D－V^N，D（RM+h^D）2δhDδλ·D=secL^DRN+h^DδvE，D－v^E，DsecL^D（RN+h^D）2δhD+

v^E，DsecL^DtanL^DRN+h^DδLDδh·D=δvU，D（20）

c.姿态误差方程

对式（15）的姿态更新方程等式两边同时求微分，可得姿态误差方程：

φ·nD=δωnie，D－ω^nin，D×φnD+δωnen，D－C^nb，Dδωbib（21）

选取SINS的姿态失准角φs、速度误差δvs、位置误差δps，陀螺/多普勒测速仪的姿态失准角φD、位置误差δpD，陀螺的常值漂移εb，加速度计的零偏Δb，多普勒测速仪的俯仰和方位安装误差角α′=［αθαφ］T，标度因数误差δk以及重力扰动及其一阶导数δg和δg·作为状态量，则状态量共30维，即

X=［φs，δvs，δps，φD，δpD，εb，Δb，α′0，δk，δg，δg·］T（22）

将式（16）～（21）所示4部分的状态模型联立，可得滤波器的状态模型：

X·（t）=F（t）X（t）+G（t）W（t）（23）

式中：F（t）为系统矩阵，F（t）=F1F206×24F3；F1与SINS误差、陀螺/LDV导航系统误差和各传感器误差有关，可以表示为F1=MsMD09×24；

Ms=Mφφ，sMφv，sMφp，s03×6－C^nb，s03×303×3

Mvφ，sMvv，sMvp，s03×603×3C^nb，s03×3

03×3Mpv，sMpp，s03×603×303×303×3，

Mφv，s=0－1RM+h^s01RM+h^s000tanL^sRM+h^s0，

Mφp，s=00v^N，s（RN+h^s）2

－ω^ie，ssinL^s0－v^E，s（RN+h^s）2

ω^ie，scosL^s+v^E，sRN+h^ssec2L^s0－v^E，stanL^s（RN+h^s）2，

Mvp，s=Av，s（Mω，s+Mφp，s），Mvv，s=Av，sMφv，s－Aω，s，

Mω，s=000

－ω^ie，ssinL^s00

ω^ie，scosL^s00，

Mpv，s=01RM+h^s0

secL^sRN+h^s00001，

Mpp，s=00－v^N，s（RM+h^s）2

v^E，ssecL^stanL^sRN+h^s0－v^E，ssecL^sRN+h^s000，

Mφφ，s为－ω^nin，s构成的反对称阵，Mvφ，s为n系下比力f^n构成的反对称阵，Av，s为v^ns构成的反对称阵，Aω，s为（2ω^nie，s+ω^nen，s）构成的反对称阵；

MD=03×9Mφφ，DMφp，D－C^nb，D03×3Mφk，D

03×9Mpφ，DMpp，D03×303×3Mpk，D，

Mφφ，D=Mφφ1，D－Av，DMφv，D，

Mφv，D=0－1RM+h^D0

1RN+h^D00

tanL^DRN+h^D00，

Mφφ1，D为v^N，DRM+h^D

－ω^ie，DcosL^D－v^E，DRN+h^D

－ω^ie，DsinL^D－v^E，DtanL^DRN+h^D所构成的反对称阵，

Mφp，D=000－ω^ie，DsinL^D00

ω^ie，DcosL^D+v^E，Dsec2L^DRN00，

Mφk，D=c23RM+hD－c21RM+hD－c22RM+hD－c13RN+hDc11RN+hDc12RN+hD

－c13tanLDRM+hDc11tanLDRM+hDc12tanLDRM+hD，

Mpφ，D=－Mpv，DAv，D，

Mpv，D=01RM+h^D0

secL^DRN+h^D00

001，

Mpp，D=000v^nE，DsecL^DtanL^DRN+h^D00000，

Mpk，D=vD－c23RM+hDc21RM+hDc22RM+hD－c13secLDRN+hDc11secLDRN+hDc12secLDRN+hD－c33c31c32，

Av，D为v^nD构成的反对称阵，cij为姿态矩阵C^nb，D的第i行、第j列元素；F2和F3与重力扰动有关，

F2=03×303×3－I3×303×3018×3018×3，F3=03×3I3×3MgMdg，

Mg=－ω20，E000－ω20，N000－ω20，U，

Mdg=－2ξω20，E000－2ξω20，N000－2ξω20，U；

W（t）为状态噪声，包含了陀螺和加速度计的测量白噪聲wg，wa，以及重力扰动马尔科夫模型的驱动噪声q，可表示为W（t）=wgwaqT；G（t）为状态噪声系数矩阵，G（t）=－C^nb，s03×303×3

03×3C^nb，s03×3

03×303×303×3

－C^nb，D03×303×3

015×3015×3015×3

03×303×3I3×3。

综合以上各式可得重力扰动补偿系统的状态模型。

3.2量测模型

（1）位置量测模型

以SINS与陀螺/LDV系统解算所得位置之差Zp作为量测值，则可建立位置误差量测模型：

Zp=p^ns－p^nD=（pn+δpns）－（pn+δpnD）=

δpns－δpnD=HpX+Vp（24）

式中：Hp=［03×6I3×303×3－I3×303×15］；Vp为虚拟位置量测噪声。

（2）速度量测模型

以SINS与陀螺/LDV系统解算所得速度之差Zv作为量测值，则可建立速度误差量测模型：

Zv=v^ns－v^nD=（vn+δvns）－（vn+δvnD）=

δvns－δvnD=HvX+Vv（25）

将式（19）陀螺/激光多普勒测速系统速度误差方程代入式（25）可得速度量测矩阵的具体形式：

Hv=［03×3I3×303×3H103×9H203×6］

式中：H1为－v^nD构成的反对称矩阵；H2的具体形式为H2=vD－c13c11－c12－c23c21－c22－c33c31－c32；Vv为激光多普勒测速仪的速度量测噪声。

（3）姿态量测模型

将SINS与陀螺/LDV系统解算所得捷联矩阵C^nb，s转置与C^nb，D相乘，构造得到姿态误差阵Cdφ：

Cdφ=C^nb，D（C^nb，s）T=（I－ΦnD）Cnb［（I－Φns）Cnb］T=

I+Φns－ΦnD（26）

式中：Φns和ΦnD分別为φns和φnD所构成的反对称阵。

以Cdφ的非对角线元素Zφ=φns－φnD作为量测值，可以建立姿态量测模型：

Zφ=HφX+Vφ（27）

式中：Hφ=［I3×303×6－I3×303×18］；Vφ为虚拟姿态量测噪声。

（4）高度量测模型

将气压高度计测量所得hBA与惯导解算所得高度h^s之差作为量测值，则可建立高度量测模型：

Zh=HhX+Vh（28）

式中：Hh=［01×2101×27］；Vh为高度计量测噪声。

联立式（24）～（28）可得组合滤波器的量测模型：

Z=HX+V（29）

式中：Z为量测值；H为观测矩阵；V为系统噪声阵，具体可写为

Z=ZpZvZφZh，H=HpHvHφHh，V=VpVvVφVh。

以上给出了重力扰动补偿系统滤波器的状态方程和量测方程，根据这两组方程再加上必要的初始条件即可进行卡尔曼滤波，从而有效地估计出惯导系统误差以及重力扰动矢量信息。

3.3系统可观性分析

系统的可观性代表由系统输出或观测得到估计状态的可能性。系统能否利用量测信息有效地对状态量进行估计，主要由系统的可观性决定。因此，可观性分析是判断系统有效性的重要环节。

由于载体的运动，重力扰动补偿系统是一个时变的系统，在滤波过程中状态变量的可观测性不恒定，给系统的可观性分析带来了困难。针对这一问题，可以利用分段定常（PWCS）的方法对系统进行处理［21-22］。此处以所建立的系统模型为基础，利用分段线性定常可观测性分析理论对典型机动方式下系统的可观性进行分析，得到系统可观性矩阵的秩如表1所示，表中“#”代表不同的时段。

由表1结果可知，载体机动方式对重力扰动补偿系统的可观性有一定的影响。当载体以匀速直线和加速直线运动时，只有线速度或线加速度发生变化，不利于状态量之间的解耦；当载体进行转弯和爬升等角机动时，线速度/线加速度发生变化的同时，方向余弦和角速度也发生变化，有利于状态量之间的解耦，因此系统的可观性较好。

此外，当载体进行俯冲与爬升后，系统完全可观，因此在实际应用过程中，可先进行俯冲与爬升机动以改善系统的可观性，再对重力扰动进行补偿。

4仿真验证

4.1仿真条件

运载体从30°N，120°E出发，初始速度与姿态角均为0，行驶过程中进行直线行驶、加减速、转弯、爬坡等机动。陀螺仪的常值漂移为0.01（°）/h，量测噪声标准差为0.005（°）/h；加速度计常值偏置为10×10-6g，零偏白噪声为5×10-6g，三向初始位置误差为10m，三向初始速度误差为0.1m/s，三向初始对准误差均为30″。LDV标度因数误差为0.01，量测噪声为0.01m/s；气压高度计量测噪声为20m。IMU数据更新率为100Hz，LDV和气压高度计测量频率均为10Hz。仿真总时间为1h，仿真步长为0.01s。

在仿真中使用的重力扰动数据由GGMplus高分地球重力数据库提供。GGMplus是美国宇航局与德国航天局联合研制的重力数据库［23］，其所包含的重力场数据及其描述性统计结果如表2所示。

该数据库综合利用了重力卫星数据、EGM重力场球谐模型与地形重力数据，能提供分辨率为7.2″的地球重力场网格信息，是现阶段分辨率和精度最高的区域重力数据库，能有效反映实际重力扰动的分布情况。基于该模型提供的重力扰动进行仿真可有效验证所设计方案的性能。由GGMplus数据库得到的载体运动轨迹上重力扰动如图7所示。

由图7可知，仿真轨迹上的重力扰动以近似于常值的长波分量为主，并且包含局部变化较大的短波分量，与实际重力扰动特性相符。

4.2仿真结果与分析

4.2.1重力扰动的估计结果

为了验证所提方案对重力扰动的估计效果，分别采用基于EGM球谐模型补偿的方案（传统方法1）、基于状态估计补偿的方案（传统方法2）与所提方法对重力扰动进行估计，得到相应垂线偏差和重力异常的估计误差曲线，如图8所示。

由图8对重力扰动的估计误差可知，直接利用EGM2008球谐模型对重力扰动的估计精度较差，对垂线偏差的估计精度约为7.17″，对重力异常的估计精度约为8.52mGal。此外，由于球谐模型仅能计算较长波段的重力扰动分量，对重力场变化的细节不敏感，因此在重力扰动变化较为剧烈时的估计精度会进一步下降。

在状态估计补偿方案中，由于重力扰动马尔科夫模型所能覆盖的波段有限，对长波扰动分量的描述精度较差，导致该方案对近似于常值的长波重力扰动分量估计效果不佳，即估计结果与实际值相比存在约5mGal的常值偏置，较难适用于全波段重力扰动的估计。该方案对垂线偏差的估计精度约为6.22″，对重力异常的估计精度约为7.83mGal，整体上估计精度有限。

所设计方案经过约700s的运动后，对重力扰动的估计误差逐渐收敛，在实际重力扰动变化较为剧烈时也能维持较好的估计效果，对垂线偏差的估计精度约为3.14″，对重力异常的估计精度约为4.20mGal。由此可知，通过利用EGM模型预先对长波重力扰动分量进行补偿，可有效改善马尔科夫模型对重力扰动残差描述的准确度，进而提高其估计精度，从而能够对长波以及中短波重力扰动分量进行有效的估计。此外，由于低阶球谐模型的计算结果接近实际重力扰动的常值分量，经过补偿后可减小重力扰动残差的常值偏置，进而减小滤波估计的初始误差，因此当载体以同样规律进行机动时，所设计方案对重力扰动估计误差的收敛速度更快。

4.2.2导航解算结果

将未补偿重力扰动的纯惯导解算结果、利用球谐模型补偿后的导航结果、利用状态估计法补偿后的导航结果以及利用所提方法补偿重力扰动后的导航结果进行對比，其位置误差曲线如图9所示。

由图9可知，当不对重力扰动进行补偿时，导航位置误差随时间迅速发散。利用球谐模型对重力扰动进行补偿后，能一定程度减缓导航误差的发散。由于状态估计法在对重力扰动估计与补偿的同时，能对惯导系统误差进行反馈矫正，因此传统方法2和所提方法的导航误差较小。不过由于SINS和LDV航位推算均具有相同的初始误差，而气压高度计只能提供高度方向的阻尼，系统没有水平方向的绝对位置量测值，导致水平方向的位置误差不收敛，东向和北向的位置误差始终在初始值附近震荡。

将20次仿真结果取平均，得到这几种方法的平均导航误差，如表3所示。

由表3结果可知，相比于单纯的状态估计法补偿，所提方法通过改善重力扰动的估计与补偿精度，可进一步减小重力扰动对惯导系统的影响，因此能够有效提升导航精度，其水平定位精度约为3m，垂直定位精度约为0.05m，实现了高精度的重力扰动补偿和组合导航。

5结论

实际重力扰动所覆盖的波长范围较宽，传统方法较难同时对各波段的重力扰动分量进行高精度补偿。针对这一问题，提出了一种球谐模型补偿与状态估计补偿相互结合的高精度重力扰动补偿方法，对长波以及中短波重力扰动分量采取两种不同途径进行估计。利用EGM2008球谐模型在长波段计算精度高的特点，对长波扰动分量进行精确的补偿；利用马尔科夫模型对中短波重力扰动分量描述效果较好的特点，对残余的重力扰动分量进行高精度的状态建模与估计，进一步提高对中短波分量的估计精度，整体上具有较好的补偿效果，可实现高精度的重力扰动补偿和组合导航。

参考文献：

［1］WestonJL，TittertonDH.ModernInertialNavigationTechnologyandItsApplication［J］.Electronics&CommunicationEngineeringJournal，2000，12（2）：49-64.

［2］陆志东，王晶.高精度惯性导航系统重力补偿方法［J］.航空科学技术，2016，27（8）：1-6.

LuZhidong，WangJing.GravityCompensationMethodsforHighAccuracyINS［J］.AeronauticalScience＆Technology，2016，27（8）：1-6.（inChinese）

［3］KwonJH.GravityCompensationMethodsforPrecisionINS［C］∥60thAnnualMeetingoftheInstituteofNavigation，2001：483-490.

［4］ChangLB，QinFJ，WuMP.GravityDisturbanceCompensationforInertialNavigationSystem［J］.IEEETransactionsonInstrumentationandMeasurement，2019，68（10）：3751-3765.

［5］王晶，杨功流，李湘云，等.重力扰动矢量对惯导系统影响误差项指标分析［J］.中国惯性技术学报，2016，24（3）：285-290.

WangJing，YangGongliu，LiXiangyun，etal.ErrorIndicatorAnalysisforGravityDisturbingVectorsInfluenceonInertialNavigationSystem［J］.JournalofChineseInertialTechnology，2016，24（3）：285-290.（inChinese）

［6］ZhuZS，TanH，JiaY，etal.ResearchontheGravityDisturbanceCompensationTerminalforHighPrecisionPositionandOrientationSystem［J］.Sensors，2020，20（17）：4932.

［7］丛琳，赵忠，杨小步.长航时捷联惯导重力扰动影响及补偿［J］.计算机工程与应用，2014，50（11）：251-255.

CongLin，ZhaoZhong，YangXiaobu.AnalysisandCompensationofDisturbingGravityEffectonLongEnduranceStrapDownInertialNavigationSystem［J］.ComputerEngineeringandApplications，2014，50（11）：251-255.（inChinese）

［8］WengJ，LiuJN，JiaoMX，etal.AnalysisandOnLineCompensationofGravityDisturbanceinaHighPrecisionInertialNavigationSystem［J］.GPSSolutions，2020，24（3）：1-8.

［9］PavlisNK，HolmesSA，KenyonSC，etal.TheDevelopmentandEvaluationoftheEarthGravitationalModel2008（EGM2008）［J］.JournalofGeophysicalResearch：SolidEarth，2012，117（B4）：531-535.

［10］刘雨生，楼立志.不同分辨率数字高程模型对EGM2008截断误差补偿效果分析［J/OL］.地球物理学进展，2022.https：∥kns-cnki-net.e2.buaa.edu.cn/kcms/detail/11.2982.P.20210723.0858.012.html.

LiuYusheng，LouLizhi.AnalysisofCompensationforTruncationErrorofEGM2008byDigitalElevationModelwithDifferentResolutions［J/OL］.ProgressinGeophysics，2022.https：∥kns-cnki-net.e2.buaa.edu.cn/kcms/detail/11.2982.P.20210723.0858.012.html.（inChinese）

［11］杨金玉，张训华，张菲菲，等.EGM2008地球重力模型数据在中国大陆地区的精度分析［J］.地球物理学进展，2012，27（4）：1298-1306.

YangJinyu，ZhangXunhua，ZhangFeifei，etal.OntheAccuracyofEGM2008EarthGravitationalModelinChineseMainland［J］.ProgressinGeophysics，2012，27（4）：1298-1306.（inChinese）

［12］WangJ，YangGL，LiXY，etal.ApplicationoftheSphericalHarmonicGravityModelinHighPrecisionInertialNavigationSystems［J］.MeasurementScienceandTechnology，2016，27（9）：095103.

［13］JekeliC.AirborneVectorGravimetryUsingPrecise，PositionAidedInertialMeasurementUnits［J］.BulletinGéodésique，1994，69（1）：1-11.

［14］AnW，XuJN，HeHY，etal.AMethodofDeflectionoftheVerticalMeasurementBasedonAttitudeDifferenceCompensation［J］.IEEESensorsJournal，2021，21（12）：13125-13136.

［15］戴东凯.基于GNSS/SINS组合的船载高精度垂線偏差测量方法研究［D］.长沙：国防科学技术大学，2016.

DaiDongkai.ResearchonHighPrecisionShipborneMeasurementMethodofDeflectionsoftheVerticalBasedonGNSS/SINSIntegration［D］.Changsha：NationalUniversityofDefenseTechnology，2016.（inChinese）

［16］赵德军，吴晓平.空间扰动引力的谱分析［J］.海洋测绘，2005，25（2）：6-8.

ZhaoDejun，WuXiaoping.SpectralAnalysisofSpatialDisturbingGravity［J］.HydrographicSurveyingandCharting，2005，25（2）：6-8.（inChinese）

［17］JekeliC.AlgorithmsandPreliminaryExperienceswiththeLN93andLN100forAirborneVectorGravimetry［R］.DefenseTechnicalInformationCenter，1998.

［18］MoritzH.GeodeticReferenceSystem1980［J］.BulletinGéodésique，1980，54（3）：395-405.

［19］铁俊波，吴美平，蔡劭琨，等.基于EGM2008重力场球谐模型的水平重力扰动计算方法［J］.中国惯性技术学报，2017，25（5）：624-629.

TieJunbo，WuMeiping，CaiShaokun，etal.GravityDisturbanceCalculationMethodBasedonEarthGravitationalModel2008［J］.JournalofChineseInertialTechnology，2017，25（5）：624-629.（inChinese）

［20］WangQ，GaoCF，ZhouJ，etal.TwoDimensionalLaserDopplerVelocimeterandItsIntegratedNavigationwithaStrapdownInertialNavigationSystem［J］.AppliedOptics，2018，57（13）：3334-3339.

［21］周卫东，蔡佳楠，孙龙，等.GPS/SINS超紧组合导航系统可观测性分析［J］.北京航空航天大学学报，2013，39（9）：1157-1162.

ZhouWeidong，CaiJianan，SunLong，etal.ObservabilityAnalysisofGPS/SINSUltraTightlyCoupledSystem［J］.JournalofBeijingUniversityofAeronauticsandAstronautics，2013，39（9）：1157-1162.（inChinese）

［22］肖佳敏，朱锋，张小红.GNSS/SINS松组合系统的可观测性分析［J］.导航定位学报，2018，6（4）：35-41.

XiaoJiamin，ZhuFeng，ZhangXiaohong.ObservabilityAnalysisonLooselyCoupledGNSS/SINSIntegrationSystem［J］.JournalofNavigationandPositioning，2018，6（4）：35-41.（inChinese）

［23］HirtC，ClaessensS，FecherT，etal.NewUltrahighResolutionPictureofEarthsGravityField［J］.GeophysicalResearchLetters，2013，40（16）：4279-4283.

[HJ*3][HJ][JZ（]AnAccurateGravityDisturbanceCompensationMethod

BasedonSphericalHarmonicModelandMultiSensorFusion

LiuYuxin1，WangXinlong1*，WangXun2，GaoWenning3，HuXiaodong4

（1.SchoolofAstronautics，BeihangUniversity，Beijing100083，China;

2.BeijingInstituteofAutomaticControlEquipment，Beijing100074，China;

3.BeijingInstituteofSatelliteInformationEngineering，Beijing100095，China；

4.AVICXianFlightAutomaticControlResearchInstitute，Xian710065，China）

Abstract：Withtheimprovementoftheaccuracyrequirementsofhighprecisionlongenduranceinertialnavigationsystem，thegravitydisturbancehasbecomethemainerrorsourceofinertialnavigationsystem.Isthekeyfactortoimprovethenavigationaccuracywhetheritcanbeeffectivelycompensated.Thetraditionalcompensationmethodbasedonsphericalharmonicmodelcannotreflectthedetailedinformationoftheearthsgravityfield，andthushaspoorcompensationresultformediumandshortwavegravitydisturbancecomponents.Inthetraditionalstateestimationmethod，theMarkovstatemodelhaspooraccuracyindescribingthelongwavegravitydisturbancecomponent，sothecompensationresultofthelongwavecomponentispoor.Theabovemethodscannotcompensatetheactualgravitydisturbancethathaswidebandwithhighprecision.Tosolvethisproblem，ahighprecisiongravitydisturbancecompensationmethodbasedonsphericalharmonicmodelandmultisensorinformationfusionisdesignedinthispaper.Ontheonehand，themethodcompensatesthelongwavegravitydisturbancecomponentbyusingthehighcalculationaccuracyofthelowordersphericalharmonicmodelinthelongwaveband.Ontheotherhand，usingstrapdowninertialnavigation/laserDopplervelocimetry/barometricaltimetertoformafullyautonomousintegratednavigationsystem，theresidualmediumandshortwavegravitydisturbancecomponentisestablishedasahighprecisionMarkovmodel，soastorealizethestateestimationandcompensationofmediumandshortwavecomponent.Thesimulationresultsshowthattheproposedgravitydisturbancecompensationmethodcaneffectivelyimprovetheestimationeffectofgravitydisturbance，realizehighprecisiongravitydisturbancecompensationandhavehighnavigationaccuracy.

Keywords：strapdowninertialnavigationsystem;gravitydisturbancecompensation;sphericalharmonicgravitymodel;stateestimation;integratednavigation

收稿日期：2022-03-01

基金項目：国家自然科学基金项目（61673040）；重点基础研究项目（2020-JCJQ-ZD-136-12）；航空科学基金项目（20170151002）；天地一体化信息技术国家重点实验室基金项目（2015-SGIIT-KFJJ-DH-01）

作者简介：刘宇鑫（1999-），男，北京人，硕士研究生。

*通信作者：王新龙（1969-），男，陕西渭南人，教授。