几乎任何递推和迭代算法都有发散的可能,卡尔曼滤波也不例外。造成卡尔曼滤波发散的原因无外乎以下两点:
状态方程描述的动力学模型不准确,或者噪声的统计模型不准确,这样会使模型和量测值不匹配,导致发散。其实模型不准确的情况很少,除非原理上错误。更多的是由于模型太过粗糙而导致的。
递推过程计算机摄入误差累积,使得方差阵失去正定性或者失去对称性,导致增益计算失去加权效果,从而导致滤波器发散。
先说模型不准确导致的发散处理技术,累积误差导致发散的处理后面有机会再提。
首先还是要先理解一下,卡尔曼滤波器的本质是线性最小方差估计,它虽然是递推形式的,可是每一步递推获得的估计值其实是根据当前时刻之前的所有量测值而获得的最优估计结果。也就是说它实际是有记忆的。
而量测值其实是完全不会发散的东西,因为它会始终在它的数学期望周围分散,我们的思路就是,在滤波器中,我们应该更加相信新鲜的量测值,而减小对陈旧量测值的信任程度,也就是增大新鲜量测值的权重,减少陈旧量测值的权重。根据这一思路,我们可以重新设计卡尔曼滤波器。其设计思路有两种,分别是衰减记忆法和限定记忆法
一、衰减记忆法
再来看看量测方程:
Z
k
=
H
k
X
k
+
V
k
Z_k=H_kX_k+V_k
Zk=HkXk+Vk 根据上面的思路,我们认为,当前时刻
Z
N
Z_N
ZN是最准的,而之前的时刻的
Z
k
(
k
<
N
)
Z_k(k<N)
Zk(k<N)要稍微不准一些,离N越远,说明量测的数据陈旧,量测值也越不准。
不准是什么概念,就是那个
R
k
R_k
Rk大,也就是说,我们希望离当前时刻越远的
R
k
R_k
Rk最好越大。这好办,把
R
k
R_k
Rk乘个系数好了
R
N
=
R
k
s
N
−
k
R_N=R_ks^{N-k}
RN=RksN−k 其中s>1,这样不就实现,在当前时刻R最小,之前的时刻R都要大一些。
基于这个思路,重新对卡尔曼滤波器进行推导,就会得出下面的滤波方程(推导过程不想写了!):
P
k
/
k
−
1
∗
=
Φ
k
,
k
−
1
(
s
P
k
−
1
∗
)
Φ
k
,
k
−
1
T
+
Q
k
−
1
P_{k/k-1}^*=\Phi_{k,k-1}(sP_{k-1}^*)\Phi_{k,k-1}^T+Q_{k-1}
Pk/k−1∗=Φk,k−1(sPk−1∗)Φk,k−1T+Qk−1
K
k
∗
=
P
k
/
k
−
1
∗
H
k
T
[
H
k
P
k
/
k
−
1
∗
H
k
T
+
R
k
]
−
1
K_k^*=P_{k/k-1}^*H_k^T[H_kP_{k/k-1}^*H_k^T+R_k]^{-1}
Kk∗=Pk/k−1∗HkT[HkPk/k−1∗HkT+Rk]−1
P
k
∗
=
(
I
−
K
k
∗
H
k
)
P
k
/
k
−
1
∗
P_k^*=(I-K_k^*H_k)P_{k/k-1}^*
Pk∗=(I−Kk∗Hk)Pk/k−1∗
X
^
k
∗
=
X
^
k
−
1
∗
+
K
k
∗
[
Z
k
−
H
k
Φ
k
,
k
−
1
X
^
k
−
1
∗
]
\hat X_k^*=\hat X_{k-1}^*+K_k^*[Z_k-H_k\Phi_{k,k-1}\hat X_{k-1}^*]
X^k∗=X^k−1∗+Kk∗[Zk−HkΦk,k−1X^k−1∗]
对比这套滤波方程跟卡尔曼滤波器的基本方程,其实只有第一个方程里多了个s。由于s>1,所以这样计算出来的
P
k
/
k
−
1
∗
P_{k/k-1}^*
Pk/k−1∗总比根据一般方程计算出来的
P
k
/
k
−
1
P_{k/k-1}
Pk/k−1要大,对应这
K
k
∗
K_k^*
Kk∗也比基本方程里的
K
k
K_k
Kk大。根据前面实例分析那篇文章里的分析结论,采用这套方程对新的量测值的利用权重要比基本方程里的权重大。也就是量测的权重增加而预测的权重减小,这样就达到了降低陈旧量测值对估计结果的影响了。而且越陈旧的量测,影响越小,所以这种方法叫做衰减记忆法
二、限定记忆法
前面提到的衰减记忆法思路是随着滤波次数的增加,陈旧量测的权重会逐渐减小,从而达到增加新量测的权重的效果。而限定记忆法则更直接,它的思路是:我只对当前时刻之前的N次量测感兴趣,之前的量测不要考虑。这话怎么说?
咱们不是说用卡尔曼滤波的基本方程时,
X
^
k
\hat X_k
X^k是由
Z
1
,
Z
2
,
…
,
Z
k
Z_1,Z_2,\dots,Z_k
Z1,Z2,…,Zk这些量测获得的估计值吗。而现在,我们要求的是由
Z
d
,
Z
d
+
1
,
…
,
Z
k
Z_d,Z_{d+1},\dots,Z_k
Zd,Zd+1,…,Zk这些量测来获得对X的估计,量测数量固定有N+1个(当前时刻之前的N个量测加上当前时刻的量测)。为了区分,把这个估计值定义成
X
^
k
N
\hat X_k^N
X^kN吧。
滤波方程推导有些复杂冗长,这里为了简洁,假定系统能够为0,也就是Q=0。可即使这样,得到的滤波方程也不简单:
X
^
k
N
=
Φ
k
,
k
−
1
X
^
k
−
1
N
+
K
k
[
Z
k
−
H
k
Φ
k
,
k
−
1
X
^
k
−
1
N
]
−
K
‾
k
[
Z
d
−
H
d
Φ
d
,
k
Φ
k
,
k
−
1
X
^
k
−
1
N
]
\hat X^N_k=\Phi_{k,k-1}\hat X_{k-1}^N+K_k[Z_k-H_k\Phi_{k,k-1}\hat X_{k-1}^N]-\overline K_k[Z_d-H_d\Phi_{d,k}\Phi_{k,k-1}\hat X_{k-1}^N]
X^kN=Φk,k−1X^k−1N+Kk[Zk−HkΦk,k−1X^k−1N]−Kk[Zd−HdΦd,kΦk,k−1X^k−1N]
K
k
=
P
k
N
H
k
T
R
K
−
1
K_k=P_k^NH_k^TR_K^{-1}
Kk=PkNHkTRK−1
K
‾
k
=
P
k
N
Φ
d
,
k
T
H
d
T
R
d
−
1
\overline K_k=P_k^N\Phi_{d,k}^TH_d^TR_d^{-1}
Kk=PkNΦd,kTHdTRd−1
(
P
k
N
)
−
1
=
Φ
k
−
1
,
k
T
(
P
k
−
1
N
)
−
1
Φ
k
−
1
,
k
+
H
k
T
R
k
−
1
H
k
−
Φ
d
,
k
T
H
d
T
R
d
−
1
H
d
Φ
d
,
k
(P_k^N)^{-1}=\Phi_{k-1,k}^T(P_{k-1}^N)^{-1}\Phi_{k-1,k}+H_k^TR_k^{-1}H_k-\Phi_{d,k}^TH_d^TR_d^{-1}H_d\Phi_{d,k}
(PkN)−1=Φk−1,kT(Pk−1N)−1Φk−1,k+HkTRk−1Hk−Φd,kTHdTRd−1HdΦd,k
这些式子的参数有些是需要说明的。
首先,这是在
k
>
N
k>N
k>N时才能使用的。道理也好理解,N是预计要限定记忆的长度,当k比N小的时候就直接用一般方程好了。
其次,这里有好几个
Φ
\Phi
Φ,也就是好几个状态转移矩阵,要注意状态转移矩阵的定义哦,
Φ
k
,
k
−
1
\Phi_{k,k-1}
Φk,k−1表示的是从k-1时刻到k时刻的状态转移矩阵,我们可从没有说过转移矩阵只能从前向后转移,从后向前也一样的。
所以,
Φ
k
−
1
,
k
\Phi_{k-1,k}
Φk−1,k是从k时刻到k-1时刻的转移矩阵。而
Φ
d
,
k
\Phi_{d,k}
Φd,k则表示k时刻到k之前的N个周期的d时刻的转移矩阵。这些倒着转移的矩阵其实可以认为是正着转移矩阵的逆矩阵。
还剩一个问题,就是这套滤波器的初值问题。一定要注意,这套算法并不是从k=0开始的,而是从k>N时才开始的。所以它的初值就是在N时刻的基本方程计算结果。而且为了体现限定记忆的精神,在k=N的瞬间,还要把初值的影响的记忆也要忘掉。所以一般按下面这种方法取值。
P
N
N
=
[
P
N
−
1
−
Φ
0
,
N
T
P
0
−
1
Φ
0
,
N
]
−
1
P_N^N=[P_N^{-1}-\Phi_{0,N}^TP_0^{-1}\Phi_{0,N}]^{-1}
PNN=[PN−1−Φ0,NTP0−1Φ0,N]−1
X
^
N
N
=
P
N
N
[
P
N
−
1
X
^
N
−
Φ
0
,
N
T
P
0
−
1
X
^
0
]
\hat X_N^N=P_N^N[P_N^{-1}\hat X_N-\Phi_{0,N}^TP_0^{-1}\hat X_0]
X^NN=PNN[PN−1X^N−Φ0,NTP0−1X^0]
这里的
P
N
P_N
PN和
X
^
N
\hat X_N
X^N就是在前N个周期内用基本方程滤波的结果。
小结
这样一来,我们基本把模型不够准确时的发撒抑制方法解释清楚了,当然,这也只是通常的做法。而且,更通常的情况是,这两种方法都不会轻易用,一是滤波的精度会下降,再有就是会增加计算量。尤其是限定记忆法,如果考虑状态噪声的影响,计算公式就会变得非常繁琐。所以最合理的方案还是尽量把系统数学模型建立的准确,让建立的状态方程能够尽可能的反应真实的状态才是正解。