1 状态观测器
在工程上，部分状态参数测量的成本过高，或是用现有仪器无法测得，此时需要引入状态观测器。状态观测器以原系统输入、输出为输入，输出估计的状态变量。
原系统状态描述：
{
x
k
=Ax
k−1
+Bu
k−1
y
k
=Cx
k
状态观测器状态描述：
{ x ~ k = A x ~ k − 1 + B u k − 1 + f [ y k − 1 − y ~ k − 1 ] y ~ k = C x ~ k
{x~k=Ax~k−1+Buk−1+f[yk−1−y~k−1]y~k=Cx~k
{x~k=Ax~k−1+Buk−1+f[yk−1−y~k−1]y~k=Cx~k
{
x
~
k
=A
x
~
k−1
+Bu
k−1
+f[y
k−1
−
y
~
k−1
]
y
~
k
=C
x
~
k
定义误差向量 e k = x k − x ~ k e_k=x_k-\tilde{x}_k e
k
=x
k
−
x
~
k
，则两系统状态方程相减得
系 统 误 差 传 递 方 程 e k = ( A − f C ) e k − 1 系统误差传递方程e_k=\left( A-fC \right) e_{k-1}
系统误差传递方程e
k
=(A−fC)e
k−1
其通解为矩阵指数函数，所以当 A − f C A-fC A−fC特征值小于0时，误差向量 e k e_k e
k
各分量均趋于0，即状态观测器能直接估计原系统状态。
本质上，状态观测器是针对状态空间方程描述的确定系统，从错误的系统状态估计值不断收敛到正确的系统状态估计值的数学模型。
2 状态滤波器
实际的系统往往是随机过程，状态描述为：
{ x k = A x k − 1 + B u k − 1 + w k − 1 y k = C x k + v k
{xk=Axk−1+Buk−1+wk−1yk=Cxk+vk
{xk=Axk−1+Buk−1+wk−1yk=Cxk+vk
{
x
k
=Ax
k−1
+Bu
k−1
+w
k−1
y
k
=Cx
k
+v
k
其中 w k w_k w
k
称为过程噪声，主要由系统运行过程中的外力导致，如旋翼飞行器运动过程中受到的阵风； v k v_k v
k
称为测量噪声，主要由系统测量过程中测量仪器误差、精度不足导致。
本质上，状态滤波器是针对随机状态空间方程描述的随机系统，从不确定观测中提取信息，得到系统状态最优估计的数学模型。当 w k w_k w
k
与 v k v_k v
k
为高斯白噪声且相互独立时，状态滤波器为卡尔曼滤波器。
状态滤波器的核心是通过贝叶斯原理不断调整滤波器增益矩阵，以减小随机干扰，使估计的系统状态趋近于真实状态；而状态观测器的核心是通过极点配置确定观测器增益矩阵，使估计的状态跟踪当前系统的状态。
一般地，状态滤波器性能优于状态观测器，但在实际系统中，若随机噪声非高斯噪声，或有若干非线性环节，按常规噪声模型建立的状态滤波器可能发散，此时则更需要状态观测器的稳定性。
3 卡尔曼滤波器
机器学习：详解贝叶斯网络+例题分析曾提到贝叶斯方法的思考模式：
参数先验信息 π ( θ ) + 样本观测数据 X = 后验分布 P ( θ ∣ X ) \text{参数先验信息}\pi \left( \theta \right) +\text{样本观测数据}X=\text{后验分布}P\left( \theta |X \right)
参数先验信息π(θ)+样本观测数据X=后验分布P(θ∣X)
即在得到新的样本信息之前，人们对模型的认知是先验分布 π ( θ ) \pi \left( \theta \right) π(θ)，在得到新的样本信息 X X X后，人们对模型的认知为后验分布 P ( θ ∣ X ) P\left( \theta |X \right) P(θ∣X)。
由此引出卡尔曼滤波器。考虑一个随机系统：
{ x k = A x k − 1 + B u k − 1 + w k − 1 z k = C x k + v k
{xk=Axk−1+Buk−1+wk−1zk=Cxk+vk
{xk=Axk−1+Buk−1+wk−1zk=Cxk+vk
{
x
k
=Ax
k−1
+Bu
k−1
+w
k−1
z
k
=Cx
k
+v
k
其中 w k N ( 0 , Q ) w_k~N\left( 0, Q \right) w
k
N(0,Q)、 v k N ( 0 , R ) v_k~N\left( 0, R \right) v
k
N(0,R)，且 w k w_k w
k
与 v k v_k v
k
相互独立。
忽略噪声可建立该系统的先验数学模型(噪声仅满足随机分布，无法建模)：
状态预测方程 { x ^ k − = A x ^ k − 1 + B u k − 1 z ^ k − = C x ^ k {\text{状态预测方程}
{x^−k=Ax^k−1+Buk−1z^−k=Cx^k
{x^k−=Ax^k−1+Buk−1z^k−=Cx^k
}
状态预测方程{
x
^
k
−
=A
x
^
k−1
+Bu
k−1
z
^
k
−
=C
x
^
k
基于上述系统，从贝叶斯方法引出卡尔曼滤波器状态更新方程：
状态更新方程 : x ^ k = x ^ k − + K k ( z k − z ^ k − ) {\text{状态更新方程}: \hat{x}_k=\hat{x}_{k}^{-}+K_k\left( z_k-\hat{z}_{k}^{-} \right) }
状态更新方程:
x
^
k
=
x
^
k
−
+K
k
(z
k
−
z
^
k
−
)
其中 x ^ k \hat{x}_k
x
^
k
是后验状态估计，即经过一次修正后，当前系统状态的最优估计值； x ^ k − \hat{x}_{k}^{-}
x
^
k
−
为系统状态的先验估计； z k z_k z
k
为当前状态的测量样本值； z ^ k − \hat{z}_{k}^{-}
z
^
k
−
为当前状态的先验测量值； z k − z ^ k − z_k-\hat{z}_{k}^{-} z
k
−
z
^
k
−
代表实际系统与估计系统间随机误差的影响； K k K_k K
k
权衡先验模型与实测数据间对后验分布的关系，称为卡尔曼增益，显然， z k − z ^ k − z_k-\hat{z}_{k}^{-} z
k
−
z
^
k
−
越小表明实际与预测越接近，则后验分布越趋近于先验分布， z k − z ^ k − z_k-\hat{z}_{k}^{-} z
k
−
z
^
k
−
越大表明预测越不可信，需要通过实测来大幅修正先验分布。
下面要确定 K k K_k K
k
使后验模型的修正效果最佳，即 x ^ k \hat{x}_k
x
^
k
越接近真实值 x k {x}_k x
k
。定义误差向量
{ 后验误差向量 e k = x k − x ^ k 先验误差向量 e k − = x k − x ^ k −
{后验误差向量ek=xk−x^k先验误差向量e−k=xk−x^−k
{后验误差向量ek=xk−x^k先验误差向量ek−=xk−x^k−
{
后验误差向量e
k
=x
k
−
x
^
k
先验误差向量e
k
−
=x
k
−
x
^
k
−
定义后验误差协方差矩阵为:
P k = E ( e k e k T ) = [ σ e 1 2 σ e 1 σ e 2 ⋯ σ e 1 σ e n σ e 2 σ e 1 σ e 2 2 ⋯ σ e 2 σ e n ⋮ ⋮ ⋱ ⋮ σ e n σ e 1 σ e n σ e 2 ⋯ σ e n 2 ] P_k=E\left( e_ke_{k}^{T} \right) =\left[
σ2e1σe2σe1⋮σenσe1σe1σe2σ2e2⋮σenσe2⋯⋯⋱⋯σe1σenσe2σen⋮σ2en
σe12σe1σe2⋯σe1σenσe2σe1σe22⋯σe2σen⋮⋮⋱⋮σenσe1σenσe2⋯σen2
\right]
P
k
=E(e
k
e
k
T
)=
⎣
⎢
⎢
⎢
⎡
σ
e1
2
σ
e2
σ
e1
⋮
σ
en
σ
e1
σ
e1
σ
e2
σ
e2
2
⋮
σ
en
σ
e2
⋯
⋯
⋱
⋯
σ
e1
σ
en
σ
e2
σ
en
⋮
σ
en
2
⎦
⎥
⎥
⎥
⎤
同理也有先验误差协方差矩阵 P k − P_{k}^{-} P
k
−
。根据最小方差估计原理，设损失函数为 t r ( P k ) \mathrm{tr}\left( P_k \right) tr(P
k
)，即要求
K k = a r g min [ t r ( P k ) ] K_k=\mathrm{arg}\min \left[ \mathrm{tr}\left( P_k \right) \right]
K
k
=argmin[tr(P
k
)]
将 P k P_k P
k
展开为：
上面运用了状态更新方程与状态预测方程，考虑到 e k − e_{k}^{-} e
k
−
与 v k v_k v
k
相互独立，则进一步：
P
k
=(I−K
k
C)E(e
k
−
e
k
−
T
)(I−K
k
C)
T
+K
k
E(v
k
v
k
T
)K
k
T
=(I−K
k
C)P
k
−
(I−K
k
C)
T
+K
k
RK
k
T
令 ∂ t r ( P k ) ∂ K k = 0 \frac{\partial \mathrm{tr}\left( P_k \right)}{\partial K_k}=0
∂K
k
∂tr(P
k
)
=0，即
卡尔曼增益调整方程 K k = P k − C T C P k − C T + R {\text{卡尔曼增益调整方程}K_k=\frac{P_{k}^{-}C^T}{CP_{k}^{-}C^T+R}}
卡尔曼增益调整方程K
k
=
CP
k
−
C
T
+R
P
k
−
C
T
将 K k K_k K
k
代入后验误差协方差矩阵表达式，即得
协方差更新方程 P k = ( I − K k H ) P k − {\text{协方差更新方程}P_k=\left( I-K_kH \right) P_{k}^{-}}
协方差更新方程P
k
=(I−K
k
H)P
k
−
要更新 K k K_k K
k
，则只需要确定先验误差协方差矩阵 P k − P_{k}^{-} P
k
−
。同样地，将 P k − P_{k}^{-} P
k
−
展开为
P k − = E ( e k − e k − T ) = E [ ( x k − x ^ k − ) ( x k − x ^ k − ) T ] = E [ ( A e k − 1 + w k − 1 ) ( A e k − 1 + w k − 1 ) T ] = A E ( e k − 1 e k − 1 T ) A T + E ( w k − 1 w k − 1 T ) P_{k}^{-}=E\left( e_{k}^{-}e_{k}^{-T} \right) \\=E\left[ \left( x_k-\hat{x}_{k}^{-} \right) \left( x_k-\hat{x}_{k}^{-} \right) ^T \right] \\=E\left[ \left( Ae_{k-1}+w_{k-1} \right) \left( Ae_{k-1}+w_{k-1} \right) ^T \right] \\=AE\left( e_{k-1}e_{k-1}^{T} \right) A^T+E\left( w_{k-1}w_{k-1}^{T} \right)
P
k
−
=E(e
k
−
e
k
−T
)
=E[(x
k
−
x
^
k
−
)(x
k
−
x
^
k
−
)
T
]
=E[(Ae
k−1
+w
k−1
)(Ae
k−1
+w
k−1
)
T
]
=AE(e
k−1
e
k−1
T
)A
T
+E(w
k−1
w
k−1
T
)
考虑到 e k e_k e
k
与 w k w_k w
k
相互独立，进一步得到
协方差预测方程 P k − = A P k − 1 A T + Q {\text{协方差预测方程}P_{k}^{-}=AP_{k-1}A^T+Q}
协方差预测方程P
k
−
=AP
k−1
A
T
+Q
至此，得到卡尔曼滤波的五大基本公式，其中核心方程为基于最小方差估计的卡尔曼增益调整方程，具体工作流程如图所示。
4 具体案例：船舶GPS定位
有一船舶出港沿某直线方向航行，辅助北斗卫星进行定位和测速。假设
① 船舶加速度 a ( t ) a\left( t \right) a(t)=机动加速度 u ( t ) u\left( t \right) u(t)+随机加速度 w ( t ) w\left( t \right) w(t)，其中 w ( t ) w\left( t \right) w(t)符合高斯分布；
② GPS观测噪声 v ( t ) v\left( t \right) v(t)符合高斯分布。
要求用卡尔曼滤波器估计真实运动轨迹。
解决方案
首先建立随机系统的真实模型。以码头出发点为原点，采样周期(雷达扫描周期)为 T T T，用 s ( k ) s(k) s(k)表示船舶在采样时刻 k T kT kT处的真实位置，用 z ( k ) z\left( k \right) z(k)表示在时刻 k T kT kT处的GPS定位观测值，则真实系统输出方程为 z ( k ) = s ( k ) + v ( k ) z(k)=s(k)+v(k) z(k)=s(k)+v(k)。
记在 k T kT kT时刻处船舶速 s ˙ ( k ) \dot{s}\left( k \right)
s
˙
(k)，加速度为 a ( k ) a\left( k \right) a(k)，由匀加速公式有：
s ( k + 1 ) = s ( k ) + s ˙ ( k ) T + 1 2 a ( k ) T 2 ⇒ s ˙ ( k + 1 ) = s ˙ ( k ) + T a ( k ) s\left( k+1 \right) =s\left( k \right) +ṡ\left( k \right) T+\frac{1}{2}a\left( k \right) T^2\\\Rightarrow \dot{s}\left( k+1 \right) =ṡ\left( k \right) +Ta\left( k \right)
s(k+1)=s(k)+
s
˙
(k)T+
2
1
a(k)T
2
⇒
s
˙
(k+1)=
s
˙
(k)+Ta(k)
其中 a ( k ) = u ( k ) + w ( k ) a(k)=u(k)+w(k) a(k)=u(k)+w(k)。
定义在采样时刻 k T kT kT系统状态 x ( k ) x\left( k \right) x(k)为船舶的位置和速度，即 x ( k ) = [ s ( k ) s ˙ ( k ) ] T x(k)=\left[
s(k)s˙(k)
s(k)s˙(k)
\right] ^T x(k)=[
s(k)
s
˙
(k)
]
T
可得到船舶运动的状态空间模型
{ [ s ( k + 1 ) s ˙ ( k + 1 ) ] = [ 1 T 0 1 ] [ s ( k ) s ˙ ( k ) ] + [ 0.5 T T ] u ( k ) + [ 0.5 T T ] w ( k ) z ( k ) = [ 1 0 ] [ s ( k ) s ˙ ( k ) ] + v ( k )
⎧⎩⎨⎪⎪⎪⎪⎪⎪[s(k+1)s˙(k+1)]=[10T1][s(k)s˙(k)]+[0.5TT]u(k)+[0.5TT]w(k)z(k)=[10][s(k)s˙(k)]+v(k)
{[s(k+1)s˙(k+1)]=[1T01][s(k)s˙(k)]+[0.5TT]u(k)+[0.5TT]w(k)z(k)=[10][s(k)s˙(k)]+v(k)
⎩
⎪
⎪
⎨
⎪
⎪
⎧
[
s(k+1)
s
˙
(k+1)
]=[
1
0
T
1
][
s(k)
s
˙
(k)
]+[
0.5T
T
]u(k)+[
0.5T
T
]w(k)
z(k)=[
1
0
][
s(k)
s
˙
(k)
]+v(k)
在不考虑机动目标的动力因素即 u ( k ) = 0 u(k)=0 u(k)=0时，将匀速直线运动的船舶系统扩展到四维，状态包含水平和纵向的位置和速度，则系统方程可化为
[ x ( k ) x ˙ ( k ) y ( k ) y ˙ ( k ) ] = [ 1 T 0 0 0 1 0 0 0 0 1 T 0 0 0 1 ] [ x ( k − 1 ) x ˙ ( k − 1 ) y ( k − 1 ) y ˙ ( k − 1 ) ] + [ 0.5 T 2 T 0.5 T 2 T ] [ w 1 ( k ) w 2 ( k ) w 3 ( k ) w 4 ( k ) ] z ( k ) = [ 1 0 0 0 0 0 1 0 ] [ x ( k ) x ˙ ( k ) y ( k ) y ˙ ( k ) ] + [ v 1 ( k ) v 2 ( k ) ] \left[
x(k)x˙(k)y(k)y˙(k)
x(k)x˙(k)y(k)y˙(k)
\right] =\left[
1000T100001000T1
1T000100001T0001
\right] \left[
x(k−1)x˙(k−1)y(k−1)y˙(k−1)
x(k−1)x˙(k−1)y(k−1)y˙(k−1)
\right] +\left[
0.5T2T0.5T2T
0.5T2T0.5T2T
\right] \left[
w1(k)w2(k)w3(k)w4(k)
w1(k)w2(k)w3(k)w4(k)
\right] \\z(k)=\left[
10000100
10000010
\right] \left[
x(k)x˙(k)y(k)y˙(k)
x(k)x˙(k)y(k)y˙(k)
\right] +\left[
v1(k)v2(k)
v1(k)v2(k)
\right]
⎣
⎢
⎢
⎡
x(k)
x
˙
(k)
y(k)
y
˙
(k)
⎦
⎥
⎥
⎤
=
⎣
⎢
⎢
⎡
1
0
0
0
T
1
0
0
0
0
1
0
0
0
T
1
⎦
⎥
⎥
⎤
⎣
⎢
⎢
⎡
x(k−1)
x
˙
(k−1)
y(k−1)
y
˙
(k−1)
⎦
⎥
⎥
⎤
+
⎣
⎢
⎢
⎡
0.5T
2
T
0.5T
2
T
⎦
⎥
⎥
⎤
⎣
⎢
⎢
⎡
w
1
(k)
w
2
(k)
w
3
(k)
w
4
(k)
⎦
⎥
⎥
⎤
z(k)=[
1
0
0
0
0
1
0
0
]
⎣
⎢
⎢
⎡
x(k)
x
˙
(k)
y(k)
y
˙
(k)
⎦
⎥
⎥
⎤
+[
v
1
(k)
v
2
(k)
]
现假设初始位置为(-100m，200m)，水平运动初速度为2m/s，垂直运动初速度为20m/s，雷达扫描周期 T = 1 s T=1s T=1s，则系统状态空间进一步化为
{ x ( k + 1 ) = A x ( k ) + B u ( k ) + Γ w ( k ) z ( k ) = H x ( k ) + v ( k )
{x(k+1)=Ax(k)+Bu(k)+Γw(k)z(k)=Hx(k)+v(k)
{x(k+1)=Ax(k)+Bu(k)+Γw(k)z(k)=Hx(k)+v(k)
{
x(k+1)=Ax(k)+Bu(k)+Γw(k)
z(k)=Hx(k)+v(k)
其中 A = [ 1 1 0 0 0 1 0 0 0 0 1 1 0 0 0 1 ] A=\left[
1000110000100011
1100010000110001
\right] A=
⎣
⎢
⎢
⎡
1
0
0
0
1
1
0
0
0
0
1
0
0
0
1
1
⎦
⎥
⎥
⎤
， B = Γ = [ 0.5 1 0.5 1 ] B=\varGamma =\left[
0.510.51
0.510.51
\right] B=Γ=
⎣
⎢
⎢
⎡
0.5
1
0.5
1
⎦
⎥
⎥
⎤
, H = [ 1 0 0 0 0 0 1 0 ] H=\left[
10000100
10000010
\right] H=[
1
0
0
0
0
1
0
0
] 。
设过程噪声 w ∼ N ( 0 , 0. 1 2 ) w\sim N(0,0.1^2) w∼N(0,0.1
2
)，其四个分量独立同分布，则
Q ∗ = E ( w 4 × 1 w 4 × 1 T ) = σ w 2 [ 1 1 1 1 ] Q^*=E(w_{4\times 1}w_{4\times 1}^{T})=\sigma _{w}^{2}\left[
1111
1111
\right]
Q
∗
=E(w
4×1
w
4×1
T
)=σ
w
2
⎣
⎢
⎢
⎡
1
1
1
1
⎦
⎥
⎥
⎤
考虑噪声矩阵 Γ \varGamma Γ的加权得到过程噪声协方差矩阵
Q = Q ∗ ⋅ Γ = σ w 2 [ 0.5 1 0.5 1 ] Q=Q^*\cdot \varGamma =\sigma _{w}^{2}\left[
0.510.51
0.510.51
\right]
Q=Q
∗
⋅Γ=σ
w
2
⎣
⎢
⎢
⎡
0.5
1
0.5
1
⎦
⎥
⎥
⎤
观测噪声 v ∼ N ( 0 , 1 0 2 ) v\sim N(0,10^2) v∼N(0,10
2
)，其两个分量独立同分布，则观测噪声协方差矩阵
R = σ v 2 [ 1 1 ] R=\sigma _{v}^{2}\left[
11
11
\right]
R=σ
v
2
[
1
1
]
依此建立卡尔曼滤波器模型。
设置最优估计初值 x ^ ( 1 ) = [ − 100 2 − 200 20 ] T \hat{x}(1)=\left[
−1002−20020
−1002−20020
\right] ^T
x
^
(1)=[
−100
2
−200
20
]
T
、先验协方差矩阵 P 1 − = [ 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 ] P_{1}^{-}=\left[
1000010000100001
1000010000100001
\right] P
1
−
=
⎣
⎢
⎢
⎡
1
0
0
0
0
1
0
0
0
0
1
0
0
0
0
1
⎦
⎥
⎥
⎤
，并以此开始迭代，得到实验图像如图所示。
在使用卡尔曼滤波器时，应至少保证系统模型或系统测量至少一个有足够的精度，否则卡尔曼滤波器无法从中提取到正确的估计信息。这是因为卡尔曼滤波调整的是对先验模型值与仪器测量值间的信任权重，若二者都不准确则卡尔曼跟踪的数值亦不准确。为验证这一结论，下面保持测量噪声不变，调整过程噪声方差，如图所示。
完整代码关注下方公众号回复“EE001”获取。
🔥 更多精彩专栏：
