CN104567930A - 一种能够估计和补偿机翼挠曲变形的传递对准方法 - Google Patents

Abstract

本发明公开了一种能够估计和补偿机翼挠曲变形的传递对准方法。包括以下步骤：将主惯导系统的导航信息即速度、姿态和位置装订给子惯导系统完成粗对准；利用主惯导以固定频率发送给子惯导系统的速度与姿态信息和子惯导系统的速度与姿态信息构造观测量；将机翼挠曲变形模型的相关变量引入状态方程中，采用“量测失准角+速度”匹配方式建立系统状态方程和观测方程；采用常规卡尔曼滤波进行迭代解算，估计出子惯导系统的姿态失准角和速度误差等信息，用其修正子惯导系统的速度和姿态，得到子惯导系统的初始导航信息，完成传递对准。本发明采用“量测失准角+速度”匹配方式进行传递对准，具有计算量小，形式简单、易于理解等优点。

Description

一种能够估计和补偿机翼挠曲变形的传递对准方法

技术领域

本发明属于惯性导航技术领域，特别是一种能够估计和补偿机翼挠曲变形的传递对准方法。

背景技术

传递对准是指载体航行时，载体上需要对准的子惯导利用已对准好的主惯导的信息进行对准的一种方法。在传递对准研究中，为了体现某一特定匹配方法的有效性，通常不考虑载体挠曲变形的影响，认为载体是一绝对的刚体，而在实际中这种假设是不成立的。

挠曲变形是指机翼本身并非刚体，在外力和外力矩、气动载荷、湍流的作用下产生变形和结构振动的现象。这些现象对传递对准的精度影响相当大，必须在方程编排中予以考虑。载体在飞行过程中的挠曲变形主要有两种，一种是在机动动作时产生的挠性慢变形，另一种是载体内部振源或阵风等使载体产生的振动变形。如果忽略载体的挠曲变形，利用卡尔曼滤波进行状态估计时就不用考虑挠曲变形角，从而减少滤波器的维数。但这种近似在载体振幅较大的时候产生的误差会很大。因此，如何估计和补偿机翼挠曲变形，是现在传递对准技术研究的重点。

文献《一种快速传递对准的方法》中提出了一种速度+姿态匹配方式的传递对准方法，但在考虑机翼挠曲变形时，这种匹配方式已不再适用。经过仿真发现，利用速度+姿态匹配方式进行卡尔曼滤波时，安装误差角的估计误差很大，显然这种方法对于安装误差角和机翼挠曲变形角的估计已经失效。

发明内容

本发明的目的在于提供一种对准速度快、精度高且能够估计和补偿机翼挠曲变形的传递对准方法。

实现本发明目的的技术解决方案为：一种能够估计和补偿机翼挠曲变形的传递对准方法，其特征在于，包括以下步骤：

步骤1，将主惯导系统的导航信息即速度、姿态和位置装订给子惯导系统完成粗对准；

步骤2，根据主惯导系统发送给子惯导系统的速度与姿态信息、子惯导系统自身的速度与姿态信息构造观测量；

步骤3，将机翼挠曲变形模型的变量引入状态方程中，采用速度加量测失准角的匹配方式建立系统状态方程和观测方程；

步骤4，采用卡尔曼滤波进行迭代解算，估计出子惯导系统的姿态失准角和速度误差，修正子惯导系统的速度和姿态，得到子惯导系统的初始导航信息，完成传递对准。

本发明与现有技术相比，其显著优点是：(1)相对于传统的速度+姿态匹配方式，本发明采用速度加量测失准角匹配方式估计和补偿机翼挠曲变形，具有对准速度快，对准精度高且可以更好的估计出机翼挠曲变形角；(2)采用速度加量测失准角匹配方式进行传递对准，具有形式简单、易于理解，可大大减小计算量的优点。

附图说明

图1为本发明能够估计和补偿机翼挠曲变形的传递对准方法工作流程图。

图2为本发明能够估计和补偿机翼挠曲变形的传递对准方法的X轴姿态角估计误差曲线。

图3为本发明能够估计和补偿机翼挠曲变形的传递对准方法的Y轴姿态角估计误差曲线。

图4为本发明能够估计和补偿机翼挠曲变形的传递对准方法的Z轴姿态角估计误差曲线。

图5为本发明能够估计和补偿机翼挠曲变形的传递对准方法的机翼挠曲变形实际值与估计值曲线。

具体实施方式

下面结合附图及具体实施例对本发明做进一步详细描述。

结合图1，本发明能够估计和补偿机翼挠曲变形的传递对准方法，包括以下步骤：

步骤1，将主惯导系统的导航信息即速度、姿态和位置装订给子惯导系统完成粗对准，即将主惯导系统的速度、姿态、位置信息传输给子惯导系统，子惯导系统利用主惯导系统传输的信息完成初始化工作；

步骤2，根据主惯导系统发送给子惯导系统的速度与姿态信息、子惯导系统自身的速度与姿态信息构造观测量；

利用主惯导以固定频率发送给子惯导系统的速度与姿态信息和子惯导系统的速度与姿态信息构造观测量

在国外文献中，使用主、子惯导姿态矩阵的乘积结果作为量测量的方法得到了广泛的使用,本文将此匹配方法称为量测失准角匹配法。此方法的起源是在Kain等提出的量测失准角(机体系m到计算弹体系s'的相对姿态失准角)的基础上发展和推导的，其中引入计算弹体系s'，s'的主要作用是将数学平台坐标系n'变换到s'系，s'系对应的理想坐标系是s系，n'系对应的理想坐标系是n系，在理想情况下s'系和s系重合，n'和n系重合；设量测失准角为小量，定义如下：

由于存在机翼挠曲变形角，则推导公式如下：

$\begin{array}{c}{C}_m^{s^{\′}}=C_n^{s^{\′}}{C}_m^n=C_n^s{C}_{n^{\′}}^n{C}_m^n\\ =C_m^s{C}_n^m{C}_{n^{\′}}^n{C}_m^n\\ =(I-(\mathrm{\λ}+\mathrm{\θ})\×)C_n^m(I+\mathrm{\φ}\×)C_m^n\\ =(I-(\mathrm{\λ}+\mathrm{\θ})\×)(I+C_n^m(\mathrm{\φ}\×)C_m^n)\end{array}$

式中，为主惯导到计算弹体系间的转换矩阵，为主惯导到实际弹体系间的转换矩阵，为计算弹体系中子惯导系统姿态矩阵的转置，为实际弹体系中子惯导系统姿态矩阵的转置，为主惯导系统姿态矩阵，为主惯导系统姿态矩阵的转置，为数学平台坐标系到导航坐标系间的转换矩阵，I为单位矩阵，λ＝[λ_x λ_y λ_z]^T为主子惯导间的安装误差角，θ＝[θ_x θ_y θ_z]^T为机翼挠曲变形角矢量，φ＝[φ_x φ_y φ_z]^T为子惯导系统的平台失准角，下标x，y，z分别表示机体坐标系下x轴，y轴，z轴的分量。并且φ×表示

$(\mathrm{\φ}\×)=\left[\begin{array}{ccc}0& -{\mathrm{\φ}}_z& {\mathrm{\φ}}_y\\ {\mathrm{\φ}}_z& 0& -{\mathrm{\φ}}_x\\ -{\mathrm{\φ}}_y& {\mathrm{\φ}}_x& 0\end{array}\right]$

将上式右边展开，省略二阶小量，得

$C_m^{s^{\′}}=I+(\left(C_n^m\mathrm{\φ}\right)\×)-(\mathrm{\λ}-\mathrm{\θ})\×$

根据反对称阵性质，得

因此采用速度加量测失准角匹配算法的姿态观测方程；速度观测方程由主、子惯导系统东、北、天三个方向的速度差作为观测量。

步骤3，将机翼挠曲变形模型的变量引入状态方程中，采用速度加量测失准角的匹配方式建立系统状态方程和观测方程，具体如下：

第1步，根据捷联惯导系统的力学编排建立惯导误差模型，得到系统的姿态误差方程、速度误差方程及惯性器件误差方程；采用二阶马尔科夫模型作为机翼挠曲变形的模型，具体如下：

(1)由于传递对准过程中存在机翼挠曲变形的影响，因此需对机翼挠曲变形进行补偿，本发明采用二阶马尔可夫过程作为机翼挠曲变形的模型，通过将机翼挠曲变形模型的相关变量引入卡尔曼滤波器状态变量中，估计和补偿机翼挠曲变形，二阶马尔可夫模型为：

${\stackrel{\·\·}{\mathrm{\θ}}}_i+2{\mathrm{\β}}_i{\stackrel{\·}{\mathrm{\θ}}}_i+{\mathrm{\β}}_i^2{\mathrm{\θ}}_i=,(i=x,y,z)$

即：

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\θ}}}_i={\mathrm{\ω}}_i\\ {\stackrel{\·}{\mathrm{\ω}}}_i=-{\mathrm{\β}}_i{\mathrm{\θ}}_i-2{\mathrm{\β}}_i{\mathrm{\ω}}_i+{\mathrm{\η}}_i\end{array}\right.(i=x,y,z)$

式中，θ＝[θ_x θ_y θ_z]^T为机翼挠曲变形角矢量，机翼挠曲变形角均方差σ＝[σ_x σ_y σ_z]^T；ω＝[ω_x ω_y ω_z]^T为动态挠曲变形角速率矢量；η＝[η_x η_y η_z]^T为白噪声，η的频谱密度 $Q_{\mathrm{\η}}={\left[\begin{array}{ccc}{Q}_{{\mathrm{\η}}_x}& Q_{{\mathrm{\η}}_y}& Q_{{\mathrm{\η}}_z}\end{array}\right]}^T,$ 即η～N(0,Q_η)；β＝[β_x β_y β_z]^T为常数；其中，Q_η、σ和β之间的关系为每个随机过程的相关时间τ_i与相应β_i之间的关系为β_i＝2.146/τ_i；下标x，y，z分别表示机体坐标系下x轴，y轴，z轴的分量；

(2)引入机翼挠曲变形模型变量，采用速度加量测失准角的匹配方式建立系统状态方程如下：

$\left\{\begin{array}{c}\mathrm{\δ}\stackrel{\·}{V^n}=f^n\×{\mathrm{\φ}}^n-2(\mathrm{\δ}{\mathrm{\ω}}_{\mathrm{ie}}^n+\mathrm{\δ}{\mathrm{\ω}}_{\mathrm{en}}^n)\×V^n\\ -2({\mathrm{\ω}}_{\mathrm{ie}}^n+{\mathrm{\ω}}_{\mathrm{en}}^n)\×\mathrm{\δ}{V}^n+{\▿}^n\\ \stackrel{\·}{\mathrm{\φ}=-{\mathrm{\ω}}_{\mathrm{in}}^n\×\mathrm{\φ}-C_b^n\mathrm{\δ}{\mathrm{\ω}}_{\mathrm{ib}}^b+\mathrm{\δ}{\mathrm{\ω}}_{\mathrm{in}}^n}\\ \stackrel{\·}{\mathrm{\ϵ}}=0\\ \stackrel{\·}{\▿}=0\\ \stackrel{\·}{\mathrm{\λ}}=0\\ \stackrel{\·}{\mathrm{\θ}}=\mathrm{\ω}\\ \stackrel{\·}{\mathrm{\ω}}=-{\mathrm{\β}}^2\mathrm{\θ}-2\mathrm{\β\ω}+\mathrm{\η}\end{array}\right.$

式中，φ＝[φ_x φ_y φ_z]为子惯导系统的失准角， 为地球自传角速率在导航坐标下的投影，为导航系相对地球系的旋转角速率，为子惯导的姿态矩阵，δV＝[δV_e δV_n δV_u]分别为子惯导系统的东、北、天方向的速度误差， $V^n=\left[\begin{array}{ccc}{V}_e^n& V_n^n& V_u^n\end{array}\right]$ 分别为子惯导系统的东、北、天方向的速度，w_ie为地球转动角速率， $f^n=\left[\begin{array}{ccc}{f}_e^n& f_n^n& f_u^n\end{array}\right]$ 是主惯导系统加速度计测量到的比力通过坐标转换后得到的导航系下东、北、天三个方向上的比力分量，ε＝[ε_x ε_y ε_z]为子惯导陀螺在体坐标系下的随机常值漂移，▽＝[▽_x ▽_y ▽_z]为子惯导加速度计在体坐标系下的随机常值偏置，λ＝[λ_x λ_y λ_z]为子惯导系统相对主惯导系统的安装误差角在体坐标系下的分量，下标x，y，z分别表示机体坐标系下x轴，y轴，z轴的分量。

第2步，根据惯导误差模型、系统误差方程以及机翼挠曲变形的模型建立捷联惯导系统的状态方程，如下式所示：

$\stackrel{\·}{X}=\mathrm{FX}+\mathrm{GW}$

式中，X为系统状态矢量，F为系统状态转移矩阵，G为系统噪声驱动矩阵，W为系统噪声矢量，具体如下：

X为系统状态矢量：

$X={\left[\begin{array}{ccccccc}\mathrm{\φ}& \mathrm{\δV}& \mathrm{\ϵ}& \▿& \mathrm{\λ}& \mathrm{\θ}& \mathrm{\ω}\end{array}\right]}_{1\×21}^T$

F为系统状态转移矩阵：

系统状态转移矩阵中矩阵块F¹与F²如下所示：

$F^1=\left[\begin{array}{ccc}0& {\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_e^n}{R_N+h}\tan L& -({\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_e^n}{R_N+h})\\ -({\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_e^n}{R_N+h}\tan L)& 0& \frac{V_n^n}{R_M+h}\\ {\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_e^n}{R_N+h}& \frac{V_n^n}{R_M+h}& 0\\ 0& -f_u^n& f_n^n\\ f_u^n& 0& -f_e^n\\ -f_n^n& f_e^n& 0\end{array}\right]$

$F^2=\left[\begin{array}{ccc}0& -\frac{1}{R_M+h}& 0\\ \frac{1}{R_N+h}& 0& 0\\ \frac{\tan L}{R_N+h}& 0& 0\\ 0& 2{\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_e^n}{R_N+h}\tan L& -(2{\mathrm{\ω}}_{\mathrm{ie}}\mathrm{cisL}+\frac{V_e^n}{R_N+j})\\ -({\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_e^n}{R_N+h}\tan L)& 0& -\frac{V_n^n}{R_M+h}\\ 2{\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_e^n}{R_N+h}& \frac{V_n^n}{R_M+h}& 0\end{array}\right]$

$F^3=\left[\begin{array}{ccc}-{\mathrm{\β}}_x^2& 0& 0\\ 0& -{\mathrm{\β}}_y^2& 0\\ 0& 0& -{\mathrm{\β}}_z^2\end{array}\right]F^4=\left[\begin{array}{ccc}-2{\mathrm{\β}}_x& 0& 0\\ 0& -2{\mathrm{\β}}_y& 0\\ 0& 0& -2{\mathrm{\β}}_z\end{array}\right]$

其中，R_M为子午圈曲率半径；R_N为卯酉圈曲率半径；h为载体当前的高度；L为载体当前所处位置的纬度；

捷联惯导系统的卡尔曼滤波器观测方程如下所示：

Z＝HX+V

式中，为观测量，分别为俯仰、滚转、偏航方向上的姿态观测量；δV_e，δV_n，δV_u为东向、北向和天向方向上的主、子惯导速度差；V是观测噪声且为均值为零的高斯白噪声，V的协方差为E[VV^T]＝R，R为观测噪声方差阵；

结合观测方程得系统的观测矩阵为：

$H={\left[\begin{array}{ccccccc}-C_n^m& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& I_{3\×3}& I_{3\×3}& 0_{3\×3}\\ 0_{3\×3}& I_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}\end{array}\right]}_{6\×21}$

其中：I为单位矩阵，为主惯导姿态矩阵的转置。

步骤4，采用卡尔曼滤波进行迭代解算，估计出子惯导系统的姿态失准角和速度误差，修正子惯导系统的速度和姿态，得到子惯导系统的初始导航信息，完成传递对准，具体为：

根据建立的系统状态方程、系统观测方程以及系统的观测量进行卡尔曼滤波迭代解算

系统离散化后的状态方程和量测方程为：

X_k＝Φ_k,k-1X_k-1+Γ_k-1W_k-1

Z_k＝H_kX_k+V_k

其中，X_k表示被估计状态，Φ_k,k-1为t_k-1时刻的一步转移矩阵，Γ_k-1为k-1时刻系统噪声驱动阵，W_k为系统激励噪声序列，Z_k为观测量，H_k为量测阵，V_k为量测噪声序列。根据适用于离散系统的卡尔曼滤波器基本方程的卡尔曼滤波算法具体如下：

状态一步预测方程：

X_k/k-1＝φ_k,k-1X_k-1

其中，X_K/K-1为k-1时刻到k时刻的系统状态一步预测值，X_k-1为k-1时刻系统状态估计值，φ_k,k-1为k-1时刻到k时刻的系统状态转移矩阵；

一步预测均方误差方程：

$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$

其中，P_k/k-1为k-1时刻到k时刻的系统状态协方差阵，P_k-1为k-1时刻的系统状态协方差阵，Q_k-1为k-1时刻系统噪声矩阵；

最优滤波增益方程：

$K_k=P_{k|k-1}{H}_k^T{[H_k{P}_{k|k-1}{H}_k^T+R_k]}^{-1}$

其中，K_k为k时刻系统增益矩阵，H_k为k时刻系统量测矩阵，R_k为k时刻系统量测噪声矩阵；

状态估计方程：

X_k＝X_k/k-1+K_k(Z_k-H_kX_k/k-1)

估计均方误差方程：

$P_k=(I-K_k{H}_k)P_{k/k-1}{(I-K_k{H}_k)}^T+K_k{R}_k{K}_{k-1}^T$

其中，Q_k是系统过程噪声方差阵，P_k为均方误差阵。

联立状态一步预测方程、一步预测均方误差方程、最优滤波增益方程、状态估计方程和估计均方误差方程可估计出状态矢量 $x={\left[\begin{array}{ccccccc}\mathrm{\φ}& \mathrm{\δV}& \mathrm{\ϵ}& \▿& \mathrm{\λ}& \mathrm{\θ}& \mathrm{\ω}\end{array}\right]}_{1\×21}^T$ 的状态估计值。利用状态矢量 $x={\left[\begin{array}{ccccccc}\mathrm{\φ}& \mathrm{\δV}& \mathrm{\ϵ}& \▿& \mathrm{\λ}& \mathrm{\θ}& \mathrm{\ω}\end{array}\right]}_{1\×21}^T$ 中的姿态失准角φ的估计值和速度误差δV的估计值对子惯导系统的导航信息进行修正，得到子惯导系统的初始导航信息，从而完成传递对准。

实施例1

为了进一步对上述方法进行说明，充分展现出该方法的可靠性与精确性，进行了数字仿真试验。

(1)仿真初始条件及参数设置

11)载体初始速度为东向速度100m/s、北向速度100m/s、天向速度0m/s；载体初始位置：纬度为30°，经度为110°，高度为1000m；载体初始姿态为：俯仰角为0°，滚转角为0°，偏航角为-45°。

12)子惯导加速度计常值偏置设为5mg；子惯导陀螺常值漂移设为10(deg/h)；子惯导在三个轴向相对于主惯导的安装角设为2°、1.5°、1°。

13)机翼挠曲变形相关时间τ_x＝60s，τ_y＝80s，τ_z＝70s，二阶马尔可夫过程方差 ${\mathrm{\σ}}_i^2=0.1{\mathrm{rad}}^2/s^4(i=x,y,z).$

14)仿真总时间为50s，滤波周期设为0.1s，在对准过程中采用爬升+摇翼机动方式。

(2)仿真试验结果及分析

图2-4给出了补偿机翼挠曲变形角和不补偿机翼挠曲变形角时子惯导的姿态角估计误差曲线，图5给出了机翼挠曲变形角的实际值与估计值。从图2-4中可以看出，不补偿机翼挠曲变形角时，姿态角估计误差很大，分别为0.32°，0.15°，0.3°，姿态角估计精度很差，不能满足传递对准高精度的要求；补偿机翼挠曲变形角时，姿态角估计误差很小，相对于不补偿机翼挠曲变形角时有明显提高，姿态角估计误差分别为-0.02°，-0.01°，-0.01°，完全满足传递对准对精度的要求。并且通过图5可以看出，本发明可以对机翼挠曲变形角进行有效估计。

Claims (6)

1.一种能够估计和补偿机翼挠曲变形的传递对准方法，其特征在于，包括以下步骤：

步骤1，将主惯导系统的导航信息即速度、姿态和位置装订给子惯导系统完成粗对准；

步骤2，根据主惯导系统发送给子惯导系统的速度与姿态信息、子惯导系统自身的速度与姿态信息构造观测量；

步骤3，将机翼挠曲变形模型的变量引入状态方程中，采用速度加量测失准角的匹配方式建立系统状态方程和观测方程；

步骤4，采用卡尔曼滤波进行迭代解算，估计出子惯导系统的姿态失准角和速度误差，修正子惯导系统的速度和姿态，得到子惯导系统的初始导航信息，完成传递对准。

2.根据权利要求1所述的能够估计和补偿机翼挠曲变形的传递对准方法，其特征在于，步骤1所述将主惯导系统的导航信息即速度、姿态和位置装订给子惯导系统完成粗对准，即将主惯导系统的速度、姿态、位置信息传输给子惯导系统，子惯导系统利用主惯导系统传输的信息完成初始化工作。

3.根据权利要求1所述的能够估计和补偿机翼挠曲变形的传递对准方法，其特征在于，步骤2所述根据主惯导系统发送给子惯导系统的速度与姿态信息、子惯导系统自身的速度与姿态信息构造观测量，具体如下：

引入计算弹体系s'，s'的主要作用是将数学平台坐标系n'变换到s'系，s'系对应的理想坐标系是s系，n'系对应的理想坐标系是n系，在理想情况下s'系和s系重合，n'和n系重合；设量测失准角为小量，定义如下：

由于存在机翼挠曲变形角，则推导公式如下：

$\begin{array}{c}{C}_m^{s\′}=C_n^{s\′}{C}_m^n=C_n^s{C}_{n\′}^s{C}_{n\′}^n{C}_m^n\\ =C_m^s{C}_n^m{C}_{n\′}^n{C}_m^n\\ =(I-(\mathrm{\λ}+\mathrm{\θ})\×)C_n^m(I+\mathrm{\φ}\×)C_m^n\\ =(I-(\mathrm{\λ}+\mathrm{\θ})\×)(I+C_n^m(\mathrm{\φ}\×)C_m^n)\end{array}$

式中，为主惯导到计算弹体系间的转换矩阵，为主惯导到实际弹体系间的转换矩阵，为计算弹体系中子惯导系统姿态矩阵的转置，为实际弹体系中子惯导系统姿态矩阵的转置，为主惯导系统姿态矩阵，为主惯导系统姿态矩阵的转置，为数学平台坐标系到导航坐标系间的转换矩阵，I为单位矩阵，λ＝[λ_x λ_y λ_z]^T为主子惯导间的安装误差角，θ＝[θ_x θ_y θ_z]^T为机翼挠曲变形角矢量，φ＝[φ_x φ_y φ_z]^T为子惯导系统的平台失准角，下标x，y，z分别表示机体坐标系下x轴，y轴，z轴的分量，并且

$(\mathrm{\φ}\×)=\left[\begin{array}{ccc}0& -{\mathrm{\φ}}_z& {\mathrm{\φ}}_y\\ {\mathrm{\φ}}_z& 0& -{\mathrm{\φ}}_x\\ -{\mathrm{\φ}}_y& {\mathrm{\φ}}_x& 0\end{array}\right]$

将上式右边展开，省略二阶小量，得

$C_m^{s^{\′}}=I+(\left(C_n^m\mathrm{\φ}\right)\×)-(\mathrm{\λ}+\mathrm{\θ})\×$

根据反对称阵性质，得

因此采用速度加量测失准角匹配算法的姿态观测方程；速度观测方程由主、子惯导系统东、北、天三个方向的速度差作为观测量。

4.根据权利要求1所述的能够估计和补偿机翼挠曲变形的传递对准方法，其特征在于，步骤3所述采用速度加量测失准角的匹配方式建立系统状态方程和观测方程，具体如下：

第1步，根据捷联惯导系统的力学编排建立惯导误差模型，得到系统的姿态误差方程、速度误差方程及惯性器件误差方程；采用二阶马尔科夫模型作为机翼挠曲变形的模型；

第2步，根据惯导误差模型、系统误差方程以及机翼挠曲变形的模型建立捷联惯导系统的状态方程。

5.根据权利要求4所述的能够估计和补偿机翼挠曲变形的传递对准方法，其特征在于，第1步所述根据捷联惯导系统的力学编排建立惯导误差模型，得到系统的姿态误差方程、速度误差方程及惯性器件误差方程；采用二阶马尔科夫模型作为机翼挠曲变形的模型，具体如下：

(1)采用二阶马尔科夫模型作为机翼挠曲变形的模型，通过将机翼挠曲变形模型的变量引入卡尔曼滤波器状态变量中，估计和补偿机翼挠曲变形，二阶马尔可夫模型如下：

${\stackrel{\·\·}{\mathrm{\θ}}}_i+2{\mathrm{\β}}_i{\stackrel{\·}{\mathrm{\θ}}}_i+{\mathrm{\β}}_i^2{\mathrm{\θ}}_i={\mathrm{\η}}_i,(i=x,y,z)$

即：

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\θ}}}_i={\mathrm{\ω}}_i\\ {\stackrel{\·}{\mathrm{\ω}}}_i=-{\mathrm{\β}}_i{\mathrm{\θ}}_i-2{\mathrm{\β}}_i{\mathrm{\ω}}_i+{\mathrm{\η}}_i\end{array}\right.,(i=x,y,z)$

式中，θ＝[θ_x θ_y θ_z]^T为机翼挠曲变形角矢量，机翼挠曲变形角均方差σ＝[σ_x σ_y σ_z]^T；ω＝[ω_x ω_y ω_z]^T为动态挠曲变形角速率矢量；η＝[η_x η_y η_z]^T为白噪声，η的频谱密度即η～N(0,Q_η)；β＝[β_x β_y β_z]^T为常数；其中，Q_η、σ和β之间的关系为每个随机过程的相关时间τ_i与相应β_i之间的关系为β_i＝2.146/τ_i；下标x，y，z分别表示机体坐标系下x轴，y轴，z轴的分量；

(2)引入机翼挠曲变形模型变量，采用速度加量测失准角的匹配方式建立系统状态方程如下：

$\left\{\begin{array}{c}\mathrm{\δ}{\stackrel{\·}{V}}^n=f^n\×{\mathrm{\φ}}^n-2(\mathrm{\δ}{\mathrm{\ω}}_{\mathrm{ie}}^n+\mathrm{\δ}{\mathrm{\ω}}_{\mathrm{en}}^n)\×V^n\\ -2({\mathrm{\ω}}_{\mathrm{ie}}^n+{\mathrm{\ω}}_{\mathrm{en}}^n)\×\mathrm{\δ}{V}^n+{\▿}^n\\ \stackrel{\·}{\mathrm{\φ}}=-{\mathrm{\ω}}_{\mathrm{in}}^n\×\mathrm{\φ}-C_b^n\mathrm{\δ}{\mathrm{\ω}}_{\mathrm{ib}}^b+\mathrm{\δ}{\mathrm{\ω}}_{\mathrm{in}}^n\\ \stackrel{\·}{\mathrm{\ϵ}}=0\\ \stackrel{\·}{\▿}=0\\ \stackrel{\·}{\mathrm{\λ}}=0\\ \stackrel{\·}{\mathrm{\θ}}=\mathrm{\ω}\\ \stackrel{\·}{\mathrm{\ω}}=-{\mathrm{\β}}^2\mathrm{\θ}-2\mathrm{\β\ω}+\mathrm{\η}\end{array}\right.$

式中，φ＝[φ_x φ_y φ_z]为子惯导系统的失准角， 为地球自传角速率在导航坐标下的投影，为导航系相对地球系的旋转角速率，为子惯导的姿态矩阵，δV＝[δV_e δV_n δV_u]分别为子惯导系统的东、北、天方向的速度误差，分别为子惯导系统的东、北、天方向的速度，w_ie为地球转动角速率，是主惯导系统加速度计测量到的比力通过坐标转换后得到的导航系下东、北、天三个方向上的比力分量，ε＝[ε_x ε_y ε_z]为子惯导陀螺在体坐标系下的随机常值漂移，为子惯导加速度计在体坐标系下的随机常值偏置，λ＝[λ_x λ_y λ_z]为子惯导系统相对主惯导系统的安装误差角在体坐标系下的分量，下标x，y，z分别表示机体坐标系下x轴，y轴，z轴的分量。

6.根据权利要求4所述的能够估计和补偿机翼挠曲变形的传递对准方法，其特征在于，第2步所述根据惯导误差模型、系统误差方程以及机翼挠曲变形的模型建立捷联惯导系统的状态方程，如下式所示：

$\stackrel{\·}{X}=\mathrm{FX}+\mathrm{GW}$

式中，X为系统状态矢量，F为系统状态转移矩阵，G为系统噪声驱动矩阵，W为系统噪声矢量，具体如下：

X为系统状态矢量：

$X={\left[\begin{array}{ccccccc}\mathrm{\φ}& \mathrm{\δV}& \mathrm{\ϵ}& \▿& \mathrm{\λ}& \mathrm{\θ}& \mathrm{\ω}\end{array}\right]}_{1\×21}^T$

F为系统状态转移矩阵：

$F={\left[\begin{array}{ccccccc}{F}_{6\×3}^1& F_{6\×3}^2& C_b^n& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}\\ & & 0_{3\×3}& C_b^n& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}\\ 0_{9\×3}& 0_{9\×3}& 0_{9\×3}& 0_{9\×3}& 0_{9\×3}& 0_{9\×3}& 0_{9\×3}\\ 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& I_{3\×3}\\ 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& F_{3\×3}^3& F_{3\×3}^4\end{array}\right]}_{21\×21}$

系统状态转移矩阵中矩阵块F¹与F²如下所示：

$F^1=\left[\begin{array}{ccc}0& {\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_e^n}{R_N+h}\tan L& -({\mathrm{\ω}}_{\mathrm{ie}}\mathrm{cps\; L}+\frac{V_e^n}{R_N+h})\\ -({\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_e^n}{R_N+h}\tan L)& 0& -\frac{V_n^n}{R_M+h}\\ {\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_e^n}{R_N+h}& \frac{V_n^n}{R_M+h}& 0\\ 0& -f_u^n& f_n^n\\ f_u^n& 0& -f_e^n\\ -f_n^n& f_e^n& 0\end{array}\right]$

$F^2=\left[\begin{array}{ccc}0& -\frac{1}{R_M+h}& 0\\ \frac{1}{R_N+h}& 0& 0\\ \frac{\tan L}{R_N+h}& 0& 0\\ 0& 2{\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_e^n}{R_N+h}\tan L& -({2\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_e^n}{R_N+h})\\ -({2\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_e^n}{R_N+h}\tan L)& 0& -\frac{V_n^n}{R_M+h}\\ {2\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{v_e^n}{R_N+h}& \frac{V_n^n}{R_M+h}& 0\end{array}\right]$

$F^3=\left[\begin{array}{ccc}-{\mathrm{\β}}_x^2& 0& 0\\ 0& -{\mathrm{\β}}_y^2& 0\\ 0& 0& -{\mathrm{\β}}_z^2\end{array}\right],F^4=\left[\begin{array}{ccc}-2{\mathrm{\β}}_x& 0& 0\\ 0& -2{\mathrm{\β}}_y& 0\\ 0& 0& -2{\mathrm{\β}}_z\end{array}\right]$

其中，R_M为子午圈曲率半径；R_N为卯酉圈曲率半径；h为载体当前的高度；L为载体当前所处位置的纬度；

捷联惯导系统的卡尔曼滤波器观测方程如下所示：

Z＝HX+V

式中，为观测量，分别为俯仰、滚转、偏航方向上的姿态观测量；δV_e，δV_n，δV_u为东向、北向和天向方向上的主、子惯导速度差；V是观测噪声且为均值为零的高斯白噪声，V的协方差为E[VV^T]＝R，R为观测噪声方差阵；

结合观测方程得系统的观测矩阵为：

$H={\left[\begin{array}{ccccccc}-C_n^m& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& I_{3\×3}& I_{3\×3}& 0_{3\×3}\\ 0_{3\×3}& I_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}\end{array}\right]}_{6\×21}$

其中：I为单位矩阵，为主惯导姿态矩阵的转置。