详细设计¶
系统功能设计¶
- GRAPES OSE：评估雷达资料有无、时空密度对GRAPES业务系统预报的总体影响效果；
- GRAPES OSSE：构建雷达反射率、径向风速度的观测模拟算子；建立Nature Run参考大气实验数据集；通过对比试验评估新增雷达资料对数值预报的影响；
- GRAPES FSO：建立基于GRAPES三维变分同化的集合FSO评估模块。
系统采用方法¶
OSE
- GRAPES OSE 以业务运行的GRAPES-MESO（V4）或已定型的GRAPES区域高分辨率版本为基础；
- 观测资料处理采用GRAPES系统中的obsproc模块，该模块可以处理常规观测和雷达回波资料；
- 同化方法为GRAPES 3DVar;
- 控制实验为同化常规观测资料的GRAPES预报；
- 对比试验为同化常规观测资料、不同区域、不同时空密度及分辨率的雷达反射率;
- 检验方法采用RMSE、TS评分。
OSSE
- GRAPES OSSE 以业务运行的GRAPES-MESO（V4.3）或已定型的GRAPES区域高分辨率版本为基础；
- Nature Run参考大气试验数据集：以ERA5为背景场，同化区域内常规观测资料；
- 同化方法为GRAPES 3DVar;
- 开发常规观测资料、3D雷达综合反射率、雷达径向风速度模拟观测算子；
- 对比试验以NCEP-GFS/GRAPES-GFS为背景场，分别同化常规观测资料、不同区域、不同时空密度及分辨率的3D雷达综合反射率、雷达径向风速度；
- 检验方法采用RMSE。
集合FSO
FSO是初值敏感性分析工具。敏感性分析是指分析和确定数值预报初始条件中造成预报误差的分析误差的特征及其所在区域，将关注的重点集中到分析与理解初始场中与预报效果密切相关的分析误差（这些误差在预报过程中增长最快，故而对预报影响最大。从另一个角度理解，也可视之为对预报事件发生发展影响最大的扰动）的特征和区域上。敏感性分析的一种实用工具是伴随模式。伴随模式是和资料变分同化概念密切相联系的。我们知道，积分数值模式需要由观测资料经分析形成初始条件，这种分析方法现在统称为资料的同化分析。具体方案从早期的多项式拟合发展到最优插值时，遇到了不能直接同化数量与日倍增的非常规观测资料的困难，资料变分同化方法由此发展出来。现今实际应用的方案已发展有三维变分分析、集合卡尔曼滤波、四维变分分析等。
目前，FSO方法的实现上主要有两种方式：
- 基于伴随模式的敏感性分析方案；
- 基于集合预报的敏感性分析方案。
两种方案的数学原理是一致的，主要不同之处实践过程中如何描述扰动的线性演变与同化分析过程中增益矩阵的处理。
伴随FSO方法
伴随模式的优点就在于它提供了一种定量确定对某一天气特征有着重要影响的初始扰动或物理因子的有效方法。其实气象研究中经常讨论敏感性的问题，也就是估计如初始条件或物理、动力因子改变后将如何影响预报的天气特征和为什么会带来这些 影响。过去常用的方法是首先确定一个基准试验，然后改变某些条件后，比如OSE试验中同化不同的观测资料，数值试验中多用的平滑、倍增或减小某些物理场等，再次向前积分数值模式，比较与基准试验之间的差异来讨论。其实这些作法在某些情况下是不太合适的，并且对于所关心的天气特征也缺乏针对性。另外，由于有些相似的改变也许会产生差异相当大的结果，为了要得到带有普遍性的结论，就必须进行大量改变试验结果的集成，带来计算量的巨增。而伴随模式针对所定义的天气特征，通过反向积分伴随模式则能够直接高效地得到对影响该天气特征有重要动力学意义的初始扰动和物理、动力过程的特征及所在区域，并且还可以对关心的不同天气特征灵活设计不同的敏感性分析。
基于伴随模式的初值敏感性分析如图所示:
基于伴随模式的敏感性分析方法示意图
基于伴随模式的敏感性分析主要有两种方法。其一是梯度敏感，Errico和Vukiceive1992年使用MM4不含线性化湿过程的伴随模式对两个实际气旋个例作了研究；Vukiceive等1995年还应用该方法讨论了阿尔卑斯山背风坡气旋的可能诱发机制。Rabier等1993、1995年对欧洲中期数值预报中心48h短期数值预报误差的敏感性研究表明主要预报误差在多数个例中都可用初始条件的分析误差来解释，指出了ECMWF分析场的一些明显系统性误差,并给出了其一般型式，并指出分析中一些细微的，但是适合的改动往往可以极大地改善预报效果，此方法可用于初始误差对预报误差作用的诊断分析。利用伴随模式讨论敏感性问题的另一种方法是奇异向量分析。梯度敏感分析得到的预报误差相对于初始条件的梯度可以解释为分析误差中快速增长分量的总和，而奇异向量则代表一段时间内增长最快的扰动，这种扰动也称为“最优扰动”，它相当于分析误差投影到奇异向量上，由此依奇异值大小降序排列对应的就是可能增长速率由快到慢的奇异向量。前面的主导奇异向量虽然只代表部分分析误差，但是它的增长却左右了预报误差，所以奇异向量能准确定量地描述在特定时段内最不稳定扰动的特征。伴随敏感性分析得到的敏感区确定了对于特定预报事件后期演变最为关键的地区，从而可以作为目标观测设计的目标区。同时，分析还能够显示初始条件中那些显著影响预报的场及其特征，也就指导了观测应该如何布局以减少分析误差，伴随敏感性分析由此成为目标观测设计的一个基础（当然，除去这里所讲的基于伴随模式的梯度敏感与奇异向量分析，还有如前所述的其它方法）。另外，敏感性分 析也可以作为制作集合预报中产生初始小扰动的方法，这种方法构造出的初始扰动考虑了初始分析误差的可能分布状态，利于提高最终预报的准确度。目前，奇异向量是业务集合预报中除孵化法外生成初始扰动常用的一种方法。
本项目技术路线参考美国国家大气研究中心开发的基于WRF和WRFDA的FSO系统。该系统旨在用于计算观测资料对区域数值准确性的贡献。通常来说，同化观测资料可以提高模式初始场的精度，进而提高数值预报的准确性，但一般很难去区分不同观测资料在模式预报结果中的贡献。FSO就是基于此目的开发的。FSO利用WRFDA中的4DVAR伴随模式（ADM）来计算预报误差对观测增量的梯度（敏感度），来评估每种观测资料对预报误差的影响。FSO不需要增加或减少观测资料既能了解观测资料对同化的影响。
WRFDA-FSO系统基于WRFDA同化系统、WRF的切线性和伴随模式（WRFPLUS）（Auligne，2010）而建立，其系统流程结构如图所示：
WRFDA-FSO流程结构图
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 x_a}{\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\) 为正值时，预报误差增大，同化观测降低预报能力。具体实现上，在四维变分同化系统中，可以利用既有的伴随模式及隐式求解增益矩阵 来评估分析增量及每种观测的对于预报误差的贡献。
在实际的同化系统中，同化资料后主要考察在背景场的基础上同化资料对预报误差的贡献。由背景场作为初始场的预报误差和分析场作为初始场的预报误差分别为:
定义 \(e_a\) 和 \(e_b\) 的差为：
(10)¶\[\delta e = e_a - e_b\]所以
(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\](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)]\](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) 中，增益矩阵及线性模式与增益矩阵的组合可以表示为 (Kalnay 等，2012):
(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\) 为分析扰动在观测空间的表示。
- 集合产生：采用GRAPES区域集合（ MEPS）方法产生；集合水平分辨率为9km，成员数为15个；
- 集合预报输出预处理：对集合预报模式输出进行格式和变量转换；
- 观测资料处理模块：采用业务GRAPES-3DVar obsproc进行资料处理；
- 影响评估：开发观测影响评估模块。
系统流程设计¶
OSE
OSE系统流程设计
OSSE
OSSE系统流程设计
EFSO
EFSO系统流程设计
系统接口设计¶
- 以业务运行的GRAPES-MESO（V4.3）业务应用为基础；
- 背景场输入格式为grib1/grib2;
- GRAPES模式输出为bin和ctl;
- 观测资料输入为cimiss库；质控处理后格式为GRAPES-3DVar格式；
- 评分、统计检验结果为txt；格点误差检验为bin或netcdf;
- 图形产品格式为png;
- 文字产品为word和pdf
系统开发环境与部署环境¶
- 系统开发环境：
- 曙光Linux系统，80核，Intel2016编译器。
- 系统部署环境：
- 曙光派HPC，Intel2017编译器。
系统运行维护设计¶
系统运行调度：
- 方法一：采用ECFlow调度（适应数值中心业务运行需求）；
- 方法二：采用cron与python结合调度（适应研究应用需求）。
系统运行维护：
- 日志维护：每一个进程、任务都产生对应的日志文件；
- 网页监控：通过ECFlow图形窗口进行运行状态监控。