Pu-Kalnay2018_NWP

数值天气预报基础：模型、数值方法和数据同化

赵晓普和尤金尼亚·卡尔奈

内容

介绍：基本概念和历史概述

NWP 模型和数值方法

2.1

基本方程式

2.2

数值天气预报模型的框架

2.3

全球和区域模型

2.4

物理参数化

2.5

陆面和海洋模型以及数值天气预报中的耦合数值模型

数据同化

3.1

最小二乘法理论

3.2

同化方法

4 最近的发展与挑战

参考文献

摘要

数值天气预报已成为全球天气预报最重要的工具。本章概述了数值天气预报的基本原理，包括模型的数值框架、数值方法、物理参数化和数据同化。介绍并讨论了历史变革、近期发展和未来方向。

Z. Pu (*) 大气科学系，犹他大学，美国犹他州盐湖城，电子邮件：Zhaoxia.Pu@utah.edu

E. Kalnay 大气与海洋科学系，马里兰大学，大学公园，马里兰州，美国电子邮件：ekalnay@atmos.umd.edu

# 施普林格出版社德国，施普林格自然的一部分 2018

Q. Duan 等（主编），《水文气象集合预报手册》

https://doi.org/10.1007/978-3-642-40457-3_11-1

引用此条目为：Pu Z., Kalnay E. (2018) 数值天气预报基础：模型、数值方法和数据同化。见：Duan Q., Pappenberger F., Thielen J., Wood A., Cloke H., Schaake J. (主编) 水文气象集合预报手册。施普林格，柏林，海德堡

关键词

数值天气预报 · 数值方法 · 物理参数化 · 数据同化

介绍：基本概念和历史概述

随着计算机技术的进步，数值天气预报（NWP）已成为天气预报的核心组成部分。例如，在美国，日常天气预报始于位于华盛顿特区的国家海洋和大气管理局（NOAA）的一台超级计算机。在欧洲，世界上最大的数值天气预报中心——欧洲中期天气预报中心（ECMWF）为所有欧盟成员国提供先进的天气指导。在全球范围内，大多数国家将 NWP 作为其操作天气预报的关键指导。

数值天气预报（NWP）的基本概念是求解一组控制大气运动和演变的偏微分方程（PDEs）（Kalnay 2003）。正如在第 2.1.2 节中将要描述的，这组偏微分方程描述了基本的守恒定律，包括动量、质量、能量和水汽的守恒。为了预测未来的大气状态，我们必须向前积分这组方程。因此，数值天气预报是一个初值问题：给定当前的大气条件（初始条件），我们积分这组偏微分方程以获得未来的大气状态。

这个初值问题在 20 世纪初就被定义了（例如，Bjerknes 1904）。Bjerknes（1904）指出，气象学中的最终问题是天气预报（预测未来的气象条件），并概述了一种解决方法。根据他的方法，成功预测未来大气状态必须满足两个条件：

I. 目前的气象条件必须尽可能准确地描述。 II. 必须了解内在规律，依据这些规律，后续状态从前一状态发展而来。

他概述了一个程序，该程序被细分为三个部分问题或组成部分：

1. 观察组件 2. 诊断或分析组件 3. 预后组件组件 1 和 2 与当前状态的特征化（状态 I）相关，而组件 3 与状态 II 相关。

在今天的术语中，Bjerknes 的方法被称为确定性，因为预测被假定为完全由当前状态决定。此外，这里的组件 1 指的是全球观测系统，尽管

Z. Pu 和 E. Kalnay

这项工作直到 1970 年代才建立。组件 2 现在被称为数据同化，它结合了观测信息和短期天气预报，以形成最佳的初始条件。组件 3 涉及使用数值方法求解偏微分方程。

在比耶克内斯之后，理查德森（1922 年）首次尝试手动进行数值天气预报。他使用了完整的原始方程和有限差分方案。他将感兴趣的区域划分为单元，就像棋盘上的方格一样。他通过手动插值从天气图中读取大气条件。尽管他的方法论无可挑剔，但他在 6 小时内获得的 145 毫巴的地面气压变化预报却失败了！理查德森预报的失败使数值天气预报的概念在理论世界中停滞了许多年。直到 1950 年代，数值天气预报才再次被尝试。查尼使用条流势涡度方程成功进行了第一次数值天气预报（查尼等，1950 年）。具体来说，在 1920 年代到 1950 年代之间，在以下领域取得了显著进展：

•

动态气象学。大气运动包括多个时间和空间尺度。因此，尺度分析方法可以用来简化数值天气预报方程，基于感兴趣的天气预报运动尺度。根据尺度分析，定义了罗斯比数以验证中纬度天气系统大气中的地转流。准地转理论被推导出来以解释中纬度天气系统流的环流。还推导了斜压不稳定性的条件。动态气象学对数值天气预报的主要好处在于根据目标预报尺度解决简化方程，而不必处理完整的原始方程（详见 Holton 2004 和 Kalnay 2003）。

•

数值分析的进展。由于气象预报方程的解析解不存在，因此必须使用数值方法在离散网格空间中进行数值积分。Courant-Friedrichs-Lewy (CFL) 条件设定了网格空间大小和时间步长之间的基本要求，以保持计算稳定性。对非线性计算稳定性的更好理解也有助于设计可以用数值方法准确高效地解决偏微分方程的方案（详见 Kalnay 2003）。

•

大气观测。无线电探空仪的发明使得探测高层大气的条件成为可能。在第二次世界大战期间，许多国家开始建立无线电探空仪网络。这对提高数值天气预报所需初始条件的准确性有很大帮助。

•

电子计算机的发明。电子计算机的发明极大地提高了科学计算的效率。因此，计算机帮助科学家在实际操作中实施数值天气预报，而运营数值天气预报中心是超级计算机的主要用户。

在 1950 年，当第一个计算机预报使用世界上第一台电子计算机电子数值积分器和计算器（ENIAC）生成时，数值天气预报（NWP）变得在操作上可行。

自 1950 年代以来，数值天气预报（NWP）领域持续快速发展

数值天气预报基础：模型、数值方法和数据……

• 1970 年代全球观测系统的成熟（http://www.wmo.int/pages/prog/www/OSY/GOS.html）以及卫星和雷达观测的使用

•

数据同化技术的进展（见 Daley 1991；Kalnay 2003；Evensen 1994，详细信息见“数据同化”部分）

•

计算机技术的重大进展和数值方法的发展，特别是全球谱模型、半拉格朗日模型和高分辨率区域模型的发展（见 Robert 1982；Williamson 2007；Lynch 2008）

•

物理参数化的快速发展（见 Arakawa 1997, 2004; Stensrud 2007）

此外，近年来发生了一些新进展：

•

集合预报：与单一的确定性预报不同，集合预报已成为今天操作性数值天气预报的主流方法（见 Kalnay 2003；Bauer 等，2015）。

•

卫星和雷达数据的高级数据同化方法（例如，Lorenc 1986；Daley 1991；Kalnay 2003；Houtekamer 和 Zhang 2016）。

•

耦合的气候-海洋-陆地模型（例如，Hodur 1997；Chen 和 Dudhia 2001；Ek 等 2003；Tolman 2014）。

今天，数值天气预报（NWP）已成为研究和操作环境中的一门多学科科学。许多国家拥有自己的 NWP 系统用于日常操作预测，包括全球和区域模型系统。此外，NWP 计算模型已成为研究和教育的有用工具。由美国机构主导的社区天气和气候模型，例如国家大气研究中心（NCAR）、宾州州立大学等，使这些天气和气候研究及预测模型可供全球许多研究机构和大学使用。

随着过去 60 年计算机性能和计算机科学的快速发展，数值天气预报（NWP）技能稳步提高。图 1 显示了 1981 年至 2014 年间欧洲中期天气预报中心（ECMWF）的预报技能提升（注意其他中心也遵循类似的改进趋势）。显然，NWP 产品在 5 天范围内相当可靠，在 7 天范围内也很有用。自 1990 年代末以来，利用卫星观测的先进应用所取得的改进是显著的，特别是在南半球。

NWP 模型和数值方法

2.1

基本方程式

有一整套七个方程和七个未知数，控制着大气的演变：牛顿第二定律或动量守恒（针对三个速度分量的三个方程），连续性方程或

Z. Pu 和 E. Kalnay

质量守恒，理想气体状态方程，热力学第一定律或能量守恒，以及水质量的守恒方程。

在笛卡尔坐标系中，这组方程可以写成如下形式：

d V

⇀

dt ¼ "

α ∇

⇀

" ∇

⇀

þ F

⇀

2Ω$ V

⇀

(1)

@ρ

@t ¼ " ∇

⇀

ρ V

⇀

(2)

pα

(3)

dt "

(4)

@ρq

@t ¼ " ∇

⇀

ρ V

⇀

ρ E

" C

(5)

500hPa 位势高度异常相关性 12 个月滑动平均 100

（居中于窗口中间）

1981 1983 1985 1987 1989 1991 1993 1995 1997 1999 2001 2003 2005 2007 2009 2011 2013

第 7 天 NHem

第 7 天 SHem

第 10 天 NHem

第 10 天 SHem

第三天 NHem

第三天 SHem

第五天 NHem

第五天 SHem

图 1 1981 年 1 月至今 500 hPa 地势高度预报与实际分析的异常相关性年运行平均时间序列。特定月份的值是该月、前 5 个月和后 6 个月的平均值。显示了 3、5、7 和 10 天的预报提前期，分数分别在北半球（粗线）和南半球（细线）外热带的平均值。阴影显示了在所示预报范围内两个半球之间分数的差异。（改编自 Simmons 和 Hollingsworth 2002。验证遵循更新的 WMO/CBS 指南，如《全球数据处理和预报系统手册》第 1 卷第 II 部分附件 II.7 表 F（2010 年版 - 2012 年更新）所规定；异常是相对于 ERA-Interim 基础气候计算的（感谢 ECMWF 网站 http://www.ecmwf.int；另见 Bauer 等人 2015 年）。

数值天气预报基础：模型、数值方法和数据……

哪里 V

¼ u,v,wð

表示空气的速度，t 是任意时间，α 是比容，ρ 是密度，p 和 T 分别是压力和温度，Φ 是地势高度，q 是水汽混合比，Q 是加热，E 和 C 分别表示蒸发和凝结，R 是气体常数，F

⇀

是摩擦力。

在球坐标中，假设 λ 和 ϕ 分别是经度和纬度，r 是地球的半径：

区域正东向

Þ ¼

r cos φ

dλ

南向正向北

Þ ¼

dφ

垂直向上

Þ ¼

由于 V

u i

v j

w kand r

þ z

;

r z；

@z.

方程（1）、（2）、（3）、（4）和（5）可以写成：

@t ¼ "

@x "

@y "

@z þ

紫外线切线 ϕ

a "

@x " 2

Ω w cos ϕ

v sin ϕ

Þ þ

(6)

@t ¼ "

@x "

@y "

@z "

正切ϕ

a "

@y " 2Ω

u sin ϕ

(7)

@t ¼ "

@x "

@y "

@z "

þ v

@z þ 2Ω

u cos ϕ

(8)

@t ¼ "

@x "

@y þ γ" γð

(9)

@ρ

@t ¼ "

@ρ

@x þ "

@ρ

@y "

@ρ

@z "

@x þ

@y þ

(10)

@t ¼ "

@x "

@y "

@z þ

(11)

pα

(12)

其中 Υ 和 Υ 分别是温度递减率和干绝热递减率。请注意，与非绝热效应 (H)、摩擦 (F)、通过相变的水的增益或损失 (Q) 相关的方程中的项必须在模型中定义。通常，上述方程组称为原始方程。这七个方程与七个未知数代表了数值天气预报的通用方程组。

Z. Pu 和 E. Kalnay

实际上，由于数值模型是为处理不同尺度的各种目的而构建的，因此预计这组方程将通过一些假设进行简化。例如，如果模型仅设计用于处理大尺度，则静水平衡方程将替代垂直运动方程。准博辛斯基或非弹性近似将用于使密度保持常数，并消除解中的高频波，以保持计算稳定性。此外，方程的形式可以根据使用的各种坐标系统（例如，笛卡尔坐标与球坐标或压力坐标），特别是垂直坐标系统（例如，西格玛坐标与压力垂直坐标）进行更改。

我们必须在大气的底部和顶部添加适当的边界条件，然后按照理查森（1922 年）所建议的积分过程来求解这些方程。

@φ

@t ¼

F φ,t

¼ φ

jþ

F φj,t

% &

Δt

(13)

2.2

数值天气预报模型的框架

2.2.1

有限差分方程 (FDEs)

上述方程组（方程 6、7、8、9、10、11 和 12）的解析解是不可行的。必须使用离散形式和数值方法来求解这些方程。因此，可以使用有限差分方程（FDEs）来寻找偏微分方程（PDEs）的近似解。

使用对流方程作为例子：

@t þ

@x ¼

(14)

我们对 x 和 t 取离散值：x = jΔx 和 t = nΔt，其中 Δx 是网格间距，Δt 是积分的时间步长。FDE 的解在离散点 (x, t) = (jΔx, nΔt) 处定义：

U jΔx, nΔt

Þ ¼

U x, t

(15)

在这里，我们用小写字母 u 表示偏微分方程（连续）的解，用大写字母 U 表示分数微分方程（离散）的解。

用于近似 PDE (14) 的 FDE 可以写成如下形式：

U j

Δt

Δx

(16)

数值天气预报基础：模型、数值方法和数据……

这被称为上游方案，如果我们假设 c > 0。请注意，这两个差异相对于点 (x, t) = (jΔx, nΔt) 都是非中心的。

由于我们使用有限差分法（FDE）来近似偏微分方程（PDE），因此应满足两个基本条件：

(i) FDE 应与 PDE 一致。

对于给定的时间 t > 0，FDE 的解应随着 Δx 收敛到 PDE 的解

0 和 Δt

为了满足这两个要求，FDE 中使用的数值方案必须尽可能准确。如果在极限Δx 中，FDE 与 PDE 一致，我们就说 FDE 与 PDE 一致。

0, Δt

0，FDE 与 PDE 重合。这要求 FDE 的解与 PDE 的解是一致的近似。PDE 和 FDE 之间的差异是离散化误差或局部（在空间和时间上）截断误差。

此外，在积分（预测）过程中保持计算稳定性非常重要。通常，在为给定网格大小 Δx 指定时间步长 Δt 时，必须满足 Courant-Friedrichs-Lewy（CFL）条件：

0 c

Δt

Δx

(17)

其中 c 是所谓的 CFL 数，c 是一个翻译速度。它通常可以通过模型中所处理的天气事件的典型最大风速来指定。

然而，这个 CFL 条件仅仅是一个必要条件，它确保了 FDE 在计算上是稳定的，以便在固定时间 t = n Δt 时，FDE 的解在 Δt 变小时保持有界

由于数值天气预报方程的非线性，即使满足 CFL 条件，计算不稳定性仍然可能发生。通常，时间步长应设置得小于满足 CFL 条件的值。然而，一些隐式数值方案允许使用较大的时间步长（例如，有限体积隐式方案）。为了实现时间和空间离散化的准确性和稳定的数值方案，计算数学领域取得了许多进展，如半拉格朗日方案、有限体积方案和高阶龙格-库塔方案（参见 Robert 1982；Durran 1999；Lin 和 Rood 1996；Lin 1997）。

2.2.2 光谱模型

除了网格空间离散化，谱模型使用 Galerkin 方法描述大气的当前和未来状态，通过一组基函数φ(x)进行空间离散化

U x,t

ð Þ ¼

k¼1

ð Þφ

ð Þ

(18)

Z. Pu 和 E. Kalnay

空间导数直接从已知的 dx 计算得出。该过程导致一组关于系数 A(t) 的常微分方程 (ODE)。例如，考虑单向平流方程：

@t þ

@x ¼

(19)

我们使用傅里叶变换对：

ξ u

ð Þ ¼

U k,t

ð Þ ¼

u x,t

ð Þ

exp "ikx

ξU

ð Þ ¼

u x,t

ð Þ ¼

2π

U k,t

ð Þ

exp ikx

(20)

其中 ξ 代表正向傅里叶变换算子，ξ 表示逆傅里叶变换算子。根据上述定义，我们可以证明：

dU k,t

ð Þ

ikU k,t

ð Þ

(21)

以便：

•

函数的时间导数的傅里叶变换等于该函数的傅里叶变换的时间导数。

•

函数的空间导数的傅里叶变换等于该函数本身的傅里叶变换乘以 ik。

我们可以使用这个，并对平流方程进行傅里叶变换，如下所示：

c ξ

(22)

这与以下内容相同：

dt þ

ickU

(23)

通过使用傅里叶变换，我们将偏微分方程转化为常微分方程。然后，我们可以将常微分方程向前积分，以找到 U(k, t) 的未来值，然后进行逆变换以找到 u(x, t)。这就是谱方法的本质。我们将实空间中的偏微分方程转换为波空间中的常微分方程，然后求解它们。

数值天气预报基础：模型、数值方法和数据……

基于谱表示的空间离散化极其准确（空间截断误差为“无限”阶），因为空间导数是通过解析计算的，而不是数值计算。鉴于这一优势，谱模型比具有相同分辨率的网格点模型更适合进行长时间范围的预报。因此，今天许多操作性全球模型都是谱模型（例如，NCEP 全球预报系统）。然而，局部强迫过程（例如，潜热释放、差异表面热通量）有时是间断的，只能在物理空间中表示。此外，当使用波的线性组合（例如，谱谐波）来表示大梯度或间断时，可能会产生伪波（吉布斯现象）。对于更高的分辨率，谱模型在计算上比网格点模型更具挑战性。此外，谱模型在质量或能量的精确度上并不保守。出于这些原因，只有少数区域性、有限区域的谱模型被开发并用于研究和操作性预报。最广泛使用的模型之一是 NCEP 区域谱模型（Juang 和 Kanamitsu 1994）。

2.2.3 网格错位方法

一旦连续偏微分方程在网格中离散化，所有模型变量都在网格中定义。即使在谱模型中，由于谱空间到网格和从网格到谱空间的变换是必要且常用的，模型变量在某种程度上也在网格空间中定义。模型变量在不同网格点上的排列成为设计数值方案时的一个考虑因素。许多数值模型采用错位网格方法，而不是将所有变量安排在同一个网格点上。

交错网格结合了位于不同几何位置的几种类型的节点，显得相当复杂。然而，交错网格允许对几个关键的偏微分方程进行自然且更准确的有限差分公式化，因此在数值模型中被广泛使用。图 2 展示了一个水平方向的交错网格示例。在垂直方向上，大多数模型采用了交错网格，例如，垂直速度在层的边界定义，而预报变量位于层的中心（图 3）。虽然可以使用非交错垂直网格，这样可以简单实现更高阶的垂直差分，但这也会在解中引入更多的计算模式。

2.2.4 边界条件

由于数值模型通常只处理宇宙的一部分，因此需要边界条件。例如，全球大气模型应给出上下边界条件。区域大气模型除了顶部和底部边界条件外，还需要侧向边界条件。

上边界条件：模型顶部的高度通常高于气象感兴趣的特征。通常，它位于平流层或更高。由于较长时间尺度的过程主导平流层，气候模型使用较高的上边界。在上边界之上，唯一感兴趣的输入是来自入射太阳辐射，通常是参数化的。表示上边界的方法有很多。例如，可以在某个指定高度用刚性盖子封闭模型。

Z. Pu 和 E. Kalnay

高度使得到达这个盖子的能量向下反射。自由表面方法将模型大气和更高的高度视为两种不同的、不混合的流体，并同样向下反射能量。由于表示上边界条件的关键问题是如何处理重力波向上和向外传递能量，因此在刚性盖和自由表面方法中都引入了一个吸收/阻尼层。这个阻尼层通常放置在模型的顶部下方，并对选定的垂直层施加扩散/阻尼算子，以减弱向上传播的能量。然而，它必须相对厚，以减轻大垂直梯度和波反射问题的发展，并可以减弱到预定义的参考状态或模型大气定义的状态。在某些模型中，还使用辐射边界条件来模拟波能量向上传播并从域的顶部逸出的效果。

图 3 垂直网格的错位（根据 Arakawa 和 Konor 1996）

图 2 Arakawa C 交错网格方法（改编自 Skamarock 等，2008 年）

数值天气预报基础：模型、数值方法和数据……

它必须相对厚，以减轻大垂直梯度和波反射问题的发展，并可以减弱到预定义的参考状态或模型大气模型所定义的状态。通常应选择一个足够高的上边界，以减轻上边界条件带来的问题。有关详细信息，请参见 Durran (1999)、Kalnay (2003) 和 Warner (2011)。

底部边界条件：数值天气预报模型的底部边界条件非常复杂，因为地表特征差异显著。因此，数值天气预报模型的底部条件通常通过参数化或采用热扩散表面模型、陆面和海洋模型或表面摩擦方案来表示（如下一节所述）。

侧边界条件（LBCs）：使用区域模型进行天气预测的原因在于希望通过提高水平分辨率来减少模型误差，而这种提高在全球模型中是无法实现的。自 1970 年代以来，操作区域模型已嵌入或“嵌套”在较粗分辨率的半球或全球模型中。例如，目前的 NCEP 北美中尺度预报系统（NAM）模型嵌套在 NCEP 全球预报系统（GFS）模型中。区域模型的嵌套需要使用从全球模型获得的更新侧边界条件。通常，如果（a）它能够传递来自“主机”模型的入射波，并在相位或幅度上没有显著变化地提供边界信息，以及（b）在出流边界处，反射波不会以显著幅度重新进入感兴趣的区域，则应满足侧边界条件。

在实践中，边界条件是根据实际情况选择的，并通过数值测试来检查其适用性。常见的侧边界条件选择包括单向和双向嵌套方案。在单向侧边界条件中，具有较粗分辨率的主模型向嵌套区域模型提供边界值的信息，但不受区域模型解的影响。在边界条件的双向交互中，即（假设更准确的）区域解也会反过来影响全局解。

此外，为了在全球模型中嵌套区域模型，许多区域模型使用嵌套域技术来实现高分辨率的模拟和预测，高分辨率域嵌套在较粗的区域模型域内。在这种情况下，侧边界条件也应以单向或双向交互的方式进行处理。此外，近年来开发了变分分辨率模型。通过使用连续拉伸的水平坐标，变分分辨率模型仅在感兴趣的区域以高分辨率进行求解。显然，采用这种方法，区域高分辨率区域的方程不需要特殊的边界条件，并且它们确实会影响较粗分辨率区域的解，因此可以视为双向交互嵌套（参见 Kalnay 2003；Warner 2011）。它们的缺点是，最小网格大小要求整个域使用短时间步长。

2.3

全球和区域模型

全球和区域模型都用于数值天气预报。全球模型通常用于中期预报（超过 3 天）的指导以及气候研究。

Z. Pu 和 E. Kalnay

模拟。在 NCEP，例如，全球模型每天运行 16 天。由于这些全球模型的水平范围是整个地球，因此通常无法以高分辨率运行。然而，随着计算能力的提高，全球模型的分辨率显著增加。例如，NCEP 全球预报系统（GFS）和 ECMWF 中期预报模型现在以近 16 公里（T1297）的水平分辨率运行，约是 20 年前水平分辨率的十倍！

为了获得更详细的预报，有必要提高分辨率，而这只能在有限的关注区域内进行。区域模型用于短期预报（通常为 1-3 天），其分辨率比全球模型高出两倍或更多。例如，在 1997 年，NCEP 全球模型运行时有 28 个垂直层，第一周的水平分辨率为 100 公里，第二周为 200 公里。区域（Eta）模型的水平分辨率为 29 公里，垂直层数为 50 层。如今，NCEP 区域北美中尺度预报系统（NAM）模型的网格间距为几公里（<10 公里），垂直层数约为 100 层。由于其更高的分辨率，区域模型能够再现小尺度现象，如锋面、阵风线，以及比全球模型更好的地形强迫。另一方面，区域模型并不是“自给自足”的，因为它们在水平域的边界处需要侧向边界条件。这些边界条件必须尽可能准确，否则区域模型的内部解会迅速恶化。因此，通常将区域模型“嵌套”在另一个分辨率较粗的模型中，该模型的预测提供不断变化的边界条件。因此，区域模型仅用于短期预测。在一定时间后，这段时间与模型的大小成正比，高分辨率初始条件中包含的信息会被边界条件的影响“冲刷掉”，区域模型仅成为粗糙模型预测在区域域中的“放大镜”。这仍然是有用的，例如，在进行长时间（季节到多年）的气候模拟时，因此往往以较粗的分辨率运行。“区域气候模型”可以提供感兴趣区域的粗糙气候模拟的更详细版本。包括英国、法国、德国在内的欧洲几个主要数值天气预报中心，以及日本、澳大利亚和加拿大也有类似的全球和区域模型，其详细信息可以在他们的网站上获取。

2.4

物理参数化

2.4.1

基本原则

NWP 模型的基本方程（例如，方程 1、2、3、4、5、6、7、8、9、10、11 和 12）包括摩擦、加热源以及蒸发和凝结过程的项。具体而言，动量方程考虑了动量的涡流通量；热力学方程包括辐射加热和冷却、显热通量以及凝结和蒸发；水汽方程包括凝结和蒸发，以及水分通量。这些物理

数值天气预报基础：模型、数值方法和数据……

以及数值模型中的湿度通量过程代表它们的贡献。因此，模型应包括表面和行星边界层过程、辐射传输和云微物理，以便代表它们的贡献。

大气运动包括广泛的时间和空间尺度。时间尺度从 1 秒到 10 秒及更长，包括小型湍流气团的生命周期到局部风暴、天气系统运动以及季节性到年际变化。空间尺度从 1 厘米到 10,000 公里，包括湍流微尺度、对流尺度、中尺度和大尺度。

由于使用数值离散化方法来求解偏微分方程，气候模型的网格分辨率总是有限的。因此，任何发生在小于网格空间尺度的过程都无法在数值模型中明确表示，尽管它们的贡献不可忽视。举个例子，这里我们将雷诺平均应用于运动方程的 u 分量（方程 6）。

假设任何变量（例如，u、v、w、T、p）可以分为可解析和不可解析的成分，即可以将所有依赖变量分别分为平均部分和湍流部分。平均值被定义为网格单元的平均值，如 Pielke（2002）所述。例如：

þ u,

T，和

þ p:

这些表达式被代入方程（6）、（7）、（8）、（9）、（10）、（11）和（12）；产生如下方程（6）的展开式：

@x ¼

þ u

Þ ¼

@ u

@x þ

@x þ u

@ u

@x þ u

0 @u @x

@x ¼

@ u

@x þ

@x þ u

@ u

@x þ u

自从

a 和 ab

¼ a

b，和

¼ 0:

因此，

Z. Pu 和 E. Kalnay

因此，

@t ¼ "

@x "

@y "

@z "

@x þ

@τ

@x þ

@τ

@y þ

@τ

哪里

@ u

@t ¼ "

@ u

@x "

@ u

@y "

@ u

@z "

@ p

@x þ

f v

" u

@x " v

@y " w

@z þ

1 ρ

@τ

@x þ

@τ

@y þ

@τ @z

@x þ

@y þ

@w @z ¼ 0

@ u

@t ¼ "

@ u

@x "

@ u

@y "

@ u

@z "

1 ρ

@ p

@x þ

f v

@uu

@x "

@uv

@y "

@uw

@z þ

1 ρ

@τ

@x þ

@τ

@y þ

@τ

¼ "ρuu,

¼ "ρuv,

¼ "ρ uw:

@ u

@t ¼ "

@ u

@x "

@ u

@y "

@ u

@z "

1 ρ

@ p

@x þ

f v

1 ρ

Þ þ

(24)

在这些方程中，五个项的第一个分量可以通过模型网格值明确表示。括号内三个项的第二个分量无法在模型网格点上明确解析，但它们的贡献不能被忽视，因为这些子网格尺度过程依赖于并反过来影响由数值模型明确解析的大尺度场和过程。因此，参数化方案是必要的，以便正确描述这些子网格尺度机制对大尺度大气流动的影响。换句话说，子网格尺度过程的集体效应必须以解析的网格尺度变量来表述。此外，预测天气参数，如 2 米温度、降水量和云量，是通过模型的物理参数化计算得出的。

总体而言，数值天气预报模型由两个主要部分组成，如图 4 所示：“模型的动力学”示意性地表示了解决的过程，而“模型物理”则是必须参数化的过程。由于所有物理过程不仅与模型的动力学相互作用，还相互影响，因此数值天气预报模型在数值上复杂且计算成本高。

2.4.2 物理参数化概述

数值模型中的物理参数化方案应设计为 (1) 表示与动力学相互作用的物理过程；以及 (2) 明确

数值天气预报基础：模型、数值方法和数据……

计算来自子网格尺度过程的贡献，这些过程被参数化为大尺度、已解析尺度的函数。数值模型中包含的常见参数化方案，以当前 ECMWF 全球模型为例，可以简要描述如下。

辐射和化学过程：辐射方案使用温度、湿度、云和气溶胶及主要痕量气体（CO、O、CH、NO 等）的月均气候值的预测值来计算短波和长波辐射通量。辐射参数化描述了辐射传输过程。通常考虑云-辐射相互作用。由于求解辐射传输方程以获得通量在计算上非常昂贵，具体取决于模型配置，通常在减少（粗糙）辐射网格和/或减少时间频率上进行完整的辐射计算。然后将结果插值回原始网格。（有关辐射过程的详细信息，请参见 Liou 1980；Stephens 1984。）

对流：湿对流方案代表深层（包括积云）、浅层和中层（升高的湿层）对流。对流方案中区分了深层和浅层对流。湿对流

图 4 大气中的物理过程及其相互作用。模型“动力学”明确计算的可解析尺度的动力过程用粗体表示。其他亚网格尺度过程则通过解析尺度场进行参数化（改编自 Arakawa 1997, 2004）。

Z. Pu 和 E. Kalnay

还解决了卷入过程和对流的昼夜变化。上升气流和下降气流的影响也被模拟（见 Arakawa 2004）。

云微物理与降水：云微物理包括所有

云过程发生在云滴和水汽颗粒的尺度上，包括云滴、雨滴、冰晶、雪花、霜冻冰粒、霰粒和冰雹，而不是在云本身的尺度上。微物理参数化旨在尽可能全面地表示微物理过程中的过程，包括凝结、聚合、蒸发、冰雪聚集、被冰冻颗粒聚合、蒸汽沉积、融化和冻结。

在大规模模型中，云和大规模降水通过一系列预报方程进行参数化，这些方程涉及云液体、云冰、雨水和雪水含量，以及子网格分数云覆盖。云方案表示由于主要的生成和破坏过程而导致的云和降水的来源和去向，包括通过积云对流的排放、凝结、沉积、蒸发、聚集、融化和冻结等过程形成的云（参见 Houze 1993；Straka 2009）。

在高分辨率模型中，特别是在允许云的区域模型中，云的微物理过程通过液体、冰和蒸气的微物理学以详细的配置和相变显式表示。在允许云的模型中，由于云被明确解析，积云对流方案可以被消除。

土壤/表面：表面参数化方案表示能量和水的表面通量及相应的地下数量。该方案应描述植被、裸土、雪和开水的不同子网格表面类型。表面能量平衡方程也包括在内。土壤层以及雪的质量和密度被表示。蒸发通量考虑了雪覆盖、湿润和干燥植被以及裸土的分数贡献。截留层收集降水和露水的水分，渗透和径流应根据土壤质地和子网格地形进行表示。可以包括碳循环，并且可以对二氧化碳的陆地-大气交换进行参数化，以响应水和能量循环中的日变化和天气变化。土壤/表面参数化可以被分类为简单的整体方案或完整的陆地表面模型。

湍流扩散和行星边界层方案：湍流扩散-

该方案表示通过亚网格尺度湍流的热量、动量和湿度的垂直交换。表面层的垂直湍流输送与其上方的处理方式不同。例如，在 2013 年 ECMWF 全球模型的表面层中，湍流通量是基于 Monin-Obukhov (MO) 相似理论使用一阶 K 扩散闭合计算的。在表面层以上，除了不稳定边界层外，处处使用 K 扩散湍流闭合，而在不稳定边界层中应用涡动扩散质量通量 (EDMF) 框架，以表示非局部边界层涡流通量。该方案以湿气守恒变量（液体静能和总水量）编写，并预测总水量方差。使用总水分布函数进行转换。

数值天气预报基础：模型、数值方法和数据……

该方案使用湿气保守变量（液态静能和总水量）编写，并预测湿气保守变量与预报云变量（液态/冰态水含量和云量）的总水方差，但仅适用于层积云的处理。对流云则由浅对流方案单独处理。（有关边界层和湍流的详细信息，请参见 Stull 1988。）

地形阻力：由于分辨率的限制，数值天气预报模型无法完全解析地形的地形特征。未解析的地形对大气流动的影响被参数化为动量的损失（阻力）。湍流扩散方案在下层大气中包含一个参数化，以表示由小尺度（<5 公里）地形引起的湍流地形阻力。此外，在稳定分层流动中，地形阻力参数化表示由于未解析地形造成的低层阻塞的影响（阻塞流动阻力）以及垂直传播的重力波对动量预算的吸收和/或反射（重力波阻力）。

非地形重力波阻力：非地形重力波阻力参数化考虑了未解决的非地形重力波的影响。这些波在自然界中由深对流、锋面扰动和剪切区等过程产生。从对流层向上传播，这些波在中层大气中破裂，包括平流层和中间层，在那里它们沉积动量并对流动施加强大的阻力（见 Teixeira 2014）。

2.5

陆面和海洋模型以及数值天气预报中的耦合数值模型

2.5.1

陆面模型

土壤、植被、雪、山脉和水体的表征是数值天气预报系统的重要组成部分。土地可以影响天气、天气效应的强度以及人类活动的演变。地表状态异常的影响可能持续数天，从而增加了正确初始条件和模型演变的重要性。对地表过程的精细表征及其准确初始化有潜力进一步改善天气预报，达到月度范围，正如可预报性研究所示。

陆面模型（LSMs）用于表示和参数化陆面过程。LSMs 之所以重要，是因为它们为数值天气预报（NWP）和气候模型提供了必要的下边界条件（Ek 等，2003）。它们还计算辐射通量、热通量和水分通量，以供 NWP 或气候模型使用。这些通量通常是中尺度环流的主要驱动机制。此外，陆面过程还会影响近地面预报，如 2 米温度、10 米风速、边界层结构和降水预报。

地表的能量和水分预算控制着基质和植被的温度和湿度，这些与大气相互作用。能量守恒方程可以写成单位质量或单位面积的表面正在经历能量增益或损失的形式：

Z. Pu 和 E. Kalnay

(25)

其中 R 是表面的净辐射，即向下的长波辐射、向下的短波辐射、向上短波辐射（由表面反射，由表面的反照率控制）和向上长波辐射（表面发射）的总和。G 是地面热通量。它也可以解释为表面下方热储存的负速率。E 是蒸发速率，LE 代表潜热通量，H 是显热通量。因此，LSM 使用表面和冠层层中的参数计算显热和潜热通量。

土壤温度传递方程是土壤水分模型（LSM）的一部分。例如，可以使用傅里叶扩散定律来控制土壤的热量和水分传递：

ρC

(26)

其中 (ρC) 是体积土壤热容 (J mK)。它是土壤质地和土壤湿度的函数。T 是土壤温度，z 是垂直坐标（距离表面的距离），t 是时间（秒），λ 是热导率。

土壤热容可以估算为其各相热容的加权和。然后

ρC

Þ¼

" θ

Þ ρCð

Þþ θ

ρC

下标 m 和 w 分别指土壤基质和水。土壤水分运动遵循“理查兹”方程：

@θ

@t ¼

@θ

@γ

@z " S

(27)

λ和γ是土壤质地和土壤湿度的函数。S 代表源（降雨）和汇（蒸发）。

表面水预算方程表示水的守恒为：

¼ P-R-

(28)

其中 DS 是土壤水分变化，P 是降水量，R 是径流，E 是蒸发。具体来说，蒸发是土壤湿度和植被类型、根系深度/密度、绿色植被覆盖等的函数。具体而言：

¼ Eþ

þ E

(29)

E 是来自土壤/植被的总表面蒸发；E 是来自土壤的直接蒸发；E 是通过植物冠层的蒸腾；E 是来自冠层截留降雨的蒸发；E 是来自雪层的升华。

在数值天气预报模型中，全面考虑表面能量、水分平衡、加热和水分输送的土地表面模型（LSMs）的纳入已被

数值天气预报基础：模型、数值方法和数据……

已被证明对近地面大气条件的数值预测以及定量降水预报（例如，陈和杜迪亚 2001）是有益的。

许多不同中心开发了多种土地表面模型（LSM），用于不同的目的。例如，目前的中尺度社区天气研究与预报（WRF）模型（Skamarock et al. 2008）有六种土地表面模型选项。由多个美国政府机构（包括 NCEP 和空军）和大学（如俄勒冈州立大学）开发的被称为 NOAH 土地表面模型（Ek et al. 2003），也用于 NCEP NAM 区域模型。

2.5.2 海洋模型

地球的海洋成分在不同时间尺度上对大气有重要影响。一个完全耦合的海洋系统模型可能包括表面波、海洋和海冰。

众所周知，波浪通过改变表面粗糙度影响大气的海洋边界层。大多数气候模型与海洋模型耦合。一些数值天气预报模型（根据其应用）也与海洋模型耦合。集合和季节预报系统使用耦合的气候-海洋模型，其中包括对海洋一般环流及其相关耦合反馈过程的模拟。

海洋波浪建模：海洋波浪建模用于预测海洋波浪的产生和演变及其相关能量。许多运营中心已经开发了独立或耦合模式的海洋波浪模型。例如，欧洲中期天气预报中心（ECMWF）开发了波浪模型（WAM），该模型与大气模型耦合或以独立模型的形式在有限区域波浪（LAW）配置中运行。自 1998 年以来，ECMWF 一直在运行一个耦合预报系统，其中综合预报系统（IFS）的大气部分通过交换查尔诺克参数与 WAM 进行通信，该参数决定了海面粗糙度（Janssen 2004）。在美国国家海洋和大气管理局（NOAA）国家环境预报中心（NCEP），开发了波浪模型 WAVEWATCH III（例如，Tolman 2014）。WAVEWATCH III 解决了波数-方向谱的随机相位谱作用密度平衡方程。该方程的隐含假设是介质（水深和水流）以及波场本身的特性在时间和空间尺度上变化，远大于单个波的变化尺度（详见 Tolman 2014）。

海洋建模：除了海洋波浪模型外，海洋模型还被开发用于计算海面高度、洋流、盐度和温度的时间演变。与大气模型类似，海洋模型设计为静水或非静水。一个斜压的原始方程海洋模型包含质量（连续性）和动量的守恒方程。热力学也用于描述盐度和温度。许多海洋模型的数值方案和参数化方法与大气模型中使用的相似，尽管近年来有许多新发展。海洋模型需要初始化。许多现有海洋模型使用三维变分数据同化方法。海洋模型已被广泛应用。

Z. Pu 和 E. Kalnay

在操作性集合预报系统和季节性预报系统中。例如，欧洲中期天气预报中心（ECMWF）目前使用分辨率为 1 的海洋模型，该模型由三维变分同化系统初始化。自 2013 年以来，集合预报从预报开始就与大气-波浪-海洋模型耦合。这一点很重要，因为它能够捕捉大气与海面温度之间的双向反馈；例如，当热带气旋缓慢移动时，它可以使海面降温。国家环境预报中心（NCEP）正在实施一种新的局部集合变换卡尔曼滤波器（LETKF）和一种混合增益，结合了三维变分 GODAS 和 LETKF，参考 Penny（2014）。

海冰建模：海冰是地球系统的重要组成部分；它具有很高的反射率，改变了吸收的太阳辐射量；它改变了形成和融化海冰的海洋的盐度；并且它作为大气与海洋之间热量和动量通量交换的障碍。然而，目前的气象预报系统通常不动态预测海冰。在耦合预报系统中，海冰建模与海洋建模耦合，以表示海冰的动态和热力学演变，主要用于季节性和年际预测。

2.5.3 数值天气预报中的耦合数值模型

随着陆地和海洋模型的快速发展，许多数值天气预报（NWP）模型已成为耦合数值模型。例如，已证明陆面模型与大气模型的耦合可以显著改善定量降水预报（QPF）。此外，为了准确预测热带气旋的形成、增强和消散，耦合的海洋和大气模型被认为是必要的。在强调飓风强度预报的区域模型中（如飓风 WRF 或 HWRF），海洋喷雾也已被纳入参数化考虑。耦合的大气、陆地和海洋模型对于中期天气预报和延长范围预报（超过 10 天）是必要的。

陆地-大气耦合可以通过陆面模型实现，输入近地表和大气条件以及地球表面的特征到陆面模型，并通过提供水分和能量通量输出到大气模型。

在耦合的海洋-大气建模系统中，最多使用三种模型：大气模型、海洋模型和波浪模型。然而，由于计算成本、关注的时间尺度和预期的应用，许多系统今天只耦合这三种模型中的两种，通常是海洋-大气耦合建模系统或波浪-大气耦合建模系统。耦合发生在空气-海洋界面。例如，当使用所有三种模型时，大气模型向波浪模型提供表面应力，波浪模型利用这些信息推导二维波能谱。波浪模型向大气模型提供波浪引起的粗糙度长度，以便用于计算表面通量，这也需要海洋模型提供的海表温度（SST）。波浪模型产生的波浪引起的应力与大气模型的表面通量和辐射一起被海洋模型用于推导海表温度（SST）。

数值天气预报基础：模型、数值方法和数据……

虽然综合观测数据通常足以启动大气模型，但海洋表面以下的观测数据稀少。这导致海洋模型在任何模拟或预测的开始时间之前运行数月或数年，以便开发出一个代表性的三维海洋状态。在同化期间，海洋模型受到全球分析或大气模型提供的表面风应力、从卫星数据得出的观测海面高度异常以及观测到的海表温度（SST）的强迫。波浪模型通常不需要在耦合模型模拟开始之前进行初始化，除非第一天的波浪特征很重要。

数据同化

如本章开头所提到的，要进行预测，我们需要了解大气和地球表面（陆地和海洋）的当前状态。现代数值天气预报广泛利用地面和卫星观测，以及常规观测（例如，地面观测、气象站的探空仪、船只、浮标和其他组件）。这些观测提供了大气、海洋和陆地表面信息。虽然卫星现在提供了大部分数据，但更多的观测仍然很重要。

气象预报中心生成的天气预报使用数据同化来根据观测数据估计预报模型的初始条件。预报的质量取决于我们如何有效利用来自全球观测系统的实时信息，该系统由众多卫星仪器、气象站、船只、浮标和其他组件组成。数据同化的目的是利用观测数据和短期预报确定最佳的大气状态。数据同化通常是一个顺序时间步进过程，在这个过程中，先前的模型预报与新接收到的观测数据进行比较，然后更新模型状态以反映观测结果，接着启动新的预报，依此类推。这个过程中的更新步骤通常称为分析；用于生成分析的短期模型预报称为背景。

3.1

最小二乘法理论

根据 Talagrand（1997）的说法，最佳的气氛状态估计（分析）是通过将关于气氛的先前信息（背景或初始猜测）与观测结果结合而获得的，但为了最佳地结合这些信息，我们还需要关于这些“信息片段”中错误的统计信息。给定两个独立的观测值（或信息片段）T 和 T，确定标量（例如，真实温度 T）真实值的最佳估计的经典例子，作为统计估计的引入：

þ e

;

þ e

(30)

Z. Pu 和 E. Kalnay

观察存在我们未知的误差。让 E( ) 代表期望值，即如果进行许多类似的测量所获得的平均值。我们假设测量 T 和 T 的仪器是无偏的：E (T

T) = E (T

T) = 0 或者，等价地，

E e

ð Þ ¼

E e

ð Þ ¼

(31)

我们知道观测误差的方差：

E e

¼ σ

1 和 E e

¼ σ

(32)

我们还假设这两个测量的误差是无关的：

E ee

Þ ¼

(33)

方程 (31)、(32) 和 (33) 表示我们所需的关于实际观察的统计信息。我们尝试通过两个观察值的线性组合来估计 T，因为它们代表了我们对 T 的真实值所拥有的所有信息：

¼ a

þ a

(34)

“分析”应该是公正的：

E T

ð Þ ¼

E T

ð Þ

(35)

这意味着

þ a¼

(36)

T 将是 T 的最佳估计，如果选择的系数能够最小化 T 的均方误差：

E T

E aT

þ a

(37)

受约束条件 (37) 的限制。代入 a = 1

" a

最小化σ与

尊重给予：

a¼

σþ

a¼

σþ

(38)

数值天气预报基础：模型、数值方法和数据……

þ σ

2 2

þ σ

2 2

(39)

即，观测值的权重与测量的“精度”或准确性成正比（定义为观测误差方差的倒数）。此外，将系数（39）代入（37），我们得到分析方差与观测方差之间的关系：

σ¼

σþ

(40)

即，如果系数是最优的，并且误差的统计是准确的，那么分析的“精度”（定义为方差的倒数）是测量精度的总和。更重要的是，方程（40）还表明“分析”的误差方差σ

% &

小于任一的误差方差

“背景”或“观察。”

根据最小二乘法理论，可以从两条不确定的信息（背景和观察）中进行分析，并使分析的准确性超过单独使用其中任何一条信息的准确性。

3.2

同化方法

在实际操作中，分析 xis 是通过将创新与背景（模型预测或初始猜测）相加，并使用基于预测和观测的估计统计误差协方差确定的权重 W

¼ x

W y

H x

% &

(41)

创新在此被定义为观测值与模型“猜测”观测值之间的差异，即 y "H(x)。具体而言，背景（模型预报）被插值到观测位置，如果它们是不同类型，则从模型变量转换为观测变量 y（例如卫星辐射或雷达反射率）。因此，观测的第一次猜测是 H(x)，其中 H 是执行必要插值和从模型变量到观测空间转换的观测算子。

不同的分析方案基于（41），但在结合背景和观测以产生分析的方式上有所不同。早期的方法，如连续修正法，形式类似于（41），权重是通过经验确定的。权重是观测与网格点之间距离的函数，分析会迭代多次。在最优插值中，权重矩阵 W 是通过最小化每个网格点的分析误差来确定的。即使是现代先进的数据同化方法也如此。

Z. Pu 和 E. Kalnay

许多数据同化系统遵循这一理论。

方法可以以类似的方式进行解释，如下所述。详细信息请参见 Kalnay（2003）。

3.2.1

三维变分（3D-VAR）数据同化方法

在 3D-VAR 方法中，定义一个与分析、背景和观测之间距离的平方成正比的代价函数。直接最小化代价函数以获得分析。Lorenc（1986）表明，如果将代价函数定义为，则最优插值（OI）和 3D-VAR 方法是等价的：

¼ 1=

2 y

H x

ð Þ

H x

ð Þ

þ 1=

2 x

" x

(42)

成本函数 J 在 (42) 中衡量了场 x 与观测值之间的距离（成本函数中的第一项）以及与第一次猜测或背景 x 之间的距离（成本函数中的第二项）。这些距离分别由观测误差协方差 R 和背景误差协方差 B 进行缩放。成本函数的最小值在 x = x 时获得，这被定义为“分析”。在 (41) 和 (42) 中获得的分析是相同的，如果 (41) 中的权重矩阵由以下给出：

BHHBH

þ R

(43)

在 3D-VAR 中，直接对（42）进行最小化，这允许额外的灵活性和数据的同时全球使用。

3.2.2

四维变分（4D-VAR）数据同化方法

4D-VAR 是 3D-VAR 的重要扩展，它允许在时间区间(t, t)内分布的观测值（例如，Courtier 和 Talagrand 1987；Derber 1989；Bouttier 和 Rabier 1997）。成本函数包括一个测量在区间开始时与背景的距离的项，以及对每个观测增量的成本函数在观测时间的模型积分上进行的时间求和：

J x t

ð Þ

x t

Þ " x

ð Þ

Bx t

Þ " x

ð Þ

i¼0

H x

Þ " x

H x

Þ "

(44)

控制变量（成本函数最小化的变量）是模型的初始状态，时间间隔为 x(t)，而在区间结束时的分析由模型积分给出，解为 x(t) = M[x(t)]。因此，模型被用作强约束，即分析。

数值天气预报基础：模型、数值方法和数据……

解法必须满足模型方程。换句话说，4D-VAR 寻求一个初始条件，使得预测在同化区间内最符合观测值。4D-VAR 方法假设模型是完美的，这一缺点在于，例如，它对区间开始时的旧观测和区间结束时的新观测给予相同的可信度。Derber（1989）提出了一种修正恒定模型误差的方法（在同化区间内的恒定形状）。

为了最小化成本函数，成本函数相对于背景和观测组件的梯度可以表示为：

@x t

ð Þ

¼ B

0 x t ð Þ " x

ð Þ

(45)

@ J

@x t

ð Þ

i¼0

L t, t

殿下 x

Þ "

(46)

方程（46）表明，4D-VAR 最小化需要计算梯度，即在前向积分过程中计算观测时间 t 的增量[H (x)"y ]，将其乘以 HR，并使用伴随模型将这些加权增量积分回初始时间。由于向后伴随积分的部分在多个时间间隔中是共同的，因此（46）中的求和可以更方便地安排。例如，假设同化的时间间隔是从 00 到 12 Z，并且每 3 小时有一次观测。我们在前向积分过程中计算加权的负观测增量 d。

殿下 x

Þ "

HRd。伴随模型 L(t, t) = Lapplied 在一个向量上“推进”它从 t 到 t。然后我们可以将 (46) 写为：

@ J

@x¼ dþ L

þ L

(47)

从（45）加（46）或（47），我们得到成本函数的梯度，最小化算法适当地修改控制变量 x(t)。在此变化之后，计算新的前向积分和新的观测增量，并重复该过程。

4D-VAR 最重要的优点是，如果我们假设(a)模型是完美的，(b)初始时刻 B 的先验误差协方差是正确的，那么可以证明，最终时刻的 4D-VAR 分析与扩展卡尔曼滤波器的分析是相同的（Lorenc 1986；Daley 1991）。这意味着，隐含地，4D-VAR 能够将预测误差协方差从 B 演变到最终时刻。不幸的是，这种隐含的协方差在周期结束时不可用，新的分析误差协方差也不可用。换句话说，4D-VAR 能够找到最佳线性无偏估计，但无法获得其误差协方差，除非以近似形式。

Z. Pu 和 E. Kalnay

与此同时，在 3D-VAR 中，背景项是统计定义的，并且随时间静态。因此，3D-VAR 和 4D-VAR 都无法表示数据同化中的随时间变化的背景（模型预测）误差协方差。

3.2.3 集合卡尔曼滤波器

在随机集合卡尔曼滤波器中，同时进行 K 个数据同化周期的集合（Evensen 1994）。所有周期同化相同的真实观测值，但必须向每个集合数据同化的成员中同化的观测值添加不同的随机扰动。另一种类型的集合卡尔曼滤波器被称为确定性或平方根滤波器（Tippett 等，2003；Ott 等，2004）。在这里，只有一个数据同化，新的分析扰动是通过平方根算法从背景（预测）集合扰动中获得的。

在时间 t 完成分析组合后，K 个预测 x(t) = M[x(t)]，可以从 K 个预测 x(t) 中获得预测误差协方差的估计，如下所示：

" 1X

k¼1

" x

(48)

其中上划线表示集合平均。这往往会低估预测误差的方差，因为存在非线性。因此，通常会在集合卡尔曼滤波器中实施协方差膨胀。

3.2.4 混合 3D-VAR/4D-VAR 和集合卡尔曼滤波器

Hamill 和 Snyder（2000）还建议了一种 3D-VAR 和集合卡尔曼滤波的混合方法，其中预测误差协方差是通过（常数）3D-VAR 协方差 B 的线性组合获得的

Pð

" α

þ αB-

(49)

其中 α 是一个可调参数，范围从 0（纯 3D-VAR 和纯集合卡尔曼滤波）到 1。由于集合卡尔曼滤波的协方差仅从有限的集合成员样本中估计，因此其秩为 K。

1，远小于模型的自由度，因此它是秩缺陷的。与 3D-VAR 的结合，从许多估计的预报误差中计算得出（例如，使用 Parrish 和 Derber 1992 的方法），可能会改善这个抽样问题并“填充”误差协方差。在 Hamill 和 Snyder（2000）的实验中，最佳结果是在α值较低时获得的，介于 0.1 和 0.4 之间，表明使用集成演变的预报误差协方差的良好影响。他们发现 25-50 个集成成员足以提供集成卡尔曼滤波的好处（但在使用比这里使用的准地转模型更复杂的模型时，这可能会有所不同）。最近，混合 3D-VAR 和 4D-VAR 以及集成卡尔曼滤波方案（参见

数值天气预报基础：模型、数值方法和数据……

最近，3dEnVar 和 4dEnVar 分别在 NCEP 实施并运行（Wang et al. 2013; Kleist and Ide 2015），与之前的 3D-VAR 分析方案，即网格点统计插值（GSI）系统相比，取得了重大改进。另一种混合方法（Penny 2014）结合了卡尔曼增益而不是协方差，并已由 ECMWF 进行了非常好的测试（Hamrud et al. 2015）。

集合卡尔曼滤波方法有几个优点：(a) K 的数量级为 10–100，因此计算成本（与 OI 或 3D-VAR 相比）增加了 10–100 倍。尽管这个增加的成本看起来很大，但与扩展卡尔曼滤波相比，它是小得多的，后者需要的成本增加与模型的自由度数量级相当。(b) 集合卡尔曼滤波不需要开发线性和伴随模型。(c) 它不需要对预报误差协方差的演变进行线性化。(d) 它可以为集合预报提供出色的初始扰动。目前，集合卡尔曼滤波似乎是未来最有前途的方法之一（Houtekamer 和 Zhang 2016）。