最小二乘估计方法是一种利用观测数据估计线性模型中未知参数的方法，其基本思路是，选择合适的估计参数使模型输出与传感器实测输出数据之差的平方和最小。
在介绍最小二乘方法之前先明确一个现象，即传感器的测量数据和实际的真实值通常不是完全一致的。测量时，由于各种因素的影响，会对测量造成误差，导致物理量的测量数据和真实值不一样，这就是通常所说的“测得不准”。
1、最小二乘估计方法描述
用 θ\pmb\thetaθ 表示待估计量，Z(k)、N(k)Z(k)、N(k)Z(k)、N(k) 分别表示第 kkk 次观测数据和观测噪声，H(k)\pmb H(k)H(k) 为观测矩阵，线性观测方程可以用如下方程描述：
Z(k)=H(k)θ+N(k)(1)
Z(k) = \pmb H(k)\pmb\theta + N(k)\tag{1}
Z(k)=H(k)θ+N(k)(1)
假设某室内温度数据（1000个室温数据）如图所示，且真实值变化规律为 y=θ1x+θ2y =\theta_1x+\theta_2y=θ1x+θ2，求 θ1、θ2\theta_1、\theta_2θ1、θ2 的最优估计。
我们将观测方程设为 Z(k)=θ1t(k)+θ2+N(k)Z(k)=\theta_1t(k) + \theta_2 + N(k)Z(k)=θ1t(k)+θ2+N(k)，即
Z(k)=[ t(k)1 ][ θ1 θ2]+N(k)(2)
Z(k) = \textbf{[}\ t(k)\quad1\ \textbf{]}
\left[
\begin{matrix}
\ \theta_1 \\[1ex]
\ \theta_2
\end{matrix}
\right]
+N(k)\tag{2}
Z(k)=[ t(k)1 ][ θ1 θ2]+N(k)(2)
将 kkk 个观测数据写成向量的形式：
Zk=[ Z(1) Z(2)⋮ Z(k)],Hk=[ H(1) H(2)⋮ H(k)],Nk=[ N(1) N(2)⋮ N(k)](3)
\pmb Z_k = \left[
\begin{matrix}
\ Z(1) \\[1ex]
\ Z(2) \\[1ex]
\vdots \\[1ex]
\ Z(k)
\end{matrix}
\right],\quad \pmb H_k= \left[
\begin{matrix}
\ \pmb H(1) \\[1ex]
\ \pmb H(2) \\[1ex]
\vdots \\[1ex]
\ \pmb H(k)
\end{matrix}
\right],\quad \pmb N_k= \left[
\begin{matrix}
\ N(1) \\[1ex]
\ N(2) \\[1ex]
\vdots \\[1ex]
\ N(k)
\end{matrix}
\right]\tag{3}
Zk= Z(1) Z(2)⋮ Z(k),Hk= H(1) H(2)⋮ H(k),Nk= N(1) N(2)⋮ N(k)(3)
则 kkk 个观测数据满足如下方程：
Zk=Hkθ+Nk(4)
\pmb Z_k = \pmb H_k\pmb\theta + \pmb N_k\tag{4}
Zk=Hkθ+Nk(4)
前面已经说过，最小二乘估计方法要求测量模型输出 H(k)θ\pmb H(k)\pmb\thetaH(k)θ 与实际测量数据输出之差 Z(k)Z(k)Z(k) 的平方和最小，即 ∑(Z(k)−H(k)θ^)2\sum\big(Z(k)-\pmb H(k)\hat{\pmb \theta}\big)^2∑(Z(k)−H(k)θ^)2 达到最小。将求和性能指标写成下面矩阵形式：
J(θ^)=(Zk−Hkθ^)T(Zk−Hkθ^)(5)
\pmb J(\hat{\pmb\theta})=(\pmb Z_k - \pmb H_k\hat{\pmb\theta})^T(\pmb Z_k - \pmb H_k\hat{\pmb\theta})\tag{5}
J(θ^)=(Zk−Hkθ^)T(Zk−Hkθ^)(5)
将式（4）代入式（5）后，可以得到
J(θ^)=(Hkθ+Nk−Hkθ^)T(Hkθ+Nk−Hkθ^)=(θ−θ^)THkTHk(θ−θ^)+(θ−θ^)THkTNk+NkTHk(θ−θ^)+NkTNk(6)
\begin{aligned}
\pmb J(\hat{\pmb\theta})&=( \pmb H_k\pmb\theta + \pmb N_k - \pmb H_k\hat{\pmb\theta})^T( \pmb H_k\pmb\theta + \pmb N_k - \pmb H_k\hat{\pmb\theta})\\[1ex]
&=(\pmb\theta-\hat{\pmb\theta})^T\pmb H_k^T\pmb H_k(\pmb\theta-\hat{\pmb\theta}) + (\pmb\theta-\hat{\pmb\theta})^T\pmb H_k^T\pmb N_k+\pmb N_k^T\pmb H_k(\pmb\theta-\hat{\pmb\theta}) + \pmb N_k^T\pmb N_k
\end{aligned}\tag{6}
J(θ^)=(Hkθ+Nk−Hkθ^)T(Hkθ+Nk−Hkθ^)=(θ−θ^)THkTHk(θ−θ^)+(θ−θ^)THkTNk+NkTHk(θ−θ^)+NkTNk(6)
如果真实的参数 θ\pmb\thetaθ 和估计 θ^\hat{\pmb \theta}θ^ 与测量噪声 Nk\pmb N_kNk 无关，上式两边取期望可以得到
J(θ^)=E[(θ−θ^)THkTHk(θ−θ^)]+E(NkTNk)(7)
\pmb J(\hat{\pmb\theta})=\Bbb E\big[(\pmb\theta-\hat{\pmb\theta})^T\pmb H_k^T\pmb H_k(\pmb\theta-\hat{\pmb\theta}) \big] + \Bbb E\big( \pmb N_k^T\pmb N_k\big)\tag{7}
J(θ^)=E[(θ−θ^)THkTHk(θ−θ^)]+E(NkTNk)(7)
仔细观测发现，发现该式由两项为正的子式加和组成，后一项表明测量噪声的方差和估计参数无关。对于前一项来说，真实的参数 θ\pmb\thetaθ 和估计参数 θ^\hat{\pmb\theta}θ^ 越接近，该项的数值就会越小，因此，只要对对性能指标函数 J(θ^)\pmb J(\hat{\pmb\theta})J(θ^) 求极小值，找到的 θ^\hat{\pmb\theta}θ^ 就是所求得估计参数。
∂∂θ^[(Zk−Hkθ^)T(Zk−Hkθ^)]=−2HkT(Zk−Hkθ^)(8) \frac{\partial}{\partial\hat{\pmb\theta}}[(\pmb Z_k - \pmb H_k\hat{\pmb\theta})^T(\pmb Z_k - \pmb H_k\hat{\pmb\theta})]=-2\pmb H_k^T(\pmb Z_k-\pmb H_k\hat{\pmb\theta})\tag{8} ∂θ^∂[(Zk−Hkθ^)T(Zk−Hkθ^)]=−2HkT(Zk−Hkθ^)(8)
令式（8）等于 0，可得
HkTZk−HkTHkθ^=0(9)
\pmb H_k^T\pmb Z_k-\pmb H_k^T\pmb H_k\hat{\pmb\theta}=0\tag{9}
HkTZk−HkTHkθ^=0(9)
进一步推导可得
HkTZk=HkTHkθ^(10)
\pmb H_k^T\pmb Z_k=\pmb H_k^T\pmb H_k\hat{\pmb\theta}\tag{10}
HkTZk=HkTHkθ^(10)
假设 HkTHk\pmb H_k^T\pmb H_kHkTHk 是满秩的，可得
θ^=(HkTHk)−1HkZk(11)
\boxed{\hat{\pmb\theta}=(\pmb H_k^T\pmb H_k)^{-1}\pmb H_k\pmb Z_k}\tag{11}
θ^=(HkTHk)−1HkZk(11)
利用上述估计方法估计室内温度的变化可得 θ^=[ θ^1 θ^2]=[ 1.9808 20.0125]\hat{\pmb\theta}=\left[ \begin{matrix} \ \hat\theta_1 \\[1ex] \ \hat\theta_2 \end{matrix} \right]=\left[ \begin{matrix} \ 1.9808 \\[1ex] \ 20.0125 \end{matrix} \right]θ^=[ θ^1 θ^2]=[ 1.9808 20.0125]
2、练习题
美国在 1946 ~ 1956 年的钢产量分别为 66.6 百万吨、84.9 百万吨、88.6 百万吨、78.0 百万吨、96.8 百万吨、105.2 百万吨、93.2 百万吨、111.6 百万吨、88.3 百万吨、117.0 百万吨、115.2 百万吨，请分别用一次线性曲线、二次曲线、三次曲线以及四次曲线来拟合这些数据，并预测1957 ~ 1960年的钢产量。
（1）一次曲线方程为
z(k)=a1k+a0+w(k)(12)
z(k) = a_1k+a_0+w(k)\tag{12}
z(k)=a1k+a0+w(k)(12)
写成矩阵形式测量方程
z(k)=[ k1 ][ a1 a0]+w(k)(13)
z(k) = \textbf{[}\ k\quad1\ \textbf{]}
\left[
\begin{matrix}
\ a_1 \\[1ex]
\ a_0
\end{matrix}
\right]
+w(k)\tag{13}
z(k)=[ k1 ][ a1 a0]+w(k)(13)
Hk\pmb H_kHk 和待估计参数 θ\pmb\thetaθ 分别为
Hk=[112131⋮⋮111],θ=[ a1 a0](14)
\pmb H_k = \left[
\begin{matrix}
1 & 1 \\[1ex]
2 & 1 \\[1ex]
3 & 1 \\[1ex]
\vdots & \vdots\\[1ex]
11 & 1
\end{matrix}
\right],\quad \pmb\theta=\left[
\begin{matrix}
\ a_1 \\[1ex]
\ a_0
\end{matrix}
\right]\tag{14}
Hk=123⋮11111⋮1,θ=[ a1 a0](14)
（2）二次曲线测量方程为
z(k)=a2k2+a1k+a0+w(k)=[ k2k1 ][ a2 a1 a0]+w(k)(15)
\begin{aligned}
z(k)&=a_2k^2 + a_1k + a_0 + w(k)\\
&=\textbf{[}\ k^2\quad k\quad 1\ \textbf{]}
\left[
\begin{matrix}
\ a_2 \\[1ex]
\ a_1 \\[1ex]
\ a_0
\end{matrix}
\right]
+w(k)
\end{aligned}\tag{15}
z(k)=a2k2+a1k+a0+w(k)=[ k2k1 ] a2 a1 a0+w(k)(15)
与前面类似，
Hk=[11122213231⋮⋮112111],θ=[ a2 a1 a0](16)
\pmb H_k = \left[
\begin{matrix}
1 & 1 & 1\\[1ex]
2^2 & 2 & 1 \\[1ex]
3^2 & 3 & 1 \\[1ex]
\vdots & \vdots\\[1ex]
11^2 & 11 & 1
\end{matrix}
\right],\quad \pmb\theta=\left[
\begin{matrix}
\ a_2 \\[1ex]
\ a_1 \\[1ex]
\ a_0
\end{matrix}
\right]\tag{16}
Hk=12232⋮112123⋮111111,θ= a2 a1 a0(16)
（3）三次曲线同理
（4）四次曲线
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
