引言

大气资料同化源于数值天气预报(Numerical Weather Prediction, NWP)初值形成的客观分析：它就是在功能上实现按照某种算法通过计算机程序(无需人工干预)将不规则分布的观测资料结合天气动力学和观测资料特征而得到某一预定规则网格点上最可能的值(Panofsky，1949)，即表示尽可能准确的大气状态的分析场，来提供NWP预报模式的初值条件。其基本的主流方法(国际上主要业务NWP中心使用的分析同化方法)先后经历了多项式函数拟合方法(Panofsky，1949)、逐步订正方法(Successive Correction Method, SCM；Bergthorsson et al, 1955；Cressman, 1959)、最优插值方法(Optimal Interpolation, OI；Gandin，1965；Bergman，1979；Lorenc，1981)、变分方法(Variational methods, Var；Lorenc，1986；Parrish et al，1992；Courtier et al，1998)和集合卡尔曼滤波方法(Ensemble Kalman Filter, EnKF；Evensen，1994；Hamill，2006)。

早期的基本主流方法中，Panofsky(1949)的多项式拟合方法的思想及其所用的三次多项式是经验假定的，Cressman(1959)的订正公式形式及其订正权重是经验给定的，所以它们虽然实现了通过计算机程序进行客观分析的功能，但是主要地表现为一项经验工艺。

随着客观分析方法的发展，OI以分析场误差方差最小为标准来确定统计最优权重，Var由分析场似然概率最大的解转化为等价于一个目标泛函最小的解来得到分析场。这期间，预报场作为客观分析的先验背景场的引入，使得客观分析是通过观测信息本身和NWP 短期预报的统计结合产生分析场：

$ {\boldsymbol{x}_a} = {\boldsymbol{x}_b} + \boldsymbol{W}\left[ {{\boldsymbol{y}_o} - H\left({{\boldsymbol{x}_b}} \right)} \right] $

(1)

式中，x_a和x_b分别表示分析场和背景场，y_o表示观测；H是观测算子(也称为向前观测算子)，它给出了大气状态物理量和观测物理量之间的具体确定联系；W是统计权重。这使得预报和分析构成循环，即：前期的NWP模式预报场作为分析某一时刻观测资料的先验背景场，得到分析场，这被称为分析步骤；之后以该分析场作为初始条件积分预报模式到下一个时刻，得到预报场，这被称为预报步骤；再由这个预报场作为分析新时刻观测资料的背景场，得到分析场，这是又一循环的分析步骤；如此循环。在这个循环过程中不同时刻的观测资料和大气的动力数值模式这两者融合来尽可能精确地确定当前时刻的大气状态。正是由此造出了“同化”这个字眼，表示这个过程(Talagrand，1997)。这样，“资料同化”除了实现客观分析的功能之外还有了“同化”字眼的自身内涵，“资料同化”也成为现在越来越普遍使用的术语。因此，自OI开始的基本主流方法成为“理论上基于估计理论”的一门独立的科学工程学科。

大气资料同化发展到今天，它的具体实现表现为各种可用信息的资料融合，来尽可能准确地确定大气状态。这些可用信息包括观测资料本身和大气运动状态的先验背景估计(如NWP短期预报场)。同时，为了使用这些信息，还应用大气运动演变的物理定律(如四维变分中的NWP 预报模式)和大气运动的统计或动力性质(如中纬度地区的近似地转平衡)等。

是什么使得大气资料同化有了理论上的数学基础和数学表示形式而从一项经验工艺发展成为一门科学工程学科？我们知道，随机变量是概率论与数理统计学的基本概念。本文将阐明随机变量是大气资料同化的根本概念。其根本性可以体现在两个方面：(1) 原理上，大气资料同化的科学实质是建立在随机变量的概念上，从而使大气资料同化有了基于估计理论的最优标准及其数学形式；(2) 实施上，基于随机变量这个概念下的数据可以理解资料同化发展进程中的可用信息在内涵和种类上的不断扩展，这个扩展是上下承接、循序渐进的，在认识上是非常清晰、简单的。

1 随机变量及其在资料同化和NWP的实际应用意义 1.1 随机变量

数学上随机变量的概念是相当简单的。如果一个变量以一定的概率取各种可能值，则该变量称为随机变量(random variable)(陆果，1997)。

一个随机变量由它的概率分布(概率函数或概率密度函数)完整地描述，即：随机变量的概率分布包含了这个随机变量的所有信息。离散型随机变量的概率函数是取各个可能值的概率，连续型随机变量的概率密度函数(Probability Density Function, PDF)反映概率在各个可能值上的“密集程度”。随机变量的数字特征是由概率函数或概率密度函数决定的常数，它描述了随机变量的某一方面的性质；如：随机变量的均值、方差和协方差分别描述了随机变量的中心位置、散布度和与其他随机变量之间的关系(陈希孺，2009)。对于连续型的随机变量X，用x表示其可能取值，其概率密度函数表示为p(x)；以均值这个数字特征为例，它被定义为X=〈X〉=∫xp(x)dx，该积分遍及X的取值范围[x₁, x_n]，且可积。

1.2 随机变量在资料同化与NWP的实际应用意义

虽然随机变量的概念是相当简单的，然而它的实际应用有着深刻的意义，可以有以下的理解：

(1) 数学上，随机变量是“各种可能值”的可能性(随机性)和“一定概率”的确定性这二者不可缺一的一体，是偶然性和必然性的一体。可能性和确定性是随机变量的一体二面：可能性是指存在多种可能的取值，确定性是指对于任一可能值其概率是一定的。

(2) 物理测量上，一个基本的真实是：测量都有误差，真值永远不可知。所以，对于任一具体物理对象，如某时某地的温度，我们不能够知道其测量值是否是真值，只能看作是真值的一个可能值。为了得到尽可能准确的值，就进行多次的测量(在2.3.3 节中阐明)。每次测量值都看作是真值的一个可能值，并假定各个可能值的概率(或是测量误差的概率)是一定的(如符合正态分布)；这样，把测量对象看作一个随机变量，于是任一测量值就成为该随机变量的一个可能值。具体地，某时某地的温度作为一个具体的测量对象，把它看作一个随机变量T，则任一测量值t是随机变量T的一个可能值，t对应的概率是P(T=t)。

(3) 因此，随机变量概念的一个实际应用是：把测量对象看作一个随机变量，于是每次测量值成为该随机变量的一个可能值。由此，具体地，把“某一观测的温度物理量”看作一个随机变量，其观测值是该随机变量的一个可能值；同样地，把“某一先验(或称背景)的温度物理量”看作一个随机变量，其先验值是该随机变量的一个可能值；“把待确定的温度物理量”作为一个随机变量，所欲求的分析值是该随机变量在某一最优标准下的一个尽可能准确值。

进一步地，这个实际应用的深刻意义是：以随机变量为研究对象的概率论与数理统计学就自然地成为资料同化的数学基础。“使用观测和先验的信息来尽可能准确地确定大气运动的状态”这一问题，可以表述为由若干已知随机变量的信息求解未知随机变量的最优值的估计问题。如柏拉图说：“数学是一切知识中的最高形式”。在第2节将看到，因为随机变量，资料同化有了基于估计理论的最优标准及其数学公式表述形式[式(5)、(6)和(9)]。

在这个意义上，资料同化的科学实质是建立在随机变量的概念上。

(4) 随机变量概念的另一个实际应用是：资料同化中可用信息的数据成为一个随机变量概念下的数据。于是，数据的内涵，不再简单地只是一个数值，而是包含值 (value)、准确度 (accuracy)、精度 (precision)与不确定性 (uncertainty)和协方差 (covariance)。其含义为：①数据的值是作为一个随机变量的一个可能的值；②数据的准确度是一个随机变量的统计均值(即其中心位置)与真值的接近程度，表征了数据的系统误差的无偏性和可靠性，可表示为偏差b，它在资料同化中直接联系着资料的偏差订正；数据的精度是一个随机变量在其中心位置附近的散布程度，表征了数据的随机误差的离散度, 可表示为标准差s，它在资料同化中直接联系着资料的方差指定；数据的不确定性是其系统误差和随机误差的综合表征，用u表示，则有$ u = \sqrt {{b^2} + {s^2}} $；③协方差表征了不同数据之间的相关联系。由于事实上测量都有误差以及真值永远都不可知，所以随机变量概念下的数据使得资料同化中的数据概念更符合真实而有了更科学的上述内涵。

在第3节将看到，基于随机变量概念下的数据可以清晰地揭示和理解在资料同化发展进程中可用信息在内涵和种类上的不断扩展。

(5) 补充上述的(1)和(2)，需要理解在含义上的一个不同：对于数学上随机变量这个概念，真值、准确值是没有意义的；有意义的“真值”或“准确值”是在物理测量上。因为“多个可能值”是数学上随机变量概念的随机性特征，这些多个可能值，虽值不同、但地位平等，并不含有哪个值是真值、准确值；而对于一个具体的物理测量对象，它的真值是客观存在的，尽管由于测量本身存在误差，这个准确的真值是不可知的。在物理测量上如何得到“尽可能准确的值”在2.3.3节中将予提及。

(6) 补充上述(1)，在NWP领域的实际应用中，随机变量的一体二面之可能性和确定性，使得对NWP的认识更接近真实而有了更科学合理的含义：其一，对于当今NWP两个新发展之“集合预报”和“观测预报交互式预报”，不确定性是这两者的一个共同实质。两者所不同的是对应在效果上。较之于传统的确定性预报，集合预报更直接地体现在NWP的输出结果更科学准确，而观测预报交互式的预报更重要地体现在NWP的做法方式更合理经济。其二，概率思维是从确定论转到了概率论、从过去概念下的确定的值转到了新的概念下确定的概率，这体现在气象预报和服务上从天气的确定性预报(定时、定点、定量)转到天气的概率预报以及概率思维意识下天气灾害的风险评估。可见，随机变量的一体二面体现了更科学地正确认识和对待事物的科学素养。由于这些超出了本文的内容，只顺便提及。

2 大气资料同化的理想化分析方程 2.1 大气运动状态和大气资料同化

一般地，按照流体连续介质假设，大气状态是某一时刻的一个三维连续统一体，但是随研究对象的不同，它可以考虑为某一物理空间点的温度，或是二维物理空间的高度场。对于NWP，模式大气状态变量代表大气状态；具体地，一个格点模式，其模式大气状态变量是离散化的所有三维物理空间格点上所有模式变量构成的一个整体，数学上表示为一个N_x维空间(这里称作模式大气状态变量空间)的向量X，维数N_x等于所有格点数与模式变量个数的乘积。

如引言中所述，大气资料同化是使用所有可用的信息尽可能准确地确定大气状态，在实施上是由观测和NWP短期预报场产生分析场，如式(1)。观测和预报场都是一个向量。预报场是模式大气状态变量空间的一个向量，用作先验背景场，表示为X_b，其维数是N_x。类似地，定义一个观测变量空间，其维数N_y等于所用观测的数据个数，则观测是该观测变量空间的一个向量，表示为Y_o。分析场是待估计的尽可能准确的大气状态，所以它是与先验估计X_b同属于模式大气状态变量空间的一个向量，表示为X，其维数也是N_x。

2.2 基于随机向量概念下的大气状态的最优估计

多维随机变量是若干个随机变量构成的一个整体。一般地，设N个随机变量X₁, X₂, …, X_N，则X=(X₁, X₂, …, X_N)^T为一个N维随机变量，这里用黑体大写字母表示。把它看作N维空间的一个点，就形成在这个N维空间的一个随机向量，通常表示为列向量；它的每个分量是一个随机变量。随机变量是随机向量的特例，是一维随机向量。

如1.2节(3)所述，把观测和背景场视为已知的测量对象，把待估计的大气状态变量视为尚未知的测量对象，于是Y_o、X_b和X可以看作为随机向量。作为随机向量的一种可能值，y_o表示观测向量的实际观测值，x和x_b分别表示模式大气状态随机向量和其先验背景场随机向量的值。对应着x的最优值是分析场的值x_a。

于是，在随机向量的概念下，可以清晰地理解：“由观测和NWP短期预报场产生分析场”就成为“由两个已知随机向量的信息求解一个未知随机向量的最优值”的估计问题。

我们知道，一个随机变量由它的概率分布(概率函数或概率密度函数)完整地描述。于是，根据估计理论的贝叶斯定理(Bayes’ theorem)，可以导出大气状态最优估计的理想化分析方程(Lorenc, 1986；Hamill, 2006)。

2.3 贝叶斯定理和理想化分析方程

贝叶斯定理指出：随机事件A在随机事件B发生的条件下的后验概率，与事件A的先验概率和事件B在事件A发生的条件下的后验概率的乘积成正比：P(A|B)∝P(B|A)P(A)。

2.3.1 待估计的大气状态随机向量X的概率是已知当前和过去所有观测的条件下的后验概率

如引言中所述，引入NWP 短期预报场作为先验背景场，使得模式预报和同化分析构成循环。不妨假定：最初的背景场是通过观测平均或是基于观测的气候值给定，那么可以理解：前期的模式预报场是来自前期之前的所有观测，于是分析场是来自当前和过去的所有观测。

设当前时刻表示为t时刻，X^t表示待估计的该时刻模式大气状态随机向量，x^t是它的一个可能取值，X^t=x^t代表随机事件A；Z^t表示已知的当前时刻t的观测向量$ \boldsymbol{Y} _o^t$以及之前时刻的所有观测向量，即$ {\boldsymbol{Z}^t} = \left[ {\boldsymbol{Y}_o^t, {\boldsymbol{Z}^{t - 1}}} \right]$(其中Z^t-1=[Y^t-1, …, Y^o])，$ {\boldsymbol{z}^t} = \left[ {\boldsymbol{y}_o^t, {\boldsymbol{z}^{t - 1}}} \right]$是已知的观测值，Z^t=z^t代表随机事件B。根据贝叶斯定理，则有：

$ \begin{array}{l} P\left( {{\mathit{\boldsymbol{X}}^t} = {\mathit{\boldsymbol{x}}^t}|{\mathit{\boldsymbol{Z}}^t} = {\mathit{\boldsymbol{z}}^t}} \right) \propto \\ P\left( {{\mathit{\boldsymbol{Z}}^t} = {\mathit{\boldsymbol{z}}^t}|{\mathit{\boldsymbol{X}}^t} = {\mathit{\boldsymbol{x}}^t}} \right)P\left( {{\mathit{\boldsymbol{X}}^t} = \mathit{\boldsymbol{x}}t} \right) \end{array} $

(2)

这里，P(X^t=x^t|Z^t=z^t)就是待估计的模式大气状态随机向量X的概率密度函数，它表述为已知当前和过去的所有观测的条件下的概率密度函数；这在含义上对应着：分析场是来自当前和过去的所有观测。

顺便说明，贝叶斯定理的表示通常在式(2)右端项的分母中还有一个标准化的常数；为简单起见，这里分母的这一项省去，且假定研发者在程序编码时将确保概率密度的积分等于1。

2.3.2 待估计的随机向量X的概率是已知先验背景场和当前观测的条件下的后验概率

式(1)右端项P(X^t=x^t)表示无条件下的当前时刻真实大气状态的概率密度函数；至少由于物理测量的误差，真实大气状态是不可知的，所以它实际上是不可能知道的，因此需要进一步地推导待估计的概率密度函数P(X^t=x^t|Z^t=z^t)。

假定不同时刻的观测误差是独立的(这是实际应用中对于独立观测的通常做法，但可能不符合卫星观测，其仪器偏差可能很难去除；此外，观测的代表性误差(Daley，1993)可能是依赖于天气演变的，在时间上相关)，有：

$ \begin{array}{l} P\left({{\mathit{\boldsymbol{Z}}^t} = {\mathit{\boldsymbol{z}}^t}|{\mathit{\boldsymbol{X}}^\mathit{t}} = {\mathit{\boldsymbol{x}}^\mathit{t}}} \right) = P\left({\mathit{\boldsymbol{Y}}_o^t = \mathit{\boldsymbol{y}}_o^y|{\mathit{\boldsymbol{X}}^t} = {\mathit{\boldsymbol{x}}^t}} \right) \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;P\left({{\mathit{\boldsymbol{Z}}^{t - 1}} = {\mathit{\boldsymbol{z}}^{t - 1}}|{\mathit{\boldsymbol{X}}^t} = {\mathit{\boldsymbol{x}}^t}} \right) \end{array} $

代入式(2)，有：

$ \begin{array}{l} P\left({{\mathit{\boldsymbol{X}}^t} = {\mathit{\boldsymbol{x}}^t}|{\mathit{\boldsymbol{Z}}^t} = {\mathit{\boldsymbol{z}}^t}} \right) \propto P\left({\mathit{\boldsymbol{Y}}_o^t = \mathit{\boldsymbol{y}}_o^t|{\mathit{\boldsymbol{X}}^t} = {\mathit{\boldsymbol{x}}^t}} \right) \times \\ \;\;\;\;\;\;\;\;\;\;P\left({{\mathit{\boldsymbol{Z}}^{t - 1}} = {\mathit{\boldsymbol{z}}^{t - 1}}|{\mathit{\boldsymbol{X}}^t} = {\mathit{\boldsymbol{x}}^t}} \right)P\left({{\mathit{\boldsymbol{X}}^{\rm{t}}} = {\mathit{\boldsymbol{x}}^t}} \right) \end{array} $

(3)

再次根据贝叶斯定理，式(3)右端中的后二项乘积写为：

$ \begin{array}{l} P\left({{\mathit{\boldsymbol{Z}}^{t - 1}} = {\mathit{\boldsymbol{z}}^{t - 1}}|{\mathit{\boldsymbol{X}}^t} = {\mathit{\boldsymbol{x}}^t}} \right)P\left({{\mathit{\boldsymbol{X}}^t} = {\mathit{\boldsymbol{x}}^t}} \right) \propto \\ \;\;\;\;\;\;\;\;\;\;\;P\left({{\mathit{\boldsymbol{X}}^t} = {\mathit{\boldsymbol{x}}^t}|{\mathit{\boldsymbol{Z}}^{t - 1}} = {\mathit{\boldsymbol{z}}^{t - 1}}} \right) \end{array} $

于是，式(3)简化成：

$ \begin{array}{l} P\left({{\mathit{\boldsymbol{X}}^t} = {\mathit{\boldsymbol{x}}^t}|{\mathit{\boldsymbol{Z}}^t} = {\mathit{\boldsymbol{z}}^t}} \right) \propto P\left({\mathit{\boldsymbol{Y}}_o^t = \mathit{\boldsymbol{y}}_o^t|{\mathit{\boldsymbol{X}}^t} = {\mathit{\boldsymbol{x}}^t}} \right) \times \\ \;\;\;\;\;\;\;\;\;\;P\left({{\mathit{\boldsymbol{X}}^t} = {\mathit{\boldsymbol{x}}^t}|{\mathit{\boldsymbol{Z}}^{t - 1}} = {\mathit{\boldsymbol{z}}^{t - 1}}} \right) \end{array} $

(4)

式中右端后一项P(X^t=x^t|Z^t-1=z^t-1)就是“先验”概率密度函数，它表述为给定直至t-1时刻所有过去观测的条件下的大气状态在t时刻的概率密度函数；这在含义上对应着：分析时刻t的先验背景场(即前期的模式预报场)是来自之前的所有观测。右端前一项$ P\left({\mathit{\boldsymbol{Y}}_o^t = \mathit{\boldsymbol{y}}_o^t|{\mathit{\boldsymbol{X}}^t} = {\mathit{\boldsymbol{x}}^t}} \right)$是当前时刻观测的概率密度函数。顺便指出，由于先验背景场来自之前的所有观测，所以前面的“不同时刻的观测误差是独立的”这个假定，也意味着假定先验背景场误差和当前观测误差是不相关的。

在一般意义上，式(4)简单而优美。它表达了一个递归关系：“后验”概率密度函数P(X^t=x^t|Z^t=z^t)是其“先验”概率密度函数P(X^t=x^t|Z^t-1=z^t-1)和当前时刻观测随机向量的概率密度函数$ P\left({\mathit{\boldsymbol{Y}}_o^t = \mathit{\boldsymbol{y}}_o^t|{\mathit{\boldsymbol{X}}^t} = {\mathit{\boldsymbol{x}}^t}} \right)$的乘积。

2.3.3 大气状态最优估计的理想化分析方程

后验概率密度函数P(X^t=x^t|Z^t=z^t)式包含了待估计的当前时刻t的模式大气状态随机向量X的所有信息。于是，这样判定：当前时刻t的模式大气状态随机向量X的最优值x_a是X的均值，即：

$ {x_a} = \int {{x^t}} P\left({{\mathit{\boldsymbol{X}}^t} = {\mathit{\boldsymbol{x}}^t}|{\mathit{\boldsymbol{Z}}^t} = {\mathit{\boldsymbol{z}}^t}} \right){\rm{d}}{\mathit{\boldsymbol{x}}^t} $

(5)

或是X的众数：

$ {\mathit{\boldsymbol{x}}_a} = {\mathit{\boldsymbol{x}}^t}, 此时满足P\left({{\mathit{\boldsymbol{X}}^t} = {\mathit{\boldsymbol{x}}^t}|{\mathit{\boldsymbol{Z}}^t} = {\mathit{\boldsymbol{z}}^t}} \right)最大 $

(6)

由此得到的x_a就是分析场的值，它是使用观测信息和先验背景信息得到的大气状态的最优估计。如果概率密度函数是正态分布的，这时均值和众数是相等的。

式(5)或(6)就是大气资料同化基于估计理论的理想化分析方程。可以理解，对于任一随机变量ξ，则E(ξ-c)²=E(ξ²)-2cE(ξ)+c²；它是c的二次函数，在c=E(ξ)处到达最小(陈希孺，2009)。可见，式(5)的均值实质上表示了大气状态随机向量X的最小方差估计。根据众数的定义，式(6)的众数表示了大气状态随机向量X的最大似然估计。因此，式(5)和(6)也是x_a的最优标准，分别是分析方差最小和似然概率最大。

自OI开始的基本主流方法正是由这二个最优标准得到x_a的数学公式形式，如：OI和EnKF都是以分析方差最小确定统计最优权重(Bergman，1979；Evensen，1994)，Var以分析似然概率最大导出所要极小化的目标泛数(Lorenc，1986)。这些同化分析方法的具体数学公式见3.1节。

为了补充上述的1.2节的(5)，顺便指出，在物理测量上，对于一个具体测量对象，正是基于式(5)和(6)，为了得到尽可能准确的值，就进行多次的测量；把该测量对象看作一个随机变量，每次测量结果都认为是真值的一个可能值；这样，如果假定测量误差是正态的，那么当测量次数趋于无穷大时，用它的均值或众数作为这个随机变量的最佳估计，就是物理测量对象的准确的真值。当然，由于在实际中不可能进行无穷次测量，所以这个真值是得不到的和不可知的。

3 大气资料同化历史进程中的可用信息

为了直接看出同化发展进程中的可用信息，首先力图简明地阐述基本主流同化方法的数学公式及其由来。

3.1 基本主流方法的数学公式

Panofsky(1949)提出的多项式拟合方法的思想是：假定气象变量和三维空间坐标存在着解析函数的客观联系，那么就能无需人工干涉地将代码报文的观测转变为预报模式的初值条件，即完成客观分析。于是分析的问题是找到一个解析函数p(x, y, z)来表示气象变量p的三维空间坐标(x, y, z)的分布。具体实施中，客观分析在二维等压面上、且限定在相对小的区域，这样，二维的三次多项式函数在该区域可以令人满意地拟合观测。其数学公式形式表示为：

$ \begin{array}{l} p\left({x, y} \right) = \sum\limits_{i, j} {{a_{ij}}{x^i}{y^i}, \left({i + j \le 3} \right)} \\ \;\;\;\;\;\;\;\;\;\;\;\;系数{a_{ij}}有10个。\end{array} $

(7)

由区域内M(>10)个观测p_o(x^k, y^k)(k=1, 2, …, M)，按照最小二乘原则确定其10个系数a_ij。

SCM引入了初猜场x_g；第一次初猜场由预报、气候场、或是两者混合提供。Cressman(1959)SCM是一个单点分析方案，其要点是经过数步“扫描”、逐个格点地根据影响半径N(随逐步扫描而减少)范围内的n个观测资料逐步地订正该格点的初猜场来得到分析值x_a。一般地，对于某步扫描下的任一格点i，其数学公式形式表示为：

$ \begin{array}{l} {x_a}\left(i \right) = {x_g}\left(i \right) + \frac{{\sum\limits_{k = 1}^n {\left\{ {W\left({k, i} \right)\left[ {{y_o}\left(k \right) - {x_g}\left(k \right)} \right]} \right\}} }}{{\sum\limits_{k = 1}^n {W\left({k, i} \right)} }}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{x_g}\left(k \right) = H\left({{x_g}} \right) \end{array} $

(8)

式中，H是将初猜场变量转换为观测变量的观测算子，这里只是空间插值算子。权重因子W由经验给定：d<N时W(k, i)=$ \frac{{{N^2} - {d^2}\left({k, i} \right)}}{{{N^2} + {d^2}\left({k, i} \right)}}$；d≥N时W=0；其中d是测站k和格点i之间的距离。

OI实质也是属于单点分析方案，其公式形式和SCM的式(8)相近。OI和SCM的主要不同是，OI引入了先验信息的背景场，在背景误差和观测误差不相关且都无偏的假定下，以分析误差方差最小为标准确定统计最优的权重，由此得到分析值x_a，而无需逐步的迭代过程。对一个受n个观测值影响的具体格点i，其数学公式形式表示为(Bergman，1979；Kalnay，2005)：

$ \begin{array}{l} {x_a}\left(i \right) = {x_b}\left(i \right) + \sum\limits_{k = 1}^n {\left\{ {W\left({k, i} \right)\left[ {{y_o}\left(k \right) - {x_b}\left(k \right)} \right]} \right\}} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{x_b}\left(k \right) = H\left({{\mathit{\boldsymbol{x}}_{\rm{b}}}} \right) \end{array} $

对于一个有限空间体积内的格点和变量，写成向量和矩阵形式(Lorenc，1981；Kalnay，2005)为：

$ {\mathit{\boldsymbol{x}}_\mathit{a}} = {\mathit{\boldsymbol{x}}_b} + \mathit{\boldsymbol{W}}\left[ {{\mathit{\boldsymbol{y}}_o} - H\left({{\mathit{\boldsymbol{x}}_b}} \right)} \right] $

(9)

它是引言中的式(1)，此时W是最优权重矩阵，也称为增益矩阵K，满足W=BH^T(HBH^T+R)^-1，其中B是背景场误差协方差，R是观测误差协方差，H是观测算子H的切线性算子；上标T和-1分别表示矩阵的转置和矩阵的逆。

Var是在背景场误差和观测误差不相关和都无偏，且满足正态分布的假设下：

$ p\left({\mathit{\boldsymbol{y}}_o^t|{\mathit{\boldsymbol{x}}^t}} \right) \propto \exp \left\{ {\frac{1}{2}{{\left[ {H\left({{\mathit{\boldsymbol{x}}^t}} \right) - {\mathit{\boldsymbol{y}}_o}} \right]}^{\rm{T}}}{\mathit{\boldsymbol{R}}^{ - 1}}\left[ {H\left({{\mathit{\boldsymbol{x}}^t}} \right) - {\mathit{\boldsymbol{y}}_o}} \right]} \right\} $

式中，R包含观测仪器误差、观测的代表性误差和观测算子误差的协方差。

$ p\left({{\mathit{\boldsymbol{X}}^t}|{\mathit{\boldsymbol{z}}^{t - 1}}} \right) \propto \exp {\left[ { - \frac{1}{2}\left({{\mathit{\boldsymbol{x}}^t} - \mathit{\boldsymbol{x}}_b^t} \right)} \right]^{\rm{T}}}{\mathit{\boldsymbol{B}}^{ - 1}}\left. {\left({{\mathit{\boldsymbol{x}}^t} - \mathit{\boldsymbol{x}}_b^t} \right)} \right] $

便可以由式(4)得到：

$ \begin{array}{l} P\left({{\mathit{\boldsymbol{X}}^t} = {\mathit{\boldsymbol{x}}^t}|{\mathit{\boldsymbol{Z}}^t} = {\mathit{\boldsymbol{z}}^t}} \right) \propto \\ \exp \left\{ {\frac{1}{2}{{\left[ {H\left({{\mathit{\boldsymbol{x}}^t}} \right) - {\mathit{\boldsymbol{y}}_o}} \right]}^{\rm{T}}}{\mathit{\boldsymbol{R}}^{ - 1}}\left[ {H\left({{\mathit{\boldsymbol{x}}^t}} \right) - {\mathit{\boldsymbol{y}}_o}} \right] + {{\left({{\mathit{\boldsymbol{x}}^t} - \mathit{\boldsymbol{x}}_b^t} \right)}^{\rm{T}}}{\mathit{\boldsymbol{B}}^{ - 1}}\left({{\mathit{\boldsymbol{x}}^t} - \mathit{\boldsymbol{x}}_b^t} \right)} \right\} \end{array} $

之后，由式(6)的分析似然概率P(X^t=x^t|Z^t=z^t)最大的解，便转化为等价于一个目标泛数J(x)最小的解，得到分析场的值x_a(Lorenc，1986)；它表示为以下数学公式形式(省去表示当前时刻t的上标)：

$ J\left({{\mathit{\boldsymbol{x}}_a}} \right) = \mathop {\min }\limits_x \mathit{J}\left(\mathit{\boldsymbol{x}} \right) $

(10)

式中，$ J\left(\mathit{\boldsymbol{x}} \right) = \frac{1}{2}{\left({\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_b}} \right)^{\rm{T}}}{\mathit{\boldsymbol{B}}^{ - 1}}\left({\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_b}} \right) + \frac{1}{2}{\left[ {H\left(\mathit{\boldsymbol{x}} \right) - {\mathit{\boldsymbol{y}}_o}} \right]^{\rm{T}}}{\mathit{\boldsymbol{R}}^{ - 1}}\left[ {H\left(\mathit{\boldsymbol{x}} \right) - {\mathit{\boldsymbol{y}}_o}} \right]$。

当观测算子H在x_b附近做线性化近似的情况下，它能够得到显式解；又在不考虑观测算子误差的情况下，该显式解也是式(9)，即在数学公式形式上和OI相同。

和OI相同，EnKF的数学公式形式也是式(9)(Evensen，1994)，且它的增益矩阵也是以分析方差最小为标准的统计最优权重。

对于正态分布的随机变量，以方差最小为最优标准的均值和以似然概率最大为最优标准的众数是相等的，由此可以理解：OI、Var和EnKF这三种方法的x_a有着相同的数学公式形式。需要指出，除了最优标准的不同，这三种方法的主要不同体现在B矩阵的处理和x_a的求解方式上；它们不仅体现了这些方法各自的显著特征，也和3.2节要阐述的各自可用信息有直接联系，在此做些阐述。

在背景场误差协方差矩阵B的处理上，由于B的巨大维数和性状而导致的数值求解困难，OI和Var把B当作不随时间变化的静态常量；EnKF使用短期预报集合来估计动态演变的B，不过，具体实施中并不是需要显式地存储访问和计算B，而是使用短期预报集合分开计算增益矩阵中的HBH^T和BH^T。

在求解方式上，首先，OI和Var一次性地使用一个同化时间窗内各时次的所有观测资料进行同化分析，EnKF是逐时次序贯地同化时间序列的观测资料。之后的求解实现上，OI将最优权重的等式重写为W(HBH^T+R)=BH^T；通过求解该线性方程组先得到权重W，再代入式(9)来得到x_a；Var是通过式(10)泛函极小的全局求解来得到x_a；EnKF是通过直接计算增益矩阵K、由式(9)来得到x_a。

基于OI、Var和EnKF经典文献的这三种方法的基本思路、做法要点和显著特征的理解评论在另一篇论文中综述。

3.2 可用信息的不断扩展

“随机变量”概念下的数据包含值、不确定性和协方差。由上述各方法数学公式的输入信息，能够清晰地看到大气资料同化的历史进程中可用信息在内涵和种类上的不断扩展。

对于多项式函数拟合方法，式(7)的输入信息只有观测，且它的可用信息只包含观测的“值”p_o(x^k, y^k)，用来确定系数a_ij。

对于SCM，式(8)的输入信息除了观测还有第一次初猜场，而它们的可用信息也只包含“值”y_o和x_g。第一次初猜场的值x_g由NWP预报场、气候场、或是两者混合提供；因为预报场是来自NWP预报模式，它意味着源自大气运动演变物理规律的信息；气候场是来自气候平均，它意味着源自大气运动之统计性质的信息。

对于OI、Var和EnKF，它们的x_a有着相同的数学公式形式式(9)，其输入信息包括观测和背景场，此时它们的可用信息不仅包含“值”y_o和x_b，还有相应的“协方差”B和R；由于任两个随机变量X和Y的协方差矩阵可以写成各自标准差和它们相关系数矩阵的乘积：Cov(X, Y)=σ_XCorr(X, Y)σ_Y，所以B和R包含了其数据的精度和变量之间交叉相关。

所用观测资料的种类上，由于OI得到x_a的过程中需要先求解关于权重的线性方程组，限制了观测算子H不能太复杂(一般是线性插值)，这限制了OI能被同化的观测资料种类；此外，还由于OI限于实质上的单点分析方式和求解线性方程组的规模，因此存在对于某一格点或一组格点分析时观测资料的区域选择问题。Var是通过泛函极小的全局求解来得到x_a；由于是通过泛函极小的方式，因此能够引入复杂的非线性观测算子，使得Var能够直接同化如卫星、雷达等大量的非常规、非模式变量的观测资料(如辐射率、GPS掩星、回波反射率、径向风等)；此外，由于是全局求解，因此不存在观测资料的区域选择问题。

对于背景场误差协方差，OI和Var只考虑了静态不变的B矩阵；EnKF使用短期预报集合来考虑动态演变的B矩阵。静态B矩阵的实质是包含真实B矩阵平均意义上的一些非常重要的性质，其变量之间交叉相关内含了大气运动之准平衡特征的动力统计性质(如运动场和质量场之间的地转平衡、质量场之气压和位势之间的静力平衡等约束关系)；对于某一时刻的大气状态，这些统计性质只是一种平均意义上的近似。由于短期预报集合来自NWP的集合预报，所以EnKF的动态B矩阵内含了模式大气运动(包含非线性)的结构特征和动力演变，能够更接近真实B矩阵。

表 1示出上述的大气资料同化历史进程中可用信息的不断扩展。

可以看到，这个扩展在历史进程中是上下承接和循序渐进的，具体地体现了同化发展史的一个内在发展逻辑，在认识上是非常清晰、简单的。

顺便指出，在具体实施中还有着多方面信息的综合使用，如：Panofsky(1949)多项式拟合方法为得到最可信的拟合还利用天气动力学知识和具体观测资料特征来考虑平滑，Cressman(1959)SCM还利用地转关系来使用风观测对高度场分析的订正，等等。

4 结 论

是什么使得大气资料同化有了理论上的数学基础和数学表示形式而发展成为一门科学工程学科？为此本文阐明随机变量是大气资料同化的根本概念。其根本性可以体现在两个方面：

(1) 原理上，大气资料同化的科学实质建立在随机变量的概念上，因为以随机变量为研究对象的概率论与数理统计学得以成为资料同化的数学基础，这使得资料同化有了基于估计理论的最优标准及其数学形式；

(2) 实施上，基于随机变量概念下的数据可以理解资料同化所使用的可用信息在其历史进程中的不断扩展，它在认识上是非常简单的，清晰地展示了上下承接的纵向联系和循序渐进的发展创新逻辑。