最小二乘估计法介绍
最小二乘估计是一种利用观测数据估计线性模型中未知参数的方法，其基本的思想是选择合适的估计参数使得模型输出与传感器实测输出数据之差的平方和最小。
对于一个线性模型，其含有 \(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)]\)
\[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} \]
我愿潇洒如鹰，远离地上宿命
