最小二乘估计

最小二乘估计法介绍

​ 最小二乘估计是一种利用观测数据估计线性模型中未知参数的方法，其基本的思想是选择合适的估计参数使得模型输出与传感器实测输出数据之差的平方和最小。

​ 对于一个线性模型，其含有 \(m+1\) 种可观测的变量 \((\Omega_0,\Omega_1,...,\Omega_m)\)，每个参数（除 \(\Omega_0\) 之外）具有一个未知的参数 \((O_1,O_2,...,O_m)\).

​ 模型方程为：

\[\Omega_0 + \sum_{j=1}^mO_j\Omega_j = 0\tag{1} \]

​ 现在需要进行 \(n\) 次观测来估计这些未知参数，要求 \(n\ge m\)。由于每次观测存在误差，所以对于每次观测，其方程应修改为：

\[\Omega_{i0} + \sum_{j=1}^mO_{ij}\Omega_{ij} = \varepsilon_i\tag{2} \]

式中，\(\varepsilon_i\) 即为每次的观测误差。

矩阵描述

​ 用 \(\theta\) 表示待估计量，\(Z(k)\)、\(N(k)\) 表示第 \(k\) 次观测数据和观测噪声，\(H(k)\) 为观测矩阵，则线性观测方程可描述为：

\[Z(k) = H(k)\theta + N(k)\tag{3} \]

使用矩阵描述可以大大简化后面的计算推导。

​ 用 \(\hat\theta\) 表示估计量，则最小二乘法要求测量模型输出 \(H(k)\hat\theta\) 与实际测量数据 \(Z(k)\) 的平方和最小，将其表示为与 \(\hat\theta\) 相关的一个函数：

\[J(\hat\theta) = (Z_k-H_k\hat\theta)^T(Z_k-H_k\hat\theta)\tag{4} \]

补充矩阵求导的相关公式

设 \(x\in F^{n\times n}\) 为常矩阵，有 \(f=x^TAx\)，则 \(\dfrac{\mathrm df}{\mathrm dx} = (A+A^T)x\)。

证明 根据乘法公式：

\[\frac{\mathrm df}{\mathrm dx} = \frac{\mathrm d}{\mathrm dx}(x^TAx) = \frac{\mathrm dx^T}{\mathrm dx}Ax+\frac{\mathrm d(Ax)^T}{\mathrm dx}x = Ax+A^Tx=(A+A^T)x \]

​ 对 \(J(\hat\theta)\) 求导得：

\[\frac{\part}{\part\hat\theta}[(Z_k-H_k\hat\theta)^T(Z_k-H_k\hat\theta)] = -2H_k^T(Z_k-H_k\hat\theta)\tag{5} \]

令导数等于0得

\[H_k^TZ_k = H_k^TH_k\hat\theta\tag{6} \]

若 \(H_k^TH_k\) 为满秩矩阵，则

\[\hat\theta = (H_k^TH_k)^{-1}H_k^TZ_k\tag{7} \]

最小二乘加权估计

​ 在上面的最小二乘估计中并没有使用测量设备的噪声方差，如果在已知噪声方差的情况下，理论上应该可以减小估计结果的偏差。方差更大的数据在估计中的作用应该越小，也就是对测量数据进行加权处理。

​ 假设测量噪声为 \(N_k\)，并有 \(E(N_k) = 0, E(N_kN_k^T)=R_k\)，取权值矩阵 \(W=R_k^{-1}\)。则函数 \(J(\hat\theta)\) 修改为：

\[J_w(\hat\theta) = (Z_k-H_k\hat\theta)^TW(Z_k-H_k\hat\theta)\tag{8} \]

求极值得：

\[\hat\theta = (H_k^TWH_k)^{-1}H_k^TWZ_k\tag{9} \]

线性最小二乘递推估计

​ 之前讲到的最小二乘法都需要等到收集到所有的测量数据之后才能进行估计，不具有实时性。递推估计的核心思想是，在获得测量数据之后及时进行处理，在估计当前时刻的参数时，利用前一次得到的结果。

​ 根据式9, 在 \(k-1\) 时刻的估计结果为：

\[\hat\theta(k-1) = [H_{k-1}^TW_{k-1}H_{k-1}]^{-1}H_{k-1}^TW_{k-1}Z_{k-1}\tag{10} \]

令 \(M_k = [H_k^TW_kH_k]^{-1}\)，则由 \(k-1\) 时刻相关数据递推 \(k\) 时刻的公式为：

\[M_k = \left\{\begin{bmatrix}H_{k-1}^T&H^T(k)\end{bmatrix}\begin{bmatrix}W_{k-1}&0\\0&W(k)\end{bmatrix}\begin{bmatrix}H_{k-1}\\H(k)\end{bmatrix}\right\}^{-1}\tag{11} \]

注：这里分清 \(H_k\) 和 \(H(k)\) 的区别，\(H_k\) 表示 0 到 \(k\) 时刻的数据所组成的向量，而 \(H(k)\) 单指 \(k\) 时刻，其他变量同理。

​ 进而可以计算得到

\[M_k = [M_{k-1}^{-1}+H^T(k)W(k)H(k)]^{-1}\tag{12} \]

由式10 得

\[\begin{aligned} \hat\theta(k) = M_kH_k^TW_kZ_k &= M_k\begin{bmatrix}H_{k-1}^T&H^T(k) \end{bmatrix} \begin{bmatrix}W_{k-1}&0\\0&W(k) \end{bmatrix}\begin{bmatrix}Z_{k-1}\\Z(k) \end{bmatrix}\\ &=M_k[H_{k-1}^TW_{k-1}Z_{k-1}+H^T(k)W(k)Z(k)] \end{aligned}\tag{13} \]

又 \(H_{k-1}^TW_{k-1}Z_{k-1} = M_{k-1}^{-1}\hat\theta(k-1)\)，代入式13 再消去 \(M_{k-1}^{-1}\) 得

\[\hat\theta(k) = \hat\theta(k-1) + M_kH^T(k)[Z(k)-H(k)\hat\theta(k-1)]\tag{14} \]

​ 现将递推最小二乘估计总结如下：

​ 设 \(M_0\)、\(\hat\theta(0)\) 为迭代初始值，\(\hat\theta(0)\) 为一合适的数，\(M_0\) 为一个大的数或正定矩阵，其维度与 \(H^T(k)W(k)H(k)\) 相同，递推最小二乘法由两个递推公式组成：

\[M_k^{-1} = M_{k-1}^{-1} + H^T(k)W(k)H(k)\\ \hat\theta(k) = \hat\theta(k-1) + M_kH^T(k)W(k)[Z(k) - H(k)\hat\theta(k-1)]\tag{16} \]

式中，\(M_kH^T(k)W(k)\) 为滤波增益\(K(k)\)。

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

​ 估计方差的定义为

\[P(k) = E[(\theta-\hat\theta(k))(\theta-\hat\theta(k))^T]\tag{17} \]

代入 \(\hat\theta(k) = \hat\theta(k-1)+K(k)[Z(k)-H(k)\hat\theta(k-1)]\) 和式3 得

\[P(k) = [I-K(k)H(k)]P(k-1)[I-K(k)H(k)]^T+K(k)R(k)K^T(k)\tag{18} \]

\(P(k)\) 对滤波增益求偏导得

\[\frac{\part P(k)}{\part K(k)} = 2[I-K(k)H(k)]P(k-1)[-H^T(k)]+2K(k)R(k)= 0\tag{19} \]

整理得

\[K(k) = P(k-1)H^T(k)[R(k)+H(k)P(k-1)H^T(k)]^{-1}\tag{20} \]