卡尔曼滤波与组合导航》课程实验报告
实验 捷联惯导 /GPS 组合导航系统静态导航实验
实验序号 3 姓名 陈星宇
系院专业
17
班级
ZY11172
学号
ZY1117212
日期
2012-5-15
指导教师
宫晓琳
成绩
、实验目的
① 掌握捷联惯导 /GPS 组合导航系统的构成和基本工作原理； ②掌握采用卡尔曼滤波方法进行捷联惯导 /GPS 组合的基本原理； ③掌握捷联惯导 /GPS 组合导航系统静态性能；
④了解捷联惯导 /GPS 组合导航静态时的系统状态可观测性；
、实验原理
（ 1）系统方程 X FX GW
系统噪声矢量由陀螺仪和加速度计的随机误差组成，表达式为：
2）量测方程
和 H 分别为捷联解算与 GPS 的东向速度、北向速度、天向速度、纬度、经度和高度之 差；量测矩阵 H H V H P T
，H P 03 6 diag R M H, （R N H ）cos L, 036 ，
H V 033 diag 1, 1, 1 039
，v v V E
v V
N
v V U v L v v H 为量测噪声。 量测噪声
v E v N
T
v
U
L h
x y z x y z
其中， E 、 N 、 U 为数学平台失准角；
v E 、 v N 、 v U 分别为载体的东向、北向和天向速度误差；
L 、 、 h 分别为纬度误差、经度误差和高度误差；
x 、 y 、 z 、 x 、 y 、 z 分别为陀螺随
机常值漂移和加速度计随机常值零偏。（下
标 系统的噪声转移矩阵 G 为：
E 、N 、 U 分别代表东、北、天）
C b
n
3 3 0
9 3
3 3 C n C b
9 3
15 6
系统的状态转移矩阵 w w w w
F 组成内容为：
w
z
F
06N
9
F S
F
M
，其中 F N 中非零元素为可由惯导误差模型获得。 F S
C b
n
3 3 0
3 3
3 3
C b
n
3 3
96
量测变量 z V E V N V U L H ，
， V E 、 V N 、 V U 、 L 、
X
U
方
差阵R 根据GPS的位置、速度噪声水平选取。
(3)卡尔曼滤波方程
状态一步预测：X?k/k 1 k,k 1X?k 1
状态估计：X?k X?k/k 1 K k(Z k H k X?k/k 1)
滤波增益：K k P k/k 1H k T(H k P k/k 1H k T R k) 1 2 一步预测均方误差：P k/k 1 k,k
三、实验内容及步骤
1、实验内容
① 捷联惯导/GPS 组合导航系统静态导航实验；
2、实验步骤
2实验准备(由实验教师操作)
① 将IMU 安装在大理石实验台上，确认IMU 的安装基准面靠在大理石实验平台上的方位基准尺上；
②将IMU 与导航计算机、导航计算机与稳压电源、导航计算机与监控计算机、GPS 接收机与导航计算机、GPS天线与GPS接收机、GPS 接收机与GPS 电池之间的连接线正确连接；
③打开GPS 信号转发器；
④打开监控计算机中的监控软件；
⑤打开稳压电源开关，确认稳压电源输出为28V ；
⑥打开捷联惯导/GPS 组合实验系统电源，实验系统开始启动，注意30s 内严禁动
IMU；
⑦打开GPS接收机电源，确认通过信号转发器可以接收到 4 颗以上卫星；
⑧将监控软件设置为“准备”状态，准备时间 5 分钟；
⑨准备过程中由实验教师介绍捷联惯导/GPS 组合实验系统的组成、工作原理；
1P k 1 k T,k 1 Q k 1 估计均方误差：P k (I K k H k)P k/k 1
2) 捷联惯导 /GPS 组合导航系统静态导航实验
① 实验系统准备 5 分钟之后， 通过监控软件， 将实验系统设置为 “组合导航” 状态；
② 记录 IMU 的原始输出，即角增量和比力信息；
③ 记录数据过程中，由实验教师介绍采用卡尔曼滤波方法进行捷联惯导 /GPS 组合导 航的基本原理；
④ 记录 IMU 数据 5 分钟后，停止记录数据，利用监控计算机中的捷联惯
导
/GPS 组 合导航软件进行静态导航解算，并显示静态导航结果； 四、实验结果
与分析 1、SINS/GPS 组合导航后得纬度曲线和 GPS 导航纬度曲线对比如下图
）度
（度34.95 34.9
34.85
34.8 34.75 34.7 34.65 34.6 34.55 34.5
34.45 0 1 2 3 4 5 6 7 8
4 x 10
2、SINS/GPS 组合导航后得经度曲线和 GPS 导航经度曲线对
比如下图
GPS 纬 度 组合导航
后纬度 )
度 (
度 经
110.1 110 109.9 109.8 109.7
109.6 109.5 109.4 109.3 109.2 109.1
GPS 经 度 组合导航
后经度
4
x 10
3、SINS/GPS 组合导航后得高度曲线和 GPS 导航高度曲线对
比如下图
3500
3000 GPS 高 度 组合导航后
高度
2500 2000
（度
高
1500
高 1000
500
8
4
x 10
4、 SINS/GPS 组合导航后得东向速度曲线和
GPS 导航东向速度曲线对比如
下图
80 60 GPS 东 向 速 度 组 合 导 航 后 东 向 速 度
40 20 0 20
40
）
s /m
（度速向东
-60
1 2 3 4 5 6 7 8 4
x 10
5、 SINS/GPS 组合导航后得北向速度曲线和
GPS 导航北向速度曲线对比如
下图
60
40
20
-20
-40
-60
GPS 北 向 速 度 组合导航后北向速度
8
4
x 10 4
6、 SINS/GPS 组合导航后得天向速度曲线和
GPS 导航天向速度曲线对比如
下图
5
67
GPS 天 向 速 度 组 合 导 航 后 天 向 速 度
6 4 2 0 ）
s /m
（度速向天
24
4
x 10
7、 SINS/GPS 组合导航后航向角曲线、俯仰角曲线和横滚角
曲线如下图
200 150
100
50
-50
-100
-150
-200
组合导航后航向角 组合导航后俯仰角 组合导航后横滚角
012 4
x 10
8、纯惯性导航轨迹、 GPS 导航轨迹和 SINS/GPS 组合导航轨迹对比如下图
110.1 110
109.9
109.8 109.7 度 109.6 经
109.5 109.4 109.3 109.2
109.1
34.45 34.5 34.55 34.6 34.65 34.7 34.75 34.8 34.85 34.9 34.95 纬度 纯惯性导航轨迹 GPS 导
航 轨 迹 组合导航导航轨迹
9、平台失准角的估计均方差曲线如下图
0.02 度 0.01 0
0.02 度 0.01 0
0 1 2
3北 向 失 准4 角 方 差 5 6
7
8
4
x 10
东向失准角方差
0 1 2 3 4 5 6 7
8 0.5
天向失准角方差
4
x 10
0 1 2 3 4 5 6 7 8 4
x 10
10、速度和位置的估计均方差曲线如下图
东向速度误差方差
0.01
/s m
0.005
0 2 4 6 8 4
x 10 天向速度误差方差 0.01
/s m 0.005
0 02468 4 x
10 经度误差方差
0.1
m 0.05
0 02468
4 x 10
北 向 速 度 误 差 方 差
0.01
/s m
0.005
0 2 4 6 8
4
x 10 纬
度误差方差
0.1
m 0.05
0 2 4 6 8
4
x 10 高
度误差方差
0.2
m 0.1
0 2 4 6 8
4
x 10
11、陀螺漂移的估计均方差曲线和加速度计偏置的估计均方差曲线如下图
结果分析：从仿真结果可以看出，滤波之后载体的位置和速度比
GPS 输出的位置和速度精度高，
载体姿态在滤波校正后结果较好， INS/GPS 组合导航有效地抑制了纯惯性导航的发散，也有效地去除 了 GPS 量测输出的噪声。 东北天方向的速度误差均能够估计出来， 天向陀螺漂移估计不出来， 在动基 座情况下，东向和北向加计零偏、天向平台失准角以及东向和北向陀螺漂移均变得可观测，收敛性变 好。可见， INS/GPS 是一种较为理想的组合导航方式。
五、源程序
clear; clc; % 载入数据
卡尔曼 \IMU.dat'); 卡尔曼 \GPS.dat');
%%%%%%%%%%定义常数 e=1/298.3; re=6378245;
/时
小度
度
0.1 0.1 /时
小度
度
0.1 东向陀螺漂移方差
0.05
/时
小
度
度
x 10
4
天向陀螺漂移方差 北向陀螺漂移方差
0.05
x 10
4
东向加计偏置方差
50
g
μ
0.05
50
g
μ
g
μ
4
x 10
4
北向加计偏置方差
4
x 10 4
4 x 10
4
北向加计偏置方差
4
x 10 4
A (6σ51)s c u ①NI G O1)① >9 G ol )s c
u ① Z
IΓ1L H ?≡≡J M s ?≡t ? 世
≡1*
沢抿&&&&&&
&&
&&&
±f +(0v (ωDU
ω?*
c?① *0+
L)a)」H
Q
±f +(0v (ωDU
ω?丄)
05
」
H X M
(O H H (S )P O E )
七
?8L *d ?E
e 6-08L *d '2o -08L *?e s
①
兰
H C -￥Pn ≡ω
EOEddUO--OFdAe=
H C-Duomsod
π>A >x >H U D A l o
o α)> <^≡t ?
鋼曰，翌?旨旨^*<吟擊型>
&&&&&&&&&&&&
PU ①
PU
①
?+Llelue6Helue6
φω-φ
?ILIeEe6Helue6
OAGL)
q u o
七
φω-φ
=IeIUe6Helue6
o ^e o
)q u o 4φs φ
PU
①
0ddelue6
φω-φ
ZAdHelUeC
OAGL)
q u o
七
O H H o o )q u o
七
PU
①
fai_x
Fn(5,6)=-vy/Ry;
kesai_kf=-pi/2;
');ylabel('经度(度) '); figure(3) plot(1:72001,GPS(:,5));hold on
plot(t,position_kf(:,3),'r');hold on legend('GPS高度','组合导航后高度');ylabel('高度( m)'); figure(4) plot(1:72001,GPS(:,6));hold on
plot(t,velocity_kf(:,1),'r');hold on legend('GPS东向速度','组合导航后东向速度');ylabel('东向速度( m/s)'); figure(5)
plot(1:72001,GPS(:,7));hold on plot(t,velocity_kf(:,2),'r');hold on legend('GPS北向速度','组合导航后北向速度');ylabel('北向速度
( m/s)'); figure(6) plot(1:72001,GPS(:,8));hold on
plot(t,velocity_kf(:,3),'r');hold on legend('GPS天向速度','组合导航后天向速度');ylabel('天向速度( m/s)') figure(7)
plot(t,attitude_kf(:,1));hold on plot(t,attitude_kf(:,2),'r');hold on
plot(t,attitude_kf(:,3),'g');hold on legend('组合导航后航向角','组合导航后俯仰角','组合导航后横滚角');ylabel(' 度'); figure(8)
plot(position(:,1),position(:,2));hold on plot(GPS(:,3),GPS(:,4),'r');hold on plot(position_kf(:,1),position_kf(:,2),'g');hold on legend('纯惯性导航轨迹','GPS导航轨迹','组合导航导航轨迹');xlabel('纬度'); ylabel('经度');
figure(9) subplot(3,1,1);
plot(t,p_kf(:,1));title(' 东向失准角方差');ylabel(' 度');
subplot(3,1,2);
plot(t,p_kf(:,2));title(' 北向失准角方差');ylabel(' 度');
subplot(3,1,3);
plot(t,p_kf(:,3));title(' 天向失准角方差');ylabel(' 度'); figure(10) subplot(3,2,1); plot(t,p_kf(:,4));title(' 东向速度误差方差');ylabel('m/s');
subplot(3,2,2); plot(t,p_kf(:,5));title(' 北向速度误差方差');ylabel('m/s');
subplot(3,2,3); plot(t,p_kf(:,6));title(' 天向速度误差方差');ylabel('m/s');
subplot(3,2,4); plot(t,p_kf(:,7));title(' 纬度误差方差');ylabel('m'); subplot(3,2,5); plot(t,p_kf(:,8));title(' 经度误差方差');ylabel('m');
subplot(3,2,1);
plot(t,p_kf(:,10));title(' 东向陀螺漂移方差');ylabel(' 度/小时'); subplot(3,2,2);
plot(t,p_kf(:,11));title(' 北向陀螺漂移方差');ylabel('度/小时'); subplot(3,2,3);
plot(t,p_kf(:,12));title(' 天向陀螺漂移方差');ylabel(' 度/小时'); subplot(3,2,4);
plot(t,p_kf(:,13));title(' 东向加计偏置方差');ylabel(' μg'); subplot(3,2,5);
plot(t,p_kf(:,14));title(' 北向加计偏置方差');ylabel(' μg'); subplot(3,2,6);
plot(t,p_kf(:,15));title(' 北向加计偏置方差');ylabel(' μg');
卡尔曼滤波计算举例 ?计算举例 ?卡尔曼滤波器特性
假设有一个标量系统，信号与观测模型为 [1][][]x k ax k n k +=+[][][] z k x k w k =+其中a 为常数，n [k ]和w [k ]是不相关的零均值白噪声，方差分别为和。 系统的起始变量x [0]为随机变量，其均值为零，方差为。2n σ2 σ[0]x P （1）求估计x [k ]的卡尔曼滤波算法；（2）当时的卡尔曼滤波增益和滤波误差方差。 22 0.9,1,10,[0]10 n x a P =σ=σ==1. 计算举例
根据卡尔曼算法，预测方程为: ??[/1][1/1]x k k ax k k -=--预测误差方差为: 2 2 [/1][1/1]x x n P k k a P k k -=--+σ 卡尔曼增益为: () 1 22 22 22 [][/1][/1][1/1][1/1]x x x n x n K k P k k P k k a P k k a P k k -=--+σ --+σ=--+σ+σ ???[/][/1][]([][/1])??[1/1][]([][1/1])?(1[])[1/1][][]x k k x k k K k z k x k k ax k k K k z k ax k k a K k x k k K k z k =-+--=--+---=---+滤波方程:
()() 2 2222222 222 22 [/](1[])[/1] [1/1]1[1/1][1/1][1/1][1/1]x x x n x n x n x n x n P k k K k P k k a P k k a P k k a P k k a P k k a P k k =--??--+σ=---+σ ?--+σ+σ??σ--+σ = --+σ+σ 滤波误差方差 起始：?[0/0]0x =[0/0][0] x x P P =
Kalman_Filter(float Gyro,float Accel) { Angle+=(Gyro - Q_bias) * dt; Pdot[0]=Q_angle - PP[0][1] - PP[1][0]; Pdot[1]= - PP[1][1]; Pdot[2]= - PP[1][1]; Pdot[3]=Q_gyro; PP[0][0] += Pdot[0] * dt; PP[0][1] += Pdot[1] * dt; PP[1][0] += Pdot[2] * dt; PP[1][1] += Pdot[3] * dt; Angle_err = Accel - Angle; PCt_0 = C_0 * PP[0][0]; PCt_1 = C_0 * PP[1][0]; E = R_angle + C_0 * PCt_0; K_0 = PCt_0 / E; K_1 = PCt_1 / E; t_0 = PCt_0; t_1 = C_0 * PP[0][1]; PP[0][0] -= K_0 * t_0; PP[0][1] -= K_0 * t_1; PP[1][0] -= K_1 * t_0; PP[1][1] -= K_1 * t_1; Angle += K_0 * Angle_err; Q_bias += K_1 * Angle_err; Gyro_x = Gyro - Q_bias; } 首先是卡尔曼滤波的5个方程： -=--+（1）先验估计 X k k AX k k Bu k (|1)(1|1)() -=--+（2）协方差矩阵的预测(|1)(1|1)' P k k AP k k A Q
卡尔曼滤波简介及其算法实现代码 卡尔曼滤波算法实现代码（C，C＋＋分别实现） 卡尔曼滤波器简介 近来发现有些问题很多人都很感兴趣。所以在这里希望能尽自己能力跟大家讨论一些力所能及的算法。现在先讨论一下卡尔曼滤波器，如果时间和能力允许，我还希望能够写写其他的算法，例如遗传算法，傅立叶变换，数字滤波，神经网络，图像处理等等。 因为这里不能写复杂的数学公式，所以也只能形象的描述。希望如果哪位是这方面的专家，欢迎讨论更正。 卡尔曼滤波器– Kalman Filter 1．什么是卡尔曼滤波器 （What is the Kalman Filter?） 在学习卡尔曼滤波器之前，首先看看为什么叫“卡尔曼”。跟其他著名的理论（例如傅立叶变换，泰勒级数等等）一样，卡尔曼也是一个人的名字，而跟他们不同的是，他是个现代人！ 卡尔曼全名Rudolf Emil Kalman，匈牙利数学家，1930年出生于匈牙利首都布达佩斯。1953，1954年于麻省理工学院分别获得电机工程学士及硕士学位。1957年于哥伦比亚大学获得博士学位。我们现在要学习的卡尔曼滤波器，正是源于他的博士论文和1960年发表的论文《A New Approach to Linear Filtering and Prediction Problems》（线性滤波与预测问题的新方法）。如果对这编论文有兴趣，可以到这里的地址下载： https://www.doczj.com/doc/3318903981.html,/~welch/media/pdf/Kalman1960.pdf。 简单来说，卡尔曼滤波器是一个“optimal recursive data processing algorithm（最优化自回归数据处理算法）”。对于解决很大部分的问题，他是最优，效率最高甚至是最有用的。他的广泛应用已经超过30年，包括机器人导航，控制，传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理，例如头脸识别，图像分割，图像边缘检测等等。 2．卡尔曼滤波器的介绍 （Introduction to the Kalman Filter） 为了可以更加容易的理解卡尔曼滤波器，这里会应用形象的描述方法来讲解，而不是像大多数参考书那样罗列一大堆的数学公式和数学符号。但是，他的5条公式是其核心内容。结合现代的计算机，其实卡尔曼的程序相当的简单，只要你理解了他的那5条公式。 在介绍他的5条公式之前，先让我们来根据下面的例子一步一步的探索。 假设我们要研究的对象是一个房间的温度。根据你的经验判断，这个房间的温度是恒定的，也就
文章编号:1008-8652(2005)01-001-004 联合卡尔曼滤波器在数据融合中的应用 胡宏灿1,2 郭 立1 朱俊株1 (1.中国科学技术大学 合肥 230026; 2.海军大连舰艇学院 大连 116018) 【摘要】 介绍多传感器数据融合中联合卡尔曼滤波器的设计步骤,并将此方法用于舰船组合导航系统,计算机仿真和理论分析表明,该滤波器可以做到全局最优,其结构遵循信息分配原则,提高了系统的数值稳定性和容错性,减小了数据传输的工作量与计算量,便于计算机实现,能够满足组合导航系统需要。 关键词:组合导航系统;数据融合;联合卡尔曼滤波 中图分类号:T P391.7 文献标识码:A The Application of Federal Kalman Filter in Data Fusion System Hu Ho ng can1,2 Guo Li1 Zhu Junzhu1 (1.University of Science and T echnology of China H ef ei230026; 2.Dalian N av al Vessels A cademy Dalian116018) Abstract:A new design o f Kalman filter based on data fusion is presented in the paper.Fistly,the fr ame Kalman filter is intr oduced.T hen,the algo rithm is given.T he simulatio n results show that the metho d is useful in integr ated navigation sy stem because it can impr ove accur acy and r eliability,and it has hig h fault-tolerant ability. Keywords:integ rated nav igatio n sy stem;data fusio n;feder al Kalman filter 1 引言 数据融合技术是近年来新兴的一门实践性较强的技术,它是对系统多个传感器的数据进行处理的过程。众所周知,由于任何传感器都有自身的不足之处,所以单一传感器具有误报风险大,可靠性和容错能力低等缺点。为了对测量环境或对象的特征有个全面、正确的认识,克服单一传感器的上述缺点,多传感器数据融合技术应运而生。简单的说,多传感器融合技术就是融合多个传感器的信息,以产生比单个传感器更可靠、更准确的信息。常用方法有贝叶斯估计法和DS证据理论法及经典推算法等,神经网络、小波分析等智能方法近年来也是研究数据融合的重要方法和手段。卡尔曼滤波器自上世纪六十年代被提出以后,作为一种新型的滤波手段在控制、跟踪、测量领域得到广泛应用。由于卡尔曼滤波器对数据的估计是无偏最优估计,滤波器结构简单等特点,使得卡尔曼滤波器在多传感器数据融合中应用极为广泛。过去使用的集中式卡尔曼滤波器要集中处理所有传感器的数据,计算量大,实时性差,并且不具备容错性。本文基于Car lson提出的联合卡尔曼滤波算法,介绍了利用信息分配原则实现多传感器信息最优融合的滤波器的设计,不仅使系统具备了一定的容错能力,实时性也有较大幅度的提高。最后给出了联合卡尔曼滤波器在舰船组合导航中的应用实例。 2 联合卡尔曼滤波器的设计步骤 联合卡尔曼滤波器的设计主要围绕两个方面,第一是对数据进行分散处理,第二是分散处理过的数据X收稿日期:2004-11-26
卡尔曼滤波实验报告 捷联惯导静基座初始对准实验 一、实验目的 ①掌握捷联惯导的构成和基本工作原理； ②掌握捷联惯导静基座对准的基本工作原理； ③了解捷联惯导静基座对准时的每个系统状态的可观测性； ④了解双位置对准时系统状态的可观测性的变化。 二、实验原理 选取状态变量为：[]T E N E N U x y x y z X V V δδεεε=ψψψ??，其
中导航坐标系选为东北天坐标系，E V δ为东向速度误差，N V δ为北向速度误差，E ψ为东向姿态误差角，N ψ为北向姿态误差角，U ψ为天向姿态误差角，x ?为东向加速度偏置，y ?为北向加速度偏置，x ε为东向陀螺漂移，y ε为北向陀螺漂移，z ε为天向陀螺漂移。则系统的状态模型为： X AX W =+ (1) 其中， 1112212211 12 1321222331323302sin 000002sin 000000000sin cos 0000sin 000000cos 0000000000000000000000000000000000000000000000000000 0L g C C L g C C L L C C C L C C C L C C C A Ω-? ? ??-Ω????Ω-Ω? ?-Ω????Ω=? ?????? ?????????? ? [00000]E N E N U T V V W W W W W W δδψψψ=，E D V W W δψ 为零均值高斯 白噪声，分别为加速度计误差和陀螺漂移的噪声成分，Ω为地球自转角速度，ij C 为姿态矩 阵n b C 中的元素，L 为当地纬度。 量测量选取两个水平速度误差：[ ]T E N Z V V δδ=，则量测方程为： 10000000000100000000E E N N V X V δηδη???? ??=+???????????? (2) 即Z HX η=+ 其中，H 为量测矩阵，[]T E N ηηη=为量测方程的随机噪声状态矢量，为零均值高 斯白噪声。 要利用基本卡尔曼滤波方程进行状态估计，需要将状态方程和量测方程进行离散化。 系统转移矩阵为： 2323/1111102!3!! n n k k k k k k n T T T I TA A A A n ∞ -----=Φ=++++=∑ (3)
/********************************************************* // 卡尔曼滤波 //********************************************************* //在程序中利用Angle+=(Gyro - Q_bias) * dt计算出陀螺仪积分出的角度，其中Q_bias是陀螺仪偏差。 //此时利用陀螺仪积分求出的Angle相当于系统的估计值，得到系统的观测方程；而加速度计检测的角度Accel相当于系统中的测量值，得到系统状态方程。 //程序中Q_angle和Q_gyro分别表示系统对加速度计及陀螺仪的信任度。根据Pdot = A*P + P*A' + Q_angle计算出先验估计协方差的微分，用于将当前估计值进行线性化处理。其中A 为雅克比矩阵。 //随后计算系统预测角度的协方差矩阵P。计算估计值Accel与预测值Angle间的误差Angle_err。 //计算卡尔曼增益K_0,K_1，K_0用于最优估计值，K_1用于计算最优估计值的偏差并更新协方差矩阵P。 //通过卡尔曼增益计算出最优估计值Angle及预测值偏差Q_bias，此时得到最优角度值Angle 及角度值。 //Kalman滤波，20MHz的处理时间约0.77ms； void Kalman_Filter(float Accel,float Gyro) { Angle+=(Gyro - Q_bias) * dt; //先验估计 Pdot[0]=Q_angle - PP[0][1] - PP[1][0]; // Pk-先验估计误差协方差的微分 Pdot[1]=- PP[1][1]; Pdot[2]=- PP[1][1]; Pdot[3]=Q_gyro; PP[0][0] += Pdot[0] * dt; // Pk-先验估计误差协方差微分的积分 PP[0][1] += Pdot[1] * dt; // =先验估计误差协方差 PP[1][0] += Pdot[2] * dt; PP[1][1] += Pdot[3] * dt; Angle_err = Accel - Angle; //zk-先验估计 PCt_0 = C_0 * PP[0][0]; PCt_1 = C_0 * PP[1][0]; E = R_angle + C_0 * PCt_0; K_0 = PCt_0 / E;
卡尔曼滤波算法实现代码 C++实现代码如下： ============================kalman.h================= =============== // kalman.h: interface for the kalman class. // ////////////////////////////////////////////////////////////////////// #if !defined(AFX_KALMAN_H__ED3D740F_01D2_4616_8B74_8BF57636F2C0__IN CLUDED_) #define AFX_KALMAN_H__ED3D740F_01D2_4616_8B74_8BF57636F2C0__INCLU DED_ #if _MSC_VER > 1000 #pragma once #endif// _MSC_VER > 1000 #include
kalman(int x=0,int xv=0,int y=0,int yv=0); //virtual ~kalman(); }; #endif// !defined(AFX_KALMAN_H__ED3D740F_01D2_4616_8B74_8BF57636F2C 0__INCLUDED_) ============================kalman.cpp=============== ================= #include "kalman.h" #include
多传感器融合中的卡尔曼滤波探讨 1 引言 目前靠单一的信息源很难保证获取环境信息的快速性和准确性的要求，会给系统对周围环境的理解及系统的决策带来影响，另外，单一传感器获得的仅仅是环境特征的局部、片面的信息，它的信息量是十分有限的。而且每个传感器采集到的信息还受到自身品质、性能噪声的影响，采集到的信息往往是不完整的，带有较大的不确定性，偶尔甚至是错误的。而且在传统方式中，各传感器采集的信息单独、孤立的进行加工处理，不仅会导致处理工作量增加，而且割断了各传感器信息的联系丢失了信息的有机组合蕴涵的信息特征，也造成信息资 源的浪费[3-7]。在运动控制系统中，传统上就往往将速度传感器测量到的速度和加速度计测量到的加速度进行单独处理，没有将两者的信息进行数据融合。由物理定律可知，加速度与速度成导数关系，所以两者的数据是存在内在联系的，完全可以根据信息融合理论对两者数据进行综合处理，从而得到更加准确的结果。卡尔曼滤波器是常用的一种数据融合技术，它利用迭代递推计算的方式，对存贮空间要求很小，适合于存贮空间和计算速度受限的场合 [1,2]。本文分析了数度传感器和加速度计各自的优缺点，给出了一种应用卡尔曼滤波器原理对两者进行数据融合的方法。 2 传感器简介 2.1 光电编码器 光电编码器通常用于角度、位移、或转速测量，通过对光脉冲的个数进行计数再经过计算而得到测量值。假设在周长为L 的圆盘上有M 个过光孔，离散系统中，在周期时间T 内对脉冲进行计数值为N ，则第k 次测量的线速度v 可表达为 MT k Le k v MT k Le MT k LN k e k N MT L k v )()()()())()(()(+=+=+=∧ (2-1) e 是随机误差，为光脉冲取整后的剩余值，取值范围为（－1， 1），可看作均匀分布。∧v 为实际的观测值，与真值v 之间相差MT Le 。可见，在固定长度的L 上，加大M 或T 的值，都可以减小误差。但是加大M 需要付出昂贵的成本，使传感器价格大幅提高，如光栅式光电传感器；而加大T 又会降低系统的动态响应性能，所以在实际应用中，这两者均难如愿。 在需要同时测量加速度的场合，理论上可以由对速度求差分方程得出，即 2 2))1(()(())1()(()1()()(MT k e k Le MT k N k N L T k v k v k a --+--=--= (2-2) 容易看出，相对误差显著提高，数据几乎不可用，所以需要专门的加速度计对加速度进行测量。 2.2 加速度计 加速度计用于测量物体的线性加速度，根据不同的测量原理，有很多种类，本文中使用的MMA7260是一款低成本、低功耗、小体积、功能完善的单芯片加速度计，主要用于运动检测、惯性导航、震动检测、交通安全等。MMA7260响应快、带宽可调整、可响应高频率输入，但是其测量数据噪声与带宽的平方根成正比，会随着带宽增加而增加。 5.1350)(?=BW g rms Noise μ (2-3) 式中BW 为传感器带宽(HZ)。因此在设计时，首先要确定被测加速度的频率范围，然后再
LMS 自适应滤波实验报告 姓名： 学号： 日期：2015.12.2 实验内容： 利用自适应滤波法研究从宽带信号中提取单频信号的方法。 设()()()()t f B t f A t s t x 212cos 2cos π?π+++=，()t s 是宽带信号，A ，B ，1f ，2f ， ?任选 （1）要求提取两个单频信号； （2）设f f f ?+=12，要求提取单频信号()t f 22cos π，研究f ?的大小对提取单频信号的影响。 1. 自适应滤波器原理 自适应滤波器理论是现代信号处理技术的重要组成部分，它对复杂信号的处理具有独特的功能。自适应滤波器在信号处理中属于随机信号处理的范畴。在一些信号和噪声特性无法预知或他们是随时间变化的情况下，自适应滤波器通过自适应滤波算法调整滤波器系数，使得滤波器的特性随信号和噪声的变化，以达到最优滤波的效果，解决了固定全系数的维纳滤器和卡尔曼滤波器的不足。 （1） 自适应横向滤波器 所谓自适应滤波，就是利用前一时刻已获得的滤波器参数等结果，自动调节现时刻的滤波器参数，以适应信号和噪声未知或随时间变化的统计特性，从而实现最优滤波。自适应滤波器由两个部分组成：滤波器结构和调节滤波器系数的自适应算法。自适应滤波器的特点是自动调节自身的冲激响应，达到最优滤波，此算法适用于平稳和非平稳随机信号，并且不要求知道信号和噪声的统计特性。 一个单输入的横向自适应滤波器的原理框图如图所示：
实际上这种单输入系统就是一个FIR 网络结构，其输出()n y 用滤波器单位脉冲响应表示成下式： ()()()∑-=-=1 N m m n x m w n y 这里()n w 称为滤波器单位脉冲响应，令：()()n i n x x i w w m i i i ,1,1,1+-=-=+=用j 表示，上式可以写成 ∑==N i ij i j x w y 1 这里i w 也称为滤波器加权系数。用上面公式表示其输出，适用于自适应线性组合器，也适用于FIR 滤波器。将上式表示成矩阵形式： X W W X j T T j j y == 式中 [][ ] T Nj j j j T N x x x w w w X W ,...,,, ,...,,2121== 误差信号表示为 X W j T j j j j d y d e -=-= （2） 最小均方（LMS ）算法 Widrow 等人提出的最小均方算法，是用梯度的估计值代替梯度的精确值，这种算法简单易行，因此获得了广泛的应用。 LMS 算法的梯度估计值用一条样本曲线进行计算，公式如下：
卡尔曼，美国数学家和电气工程师。1930年5月 19日生于匈牙利首都布达佩斯。1953年在美国麻省理工学院毕业获理学士学位,1954年获理学硕士学位,1957年在哥伦比亚大学获科学博士学位。1957～1958年在国际商业机器公司(IBM)研究大系统计算机控制的数学问题。1958～1964年在巴尔的摩高级研究院研究控制和数学问题。1964～1971年到斯坦福大学任教授。1971年任佛罗里达大学数学系统理论研究中心主任，并兼任苏黎世的瑞士联邦高等工业学校教授。1960年卡尔曼因提出著名的卡尔曼滤波器而闻名于世。卡尔曼滤波器在随机序列估计、空间技术、工程系统辨识和经济系统建模等方面有许多重要应用。1960年卡尔曼还提出能控性的概念。能控性是控制系统的研究和实现的基本概念，在最优控制理论、稳定性理论和网络理论中起着重要作用。卡尔曼还利用对偶原理导出能观测性概念，并在数学上证明了卡尔曼滤波理论与最优控制理论对偶。为此获电气与电子工程师学会(IEEE)的最高奖──荣誉奖章。卡尔曼著有《数学系统概论》(1968)等书。 什么是卡尔曼滤波 最佳线性滤波理论起源于40年代美国科学家Wiener和前苏联科学家Ｋолмогоров等人的研究工作，后人统称为维纳滤波理论。从理论上说，维纳滤波的最大缺点是必须用到无限过去的数据，不适用于实时处理。为了克服这一缺点，60年代Kalman把状态空间模型引入滤波理论，并导出了一套递推估计算法，后人称之为卡尔曼
滤波理论。卡尔曼滤波是以最小均方误差为估计的最佳准则，来寻求一套递推估计的算法，其基本思想是：采用信号与噪声的状态空间模型，利用前一时刻地估计值和现时刻的观测值来更新对状态变量的估计，求出现时刻的估计值。它适合于实时处理和计算机运算。 卡尔曼滤波的实质是由量测值重构系统的状态向量。它以“预测—实测—修正”的顺序递推，根据系统的量测值来消除随机干扰，再现系统的状态，或根据系统的量测值从被污染的系统中恢复系统的本来面目。 释文：卡尔曼滤波器是一种由卡尔曼（Kalman）提出的用于时变线性系统的递归滤波器。这个系统可用包含正交状态变量的微分方程模型来描述，这种滤波器是将过去的测量估计误差合并到新的测量误差中来估计将来的误差。 卡尔曼滤波的应用 斯坦利.施密特(Stanley Schmidt)首次实现了卡尔曼滤波器.卡尔曼在NASA埃姆斯研究中心访问时,发现他的方法对于解决阿波罗计划的轨道预测很有用,后来阿波罗飞船的导航电脑使用了这种滤波器. 关于这种滤波器的论文由Swerling (1958), Kalman (1960)与 Kalman and Bucy (1961)发表.
LMS 自适应滤波实验报告 ： 学号： 日期：2015.12.2 实验容： 利用自适应滤波法研究从宽带信号中提取单频信号的方法。 设()()()()t f B t f A t s t x 212cos 2cos π?π+++=，()t s 是宽带信号，A ，B ，1f ，2f ， ?任选 （1）要求提取两个单频信号； （2）设f f f ?+=12，要求提取单频信号()t f 22cos π，研究f ?的大小对提取单频信号的影响。 1. 自适应滤波器原理 自适应滤波器理论是现代信号处理技术的重要组成部分，它对复杂信号的处理具有独特的功能。自适应滤波器在信号处理中属于随机信号处理的畴。在一些信号和噪声特性无法预知或他们是随时间变化的情况下，自适应滤波器通过自适应滤波算法调整滤波器系数，使得滤波器的特性随信号和噪声的变化，以达到最优滤波的效果，解决了固定全系数的维纳滤器和卡尔曼滤波器的不足。 （1） 自适应横向滤波器 所谓自适应滤波，就是利用前一时刻已获得的滤波器参数等结果，自动调节现时刻的滤波器参数，以适应信号和噪声未知或随时间变化的统计特性，从而实现最优滤波。自适应滤波器由两个部分组成：滤波器结构和调节滤波器系数的自适应算法。自适应滤波器的特点是自动调节自身的冲激响应，达到最优滤波，此算法适用于平稳和非平稳随机信号，并且不要求知道信号和噪声的统计特性。
一个单输入的横向自适应滤波器的原理框图如图所示： 实际上这种单输入系统就是一个FIR 网络结构，其输出()n y 用滤波器单位脉冲响应表示成下式： ()()()∑-=-=1 N m m n x m w n y 这里()n w 称为滤波器单位脉冲响应，令： ()()n i n x x i w w m i i i ,1,1,1+-=-=+=用j 表示，上式可以写成 ∑==N i ij i j x w y 1 这里i w 也称为滤波器加权系数。用上面公式表示其输出，适用于自适应线性组合器，也适用于FIR 滤波器。将上式表示成矩阵形式： X W W X j T T j j y == 式中 [][ ] T Nj j j j T N x x x w w w X W ,...,,, ,...,,2121== 误差信号表示为 X W j T j j j j d y d e -=-= （2） 最小均方（LMS ）算法 Widrow 等人提出的最小均方算法，是用梯度的估计值代替梯度的精确值，这种算法简单易行，因此获得了广泛的应用。
Kalman 滤波算法 姓名：刘金强 专业：控制理论与控制工程 学号：2007255 ◆实验目的： （1）、掌握klman 滤波实现的原理和方法 （2）、掌握状态向量预测公式的实现过程 （3）、了解Riccati 差分方程实现的过程和新息的基本性质和过程的计算 ◆实验要求： 问题： F=[a1,a2,a3],其中a1=[1.0 0 0]的转置，a2=[0.3 1.0 0]的转置，a3=[0.1 0.2 0.4]的转置，x(0)=[3,-1,2]的转置；C=[b1,b2,b3],其中b1=[0.3 0.5]的转置，b2=[1,0.4]的转置，b3=[0.8 -0.7]的转置；V1(n)=[0 0 n1(n)sin(0.1n)]的转置，V2(n)=[n2(n) n3(n)];n1(n)为均值为零，方差为1的均匀分布白噪声;n2(n),n3(n)为均值为0，方差为0.1的均匀分布白噪声,n1(n),n2(n),n3(n)相互独立，试用卡尔曼滤波器算法估计x^(n). ◆实验原理： 初始条件： 1?(1)x =E{x(1)} K(1,0)=E{[x(1)- (1)x ][x(1)- (1)H x ]},其中(1)x =E{x(1)} 输入观测向量过程： 观测向量序列={y(1),…………y(n)} 已知参数： 状态转移矩阵F(n+1,n) 观测矩阵C(n) 过程噪声向量的相关矩阵1()Q n 观测噪声向量的相关矩阵2()Q n 计算：n=1,2,3,………………. G(n)=F(n+1,n)K(n,n+1) ()H C n 12[()(,1)()()]H C n K n n C n Q n --+ Kalman 滤波器是一种线性的离散时间有限维系统。Kalman 滤波器的估计性能是：它使滤波后的状态估计误差的相关矩阵P(n)的迹最小化。这意味着，kalman 滤波器是状态向量x(n)的线性最小方差估计。 ◆实验结果： ◆程序代码： （1）主程序
卡尔曼滤波实验报告
2014 年 4 月 GPS 静/动态滤波实验 一、实验要求 1、分别建立GPS 静态及动态卡尔曼滤波模型，编写程序对静态和动态GPS 数据进行Kalman 滤波。 2、对比滤波前后导航轨迹图。 3、画出滤波过程中估计均方差(P 阵对角线元素开根号)的变化趋势。 4、思考：① 简述动态模型与静态模型的区别与联系；② R 阵、Q 阵，P0阵的选取对滤波精度及收敛速度有何影响，取值时应注意什么；③ 本滤波问题是否可以用最小二乘方法解决，如果可以，请阐述最小二乘方法与Kalman 滤波方法的优劣对比。 二、实验原理 2.1 GPS 静态滤波 选取系统的状态变量为[ ]T L h λ=X ,其中L 为纬度(deg),λ为经度(deg)，h 为高度 (m)。设()w t 为零均值高斯白噪声，则系统的状态方程为： 310()w t ?=+X (1) 所以离散化的状态模型为： ,111k k k k k W ---=+X X Φ (2) 式中，,1k k -Φ为33?单位阵，k W 为系统噪声序列。 测量数据包括：纬度静态量测值、经度静态量测值和高度构成31?矩阵Z ，量测方程
可以表示为： k k k Z HX V =+ (3) 式中，H 为33?单位阵，k V 为量测噪声序列。 系统的状态模型是十分准确的，所以系统模型噪声方差阵可以取得十分小，取Q 阵零矩阵。 系统测量噪声方差阵R 由测量确定，由于位置量测精度为5m ，采用克拉索夫斯基地球椭球模型，长半径e R 为6378245m ，短半径p R 为6356863m 。所以R 阵为： 2 2 25180()0 05180 ( )0cos()00 5p e R R L ππ ??? ?? ? ??= ??? ? ? ?? ? R (4) 2.2 GPS 动态滤波 动态滤波基于当前统计模型，在地球坐标系下解算。选取系统的状态变量为 T x x x y y y z z z X x v a y v a z v a εεε??=??,其中,,,x x x x v a ε依次为地球坐标系下x 轴上的位置、速度、加速度和位置误差分量，,y z 轴同理。系统的状态模型可以表示为： ()()()()t t t t =++X AX U W (5) 式中，位置误差视为有色噪声，为一阶马尔科夫过程，可表示为： x x x x y y y y z z z z w w w εετεετεετ?=-+????=-+????=-+?? 1 11 (6) 其中，i τ(,,i x y z =)为对应马尔科夫过程的相关时间常数，(,,)i w i x y z =为零均值高斯白噪声。
自适应卡尔曼滤波 卡尔曼滤波发散的原因 如果卡尔曼滤波是稳定的，随着滤波的推进，卡尔曼滤波估计的精度应该越来越高，滤波误差方差阵也应趋于稳定值或有界值。但在实际应用中，随着量测值数目的增加，由于估计误差的均值和估计误差协方差可能越来越大，使滤波逐渐失去准确估计的作用，这种现象称为卡尔曼滤波发散。 引起滤波器发散的主要原因有两点： （1）描述系统动力学特性的数学模型和噪声估计模型不准确，不能直接真实地反映物理过程，使得模型与获得的量测值不匹配而导致滤波发散。这种由于模型建立过于粗糙或失真所引起的发散称为滤波发散。 （2）由于卡尔曼滤波是递推过程，随着滤波步数的增加，舍入误差将逐渐积累。如果计算机字长不够长，这种积累误差很有可能使估计误差方差阵失去非负定性甚至失去对称性，使滤波增益矩阵逐渐失去合适的加权作用而导致发散。这种由于计算舍入误差所引起的发散称为计算发散。 针对上述卡尔曼滤波发散的原因，目前已经出现了几种有效抑制滤波发散的方法，常用的有衰减记忆滤波、限定记忆滤波、扩充状态滤波、有限下界滤波、平方根滤波、和自适应滤波等。这些方法本质上都是以牺牲滤波器的最优性为代价来抑制滤波发散，也就是说，多数都是次优滤波方法。 自适应滤波 在很多实际系统中，系统过程噪声方差矩阵Q和量测误差方差阵R事先是不知道的，有时甚至连状态转移矩阵 或量测矩阵H也不能确切建立。如果所建立的模型与实际模型不符可能回引起滤波发散。自适应滤波就是这样一种具有抑制滤波发散作用的滤波方法。在滤波过程中，自适应滤波一方面利用量测值修正预测值，同时也对未知的或不确切的系统模型参数和噪声统计参数进行估计修正。自适应滤波的方法很多，包括贝叶斯法、极大似然法、相关法与协方差匹配法，其中最基本也是最重要的是相关法，而相关法可分为输出相关法和新息相关法。 在这里只讨论系统模型参数已知，而噪声统计参数Q和R未知情况下的自适应滤波。由于Q和R等参数最终是通过增益矩阵K影响滤波值的，因此进行自适应滤波时，也可以不去估计Q和R等参数而直接根据量测数据调整K就可以了。
卡尔曼滤波实验报告 2014 年 4 月 GPS静/动态滤波实验 一、实验要求 1、分别建立GPS静态及动态卡尔曼滤波模型，编写程序对静态和动态GPS数据进行Kalman滤波。 2、对比滤波前后导航轨迹图。
3、画出滤波过程中估计均方差(P 阵对角线元素开根号)的变化趋势。 4、思考：① 简述动态模型与静态模型的区别与联系；② R 阵、Q 阵，P0阵的选取对滤波精度及收敛速度有何影响，取值时应注意什么；③ 本滤波问题是否可以用最小二乘方法解决，如果可以，请阐述最小二乘方法与Kalman 滤波方法的优劣对比。 二、实验原理 2.1 GPS 静态滤波 (deg) 度(m) (1) 所以离散化的状态模型为： (2) 可以表示为： (3) 矩阵。 5m ，采用克拉索夫斯基地球 6378245m 6356863m (4) 2.2 GPS 动态滤波 动态滤波基于当前 统计模型，在地球坐标系下解算。选取系统的状态变量为 (5)
式中，位置误差视为有色噪声，为一阶马尔科夫过程，可表示为： ε τεετεετ-=- =-1 1 (6) 白噪声。 (7) (8) 系统噪声为： (9) 量测量为纬度动态量测值、经度动态量测值、高度和三向速度量测值。由于滤波在地球 坐标系下进行，为了简便首先将纬度、经度和高度转化为三轴位置坐标值，转化方式如下： (10) 量测方程为： (11)
综上，离散化的Kalman滤波方程为： (12) 离散化的系统噪声协方差阵为： 2 [ π ?] ? k x = +<0 “当前”加速度 (13) 离散化量测噪声协方差阵为：diag = R 三、实验结果 3.1 GPS静态滤波
卡尔曼滤波与组合导航》课程实验报告 实验 捷联惯导 /GPS 组合导航系统静态导航实验 实验序号 3 姓名 陈星宇 系院专业 17 班级 ZY11172 学号 ZY1117212 日期 2012-5-15 指导教师 宫晓琳 成绩 、实验目的 ① 掌握捷联惯导 /GPS 组合导航系统的构成和基本工作原理； ②掌握采用卡尔曼滤波方法进行捷联惯导 /GPS 组合的基本原理； ③掌握捷联惯导 /GPS 组合导航系统静态性能； ④了解捷联惯导 /GPS 组合导航静态时的系统状态可观测性； 、实验原理 （ 1）系统方程 X FX GW 系统噪声矢量由陀螺仪和加速度计的随机误差组成，表达式为： 2）量测方程 和 H 分别为捷联解算与 GPS 的东向速度、北向速度、天向速度、纬度、经度和高度之 差；量测矩阵 H H V H P T ，H P 03 6 diag R M H, （R N H ）cos L, 036 ， H V 033 diag 1, 1, 1 039 ，v v V E v V N v V U v L v v H 为量测噪声。 量测噪声 v E v N T v U L h x y z x y z 其中， E 、 N 、 U 为数学平台失准角； v E 、 v N 、 v U 分别为载体的东向、北向和天向速度误差； L 、 、 h 分别为纬度误差、经度误差和高度误差； x 、 y 、 z 、 x 、 y 、 z 分别为陀螺随 机常值漂移和加速度计随机常值零偏。（下 标 系统的噪声转移矩阵 G 为： E 、N 、 U 分别代表东、北、天） C b n 3 3 0 9 3 3 3 C n C b 9 3 15 6 系统的状态转移矩阵 w w w w F 组成内容为： w z F 06N 9 F S F M ，其中 F N 中非零元素为可由惯导误差模型获得。 F S C b n 3 3 0 3 3 3 3 C b n 3 3 96 量测变量 z V E V N V U L H ， ， V E 、 V N 、 V U 、 L 、 X U
两传感器分布式kalman滤波融合算法及其仿真分析 摘要：讨论了基于两传感器kalman滤波的数据融合算法，对FAFSS算法机理进行了描述并融合算法进行了仿真，分析了融合结果。 关键字：kalman滤波；分布式传感器信息融合；分布式滤波数据融合算法
卡尔曼滤波实验报告 2014 年4 月 GPS静/动态滤波实验 一、实验要求 1、分别建立GPS静态及动态卡尔曼滤波模型，编写程序对静态和动态GPS数据进行Kalman滤波。 2、对比滤波前后导航轨迹图。
3、画出滤波过程中估计均方差(P 阵对角线元素开根号)的变化趋势。 4、思考：① 简述动态模型与静态模型的区别与联系；② R 阵、Q 阵，P0阵的选取对滤波精度及收敛速度有何影响，取值时应注意什么；③ 本滤波问题是否可以用最小二乘方法解决，如果可以，请阐述最小二乘方法与Kalman 滤波方法的优劣对比。 二、实验原理 GPS 静态滤波 选取系统的状态变量为[ ]T L h λ=X ,其中L 为纬度(deg),λ为经度(deg)，h 为高度(m)。设()w t 为零均值高斯白噪声，则系统的状态方程为： 31 0()w t ?=+&X (1) 所以离散化的状态模型为： ,111k k k k k W ---=+X X Φ (2) 式中，,1k k -Φ为33?单位阵，k W 为系统噪声序列。 测量数据包括：纬度静态量测值、经度静态量测值和高度构成31?矩阵Z ，量测方程 可以表示为： k k k Z HX V =+ (3) 式中，H 为33?单位阵，k V 为量测噪声序列。 系统的状态模型是十分准确的，所以系统模型噪声方差阵可以取得十分小，取Q 阵零矩阵。 系统测量噪声方差阵R 由测量确定，由于位置量测精度为5m ，采用克拉索夫斯基地球椭球模型，长半径e R 为6378245m ，短半径p R 为6356863m 。所以R 阵为： 22 25180()0 05180 ( )0cos()00 5p e R R L ππ ??? ?? ? ??= ??? ? ? ?? ? R (4) GPS 动态滤波 动态滤波基于当前统计模型，在地球坐标系下解算。选取系统的状态变量为 T x x x y y y z z z X x v a y v a z v a εεε??=??,其中,,,x x x x v a ε依次为地球坐标系下x 轴上的位置、速度、加速度和位置误差分量，,y z 轴同理。系统的状态模型可以表示为： ()()()()t t t t =++X AX U W & (5) 式中，位置误差视为有色噪声，为一阶马尔科夫过程，可表示为：