CN102506870A - 基于角速度的欧拉角埃米特指数近似输出方法 - Google Patents

Abstract

本发明公开了一种基于角速度的欧拉角埃米特指数近似输出方法，用于解决现有的飞行器机动飞行时欧拉角输出精度差的技术问题。技术方案是通过引入多个参数并将滚转、俯仰、偏航角速度按照埃米特正交多项式展开，通过按照依次求解俯仰角、滚转角、偏航角，直接对欧拉角的表达式进行高阶逼近积分，使得欧拉角的求解按照超线性逼近，从而保证了确定欧拉角的时间更新迭代计算精度和惯性单元的输出精度。

Description

基于角速度的欧拉角埃米特指数近似输出方法

技术领域

本发明涉及一种飞行器机动飞行姿态确定方法，特别是涉及一种基于角速度的欧拉角埃米特指数近似输出方法。 

背景技术

惯性设备在运动体导航和控制中具有重要作用；刚体运动的加速度、角速度和姿态等通常都依赖于惯性设备输出，因此提高惯性设备的输出精度具有明确的实际意义；在惯性设备中，加速度采用加速度计、角速度采用角速率陀螺直接测量方式，刚体的姿态精度要求很高时如飞行试验等采用姿态陀螺测量，但在很多应用领域都有角速度等测量直接解算输出；主要原因是由于动态姿态传感器价格昂贵、体积大，导致很多飞行器采用角速率陀螺等解算三个欧拉角，使得姿态时间更新输出成为导航等核心内容，也使其成为影响惯导系统精度的主要因素之一，因此设计和采用合理的姿态时间更新输出方法就成为研究的热点课题；从公开发表的文献中对姿态输出主要基于角速度采用欧拉方程直接近似法或采用近似龙格库塔方法解算(孙丽、秦永元，捷联惯导系统姿态算法比较，中国惯性技术学报，2006，Vol.14(3)：6-10；Pu Li，Wang TianMiao，Liang JianHong，Wang Song，An Attitude Estimate Approach using MEMS Sensors for Small UAVs，2006，IEEE International Conference on Industrial Informatics，1113-1117)：由于欧拉方程中三个欧拉角互相耦合，属于非线性微分方程，在不同初始条件和不同飞行状态下的误差范围不同，难以保证实际工程要求的精度。 

发明内容

为了克服现有的飞行器机动飞行时欧拉角输出精度差的问题，本发明提供一种基于角速度的欧拉角埃米特指数近似输出方法。该方法通过引入多个参数并将滚转、俯仰、偏航角速度按照埃米特正交多项式展开，通过按照依次求解俯仰角、滚转角、偏航角，直接对欧拉角的表达式进行高阶逼近积分，使得欧拉角的求解按照超线性逼近，从而可以保证确定欧拉角的时间更新迭代计算精度和惯性单元的输出精度。 

本发明解决其技术问题所采用的技术方案是：一种基于角速度的欧拉角埃米特指数近似输出方法，其特点是包括以下步骤： 

1、(a)根据欧拉方程： 

式中： 分别指滚转、俯仰、偏航角；p，q，r分别为滚转、俯仰、偏航角速度；全文参数定义相同；这三个欧拉角的计算按照依次求解俯仰角、滚转角、偏航角的步骤进行；全文符号定义相同；滚转、俯仰、偏航角速度p，q，r的展开式分别为 

p(t)＝pξ，q(t)＝qξ，r(t)＝rξ 

其中 

p＝[p₀p₁L  p_n-1p_n]q＝[q₀q₁L  q_n-1q_n] 

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

$\left.\begin{array}{c}{\mathrm{\ξ}}_0\left(t\right)=1\\ {\mathrm{\ξ}}_1\left(t\right)=2t\\ {\mathrm{\ξ}}_2\left(t\right)={4t}^2-2\\ {\mathrm{\ξ}}_3\left(t\right)={8t}^3-12t\\ {\mathrm{\ξ}}_4\left(t\right)={16t}^4-{48t}^2+12\\ {\mathrm{\ξ}}_5\left(t\right)={32t}^5-{160t}^3+120t\\ {\mathrm{\ξ}}_6\left(t\right)={64t}^6-{480t}^4+720t-120\\ M\\ {\mathrm{\ξ}}_{i+1}\left(t\right)={2\mathrm{t\ξ}}_i\left(t\right)-{2\mathrm{i\ξ}}_{i-1}\left(t\right)\end{array}\right\},i=\mathrm{2,3},L,n-1$

为埃米特正交多项式的递推形式； 

(b)俯仰角的时间更新求解式为： 

式中： 

$a_1={\left(\mathrm{qH}{\mathrm{\ξ}|}_{\mathrm{kT}}^{(k+1)T}\right)}^2+{\left({\mathrm{rH\ξ}|}_{\mathrm{kT}}^{(k+1)T}\right)}^2-{\left({\mathrm{pH\ξ}|}_{\mathrm{kT}}^{(k+1)T}\right)}^2$

$a_2=p{\∫}_{\mathrm{kT}}^{(k+1)T}\left[\mathrm{\ξ}\left(t\right){\mathrm{\ξ}}^T\left(t\right)\right]\mathrm{dt}{H}^T{r}^T-p{\mathrm{H\ξ}|}_{\mathrm{kT}}^{(k+1)T}r{\mathrm{H\ξ}|}_{\mathrm{kT}}$

$a_3=p{\∫}_{\mathrm{kT}}^{(k+1)T}\left[\mathrm{\ξ}\left(t\right){\mathrm{\ξ}}^T\left(t\right)\right]\mathrm{dt}{H}^T{q}^T-p{\mathrm{H\ξ}|}_{\mathrm{kT}}^{(k+1)T}q{\mathrm{H\ξ}|}_{\mathrm{kT}}$

$\left|\mathrm{\λ}\right|=\left\{p{\∫}_{\mathrm{kT}}^{(k+1)T}\right[\mathrm{\ξ}\left(t\right){\mathrm{\ξ}}^T\left(t\right)]\mathrm{dt}{H}^T{p}^T-p{\mathrm{H\ξ}|}_{\mathrm{kT}}^{(k+1)T}p{\mathrm{H\ξ}|}_{\mathrm{kT}}$

$+q{\∫}_{\mathrm{kT}}^{(k+1)T}\left[\mathrm{\ξ}\left(t\right){\mathrm{\ξ}}^T\left(t\right)\right]\mathrm{dt}{H}^T{q}^T-q{\mathrm{H\ξ}|}_{\mathrm{kT}}^{(k+1)T}q{\mathrm{H\ξ}|}_{\mathrm{kT}}$

$+r{\∫}_{\mathrm{kT}}^{(k+1)T}\left[\mathrm{\ξ}\left(t\right){\mathrm{\ξ}}^T\left(t\right)\right]\mathrm{dt}{H}^T{r}^T-r{\mathrm{H\ξ}|}_{\mathrm{kT}}^{(k+1)T}r{\mathrm{H\ξ}|}_{\mathrm{kT}}{\}}^{\frac{1}{2}}$

T为采样周期， 

$H={\left\{h_{\mathrm{ij}}\right\}}_{(n+1)(n+2)}=\left[\begin{array}{ccccccc}0& 0.5& 0& 0& L& 0& 0\\ 0.5& 0& 0.25& 0& L& 0& 0\\ 0& 0& 0& \frac{1}{6}& L& 0& 0\\ -1.5& 0& 0& 0& L& 0& 0\\ M& M& M& O& M& M& M\\ h_{n1}& 0& 0& 0& L& h_{n(n+1)}& 0\\ h_{(n+1)1}& 0& 0& 0& L& 0& h_{(n+1)(n+2)}\end{array}\right]$

h₁₂＝0.5，h₂₃＝0.25， $h_{i(i+1)}=\frac{h_{(i-1)i}}{1+2h_{(i-1)i}},L$

h₂₁＝0.5，h₄₁＝-1.5， $h_{\left(2i\right)1}=\frac{4i{h}_{(2i-2)1}}{1+2h_{\left(2i\right)(2i+1)}},L$

其余h_ij＝0； 

2、在已知俯仰角的情况下，滚转角的时间更新求解式为： 

其中 

$a_4={\left(\mathrm{pH}{\mathrm{\ξ}|}_{\mathrm{kT}}^{(k+1)T}\right)}^2+{\left({\mathrm{rH\ξ}|}_{\mathrm{kT}}^{(k+1)T}\right)}^2-{\left({\mathrm{qH\ξ}|}_{\mathrm{kT}}^{(k+1)T}\right)}^2$

$a_5=q{\∫}_{\mathrm{kT}}^{(k+1)T}\left[\mathrm{\ξ}\left(t\right){\mathrm{\ξ}}^T\left(t\right)\right]\mathrm{dt}{H}^T{p}^T-q{\mathrm{H\ξ}|}_{\mathrm{kT}}^{(k+1)T}p{\mathrm{H\ξ}|}_{\mathrm{kT}}$

$a_6=q{\∫}_{\mathrm{kT}}^{(k+1)T}\left[\mathrm{\ξ}\left(t\right){\mathrm{\ξ}}^T\left(t\right)\right]\mathrm{dt}{H}^T{r}^T-q{\mathrm{H\ξ}|}_{\mathrm{kT}}^{(k+1)T}r{\mathrm{H\ξ}|}_{\mathrm{kT}}$

3、在俯仰角、滚转角已知情况下，偏航角的求解式为： 

$\mathrm{\ψ}\left(t\right)=\mathrm{\ψ}\left(\mathrm{kT}\right)+{\∫}_{\mathrm{kT}}^t[b_1\left(t\right)+b_2\left(t\right)]\mathrm{dt}$

式中： 

本发明的有益效果是：由于通过引入多个参数并将滚转、俯仰、偏航角速度按照埃米特正交多项式展开，通过按照依次求解俯仰角、滚转角、偏航角，直接对欧拉角的表达式进行高阶逼近积分，使得欧拉角的求解按照超线性逼近，从而保证了确定欧拉角的时间更新迭代计算精度和惯性单元的输出精度。 

下面结合具体实施方式对本发明作详细说明。 

具体实施方式

1、(a)根据刚体姿态方程(欧拉方程)： 

其中： 

分别指滚转、俯仰、偏航角；p，q，r分别为滚转、俯仰、偏航角速度；全文参数定义相同；这三个欧拉角的计算按照依次求解俯仰角、滚转角、偏航角的步骤进行；滚转、俯仰、偏航角速度p，q，r的展开式分别为 

p(t)＝pξ，q(t)＝qξ，r(t)＝rξ 

p＝[p₀p₁L  p_n-1p_n]q＝[q₀q₁L  q_n-1q_n] 

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

其中 

$\left.\begin{array}{c}{\mathrm{\ξ}}_0\left(t\right)=1\\ {\mathrm{\ξ}}_1\left(t\right)=2t\\ {\mathrm{\ξ}}_2\left(t\right)={4t}^2-2\\ {\mathrm{\ξ}}_3\left(t\right)={8t}^3-12t\\ {\mathrm{\ξ}}_4\left(t\right)={16t}^4-{48t}^2+12\\ {\mathrm{\ξ}}_5\left(t\right)={32t}^5-{160t}^3+120t\\ {\mathrm{\ξ}}_6\left(t\right)={64t}^6-{480t}^4+720t-120\\ M\\ {\mathrm{\ξ}}_{i+1}\left(t\right)={2\mathrm{t\ξ}}_i\left(t\right)-{2\mathrm{i\ξ}}_{i-1}\left(t\right)\end{array}\right\},i=\mathrm{2,3},L,n-1$

为埃米特正交多项式的递推形式，T为采样周期； 

(b)俯仰角的时间更新求解式为： 

式中： 

$a_1={\left(\mathrm{qH}{\mathrm{\ξ}|}_{\mathrm{kT}}^{(k+1)T}\right)}^2+{\left({\mathrm{rH\ξ}|}_{\mathrm{kT}}^{(k+1)T}\right)}^2-{\left({\mathrm{pH\ξ}|}_{\mathrm{kT}}^{(k+1)T}\right)}^2$

$a_2=p{\∫}_{\mathrm{kT}}^{(k+1)T}\left[\mathrm{\ξ}\left(t\right){\mathrm{\ξ}}^T\left(t\right)\right]\mathrm{dt}{H}^T{r}^T-p{\mathrm{H\ξ}|}_{\mathrm{kT}}^{(k+1)T}r{\mathrm{H\ξ}|}_{\mathrm{kT}}$

$a_3=p{\∫}_{\mathrm{kT}}^{(k+1)T}\left[\mathrm{\ξ}\left(t\right){\mathrm{\ξ}}^T\left(t\right)\right]\mathrm{dt}{H}^T{q}^T-p{\mathrm{H\ξ}|}_{\mathrm{kT}}^{(k+1)T}q{\mathrm{H\ξ}|}_{\mathrm{kT}}$

$\left|\mathrm{\λ}\right|=\left\{p{\∫}_{\mathrm{kT}}^{(k+1)T}\right[\mathrm{\ξ}\left(t\right){\mathrm{\ξ}}^T\left(t\right)]\mathrm{dt}{H}^T{p}^T-p{\mathrm{H\ξ}|}_{\mathrm{kT}}^{(k+1)T}p{\mathrm{H\ξ}|}_{\mathrm{kT}}$

$+q{\∫}_{\mathrm{kT}}^{(k+1)T}\left[\mathrm{\ξ}\left(t\right){\mathrm{\ξ}}^T\left(t\right)\right]\mathrm{dt}{H}^T{q}^T-q{\mathrm{H\ξ}|}_{\mathrm{kT}}^{(k+1)T}q{\mathrm{H\ξ}|}_{\mathrm{kT}}$

$+r{\∫}_{\mathrm{kT}}^{(k+1)T}\left[\mathrm{\ξ}\left(t\right){\mathrm{\ξ}}^T\left(t\right)\right]\mathrm{dt}{H}^T{r}^T-r{\mathrm{H\ξ}|}_{\mathrm{kT}}^{(k+1)T}r{\mathrm{H\ξ}|}_{\mathrm{kT}}{\}}^{\frac{1}{2}}$

$H={\left\{h_{\mathrm{ij}}\right\}}_{(n+1)(n+2)}=\left[\begin{array}{ccccccc}0& 0.5& 0& 0& L& 0& 0\\ 0.5& 0& 0.25& 0& L& 0& 0\\ 0& 0& 0& \frac{1}{6}& L& 0& 0\\ -1.5& 0& 0& 0& L& 0& 0\\ M& M& M& O& M& M& M\\ h_{n1}& 0& 0& 0& L& h_{n(n+1)}& 0\\ h_{(n+1)1}& 0& 0& 0& L& 0& h_{(n+1)(n+2)}\end{array}\right]$

h₁₂＝0.5，h₂₃＝0.25， $h_{i(i+1)}=\frac{h_{(i-1)i}}{1+2h_{(i-1)i}},L$

h₂₁＝0.5，h₄₁＝-1.5， $h_{\left(2i\right)1}=\frac{4i{h}_{(2i-2)1}}{1+2h_{\left(2i\right)(2i+1)}},L$

其余h_ij＝0； 

2、在已知俯仰角的情况下，滚转角的时间更新求解式为： 

其中 

$a_4={\left(\mathrm{pH}{\mathrm{\ξ}|}_{\mathrm{kT}}^{(k+1)T}\right)}^2+{\left({\mathrm{rH\ξ}|}_{\mathrm{kT}}^{(k+1)T}\right)}^2-{\left({\mathrm{qH\ξ}|}_{\mathrm{kT}}^{(k+1)T}\right)}^2$

$a_5=q{\∫}_{\mathrm{kT}}^{(k+1)T}\left[\mathrm{\ξ}\left(t\right){\mathrm{\ξ}}^T\left(t\right)\right]\mathrm{dt}{H}^T{p}^T-q{\mathrm{H\ξ}|}_{\mathrm{kT}}^{(k+1)T}p{\mathrm{H\ξ}|}_{\mathrm{kT}}$

$a_6=q{\∫}_{\mathrm{kT}}^{(k+1)T}\left[\mathrm{\ξ}\left(t\right){\mathrm{\ξ}}^T\left(t\right)\right]\mathrm{dt}{H}^T{r}^T-q{\mathrm{H\ξ}|}_{\mathrm{kT}}^{(k+1)T}r{\mathrm{H\ξ}|}_{\mathrm{kT}}$

3、在俯仰角、滚转角已知情况下，偏航角的求解式为： 

$\mathrm{\ψ}\left(t\right)=\mathrm{\ψ}\left(\mathrm{kT}\right)+{\∫}_{\mathrm{kT}}^t[b_1\left(t\right)+b_2\left(t\right)]\mathrm{dt}$

式中： 

当对惯性设备直接输出滚转、俯仰、偏航角速度p，q，r采用三阶逼近描述时，所得结果也接近O(T³)，相比欧拉方程直接近似法或采用近似龙格库塔方法解算等方法的O(T²)精度要高。 

Claims (1)

1.一种基于角速度的欧拉角埃米特指数近似输出方法，其特征在于包括以下步骤：

步骤1、(a)根据欧拉方程：

式中： 分别指滚转、俯仰、偏航角；p，q，r分别为滚转、俯仰、偏航角速度；全文参数定义相同；这三个欧拉角的计算按照依次求解俯仰角、滚转角、偏航角的步骤进行；全文符号定义相同；滚转、俯仰、偏航角速度p，q，r的展开式分别为

p(t)＝pξ，q(t)＝qξ，r(t)＝rξ

其中

p＝[p₀ p₁ L  p_n-1 p_n]q＝[q₀ q₁ L  q_n-1 q_n]

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

为埃米特正交多项式的递推形式；

(b)俯仰角的时间更新求解式为：

式中： 

T为采样周期，

h₁₂＝0.5，h₂₃＝0.25，

h₂₁＝0.5，h₄₁＝-1.5，

其余h_ij＝0；

步骤2、在已知俯仰角的情况下，滚转角的时间更新求解式为：

其中

步骤3、在俯仰角、滚转角已知情况下，偏航角的求解式为：

式中：