0 引　言

伴随构建人类命运共同体进程的不断推进，北极的气候变化、环境、科研、航道利用、资源勘探与开发、安全、国际治理等问题关系到世界各国和人类的共同生存与发展。当下中国的北极活动也已经由单纯的科学研究逐步拓展至涉及全球治理、区域合作、多边和双边机制等多个层面^[1]。作为确保极区科研探索活动安全有序进行的重要技术保障，高精度、稳定可靠的导航系统必不可少。

惯性导航系统（Inertial Navigation System，INS）当地水平固定指北机械编排由于极区经线快速收敛存在北向难以确定、极点北向缺失等问题^{[2 – 3]}，使用横向惯导力学编排^{[4 – 5]}或格网惯导力学编排^{[6 – 7]}能够有效地解决上述问题，但INS的导航误差会随着时间不断累积，通常需要其他导航系统提供外部信息辅助误差修正。面向极区水下环境的导航定位需求，来自中高轨道的全球卫星导航系统（Global Navigation Satellite System，GNSS）的电磁波信号在水体中严重衰减^[8]，无法支撑水下导航定位。水声导航采用的声信号在水中衰减率小、传播性能好，更能满足水下导航对信号传输的实际需求。超短基线（Ultra-Short Baseline，USBL）水声导航采用的基阵规模小、便于载体携带，可稳定提供载体位置信息，故而被广泛应用与INS形成组合导航系统，抑制INS误差发散^{[9 – 10]}。常规的USBL定位方式^{[11 – 12]}需要提前在海底布置应答器，并预先准确测量出应答器的位置。这使得USBL的适用海域十分受限，对于不常涉足且无应答器布设的海域，仅为一次航行布设海底应答器成本过高，实际应用价值大打折扣。因此，基于GNSS浮标的USBL定位方法可用于解决此类陌生海域无应答器导致USBL不可用的问题。

汤均博等^[13]提出在海面安置GNSS浮标实时测定瞬时海面的高度，将实时水位信息发送给用户并结合电子海图实现导航的方法，最高水深测量精度能达到标准偏差3 cm。赵韵^[14]针对USBL声学观测不稳定的问题，采用新息匹配的方法检测异常点，利用USBL预测模型来修正INS误差，东向、北向位置精度分别提高45.0%和55.0%。WANG等^[15]提出一种将USBL方位角、倾斜距离和高度的原始信息作为测量信息引入的组合方式，与传统方法相比水平位置均方根误差分别降低了21.6%和13.2%。CHENG等^[16]提出一种考虑极地声通信延迟的横向SINS/USBL组合导航算法，通过设计自适应滤波的方法同时减弱水声通信延迟和极区水下干扰对组合导航性能的影响，在横北、横东、地3个方向上定位精度分别提升了58.9%、69.2%和87.2%。

上述基于GNSS浮标的INS/USBL组合导航模型将GNSS解算得到的浮标位置直接应用于组合导航系统，均未考虑浮标位置误差对组合导航结果的影响，而在实际应用中，极区GNSS卫星低高度角、电离层闪烁和磁暴频发等弱观测环境严重降低浮标定位准确性，而浮标定位的误差又会通过USBL测量最终传递至载体，最终导致组合导航结果受浮标的GNSS定位误差影响。

针对浮标定位误差对基于GNSS浮标的INS/USBL组合导航模型的影响，本文首先介绍了基于GNSS浮标的INS/USBL组合导航系统组成，构建横向INS/USBL组合导航模型，在此基础上分析浮标GNSS位置误差对组合导航结果的影响，提出一种顾及浮标GNSS位置误差的横向INS/USBL组合模型，并采用极区仿真实验验证其有效性。

1 系统组成

本文基于浮标的极区横向INS/USBL组合导航系统如图1所示，由一个水面浮标和自主无人水下潜航器（Autonomous Underwater Vehicles，AUV）2部分组成，其中浮标搭载GNSS定位系统和应答器，AUV搭载换能器、声学基阵和捷联惯导系统。

图1所示系统中，水面浮标接收GNSS信号，实时解算出浮标位置信息，通过水声通信将位置信息传送至AUV。USBL应答器放置于浮标下端，与浮标之间刚性连接，应答器与GNSS接收端相对位置可准确测量，由此得到应答器位置。当水下AUV需要导航辅助信息时，通过AUV上端的换能器发出信号，声学基阵接收应答器的声波，并由声波信息解算出AUV与浮标的相对位置，作为修正惯性导航系统误差累积的辅助信息，如图2所示。

2 横向INS/USBL组合导航模型

本文使用了基于WGS-84参考旋转圆球模型的横向坐标系^[4]，利用横向惯导的基本方程推导出误差方程，进而建立卡尔曼滤波的状态方程及观测方程。

2.1 横向惯导误差方程

设传统地球坐标系为e系、横向地理坐标系为t系、横向地球坐标系为E系、横向地理坐标系为H系，载体系为b系，地球半径为R。将横向坐标系作为导航坐标系的捷联惯性导航系统的方向余弦矩阵、速度微分方程和位置微分方程分别为^[3]：

$ \left\{ {\begin{array}{*{20}{l}} {{{\dot {\boldsymbol C}}}_{{b}}^{{H}} = {\boldsymbol{C}}_{{b}}^{{H}}({\boldsymbol{\omega }}_{{{Hb}}}^{{b}} \times )} ，\\ {{{{\boldsymbol{\dot v}}}^{{H}}} = {\boldsymbol{C}}_{{b}}^{{H}}{{\boldsymbol{f}}^{{b}}} - (2{\boldsymbol{C}}_{{E}}^{{H}}{\boldsymbol{\omega }}_{{{iE}}}^{{E}} + {\boldsymbol{\omega }}_{{{EH}}}^{{H}}) \times {{\boldsymbol{v}}^{{H}}} + {{\boldsymbol{g}}^{{H}}}}，\\ {{{{\boldsymbol{\dot r}}}^{{H}}} = {{\boldsymbol{v}}^{{H}}} - {\boldsymbol{\omega }}_{{{EH}}}^{{H}} \times {{\boldsymbol{r}}^{{H}}}} 。\end{array}} \right. $

(1)

式中：$ {\boldsymbol{C}}_{{b}}^{{H}} $为b系到H系的方向余弦矩阵；$ {\boldsymbol{C}}_{{E}}^{{H}} $为E系到H系的方向余弦矩阵；$ {\boldsymbol{\omega }}_{{{Hb}}}^{{b}} \times $为b系相对于H系的角速度在b系下的投影；$ {{\boldsymbol{v}}^{{H}}} $为船速在H系的投影；$ {{\boldsymbol{f}}^{{b}}} $为惯导的比力输出；$ {\boldsymbol{\omega }}_{{{iE}}}^{{E}} $为地球自转角速度在E系下的投影；$ {\boldsymbol{\omega }}_{{{EH}}}^{{H}} $为H系相对E系的角速度在H系下的投影；$ {{\boldsymbol{g}}^{{H}}} $为重力矢量在H系下的投影；$ {{\boldsymbol{r}}^{{H}}} $为H系下的位置矢量。

由姿态、速度、位置微分方程推出对应的误差方程为^[17]：

$ \left\{ {\begin{array}{*{20}{l}} {\dot \phi = - ({\boldsymbol{\omega }}_{{{iE}}}^{{H}} + {\boldsymbol{\omega }}_{{{EH}}}^{{H}}) \times \phi + \delta {\boldsymbol{\omega }}_{{{iE}}}^{{H}} + \delta {\boldsymbol{\omega }}_{{{EH}}}^{{H}} - {\boldsymbol{C}}_{{b}}^{{H}}\delta {\boldsymbol{\omega }}_{{{ib}}}^{{b}}} ，\\ \delta {{{\boldsymbol{\dot v}}}^{{H}}} = {{\boldsymbol{f}}^{{H}}} \times \phi - (2{\boldsymbol{\omega }}_{{{iE}}}^{{H}} + {\boldsymbol{\omega }}_{{{EH}}}^{{H}})\times\\\quad\quad\quad \delta {{\boldsymbol{v}}^{{H}}} + {{\boldsymbol{v}}^{{H}}} \times (2\delta {\boldsymbol{\omega }}_{{{iE}}}^{{H}} + \delta {\boldsymbol{\omega }}_{{{EH}}}^{{H}}) + {\boldsymbol{C}}_{{b}}^{{H}}{{\boldsymbol{\xi }}^{{b}}} ，\\ {\delta {\boldsymbol{\dot p}} = {{\boldsymbol{F}}_{{{pv}}}}\delta {{\boldsymbol{v}}^{{H}}} + {{\boldsymbol{F}}_{{{pp}}}}\delta {\boldsymbol{p}}} 。\end{array}} \right. $

(2)

式中：$ {\boldsymbol{C}}_{{b}}^{{H}} $为b系到H系的方向余弦矩阵；$ \delta {\boldsymbol{\omega }}_{{{ib}}}^{{b}} $为陀螺的输出误差；$ {{\boldsymbol{f}}^{{H}}} $为H系下的比力向量；$ {{{\xi }}^{\text{b}}} $为加速度计的输出误差。

2.2 状态方程

横向INS/USBL组合导航模型以横向导航系统姿态误差、横向速度误差、横向经纬高位置误差、及惯性器件漂移误差为状态量，构成15维状态向量：

$ {{\boldsymbol{X}} = {\left[ {\begin{array}{*{20}{l}} {{\phi _{\mathit{x}}}}&{{\phi _{\mathit{y}}}}&{{\phi _{\mathit{z}}}}&{\delta {v_{\mathit{x}}}}&{\delta {v_{\mathit{y}}}}&{\delta {v_{\mathit{z}}}}&{\delta \varphi }&{\delta \lambda }&\\ {\delta h}& {{\varepsilon _{\mathit{x}}}}&{{\varepsilon _{\mathit{y}}}}&{{\varepsilon _{\mathit{z}}}}&{{\xi _{\mathit{x}}}}&{{\xi _{\mathit{y}}}}&{{\xi _{\mathit{z}}}} \end{array}} \right]^{\text{T}}}。}$

(3)

式中：$ {\phi _i}(i = {x,y,z}) $为姿态误差；$ \delta {v_i}(i = {x,y,z}) $为横向速度误差；$ \delta \varphi $、$ \delta \lambda $、$ \delta h $为横向经纬度和高度位置误差；$ {\varepsilon _i}(i = {{x,y,z}}) $为陀螺常值漂移；$ {\xi _i}(i = {{x,y,z}}) $为加速度计零偏。

状态空间模型可以表示为：

$ {{\boldsymbol{\dot X}}_k} = {{\boldsymbol{F}}_k}{{\boldsymbol{X}}_k} + {{\boldsymbol{G}}_k}{{\boldsymbol{W}}_k} 。$

(4)

式中：$ {{\boldsymbol{F}}_k} $为状态转移矩阵；$ {{\boldsymbol{G}}_k} $为系统噪声矩阵；$ {{\boldsymbol{W}}_k} $为系统噪声^[17]。

2.3 观测方程

横向INS/USBL组合导航模型采用载体和信标之间的方向角与斜距为位置量测：

$ {{\boldsymbol{Z}}_k} = \left[ {\begin{array}{*{20}{c}} {{\alpha _{{\text{S,}}k}} - {\alpha _{{\text{U,}}k}}} \\ {{\beta _{{\text{S,}}k}} - {\beta _{{\text{U,}}k}}} \\ {{R_{{\text{S,}}k}} - {R_{{\text{U,}}k}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\delta {\alpha _{{\text{S,}}k}} - \delta {\alpha _{{\text{U,}}k}}} \\ {\delta {\beta _{{\text{S,}}k}} - \delta {\beta _{{\text{U,}}k}}} \\ {\delta {R_{{\text{S,}}k}} - \delta {R_{{\text{U,}}k}}} \end{array}} \right] 。$

(5)

式中：$ {\left[ {{\alpha _{{\text{U,}}k}}}\;\;\;{{\beta _{{\text{U,}}k}}}\;\;\;{{R_{{\text{U,}}k}}} \right]^{\text{T}}} $和$ {\left[ {{\alpha _{{\text{S,}}k}}}\;\;\;{{\beta _{{\text{S,}}k}}}\;\;\;{{R_{{\text{S,}}k}}} \right]^{\text{T}}} $分别为第k个历元的超短基线和横向惯性导航提供的相对位置信息；$ {\left[ {\delta {\alpha _{{\text{U,}}k}}}\;\;\;{\delta {\beta _{{\text{U,}}k}}}\;\;\;{\delta {R_{{\text{U,}}k}}} \right]^{\text{T}}} $和$ \left[ {\delta {\alpha _{{\text{S,}}k}}}\;\;\;{\delta {\beta _{{\text{S,}}k}}} \;\;\; {\delta {R_{{\text{S,}}k}}} \right]^{\text{T}} $分别为第k个历元的超短基线与横向惯性导航相对位置误差；$ \alpha $为载体系下斜距方向与oxy平面的夹角；$ \beta $为载体系下斜距在oxy平面的投影与y轴的夹角；$ R $为斜距。建立量测模型如下：

$ {{\boldsymbol{Z}}_k} = {{\boldsymbol{H}}_k}{{\boldsymbol{X}}_k} + {{\boldsymbol{V}}_k} 。$

(6)

式中：$ {{\boldsymbol{H}}_k} $为设计矩阵^[18]：

$ {{\boldsymbol{H}}_k} = {\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{H}}_{\phi {\text{,}}k}}}&{{0_{3 \times 3}}}&{{{\boldsymbol{H}}_{{\delta p}{\text{,}}k}}}&{{0_{3 \times 3}}}&{{0_{3 \times 3}}} \end{array}} \right]^{\text{T}}} ，$

(7)

$ {{\boldsymbol{H}}_{\phi {\text{,}}k}} = - {{\boldsymbol{A}}_k}{\boldsymbol{C}}_{{{H,}}{k}}^{{u}}{{\boldsymbol{H}}_{1{\text{,}}k}} ，$

(8)

$ {{\boldsymbol{H}}_{{\delta p}{\text{,}}k}} = - {{\boldsymbol{A}}_k}{\boldsymbol{C}}_{{{H,}}{k}}^{{u}}({{\boldsymbol{H}}_{2{\text{,}}k}} - {\boldsymbol{C}}_{{{E,}}{k}}^{{H}}{{\boldsymbol{H}}_{3{\text{,}}k}})，$

(9)

$ {{\boldsymbol{H}}_{{\text{1,}}k}} = (({\boldsymbol{C}}_{{{E,}}{k}}^{{H}}{\boldsymbol{P}}_{{{br,}}{k}}^{{E}}) \times ) 。$

(10)

式中：$ {\boldsymbol{P}}_{{{br,}}{k}}^{{E}} $为第k个历元载体和浮标在E系下的坐标差，$ {\boldsymbol{P}}_{{{br,}}{k}}^{{E}} = {\left[ {x_{{{br,}}{k}}^{{E}}}\;\;\;{y_{{{br,}}{k}}^{{E}}}\;\;\;{z_{{{br,}}{k}}^{{E}}} \right]^{\text{T}}} $；$ {{\boldsymbol{H}}_{2{\text{,}}k}} $及$ {{\boldsymbol{H}}_{3{\text{,}}k}} $分别为第k个历元横向速度误差和横向经纬高位置误差向方位角及斜距$ \alpha $、$ \beta $、$ R $转换的矩阵，可写为如下形式：

$ {{\boldsymbol{H}}_{2{\text{,}}k}} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{M}}_{{\text{1,}}k}}}&{{{\boldsymbol{M}}_{{\text{2,}}k}}}&{{0_{3 \times 1}}} \end{array}} \right] ，$

(11)

$ {{{\boldsymbol{M}}_{{\text{1,}}k}} = \left[ {\begin{array}{*{20}{c}} 0 \\ { - x_{{{br,}}{k}}^{{E}}\cos {{\varphi }_k}\cos {{\lambda }_k} - y_{{{br,}}{k}}^{{E}}\cos {{\varphi }_k}\sin {{\lambda }_k} - z_{{{br,}}{k}}^{{E}}\sin {{\varphi }_k}} \\ { - x_{{{br,}}{k}}^{{E}}\sin {{\varphi }_k}\cos {{\lambda }_k} - y_{{{br,}}{k}}^{{E}}\sin {{\varphi }_k}\sin {{\lambda }_k} + z_{{{br,}}{k}}^{{E}}\cos {{\varphi }_k}} \end{array}} \right] ，}$

(12)

$ {{\boldsymbol{M}}_{{\text{2,}}k}} = \left[ {\begin{array}{*{20}{c}} { - x_{{{br,}}{k}}^{{E}}\cos {{\lambda }_k} - y_{{{br,}}{k}}^{{E}}\sin {{\varphi }_k}} \\ {x_{{{br,}}{k}}^{{E}}\sin {{\varphi }_k}\sin {{\lambda }_k} - y_{{{br,}}{k}}^{{E}}\sin {{\varphi }_k}\cos {{\lambda }_k}} \\ { - x_{{{br,}}{k}}^{{E}}\cos {{\varphi }_k}\sin {{\lambda }_k} + y_{{{br,}}{k}}^{{E}}\cos {{\varphi }_k}\cos {{\lambda }_k}} \end{array}} \right] ，$

(13)

${ {{\boldsymbol{H}}_{{\text{3,}}k}} = \left[ {\begin{array}{*{20}{c}} { - {\text{(R}} + {h_k})\sin {{\varphi }_k}\cos {{\lambda }_k}}&{ - {\text{(R}} + {h_k})\cos {{\varphi }_k}\sin {{\lambda }_k}}&{\cos {{\varphi }_k}\cos {{\lambda }_k}} \\ { - {\text{(R}} + {h_k})\sin {{\varphi }_k}\sin {{\lambda }_k}}&{{\text{(R}} + {h_k})\cos {{\varphi }_k}\cos {{\lambda }_k}}&{\cos {{\varphi }_k}\sin {{\lambda }_k}} \\ {{\text{(R}} + {h_k})\cos {{\varphi }_k}}&0&{\sin {{\varphi }_k}} \end{array}} \right] }。$

(14)

设载体在基阵坐标系下的位置坐标为$ \left[ {{x_{{\text{US,}}k}}} \; {{y_{{\text{US,}}k}}} {{z_{{\text{US,}}k}}} \right] $，$ {{R}_{{\text{US,}}k}} = \sqrt {x_{{\text{US,}}k}^2 + y_{{\text{US,}}k}^2 + z_{{\text{US,}}k}^2} $。

$ {{{\boldsymbol{A}}_k} = \left[ {\begin{array}{*{20}{c}} {\dfrac{{ - {x_{{\text{US,}}k}}{z_{{\text{US,}}k}}}}{{{R}_{{\text{US,}}k}^2\sqrt {x_{{\text{US,}}k}^2 + y_{{\text{US,}}k}^2} }}}&{\dfrac{{ - {x_{{\text{US,}}k}}{y_{{\text{US,}}k}}}}{{{R}_{{\text{US,}}k}^2\sqrt {x_{{\text{US,}}k}^2 + y_{{\text{US,}}k}^2} }}}&{\dfrac{{\sqrt {x_{{\text{US,}}k}^2 + y_{{\text{US,}}k}^2} }}{{{R}_{{\text{US,}}k}^2}}} \\ {\dfrac{{{y_{{\text{US,}}k}}}}{{\sqrt {x_{{\text{US,}}k}^2 + y_{{\text{US,}}k}^2} }}}&{\dfrac{{ - {x_{{\text{US,}}k}}}}{{\sqrt {x_{{\text{US,}}k}^2 + y_{{\text{US,}}k}^2} }}}&0 \\ {\dfrac{{{x_{{\text{US,}}k}}}}{{{{R}_{{\text{US,}}k}}}}}&{\dfrac{{{y_{{\text{US,}}k}}}}{{{{R}_{{\text{US,}}k}}}}}&{\dfrac{{{z_{{\text{US,}}k}}}}{{{{R}_{{\text{US,}}k}}}}} \end{array}} \right]} 。$

(15)

$ {{\boldsymbol{V}}_k} $为量测噪声：

$ {{\boldsymbol{V}}_k} = {{\boldsymbol{V}}_{{\text{U}},k}} + {{\boldsymbol{V}}_{{\text{G,}}k}}。$

(16)

式中：$ {{\boldsymbol{V}}_{{\text{U}},k}} $为第k个历元USBL量测噪声；$ {{\boldsymbol{V}}_{{\text{G,}}k}} $为GNSS位置误差产生的量测噪声。

$ {{\boldsymbol{V}}_{{\text{G,}}k}} = {{\boldsymbol{A}}_k}{\boldsymbol{C}}_{{\text{e,}}k}^{\text{u}}{\boldsymbol{\delta }}{{\boldsymbol{r}}_{{\text{GNSS,}}k}}，$

(17)

$ {\boldsymbol{\delta }}{{\boldsymbol{r}}_{{\text{GNSS,}}k}} = \left[ {\begin{array}{*{20}{c}} {\delta {x_{{\text{GNSS,}}k}}} \\ {\delta {y_{{\text{GNSS,}}k}}} \\ {\delta {z_{{\text{GNSS,}}k}}} \end{array}} \right]。$

(18)

式中：$ {{\boldsymbol{A}}_k} $为U系坐标与方向角、斜距的转换矩阵；$ {\boldsymbol{C}}_{{{e,}}k}^{{u}} $为地心地固系与载体系之间的坐标转换矩阵$ {\boldsymbol{C}}_{{{e,}}k}^{{u}} = {\boldsymbol{C}}_{{{b,}}k}^{{u}}{\boldsymbol{C}}_{{{n,}}k}^{{b}}{\boldsymbol{C}}_{{{e,}}k}^{{n}} $；$ {\boldsymbol{\delta }}{{\boldsymbol{r}}_{{\text{GNSS,}}k}} $为GNSS定位结果在地心地固坐标系下的误差。

3 浮标GNSS位置误差对组合导航解的影响

将GNSS定位误差引入USBL观测模型中，进而分析GNSS误差对INS/USBL组合导航解的影响。将GNSS观测误差$ {{\boldsymbol{V}}_{{{G,}}k}} $引入USBL观测误差方程中可得$ {{\boldsymbol{V}}_k} = {{\boldsymbol{V}}_{{{U}},k}} + {{\boldsymbol{V}}_{{{G,}}k}} $，其方差$ {{\boldsymbol{R}}_k} $可写为：

$ {{\boldsymbol{R}}_k} = {{\boldsymbol{R}}_{{{U,}}k}} + {{\boldsymbol{R}}_{{{G,}}k}} 。$

(19)

式中：$ {{\boldsymbol{R}}_{{{U,}}k}} $为USBL自身量测噪声方差；$ {{\boldsymbol{R}}_{{{G,}}k}} $为由GNSS位置误差$ {{\boldsymbol{V}}_{{{G,}}k}} $产生的量测噪声方差：

$ {{\boldsymbol{R}}_{{{G,}}k}} = {{\boldsymbol{A}}_k}{\boldsymbol{C}}_{{{e,}}k}^{{u}}{{\boldsymbol{R}}_{{\text{GNSS,}}k}}{\boldsymbol{C}}_{{{u,}}k}^{{e}}{\boldsymbol{A}}_k^{\text{T}}。$

(20)

式中：$ {{\boldsymbol{R}}_{{\text{GNSS,}}k}} = {\text{diag}}( {\sigma _{{{x\_\text GNSS,}}k}^2}\;\;{\sigma _{{{y\_\text GNSS,}}k}^2} \;\; {\sigma _{{{z\_ GNSS,}}k}^2} ) $为$ {\boldsymbol{\delta }}{{\boldsymbol{r}}_{{\text{GNSS,}}k}} $的方差矩阵。$ \sigma _{{{x\_\text GNSS,}}k}^2 $、$ \sigma _{{{y\_\text GNSS,}}k}^2 $、$ \sigma _{{{z\_\text GNSS,}}k}^2 $分别为GNSS定位结果在e系下x、y、z三轴方向上误差的方差。$ {{\boldsymbol{R}}_{{\text{GNSS,}}k}} $可由经典高度角随机模型结合抗差自适应滤波的后验估计方差阵给出^[19]，可以准确描述当前GNSS定位精度。考虑GNSS误差$ {{\boldsymbol{R}}_{{\text{GNSS,}}k}} $，得到带有浮标误差的第k个历元横向INS/USBL组合导航参数估计结果和方差分别为：

$\begin{split} { {\boldsymbol{X}}_k} = \;& {\boldsymbol{X}}_k^ - + {\boldsymbol{P}}_k^ - {\boldsymbol{H}}_k^{\text{T}}({{\boldsymbol{H}}_k}{\boldsymbol{P}}_k^ - {\boldsymbol{H}}_k^{\text{T}} + {{\boldsymbol{R}}_{{{U,}}k}}+ \\ & {{\boldsymbol{A}}_k}{\boldsymbol{C}}_{{{e,}}k}^{{u}}{{\boldsymbol{R}}_{{\text{GNSS,}}k}}{\boldsymbol{C}}_{{{U,}}k}^{{e}}{\boldsymbol{A}}_k^{\text{T}})^{ - 1}({{\boldsymbol{Z}}_k} - {{\boldsymbol{H}}_k} {\boldsymbol{X}}_k^ - ) \end{split}，$

(21)

$\begin{split} {{\boldsymbol{P}}_k} = \;&{\boldsymbol{P}}_k^ - - {\boldsymbol{P}}_k^ - {\boldsymbol{H}}_k^{\text{T}}({{\boldsymbol{H}}_k}{\boldsymbol{P}}_k^ - {\boldsymbol{H}}_k^{\text{T}} + {{\boldsymbol{R}}_{{{U,}}k}}+ \\& {{\boldsymbol{A}}_k}{\boldsymbol{C}}_{{{e,}}k}^{{u}}{{\boldsymbol{R}}_{{\text{GNSS,}}k}}{\boldsymbol{C}}_{{{U,}}k}^{{e}}{\boldsymbol{A}}_k^{\text{T}})^{ - 1}{{\boldsymbol{H}}_k}{\boldsymbol{P}}_k^ - \end{split}。$

(22)

式中：$ {\boldsymbol{X}}_k^ - $和$ {\boldsymbol{P}}_k^ - $分别为状态量和均方差的先验估计值；$ {{\boldsymbol{A}}_k} $为坐标与方向角、斜距的转换矩阵；$ {{\boldsymbol{R}}_{{\text{GNSS,}}k}} $为一个正定对角矩阵；$ {{\boldsymbol{R}}_{{{G,}}k}} $为$ {{\boldsymbol{R}}_{{\text{GNSS,}}k}} $经过$ {{\boldsymbol{A}}_k} $和$ {\boldsymbol{C}}_{{{e,}}k}^{{u}} $两个矩阵转换后得到的一个非负定阵。由于$ {{\boldsymbol{H}}_k}{\boldsymbol{P}}_k^ - {\boldsymbol{H}}_k^{\text{T}} $和$ {{\boldsymbol{R}}_{{{U,}}k}} $也均为非负定矩阵，随着$ {{\boldsymbol{R}}_{{\text{GNSS,}}k}} $主对角线元素增大，$ {{\boldsymbol{H}}_k}{\boldsymbol{P}}_k^ - {\boldsymbol{H}}_k^{\text{T}} + {{\boldsymbol{R}}_{{{U,}}k}} + {{\boldsymbol{A}}_k}{\boldsymbol{C}}_{{{e,}}k}^{{u}}{{\boldsymbol{R}}_{{\text{GNSS,}}k}}{\boldsymbol{C}}_{{{U,}}k}^{{e}} {\boldsymbol{A}}_k^{\text{T}} $和的主对角线元素增大，$ {\boldsymbol{P}}_k^ - {\boldsymbol{H}}_k^{\text{T}}({{\boldsymbol{H}}_k}{\boldsymbol{P}}_k^ - {\boldsymbol{H}}_k^{\text{T}} + {{\boldsymbol{R}}_{{{U,}}k}} + {{\boldsymbol{A}}_k}{\boldsymbol{C}}_{{{e,}}k}^{{u}} {{\boldsymbol{R}}_{{\text{GNSS,}}k}} {\boldsymbol{C}}_{{{U,}}k}^{{e}} {\boldsymbol{A}}_k^{\text{T}})^{ - 1}{{\boldsymbol{H}}_k}{\boldsymbol{P}}_k^ - $趋向于减小。因此，GNSS定位误差增大时，$ {{\boldsymbol{P}}_k} $主对角线元素趋向于增大，即组合导航解的精度降低。

$ {{\boldsymbol{R}}_{{\text{GNSS,}}k}} $的主对角线元素均为0时，即不考虑GNSS定位误差，$ {{\boldsymbol{P}}_k} $达到最小值，即INS/USBL组合导航解精度最高，此时观测误差来源仅为USBL测量误差。$ {{\boldsymbol{R}}_{{\text{GNSS,}}k}} $的主对角线元素均为∞时，即GNSS定位解完全不可信，$ {{\boldsymbol{P}}_k} = {\boldsymbol{P}}_k^ - $，即USBL的测量结果完全失效，INS/USBL系统解完全依赖于INS自身。在$ {{\boldsymbol{R}}_{{\text{GNSS,}}k}} $由小变大的过程中，$ {{\boldsymbol{P}}_k} $单调递增，INS/USBL组合导航解的定位解的精度逐渐降低。在GNSS误差较小时，USBL定位误差是影响INS/USBL解的主因，此时忽略GNSS定位误差对INS/USBL解的影响。然而，当GNSS定位误差因极区电离层闪烁、地磁暴等弱观测环境因素急剧放大时，忽略GNSS定位误差会导致INS/USBL随机模型严重有偏，进而导致组合导航解误差增大，严重时会导致系统误差发散。因此，需要采用式(19)所示INS/USBL浮标位置误差（GNSS Buoy Position Error，GBPE）组合导航模型，实现的系统对GNSS定位误差的补偿。

4 实验结果及分析

我国地处中低纬度，开展极区实验难度大、成本高，因此对INS/USBL-GBPE组合导航模型采取仿真实验并分析结果。在实验中，各传感器误差如表1所示。

横向惯性导航系统初始姿态误差分别为：纵摇角误差3′、横摇角误差3′和航向角误差6′。载体位置为北纬75.0000°，东经126.6705°，高度−100 m，信标位置为北纬75.0017°，东经126.6767°，高度0 m。

实验中组合解算周期为1 s，实验全程1 h，在载体和浮标均静止的情况下，对GNSS定位结果分别加入e系三轴方向上方差分别为1、25、100 m²的零均值高斯噪声，组合导航的姿态、速度和位置误差输出如图3所示，均方根（Root Mean Square，RMS）误差统计于表3。GNSS定位误差方差为1 m²时INS/USBL组合导航的定位误差相对于方差为25 m²和100 m²时分别下降了61.5%、68.2%、41.7%和81.0%、86.0%、72.5%。综上，INS/USBL组合导航性能与GNSS定位精度直接相关，GNSS定位精度越高，INS/USBL组合导航性能越好。

模拟在航行中由于GNSS观测质量下降导致浮标定位精度急剧下降的情况，令浮标静止，载体匀速直线行驶，惯导解算周期及各传感器误差不变。组合导航的仿真轨迹如图4所示，其中实线段GNSS定位误差方差为1 m²，虚线段为引入GNSS仿真干扰的时段，具体设置如表4所示。图5、图6分别为实验1和实验2设置下，INS/USBL组合导航误差曲线。

图5、图6为在浮标定位精度急剧下降时，所提出INS/USBL-GBPE模型和传统INS/USBL模型性能差异。从图5可知在实验1的情况下，INS/USBL-GBPE模型航向误差不发散，速度误差和位置误差峰值及波动范围在浮标定位精度下降段均明显优于传统模型，精度分别提升80.9%、86.5%、52.4%和41.3%、61.2%、71.6%。从图6可知在实验2的情况下，INS/USBL-GBPE模型的姿态误差不发散，速度和位置误差均收敛于0，传统模型的误差曲线在GNSS受干扰时姿态、速度、位置误差均发散，组合导航结果失效。

5 结　语

针对单浮标极区INS/USBL组合导航系统中未考虑极区恶劣观测环境下浮标GNSS定位误差的问题，通过分析浮标GNSS定位误差对INS/USBL组合导航解的影响，构建INS/USBL-GBPE模型并利用仿真实验验证所提模型的有效性，得到以下结论：

1）在GNSS观测误差$ {{\boldsymbol{V}}_{{\text{G,}}k}} $由小变大的过程中，$ {{\boldsymbol{P}}_k} $主对角线元素单调递增，INS/USBL组合导航解的定位解的精度逐渐降低。INS/USBL组合导航结果显著受GNSS定位误差影响，且误差越大，对INS/USBL组合导航影响越明显。

2）实验表明，GNSS浮标定位精度对INS/USBL组合导航解的影响主要表现在，随着浮标定位精度逐渐下降，纵摇、横摇误差增幅小于1′，相对于GNSS定位误差方差为1 m²时，定位误差方差为25 m²和100 m²时的INS/USBL组合导航的速度误差收敛时间分别增加32.0%、19.7%、61.1%和50.7%、43.3%、61.4%，定位精度分别下降61.5%、68.2%、41.7%和81.0%、86.0%、72.5%。

3）面对浮标定位精度急剧下降的情况，若不及时调节观测误差方差阵，会导致航向误差在浮标定位精度恢复正常后，最终无法重新收敛，当浮标定位误差方差过大时，组合导航定位误差发散，结果失效。