刘柏年皇群博张卫民曹小群赵 军赵延来

（1 国防科学技术大学海洋科学与工程研究院，长沙 410073；2 国防科学技术大学计算机学院，长沙 410073）

集合四维变分资料同化研究进展

刘柏年1，2皇群博1，2张卫民1曹小群1赵 军1赵延来1

（1 国防科学技术大学海洋科学与工程研究院，长沙 410073；2 国防科学技术大学计算机学院，长沙 410073）

背景误差协方差矩阵的精确定义是构建高水平资料同化系统的先决条件。传统四维变分资料同化（4D-Var）方法将观测资料处理转化成以动力模式为约束的泛函极小化问题，通过调整控制变量，使指定时间窗口内由控制变量得到的模式预报结果与实际观测资料之间的偏差达到最小。该方法在同化窗口内可以利用模式的切线性和伴随隐式地改变背景误差协方差，能够在某种程度上满足快速发展的天气过程。但是大部分业务中心的四维变分资料同化系统仍采用静态化的背景误差协方差矩阵模型来缓解背景误差协方差矩阵的维度问题，即矩阵维数远大于可用信息量。随着计算机科学的迅猛发展，维度问题可以进一步通过集合的方法缓解。集合四维变分资料同化就是基于这一目标通过构造多个能反映出背景误差协方差分布特征的样本集合来弥补可用信息量的不足。该方法目前已在ECMWF、Mete-France等业务中心实现业务化，为确定性四维变分资料同化系统提供流依赖背景误差协方差估计。简要介绍了集合四维变分资料同化方法的基本原理；其次以ECMWF为例，概述了四维变分资料同化系统的业务现状，重点阐述了系统在开发过程中需要解决的扰动、滤波、校正等一些关键技术；最后探讨集合四维变分资料同化系统目前存在的问题和未来可能的研究方向。

背景误差协方差矩阵，集合四维变分资料同化，扰动，流依赖

0 引言

大气作为一种典型的混沌系统，其发展和变化对初值非常敏感。Bjerknes等［1］曾把数值天气预报归结为一个典型的偏微分方程“初/边值问题”，指出要得到一个准确的天气预报必须要满足两个条件：一是数值预报模式能足够准确地模拟天气系统的演变过程，即对控制大气运动的物理规律有较好地描述；二是模式初值能足够精确地反映初始时刻的大气状态。所以好的初始场越来越被认为是整个数值预报领域的一个重要方面，初始场的精确性直接影响着数值天气预报的成败。资料同化方法应运而生成为估计大气状态的主要工具，它定义为利用一切可用的信息估计出尽可能精确的大气或海洋流体状态的过程［2］，这里的一切可用信息除了各种仪器设备观测到的气象观测数据以外，还包括模式对当前大气状态的短时预报（背景场）以及各自的误差信息［3］。1949年，Panofsky［4］第一次提出多项式插值的资料同化方法，即用一个低阶的多项式拟合分析区域内一小范围观测值，并尝试引入动力约束来协调风场和质量场，但由于区域多项式拟合不稳定，导致了分析场在拟合区域之间存在不连续。随后，资料同化方法经历了逐步订正［5］、最优插值［6-7］、三/四维变分［8-14］、卡尔曼滤波和集合卡尔曼滤波［15-16］等发展阶段。与其他资料同化方法相比，四维变分资料同化具有以下几个优势：1）在同化周期内可以同化不同时刻、不同地区、不同性质的各种气象资料而不需要对时间做近似［17］，通过设计非线性观测算子同化日益丰富的非常规资料［18-19］；2）将动力约束和资料约束纳入同一方程，由动力模式同化观测资料，反过来观测资料又优化动力模式的要素场和某些参数，这样在动力模式和观测资料之间建立起一种可信又可观的联系，不但同化效果好，而且同化的结果满足动力约束可直接输入模式进行预报［20］，可以在同化框架中比较自然地引入重力波控制和观测资料质量控制等技术；3）通过预报模式协方差矩阵、背景误差协方差矩阵、观测误差协方差矩阵有效控制各种误差，利用切线性和伴随模式隐式演变同化窗口内的预报误差协方差，这对于快速发展天气过程使用观测资料十分重要。基于以上因素，目前国际上大部分先进业务中心仍采用四维变分资料同化方法［21］。

背景误差协方差矩阵的精确定义是构建高水平资料同化系统的先决条件［22］。在四维变分资料同化框架中，背景误差协方差矩阵［23-25］的精度很大程度上决定了分析场的准确度，也就决定了数值预报的整体水平。一方面，背景误差协方差决定了背景信息在分析场中所占的权重，ECMWF评估各种信息在分析场中所占的贡献，同化系统得到的分析场信息大约只有15%来自同化窗口内的观测资料，其余85%是来自背景场［26］；另一方面，背景误差协方差对信息传播、信息光滑、平衡关系和流型结构建立等具有十分重要的作用，它决定了观测信息如何在模式空间上传递，并通过平衡关系将观测信息从一个变量传递给其他变量［23-24］。由于背景误差协方差矩阵对分析场的质量和预报技巧起关键作用，世界上各数值天气预报中心在建立变分同化系统的过程中都十分注重背景误差协方差模型的设计与统计参数的估计［11-13，27-30］。但背景误差协方差估计同时也是一项十分困难的工作。首先是背景误差协方差矩阵的维数特别巨大，如分辨率为25km的全球分析预报系统所需的背景误差协方差矩阵元素个数达到了1016以上；其次是我们无法获取大气环境的真实状态，也就是不可能直接得到完整的背景误差。在当前大多数的业务变分同化系统，尤其是全球变分同化系统中为了避免直接表示，通常是从物理、统计和计算效率等方面考虑构造简化和近似的背景误差协方差模型，如假定误差相关函数满足均匀性、各向同性和静态假设［23-24］。然后，在此基础上利用NMC、更新向量法、集合方法等统计出平衡关系的回归系数和误差相关函数［31-35］。但是这种处理忽略了背景误差协方差中的非均匀、各向异性和时变等特性，在处理锋面、台风等快速发展系统时尤为明显。

Fisher［22］和Fisher等［10］在全球四维变分资料同化系统中引入非正交小波基，有效解决了背景误差协方差的非均匀、各项异性表征问题。Zhang等［35］也在YH4DVAR系统中构建了基于小波的背景误差协方差模型，使该模型具备了各项异性、垂直相关和水平相关不可分，背景误差相关具有空间变化等特征。但是模型中的基于历史样本统计的方差项和水平相关函数仍然是一个时间常量，无法反映快速发展的天气系统。针对这一问题，Fisher［22］基于集合卡曼滤波的思想提出了集合同化概念，即同时运行多个相互独立的同化循环，每个循环的背景场、观测资料及SST输入是在控制的基础上叠加了服从各自误差分布的扰动，利用集合同化的样本就可以实时估计出背景误差方差，这即为集合四维变分资料同化方法（En4DVar）的雏形。En4DVar集合资料同化能够精确估计流依赖背景误差协方差，被认为是解决变分资料同化主要问题的有效方法，目前已经在国际上少数业务中心（如ECMWF［22，36-37］和法国气象局（Mete-France）［38-39］）实现了业务化，其业务实现对确定性四维变分系统的水平提升起到了至关重要的作用［36-37］。

本文简要介绍了En4DVar的基本原理及其在ECMWF和Mete-France的业务应用情况，归纳了En4DVar的业务化实现过程中面临的主要问题及相应的解决方法，旨在为其他业务中En4DVar的实现提供理论参考。随着研究的深入，En4DVar在不断被完善，应用也在不断被拓展，En4DVar后续的发展将在文章最后一部分讨论。

1 En4DVar的基本原理

En4DVar的目的是统计具有流型特征的背景误差协方差，即从前一时刻起报到分析时刻的短时预报与大气真值之间的误差协方差。背景误差的主要来源包括大气观测误差和积分模式误差。En4Dvar的思想源自扰动观测的集合卡曼滤波，在观测资料、海表温度（SST）和预报模式上叠加服从各自误差分布的扰动，构造多个能反映背景误差分布特性的分析—预报样本集合。利用数据同化集合（EDA）成员之间的差别可以统计所需的背景误差协方差。图1为En4DVar的扰动分析—预报系统示意图，相对于控制系统，扰动系统采用相同的资料同化模块，但是通过扰动预报模式的随机物理过程来隐式地扰动背景场，而对观测资料和侧边界的扰动则是直接进行。事实上，背景场也可以采用直接扰动的方式，后面2.2节将会详细讨论。

图1 En4DVar的扰动分析—预报系统示意图Fig. 1 Schematic illustration showing the analysis-forecast employed in En4DVar

短时的分析—预报可以被看成是一个弱非线性系统，对图1中的控制系统可以用下式表示：

其中，K表示增益矩阵，k代表分析—预报循环步序号，yk是观测矢量，表示第k时间步的分析场，是背景场，观测算子Hk将模式空间的变量转换到观测空间。M是预报模式，第k步的分析误差协方差矩阵和背景误差协方差矩阵满足下式：

式中，Rk和Qk分别表示观测、模式误差协方差矩阵，I是单位矩阵。

扰动系统在控制的基础上考虑了观测和模式的不完美性，并采用扰动的方式来构造一个可能输入状态。对于扰动系统，数学表达形式为：

式（3）与式（1）相减可得出添加的扰动满足：

对比控制和扰动系统的协方差矩阵演化式可以看出，如果在分析—预报循环步m时满足，则对于所有k≥m步都满足

上述理论表明，如果在控制系统的基础上叠加的观测和模式扰动分别是从真实的观测误差、模式协方差中独立取样得到，则基于扰动构造的分析—预报集合，利用成员间的差别可以统计出流依赖的背景误差协方差。事实上，在系统启动时刻，式（3）中的背景场还没加入扰动信息，即但是这并不影响收敛于的结论。可以将两个不同的经扰动系统演化后的差值定义为δ，两个不同的经控制系统演化后的差值定义为Δ，很容易证明δ与Δ是等价的，即收敛性不受初值状态和的影响。文献［40］研究了集合资料同化系统冷启动的收敛性，结果表明，对于500hPa的位势高度，全球集合离散度能在一周左右达到渐进值。

综上分析，如果观测、模式及其他输入量的扰动服从各自误差协方差的真实分布，则一周后的集合离散度能正确表征出分析、背景误差方差。这即为集合四维变分资料同化的基本理论依据。业务中为了使集合离散度处于收敛状态，En4DVar往往与业务高分辨率4D-Var系统同时不间断地往前循环推进。

2 En4DVar的发展及业务现状

目前，ECMWF、Mete-France等中心已成功实现了En4DVar的业务化运营。Fisher［22］首次提出了集合资料同化的思想，扰动短期预报场在生成过程中的各个环节，来表征短期预报场的主要误差来源，短期预报场之间的差别就具有背景误差的统计特征。ECMWF在文献［37，40］试验和研究基础上，于2010年6月在c36r2版本集成预报系统中正式引进En4DVar，一方面为业务高分辨率4D-Var系统提供流依赖的背景误差协方差的估计值，另一方面结合奇异向量技术为集合预报系统提供高质量的扰动初值。业务运营的En4DVar包含了10个扰动的成员和1个控制成员，为了避免引入垂直方向的插值误差，分辨率设为T399L91，在T399/ T95/T159的三重嵌套的最小化循环迭代的框架设计下，总计算量与T1279分辨率的业务确定性四维变分资料同化系统持平。En4DVar每天运行两次，选用相同的时间滑窗与观测资料集，模式的不确定性表征采用了集合预报系统中随机扰动参数化倾向方案。

整个业务流程如图2所示，En4DVar分析使用与高分辨率4D-Var相同的观测资料。控制成员不需要扰动观测资料，而其他En4DVar成员则需要对观测资料叠加一个随机方法抽样得到的满足均值为零、标准偏差等于观测误差的高斯分布扰动。由于大气运动矢量观测［41］和海表温度场［42］等具有高相关性的观测资料，在扰动中则需要考虑其相同性。在背景场扰动方面，ECMWF采用了SPPT模式扰动方案以提高En4DVar的离散度。统计得到的背景误差在应用到4D-Var之前需要经过滤波和校正处理来消除估计值中存在的随机误差和系统偏差。

图2 En4DVar流依赖背景场误差协方差的模拟流程Fig. 2 Process in generation of flow-dependent background error covariance

2012年，为了得到流依赖背景误差协方差，EDA的成员增至25个，并改进和引入了许多新的技术，如能量后向散射（SKEB）的模式扰动方法，方差估计值的小样本随机噪声滤波方法，方差估计值的系统偏差校正技术，以及流依赖背景相关函数的估计方法等等。

Mete-France业务实现的En4DVar系统与ECMWF相似，早期考虑到现有的计算条件，集合中只包含了6个成员，因此Mete-France在随机误差处理方面进行了大量研究。随着计算条件的改善，目前已经将集合成员数扩展到了10个。

2.1 观测资料和SSSTT扰动

根据观测之间相关性的强弱，将观测资料划分为强相关观测资料、弱相关观测资料。对于大部分不相关或相关性较小的观测资料，叠加的扰动服从均值为0，标准差为观测误差协方差的高斯分布。但是，扰动云导风资料时需要考虑资料之间的相关特性。由于SST是在原始观测基础上经过最优插值和二维变分等处理过程得到的二维格点场，在扰动时，除了原始观测误差外，还需要考虑处理过程中的误差特性，且扰动也需要遵循SST的遥相关特性。

受高度指定、相似云结构追踪、质量控制过程［43］等因素影响，同步卫星的云导风观测具有很强的时间、空间相关性［44-45］，目前大部分风产品的分辨率在160km或更高，因此只能诊断该分辨率以下的相关性，而很难表征出更细的相关特征。Bormann等［41］基于一年的云导风（AMVs）-无线电探空数据集，在假定探空仪观测误差不相关的条件下，利用密集的探空观测网络研究了AMVs随机误差的空间相关性。其结果表明，云导风的显著相关距离约为800km，且相关距离对于不同卫星、不同通道和垂直层相关距离变化不大。其中热带地区的相关性大于其他地区，相关性具有各向异性。北半球高空的风分量误差年平均约为2.7～3.5m/s，冬季的误差最大。Bormann等［41］分别采用各向同性（利用站点距离分组）或各向异性（利用N-S和E-W分离进行分组）相关函数将相关的数据以一种统计的方式外插到0来估计空间相关的AMV误差大小，外插的相关函数在0点位置将探空资料方差划分为空间相关和空间不相关部分。相关部分对应AMV的观测误差，而空间不相关部分则由不相关的AMV误差、探空观测误差、探空（单点测值）与AMVs（区域平均）之间的失配误差构成。其中各向同性的相关函数可采用以下形式：

式中，L表示R的长度尺度，r表示两点之间的大圆距离，R0表示r=0时的R的大小。

海洋具有多种天气气候意义的特性，它在地-气系统热量平衡及水分平衡中具有重要作用。其中，SST不仅是海洋表面物理状态的重要参数，还是影响大气环流及长期天气变化的重要因素。在4D-Var中，SST作为侧边界输入信息，其精度或者误差分布对同化效果影响较大，同时SST又具有遥相关特性［46］，因此对SST扰动的过程也必须考虑其相关性。目前，ECMWF业务中采用的是集合预报中SST扰动方案［42］，该扰动方案通过构造两组扰动量来表征SST的不同误差来源：第一组扰动量的构造方法是，根据Reynolds等［46］的最优插值方法得到的SST分析场周平均和二维变分SST分析场的周平均差在1985—1999年的统计量来构造表征SST产品中的典型误差分布的扰动量；第二组SST扰动构造则是利用Reynolds等［46］二维变分SST分析场与其一周均值的差来构造扰动。第一组扰动代表了SST分析场中的不确定性，而第二种扰动则主要考虑到了NCEP的SST是周平均产品。这两组扰动采用随机选取的方法叠加到SST分析场上。扰动系数由海表的1线性衰减到40m海深处的0。

2.2 预报模式扰动

En4DVar采用的预报模式并不是真正的“完美模式”，而是用到了很多次网格参数化过程和随机方法，因此也需要开展模式扰动方案来表征模式中的误差特征。模式误差来源之一是参数化方案中缺少对次网格物理过程变率的描述。对于数值模式中各物理过程的参数化方案，它的整体作用表现为在控制方程中采用某一倾向项描述次网格物理过程的贡献，该描述是一确定性的结果并且依赖于网格尺度的物理量，因此忽略了物理量通量具有统计意义振荡的特性以及网格尺度运动与次网格尺度运动之间的相互作用。

模式误差的另一个来源则是参数化方案及模式积分方案本身的原理导致了系统性的动能缺失，从而使得模式大气的动能谱与实际大气不符。例如，模式积分中为保证计算稳定性而采用半拉格朗日平流方案并引入水平耗散项，这往往导致过强的能量耗散；在深对流参数化方案中，没有合理描述出对流产生的动能向平衡流场传输并激发重力波生成这一物理过程。所以，这一模式误差来源的主要影响在于没有描述大气动能的升尺度传播特性，导致系统性的动能谱偏差。ECMWF的En4DVar业务系统中根据如上的模式误差来源，同时引用了集合预报系统中的两套扰动方案：参数化倾向随机扰动法（SPPT）表征已有参数化方案中存在的不确定性；随机后向散射法（SKEB）表征模式中未被参数化方案描述而缺失的物理过程。

BMP方法［47］是SPPT的原始版本，最早应用在ECMWF的集合预报系统中，用来表征数值模式已有物理过程参数化方案中存在的不确定性。对于模式中的任一预报量，其预报方程如下：

式中，A表示模式网格尺度运动（非参数化部分）对预报量倾向的贡献。P表示次网格物理过程参数化对预报量倾向的贡献，为描述这一倾向分量的不确定性，在等式右边叠加与P有关的随机强迫项，故上式可改写成：

SPPT的原理与BMP方法一致，区别在于对随机数r的处理，即对所用变量采用相同的r，且在近地表和平流层引入调整因子来调整扰动的幅度。另外，为了解决BMP方法出现的不连续现象，SPPT方法在谱空间中利用谱系数构造随机扰动，使得扰动随时间和空间的变化都非常平滑。SPPT实现后改进了集合预报中降水分布的预报效果［48］。

SKEB［47，49］方法认为与随机中小尺度系统运动联系的物理过程包含了局地深对流产生的动能向平衡流场传输并激发重力波、水平耗散等过程，它将影响尺度位于模式截断波数附近的动能传输，最终使得网格尺度运动的动能谱发生改变，故最终的影响效果为风场分量的倾向。因此，随机物理过程扰动作为一强迫项作用于速度场，引入有效流函数强迫项：

SPPT和SKEB都是通过扰动模式来隐式表征背景场的不确定性，另一种显式表征方法是直接将扰动叠加到背景场上，这种方法称为XB方法［51］。扰动定义为，其中，扰动幅度调整因子是一个与纬度模式层l和模式变量x有关的函数，则是从高维高斯分布（0，B）中取样得到的随机扰动，B表示背景误差协方差矩阵。的大小为高分辨率业务信息向量方差和集合方差（这里的集合仅扰动观测）的差值的三周统计平均。ECMWF对离散度诊断后得出XB方法能得到比物理过程参数化方法更大的离散度，因此不再需要进行集合离散度的校正处理，这一改进在中高纬地区更为明显。此外，应用XB方法的En4DVar来为ECMWF的集合预报系统提供扰动初值，能够得到最大的离散度，使得中期集合预报系统更为稳定。相对于原有的集合预报系统3和7d的预报技巧有了小幅度的提升。

2.3 方差系统偏差校正

En4DVar通过使用扰动观测，SST和背景场分别表征输入到分析—预报系统的初值、边界值、背景场等主要误差来源。但是在观测误差统计，SST扰动和模式误差参数化中的任何缺陷或近似，以及其他未知的不确定项（如陆表过程）都会引起En4Dvar取样方差为分析和背景场预报的次最优估计值。这种类型的估计误差无法通过增加集合成员个数方法来消减，且转化为En4DVar取样方差和真实分析/预报误差的系统差别。

一种简单的处理方式是将En4DVar平均离散度值乘以一个定值缩放因子来减少这种系统偏差。但是En4DVar方差的系统误差具有复杂的空间和时间结构，并不能简单地通过乘以一个全局的缩放因子来表征。ECMWF采用Leutbecher［49］的方法，即利用集合离散度-误差关系定量诊断系统偏差。对于一个理想的En4Dvar系统，离散度-误差曲线应当位于对角线位置，即：

其中，N是集合成员个数，y是真实值，实际中以高分辨率业务分析场代替。等式左边表示集合方差，右边为集合误差方差（RMSE）。而实际的En4Dvar系统并不满足上式，即离散度—误差曲线斜率并不为1，其大小反映了集合离散度的条件偏移程度。处理方法按照离散度—误差差别的带状分布特征，较为精细的将全球划分为N（30°—90°N），S（30°—90°E），T（30°E—30°N）三个区域，对每个区域应用不同的膨胀系数。考虑到离散度散度-误差诊断系数的日变化和周变化不大，但有一明显的季节性漂移，大小可利用最近5d的集合样本通过离散度-误差曲线以在线的形式来确定。

2.4 方差随机误差滤波

受计算资源的限制，集合成员个数限制在o（10）～o（102），由这种小集合统计得到的B矩阵，包含了样本噪声，对B矩阵的统计精度有很大的影响，因此在使用前需要引入滤波工具消除样本噪声。文献［32， 39， 51-52］研究了长距离的假非零值。Raynaud等［53］对相关性进行了试验研究，结果表明，最优局域空间平均滤波能提高B的计算精度但效果并不明显，且低通滤波器的截断波数与变量场、垂直层、背景误差特征尺度有关，需要花费大量的工作进行最优化设置。Raynaud等［54］在此基础上进行了改进，通过计算信号、噪声的特征尺度长度，能自动计算最佳截断波数。由此发展成的客观滤波器被应用到ECMWF的集合资料同化系统中。

形式上可将EDA的计算误差分为随机项和系统项：

对于一组相互独立、服从高斯分布的N成员扰

可以推导出样本噪声协方差矩阵的矩阵形式为：

式中，n是总谱波数，Ntrunc是滤波器的截断波数。在谱空间应用此低通滤波，相当于在格点空间的误差方差添加了一个加权平均，使得较大尺度信号通过，而较小尺度样本噪声被过滤。但该方法也存有以下局限：1）它是基于样本估计的误差协方差矩阵，而真实的误差协方差矩阵是未知的；2）所采用的参数不能适用于所有变量/模式层；3）不适用于处理EDA方差能量谱的小变化。为了克服以上方法的局限性，Montmerle等［34］对以上方法进行了修正，利用同一个EDA系统的两个不同集合来计算EDA背景场标准偏差噪声能量谱：

谱形式的噪声滤波可以简单的拓展到小波形式。小波转化的主要特征是，信号与一组径向基函数的卷积、径向基函数的选择受谱带限制，对于一些连续的截断波数N1，…，NK，Legendre 转化满足：

从一组10成员EDA事件，可以计算小波空间中相关系数样本，如对于每个小波分解的波带计算空间相关系数场。小波滤波标准偏差就能通过原标准偏差和取样相关系数在小波空间中的卷积导出：

不同尺度上的截断波数可以采用Donoho方法来自动最优确定［55］，图3对应2013年8月2日09：00 UTC第91模式层上涡度场背景误差的集合估计值、应用Donoho方法经小波滤波和谱滤波方法得到结果。其中心最大值对应于第九号台风飞燕的中心位置，大小分别为8.31×10-5、8.21×10-5和6.59×10-5。对比可以看出，应用小波方法能够改善谱滤波在空间平滑过程导致部分关键信息被抹平的局限性，且由于Donoho方法无需应用式15统计噪声能量谱，因此可以极大提高滤波效率。

2.5 流依赖相关函数统计

图3 2013年8月2日09:00 UTC第91模式层上涡度场背景误差（单位：10-5s-1）（a）、（b）、（c）分别为10个样本得到的集合统计值、小波和谱方法滤波结果Fig. 3 Standard deviations of voracity at model level 91, for 0900 UTC on 2 August 2013 (unit: 10-5s-1). (a) Corresponding to raw estimates from 10 member ensemble; (b) to the filtered results with wavelet and (c) with spectral method

相关函数C对应的是背景误差协方差矩阵B中的非对角元素，表示变量内部和变量之间在不同位置、不同高度的相关性。由于C的个数是方差元素个数的二次方倍，因此对C的精确估计远大于方差估计。Fisher［22］利用不同时次的样本估计得到具有尺度相关、位置相关的各向异性相关函数，但是其统计是基于历史不同时次的集合样本，因此得到的是一个气候态的定常量，忽略了其流依赖的特性。Varella等［56］研究表明，具有弱流依赖特性的时变相关函数能够更加精确地表征出隐含的相关结构，这一优势在活跃天气系统的影响区域中更为明显。

根据样本个数的不同，对流依赖相关函数的处理也不同。当样本个数十分有限时，为了消除因样本不足导致的噪声或虚假相关，需要构造合适的空间局地化函数来消除或减少远距离虚假相关；另一种方法是构造出足够多的样本个数，减少小样本估计值中含带的远距离虚假相关。

Houtekamer等［52］于2001年首次引入局地化矩阵Cloc与原估计值进行相乘来实现局地化，对流依赖相关函数即：

式中，Cloc是变化的相关性矩阵，大小与相同，Cloc的形式直接决定了流依赖相关函的局地化的效果。一种简单局地化处理方式是，假定变量的自相关强度与距离成反比，Cloc为一个高维的高斯函数矩阵。另一种处理思路将低分辨率的局地化矩阵通过谱变换作用到高维的流依赖相关函数估计值矩阵上。其优点是，如果局地相关性在谱空间中是水平均匀的，则谱形式的局地化矩阵Clocs即为简单的对角矩阵。

实践表明，在没有空间局地化和滤波处理的前提下，至少需要500～600个扰动样本才能粗略统计得到小波B模型。如此庞大的样本个数是现阶段、以及未来几年都无法承受的。对此有两种解决思路：1）通过滞后方法实时估计B；2）采用混合的形式构造B。滞后法实时估计B的思想是结合现有的和过去12d滑动窗口内的EDA背景预报成员构成一个包含600个样本的超大EDA集合，用于统计流依赖相关函数C。该方法尽管得到足够的统计样本，但是由于滑动窗口过长，统计得到的相关函数C中的流依赖信息量过少，很难及时反映快速发展的天气系统。且如果业务EDA成员个数少，这一情况将更加严重。混合方法则是利用最新的200个集合成员表征出误差逐日变化的特征，另外再从全年的EDA样本取样出400个成员表示误差的气候态特性。ECMWF试验表明上述方法都能有效改进现有的静态相关函数C，且对同化和预报具有正效果［58］。

3 存在问题及展望

集合四维变分资料同化系统业务化以后，显著提升了确定性和集合预报水平。其成功主要基于四个方面的因素：拥有一个非常可靠的成熟的集合四维变分资料同化系统；先进的小波背景误差协方差模型；精确的切线性、伴随模式以及高精度的观测算子。这些因素使得背景误差能在同化窗内的表征非常精确，而避免采用空间局地化技术来消除协方差的噪声。但是集合资料同化系统仍存在一些尚未解决的问题：如忽略了中高层质量-风平衡约束对预报的影响。同时出于实际业务化考量，ECMWF和Mete-France简化了En4Dvar中误差表征、随机误差处理、系统误差校正等问题，随着计算条件的逐步提升，对这些问题的处理将会更加精细化。

常规观测资料扰动采用均值为0、标准偏差等于观测误差的高斯型分布来模拟。但是对于占据同化系统观测资料来源主体（约占90%）的卫星观测资料而言，部分类型的观测误差具有空间或时间相关性，扰动时必须遵循这一相关性。目前上述业务中心仅考虑了具有强相关的云导风卫星观测资料，并假设其他卫星观测资料是不相关的。

随机误差的滤波处理能减少背景误差方差估计值中的采样噪声，提高精确度，其作用十分重要。业务中采用的统计滤波都属于均匀滤波，即作用在谱系数上的滤波系数仅与波数（或尺度）有关。在格点空间相当于全球相同尺度的平滑范围和平滑强度是一样的。这种处理有可能会抹去尺度较小的局地特征信号，降低集合方差估计值的有效分辨率，如T399分辨率下，10成员的有效分辨率仅为T70。因此有必要开展更精细化的非均一滤波方法。非均一滤波系数具有各向异性，大小不仅与尺度有关，还受方向和位置影响，符合背景误差协方差各向异性的这一既定事实。文献［57］利用扩散方程构造了一种非均一滤波器，通过建立起滤波系数与尺度函数之间的关系来局地调整滤波的强度。由于ECMWF的4D-Var系统采用的背景误差协方差模型是基于小波构造的，另一种高效的方法是利用小波的多尺度、多分辨率的刻画优势，同时精确表征全局的尺度信息和局地的信息。同时将信号和与位置x有关长度尺度信息L（x）转换到小波空间，构造非均一小波滤波器

集合方法统计得到的离散度一般偏小，一种处理的方法是乘上一个膨胀系数来调整离散度，这个膨胀系数对全球而言是一致且不变的。改进后可以按照离散度的带状分布特征，较为精细的将全球划分为N（30°—90°N），S（30°—90°E），T（30°E—30°N）三个区域，对每个区域应用不同的膨胀系数，考虑到离散度的时变特性，膨胀系数可利用最近5d的集合样本通过离散度-误差曲线以在线的形式来确定。划分区域能够较为精确的调整每个区域的离散度，但是在区域之间的边界区域会存在不连续情况。

2013年，ECMWF将En4DVar的成员个数由原来的10个增加到了25个，使得到的背景误差方差更加精细化。其另一方面的作用是可以利用25个集合成员估计流依赖的背景误差相关函数，结合背景误差方差组成一个全流依赖的背景误差协方差估计值。ECMWF计划将估计值与静态的背景误差协方差矩阵按照30%和70%的权重线性组合。另一种解决思路是采用滞后法构造足够多的统计样本统计流依赖的背景误差协方差。

目前，ECMWF的业务资料同化系统已经从原来的T799L91提升到了T1279L139。随着模式和资料同化的空间分辨率进一步提升，原有相关长度尺度的各向同性假设已不再适用，现有背景误差协方差模型面临的流依赖平衡关系估计问题将更加严峻。无论是集合四维变分资料同化（En4DVar），还是目前盛行的四维集合变分资料同化（4DEnVar），都需要进一步研究集合信息引入到4D-Var中的效率。另一个面临的关键问题是En4DVar本质上属于蒙特卡罗方法，其分辨率与扰动是否能表征高分辨率资料同化系统的误差密切相关。即使En4Dvar的外层循环分辨率从T399升到T699，仍只有目前4Dvar分辨率的一半。它能够考虑大部分不确定性信息，降低参数化表征模式误差的需求，但是以4D-Var为基本模块的En4Dvar是一个计算量非常大的系统，未来研究将主要集中在提高En4Dvar的效率，如改进估计值中随机误差、系统误差的处理方法，利用流依赖集合离散度优化观测资料质量筛选和质量控制决策等。

［1］ Bjerknes V， Hesselberg T. Dynamic meteorology and hydrographs. PartⅡ： Kinematics. New York： Carnegie Institute Press， 1911.

［2］ Talagrand O. Assimilation of observations， an introduction. J Meteor Soc Japan， 1997， 75： 191-209.

［3］ Ham J. Mesoscale predictability and background error covariance estimation through ensemble forecasting： Texas： Texas A&M University College Station Press， 2002.

［4］ Panofsky H. Objective weather map analysis. J Met， 1949， 6：386-392.

［5］ Gilchrist B， Cressman G P. An experiment in objective analysis. Tellus A， 1954， 6（4）： 309-318.

［6］ Bergthorsson P， Does B. Numerical weather map analysis. Tellus A，1955， 7（3）： 329-340.

［7］ Rutherford I D. Data assimilation by statistical interpolation of foercast erorr fi elds. J Atmos Sci， 1972， 29（5）： 809-815.

［8］ Sasaki Y. Numerical variational analysis with weak constraint and application to surface analysis of severe storm gust. Mon Wea Rev，1970， 98（12）： 899-910.

［9］ Laroche S， Gauthier P. A validation of the incremental formulation of 4D variational data assimilation in a nonlinear barotropic flow. Tellus A， 1998， 50（5）： 557-572.

［10］ Fisher M， Andersson E. Developments in 4D-Var and Kalman filtering， ECMWF Tech Memo No 347， 2001. （available from European Centre for Medium-Range Weather Forecasts， Shinfi eld Park， Reading， Berkshire RG2 9AX， UK）.

［11］ Zhao J， Song J Q， Li Z J. Distributed parallelization of a global atomospheric data objective analysis system. Adv Atmos Sci， 2003，20（1）： 159-163.

［12］ 曹小群， 黄思训， 杜华栋. 变分同化中水平误差函数的正交小波模拟新方法. 物理学报， 2008， 57（3）： 1984-1989.

［13］ 曹小群， 黄思训， 张卫民， 等. 区域三维变分同化中背景误差协方差的模拟. 气象科学， 2008， 28（1）： 8-14.

［14］ 卢风顺， 宋君强， 朱小谦. WRF三维变分同化并行程序性能分析.计算机工程与科学， 2008， 29（11）： 149-151.

［15］ Daley R. Atmospheric data analysis. London： Cambridge university Press， 1993.

［16］ Anderson J L. An ensemble adjustment Kalman filter for data assimilation. Mon Wea Rev， 2001， 129： 2884-2903.

［17］ 廖洞贤， 王两铭. 数值天气预报中的若干新技术. 北京： 气象出版社， 1995.

［18］ 潘宁， 董超华， 张文建. 变分同化及卫星资料同化. 气象科技，2001， 29（2）： 29-36.

［19］ Lindskog M， Salonen K， Järvinen H， et al. Doppler radar wind data assimilation with HIRLAM 3DVAR. Mon Wea Rev， 2004， 132（5）：1081-1092.

［20］ Li Z， Navon I M， Zhu Y. Performance of 4D-Var with diff erent strategies for the use of adjoint physics with the FSU global spectral model. Mon Wea Rev， 2000， 128（3）： 668-688.

［21］ Fisher M， Courtier P. Estimating the covariance matrices of analysis and forecast error in variational data assimilation. ECMWF Tech Memo No 220， 1995. （available from European Centre for Medium-Range Weather Forecasts， Shinfield Park，Reading， Berkshire RG2 9AX， UK）.

［22］ Fisher M. Background error covariance modelling. Seminar on Recent Development in Data Assimilation for Atmosphere and Ocean. ECMWF， 2003.

［23］ Bannister R N. A review of forecast error covariance statistics in atmospheric variational data assimilation. I： Characteristics and measurements of forecast error covariances. Q J R Meteorol Soc，2008， 134（637）： 1951-1970.

［24］ Bannister R N. A review of forecast error covariance statistics in atmospheric variational data assimilation. II： Modelling the forecast error covariance statistics. Q J R Meteorol Soc， 2008，134（637）： 1971-1996.

［25］ Lorenc A C. Analysis methods for numerical weather prediction. Q J R Meteorol Soc， 1986， 112（474）： 1177-1194.

［26］ Cardinali C， Rukhovets L， Tenenbaum J. Jet stream analysis and forecast errors using GADS aircraft observations in the DAO，ECMWF， and NCEP models. Mon Wea Rev， 2004， 132（3）： 764-779.

［27］ 曹小群， 宋军强， 张卫民， 等. 一种基于复数域微分的资料同化新方法. 物理学报， 2013， 62（17）： 170504-170504.

［28］ 王舒畅， 李毅， 张卫民， 等. 资料同化中的数字滤波弱约束试验及分析. 物理学报， 2011， 60（9）： 99203-099203.

［29］ 赵延来， 黄思训， 杜华栋. 基于变分方法的有限区域风场分解与重构I： 理论框架和仿真实验. 物理学报， 2013， 62（3）： 39204-039204.

［30］ 朱孟斌， 张卫民， 曹小群. GPS掩星一维弯曲角算子在四维变分资料同化系统中的实现方法研究. 物理学报， 2013， 62（18）：189203-189203.

［31］ Hamill T M， Snyder C. Using improved background-error covariances from an Ensemble Kalman Filter for adaptive observations. Mon Wea Rev， 2002， 130（6）： 1552-1572.

［32］ Buehner M， Charron M. Spectral and spatial localization of background-error correlations for data assimilation. Q J R Meteorol Soc， 2007， 133（624）： 615-30

［33］ Carrier M J， Ngodock H. Background-error correlation model based on the implicit solution of a diffusion equation. Ocean Modelling， 2010， 35（1）： 45-53.

［34］ Montmerle T， Berre L. Diagnosis and formulation of heterogeneous background-error covariances at the mesoscale. Q J R Meteorol Soc， 2010， 136： 1408-1420.

［35］ Zhang W M， Cao X Q， Xiao Q N. Variational data assimilation using wavelet background error covariance： initialization of typhoon KAEMI （2006）. J Trop Meteor， 2010， 16： 333-340.

［36］ Bonavita M， Isaksen L， Hólm E. On the use of EDA background error variances in the ECMWF 4D-Var. Q J R Meteorol Soc，2012， 138（667）： 1540-1559.

［37］ Bonavita M， Raynaud L， Isaksen L. Estimating background-error variances with the ECMWF Ensemble of Data Assimilations system： the eff ect of ensemble size and day-to-day variability. Q J R Meteorol Soc， 2010， 137： 423-434.

［38］ Ştefănescu S E， Berre L， Pereira M B. T e evolution of dispersion spectra and the evaluation of model differences in an ensemble estimation of error statistics for a limited-area analysis. Mon Wea Rev， 2006， 134（11）： 3456-3478.

［39］ Berre L， Pannekoucke O， Desroziers G， et al. A variational assimilation ensemble and the spatial filtering of its errorcovariances： increase of sample size by local spatial averaging. Proc ECMWF Workshop on Flow-Dependent Aspects of Data Assimilation： 2006， 151-168.

［40］ Isaksen L， Fisher M， Berner J. Use of analysis ensembles in estimating fl ow-dependent background error variances. ECMWF Tech Memo No 492， 2006. （available from European Centre for Medium-Range Weather Forecasts， Shinfield Park， Reading，Berkshire RG2 9AX， UK）.

［41］ Bormann N， Saarinen S， Thepaut J， et al. The spatial structure of observation errors in atmospheric motion vectors from geostationary satellite data. Mon Wea Rev， 2003， 131： 706-718.

［42］ Vialard J. Vitartf F， Balmaseda M A. el al. An ensemble generation method for seasonal forecasting with an ocean-atmosphere coupled model. Mon Wea Rev， 2005， 133： 441-453.

［43］ Holmlund K. The utilization of statistical properties of satellitederived atmospheric motion vectors to derive quality indicators. Wea Forecasting， 1998， 13： 1093-1104.

［44］ Ingleby N B. The statistical structure of forecast errors and its representation in the Met. Office Global 3-D Variational Data Assimilation Scheme. Q J R Meteorol Soc， 2001， 127（571）：209-231.

［45］ Bouttier F， Kelly G. Observing-system experiments in the ECMWF 4D-Var data assimilation system. Q J R Meteorol Soc，2001， 127（574）： 1469-1488.

［46］ Reynolds R W， Rayner N A， Smith T M， et al. An improved in situ and satellite SST analysis for climate. J Climate， 2002， 15（13）：1609-1625.

［47］ Berner J， Ha S Y， Hacker J P， et al. Model uncertainty in a mesoscale ensemble prediction system： stochastic versus multiphysics representations. Mon Wea Rev， 2010， 139（6）： 1972-1995.

［48］ Palmer T N， Buizza R， Doblas-Reyes F， et al. Stochastic parametrization and model uncertainty. ECMWF Tech Memo No 598， 2009. （available from European Centre for Medium-Range Weather Forecasts， Shinfi eld Park， Reading， Berkshire RG2 9AX， UK）.

［49］ Berner J， Shutts G J， Leutbecher M， et al. A spectral stochastic kinetic energy backscatter scheme and its impact on flowdependent predictability in the ECMWF ensemble prediction system. J Atmos Sci， 2009， 66（3）： 603-626.

［50］ Cardinali C， Zagar N， Radnoti G， et al. Representing model error in ensemble data assimilation. ECMWF Tech Memo No 726， 2014.

［51］ Houtekamer P L， Mitchell H L. Data assimilation using an ensemble Kalman fi lter technique. Mon Wea Rev， 1998， 126（3）： 796-811.

［52］ Houtekamer P L， Mitchell H L. A sequential ensemble Kalman filter for atmospheric data assimilation. Mon Wea Rev， 2001，129（1）： 123-137.

［53］ Raynaud L， Berre L， Desroziers G. Spatial averaging of ensemblebased background-error variances. Q J R Meteorol Soc， 2008，134（633）： 1003-1014.

［54］ Raynaud L， Berre L， Desroziers G. Objective fi ltering of ensemblebased background-error variances. Q J R Meteorol Soc， 2009，135（642）： 1177-1199.

［55］ Donoho D， Johnstone I. Ideal spatial adaptation via wavelet shrinkage. Biometrika， 1994， 81： 425-455.

［56］ Varella H， Berre L， Desroziers G. Diagnostic and impact studies of a wavelet formulation of background-error correlations in a global model. Q J R Meteorol Soc， 2011， 137： 1369-1379.

［57］ Raynaud L， Pannekoucke O. Heterogeneous fi ltering of ensemblebased background-error variances. Q J R Meteorol Soc， 2012， 138：1589-1598.

Research Progress in Ensemble Four-Dimensions Variational Data Assimilation

Liu Bainian1，2， Huang Qunbo1， Zhang Weimin1， Cao Xiaoqun1， Zhao Jun1， Zhao Yanlai1
（1 Academy of Ocean Science and Engineering， National University of Defense Technology， Changsha 410073 2 College of Computer， National University of Defense Technology， Changsha 410073）

Accurate background error covariance is the foundation for all advanced data assimilation systems. For four dimensions data assimilation （4D-Var）， assimilating the observation data is converted to a question of cost function minimization which is constricted by atmosphere dynamic model. By adjusting the control vectors， the distance between model trajectory and real time observations reached its minimal value over whole assimilation time window. As background error covariance evolves according to the adjoint and tangent linear model， it can adapt to rapid development weather. However， most of operational 4D-Var systems still adopt simi-climatic background error covariance model compromised by huge dimensionality， which can't be exactly defi ned with all available information. As the rapid development of computer science， the problem of dimensionality can be released by ensemble method. Ensemble four dimensionality data assimilation （En4DVar） employed several independent perturbed analysis forecast cycles to remedy the limited information synchronously. In this scheme， fl ow-dependent background error covariance can be estimated from the differences between ensemble members. Several famous numeric prediction centers， such as ECMWF，Mete-France， adopted it to provide fl ow-depended background error covariance for the high-resolution determined 4D-Var system. In this thesis， the basic theory of the En4DVar method is demonstrated briefl y， followed by a description of currently application at ECMWF， and focusing on the disturbing， filtering， calibration as well as other key techniques for helping to improve the precision of estimates. The last part presents an investigation of some issues in current operation and possibly future research fi elds in the En4DVar.

background error covariance， ensemble， four-dimensional data assimilation， disturb， fl ow-dependent

10.3969/j.issn.2095-1973.2016.05.002

2015年5月21日；

2015年7月11日

刘柏年（1985—），Email： bnliu@nudt.edu.cn

资助信息： 国家自然科学基金项目（41375113，41105063，41475094，41305101）