文章目录
- 一、连续时间系统方程离散化
- 1、连续时间模型
- 2、状态转移矩阵计算
- 3、激励噪声的等效计算
- 4、最终离散化结论
- 5、常见简单随机过程离散化
- 6、实际物理信号的噪声单位
- 二、连续时间量测方程离散化
- 三、连续时间Kalman滤波
- 1、连续状态空间模型
- 2、离散时间Kalman滤波
- 3、增益矩阵的连续化
- 4、状态估计的连续化
- 5、均方差阵的连续化
- 6、连续时间Kalman滤波方程汇总
实际建模中给的可能是连续的状态方程、连续的量测方程;要用Kalman滤波得先进行离散化。
一、连续时间系统方程离散化
1、连续时间模型
X ˙ ( t ) = F ( t ) X ( t ) + G ( t ) w ( t ) E [ w ( t ) ] = 0 , E [ w ( t ) w T ( τ ) ] = q ( t ) δ ( t − τ ) \dot{\boldsymbol{X}}(t)=\boldsymbol{F}(t) \boldsymbol{X}(t)+\boldsymbol{G}(t) \boldsymbol{w}(t) \quad \mathrm{E}[\boldsymbol{w}(t)]=0, \quad \mathrm{E}\left[\boldsymbol{w}(t) \boldsymbol{w}^{\mathrm{T}}(\tau)\right]=\boldsymbol{q}(t) \delta(t-\tau) X˙(t)=F(t)X(t)+G(t)w(t)E[w(t)]=0,E[w(t)wT(τ)]=q(t)δ(t−τ)
连续时间的白噪声不太好理解,它处处连续,处处不可导,幅值无穷大,任何两个时间不相关,且是高斯白噪声
离散化形式
X
k
=
Φ
k
/
k
−
1
X
k
−
1
+
η
k
−
1
\boldsymbol{X}_{k}=\boldsymbol{\Phi}_{k / k-1} \boldsymbol{X}_{k-1}+\boldsymbol{\eta}_{k-1}
Xk=Φk/k−1Xk−1+ηk−1
其中:
X
k
=
X
(
t
k
)
Φ
k
/
k
−
1
=
Φ
(
t
k
,
t
k
−
1
)
≈
e
∫
t
k
−
1
t
k
F
(
τ
)
d
τ
η
k
−
1
=
∫
t
k
−
1
t
k
Φ
(
t
k
,
τ
)
G
(
τ
)
w
(
τ
)
d
τ
\begin{array}{l}\boldsymbol{X}_{k}=\boldsymbol{X}\left(t_{k}\right) \\ \boldsymbol{\Phi}_{k / k-1}=\boldsymbol{\Phi}\left(t_{k}, t_{k-1}\right) \approx \mathrm{e}^{\int_{t_{k-1}}^{t_{k}} F(\tau) \mathrm{d} \tau} \\ \boldsymbol{\eta}_{k-1}=\int_{t_{k-1}}^{t_{k}} \boldsymbol{\Phi}\left(t_{k}, \tau\right) \boldsymbol{G}(\tau) \boldsymbol{w}(\tau) \mathrm{d} \tau\end{array}
Xk=X(tk)Φk/k−1=Φ(tk,tk−1)≈e∫tk−1tkF(τ)dτηk−1=∫tk−1tkΦ(tk,τ)G(τ)w(τ)dτ
2、状态转移矩阵计算
泰勒展开之后取前两项;后面高阶项计算量大,收益小:
Φ
k
/
k
−
1
≈
e
F
(
t
k
−
1
)
T
s
=
I
+
F
(
t
k
−
1
)
T
s
+
F
2
(
t
k
−
1
)
I
s
2
∂
!
+
F
3
(
t
k
−
1
)
I
s
∂
!
+
⋯
≈
I
+
F
(
t
k
−
1
)
T
s
\boldsymbol{\Phi}_{k / k-1} \approx \mathrm{e}^{F\left(t_{k-1}\right) T_{s}}=\boldsymbol{I}+\boldsymbol{F}\left(t_{k-1}\right) T_{s}+\boldsymbol{F}^{2}\left(t_{k-1}\right) \frac{I_{s}^{2}}{\partial !}+\boldsymbol{F}^{3}\left(t_{k-1}\right) \frac{I_{s}}{\partial !}+\cdots \approx \boldsymbol{I}+\boldsymbol{F}\left(t_{k-1}\right) T_{s}
Φk/k−1≈eF(tk−1)Ts=I+F(tk−1)Ts+F2(tk−1)∂!Is2+F3(tk−1)∂!Is+⋯≈I+F(tk−1)Ts
3、激励噪声的等效计算
η
k
−
1
=
∫
t
k
−
1
t
k
Φ
(
t
k
,
τ
)
G
(
τ
)
w
(
τ
)
d
τ
\boldsymbol{\eta}_{k-1}=\int_{t_{k-1}}^{t_{k}} \boldsymbol{\Phi}\left(t_{k}, \tau\right) \boldsymbol{G}(\tau) \boldsymbol{w}(\tau) \mathrm{d} \tau
ηk−1=∫tk−1tkΦ(tk,τ)G(τ)w(τ)dτ
连续时间系统的噪声是高斯白噪声,经过线性变换(微积分也是线性变换),得到的还是高斯白噪声。想了解它,只要知道它的一阶矩和二阶矩就行了
激励噪声均值:
E
[
η
k
−
1
]
=
E
[
∫
t
k
−
1
t
k
Φ
(
t
k
,
τ
)
G
(
τ
)
w
(
τ
)
d
τ
]
=
∫
t
k
−
1
t
k
Φ
(
t
k
,
τ
)
G
(
τ
)
E
[
w
(
τ
)
]
d
τ
=
0
\mathrm{E}\left[\boldsymbol{\eta}_{k-1}\right]=\mathrm{E}\left[\int_{t_{k-1}}^{t_{k}} \boldsymbol{\Phi}\left(t_{k}, \tau\right) \boldsymbol{G}(\tau) \boldsymbol{w}(\tau) \mathrm{d} \tau\right]=\int_{t_{k-1}}^{t_{k}} \boldsymbol{\Phi}\left(t_{k}, \tau\right) \boldsymbol{G}(\tau) \mathrm{E}[\boldsymbol{w}(\tau)] \mathrm{d} \tau=\mathbf{0}
E[ηk−1]=E[∫tk−1tkΦ(tk,τ)G(τ)w(τ)dτ]=∫tk−1tkΦ(tk,τ)G(τ)E[w(τ)]dτ=0
激励噪声协方差:
E
[
η
k
−
1
η
j
−
1
T
]
=
0
k
≠
j
\mathrm{E}\left[\boldsymbol{\eta}_{k-1} \boldsymbol{\eta}_{j-1}^{\mathrm{T}}\right]=\mathbf{0} \quad k \neq j
E[ηk−1ηj−1T]=0k=j
E [ η k − 1 η k − 1 T ] = E { ∫ t k − 1 t k Φ ( t k , τ ) G ( τ ) w ( τ ) d τ ⋅ [ ∫ t k − 1 t k Φ ( t k , s ) G ( s ) w ( s ) d s ] T } k = j = E [ ∫ t k − 1 t k Φ ( t k , τ ) G ( τ ) w ( τ ) ∫ t k − 1 t k w T ( s ) G T ( s ) Φ T ( t k , s ) d s d τ ] = ∫ t k − 1 t k Φ ( t k , τ ) G ( τ ) ∫ t k − 1 t k E [ w ( τ ) w T ( s ) ] G T ( s ) Φ T ( t k , s ) d s d τ = ∫ t k − 1 t k Φ ( t k , τ ) G ( τ ) ∫ t k − 1 t k q ( τ ) δ ( τ − s ) G T ( s ) Φ T ( t k , s ) d s d τ = ∫ t k − 1 t k Φ ( t k , τ ) G ( τ ) q ( τ ) G T ( τ ) Φ T ( t k , τ ) d τ \begin{aligned} \mathrm{E}\left[\boldsymbol{\eta}_{k-1} \boldsymbol{\eta}_{k-1}^{\mathrm{T}}\right] & =\mathrm{E}\left\{\int_{t_{k-1}}^{t_{k}} \boldsymbol{\Phi}\left(t_{k}, \tau\right) \boldsymbol{G}(\tau) \boldsymbol{w}(\tau) \mathrm{d} \tau \cdot\left[\int_{t_{k-1}}^{t_{k}} \boldsymbol{\Phi}\left(t_{k}, s\right) \boldsymbol{G}(s) \boldsymbol{w}(s) \mathrm{d} s\right]^{\mathrm{T}}\right\} \quad k=j \\ & =\mathrm{E}\left[\int_{t_{k-1}}^{t_{k}} \boldsymbol{\Phi}\left(t_{k}, \tau\right) \boldsymbol{G}(\tau) \boldsymbol{w}(\tau) \int_{t_{k-1}}^{t_{k}} \boldsymbol{w}^{\mathrm{T}}(s) \boldsymbol{G}^{\mathrm{T}}(s) \boldsymbol{\Phi}^{\mathrm{T}}\left(t_{k}, s\right) \mathrm{d} s \mathrm{~d} \tau\right] \\ & =\int_{t_{k-1}}^{t_{k}} \boldsymbol{\Phi}\left(t_{k}, \tau\right) \boldsymbol{G}(\tau) \int_{t_{k-1}}^{t_{k}} \mathrm{E}\left[\boldsymbol{w}(\tau) \boldsymbol{w}^{\mathrm{T}}(s)\right] \boldsymbol{G}^{\mathrm{T}}(s) \boldsymbol{\Phi}^{\mathrm{T}}\left(t_{k}, s\right) \mathrm{d} s \mathrm{~d} \tau \\ & =\int_{t_{k-1}}^{t_{k}} \boldsymbol{\Phi}\left(t_{k}, \tau\right) \boldsymbol{G}(\tau) \int_{t_{k-1}}^{t_{k}} \boldsymbol{q}(\tau) \delta(\tau-s) \boldsymbol{G}^{\mathrm{T}}(s) \boldsymbol{\Phi}^{\mathrm{T}}\left(t_{k}, s\right) \mathrm{d} s \mathrm{~d} \tau \\ & =\int_{t_{k-1}}^{t_{k}} \boldsymbol{\Phi}\left(t_{k}, \tau\right) \boldsymbol{G}(\tau) \boldsymbol{q}(\tau) \boldsymbol{G}^{\mathrm{T}}(\tau) \boldsymbol{\Phi}^{\mathrm{T}}\left(t_{k}, \tau\right) \mathrm{d} \tau\end{aligned} E[ηk−1ηk−1T]=E{∫tk−1tkΦ(tk,τ)G(τ)w(τ)dτ⋅[∫tk−1tkΦ(tk,s)G(s)w(s)ds]T}k=j=E[∫tk−1tkΦ(tk,τ)G(τ)w(τ)∫tk−1tkwT(s)GT(s)ΦT(tk,s)ds dτ]=∫tk−1tkΦ(tk,τ)G(τ)∫tk−1tkE[w(τ)wT(s)]GT(s)ΦT(tk,s)ds dτ=∫tk−1tkΦ(tk,τ)G(τ)∫tk−1tkq(τ)δ(τ−s)GT(s)ΦT(tk,s)ds dτ=∫tk−1tkΦ(tk,τ)G(τ)q(τ)GT(τ)ΦT(tk,τ)dτ
想往下算,得做一些假设,在
t
k
−
1
t_{k-1}
tk−1 到
t
k
t_k
tk 这一小段时间内,
Φ
(
t
k
,
τ
)
G
(
τ
)
\boldsymbol{\Phi}\left(t_{k}, \tau\right) \boldsymbol{G}(\tau)
Φ(tk,τ)G(τ) 都为常值,时不变;不是常值可以取一阶项
E
[
η
k
−
1
η
k
−
1
T
]
=
∫
t
k
−
1
t
k
Φ
(
t
k
,
τ
)
G
(
τ
)
q
(
τ
)
G
T
(
τ
)
Φ
T
(
t
k
,
τ
)
d
τ
≈
∫
t
k
−
1
t
k
[
I
+
F
(
t
k
−
1
)
(
t
k
−
τ
)
]
G
(
t
k
−
1
)
q
(
t
k
−
1
)
G
T
(
t
k
−
1
)
[
I
+
F
(
t
k
−
1
)
(
t
k
−
τ
)
]
T
d
τ
=
∫
t
k
−
1
t
k
G
(
t
k
−
1
)
q
G
T
(
t
k
−
1
)
d
τ
+
∫
t
k
−
1
t
k
G
(
t
k
−
1
)
q
(
t
k
−
1
)
G
T
(
t
k
−
1
)
F
T
(
t
k
−
1
)
(
t
k
−
τ
)
d
τ
+
∫
t
k
−
1
t
k
F
(
t
k
−
1
)
G
(
t
k
−
1
)
q
(
t
k
−
1
)
G
T
(
t
k
−
1
)
(
t
k
−
τ
)
d
τ
+
∫
t
k
−
1
t
k
F
(
t
k
−
1
)
G
(
t
k
−
1
)
q
(
t
k
−
1
)
G
T
(
t
k
−
1
)
F
T
(
t
k
−
1
)
(
t
k
−
τ
)
2
d
τ
=
G
(
t
k
−
1
)
q
(
t
k
−
1
)
G
T
(
t
k
−
1
)
T
s
+
1
2
G
(
t
k
−
1
)
q
(
t
k
−
1
)
G
T
(
t
k
−
1
)
F
T
(
t
k
−
1
)
T
s
2
+
1
2
F
(
t
k
−
1
)
G
(
t
k
−
1
)
q
(
t
k
−
1
)
G
T
(
t
k
−
1
)
T
s
2
+
1
3
F
(
t
k
−
1
)
G
(
t
k
−
1
)
q
(
t
k
−
1
)
G
T
(
t
k
−
1
)
F
T
(
t
k
−
1
)
T
s
3
=
[
I
+
1
2
F
(
t
k
−
1
)
T
s
]
⋅
[
G
(
t
k
−
1
)
q
(
t
k
−
1
)
G
T
(
t
k
−
1
)
T
T
s
]
⋅
[
I
+
1
2
F
(
t
k
−
1
)
T
s
]
T
+
1
12
F
(
t
k
−
1
)
G
(
t
k
−
1
)
q
(
t
k
−
1
)
G
T
(
t
k
−
1
)
F
T
(
t
k
−
1
)
T
s
3
≈
{
[
I
+
1
2
F
(
t
k
−
1
)
T
s
]
G
(
t
k
−
1
)
}
⋅
[
q
(
t
k
−
1
)
T
s
]
⋅
{
[
I
+
1
2
F
(
t
k
−
1
)
T
s
]
G
(
t
k
−
1
)
}
T
\begin{aligned} \mathrm{E} & {\left[\boldsymbol{\eta}_{k-1} \boldsymbol{\eta}_{k-1}^{\mathrm{T}}\right]=\int_{t_{k-1}}^{t_{k}} \boldsymbol{\Phi}\left(t_{k}, \tau\right) \boldsymbol{G}(\tau) \boldsymbol{q}(\tau) \boldsymbol{G}^{\mathrm{T}}(\tau) \boldsymbol{\Phi}^{\mathrm{T}}\left(t_{k}, \tau\right) \mathrm{d} \tau } \\ \approx & \int_{t_{k-1}}^{t_{k}}\left[\boldsymbol{I}+\boldsymbol{F}\left(t_{k-1}\right)\left(t_{k}-\tau\right)\right] \boldsymbol{G}\left(t_{k-1}\right) \boldsymbol{q}\left(t_{k-1}\right) \boldsymbol{G}^{\mathrm{T}}\left(t_{k-1}\right)\left[\boldsymbol{I}+\boldsymbol{F}\left(t_{k-1}\right)\left(t_{k}-\tau\right)\right]^{\mathrm{T}} \mathrm{d} \tau \\ = & \int_{t_{k-1}}^{t_{k}} \boldsymbol{G}\left(t_{k-1}\right) \boldsymbol{q} \boldsymbol{G}^{\mathrm{T}}\left(t_{k-1}\right) \mathrm{d} \tau+\int_{t_{k-1}}^{t_{k}} \boldsymbol{G}\left(t_{k-1}\right) \boldsymbol{q}\left(t_{k-1}\right) \boldsymbol{G}^{\mathrm{T}}\left(t_{k-1}\right) \boldsymbol{F}^{\mathrm{T}}\left(t_{k-1}\right)\left(t_{k}-\tau\right) \mathrm{d} \tau \\ & \quad+\int_{t_{k-1}}^{t_{k}} \boldsymbol{F}\left(t_{k-1}\right) \boldsymbol{G}\left(t_{k-1}\right) \boldsymbol{q}\left(t_{k-1}\right) \boldsymbol{G}^{\mathrm{T}}\left(t_{k-1}\right)\left(t_{k}-\tau\right) \mathrm{d} \tau+\int_{t_{k-1}}^{t_{k}} \boldsymbol{F}\left(t_{k-1}\right) \boldsymbol{G}\left(t_{k-1}\right) \boldsymbol{q}\left(t_{k-1}\right) \boldsymbol{G}^{\mathrm{T}}\left(t_{k-1}\right) \boldsymbol{F}^{\mathrm{T}}\left(t_{k-1}\right)\left(t_{k}-\tau\right)^{2} \mathrm{~d} \tau \\ = & \boldsymbol{G}\left(t_{k-1}\right) \boldsymbol{q}\left(t_{k-1}\right) \boldsymbol{G}^{\mathrm{T}}\left(t_{k-1}\right) T_{s}+\frac{1}{2} \boldsymbol{G}\left(t_{k-1}\right) \boldsymbol{q}\left(t_{k-1}\right) \boldsymbol{G}^{\mathrm{T}}\left(t_{k-1}\right) \boldsymbol{F}^{\mathrm{T}}\left(t_{k-1}\right) T_{s}^{2} \\ & \quad+\frac{1}{2} \boldsymbol{F}\left(t_{k-1}\right) \boldsymbol{G}\left(t_{k-1}\right) \boldsymbol{q}\left(t_{k-1}\right) \boldsymbol{G}^{\mathrm{T}}\left(t_{k-1}\right) T_{s}^{2}+\frac{1}{3} \boldsymbol{F}\left(t_{k-1}\right) \boldsymbol{G}\left(t_{k-1}\right) \boldsymbol{q}\left(t_{k-1}\right) \boldsymbol{G}^{\mathrm{T}}\left(t_{k-1}\right) \boldsymbol{F}^{\mathrm{T}}\left(t_{k-1}\right) T_{s}^{3} \\ = & {\left[\boldsymbol{I}+\frac{1}{2} \boldsymbol{F}\left(t_{k-1}\right) T_{s}\right] \cdot\left[\boldsymbol{G}\left(t_{k-1}\right) \boldsymbol{q}\left(t_{k-1}\right) \boldsymbol{G}^{\mathrm{T}}\left(t_{k-1}\right) T T_{s}\right] \cdot\left[\boldsymbol{I}+\frac{1}{2} \boldsymbol{F}\left(t_{k-1}\right) T_{s}\right]^{\mathrm{T}}+\frac{1}{12} \boldsymbol{F}\left(t_{k-1}\right) \boldsymbol{G}\left(t_{k-1}\right) \boldsymbol{q}\left(t_{k-1}\right) \boldsymbol{G}^{\mathrm{T}}\left(t_{k-1}\right) \boldsymbol{F}^{\mathrm{T}}\left(t_{k-1}\right) T_{s}^{3} } \\ \approx & \left\{\left[\boldsymbol{I}+\frac{1}{2} \boldsymbol{F}\left(t_{k-1}\right) T_{s}\right] \boldsymbol{G}\left(t_{k-1}\right)\right\} \cdot\left[\boldsymbol{q}\left(t_{k-1}\right) T_{s}\right] \cdot\left\{\left[\boldsymbol{I}+\frac{1}{2} \boldsymbol{F}\left(t_{k-1}\right) T_{s}\right] \boldsymbol{G}\left(t_{k-1}\right)\right\}^{\mathrm{T}}\end{aligned}
E≈===≈[ηk−1ηk−1T]=∫tk−1tkΦ(tk,τ)G(τ)q(τ)GT(τ)ΦT(tk,τ)dτ∫tk−1tk[I+F(tk−1)(tk−τ)]G(tk−1)q(tk−1)GT(tk−1)[I+F(tk−1)(tk−τ)]Tdτ∫tk−1tkG(tk−1)qGT(tk−1)dτ+∫tk−1tkG(tk−1)q(tk−1)GT(tk−1)FT(tk−1)(tk−τ)dτ+∫tk−1tkF(tk−1)G(tk−1)q(tk−1)GT(tk−1)(tk−τ)dτ+∫tk−1tkF(tk−1)G(tk−1)q(tk−1)GT(tk−1)FT(tk−1)(tk−τ)2 dτG(tk−1)q(tk−1)GT(tk−1)Ts+21G(tk−1)q(tk−1)GT(tk−1)FT(tk−1)Ts2+21F(tk−1)G(tk−1)q(tk−1)GT(tk−1)Ts2+31F(tk−1)G(tk−1)q(tk−1)GT(tk−1)FT(tk−1)Ts3[I+21F(tk−1)Ts]⋅[G(tk−1)q(tk−1)GT(tk−1)TTs]⋅[I+21F(tk−1)Ts]T+121F(tk−1)G(tk−1)q(tk−1)GT(tk−1)FT(tk−1)Ts3{[I+21F(tk−1)Ts]G(tk−1)}⋅[q(tk−1)Ts]⋅{[I+21F(tk−1)Ts]G(tk−1)}T
如果时间极短,还可以再简化:
E
[
η
k
−
1
η
k
−
1
T
]
≈
G
(
t
k
−
1
)
⋅
[
q
(
t
k
−
1
)
T
s
]
⋅
G
T
(
t
k
−
1
)
F
(
t
k
−
1
)
T
s
≪
I
\mathrm{E}\left[\boldsymbol{\eta}_{k-1} \boldsymbol{\eta}_{k-1}^{\mathrm{T}}\right] \approx \boldsymbol{G}\left(t_{k-1}\right) \cdot\left[\boldsymbol{q}\left(t_{k-1}\right) T_{s}\right] \cdot \boldsymbol{G}^{\mathrm{T}}\left(t_{k-1}\right) \quad \boldsymbol{F}\left(t_{k-1}\right) T_{s} \ll \boldsymbol{I}
E[ηk−1ηk−1T]≈G(tk−1)⋅[q(tk−1)Ts]⋅GT(tk−1)F(tk−1)Ts≪I
4、最终离散化结论
或等价表示为:
两者都满足:
E
[
(
Γ
k
−
1
W
k
−
1
)
(
Γ
k
−
1
W
k
−
1
)
T
]
=
G
(
t
k
−
1
)
⋅
[
q
(
t
k
−
1
)
T
s
]
⋅
G
T
(
t
k
−
1
)
=
E
[
(
Γ
k
−
1
′
W
k
−
1
′
)
(
Γ
k
−
1
′
W
k
−
1
′
)
T
]
\mathrm{E}\left[\left(\boldsymbol{\Gamma}_{k-1} \boldsymbol{W}_{k-1}\right)\left(\boldsymbol{\Gamma}_{k-1} \boldsymbol{W}_{k-1}\right)^{\mathrm{T}}\right]=\boldsymbol{G}\left(t_{k-1}\right) \cdot\left[\boldsymbol{q}\left(t_{k-1}\right) T_{s}\right] \cdot \boldsymbol{G}^{\mathrm{T}}\left(t_{k-1}\right)=\mathrm{E}\left[\left(\boldsymbol{\Gamma}_{k-1}^{\prime} \boldsymbol{W}_{k-1}^{\prime}\right)\left(\boldsymbol{\Gamma}_{k-1}^{\prime} \boldsymbol{W}_{k-1}^{\prime}\right)^{\mathrm{T}}\right]
E[(Γk−1Wk−1)(Γk−1Wk−1)T]=G(tk−1)⋅[q(tk−1)Ts]⋅GT(tk−1)=E[(Γk−1′Wk−1′)(Γk−1′Wk−1′)T]
5、常见简单随机过程离散化
1.一阶马尔可夫: X ˙ ( t ) = − β X ( t ) + w ( t ) X k = a 1 X k − 1 + W k − 1 \dot{X}(t)=-\beta X(t)+w(t) \quad X_{k}=a_{1} X_{k-1}+W_{k-1} X˙(t)=−βX(t)+w(t)Xk=a1Xk−1+Wk−1
2.二阶马尔可夫: X ¨ ( t ) = − 2 β X ˙ ( t ) − β 2 X ( t ) + w ( t ) X k = a 1 X k − 1 + a 2 X k − 2 + W k − 1 \ddot{X}(t)=-2 \beta \dot{X}(t)-\beta^{2} X(t)+w(t) \quad X_{k}=a_{1} X_{k-1}+a_{2} X_{k-2}+W_{k-1} X¨(t)=−2βX˙(t)−β2X(t)+w(t)Xk=a1Xk−1+a2Xk−2+Wk−1
3.随机游走: X ˙ ( t ) = w ( t ) X k = X k − 1 + W k − 1 \dot{X}(t)=w(t) \quad X_{k}=X_{k-1}+W_{k-1} X˙(t)=w(t)Xk=Xk−1+Wk−1
4.随机常值: X ˙ ( t ) = 0 X k = X 0 \dot{X}(t)=0 \quad X_{k}=X_{0} X˙(t)=0Xk=X0
6、实际物理信号的噪声单位
- τ \tau τ 表示相关时间,很多时候等于 1 1 1,所以 − 1 / τ -1 / \tau −1/τ 可能被省略
- q − U / s / H z \sqrt{q}-\mathrm{U} / \mathrm{s} / \sqrt{Hz } \mathrm{} q−U/s/Hz 激励噪声密度,描述激励噪声本身。 U / S \mathrm{U} / \sqrt{ S} U/S 随机游走噪声系数,是想描述 X X X 的噪声。描述的主体不同。比如惯导中有:
二、连续时间量测方程离散化
连续量测方程的只是理论上的抽象,没有实际意义
连续量测模型:
Z
(
t
)
=
H
(
t
)
X
(
t
)
+
v
(
t
)
E
[
v
(
t
)
]
=
0
,
E
[
v
(
t
)
v
T
(
τ
)
]
=
r
(
t
)
δ
(
t
−
τ
)
\boldsymbol{Z}(t)=\boldsymbol{H}(t) \boldsymbol{X}(t)+\boldsymbol{v}(t) \quad \mathrm{E}[\boldsymbol{v}(t)]=\mathbf{0}, \quad \mathrm{E}\left[\boldsymbol{v}(t) \boldsymbol{v}^{\mathrm{T}}(\tau)\right]=\boldsymbol{r}(t) \delta(t-\tau)
Z(t)=H(t)X(t)+v(t)E[v(t)]=0,E[v(t)vT(τ)]=r(t)δ(t−τ)
在 $ t_{k-1} $ 到
t
k
t_{k}
tk 时间段,两边积分平均:
1
T
s
∫
t
k
−
1
t
k
Z
(
τ
)
d
τ
=
1
T
v
∫
t
k
−
1
t
k
H
(
τ
)
X
(
τ
)
+
v
(
τ
)
d
τ
=
1
T
v
∫
t
k
−
1
t
k
H
(
τ
)
X
(
τ
)
d
τ
+
1
T
s
∫
t
k
−
1
t
k
v
(
τ
)
d
τ
\frac{1}{T_{s}} \int_{t_{k-1}}^{t_{k}} \boldsymbol{Z}(\tau) \mathrm{d} \tau=\frac{1}{T_{v}} \int_{t_{k-1}}^{t_{k}} \boldsymbol{H}(\tau) \boldsymbol{X}(\tau)+\boldsymbol{v}(\tau) \mathrm{d} \tau=\frac{1}{T_{v}} \int_{t_{k-1}}^{t_{k}} \boldsymbol{H}(\tau) \boldsymbol{X}(\tau) \mathrm{d} \tau+\frac{1}{T_{s}} \int_{t_{k-1}}^{t_{k}} \boldsymbol{v}(\tau) \mathrm{d} \tau
Ts1∫tk−1tkZ(τ)dτ=Tv1∫tk−1tkH(τ)X(τ)+v(τ)dτ=Tv1∫tk−1tkH(τ)X(τ)dτ+Ts1∫tk−1tkv(τ)dτ
离散化记为:
Z
k
≈
H
k
X
k
+
V
k
\boldsymbol{Z}_{k} \approx \boldsymbol{H}_{k} \boldsymbol{X}_{k}+\boldsymbol{V}_{k}
Zk≈HkXk+Vk
观测噪声均值:
E
[
V
k
]
=
E
[
1
T
s
∫
t
k
−
1
t
k
v
(
τ
)
d
τ
]
=
1
T
s
∫
t
k
−
1
t
k
E
[
v
(
τ
)
]
d
τ
=
0
\mathrm{E}\left[\boldsymbol{V}_{k}\right]=\mathrm{E}\left[\frac{1}{T_{s}} \int_{t_{k-1}}^{t_{k}} \boldsymbol{v}(\tau) \mathrm{d} \tau\right]=\frac{1}{T_{s}} \int_{t_{k-1}}^{t_{k}} \mathrm{E}[\boldsymbol{v}(\tau)] \mathrm{d} \tau=\mathbf{0}
E[Vk]=E[Ts1∫tk−1tkv(τ)dτ]=Ts1∫tk−1tkE[v(τ)]dτ=0
观测噪声协方差:
E
[
V
k
V
j
T
]
=
E
[
(
1
T
s
∫
t
k
−
1
t
k
v
(
τ
)
d
τ
)
(
1
T
s
∫
t
j
−
1
t
j
v
(
s
)
d
s
)
T
]
=
1
T
s
2
∫
t
k
−
1
t
k
∫
t
j
−
1
t
j
E
[
v
(
τ
)
v
T
(
s
)
]
d
s
d
τ
=
1
T
s
2
∫
t
k
−
1
t
k
∫
t
j
−
1
t
j
r
(
τ
)
δ
(
s
−
τ
)
d
s
d
=
1
T
s
2
∫
t
k
−
1
t
k
r
(
τ
)
δ
k
j
d
τ
≈
r
(
t
k
)
T
s
δ
k
j
≜
R
k
δ
k
j
\begin{aligned} \mathrm{E}\left[\boldsymbol{V}_{k} \boldsymbol{V}_{j}^{\mathrm{T}}\right] & =\mathrm{E}\left[\left(\frac{1}{T_{s}} \int_{t_{k-1}}^{t_{k}} \boldsymbol{v}(\tau) \mathrm{d} \tau\right)\left(\frac{1}{T_{s}} \int_{t_{j-1}}^{t_{j}} \boldsymbol{v}(s) \mathrm{d} s\right)^{\mathrm{T}}\right] \\ & =\frac{1}{T_{s}^{2}} \int_{t_{k-1}}^{t_{k}} \int_{t_{j-1}}^{t_{j}} \mathrm{E}\left[\boldsymbol{v}(\tau) \boldsymbol{v}^{\mathrm{T}}(s)\right] \mathrm{d} s \mathrm{~d} \tau=\frac{1}{T_{s}^{2}} \int_{t_{k-1}}^{t_{k}} \int_{t_{j-1}}^{t_{j}} \boldsymbol{r}(\tau) \delta(s-\tau) \mathrm{d} s \mathrm{~d} \\ & =\frac{1}{T_{s}^{2}} \int_{t_{k-1}}^{t_{k}} \boldsymbol{r}(\tau) \delta_{k j} \mathrm{~d} \tau \approx \frac{\boldsymbol{r}\left(t_{k}\right)}{T_{s}} \delta_{k j} \triangleq \boldsymbol{R}_{k} \delta_{k j}\end{aligned}
E[VkVjT]=E
(Ts1∫tk−1tkv(τ)dτ)(Ts1∫tj−1tjv(s)ds)T
=Ts21∫tk−1tk∫tj−1tjE[v(τ)vT(s)]ds dτ=Ts21∫tk−1tk∫tj−1tjr(τ)δ(s−τ)ds d=Ts21∫tk−1tkr(τ)δkj dτ≈Tsr(tk)δkj≜Rkδkj
离散化量测间隔越小,等效出来的量测方差 $ \boldsymbol{R}_{k}$ 越大
三、连续时间Kalman滤波
连续的函数可以认为是离散函数,但离散间隔趋于 0 0 0 ;推导不是太严格,但比较直观。
1、连续状态空间模型
函数模型:
{
X
˙
(
t
)
=
F
(
t
)
X
(
t
)
+
G
(
t
)
w
(
t
)
Z
(
t
)
=
H
(
t
)
X
(
t
)
+
v
(
t
)
\left\{\begin{array}{l}\dot{\boldsymbol{X}}(t)=\boldsymbol{F}(t) \boldsymbol{X}(t)+\boldsymbol{G}(t) \boldsymbol{w}(t) \\ \boldsymbol{Z}(t)=\boldsymbol{H}(t) \boldsymbol{X}(t)+\boldsymbol{v}(t)\end{array}\right.
{X˙(t)=F(t)X(t)+G(t)w(t)Z(t)=H(t)X(t)+v(t)
随机模型:
{
E
[
w
(
t
)
]
=
0
,
E
[
w
(
t
)
w
T
(
τ
)
]
=
q
(
t
)
δ
(
t
−
τ
)
E
[
v
(
t
)
]
=
0
,
E
[
v
(
t
)
v
T
(
τ
)
]
=
r
(
t
)
δ
(
t
−
τ
)
E
[
w
(
t
)
v
T
(
τ
)
]
=
0
\left\{\begin{array}{l}\mathrm{E}[\boldsymbol{w}(t)]=\mathbf{0}, \quad \mathrm{E}\left[\boldsymbol{w}(t) \boldsymbol{w}^{\mathrm{T}}(\tau)\right]=\boldsymbol{q}(t) \delta(t-\tau) \\ \mathrm{E}[\boldsymbol{v}(t)]=\mathbf{0}, \quad \mathrm{E}\left[\boldsymbol{v}(t) \boldsymbol{v}^{\mathrm{T}}(\tau)\right]=\boldsymbol{r}(t) \delta(t-\tau) \\ \mathrm{E}\left[\boldsymbol{w}(t) \boldsymbol{v}^{\mathrm{T}}(\tau)\right]=\mathbf{0}\end{array}\right.
⎩
⎨
⎧E[w(t)]=0,E[w(t)wT(τ)]=q(t)δ(t−τ)E[v(t)]=0,E[v(t)vT(τ)]=r(t)δ(t−τ)E[w(t)vT(τ)]=0
2、离散时间Kalman滤波
- 状态一步预测: X ^ k ∣ k − 1 = Φ k ∣ k − 1 X ^ k − 1 \hat{\boldsymbol{X}}_{k \mid k-1}=\boldsymbol{\Phi}_{k \mid k-1} \hat{\boldsymbol{X}}_{k-1} X^k∣k−1=Φk∣k−1X^k−1
- 状态一步预测均方误差: P k / k − 1 = Φ k / k − 1 P k − 1 Φ k / k − 1 T + Γ k − 1 Q k − 1 Γ k − 1 T \boldsymbol{P}_{k / k-1}=\boldsymbol{\Phi}_{k / k-1} \boldsymbol{P}_{k-1} \boldsymbol{\Phi}_{k / k-1}^{\mathrm{T}}+\boldsymbol{\Gamma}_{k-1} \boldsymbol{Q}_{k-1} \boldsymbol{\Gamma}_{k-1}^{\mathrm{T}} Pk/k−1=Φk/k−1Pk−1Φk/k−1T+Γk−1Qk−1Γk−1T
- 滤波增益: K k = P k / k − 1 H k T ( H k P k / k − 1 H k T + R k ) − 1 \boldsymbol{K}_{k}=\boldsymbol{P}_{k / k-1} \boldsymbol{H}_{k}^{\mathrm{T}}\left(\boldsymbol{H}_{k} \boldsymbol{P}_{k / k-1} \boldsymbol{H}_{k}^{\mathrm{T}}+\boldsymbol{R}_{k}\right)^{-1} Kk=Pk/k−1HkT(HkPk/k−1HkT+Rk)−1
- 状态估计: X ^ k = X ^ k / k − 1 + K k ( Z k − H k X ^ k / k − 1 ) \hat{\boldsymbol{X}}_{k}=\hat{\boldsymbol{X}}_{k / k-1}+\boldsymbol{K}_{k}\left(\boldsymbol{Z}_{k}-\boldsymbol{H}_{k} \hat{\boldsymbol{X}}_{k / k-1}\right) X^k=X^k/k−1+Kk(Zk−HkX^k/k−1)
- 状态估计均方误差: P k = ( I − K k H k ) P k / k − 1 \boldsymbol{P}_{k}=\left(\boldsymbol{I}-\boldsymbol{K}_{k} \boldsymbol{H}_{k}\right) \boldsymbol{P}_{k / k-1} Pk=(I−KkHk)Pk/k−1
$ t_{k-1} $ 到 t k t_{k} tk 时间段趋于 0 0 0 就是连续。
3、增益矩阵的连续化
将离散方程里的量测噪声
R
k
\boldsymbol{R}_{k}
Rk 换成连续方程里的
r
(
t
k
)
T
s
\frac{\boldsymbol{r}\left(t_{k}\right)}{T_{s}}
Tsr(tk) 得:
K
k
=
P
k
H
k
T
R
k
−
1
=
P
k
H
k
T
[
r
(
t
k
)
T
s
]
−
1
=
T
s
P
k
H
k
T
r
−
1
(
t
k
)
\boldsymbol{K}_{k}=\boldsymbol{P}_{k} \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{R}_{k}^{-1}=\boldsymbol{P}_{k} \boldsymbol{H}_{k}^{\mathrm{T}}\left[\frac{\boldsymbol{r}\left(t_{k}\right)}{T_{s}}\right]^{-1}=T_{s} \boldsymbol{P}_{k} \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{r}^{-1}\left(t_{k}\right)
Kk=PkHkTRk−1=PkHkT[Tsr(tk)]−1=TsPkHkTr−1(tk)
将
T
s
T_s
Ts 除到左边,求
T
s
T_s
Ts 趋于
0
0
0 的极限:
K
(
t
)
=
lim
T
s
→
0
K
k
T
s
=
lim
T
s
→
0
P
k
H
k
T
r
−
1
(
t
k
)
=
P
(
t
)
H
T
(
t
)
r
−
1
(
t
)
\boldsymbol{K}(t)=\lim \limits_{T_{s} \rightarrow 0} \frac{\boldsymbol{K}_{k}}{T_{s}}=\lim \limits_{T_{s} \rightarrow 0} \boldsymbol{P}_{k} \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{r}^{-1}\left(t_{k}\right)=\boldsymbol{P}(t) \boldsymbol{H}^{\mathrm{T}}(t) \boldsymbol{r}^{-1}(t)
K(t)=Ts→0limTsKk=Ts→0limPkHkTr−1(tk)=P(t)HT(t)r−1(t)
即离散方程的增益
K
k
{K}_{k}
Kk 除以时间间隔
T
s
T_s
Ts 记成连续方程的增益;如果
T
s
T_s
Ts 趋于
0
0
0 ,
K
k
{K}_{k}
Kk 也趋于 0,时间很短,预测很小,量测的修正自然也很小。
4、状态估计的连续化
将一阶的状态转移矩阵带入:
X
^
k
=
Φ
k
/
k
−
1
X
^
k
−
1
+
K
k
(
Z
k
−
H
k
Φ
k
/
k
−
1
X
^
k
−
1
)
=
[
I
+
F
(
t
k
−
1
)
T
s
]
X
^
k
−
1
+
K
k
{
Z
k
−
H
k
[
I
+
F
(
t
k
−
1
)
T
s
]
X
^
k
−
1
}
=
X
^
k
−
1
+
F
(
t
k
−
1
)
X
^
k
−
1
T
s
+
K
k
[
Z
k
−
H
k
X
^
k
−
1
−
H
k
F
(
t
k
−
1
)
X
^
k
−
1
T
s
]
\begin{aligned} \hat{\boldsymbol{X}}_{k} & =\boldsymbol{\Phi}_{k / k-1} \hat{\boldsymbol{X}}_{k-1}+\boldsymbol{K}_{k}\left(\boldsymbol{Z}_{k}-\boldsymbol{H}_{k} \boldsymbol{\Phi}_{k / k-1} \hat{\boldsymbol{X}}_{k-1}\right) \\ & =\left[\boldsymbol{I}+\boldsymbol{F}\left(t_{k-1}\right) T_{s}\right] \hat{\boldsymbol{X}}_{k-1}+\boldsymbol{K}_{k}\left\{\boldsymbol{Z}_{k}-\boldsymbol{H}_{k}\left[\boldsymbol{I}+\boldsymbol{F}\left(t_{k-1}\right) T_{s}\right] \hat{\boldsymbol{X}}_{k-1}\right\} \\ & ={\color{red}\hat{\boldsymbol{X}}_{k-1}}+\boldsymbol{F}\left(t_{k-1}\right) \hat{\boldsymbol{X}}_{k-1} T_{s}+\boldsymbol{K}_{k}\left[\boldsymbol{Z}_{k}-\boldsymbol{H}_{k} \hat{\boldsymbol{X}}_{k-1}-\boldsymbol{H}_{k} \boldsymbol{F}\left(t_{k-1}\right) \hat{\boldsymbol{X}}_{k-1} T_{s}\right]\end{aligned}
X^k=Φk/k−1X^k−1+Kk(Zk−HkΦk/k−1X^k−1)=[I+F(tk−1)Ts]X^k−1+Kk{Zk−Hk[I+F(tk−1)Ts]X^k−1}=X^k−1+F(tk−1)X^k−1Ts+Kk[Zk−HkX^k−1−HkF(tk−1)X^k−1Ts]
将
X
^
k
−
1
{\color{red}\hat{\boldsymbol{X}}_{k-1}}
X^k−1 减到左边,
k
k
k 时刻的值,减去
k
−
1
k-1
k−1 时刻的值,除以
T
s
T_s
Ts,求
T
s
T_s
Ts 趋于
0
0
0 的极限,就是连续的方程:
X
^
˙
(
t
)
=
lim
T
s
→
0
X
k
−
X
k
−
1
T
s
=
lim
T
s
→
0
F
(
t
k
−
1
)
X
^
k
−
1
+
K
k
T
s
[
Z
k
−
H
k
X
^
k
−
1
−
H
k
F
(
t
k
−
1
)
X
^
k
−
1
T
s
]
=
F
(
t
)
X
^
(
t
)
+
K
(
t
)
[
Z
(
t
)
−
H
(
t
)
X
^
(
t
)
]
\begin{aligned} \dot{\hat{\boldsymbol{X}}}(t) & =\lim \limits_{T_{s} \rightarrow 0} \frac{\boldsymbol{X}_{k}-\boldsymbol{X}_{k-1}}{T_{s}}=\lim \limits_{T_{s} \rightarrow 0} \boldsymbol{F}\left(t_{k-1}\right) \hat{\boldsymbol{X}}_{k-1}+\frac{\boldsymbol{K}_{k}}{T_{s}}\left[\boldsymbol{Z}_{k}-\boldsymbol{H}_{k} \hat{\boldsymbol{X}}_{k-1}-\boldsymbol{H}_{k} \boldsymbol{F}\left(t_{k-1}\right) \hat{\boldsymbol{X}}_{k-1} T_{s}\right] \\ & =\boldsymbol{F}(t) \hat{\boldsymbol{X}}(t)+\boldsymbol{K}(t)[\boldsymbol{Z}(t)-\boldsymbol{H}(t) \hat{\boldsymbol{X}}(t)]\end{aligned}
X^˙(t)=Ts→0limTsXk−Xk−1=Ts→0limF(tk−1)X^k−1+TsKk[Zk−HkX^k−1−HkF(tk−1)X^k−1Ts]=F(t)X^(t)+K(t)[Z(t)−H(t)X^(t)]
5、均方差阵的连续化
将一阶的状态转移矩阵带入:
P
k
=
(
I
−
K
k
H
k
)
(
Φ
k
/
k
−
1
P
k
−
1
Φ
k
/
k
−
1
T
+
Γ
k
−
1
Q
k
−
1
Γ
k
−
1
T
)
=
(
I
−
K
k
H
k
)
{
[
I
+
F
(
t
k
−
1
)
T
s
]
P
k
−
1
[
I
+
F
(
t
k
−
1
)
T
s
]
T
+
G
(
t
k
−
1
)
⋅
q
(
t
k
−
1
)
T
s
⋅
G
T
(
t
k
−
1
)
}
=
P
k
−
1
+
F
(
t
k
−
1
)
P
k
−
1
T
s
+
P
k
−
1
F
T
(
t
k
−
1
)
T
s
+
G
(
t
k
−
1
)
q
(
t
k
−
1
)
G
T
(
t
k
−
1
)
T
s
+
O
(
T
s
2
)
−
K
k
H
k
[
P
k
−
1
+
O
(
T
s
)
]
\begin{aligned} \boldsymbol{P}_{k} & =\left(\boldsymbol{I}-\boldsymbol{K}_{k} \boldsymbol{H}_{k}\right)\left(\boldsymbol{\Phi}_{k / k-1} \boldsymbol{P}_{k-1} \boldsymbol{\Phi}_{k / k-1}^{\mathrm{T}}+\boldsymbol{\Gamma}_{k-1} \boldsymbol{Q}_{k-1} \boldsymbol{\Gamma}_{k-1}^{\mathrm{T}}\right) \\ & =\left(\boldsymbol{I}-\boldsymbol{K}_{k} \boldsymbol{H}_{k}\right)\left\{\left[\boldsymbol{I}+\boldsymbol{F}\left(t_{k-1}\right) T_{s}\right] \boldsymbol{P}_{k-1}\left[\boldsymbol{I}+\boldsymbol{F}\left(t_{k-1}\right) T_{s}\right]^{\mathrm{T}}+\boldsymbol{G}\left(t_{k-1}\right) \cdot \boldsymbol{q}\left(t_{k-1}\right) T_{s} \cdot \boldsymbol{G}^{\mathrm{T}}\left(t_{k-1}\right)\right\} \\ & ={\color{red}\boldsymbol{P}_{k-1}}+\boldsymbol{F}\left(t_{k-1}\right) \boldsymbol{P}_{k-1} T_{s}+\boldsymbol{P}_{k-1} \boldsymbol{F}^{\mathrm{T}}\left(t_{k-1}\right) T_{s}+\boldsymbol{G}\left(t_{k-1}\right) \boldsymbol{q}\left(t_{k-1}\right) \boldsymbol{G}^{\mathrm{T}}\left(t_{k-1}\right) T_{s}+O\left(T_{s}^{2}\right)-\boldsymbol{K}_{k} \boldsymbol{H}_{k}\left[{\boldsymbol{P}_{k-1}}+O\left(T_{s}\right)\right]\end{aligned}
Pk=(I−KkHk)(Φk/k−1Pk−1Φk/k−1T+Γk−1Qk−1Γk−1T)=(I−KkHk){[I+F(tk−1)Ts]Pk−1[I+F(tk−1)Ts]T+G(tk−1)⋅q(tk−1)Ts⋅GT(tk−1)}=Pk−1+F(tk−1)Pk−1Ts+Pk−1FT(tk−1)Ts+G(tk−1)q(tk−1)GT(tk−1)Ts+O(Ts2)−KkHk[Pk−1+O(Ts)]
同理,将
P
k
−
1
{\color{red}P_{k-1}}
Pk−1 减到左边,
k
k
k 时刻的值,减去
k
−
1
k-1
k−1 时刻的值,除以
T
s
T_s
Ts,求
T
s
T_s
Ts 趋于
0
0
0 的极限,就是连续的方程:
P
˙
(
t
)
=
lim
T
s
→
0
P
k
−
P
k
−
1
T
s
=
lim
T
s
→
0
F
(
t
k
−
1
)
P
k
−
1
+
P
k
−
1
F
T
(
t
k
−
1
)
+
G
(
t
k
−
1
)
q
(
t
k
−
1
)
G
T
(
t
k
−
1
)
−
K
k
T
s
H
k
P
k
−
1
=
F
(
t
)
P
(
t
)
+
P
(
t
)
F
T
(
t
)
+
G
(
t
)
q
(
t
)
G
T
(
t
)
−
K
(
t
)
H
(
t
)
P
(
t
)
=
F
(
t
)
P
(
t
)
+
P
(
t
)
F
T
(
t
)
−
K
(
t
)
r
(
t
)
K
T
(
t
)
+
G
(
t
)
q
(
t
)
G
T
(
t
)
\begin{aligned} \dot{\boldsymbol{P}}(t) & =\lim \limits_{T_{s} \rightarrow 0} \frac{\boldsymbol{P}_{k}-\boldsymbol{P}_{k-1}}{T_{s}}=\lim \limits_{T_{s} \rightarrow 0} \boldsymbol{F}\left(t_{k-1}\right) \boldsymbol{P}_{k-1}+\boldsymbol{P}_{k-1} \boldsymbol{F}^{\mathrm{T}}\left(t_{k-1}\right)+\boldsymbol{G}\left(t_{k-1}\right) \boldsymbol{q}\left(t_{k-1}\right) \boldsymbol{G}^{\mathrm{T}}\left(t_{k-1}\right)-\frac{\boldsymbol{K}_{k}}{T_{s}} \boldsymbol{H}_{k} \boldsymbol{P}_{k-1} \\ & =\boldsymbol{F}(t) \boldsymbol{P}(t)+\boldsymbol{P}(t) \boldsymbol{F}^{\mathrm{T}}(t)+\boldsymbol{G}(t) \boldsymbol{q}(t) \boldsymbol{G}^{\mathrm{T}}(t)-\boldsymbol{K}(t) \boldsymbol{H}(t) \boldsymbol{P}(t) \\ & =\boldsymbol{F}(t) \boldsymbol{P}(t)+\boldsymbol{P}(t) \boldsymbol{F}^{\mathrm{T}}(t)-\boldsymbol{K}(t) \boldsymbol{r}(t) \boldsymbol{K}^{\mathrm{T}}(t)+\boldsymbol{G}(t) \boldsymbol{q}(t) \boldsymbol{G}^{\mathrm{T}}(t)\end{aligned}
P˙(t)=Ts→0limTsPk−Pk−1=Ts→0limF(tk−1)Pk−1+Pk−1FT(tk−1)+G(tk−1)q(tk−1)GT(tk−1)−TsKkHkPk−1=F(t)P(t)+P(t)FT(t)+G(t)q(t)GT(t)−K(t)H(t)P(t)=F(t)P(t)+P(t)FT(t)−K(t)r(t)KT(t)+G(t)q(t)GT(t)
此方程也称:矩阵黎卡蒂(Riccati)方程,下面简单介绍一下Riccati方程:
标量Riccati方程:
矩阵Riccati方程:
P
˙
(
t
)
=
F
(
t
)
P
(
t
)
+
P
(
t
)
F
T
(
t
)
−
P
(
t
)
H
T
(
t
)
r
−
1
(
t
)
H
(
t
)
P
(
t
)
+
G
(
t
)
q
(
t
)
G
T
(
t
)
\dot{\boldsymbol{P}}(t)=\boldsymbol{F}(t) \boldsymbol{P}(t)+\boldsymbol{P}(t) \boldsymbol{F}^{\mathrm{T}}(t)-\boldsymbol{P}(t) \boldsymbol{H}^{\mathrm{T}}(t) \boldsymbol{r}^{-1}(t) \boldsymbol{H}(t) \boldsymbol{P}(t)+\boldsymbol{G}(t) \boldsymbol{q}(t) \boldsymbol{G}^{\mathrm{T}}(t)
P˙(t)=F(t)P(t)+P(t)FT(t)−P(t)HT(t)r−1(t)H(t)P(t)+G(t)q(t)GT(t)
Kalman滤波里的
P
P
P 一定是对称阵,转置还是本身,所以都没写
简记为: P ˙ = F P + P F F T − P R P + Q \dot{\boldsymbol{P}}=\boldsymbol{F P}+\boldsymbol{P F} \boldsymbol{F}^{\mathrm{T}}-\boldsymbol{P} \boldsymbol{R P}+\boldsymbol{Q} P˙=FP+PFFT−PRP+Q
它的解等价于 (化非线性问题为线性求解):
P
=
Y
D
−
1
其中
[
Y
˙
D
˙
]
=
[
F
Q
R
−
F
T
]
[
Y
D
]
\boldsymbol{P}=\boldsymbol{Y} \boldsymbol{D}^{-1} \quad \text { 其中 } \quad\left[\begin{array}{c} \dot{\boldsymbol{Y}} \\ \dot{\boldsymbol{D}} \end{array}\right]=\left[\begin{array}{cc} \boldsymbol{F} & \boldsymbol{Q} \\ \boldsymbol{R} & -\boldsymbol{F}^{\mathrm{T}} \end{array}\right]\left[\begin{array}{l} \boldsymbol{Y} \\ \boldsymbol{D} \end{array}\right]
P=YD−1 其中 [Y˙D˙]=[FRQ−FT][YD]
验证如下:
Y = P D Y ˙ = P ˙ D + P D ˙ P ˙ = Y ˙ D − 1 − P D ˙ D − 1 = ( F Y + Q D ) D − 1 − P ( R Y − F T D ) D − 1 = F Y D − 1 + Q − P R Y D − 1 + P F T = F P + Q − P R P + P F T \begin{aligned} \boldsymbol{Y} & =\boldsymbol{P D} \\ \dot{\boldsymbol{Y}} & =\dot{\boldsymbol{P}} \boldsymbol{D}+\boldsymbol{P} \dot{\boldsymbol{D}} \\ \dot{\boldsymbol{P}} & =\dot{\boldsymbol{Y}} \boldsymbol{D}^{-1}-\boldsymbol{P} \dot{\boldsymbol{D}} \boldsymbol{D}^{-1}=(\boldsymbol{F} \boldsymbol{Y}+\boldsymbol{Q D}) \boldsymbol{D}^{-1}-\boldsymbol{P}\left(\boldsymbol{R} \boldsymbol{Y}-\boldsymbol{F}^{\mathrm{T}} \boldsymbol{D}\right) \boldsymbol{D}^{-1} \\ & =\boldsymbol{F} \boldsymbol{Y} \boldsymbol{D}^{-1}+\boldsymbol{Q}-\boldsymbol{P} \boldsymbol{R} \boldsymbol{Y} \boldsymbol{D}^{-1}+\boldsymbol{P} \boldsymbol{F}^{\mathrm{T}}=\boldsymbol{F P}+\boldsymbol{Q}-\boldsymbol{P} \boldsymbol{R} \boldsymbol{P}+\boldsymbol{P} \boldsymbol{F}^{\mathrm{T}}\end{aligned} YY˙P˙=PD=P˙D+PD˙=Y˙D−1−PD˙D−1=(FY+QD)D−1−P(RY−FTD)D−1=FYD−1+Q−PRYD−1+PFT=FP+Q−PRP+PFT
但线性化后此矩阵微分方程没有初等解,只有毕卡积分级数解:
时变状态/矩阵方程 X ˙ ( t ) = F ( t ) X ( t ) \dot{\boldsymbol{X}}(t)=\boldsymbol{F}(t) \boldsymbol{X}(t) X˙(t)=F(t)X(t) 的毕卡积分级数解:
上面方程两边积分得:
X ( t ) = X ( 0 ) + ∫ 0 t F ( τ ) X ( τ ) d τ \boldsymbol{X}(t) =\boldsymbol{X}(0)+\int_{0}^{t} \boldsymbol{F}(\tau) {\color{red}\boldsymbol{X}(\tau)} \mathrm{d} \tau X(t)=X(0)+∫0tF(τ)X(τ)dτ
但此方程右边还含有未知量,迭代一次,将 X ( τ ) {\color{red}\boldsymbol{X}(\tau)} X(τ) 替换成上式
X ( t ) = X ( 0 ) + ∫ 0 t F ( τ ) [ X ( 0 ) + ∫ 0 τ F ( τ 1 ) X ( τ 1 ) d τ 1 ] d τ \boldsymbol{X}(t) =\boldsymbol{X}(0)+\int_{0}^{t} \boldsymbol{F}(\tau)\left[\boldsymbol{X}(0)+\int_{0}^{\tau} \boldsymbol{F}\left(\tau_{1}\right) {\color{red}\boldsymbol{X}\left(\tau_{1}\right) }\mathrm{d} \tau_{1}\right] \mathrm{d} \tau X(t)=X(0)+∫0tF(τ)[X(0)+∫0τF(τ1)X(τ1)dτ1]dτ
里面又含未知变量 X ( τ 1 ) {\color{red}\boldsymbol{X}\left(\tau_{1}\right) } X(τ1) ,再进行替换。反复迭代,一直存在未知变量,最后就变成无穷级数,而且还是无穷重积分
= X ( 0 ) + ∫ 0 t F ( τ ) d τ X ( 0 ) + ∫ 0 t F ( τ ) ∫ 0 τ F ( τ 1 ) X ( τ 1 ) d τ 1 d τ = ⋯ = [ I + ∫ 0 t F ( τ ) d τ + ∫ 0 t F ( τ ) ∫ 0 τ F ( τ 1 ) d τ 1 d τ + ∫ 0 t F ( τ ) ∫ 0 τ F ( τ 1 ) ∫ 0 τ 1 F ( τ 2 ) d τ 2 d τ 1 d τ + ⋯ ] X ( 0 ) \begin{aligned} & =\boldsymbol{X}(0)+\int_{0}^{t} \boldsymbol{F}(\tau) \mathrm{d} \tau \boldsymbol{X}(0)+\int_{0}^{t} \boldsymbol{F}(\tau) \int_{0}^{\tau} \boldsymbol{F}\left(\tau_{1}\right) \boldsymbol{X}\left(\tau_{1}\right) \mathrm{d} \tau_{1} \mathrm{~d} \tau=\cdots \\ & =\left[\boldsymbol{I}+\int_{0}^{t} \boldsymbol{F}(\tau) \mathrm{d} \tau+\int_{0}^{t} \boldsymbol{F}(\tau) \int_{0}^{\tau} \boldsymbol{F}\left(\tau_{1}\right) \mathrm{d} \tau_{1} \mathrm{~d} \tau+\int_{0}^{t} \boldsymbol{F}(\tau) \int_{0}^{\tau} \boldsymbol{F}\left(\tau_{1}\right) \int_{0}^{\tau_{1}} \boldsymbol{F}\left(\tau_{2}\right) \mathrm{d} \tau_{2} \mathrm{~d} \tau_{1} \mathrm{~d} \tau+\cdots\right] \boldsymbol{X}(0) \\ \end{aligned} =X(0)+∫0tF(τ)dτX(0)+∫0tF(τ)∫0τF(τ1)X(τ1)dτ1 dτ=⋯=[I+∫0tF(τ)dτ+∫0tF(τ)∫0τF(τ1)dτ1 dτ+∫0tF(τ)∫0τF(τ1)∫0τ1F(τ2)dτ2 dτ1 dτ+⋯]X(0)
前面括号内的矩阵就是从 0 0 0 时刻转移到 t t t 时刻的状态转移矩阵:
Φ ( t , 0 ) = I + ∫ 0 t F ( τ ) d τ + ∫ 0 t F ( τ ) ∫ 0 τ F ( τ 1 ) d τ 1 d τ + ∫ 0 t F ( τ ) ∫ 0 τ F ( τ 1 ) ∫ 0 τ 1 F ( τ 2 ) d τ 2 d τ 1 d τ + ⋯ \boldsymbol{\Phi}(t, 0) =\boldsymbol{I}+\int_{0}^{t} \boldsymbol{F}(\tau) \mathrm{d} \tau+\int_{0}^{t} \boldsymbol{F}(\tau) \int_{0}^{\tau} \boldsymbol{F}\left(\tau_{1}\right) \mathrm{d} \tau_{1} \mathrm{~d} \tau+\int_{0}^{t} \boldsymbol{F}(\tau) \int_{0}^{\tau} \boldsymbol{F}\left(\tau_{1}\right) \int_{0}^{\tau_{1}} \boldsymbol{F}\left(\tau_{2}\right) \mathrm{d} \tau_{2} \mathrm{~d} \tau_{1} \mathrm{~d} \tau+\cdots Φ(t,0)=I+∫0tF(τ)dτ+∫0tF(τ)∫0τF(τ1)dτ1 dτ+∫0tF(τ)∫0τF(τ1)∫0τ1F(τ2)dτ2 dτ1 dτ+⋯可以省略高阶项进行近似,或者有时候高阶项就为 0 0 0
一般情况下不可交换 F ( τ ) F ( τ 1 ) ≠ F ( τ 1 ) F ( τ ) \boldsymbol{F}(\tau) \boldsymbol{F}\left(\tau_{1}\right) \neq \boldsymbol{F}\left(\tau_{1}\right) \boldsymbol{F}(\tau) F(τ)F(τ1)=F(τ1)F(τ), 上式即为最终毕卡级数解;
可交换时 F ( τ ) F ( τ 1 ) = F ( τ 1 ) F ( τ ) \boldsymbol{F}(\tau) \boldsymbol{F}\left(\tau_{1}\right)=\boldsymbol{F}\left(\tau_{1}\right) \boldsymbol{F}(\tau) F(τ)F(τ1)=F(τ1)F(τ) (特殊如常值 F ( τ ) = F ) \left.\boldsymbol{F}(\tau)=\boldsymbol{F}\right) F(τ)=F), 才有闭合解:
可交换时,多重积分等于单重积分的 n n n 次方,除以 n n n 的阶乘分之一
∫ 0 τ 1 F ( τ 1 ) ⋯ ∫ 0 τ n F ( τ n ) d τ n ⋯ d τ 1 = 1 n ! [ ∫ 0 τ 1 F ( τ ) d τ 1 ] n Φ ( t , 0 ) = I + 1 1 ! ∫ 0 t F ( τ ) d τ + 1 2 ! [ ∫ 0 t F ( τ ) d τ ] 2 + 1 3 ! [ ∫ 0 t F ( τ ) d τ ] 3 + ⋯ = e ∫ 0 t F ( τ ) d τ 40 \begin{array}{l} \int_{0}^{\tau_{1}} \boldsymbol{F}\left(\tau_{1}\right) \cdots \int_{0}^{\tau_{n}} \boldsymbol{F}\left(\tau_{n}\right) \mathrm{d} \tau_{n} \cdots \mathrm{d} \tau_{1}=\frac{1}{n !}\left[\int_{0}^{\tau_{1}} \boldsymbol{F}(\tau) \mathrm{d} \tau_{1}\right]^{n} \\ \boldsymbol{\Phi}(t, 0)=\boldsymbol{I}+\frac{1}{1 !} \int_{0}^{t} \boldsymbol{F}(\tau) \mathrm{d} \tau+\frac{1}{2 !}\left[\int_{0}^{t} \boldsymbol{F}(\tau) \mathrm{d} \tau\right]^{2}+\frac{1}{3 !}\left[\int_{0}^{t} \boldsymbol{F}(\tau) \mathrm{d} \tau\right]^{3}+\cdots=\mathrm{e}^{\int_{0}^{t} F(\tau) \mathrm{d} \tau} 40 \end{array} ∫0τ1F(τ1)⋯∫0τnF(τn)dτn⋯dτ1=n!1[∫0τ1F(τ)dτ1]nΦ(t,0)=I+1!1∫0tF(τ)dτ+2!1[∫0tF(τ)dτ]2+3!1[∫0tF(τ)dτ]3+⋯=e∫0tF(τ)dτ40
6、连续时间Kalman滤波方程汇总
K ( t ) = P ( t ) H T ( t ) r − 1 ( t ) X ^ ˙ ( t ) = F ( t ) X ^ ( t ) + K ( t ) [ Z ( t ) − H ( t ) X ^ ( t ) ] P ˙ ( t ) = F ( t ) P ( t ) + P ( t ) F T ( t ) − K ( t ) r ( t ) K T ( t ) + G ( t ) q ( t ) G T ( t ) \begin{array}{l}\boldsymbol{K}(t)=\boldsymbol{P}(t) \boldsymbol{H}^{\mathrm{T}}(t) \boldsymbol{r}^{-1}(t) \\ \dot{\hat{\boldsymbol{X}}}(t)=\boldsymbol{F}(t) \hat{\boldsymbol{X}}(t)+\boldsymbol{K}(t)[\boldsymbol{Z}(t)-\boldsymbol{H}(t) \hat{\boldsymbol{X}}(t)] \\ \dot{\boldsymbol{P}}(t)=\boldsymbol{F}(t) \boldsymbol{P}(t)+\boldsymbol{P}(t) \boldsymbol{F}^{\mathrm{T}}(t)-\boldsymbol{K}(t) \boldsymbol{r}(t) \boldsymbol{K}^{\mathrm{T}}(t)+\boldsymbol{G}(t) \boldsymbol{q}(t) \boldsymbol{G}^{\mathrm{T}}(t) \\\end{array} K(t)=P(t)HT(t)r−1(t)X^˙(t)=F(t)X^(t)+K(t)[Z(t)−H(t)X^(t)]P˙(t)=F(t)P(t)+P(t)FT(t)−K(t)r(t)KT(t)+G(t)q(t)GT(t)
连续时间Kalman滤波方程中的状态估计和方差都是微分方程,相求结果得解微分方程。
离散时间Kalman滤波方程中的状态估计和方差可以认为是差分方程。
与确定性系统的状态观测器非常像