CN102621565A - 一种机载分布式pos的传递对准方法 - Google Patents

Abstract

一种机载分布式POS的传递对准方法，利用力学的方法对载机机翼的挠曲运动进行建模，并将由机翼挠曲运动产生的挠曲变形角和挠曲变形角速度增广为卡尔曼滤波的状态变量。在此基础上，采用“速度+姿态”的匹配方式设计卡尔曼滤波器，利用主惯性测量单元与GPS进行信息融合后得到的速度、姿态基准信息对子惯性测量单元进行传递对准，最终获得每个子惯性测量单元安装点的经校正后的速度和姿态信息。本发明具有自主性强、精度高的特点，可用于提高载机存在挠曲变形时分布式POS的传递对准精度。

Description

一种机载分布式POS的传递对准方法

技术领域

本发明涉及一种机载分布式POS的传递对准方法，可用于提高载机存在挠曲变形时分布式POS的传递对准精度。

背景技术

航空遥感是以飞机为观测平台，利用各种成像载荷获取地球表面或表层的大范围、高精度图像的战略技术，对国家经济建设和国家安全具有重大意义。随着我国经济的发展，航空遥感不断向着高分辨率、高精度、多遥感载荷集成及阵列载荷方向发展。

为实现遥感载荷的高精度成像，飞机需做匀速直线运动，但飞机在实际飞行中受气流扰动、飞控系统误差等因素影响，必然偏离理想的匀速直线运动，导致遥感成像分辨率和精度严重下降。因此，高精度航空遥感必须进行运动补偿。位置姿态测量系统(Position and Orientation System，POS)是一种特殊的惯性/卫星组合测量系统，可精确测量遥感载荷中心的位置、速度和姿态等运动参数，是高分辨率航空遥感系统实现运动补偿的关键设备。

但是对于装备了多个观测载荷的高性能航空遥感系统，由于多个观测载荷安装在飞机的不同位置，飞机弹性变形导致载荷间的空间相对关系发生变化。此时，采用传统的单一POS无法实现多点的高精度位置姿态测量。因此，必须建立起高精度分布式时空基准系统(分布式POS)为高性能航空遥感系统中所有载荷提供高精度的时间、空间信息。分布式POS由主导航系统和子惯性测量单元组成。其中主导航系统即为采用高精度惯性测量单元的POS，可精确测量载体的运动参数；子惯性测量单元采用低精度的惯性测量单元。主惯性测量单元(导航解算后称为主惯导)通常安装在机腹中央，子惯性测量单元(导航解算后称为子惯导)安装在机翼两侧的不同载荷附近，用于测量载荷中心的运动参数。

但是，子惯导精度不高，测量误差随时间累积，需要主惯导不断对其进行传递对准，以达到高精度测量。在传递对准过程中，主惯导与GPS信息进行融合后所提供的信息可视为高精度的参考信息。利用精度较高的主惯导来校准子惯导的传递对准技术就是设法估计出失准角并消除其影响，提高对准精度。因此传递对准的核心是确定主惯性测量单元与子惯性测量单元之间的姿态失准角。此时，在飞机飞行过程中，主惯性测量单元与子惯性测量单元之间的姿态失准角包括两部分：固定安装误差角和机翼弹性运动引起的挠曲变形角。其中，固定安装误差角是常值，并且不受外界干扰的影响，容易确定。而挠曲变形角是由载机机动、载机内部振源或阵风等使载机机翼产生挠曲变形而产生的动态变形角，不易确定，对传递对准的精度影响较大。现有的对挠曲变形角进行补偿的方法有两种：一种是用经验建模的方法，将机翼的挠曲变形考虑为二阶或三阶的马尔科夫过程，并将挠曲变形角增广为卡尔曼滤波的状态变量，通过卡尔曼滤波估计出该挠曲变形角并进行补偿。用此方法时，模型的有关参数全凭经验设定。并且，对于安装在机翼不同位置的子惯性测量单元，经验模型和模型参数也完全相同。而对于安装在机翼不同位置的子惯性测量单元来说，挠曲运动变形是有明显差别的。靠近机腹处的挠曲变形程度较低，而靠近机翼尖端处的挠曲变形程度较高。这种不加区分的对整个机翼建立统一的挠曲变形运动模型是不够精确的，必然影响了传递对准的精度；另一种补偿方法是采用加大卡尔曼滤波器过程噪声矩阵来减小机翼变形和振动的影响。此方法是依据真实模型的协方差分析结果来确定注入白噪声的强度，即通过加大过程噪声来补偿建模的挠曲变形。该方法除了可以补偿挠曲运动外，还可以增加滤波器的鲁棒性，但是以降低传递对准的精度为代价的。对于机载分布式POS来说，其用途是运动补偿，因此对传递对准精度要求较高，现有的对机翼挠曲变形进行补偿的两种方法都很难保证对准精度。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提出一种机载分布式POS的传递对准方法，该方法可提高载机机翼存在挠曲变形时分布式POS的传递对准精度。

本发明的技术解决方案为：一种机载分布式POS的传递对准方法，具体步骤如下：

机载分布式POS的传递对准是利用卡尔曼滤波技术估计出主、子惯性测量单元间的失准角，包括：固定安装误差角λ和挠曲变形角μ。利用主惯性测量单元与GPS进行信息融合后得到的基准信息，对主、子惯性测量单元的失准角进行校正，最终获得每个子惯性测量单元安装点的高精度速度、姿态信息，完成传递对准。

(1)利用力学的方法建立载机机翼挠曲变形的运动模型，得到挠曲变形角μ和挠曲变形角速度

的运动方程；

(2)采用卡尔曼滤波技术进行主惯性测量单元与GPS的信息融合，获得主惯性测量单元安装点的速度和姿态基准信息；

(3)采用“速度+姿态”匹配方式建立卡尔曼滤波模型，将机翼由于挠曲变形运动产生的挠曲变形角μ和挠曲变形角速度

，增广为卡尔曼滤波的状态变量。通过卡尔曼滤波估计出主、子惯性测量单元间速度误差，固定安装误差角λ和挠曲变形角μ；

(4)利用步骤(2)得到的基准信息和步骤(3)估计出的主、子惯性测量单元间的姿态失准角，对每个子惯性测量单元的速度和姿态进行校正，最终获得每个子惯性测量单元安装点的高精度的速度和姿态信息。

上述所采用的对机翼挠曲变形进行建模的方法为一种基于力学的建模方法，该方法的具体步骤为：

(1)利用ANSYS软件对机翼结构进行模态分析，得到机翼的第一阶弯曲模态频率ω_m和扭转模态频率ω_n；

(2)建立机翼弯曲和扭转运动与时间相关项的数学方程为：

$\left\{\begin{array}{c}\stackrel{\·\·}{m}\left(t\right)+2{\mathrm{\ξ}}_m{\mathrm{\ω}}_m\stackrel{\·}{m}\left(t\right)+{{\mathrm{\ω}}_m}^{2}m\left(t\right)=f_m\left(t\right)\\ \stackrel{\·\·}{n}\left(t\right)+2{\mathrm{\ξ}}_n{\mathrm{\ω}}_n\stackrel{\·}{n}\left(t\right)+{{\mathrm{\ω}}_n}^{2}n\left(t\right)=f_n\left(t\right)\end{array}\right.$

式中，m(t)是弯曲运动与时间相关项；

是对m(t)时间t的一阶导数；

是m(t)对时间t的二阶导数；n(t)是扭转运动与时间相关项；

是n(t)对时间t的一阶导数；

是n(t)对时间t的二阶导数；ξ_m、ξ_n分别是弯曲和扭转的模态阻尼系数；ω_m、ω_n分别是弯曲和扭转的模态频率；f_m(t)、f_n(t)是白噪声；

(3)由于用于对地观测成像的飞机一般为平直翼飞机，且处于亚音速飞行状态，建立机翼弯曲和扭转运动与位置相关项的数学方程为：

$\left\{\begin{array}{c}p\left(x\right)=\cosh \left(\mathrm{\βx}\right)-\cos \left(\mathrm{\βx}\right)+C_1(\sinh \left(\mathrm{\βx}\right)-\sin \left(\mathrm{\βx}\right))\\ q\left(x\right)=C_2\sin \frac{\mathrm{\πx}}{2L}\end{array}\right.$

式中：p(x)是机翼弯曲运动与空间位置相关项，q(x)是机翼扭转运动与空间位置相关项，自变量x为机翼上的坐标点，L为机翼长度，C₁，C₂为比例常数，

是机翼结构特性固有参数，ω是机翼自由弯曲振动时的振动圆频率，

是与机翼刚度和线质量相关的固有参数，EI是机翼的抗弯刚度，m代表机翼的线质量；

运用麦夸特曲线拟合的方法拟合得到p(x)，q(x)，即确定上式中的系数C₁，C₂，β；

(4)结合步骤(2)得到的m(t)，n(t)以及步骤(3)得到的p(x)，q(x)，可建立载机机翼弹性运动方程为：

$\left\{\begin{array}{c}{\mathrm{\γ}}_c(x,t)=p\left(x\right)m\left(t\right)\\ {\mathrm{\γ}}_s(x,t)=q\left(x\right)n\left(t\right)\end{array}\right.$

式中，γ_c(x，t)是机翼沿z方向的弯曲位移，γ_s(x，t)是机翼绕ox轴的扭转角位移，γ_c(x，t)和γ_s(x，t)均是随时间变化的；

(5)建立机翼挠曲变形角μ(x，t)＝[μ_x(x，t) μ_y(x，t) μ_z(x，t)]^T与扭转和弯曲位移关系为：

$\left\{\begin{array}{c}{\mathrm{\μ}}_x(x,t)={\mathrm{\γ}}_s(x,t)\\ {\mathrm{\μ}}_y(x,t)=\frac{{\∂\mathrm{\γ}}_c(x,t)}{\∂x}\\ {\mathrm{\μ}}_z(x,t)=0\end{array}\right.$

式中，μ_x(x，t)是机翼沿x方向的挠曲变形角，μ_y(x，t)是机翼沿y方向的挠曲变形角，μ_z(x，t)是机翼沿z方向的挠曲变形角。μ_x(x，t)，μ_y(x，t)，μ_z(x，t)均是随时间变化的。μ(x，t)的具体表达式为：

$\left\{\begin{array}{c}{\mathrm{\μ}}_x(x,t)=C_2\sin \frac{\mathrm{\πx}}{2L}n\left(t\right)\\ {\mathrm{\μ}}_y(x,t)=p^{\′}\left(x\right)m\left(t\right)=(\frac{\mathrm{\β}}{2}\sinh \left(\mathrm{\βx}\right)+\mathrm{\β}\sin \left(\mathrm{\βx}\right)+C_1(\frac{\mathrm{\β}}{2}\cosh \left(\mathrm{\βx}\right)-\mathrm{\β}\cos \left(\mathrm{\βx}\right)))m\left(t\right)\\ {\mathrm{\μ}}_z(x,t)=0\end{array}\right.$

机翼挠曲变形角速度 $\stackrel{\·}{\mathrm{\μ}}(x,t)={\left[\begin{array}{ccc}{\stackrel{\·}{\mathrm{\μ}}}_x(x,t)& {\stackrel{\·}{\mathrm{\μ}}}_y(x,t)& {\stackrel{\·}{\mathrm{\μ}}}_z(x,t)\end{array}\right]}^T$ 为：

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\μ}}}_x(x,t)=C_2\sin \frac{\mathrm{\πx}}{2L}\stackrel{\·}{n}\left(t\right)\\ {\stackrel{\·}{\mathrm{\μ}}}_y(x,t)=p^{\′}\left(x\right)m\left(t\right)=(\frac{\mathrm{\β}}{2}\sinh \left(\mathrm{\βx}\right)+\mathrm{\β}\sin \left(\mathrm{\βx}\right)+C_1(\frac{\mathrm{\β}}{2}\cosh \left(\mathrm{\βx}\right)-\mathrm{\β}\cos \left(\mathrm{\βx}\right)))\stackrel{\·}{m}\left(t\right)\\ {\stackrel{\·}{\mathrm{\μ}}}_z(x,t)=0\end{array}\right.$

式中，

是机翼沿x方向的挠曲变形角速度，

是机翼沿y方向的挠曲变形角速度，

机翼沿z方向的挠曲变形角速度。

均是随时间变化的。

本发明的原理是：如图4所示，通过对机翼挠曲变形运动进行动力学建模，将由机翼挠曲运动引起的挠曲变形角和挠曲变形角速度增广为卡尔曼滤波的状态变量。利用卡尔曼滤波估计出主、子惯性测量单元间的固定安装误差角λ和挠曲变形角μ。利用主惯性测量单元与GPS进行信息融合后得到的基准信息，对子惯性测量单元速度和姿态进行校正，最终获得每个子惯性测量单元安装点的高精度速度、姿态信息，完成传递对准。

如图1所示，在对机翼挠曲变形运动进行建模时，充分考虑了机翼挠曲运动受时间因素和空间因素两方面影响。时间因素是指在不同时刻，机翼的挠曲变形程度不同。因此在建模时建立了机翼挠曲运动(包括弯曲和扭转)与时间相关项的数学方程；空间因素是指机翼不同位置的挠曲变形程度不同。比如，靠近机腹处的挠曲变形程度较低，而靠近机翼尖端处的挠曲变形程度较高。因此在建模时建立了机翼挠曲运动(包括弯曲和扭转)与位置相关项的数学方程。最后，将两方面因素综合到一起，建立了机翼挠曲变形运动整体的数学模型。数学模型中包括了机翼挠曲运动引起的挠曲变形角和挠曲变形角速度的模型，基于该模型可建立相应的卡尔曼滤波方程。通过卡尔曼滤波可估计出挠曲变形角。在此基础上，利用基准信息完成姿态校正，从而完成传递对准。

本发明与现有技术相比的优点在于：充分考虑了机翼挠曲变形对传递对准精度的影响，并通过将挠曲变形角及变形角速度增广为卡尔曼滤波状态变量的方法对其进行补偿。在建立卡尔曼滤波模型时，利用力学的方法对载机机翼挠曲变形进行建模，是在充分考虑载机机翼的结构特性的基础上，从力学的角度精确的建立机翼挠曲变形运动的数学模型，模型中考虑了机翼不同位置处挠曲变形的不同情况。在此基础上，利用卡尔曼滤波可精确的估计出安装误差角和挠曲变形角，从而完成子惯性测量单元的姿态校正，实现高精度的传递对准。

附图说明

图1为本发明的流程图；

图2为本发明的机翼挠曲变形运动建模的流程图；

图3卡尔曼滤波算法基本编排图；

图4为本发明的“速度+姿态”匹配传递对准原理图；

具体实施方式

由于分布式POS的主惯性测量单元通常安装在机腹中央，子惯性测量单元安装在机翼两侧的不同载荷附近，在大气湍流的影响下，机翼较机腹更容易发生挠曲变形，因此本发明对机翼的挠曲变形进行建模。

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

1、利用力学的方法建立载机机翼挠曲变形的运动模型，得到挠曲变形角μ和挠曲变形角速度

的运动方程；

(1)利用ANSYS软件对机翼结构进行模态分析，得到机翼的第一阶弯曲模态频率ω_m和扭转模态频率ω_n；

(2)机翼弯曲和扭转运动与时间相关项的数学方程为：

$\left\{\begin{array}{c}\stackrel{\·\·}{m}\left(t\right)+2{\mathrm{\ξ}}_m{\mathrm{\ω}}_m\stackrel{\·}{m}\left(t\right)+{{\mathrm{\ω}}_m}^{2}m\left(t\right)=f_m\left(t\right)\\ \stackrel{\·\·}{n}\left(t\right)+2{\mathrm{\ξ}}_n{\mathrm{\ω}}_n\stackrel{\·}{n}\left(t\right)+{{\mathrm{\ω}}_n}^{2}n\left(t\right)=f_n\left(t\right)\end{array}\right.---\left(1\right)$

式中，m(t)是弯曲运动与时间相关项；

是对m(t)时间t的一阶导数；

是m(t)对时间t的二阶导数；n(t)是扭转运动与时间相关项；

是n(t)对时间t的一阶导数；是n(t)对时间t的二阶导数；ξ_m、ξ_n分别是弯曲和扭转的模态阻尼系数；ω_m、ω_n分别是弯曲和扭转的模态频率；f_m(t)、f_n(t)是白噪声；

(3)由于用于对地观测成像的飞机一般为平直翼飞机，且处于亚音速飞行状态，建立机翼弯曲和扭转运动与位置相关项的数学方程为：

$\left\{\begin{array}{c}p\left(x\right)=\cosh \left(\mathrm{\βx}\right)-\cos \left(\mathrm{\βx}\right)+C_1(\sinh \left(\mathrm{\βx}\right)-\sin \left(\mathrm{\βx}\right))\\ q\left(x\right)=C_2\sin \frac{\mathrm{\πx}}{2L}\end{array}\right.---\left(2\right)$

式中：p(x)是机翼弯曲运动与空间位置相关项，q(x)是机翼扭转运动与空间位置相关项，自变量x为机翼上的坐标点，L为机翼长度，C₁，C₂为比例常数，

是机翼结构特性固有参数，ω是机翼自由弯曲振动时的振动圆频率，是与机翼刚度和线质量相关的固有参数，EI是机翼的抗弯刚度，m代表机翼的线质量。

通过ANSYS对机翼结构的分析，可以得到x-p(x)，x-q(x)曲线图。运用麦夸特曲线拟合的方法拟合得到p(x)，q(x)即确定(2)式中的系数C₁，C₂，β；

(4)结合步骤(2)得到的m(t)和n(t)以及步骤(3)得到的p(x)和q(x)，可建立载机机翼弹性运动方程为：

$\left\{\begin{array}{c}{\mathrm{\γ}}_c(x,t)=p\left(x\right)m\left(t\right)\\ {\mathrm{\γ}}_s(x,t)=q\left(x\right)n\left(t\right)\end{array}\right.---\left(3\right)$

式中，γ_c(x，t)是机翼沿z方向的弯曲位移，γ_s(x，t)是机翼绕ox轴的扭转角位移，γ_c(x，t)和γ_s(x，t)均是随时间变化的；

(5)建立机翼挠曲变形角μ(x，t)＝[μ_x(x，t) μ_y(x，t) μ_z(x，t)]^T与扭转和弯曲位移关系为：

$\left\{\begin{array}{c}{\mathrm{\μ}}_x(x,t)={\mathrm{\γ}}_s(x,t)\\ {\mathrm{\μ}}_y(x,t)=\frac{{\∂\mathrm{\γ}}_c(x,t)}{\∂x}\\ {\mathrm{\μ}}_z(x,t)=0\end{array}\right.---\left(4\right)$

将式(2)，式(3)代入式(4)，得到：

$\left\{\begin{array}{c}{\mathrm{\μ}}_x(x,t)=C_2\sin \frac{\mathrm{\πx}}{2L}n\left(t\right)\\ {\mathrm{\μ}}_y(x,t)=p^{\′}\left(x\right)m\left(t\right)=(\frac{\mathrm{\β}}{2}\sinh \left(\mathrm{\βx}\right)+\mathrm{\β}\sin \left(\mathrm{\βx}\right)+C_1(\frac{\mathrm{\β}}{2}\cosh \left(\mathrm{\βx}\right)-\mathrm{\β}\cos \left(\mathrm{\βx}\right)))m\left(t\right)\\ {\mathrm{\μ}}_z(x,t)=0\end{array}\right.---\left(5\right)$

式中，μ_x(x，t)是机翼沿x方向的挠曲变形角，μ_y(x，y)是机翼沿y方向的挠曲变形角，μ_z(x，t)是机翼沿z方向的挠曲变形角。μ_x(x，t)，μ_y(x，t)，μ_z(x，t)均是随时间变化的。p′(x)是机翼弯曲运动与空间位置相关项p(x)的一阶导数。

由于挠曲变形角速度

即为挠曲变形角μ对时间的导数，因此对(5)式求导有：

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\μ}}}_x(x,t)=C_2\sin \frac{\mathrm{\πx}}{2L}\stackrel{\·}{n}\left(t\right)\\ {\stackrel{\·}{\mathrm{\μ}}}_y(x,t)=p^{\′}\left(x\right)m\left(t\right)=(\frac{\mathrm{\β}}{2}\sinh \left(\mathrm{\βx}\right)+\mathrm{\β}\sin \left(\mathrm{\βx}\right)+C_1(\frac{\mathrm{\β}}{2}\cosh \left(\mathrm{\βx}\right)-\mathrm{\β}\cos \left(\mathrm{\βx}\right)))\stackrel{\·}{m}\left(t\right)\\ {\stackrel{\·}{\mathrm{\μ}}}_z(x,t)=0\end{array}\right.---\left(6\right)$

式中，是机翼沿x方向的挠曲变形角速度，

是机翼沿y方向的挠曲变形角速度，

机翼沿z方向的挠曲变形角速度。

均是随时间变化的。

2、采用卡尔曼滤波进行主惯性测量单元与GPS的信息融合，获得主惯性测量单元安装点的速度和姿态基准信息，具体公式如下：

(1)建立系统状态方程：

$\stackrel{\·}{X}=\mathrm{FX}+\mathrm{GW}---\left(7\right)$

式中，X为系统状态矢量，W为系统噪声矢量，F为系统转移矩阵，G为噪声转换矩阵：

$W={\left[\begin{array}{cccccc}{w}_{{\mathrm{\ϵ}}_x}& w_{{\mathrm{\ϵ}}_y}& w_{{\mathrm{\ϵ}}_z}& w_{{\▿}_x}& w_{{\▿}_y}& w_{{\▿}_z}\end{array}\right]}^T$

$F=\left[\begin{array}{cc}{F}_{\mathrm{INS}}& F_S\\ o_{6\×6}& F_M\end{array}\right]$ $F_S=\left[\begin{array}{cc}{C}_b^n& 0_{3\×3}\\ 0_{3\×3}& C_b^n\\ 0_{3\×3}& 0_{3\×3}\end{array}\right]$ F_M＝[0_6×15]

$G=\left[\begin{array}{cc}{C}_b^n& 0_{3\×3}\\ 0_{3\×3}& C_b^n\\ 0_{9\×3}& 0_{9\×3}\end{array}\right]$

其中，φ_x、φ_y和φ_z为平台失准角，δv_x、δv_y和δv_z为速度误差，δL、δλ和δh分别为纬度误差、经度误差和高度误差，ε_x、ε_y和ε_z为陀螺仪随机常值漂移误差，

和

为加速度计随机常值偏置误差，F_INS为惯性系统矩阵，

为系统姿态转换矩阵；

(2)建立系统的量测方程：

Z＝HX+η    (8)

其中，Z为观测矢量，H为观测矩阵，η为量测噪声，I为单位矩阵：

Z＝[δ_L δ_λ δ_h δV_E δV_N δV_U]^T

$H=\left[\begin{array}{ccc}{0}_{3\×6}& I_{3\×3}& 0_{3\×6}\\ 0_{3\×3}& I_{3\×3}& 0_{3\×9}\end{array}\right]$

$\mathrm{\η}={\left[\begin{array}{cccccc}{\mathrm{\η}}_L& {\mathrm{\η}}_{\mathrm{\λ}}& {\mathrm{\η}}_h& {\mathrm{\η}}_{V_E}& {\mathrm{\η}}_{V_N}& {\mathrm{\η}}_{V_U}\end{array}\right]}^T$

(3)卡尔曼滤波基本算法编排，该算法的流程图如图3所示。

状态一步预测方程：

${\stackrel{\mathrm{\Λ}}{X}}_{k/k-1}={\mathrm{\φ}}_{k,k-1}{\stackrel{\mathrm{\Λ}}{X}}_{k-1}---\left(9\right)$

其中，

φ_k，k-1分别为k时刻系统状态一步预测值、k-1时刻系统状态估值、k-1时刻到k时刻的系统状态转移矩阵；

状态估值计算方程；

${\stackrel{\mathrm{\Λ}}{X}}_k={\stackrel{\mathrm{\Λ}}{X}}_{k/k-1}+K_k(Z_k-H_k{\stackrel{\mathrm{\Λ}}{X}}_{k/k-1})---\left(10\right)$

其中，

K_k、Z_k、H_k分别为k时刻系统状态估值、系统增益矩阵、量测向量和量测矩阵；

滤波增量方程：

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

其中，P_k/k-1、R_k分别为k时刻系统协方差阵的一步预测、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---\left(12\right)$

其中，P_k-1、Q_k-1、Г_k-1分别为k-1时刻系统协方差阵、k-1时刻系统噪声矩阵、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^T---\left(13\right)$

其中，P_k为k时刻系统状态协方差阵，I为单位阵；

3、采用“速度+姿态”的匹配方式建立卡尔曼滤波模型，在考虑机翼挠曲变形的情况下，结合步骤1的建模，将机翼因挠曲变形运动产生的挠曲变形角μ和挠曲变形角速度增广为卡尔曼滤波的状态变量，即在确定出机翼挠曲变形角和挠曲变形角速度的运动模型后，将x，y方向的挠曲变形角μ_x，μ_y和将x，y方向的挠曲变形角速度

增广为卡尔曼滤波的状态变量(z方向的挠曲变形角μ_z和挠曲变形角速度

均为零)；结合步骤2得到的主惯性测量单元安装点的速度和姿态的基准信息，通过卡尔曼滤波估计出主、子惯性测量单元间的姿态失准角；其中卡尔曼滤波器的模型包括状态方程和量测方程，分别如式(14)和式(15)所示。

(1)建立系统状态方程：

$\stackrel{\·}{X}=\mathrm{FX}+\mathrm{GW}---\left(14\right)$

式中：X为系统状态矢量，W为系统噪声矢量，F为系统转移矩阵，G为噪声转换矩阵：

$X=[{\mathrm{\φ}}_E,{\mathrm{\φ}}_N,{\mathrm{\φ}}_U,\mathrm{\δ}{V}_E,\mathrm{\δ}{V}_N,\mathrm{\δ}{V}_U,{\mathrm{\ϵ}}_x,{\mathrm{\ϵ}}_y,{\mathrm{\ϵ}}_z,{\▿}_x,{\▿}_y,{\▿}_z,$

${{\mathrm{\λ}}_x,{\mathrm{\λ}}_y,{\mathrm{\λ}}_z,{\mathrm{\μ}}_x,{\mathrm{\μ}}_y,{\stackrel{\·}{\mathrm{\μ}}}_x,{\stackrel{\·}{\mathrm{\μ}}}_y]}^T$

其中，φ_E是子惯性测量单元东向平台失准角，φ_N是子惯性测量单元东向平台失准角，φ_U是子惯性测量单元天向平台失准角；δV_E是主、子惯性测量单元间的东向速度误差，δV_N是主、子惯性测量单元间的北向速度误差，δV_U是主、子惯性测量单元间的天向速度误差；ε_x是子惯性测量单元载体坐标系x轴陀螺常值漂移，ε_y是子惯性测量单元载体坐标系y轴陀螺常值漂移，ε_z是子惯性测量单元载体坐标系z轴陀螺常值漂移；

是子惯性测量单元载体坐标系x轴加计常值偏置，

是子惯性测量单元载体坐标系y轴加计常值偏置，

是子惯性测量单元载体坐标系z轴加计常值偏置；λ_x是子惯性测量单元载体坐标系x轴固定安装误差角，λ_y是子惯性测量单元载体坐标系y轴固定安装误差角，λ_z是子惯性测量单元载体坐标系z轴固定安装误差角；μ_x是子惯性测量单元载体坐标系x轴挠曲变形角，μ_y是子惯性测量单元载体坐标系y轴挠曲变形角；

是子惯性测量单元载体坐标系x轴挠曲变形角速度，

是子惯性测量单元载体坐标系y轴挠曲变形角速度。

$W={\left[\begin{array}{cccccccc}{w}_{\mathrm{\ϵx}}& w_{\mathrm{\ϵy}}& w_{\mathrm{\ϵz}}& w_{\▿x}& w_{\▿y}& w_{\▿z}& f_m\left(t\right)& f_n\left(t\right)\end{array}\right]}^T$

w_εx是子惯性测量单元载体坐标系x轴陀螺随机误差，w_εy是子惯性测量单元载体坐标系y轴陀螺随机误差，w_εz是子惯性测量单元载体坐标系z轴陀螺随机误差；

是子惯性测量单元载体坐标系x轴加计随机误差，是子惯性测量单元载体坐标系y轴加计随机误差，是子惯性测量单元载体坐标系z轴加计随机误差；f_m(t)是机翼弯曲运动白噪声，f_n(t)是机翼扭转运动白噪声。

$F=\left[\begin{array}{ccccccccc}-[{\mathrm{\ω}}_{\mathrm{in}}^n\×]& F_{12}& C_s^n& 0_{3\×3}& 0_{3\×3}& 0_{3\×1}& 0_{3\×1}& 0_{3\×1}& 0_{3\×1}\\ [f^n\×]& F_{22}& 0_{3\×3}& C_s^n& 0_{3\×3}& 0_{3\×1}& 0_{3\×1}& 0_{3\×1}& 0_{3\×1}\\ 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×1}& 0_{3\×1}& 0_{3\×1}& 0_{3\×1}\\ 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×1}& 0_{3\×1}& 0_{3\×1}& 0_{3\×1}\\ 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×1}& 0_{3\×1}& 0_{3\×1}& 0_{3\×1}\\ 0_{n\×3}& 0_{n\×3}& 0_{n\×3}& 0_{n\×3}& 0_{n\×3}& 0& 0& 1& 0\\ 0_{n\×3}& 0_{n\×3}& 0_{n\×3}& 0_{n\×3}& 0_{n\×3}& 0& 0& 0& 1\\ 0_{n\×3}& 0_{n\×3}& 0_{n\×3}& 0_{n\×3}& 0_{n\×3}& -{\mathrm{\ω}}_n^2& 0& -2{\mathrm{\ξ}}_n{\mathrm{\ω}}_n& 0\\ 0_{n\×3}& 0_{n\×3}& 0_{n\×3}& 0_{n\×3}& 0_{n\×3}& 0& -{\mathrm{\ω}}_m^2& 0& -2{\mathrm{\ξ}}_m{\mathrm{\ω}}_m\end{array}\right]$

其中，

是主惯性测量单元的导航坐标系相对惯性坐标系的旋转角速度在导航坐标系上的投影分量组成的反对称阵；[fⁿ×]是子惯性测量单元加计测得的比例在导航坐标系上的投影分量组成的反对称阵；是子惯性测量单元载体坐标系到导航坐标系的转换矩阵；ξ_m、ξ_n分别是弯曲和扭转的模态阻尼系数；ω_m、ω_n分别是弯曲和扭转的模态频率。

$F_{12}=\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\end{array}\right]$

其中，R_M是地球子午圈主曲率半径，R_N是地球卯酉圈主曲率半径，H是高度。

$F_{22}=\left[\begin{array}{ccc}\frac{V_N\tan L}{R_N+H}-\frac{V_U}{R_N+H}& 2{\mathrm{\ω}}_{\mathrm{ie}}\sin L+\frac{V_E\tan L}{R_N+H}& -2{\mathrm{\ω}}_{\mathrm{ie}}\cos L-\frac{V_E}{R_N+H}\\ -2{\mathrm{\ω}}_{\mathrm{ie}}\sin L-\frac{{2V}_E\tan L}{R_N+H}& -\frac{V_U}{R_M+H}& -\frac{V_N}{R_M+H}\\ 2({\mathrm{\ω}}_{\mathrm{ie}}\cos L+\frac{V_E}{R_N+H})& \frac{2V_N}{R_M+H}& 0\end{array}\right]$

其中，V_E，V_N，V_U分别是惯性测量单元与GPS信息进行融合后东、北、天方向的基准速度，L是基准纬度，ω_ie是地球自转角速度。

$G^w=\left[\begin{array}{cccc}{C}_s^n& 0_{3\×3}& 0_{3\×1}& 0_{3\×1}\\ 0_{3\×3}& C_s^n& 0_{3\×1}& 0_{3\×1}\\ 0_{11\×3}& 0_{11\×3}& 0_{11\×1}& 0_{11\×1}\\ 0_{1\×3}& 0_{1\×3}& 1& 0\\ 0_{1\×3}& 0_{1\×3}& 0& 1\end{array}\right]$

(2)建立系统量测方程：

当采用“速度+姿态”匹配方法，即量测量选为主、子惯性测量单元的速度误差和姿态误差进行分布式POS的传递对准时，量测值为主、子惯性测量单元的姿态角误差和速度误差。系统量测方程为：

Z＝HX+η    (15)

式中：Z为量测变量，H为量测矩阵，η是量测噪声，假定其为零均值的高斯白噪声；

Z＝[δψ δθ δγ δV_E δV_N δV_U]^T

其中，δψ、δθ、δγ分别是是主、子惯性测量单元的航向角误差、俯仰角误差和横滚角误差，即三个方向的姿态误差。δV_E、δV_N、δV_U分别是是主、子惯性测量单元的东向、北向、天向速度误差，即三个方向的速度误差。

$H=\left[\begin{array}{ccccccccc}{H}_{11}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& H_{15}& H_{16}& H_{17}& 0_{3\×1}& 0_{3\×1}\\ 0_{3\×3}& I_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×3}& 0_{3\×1}& 0_{3\×1}& 0_{3\×1}& 0_{3\×1}\end{array}\right]$

记为主惯性测量单元载体坐标系到导航坐标系的转换矩阵，令

$T_a=\left[\begin{array}{ccc}{T}_a^{\left(11\right)}& T_a^{\left(12\right)}& T_a^{\left(13\right)}\\ T_a^{\left(21\right)}& T_a^{\left(22\right)}& T_a^{\left(23\right)}\\ T_a^{\left(31\right)}& T_a^{\left(32\right)}& T_a^{\left(33\right)}\end{array}\right]$

$H_{11}=\left[\begin{array}{ccc}\frac{T_a^{12}{T}_a^{32}}{{\left(T_a^{12}\right)}^2+{\left(T_a^{22}\right)}^2}& \frac{T_a^{22}{T}_a^{32}}{{\left(T_a^{12}\right)}^2+{\left(T_a^{22}\right)}^2}& -1\\ -\frac{T_a^{22}}{\sqrt{1-{\left(T_a^{32}\right)}^2}}& \frac{T_a^{12}}{\sqrt{1-{\left(T_a^{32}\right)}^2}}& 0\\ \frac{T_a^{21}{T}_a^{33}-T_a^{31}{T}_a^{23}}{{\left(T_a^{33}\right)}^2+{\left(T_a^{31}\right)}^2}& \frac{T_a^{31}{T}_a^{13}-T_a^{11}{T}_a^{33}}{{\left(T_a^{33}\right)}^2+{\left(T_a^{31}\right)}^2}& 0\end{array}\right]$

$H_{15}=\left[\begin{array}{ccc}\frac{T_a^{12}{T}_a^{23}-T_a^{13}{T}_a^{22}}{{\left(T_a^{12}\right)}^2+{\left(T_a^{22}\right)}^2}& 0& \frac{T_a^{11}{T}_a^{22}-T_a^{12}{T}_a^{21}}{{\left(T_a^{12}\right)}^2+{\left(T_a^{22}\right)}^2}\\ \frac{T_a^{33}}{\sqrt{1-{\left(T_a^{32}\right)}^2}}& 0& -\frac{T_a^{31}}{\sqrt{1-{\left(T_a^{32}\right)}^2}}\\ \frac{-T_a^{31}{T}_a^{32}}{{\left(T_a^{33}\right)}^2+{\left(T_a^{31}\right)}^2}& 1& \frac{-T_a^{32}{T}_a^{33}}{{\left(T_a^{33}\right)}^2+{\left(T_a^{31}\right)}^2}\end{array}\right]$

$H_{16}=\left[\begin{array}{c}0\\ 0\\ p^{\′}\left(x\right)\end{array}\right]$ $H_{17}=\left[\begin{array}{c}\frac{(T_a^{12}{T}_a^{23}-T_a^{13}{T}_a^{22})q\left(x\right)}{{\left(T_a^{12}\right)}^2+{\left(T_a^{22}\right)}^2}\\ \frac{T_a^{33}q\left(x\right)}{\sqrt{1-{\left(T_a^{32}\right)}^2}}\\ \frac{-T_a^{31}{T}_a^{32}q\left(x\right)}{{\left(T_a^{33}\right)}^2+{\left(T_a^{31}\right)}^2}\end{array}\right]$

其中，p′(x)是机翼弯曲运动与空间位置相关项p(x)的一阶导数。q(x)是机翼扭转运动与空间位置相关项。

(3)卡尔曼滤波基本算法编排，根据上述系统状态方程和量测方程，可建立离散卡尔曼滤波方程如下：

状态一步预测方程：

${\stackrel{\mathrm{\Λ}}{X}}_{k/k-1}={\mathrm{\φ}}_{k,k-1}{\stackrel{\mathrm{\Λ}}{X}}_{k-1}---\left(16\right)$

其中，

φ_k，k-1分别为k时刻系统状态一步预测值、k-1时刻系统状态估值、k-1时刻到k时刻的系统状态转移矩阵；

状态估值计算方程：

${\stackrel{\mathrm{\Λ}}{X}}_k={\stackrel{\mathrm{\Λ}}{X}}_{k/k-1}+K_k(Z_k-H_k{\stackrel{\mathrm{\Λ}}{X}}_{k/k-1})---\left(17\right)$

其中，K_k、Z_k、H_k分别为k时刻系统状态估值、系统增益矩阵、量测向量和量测矩阵；

滤波增量方程：

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

其中，P_k/k-1、R_k分别为k时刻系统协方差阵的一步预测、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---\left(19\right)$

其中，P_k-1、Q_k-1、Г_k-1分别为k-1时刻系统协方差阵、k-1时刻系统噪声矩阵、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^T---\left(20\right)$

其中，P_k为k时刻系统状态协方差阵，I是单位阵；

4、利用步骤2得到的基准信息和步骤3估计出的主、子惯性测量单元间的姿态失准角，对每个子惯性测量单元计算出速度、姿态进行校正，最终获得每个子惯性测量单元安装点的经校正后的速度和姿态信息。

(1)速度校正：

速度校正是直接利用子惯性测量单元量测得的速度减去在步骤3中利用卡尔曼滤波估计出的主、子惯性测量单元间的速度误差，即可得到子惯性测量单元被校正后的速度信息。

V_{s_new}＝V_{s_old}-δV    (21)

其中，V_{s_new}是子惯性测量单元被更新后的速度，V_{s_old}是子惯性测量单元测得到的速度，δV是用卡尔曼滤波估计出的主、子惯性测量单元间的速度误差。

(2)姿态校正：

姿态校正是指对每个子惯性测量单元计算得到的姿态阵

进行校正。再由解算出相应的姿态角。

利用估计出的固定安装误差角λ和挠曲变形角μ，计算主、子惯性测量单元的相对姿态矩阵为：

$C_s^m=I_{3\×3}+(\mathrm{\λ}+\mathrm{\μ})\×---\left(22\right)$

式中：(λ+μ)×为子系统相对主系统的失准角组成的反对称阵。

经展开计算后得到：

$C_s^m=\left[\begin{array}{ccc}1& -({\mathrm{\λ}}_z+{\mathrm{\μ}}_z)& {\mathrm{\λ}}_y+{\mathrm{\μ}}_y\\ {\mathrm{\λ}}_z+{\mathrm{\μ}}_z& 1& -({\mathrm{\λ}}_x+{\mathrm{\μ}}_x)\\ -({\mathrm{\λ}}_y+{\mathrm{\μ}}_y)& {\mathrm{\λ}}_x+{\mathrm{\μ}}_x& 1\end{array}\right]$

并结合主系统的姿态阵

更新子惯性测量单元的姿态阵

为：

$C_s^n=C_m^n{C}_s^m---\left(23\right)$

由被更新后的子惯性测量单元的姿态阵

计算子惯性测量单元安装点的航向角ψ_s，俯仰角θ_s和横滚角γ_s。将(23)式计算出的

记为

$C_s^n=\left[\begin{array}{ccc}{T}_{11}& T_{12}& T_{13}\\ T_{21}& T_{22}& T_{23}\\ T_{31}& T_{32}& T_{33}\end{array}\right]---\left(24\right)$

又因为

因此，由(23)式和(24)式，可以确定出航向角

俯仰角θ_s和横滚角γ_s的主值，即

θ_s主＝arcsin(T₃₂)    (26)

若航向角

俯仰角θ_s和横滚角γ_s的取值范围分别定义为[0，2π]、

[-π，+π]。那么，

θ_s和γ_s的真值可由下式确定：

θ_s＝θ_s主    (28)

通过对子惯性测量单元的速度、姿态进行校正，能够得到子惯性测量单元安装点的高精度速度、姿态信息，完成传递对准。

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

Claims (2)

1.一种机载分布式POS的传递对准方法，具体步骤为：

(1)利用力学的方法建立载机机翼挠曲变形的运动模型，得到挠曲变形角μ和挠曲变形角速度

的运动方程；

(2)采用卡尔曼滤波技术进行主惯性测量单元与GPS的信息融合，获得主惯性测量单元安装点的速度和姿态基准信息；

(3)采用“速度+姿态”匹配方式建立卡尔曼滤波模型，将机翼因挠曲变形运动产生的挠曲变形角μ和挠曲变形角速度

增广为卡尔曼滤波的状态变量，通过卡尔曼滤波估计出主、子惯性测量单元间的姿态失准角；

(4)利用步骤(2)得到的基准信息和步骤(3)估计出的主、子惯性测量单元间的姿态失准角，对每个子惯性测量单元计算出速度、姿态进行校正，最终获得每个子惯性测量单元安装点的经校正后的速度和姿态信息。

2.根据权利要求1所述的一种机载分布式POS的传递对准方法，所述的步骤(1)利用基于力学的方法建立载机机翼挠曲变形运动模型，其具体步骤为：

(1.1)利用有限元软件对机翼结构进行模态分析，得到机翼的第一阶弯曲模态频率ω_m和扭转模态频率ω_n；

(1.2)建立机翼弯曲和扭转运动与时间相关项的数学方程为：

$\left\{\begin{array}{c}\stackrel{\·\·}{m}\left(t\right)+2{\mathrm{\ξ}}_m{\mathrm{\ω}}_m\stackrel{\·}{m}\left(t\right)+{{\mathrm{\ω}}_m}^{2}m\left(t\right)=f_m\left(t\right)\\ \stackrel{\·\·}{n}\left(t\right)+2{\mathrm{\ξ}}_n{\mathrm{\ω}}_n\stackrel{\·}{n}\left(t\right)+{{\mathrm{\ω}}_n}^{2}n\left(t\right)=f_n\left(t\right)\end{array}\right.---\left(1\right)$

式中，m(t)是弯曲运动与时间相关项；

是对m(t)时间t的一阶导数；

是m(t)对时间t的二阶导数；n(t)是扭转运动与时间相关项；

是n(t)对时间t的一阶导数；

是n(t)对时间t的二阶导数；ξ_m、ξ_n分别是弯曲和扭转的模态阻尼系数；ω_m、ω_n分别是弯曲和扭转的模态频率；f_m(t)、f_n(t)是白噪声；

(1.3)机翼弯曲和扭转运动与位置坐标相关项的数学方程为：

$\left\{\begin{array}{c}p\left(x\right)=\cosh \left(\mathrm{\βx}\right)-\cos \left(\mathrm{\βx}\right)+C_1(\sinh \left(\mathrm{\βx}\right)-\sin \left(\mathrm{\βx}\right))\\ q\left(x\right)=C_2\sin \frac{\mathrm{\πx}}{2L}\end{array}\right.---\left(2\right)$

式中：p(x)是机翼弯曲运动与空间位置相关项，q(x)是机翼扭转运动与空间位置相关项，自变量x为机翼上的坐标点，L为机翼长度，C₁，C₂为比例常数，

是机翼结构特性固有参数，ω是机翼自由弯曲振动时的振动圆频率，

是与机翼刚度和线质量相关的固有参数，EI是机翼的抗弯刚度，m代表机翼的线质量；

运用麦夸特曲线拟合的方法拟合得到p(x)，q(x)即确定(2)式中的系数C₁，C₂，β；

(1.4)结合步骤(2)得到的m(t)，n(t)以及步骤(3)得到的p(x)，q(x)，建立载机机翼完整弹性运动方程为：

$\left\{\begin{array}{c}{\mathrm{\γ}}_c(x,t)=p\left(x\right)m\left(t\right)\\ {\mathrm{\γ}}_s(x,t)=q\left(x\right)n\left(t\right)\end{array}\right.---\left(3\right)$

式中，γ_c(x，t)是机翼沿z方向的弯曲位移，γ_s(x，t)是机翼绕ox轴的扭转角位移，γ_c(x，t)和γ_s(x，t)均是随时间变化的；

(1.5)建立机翼挠曲变形角μ(x，t)＝[μ_x(x，t) μ_y(x，t) μ_z(x，t)]^T与扭转和弯曲位移关系为：

$\left\{\begin{array}{c}{\mathrm{\μ}}_x(x,t)={\mathrm{\γ}}_s(x,t)\\ {\mathrm{\μ}}_y(x,t)=\frac{{\∂\mathrm{\γ}}_c(x,t)}{\∂x}\\ {\mathrm{\μ}}_z(x,t)=0\end{array}\right.---\left(4\right)$

式中，μ_x(x，t)是机翼沿x方向的挠曲变形角，μ_y(x，t)是机翼沿y方向的挠曲变形角，μ_z(x，t)是机翼沿z方向的挠曲变形角。μ_x(x，t)，μ_y(x，t)，μ_z(x，t)均是随时间变化的。μ(x，t)的具体表达式为：

$\left\{\begin{array}{c}{\mathrm{\μ}}_x(x,t)=C_2\sin \frac{\mathrm{\πx}}{2L}n\left(t\right)\\ {\mathrm{\μ}}_y(x,t)=p^{\′}\left(x\right)m\left(t\right)=(\frac{\mathrm{\β}}{2}\sinh \left(\mathrm{\βx}\right)+\mathrm{\β}\sin \left(\mathrm{\βx}\right)+C_1(\frac{\mathrm{\β}}{2}\cosh \left(\mathrm{\βx}\right)-\mathrm{\β}\cos \left(\mathrm{\βx}\right)))m\left(t\right)\\ {\mathrm{\μ}}_z(x,t)=0\end{array}\right.---\left(5\right)$

机翼挠曲变形角速度 $\stackrel{\·}{\mathrm{\μ}}(x,t)={\left[\begin{array}{ccc}{\stackrel{\·}{\mathrm{\μ}}}_x(x,t)& {\stackrel{\·}{\mathrm{\μ}}}_y(x,t)& {\stackrel{\·}{\mathrm{\μ}}}_z(x,t)\end{array}\right]}^T$ 为：

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\μ}}}_x(x,t)=C_2\sin \frac{\mathrm{\πx}}{2L}\stackrel{\·}{n}\left(t\right)\\ {\stackrel{\·}{\mathrm{\μ}}}_y(x,t)=p^{\′}\left(x\right)m\left(t\right)=(\frac{\mathrm{\β}}{2}\sinh \left(\mathrm{\βx}\right)+\mathrm{\β}\sin \left(\mathrm{\βx}\right)+C_1(\frac{\mathrm{\β}}{2}\cosh \left(\mathrm{\βx}\right)-\mathrm{\β}\cos \left(\mathrm{\βx}\right)))\stackrel{\·}{m}\left(t\right)\\ {\stackrel{\·}{\mathrm{\μ}}}_z(x,t)=0\end{array}\right.---\left(6\right)$

式中，是机翼沿x方向的挠曲变形角速度，

是机翼沿y方向的挠曲变形角速度，

机翼沿z方向的挠曲变形角速度。

均随时间变化。

Images