在工程实际中，特别是用数值法解运动学方程时，四元数运动学方程得到了广泛的应用，这是因为四元数运动学方程对任何参数值都不退化(无奇点)，而且参数数目少(只有4个)，联系式数目最少(只有1个，即  。假设运载火箭的角运动是刚体角运动，可得四元数运动学方程：

  在ωx1，ωy1，ωz1已知的情况下，由式(4)进行差分法计算可解算出q0，q1，q2，q3。

1   四元数变换

设A为箭体坐标系到惯性坐标系的转换矩阵，则有

式中

A 即为所说的由箭体到惯性系的数学平台转换矩阵，它是通过解算四元数运动方程而得到的。

2  四元数方程的解算

       由式(4)用数值迭代法可解算出q0，q1，q2，q3。捷联系统中激光陀螺输出脉冲数代表箭体在采样时间内的转角增量，若已知箭体的初始转动四元数，并测得各采样时刻箭体的转角增量，则可用下述的一种数值迭代法求解四元数微分方程。

       设陀螺仪测出箭体从第(k-1)次采样时间到第k次采样时间内的转角增量为

式中 Δθx1(k)，Δθy1(k)，Δθz1(k)分别为从第k-1次到第k次采样时间内箭体绕3个坐标轴的转角增量；tk-1为第k-1次采样时的时间；tk为第k次采样时的时间(k=1,2,3, …)。

从第k-1次到第k次采样时间内，箭体转角增量的模为

则四元数的差分解为：

3  四元数的初值

       按式(7)求解转动四元数时，需要知道第1次采样时刻（k=1）的初始四元数q(0)，通常q(0)值可在箭体进行初始对准时求得。选用运载火箭起飞时刻发射点的惯性坐标系作为导航坐标系oxyz，且箭体坐标系ox1y1z1相对导航坐标系的初始不对准角度以图1所示的克雷洛夫角表示。按图1所示，箭体坐标系相对导航坐标系的初始位置可由与克雷洛夫角转动顺序有关的四元数表示，这些四元数可写为

         克雷洛夫角的转动顺序为φ0→ψ0→γ0，因此可得合成转动四元数q(0)=qφ0qψ0qγ0。

将式(8)代入上式，并按四元数乘法规则展开，可得出以矩阵形式给出火箭起飞时的初始四元数为

图1 坐标变换关系图

function Ztaij=Angle()
global Q
Pi=3.1415926;
q0=Q(1);
q1=Q(2);
q2=Q(3);
q3=Q(4);
T11=q1^2+q0^2-q3^2-q2^2;%初始姿态矩阵元素确定
T12=2*(q1*q2-q0*q3);
T13=2*(q1*q3+q0*q2);
T21=2*(q1*q2+q0*q3);
T22=q2^2-q3^2+q0^2-q1^2;
T23=2*(q2*q3-q0*q1);
T31=2*(q1*q3-q0*q2);
T32=2*(q2*q3+q0*q1);
T33=q3^2-q2^2-q1^2+q0^2;
a=asin(T32);
b_1=atan(-T31/T33);
if abs(T22)<0.0001
    if T12<0
        b=-Pi/2;
    else b=Pi/2;
    end
elseif T22>0
    b=b_1;
    elseif T12>0
        b=b_1+Pi;
else  b=b_1-Pi;
end
c_1=atan(T12/T22);
if T33>=0
    c=c_1;
elseif c_1>=0
    c=c_1-Pi;
else  c=c_1+Pi;
end
Ztaij=[a c b]'