CN105116431A - 一种惯性导航平台和北斗卫星的高精度超紧耦合导航方法 - Google Patents

Abstract

本发明涉及一种惯性导航平台和北斗卫星的高精度超紧耦合导航方法。根据运动物体状态方程及测量方程建立了随机线性定常离散系统力学模型；利用对系统的冗余测量值的一阶、二阶差分序列进行统计分析，实现了基于冗余测量的自适应卡尔曼滤波算法，准确的估计了系统量测噪声，进而自适应调节噪声方差R；根据超紧耦合传统结构，引入神经网络预测模块，用于辅助失锁状态下的北斗卫星载波跟踪环和码跟踪环。采用时滞滤波器用于北斗卫星接收机，显著改善了超紧耦合系统中的环路跟踪误差与更新时间的相关性；本发明属于导航控制技术领域，可应用于各类搭载捷联惯性导航平台的运载体的高精度导航。

Description

一种惯性导航平台和北斗卫星的高精度超紧耦合导航方法

技术领域

本发明涉及一种惯性导航平台和北斗卫星的高精度超紧耦合导航方法，适用于各类搭载捷联惯性导航平台的运载体的高精度导航。

技术背景

当今导航领域的两大基础导航方式是惯性导航系统(INS)和以GPS、BDS等为代表的星基定位系统。同时，这两种导航系统自身又存在较强的互补性，一直都是近年来组合导航研究的热点；惯性导航具有强自主性、不受环境影响、短期内导航精度高等优点；卫星导航具有定位精度高，时间稳定性较好的特点。两者的工作特性恰好在性能上互补，卫星导航可以加快惯性导航的初始对准时间，改善精度的时间漂移性；而惯性导航可以在卫星导航失锁的时候进行独立的导航，继续提供定位信息。当前的典型组合方式有松耦合、紧耦合和超紧耦合三种方式，其中以超紧耦合对于导航性能的提升最为显著。松耦合组合方式下，惯导和卫星独立工作，并将两者导航数据进行处理用于反馈惯导进行修正；紧耦合组合方式下，卫星提供的用于对惯导惯性测量单元(IMU)进行校正的是信息是伪距、载波相位等用于定位解算的原始信息；超紧耦合方式相对于紧耦合还要完成将INS的量测信息反馈给接收机或者直接利用MU量测信息辅助卫星的码跟踪环和载波跟踪环，提高了接收机的跟踪锁定能力。

但是，在整体方向上，超紧耦合多数只局限于基于GPS和SINS相结合的研究；而由于GPS的信号控制并非我国所掌握，极容易出现导航信号限制甚至丧失导航能力的情况；因而将组合导航应用于我国自主研制的北斗系列导航卫星已成为当前的一大主要研究方向，也是本专利的主要研究方面。

组合导航的相关研究，现已有诸多成果。周启帆设计了一种基于冗余测量的自适应卡尔曼滤波算法，并成功应用于GPS/SINS的组合导航，取得了较为显著的性能提升，但只应用于了两系统独立工作的松耦合系统中；张静娴提出了一种基于神经网络的分布式卡尔曼滤波算法，但只用于松耦合，同时也只能对SINS进行位置和速度的校正，具有较强的局限性；于洁、孙熙、宋高顺等都针对超紧耦合进行了不同方向的研究，但最终均存在非线性计算量巨大或者时间相关性强的性能问题。

超紧耦合对于组合导航而言，是一种新概念导航方式，具有较高的实用意义和性能参数，但同样的，其系统复杂度也是最高的。因此，如何在保证超紧耦合的工作效能的基础上，简化计算，提高系统鲁棒性和普适性，就成为当前研究的主要方向和重点难点所在。

发明内容

本发明的技术解决问题是：为了避免在超紧耦合实现过程中的非线性计算量巨大、接收机信号时间相关性以及系统测量误差等存在性能限制的问题，提出了一种基于惯性导航平台和北斗卫星的高精度超紧耦合导航方法。该方法不仅实现了超紧耦合的双向调节功能，同时降低了环路跟踪误差与更新时间的相关性，通过冗余测量修正了SINS/BDS的测量系统误差，并采用神经网络实现系统对于BDS接收机性能的调整，为搭载捷联惯导的运载体导航提供了一种全新的技术途径。

本发明的技术解决方案是：根据运动物体状态方程及测量方程建立随机线性定常离散系统模型；利用对系统的冗余测量值的差分序列进行统计分析，实现基于冗余测量的卡尔曼滤波算法，准确的估计了系统量测噪声并自适应调节捷联惯导(SINS)和北斗卫星(BDS)噪声方差R；根据超紧耦合传统结构，引入神经网络预测模块，用于辅助失锁状态下的北斗卫星载波跟踪环和码跟踪环；在北斗卫星接收机中添加时滞滤波模块，显著改善了超紧耦合系统中的环路跟踪误差与更新时间的相关性。具体包括以下步骤：

(1)根据运动物体状态方程及测量方程建立随机线性定常离散系统模型：

首先构造SINS和BDS对应的状态方程，BDS/SINS的状态分量的转移过程满足如下等式：

X_k，BDS＝Φ_{k，k-1，BDS}X_k-1，BDS+Г_{k，k-1，BDS}W_k-1，BDS

X_k，SINS＝Φ_{k，k-1，SINS}X_k-1，SINS+Г_{k，k-1，SINS}W_k-1，SINS

于是，归结上两式为：

X_k＝Φ_k，k-1X_k-1+Г_k，k-1W_k-1

而后进行测量方程的分析，BDS/SINS的观测值变化规律满足如下等式：

Z_k，SINS＝H_k，SINSX_k，SINS+V_k，SINS

Z_k，BDS＝H_k，BDSX_k，BDS+V_k，BDS

于是，同样归结上面两式为：

Z_k＝H_kX_k+V_k

则最终得到组合导航系统的状态方程和测量方程，以此表征系统的状态变化和观测值变化规律，作为系统模型的构建规律：

$\left\{\begin{array}{c}{X}_k={\mathrm{\Φ}}_{k,k-1}{X}_{k-1}+{\mathrm{\Γ}}_{k,k-1}{W}_{k-1}\\ Z_k=H_k{X}_k+V_k\end{array}\right.$

状态方程式中，X_k为SINS/BDS系统的n维状态变量；Φ_k，k-1为系统的nⅹn维状态转移矩阵；Г_k，k-1为nⅹp维系统过程噪声输入矩阵；W_k-1为p维系统随机过程噪声序列；

测量方程式中，Z_k为测量系统的m维观测序列，即测量值；H_k为mⅹn维观测矩阵；V_k为m维系统随机观测噪声序列；

上面两式中X_k包含的状态分量有：

$X_k={\left[\begin{array}{ccccccccccccccc}{\mathrm{\φ}}_E& {\mathrm{\φ}}_N& {\mathrm{\φ}}_U& {\mathrm{\δv}}_E& {\mathrm{\δv}}_N& {\mathrm{\δv}}_U& \mathrm{\δ}L& \mathrm{\δ}\mathrm{\λ}& \mathrm{\δ}h& {\mathrm{\ϵ}}_{bx}& {\mathrm{\ϵ}}_{by}& {\mathrm{\ϵ}}_{bz}& {\▿}_x& {\▿}_y& {\▿}_z\end{array}\right]}^T$

其中，下标E、N、U分别代表东、北、天三个方向；下标x、y、z分别代表待导航机体的右、前和上轴向；φ为惯导平台误差角；δv_E等表示速度误差；δL、δλ和δh分别表示三个方向的位置误差；ε_bx等表示陀螺漂移；等表示加速度计漂移；

在自适应卡尔曼滤波的关键步骤就是针对测量噪声的准确估计，即需要对V_k的性质和参数进行准确的分析。

(2)基于冗余测量的卡尔曼滤波算法，实现噪声方差的自适应估计与调节：

在系统量测噪声未知的条件下，通过针对系统冗余值的一阶二阶差分分析，获得统计特性，进而自适应调节方差；

首先，两个测量系统的差分序列表达式为：

$\left\{\begin{array}{c}\mathrm{\Δ}{Z}_1=Z_1\left(k\right)-Z_1(k-1)\\ \mathrm{\Δ}{Z}_2=Z_2\left(k\right)-Z_2(k-1)\end{array}\right.$

于是可以得到系统的二阶差分序列：

ΔZ₁-ΔZ₂

通过计算二阶差分序列的自相关可得：

$\begin{array}{c}E\[({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2){({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2)}^T\]\\ =E\{\[\left(V_1\right(k+1)-V_1(k\left)\right)-\left(V_2\right(k+1)-V_2(k\left)\right)\]\\ \·{\[\left(V_1\right(k+1)-V_1(k\left)\right)-\left(V_2\right(k+1)-V_2(k\left)\right)\]}^T\}\end{array}$

其中V₁(k)和V₂(k)为不同测量系统下的观测噪声；

由于V₁(k)和V₂(k)均为零均值白噪声且不具有时间相关性：

$\left\{\begin{array}{c}E\left(V_1\left(k+1\right)V_1^T\left(k\right)\right)=0\\ E\left(V_1\left(k\right)V_2^T\left(k\right)\right)=0\\ E\left(V_2\left(k+1\right)V_2^T\left(k\right)\right)=0\end{array}\right.$

于是，二阶差分序列自相关可以简化为：

$E\[({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2){({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2)}^T\]=2E\left(V_1\right(k\left)V_1^T\right(k\left)\right)+2E\left(V_2\right(k\left)V_2^T\right(k\left)\right)$

即表达为下式：

E[(ΔZ₁-ΔZ₂)(ΔZ₁-ΔZ₂)^T]＝2(R₁+R₂)

计算两个一阶差分序列自相关：

$\begin{array}{c}E\left({\mathrm{\ΔZ}}_1{\mathrm{\ΔZ}}_1^T\right)=E\[\left(Z\right(k)-Z(k-1\left)\right)\left(Z\right(k)-Z^T(k-1\left)\right)\]+E\left(V_1\right(k\left)V_1^T\right(k\left)\right)\\ +E\left(V_1\right(k-1\left)V_1^T\right(k-1\left)\right)\end{array}$

$\begin{array}{c}E\left({\mathrm{\ΔZ}}_2{\mathrm{\ΔZ}}_2^T\right)=E\[\left(Z\right(k)-Z(k-1\left)\right)\left(Z\right(k)-Z^T(k-1\left)\right)\]+E\left(V_2\right(k\left)V_2^T\right(k\left)\right)\\ +E\left(V_2\right(k-1\left)V_2^T\right(k-1\left)\right)\end{array}$

由于V₁(k)和V₂(k)均为零均值白噪声，上两式化简并相减得到：

$E\left({\mathrm{\ΔZ}}_1{\mathrm{\ΔZ}}_1^T\right)-E\left({\mathrm{\ΔZ}}_2{\mathrm{\ΔZ}}_2^T\right)=2(R_1-R_2)$

于是得到噪声方差矩阵：

$\left\{\begin{array}{c}{R}_1\left(k\right)=\frac{E\[({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2){({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2)}^T\]+\left(E\right({\mathrm{\ΔZ}}_1{\mathrm{\ΔZ}}_1^T)-E({\mathrm{\ΔZ}}_2{\mathrm{\ΔZ}}_2^T\left)\right)}{4}\\ R_2\left(k\right)=\frac{E\[({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2){({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2)}^T\]-\left(E\right({\mathrm{\ΔZ}}_1{\mathrm{\ΔZ}}_1^T)-E({\mathrm{\ΔZ}}_2{\mathrm{\ΔZ}}_2^T\left)\right)}{4}\end{array}\right.$

将对应的SINS/BDS系统差分序列代入上式：

$\left\{\begin{array}{c}{R}_{SINS}\left(k\right)=\frac{E{\left({\mathrm{\ΔZ}}_{SB}{{\mathrm{\ΔZ}}_{SB}}^T\right)}_{k-M:k}+(E{\left({\mathrm{\ΔZ}}_{SINS}{\mathrm{\ΔZ}}_{SINS}^T\right)}_{k-M:k}-E{\left({\mathrm{\ΔZ}}_{BDS}{\mathrm{\ΔZ}}_{BDS}^T\right)}_{k-M:k})}{4}\\ R_{BDS}\left(k\right)=\frac{E{\left({\mathrm{\ΔZ}}_{SB}{{\mathrm{\ΔZ}}_{SB}}^T\right)}_{k-M:k}-(E{\left({\mathrm{\ΔZ}}_{SINS}{\mathrm{\ΔZ}}_{SINS}^T\right)}_{k-M:k}-E{\left({\mathrm{\ΔZ}}_{BDS}{\mathrm{\ΔZ}}_{BDS}^T\right)}_{k-M:k})}{4}\end{array}\right.$

其中：

ΔZ_SB＝ΔZ_SINS-ΔZ_BDS

式中，ΔZ_SINS为SINS的测量系统的一阶差分序列；ΔZ_BDS为BDS测量系统的；ΔZ_BDS为SINS和BDS的二阶差分序列；

在实际工程中采用划窗方式保证计算结果的实时性，上式中，k表示当前测量时刻，M为划窗长度；实际选取k-M:k范围内的信息作为参考值，以此保证自适应的实时效应；

(3)引入神经网络预测模块，实现对于BDS的载波跟踪环和码跟踪环的校正：

最基本的RBF神经网络分为输入层输出层和隐层，隐层用于进行从输入空间到隐藏空间的变换，一般情况下隐层空间具有较高维度；输出层作用在输入层的激活模式并提供响应；

依据超紧耦合的经典模式，系统需要在输出导航变量的同时，输出对卫星接收机的辅助变量，通常为伪距和多普勒频移的估计值；本方法中采用的就是将伪距和多普勒频移的辅助变量估计通过神经网络来实现；

神经网络的工作过程分为两个过程，分别是学习过程和输出过程；其中学习过程又分为自组织学习过程和监督学习过程，分别确定隐层函数和输出层权值；

本方法的工作方式是，当BDS和SINS正常工作时，对神经网络进行训练，BDS和SINS的组合导航输出作为目标神经网络的输出，BDS本身的伪距和多普勒频移测量结果作为目标神经网络的输入；训练算法迭代修正网络参数，从而达到误差的最小均方值；当SINS性能恶化，无法针对BDS完成直接辅助的时候，将BDS的伪距和多普勒结果作为训练完成的神经网络输入，从而预测伪距和多普勒频移的误差，通过矢量相加的方式实现辅助功能；

神经网络的主要模块在于通过学习构建隐层系统参数，前馈型神经网络的结构中，迭代过程为求解最优化问题：

构造误差信号能量函数为：

$J\left(\mathrm{\ω}\right)=E\[\frac{1}{2}\underset{k}{\Σ}{e}_k^2\left(n\right)\]$

得到能量信号之后解决最优化问题：

$\underset{\mathrm{\ω}}{min}J\left(\mathrm{\ω}\right)=E\[\frac{1}{2}\underset{k}{\Σ}{e}_k^2\left(n\right)\]$

最终得到满足上式的最优化系统参数ω；

通过数值迭代算法可得系统参数修正：

${\mathrm{\Δ\ω}}_{kj}=-\mathrm{\η}\frac{\∂J\left(\mathrm{\ω}\right)}{\∂\mathrm{\ω}}$

(4)在接收机中添加时滞滤波模块，改善环路跟踪误差与更新时间的相关性：

将脉冲序列与参考命令卷积形成的整形命令作为控制信号，滤掉参考命令中引起系统振动的频率分量，消除二阶震荡系统的震荡残留；

通常接收机中采用典型的二阶PLL(锁相环)的结构，接收机环路的传递函数为：

而这里的相位误差为：

于是，针对上面等式的相位误差传递函数为：

在本设计之中，存在SINS提供的辅助信息用于调整跟踪环性能，所以需要在原有的载体经历动态中考虑额外的影响，加入了前馈环节之后得到如下传递函数：

$H\left(s\right)=\frac{2{\mathrm{\ξ\ω}}_{n}s+{\mathrm{\ω}}_n^2+\frac{{\mathrm{\α}}_{IMU}{s}^2}{s+{\mathrm{\α}}_{IMU}}}{s^2+2{\mathrm{\ξ\ω}}_{n}s+{\mathrm{\ω}}_n^2}$

其中，ξ为二阶系统的阻尼系数，ω_n为自然圆频率，α_IMU为SINS的辅助数据带宽；

此处采用双脉冲零振时滞滤波器(ZV)，当滤波器模型参数与系统本身参数相同，且不存在参数不确定性的时候，滤波器传递函数如下：

$F\left(s\right)=\underset{i=1}{\overset{2}{\Σ}}{A}_i{e}^{-t_{i}s}$

其中，脉冲幅度和时延表达式如下：

$A_1=\frac{1}{1+K}$

$A_2=\frac{K}{1+K}$

t₁＝0

$t_2=\frac{\mathrm{\π}}{{\mathrm{\ω}}_n\sqrt{1-{\mathrm{\ξ}}^2}}$

$K=\frac{-\mathrm{\ξ}\mathrm{\π}}{e_e\sqrt{1-{\mathrm{\ξ}}^2}}$

于是，最终载波环路误差的最终传递函数为：

H_F(s)＝F(s)·H(s)

本发明的原理是：根据运动物体状态方程及测量方程建立了随机线性定常离散系统力学模型，将运载体的机体信息与SINS的导航信息以及BDS提供的导航信息进行综合分析，实现超紧耦合的高精度导航方式；将SINS和BDS的测量系统结果(速度、位置)直接送到滤波器中，通过对测量结果的冗余测量值的一阶、二阶差分序列进行统计分析，实现了基于冗余测量的自适应卡尔曼滤波算法，准确的估计了系统量测噪声的信息，并通过前馈机制将自适应调节噪声方差R传送到测量端，实现测量校正；根据超紧耦合传统结构，引入神经网络预测模块，神经网络与卡尔曼滤波器并立，在SINS正常工作时进行训练，输入信息为BDS本身的伪距和多普勒频移测量结果，同时还有SINS的导航信号，输出为伪距和多普勒频移的误差，用于辅助失锁状态下的北斗卫星载波跟踪环和码跟踪环。采用时滞滤波器用于北斗卫星接收机，对正常接收的I-Q正交信号直接进行时滞滤波，将滤波处理后的信号输入到鉴频鉴相器，进行多普勒等参数的计算，可以改善超紧耦合系统中的环路跟踪误差与更新时间的相关性。

该超紧耦合系统的示意图如图1所示，SINS信号经过惯导传感器获得运载体自身的运动特性信息(陀螺偏移、加速度计等)，将信号输入到SINS自带的定位算法中，确定运载体的位置和速度等导航数据；得到的信号一方面用于确定最终的导航参数，一方面用于输入到综合算法模块，与BDS提供的运动信息共同输入到卡尔曼滤波器，经过非线性的运算和最优估计之后得到针对SINS信号的速度和位置修正，将之前得到的导航数据修正后，输出最终的导航参数；而在BDS的主要信息环路中，接收机获得卫星定位信息以及卫星的自身运转情况，将卫星星历和解算得到的伪距、多普勒、速度位置等运载体信息和卫星信息输入到卡尔曼滤波器的综合计算单元，用于输出导航信息的误差估计，完成最终导航信号的输出；在神经网络回路中，输入到卡尔曼滤波器的信号也同样输入到神经网络之中，用于完成对于BDS跟踪环性能的调整，在BDS和SINS都正常工作的时候，对神经网络进行参数的迭代训练，训练结束后用于在BDS卫星信号进行跟踪性能的调整，输出为伪距与多普勒频移的估计值，作为参考信号输入到BDS的接收机中用于数据解算；此外，基于冗余测量的卡尔曼滤波器对于SINS和BDS的接收传感器(测量系统)同时输出辅助信号，具体信号内容为测量系统的噪声方差，用于对于收集到的信号进行预先的噪声处理，提高数据正确性。

上述过程中，定义随机线性定常离散系统，基本参考方程为：

$\left\{\begin{array}{c}{X}_k={\mathrm{\Φ}}_{k,k-1}{X}_{k-1}+{\mathrm{\Γ}}_{k,k-1}{W}_{k-1}\\ Z_k=H_k{X}_k+V_k\end{array}\right.$

其中第一式为状态方程，定义了运载体自身物理运动状态的变化方程；第二式为测量方程，用于定义测量结果与系统实际运动状态之间的转化过程。

上面两式中X_k包含的状态分量有：

$X_k={\left[\begin{array}{ccccccccccccccc}{\mathrm{\φ}}_E& {\mathrm{\φ}}_N& {\mathrm{\φ}}_U& {\mathrm{\δv}}_E& {\mathrm{\δv}}_N& {\mathrm{\δv}}_U& \mathrm{\δ}L& \mathrm{\δ}\mathrm{\λ}& \mathrm{\δ}h& {\mathrm{\ϵ}}_{bx}& {\mathrm{\ϵ}}_{by}& {\mathrm{\ϵ}}_{bz}& {\▿}_x& {\▿}_y& {\▿}_z\end{array}\right]}^T$

其中，下标E、N、U分别代表东、北、天三个方向；下标x、y、z分别代表待导航机体的右、前和上轴向；φ为惯导平台误差角；δv_E等表示速度误差；δL、δλ和δh分别表示三个方向的位置误差；ε_bx等表示陀螺漂移；等表示加速度计漂移。以上信息均可通过IMU模块获得，用于导航信息的解算。

BDS卫星接收机的I-Q信号首先进行时滞滤波操作，时滞滤波为频域操作，原有二阶PLL的传递函数模型为：

$H\left(s\right)=\frac{2{\mathrm{\ξ\ω}}_{n}s+{\mathrm{\ω}}_n^2+\frac{{\mathrm{\α}}_{IMU}{s}^2}{s+{\mathrm{\α}}_{IMU}}}{s^2+2{\mathrm{\ξ\ω}}_{n}s+{\mathrm{\ω}}_n^2}$

上式已经考虑IMU模块针对BDS信号进行的辅助，具体表现为在传递函数中添加IMU辅助信号带宽α_IMU。由于辅助带宽通常较大，此时的数值会极大趋向于1，此时H(s)可以保证输出信号较好的体现输入相位，此时相位误差基本为0。

通常在实际系统中采用双脉冲零振时滞滤波器，滤波器函数为：

$F\left(s\right)=\underset{i=1}{\overset{2}{\Σ}}{A}_i{e}^{-t_{i}s}$

式中i的取值为1或2，分别对应两个脉冲，用于实现时滞效果，而最终的环路传递函数为：

H_F(s)＝F(s)·H(s)

上式描述了接收机内部的环路基本特征；

通过SINS的计算结果和BDS的结算结果输入到卡尔曼滤波器后，除了基本的非线性运算输出导航误差修正信息之外，还进行了针对惯导和北斗接收机传感器模块的误差修正计算：

通过针对系统冗余值的一阶二阶差分分析，获得统计特性，进而自适应调节方差；

通过计算二阶差分序列的自相关可得：

E[(ΔZ₁-ΔZ₂)(ΔZ₁-ΔZ₂)^T]＝2(R₁+R₂)

计算两个一阶差分序列自相关，并相减：

$E\left({\mathrm{\ΔZ}}_1{\mathrm{\ΔZ}}_1^T\right)-E\left({\mathrm{\ΔZ}}_2{\mathrm{\ΔZ}}_2^T\right)=2(R_1-R_2)$

于是得到噪声方差矩阵：

$\left\{\begin{array}{c}{R}_1\left(k\right)=\frac{E\[({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2){({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2)}^T\]+\left(E\right({\mathrm{\ΔZ}}_1{\mathrm{\ΔZ}}_1^T)-E({\mathrm{\ΔZ}}_2{\mathrm{\ΔZ}}_2^T\left)\right)}{4}\\ R_2\left(k\right)=\frac{E\[({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2){({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2)}^T\]-\left(E\right({\mathrm{\ΔZ}}_1{\mathrm{\ΔZ}}_1^T)-E({\mathrm{\ΔZ}}_2{\mathrm{\ΔZ}}_2^T\left)\right)}{4}\end{array}\right.$

将对应的SINS/BDS系统差分序列代入上式：

$\left\{\begin{array}{c}{R}_{SINS}\left(k\right)=\frac{E{\left({\mathrm{\ΔZ}}_{SB}{{\mathrm{\ΔZ}}_{SB}}^T\right)}_{k-M:k}+(E{\left({\mathrm{\ΔZ}}_{SINS}{\mathrm{\ΔZ}}_{SINS}^T\right)}_{k-M:k}-E{\left({\mathrm{\ΔZ}}_{BDS}{\mathrm{\ΔZ}}_{BDS}^T\right)}_{k-M:k})}{4}\\ R_{BDS}\left(k\right)=\frac{E{\left({\mathrm{\ΔZ}}_{SB}{{\mathrm{\ΔZ}}_{SB}}^T\right)}_{k-M:k}-(E{\left({\mathrm{\ΔZ}}_{SINS}{\mathrm{\ΔZ}}_{SINS}^T\right)}_{k-M:k}-E{\left({\mathrm{\ΔZ}}_{BDS}{\mathrm{\ΔZ}}_{BDS}^T\right)}_{k-M:k})}{4}\end{array}\right.$

将上面的两个噪声方差反馈给测量单元即完成了针对测量系统的误差校正操作。

于是，便实现了结合冗余测量、神经网络和时滞滤波的超紧耦合导航算法，实现了针对使用SINS的运载体的高精度组合导航。

本发明的方案与现有方案相比，主要优点在于：除了超紧耦合自身带来的BDS性能辅助和SINS功能校正外，时滞滤波引入明显改善超紧耦合系统之中跟踪环路的时间常数，提高跟踪环路的性能，降低环路跟踪误差与更新时间的相关性；神经网络的引入一方面免除了在卡尔曼滤波器中进行BDS辅助信号计算的非线性矩阵维度扩大，优化了计算复杂度，同时神经网络的隐层参数适应性保证了可以动态高效地对BDS跟踪性能提供辅助；采用基于冗余测量的卡尔曼滤波器，实现了针对测量系统的误差校正，使得初始测量数据的正确性得到提高；为使用捷联惯导的运载机体提供了一种新的高精度组合导航技术途径。

附图说明

图1为基于捷联惯导和北斗的超紧耦合导航系统结构示意图；

图2为添加时滞滤波模块后的BDS接收机原理框图；

图3为基于冗余测量的自适应卡尔曼滤波器的误差计算流程；

图4为RBF神经网络的基本结构；

图5为PVT解算并添加误差的运载体运动轨迹图(BDS)；

图6为添加误差的运载体运动轨迹图(SINS)；

图7为最终得到的位置噪声方差曲线(东经方向)。

具体实施方案

本发明的具体实施方案如图1所示，具体实施步骤如下：

(1)根据运动物体状态方程及测量方程建立随机线性定常离散系统模型：

$\left\{\begin{array}{c}{X}_k={\mathrm{\Φ}}_{k,k-1}{X}_{k-1}+{\mathrm{\Γ}}_{k,k-1}{W}_{k-1}\\ Z_k=H_k{X}_k+V_k\end{array}\right.$

状态方程式中，X_k为SINS/BDS系统的n维状态变量；Φ_k，k-1为系统的nⅹn维状态转移矩阵；Г_k，k-1为nⅹp维系统过程噪声输入矩阵；W_k-1为p维系统随机过程噪声序列；

测量方程式中，Z_k为测量系统的m维观测序列，即测量值；H_k为mⅹn维观测矩阵；V_k为m维系统随机观测噪声序列；

上面两式中X_k包含的状态分量有：

$X_k={\left[\begin{array}{ccccccccccccccc}{\mathrm{\φ}}_E& {\mathrm{\φ}}_N& {\mathrm{\φ}}_U& {\mathrm{\δv}}_E& {\mathrm{\δv}}_N& {\mathrm{\δv}}_U& \mathrm{\δ}L& \mathrm{\δ}\mathrm{\λ}& \mathrm{\δ}h& {\mathrm{\ϵ}}_{bx}& {\mathrm{\ϵ}}_{by}& {\mathrm{\ϵ}}_{bz}& {\▿}_x& {\▿}_y& {\▿}_z\end{array}\right]}^T$

其中，下标E、N、U分别代表东、北、天三个方向；下标x、y、z分别代表待导航机体的右、前和上轴向；φ为惯导平台误差角；δv_E等表示速度误差；δL、δλ和δh分别表示三个方向的位置误差；ε_bx等表示陀螺漂移；等表示加速度计漂移；

在自适应卡尔曼滤波的关键步骤就是针对测量噪声的准确估计，即需要对V_k的性质和参数进行准确的分析。

(2)基于冗余测量的卡尔曼滤波算法，实现噪声方差的自适应估计与调节：

在系统量测噪声未知的条件下，通过针对系统冗余值的一阶二阶差分分析，获得统计特性，进而自适应调节方差；

通过计算二阶差分序列的自相关可得：

E[(ΔZ₁-ΔZ₂)(ΔZ₁-ΔZ₂)^T]＝2(R₁+R₂)

计算两个一阶差分序列自相关，并相减：

$E\left({\mathrm{\ΔZ}}_1{\mathrm{\ΔZ}}_1^T\right)-E\left({\mathrm{\ΔZ}}_2{\mathrm{\ΔZ}}_2^T\right)=2(R_1-R_2)$

于是得到噪声方差矩阵：

$\left\{\begin{array}{c}{R}_1\left(k\right)=\frac{E\[({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2){({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2)}^T\]+\left(E\right({\mathrm{\ΔZ}}_1{\mathrm{\ΔZ}}_1^T)-E({\mathrm{\ΔZ}}_2{\mathrm{\ΔZ}}_2^T\left)\right)}{4}\\ R_2\left(k\right)=\frac{E\[({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2){({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2)}^T\]-\left(E\right({\mathrm{\ΔZ}}_1{\mathrm{\ΔZ}}_1^T)-E({\mathrm{\ΔZ}}_2{\mathrm{\ΔZ}}_2^T\left)\right)}{4}\end{array}\right.$

将对应的SINS/BDS系统差分序列代入上式：

$\left\{\begin{array}{c}{R}_{SINS}\left(k\right)=\frac{E{\left({\mathrm{\ΔZ}}_{SB}{{\mathrm{\ΔZ}}_{SB}}^T\right)}_{k-M:k}+(E{\left({\mathrm{\ΔZ}}_{SINS}{\mathrm{\ΔZ}}_{SINS}^T\right)}_{k-M:k}-E{\left({\mathrm{\ΔZ}}_{BDS}{\mathrm{\ΔZ}}_{BDS}^T\right)}_{k-M:k})}{4}\\ R_{BDS}\left(k\right)=\frac{E{\left({\mathrm{\ΔZ}}_{SB}{{\mathrm{\ΔZ}}_{SB}}^T\right)}_{k-M:k}-(E{\left({\mathrm{\ΔZ}}_{SINS}{\mathrm{\ΔZ}}_{SINS}^T\right)}_{k-M:k}-E{\left({\mathrm{\ΔZ}}_{BDS}{\mathrm{\ΔZ}}_{BDS}^T\right)}_{k-M:k})}{4}\end{array}\right.$

其中：

ΔZ_SB＝ΔZ_SINS-ΔZ_BDS

式中，SZ_SINS为SINS的测量系统的一阶差分序列；ΔZ_BDS为BDS测量系统的；ΔZ_BDS为SINS和BDS的二阶差分序列；

在实际工程中采用划窗方式保证计算结果的实时性，上式中，k表示当前测量时刻，M为划窗长度；实际选取k-M:k范围内的信息作为参考值，以此保证自适应的实时效应；

在本专利中，针对冗余测量的滤波器进行了相关的数据仿真:

采用的初始数据为卫星的三维坐标，为保证正常运行提供七颗卫星的记录数据，同时提供伪距信息用于解算运载体位置信息；

解算中主要应用到PVT解算函数，其原理主要应用到了两个数学工具：牛顿迭代法和最小二乘法；

首先，位置求解遵循如下线性方程组：

Δρ＝H·Δp

其中，ρ是已给出的卫星伪距测量值；H为假定运载体位置与卫星位置、钟差的观测矩阵；Δp为对假定运载体位置、本地钟差的校正量；

由于给出七颗卫星的数据，方程组过定，故而需要用到最小二乘法计算上式，转化为：

Δp＝(H^TH)^-1H^TΔρ

针对上式采用迭代法求解，遵循如下关系：

$x_k=x_{k-1}+\mathrm{\Δ}x=x_k+\left[\begin{array}{c}\mathrm{\Δ}x\\ \mathrm{\Δ}y\\ \mathrm{\Δ}z\end{array}\right]$

其中Δx等三个量为针对位置的迭代修正值，设置初始位置为原定(0,0,0)；采用牛顿迭代法计算，当迭代收敛到所需要精度的时候，迭代终止，即求出用户位置；

计算位置结束后，由于仿真设定相邻两个时刻之间为1秒，所以速度可以直接采用差分序列进行速度的求解；

通过在已有位置的基础上叠加随机误差，完成测量序列的仿真构建，并将得到的仿真序列用于最终的冗余测量滤波器噪声方差计算；

(3)引入神经网络预测模块，实现对于BDS的载波跟踪环和码跟踪环的校正：

最基本的RBF神经网络分为输入层输出层和隐层，隐层用于进行从输入空间到隐藏空间的变换，一般情况下隐层空间具有较高维度；输出层作用在输入层的激活模式并提供响应；

依据超紧耦合的经典模式，系统需要在输出导航变量的同时，输出对卫星接收机的辅助变量，通常为伪距和多普勒频移的估计值；本方法中采用的就是将伪距和多普勒频移的辅助变量估计通过神经网络来实现；

神经网络的工作过程分为两个过程，分别是学习过程和输出过程；其中学习过程又分为自组织学习过程和监督学习过程，分别确定隐层函数和输出层权值；

本方法的工作方式是，当BDS和SINS正常工作时，对神经网络进行训练，BDS和SINS的组合导航输出作为目标神经网络的输出，BDS本身的伪距和多普勒频移测量结果作为目标神经网络的输入；训练算法迭代修正网络参数，从而达到误差的最小均方值；当SINS性能恶化，无法针对BDS完成直接辅助的时候，将BDS的伪距和多普勒结果作为训练完成的神经网络输入，从而预测伪距和多普勒频移的误差，通过矢量相加的方式实现辅助功能；

神经网络的主要模块在于通过学习构建隐层系统参数，前馈型神经网络的结构中，迭代过程为求解最优化问题：

构造误差信号能量函数为：

$J\left(\mathrm{\ω}\right)=E\[\frac{1}{2}\underset{k}{\Σ}{e}_k^2\left(n\right)\]$

其中，为误差信号的平方；E(·)为求期望算子；

得到能量信号之后解决最优化问题：

$\underset{\mathrm{\ω}}{min}J\left(\mathrm{\ω}\right)=E\[\frac{1}{2}\underset{k}{\Σ}{e}_k^2\left(n\right)\]$

最终得到满足上式的最优化系统参数ω；

通过数值迭代算法可得：

${\mathrm{\Δ\ω}}_{kj}=-\mathrm{\η}\frac{\∂J\left(\mathrm{\ω}\right)}{\∂\mathrm{\ω}}$

其中，Δω_kj为系统参数修正值，η为学习步长；

(4)在接收机中添加时滞滤波模块，改善环路跟踪误差与更新时间的相关性：

将脉冲序列与参考命令卷积形成的整形命令作为控制信号，滤掉参考命令中引起系统振动的频率分量，消除二阶震荡系统的震荡残留；

通常接收机中采用典型的二阶PLL(锁相环)的结构，接收机环路的传递函数为：

$H\left(s\right)=\frac{2{\mathrm{\ξ\ω}}_{n}s+{\mathrm{\ω}}_n^2+\frac{{\mathrm{\α}}_{IMU}{s}^2}{s+{\mathrm{\α}}_{IMU}}}{s^2+2{\mathrm{\ξ\ω}}_{n}s+{\mathrm{\ω}}_n^2}$

其中，ξ为二阶系统的阻尼系数，ω_n为自然圆频率，α_IMU为SINS的辅助数据带宽；

此处采用双脉冲零振时滞滤波器(ZV)，当滤波器模型参数与系统本身参数相同，且不存在参数不确定性的时候，滤波器传递函数如下：

$F\left(s\right)=\underset{i=1}{\overset{2}{\Σ}}{A}_i{e}^{-t_{i}s}$

其中，脉冲幅度和时延表达式如下：

$A_1=\frac{1}{1+K}$

$A_2=\frac{K}{1+K}$

t₁＝0

$t_2=\frac{\mathrm{\π}}{{\mathrm{\ω}}_n\sqrt{1-{\mathrm{\ξ}}^2}}$

$K=\frac{-\mathrm{\ξ}\mathrm{\π}}{e_e\sqrt{1-{\mathrm{\ξ}}^2}}$

于是，最终载波环路误差的最终传递函数为：

H_F(s)＝F(s)·H(s)

对比其他实际应用的导航方式，可以看出本设计具有较为完备的数据计算体系和性能调节模块，说明采用本发明方法可以很好地实现基于惯性导航的高精度组合导航，且各模块计算及硬件实现过程较为简洁明了，工程实现性较强。

本发明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (1)

1.一种惯性导航平台和北斗卫星的高精度超紧耦合导航方法，其特征在于：根据运动物体状态方程及测量方程建立随机线性定常离散系统模型；利用对系统的冗余测量值的差分序列进行统计分析，实现基于冗余测量的卡尔曼滤波算法，准确的估计了系统量测噪声并自适应调节捷联惯导(SINS)和北斗卫星(BDS)噪声方差R；根据超紧耦合传统结构，引入神经网络预测模块，用于辅助失锁状态下的北斗卫星载波跟踪环和码跟踪环；在北斗卫星接收机中添加时滞滤波模块，显著改善了超紧耦合系统中的环路跟踪误差与更新时间的相关性，具体包括以下步骤：

(1)根据运动物体状态方程及测量方程建立随机线性定常离散系统模型：

$\left\{\begin{array}{c}{X}_k={\mathrm{\Φ}}_{k,k-1}{X}_{k-1}+{\mathrm{\Γ}}_{k,k-1}{W}_{k-1}\\ Z_k=H_k{X}_k+V_k\end{array}\right.$

状态方程式中，X_k为SINS/BDS系统的n维状态变量；Φ_k，k-1为系统的nⅹn维状态转移矩阵；Γ_k，k-1为nⅹp维系统过程噪声输入矩阵；W_k-1为p维系统随机过程噪声序列；

测量方程式中，Z_k为测量系统的m维观测序列，即测量值；H_k为mⅹn维观测矩阵；V_k为m维系统随机观测噪声序列；

上面两式中X_k包含的状态分量有：

$X_k={\left[\begin{array}{ccccccccccccccc}{\mathrm{\φ}}_E& {\mathrm{\φ}}_N& {\mathrm{\φ}}_U& {\mathrm{\δv}}_E& {\mathrm{\δv}}_N& {\mathrm{\δv}}_U& \mathrm{\δ}L& \mathrm{\δ}\mathrm{\λ}& \mathrm{\δ}h& {\mathrm{\ϵ}}_{bx}& {\mathrm{\ϵ}}_{by}& {\mathrm{\ϵ}}_{bz}& {\▿}_x& {\▿}_y& {\▿}_z\end{array}\right]}^T$

其中，下标E、N、U分别代表东、北、天三个方向；下标x、y、z分别代表待导航机体的右、前和上轴向；φ为惯导平台误差角；δv_E等表示速度误差；δL、δλ和δh分别表示三个方向的位置误差；ε_bx等表示陀螺漂移；等表示加速度计漂移；

在自适应卡尔曼滤波的关键步骤就是针对测量噪声的准确估计，即需要对V_k的性质和参数进行准确的分析；

(2)基于冗余测量的卡尔曼滤波算法，实现噪声方差的自适应估计与调节：

在系统量测噪声未知的条件下，通过针对系统冗余值的一阶二阶差分分析，获得统计特性，进而自适应调节方差；

通过计算二阶差分序列的自相关可得：

E[(ΔZ₁-ΔZ₂)(ΔZ₁-ΔZ₂)^T]＝2(R₁+R₂)

计算两个一阶差分序列自相关，并相减：

$E\left({\mathrm{\ΔZ}}_1{\mathrm{\ΔZ}}_1^T\right)-E\left({\mathrm{\ΔZ}}_2{\mathrm{\ΔZ}}_2^T\right)=2(R_1-R_2)$

于是得到噪声方差矩阵：

$\left\{\begin{array}{c}{R}_1\left(k\right)=\frac{E\[({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2){({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2)}^T\]+\left(E\right({\mathrm{\ΔZ}}_1{\mathrm{\ΔZ}}_1^T)-E({\mathrm{\ΔZ}}_2{\mathrm{\ΔZ}}_2^T\left)\right)}{4}\\ R_2\left(k\right)=\frac{E\[({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2){({\mathrm{\ΔZ}}_1-{\mathrm{\ΔZ}}_2)}^T\]-\left(E\right({\mathrm{\ΔZ}}_1{\mathrm{\ΔZ}}_1^T)-E({\mathrm{\ΔZ}}_2{\mathrm{\ΔZ}}_2^T\left)\right)}{4}\end{array}\right.$

将对应的SINS/BDS系统差分序列代入上式：

$\left\{\begin{array}{c}{R}_{SINS}\left(k\right)=\frac{E{\left({\mathrm{\ΔZ}}_{SB}{{\mathrm{\ΔZ}}_{SB}}^T\right)}_{k-M:k}+(E{\left({\mathrm{\ΔZ}}_{SINS}{\mathrm{\ΔZ}}_{SINS}^T\right)}_{k-M:k}-E{\left({\mathrm{\ΔZ}}_{BDS}{\mathrm{\ΔZ}}_{BDS}^T\right)}_{k-M:k})}{4}\\ R_{BDS}\left(k\right)=\frac{E{\left({\mathrm{\ΔZ}}_{SB}{{\mathrm{\ΔZ}}_{SB}}^T\right)}_{k-M:k}-(E{\left({\mathrm{\ΔZ}}_{SINS}{\mathrm{\ΔZ}}_{SINS}^T\right)}_{k-M:k}-E{\left({\mathrm{\ΔZ}}_{BDS}{\mathrm{\ΔZ}}_{BDS}^T\right)}_{k-M:k})}{4}\end{array}\right.$

其中：

ΔZ_SB＝ΔZ_SINS-ΔZ_BDS

式中，ΔZ_SINS为SINS的测量系统的一阶差分序列；ΔZ_BDS为BDS测量系统的；ΔZ_BDS为SINS和BDS的二阶差分序列；

在实际工程中采用划窗方式保证计算结果的实时性，上式中，k表示当前测量时刻，M为划窗长度；实际选取k-M:k范围内的信息作为参考值，以此保证自适应的实时效应；

(3)引入神经网络预测模块，实现对于BDS的载波跟踪环和码跟踪环的校正：

最基本的RBF神经网络分为输入层输出层和隐层，隐层用于进行从输入空间到隐藏空间的变换，一般情况下隐层空间具有较高维度；输出层作用在输入层的激活模式并提供响应；

依据超紧耦合的经典模式，系统需要在输出导航变量的同时，输出对卫星接收机的辅助变量，通常为伪距和多普勒频移的估计值；本方法中采用的就是将伪距和多普勒频移的辅助变量估计通过神经网络来实现；

神经网络的工作过程分为两个过程，分别是学习过程和输出过程；其中学习过程又分为自组织学习过程和监督学习过程，分别确定隐层函数和输出层权值；

本方法的工作方式是，当BDS和SINS正常工作时，对神经网络进行训练，BDS和SINS的组合导航输出作为目标神经网络的输出，BDS本身的伪距和多普勒频移测量结果作为目标神经网络的输入；训练算法迭代修正网络参数，从而达到误差的最小均方值；当SINS性能恶化，无法针对BDS完成直接辅助的时候，将BDS的伪距和多普勒结果作为训练完成的神经网络输入，从而预测伪距和多普勒频移的误差，通过矢量相加的方式实现辅助功能；

神经网络的主要模块在于通过学习构建隐层系统参数，前馈型神经网络的结构中，迭代过程为求解最优化问题：

构造误差信号能量函数为：

$J\left(\mathrm{\ω}\right)=E\[\frac{1}{2}\underset{k}{\Σ}{e}_k^2\left(n\right)\]$

其中，为误差信号的平方；e(·)为求期望算子；

得到能量信号之后解决最优化问题：

$\underset{\mathrm{\ω}}{min}J\left(\mathrm{\ω}\right)=E\[\frac{1}{2}\underset{k}{\Σ}{e}_k^2\left(n\right)\]$

最终得到满足上式的最优化系统参数ω；

通过数值迭代算法可得：

${\mathrm{\Δ\ω}}_{kj}=-\mathrm{\η}\frac{\∂J\left(\mathrm{\ω}\right)}{\∂\mathrm{\ω}}$

其中，Δω_kj为系统参数修正值，η为学习步长；

(4)在接收机中添加时滞滤波模块，改善环路跟踪误差与更新时间的相关性：

将脉冲序列与参考命令卷积形成的整形命令作为控制信号，滤掉参考命令中引起系统振动的频率分量，消除二阶震荡系统的震荡残留；

通常接收机中采用典型的二阶PLL(锁相环)的结构，接收机环路的传递函数为：

$H\left(s\right)=\frac{2{\mathrm{\ξ\ω}}_{n}s+{\mathrm{\ω}}_n^2+\frac{{\mathrm{\α}}_{IMU}{s}^2}{s+{\mathrm{\α}}_{IMU}}}{s^2+2{\mathrm{\ξ\ω}}_{n}s+{\mathrm{\ω}}_n^2}$

其中，ξ为二阶系统的阻尼系数，ω_n为自然圆频率，α_IMU为SINS的辅助数据带宽；

此处采用双脉冲零振时滞滤波器(ZV)，当滤波器模型参数与系统本身参数相同，且不存在参数不确定性的时候，滤波器传递函数如下：

$F\left(s\right)=\underset{i=1}{\overset{2}{\Σ}}{A}_i{e}^{-t_{i}s}$

其中，脉冲幅度和时延表达式如下：

$A_1=\frac{1}{1+K}$

$A_2=\frac{K}{1+K}$

t₁＝0

$t_2=\frac{\mathrm{\π}}{{\mathrm{\ω}}_n\sqrt{1-{\mathrm{\ξ}}^2}}$

$K=\frac{-\mathrm{\ξ}\mathrm{\π}}{e_e\sqrt{1-{\mathrm{\ξ}}^2}}$

于是，最终载波环路误差的最终传递函数为：

H_F(s)＝F(s)·H(s)。