第三讲 序贯蒙特卡罗滤波

非线性/非高斯滤波讲义
i i = xk Pr ( zk | xk x= | z1:k −1 ) k , z1:k −1 ) Pr ( xk = Pr ( zk | z1:k −1 )
=
i i = xk Pr | z1:k −1 ) ( zk | xk x= k ) Pr ( xk
Pr ( zk | z1:k −1 )
i ˆ ( x0:k | z1:k ) = ∑ ω kjδ ( x0:k − x0 p :k ) j =1 N
i 来 逼 近 真 实 的 系 统 S ∈ S SSM 的 条 件 分 布 p ( x0:k | z1:k ) 。 其 中 x0:k 为 第 i 条 样 本 路 径 ，
i = 1, , N 。我们形象地称在 k 时刻的 xk 的一个随机抽样为“粒子” ，所有 k 时刻由蒙特卡
Pr ( x = ∑
j =1 Ns
= ∑ aij ωkj−1|k −1
j =1
ωki |k −1
故 xk 的预测概率分布密度为
= p ( xk | z1:k −1 )
∑ω
i =1
Ns
i k | k −1
i δ ( xk − xk )
当 zk 给定时，利用贝叶斯公式，我们有
张永安
第三讲
i Pr ( xk = xk | z1:k )
x∈=
g ( x ) p ( x ) dx E p ( x ) g ( x )
ˆ ( x) 近似，其数学期望可以表示为 若 p( x) 形式较复杂导致以上积分计算困难时，我们采用 p E g ( x ) ≈∫
x∈=
ˆ ( x ) dx = g ( x) p Ep ˆ ( x) g ( x )
类似我们同样可以求得最大验后估计。 对于非线性系统，从定理 2.1 的递推贝叶斯估计公式可以看到，归一化常数 p ( z k | z1:k ) 一般来说是未知的，因此状态的后验概率密度 p ( x0:k | z1:k ) 是不可直接抽样的，但我们可以 通过对另外一个与其具有相同或者更大支撑集合的概率密度函数 q ( x0:k | z1:k ) (称为重要性 抽样函数或重要性函数)进行抽样，也就是说
-2-
张永安
第三讲
非线性/非高斯滤波讲义
3.2 抽样原理
ˆ ( x) 为 p ( x) 的一个近似， ˆ ( x) 也是一个 pdf， 设 p( x) 是一个概率分布密度函数 （pdf） ，p 且p
g ( x ) 为 上的一个有界函数，则 g ( x ) 的数学期望为
E g ( x ) =∫
j =1 Ns
i x∈xk
p ( x | xkj−1 ) dx
ω
i k |k

ωki |k −1 ∫
i x∈xk
p ( zk | x ) dx p ( zk | x ) dx
∑ω
j =1
Ns
j k | k −1 x∈xi k
∫
这里 xkj−1 表示 k − 1 时刻第 j 个网格单元的中心。为了简化计算，我们通常用网格中心点的概 率密度来代替整个网格单元的概率加权值，这样
ω ki ∝
因而我们可以将 p ( x 0:k | z1:k ) 近似估计为
i p ( x0 :k | z1:k ) i q ( x0:k | z1:k )
ˆ ( x0:k | z1:k ) = ∑ ω kjδ ( x0:k − x0j:k ) p
j =1
N
若选取 q ( x0:k | z1:k ) 满足如下形式：
ˆ ( x) 。 下面， 我们利用抽样来构造 p 我们从 p( x) 抽取 N s 个相互独立的样本 { xi | i = 1, 2, , N s } ，
我们用这组样本来代替整个总体，样本的概率分布密度函数为
ˆ ( x) = p 1 Ns
∑δ ( x − x )
i i =1
Ns
它是 p( x) 的一个近似。根据强大数定律，
i
{
}
q ( x0:k | z1:k ) 的一组样本，则 q ( x0:k | z1:k ) 可以近似为
i q ( x0:k | z1:k ) ≈ ∑ δ ( x0:k − x0 :k ) k =1 N
因此，我们有
p ( x0:k | z1:k ) =
其中 ω k 为权函数，满足
i
N p ( x0:k | z1:k ) q ( x0:k | z1:k ) ≈ ∑ ω kjδ ( x0:k − x0j:k ) q ( x0:k | z1:k ) j =1
张永安 非线性/非高斯滤波讲义
第三讲 序贯蒙特卡罗滤波
3.1 近似网格滤波
假定系统是隐马尔科夫链，且状态空间是有限集，也就是系统状态空间 可以表示为
N 2 = = { x1 k , xk , , xk
s
}
这样，在任意时刻，其概率分布是一个离散型分布。设 k − 1 时刻，系统状态的后验分布为
i = Pr ( xk −1 x= ωki −1|k −1 ， i = 1, 2, , N s k −1 | z1:k −1 )
这里
ω(x ) =
i
q ( xi )
p ( xi )
k 为归一化常数，满足
-3-
张永安
第三讲
Ns
非线性/非高斯滤波讲义
ˆ ( x ) dx =1 ⇒ ∑ ω ( x ) =1 ⇒ k = ∫p N
i s i =1
k
Ns
∑ω ( x )
i i =1
Ns
故有
ˆ ( x) = p
∑ω δ ( x − x )
i δ ( xk − xk )
i 若系统状态空间是连续的，但可以划分为 N s 个网格单元 { xk 则其后验分 | i = 1, 2, , N s } ，
布密度可以近似为一个离散型分布：
i p ( xk −1 | z1:k −1 ) ≈ ∑ ωki −1|k −1δ ( xk −1 − xk −1 ) i =1 Ns
ˆ ( x) q = 1 Ns
∑δ ( x − x )
i i =1
Ns
从而
i p ( x) k Ns p ( x ) k Ns i ˆ ( x) k= ˆ ( x) p q x x = δ − ω ( xi )δ ( x − xi ) ( ) ∑ q xi ∑ q ( x) N s i 1= N 1 i = ( ) s
来表征。 根据全概率公式，我们有
i Pr ( xk = xk | z1:k −1 )
= =
Pr ( x = ∑
j =1 Ns
Ns
i x= xkj−1 , z1:k −1 ) Pr = ( xk −1 xkj−1 | z1:k −1 ) k | xk −1 i x= xkj−1 ) Pr = ( xk −1 xkj−1 | z1:k −1 ) k | xk −1
则系统在 k 时刻状态的预测概率分布密度和修正后后验分布密度可以近似为
= p ( xk | z1:k −1 ) = p ( xk | z1:k )
∑ω
i =1 Ns
Ns
i k | k −1
i δ ( xk − xk )
∑ω
i =1
i k |k
i δ ( xk − xk )
其中
ωki |k −1 ∑ ωkj−1|k −1 ∫
存在，我们就可以求出其近似最小均方估计
-4-
张永安
第三讲
非线性/非高斯滤波讲义
ˆ k ( x0:k ) = ∫ g k ( x0:k ) p ˆ ( x0:k | z1:k )dx1:k = g
则根据强大数定律，有
1 N
∑g
j =1
N
k
( x0j:k )
.s. ˆ k ( x0:k ) a g → E{g k ( x0:k )} N →∞
。为了便于后面问题的叙述，我们称随机 罗随机抽样所获得的 xk 的样本集合称为“粒子云” 过程蒙特卡罗抽样样本路径中非零概率权值、且取值互不相等的一组样本路径为有效样本； 类似地，粒子云中非零概率权值，且取值互不相等的一组粒子为有效粒子。对于任意的非线 性函数 g k ( x0:k ) ，若其期望
E{g k ( x0:k )} = ∫ g k ( x0:k ) p( x0:k | z1:k )dx0:k
i i i =1
Ns
其中
ω =
i
ω ( xi )
∑ω ( x )
j j =1
Ns
3.3 序贯蒙特卡罗滤波（粒子滤波）
我们注意到， 前面的近似高斯滤波一个共同点是， 在对状态的条件概率分布的递推计算 过程中， 仅对条件分布的一阶矩和二阶矩做近似递推计算。 因而它们的适应范围只是随机非 线性系统的特殊情况： S ∈ S GSSM 。并且，到现在为止，关于高斯滤波还没有宽松的收敛性保 证。这是因为，在非线性、非高斯系统中，状态的条件分布实际上都有无穷个参数，这时， 仅在递推参数中传递两个低阶矩参数或者高斯分布是不够的。也就是说，在非线性/非高斯 系统模型，采用高斯假设往往是难以得出满意的结果。一种思路是采用非参数的方法，完全 放弃对状态的条件分布所作的高斯假设， 这样， 滤波的适应范围才能扩大到一般的状态空间 模型集合 S SSM 。粒子滤波就是这样一种方法。不同于前面几种高斯滤波方法，粒子滤波利 用序贯蒙特卡罗模拟的办法来近似整个条件概率分布。 因此， 粒子滤波是一种全局近似最优 滤波。 其基本思想是用一组简单随机样本来近似随机量的概率分布， 也就是用随机离散型概 率分布列
p ( x0:k | z1:k ) > 0 ⇒ q ( x0:k | z1:k ) > 0
(3-1)
条件(3-1)是为了保证 p ( x0:k | z1:k ) 的每一样本在对重要性函数 q ( x0:k | z1:k ) 抽样过程中都有 可能出现。在该条件下，我们同样可以对后验分布进行近似。设 x0:k : k = 1,2, , N 为
i 的概率： 素 aij , i, j = 1, 2, , N s 定义为状态从 xkj−1 转移到 xk
j Pr= xk −1 x= aij ,= i, j 1, 2, , N s ( xk xki | = k −1 )
观测过程则由似然概率
Pr ( zk | xk = xkj ) ， j = 1, 2, , N s
q ( x0:k | z1:k ) = q ( xk | x0:k −1 , z1:k )q ( x0:k −1 | z1:k −1 )
则根据贝叶斯递推估计公式，我们有-Biblioteka -(3-2)张永安
第三讲
非线性/非高斯滤波讲义
则系统状态的后验分布密度可以写成：
= p ( xk −1 | z1:k −1 )
∑ω
i =1
Ns
i k −1| k −1
i δ ( xk −1 − xk −1 )
设系统的状态过程是一 Markov 链，则
Pr ( xk | x1:k −1 ) = Pr ( xk | xk −1 )
由于系统状态空间是有限的， 故系统的转移概率可以用一个矩阵来描述， 其中矩阵的每个元
q ( x) > 0 ⇒ p ( x) > 0 ；
(2) 从重要性函数 q ( x ) 抽取 N s 个相互独立的样本 { xi | i = 1, 2, , N s } ；
ˆ ( x ) 计算 p ( x) 的近似函数 p ˆ ( x) (3) 由 q ( x ) 的近似 q ˆ ( x ) 的近似函数可以由样本集给出： q
N s →∞ a.s. lim E p ˆ ( x) g ( x ) → E p ( x ) g ( x )
若 p( x) 不能直接抽样，我们这时需要构造新的抽样函数，其方法如下： (1) 构 造 一 个 与 p( x) 相 近 的 概 率 分 布 密 度 函 数 q ( x ) （ 称 为 重 要 性 函 数 ） ，满足
ω
i k |k
这里 Pr ( zk | z1:k −1 ) 为归一化常数
= Pr ( z k | z1:k −1 ) = (z | x ∑ Pr
i =1 k k Ns i i x= xk | z1:k −1 ) k ) Pr ( xk
因而
= p ( xk | z1:k )
∑ω
i =1
Ns
i k |k
ωki |k −1 ≈ ∑ ωkj−1|k −1 p ( xki | xkj−1 )
j =1 Ns
ωki |k ≈
ωki |k −1 p ( zk | xki )
∑ω
j =1
Ns
j k | k −1
p ( zk | xki )
当网格无数增加时， 近似网格滤波为一种最优滤波的解决方法。 但实际上网格数随着维数的 增加而指数增长，因此在实际的最优滤波中是很难实现的。为此，我们以下利用重要性抽样 算法来克服该困难。