1、最小二乘加权估计
最小二乘估计(1)中并没有使用测量噪声方差,加权的依据来源于测量的准确性。假设测量噪声为高斯白噪声
N
k
\pmb N_k
Nk,其均值和协方差分别为
E
(
N
k
)
=
0
,
E
(
N
k
N
k
T
)
=
R
k
E(\pmb N_k)=0,\ E(\pmb N_k\pmb N_k^T)=\pmb R_k
E(Nk)=0, E(NkNkT)=Rk,若每次测量噪声不相关,则
R
k
\pmb R_k
Rk 为对角阵,即
R
k
=
[
R
(
1
)
R
(
2
)
⋱
R
(
k
)
]
(1)
\pmb R_k=\left[ \begin{matrix} R(1) & & & \\ & R(2) & & \\ & & \ddots & \\ & & & R(k)\\ \end{matrix} \right]\tag{1}
Rk=
R(1)R(2)⋱R(k)
(1)
式中, R ( 1 ) , R ( 2 ) , ⋯ , R ( k ) R(1),R(2),\cdots,R(k) R(1),R(2),⋯,R(k) 是测量噪声在 1 , 2 , ⋯ , k 1,2,\cdots,k 1,2,⋯,k 时刻的方差。设权值用 W \pmb W W 表示,可取 W = R k − 1 \pmb W = \pmb R_k^{-1} W=Rk−1。这样做的原因在于当数据测量方差大时,其权重应设置小些。当然,也可以使用其他函数实现这种关系,但是一般情况下都将测量噪声方差的“逆”取为权值。
考虑每个测量时刻的测量噪声方差,设线性最小二乘加权估计的性能指标为
J
w
(
θ
^
)
=
(
Z
k
−
H
k
θ
^
)
W
(
Z
k
−
H
k
θ
^
)
(2)
\pmb J_w(\hat{\pmb\theta}) = (\pmb Z_k - \pmb H_k\hat{\pmb\theta})\pmb W(\pmb Z_k - \pmb H_k\hat{\pmb\theta})\tag{2}
Jw(θ^)=(Zk−Hkθ^)W(Zk−Hkθ^)(2)
接下来的任务是选择
θ
^
\hat{\pmb\theta}
θ^ 使之达到最小。与 最小二乘估计(1)类似,使用求极值的方式来解决问题。令
∂
J
w
(
θ
^
)
∂
θ
^
=
0
\dfrac{\partial\pmb J_w(\hat{\pmb\theta})}{\partial\hat{\pmb\theta}}=0
∂θ^∂Jw(θ^)=0,则
∂
J
w
(
θ
^
)
∂
θ
^
=
∂
∂
θ
^
[
(
Z
k
−
H
k
θ
^
)
W
(
Z
k
−
H
k
θ
^
)
]
=
−
2
H
k
T
W
(
Z
k
−
H
k
θ
^
)
=
0
(3)
\dfrac{\partial\pmb J_w(\hat{\pmb\theta})}{\partial\hat{\pmb\theta}}=\dfrac{\partial}{\partial\hat{\pmb\theta}}[(\pmb Z_k - \pmb H_k\hat{\pmb\theta})\pmb W(\pmb Z_k - \pmb H_k\hat{\pmb\theta})]=-2\pmb H_k^T\pmb W(\pmb Z_k-\pmb H_k\hat{\pmb\theta})=0\tag{3}
∂θ^∂Jw(θ^)=∂θ^∂[(Zk−Hkθ^)W(Zk−Hkθ^)]=−2HkTW(Zk−Hkθ^)=0(3)
解上述方程得到
θ
^
=
(
H
k
T
W
H
k
)
−
1
H
k
T
W
Z
k
(4)
\hat{\pmb\theta}=(\pmb H_k^T\pmb W\pmb H_k)^{-1}\pmb H_k^T\pmb W\pmb Z_k\tag{4}
θ^=(HkTWHk)−1HkTWZk(4)
2、线性最小二乘递推估计
递推估计的核心思想是,在获得测量数据以后及时地进行处理,而不必等到所有测量数据都获得后再进行估计,而且在估计当前时刻的参数时,可以利用前一次得到的结果,而不必将数据全部重新计算一遍。也就是说,线性最小二乘递推估计方法关注的是:如果已经得到第 k − 1 k-1 k−1 步的估计 θ ^ ( k − 1 ) \hat{\pmb\theta}(k-1) θ^(k−1),当第 k k k 步的测量 Z ( k ) Z(k) Z(k) 到达后,在 H ( k ) \pmb H(k) H(k) 已知情况下,如何利用 θ ^ ( k − 1 ) 、 Z ( k ) \hat{\pmb\theta}(k-1)、Z(k) θ^(k−1)、Z(k) 和 H ( k ) \pmb H(k) H(k) 构造估计量 θ ^ ( k ) \hat{\pmb\theta}(k) θ^(k) 使估计结果越来越接近估计量的真实值。
设加权矩阵为
W
k
−
1
=
d
i
a
g
[
W
(
1
)
,
W
(
2
)
,
⋯
,
W
(
k
−
1
)
]
\pmb W_{k-1}=\rm{diag}[\pmb W(1),\pmb W(2),\cdots,\pmb W(k-1)]
Wk−1=diag[W(1),W(2),⋯,W(k−1)],根据线性最小二乘加权估计方法式(4),可知估计矢量的计算方法
θ
^
(
k
−
1
)
\hat{\pmb\theta}(k-1)
θ^(k−1) 为
θ
^
(
k
−
1
)
=
(
H
k
−
1
T
W
k
−
1
H
k
−
1
)
−
1
H
k
−
1
T
W
k
−
1
Z
k
−
1
(5)
\hat{\pmb\theta}(k-1)=(\pmb H_{k-1}^T\pmb W_{k-1}\pmb H_{k-1})^{-1}\pmb H_{k-1}^T\pmb W_{k-1}\pmb Z_{k-1}\tag{5}
θ^(k−1)=(Hk−1TWk−1Hk−1)−1Hk−1TWk−1Zk−1(5)
式中,测量数据
Z
k
−
1
\pmb Z_{k-1}
Zk−1、测量矩阵
H
k
−
1
\pmb H_{k-1}
Hk−1 以及相应的测量噪声
N
k
−
1
\pmb N_{k-1}
Nk−1 分别为
Z
k
−
1
=
[
z
(
1
)
z
(
2
)
⋮
z
(
k
−
1
)
]
,
H
k
−
1
=
[
H
(
1
)
H
(
2
)
⋮
H
(
k
−
1
)
]
,
N
k
−
1
=
[
N
(
1
)
N
(
2
)
⋮
N
(
k
−
1
)
]
(6)
\pmb Z_{k-1}=\left[ \begin{matrix} \ z(1)\\ \ z(2)\\ \vdots\\ \ z(k-1)\\ \end{matrix} \right],\quad \pmb H_{k-1}=\left[ \begin{matrix} \ \pmb H(1)\\ \ \pmb H(2)\\ \vdots\\ \ \pmb H(k-1)\\ \end{matrix} \right],\quad \pmb N_{k-1}=\left[ \begin{matrix} \ N(1)\\ \ N(2)\\ \vdots\\ \ N(k-1)\\ \end{matrix} \right]\tag{6}
Zk−1=
z(1) z(2)⋮ z(k−1)
,Hk−1=
H(1) H(2)⋮ H(k−1)
,Nk−1=
N(1) N(2)⋮ N(k−1)
(6)
令
M
k
−
1
=
(
H
k
−
1
T
W
k
−
1
H
k
−
1
)
−
1
\pmb M_{k-1}=(\pmb H_{k-1}^T\pmb W_{k-1}\pmb H_{k-1})^{-1}
Mk−1=(Hk−1TWk−1Hk−1)−1,则
θ
^
(
k
−
1
)
=
M
k
−
1
H
k
−
1
T
W
k
−
1
Z
k
−
1
(7)
\hat{\pmb\theta}(k-1)=\pmb M_{k-1}\pmb H_{k-1}^T\pmb W_{k-1}\pmb Z_{k-1}\tag{7}
θ^(k−1)=Mk−1Hk−1TWk−1Zk−1(7)
在
H
(
k
)
\pmb H(k)
H(k) 已知的情况下,当得到第
k
−
1
k-1
k−1 步的估计
θ
^
(
k
−
1
)
\hat{\pmb\theta}(k-1)
θ^(k−1),并且测量
Z
(
k
)
Z(k)
Z(k) 到达后,根据最小二乘加权方法式(4),可以得到
θ
^
(
k
)
\hat{\pmb\theta}(k)
θ^(k):
θ
^
(
k
)
=
(
H
k
T
W
k
H
k
)
−
1
H
k
T
W
k
Z
k
(8)
\hat{\pmb\theta}(k)=(\pmb H_{k}^T\pmb W_{k}\pmb H_{k})^{-1}\pmb H_{k}^T\pmb W_{k}\pmb Z_{k}\tag{8}
θ^(k)=(HkTWkHk)−1HkTWkZk(8)
前面提过,测量数据
Z
k
\pmb Z_{k}
Zk、测量矩阵
H
k
\pmb H_{k}
Hk 以及相应的测量噪声
N
k
\pmb N_{k}
Nk 分别为
Z
k
=
[
z
(
1
)
z
(
2
)
⋮
z
(
k
)
]
,
H
k
−
1
=
[
H
(
1
)
H
(
2
)
⋮
H
(
k
)
]
,
N
k
−
1
=
[
N
(
1
)
N
(
2
)
⋮
N
(
k
)
]
(9)
\pmb Z_{k}=\left[ \begin{matrix} \ z(1)\\ \ z(2)\\ \vdots\\ \ z(k)\\ \end{matrix} \right],\quad \pmb H_{k-1}=\left[ \begin{matrix} \ \pmb H(1)\\ \ \pmb H(2)\\ \vdots\\ \ \pmb H(k)\\ \end{matrix} \right],\quad \pmb N_{k-1}=\left[ \begin{matrix} \ N(1)\\ \ N(2)\\ \vdots\\ \ N(k)\\ \end{matrix} \right]\tag{9}
Zk=
z(1) z(2)⋮ z(k)
,Hk−1=
H(1) H(2)⋮ H(k)
,Nk−1=
N(1) N(2)⋮ N(k)
(9)
同样令
M
k
=
(
H
k
T
W
k
H
k
)
−
1
\pmb M_{k}=(\pmb H_{k}^T\pmb W_{k}\pmb H_{k})^{-1}
Mk=(HkTWkHk)−1,则
θ
^
(
k
)
=
M
k
H
k
T
W
k
Z
k
(10)
\hat{\pmb\theta}(k)=\pmb M_{k}\pmb H_{k}^T\pmb W_{k}\pmb Z_{k}\tag{10}
θ^(k)=MkHkTWkZk(10)
现在需要建立
θ
^
(
k
)
\hat{\pmb\theta}(k)
θ^(k) 和
θ
^
(
k
−
1
)
\hat{\pmb\theta}(k-1)
θ^(k−1) 之间的关系,首先
Z
k
=
[
Z
k
−
1
Z
(
k
)
]
\pmb Z_k=\left[ \begin{matrix} \ \pmb Z_{k-1}\\[1ex] \ Z(k)\\ \end{matrix} \right]
Zk=[ Zk−1 Z(k)],同样,
H
k
=
[
H
k
−
1
H
(
k
)
]
\pmb H_k=\left[ \begin{matrix} \ \pmb H_{k-1}\\[1ex] \ H(k)\\ \end{matrix} \right]
Hk=[ Hk−1 H(k)],对于测量噪声来说,噪声的关系矩阵
N
k
=
[
N
k
−
1
N
(
k
)
]
\pmb N_k=\left[ \begin{matrix} \ \pmb N_{k-1}\\[1ex] \ N(k)\\ \end{matrix} \right]
Nk=[ Nk−1 N(k)] 意味着噪声方差的关系满足
W
k
=
[
W
k
−
1
0
0
W
(
k
)
]
\pmb W_k=\left[ \begin{matrix} \ \pmb W_{k-1} & 0\\[1ex] \ 0 & \pmb W(k)\\ \end{matrix} \right]
Wk=[ Wk−1 00W(k)],很显然,
M
k
\pmb M_k
Mk 与
M
k
−
1
\pmb M_{k-1}
Mk−1 也是有关系的,将
H
k
\pmb H_k
Hk 和
W
k
\pmb W_k
Wk 代入
M
k
\pmb M_k
Mk 得
M
k
=
{
[
H
k
−
1
T
H
T
(
k
)
]
[
W
k
−
1
0
0
W
(
k
)
]
[
H
k
−
1
H
(
k
)
]
}
−
1
(11)
\pmb M_k=\left\{ \left[\ \pmb H_{k-1}^T\quad \pmb H^T(k)\ \right]\left[ \begin{matrix} \ \pmb W_{k-1} & 0\\[1ex] \ 0 & \pmb W(k)\\ \end{matrix} \right]\left[ \begin{matrix} \ \pmb H_{k-1}\\[1ex] \ H(k)\\ \end{matrix} \right]\right\}^{-1}\tag{11}
Mk={[ Hk−1THT(k) ][ Wk−1 00W(k)][ Hk−1 H(k)]}−1(11)
进行矩阵相乘运算得到
M
k
=
[
H
k
−
1
T
W
k
−
1
H
k
−
1
+
H
T
(
k
)
W
(
k
)
H
(
k
)
]
−
1
=
[
M
k
−
1
−
1
+
H
T
(
k
)
W
(
k
)
H
(
k
)
]
−
1
(12)
\begin{aligned} \pmb M_k&=[\ \pmb H_{k-1}^T\pmb W_{k-1}\pmb H_{k-1}+\pmb H^T(k)\pmb W(k)\pmb H(k)\ ]^{-1}\\[1ex] &=[\ \pmb M_{k-1}^{-1}+\pmb H^T(k)\pmb W(k)\pmb H(k)\ ]^{-1} \end{aligned}\tag{12}
Mk=[ Hk−1TWk−1Hk−1+HT(k)W(k)H(k) ]−1=[ Mk−1−1+HT(k)W(k)H(k) ]−1(12)
两边取逆后得到
M
k
−
1
\pmb M_{k-1}
Mk−1 与
M
k
\pmb M_{k}
Mk 的关系如下:
M
k
−
1
−
1
=
M
k
−
1
−
H
T
(
k
)
W
(
k
)
H
(
k
)
(13)
\boxed{\pmb M_{k-1}^{-1}=\pmb M_{k}^{-1} - \pmb H^T(k)\pmb W(k)\pmb H(k)\tag{13}}
Mk−1−1=Mk−1−HT(k)W(k)H(k)(13)
接下来我们将
Z
k
、
H
k
、
W
k
、
M
k
\pmb Z_k、\pmb H_k、\pmb W_k、\pmb M_k
Zk、Hk、Wk、Mk 代入式(10),有
θ
^
(
k
)
=
M
k
H
k
T
W
k
Z
k
=
M
k
[
H
k
−
1
T
H
T
(
k
)
]
[
W
k
−
1
0
0
W
(
k
)
]
[
Z
k
−
1
Z
(
k
)
]
=
M
k
[
H
k
−
1
T
W
k
−
1
Z
k
−
1
+
H
T
(
k
)
W
(
k
)
Z
(
k
)
]
=
M
k
[
M
k
−
1
−
1
θ
^
(
k
−
1
)
+
H
T
(
k
)
W
(
k
)
Z
(
k
)
]
=
M
k
[
(
M
k
−
1
−
H
T
(
k
)
W
(
k
)
H
(
k
)
)
θ
^
(
k
−
1
)
+
H
T
(
k
)
W
(
k
)
Z
(
k
)
]
=
θ
^
(
k
−
1
)
−
M
k
H
T
(
k
)
W
(
k
)
H
(
k
)
θ
^
(
k
−
1
)
+
M
k
H
T
(
k
)
W
(
k
)
Z
(
k
)
=
θ
^
(
k
−
1
)
+
M
k
H
T
(
k
)
W
(
k
)
[
Z
(
k
)
−
H
(
k
)
θ
^
(
k
−
1
)
]
(14)
\boxed{\begin{aligned} \hat{\pmb\theta}(k)&=\pmb M_{k}\pmb H_{k}^T\pmb W_{k}\pmb Z_{k}\\[1ex] &=\pmb M_{k}\left[\ \pmb H_{k-1}^T\quad \pmb H^T(k)\ \right]\left[ \begin{matrix} \ \pmb W_{k-1} & 0\\[1ex] \ 0 & \pmb W(k)\\ \end{matrix} \right]\left[ \begin{matrix} \ \pmb Z_{k-1}\\[1ex] \ Z(k)\\ \end{matrix} \right]\\[1ex] &=\pmb M_k\left[\ \pmb H_{k-1}^T\pmb W_{k-1}\pmb Z_{k-1}+\pmb H^T(k)\pmb W(k)Z(k)\ \right]\\[1ex] &=\pmb M_k\left[\ \pmb M_{k-1}^{-1}\hat{\pmb\theta}(k-1)+\pmb H^T(k)\pmb W(k)Z(k)\ \right]\\[1ex] &=\pmb M_k\left[\ \left(\pmb M_{k}^{-1} - \pmb H^T(k)\pmb W(k)\pmb H(k)\right)\hat{\pmb\theta}(k-1)+\pmb H^T(k)\pmb W(k)Z(k)\ \right]\\[1ex] &=\hat{\pmb\theta}(k-1)-\pmb M_k\pmb H^T(k)\pmb W(k)\pmb H(k)\hat{\pmb\theta}(k-1) + \pmb M_k\pmb H^T(k)\pmb W(k)Z(k)\\[1ex] &=\hat{\pmb\theta}(k-1)+\pmb M_k\pmb H^T(k)\pmb W(k)\left[Z(k)-\pmb H(k)\hat{\pmb\theta}(k-1) \right] \end{aligned}\tag{14}}
θ^(k)=MkHkTWkZk=Mk[ Hk−1THT(k) ][ Wk−1 00W(k)][ Zk−1 Z(k)]=Mk[ Hk−1TWk−1Zk−1+HT(k)W(k)Z(k) ]=Mk[ Mk−1−1θ^(k−1)+HT(k)W(k)Z(k) ]=Mk[ (Mk−1−HT(k)W(k)H(k))θ^(k−1)+HT(k)W(k)Z(k) ]=θ^(k−1)−MkHT(k)W(k)H(k)θ^(k−1)+MkHT(k)W(k)Z(k)=θ^(k−1)+MkHT(k)W(k)[Z(k)−H(k)θ^(k−1)](14)