Download presentation
Presentation is loading. Please wait.
1
第四章 维纳滤波和卡尔曼滤波 引言 维纳滤波器的离散形式——时域解 离散维纳滤波器的z域解 维纳预测 卡尔曼(Kalman)滤波
2
4.1 引 言 问题提出： x(n) 信道噪声（测量噪声） s(n) 信道 s(n)：原始输入（发射）信号，随机平稳
4.1 引 言 问题提出： 信道噪声（测量噪声） s(n) 信道 x(n) s(n)：原始输入（发射）信号，随机平稳 x(n)：接收（测量）信号， 随机平稳
3
s(n) 已知：s(n), (n) 的统计特性，要求： 设计线性移不变滤波器 h(n), 从x(n)中恢复s(n) FIR,IIR
逼近（准则） 准则：最大后验准则。均方准则，最大似然准则， 线性均方准则（最小二乘滤波）
4
2.横向滤波器结构 M个权系数（抽头）的横向滤波器 定义： ：输入信号 ：输入向量
5
:滤波器的权系数 :滤波器权向量 :期望响应 :对期望响应的估计 :估计误差
6
假设 由信号 与噪声 组成 ⑴ 如果 ，上图的系统称为滤波（filtering）； ⑵ 如果 ，上图的系统称为预测（prediction）； ⑶ 如果 ，上图的系统称为平滑 （smoothing）。
7
4.2 维纳滤波器的离散形式—时域解 （加性干扰） d(n)=s(n) 最小均方误差准则：
8
一、维纳—霍夫方程 为此令 正交性原理：使代价函数最小化的充要条件是 n 时刻的最优估计误差正交于n 时刻滤波器的每个输入值，或者说正交于n时刻的输入信号空间。
9
一、维纳—霍夫方程 推论： n时刻的最优估计误差正交于n时刻滤波器的最优输出值
10
一、维纳—霍夫方程 由正交方程可得：
11
一、维纳—霍夫方程 定义 可得 维纳—霍夫（Wiener-Holf)方程或标准方程 求和范围（i）随滤波器的不同取不同区间
12
FIR 维纳滤波器 对FIR结构，假设其长度为N,期望信号为s(n) 令
14
令 Toeplitz矩阵,NXN对称半正定
15
令 代入可得：
16
则 令
17
可以看出，均方误差与滤波器的单位脉冲响应是一个二次函数关系。由于单位脉冲响应h(n) 为M维向量，因此均方误差是一个超椭圆抛物形曲面，该曲面有极小点存在。当滤波器工作于最佳状态时，均方误差取得最小值。
18
即 上式表明已知期望信号与观测数据的互相关函数及观测数据的自相关函数时，可以通过矩阵求逆运算，得到维纳滤波器的最佳解。同时可以看到，直接从时域求解因果的维纳滤波器，当选择的滤波器的长度M较大时，计算工作量很大，并且需要计算Rxx 的逆矩阵，从而要求的存贮量也很大。此外，在具体实现时，滤波器的长度是由实验来确定的，如果想通过增加长度提高逼近的精度，就需要在新N基础上重新进行计算。因此，从时域求解维纳滤波器，并不是一个有效的方法。
19
由正交方程可知：误差与输入信号矢量正交，可推得其与估计值也正交，用下图表示。
几何解释： 图 期望信号、 估计值与误差信号的几何关系
20
图表明在滤波器处于最佳工作状态时， 估计值加上估计偏差等于期望信号， 即
注意我们所研究的是随机信号，图中各矢量的几何表示应理解为相应量的统计平均或者是数学期望。再从能量的角度来看，假定输入信号和期望信号都是零均值, 应用正交性原理，则 , 因此在滤波器处于最佳状态时，估计值的能量总是小于等于期望信号的能量。
21
非因果性考虑 可以证明：非因果Wiener滤波器的性能(误差方差性能)要优于因果Wiener滤波器(参见郑南宁编《数字信号处理》)。所以，在实际FIR滤波器中，常用时延方法用可实現的因果系统逼近非因果系统。 令
22
因果系统n时刻的输出可以逼近非因果系统(n-M)的输出。
23
例1： 已知： 设计N=4的FIR最佳滤波器
24
解：
25
同样：
26
例 设y(n)=x(n)+v2(n)，v2(n)是一白噪声，方差σ22=0
例 设y(n)=x(n)+v2(n)，v2(n)是一白噪声，方差σ22=0.1。期望信号x1(n)的信号模型如图(a)所示，其中白噪声v1(n)的方差σ21=0.27，且b0=0.8458。x(n)的信号模型如图（b）所示，b1=0.9458。假定v1(n)与v2(n)、x1(n)与y(n)不相关，并都是实信号。设计一个维纳滤波器，得到该信号的最佳估计，要求滤波器是一长度为 2 的FIR滤波器。
27
图 输入信号与观测数据的模型 解 这个问题属于直接应用维纳-霍夫方程的典型问题，其关键在于求出观测信号的自相关函数和观测信号与期望信号的互相关函数。 图 维纳滤波器的框图
28
根据题意，画出这个维纳滤波器的框图，如图所示。 用H1(z)和H2(z)分别表示x1(n)和x(n)的信号模型，那么滤波器的输入信号x(n)可以看作是v1(n)通过H1(z)和H２(z)级联后的输出， H1(z)和H２(z)级联后的等效系统用H(z)表示，输出信号y(n)就等于x(n)和v2(n)之和。因此求出输出信号的自相关函数矩阵Ryy和输出信号与期望信号的互相关矩阵Ryd是解决问题的关键。相关函数矩阵由相关函数值组成，已知x(n)与v2(n)不相关，那么
29
(1) 求出期望信号的方差。根据图 (a)，期望信号的时间序列模型所对应的差分方程为
x1(n)=v1(n)-b0x1(n-1) 这里，b0=0.8458, 由于x1(n)的均值为零，其方差与自相关函数在零点的值相等。
30
(2) 计算输入信号和输出信号的自相关函数矩阵。根据自相关函数、功率谱密度和时间序列信号模型的等价关系，已知时间序列信号模型，就可以求出自相关函数。这里，信号的模型H(z)可以通过计算得到。
这是一个二阶系统，所对应的差分方程为 x(n)+a1x(n-1)+a2x(n-2)=v1(n) 式中，a1=-0.1，a2=-0.8。由于v1(n)、v2(n)的均值为零，因此， x(n)的均值为0。将方程两边同乘以x*(n-m)，并取数学期望，得 rxx(m)+a1rxx(m-1)+a2rxx (m-2)=0 m＞0 (1) rxx(0)+ a1rxx(1)+a2rxx(2)=σ21 m=0 (2)
31
对方程（1）取m=1, 2，得到 rxx(1)+a1rxx(0)+a2rxx(1)=0 (3) rxx(2)+a1rxx(1)+a2rxx(0)=0 (4) 方程(2)、(3)、(4)联立求解，得 至此， 输入信号的自相关矩阵Rxx可以写出：
32
v2(n)是一个零均值的白噪声，它的自相关函数矩阵呈对角形， 且 ，
因此，输出信号的自相关Ryy为
33
(3) 计算输出信号与期望信号的互相关函数矩阵。 由于两个信号都是实信号，故
ryd(m)=E［y(n)d(n-m)］=E［y(n)x1(n-m)］ =E［(x(n)+v2(n))x1(n-m)］ =E［x(n)x1(n-m)］ m=0, 1 根据系统H2(z)的输入与输出的关系，有 x1(n)-b1x(n-1)=x(n) 推出 x1(n)=x(n)+b1x(n-1) ryd(m) =E［x(n)x1(n-m)］ =E［x(n)(x(n-m)+b1x(n-1-m))］ =rxx(m)+b1rxx(m-1) 这样
34
将m=0, m=1代入上式, 得 ryd(0)=rxx(0)+b1rxx(-1)= ×0.5=0.5272 ryd(1)=rxx(1)+ b1rxx(0)= ×1= 因此，输出信号与期望信号的互相关Ryd为 求出输出信号自相关的逆矩阵, 并乘以Ryd, 就可以得到维纳滤波器的最佳解Wopt:
35
可以计算出该维纳滤波达到最佳状态最小值E［|e(n)|2］min：
36
2.3 离散维纳滤波器的 z 域解 一、非因果IIR 维纳滤波器 标准方程： 双边Z变换
37
Sxx(z)=Sss(z)+Svv(z)
假设信号和噪声不相关，即rsv(m)=0，则 Sxs(z)=Sss(z) Sxx(z)=Sss(z)+Svv(z) 则有 当噪声为0时，信号全部通过；当信号为0时， 噪声全部被抑制掉，因此维纳滤波确有滤除噪声的能力。把信号的频谱用Pss(ejω)表示，噪声的频谱用Pvv(ejω)表示，那么非因果的维纳滤波器的传输函数Hopt(ejω)的幅频特性如图所示。
38
Pss(ejω)≠0, Pvv(ejω)=0 Pss(ejω)≠0, Pvv(ejω) ≠ 0
图 非因果维纳滤波器的传输函数的幅频特性
39
因果情况处理思路： 然而实际的系统都是因果的。对于一个因果系统，不能直接转入频域求解的原因是由于输入信号与期望信号的互相关序列是一个因果序列，如果能够把因果维纳滤波器的求解问题转化为非因果问题，求解方法将大大简化。那么怎样把一个因果序列转化为一个非因果序列呢？
40
二、因果IIR 维纳滤波器 标准方程： 一般情况下直接求解比较困难。但如果滤波器输入 是白噪声，则维纳-霍夫方程容易求解；而任何平稳
v(n) x(n) S^(n) 信道 s(n) h(n) IIR 标准方程： 一般情况下直接求解比较困难。但如果滤波器输入 是白噪声，则维纳-霍夫方程容易求解；而任何平稳 随机信号可变换为等效的白噪声过程，故借助谱分解 定理可找到一种简单解决方法。
41
因果维纳滤波器的求解方法1 回顾前面讲到的时间序列信号模型，假设x(n)的信号模型B(z)已知（如图 (a)所示），求出信号模型的逆系统B-1(z)， 并将x(n)作为输入，那么逆系统B-1(z)的输出ω(n)为白噪声，白化滤波器（如图 (b)所示）。 (a) (b) 图 x(n)的时间序列信号模型及其白化滤波器
42
（1）若滤波器的输入是白噪声时 对应传递函数：
43
（2）若滤波器的输入是平稳随机信号时 则x(n)可看作是由白噪声w(n)激励一个线性移不变系统的输出 B(z) 为有理分式，
b(n) w(n) w(n) x(n) y(n) 信号模型 白化 滤波器 s(n) B(z) 为有理分式， N(z),D(z)为最小 相位多项式 问题转化为求 的问题
44
由（1） 其中
45
2、因果维纳滤波器的求解方法2 具体思路如图所示。用白噪声作为待求的维纳滤波器的输入，设定1/B(z)为信号x(n)的白化滤波器的传输函数，那么维纳滤波器的传输函数G(z)的关系为 因此，维纳滤波器的传输函数H(z)的求解转化为G(z) 的求解。 图 维纳滤波解题思路
46
假设待求维纳滤波器的单位脉冲响应为g(n)，期望信号d(n)=s(n)，系统的输出信号y(n)=s(n)，g(n)是G(z)的逆Z变换， 则输出信号可表示为
48
可以看出，均方误差的第一项和第三项都是非负数, 要使均方误差为最小，当且仅当
-∞＜k＜∞ 因此g(n)的最佳值为 -∞＜k＜∞ 对上式两边同时做Z变换，得到
49
这样，非因果维纳滤波器的最佳解为 因为s(n)=s(n)*δ(n)，且x(n)=ω(n)*b(n)，根据相关卷积定理，得到 rxs(m)=rωs(m)*b(-m) 对上式两边做Z变换，得到 Sxs (z)=Sωs(z)B(z-1) 因此
50
根据x(n)的信号模型，得到非因果的维纳滤波器的复频域最佳解的一般表达式
假定信号与噪声不相关，即当E［s(n)v(n)］=0时，有 rxs(m)=E［(s(n)+v(n))*s(n+m)］=rss(m) rxx(m)=E［(s(n)+v(n))*(s(n+m)+v(n+m))］ =rss(m)+rvv(m) 对上边两式做Z变换，得到 Sxs(z)=Sss(z) Sxx(z)=Sss(z)+Svv(z)
51
信号和噪声不相关时，非因果维纳滤波器的复频域最佳解和频率响应分别为
52
滤波器的最小均方误差E［|e(n)|2］min的计算，
根据围线积分法求逆Z变换的公式，rss(m)用下式表示: 得出
53
由复卷积定理 取y(n)=x(n)，有 因此
54
因为实信号的自相关函数是偶函数，即rss(m)=rss(-m)，
因此 Sss(z)=Sss(z-1) 假定信号与噪声不相关，E［s(n)v(n)］=0， 则
55
因果维纳滤波器的求解 g(n)=0 n＜0 则滤波器的输出信号 估计误差的均方值 E［|e(n)|2］=E［|s(n)-y(n)|2］ 得到
若维纳滤波器是一个因果滤波器， 要求 g(n)=0 n＜0 则滤波器的输出信号 估计误差的均方值 E［|e(n)|2］=E［|s(n)-y(n)|2］ 得到
56
要使均方误差取得最小值， 当且仅当 令
57
所以因果维纳滤波器的复频域最佳解为
58
维纳滤波的最小均方误差为
59
比较可以看出因果维纳滤波器的最小均方误差与非因果维纳滤波器的最小均方误差的形式相同，但公式中的Hopt(z)的表达式不同。 前面已经导出， 对于非因果情况，
对于因果情况， 比较两式，它们的第二项求和域不同，因为因果情况下,k=0~+∞， 因此可以说明非因果情况的E［|e(n)|2］min一定小于等于因果情况E［|e(n)|2］min。在具体计算时,可以选择单位圆作为积分曲线， 应用留数定理， 计算积分函数在单位圆内的极点的留数来得到。
60
因果维纳滤波器设计的一般步骤： (1) 根据观测信号x(n)的功率谱求出它所对应的信号模型的传输函数，即采用谱分解的方法得到B(z)。具体方法为Sxx(z)=σ2ωB(z)B(z-1)，把单位圆内的零极点分配给B(z)，单位圆外的零极点分配给B(z-1)，系数分配给σ2ω。 (2)求 的Z反变换，取其因果部分再做Z变换，即舍掉单位圆外的极点，得 (3) 计算Hopt(z)， E［|e(n)|2］min。
61
例 已知 信号和噪声不相关，即rsv(m)=0，噪声v(n)是零均值、单位功率的白噪声（σ2v=1，mv=0)，求Hopt(z)和E［|e(n)|2］min。 解 根据白噪声的特点得出Svv(z)=1，由噪声和信号不相关， 得到 rxx(m)=rss(m)+rvv(m)。 对上式两边做Z变换，并代入已知条件，对x(n)进行功率谱分解:
62
考虑到B(z)必须是因果稳定的系统，得到连
(1) 首先分析物理可实现情况，应用公式（2.3.38）: 令 F(z)的极点为0.8和2，考虑到因果性、稳定性，仅取单位圆内的极点zi=0.8，f(n)为F(z)的Z反变换。用Res表示留数，应用留数定理，有
63
取因果部分， f+(n)=0.6×0.8n×u(n)
64
令
65
单位圆内只有极点zi=0.5， 未经滤波器的均方误差
66
(2) 对于非物理可实现情况，有
67
令 单位圆内有两个极点0.8和0.5，应用留数定理，有 比较两种情况下的最小均方误差,可以看出非物理可实现情况的最小均方误差小于物理可实现情况的均方误差。
68
4.4 维 纳 预 测 一、 维纳预测的计算 在维纳滤波中，期望的输出信号 yd(n)=s(n)，实际的输出为y(n)=s(n)。在维纳预测中，期望的输出信号yd(n)=s(n+N)，实际的输出y(n)=s(n+N)。前面已经推导得到维纳滤波的最佳解为 ^ 其中,Sxx(z)是观测数据的功率谱;Sxyd(z)是观测数据与期望信号的互功率谱，即互相关函数rxyd(k)的傅里叶变换
69
对应于维纳预测器，其输出信号y(n)和预测误差信号e(n+N)分别为
同理，要使预测误差的均方值为最小，须满足 其中，hk表示h(k)。
70
观测数据与期望的输出的互相关函数rxyd(k)和互谱密度Sxyd(z)分别为
这样，非因果维纳预测器的最佳解为 因果维纳预测器的最佳解为
71
维纳预测的最小均方误差为 从上面分析可以看出， 维纳预测的求解和维纳滤波器的求解方法是一致的。
72
二、 纯预测 假定维纳预测器是因果的，仍设s(n)与v(n)不相关，纯预测情况下的输入信号的功率谱及维纳预测器的最佳解分别为
假设x(n)=s(n)+v(n)，式中v(n)是噪声，且v(n)=0，期望信号为s(n+N)，N＞0，此种情况称为纯预测。 假定维纳预测器是因果的，仍设s(n)与v(n)不相关，纯预测情况下的输入信号的功率谱及维纳预测器的最佳解分别为
73
纯预测器的最小均方误差为 应用复卷积定理
74
取y(n)=x(n) 并考虑到b(n)是因果系统，得到
可以看到，随着N增加，E［|e(n+N)|2］min也增加。这一点也容易理解，当预测的距离越远，预测的效果越差，偏差越大，因而E［|e(n+N)|2］min越大。
75
例 已知 其中-1＜a＜1， 求: (1) 最小均方误差下的s(n+N)； (2) E［|e(n+N)|2］min。 ^ 解 首先对Sxx(z)进行功率谱分解。因为 所以
76
其次，求出B(z)的Z反变换 然后，应用Z变换的性质，得到 图 纯预测维纳滤波器
77
由Hopt(z)=aN，此时可以把纯预测的维纳滤波器看作是一个线性比例放大器（如图4.4.1所示）。根据x(n)的信号模型
x(n)=ω(n)+ax(n-1) 将信号x(n)通过纯预测维纳滤波器，随着时间的递增，可以得到 当N=1时，x(n+1)=ax(n)=as(n) 当N=2时，x(n+2)=ax(n+1)=a2s(n) 当N=N时，x(n+N)=ax(N+n-1)=aNs(n) …
78
以上推导结果相当于在n+N时刻，ω(n+N)=0，即去掉噪声时的结果。设N＞0时，ω(n+N)=0，则
x(n+N)=ax(n+N-1) 此时，从统计意义上讲，当N＞0时，白噪声信号ω(n+N)对x(n)无影响。 这一结论还可以推广，对于任何均值为零的x(n)，要估计s(n+N)时，只需要考虑B(z)的惯性，即可认为ω(n+N)=0，N＞0，这样估计出来的结果将有最小均方误差。 ^
79
终值定理 表明一个信号的功率谱在单位圆上没有极点与信号均值等于0等价，因此对于功率谱在单位圆上没有极点的信号，要估计s(n+N)时，可认为ω(n+N)=0, N＞0，即仅需要考虑B(z)的惯性，这样估计出来的结果将有最小均方误差。 ^
80
三、一步线性预测的时域解 令apk=-h(k)，则 预测误差
已知x(n-1), x(n-2)，…，x(n-p)，预测x(n)，假设噪声v(n)=0，这样的预测称为一步线性预测。设定系统的单位脉冲响应为h(n)，根据线性系统的基本理论，输出信号 令apk=-h(k)，则 预测误差
81
其中，ap0=1, 要使均方误差为最小值，要求 同维纳滤波的推导过程一样, 可以得到 E［e*(n)x(n-l)］=0 l=1, 2, …, p l=1, 2, …, p
82
由于预测器的输出 是输入信号的线性组合， 得到
说明误差信号与输入信号满足正交性原理预测误差与预测的信号值同样满足正交性原理。 预测误差的最小均方值
83
得到下面的方程组： 将方程组写成矩阵形式
84
这就是有名的Yule-Walker方程，可以看出Yule-Walker方程具有以下特点：
(1) 除了第一个方程外，其余都是齐次方程; (2) 与维纳-霍夫方程相比，不需要知道观测数据x(n)与期望信号s(n)的互相关函数。 该方程组有 p+1个方程，对应地，可以确定apk，k=1, 2, …， p 和E［e2(n)］min，共计 p+1个未知数，因此可用来求解AR模型参数。这就是后面要介绍的AR模型法进行功率谱估计的原理，它再一次揭示了时间序列信号模型、功率谱和自相关函数描述一个随机信号的等价性。
85
例 已知 x(n)为AR模型，求AR模型参数。 解 求解AR模型参数包括确定AR模型的阶数 p 及系数ap1，ap2，…，app。 首先对Sxx(z)做傅里叶反变换，得到x(n)的自相关函数rxx(m), rxx(m)=0.8|m| 采用试验的方法确定模型阶数 p。首先取 p=2，各相关函数值由上式计算，并代入(3.4.29) 式，可得
86
计算得到 a1=-0.8， a2=0, σ2ω=0.36 如果取p=3，可计算出a1=-0.8, a2=a3=0, σ2ω=0.36，说明AR模型的阶数只能是一阶的。采用谱分解的方法，即对Sxx(z)进行谱分解，得到的模型也是一阶的，其时间序列模型和差分方程为
87
4.5 卡尔曼(Kalman)滤波 卡尔曼滤波是用状态空间法描述系统的，由状态方程和量测方程所组成。卡尔曼滤波用前一个状态的估计值和最近一个观测数据来估计状态变量的当前值，并以状态变量的估计值的形式给出。卡尔曼滤波具有以下的特点： (1) 算法是递推的，且状态空间法采用在时域内设计滤波器的方法，因而适用于多维随机过程的估计；离散型卡尔曼算法适用于计算机处理。 (2) 用递推法计算，不需要知道全部过去的值，用状态方程描述状态变量的动态变化规律，因此信号可以是平稳的，也可以是非平稳的， 即卡尔曼滤波适用于非平稳过程。 (3) 卡尔曼滤波采取的误差准则仍为均方误差最小准则。
88
一、卡尔曼滤波的状态方程和量测方程 假设某系统k时刻的状态变量为xk，状态方程和量测 方程（也称为输出方程）表示为 (4.5.1a)
(4.5.1b) 其中，k表示时间，这里指第k步迭代时，相应信号的取值；输入信号ωk是一白噪声，输出信号的观测噪声 vk 也是一个白噪声，输入信号到状态变量的支路增益等于1，即B=1；
89
A表示状态变量之间的增益矩阵，可以随时间发生变化，用Ak表示第k步迭代时，增益矩阵A的取值；C表示状态变量与输出信号之间的增益矩阵，可以随时间变化，第 k 步迭代时，取值用Ck表示，其信号模型如图4.5.1所示。将状态方程中时间变量 k 用 k-1 代替，得到的状态方程和量测方程如下所示： xk=Ak-1xk-1+ωk-1 (4.5.2) yk=Ckxk+vk (4.5.3) 其中，xk是状态变量;ωk-1表示输入信号是白噪声; vk是观测噪声; yk是观测数据。
90
图 卡尔曼滤波器的信号模型
91
为了后面的推导简单起见，假设状态变量的增益矩阵A不随时间发生变化，ωk，vk都是均值为零的正态白噪声，方差分别是Qk和Rk，并且初始状态与ωk，vk都不相关， 表示相关系数。即
其中
92
二、标量新息过程及其性质 新息过程 为线性预测器在最小均方误差意义下的预测误差。
根据维纳滤波的正交原理（估计误差与输入信号向量正交）， 与输入信号向量 正交，因此， 包含了存在于当前观测样本 中的新的信息，“新息”的含义即在于此。
94
三 卡尔曼滤波的递推算法 当不考虑观测噪声和输入信号时，状态方程和量测方程为 (4.5.4) (4.5.5)
卡尔曼滤波是采用递推的算法实现的，其基本思想是先不考虑输入信号ωk和观测噪声vk的影响，得到状态变量和输出信号（即观测数据）的估计值，再用输出信号的估计误差加权后校正状态变量的估计值，使状态变量估计误差的均方值最小。 因此, 卡尔曼滤波的关键是计算出加权矩阵的最佳值。 当不考虑观测噪声和输入信号时，状态方程和量测方程为 (4.5.4) (4.5.5)
95
显然，由于不考虑观测噪声的影响，输出信号的估计值与实际值是有误差的，用 表示
(4.5.6) 为了提高状态估计的质量，用输出信号的估计误差 来校正状态变量 (4.5.7) 其中，Hk为增益矩阵，实质是一加权矩阵。经过校正后的状态变量的估计误差及其均方值分别用 和Pk表示，把未经校正的状态变量的估计误差的均方值用 表示
96
(4.5.8) (4.5.9) (4.5.10) 卡尔曼滤波要求状态变量的估计误差的均方值Pk为最小，因此卡尔曼滤波的关键就是要得到Pk与Hk的关系式，即通过选择合适的Hk,使Pk取得最小值。 首先推导状态变量的估计值 和状态变量的估计误差 ，然后计算 的均方值Pk ，并通过化简Pk，得到一组卡尔曼滤波的递推公式。
97
将（4.5.3）、 (4.5.5)式代入（4.5.7）式 同理，状态变量的估计误差 为 (4.5.11) (4.5.12)
98
由上式可以看出，状态变量的估计误差 由三部分组成， 可记为
其中 (4.5.13b) (4.5.13c) (4.5.13d) 那么，状态变量的估计误差的均方值Pk就由9项组成: (4.5.14a)
99
其中 (4.5.14b) (4.5.14d) (4.5.14c) 下面化简Pk的表达式，根据假设的条件，状态变量的增益矩阵A不随时间发生变化，起始时刻为k0，则（4.5.2）式经过迭代，得到 令l=k-k0-j，得到
100
取k0=0，k=k-1，得到 (4.5.15) 所以xk-1仅依赖于x0,ω0, ω1,…,ωk-2 , 与ωk-1不相关，即 (4.5.16) 又据(4.5.7)式和(4.5.3)式，得 (4.5.17)
101
所以 仅依赖于xk-1, vk-1，而与vk不相关，即
(4.5.18) (4.5.19) 把（4.5.15）～（4.5.19）式代入(4.5.14)式, Pk中的9项可以分别化简为 (4.5.20a) (4.5.20b)
102
(4.5.20c) (4.5.20d) (4.5.20e) (4.5.20f) (4.5.20g) (4.5.20h) (4.5.20j)
103
也就是说，Pk仅有其中的三项不为零， 化简成
(4.5.21)
104
为了进一步化简Pk，推导未经误差校正的状态估计误差的均方值Pk′，由下面推导结果可以看出，Pk′是一对称矩阵，满足Pk′=(Pk′)T。
(4.5.22)
105
将(4.5.22)式代入(4.5.21)式，即把Pk′代入Pk, (4.5.23) 其中, 是正定阵，记 (4.5.24) 令 (4.5.25) 将上式代入（4.5.23）式，得 (4.5.26)
106
将(4.5.26)式后三项配对 (4.5.27) 第二项和第三项均与Hk无关，第一项为一半正定阵，因此使Pk最小的Hk应满足 (4.5.28) (4.5.29)
107
将Hopt代入Pk，得到最小均方误差阵 (4.5.30) 将(4.5.7)、 (4.5.22)、 (4.5.29)式和(4.5.30)式联立， 得到一组卡尔曼递推公式 (4.5.31a) (4.5.31b) (4.5.31c) (4.5.31d)
108
假设初始条件Ak,Ck,Qk,Rk,yk, xk-1, Pk-1已知，其中 x0=E［x0］, P0=var［x0］, 那么，递推流程见图4
^ 图 卡尔曼滤波递推流程
109
例 已知 信号与噪声不相关，yk=xk+vk,求卡尔曼信号模型中的Ak和Ck。 解 由yk=xk+vk知道，Ck=1。 对Sxx(z)进行谱分解，确定x(n)的信号模型B(z)，从而确定Ak。 根据Sxx(z) =σ2ωB(z)B(z-1)，得出
110
上式与卡尔曼状态方程相比，不同之处在于输入信号ω(n)的时间不同，因此将Sxx(z)改写为
(解毕)
111
卡尔曼滤波和维纳滤波都是采用均方误差最小的准则来实现信号滤波的，但维纳滤波是在信号进入了稳态后的分析，卡尔曼滤波是从初始状态采用递推的方法进行滤波。对于平稳随机信号，当过渡过程结束以后，卡尔曼滤波与维纳滤波的结果间存在什么关系呢? 下面举一例说明。
112
例 已知 在k=0时开始观察yk, yk=xk+vk，用卡尔曼过滤的计算公式求xk， 并与维纳过滤的方法进行比较。 解 (1) 由x(n)功率谱及量测方程，确定卡尔曼递推算法。 首先对Sxx(z)进行功率谱分解，由例4.5.1的结果，得到卡尔曼滤波的状态方程为 xk=0.8xk-1+ωk-1, 确定Ak=0.8
113
由量测方程yk=xk+vk, 确定Ck=1， 将参数矩阵Ak,Ck,Rk代入卡尔曼递推公式（4.5.30）， 得到 (4.5.32a) (4.5.32b) (4.5.32c) (4.5.32d)
114
(2) 求出卡尔曼滤波的输出。由卡尔曼递推公式，以及
(2) 求出卡尔曼滤波的输出。由卡尔曼递推公式，以及 ，P0=var［x0］=1, 可得到Pk′,Hk,Pk及xk（k表示迭代次数），迭代流程为： 由具体迭代结果可以看出，原先的增益矩阵Ak,由于只选择了一个状态变量， 变成了加权系数。见表4.5.1。
115
表4.5.1 Kalman滤波迭代结果
116
(3) 求出卡尔曼滤波的稳态解。将(4.5.32b)式代入方程(4.5.32d), 得到第5个方程
(4.5.32e) 将方程(4.5.32c) 、(4.5.32e) 代入方程(4.5.32d), 消去Pk′， 可以得到Pk的递推关系: Pk=(1-Pk)［0.64 Pk ］ =0.64 Pk Pk -1 Pk Pk
117
化简上式，得到 1.36Pk+0.64Pk-1Pk=0.64Pk 要求的是稳态解，因此将Pk, Pk-1都用P∞代替, 得到
118
根据P∞，可以确定达到稳态后的卡尔曼滤波的状
态方程: (4.5.33)
119
(4) 用维纳滤波的方法分析。采用功率谱分解的方法，得到x(n)的时间序列信号模型的传输函数H(z):
120
上式说明x是一阶AR模型，对H(z)做Z反变换得到
121
(4.5.34) 比较（4.5.33）式和（4.5.34）式，可以看出卡尔曼滤波的稳态解与维纳解是相等的。 (解毕)
122
通过上面的例题，可以看出维纳滤波是已知前p个观测数据及信号与噪声的相关函数，通过建立模型的方法分析的。卡尔曼滤波要求已知前一个时刻的状态估计值x(k-1)和当前的观测值yk，由状态方程和量测方程递推得到结果。 维纳滤波的解以H(z)的形式给出，卡尔曼滤波是以状态变量的估计值给出解的形式。它们都采用均方误差最小的准则, 但卡尔曼滤波有一个过渡过程，其结果与维纳滤波不完全相同，但到达稳态后， 结果相同。 ^
123
三、 应用举例 下面举一个雷达跟踪目标物的例子说明卡尔曼滤波的应用。
雷达跟踪目标的基本原理是通过发射脉冲，根据接收到的脉冲与发射脉冲的时间间隔，来确定目标物的距离和速度。由于干扰的影响，接收到的脉冲波形变化很大，那么一次的测量结果可能存在很大的误差。为了减小误差，往往采取发射一串脉冲的方法进行测量。
124
我们知道，空间中的一点需要由径向距离和方位角来确定。 假设雷达跟踪的目标为飞行器，发射的脉冲时间间隔为T，在时间k，径向距离为R+ρ(k)，在时间k+1，距离为R+ρ(k+1)， 两者之间有秒的延时，因此T表示空间一次扫描的时间间隔。 平均距离用R表示，ρ(k)和ρ(k+1)表示对平均值的偏差。 假定偏差是统计随机的，均值为零, 那么可以写出距离方程 式中， 表示速度。用u表示加速度，则可以得到加 速度方程
125
假定加速度u(k)是零均值的平稳白噪声，即满足
为了后面叙述方便，定义x(k)表示第k个雷达回波脉冲获得的目标距离，z(k)表示在第k个雷达脉冲进行数据处理之后的目标距离估计，z(k)表示在第k个雷达脉冲进行数据处理之后的目标速度估计。设定状态变量为x(k)，选择的状态变量有4个，分别表示径向距离、径向速度、方位角和方位角速度，即 . （4.5.35）
126
根据状态变量的物理含义，得到以下方程： 式中，u1(k)和u2(k)表示在区间T径向速度和方向的变化。 将上式写成矩阵形式
127
由此得到卡尔曼滤波信号模型的状态方程 （4.5.36） 再来看量测方程， 距离和方向的估计值为 这里v1(k), v2(k)为观测偏差，将上两式分别写成向量形式和矩阵形式: Z(k)=CX(k)+V(k) （4.5.37a） （4.5.37b）
128
观测噪声V(k)假定为高斯噪声，均值为0，方差为
和 。状态方程激励信号的协方差阵为 (4.5.38) 其中， 为径向加速度在T时刻的方差； 为角加速度在T时刻的方差。
129
量测方程的噪声协方差阵 (4.5.39)
130
为了简单起见，假定在各个方向，加速度服从均匀分布，其概率密度函数 ，并将u的值限制在±M之间，那么，加速度的方差为
根据误差信号协方差阵P(k)的定义 P(k)=E［e(k)eT(k)］ 可以计算出，单个信号的均方误差和两个信号的协方差矩阵分别为
131
根据接收到的相邻两个回波脉冲，可以测量出飞行器的距离z1(1)和z2(2)，方位角z2(1)和z2(2)。根据这四个数据，用(4. 5
因此，k时刻的状态向量 可以写为
132
取k=2,将上式中的观测信号z1(1),z1(2),z2(1),z2(2)用(4.5.37)式代入，得到状态变量
（4.5.40） 根据(4.5.35)式, k=2时的状态向量为 （4.5.41）
133
计算k=2时刻的协方差阵 式中
134
因此，误差协方差阵是4×4阶矩阵，假设噪声源u和v是独立的， 则协方差阵为
其中
135
4.5.4 发散问题及其抑制 从理论上讲，卡尔曼滤波的递推算法可以无限地继续下去。然而在实际问题中的某些条件下，可能产生发散问题。也就是说，实际应用中发现估计误差大大地超过了理论误差的预测值，而且误差不但不减小，反而越来越大，即不收敛。 导致发散的一个原因是舍入误差的影响以及递推算法使得舍入误差积累的影响。计算机存贮单元的长度有限，使得舍入误差不可避免地存在，它相当于在状态方程和量测方程中又加入了噪声，带来的后果是有可能改变某些矩阵的性质，引起误差矩阵失去正定性和对称性。如果均方误差阵受到扰动而离开稳定解，只要它没有失去正定性，那么仍可能返回稳定解。
136
舍入误差引起的发散现象可以采用双精度运算得以改善，但运算量要增加许多。目前多采用平方根法，即把递推公式中的均方误差阵P改用其平方根P1/2实现。
另一种类型的发散问题是由于待估计过程模型的不精确引起的。人们在设计卡尔曼滤波时，认为分析过程是按某一规律发展的，但实际上是按另一规律演变的。如假定待分析过程的模型是一随机数，而实际过程是一个随机斜面，这样滤波器将连续地试着用错误曲线去拟合观测数据，结果导致发散。
137
当选择系统模型不准确时，由于新观测值对估计值的修正作用下降，陈旧观测值的修正作用相对上升，是引发滤波发散的一个重要因素。因此逐渐减小陈旧观测值的权重，相应地增大新观测值的权重，是抑制这类发散的一个可行途径。常用的方法有衰减记忆法、限定记忆法、限定下界法等。另外，通过人为地增加８模型输入噪声方差，用扩大了的系统噪声来补偿模型误差，抑制模型不准确所造成的发散现象，也是一种常见的策略，常用的方法有伪随机噪声法等。 还存在第三种发散问题，它是由于系统不可观察引起的。所谓不可观察，是指系统有一个或几个状态变量是隐含的，现有的观测数据不能提供足够的信息来估计所有的状态变量。这种发散问题表现为估计值误差不稳定或者均方误差阵的主对角线上有一项或几项无限增长。
Similar presentations
