CN104344837A - 一种基于速度观测的冗余惯导系统加速度计系统级标定方法 - Google Patents

Abstract

一种基于速度观测的冗余惯导系统加速度计系统级标定方法，一、将惯导安装转台上，确定载体的初始位置参数；二、确定加速度计轴向与惯导本体坐标系安装关系；三、惯导预热，采集加速度计输出数据进行精标定；四、使惯组位于东北天位置静止不动，第一次标定参数修正；五、使惯组绕X轴转90°至东天南位置静止不动，第二次标定参数修正；六、使惯组绕Z轴转90°至天西南位置静止不动，第三次标定参数修正；七、使惯组绕Y轴转-90°至南西地位置静止不动，第四次标定参数修正；八、使惯组绕X轴转180°至南东天位置静止不动，第五次标定参数修正；九、对第五次标定参数修正，得到斜置加速度计高精度的标定参数零偏、标度因数、失准角结果。

Description

一种基于速度观测的冗余惯导系统加速度计系统级标定方法

技术领域

本发明属于惯性导航技术领域，特别是涉及一种基于速度观测的冗余惯导系统加速度计系统级标定方法。

背景技术

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

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

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

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

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

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

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

发明内容

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

本发明一种基于速度观测的冗余惯导系统加速度计系统级标定方法，该方法具体步骤如下：

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

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

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

加速度计输出的数据为载体相对于惯性参考系的比力f^b。

步骤四：使光纤捷联惯组位于东北天位置静止不动，进行第一次标定参数修正；

步骤五：使光纤捷联惯组绕X轴转90°至东天南位置静止不动，进行第二次标定参数修正；

步骤六：使光纤捷联惯组绕Z轴转90°至天西南位置静止不动，进行第三次标定参数修正；

步骤七：使光纤捷联惯组绕Y轴转-90°至南西地位置静止不动，进行第四次标定参数修正；

步骤八：使光纤捷联惯组绕X轴转180°至南东天位置静止不动，进行第五次标定参数修正；

步骤九：通过步骤八的第五次标定参数修正，得到斜置加速度计高精度的标定参数(零偏、标度因数、失准角)结果。

其中，步骤四至步骤八中采用基于卡尔曼滤波技术的误差标定方法，利用速度误差作为观测量，通过卡尔曼滤波迭代，估计加速度计零偏误差、标度因数误差及失准角，对加速度计粗标定结果进行修正。具体步骤如下：

步骤一：建立加速度计标定的系统状态方程和观测方程。

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

子惯导系统中加速度计误差项系统级精标定卡尔曼滤波器的状态方程为：

${\stackrel{.}{X}}_f=A_f{X}_f+W_f---\left(1\right)$

其中15维状态矢量包括东、北、天向速度误差δv^T；加速度计零偏残差矢量：ΔB_f＝[ΔB_f1,…,ΔB_fn]^T，加速度计标度因数误差残差矢量：加速度计安装失准角残差矢量： $\mathrm{\Δ}{\stackrel{\‾}{\mathrm{\α}}}_f={[{{\mathrm{\δ\α}}^{\′}}_{f1},...,{{\mathrm{\δ\α}}^{\′}}_{\mathrm{fn}}]}^T,\mathrm{\Δ}{\stackrel{\‾}{\mathrm{\β}}}_f={[{{\mathrm{\δ\β}}^{\′}}_{f1},...,{{\mathrm{\δ\β}}^{\′}}_{\mathrm{fn}}]}^T.$ 表示系统状态的微分，状态方程中15阶方阵(状态矩阵)可表示为如下形式：

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

其中A_f1A_f2可表示为如下形式：

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

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

表示地球自转角速度在导航坐标系n系下的投影(下标的x,y及z表示沿导航坐标系的三个坐标轴)。表示导航坐标系n系相对地球坐标系e系的角速度在n系下的投影。为子惯导系统的配置矩阵其中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)。

p_i＝[sin(α_i)cos(β_i) -cos(α_i)cos(β_i) 0]q_i＝[-cos(α_i)sin(β_i) -sin(α_i)sin(β_i) cos(β_i)]，

与h_i共同表征加速度计轴向与系统本体坐标系的安装关系(如图1、图2)。为转台所示捷联姿态矩阵。

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

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

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

Z_f＝F_fX_f+V_f    (6)

上式中状态矢量X_f的定义与式(1)相同，观测量Z_f＝[v_x,v_y,v_z]^T。量测矩阵F_f为15阶方阵，可表示为如下形式：

$F_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_f为量测噪声方差阵，式(6)中V_f为服从正态分布N(0,R_f)的系统噪声，满足如下条件：

$\left\{\begin{array}{c}E\left[V_f\left(i\right)\right]=0\\ E\left[V_f\left(i\right)V_f{\left(i\right)}^T\right]={\mathrm{R\Θ}}_f\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为滤波周期。

系统模型噪声的方差为：

$\mathrm{\Φ}\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，本发明是一种基于速度观测的冗余惯导系统加速度计系统级标定方法，它包括以下几个步骤：

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

步骤二：确定加速度计轴向与系统本体坐标系安装关系即安装角，计算安装矩阵；

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

加速度计输出的数据为载体相对于惯性参考系的比力f^b。

步骤四：使光纤捷联惯组分别位于东北天位置静止不动，进行第一次标定参数修正；

步骤五：使光纤捷联惯组绕X轴转90°至东天南位置静止不动，进行第二次标定参数修正；

步骤六：使光纤捷联惯组绕Z轴转90°至天西南位置静止不动，进行第三次标定参数修正；

步骤七：使光纤捷联惯组绕Y轴转-90°至南西地位置静止不动，进行第四次标定参数修正；

步骤八：使光纤捷联惯组绕X轴转180°至南东天位置静止不动，进行第五次标定参数修正；

每次转动过程持续5s，转动完成后静止1min。

步骤九：通过步骤八的第五次标定参数修正，得到斜置加速度计高精度的标定参数(零偏、标度因数、失准角)结果。

步骤四至步骤八中采用基于卡尔曼滤波技术的误差标定方法，利用速度误差作为观测量，通过卡尔曼滤波迭代，估计加速度计零偏误差、标度因数误差及失准角，对加速度计粗标定结果进行修正。具体步骤如下：

步骤一：建立加速度计标定的系统状态方程和观测方程。

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

子惯导系统中加速度计误差项系统级精标定卡尔曼滤波器的状态方程为：

${\stackrel{.}{X}}_f=A_f{X}_f+W_f---\left(1\right)$

其中15维状态矢量包括东、北、天向速度误差δv^T；加速度计零偏残差矢量：ΔB_f＝[ΔB_f1,…,ΔB_fn]^T，加速度计标度因数误差残差矢量：加速度计安装失准角残差矢量： $\mathrm{\Δ}{\stackrel{\‾}{\mathrm{\α}}}_f={[{{\mathrm{\δ\α}}^{\′}}_{f1},...,{{\mathrm{\δ\α}}^{\′}}_{\mathrm{fn}}]}^T,\mathrm{\Δ}{\stackrel{\‾}{\mathrm{\β}}}_f={[{{\mathrm{\δ\β}}^{\′}}_{f1},...,{{\mathrm{\δ\β}}^{\′}}_{\mathrm{fn}}]}^T.$ 表示系统状态的微分，状态方程中15阶方阵(状态矩阵)可表示为如下形式：

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

其中A_f1A_f2可表示为如下形式：

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

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

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

其中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)。

p_i＝[sin(α_i)cos(β_i) -cos(α_i)cos(β_i) 0]q_i＝[-cos(α_i)sin(β_i) -sin(α_i)sin(β_i) cos(β_i)]，与h_i共同表征加速度计轴向与系统本体坐标系的安装关系(如图1、图2所示)。为转台所示捷联姿态矩阵。

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

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

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

Z_f＝F_fX_f+V_f    (6)

上式中状态矢量X_f的定义与式(1)相同，观测量Z_f＝[v_x,v_y,v_z]^T。量测矩阵F_f为15阶方阵，可表示为如下形式：

$F_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_f为量测噪声方差阵，式(6)中V_f为服从正态分布N(0,R_f)的系统噪声，满足如下条件：

$\left\{\begin{array}{c}E\left[V_f\left(i\right)\right]=0\\ E\left[V_f\left(i\right)V_f{\left(i\right)}^T\right]=R_f\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为滤波周期。

系统模型噪声的方差为：

$\mathrm{\Φ}\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)系统级标定仿真条件及结论

仿真过程中的子惯导系统由图4所示六冗余斜置RFINS系统结构中的ABC光纤陀螺与加速度计构成，如图5所示。其中AB轴位于面M2，C轴位于面M1，其与坐标轴的夹角为α＝31°43'2.9″。

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

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

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

表1冗余加速度计标定参数

	加速度计
		零偏	300ug
标度因数(ppm)	100
		失准角A(角分)	1
失准角B(角分)	1
		噪声方差	10ug

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

从图6a-d中可以看出，当以系统速度误差作为观测量时(每次转动完成后系统姿态角替换为转台姿态角)，在前三个转动序列中加速度计误差逐渐收敛，当转至第四个位置后加速度计的各项误差逐渐逼近误差真值。

Claims (2)

1.一种基于速度观测的冗余惯导系统加速度计系统级标定方法，其特征在于：该方法具体步骤如下：

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

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

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

步骤四：使光纤捷联惯组位于东北天位置静止不动，进行第一次标定参数修正；

步骤五：使光纤捷联惯组绕X轴转90°至东天南位置静止不动，进行第二次标定参数修正；

步骤六：使光纤捷联惯组绕Z轴转90°至天西南位置静止不动，进行第三次标定参数修正；

步骤七：使光纤捷联惯组绕Y轴转-90°至南西地位置静止不动，进行第四次标定参数修正；

步骤八：使光纤捷联惯组绕X轴转180°至南东天位置静止不动，进行第五次标定参数修正；

步骤九：通过步骤八的第五次标定参数修正，得到斜置加速度计高精度的标定参数零偏、标度因数、失准角结果。

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

步骤一：建立加速度计标定的系统状态方程和观测方程；

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

子惯导系统中加速度计误差项系统级精标定卡尔曼滤波器的状态方程为：

${\stackrel{.}{X}}_f=A_f{X}_f+W_f---\left(1\right)$

其中15维状态矢量包括东、北、天向速度误差δv^T；加速度计零偏残差矢量：ΔB_f＝[ΔB_f1,…,ΔB_fn]^T，加速度计标度因数误差残差矢量：

$\mathrm{\Δ}{\stackrel{\‾}{K}}_f={[{\mathrm{\Δk}}_{f1},\·\·\·,{\mathrm{\Δk}}_{\mathrm{fn}}]}^T,$ 加速度计安装失准角残差矢量： $\mathrm{\Δ}{\stackrel{\‾}{\mathrm{\α}}}_f={[{{\mathrm{\δ\α}}^{\′}}_{f1},\·\·\·,{{\mathrm{\δ\α}}^{\′}}_{\mathrm{fn}}]}^T,$

表示系统状态的微分，状态方程中15阶方阵即状态矩阵表示为如下形式：

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

其中A_f1 A_f2表示为如下形式：

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

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

表示地球自转角速度在导航坐标系n系下的投影，下标的x,y及z表示沿导航坐标系的三个坐标轴，表示导航坐标系n系相对地球坐标系e系的角速度在n系下的投影，为子惯导系统的配置矩阵其中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平面的夹角p_i＝[sin(α_i)cos(β_i) -cos(α_i)cos(β_i) 0]q_i＝[-cos(α_i)sin(β_i) -sin(α_i)sin(β_i) cos(β_i)]，与h_i共同表征加速度计轴向与系统本体坐标系的安装关系，为转台所示捷联姿态矩阵；假设Θ_f为系统噪声方差阵，式(1)中W_f为服从正态分布N(0,Θ_f)的系统噪声，满足如下条件：

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

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

Z_f＝F_fX_f+V_f        (6)

上式中状态矢量X_f的定义与式(1)相同，观测量Z_f＝[v_x,v_y,v_z]^T，量测矩阵F_f为15阶方阵，表示为如下形式：

$F_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_f为量测噪声方差阵，式(6)中V_f为服从正态分布N(0,R_f)的系统噪声，满足如下条件：

$\left\{\begin{array}{c}E\left[V_f\left(i\right)\right]=0\\ E\left[V_f\left(i\right)V_f{\left(i\right)}^T\right]=R_f\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)$

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