方法和技术路线¶
- 采用动态脚本语言建立数据收集与处理、模式编译与运行系统，实现“一键”运行，简化试验准备环节，建立循环同化试验的流程管理平台，提供网页监控界面，使得可以远程监控和管理运行试验。
- 动态脚本语言采用常用的Python语言，所有程序需要采用模块化设计，并使用GIT版本管理软件管理，流程管理采用Python Airflow实现，输出数据格式采用标准的NetCDF格式，提高易用性，绘图程序采用NCL语言编写。
FSO的数学原理简要描述如下（Baker，2000）。在一个循环资料同化预报系统中，当前分析时刻的分析场一般由前一个时刻分析场产生的短时预报场经过当前观测修正而得到。令 \(x_a\) 表示分析时刻 \(t=0\) 的分析场，而 \(x_b\) 表示前一个同化预报循环产生的短期预报背景场。则观测对分析的影响可以表示为:
\[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) 可以推导出分析场对观测和背景场的敏感梯度可写为：
\[\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\) 代表了同化系统的伴随系统。
考虑一个预报，其预报方程可以写为:
\[x^f = M(x^0)\]
式 (3) 中 \(x^f\) 为模式预报场，\(x^0\) 为模式初始场，\(M\) 为非线性数值预报模式。
定义在某时刻的一个模式预报误差函数:
\[J = \frac{1}{2} ( x^f - x^t )^\intercal \mathrm{C} ( x^f - x^t )\]
式 (4) 中，\(\mathrm{C}\) 为权重矩阵，一般为对角矩阵。从 (4) 中可以推导出预报误差函数 \(J\) 对初始场（如果经过了同化分析，初始场就是同化系统中的分析场，\(x^0\) 即为 \(x_a\) ）的敏感性梯度：
\[\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\) 对观测的敏感梯度：
\[\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）一级近似为：
\[\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_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\]\[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\) 的差为：
\[\delta e = e_a - e_b\]
\[\frac{\partial e_a}{\partial x_a^f} = \mathrm{C} (x_a^f - x^t)\]\[\frac{\partial e_a}{\partial x_a^f} = \mathrm{C} (x_a^f - x^t)\]
所以
\[\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\]
\[\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) 可以推导出右边的表达式
\[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种不同的表达形式：
\[\delta e_1 = 2(x_a - x_b^\intercal \mathrm{M}_b^\intercal \mathrm{C} (x_a^f - x^t)\]\[\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)]\]\[\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)]\]\[\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)]\]\[\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)]\]
\[\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):
\[\mathrm{K} = \frac{1}{K-1} \mathrm{X}_0^a {\mathrm{X}_0^a}^\intercal \mathrm{H}^\intercal \mathrm{R}^{-1}\]\[\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}\]\[\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\) 为分析扰动在观测空间的表示。