线性估计

线性估计

设信号模型为\(x(t)=a+bt+N(t)\)​​​​​​​，其中\(a,b\)​​​​​​​为待估计参数，\(N(t)\)​​​​​​​表示噪声。我们首先需要对上述信号进行采样，设采样间隔为\(\Delta t\)​​​​​​​，则第\(k\)​​​​​​​个采样点可以表示为\(x_k=a+b k\Delta t+N_k\)​​​​​​​​，现在需要利用采样点对\(a,b\)​​​进行估计。定义如下性能指标：

\[g(a,b)=\sum_{k=1}^n[(x_k-a-b k\Delta t)^2]\tag{1} \]

(1)表示所有采样点与信号模型之间的垂直距离的平方和最小。为描述简便，我们将上述问题用矩阵形式表示出来，为此定义：\(\boldsymbol{x}=[x_1,...,x_n]^{\rm T}\)​表示采样序列，\(\boldsymbol{\theta}=[a,b]^{\rm T}\)​​表示待估计参数，\(\boldsymbol{H}=\begin{bmatrix}1,\Delta t\\ \vdots ,\vdots\\ 1,n\Delta t\end{bmatrix}\)表示系数矩阵，\(\boldsymbol{N}=[N_1,...,N_n]^{\rm T}\)​表示噪声向量，则信号模型和性能指标可以表示为

\[\begin{equation} \begin{aligned} \boldsymbol{x}&=\boldsymbol{H \theta}+\boldsymbol{N}\\ g(\boldsymbol{\theta})&=(\boldsymbol{x}-\boldsymbol{H\theta})^{\rm T}(\boldsymbol{x}-\boldsymbol{H\theta}) \end{aligned} \end{equation}\tag{3} \]

显然，为求\(\boldsymbol{\theta}\)，一个直观的思路是对\(g(\boldsymbol{\theta})\)求梯度，为此，先对\(g(\boldsymbol{\theta})\)展开整理：

\[\begin{equation} \begin{aligned} g(\boldsymbol{\theta})&=\boldsymbol{x}^{\rm T}\boldsymbol{x}-\boldsymbol{x}^{\rm T}\boldsymbol{H \theta}-\boldsymbol{\theta}^{\rm T}\boldsymbol{H}^{\rm T}\boldsymbol{x}+\boldsymbol{\theta}^{\rm T}\boldsymbol{H}^{\rm T}\boldsymbol{H \theta} \end{aligned} \end{equation}\tag{3} \]

因此

\[\nabla_{\boldsymbol{\theta}}g(\boldsymbol{\theta})=-2\boldsymbol{H}^{\rm T}\boldsymbol{x}+2\boldsymbol{H}^{\rm T}\boldsymbol{H}\boldsymbol{\theta}\tag{4} \]

在上面的求导过程中，可以应用一些现有的求导公式：\(\nabla_{\boldsymbol{\theta}}(\boldsymbol{A \theta})=\boldsymbol{A}^{\rm T}\)，\(\nabla_{\boldsymbol{\theta}}(\boldsymbol{\theta^{\rm T} A})=\boldsymbol{A}\)，\(\nabla_{\boldsymbol{\theta}}(\boldsymbol{\theta^{\rm T} A \boldsymbol{\theta}})=(\boldsymbol{A}^{\rm T}+\boldsymbol{A})\boldsymbol{\theta}\)​。

在(4)中，令\(\nabla_{\boldsymbol{\theta}}g(\boldsymbol{\theta})=0\)可以得到\(\boldsymbol{\theta}=(\boldsymbol{H}^{\rm T}\boldsymbol{H})^{-1}\boldsymbol{H}^{\rm T}\boldsymbol{x}\)。该结果称为最小二乘解。

---

上面从线性拟合的角度求出了线性模型中的参数，整个过程中并没有考虑到噪声的分布，相当于忽略了噪声的影响。下面从统计的角度来重新求一下上述线性模型中的参数。设\(N_k \backsim N(0,\sigma^2)\)，则要求\(\boldsymbol{\theta}\)的最优估计，可从CRLB入手，因此进行如下操作:

(1) 列出所有采样样本\(\boldsymbol{x}=[x_1,...,x_n]^{\rm T}\)的联合概率密度函数:

\[f(\boldsymbol{x},\boldsymbol{\theta})=(\frac{1}{\sqrt{2\pi}\sigma})^n {\rm exp}[-\frac{(\boldsymbol{x-\boldsymbol{H \theta}})^{\rm T}(\boldsymbol{x-\boldsymbol{H \theta}})}{2\sigma^2}]\tag{5} \]

(2) 对上述联合概率密度求对数:

\[{\rm ln}[f(\boldsymbol{x},\boldsymbol{\theta})]=-n{\rm ln}(\sqrt{2\pi}\sigma)-\frac{1}{2\sigma^2}(\boldsymbol{x-\boldsymbol{H \theta}})^{\rm T}(\boldsymbol{x-\boldsymbol{H \theta}})\tag{6} \]

(3) 对对数联合概率密度函数求梯度：

\[\begin{equation} \begin{aligned} \nabla_{\boldsymbol{\theta}}({\rm ln}[f(\boldsymbol{x},\boldsymbol{\theta})])&=\frac{1}{\sigma^2}(\boldsymbol{H}^{\rm T}\boldsymbol{x}-\boldsymbol{H}^{\rm T}\boldsymbol{H}\boldsymbol{\theta})\\ &=\frac{\boldsymbol{H}^{\rm T}\boldsymbol{H}}{\sigma^2}[(\boldsymbol{H}^{\rm T}\boldsymbol{H})^{-1}\boldsymbol{H}^{\rm T}\boldsymbol{x}-\boldsymbol{\theta}] \end{aligned} \end{equation}\tag{7} \]

到这边其实已经没有必要继续往下算了，从(7)式就可以得到待估参数\(\boldsymbol{\theta}\)对应的Fisher矩阵\(\boldsymbol{I}(\boldsymbol{\theta})=\frac{\boldsymbol{H}^{\rm T}\boldsymbol{H}}{\sigma^2}\)，而最优估计量为\(\hat{\boldsymbol{\theta}}=\boldsymbol{H}^{\rm T}\boldsymbol{H})^{-1}\boldsymbol{H}^{\rm T}\boldsymbol{x}\)，显然与上面前面从线性拟合角度得到的结果是一致的。

从上面的分析可以知道，最小二乘估计结果是线性估计中的最优结果(此时不管噪声的具体分布)，同时也是高斯噪声下任意估计(不局限于线性估计)中的最优结果(因为此时最小二乘结果达到了CRLB下界)。

---

接下来进一步考虑线性估计问题，设此时信号模型为\(x_k=bs_k+N_k\)，其中\(b\)为待估参数。\(s_k\)在不同系统中可能具有不同的含义，比如在雷达系统中可能表示波形、在通信系统中可能表示码元。我们希望从采样序列\(\boldsymbol{x}=[x_1,...,x_n]^{\rm T}\)中估计出\(b\)，且这个估计值来源与采样序列的线性组合，即\(\hat{\theta}(\boldsymbol{x})=\sum_{k=1}^n{\alpha_k x_k}=\boldsymbol{\alpha}^{\rm T}\boldsymbol{x}\)，同时我们要求这个估计是无偏的，所以\({\rm E}[\hat{\theta}(\boldsymbol{x})]=\boldsymbol{\alpha}^{\rm T}{\rm E}[\boldsymbol{x}]=b\)​。据此，我们得到我们的估计问题如下：

\[\underset{\boldsymbol{\alpha}}{{\rm min}}{\rm E}[(b-\boldsymbol{\alpha}^{\rm T} \boldsymbol{x})^2] \quad \quad s.t. \quad {\rm E}[\hat{\theta}(\boldsymbol{x})]=\boldsymbol{\alpha}^{\rm T}{\rm E}[\boldsymbol{x}]=b \tag{8} \]

对上面的式子进行整理:

\[\begin{equation} \begin{aligned} {\rm E}[(b-\boldsymbol{\alpha}^{\rm T} \boldsymbol{x})^2]&={\rm E}[b^2-2b\boldsymbol{\alpha}^{\rm T} \boldsymbol{x}+\boldsymbol{\alpha}^{\rm T}\boldsymbol{x}\boldsymbol{x}^{\rm T}\boldsymbol{\alpha}]\\ &=b^2-2b\boldsymbol{\alpha}^{\rm T}{\rm E}[\boldsymbol{x}]+\boldsymbol{\alpha}^{\rm T}{\rm E}[\boldsymbol{x}\boldsymbol{x}^{\rm T}]\boldsymbol{\alpha}\\ &=-b^2+\boldsymbol{\alpha}^{\rm T}{\rm E}[\boldsymbol{x}\boldsymbol{x}^{\rm T}]\boldsymbol{\alpha} \end{aligned} \end{equation}\tag{9} \]

(9)中含有待估参数\(b\)，这显然是我们不希望的，因此要想办法将(9)中的\(b\)​去掉。此时可以利用(8)中约束条件，将其代入目标函数，可以得到：

\[\begin{equation} \begin{aligned} {\rm E}[(b-\boldsymbol{\alpha}^{\rm T} \boldsymbol{x})^2]&={\rm E}[(\boldsymbol{\alpha}^{\rm T}{\rm E}[\boldsymbol{x}]-\boldsymbol{\alpha}^{\rm T} \boldsymbol{x})^2]\\ &={\rm E}[\boldsymbol{\alpha}^{\rm T}({\rm E}[\boldsymbol{x}]-\boldsymbol{x})({\rm E}[\boldsymbol{x}]-\boldsymbol{x})^{\rm T}\boldsymbol{\alpha}]\\ &=\boldsymbol{\alpha}^{\rm T}{\rm E}[{\rm E}[\boldsymbol{x}]-\boldsymbol{x})({\rm E}[\boldsymbol{x}]-\boldsymbol{x})^{\rm T}]\boldsymbol{\alpha}\\ &=\boldsymbol{\alpha}^{\rm T}{\rm Cov}(\boldsymbol{x})\boldsymbol{\alpha} \end{aligned} \end{equation}\tag{9} \]

而\({\rm E}[\boldsymbol{x}]-\boldsymbol{x})({\rm E}[\boldsymbol{x}]-\boldsymbol{x})^{\rm T}=\boldsymbol{N}\boldsymbol{N}^{\rm T}\)，所以\({\rm Cov}(\boldsymbol{x})={\rm E}[\boldsymbol{N}\boldsymbol{N}^{\rm T}]=\boldsymbol{C}_N\)，此外进一步对(8)中的约束条件进行整理(目的是消除待估参数\(b\))，可以得到\(\boldsymbol{\alpha}^{\rm T}{\rm E}[\boldsymbol{x}]=b=\boldsymbol{\alpha}^{\rm T}b\boldsymbol{s} \rightarrow \boldsymbol{\alpha}^{\rm T}\boldsymbol{s}=1\)，所以经过上面的化简(8)中的约束最优化问题转化为：

\[\underset{\boldsymbol{\alpha}}{{\rm min}}\boldsymbol{\alpha}^{\rm T}\boldsymbol{C}_N\boldsymbol{\alpha} \quad \quad s.t. \quad \boldsymbol{\alpha}^{\rm T}\boldsymbol{s}=1 \tag{10} \]

对于约束最优化问题，通常可以采用拉格朗日乘子法，为此构造下述目标函数:

\[L(\boldsymbol{\alpha},\lambda)=\frac{1}{2}\boldsymbol{\alpha}^{\rm T}\boldsymbol{C}_N\boldsymbol{\alpha}-\lambda(\boldsymbol{\alpha}^{\rm T}\boldsymbol{s}-1)\tag{11} \]

所以

\[\nabla_{\boldsymbol{\alpha}}L(\boldsymbol{\alpha},\lambda)=\frac{1}{2}(\boldsymbol{C}_N+\boldsymbol{C}_N^{\rm T})\boldsymbol{\alpha}-\lambda \boldsymbol{s}=\boldsymbol{C}_N \boldsymbol{\alpha}-\lambda \boldsymbol{s}\tag{12} \]

令\(\nabla_{\boldsymbol{\alpha}}L(\boldsymbol{\alpha},\lambda)=0 \rightarrow \boldsymbol{\alpha}=\lambda\boldsymbol{C}_N^{-1}\boldsymbol{s}\)​，所以\(\boldsymbol{\alpha}^{\rm T}=\lambda \boldsymbol{s}^{\rm T}\boldsymbol{C}_N^{-1}\)​，所以\(\boldsymbol{\alpha}^{\rm T}\boldsymbol{s}=\lambda \boldsymbol{s}^{\rm T}\boldsymbol{C}_N^{-1}\boldsymbol{s}=1\)​，所以\(\lambda=\frac{1}{\boldsymbol{s}^{\rm T}\boldsymbol{C}_N^{-1}\boldsymbol{s}}\)​，所以\(\boldsymbol{\alpha}=\frac{\boldsymbol{C}_N^{-1}\boldsymbol{s}}{\boldsymbol{s}^{\rm T}\boldsymbol{C}_N^{-1}\boldsymbol{s}}\)​，所以估计结果可以表示为\(\hat{\theta}(\boldsymbol{x})=\boldsymbol{\alpha}^{\rm T}\boldsymbol{x}=\frac{\boldsymbol{s}^{\rm T}\boldsymbol{C}_N^{-1}}{\boldsymbol{s}^{\rm T}\boldsymbol{C}_N^{-1}\boldsymbol{s}}\boldsymbol{x}\)​。该结果即为最优线性无偏估计。

需要注意该结果与最小方差无偏估计(MVUE)的差别，MVUE需要知道信号的详细分布，然后利用分布算CRLB，再寻求最优估计，它对任何估计(不局限于线性估计)都是最优的。而上面的最优线性无偏估计(BLUE)则不要求信号的具体分布，不过它只在线性估计范围内寻找最优估计结果。所以两者的前提和适用范围都是不一样的。