捷联式姿态解算过程

第三部分：基于“存储数据与迭代计算对准”罗经法对准
3.2 罗经法对准过程中的调整策略（以北向通道为例）
g
y
f
p N
1
VN
1
s
R
-
K1
K2 R
K3 s
Control algorithm
cE -
x
1
s
x
-
ie cos L
z
实线所示的北向通道：本质上是一个休拉回路，失准角作无阻尼振荡。
采取的策略：1）引入内反馈环节（虚线）实现衰减振荡；2）引入前馈环节（点画 线）缩短振荡周期；3）引入积分环节（双点画线）消除罗经项的影响。
3.5 SINS罗经法对准如何实现迭代计算？
fˆNn -
b ib
fb cU
Cˆbn
Cˆbn
b ib
Cnbine
Cnbc
cN
fˆ n Cˆbn f b
1
VN
s
1
cE
R
K1
K2 R
K3 s
Control algorithm
上述过程中，可以实现迭代计算。
Page 15
第三部分：基于“存储数据与迭代计算对准”罗经法对准
导航坐标系 n （b）SINS
GINS中的测量数据直接反映失准角的大小； SINS中的测量数据不直接反映失准角；只有投影数据能够反映失准角的大小；相同 的测量数据经过不同的姿态矩阵进行投影，可以获取不同的投影数据。 注：上述均不考虑仪表误差。
对于SINS而言，分析一种理想的情况：仪表无误差，载体无机动，此时在整个对准 过程中，仪表测量数据均相等。整个对准过程，其实只用了一组仪表参数。
3.6 SINS罗经法对准中存储数据如何使用？

Pixhawk之姿态解算篇（1）_⼊门篇（DCMNomalize）⼀、开篇慢慢的、慢慢的、慢慢的就快要到飞控的主要部分了，飞控飞控就是所谓的飞⾏控制呗，⼀个是姿态解算⼀个是姿态控制，解算是解算，控制是控制，各⾃负责各⾃的任务。

我也不懂。

还在学习中~~~~近期看姿态预计部分看的太累了，明显发现基础知识太薄弱，什么欧拉⾓、DCM、四元数、gyro误差、矫正、正交化等各个概念。

然后就是各种转换公式。

接下来结合代码介绍⼀些主要的东西。

太深⼊的还不了解~~~⼀定要多看论⽂啊，英⽂版的论⽂（也没有中⽂的。

国⼈的悲哀啊）。

尽管看着头疼，看是看完了以后就会发现很多不了解的迷惑的就⾃然都解开了。

三、实验平台Software Version：ArduCopter（Ver_3.3）Hardware Version：pixhawkIDE：eclipse Juno （Windows）四、基本知识1、怎样实现控制⼀个⽆⼈机系统的算法主要有两类：姿态检測算法、姿态控制算法。

姿态控制、被控对象、姿态检測三个部分构成⼀个闭环控制系统。

被控对象的模型是由其物理系统决定，设计⽆⼈机的算法就是设计姿态控制算法、姿态检測算法。

1）姿态检測算法：姿态的表⽰能够⽤欧拉⾓，也能够⽤四元数。

姿态检測算法的作⽤就是将加速度计、陀螺仪等传感器的測量值解算成姿态，进⽽作为系统的反馈量。

在获取sensors值之前须要对数据进⾏滤波，滤波算法主要是将获取到的陀螺仪和加速度计的数据进⾏去噪声及融合，得出正确的⾓度数据（欧拉⾓或四元数），主要採⽤互补滤波或者⾼⼤上的卡尔曼滤波。

2）姿态控制算法：控制⽆⼈机姿态的三个⾃由度，⽤给定姿态与姿态检測算法得出的姿态偏差作为输⼊，被控对象模型的输⼊量作为输出（⽐如姿态增量），从⽽达到控制⽆⼈机姿态的作⽤。

最经常使⽤的就是PID控制及其各种PID扩展（分段、模糊等，我的毕设就是模糊PID控制算法，慘了，啥都不懂，还能毕业不，哎~~~），⾼端点的有⾃适应控制（⾃⼰主动壁障应该就⽤这个）。

捷联惯性导航系统的解算方法捷联惯性导航系统（Inertial Navigation System，简称INS）是一种利用陀螺仪和加速度计等惯性测量单元测量物体的加速度和角速度，然后通过对这些测量值的积分计算出物体的速度和位置的导航系统。

INS广泛应用于航空航天、无人驾驶车辆和船舶等领域，具有高精度和自主性等特点。

INS的解算方法一般分为初始对准、运动状态估计和航位推算三个主要过程。

初始对准是指在启动导航系统时，通过利用外部辅助传感器（如GPS）或静态校准等方法将惯性传感器的输出与真实姿态和位置进行初次校准。

在初始对准过程中，需要获取传感器的初始偏差和初始姿态，一般采用标定或矩阵运算等方法进行。

运动状态估计是指根据惯性传感器的测量值，使用滤波算法对物体的加速度和角速度进行实时估计。

常用的滤波算法包括卡尔曼滤波、扩展卡尔曼滤波和粒子滤波等。

其中，卡尔曼滤波是一种最优估计算法，通过对观测值和状态进行线性组合，得到对真实状态的最佳估计。

扩展卡尔曼滤波则是基于卡尔曼滤波的非线性扩展，可以应用于非线性INS系统。

粒子滤波是一种利用蒙特卡洛采样技术进行状态估计的方法，适用于非高斯分布的状态估计问题。

航位推算是指根据运动状态估计的结果，对物体的速度和位置进行推算。

INS最基本的航位推算方法是利用加速度值对速度进行积分，然后再对速度进行积分得到位置。

但是，在实际应用中，由于传感器本身存在噪声和漂移等误差，导致航位推算过程会出现积分漂移现象。

为了解决这个问题，通常采用辅助传感器（如GPS）和地图等数据对INS的输出进行校正和修正。

当前，还有一些先进的INS解算方法被提出，如基于深度学习的INS 解算方法。

这些方法利用神经网络等深度学习模型，结合原始传感器数据进行端到端的学习和预测，以实现更高精度的位置和姿态估计。

综上所述，捷联惯性导航系统的解算方法主要包括初始对准、运动状态估计和航位推算三个过程。

其中，运动状态估计过程利用滤波算法对传感器的测量值进行处理，得到物体的加速度和角速度的估计。

收稿日期 :2008208207 ; 修订日期 :2008209205. 基金项目 : 吉林省杰出青年基金资助项目 ( No . 20060115) ; 中国科学院 “三期创新” 平台资助项目
© 1994-2009 China Academic Journal Electronic Publishing House. All rights reserved.
摘要 : 提出了一种实现捷联惯性导航系统姿态解算模块的方法 。基于 DSP 的硬件平台和四元数的数学平台 ,设计了捷 联惯性导航系统的姿态解算模块 。介绍了捷联惯性导航系统的工作原理和姿态解算的基本算法 , 并给出了四元数法的 四阶龙格2库塔数值解法 。设计了姿态解算模块的硬件电路和软件实现程序 。实验测试结果表明 ,在增量角 < 5° 的情况 下 ,用四阶龙格2库塔法进行姿态解算 , 误差 < 0. 005 3 % ; 应用 TMS320C6713B 进行硬件电路设计 , 每次解算时间 <

和调试 。综上所述 ,DSP 系统具有速度快 、 精度 高、 体积小 、 成本低 、 可靠性高以及抗干扰能力强 等优点 ,可以满足捷联式惯性导航系统对硬件平 台的要求 [ 122 ] 。
© 1994-2009 China Academic Journal Electronic Publishing House. All rights reserved.
,
( 5)
2
K23 T
2
如下 : ψ θ γ ψ θ γ 0 0 0 0 0 0 cos cos cos + sin sin sin 2 2 2 2 2 2 λ (0) ψ θ γ ψ θ γ 0 0 0 0 0 0 cos cos sin - sin sin cos p1 (0) 2 2 2 2 2 2 = , ψ θ γ ψ θ γ p2 (0) 0 0 0 0 0 0 cos sin cos + sin cos sin 2 2 2 2 2 2 p3 (0) ψ θ γ ψ θ γ 0 0 0 0 0 0 sin cos cos - cos sin sin 2 2 2 2 2 2

02
CATALOGUE
捷联惯性导航系统组成及工作 原理
主要组成部分介绍
惯性测量单元
包括加速度计和陀螺仪，用于测量载体在三个正交轴上的加速度 和角速度。
导航计算机
用于处理惯性测量单元的测量数据，解算出载体的姿态、速度和 位置信息。
控制与显示单元
用于实现人机交互，包括设置导航参数、显示导航信息等。
工作原理简述
学生自我评价报告
知识掌握情况
学生对捷联惯性导航系统的基本原理、解算 方法和实现技术有了深入的理解和掌握。
实践能力提升
通过实验和仿真，学生的动手实践能力得到了提升 ，能够独立完成相关的实验和仿真验证。
团队协作能力
在课程项目中，学生之间的团队协作能力得 到了锻炼和提升，能够相互协作完成项目任 务。
对未来发展趋势的预测和建议
捷联惯性导航系统的解算 方法课件
CATALOGUE
目 录
• 捷联惯性导航系统概述 • 捷联惯性导航系统组成及工作原理 • 捷联惯性导航系统解算方法 • 误差分析及补偿策略 • 实验验证与结果展示 • 总结与展望
01
CATALOGUE
捷联惯性导航系统概述
定义与基本原理
定义
捷联惯性导航系统是一种基于惯性测量元件（加速度计和陀螺仪）来测量载体（如飞机、导弹等）的加速度和角 速度，并通过积分运算得到载体位置、速度和姿态信息的自主导航系统。
01
高精度、高可靠性
02
多传感器融合技术
随着科技的发展和应用需求的提高， 捷联惯性导航系统需要进一步提高精 度和可靠性，以满足更高层次的应用 需求。
为了克服单一传感器的局限性，可以 采用多传感器融合技术，将捷联惯性 导航系统与其他传感器进行融合，提 高导航系统的性能和鲁棒性。

们可以利用前一时刻的角速度采样值对旋
转矢量双子样法中的角增量进行修正 ,由此 4 数字仿真
提出双子样改进算法 。
我们用 MATLAB 编写了仿真程序 ,采用
修正后 ,双子样法的角增量由以下公式 角速度输出的陀螺 ,在圆锥运动下进行了数
提取
字仿真 。不同算法的时间复杂度取捷联矩
Δθ1 = [ 8ω(t - T2/ 2) - ω(t) + 5ω(t + T2/ 阵一次递推所需运算次数 ,如下表所示
<20 8
+ζ·3<8404 ,d2≈0. 5 -
<20 48
Δ; (10) Δθ3 = [ω(t + 2T3/ 3) +ω(t + T3) ]·T3/ 6 +Δ;
式中 ζ, = 0. 8 ,目的是为了减小漂移误差 。 若认为在式 (5) 中 ,一个计算周期内ω
= a + b·t ,并作一定的近似 ,可推导出等效 旋转矢量双子样计算公式为
自动化技术与应用 Vol. 19 ,No. 1 2000 年第 1 期 7
q (t + T1) = q (t) + (k1 + k2 + k3 + k4) / 6
若认为在一个计算周期内ω= a + b·t +
(7) c·t2 ,可导出等效旋转矢量三子样计算公式
式中 ,T1 为计算周期 ; k1 = 0. 5·T1[ωb (t) ]·q (t) ; k2 = 0. 5·T1 [ωb (t + T1/ 2) ] ·( q (t) + k1/
对于双子样法 Δθ1 = [ω(t) +ω(t + T2/ 2) ]·T2/ 4 ;

b EbX cos H sin C sin H cos C sin H X b 0 cos H cosC - sin C cos H X EbY .
0
sin C
cos C
b X EbZ
( 1) 根据欧拉角微分方程, 由角速度可 以求解 3 个姿态角。欧拉角微分方程式只有 3 个, 但每个 方程 Û x = f ( co s x , sin x ) X 都含有三角函数的运 算, 计算速度慢 , 且方程会出现/ 奇点0 , 方程式退 化, 故不能全姿态工作。 2. 2 方向余弦矩阵微分方程式 当一个坐标系相对另一个坐标系做一次或多 次旋转后可得到另外一个新的坐标系 , 前者往往 被称为参考坐标系或固定坐标系, 后者被称为动 坐标系, 他们之间的相互关系可用方向余弦表来 表示。方向余弦矩阵微分方程式可写为载体坐标 系相对导航坐标系旋转角速度的斜对称矩阵表达 式, 方向余弦表是对这两种坐标系相对转动的一 种数学描述。
E b C ÛE b = Cb 8 Eb , b Eb
2
姿态矩阵的计算方法
由于载体的姿态方位角速率较大 , 所以针对
姿态矩阵的实时计算提出了更高的要求。通常假 定捷联系统/ 数学平台0 模拟地理坐标系 , 即导航 坐标系; 而确定载体的姿态矩阵即为研究载体坐 标系 ( b) 和导航坐标系 ( E ) 的空间转动关系 , 一般
第 16 卷
第 10 期
光学 精密工程
O pt ics and Precision Eng ineer ing
2008 年 10 月 文章编号 1004 -924X( 2008) 10 - 1963 - 08
V ol. 16 N o. 10 Oct. 2008
பைடு நூலகம்

%====本程序为捷联惯导的解算程序(由惯性器件的输出解算出飞行器的位置、速度、姿态信息)======clear all;close all;clc;deg_rad=pi/180; %由度转化成弧度rad_deg=180/pi; %由弧度转化成度%-------------------------------从源文件中读入数据----------------------------------fid_read=fopen('IMUout.txt','r'); %path1_Den.dat 是由轨迹发生器产生的数据[AllDataNumofAllData]=fscanf(fid_read,'%g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g',[17 inf]);AllData=AllData';NumofEachData=round(NumofAllData/17);Time=AllData(:,1);longitude=AllData(:,2); %经度单位：弧度latitude=AllData(:,3); %纬度单位：弧度High=AllData(:,4); %高度单位：米Ve=-AllData(:,6); % 东向、北向、天向速度单位：米/妙Vn=AllData(:,5);Vu=AllData(:,7);fb_x=AllData(:,9); %比力（fx，fy，fz）fb_y=AllData(:,8); %指向右机翼方向为x正方向，指向机头方向为y正向，z轴与x轴和y轴构成右手坐标系单位：米/秒2fb_z=-AllData(:,10); %右前上pitch=AllData(:,11); %俯仰角（向上为正）单位:弧度head=-AllData(:,13); %偏航角（偏西为正）roll=AllData(:,12); %滚转角（向右为正）omigax=AllData(:,15); %陀螺输出（单位：弧度/秒，坐标轴的定义与比力的相同）omigay=AllData(:,14);omigaz=-AllData(:,16);%-------------------------------程序初始化--------------------------------------latitude0=latitude(1);longitude0=longitude(1); %初始位置High0=High(1);Ve0=Ve(1);Vn0=Vn(1); %初始速度Vu0=Vu(1);pitch0=pitch(1);head0=head(1); %初始姿态roll0=roll(1);TimeEach=0.005; %周期和仿真总时间TimeAll=(NumofEachData-1)*TimeEach;Omega_ie=0.7292115147E-4;%0.00007272205216643040; %地球自转角速度单位：弧度每妙g0=9.78;%------------------------------导航解算开始--------------------------------------%假设没有初始对准误差pitch_err0=pitch0+0*deg_rad;head_err0=head0+0*deg_rad;roll_err0=roll0+0*deg_rad;%初始捷联矩阵的计算《捷联惯导系统》P63 旋转顺序 head - pitch - roll %导航坐标系n为东北天方向载体坐标系b为右前上偏航角北偏西为正Tbn(1,1)=cos(roll_err0)*cos(head_err0)-sin(roll_err0)*sin(pitch_err0) *sin(head_err0);Tbn(1,2)=cos(roll_err0)*sin(head_err0)+sin(roll_err0)*sin(pitch_err0) *cos(head_err0);Tbn(1,3)=-sin(roll_err0)*cos(pitch_err0);Tbn(2,1)=-cos(pitch_err0)*sin(head_err0);Tbn(2,2)=cos(pitch_err0)*cos(head_err0);Tbn(2,3)=sin(pitch_err0);Tbn(3,1)=sin(roll_err0)*cos(head_err0)+cos(roll_err0)*sin(pitch_err0) *sin(head_err0);Tbn(3,2)=sin(roll_err0)*sin(head_err0)-cos(roll_err0)*sin(pitch_err0) *cos(head_err0);Tbn(3,3)=cos(roll_err0)*cos(pitch_err0);Tnb=Tbn';%位置矩阵的初始化《捷联惯导系统》P46 其中游动方位角 a＝0 假使初始经纬度确知Cne(1,1) = - sin(longitude0);Cne(1,2) = cos(longitude0);Cne(1,3) = 0;Cne(2,1) = - sin(latitude0) * cos(longitude0);Cne(2,2) = - sin(latitude0) * sin(longitude0);Cne(2,3) = cos(latitude0);Cne(3,1) = cos(latitude0) * cos(longitude0);Cne(3,2) = cos(latitude0) * sin(longitude0);Cne(3,3) = sin(latitude0);Cen=Cne';%初始四元数的确定《捷联惯导系统》 P151-152 方法本身保证了q1^2+q2^2+q3^2+q4^2=1q(2,1) = sqrt(abs(1.0 + Tnb(1,1) - Tnb(2,2) - Tnb(3,3))) / 2.0; q(3,1) = sqrt(abs(1.0 - Tnb(1,1) + Tnb(2,2) - Tnb(3,3))) / 2.0; q(4,1) = sqrt(abs(1.0 - Tnb(1,1) - Tnb(2,2) + Tnb(3,3))) / 2.0;q(1,1) = sqrt(abs(1.0 - q(2,1) ^2 - q(3,1) ^2 - q(4,1) ^2));% 判断q(1,1)的符号flag_q11=cos(head_err0/2.0)*cos(pitch_err0/2.0)*cos(roll_err0/2.0)-sin(head_err0/2.0)*sin(pitch_err0/2.0)*sin(roll_err0/2.0);if (flag_q11 >0) %此时q(1,1)取正if (Tnb(3,2) < Tnb(2,3))q(2,1) = - q(2,1);endif (Tnb(1,3) < Tnb(3,1))q(3,1) = - q(3,1);endif (Tnb(2,1) < Tnb(1,2))q(4,1) = - q(4,1);endelse %此时q(1,1)取负或0q(1,1) = - q(1,1);if (Tnb(3,2) > Tnb(2,3))q(2,1) = - q(2,1);endif (Tnb(1,3) > Tnb(3,1))q(3,1) = - q(3,1);endif (Tnb(2,1) > Tnb(1,2))q(4,1) = - q(4,1);endend%-------------------------迭代推算用到的参数的初始化------------------------Wiee_e = 0;Wiee_n = 0;Wiee_u = Omega_ie;Wiee = [Wiee_e Wiee_n Wiee_u]'; %地球速率在地球系中的投影东－北－天Lat_err(1)=latitude0;Lon_err(1)=longitude0;High_err(1)=High0;Ve_err(1)=Ve0;Vn_err(1)=Vn0;Vu_err(1)=Vu0;pitch_err(1)=pitch_err0;head_err(1)=head_err0;roll_err(1)=roll_err0;Re=6378137.0;%6378245.0; %地球长轴《惯性导航系统》 P28e=0.0033528106647474807198455286185206; %地球扁率精确值ee=0.00669437999014131699614;%----------------------------迭代推算开始-----------------------------------for i=1:NumofEachData%----------------------------惯性仪表数据的获得------------------------Wibb(1,1)=omigax(i); %指向右机翼方向为x正方向，指向机头方向为y正向，z轴与x轴和y轴构成右手坐标系Wibb(2,1)=omigay(i); %单位：弧度/妙Wibb(3,1)=omigaz(i); %右前上fb(1,1)=fb_x(i); %指向右机翼方向为x正方向，指向机头方向为y正向，z轴与x轴和y轴构成右手坐标系fb(2,1)=fb_y(i); %单位：米/秒2fb(3,1)=fb_z(i); %右前上%--------计算在姿态矩阵和位置矩阵更新时用到的参数------------------ RM=Re*(1.0-2.0*e+3.0*e*Cne(3,3)^2)+High_err(i); %《捷联惯导系统》P233 P235RN=Re*(1.0+e*Cne(3,3)^2)+High_err(i);% RN=Re*(1-ee)/(sqrt(1-ee*sin(Lat_err(i))))^3+High_err(i);% RM=Re/sqrt(1-ee*sin(Lat_err(i)))+High_err(i);%实验当地重力加速度计算《捷联惯导系统》P150 《惯性导航系统》 P35g=g0*((1.0+0.0052884*Cne(3,3)^2)-0.0000059*(1-(1-2*Cne(3,3)^2)^2))*(1 .0-2.0*High_err(i)/Re);tmp_slat=sin(Lat_err(i))*sin(Lat_err(i));Wien = Cne * Wiee; %地球速率在导航系中的投影Wenn(1,1) = -Vn_err(i)/RM;Wenn(2,1) = Ve_err(i)/RN; % <<惯性导航系统>> P45 考虑了地球转动的影响.Wenn(3,1) = Ve_err(i)*tan(Lat_err(i))/RN; %计算Wenn(不太精确)，更新速度和位置矩阵时用Winn=Wien+Wenn;Winb=Tbn*Winn;Wnbb=Wibb-Winb; %姿态速率在姿态更新时用到fn=Tnb*fb; % x-y-z 东－北－天% 速度的更新《捷联惯导系统》 P30 33 东－北－天difVe_err=fn(1,1)+(2*Wien(3,1)+Wenn(3,1))*Vn_err(i)-(2*Wien(2,1)+Wenn (2,1))*Vu_err(i);difVn_err=fn(2,1)-(2*Wien(3,1)+Wenn(3,1))*Ve_err(i)+(2*Wien(1,1)+Wenn (1,1))*Vu_err(i);difVu_err=fn(3,1)+(2*Wien(2,1)+Wenn(2,1))*Ve_err(i)-(2*Wien(1,1)+Wenn (1,1))*Vn_err(i)-g;Ve_err(i+1)=Ve_err(i)+difVe_err*TimeEach;Vn_err(i+1)=Vn_err(i)+difVn_err*TimeEach;Vu_err(i+1)=Vu_err(i)+difVu_err*TimeEach;High_err(i+1)=High_err(i)+Vu_err(i)*TimeEach;% 位置矩阵的实时更新《惯性导航系统》 P190Cne(1,1)=Cne(1,1)+TimeEach*(Wenn(3,1)*Cne(2,1)-Wenn(2,1)*Cne(3,1)); Cne(1,2)=Cne(1,2)+TimeEach*(Wenn(3,1)*Cne(2,2)-Wenn(2,1)*Cne(3,2)); Cne(1,3)=Cne(1,3)+TimeEach*(Wenn(3,1)*Cne(2,3)-Wenn(2,1)*Cne(3,3)); Cne(2,1)=Cne(2,1)+TimeEach*(-Wenn(3,1)*Cne(1,1)+Wenn(1,1)*Cne(3,1)); Cne(2,2)=Cne(2,2)+TimeEach*(-Wenn(3,1)*Cne(1,2)+Wenn(1,1)*Cne(3,2)); Cne(2,3)=Cne(2,3)+TimeEach*(-Wenn(3,1)*Cne(1,3)+Wenn(1,1)*Cne(3,3)); Cne(3,1)=Cne(3,1)+TimeEach*(Wenn(2,1)*Cne(1,1)-Wenn(1,1)*Cne(2,1)); Cne(3,2)=Cne(3,2)+TimeEach*(Wenn(2,1)*Cne(1,2)-Wenn(1,1)*Cne(2,2)); Cne(3,3)=Cne(3,3)+TimeEach*(Wenn(2,1)*Cne(1,3)-Wenn(1,1)*Cne(2,3));% Mat_Wenn(1,1)=0;% Mat_Wenn(1,2)=Wenn(3,1);% Mat_Wenn(1,3)=-Wenn(2,1); %Wenn的反对阵矩阵取负% Mat_Wenn(2,1)=-Wenn(3,1); %这里位置矩阵的及时修正为:dCne/dt=Mat_Wenn*Cne% Mat_Wenn(2,2)=0;% Mat_Wenn(2,3)=Wenn(1,1);% Mat_Wenn(3,1)=Wenn(2,1);% Mat_Wenn(3,2)=-Wenn(1,1);% Mat_Wenn(3,3)=0;%% Mat_Wenn=Mat_Wenn*Cne*TimeEach;% Cne=Cne+Mat_Wenn;Cen=Cne';% 计算经纬度Lat_err(i+1)=asin(Cne(3,3));Lon_err(i+1)=atan(Cne(3,2)/Cne(3,1)); %这是经度的主值if (Cne(3,1) < 0)if (Lon_err(i+1) > 0)Lon_err(i+1) = Lon_err(i+1) - pi;elseLon_err(i+1) = Lon_err(i+1) + pi;endend% 四元数的及时修正《惯性导航系统》 P194% Mat_Wnbb=[ 0, -Wnbb(1,1), -Wnbb(2,1), -Wnbb(3,1); % Wnbb(1,1), 0, Wnbb(3,1), -Wnbb(2,1); % Wnbb(2,1), -Wnbb(3,1), 0, Wnbb(1,1); % Wnbb(3,1), Wnbb(2,1), -Wnbb(1,1), 0];% q=q+Mat_Wnbb*q*TimeEach/2.0;q(1,1)=q(1,1)+TimeEach*(-Wnbb(1,1)*q(2,1)-Wnbb(2,1)*q(3,1)-Wnbb(3,1)* q(4,1))/2.0;q(2,1)=q(2,1)+TimeEach*(Wnbb(1,1)*q(1,1)+Wnbb(3,1)*q(3,1)-Wnbb(2,1)*q (4,1))/2.0;q(3,1)=q(3,1)+TimeEach*(Wnbb(2,1)*q(1,1)-Wnbb(3,1)*q(2,1)+Wnbb(1,1)*q (4,1))/2.0;q(4,1)=q(4,1)+TimeEach*(Wnbb(3,1)*q(1,1)+Wnbb(2,1)*q(2,1)-Wnbb(1,1)*q (3,1))/2.0;% 四元数归一化处理q_norm=sqrt(sum(q.*q));q=q/q_norm;% 计算姿态矩阵 TnbTnb(1,1) = q(1,1) ^2 + q(2,1) ^2 - q(3,1)^2 - q(4,1)^2;Tnb(1,2) = 2.0 * (q(2,1) * q(3,1) - q(1,1) * q(4,1));Tnb(1,3) = 2.0 * (q(2,1) * q(4,1) + q(1,1) * q(3,1));Tnb(2,1) = 2.0 * (q(2,1) * q(3,1) + q(1,1) * q(4,1));Tnb(2,2) = q(1,1)^2 - q(2,1)^2 + q(3,1)^2 - q(4,1)^2;Tnb(2,3) = 2.0 * (q(3,1) * q(4,1) - q(1,1) * q(2,1));Tnb(3,1) = 2.0 * (q(2,1) * q(4,1) - q(1,1) * q(3,1));Tnb(3,2) = 2.0 * (q(3,1) * q(4,1) + q(1,1) * q(2,1));Tnb(3,3) = q(1,1)^2 - q(2,1)^2 - q(3,1)^2 + q(4,1)^2;Tbn=Tnb';flag_pitch=asin(Tnb(3,2));flag_roll=atan(-Tnb(3,1)/Tnb(3,3));flag_head=atan(-Tnb(1,2)/Tnb(2,2));if(Tnb(3,3)<0)if(flag_roll<0)flag_roll=flag_roll+pi;endif(flag_roll>0)flag_roll=flag_roll-pi;endend% 偏航角范围 -180度——180度北偏西为正if(Tnb(2,2)<0)if(flag_head<0)flag_head=flag_head+pi;endif(flag_head>0)flag_head=flag_head-pi;endend% 姿态角更新pitch_err(i+1)=flag_pitch;head_err(i+1)=flag_head;roll_err(i+1)=flag_roll;% 解算完毕由对准结果、陀螺、加表的输出解算出载体的位置、速度、姿态%----------------------计算解算误差------------------ddLat(i)=(Lat_err(i)-latitude(i))*rad_deg; %纬度误差单位：度ddLog(i)=(Lon_err(i)-longitude(i))*rad_deg; %经度误差单位：度ddHigh(i)=High_err(i)-High(i); %高度误差单位：米ddVe(i)=Ve_err(i)-Ve(i);ddVn(i)=Vn_err(i)-Vn(i); % 速度误差单位：米/妙2ddVu(i)=Vu_err(i)-Vu(i);ddpitch(i)=(pitch_err(i)-pitch(i))*rad_deg*3600; %姿态误差单位：度ddhead(i)=(head_err(i)-head(i))*rad_deg*3600;ddroll(i)=(roll_err(i)-roll(i))*rad_deg*3600;endfclose(fid_read);%---------------------------绘图开始--------------------------------- figure(1)plot(Time,ddLog)ylabel('经度误差(度)'),xlabel('时间(秒)');figure(2)plot(Time,ddLat)ylabel('纬度误差(度)'),xlabel('时间(秒)');figure(3)plot(Time,ddHigh);ylabel('高度误差(米)'),xlabel('时间(秒)');figure(4)plot(Time,ddhead)ylabel('偏航角误差(角妙)'),xlabel('时间(秒)'); figure(5)plot(Time,ddpitch)ylabel('俯仰角误差(角妙)'),xlabel('时间(秒'); figure(6)plot(Time,ddroll);ylabel('滚转角误差(角妙)'),xlabel('时间(秒)'); figure(7)plot(Time,ddVe);ylabel('东向速度误差(米/秒)'),xlabel('时间(秒)'); figure(8)plot(Time,ddVn)ylabel('北向速度误差(米/秒)'),xlabel('时间(秒)'); figure(9)plot(Time,ddVu)ylabel('天向速度误差(米/秒)'),xlabel('时间(秒)');%------------------------------绘图结束-------------------------------。

0 cos sin ， Rz sin 0 cos
sin cos 0
0 0 1
cos cos sin sin sin cos cos sin sin cos sin cos T11 T12 T13 Ry Rx Rz cos sin cos cos sin T21 T22 T23 sin cos cos sin sin sin sin cos sin cos cos cos T T T 31 32 33 b 由姿态矩阵 C n 反解飞行器姿态欧拉角：
（5） 速度的计算
t t t t t 0 2iez etz ety 2iey Vxt Vx 0 t t b t t t t 0 2iex etx Vyt 0 Vy Cb f 2iez etz t Vz g Vzt 2 t t 2 t t 0 iey ety iex etx
o o sin 1 T23 , 90 , 90
tg 1
T13 180o , 180o , T33
tg 1
T21 o o , 180 , 180 T 22
图 6 东向北向速度变化曲线
阶段总结：1.学习了平台式和捷联式惯导的惯导解算方法并进行了仿真计算。 2.平台式惯导物理平台时刻跟踪当地水平东北天地理系， 加速计的比 力信息直接投影在导航系中，可直接进行导航速度和位置解算。载体的姿态可直 接从平台框架直接得出；而捷联式惯导用数学平台取代实体的物理平台，通过求

设两个四元数为
Λ = λ0 + λ1i + λ2 j + λ3k M = m0 + m1i + m2 j + m3k
两者和
Λ + M = ( λ0 + m0 ) + ( λ1 + m1 ) i + ( λ2 + m2 ) j + ( λ3 + m3 ) k Λ − M = ( λ0 − m0 ) + ( λ1 − m1 ) i + ( λ2 − m2 ) j + ( λ3 − m3 ) k
弊端：数值积分时进行超越函数运算 纵摇角为90o度，微分方程出现奇点
捷联姿态算法
方向余弦法
方向余弦姿态矩阵微分方程的形式为
b b &b Cn = Cn ωnb
写成矩阵的形式有
& & C11 C12 & & C21 C22 & & C31 C32 & C13 C11 C12 & C23 = C21 C22 & C33 C31 C32 C13 0 b C23 ωnbz b C33 −ωnby
则两者的乘积
N = M ⊗ Λ = n0 + n1i + n2 j + n3 k
n0 λ0 n λ 1 = 1 n2 λ2 n3 λ3
−λ1
−λ2
λ0 −λ3 λ2
λ3 λ0 −λ1
−λ3 m0 −λ2 m1 λ1 m2 λ0 m3
惯 性 导 航 原 理
第五讲

;no dirg ;)')s/m( /yticoleV'(lebaly ;)')sm01( / emiT'(lebalx ;)'noitaluclaC yticoleV'(eltit ;)'Z :etanidrooC noitagivaN','Y :etanidrooC noitagivaN','X :etanidrooC noitagivaN'(dnegel )yticolev,mt(tolp erugif ;no dirg ;)')eerged( /legnA'(lebaly ;)')sm01( / emiT'(lebalx ;)'noitaluclaC erutseG'(eltit ;)'legnA waY','legnA lloR','legnA hctiP'(dnegel ;)'.-b',):,3(ELGNA,mt,'.g',):,2(ELGNA,mt,'-r',):,1(ELGNA,mt(tolp ;ip/081*ELGNA=ELGNA erugif dne
值比度速加/压电的向轴个三计度速加% ?=za_orez% ?=ya_orez% ?=xa_orez% 压电点零的向轴个三计度速加%
*********************************************************% % 值度速加的应相为换转据数将% 3*4201/cca_atad=cca_atad ;]z_cca;y_cca;x_cca[=cca_atad 压电拟模为换转字数度速加% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ;)'txt.zetarelecca'(daol=z_cca ;)'txt.yetarelecca'(daol= y_cca ;)'txt.xetarelecca'(daol= x_cca ;4201/3=etar_cca% **************据数的计度速加入读****************% ;)'txt.zoryg'(daol=z_oryg ;)'txt.yoryg'(daol=y_oryg ;)'txt.xoryg'(daol=x_oryg 据数的仪螺陀入读**********% 据数入读****************************% ;']g-,0,0[=G 2^s/m8.9 位单%;24108.9=g 量矢度速加的生产力重% ;lla raelc ;lla esolc

用四元数法的捷联惯性导航姿态解算程序close all;clear all;%重力产生的加速度矢量g=9.80142;%单位9.8m/s^2G=[0,0,-g]';%****************************读入数据%**********读入陀螺仪的数据gyro_x=load('gyrox.txt');gyro_y=load('gyroy.txt');gyro_z=load('gyroz.txt');%****************读入加速度计的数据**************%acc_rate=3/1024;acc_x =load('acceleratex.txt');acc_y =load('acceleratey.txt');acc_z=load('acceleratez.txt'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %加速度数字转换为模拟电压data_acc=[acc_x;acc_y;acc_z];data_acc=data_acc/1024*3%将数据转换为相应的加速度值%%*********************************************************%加速度计三个轴向的零点电压%zero_ax=?%zero_ay=?%zero_az=?%加速度计三个轴向的电压/加速度比值%rate_ax=? %单位是m/s^2/V%rate_ay=?%rate_az=?%acc_x=acc_x*acc_rate;%acc_y=acc_y*acc_rate;%acc_z=acc_z*acc_rate;aver_acc_x=mean(acc_x)aver_acc_y=mean(acc_y)aver_acc_z=mean(acc_z)%采样时间dtime=0.01;tm=0:dtime:0.01* (size(gyro_x,2)-1);%个数numn_point=size(gyro_x,2);%图1figureplot(tm,data_acc(1,:),'-',tm,data_acc(2,:),'.',tm,data_acc(3,:),'-.'); title('加速度计的采样曲线');legend('x_ACC','Y_ACC','Z_ACC');xlabel('Time / (10ms)');ylabel('Accelerate/ (m/s'')');grid on;%plot(tm,acc_x,'-',tm,acc_y,'.',tm,acc_z,'-.');%title('加速度的计的采样曲线'):%对采样曲线进行低通滤波a=[1,2,4,2,1];%gyro_x=filter(a/sum(a),1,gyro_x);%gyro_y=filter(a/sum(a),1,gyro_y);%gyro_z=filter(a/sum(a),1,gyro_z);%比例变换gyro_x=gyro_x/1024*3/0.6;gyro_y=gyro_y/1024*3/0.6;gyro_z=gyro_z/1024*3/0.6;%零点电压--陀螺仪，取前80个数的平均电压zero_gx=sum(gyro_x(1:80))/80zero_gy=sum(gyro_y(1:80))/80zero_gz=sum(gyro_z(1:80))/80%减去零点gyro_x=(gyro_x-zero_gx)/0.0125/180*pi; gyro_y=(gyro_y-zero_gy)/0.0125/180*pi; gyro_z=(gyro_z-zero_gz)/0.0125/180*pi;%gyro_x=(gyro_x-2.5)/0.0125/180*pi;%gyro_y=(gyro_y-2.5)/0.0125/180*pi;%gyro_z=(gyro_z-2.5)/0.0125/180*pi;%测试数据accelerate=zeros(3,n_point);accelerate(1,1:100)=10;accelerate(1,101:200)=-10;accelerate(1,201:300)=0;%陀螺仪数据gyro_x=zeros(1,n_point);gyro_y=zeros(1,n_point);gyro_z=zeros(1,n_point);gyro_z(1:100)=pi/3;gyro_z(101:200)=-pi/3;%重力轴始终有加速度accelerate(3,:)=accelerate(3,:)+9.8;figureplot(tm,accelerate(1,:),'-',tm,accelerate(2,:),'.',tm,accelerate(3,:),'-.');title('加速度计的采样曲线');legend('x_ACC','Y_ACC','Z_ACC');xlabel('Time / (10ms)');ylabel('Accelerate/ (m/s'')');grid on;%画出陀螺仪的采样曲线figureplot(tm,gyro_x,'r-',tm,gyro_y,'g.',tm,gyro_z,'b-.');title('陀螺仪的采样曲线');legend('x_Gyro','Y_Gyro','Z_Gyro');xlabel('Time / (10ms)');ylabel('Angel_rate/ (degree/s)');grid on;%size(gyro_x)%size(gyro_y)%size(gyro_z)data_gyro=[gyro_x;gyro_y;gyro_z];%转移矩阵--即方向余弦矩阵T=eye(3); %T是3*3的单位矩阵,初始转移矩阵%四元数矩阵，存储每步更新之后的四元数，方便以后绘图Q=zeros(4,n_point);%四元数的初始值确定，假定一开始导航坐标系与载体坐标系是重合的，因此方向余弦矩阵，是单位矩阵，利用它们之间的关系确定四元数的初始值。

b nbz
CC
RP
0 0
H.
C
R
P0. 0
R0.
0
2010-03-19
整理课件
Yb'' Yb Yb' H
N
17
b nbx b nby
cos
0
R
b nbz
sin R
P. cos R 0
R.
0
1
H. sin R 0
0 sin R cos P P.
1
sin P R.
kb jn
ib j
kn k
b n kb kn
2010-03-19
21
整理课件
固定矢量的坐标变换
矢量的坐标变换
固定矢量的坐标变换是一个在空间大小和方向都不 变的矢量在两个不同方位的坐标系轴向分量之间的变 换关系，也即同一个矢量在两个不同的坐标系轴向投 影之间的变换关系。
旋转矢量的坐标变换
是指一个矢量大小不变，但在方向上转动了一个位 置，这个矢量转动前和转动后在同一个坐标系轴向 分量之间的变换关系。
2010-03-19
22
整理课件
固定矢量的坐标变换
一个矢量r，写成载体坐标系轴向分量形式：
r X bib Yb jb Z b kb r bT b
rb
X Y
b
b
ib j
b
Zb
kbb
b:载体坐标系 n:地理坐标系
同一个矢量r，如果写成地理坐标系轴向分量形式：
r X nin Yn jn Z n k n r nT n 由于r是同一个矢量，故
惯性导航系统原理
1
整理课件
3 捷联式惯导系统

z
0 x
y
姿态矩阵计算方法
3 四元数微分方程及其解
1 qk 1 exp( )qk Mq k qk rk 2
rk ac as x i as y j as z k
(0.5 )2 (0.5 )4 ac cos( ) 1 2 2! 4!
tk b nb
t
v f dt
b tk
t
速度微分方程求解
2 速度积分算法
v f dt [2ω ωen ] v dt g1 dt
n n 0 0 n ie n n n 0
t
t
t
v
n k 1
v u gt
n k n
b
v
n l 1
[ I 2Ωie tl Ωen tl ]v
v f dt
b
t
1 tk 1 αk 1 α ωdt 2 tk
α ωdt
tk
tk
t
捷联导航算法
捷联惯导误差分析
主要误差源

1 初始对准误差
2 惯性器件误差
3 计算误差
误差传播方程：误差源与系统性能的相互关系的 方程，主要用来阐述各种类型误差随时间的传播 过程。
捷联惯导误差分析
sin( / 2) (0.5 ) (0.5 ) as 0.5(1 ) 3! 5!
2 4
姿态矩阵计算方法
4 圆锥误差
σ αk 1 αk 1
α ωdt
tk t
1 tk 1 αk 1 α ωdt 2 tk
姿态矩阵计算方法
5 姿态角计算
b nbz

捷联系统的四元数法姿态算法算法输入：物体的初始姿态，3轴陀螺仪不同时刻的Yaw，Pitch，Roll的角速度；算法输出：物体的当前姿态。

具体算法：1. 初始姿态的四元数(w,x,y,z)=(1,0,0,0) 命名为A2. 读取3轴陀螺仪当前时刻的Yaw,Pitch,Roll角速度，乘以上次计算以来的间隔时间，得到上一时刻以来(Yaw,Pitch,Roll)的变化量，命名为欧拉角b3. b是Tait-Bryan angle定义的欧拉角，将其转为四元数B4. A=A×B，做四元数乘法，即可得到当前姿态对应的新的四元数A5.重复2～4部，即可连续更新姿态6.将四元数A重新转换为Tait-Bryan angle形式的欧拉角a，就可以以直观的形式查看当前姿态核心算法1，欧拉角转四元数void Quaternion::FromEulerAngle(const EulerAngle &ea){float fCosHRoll = cos(ea.fRoll * .5f);float fSinHRoll = sin(ea.fRoll * .5f);float fCosHPitch = cos(ea.fPitch * .5f);float fSinHPitch = sin(ea.fPitch * .5f);float fCosHYaw = cos(ea.fYaw * .5f);float fSinHYaw = sin(ea.fYaw * .5f);w = fCosHRoll * fCosHPitch * fCosHYaw + fSinHRoll * fSinHPitch * fSinHYaw;x = fCosHRoll * fSinHPitch * fCosHYaw + fSinHRoll * fCosHPitch * fSinHYaw;y = fCosHRoll * fCosHPitch * fSinHYaw - fSinHRoll * fSinHPitch * fCosHYaw;z = fSinHRoll * fCosHPitch * fCosHYaw - fCosHRoll * fSinHPitch * fSinHYaw;}核心算法2，四元数转欧拉角EulerAngle Quaternion::ToEulerAngle() const{EulerAngle ea;ea.fRoll = atan2(2 * (w * z + x * y) , 1 - 2 * (z * z + x * x));ea.fPitch = asin(CLAMP(2 * (w * x - y * z) , -1.0f , 1.0f));ea.fYaw = atan2(2 * (w * y + z * x) , 1 - 2 * (x * x + y * y));return ea;}核心算法3，四元数乘法Quaternion Quaternion::Multiply(const Quaternion &b){Quaternion c;c.w=w*b.w -x*b.x -y*b.y -z*b.z;c.x=w*b.x +x*b.w +y*b.z -z*b.y;c.y=w*b.y -x*b.z +y*b.w +z*b.x;c.z=w*b.z +x*b.y -y*b.x +z*b.w;c.Normalize();return c;}次要的规范化算法：void Quaternion::Normalize(){float s=getS();w/=s;x/=s;y/=s;z/=s;}float Quaternion::getS(){return sqrt(w*w+x*x+y*y+z*z);}我的loop函数，算法的集成部分：Quaternion nowQ;void loop() {int intx, inty,intz;float pitch,roll,yaw;gyro.ReadGyroOutCalibrated_Radian(&pitch, &roll, &yaw); EulerAngle dt;dt.fRoll=roll;dt.fPitch=pitch;dt.fYaw=-yaw;Quaternion dQ;dQ.FromEulerAngle(dt);nowQ=nowQ.Multiply(dQ);count++;if (count>1000){EulerAngle nowG=nowQ.ToEulerAngle();Serial.print(nowG.fRoll/3.1415926535*180,11);//横滚Serial.print(",");Serial.print(nowG.fPitch/3.1415926535*180,11);//俯仰Serial.print(",");Serial.print(nowG.fYaw/3.1415926535*180,11);//偏航Serial.print(",");Serial.print(nowQ.getS(),11);//偏航Serial.println();count=0;}}四元数与欧拉角之间的转换在3D图形学中，最常用的旋转表示方法便是四元数和欧拉角，比起矩阵来具有节省存储空间和方便插值的优点。

捷联式惯性导航积分算法设计上篇：姿态算法Paul G. SavageStrapdown Associates, Inc., Maple Plain, Minnesota55359摘要：本论文分上下两篇，用于给现代捷联惯导系统的主要软件算法设计提供一个严密的综合方法：将角速率积分成姿态角，将加速度变换或积分成速度以及将速度积分成位置。

该算法是用两速修正法构成的，而两速修正法是具有一定创新程度的新颖算法，是为姿态修正而开发出来的，在姿态修正中，以中速运用精密解读方程去校正积分参数(姿态、速度或位置>，其输入是由在参数修正(姿态锥化修正、速度摇橹修正以及高分辨率位置螺旋修正>时间间隔内计算运动角速度和加速度的高速算法提供的。

该设计方法考虑了通过捷联系统惯性传感器对角速度或比力加速度所进行的测量以及用于姿态基准和矢量速度积分的导航系旋转问题。

本论文上篇定义了捷联惯导积分函数的总体设计要求，并开发出了用于姿态修正算法的方向余弦法和四元数法；下篇着重讨论速度和位置积分算法的设计。

尽管上下两篇讨论中常常涉及到基本的惯性导航概念，然而，本论文是为那些已对基础惯导概念很熟悉的实际工作者而写的。

专门用语：123,,,A A A A =任意坐标系12A A C =将矢量从2A 坐标系投影到1A 坐标系的方向余弦矩阵 I =单位矩阵12A A q =从2A 坐标系投影到1A 坐标系的旋转矢量所构成的姿态变化四元数 1*2A A q =12A A q 的共轭四元数，它的第1项与12A A q 的首项相同，余下的2~4项与 12A A q 的互为相反数 1q =单位四元数，它的第1项为1，其余3项为0V =无具体坐标系定义的矢量A V =列向量，它的各项元素等于矢量V 在坐标系A 的各轴上的投影A V ⨯（） =向量AV 的反对称<或交叉积）形式，代表如下矩阵：000ZA YA ZA XA YAXAV V V V V V -⎡⎤⎢⎥-⎢⎥⎢⎥-⎣⎦其中：XA V ，YA V ，ZA V 是AV 的分量，AV ⨯（）与A 系矢量的矩阵乘积等于AV 与该矢量的叉积A q V =与A V 等量的四元数矢量，0A V ⎡⎤⎢⎥⎣⎦2A ω1A =2A 坐标系相对于1A 坐标系的角速率，当1A 为惯性系<I 系）时，2A ω1A 是由安装在2A 坐标系上的角速率传感器所测到的角速率1.概论惯性导航是通过对速度积分得到位置并对总加速度积分得到速度的过程。

收稿日期 : 2003207225 基金项目 : 国防预研资助项目 (18 YXGFK D118) 作者简介 : 刘 危 (1974 - ) ,男 ,江西都昌人 ,博士生 , bullkill @21cn. com.
© 1994-2007 China Academic Journal Electronic Publishing House. All rights reserved.
2005 年 1 月 第 31 卷 第 1 期
北京航空航天大学学报 Journal of Beijing University of Aeronautics and Astronautics
January 2005 Vol. 31 No11
捷联惯性导航系统的姿态算法
刘 危 解旭辉 李圣怡
( 19)
Δ T ,则近似中的算法误差为 若令 λ= Ω Δ δ ^Φ - δ Φ ≈δ Φ =
N
- 2 ( p + N)
2
p + N +1
∏
( p + N) ! 2 ( p + N ) +1 λ ab p + N +1 ( 2 k - 1) k =1
( 13)
A =
… … … … … … … A i ,1 … A i , N - 1 A i , N … A i , N + p - 1 A i , N + p … … … … …

第 1 期 刘 危等 : 捷联惯性导航系统的姿态算法
47
k1 K( N + p - 1) × 1 = k2
d1
其中 , A ij ( j ≠( N + p) ) 的定义如式 ( 11) ; 而 A i , N + P 的定义如下 :