概要设计

OSSE

OSSE是一种建模实验,用于在没有实际观测数据时评估新观测系统对运行预测的影响。OSSE是用一个性能良好的模式,模式自由运行一段时间模拟大气作为大气的真值,称为Nature run, 并且用此真值加上随机误差构造观测,然后将该构造的资料应用在另一个数值模式系统,做同化预报,结果和Nature run比较,以此来评估新增观测资料的影响。其流程结构如下图所示:

../_images/osse.png

OSSE试验流程概念图

OSSE试验一般由以下功能模块构成:

  • 大气三维综合仿真观测构建模块:建立Nature Run参考大气,构建大气三维综合仿真观测系统,模拟实现天地空各类观测。可以使用除GRAPES之外的其他模式,如WRF,GFS等;可以使用资料同化系统的Observer功能,利用“真实”大气模拟产生所需要的模拟(合成)观测资料

  • 同化预报模块:同化模拟(合成)观测资料

  • 开发常规观测资料、3D雷达综合反射率、雷达径向风速度模拟观测算子;

  • GRAPES模式预报模块:

  • 试验结果评估模块:利用GRAPES模式后处理模块以及MET模式预报验证模块计算预报结果的预报偏差、降水TS评分、预报与实况对比、地面要素检验等,再进行统计对比分析,对模拟观测的贡献及影响进行评估。

基于集合预报的FSO方法

FSO是评估预报误差对观测敏感性的一种分析方法。它不但可以评估观测系统中具体组成部分(如探空观测系统,地面观测系统)对改进预报技巧的有效性,而且能够定量给出各观测变量(如风场,温度)在当前预报中的相对影响,同时给出这种影响的空间分布。FSO应用的一个重要假设位观测作用可以线性分解与叠加的,即预报误差中的主要部分是可以用非线性模式的切线性模式或是一组集合预报非线性扰动的线性组合表示。

目前,FSO方法的实现上主要有两种方式:1)基于伴随模式的敏感性分析方案;2)基于集合预报的敏感性分析方案。两种方案的数学原理是一致的,主要不同之处实践过程中如何描述扰动的线性演变与同化分析过程中增益矩阵的处理。

FSO的数学原理简要描述如下(Baker,2000)。在一个循环资料同化预报系统中,当前分析时刻的分析场一般由前一个时刻分析场产生的短时预报场经过当前观测修正而得到。令 \(x_a\) 表示分析时刻 \(t=0\) 的分析场,而 \(x_b\) 表示前一个同化预报循环产生的短期预报背景场。则观测对分析的影响可以表示为:

(1)\[x_a = x_b + \mathrm{K} \delta y \; ;\; \delta y = y - H(x_b)\]

(1)\(x_a\) 为分析场,\(x_b\) 为背景场,\(y\) 为观测值,\(\mathrm{K}\)Kalman 增益矩阵,\(H\) 为观测算子,\(\delta y\) 为更新变量。从 (1) 可以推导出分析场对观测和背景场的敏感梯度可写为:

(2)\[\frac{\partial x_a}{\partial y} = \mathrm{K}^\intercal \; ; \; \frac{\partial x_a}{\partial x_b} = \mathbf{I} - \mathrm{H}^\intercal \mathrm{K}^\intercal \; ; \; \mathrm{H} = \frac{\partial H}{\partial x_b}\]

(2)\(\mathrm{H}\) 为切线性观测算子,上标 \(\intercal\) 代表矩阵转置,\(\mathrm{K}^\intercal\) 代表了同化系统的伴随系统。

考虑一个预报,其预报方程可以写为:

(3)\[x^f = M(x^0)\]

(3)\(x^f\) 为模式预报场,\(x^0\) 为模式初始场,\(M\) 为非线性数值预报模式。

定义在某时刻的一个模式预报误差函数:

(4)\[J = \frac{1}{2} ( x^f - x^t )^\intercal \mathrm{C} ( x^f - x^t )\]

(4) 中,\(\mathrm{C}\) 为权重矩阵,一般为对角矩阵。从 (4) 中可以推导出预报误差函数 \(J\) 对初始场(如果经过了同化分析,初始场就是同化系统中的分析场,\(x^0\) 即为 \(x_a\) )的敏感性梯度:

(5)\[\frac{\partial J}{\partial x_a} = \frac{\partial x^f}{\partial x_a} \mathrm{C} ( x^f - x^t ) = \frac{\partial M(x_a)}{\partial x_a} \mathrm{C} ( x^f - x^t ) = \mathrm{M}^\intercal \mathrm{C} ( x^f - x^t )\]

(5)\(x^t\) 为真值,\(\mathrm{M}\) 为非线性模式的切线性模式,\(\mathrm{M}^\intercal\)\(\mathrm{M}\) 模式的伴随模式。联合式 (1)(4), (5) 可推导出模式预报误差函数 \(J\) 对观测的敏感梯度:

(6)\[\frac{\partial J}{\partial y} = \frac{\partial x_a}{\partial y} \frac{\partial J}{\partial x_a} = \mathrm{K}^\intercal \frac{\partial J}{\partial x_a} = \mathrm{K}^\intercal \mathrm{M}^\intercal \mathrm{C} ( x^f - x^t )\]

因此,从式中可以看到,预报误差对观测的敏感性梯度需要计算模式的伴随和同化系统的伴随。在同化系统中,观测资料对预报误差的贡献(Observation impacts)一级近似为:

(7)\[\delta J = \langle \frac{\partial J}{\partial y}\, , \delta y \rangle = \langle \mathrm{K}^\intercal \frac{\partial J}{\partial x_a }\, , y - H(x_b) \rangle= \langle \frac{\partial J}{\partial x_a} \, , \mathrm{K} ( y - H(x_b)) \rangle = \langle \frac{\partial J}{\partial x_a} \, , \delta x_a \rangle\]

(7) 中,\(\langle \, , \rangle\) 为内积算符。

由式 (7) 可知,观测对于预报的影响可以解释为:\(\delta y = y- H(x_b)\) 在预报误差对于观测敏感性梯度上的投影;或是,分析增量在预报对于初值的敏感性的投影。从式 (7) 可以看出,当误差改变为负值时,预报误差减少,同化观测提高预报能力;当 \(\delta e\) 为正值时,预报误差增大,同化观测降低预报能力。具体实现上,在四维变分同化系统中,可以利用既有的伴随模式及隐式求解增益矩阵 来评估分析增量及每种观测的对于预报误差的贡献。

在实际的同化系统中,同化资料后主要考察在背景场的基础上同化资料对预报误差的贡献。由背景场作为初始场的预报误差和分析场作为初始场的预报误差分别为:

(8)\[e_b = \frac{1}{2} (x_b^f - x^t)^\intercal \mathrm{C} (x_b^f - x^t) = \frac{1}{2} \langle (x_b^f - x^t) , (x_b^f - x^t) \rangle\]
(9)\[e_a = \frac{1}{2} (x_a^f - x^t)^\intercal \mathrm{C} (x_a^f - x^t) = \frac{1}{2} \langle (x_a^f - x^t) , (x_a^f - x^t) \rangle\]

定义 \(e_a\)\(e_b\) 的差为:

(10)\[\delta e = e_a - e_b\]

从式 (8)(9) 可以得到一阶偏导数:

(11)\[\frac{\partial e_a}{\partial x_a^f} = \mathrm{C} (x_a^f - x^t)\]
(12)\[\frac{\partial e_a}{\partial x_a^f} = \mathrm{C} (x_a^f - x^t)\]

所以

(13)\[\delta e = \langle \frac{\partial e_a}{\partial x_a^f} , (x_a^f - x^t) \rangle - \langle \frac{\partial e_b}{\partial x_b^f} , (x_b^f - x^t) \rangle\]

利用 (11)(12)(13) 可得:

(14)\[\delta e = \langle (x_a^f - x_b^f) , \frac{\partial e_a}{\partial x_a^f} + \frac{\partial e_b}{\partial x_b^f} \rangle = \langle (x_a - x_b) , \frac{\partial e_a}{\partial x_a} + \frac{\partial e_b}{\partial x_b} \rangle\]

(14) 利用式 (5) \(x_a^f-x_b^f\) 可以近似展开 (15) 可以推导出右边的表达式

(15)\[x_a^f - x_b^f = \mathrm{M}_b^\intercal (x_a -x_b) = \mathrm{M}_a^\intercal (x_a - x_b)\]

因此采用不同的近似,(14) 有不同的表达形式。Gelaro等(2007)给出了5种不同的表达形式:

(16)\[\delta e_1 = 2(x_a - x_b^\intercal \mathrm{M}_b^\intercal \mathrm{C} (x_a^f - x^t)\]
(17)\[\delta e_2 = (x_a - x_b)^\intercal [\mathrm{M}_b^\intercal \mathrm{C} (x_a^f - x^t) + \mathrm{M}_a^\intercal \mathrm{C} (x_b^f - x^t)]\]
(18)\[\delta e_3 = (x_a - x_b)^\intercal [\mathrm{M}_b^\intercal \mathrm{C} (x_b^f - x^t) + \mathrm{M}_a^\intercal \mathrm{C} (x_a^f - x^t)]\]
(19)\[\delta e_4 = (x_a - x_b)^\intercal [\mathrm{M}_a^\intercal \mathrm{C} (x_a^f - x^t) + \mathrm{M}_a^\intercal \mathrm{C} (x_b^f - x^t)]\]
(20)\[\delta e_5 = (x_a - x_b)^\intercal [\mathrm{M}_b^\intercal \mathrm{C} (x_b^f - x^t) + \mathrm{M}_b^\intercal \mathrm{C} (x_a^f - x^t)]\]

以式 (16) 为例,将 (1) 代入 (16) :

(21)\[\delta e_1 = 2 \delta y^\intercal \mathrm{K}^\intercal \mathrm{M}_b^\intercal \mathrm{C} (x_a^f - x^t)\]

基于集合的FSO,关键在于利用一组集合扰动的线性组合来近似表示扰动的时间演变,进而省略非线性模式的切线性与伴随模式。在集合预报系统中,已知分析扰动与预报扰动的前提下,式 (6) - (7) 中,增益矩阵及线性模式与增益矩阵的组合可以表示为 ( [KOML12] ):

(22)\[\mathrm{K} = \frac{1}{K-1} \mathrm{X}_0^a {\mathrm{X}_0^a}^\intercal \mathrm{H}^\intercal \mathrm{R}^{-1}\]
(23)\[\mathrm{M} \mathrm{K} = \frac{1}{K-1} \mathrm{M} \mathrm{X}_0^a {\mathrm{X}_0^a}^\intercal \mathrm{H}^\intercal \mathrm{R}^{-1} \approx \frac{1}{K-1} \mathrm{X}_{t|0}^f {\mathrm{Y}_0^a}^\intercal \mathrm{R}^{-1}\]
(24)\[\bigtriangleup e^2 \approx \frac{1}{K-1} \delta y_0^\intercal \mathrm{R}^{-1} \mathrm{Y}_0^a \mathrm{X}_{t|0}^{f^\intercal} \mathrm{C} (e_{t|0} + e_{t|-6})\]

其中:\(\mathrm{Y}_0^a = \mathrm{H} \mathrm{X}_0^a\) 为分析扰动在观测空间的表示。

本项目建立基于集合的FSO( [KOML12] )。具体实现表示如下:我们将应用一个时间序列的集合预报的扰动,来近似表示非线性模式的线性模式进而获得预报误差对初值的敏感梯度;以及观测集合的扰动观测来简化线性模式 \(\mathrm{M}\) 与增益矩阵 \(\mathrm{K}\) 的求解,进而实现观测对预报贡献的定量评估。

二期项目的主要任务为Refactoring EFSO:

  • 以OOPS为基础,重构EFSO

  • 利用eckit,fckit,Atlas,Eigen3等工具,软件工程现代化

  • 面向对象(可扩展为适应于各种集合预报系统)

  • 并行化处理

  • 提供考虑localization的接口

  • 业务化功能

MET系统架构和数据流

模式检验评估不仅是发展预报系统的重要组成部分,而且可以用来评价模式预报的准确性,从而为用户提供客观的预报依据。事实上,由于模式预报系统的底层物理过程十分复杂,其预报结果存在很大的不确定性,因此,通过评估检验来获得模式的特定属性和外在表现就成为用户研究和应用模式的重要手段。模式检验中最古老的方法是通过目视进行误差主观对比分析,尽管主观的目视误差分析或天气个例分析能够准确给出模式对降水系统诸如锋面、雨带预报性能的详细描述,然而目视分析存在显著缺陷,它一方面无法客观分析海量数据,另一方面通过自视检验获得的结论也具有主观性和非定量性,因此很难为用户提供有效的判别标准。随着模式检验需求的不断扩大和气象科研工作者的持续努力,一系列客观的模式评估检验方法得到了快速发展。

MET模式评估工具库(Model Evaluation Tools)是由美国国家大气研究中心(NCAR)研究应用实验室(RAL)开发的数值预报检验系统。MET除了支持传统的客观站点检验外,还可针对格点化实况对高分辨率数值预报产品检验。

其主要目的是为数值预报使用者、开发者搭建分析、评估数值模式预报产品的桥梁,它不仅为模式开发者提供实时的模式产品测试环境,为新的模式预报产品投入预报业务提供一揽子评估方案,而且为模式使用者客观把握模式的预报能力提供了有效手段。MET最初版本发行于2008年,主要针对于WRF数值模式产品的检验和评估,其最新版本MET8.1发行于2019年9月,提供各种通用数据格式的模式产品检验评估接口。本系统对MET的功能、数据流、算法作了进一步的模块化包装,尤其时候针对中国的数据环境,开发了一系列数据接口,方便用户的安装,使用,维护和升级。

MET的初始研发理念是为短期天气预报提供最为先进的数值天气预报模式产品评估检验方案,这种理念使得MET系统不仅囊括了传统的经典模式检验方法,同时不断地吸纳包括空间分析、诊断分析等在内的一系列新的检验方法。此外,MET增加了模式评估结果置信度检验等功能,这也使得MET快速取代传统的模式评估系统,比如MET已经取代了NCEP的模式实时校验系统。

MET是一种开放系统,可以在DTC网站自由下载。MET4.1以前的版本不包括台风预报检验,MET4.1及以后的版本增加了模式预报台风路径检验模块,但MET开发组认为台风路径检验本质上和模式其他要素检验是不同的,DTC给出的无论是台风检验程序包,还是说明文档,均为单独的一个部分:MET-TC。 MET采用模块化设计,具有很强的适应性,一种方法的检验模块可以独立运行而不依赖于其他组件,这种设计使得MET不仅能够在具有复杂数据输入、输出接口的各种大型数据库系统上运行,而且也方便运行于个人计算机。按照MET系统组件的功能的差异,可以将其分为格式转换模块、核心算法模块、统计输出模块以及绘图输出模块四个部分。

../_images/fig11.png

MET系统的基本架构及数据流(卷角方框表示数据的输入或输出,带阴影的椭圆表示可执行程序,箭头表示数据流)

上图给出了MET系统的基本架构及数据流。图中带阴影的椭圆表示MET的实际可执行程序,卷角方框表示数据输入或输出。可以看出,MET主要接收五类观测数据:格点数据、MODIS卫星观测数据、WWMCA格式云分析数据、ASCII格式站点观测数据、PrepBufr格式观测数据以及MADIS格式观测数据,其中PrepBufr是NCEP的一种站点格式的观测数据,而MADIS是NOAA所收集和整理的另一类站点格式数据文件。这些数据格式均被转换成NetCDF格式进行计算。在新版本中也可以直接输入Grib或Grib2格式数据。数据绘图程序可以显示转码后的格点或站点数据,生成PS格式图表。

MET的核心模块包括格点检验、站点检验、集合预报检验、小波检验、面向对象(基于目标)的MODE检验和时序检验六个部分,其中概率预报检验分别包含于站点检验和格点检验之中。MET的输出主要包括三类:NetCDF、ASCII格式计算结果以及PS格式或PNG格式图片文件。为了清楚反映模式的预报性能,通常还需要对输出结果进行再分析。表1列出了MET5.0部分程序及其功能介绍。 表1 MET主要的模块程序及功能

项目组专家建议

  • 注意使用的软件与数值中心软件环境的兼容

  • OSSE和EFSO建成之后,需要进行一系列的个例试验

  • 借鉴 http://jointisse.org 的数据

  • EFSO的localization的接口

  • EFSO的后处理图形生成功能