最小二乘估计（3）

1、最小二乘的性能——估计方差

定量描述估计的误差，估计方差是一个很关键的评价指标。本节给出几种最小二乘方法滤波增益和估计方差的求解方法，这些也是卡尔曼滤波器的理论基础。

估计方差的定义为
$$P(k)=E\left[(\theta -\hat{\theta }(k))(\theta -\hat{\theta }(k){)}^T\right]$$

式中，$\theta$ 为待估计参数的真值，$\hat{\theta }(k)$ 为第 $k$ 步的估计值，$E(\cdot )$ 表示求均值。一般可以认为，所得的估计方差越小，估计方法就越好。

接下来以递推最小二乘方法为例，研究一下最小二乘方法的估计方差。前面已经讲过,权值 $W(k)$ 一般取为测量噪声方差 $R(k)$ 的倒数，如果权值 $W(k)$ 为一个矩阵的话则取为 $R(k)$ 的逆。从现在开始都用 $R^{-1}(k)$ 来代替 $W(k)$。
$$\boxed{\begin{array}{rl}\hat{\theta }(k)& =\hat{\theta }(k-1)+M_k{H}^T(k)W(k)\left[Z(k)-H(k)\hat{\theta }(k-1)\right]\\ M_k^{-1}& =M_{k-1}^{-1}+H^T(k)W(k)H(k)\end{array}}\text{\hspace{1em}(1)}$$

我们令 $K(k)=M_k{H}^T(k)R^{-1}(k)$，称其为滤波增益，将上式代入估计方差的定义，可以得到
$$\begin{array}{rl}P(k)& =E\left[(\theta -\hat{\theta }(k))(\theta -\hat{\theta }(k){)}^T\right]\\ & =E\left[\right(\theta -\hat{\theta }(k-1)-K(k)\left(Z(k)-H(k)\hat{\theta }(k-1)\right))\\ & (\theta -\hat{\theta }(k-1)-K(k)\left(Z(k)-H(k)\hat{\theta }(k-1)\right){)}^T]\end{array}\text{\hspace{1em}(2)}$$

再将测量方程 $Z(k)=H(k)\theta +N(k)$ 代入上式，考虑到待估计值与 $N(k)$ 不相关，整理后可以得到
$$P(k)=[I-K(k)H(k)]P(k-1)[I-K(k)H(k){]}^T+K(k)R(k)K^T(k)\text{\hspace{1em}(3)}$$

接下来讨论递推最小二乘方法的估计性能问题，其结果是最优的吗？

令式（3）$P(k)$ 求偏导得到的子式等于 0，即
$$\frac{\partial P(k)}{\partial K(k)}=2[I-K(k)H(k)]P(k-1)[-H^T(k)]+2K(k)R(k)=0\text{\hspace{1em}(4)}$$

整理可得
$$K(k)=P(k-1)H^T(k)[R(k)+H(k)P(k-1)H^T(k){]}^{-1}\text{\hspace{1em}(5)}$$

这个答案初看起来令人惊讶，和已经得到的滤波增益 $K(k)=M_k{H}^T(k)R^{-1}(k)$ 是不同的，一个包含 $M_k$，另一个包含 $P(k-1)$。但下面2个结论（证明方法见附录（有空再补充吧））却表明这两个公式是等价的。

1. 若考虑 $W(k)=R^{-1}(k)$，则 $M_k$ 实际上是估计方差 $P(k)$，按照式（1），估计方差也满足 $P^{-1}(k)=P^{-1}(k-1)+H^T(k)R^{-1}(k)H(k)$。
2. 若用 $P(k)$ 来代替式（1）中的 $M_k$，则式（5）所表达的滤波增益是等价的。

到这里，可以得到如下结论，利用最小二乘性能指标得到的估计方法:
$$\hat{\theta }(k)=\hat{\theta }(k-1)+K(k)\left[Z(k)-H(k)\hat{\theta }(k-1)\right]\text{\hspace{1em}(6)}$$

可以得到最小的估计方差 $P(k)$，其中滤波增益 $K(k)$ 有以下两种不同的计算方法：
$$\begin{array}{rl}K(k)& =M_k{H}^T(k)R^{-1}(k)\\ K(k)=P(k-1)H^T& (k)[R(k)+H(k)P(k-1)H^T(k){]}^{-1}\end{array}\text{\hspace{1em}(7)}$$

估计方差 $P(k)$ 可以使用如下两种求法来计算：
$$P(k)=[P^{-1}(k-1)+H^T(k)R^{-1}(k)H(k){]}^{-1}\text{\hspace{1em}(8)}$$

$$P(k)=[I-K(k)H(k)]P(k-1)[I-K(k)H(k){]}^T+K(k)R(k)K^T(k)\text{\hspace{1em}(9)}$$

在附录中还提到了下面的估计方差 $P(k)$ 计算方法：
$$P(k)=[I-K(k)H(k)]P(k-1)\text{\hspace{1em}(10)}$$

在上面的三个计算 $P(k)$ 的公式中，式（10）是一个比较简单的表达式，但是数值计算问题可能导致 $P(k)$ 不是正定的，即使 $P(k-1)$ 和 $R(k)$ 都是正定的。式（9）虽然复杂，但是却可以保证 $P(k)$ 是正定的，与理论解一致，而式（8）需要三次矩阵求逆计算，因此计算复杂，需要占用更多的计算资源，但在多传感器状态融合估计中可以使融合公式具有统一形式，所以经常在状态融合估计方法中使用。