从概率统计角度对资料同化的简单理解.md

从概率统计角度对资料同化的简单理解

之前我已经在"先验-后验-似然-极大似然估计-最大后验估计"一文中阐述了概率统计的一些基本概念，最后我们可以发现，先验对贝叶斯派的估计影响很大。顺着前面这篇文章，这里我们简单讨论下概率统计角度如何看待资料同化问题。

资料同化的问题其实可以简单考虑为下面这样的一个系统：

$v_{j+1} = \Phi(v_j) + \xi_j$

$y_{j+1} = h(v_{j+1}) + \eta_{j+1}$

这里$j$代表某个时间步，$v_{j+1}$代表$j+1$时刻的大气变量（即信号，这里我们假设信号是马尔可夫的），$\xi_j$和$\eta_{j+1}$则是高斯噪声，$h$代表的是观测算子，$y_{j+1}$代表的是$j+1$时刻的观测值。资料同化的根本目的，实际上是根据观测$y$来估计信号$v$。

现在我们在上帝视角来看待上面那两个公式。作为上帝，我们已经知道大气的演变过程服从第一个公式（随机动力系统），由于人类的观测是有限的、不完美的，通过有限的方式，也只有有限的数量，所以利用第二个公式来描述这个情况：观测算子联系了大气状态空间和观测空间，同时需要考虑到观测这个过程中噪声的引入。OK，虽然上帝是能够看穿这一切的，但是我们渺小的人类明显不行。我们只能得到一系列的观测数据$y$，但是我们想通过这些数据来估计大气实际的状态是什么样的（$v$）。

到这里，我们前面提到的贝叶斯公式又要登场了。

$P(v|y)=\dfrac{P(y|v)P(v)}{P(y)}$

和前面那篇文章中提到的一样，由于我们已经得到了观测，$P(y)$仅仅是一个无关紧要的常数，因此，

$P(v|y) \varpropto P(y|v)P(v)$

这里我们需要确定先验和似然。为了清晰，我们这里限定时间步的下标从$0$开始到$j$。

$v_{j} = \Phi(v_{j-1}) + \xi_{j-1}$

$y_{j} = h(v_{j}) + \eta_{j}$

对于先验，我们可以有如下表述：

$P(v) = P(v_j, v_{j-1},...,v_0) = P(v_j|v_{j-1},...,v_0)P(v_{j-1},...,v_0)$

这一点根据联合概率分布是很容易得到的。由于我们假定了信号具有马尔可夫性，即对于一个过程，$j$时刻的状态仅仅由$j-1$时刻的状态决定。所以，我们可以得到：

$P(v) = P(v_j|v_{j-1})P(v_{j-1},...,v_0) = \prod_{k=0}^{j}P(v_k|v_{k-1})P(v_0)$

为了得到更具体表达式，我们假设初始时刻$v$服从高斯分布，那么，

$P(v_0) \propto exp(-\frac{1}{2}|C_0^{-\frac{1}{2}}(v_0-m_0)|^2)$

根据资料同化的第一个公式：

$v_{j} = \Phi(v_{j-1}) + \xi_{j-1}$

我们可以得到：

$P(v_k|v_{k-1}) \propto exp(-\frac{1}{2}|\Sigma^{-\frac{1}{2}}(v_k-\Phi(v_{k-1}))|^2)$

其中，$C_0$，$\Sigma$代表协方差。这样，我们就可以得到$P(v)$的表达式：

$P(v) \propto exp(-J(v))$,

$J(v) = \frac{1}{2}|C_0^{-\frac{1}{2}}(v_0-m_0)|^2+\Sigma_{k=0}^{j}\frac{1}{2}|\Sigma^{-\frac{1}{2}}(v_k-\Phi(v_{k-1}))|^2$

注意第一个$\Sigma$是求和，第二个是协方差。

下面来看似然。

$P(y|v) = \prod_{k=0}^{j-1}P(y_{k+1}|v) = \prod_{k=0}^{j-1}P(y_{k+1}|v_{k+1})$

根据资料同化的第二个公式：

$y_{j} = h(v_{j}) + \eta_{j}$

我们可以发现，

$P(y|v) \propto \prod_{k=0}^{j-1}exp(-\frac{1}{2}|\Gamma^{-\frac{1}{2}}(y_{k+1}-h(v_{k+1}))|^2)$

这样我们就得到了先验和似然的表达，从而能够根据$y$来估计$v$。一般的文献中，$m_0$被称为背景场均值，而$C_0$被称为背景协方差矩阵。

下面我们先把理论结合实际来看一看。有WRF同化经历的人都知道，背景误差协方差可以通过NMC方法来获得。NMC方法对于区域应用，一般是T+24预报减去T+12预报，对于全球应用，一般是T+48预报减去T+24预报。

为了便于理解，我们拿出资料同化的优化目标函数：

$J(x)=(x-x_b)^TB^{-1}(x-x_b)+(y-H(x))^TR^{-1}(y-H(x))$

其中，$x_b$是背景场的值，$B$是背景场误差协方差，$H$是观测算子，$y$是观测值，$R$是观测误差协方差。在同化中，背景场协方差可以把观测信息从观测点传播到周围的格点和垂直层上。而背景场往往是上一个时刻模式的预报结果。

回顾下刚刚所述的资料同化的第一个公式，

$v_{j} = \Phi(v_{j-1}) + \xi_{j-1}$

这里的$v$对应着目标函数中的$x$，背景场通过模式预报而来，即对应$\Phi(v_{j-1})$，我们可以联想到NMC方法其实就是在估计高斯噪声$\xi$的情况，只是这个噪声被称为误差。

另外，目标函数的结构，和我们推导的后验概率分布有着非常好的对应关系。目标函数右侧一共有两项，第一项代表的是模式信息，第二项代表的是观测信息。所以当我们在优化目标函数得到分析值时，我们实际上是在利用极大后验方法在估计信号$v$。

现在我们重新看一下数值模式和观测资料在同化中扮演的作用。模式这部分东东对应的其实是资料同化的第一个公式，实际上模式提供了先验信息；观测对应的是资料同化的第二个公式，它们提供了似然。换种角度说，模式实际上是人类根据目前的认知，对大气运动状态的一种近似，换言之，我们先入为主地认为大气的运动是必须符合这样的状态的，这部分知识是我们已经知道的，所以是“先验”；而观测得到的是一堆数据，我们需要用数据来估计大气真实的状态，这就是“似然”，“看起来好像是这个样子”。

在之前关于先验后验那篇文章的末尾，我们注意到，在极大后验方法中，如果试验次数较少（即观测较少）的情况下，先验对估计的结果影响很大，换句话说，如果观测不够多，那么估计的结果主要被模式影响，而如果观测特别特别特别多，那么模式这一“先验”的信息作用就不强了。