数据同化 最优插值 三维变分 原理 简明教程

什么是数据同化？简单来说，对于你要研究的变量，你有模型模拟的结果，存在偏差，也有测量仪器观测的结果，同样存在误差，那么数据同化就是研究如何将两者结合起来，得到该变量更加精准的值

问题一：简述陆面过程模型和卫星观测的特点和误差来源

陆面过程模型

特点

1. 陆面数据同化是将数据同化方法应用于地球表层科学和水文学中。
2. 陆面过程模拟主要研究采用不同的数据同化算法，同化地表观测数据，卫星和雷达数据，优化地表和根区水分，温度，地表能量通量等的估算。
3. 陆面过程模型初始化和运行的许多参数都通过遥感数据获得。
4. 陆面过程模型的优势在于其可以通过内在的物理过程和动力学机制，模拟对象在时间和空间上的连续演进。
5. 陆面过程模型的初值，驱动和参数都具有不确定性。

误差源

1. 模型误差：主要指模型结构误差，参数误差和驱动数据带来的误差和模型计算误差等。其中模型计算误差是指空间，时间差分、截断误差以及其他的计算误差带来的误差
2. 观测误差：由观测算子误差、仪器误差和代表性误差组成
3. 同化算法误差：不同的数据同化方法带来的误差。

卫星观测

特点

1. 通过卫星观测可以得到测量对象在观测时刻和所代表空间上的真值
2. 观测为瞬间值，在时间上不连续
3. 卫星观测难以获取地表以下信息
4. 卫星观测不是模式变量，需要借助反演算法和观测算子

误差源

1. 与卫星有关的误差：卫星仪器误差，地球自转影响，相对论效应
2. 与信号传播有关的误差：大气效应
3. 代表性误差

问题二：简述逐步订正法，最优插值以及三维变分的异同点

相同点

1. 这三种方法都是基于一定的数学原理，将观测值和模式值结合起来，为模型提供精确的初值，都属于客观分析方法。
2. 逐步订正法和最优插值法本质上都是一种插值方法。

不同点

1. 模型原理不同：逐步订正法和最优插值是基于最小方差的原理，三维变分法是基于最大似然估计和贝叶斯理论。逐步订正法和最优插值法本质上都是一种插值方法，逐步订正法的权重函数为经验函数，而最优插值的权重函数是基于最小二乘导出，因此最优插值法比逐步订正法的权重函数更为合理。三维变分方法利用变分思想将数据同化问题转化为一个极值求解问题，在满足动态约束的条件下，最小化状态预测值与观测值之间的距离，使这种“距离”最小的状态量即为最优状态估计量。
2. 模型假设不同：逐步订正法较为简单假设较少，最优插值法假设背景场和观测场误差都是无偏的，非线性观测算子可以线性化，以及观测场和背景场误差不相关。三维变分法假设观测误差和模式结果误差均服从均值为0的高斯分布。
3. 模型求解方法不同：逐步订正法先给定第一猜场，然后用实际观测场逐步修正第一猜测场，直到订正后的场逼近观测记录，其中权重函数为与距离相关的经验函数。最优插值法与逐步订正法类似，只不过在计算权重时不再使用经验函数，而是引入误差协方差矩阵。三维变分法通过最小化目标函数来获得最优估计。其中，如果观测算法算子为线性，可使用正则方程求解，如果观测算子为非线性，可使用迭代法求解，比如拟牛顿法，共轭梯度法，牛顿下降法，最速下降法等。

问题三：简述最优插值和三维变分的基本原理

最优插值原理

最优插值法是在逐步订正法的基础上利用最小二乘原理求解权重矩阵，得到最优估计。假设$x_i^b,x_i^a,x_i^t$分别为第i个分析格点上的背景场、分析场和真值。则有

$x_i^a=x_i^b+{\sum }_{k=1}^{K_i}{W}_{ki}(x_k^{obs}-x_k^b)=x_i^b+{\mathbf{w}}_{\mathbf{i}}^{\mathbf{T}}({\mathbf{x}}^{\mathbf{o}\mathbf{b}\mathbf{s}}-{\mathbf{x}}^{\mathbf{b}})$

为了使得分析场误差方差最小，即求解下面的极小化问题：

$mi{n}_{w_i}=mi{n}_{w_i}(x_i^a-x_i^t)(x_i^a-x_i^t)$

假设背景场合观测误差都无偏合无相关，观测误差合背景误差无相关，为了使得分析场的误差方差最小，求一阶导数为0，可得：

${\mathbf{w}}_{\mathbf{i}}=(\mathbf{B}+\mathbf{O}{)}^{-1}{\mathbf{b}}_{\mathbf{i}}$

其中B矩阵O矩阵和bi矩阵如下所示

$B=\left(\begin{array}{cccc}\overline{(x_1^b-x_1^t)(x_1^b-x_1^t)}& \overline{(x_1^b-x_1^t)(x_2^b-x_2^t)}& \cdots & \overline{(x_1^b-x_1^t)(x_{Ki}^b-x_{Ki}^t)}\\ \overline{(x_2^b-x_2^t)(x_1^b-x_1^t)}& \overline{(x_2^b-x_2^t)(x_2^b-x_2^t)}& \cdots & \overline{(x_2^b-x_2^t)(x_{Ki}^b-x_i^t)}\\ ⋮& ⋮& ⋮& ⋮\\ \overline{(x_{Ki}^b-x_{Ki}^t)(x_1^b-x_1^t)}& \overline{(x_{Ki}^b-x_{Ki}^t)(x_2^b-x_2^t)}& \cdots & \overline{(x_{Ki}^b-x_{Ki}^t)(x_{ki}^b-x_{ki}^t)}\end{array}\right)$

$O=\left(\begin{array}{cccc}\overline{(x_1^{obs}-x_1^t)(x_1^{obs}-x_1^t)}& \overline{(x_1^{obs}-x_1^t)(x_2^{obs}-x_2^t)}& \cdots & \overline{(x_1^{obs}-x_1^t)(x_{Ki}^{obs}-x_{Ki}^t)}\\ \overline{(x_2^{obs}-x_2^t)(x_1^{obs}-x_1^t)}& \overline{(x_2^{obs}-x_2^t)(x_2^{obs}-x_2^t)}& \cdots & \overline{(x_2^{obs}-x_2^t)(x_{Ki}^{obs}-x_i^t)}\\ ⋮& ⋮& ⋮& ⋮\\ \overline{(x_{Ki}^{obs}-x_{Ki}^t)(x_1^{obs}-x_1^t)}& \overline{(x_{Ki}^{obs}-x_{Ki}^t)(x_2^{obs}-x_2^t)}& \cdots & \overline{(x_{Ki}^{obs}-x_{Ki}^t)(x_{ki}^{obs}-x_{ki}^t)}\end{array}\right)$

$b_i=\left(\begin{array}{c}(x_1^b-x_1^t)(x_i^b-x_i^t)\\ (x_2^b-x_2^t)(x_i^b-x_i^t)\\ ⋮\\ (x_{Ki}^b-x_{Ki}^t)(x_i^b-x_i^t)\end{array}\right)$

其中，B为背景协方差矩阵，O为观测场协方差矩阵。两者均需要真值来推导，但是由于不知道所谓的均值，因此无法准确快速求解协方差均值。因此最优插值法的难点在于协方差矩阵的推导。

三维变分法原理

由上面的推导可知，要想求解出最优差值法的权重函数，需要$(B+O)$矩阵可逆，而在实际的应用中，由于矩阵的过大，求逆计算量较大，这种情况则可以使用三维差分方法来解决。变分方法构建代价函数来描述状态量分析值和真值之间的差异，利用变分思想把数据同化问题转化为极值求解问题。下面以单点的三维变分为例，讲述三维变分的原理：

设模式模拟结果为$x_b$，标准差为${\sigma }_b$,，观测值为$x_{obs}$ ，标准差为${\sigma }_{obs}$，需要求得的分析值为$y$。

由Bayes理论可知:$P(y|x_{obs})\ast P(x_{obs})=P(x_{obs}|x_b)\ast P(x_b)$

又由于$x_{obs}$是给定的，所以$P(x_{obs})=const$

综上：$P(y|x_{obs})\sim P(x_{obs}|x_b)\ast P(x_b)$

对于三维变分法，其假设模拟误差与观测误差均服从正态分布。

则$P(x_{obs}|x_b)=\frac{1}{\sqrt{2\pi }{\sigma }_{obs}}exp[-\frac{(x-x_{obs}{)}^2}{2{\sigma }_{obs}^2}]$

$P(x_b)=\frac{1}{\sqrt{2\pi }{\sigma }_b}exp[-\frac{(x-x_b{)}^2}{2{\sigma }_b^2}]$

所以$P(y|x_{obs})=\frac{1}{2\pi {\sigma }_b{\sigma }_{obs}}exp[-\frac{(x-x_b{)}^2}{2{\sigma }_b^2}-\frac{(x-x_{obs}{)}^2}{2{\sigma }_{obs}^2}]$

由此，我们可以得到分析值的概率分布，要求解最优分析值，实际上是求解概率最大值。

对上述概率取对数得到代价函数$ J(x) = - ln[P(y|x_{obs})] = -ln({{2\pi}{\sigma_{b}}{\sigma_{obs}}}) + {\frac{(x - x_b)^{2}{2\sigma_{b}}2}}{+\frac{(x - x_{obs})^{2}{2\sigma_{obs}}2}} $

所以$y=argminJ(x)=\frac{{\sigma }_b^2{x}_{obs}+{\sigma }_{obs}^2{x}_b}{{\sigma }_b^2+{\sigma }_{obs}^2}(1.1)$

对于此问题而言，代价函数为凸函数，只需求导即可，而在实际问题中，代价函数形式较为复杂，无解析最优值，通常需要迭代法求解，比如拟牛顿法，共轭梯度法，牛顿下降法，最速下降法等。