CN104344836A - 一种基于姿态观测的冗余惯导系统光纤陀螺系统级标定方法 - Google Patents

Abstract

一种基于姿态观测的冗余惯导系统光纤陀螺系统级标定方法，有八大步骤：一、将惯导安装转台上，确定载体的初始位置参数；二、确定光纤陀螺轴向与惯导本体坐标系安装关系；三、惯导预热，采集光纤陀螺输出数据进行精标定；四、使惯导位于固定位置东、北、天位置，采集光纤陀螺输出数据，标定光纤陀螺零偏误差并实时补偿；五、在东、北、天坐标系下绕X轴以角速度ω旋转一周，进行第一次标定参数修正；六、在东、北、天坐标系下绕Y轴以角速度ω旋转一周，进行第二次标定参数修正；七、在东、北、天坐标系下绕Z轴以角速度ω旋转一周，进行第三次标定参数修正；八、通过第三次标定参数修正，得到斜置光纤陀螺高精度的标定参数零偏、标度因数、失准角结果。

Description

一种基于姿态观测的冗余惯导系统光纤陀螺系统级标定方法

技术领域

本发明属于惯性导航技术领域，特别是涉及一种基于姿态观测的冗余惯导系统光纤陀螺系统级标定方法。

背景技术

惯性技术是一项涉及多学科的综合技术，它是惯性导航和惯性制导技术、惯性仪表技术、惯性测量技术以及有关系统和装置技术的统称。惯性导航系统依靠自身的惯性敏感元件，不依赖任何外界信息测量导航参数，因此它不受天然的或人为的干扰，具有很好的隐蔽性，是一种完全自主式的导航系统。

惯性导航系统标定是通过比较系统中惯性器件的输出和已知参考输入，确定一组参数使惯导系统输出与输入相吻合的过程，惯导系统标定的理论基础是系统辨识和参数估计，其目的是确定惯性器件组合的数学误差模型或误差数学的模型参数。惯导系统使用之前必须进行标定，对器件零偏、标度因数、安装失准角等参数进行补偿。冗余惯导系统中器件安装方式与三轴惯导系统有较大差异，传统的标定方法在冗余惯导系统中实现起来特别繁琐，而且精度较低。因此，新的适用于冗余系统的标定方法已成为必然需求。

惯导系统常用标定方法主要有：分立式标定法、模观测标定方法、系统级标定方法等。

分立式标定方法也称为基于转台标定方法，需要转台为系统提供标准输入信息，对转台精度要求较高，同时，分立式标定过程依赖转台，一般只能在实验室进行。

模观测标定方法是指基于惯导系统输入加速度、角速度激励的模分别和加速度计比力测量、陀螺角速度测量的模相等的原理，以输入加速度、角速度的模作为观测，计算惯导系统参数的方法。目前模观测标定计算采用迭代算法，其收敛性严重依赖标定参数初值。

系统级标定方法主要基于导航解算误差的原理：惯导系统进入导航状态之后，其参数误差(惯性器件参数误差、初始对准姿态误差，初始位置误差等)经由导航解算会传递到导航结果(位置、速度、姿态等)中去，表现为导航误差，如能获取导航误差的全部或部分信息，就可能对惯导系统参数做出估计。系统级标定方法降低了对转台的精度要求，利用低精度转台就可以达到较高的标定精度，因此是现场标定的理想方法。

系统级标定方法相对于其他标定方法拥有较大优势，在现场标定和高精度标定的场合，系统级标定将占据重要地位。在冗余系统或者特定场合(如空间用单表陀螺等)中，光纤陀螺并不一定按照笛卡尔坐标系正交安装，而是采用特定的斜置安装方式以满足特定需求、提高系统可靠性和精度，传统意义上标定方法相对繁琐、精度低甚至不再适用，因此，斜置光纤陀螺高精度系统级标定方法拥有重大需求。

发明内容

本发明的目的在于提高冗余型光纤捷联惯性导航系统斜置光纤陀螺初始标定精度，提供了一种基于姿态观测的冗余惯导系统光纤陀螺系统级标定方法，它是一种适用于冗余光纤惯性导航系统斜置光纤陀螺系统级精标定方法。

本发明一种基于姿态观测的冗余惯导系统光纤陀螺系统级标定方法，该方法体步骤如下：

步骤一：将惯导系统安装在转台上，确定载体的初始位置参数，包括经度、纬度等；

步骤二：确定光纤陀螺轴向与惯导系统本体坐标系安装关系即安装角，计算安装矩阵；

步骤三：惯导系统预热，在已有光纤陀螺粗略标定参数(零偏、标度因数、失准角等)基础上(粗标定完成)，准备采集光纤陀螺输出数据进行精标定；

光纤陀螺输出的数据为载体相对于惯性参考系的角速度

步骤四：使惯导系统位于某一固定位置(如东、北、天位置)，采集光纤陀螺输出数据，以转台姿态角与惯导系统姿态角之间误差作为观测量，标定光纤陀螺零偏误差并实时补偿；

步骤五：在东、北、天坐标系下绕X轴以角速度ω(角速度ω应远大于地球转速，忽略地球转速影响)旋转一周，以转台姿态角与惯导系统姿态角之间的误差作为观测量，进行第一次标定参数修正；

步骤六：在东、北、天坐标系下绕Y轴以角速度ω(角速度ω应远大于地球转速，忽略地球转速影响)旋转一周，以转台姿态角与惯导系统姿态角之间的误差作为观测量，进行第二次标定参数修正；

步骤七：在东、北、天坐标系下绕Z轴以角速度ω(角速度ω应远大于地球转速，忽略地球转速影响)旋转一周，以转台姿态角与惯导系统姿态角之间的误差作为观测量，进行第三次标定参数修正；

步骤八：通过步骤七的第三次标定参数修正，得到斜置光纤陀螺高精度的标定参数(零偏、标度因数、失准角)结果。

其中，步骤四至步骤七中采用基于卡尔曼滤波技术的误差标定方法，利用转台姿态角与系统姿态角之间的误差作为观测量，通过卡尔曼滤波迭代，估计光纤陀螺零偏误差、标度因数误差及失准角，对光纤陀螺粗标定结果进行修正。具体步骤如下：

步骤一：建立光纤陀螺标定的系统状态方程和观测方程。

若以冗余系统中所有光纤陀螺标定参数构建卡尔曼滤波器，最终只能得到标定参数最小二乘解，而非真实解。此处采用任意三个编号为a,b,c的光纤陀螺和加速度计构建一套子惯导系统。

子惯导系统中光纤陀螺误差项系统级精标定卡尔曼滤波器的状态方程为：

${\stackrel{\·}{X}}_g=A_g{X}_g+W_g---\left(1\right)$

其中15维状态矢量包括俯仰、横滚及航向姿态角误差ψ^T；光纤陀螺零偏残差矢量：光纤陀螺标度因数误差残差矢量：光纤陀螺安装失准角残差矢量：

$\mathrm{\Δ}{\stackrel{\‾}{\mathrm{\α}}}_g={[{{\mathrm{\δ\α}}^{\′}}_{\mathrm{ga}},{{\mathrm{\δ\α}}^{\′}}_{\mathrm{gb}},{{\mathrm{\δ\α}}^{\′}}_{\mathrm{gc}}]}^T,\mathrm{\Δ}{\stackrel{\‾}{\mathrm{\β}}}_g={[{{\mathrm{\δ\β}}^{\′}}_{\mathrm{ga}},{{\mathrm{\δ\β}}^{\′}}_{\mathrm{gb}},{{\mathrm{\δ\β}}^{\′}}_{\mathrm{gc}}]}^T.$ 表示系统状态的微分，状态方程中15阶方阵(状态矩阵)可表示为如下形式：

$A_g={\left[\begin{array}{c}\begin{array}{cc}{A}_{g1}& A_{g2}\end{array}\\ 0_{12\×15}\end{array}\right]}_{15\×15}---\left(2\right)$

其中A_g1和A_g2可表示为如下形式：

$A_{g1}=\left[\begin{array}{ccc}0& {\mathrm{\ω}}_{\mathrm{iez}}^n& {-\mathrm{\ω}}_{\mathrm{iey}}^n\\ {-\mathrm{\ω}}_{\mathrm{iez}}^n& 0& {\mathrm{\ω}}_{\mathrm{iex}}^n\\ {\mathrm{\ω}}_{\mathrm{iey}}^n& {-\mathrm{\ω}}_{\mathrm{iex}}^n& 0\end{array}\right]---\left(3\right)$

$A_{g2}=\left[\begin{array}{cccc}{\stackrel{\‾}{C}}_b^n\·{\stackrel{\‾}{H}}^{-1}& {\stackrel{\‾}{C}}_b^n\·{\stackrel{\‾}{H}}^{-1}\·{\stackrel{\‾}{H}}_{\mathrm{\ω}}& {\stackrel{\‾}{C}}_b^n\·{\stackrel{\‾}{H}}^{-1}\·{\stackrel{\‾}{P}}_{\mathrm{\ω}}& {\stackrel{\‾}{C}}_b^n\·{\stackrel{\‾}{H}}^{-1}\·{\stackrel{\‾}{Q}}_{\mathrm{\ω}}\end{array}\right]---\left(4\right)$

表示地球自转角速度在导航坐标系n系下的投影(下标的x,y及z表示沿导航坐标系的三个坐标轴)。为子惯导系统的配置矩阵

其中h_i＝[cos(α_i)cos(β_i)]·i+[sin(α_i)cos(β_i)]·j+[sin(β_i)]·k(i＝a,b,c)，这里h_i,i,j和k表示轴H_i,X_b,Y_b和Z_b上的单位矢量，α_i表示h_i在X_b-Y_b平面上的投影向量与轴的夹角，β_i表示h_i与X_b-Y_b平面的夹角(如图1)，对角线矩阵 表示载体系(b系)相对惯性系(i系)角速度在b系下的投影(即是陀螺输出角速度值)。 ${\stackrel{\‾}{P}}_{\mathrm{\ω}}=\mathrm{diag}\{p_a\·{\mathrm{\ω}}_{\mathrm{ib}}^b,p_b\·{\mathrm{\ω}}_{\mathrm{ib}}^b,p_c\·{\mathrm{\ω}}_{\mathrm{ib}}^b\},$

其中p_i＝[sin(α_i)cos(β_i)-cos(α_i)cos(β_i)0]。

其中q_i＝[-cos(α_i)sin(β_i)-sin(α_i)sin(β_i)cos(β_i)](如图1、图2所示)。为转台所示捷联姿态矩阵。

假设Θ_g为系统噪声方差阵，式(1)中W_g为服从正态分布N(0,Θ_g)的系统噪声，满足如下条件：

$\left\{\begin{array}{c}E\left[W_g\left(i\right)\right]=0\\ E\left[W_g\left(i\right)W_g{\left(i\right)}^T\right]={\mathrm{\Θ}}_g\end{array}\right.---\left(5\right)$

以姿态误差作为观测量构建卡尔曼滤波器的量测方程，形式如下：

Z_g＝F_gX_g+         (6)

上式中状态矢量X_g的定义与式(1)相同，观测量其中θ_x,θ_y,θ_z为系统解算姿态，为转台姿态(以姿态作为观测量)。量测矩阵F_g为15阶方阵，可表示为如下形式：

$F={\left[\begin{array}{cccc}1& 0& 0& \\ 0& 1& 0& 0_{3\×12}\\ 0& 0& 1& \end{array}\right]}_{3\×15}---\left(7\right)$

假设R_g为量测噪声方差阵，式(11)中V_g为服从正态分布N(0,R_g)的系统噪声，满足如下条件：

$\left\{\begin{array}{c}E\left[V_g\left(i\right)\right]=0\\ E\left[V_g\left(i\right)V_g{\left(i\right)}^T\right]=R_g\end{array}\right.---\left(8\right)$

步骤二：对系统状态方程进行离散化。

对步骤一建立的系统状态变量进行估计，需要对系统状态方程进行离散化。离散化采用泰勒级数展开，则：

$\mathrm{\Φ}(k+1,k)=I+\mathrm{TA}\left(k\right)+\frac{T^2}{2!}{A}^2\left(k\right)+\frac{T^3}{3!}{A}^3\left(k\right)+\·\·\·---\left(9\right)$

其中：Φ(k+1,k)为状态一步转移矩阵、I为十五阶单位阵、A(k)为状态转移矩阵，T为滤波周期。

系统模型噪声的方差为：

$Q\left(k\right)=\mathrm{QT}+[\mathrm{AQ}+{\left(\mathrm{AQ}\right)}^T]\frac{T^2}{2!}+\left\{A\right[\mathrm{AQ}+{\left(\mathrm{AQ}\right)}^T]+A{[\mathrm{AQ}+{\mathrm{QA}}^T]}^T\}\frac{T^3}{3!}+\·\·\·---\left(10\right)$

其中：Q(k)为离散系统噪声方差阵、Q连续系统噪声方程强度阵、A为状态转移矩阵。

步骤三：进行卡尔曼滤波状态估计。

对卡尔曼滤波器进行迭代，状态预测估计方程、方差预测方程、状态预测估计方程、方差迭代方程以及滤波增益方程可表示为如下形式：

$P_{k/k-1}={\mathrm{\Φ}}_{k,k-1}{P}_{k-1}{\mathrm{\Φ}}_{k,k-1}^T+{\mathrm{\Γ}}_{k-1}{Q}_{k-1}{\mathrm{\Γ}}_{k-1}^T---\left(12\right)$

${\hat{X}}_k={\hat{X}}_{k/k-1}+K_k(Z_k-H_k{\hat{X}}_{k/k-1})---\left(13\right)$

P_k＝(I-K_kH_k)P_k/k-1     (14)

$K_k=P_{k/k-1}{H}_k^T{(H_k{P}_{k/k-1}{H}_k^T+R_k)}^{-1}---\left(15\right)$

最后估计得到冗余系统中斜置光纤陀螺零偏误差、标度因数误差以及失准角的标定结果，并对粗标定结果进行修正。

本发明的优点在于：

本发明所述方法在分立式标定(简称粗标定，粗标定之后光纤陀螺标定参数尚有残差)基础上，不增加传统标定程式，通过转动序列、滤波器设计的参数等的合理设计，对光纤陀螺零偏误差、标度因数误差及失准角进行进一步补偿，能够很大程度上提高斜置光纤陀螺标定精度。

附图说明

图1为理想传感器轴向与系统本体坐标系的安装关系示意图。

图2为实际传感器轴向与系统本体坐标系的安装关系示意图。

图3为系统级标定仿真平台示意图。

图4为六冗余RFINS的正十二面体安装方式示意图。

图5为系统级标定仿真过程中的子惯导系统示意图。

图6(a)为陀螺仪零偏误差随时间变化曲线示意图。

图6(b)为陀螺仪标度因数误差随时间变化曲线示意图。

图6(c)为陀螺仪失准角误差A随时间变化曲线示意图。

图6(d)为陀螺仪失准角误差B随时间变化曲线示意图。

图7为本发明流程框图。

图中符号说明如下：

OX_bY_bZ_b为系统本体坐标系；

H_i为理想传感器轴向，α_i为H_i在X_bOY_b平面投影与X_b轴的夹角，β为H_i与Z_b轴的夹角；

H_i'为实际传感器轴向，δα_i为H_i'与H_i夹角在X_bOY_b平面的投影，δβ_i为H_i'、H_i与Z_b轴夹角之差；

M1、M2、M3分别为X_bOY_b平面、X_bOZ_b平面、Y_bOZ_b平面，α角为各平面内(M1、M2、M3)传感器轴向与最近系统本体坐标轴夹角，ABCDEF分别为六个传感器轴向；

ω表示系统绕坐标轴旋转角速度(角速度ω应远大于地球转速)。

具体实施方式

见图7，本发明一种基于姿态观测的冗余惯导系统光纤陀螺系统级标定方法，包括以下几个步骤：

步骤一：将光纤捷联惯组安装在转台上，确定载体的初始位置参数，包括经度、纬度等；

步骤二：确定光纤陀螺轴向与系统本体坐标系安装关系即安装角，计算安装矩阵；

步骤三：光纤捷联惯组预热，在已有光纤陀螺粗略标定参数(零偏、标度因数、失准角等)基础上(粗标定完成)，准备采集光纤陀螺输出数据进行精标定；

光纤陀螺输出的数据为载体相对于惯性参考系的角速度

步骤四：使系统位于东、北、天位置，采集光纤陀螺输出数据，以转台姿态角与系统姿态角之间误差作为观测量，标定光纤陀螺零偏误差并实时补偿；

步骤五：在东、北、天坐标系下绕X轴以角速度ω(角速度ω应远大于地球转速，忽略地球转速影响)旋转一周，以转台姿态角与系统姿态角之间的误差作为观测量，进行第一次标定参数修正；

步骤六：在东、北、天坐标系下绕Y轴以角速度ω(角速度ω应远大于地球转速，忽略地球转速影响)旋转一周，以转台姿态角与系统姿态角之间的误差作为观测量，进行第二次标定参数修正；

步骤七：在东、北、天坐标系下绕Z轴以角速度ω(角速度ω应远大于地球转速，忽略地球转速影响)旋转一周，以转台姿态角与系统姿态角之间的误差作为观测量，进行第三次标定参数修正；

步骤八：通过步骤七的第三次标定参数修正，得到斜置光纤陀螺高精度的标定参数(零偏、标度因数、失准角)结果。

步骤四至步骤七中采用基于卡尔曼滤波技术的误差标定方法，利用转台姿态角与系统姿态角之间的误差作为观测量，通过卡尔曼滤波迭代，估计光纤陀螺零偏误差、标度因数误差及失准角，对光纤陀螺粗标定结果进行修正。具体步骤如下：

步骤一：建立光纤陀螺标定的系统状态方程和观测方程。

若以冗余系统中所有光纤陀螺标定参数构建卡尔曼滤波器，最终只能得到标定参数最小二乘解，而非真实解。此处采用任意三个编号为a,b,c的光纤陀螺和加速度计构建一套子惯导系统。

子惯导系统中光纤陀螺误差项系统级精标定卡尔曼滤波器的状态方程为：

${\stackrel{\·}{X}}_g=A_g{X}_g+W_g---\left(1\right)$

其中15维状态矢量包括俯仰、横滚及航向姿态角误差ψ^T；光纤陀螺零偏残差矢量：ΔB_g＝[ΔB_ga,ΔB_gb,ΔB_gc]^T，光纤陀螺标度因数误差残差矢量：光纤陀螺安装失准角残差矢量：

$\mathrm{\Δ}{\stackrel{\‾}{\mathrm{\α}}}_g={[{{\mathrm{\δ\α}}^{\′}}_{\mathrm{ga}},{{\mathrm{\δ\α}}^{\′}}_{\mathrm{gb}},{{\mathrm{\δ\α}}^{\′}}_{\mathrm{gc}}]}^T,\mathrm{\Δ}{\stackrel{\‾}{\mathrm{\β}}}_g={[{{\mathrm{\δ\β}}^{\′}}_{\mathrm{ga}},{{\mathrm{\δ\β}}^{\′}}_{\mathrm{gb}},{{\mathrm{\δ\β}}^{\′}}_{\mathrm{gc}}]}^T.$ ${\stackrel{\·}{X}}_g$ 表示系统状态的微分，状态方程中15阶方阵(状态矩阵)可表示为如下形式：

$A_g={\left[\begin{array}{c}\begin{array}{cc}{A}_{g1}& A_{g2}\end{array}\\ 0_{12\×15}\end{array}\right]}_{15\×15}---\left(2\right)$

其中A_g1和A_g2可表示为如下形式：

$A_{g1}=\left[\begin{array}{ccc}0& {\mathrm{\ω}}_{\mathrm{iez}}^n& {-\mathrm{\ω}}_{\mathrm{iey}}^n\\ {-\mathrm{\ω}}_{\mathrm{iez}}^n& 0& {\mathrm{\ω}}_{\mathrm{iex}}^n\\ {\mathrm{\ω}}_{\mathrm{iey}}^n& {-\mathrm{\ω}}_{\mathrm{iex}}^n& 0\end{array}\right]---\left(3\right)$

$A_{g2}=\left[\begin{array}{cccc}{\stackrel{\‾}{C}}_b^n\·{\stackrel{\‾}{H}}^{-1}& {\stackrel{\‾}{C}}_b^n\·{\stackrel{\‾}{H}}^{-1}\·{\stackrel{\‾}{H}}_{\mathrm{\ω}}& {\stackrel{\‾}{C}}_b^n\·{\stackrel{\‾}{H}}^{-1}\·{\stackrel{\‾}{P}}_{\mathrm{\ω}}& {\stackrel{\‾}{C}}_b^n\·{\stackrel{\‾}{H}}^{-1}\·{\stackrel{\‾}{Q}}_{\mathrm{\ω}}\end{array}\right]---\left(4\right)$

表示地球自转角速度在导航坐标系n系下的投影(下标的x,y及z表示沿导航坐标系的三个坐标轴)。为子惯导系统的配置矩阵

其中h_i＝[cos(α_i)cos(β_i)]·i+[sin(α_i)cos(β_i)]·j+[sin(β_i)]·k(i＝a,b,c)，这里h_i,i,j和k表示轴H_i,X_b,Y_b和Z_b上的单位矢量，α_i表示h_i在X_b-Y_b平面上的投影向量与轴的夹角，β_i表示h_i与X_b-Y_b平面的夹角(如图1)，对角线矩阵 表示载体系(b系)相对惯性系(i系)角速度在b系下的投影(即是陀螺输出角速度值)。 ${\stackrel{\‾}{P}}_{\mathrm{\ω}}=\mathrm{diag}\{p_a\·{\mathrm{\ω}}_{\mathrm{ib}}^b,p_b\·{\mathrm{\ω}}_{\mathrm{ib}}^b,p_c\·{\mathrm{\ω}}_{\mathrm{ib}}^b\},$

其中p_i＝[sin(α_i)cos(β_i)-cos(α_i)cos(β_i)0]。

其中q_i＝[-cos(α_i)sin(β_i)-sin(α_i)sin(β_i)cos(β_i)](如图1、图2所示)。为转台所示捷联姿态矩阵。

假设Θ_g为系统噪声方差阵，式(1)中W_g为服从正态分布N(0,Θ_g)的系统噪声，满足如下条件：

$\left\{\begin{array}{c}E\left[W_g\left(i\right)\right]=0\\ E\left[W_g\left(i\right)W_g{\left(i\right)}^T\right]={\mathrm{\Θ}}_g\end{array}\right.---\left(5\right)$

以姿态误差作为观测量构建卡尔曼滤波器的量测方程，形式如下：

Z_g＝F_gX_g+        (6)

上式中状态矢量X_g的定义与式(1)相同，观测量其中θ_x,θ_y,θ_z为系统解算姿态，为转台姿态(以姿态作为观测量)。量测矩阵F_g为15阶方阵，可表示为如下形式：

$F={\left[\begin{array}{cccc}1& 0& 0& \\ 0& 1& 0& 0_{3\×12}\\ 0& 0& 1& \end{array}\right]}_{3\×15}---\left(7\right)$

假设R_g为量测噪声方差阵，式(11)中V_g为服从正态分布N(0,R_g)的系统噪声，满足如下条件：

$\left\{\begin{array}{c}E\left[V_g\left(i\right)\right]=0\\ E\left[V_g\left(i\right)V_g{\left(i\right)}^T\right]=R_g\end{array}\right.---\left(8\right)$

步骤二：对系统状态方程进行离散化。

对步骤一建立的系统状态变量进行估计，需要对系统状态方程进行离散化。离散化采用泰勒级数展开，则：

$\mathrm{\Φ}(k+1,k)=I+\mathrm{TA}\left(k\right)+\frac{T^2}{2!}{A}^2\left(k\right)+\frac{T^3}{3!}{A}^3\left(k\right)+\·\·\·---\left(9\right)$

其中：Φ(k+1,k)为状态一步转移矩阵、I为十五阶单位阵、A(k)为状态转移矩阵，T为滤波周期。

系统模型噪声的方差为：

$Q\left(k\right)=\mathrm{QT}+[\mathrm{AQ}+{\left(\mathrm{AQ}\right)}^T]\frac{T^2}{2!}+\left\{A\right[\mathrm{AQ}+{\left(\mathrm{AQ}\right)}^T]+A{[\mathrm{AQ}+{\mathrm{QA}}^T]}^T\}\frac{T^3}{3!}+\·\·\·---\left(10\right)$

其中：Q(k)为离散系统噪声方差阵、Q连续系统噪声方程强度阵、A为状态转移矩阵。

步骤三：进行卡尔曼滤波状态估计。

对卡尔曼滤波器进行迭代，状态预测估计方程、方差预测方程、状态预测估计方程、方差迭代方程以及滤波增益方程可表示为如下形式：

$P_{k/k-1}={\mathrm{\Φ}}_{k,k-1}{P}_{k-1}{\mathrm{\Φ}}_{k,k-1}^T+{\mathrm{\Γ}}_{k-1}{Q}_{k-1}{\mathrm{\Γ}}_{k-1}^T---\left(12\right)$

${\hat{X}}_k={\hat{X}}_{k/k-1}+K_k(Z_k-H_k{\hat{X}}_{k/k-1})---\left(13\right)$

P_k＝(I-K_kH_k)P_k/k-1      (14)

$K_k=P_{k/k-1}{H}_k^T{(H_k{P}_{k/k-1}{H}_k^T+R_k)}^{-1}---\left(15\right)$

最后估计得到冗余系统中斜置光纤陀螺零偏误差、标度因数误差以及失准角的标定结果，并对粗标定结果进行修正。

实施例仿真：

下面结合实施例及标定仿真对本发明做进一步说明，本发明实施例以某六冗余型光纤捷联惯导系统的标定为例。

(1)系统级标定仿真平台

系统级标定仿真平台是在分立式标定仿真平台的基础上加入光纤陀螺误差精标定模块，包括光纤陀螺误差精标定卡尔曼滤波器。标定仿真平台结构框图如图3所示。

(2)系统级标定仿真条件及结论

六冗余RFINS的正十二面体安装方式由图4所示，六冗余斜置RFINS系统结构中的ABC光纤陀螺与加速度计构成，如图5所示。其中AB轴位于面M2，C轴位于面M1，其与坐标轴的夹角为α＝31°43'2.9″。

根据图5的配置方式，可得系统的安装矩阵为：

$H={\left[\begin{array}{ccc}\sin \mathrm{\α}& -\sin & \cos \mathrm{\α}\\ 0& 0& \sin \mathrm{\α}\\ \cos \mathrm{\α}& \cos \mathrm{\α}& 0\end{array}\right]}^T---\left(16\right)$

子惯导系统中光纤陀螺的零偏、标度及安装失准角等各项标定参数如下表所示，仿真过程中采样周期为20ms。

表1冗余光纤陀螺标定参数

	光纤陀螺
		零偏	0.1°/h
标度因数(ppm)	100
		失准角A(角分)	1
失准角B(角分)	1
		噪声方差	0.005°/h

在光纤陀螺误差系统级标定过程中，利用前述精标定卡尔曼滤波器以及转台转动序列。系统每次转动过程持续5s，转动完成后静止1min。下图描述了子惯导系统中光纤陀螺的零偏误差、标度因数误差和两类安装失准角误差的估计曲线。

从图6a-d可以看出，当以系统姿态与转台姿态的差值作为观测量时(系统速度置0)，由于系统静止过程中光纤陀螺标度因数及安装失准角不被激发，因此在第一个位置可以较为准确地标定陀螺零偏误差；当系统分别绕本体坐标系转动一周后，光纤陀螺的标度因数误差和安装失准角误差可以准确标定；当系统转至第三个位置时A类安装失准角即可准确标定，是由于系统从绕天向转动时A类失准角未影响系统姿态误差，而到绕东向转动时A类失准角将作用于系统姿态误差，而此时标度因数误差和B类安装失准角误差已经到达非真值的稳态，所以A类安装失准角可以准确求得。

Claims (2)

1.一种基于姿态观测的冗余惯导系统光纤陀螺系统级标定方法，其特征在于：该方法体步骤如下：

步骤一：将惯导系统安装在转台上，确定载体的初始位置参数，包括经度、纬度；

步骤二：确定光纤陀螺轴向与惯导系统本体坐标系安装关系即安装角，计算安装矩阵；

步骤三：惯导系统预热，在已有光纤陀螺粗略标定参数零偏、标度因数、失准角基础上，准备采集光纤陀螺输出数据进行精标定；光纤陀螺输出的数据为载体相对于惯性参考系的角速度

步骤四：使惯导系统位于某一固定位置东、北、天位置，采集光纤陀螺输出数据，以转台姿态角与惯导系统姿态角之间误差作为观测量，标定光纤陀螺零偏误差并实时补偿；

步骤五：在东、北、天坐标系下绕X轴以角速度ω旋转一周，以转台姿态角与惯导系统姿态角之间的误差作为观测量，进行第一次标定参数修正；

步骤六：在东、北、天坐标系下绕Y轴以角速度ω旋转一周，以转台姿态角与惯导系统姿态角之间的误差作为观测量，进行第二次标定参数修正；

步骤七：在东、北、天坐标系下绕Z轴以角速度ω旋转一周，以转台姿态角与惯导系统姿态角之间的误差作为观测量，进行第三次标定参数修正；

步骤八：通过步骤七的第三次标定参数修正，得到斜置光纤陀螺高精度的标定参数零偏、标度因数、失准角结果。

2.根据权利要求1所述的一种基于姿态观测的冗余惯导系统光纤陀螺系统级标定方法，其特征在于：步骤四至步骤七中采用基于卡尔曼滤波技术的误差标定方法，利用转台姿态角与系统姿态角之间的误差作为观测量，通过卡尔曼滤波迭代，估计光纤陀螺零偏误差、标度因数误差及失准角，对光纤陀螺粗标定结果进行修正，具体步骤如下：

步骤一：建立光纤陀螺标定的系统状态方程和观测方程；

若以冗余系统中所有光纤陀螺标定参数构建卡尔曼滤波器，最终只能得到标定参数最小二乘解，而非真实解，此处采用任意三个编号为a,b,c的光纤陀螺和加速度计构建一套子惯导系统；

子惯导系统中光纤陀螺误差项系统级精标定卡尔曼滤波器的状态方程为：

${\stackrel{\·}{X}}_g=A_g{X}_g+W_g---\left(1\right)$

其中15维状态矢量包括俯仰、横滚及航向姿态角误差ψ^T；光纤陀螺零偏残差矢量：ΔB_g＝[ΔB_ga,ΔB_gb,ΔB_gc]^T，光纤陀螺标度因数误差残差矢量：光纤陀螺安装失准角残差矢量：

$\mathrm{\Δ}{\stackrel{\‾}{\mathrm{\α}}}_g={[{{\mathrm{\δ\α}}^{\′}}_{\mathrm{ga}},{{\mathrm{\δ\α}}^{\′}}_{\mathrm{gb}},{{\mathrm{\δ\α}}^{\′}}_{\mathrm{gc}}]}^T,\mathrm{\Δ}{\stackrel{\‾}{\mathrm{\β}}}_g={[{{\mathrm{\δ\β}}^{\′}}_{\mathrm{ga}},{{\mathrm{\δ\β}}^{\′}}_{\mathrm{gb}},{{\mathrm{\δ\β}}^{\′}}_{\mathrm{gc}}]}^T;$ 表示系统状态的微分，状态方程中15阶方阵即状态矩阵，表示为如下形式：

$A_g={\left[\begin{array}{cc}{A}_{g1}& A_{g2}\\ 0_{12\×15}& \end{array}\right]}_{15\×15}---\left(2\right)$

其中A_g1和A_g2表示为如下形式：

$A_{g1}=\left[\begin{array}{ccc}0& {\mathrm{\ω}}_{\mathrm{iez}}^n& -{\mathrm{\ω}}_{\mathrm{iey}}^n\\ -{\mathrm{\ω}}_{\mathrm{iez}}^n& 0& {\mathrm{\ω}}_{\mathrm{iex}}^n\\ {\mathrm{\ω}}_{\mathrm{iey}}^n& -{\mathrm{\ω}}_{\mathrm{iex}}^n& 0\end{array}\right]---\left(3\right)$

$A_{g2}=\left[\begin{array}{cccc}{\stackrel{\‾}{C}}_b^n\·{\stackrel{\‾}{H}}^{-1}& {\stackrel{\‾}{C}}_b^n\·{\stackrel{\‾}{H}}^{-1}\·{\stackrel{\‾}{H}}_{\mathrm{\ω}}& {\stackrel{\‾}{C}}_b^n\·{\stackrel{\‾}{H}}^{-1}\·{\stackrel{\‾}{P}}_{\mathrm{\ω}}& {\stackrel{\‾}{C}}_b^n\·{\stackrel{\‾}{H}}^{-1}\·{\stackrel{\‾}{Q}}_{\mathrm{\ω}}\end{array}\right]---\left(4\right)$

表示地球自转角速度在导航坐标系n系下的投影，下标的x,y及z表示沿导航坐标系的三个坐标轴；为子惯导系统的配置矩阵

其中h_i＝[cos(α_i)cos(β_i)]·i+[sin(α_i)cos(β_i)]·j+[sin(β_i)]·k(i＝a,b,c)，这里h_i,i,j和k表示轴H_i,X_b,Y_b和Z_b上的单位矢量，α_i表示h_i在X_b-Y_b平面上的投影向量与轴的夹角，β_i表示h_i与X_b-Y_b平面的夹角，对角线矩阵 表示载体系即b系相对惯性系即i系角速度在b系下的投影，即是陀螺输出角速度值； ${\stackrel{\‾}{P}}_{\mathrm{\ω}}=\mathrm{diag}\{p_a\·{\mathrm{\ω}}_{\mathrm{ib}}^b,p_b\·{\mathrm{\ω}}_{\mathrm{ib}}^b,p_c\·{\mathrm{\ω}}_{\mathrm{ib}}^b\},$

其中p_i＝[sin(α_i)cos(β_i) -cos(α_i)cos(β_i) 0]；

其中q_i＝[-cos(α_i)sin(β_i) -sin(α_i)sin(β_i) cos(β_i)]，为转台所示捷联姿态矩阵；

假设Θ_g为系统噪声方差阵，式(1)中W_g为服从正态分布N(0,Θ_g)的系统噪声，满足如下条件：

$\left\{\begin{array}{c}E\left[W_g\left(i\right)\right]=0\\ E\left[W_g\left(i\right)W_g{\left(i\right)}^T\right]={\mathrm{\Θ}}_g\end{array}\right.---\left(5\right)$

以姿态误差作为观测量构建卡尔曼滤波器的量测方程，形式如下：

Z_g＝F_gX_g+gV    (6)

上式中状态矢量X_g的定义与式(1)相同，观测量其中θ_x,θ_y,θ_z为系统解算姿态，为转台姿态；量测矩阵F_g为15阶方阵，表示为如下形式：

$F={\left[\begin{array}{cccc}1& 0& 0& \\ 0& 1& 0& 0_{3\×12}\\ 0& 0& 1& \end{array}\right]}_{3\×15}---\left(7\right)$

假设R_g为量测噪声方差阵，式(11)中V_g为服从正态分布N(0,R_g)的系统噪声，满足如下条件：

$\left\{\begin{array}{c}E\left[V_g\left(i\right)\right]=0\\ E\left[V_g\left(i\right)V_g{\left(i\right)}^T\right]=R_g\end{array}\right.---\left(8\right)$

步骤二：对系统状态方程进行离散化；

对步骤一建立的系统状态变量进行估计，需要对系统状态方程进行离散化；离散化采用泰勒级数展开，则：

$\mathrm{\Φ}(k+1,k)=I+\mathrm{TA}\left(k\right)+\frac{T^2}{2!}{A}^2\left(k\right)+\frac{T^3}{3!}{A}^3\left(k\right)+...---\left(9\right)$

其中：Φ(k+1,k)为状态一步转移矩阵、I为十五阶单位阵、A(k)为状态转移矩阵，T为滤波周期；

系统模型噪声的方差为：

$Q\left(k\right)=\mathrm{QT}+[\mathrm{AQ}+{\left(\mathrm{AQ}\right)}^T]\frac{T^2}{2!}+\left\{A\right[\mathrm{AQ}+{\left(\mathrm{AQ}\right)}^T]+A{[\mathrm{AQ}+{\mathrm{QA}}^T]}^T\}\frac{T^3}{3!}+...---\left(10\right)$

其中：Q(k)为离散系统噪声方差阵、Q连续系统噪声方程强度阵、A为状态转移矩阵；

步骤三：进行卡尔曼滤波状态估计；

对卡尔曼滤波器进行迭代，状态预测估计方程、方差预测方程、状态预测估计方程、方差迭代方程以及滤波增益方程表示为如下形式：

$P_{k/k-1}={\mathrm{\Φ}}_{k,k-1}{P}_{k-1}{\mathrm{\Φ}}_{k,k-1}^T+{\mathrm{\Γ}}_{k-1}{Q}_{k-1}{\mathrm{\Γ}}_{k-1}^T---\left(12\right)$

${\hat{X}}_k={\hat{X}}_{k/k-1}+K_k(Z_k-H_k{\hat{X}}_{k/k-1})---\left(13\right)$

P_k＝(I-K_kH_k)P_k/k-1    (14)

$K_k=P_{k/k-1}{H}_k^T{(H_k{P}_{k/k-1}{H}_k^T+R_k)}^{-1}---\left(15\right)$

最后估计得到冗余系统中斜置光纤陀螺零偏误差、标度因数误差以及失准角的标定结果，并对粗标定结果进行修正。