文章目录
- 一、平方根滤波基本形式
- 二、Potter平方根滤波
- 1、方差阵的量测更新
- 2、方差阵的时间更新
- 3、Potter平方根滤波流程
- 4、向量量测情况下的方差阵量测更新
- 三、奇异值(SVD)分解滤波
- 1、时间更新方差方程的SVD分解
- 2、量测更新方差方程的SVD分解
- 3、SVD分解滤波流程
- 四、UD分解滤波
- 1、量测更新方差方程UD分解
- 2、时间更新方差方程UD分解
- 1.朴素算法
- 2.快速算法
- 五、平方根信息滤波SRIKF
常见矩阵分解简介:
- LU分解: M = L U M=LU M=LU , L L L 是下三角矩阵, U U U 是上三角矩阵。方便求矩阵的行列式和逆,解线性方程组。
- Cholesky分解: M = L T L M=L^TL M=LTL , L L L 是上三角矩阵。是LU分解的进阶版。但是不同于LU分解的是,Cholesky分解只适用于正定的对称阵。在复数域,Cholesky分解要求矩阵是正定的共轭对称矩阵。要求这么多,就是因为Cholesky分解可以看成是给矩阵开平方根。
- QR分解: M = Q R M=QR M=QR ,其中 Q Q Q 是正交矩阵, R R R 是上三角矩阵。主要有三种方法:Gram-Schmidt正交化法(QR分解中的Q本身就可以看作是正交化构造出来的),Household变换法,Givens变换法。 可以用来求矩阵的特征值以及求解最小二乘法。
- 特征值分解: M = V D V T M=VDV^T M=VDVT ,其中 V V V 是正交阵, D D D 是由 M M M 的特征值构成的对角矩阵。只适用于方阵,目的是提取出矩阵最重要的特征。
- 奇异值分解(SVD分解): M = U S V T M=USV^T M=USVT ,其中 U U U、 V V V 是正交矩阵, S S S 是由 M M M 的奇异值构成的对角矩阵。特征值分解只适用于方阵,对于普通矩阵,则可以采用奇异值分解,提取出奇异值。
- UD分解: M = U D U T M=UDU^T M=UDUT, U U U 为上三角且对角线元素为 1 1 1, D D D 为对角阵
一、平方根滤波基本形式
在计算机中,单精度浮点数(float)有效数字为7位,双精度 (double)为16位,为了提高前者环境下的均方差阵计算精度,须采用平方根滤波。Kalman滤波计算过程中,有三个方差阵,分别可以记为:
P
k
−
1
=
U
k
−
1
U
k
−
1
T
P
k
/
k
−
1
=
Δ
k
/
k
−
1
Δ
k
/
k
−
1
T
P
k
=
U
k
Δ
k
T
\begin{array}{l} \boldsymbol{{P}_{k-1}}=\boldsymbol{U}_{k-1} \boldsymbol{U}_{k-1}^{\mathrm{T}}\\\boldsymbol{P}_{k / k-1}=\boldsymbol{\Delta}_{k / k-1} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \\\boldsymbol{P}_{k}=\boldsymbol{U}_{k} \boldsymbol{\Delta}_{k}^{\mathrm{T}} \end{array}
Pk−1=Uk−1Uk−1TPk/k−1=Δk/k−1Δk/k−1TPk=UkΔkT
方差都可以表示成平方的形式,标量的平方直接拆成中误差的平方就行,矩阵的平方可以拆成一个矩阵乘以它的转置,最常见是表示成下三角阵乘上三角阵形式(Cholesky分解)。
标准Kalman滤波的流程图可以表示成:
将图中的 P 、 Q 、 R P、Q、R P、Q、R 都用平方根形式替换 ,得:
- 更新都用平方根,损失精度小。
- 每一次滤波只计算一次增益矩阵 K K K ,用于左边的滤波计算回路,而对右边的增益计算回路无影响,可以损失一点精度。
二、Potter平方根滤波
1、方差阵的量测更新
考虑量测为标量的情况,由标准Kalman滤波方差阵递推公式:
如果不是标量量测,就用序贯滤波的方法,一个个做滤波
K k = P k / k − 1 H k T ( H k P k / k − 1 H k T + R k ) − 1 P k = ( I − K k H k ) P k / k − 1 \begin{array}{l}\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}}+R_{k}\right)^{-1} \\ \boldsymbol{P}_{k}=\left(\boldsymbol{I}-\boldsymbol{K}_{k} \boldsymbol{H}_{k}\right) \boldsymbol{P}_{k / k-1}\end{array} Kk=Pk/k−1HkT(HkPk/k−1HkT+Rk)−1Pk=(I−KkHk)Pk/k−1
将
K
K
K 带入第二个公式:
P
k
=
P
k
/
k
−
1
−
P
k
/
k
−
1
H
k
T
(
H
k
P
k
/
k
−
1
H
k
T
+
R
k
)
−
1
H
k
P
k
/
k
−
1
\boldsymbol{P}_{k}=\boldsymbol{P}_{k / k-1}-\boldsymbol{P}_{k / k-1} \boldsymbol{H}_{k}^{\mathrm{T}}\left(\boldsymbol{H}_{k} \boldsymbol{P}_{k / k-1} \boldsymbol{H}_{k}^{\mathrm{T}}+R_{k}\right)^{-1} \boldsymbol{H}_{k} \boldsymbol{P}_{k / k-1}
Pk=Pk/k−1−Pk/k−1HkT(HkPk/k−1HkT+Rk)−1HkPk/k−1
将方差阵Cholesky分解,表示成平方根形式:
Δ
k
Δ
k
T
=
Δ
k
/
k
−
1
Δ
k
/
k
−
1
T
−
Δ
k
/
k
−
1
Δ
k
/
k
−
1
T
H
k
T
(
H
k
Δ
k
/
k
−
1
Δ
k
/
k
−
1
T
H
k
T
+
R
k
)
−
1
H
k
Δ
k
/
k
−
1
Δ
k
/
k
−
1
T
=
Δ
k
/
k
−
1
[
I
−
Δ
k
/
k
−
1
T
H
k
T
(
H
k
Δ
k
/
k
−
1
Δ
k
/
k
−
1
T
H
k
T
+
R
k
)
−
1
H
k
Δ
k
/
k
−
1
]
Δ
k
/
k
−
1
T
\begin{array}{l} \boldsymbol{\Delta}_{k} \boldsymbol{\Delta}_{k}^{\mathrm{T}}=\boldsymbol{\Delta}_{k / k-1} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}}-\boldsymbol{\Delta}_{k / k-1} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}}\left(\boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}}+R_{k}\right)^{-1} \boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \\ =\boldsymbol{\Delta}_{k / k-1}\left[\boldsymbol{I}-\boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}}{\color{red}\left(\boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}}+R_{k}\right)}^{-1} \boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1}\right] \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \\\end{array}
ΔkΔkT=Δk/k−1Δk/k−1T−Δk/k−1Δk/k−1THkT(HkΔk/k−1Δk/k−1THkT+Rk)−1HkΔk/k−1Δk/k−1T=Δk/k−1[I−Δk/k−1THkT(HkΔk/k−1Δk/k−1THkT+Rk)−1HkΔk/k−1]Δk/k−1T
其中,
H
k
Δ
k
/
k
−
1
Δ
k
/
k
−
1
T
H
k
T
+
R
k
\boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}}+R_{k}
HkΔk/k−1Δk/k−1THkT+Rk 是标量,记为
ρ
k
2
\rho_{k}^{2}
ρk2 ,可以在矩阵的前后随意移动,上式变为:
Δ
k
/
k
−
1
(
I
−
ρ
k
−
2
Δ
k
/
k
−
1
T
H
k
T
H
k
Δ
k
/
k
−
1
)
Δ
k
/
k
−
1
T
\boldsymbol{\Delta}_{k / k-1}{\color{green}\left(\boldsymbol{I}-\rho_{k}^{-2} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1}\right)} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \quad
Δk/k−1(I−ρk−2Δk/k−1THkTHkΔk/k−1)Δk/k−1T
其中,绿色部分 ( I − ρ k − 2 Δ k / k − 1 T H k T H k Δ k / k − 1 ) \left(\boldsymbol{I}-{\color{brown}\rho_{k}^{-2}} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1}\right) (I−ρk−2Δk/k−1THkTHkΔk/k−1) 的平方根分解:
取 γ k − 1 \gamma_{k}^{-1} γk−1 为待 定系数,写成分解的形式:
( I − γ k − 1 Δ k / k − 1 T H k T H k Δ k / k − 1 ) ( I − γ k − 1 Δ k / k − 1 T H k T H k Δ k / k − 1 ) T = I − 2 γ k − 1 Δ k / k − 1 T H k T H k Δ k / k − 1 + γ k − 2 Δ k / k − 1 T H k T H k Δ k / k − 1 Δ k / k − 1 T H k T H k Δ k / k − 1 = I − 2 γ k − 1 Δ k / k − 1 T H k T H k Δ k / k − 1 + γ k − 2 Δ k / k − 1 T H k T ( ρ k 2 − R k ) H k Δ k / k − 1 = I − [ 2 γ k − 1 − ( ρ k 2 − R k ) γ k − 2 ] Δ k / k − 1 T H k T H k Δ k / k − 1 \begin{aligned}(\boldsymbol{I}- & \left.\gamma_{k}^{-1} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1}\right)\left(\boldsymbol{I}-\gamma_{k}^{-1} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1}\right)^{\mathrm{T}} \\ & =\boldsymbol{I}-2 \gamma_{k}^{-1} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1}+\gamma_{k}^{-2} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1} \\ & =\boldsymbol{I}-2 \gamma_{k}^{-1} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1}+\gamma_{k}^{-2} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}}\left(\rho_{k}^{2}-R_{k}\right) \boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1} \\ & =\boldsymbol{I}-{\color{brown}\left[2 \gamma_{k}^{-1}-\left(\rho_{k}^{2}-R_{k}\right) \gamma_{k}^{-2}\right]} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1}\end{aligned} (I−γk−1Δk/k−1THkTHkΔk/k−1)(I−γk−1Δk/k−1THkTHkΔk/k−1)T=I−2γk−1Δk/k−1THkTHkΔk/k−1+γk−2Δk/k−1THkTHkΔk/k−1Δk/k−1THkTHkΔk/k−1=I−2γk−1Δk/k−1THkTHkΔk/k−1+γk−2Δk/k−1THkT(ρk2−Rk)HkΔk/k−1=I−[2γk−1−(ρk2−Rk)γk−2]Δk/k−1THkTHkΔk/k−1
最后的式子和原式很相似,只差了棕色的部分,令棕色部分相等,可解待定系数 γ k − 1 \gamma_{k}^{-1} γk−1 :
ρ k − 2 = 2 γ k − 1 − ( ρ k 2 − R k ) γ k − 2 \rho_{k}^{-2}=2 \gamma_{k}^{-1}-\left(\rho_{k}^{2}-R_{k}\right) \gamma_{k}^{-2} ρk−2=2γk−1−(ρk2−Rk)γk−2
由:
ρ k − 2 = 2 γ k − 1 − ( ρ k 2 − R k ) γ k − 2 γ k 2 − 2 ρ k 2 γ k + ρ k 2 ( ρ k 2 − R k ) = 0 ρ k 2 = H k Δ k ∣ k − 1 Δ k ∣ k − 1 T H k T + R k \begin{array}{l} \rho_{k}^{-2}=2 \gamma_{k}^{-1}-\left(\rho_{k}^{2}-R_{k}\right) \gamma_{k}^{-2} \\ \gamma_{k}^{2}-2 \rho_{k}^{2} \gamma_{k}+\rho_{k}^{2}\left(\rho_{k}^{2}-R_{k}\right)=0 \quad \rho_{k}^{2}=\boldsymbol{H}_{k} \boldsymbol{\Delta}_{k \mid k-1} \boldsymbol{\Delta}_{k \mid k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}}+R_{k} \end{array} ρk−2=2γk−1−(ρk2−Rk)γk−2γk2−2ρk2γk+ρk2(ρk2−Rk)=0ρk2=HkΔk∣k−1Δk∣k−1THkT+Rk
由于是二次,所以 γ k − 1 \gamma_{k}^{-1} γk−1 有两个解:
γ k = 2 ρ k 2 ± 4 ρ k 4 − 4 ρ k 2 ( ρ k 2 − R k ) 2 = ρ k ( ρ k ± R k ) \gamma_{k}=\frac{2 \rho_{k}^{2} \pm \sqrt{4 \rho_{k}^{4}-4 \rho_{k}^{2}\left(\rho_{k}^{2}-R_{k}\right)}}{2}=\rho_{k}\left(\rho_{k} \pm \sqrt{R_{k}}\right) γk=22ρk2±4ρk4−4ρk2(ρk2−Rk)=ρk(ρk±Rk)
将改部分进行分解:
Δ
k
Δ
k
T
=
Δ
k
/
k
−
1
(
I
−
ρ
k
−
2
Δ
k
/
k
−
1
T
H
k
T
H
k
Δ
k
/
k
−
1
)
Δ
k
/
k
−
1
T
=
Δ
k
/
k
−
1
(
I
−
γ
k
−
1
Δ
k
/
k
−
1
T
H
k
T
H
k
Δ
k
/
k
−
1
)
(
I
−
γ
k
−
1
Δ
k
/
k
−
1
T
H
k
T
H
k
Δ
k
/
k
−
1
)
T
Δ
k
/
k
−
1
T
=
[
Δ
k
/
k
−
1
(
I
−
γ
k
−
1
Δ
k
/
k
−
1
T
H
k
T
H
k
Δ
k
/
k
−
1
)
]
[
Δ
k
/
k
−
1
(
I
−
γ
k
−
1
Δ
k
/
k
−
1
T
H
k
T
H
k
Δ
k
/
k
−
1
)
]
T
\begin{aligned} \boldsymbol{\Delta}_{k} \boldsymbol{\Delta}_{k}^{\mathrm{T}} & =\boldsymbol{\Delta}_{k / k-1}\left(\boldsymbol{I}-\rho_{k}^{-2} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1}\right) \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \\ & =\boldsymbol{\Delta}_{k / k-1}\left(\boldsymbol{I}-\gamma_{k}^{-1} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1}\right)\left(\boldsymbol{I}-\gamma_{k}^{-1} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1}\right)^{\mathrm{T}} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \\ & =\left[\boldsymbol{\Delta}_{k / k-1}\left(\boldsymbol{I}-\gamma_{k}^{-1} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1}\right)\right]\left[\boldsymbol{\Delta}_{k / k-1}\left(\boldsymbol{I}-\gamma_{k}^{-1} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1}\right)\right]^{\mathrm{T}} \end{aligned}
ΔkΔkT=Δk/k−1(I−ρk−2Δk/k−1THkTHkΔk/k−1)Δk/k−1T=Δk/k−1(I−γk−1Δk/k−1THkTHkΔk/k−1)(I−γk−1Δk/k−1THkTHkΔk/k−1)TΔk/k−1T=[Δk/k−1(I−γk−1Δk/k−1THkTHkΔk/k−1)][Δk/k−1(I−γk−1Δk/k−1THkTHkΔk/k−1)]T
得到方差阵的分解:
Δ
k
=
U
k
/
k
−
1
(
I
−
γ
k
−
1
Δ
k
/
k
−
1
T
H
k
T
H
k
Δ
k
/
k
−
1
)
\boldsymbol{\Delta}_{k}=\boldsymbol{U}_{k / k-1}\left(\boldsymbol{I}-\gamma_{k}^{-1} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{H}_{k} \boldsymbol{\Delta}_{k / k-1}\right)
Δk=Uk/k−1(I−γk−1Δk/k−1THkTHkΔk/k−1)
实现方差阵的量测更新:
Δ
k
/
k
−
1
,
R
k
1
/
2
⟶
Δ
k
\boldsymbol{\Delta}_{k / k-1}, \boldsymbol{R}_{k}^{1 / 2} \longrightarrow \boldsymbol{\Delta}_{k}
Δk/k−1,Rk1/2⟶Δk
2、方差阵的时间更新
由标准Kalman滤波的时间更新:
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
−
1
Δ
k
/
k
−
1
T
=
Φ
k
/
k
−
1
Δ
k
−
1
Δ
k
−
1
T
Φ
k
/
k
−
1
T
+
Γ
k
−
1
Q
k
−
1
1
2
(
Q
k
−
1
1
2
)
T
Γ
k
−
1
T
=
[
Φ
k
/
k
−
1
Δ
k
−
1
Γ
k
−
1
Q
k
−
1
1
2
]
[
Δ
k
−
1
T
Φ
k
/
k
−
1
T
(
Q
k
−
1
1
2
)
T
Γ
k
−
1
T
]
\begin{array}{c}\boldsymbol{\Delta}_{k / k-1} \boldsymbol{\Delta}_{k / k-1}^{\mathrm{T}}=\boldsymbol{\Phi}_{k / k-1} \boldsymbol{\Delta}_{k-1} \boldsymbol{\Delta}_{k-1}^{\mathrm{T}} \boldsymbol{\Phi}_{k / k-1}^{\mathrm{T}}+\boldsymbol{\Gamma}_{k-1} \boldsymbol{Q}_{k-1}^{\frac{1}{2}}\left(\boldsymbol{Q}_{k-1}^{\frac{1}{2}}\right)^{\mathrm{T}} \boldsymbol{\Gamma}_{k-1}^{\mathrm{T}} \\ =\left[\begin{array}{ll}\boldsymbol{\Phi}_{k / k-1} \boldsymbol{\Delta}_{k-1} & \boldsymbol{\Gamma}_{k-1} \boldsymbol{Q}_{k-1}^{\frac{1}{2}}\end{array}\right]\left[\begin{array}{c}\boldsymbol{\Delta}_{k-1}^{\mathrm{T}} \boldsymbol{\Phi}_{k / k-1}^{\mathrm{T}} \\ \left(\boldsymbol{Q}_{k-1}^{\frac{1}{2}}\right)^{\mathrm{T}} \boldsymbol{\Gamma}_{k-1}^{\mathrm{T}}\end{array}\right]\end{array}
Δk/k−1Δk/k−1T=Φk/k−1Δk−1Δk−1TΦk/k−1T+Γk−1Qk−121(Qk−121)TΓk−1T=[Φk/k−1Δk−1Γk−1Qk−121]
Δk−1TΦk/k−1T(Qk−121)TΓk−1T
不能直接将 Φ k / k − 1 Δ k − 1 Γ k − 1 Q k − 1 1 2 \boldsymbol{\Phi}_{k / k-1} \boldsymbol{\Delta}_{k-1} \boldsymbol{\Gamma}_{k-1} \boldsymbol{Q}_{k-1}^{\frac{1}{2}} Φk/k−1Δk−1Γk−1Qk−121 作为 Δ k / k − 1 \boldsymbol{\Delta}_{k / k-1} Δk/k−1 ,因为 A A T = B B T AA^T=BB^T AAT=BBT 时,不能说 A = B A=B A=B
QR分解:列满秩矩阵 A m ∗ n A_{m*n} Am∗n 总可以做QR分解: A m × n = Q ^ m × n R ^ n × n \boldsymbol{A}_{m \times n}=\hat{\boldsymbol{Q}}_{m \times n} \hat{\boldsymbol{R}}_{n \times n} Am×n=Q^m×nR^n×n 且有 Q ^ m × n T Q ^ m × n = I \hat{\boldsymbol{Q}}_{m \times n}^{\mathrm{T}} \hat{\boldsymbol{Q}}_{m \times n}=\boldsymbol{I} Q^m×nTQ^m×n=I , R ^ n × n \hat{\boldsymbol{R}}_{n \times n} R^n×n 是上三角可逆。
对上式QR分解:
(
Q
^
2
n
×
n
R
^
n
×
n
)
T
(
Q
^
2
n
×
n
R
^
n
×
n
)
=
R
n
×
n
T
R
^
n
×
n
\left(\hat{\boldsymbol{Q}}_{2 n \times n} \hat{\boldsymbol{R}}_{n \times n}\right)^{\mathrm{T}}\left({\hat{\boldsymbol{Q}}_{2 n \times n} \hat{\boldsymbol{R}}_{n \times n}}\right)={\boldsymbol{R}}_{n \times n}^{\mathrm{T}} \hat{\boldsymbol{R}}_{n \times n}
(Q^2n×nR^n×n)T(Q^2n×nR^n×n)=Rn×nTR^n×n
QR分解的改进:施密特正交化法,伪代码如下:
核心就是每个新的矢量都减去它在已经正交化的矢量方向的投影,进而每次新增一个新的正交矢量。新的矢量只和之前的矢量有关,而与后面的矢量无关。
先把 [ Δ k − 1 T Φ k / k − 1 T ( Q k − 1 1 2 ) T Γ k − 1 T ] \left[\begin{array}{c}\boldsymbol{\Delta}_{k-1}^{\mathrm{T}} \boldsymbol{\Phi}_{k / k-1}^{\mathrm{T}} \\ \left(\boldsymbol{Q}_{k-1}^{\frac{1}{2}}\right)^{\mathrm{T}} \boldsymbol{\Gamma}_{k-1}^{\mathrm{T}}\end{array}\right] Δk−1TΦk/k−1T(Qk−121)TΓk−1T 求出来,再用Gram-Schmidt法即可得到方差阵的时间更新 Δ k / k − 1 \boldsymbol{\Delta}_{k / k-1} Δk/k−1
3、Potter平方根滤波流程
4、向量量测情况下的方差阵量测更新
改量测更新就行,时间更新没必要改
前述标量量测情形
同理,向量量测情形:
全流程:
- 时间更新: Δ k − 1 ⟶ Q R Δ k / k − 1 \Delta_{k-1} \stackrel{\mathrm{QR}}{\longrightarrow} \Delta_{k / k-1} Δk−1⟶QRΔk/k−1 ,要做一次QR分解
- 量测更新: Δ k / k − 1 ⟶ Q R γ k ⟶ 求逆 Δ k \boldsymbol{\Delta}_{k / k-1} \stackrel{\mathrm{QR}}{\longrightarrow} \gamma_{k} \stackrel{\text { 求逆 }}{\longrightarrow} \boldsymbol{\Delta}_{k} Δk/k−1⟶QRγk⟶ 求逆 Δk ,也要做一次QR分解(其实前面标量的开方也是QR分解)
三、奇异值(SVD)分解滤波
奇异值分解可以参考博客
M = U S V T M=USV^T M=USVT,其中 U U U、 V V V 是正交矩阵, S S S 是由 M M M 的奇异值构成的对角矩阵。我们用奇异值分解的是方差阵,它对称正定。对称所以 U U U 和 V V V 相等,正定所以奇异值都大于 0 0 0
朴素分解方式:
上式每次更新要做两次QR分解,两次三角阵求逆;考虑用SVD分解,这样就只用做对角阵求逆。
1、时间更新方差方程的SVD分解
将方程的
P
k
/
k
−
1
P_{k/k-1}
Pk/k−1 分解为
U
k
/
k
−
1
Λ
k
/
k
−
1
U
k
/
k
−
1
T
\boldsymbol{U}_{k / k-1} \boldsymbol{\Lambda}_{k / k-1} \boldsymbol{U}_{k / k-1}^{\mathrm{T}}
Uk/k−1Λk/k−1Uk/k−1T ,
P
k
−
1
P_{k-1}
Pk−1 分解为
U
k
−
1
Λ
k
−
1
U
k
−
1
T
\boldsymbol{U}_{k-1} \boldsymbol{\Lambda}_{ k-1} \boldsymbol{U}_{ k-1}^{\mathrm{T}}
Uk−1Λk−1Uk−1T ,选择合适的
Γ
\boldsymbol{\Gamma}
Γ 噪声系数分配矩阵,可以认为
Q
Q
Q 是对角阵。
l
U
k
/
k
−
1
Λ
k
/
k
−
1
U
k
/
k
−
1
T
=
Φ
k
/
k
−
1
U
k
−
1
Λ
k
−
1
U
k
−
1
T
Φ
k
/
k
−
1
T
+
Γ
k
−
1
Q
k
−
1
Γ
k
−
1
T
=
[
Φ
k
/
k
−
1
U
k
−
1
Λ
k
−
1
1
2
Γ
k
−
1
Q
k
−
1
1
2
]
[
(
Λ
k
−
1
1
2
)
T
U
k
−
1
T
Φ
k
/
k
−
1
T
(
Q
k
−
1
1
2
)
T
Γ
k
−
1
T
]
\begin{aligned}{l} \boldsymbol{U}_{k / k-1} {\color{brown}\boldsymbol{\Lambda}_{k / k-1}} \boldsymbol{U}_{k / k-1}^{\mathrm{T}} & =\boldsymbol{\Phi}_{k / k-1} \boldsymbol{U}_{k-1} \boldsymbol{\Lambda}_{k-1} \boldsymbol{U}_{k-1}^{\mathrm{T}} \boldsymbol{\Phi}_{k / k-1}^{\mathrm{T}}+\boldsymbol{\Gamma}_{k-1} \boldsymbol{Q}_{k-1} \boldsymbol{\Gamma}_{k-1}^{\mathrm{T}} \\ & =\left[\begin{array}{lll}\boldsymbol{\Phi}_{k / k-1} \boldsymbol{U}_{k-1} \boldsymbol{\Lambda}_{k-1}^{\frac{1}{2}} & \boldsymbol{\Gamma}_{k-1} \boldsymbol{Q}_{k-1}^{\frac{1}{2}}\end{array}\right]\left[\begin{array}{c}\left(\boldsymbol{\Lambda}_{k-1}^{\frac{1}{2}}\right)^{\mathrm{T}} \boldsymbol{U}_{k-1}^{\mathrm{T}} \boldsymbol{\Phi}_{k / k-1}^{\mathrm{T}} \\ \left(\boldsymbol{Q}_{k-1}^{\frac{1}{2}}\right)^{\mathrm{T}} \boldsymbol{\Gamma}_{k-1}^{\mathrm{T}}\end{array}\right] \end{aligned}
lUk/k−1Λk/k−1Uk/k−1T=Φk/k−1Uk−1Λk−1Uk−1TΦk/k−1T+Γk−1Qk−1Γk−1T=[Φk/k−1Uk−1Λk−121Γk−1Qk−121]
(Λk−121)TUk−1TΦk/k−1T(Qk−121)TΓk−1T
再对矩阵做奇异值分解,得:
=
(
S
k
/
k
−
1
D
k
/
k
−
1
V
k
/
k
−
1
T
)
(
S
k
/
k
−
1
D
k
/
k
−
1
V
k
/
k
−
1
T
)
T
=
S
k
/
k
−
1
D
k
/
k
−
1
D
k
/
k
−
1
T
S
k
/
k
−
1
T
\begin{array}{l} =\left(\boldsymbol{S}_{k / k-1} \boldsymbol{D}_{k / k-1} \boldsymbol{V}_{k / k-1}^{\mathrm{T}}\right)\left(\boldsymbol{S}_{k / k-1} \boldsymbol{D}_{k / k-1} \boldsymbol{V}_{k / k-1}^{\mathrm{T}}\right)^{\mathrm{T}} \\ =\boldsymbol{S}_{k / k-1} {\color{brown}\boldsymbol{D}_{k / k-1} \boldsymbol{D}_{k / k-1}}^{\mathrm{T}} \boldsymbol{S}_{k / k-1}^{\mathrm{T}} \end{array}
=(Sk/k−1Dk/k−1Vk/k−1T)(Sk/k−1Dk/k−1Vk/k−1T)T=Sk/k−1Dk/k−1Dk/k−1TSk/k−1T
此式与原式形式一致,令:
U
k
/
k
−
1
=
S
k
/
k
−
1
,
Λ
k
/
k
−
1
1
2
=
D
k
/
k
−
1
\boldsymbol{U}_{k / k-1}=\boldsymbol{S}_{k / k-1}, \boldsymbol{\Lambda}_{k / k-1}^{\frac{1}{2}}=\boldsymbol{D}_{k / k-1}
Uk/k−1=Sk/k−1,Λk/k−121=Dk/k−1
2、量测更新方差方程的SVD分解
同样,将方程的
P
k
/
k
−
1
P_{k/k-1}
Pk/k−1 分解为
U
k
/
k
−
1
Λ
k
/
k
−
1
U
k
/
k
−
1
T
\boldsymbol{U}_{k / k-1} \boldsymbol{\Lambda}_{k / k-1} \boldsymbol{U}_{k / k-1}^{\mathrm{T}}
Uk/k−1Λk/k−1Uk/k−1T ,
P
k
P_{k}
Pk 分解为
U
k
Λ
k
U
k
T
\boldsymbol{U}_{k} \boldsymbol{\Lambda}_{ k} \boldsymbol{U}_{ k}^{\mathrm{T}}
UkΛkUkT ,得:
U
k
Λ
k
−
1
U
k
T
=
U
k
/
k
−
1
Λ
k
/
k
−
1
−
1
U
k
/
k
−
1
T
+
H
k
T
R
k
−
1
H
k
=
[
U
k
/
k
−
1
Λ
k
/
k
−
1
−
1
2
H
k
T
R
k
−
1
2
]
[
(
Λ
k
/
k
−
1
−
1
2
)
T
U
k
/
k
−
1
T
(
R
k
−
1
2
)
T
H
k
]
\boldsymbol{U}_{k} \boldsymbol{\Lambda}_{k}^{-1} \boldsymbol{U}_{k}^{\mathrm{T}} =\boldsymbol{U}_{k / k-1} \boldsymbol{\Lambda}_{k / k-1}^{-1} \boldsymbol{U}_{k / k-1}^{\mathrm{T}}+\boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{R}_{k}^{-1} \boldsymbol{H}_{k} \\ =\left[\begin{array}{ll}\boldsymbol{U}_{k / k-1} \boldsymbol{\Lambda}_{k / k-1}^{-\frac{1}{2}} & \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{R}_{k}^{-\frac{1}{2}}\end{array}\right]\left[\begin{array}{c}\left(\boldsymbol{\Lambda}_{k / k-1}^{-\frac{1}{2}}\right)^{\mathrm{T}} \boldsymbol{U}_{k / k-1}^{\mathrm{T}} \\ \left(\boldsymbol{R}_{k}^{-\frac{1}{2}}\right)^{\mathrm{T}} \boldsymbol{H}_{k}\end{array}\right]
UkΛk−1UkT=Uk/k−1Λk/k−1−1Uk/k−1T+HkTRk−1Hk=[Uk/k−1Λk/k−1−21HkTRk−21]
(Λk/k−1−21)TUk/k−1T(Rk−21)THk
对矩阵做SVD分解,得:
=
(
S
k
D
k
V
k
T
)
(
S
k
D
k
V
k
T
)
T
=
S
k
D
k
D
k
T
S
k
T
\begin{array}{l} =\left(\boldsymbol{S}_{k} \boldsymbol{D}_{k} \boldsymbol{V}_{k}^{\mathrm{T}}\right)\left(\boldsymbol{S}_{k} \boldsymbol{D}_{k} \boldsymbol{V}_{k}^{\mathrm{T}}\right)^{\mathrm{T}} =\boldsymbol{S}_{k} \boldsymbol{D}_{k} \boldsymbol{D}_{k}^{\mathrm{T}} \boldsymbol{S}_{k}^{\mathrm{T}} \end{array}
=(SkDkVkT)(SkDkVkT)T=SkDkDkTSkT
也令:
U
k
=
S
k
,
Λ
k
−
1
2
=
D
k
\boldsymbol{U}_{k}=\boldsymbol{S}_{k}, \boldsymbol{\Lambda}_{k}^{-\frac{1}{2}}=\boldsymbol{D}_{k}
Uk=Sk,Λk−21=Dk
3、SVD分解滤波流程
( U k − 1 , Λ k − 1 1 2 , Q k − 1 1 2 ) → [ Φ k / k − 1 U k − 1 Λ k − 1 1 2 Γ k − 1 Q k − 1 1 2 ] ⟶ SVD ( U k / k − 1 , Λ k / k − 1 1 2 ) → ( U k / k − 1 , Λ k / k − 1 − 1 2 , R k − 1 2 ) → [ U k / k − 1 Λ k / k − 1 − 1 2 H k T R k − 1 2 ] ⟶ S V D ( U k , Λ k − 1 2 ) → ⋯ \begin{array}{l}\left(\boldsymbol{U}_{k-1}, \boldsymbol{\Lambda}_{k-1}^{\frac{1}{2}}, \boldsymbol{Q}_{k-1}^{\frac{1}{2}}\right) \rightarrow\left[\boldsymbol{\Phi}_{k / k-1} \boldsymbol{U}_{k-1} \boldsymbol{\Lambda}_{k-1}^{\frac{1}{2}} \quad \boldsymbol{\Gamma}_{k-1} \boldsymbol{Q}_{k-1}^{\frac{1}{2}}\right] \stackrel{\operatorname{SVD}}{\longrightarrow}\left(\boldsymbol{U}_{k / k-1}, \boldsymbol{\Lambda}_{k / k-1}^{\frac{1}{2}}\right) \\ \rightarrow\left(\boldsymbol{U}_{k / k-1}, \boldsymbol{\Lambda}_{k / k-1}^{-\frac{1}{2}}, \boldsymbol{R}_{k}^{-\frac{1}{2}}\right) \rightarrow\left[\boldsymbol{U}_{k / k-1} \boldsymbol{\Lambda}_{k / k-1}^{-\frac{1}{2}} \quad \boldsymbol{H}_{k}^{\mathrm{T}} \boldsymbol{R}_{k}^{-\frac{1}{2}}\right] \stackrel{\mathrm{SVD}}{\longrightarrow}\left(\boldsymbol{U}_{k}, \boldsymbol{\Lambda}_{k}^{-\frac{1}{2}}\right) \rightarrow \cdots \\\end{array} (Uk−1,Λk−121,Qk−121)→[Φk/k−1Uk−1Λk−121Γk−1Qk−121]⟶SVD(Uk/k−1,Λk/k−121)→(Uk/k−1,Λk/k−1−21,Rk−21)→[Uk/k−1Λk/k−1−21HkTRk−21]⟶SVD(Uk,Λk−21)→⋯
每次更新含2次SVD分解、2次对角阵求逆,运算量与“朴素”方法相比没有明显优势(SVD分解计算量远大于QR),没啥用。
四、UD分解滤波
UD分解: M = U D U T M=UDU^T M=UDUT, U U U 为上三角且对角线元素为 1 1 1, D D D 为对角阵。存储的时候可以让 U U U 占据上三角(对角线都是1),对角线上都为 D D D 。
1、量测更新方差方程UD分解
必须是标量量测,向量的还没有人推导过,只能改成标量再分解
P k = P k / k − 1 − P k / k − 1 H k T ( H k P k / k − 1 H k T + R k ) − 1 H k P k / k − 1 \boldsymbol{P}_{k}=\boldsymbol{P}_{k / k-1}-\boldsymbol{P}_{k / k-1} \boldsymbol{H}_{k}^{\mathrm{T}}\left(\boldsymbol{H}_{k} \boldsymbol{P}_{k / k-1} \boldsymbol{H}_{k}^{\mathrm{T}}+R_{k}\right)^{-1} \boldsymbol{H}_{k} \boldsymbol{P}_{k / k-1} Pk=Pk/k−1−Pk/k−1HkT(HkPk/k−1HkT+Rk)−1HkPk/k−1
将方程的
P
k
/
k
−
1
P_{k/k-1}
Pk/k−1 分解为
U
k
/
k
−
1
D
k
/
k
−
1
U
k
/
k
−
1
T
\boldsymbol{U}_{k / k-1} \boldsymbol{D}_{k / k-1} \boldsymbol{U}_{k / k-1}^{\mathrm{T}}
Uk/k−1Dk/k−1Uk/k−1T ,
P
k
P_{k}
Pk 分解为
U
k
D
k
U
k
T
\boldsymbol{U}_{k} \boldsymbol{D}_{ k} \boldsymbol{U}_{ k}^{\mathrm{T}}
UkDkUkT,
R
R
R 为标量可以前后移动,得:
U
k
D
k
U
k
T
=
U
k
/
k
−
1
D
k
/
k
−
1
U
k
/
k
−
1
T
−
U
k
/
k
−
1
D
k
/
k
−
1
U
k
/
k
−
1
T
H
k
T
(
H
k
U
k
/
k
−
1
D
k
/
k
−
1
U
k
/
k
−
1
T
H
k
T
+
R
k
)
−
1
H
k
U
k
/
k
−
1
D
k
/
k
−
1
U
k
/
k
−
1
T
=
U
k
/
k
−
1
[
D
k
/
k
−
1
−
D
k
/
k
−
1
U
k
/
k
−
1
T
H
k
T
(
H
k
U
k
/
k
−
1
D
k
/
k
−
1
U
k
/
k
−
1
T
H
k
T
+
R
k
)
−
1
H
k
U
k
/
k
−
1
D
k
/
k
−
1
]
U
k
/
k
−
1
T
=
U
k
/
k
−
1
(
D
k
/
k
−
1
−
α
−
1
g
g
T
)
U
k
/
k
−
1
T
记
{
α
=
H
k
U
k
/
k
−
1
⋅
D
k
/
k
−
1
U
k
/
k
−
1
T
H
k
T
+
R
k
=
f
T
g
+
R
k
f
=
(
H
k
U
k
/
k
−
1
)
T
g
=
D
k
/
k
−
1
U
k
/
k
−
1
T
H
k
T
=
D
k
/
k
−
1
f
\begin{array}{l} \boldsymbol{U}_{k} \boldsymbol{D}_{k} \boldsymbol{U}_{k}^{\mathrm{T}}= \boldsymbol{U}_{k / k-1} \boldsymbol{D}_{k / k-1} \boldsymbol{U}_{k / k-1}^{\mathrm{T}}- \\ \boldsymbol{U}_{k / k-1} \boldsymbol{D}_{k / k-1} \boldsymbol{U}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}}\left(\boldsymbol{H}_{k} \boldsymbol{U}_{k / k-1} \boldsymbol{D}_{k / k-1} \boldsymbol{U}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}}+R_{k}\right)^{-1} \boldsymbol{H}_{k} \boldsymbol{U}_{k / k-1} \boldsymbol{D}_{k / k-1} \boldsymbol{U}_{k / k-1}^{\mathrm{T}} \\ = \boldsymbol{U}_{k / k-1}\left[\boldsymbol{D}_{k / k-1}-\boldsymbol{D}_{k / k-1} \boldsymbol{U}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}}\left(\boldsymbol{H}_{k} \boldsymbol{U}_{k / k-1} \boldsymbol{D}_{k / k-1} \boldsymbol{U}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}}+R_{k}\right)^{-1} \boldsymbol{H}_{k} \boldsymbol{U}_{k / k-1} \boldsymbol{D}_{k / k-1}\right] \boldsymbol{U}_{k / k-1}^{\mathrm{T}} \\ = \boldsymbol{U}_{k / k-1}\left(\boldsymbol{D}_{k / k-1}-\alpha^{-1} \boldsymbol{g} \boldsymbol{g}^{\mathrm{T}}\right) \boldsymbol{U}_{k / k-1}^{\mathrm{T}} \\ \text { 记 }\left\{\begin{array}{l}\alpha=\boldsymbol{H}_{k} \boldsymbol{U}_{k / k-1} \cdot \boldsymbol{D}_{k / k-1} \boldsymbol{U}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}}+R_{k}=\boldsymbol{f}^{\mathrm{T}} \boldsymbol{g}+R_{k} \\ \boldsymbol{f}=\left(\boldsymbol{H}_{k} \boldsymbol{U}_{k / k-1}\right)^{\mathrm{T}} \\ \boldsymbol{g}=\boldsymbol{D}_{k / k-1} \boldsymbol{U}_{k / k-1}^{\mathrm{T}} \boldsymbol{H}_{k}^{\mathrm{T}}=\boldsymbol{D}_{k / k-1} \boldsymbol{f}\end{array}\right.\end{array}
UkDkUkT=Uk/k−1Dk/k−1Uk/k−1T−Uk/k−1Dk/k−1Uk/k−1THkT(HkUk/k−1Dk/k−1Uk/k−1THkT+Rk)−1HkUk/k−1Dk/k−1Uk/k−1T=Uk/k−1[Dk/k−1−Dk/k−1Uk/k−1THkT(HkUk/k−1Dk/k−1Uk/k−1THkT+Rk)−1HkUk/k−1Dk/k−1]Uk/k−1T=Uk/k−1(Dk/k−1−α−1ggT)Uk/k−1T 记 ⎩
⎨
⎧α=HkUk/k−1⋅Dk/k−1Uk/k−1THkT+Rk=fTg+Rkf=(HkUk/k−1)Tg=Dk/k−1Uk/k−1THkT=Dk/k−1f
从 D n n D_{nn} Dnn 开始先计算最后一列 ,再计算倒数第二列,直到 D 11 D_{11} D11 。
2、时间更新方差方程UD分解
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
将方程的
P
k
/
k
−
1
P_{k/k-1}
Pk/k−1 分解为
U
k
/
k
−
1
D
k
/
k
−
1
U
k
/
k
−
1
T
\boldsymbol{U}_{k / k-1} \boldsymbol{D}_{k / k-1} \boldsymbol{U}_{k / k-1}^{\mathrm{T}}
Uk/k−1Dk/k−1Uk/k−1T ,
P
k
−
1
P_{k-1}
Pk−1 分解为
U
k
−
1
D
k
−
1
U
k
−
1
T
\boldsymbol{U}_{k-1} \boldsymbol{D}_{ k-1} \boldsymbol{U}_{ k-1}^{\mathrm{T}}
Uk−1Dk−1Uk−1T ,最后两边的矩阵记为
W
k
∣
k
−
1
\boldsymbol{W}_{k \mid k-1}
Wk∣k−1
U
k
∣
k
−
1
D
k
/
k
−
1
U
k
/
k
−
1
T
=
Φ
k
l
k
−
1
U
k
−
1
D
k
−
1
U
k
−
1
T
Φ
k
/
k
−
1
T
+
Γ
k
−
1
Q
k
−
1
Γ
k
−
1
T
=
[
Φ
k
l
k
−
1
U
k
−
1
Γ
k
−
1
]
[
D
k
−
1
0
0
Q
k
−
1
]
[
U
k
−
1
T
Φ
k
l
k
−
1
T
Γ
k
−
1
T
]
≜
W
k
∣
k
−
1
D
~
k
−
1
W
k
∣
k
−
1
T
\begin{aligned} \boldsymbol{U}_{k \mid k-1} \boldsymbol{D}_{k / k-1} \boldsymbol{U}_{k / k-1}^{\mathrm{T}} & =\boldsymbol{\Phi}_{k l k-1} \boldsymbol{U}_{k-1} \boldsymbol{D}_{k-1} \boldsymbol{U}_{k-1}^{\mathrm{T}} \boldsymbol{\Phi}_{k / k-1}^{\mathrm{T}}+\boldsymbol{\Gamma}_{k-1} \boldsymbol{Q}_{k-1} \boldsymbol{\Gamma}_{k-1}^{\mathrm{T}} \\ & =\left[\begin{array}{ll} \boldsymbol{\Phi}_{k l k-1} \boldsymbol{U}_{k-1} & \boldsymbol{\Gamma}_{k-1} \end{array}\right]\left[\begin{array}{cc} \boldsymbol{D}_{k-1} & \mathbf{0} \\ \mathbf{0} & \boldsymbol{Q}_{k-1} \end{array}\right]\left[\begin{array}{c} \boldsymbol{U}_{k-1}^{\mathrm{T}} \boldsymbol{\Phi}_{k l k-1}^{\mathrm{T}} \\ \boldsymbol{\Gamma}_{k-1}^{\mathrm{T}} \end{array}\right] \triangleq \boldsymbol{W}_{k \mid k-1} \tilde{\boldsymbol{D}}_{k-1} \boldsymbol{W}_{k \mid k-1}^{\mathrm{T}} \end{aligned}
Uk∣k−1Dk/k−1Uk/k−1T=Φklk−1Uk−1Dk−1Uk−1TΦk/k−1T+Γk−1Qk−1Γk−1T=[Φklk−1Uk−1Γk−1][Dk−100Qk−1][Uk−1TΦklk−1TΓk−1T]≜Wk∣k−1D~k−1Wk∣k−1T
1.朴素算法
对
W
k
∣
k
−
1
\boldsymbol{W}_{k \mid k-1}
Wk∣k−1 再做一次QR分解:
W
k
l
k
−
1
D
~
k
−
1
W
k
/
k
−
1
T
=
R
^
T
Q
^
T
D
~
k
−
1
Q
R
^
=
R
^
T
A
R
^
\boldsymbol{W}_{k l k-1} \tilde{\boldsymbol{D}}_{k-1} \boldsymbol{W}_{k / k-1}^{\mathrm{T}}=\hat{\boldsymbol{R}}^{\mathrm{T}} {\color{green}\hat{\boldsymbol{Q}}^{\mathrm{T}} \tilde{\boldsymbol{D}}_{k-1} {\boldsymbol{Q}}} \hat{\boldsymbol{R}}=\hat{\boldsymbol{R}}^{\mathrm{T}} {\color{green}\boldsymbol{A}} \hat{\boldsymbol{R}}
Wklk−1D~k−1Wk/k−1T=R^TQ^TD~k−1QR^=R^TAR^
绿色部分正定对称,记为
A
A
A ,对其进行UD分解:
=
R
^
T
U
^
D
^
U
^
T
R
‾
=
(
R
^
T
U
^
)
D
^
(
R
^
T
U
^
)
T
=\hat{\boldsymbol{R}}^{\mathrm{T}} \hat{\boldsymbol{U}} \hat{\boldsymbol{D}} \hat{\boldsymbol{U}}^{\mathrm{T}} \overline{\boldsymbol{R}}=\left(\hat{\boldsymbol{R}}^{\mathrm{T}} \hat{\boldsymbol{U}}\right) \hat{\boldsymbol{D}}\left(\hat{\boldsymbol{R}}^{\mathrm{T}} \hat{\boldsymbol{U}}\right)^{\mathrm{T}}
=R^TU^D^U^TR=(R^TU^)D^(R^TU^)T
R
R
R 是QR分解出的三角阵,
U
U
U 是UD分解出的三角阵,乘出来的
R
T
U
R^TU
RTU 还是三角阵,还有保证对角线是
1
1
1,把这部分赋给
U
k
∣
k
−
1
\boldsymbol{U}_{k \mid k-1}
Uk∣k−1 ,不等于
1
1
1 的部分归到
D
D
D 里面。此方法要一次QR分解和一个UD分解,计算量大。
2.快速算法
用左边和右边元素一一对应得:
D
k
∣
k
−
1
,
j
j
=
∑
s
=
1
n
+
1
D
~
k
−
1
,
s
s
W
j
,
s
(
n
−
j
)
W
j
,
s
(
n
−
j
)
U
k
l
k
−
1
,
i
j
=
∑
s
=
1
n
+
1
D
~
k
−
1
,
s
s
W
i
,
s
(
n
−
j
)
W
j
,
s
(
n
−
j
)
D
k
/
k
−
1
,
j
j
D_{k \mid k-1, j j}=\sum_{s=1}^{n+1} \tilde{D}_{k-1, s s} W_{j, s}^{(n-j)} W_{j, s}^{(n-j)} \quad U_{k l k-1, i j}=\frac{\sum_{s=1}^{n+1} \tilde{D}_{k-1, s s} W_{i, s}^{(n-j)} W_{j, s}^{(n-j)}}{D_{k / k-1, j j}}
Dk∣k−1,jj=s=1∑n+1D~k−1,ssWj,s(n−j)Wj,s(n−j)Uklk−1,ij=Dk/k−1,jj∑s=1n+1D~k−1,ssWi,s(n−j)Wj,s(n−j)
滤波流程: ( U k − 1 , j , , D k − 1 , j ) → ( U k ∣ k − 1 , j j , D k l k − 1 , j ) → ( U k , j , D k , j , j ) \left(U_{k-1, j,}, D_{k-1, j}\right) \rightarrow\left(U_{k \mid k-1, j j}, D_{k l k-1, j}\right) \rightarrow\left(U_{k, j}, D_{k, j, j}\right) (Uk−1,j,,Dk−1,j)→(Uk∣k−1,jj,Dklk−1,j)→(Uk,j,Dk,j,j)
运算量比较小,计算比较紧凑,已经推导好了一个一个元素怎么计算,比较实用。
五、平方根信息滤波SRIKF
将信息矩阵 I I I 分解,与Potter平方根滤波很相似,计算量几乎一模一样,可以类比来看。