一种刚体空间运动状态的任意步长正交级数输出方法

摘要

本发明公开了一种刚体空间运动状态的任意步长正交级数输出方法，该方法通过定义三元数，使得机体轴系三个速度分量和三元数构成线性微分方程组，并采用任意步长正交多项式对滚转、俯仰、偏航角速度p，q，r进行近似逼近描述，可以按照任意阶保持器的方式求解系统的状态转移矩阵，进而得到刚体运动离散状态方程的表达式，避免了姿态方程奇异问题，从而得到刚体主要运动状态；本发明通过引入三元数使得状态转移矩阵为分块上三角形式，可以降阶求解状态转移矩阵，大大简化了计算复杂度，便于工程使用。

说明书

技术领域

本发明涉及空间运动刚体模型，特别涉及飞行器大机动飞行状态输出问题。

背景技术

机体轴系刚体运动微分方程是描述飞行器、鱼雷、航天器等空间运动的基本方程。通常， 在数据处理等应用中，体轴系的状态变量主要包含3个速度分量、三个欧拉角、以及地面坐标 系的X_E，Y_E，Z_E等，由于ZE定义为垂直地面指向地球中心，因此Z_E实际为负的飞行高度； X_E，Y_E通常主要依赖GPS、GNSS、北斗等直接给出；欧拉角表示刚体空间运动姿态，而刻画 刚体姿态的微分方程又是其中的核心，通常以三个欧拉角即俯仰、滚转和偏航角来描述。当 刚体的俯仰角为±90°时，滚转角和偏航角无法定值，同时临近该奇点的区域求解误差过大， 导致工程上不可容忍的误差而不能使用；为了避免这一问题，人们首先采用限制俯仰角取值 范围的方法，这使得方程式退化，不能全姿态工作，因而难以广泛用于工程实践。随着对飞 行器极限飞行的研究，人们又相继采用了方向余弦法、等效转动矢量法、四元数法等推算刚 体运动姿态。

方向余弦法避免了欧拉角描述方法的“奇异”现象，用方向余弦法计算姿态矩阵没有方 程退化问题，可以全姿态工作，但需要求解9个微分方程，计算量较大，实时性较差，无法满 足工程实践要求。等效转动矢量法如单子样递推、双子样转动矢量、三子样转动矢量和四子 样旋转矢量法以及在此基础上的各种修正算法和递推算法等。文献中研究旋转矢量时，都是 基于速率陀螺输出为角增量的算法。然而在实际工程中，一些陀螺的输出是角速率信号，如 光纤陀螺、动力调谐陀螺等。当速率陀螺输出为角速率信号时，旋转矢量法的算法误差明显 增大。四元数法是定义4个欧拉角的函数来计算航姿，能够有效弥补欧拉角描述方法的奇异性， 只要解4个一阶微分方程式组即可，比方向余弦姿态矩阵微分方程式计算量有明显的减少，能 满足工程实践中对实时性的要求。其常用的计算方法有毕卡逼近法、二阶、四阶龙格-库塔法 和三阶泰勒展开法等。毕卡逼近法实质是单子样算法，对有限转动引起的不可交换误差没有 补偿，在高动态情况下姿态解算中的算法漂移会十分严重。采用四阶龙格-库塔法求解四元数 微分方程时，随着积分误差的不断积累，会出现三角函数取值超出±1的现象，从而导致计算 发散；泰勒展开法也因计算精度的不足而受到制约。当刚体大机动时，角速率较大导致上述 方法的误差更大；不仅如此，姿态估计的误差常常会导致速度4个分量、高度输出的误差急剧 增大。

发明内容

为了克服现有刚体运动模型输出误差大的问题，本发明提供一种刚体空间运动状态的任 意步长正交级数输出方法，该方法通过定义三元数，使得机体轴系三个速度分量和三元数构 成线性微分方程组，并采用任意步长正交多项式对滚转、俯仰、偏航角速度p，q，r进行近似 逼近描述，可以按照任意阶保持器的方式求解系统的状态转移矩阵，进而得到刚体运动离散 状态方程的表达式，避免了姿态方程奇异问题，从而得到刚体主要运动状态。

本发明解决其技术问题采用的技术方案是，一种刚体空间运动状态的任意步长正交级数 输出方法，其特征包括以下步骤：

1、机体轴系三个速度分量输出为：

${\left(\begin{array}{c}u\left(t\right)\\ v\left(t\right)\\ w\left(t\right)\end{array}\right)}_{t=(k+1)T}={\Phi }_v[(k+1)T,\mathrm{kT}]{\left(\begin{array}{c}u\left(t\right)\\ v\left(t\right)\\ w\left(t\right)\end{array}\right)}_{t=\mathrm{kT}}+g{\Phi }_v[(k+1)T,\mathrm{kT}]{\Phi }_s[(k+1)T,\mathrm{kT}]{\left(\begin{array}{c}{s}_1\left(t\right)\\ s_2\left(t\right)\\ s_3\left(t\right)\end{array}\right)}_{t=\mathrm{kT}}$

$+g{\int }_{\mathrm{kT}}^{(k+1)T}{\Phi }_v[(k+1)T,\tau ]\left(\begin{array}{c}{n}_x\\ n_y\\ n_z\end{array}\right)\mathrm{d\tau }$

其中：u，v，w分别为沿刚体机体轴系x，y，z轴的速度分量，n_x，n_y，n_z分别为沿x，y，z轴的过载， g为重力加速度，s₁、s₂、s₃为定义的三元数，且

${\left(\begin{array}{c}{s}_1\left(t\right)\\ s_2\left(t\right)\\ s_3\left(t\right)\end{array}\right)}_{t=(k+1)T}={\Phi }_s[(k+1)T,\mathrm{kT}]{\left(\begin{array}{c}{s}_1\left(t\right)\\ s_2\left(t\right)\\ s_3\left(t\right)\end{array}\right)}_{t=\mathrm{kT}}$

${\Phi }_v[(k+1)T,\mathrm{kT}]\approx I+{\Pi }_v\zeta +{\Pi }_v\Omega {\Pi }_v^T$

${\Phi }_s[(k+1)T,\mathrm{kT}]\approx I+{\Pi }_s\zeta +{\Pi }_s\Omega {\Pi }_s^T$

p，q，r分别为滚转、俯仰、偏航角速度，T为采样周期；全文参数定义相同；

$I=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right),$ξ(t)＝[ξ₀(t)ξ₁(t)…ξ_n-1(t)ξ_n(t)]^T，

$\left(\begin{array}{c}{\xi }_0\left(t\right)=1\\ {\xi }_1\left(t\right)=\cos \left[a{\cos }^{-1}(1-2t/b)\right]\\ {\xi }_2\left(t\right)={2\xi }_1\left(t\right)·{\xi }_1\left(t\right)-1\\ ·\\ ·\\ ·\\ {\xi }_{i+1}\left(t\right)={2\xi }_1\left(t\right){\xi }_i\left(t\right)-{\xi }_{i-1}\left(t\right)\end{array}\right)i=\mathrm{2,3},···,n-\mathrm{1,0}\le t\le \mathrm{NT},b=\mathrm{NT}$

为基于Chebyshev(切比雪夫)正交多项式得到的任意步长正交多项式的递推形式，a为任意 实数，滚转、俯仰、偏航角速度p，q，r的展开式分别为

p(t)＝[p₀ p₁…p_n-1 p_n][ξ₀(t)ξ₁(t)…ξ_n-1(t)ξ_n(t)]^T

q(t)＝[q₀ q₁…q_n-1 q_n][ξ₀(t)ξ₁(t)…ξ_n-1(t)ξ_n(t)]^T

r(t)＝[r₀ r₁…r_n-1 r_n][ξ₀(t)ξ₁(t)…ξ_n-1(t)ξ_n(t)]^T

${\Pi }_v=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 1\\ 0& -1& \text{0}\end{array}\right)\left(\begin{array}{ccccc}{p}_0& p_1& ···& p_{n-1}& p_n\end{array}\right)$

$+\left(\begin{array}{ccc}0& 0& -1\\ 0& 0& 0\\ 1& 0& 0\end{array}\right)\left(\begin{array}{ccccc}{q}_0& q_1& ···& q_{n-1}& q_n\end{array}\right)+\left(\begin{array}{ccc}0& 1& 0\\ -1& 0& 0\\ 0& 0& 0\end{array}\right)\left(\begin{array}{ccccc}{r}_0& r_1& ···& r_{n-1}& r_n\end{array}\right)$

${\Pi }_s=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 1\\ 0& -1& 0\end{array}\right)\left(\begin{array}{ccccc}{p}_0& p_1& ···& p_{n-1}& p_n\end{array}\right)$

$+\left(\begin{array}{ccc}0& 0& 1\\ 0& 0& 0\\ -1& 0& 0\end{array}\right)\left(\begin{array}{ccccc}{q}_0& q_1& ···& q_{n-1}& q_n\end{array}\right)+\left(\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right)\left(\begin{array}{ccccc}{r}_0& r_1& ···& r_{n-1}& r_n\end{array}\right)$

$\zeta ={\left(\begin{array}{cccc}{\zeta }_0& {\zeta }_1& ···& {\zeta }_n\end{array}\right)}^T={\int }_{\mathrm{kT}}^{(k+1)T}\xi \left(t\right)\mathrm{dt}$

${\zeta }_i={\int }_{\mathrm{kT}}^{(k+1)T}{\xi }_i\left(t\right)\mathrm{dt}={\int }_{\mathrm{kT}}^{(k+1)T}\cos \left[\mathrm{ai}{\cos }^{-1}(1-2t/b)\right]\mathrm{dt}$

$=\frac{b}{4}\left\{\frac{1}{\mathrm{ai}-1}\cos \right[(\mathrm{ai}-1){\cos }^{-1}(1-2t/b)]$

${-\frac{1}{\mathrm{ai}+1}\cos \left[(\mathrm{ai}+1){\cos }^{-1}(1-2t/b)\right]\left\}\right|}_{\mathrm{kT}}^{(k+1)T}$

$\Omega ={\left\{{\Omega }_{\mathrm{ji}}\right\}}_{j=\mathrm{0,1},···,n;i=\mathrm{1,2},···,n}={\int }_{\mathrm{kT}}^{(k+1)T}\xi \left(t\right){\int }_{\mathrm{kT}}^T\xi \left(\tau \right)\mathrm{d\tau dt}$

${\Omega }_{\mathrm{ji}}={\int }_{\mathrm{kT}}^{(k+1)T}{\xi }_j\left(t\right){\int }_{\mathrm{kT}}^t{\xi }_i\left(\tau \right)\mathrm{d\tau dt}$

$=\frac{b}{8}\left\{\frac{1}{\mathrm{ai}-1}{\int }_{\mathrm{kT}}^{(k+1)T}\right\{\cos \left[(\mathrm{aj}-\mathrm{ai}+1){\cos }^{-1}(1-2t/b)\right]+\cos \left[(\mathrm{aj}+\mathrm{ai}-1){\cos }^{-1}(1-2t/b)\right]\}\mathrm{dt}$

$-\frac{1}{\mathrm{ai}+1}{\int }_{\mathrm{kT}}^{(k+1)T}\left\{\cos \right[(\mathrm{aj}-\mathrm{ai}-1){\cos }^{-1}(1-2t/b)]+\cos [(\mathrm{aj}+\mathrm{ai}+1){\cos }^{-1}(1-2t/b)\left]\right\}\mathrm{dt}\};$

$-\frac{b}{4}\left\{\frac{1}{\mathrm{ai}-1}\cos \right[(\mathrm{ai}-1){\cos }^{-1}(1-2\mathrm{kT}/b)]$

$-\frac{1}{\mathrm{ai}+1}\cos \left[(\mathrm{ai}+1){\cos }^{-1}(1-2\mathrm{kT}/b)\right]\left\}{\int }_{\mathrm{kT}}^{(k+1)T}\cos \right[{\mathrm{aj}\cos }^{-1}(1-2t/b)]\mathrm{dt}$

2、高度输出为：

$\stackrel{·}{h}=\left(\begin{array}{ccc}u& v& w\end{array}\right)\left(\begin{array}{c}{s}_1\\ s_2\\ s_3\end{array}\right)$

其中：h为高度；

3、姿态角的输出为：

$\psi \left(t\right)=\psi \left(\mathrm{kT}\right)+{\int }_{\mathrm{kT}}^t\frac{{\mathrm{qs}}_2\left(t\right)+{\mathrm{rs}}_3\left(t\right)}{s_2^2\left(t\right)+s_3^2\left(t\right)}\mathrm{dt}$

其中：θ，ψ分别表示滚转、俯仰、偏航角，$\left(\begin{array}{c}{s}_1\left(t\right)\\ s_2\left(t\right)\\ s_3\left(t\right)\end{array}\right)={\Phi }_s(t,\mathrm{kT}){\left(\begin{array}{c}{s}_1\left(t\right)\\ s_2\left(t\right)\\ s_3\left(t\right)\end{array}\right)}_{t=\mathrm{kT}}.$

本发明的有益效果是：通过引入三元数使得状态转移矩阵为分块上三角形式，可以降阶求 解状态转移矩阵，大大简化了计算复杂度，便于工程使用。

下面结合实施例对本发明作详细说明。

具体实施方式

1、机体轴系三个速度分量输出为：

${\left(\begin{array}{c}u\left(t\right)\\ v\left(t\right)\\ w\left(t\right)\end{array}\right)}_{t=(k+1)T}={\Phi }_v[(k+1)T,\mathrm{kT}]{\left(\begin{array}{c}u\left(t\right)\\ v\left(t\right)\\ w\left(t\right)\end{array}\right)}_{t=\mathrm{kT}}+g{\Phi }_v[(k+1)T,\mathrm{kT}]{\Phi }_s[(k+1)T,\mathrm{kT}]{\left(\begin{array}{c}{s}_1\left(t\right)\\ s_2\left(t\right)\\ s_3\left(t\right)\end{array}\right)}_{t=\mathrm{kT}}$

$+g{\int }_{\mathrm{kT}}^{(k+1)T}{\Phi }_v[(k+1)T,\tau ]\left(\begin{array}{c}{n}_x\\ n_y\\ n_z\end{array}\right)\mathrm{d\tau }$

其中：u，v，w分别为沿刚体机体轴系x，y，z轴的速度分量，n_x，n_y，n_z分别为沿x，y，z轴的过载， g为重力加速度，s₁、s₂、s₃为定义的三元数，且

${\left(\begin{array}{c}{s}_1\left(t\right)\\ s_2\left(t\right)\\ s_3\left(t\right)\end{array}\right)}_{t=(k+1)T}={\Phi }_s[(k+1)T,\mathrm{kT}]{\left(\begin{array}{c}{s}_1\left(t\right)\\ s_2\left(t\right)\\ s_3\left(t\right)\end{array}\right)}_{t=\mathrm{kT}}$

${\Phi }_v[(k+1)T,\mathrm{kT}]\approx I+{\Pi }_v\zeta +{\Pi }_v\Omega {\Pi }_v^T$

${\Phi }_s[(k+1)T,\mathrm{kT}]\approx I+{\Pi }_s\zeta +{\Pi }_s\Omega {\Pi }_s^T$

p，q，r分别为滚转、俯仰、偏航角速度，T为采样周期；全文参数定义相同；

$I=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right),$ξ(t)＝[ξ₀(t)ξ₁(t)…ξ_n-1(t)ξ_n(t)]^T，

$\left(\begin{array}{c}{\xi }_0\left(t\right)=1\\ {\xi }_1\left(t\right)=\cos \left[a{\cos }^{-1}(1-2t/b)\right]\\ {\xi }_2\left(t\right)={2\xi }_1\left(t\right)·{\xi }_1\left(t\right)-1\\ ·\\ ·\\ ·\\ {\xi }_{i+1}\left(t\right)={2\xi }_1\left(t\right){\xi }_i\left(t\right)-{\xi }_{i-1}\left(t\right)\end{array}\right)i=\mathrm{2,3},···,n-\mathrm{1,0}\le t\le \mathrm{NT},b=\mathrm{NT}$ 为基于Chebyshev(切比雪夫)正交多项式得到的任意步长正交多项式的递推形式，a为任 意实数，滚转、俯仰、偏航角速度p，q，r的展开式分别为

p(t)＝[p₀ p₁…p_n-1 p_n][ξ₀(t)ξ₁(t)…ξ_n-1(t)ξ_n(t)]^T

q(t)＝[q₀ q₁…q_n-1 q_n][ξ₀(t)ξ₁(t)…ξ_n-1(t)ξ_n(t)]^T

r(t)＝[r₀ r₁…r_n-1 r_n][ξ₀(t)ξ₁(t)…ξ_n-1(t)ξ_n(t)]^T

${\Pi }_v=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 1\\ 0& -1& \text{0}\end{array}\right)\left(\begin{array}{ccccc}{p}_0& p_1& ···& p_{n-1}& p_n\end{array}\right)$

$+\left(\begin{array}{ccc}0& 0& -1\\ 0& 0& 0\\ 1& 0& 0\end{array}\right)\left(\begin{array}{ccccc}{q}_0& q_1& ···& q_{n-1}& q_n\end{array}\right)+\left(\begin{array}{ccc}0& 1& 0\\ -1& 0& 0\\ 0& 0& 0\end{array}\right)\left(\begin{array}{ccccc}{r}_0& r_1& ···& r_{n-1}& r_n\end{array}\right)$

${\Pi }_s=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 1\\ 0& -1& 0\end{array}\right)\left(\begin{array}{ccccc}{p}_0& p_1& ···& p_{n-1}& p_n\end{array}\right)$

$+\left(\begin{array}{ccc}0& 0& 1\\ 0& 0& 0\\ -1& 0& 0\end{array}\right)\left(\begin{array}{ccccc}{q}_0& q_1& ···& q_{n-1}& q_n\end{array}\right)+\left(\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right)\left(\begin{array}{ccccc}{r}_0& r_1& ···& r_{n-1}& r_n\end{array}\right)$

$\zeta ={\left(\begin{array}{cccc}{\zeta }_0& {\zeta }_1& ···& {\zeta }_n\end{array}\right)}^T={\int }_{\mathrm{kT}}^{(k+1)T}\xi \left(t\right)\mathrm{dt}$

${\zeta }_i={\int }_{\mathrm{kT}}^{(k+1)T}{\xi }_i\left(t\right)\mathrm{dt}={\int }_{\mathrm{kT}}^{(k+1)T}\cos \left[\mathrm{ai}{\cos }^{-1}(1-2t/b)\right]\mathrm{dt}$

$=\frac{b}{4}\left\{\frac{1}{\mathrm{ai}-1}\cos \right[(\mathrm{ai}-1){\cos }^{-1}(1-2t/b)]$

${-\frac{1}{\mathrm{ai}+1}\cos \left[(\mathrm{ai}+1){\cos }^{-1}(1-2t/b)\right]\left\}\right|}_{\mathrm{kT}}^{(k+1)T}$

$\Omega ={\left\{{\Omega }_{\mathrm{ji}}\right\}}_{j=\mathrm{0,1},···,n;i=\mathrm{1,2},···,n}={\int }_{\mathrm{kT}}^{(k+1)T}\xi \left(t\right){\int }_{\mathrm{kT}}^T\xi \left(\tau \right)\mathrm{d\tau dt}$

${\Omega }_{\mathrm{ji}}={\int }_{\mathrm{kT}}^{(k+1)T}{\xi }_j\left(t\right){\int }_{\mathrm{kT}}^t{\xi }_i\left(\tau \right)\mathrm{d\tau dt}$

$=\frac{b}{8}\left\{\frac{1}{\mathrm{ai}-1}{\int }_{\mathrm{kT}}^{(k+1)T}\right\{\cos \left[(\mathrm{aj}-\mathrm{ai}+1){\cos }^{-1}(1-2t/b)\right]+\cos \left[(\mathrm{aj}+\mathrm{ai}-1){\cos }^{-1}(1-2t/b)\right]\}\mathrm{dt}$

$-\frac{1}{\mathrm{ai}+1}{\int }_{\mathrm{kT}}^{(k+1)T}\left\{\cos \right[(\mathrm{aj}-\mathrm{ai}-1){\cos }^{-1}(1-2t/b)]+\cos [(\mathrm{aj}+\mathrm{ai}+1){\cos }^{-1}(1-2t/b)\left]\right\}\mathrm{dt}\};$

$-\frac{b}{4}\left\{\frac{1}{\mathrm{ai}-1}\cos \right[(\mathrm{ai}-1){\cos }^{-1}(1-2\mathrm{kT}/b)]$

$-\frac{1}{\mathrm{ai}+1}\cos \left[(\mathrm{ai}+1){\cos }^{-1}(1-2\mathrm{kT}/b)\right]\left\}{\int }_{\mathrm{kT}}^{(k+1)T}\cos \right[{\mathrm{aj}\cos }^{-1}(1-2t/b)]\mathrm{dt}$

2、高度输出为：

$\stackrel{·}{h}=\left(\begin{array}{ccc}u& v& w\end{array}\right)\left(\begin{array}{c}{s}_1\\ s_2\\ s_3\end{array}\right)$

其中：h为高度；

3、姿态角的输出为：。

$\psi \left[(k+1)T\right]=\psi \left(\mathrm{kT}\right)+{\int }_{\mathrm{kT}}^{(k+1)T}\frac{{\mathrm{qs}}_2\left(t\right)+{\mathrm{rs}}_3\left(t\right)}{s_2^2\left(t\right)+s_3^2\left(t\right)}\mathrm{dt}$

其中：θ，ψ分别表示滚转、俯仰、偏航角，${\left(\begin{array}{c}{s}_1\left(t\right)\\ s_2\left(t\right)\\ s_3\left(t\right)\end{array}\right)}_{t=(k+1)T}={\Phi }_s[{(k+1)T,\mathrm{kT})\left(\begin{array}{c}{s}_1\left(t\right)\\ s_2\left(t\right)\\ s_3\left(t\right)\end{array}\right)}_{t=\mathrm{kT}}.$