附录：kalman滤波（起源、发展、原理、应用）
1、Kalman滤波起源及发展
1960年，匈牙利数学家卡尔曼发表了一篇关于离散数据线性滤波递推算法的论文，这意味着卡尔曼滤波的诞生。斯坦利.施密特(Stanley Schmidt)首次实现了卡尔曼滤波器，卡尔曼在NASA埃姆斯研究中心访问时，发现他的方法对于解决阿波罗计划的轨道预测很有用，后来阿波罗飞船的导航电脑使用了这种滤波器。关于这种滤波器的论文由Swerling (1958)、Kalman (1960)与Kalman and Bucy (1961)发表.
卡尔曼滤波是一种有着相当广泛应用的滤波方法，但它既需要假定系统是线性的，又需要认为系统中的各个噪声与状态变量均呈高斯分布，而这两条并不总是确切的假设限制了卡尔曼滤波器在现实生活中的应用。扩展卡尔曼滤波器（EKF）极大地拓宽了卡尔曼滤波的适用范围。EKF的基本思路是，假定卡尔曼滤滤对当前系统状态估计值非常接近于其真实值，于是将非线性函数在当前状态估计值处进行台劳展开并实现线性化。另一种非线性卡尔曼滤波叫线性化卡尔曼滤波。它与EKF的主要区别是前者将非线函数在滤波器对当前系统状态的最优估计值处线性化，而后者因为预先知道非线性系统的实际运行状态大致按照所要求、希望的轨迹变化，所以这些非线性化函数在实际状态处的值可以表达为在希望的轨迹处的台劳展开式，从而完成线性化。
不敏卡尔曼滤波器（UKF）是针对非线性系统的一种改进型卡尔曼滤波器。UKF处理非线性系统的基本思路在于不敏变换，而不敏变换从根本上讲是一种描述高斯随机变量在非线性化变换后的概率分布情况的方法。不敏卡尔曼滤波认为，与其将一个非线性化变换线性化、近似化，还不如将高斯随机变量经非线性变换后的概率分布情况用高斯分布来近似那样简单，因而不敏卡尔曼滤波算法没有非线性化这一步骤。在每一定位历元，不敏卡尔曼滤波器按照一套公式产生一系列样点，每一样点均配有一个相应的权重，而这些带权的样点被用来完整地描述系统状态向量估计值的分布情况，它们替代了原先卡尔曼滤波器中的状态向量估计值及协方差。不敏卡尔曼滤器让这些样点一一经历非线性状态方程与测量方程，然后再将这些经非线性变换后的样点按照它们的权重而综合出对当前时刻的系统状态向量估计值。
多态自适应（MMA）卡尔曼滤波器是一种受到广泛关注的滤波器，它由好多个并联、同时运行的卡尔曼滤波器组成。在这组卡尔曼滤波器中，每一个滤波器对未知的滤波参数分别做出相互不同的假设，然后各自按照自己的模型假设进行滤波计算，而多态自适应滤波器最后将它们对系统状态的各个估计值进行加权，并以此作为最优估计值输出。
2、Kalman 原理
卡尔曼滤波是基于状态空间方法的一套递推滤波算法，在状态空间方法中，引入了状态变量的概念。实际应用中，可以通过选取合适的状态变量来体现系统的特征、特点和状况的变化。卡尔曼滤波的模型包括状态空间模型和观测模型。状态模型是反映状态变化规律的模型，通过状态方程来描写相邻时刻的状态转移变化规律；观测模型反映了实际观测量与状态变量之间的关系。Kalman 滤波问题就是联合观测信息及状态转移规律来得到系统状态的最优估计。
假设动态系统的状态空间模型为
)()()1(t W t X t X Γ+Φ=+
（2-1） )()()(t V t HX t Y +=
（2-2） 其中，
X(t) 系统在时刻t 的状态
Y(t) 对状态的观测值
W(t) 系统噪声，方差阵为Q
V(t) 观测噪声，方差阵为R
Φ 状态转移矩阵
H 观测矩阵
Γ 系统噪声驱动矩阵
卡尔曼滤波的计算流程为：
计算状态估计值：
)1()1()|1(ˆ)1|1(ˆ++++=++t t K t t X t t X ε
（2-3） 计算状态一步预测：
)|(ˆ)|1(ˆt t X t t X Φ=+
（2-4） 计算新息：
)|1(ˆ)1()1(t t X H t Y t +-+=+ε
（2-5） 计算卡尔曼滤波增益：
1])|1([)|1()1(-+++=+R H t t HP H t t P t K T T
（
2-6） 计算一步预测均方误差：
T T Q t t P t t P ΓΓ+ΦΦ=+)|()|1( （2-7）
计算一步预测估计均方误差：
)|1(])1([)1|1(t t P H t K I t t P n ++-=++ （2-8）
为了更形象地说明卡尔曼滤波的原理，下面给出卡尔曼滤波的系统模型框图：
图1 卡尔曼滤波的系统模型
3、Kalman 滤波应用之一：基于VDLL 的GPS 信号跟踪算法[2]
图2 VDLL 的基本结构
系统的状态量为用户状态( 位置、速度、用户钟差、用户钟漂移) ，观测量为所有卫星通道的码环鉴相器输出组成的矢量D(X)
X =(x y z b x • y • z • b •
)T
D(X)=(D(Δτ1) D(Δτ2) … D(ΔτN ))T （3-1）
其中x ，y ，z ，b ，x •，y •，z •，b •分别为在ECEF( Earth Centered Earth Fixed, 地心坐标系) 下的用户接收机三维位置和速度， b , b •分别为接收机钟差和钟差变化率。预检测积分时间T 选择为1ms ，即滤波器迭代的k 时刻与k+1时刻间隔1ms 。假设用户为匀加速运动，可写出状态方程如下 1k 1k 4444k w X I 0I I X --+⎥⎦
⎤⎢⎣⎡=T （3-2） 其中w k-1 为状态扰动噪声。 导航滤波器把k- 1 时刻的状态和对应的卫星位置信息反馈给本地信号生成模块，产生超前、即时与滞后码序列，分别与k 时刻的1ms 卫星信号的同相支路与正交支路做相关运算，即图2 中的“积分存储”模块。为了去除对信号幅度变化的敏感性和对载波环路锁定的依赖性，码环鉴相器选择归一化的超前码与滞后码相关能量差( 第i 颗卫星)
2
2222
222,)(QL IL QE IE QL IL QE IE D k i +++--+=∆τ （3-3） 该码环鉴相器在没有噪声的理想情况下，当超前码与滞后码间距为1 chip( 码片) 时，在±0.5 chip 输入误差范围内，鉴相器输出等于真实的跟踪误差。在实际应用中，超前滞后码的移位通常以采样点为单位，1 个chip 的长度不一定为整数个采样点，但鉴相器仍然与输入误差成线性关系，所以式可写为
k i k i i k i v g D ,,,)(+∆=∆ττ （3-4）
其中g i 为待求的鉴相器特性在线性区的斜率， v i, k 为观测噪声。鉴相器输入误差k i ,τ∆是k 时刻( 接收机时刻) 第i 颗卫星1ms 长度的采样信号与本地即时码之间的相位偏移，而本地即时码的初始相位值等于k- 1时刻对卫星信号的相位估计。所以Δτi, k 实质上是k 时刻卫星信号的相位与k - 1时刻的相位估计之差，等于k 时刻信号传输延迟的真实值与k- 1 时刻信号传输延迟的估计值。 k 时刻接收到第i 颗卫星信号的理论测量延迟为
