维纳滤波和卡尔曼滤WienerandKalmanFiltering课件

第七章维纳滤波和卡尔曼滤

（WienerandKalmanFiltering）

随机信号或随机过程(randomprocess)是普遍存在的。一方面，任何确定性信号经过测量后往往就会引入随机性误差而使该信号随机化；另一方面，任何信号本身都存在随机干扰，通常把对信号或系统功能起干扰作用的随机信号称之为噪声。噪声按功率谱密度划分可以分为白噪声（whitenoise）和色噪声（colornoise），我们把均值为0的白噪声叫纯随机信号（purerandomsignal）。因此，任何其它随机信号都可看成是纯随机信号与确定性信号并存的混合随机信号或简称为随机信号。要区别干扰（interference）和噪声(noise)两种事实和两个概念。非目标信号（nonobjectivesignal）都可叫干扰。

干扰可以是确定信号，如国内的50Hz工频干扰。干扰也可以是噪声，纯随机信号（白噪声）加上一个直流成分（确定性信号），就成了最简单的混合随机信号。医学数字信号处理的目的是要提取包含在随机信号中的确定成分，并探求它与生理、病理过程的关系，为医学决策提供一定的依据。例如从自发脑电中提取诱发脑电信号，就是把自发脑电看成是干扰信号，从中提取出需要的信息成分。因此我们需要寻找一种最佳线性滤波器，当信号和干扰以及随机噪声同时输入该滤波器时，在输出端能将信号尽可能精确地表现出来。维纳滤波和卡尔曼滤波就是用来解决这样一类问题的方法：从噪声中提取出有用的信号。实际上，这种线性滤波方法也被看成是一种估计问题或者线性预测问题。设有一个线性系统，它的单位脉冲响应是，当输入一个观测到的随机信号，简称观测值，且该信号包含噪声和有用信号，简称信号，也即

(7-1)则输出为(7-2)如果该系统是因果系统，式(7-2)的m＝0，1，2，…，则输出的可以看成是由当前时刻的观测值和过去时刻的观测值、、…的估计值。用当前的和过去的观测值来估计当前的信号称为滤波；用过去的观测值来估计当前的或将来的信号，N，称为预测；用过去的观测值来估计过去的信号，N，称为平滑或者内插。本章将讨论滤波和预测问题。从图7-1的系统框图中估计到的信号和我们期望得到的有用信号可能不完全相同，这里用来表示真值和估计值之间的误差(7-3)显然是随机变量，维纳滤波和卡尔曼滤波的误差准则就是最小均方误差准则：(7-4)第一节维纳滤波器的时域解

（TimedomainsolutionoftheWienerfilter）

设计维纳滤波器的过程就是寻求在最小均方误差下滤波器的单位脉冲响应或传递函数的表达式，其实质就是解维纳－霍夫（Wiener－Hopf）方程。我们从时域入手求最小均方误差下的，用表示最佳线性滤波器。这里只讨论因果可实现滤波器的设计。一、因果维纳滤波器

设是物理可实现的，也即是因果序列：因此，从式(7-1)、(7-2)、(7-3)、(7-4)推导：

(7-5)要使得均方误差最小，则将上式对各，m＝0，1，…，求偏导，并且令其等于零，得：(7-6)(7-7)由式（7-9）进一步化简得：(7-10)二、有限脉冲响应法求解维纳－霍夫方程

如何去求解维纳－霍夫方程，即式（7-9）中解的问题，设是一个因果序列且可以用有限长（N点长）的序列去逼近它，则式(7-5)－(7-10)分别发生变化(7-11)

(7-12)

(7-13)………..(7-14)(7-15)简化形式：RxxH=Rxs

(7-17)式中，H＝[h(0)h(1)…h(N-1)]′，是待求的单位脉冲响应；Rxs＝′，是互相关序列；

Rxx＝

，是自相关矩阵

只要Rxx是非奇异的，就可以求到H：H=Rxx－1Rxs

(7-18)求得后，这时的均方误差为最小：

由式（7-15）进一步化简得：

(7-19)用有限长的来实现维纳滤波时，当已知观测值的自相关和观测值与信号的互相关时就可以按照式（7-15）在时域里求解但是当N比较大时，计算量很大，并

且涉及到求自相关矩阵的逆矩阵问题。

则式（7-15）和式（7-19）化为：(7-20)

(7-21)【例7-1】已知图7-1中

且与统计独立，其中

的自相关序列为

，是方差为1的单位白噪声，维纳滤波器来估计

，并求最小均方误差。

〖解〗依题意，已知信号的自相关和噪声的自相关为：

，代入式（7-20）得解得：

＝0.451，

＝0.165。

将上述结果代入式（7-21），求得最小均方误差：若要进一步减小误差可以适当增加维纳滤波的阶数，但相应的计算量也会增加。三、预白化法求解维纳－霍夫方程

从上面分析知求解维纳－霍夫方程比较复杂，本节用波德（Bode）和香（Shannon）

提出的白化的方法求解维纳－霍夫方程，得到系统函数

H(z)。由第三章的知识，我们知道随机信号都可以看成是由一白色噪声

w1(n)

激励一个物

理可实现

的系统

或模型的响应，如图7-2所示，

其中A(z)表示系统

的传递函数。由于

x(n)=s(n)+w(n)，

在图7-2的基础上给出x(n)的信号模型

白噪声的自相关函数为，它的z变换就等于。图7-2中输出信号的自相关函数为，根据卷积性质有令

上式

令

对式（7-22）进行Z变换得到系统函数和相关函数的z变换之间的关系：

(7-23)同样，对图7-4进行z变换得(7-24)图7-4中利用卷积性质还可以找到互相关函数之间的关系：两边z变换得到(7-25)如果已知观测信号的自相关函数，求它的z变换，

然后找到该函数的成对零点、极点，取其中在单位圆内的那一半零点、极点构成，另外在单位圆外的零、极点构成，这样就保证了

是因果的，并且是最小相位系统。将图7-1重新给出，待求的问题就是最小均方误差下的最佳，如图7-5(a)所示，为了便于求这个，将图7-5(a)的滤波器分解成两个级联的滤波器：和G（z），

如图7-5(b)所示，则

(7-27)（a）（b）图7-5利用白化方法求解模型有了上述的模型后，白化法求解维纳－霍夫方程步骤如下：1）对观测信号的自相关函数求z变换得到；2）利用等式，找到最小相位系统；3）利用均方误差最小原则求解因果的G（z）；4），即得到维纳－霍夫方程的系统函数解。在上述步骤中，可以通过已知的观测信号的自相关函数来求得，因而求解的问题就归结为求解G（z）的问题了。由于G（z）的激励源是白噪声，求解变得容易多了，下面我们分析步骤3的求解过程。按图7-5(b)有:

(7-28)均方误差为：由于，代入上式，并且进行配方得：均方误差最小也就是上式的中间一项最小，所以

(7-30)注意，这里的是因果的。对该式求z变换，得到

(7-31)表示对求单边z变换。所以维纳－霍夫方程的系统函数解表示为：由式（7-25）上式可以表示为：因果的维纳滤波器的最小均方误差为：

(7-33)利用帕塞伐尔(Parseval)定理，上式可用z域来表示：(7-34)利用帕塞伐尔(Parseval)定理，上式可用z域来表示：

围线积分可以取单位圆。(7-34)【例7-2】已知图7-1中，且统计独立，其中的自相关序列为，是方差为1的单位白噪声，试，并求最小均方误差。与设计一个物理可实现的维纳滤波器来估计〖解〗依题意，已知，，，步骤1求z变换步骤2由于，容易找到最小相位系统和白噪声方差步骤3利用式（7－32）

对括号里面求反变换，注意括号内的收敛域为取因果部分，也就是第一项，所以步骤4最小均方误差为取单位圆为积分围线，有两个单位圆内的极点，0.8和0.5，求它们的留数和，所以第二节维纳预测器（Wiener’sPredictor）

上节讨论的维纳滤波器是一种估计器，是用观测到的当前和全部过去的数据、、…来估计当前的信号值。本节将讨论维纳预测器，它同样也是一种估计器，是用过去的观测值来估计当前的或将来的信号，N，也是用真值和估计值的均方误差最小为估计准则。一、因果的维纳预测器

图7-6就是维纳预测器的模型，N>0，是希望得到的输出，而表示实际的估计值。图7-6维纳预测器本节和上节一样着重讨论预测器的系统函数以及预测的均方误差，维纳预测器和维纳滤波器比较类似，因而分析方法也都可以借鉴前面的内容。对于图7-6模型，设是物理可实现的，也即，则有是因果序列：

(7-35)

(7-36)要使得均方误差最小，则将上式对各，m＝0，1，…，求偏导，并且等于零，得：

(7-37)即用相关函数R来表达上式：(7-38)(7-39)由于，则，z变换得(7-40)借鉴维纳滤波器的结果类似给出维纳预测器的最佳传递函数，对应维纳预测器，对应维纳滤波器，故因果的预测器的传递函数为：

(7-41)最小均方误差为(7-42)利用帕塞伐尔(Parseval)定理，上式可用z域来表示【例7-3】已知图7-6中，且与统计独立，其中的自相关序列为，是方差为1的单位白噪声，，并求最小均方误差。试设计一个物理可实现的维纳预测器估计〖解〗依题意已知，，，

求z变换：由于，容易找到最小相位系统和白噪声方差：由式（7-41），N＝1，求z变换：由于，容易找到最小相位系统和白噪声方差：由式（7-41），N＝1，对括号里面求z反变换，注意括号内的收敛域为：，取因果部分，也就是第一项，所以把上式写成差分方程形式有：最小均方误差为：二、纯预测器（N步）

纯预测器指的是＝0的情况下，对的预测。如图7-7所示。图7-7N步纯预测器这时，用白化法来求解预测器的系统函数。因为，从而有：

(7-44)将上式代入式（7-41）、（7-43）得：假设B(z）是b（n）的z变换，且b（n）是实序列，则上式可以利用帕塞伐尔定(Parseval)理进一步化简：

(7-46)又因为B(z）是最小相位系统，一定是因果的，上式可以简化

(7-47)上式说明最小均方误差随着N的增加而增加，也即预测距离越远误差越大。【例7-4】已知图7-7中，其中的自相关序列为，试设计一个物理可，并求最小均方误差。，则实现的维纳预测器来估计〖解〗依题意，已知因为

容易找到最小相位系统和白噪声方差:利用式（7-45）:因为

，只取的部分，有:回到z域有：，代入得:最小均方误差为：它说明当N越大，误差越大，当N＝0时，没有误差。把上述结果用模型表示如图7-8所示。图7-8例题7-3的纯预测模型三、一步线性预测器

对于纯预测问题，有，然而预测的问题常常是要求在过去的p个观测值的基础上来预测当前值，也就是这就是一步线性预测公式，常常用下列符合表示

(7-48)式中p为阶数，。预测的均方误差为:(7-49)要使得均方误差最小，将上式右边对求偏导并且等于零，得到p个等式

(7-50)最小均方误差：

(7-51)式（7-50）就是Yule－Walker（Y-W）方程，和第三章AR模型参数估计的方程一致，如何去求解该方程在第三章有详细介绍。把Yule－Walker（Y-W）方程和维纳－霍夫方程进行比较，维纳－霍夫方程要估计的量是s(n)，Y-W方程要估计的量是x(n)本身，因而解维纳－霍夫方程要已知x(n)、y(n)的互相关函数，实际中这个互相关函数往往是未知的，而解Y-W方程只需要知道观测信号的自相关函数。因此Y-W方程比W-H方程更具有实用价值。例7-5】已知图7-7中x(n)=s(n)，其中的自相关序列为的可实现的一步线性预测器，并求最小均方误差。〖解〗，

，试设计一个p＝2利用Y-W方程，可以列出2个方程式解得：，也即结果和例（7-4）N=1时一致。第三节维纳滤波器的应用

（ApplicationofWienerFilter）要设计维纳滤波器必须知道观测信号和估计信号之间的相关函数，即先验知识。如果我们不知道它们之间的相关函数，就必须先对它们的统计特性做估计，然后才能设计出维纳滤波器，这样设计出的滤波器被称为“后验维纳滤波器”。在生物医学信号处理中比较典型的应用就是关于诱发脑电信号的提取。大脑诱发电位（EvokedPotential，EP）指在外界刺激下，从头皮上记录到的特异电位，它反映了外周感觉神经、感觉通路及中枢神经系统中相关结构在特定刺激情况下的状态反应。在神经学研究以及临床诊断、手术监护中有重要意义。EP信号十分微弱，一般都淹没在自发脑电（EEG）之中，从EEG背景中提取诱发电位一直是个难题：EP的幅度比自发脑电低一个数量级，无法从一次观察中直接得到；EP的频谱与自发脑电频谱完全重迭，使得频率滤波失效；在统计上EP是非平稳的、时变的脑诱发电位。通过多次刺激得到的脑电信号进行叠加来提取EP，这是现今最为广泛使用的EP提取方法。为了解决诱发电位提取问题，研究者利用维纳滤波来提高信噪比，先后有Walter、Doyle、Weerd等对维纳滤波方法进行了改进。在频域应用后验维纳滤波的核心就是由各次观察信号中分解出信号的谱估计和噪声的谱估计，通过设计出的滤波器来提高信噪比。本节将介绍时－频平面的维纳滤波（time－frequencyplanewienerfiltering，简称TFPW）在高分辨心电图（HRECG）中的应用。方法如下：一、观测样本设共有N次观测样本：xi(t)=s(t)+wi(t),i=1,2,…N。其中s(t)是周期确定的心电信号；wi(t)是第i次记录时的噪声，包括肌电、测量仪器噪声等，假设每次记录的噪声之间互不相关；xi(t)是观测信号；信号和噪声相互独立。对每次观测用短时傅立叶变换求时频表示（TFR）：对N次观测的时频表示（TFR）求平均：，样本平均的时频表示（TFR）为：

(1)样本平均为：

(2)从式（2）可以得到一个基于样本平均的简单时－频平面后验维纳滤波器：

(3)二、公式修正在时－频域上对式（1）（2）进行修正，给出更实际的表示：

………………..……(4)

(5)式中COV表示信号和噪声之间的方差，也就是考虑了信号和噪声并非相互独立；IF是干扰项；表示样本平均的噪声功率；表示样本噪声功率的平均。三、TFPW的计算过程TFPW的计算过程如图7-9所示。图7-10TFPW的模拟实验结果

注：(上图)原信号是两个正弦波，观测信号混有白噪声；(下图)原信号是线性调频信号，观测信号混有白噪声。

在图7-10中每一个图中从上至下分别表示：测量的单个样本，样本平均，TFPW滤波器估计的信号，原始信号。图7-10的初始信噪比设为－12dB，TFPW与叠加平均法相比，信噪比有5个dB左右的改善。五、需要进一步研究的问题FPW滤波中由于有二次TFR中的相关噪声以及IF项，滤波器可能包含虚部，也就是包含信号的相位信息，直接在时－频平面上考虑相位问题还需要进一步研究。第四节卡尔曼滤波的信号模型

（SignalModelofKalmanFiltering）

通过前面几节内容的学习，我们知道维纳滤波是根据当前和过去全部的观测值来估计信号的当前值，它的解形式是以均方误差最小为原则下的系统的传递函数或单位脉冲响应。而卡尔曼滤波不需要过去全部的观测它是根据前一个估计值和最近一个观测值来估计信号的当前推方法进行估计的，因而卡尔曼滤波对信号的平稳性和时不变性不做要求。我们利用维纳滤波的模型引入到卡尔曼滤波的信号模型。，它是用状态方程和递一、状态方程和量测方程要给出卡尔曼滤波的信号模型，先来讨论状态方程和量测方程。图7-11是维纳滤波的模型，信号可以认为是由白噪声激励一个线性系统的响应，假设响应和激励的时域关系可以用下式表示：

(7-52)上式也就是一阶AR模型。在卡尔曼滤波中信号被称为是状态变量，用矢量的

形式表示为，在k时刻的状态用表示，

在k－1时刻的状态用表示。

激励信号也用矢量表示为，激励和响应之间的关系用传递矩阵来表示，它是由系统的有一定关系。有了这些假设后结构确定的，与我们给出状态方程：(7-53)上式表示的含义就是在k时刻的状态可以由它的前一个时刻的状态来求得，即认为k－1时刻以前的各状态都已记忆在状态中了

图7-11维纳滤波的信号模型和观测信号模型卡尔曼滤波是根据系统的量测数据（即观测数据）对系统的运动进行估计的，所以除了状态方程之外，还需要量测方程。还是从维纳滤波的观测信号模型入手，图7-11的右图，观测数据和信号的关系为：，一般是均值为零的高斯白误差矢量，则量测矢量与状态矢量(7-54)上式和维纳滤波的概念上是一致的，也就是说卡尔曼滤波的一维信号模型和维纳滤波的信号模型是一致的。把式(7-54)推广就得到更普遍的多维量测方程(7-55)上式中的称为量测矩阵，它的引入原因是，的维数不一定与状态矢量的维数相同，因为我们不一定能观测到所有需要的状态参数

量测矢量假如是的矢量，是的矢量，就是的矩阵，是的矢量。二、信号模型

有了状态方程和量测方程后我们就能给出卡尔曼滤波的信号模型，如图7-12所示。图7-12卡尔曼滤波的信号模型【例7-6】设卡尔曼滤波中量测方程为，已知信号的自相关函数的z变换为，噪声的自相关函数为：，信号和噪声统计独立。求卡尔曼滤波和。信号模型中的〖解〗根据等式：可以求得：变换到时域得：因此又因为，所以＝1。第五节卡尔曼滤波方法

（MethodofKalmanFiltering）

建立好了卡尔曼滤波的信号模型以及状态方程、量测方程后，要解决的问题就是要寻找在最小均方误差下信号的估计值。一、卡尔曼滤波的一步递推法模型

把状态方程和量测方程重新给出：(7-56)(7-57)上式中和是已知的，已知，现在的问题就是如何来求当前时刻。

是观测到的数据，也是已知的，假设信号的上一个估计值的估计值上两式中如果没有与，可以立即求得，估计问题的出现就是因为信号与噪声的与，用上两式和分别用和表示，得：叠加。假设暂不考虑得到的(7-58)

(7-59)必然，观测值和估计值之间有误差

，它们之间的差称为新息（innovation）：

(7-60)显然，新息的产生是由于我们前面忽略了与所引起的，也就是说新息里面

包含了与的信息成分。因而我们用新息乘以一个修正矩阵，用它来代替式（7-56）的来对进行估计：(7-61)由（7-56）～（7-61）可以画出卡尔曼滤波对进行估计的递推模型，如图7-13所示，输入为观测值，输出为信号估计值。图7-13卡尔曼滤波的一步递推法模型二、卡尔曼滤波的递推公式从图7-13容易看出，要估计出就必须要先找到最小均方误差下的修正矩阵，结合式（7-61）、（7-56）、（7-57）得：……….(7-62)根据上式来求最小均方误差下的，然后把求到的代入（7-61）则可以得

到估计值。设真值和估计值之间的误差为：，误差是个矢量，因而均方误差是一个矩阵，用表示。把式（7-62）代入得：….(7-63)均方误差矩阵：(7-64)表示对向量取共轭转置。为了计算方便，令(7-65)找到和均方误差矩阵的关系：(7-66)把式（7-63）代入式（7-64），并且利用条件：与都是零均值的高斯白噪声，且它们和

互不相关，协方差矩阵分别为与不相关；与及不相关。最后化简得：………………….(7-67)把式（7-66）代入（7-67）得令，，代入上式化简：

(7-68)上式第一项和第二项与修正矩阵无关，第三项是，于是可以求得最小均方误差下的修正矩阵为：半正定矩阵，要使得均方误差最小，则必须(7-69)把上式代入(7-61)即可得均方误差最小条件下的递推公式。相应的式（7-68）的第三项为零，得最小均方误差为：(7-70)综上所述，得到卡尔曼滤波的一步递推公式：(7-71)(7-72)(7-73)(7-74)有了上面四个递推公式后我们就可以得到和。如果初始状态的统计特性已知，并且令且矩阵

都是已知的，以及观测量也是已知的，就能用递推

计算法得到所有的和：将初始条件代入式（7-71）求得；

将代入式（7-72）求得和代入式（7-73）求得；将初始条件和代入式（7-74）求得；…依此类推。这样递推用计算机实现

；将非常方便。和维纳滤波一样，卡尔曼滤波也可以推广到卡尔曼预测，推导过程和维纳滤波到维纳预测类似，也同样有纯卡尔曼预测，这里不再推导。【例7-7】设卡尔曼滤波中量测方程为，已知信号的自相关函数的z变换为，噪声的自相关函数为，，信号和噪声统计独立，已知，在k＝0时刻开始观测信号。试用卡尔曼滤波的公式求和，k＝0,1,2,3,4,5,6,7；以及和。稳态时的〖解〗由例7-6的结果知，＝1，，，把它们代入式(7-71)～(7-74)得（1）（2）

（3）（4）由于是一维情况，求逆，把（1）代入（2）、（3）式，消去，再把（2）和（3）联立，得到

（5）初始条件为，k＝0开始观测，利用等式（4），（5）进行递推得：k＝0，1.0000，1.0000，；k＝1，0.5000，0.5000，k＝2，0.4048，0.4048，

k＝3，0.3824，0.3824，

k＝4，0.3768，0.3768，

k＝5，0.3755，0.3755，

k＝6，0.3751，0.3751，；k＝7，0.3750，0.3750，

如果给定每个时刻的观察值就可以得到每一时刻的信号估计值，上面是递推过程，还没有达到稳态的情况。假设到了某一时刻k－1，前后时刻的均方误差相等，也就是误差不再随着递推增加而下降，达到最小的均方误差了，即稳态情况，式（5）中的误差，代入（5）式可以计算到稳态时的均方误差为

即稳态时的修正矩阵，代入式（4）得稳态时的信号估计：

化到z域有：。将上述结果和维纳滤波的例题7-2的结果相比较：，，，发现当卡尔曼滤波达到稳态时和维纳滤波的结果一致，原因就是它们两种滤波都是用的同样的估计原则：最小均方误差准则。然而在卡尔曼滤波的过渡期间的信号估计结果和维纳滤波当然完全不同。第六节卡尔曼滤波器的应用

(ApplicationKalmanFilter)

最优估计指从带有随机干扰的观测数据中估计出信号来，其中的线性最小均方误差的卡尔曼滤波占有重要的地位，自动控制系统中应用非常广泛。前面我们已经推导出卡尔曼滤波的公式，也有了卡尔曼滤波器设计的直接调用程序。应用卡尔曼滤波时，核心是把问题如何纳入卡尔曼滤波的框架里面去，往往很难获得准确可靠的噪声数据，如前面的和加上干扰和噪声的协方差矩阵不一定为零，有，一旦确定了这几个，H，]＝kalman表示均方误差矩阵，表示另外一种均方误差矩阵。矩阵，可以直接调用Matlab中的控制工具箱中的状态空间设计函数：[S，L，（sys，Q，R，N）。输入变量的含义与上面提到的相同，sys表示状态空间模型，可以用函数ss（a，b，c，…）来生成。输出变量S表示卡尔曼滤波器的状态方程的模型，L表示滤波器增益矩阵（是为了计算而定义的），H表示修正矩阵，下面用例题来看实际的计算过程。【例7-8】已知条件和例7-7一样，状态方程和量测方程为：其中，，，，信号和噪声统计独立。求卡尔曼滤波器的稳态和。〖解〗根据函数调用sys＝ss（A，B，C，D，1），得到离散卡尔曼状态模型，采样周期这里设为1。A，C已知，由于函数调用中是设计了两个观测信号的，我们这里只有一个观测信号，所以B取[01]，后一个1表示噪声的