背景误差水平相关结构对四维变分资料同化的影响研究

奥鹏网院作业

背景误差水平相关结构对四维变分资料同化的影响研究*
龚建东 张 林 王金成

国家气象中心/数值预报中心，北京，100081

摘 要 为考察GRAPES 全球四维变分同化（4DVar）的分析增量在谱空间的时间演变特征，分析当同化时间窗起始时刻与终止时刻背景误差水平相关特征明显不一致时对分析与预报造成的影响，对GRAPES 全球4DVar 的背景误差水平相关采用二阶自回归模型（SOAR）、集合资料同化生成扰动样本估计的水平相关模型以及基于这两者的背景误差谱空间融合模型进行比较。结果表明，SOAR 的分析增量在20 波以上的天气尺度波动的分析信息明显不足，而将集合资料同化样本所计算的水平相关的功率谱方差与SOAR 功率谱方差进行融合，水平相关特征呈现出多尺度水平相关的特点，可以更好地吸纳观测信息，显著改善北半球形势场、温度与风场预报效果，南半球也有改善，对赤道地区的影响中性。表明研究发展的融合水平相关方案合理、实用。

关键词 GRAPES 四维变分同化，背景误差水平相关，谱空间方差融合，数值试验

1 引言
四维变分资料同化（4DVar）是数值天气预报中的重要技术方法（Talagrand，et al，1987），在世界主要业务中心得到广泛应用（Rabier，et al，2000；Gauthier，et al，2007；Rawlins，et al，2007），中国自主发展的全球/区域分析与预报系统（GRAPES）中全球资料同化系统采用了增量四维变分同化方法（Zhang，et al，2019）。不断改善4DVar 分析质量是资料同化技术研究的重要方向，包括观测资料质量控制和资料的合理同化应用（薛纪善，2009）、背景误差协方差的合理估计（Bannister，2008a，2008b）、复杂物理过程参数化方案的线性简化处理（Lopez，et al，2005；龚建东等，2019）、切线性和伴随模式调优与下降算法的计算效率（刘永柱等，2017；张林等，2017）、分析场的动力平衡特征维持与重力波噪音控制（Gauthier，et al，2001；刘艳等，2019）等方面。

不同于早期业务使用的不考虑观测与背景场时间对应的三维变分同化（3DVar）方法，4DVar 方法则考虑了每个观测与背景场在时间上的多时刻对应关系，4DVar 寻求在整个有观测资料的时间期间（称为同化时间窗），观测与分析场的总距离以及同化时间窗起始时刻分析场与背景场的距离之和达到最小。相较于以背景场为初值的积分轨迹，起始时刻分析场的时间积分轨迹在同化时间窗内与观测值拟合得更好（Lorenc，et al，2005），分析场积分轨迹与观测拟合的改善成为4DVar 的一个重要特征。此外，背景误差协方差在资料同化中起关键作用，3DVar 使用静态的、不随时间变化的气候背景误差协方差，而4DVar 在极小化循环迭代过程中通过切线与伴随模式积分在同化时间窗内隐式演变背景误差协方差（Lorenc，et al，2005）。合理表达背景误差协方差的时间演变，对4DVar 分析效果的改善也很重要。

变分同化系统按照统计模型构造气候背景误差协方差，其可以表述为背景误差方差、背景误差水平与垂直相关等要素。水平相关一般采用二阶自回归模型（SOAR），垂直相关多采用预报样本来估计，这方面的研究工作开展较多（Bannister，2008a，2008b；王瑞春等，2015）。在估计背景误差协方差方面，Fisher（2003）、Isaksen 等（2010）提出使用扰动观测和扰动海温的方法建立集合资料同化系统（Ensemble of Data Assimilation，EDA）生成统计样本，来更为准确地估计背景误差协方差的特征。

在GRAPES 全球4DVar 业务系统中，起始时刻背景误差水平相关模型采用二阶自回归模型（张华等，2004），SOAR 仅能表述水平相关结构的粗略特征，并不能准确表述背景误差水平相关的实际特征。当起始与终止时刻的水平相关特征不一致时会对4DVar 分析和预报带来怎样的影响，现有研究工作给出清晰结论的还不多。本研究的重点是，在分析场积分轨迹与同化时间窗内观测资料拟合更好的基础上，进一步考察分析增量在同化时间窗的功率谱演变特征，并利用GRAPES 全球4DVar 系统来分析当起始与终止时刻背景误差水平相关特征差异对分析与预报造成的影响。

2 4DVar 分析增量的功率谱分布特征
2.1 公式表达及误差特征分析
GRAPES 全球4DVar 对模式初值的最优估计可以表述为求解以下目标函数的极小值（Zhang，et al，2019）

pagenumber_ebook=116,pagenumber_book=989
pagenumber_ebook=116,pagenumber_book=989
pagenumber_ebook=117,pagenumber_book=990
由式（6）可见，在 t0时通过统计模型给出背景误差协方差 B0，而在同化时间窗的其他时刻，背景误差协方差 Bi由切线模式的动力约束来控制。

Daley（1991）指出，背景误差水平相关变换的具体定义将直接影响分析增量的谱响应特征。式（2）表明，t0时分析增量（δXa）的谱响应特征由物理变换（P）、非平衡部分的背景误差方差（Σu）以及背景误差水平相关变换（U）决定。首先分析这三者对背景误差水平相关的影响。由于 Σu为纬向平均值，代表大尺度特征，其对天气尺度分析增量的谱响应的影响可以忽略，对于物理变换（P），非平衡物理量与物理全量的功率谱差异较大，以地面气压为例，旋转风解释的气压平衡部分占到全气压的85%—90%，而旋转风解释辐散风的约10%（Derber，et al，1999，图2）。龚建东（2007）研究表明，特征长度的变化与背景误差功率谱沿总波数的分布特征相关，主要由次天气尺度（20—60 波）到中尺度波（大于60 波）的误差功率谱的斜率特征决定。非平衡部分的功率谱与物理全量的误差功率谱斜率特征相近，意味着水平相关特征相近，Derber 等的研究也表明非平衡气压与气压的水平相关特征接近（Derber，et al，1999，图8）。这表明物理变换（P）不是主导影响因素，最主要的影响来自背景误差水平相关变换（U）。而ti时分析增量δXia由切线模式M0,i的动力演变特征和背景误差水平相关变换（U）特征共同决定，M0,i也将影响分析增量的水平谱响应特征。式（5）和（6）同样表明，t0时背景误差 εb的谱响应特征也由背景误差水平相关变换（U）决定，ti时背景误差 εbi由切线模式M0,i的动力学演变和背景误差水平相关变换（U）特征共同决定。Järvinen（2001）采用观测与背景之差的新息向量（di），利用拟合方法估计背景误差方差与水平相关时发现，在6 h 同化时间窗内，随着预报时效延长背景误差水平相关尺度呈现出逐渐增大的特征，其原因与背景误差随切线模式M0,i的时间积分出现误差能量跨尺度的传播，从小尺度向更大尺度波传播有关。Wee 等（2004）利用观测系统模拟试验（OSSE）对区域MM5 的4DVar 同化系统的背景与分析误差功率谱进行诊断时，也发现从t0时到 tn时误差功率谱能量从小尺度向更大尺度波传播特征。此外，他们对同化时间窗的分析增量功率谱特征诊断时，发现 tn时的分析增量在短波部分明显小于 t0时，分析增量功率谱也呈现出从小尺度向更大尺度波传播特征，分析增量的功率谱变化特征与背景误差功率谱变化特征非常类似。

以上分析表明，在6 h 同化时间窗，切线模式M0,i的动力演变特征和背景误差水平相关变换（U）的特征共同决定背景误差水平相关随时间的演变。GRAPES 全球4DVar 中背景误差水平相关特征随时间的演变特征还不清楚，需要通过诊断来分析。式（2）与（5）表明，分析增量的功率谱演变与背景误差的功率谱演变的约束机制相同，可以用分析增量的水平相关特征在同化时间窗的变化来间接诊断背景误差水平相关特征的变化。

2.2 分析增量的功率谱分布特征
4DVar 同化时间窗计算不同时刻分析增量pagenumber_ebook=117,pagenumber_book=990在球谐谱空间的功率谱特征。按照球谐谱空间的误差协相关表达，球面任意两点相对于球心夹角为θ的误差协相关可以定义为

pagenumber_ebook=117,pagenumber_book=990
式中，m为 纬向波数，n为 二维总波数，而pagenumber_ebook=117,pagenumber_book=990是球谐谱空间的正交基底，pagenumber_ebook=117,pagenumber_book=990为二维总波数为 n的方差，pagenumber_ebook=118,pagenumber_book=991为纬向平均的球谐谱函数，相当于0 阶勒让德函数pagenumber_ebook=118,pagenumber_book=991误差方差可以导出为

pagenumber_ebook=118,pagenumber_book=991
GRAPES 全球4DVar 的分析增量定义在1.0°分辨率的经纬度网格上，利用误差方差公式（8）来计算分析增量功率谱。图1 给出2018 年7 月9 日GRAPES 全球4DVar 风场分析增量与无量纲气压分析增量在6 h 同化时间窗的功率谱变化情况。由图1 可见，对同化时间窗的不同时刻，分析增量的功率谱差异较大。在起始时刻分析增量在40 波以上的短波部分明显偏小，分析增量随波数呈现出较大的斜率特征，而随着时间演变短波部分分析增量逐步增大，斜率有所减小。相对于起始时刻，在终止时刻分析增量的短波部分明显增大，且可以发现在同化时间窗的后期，分析增量在短波部分变化减小，趋于稳定，斜率变化不大。对30 波以下的天气尺度到长波部分，分析增量功率谱呈现出减小的趋势。对无量纲气压同样可以发现在起始时分析增量的40 波以上的短波部分明显偏小，而在终止时刻分析增量在短波部分明显增大。这种特征与分析变量无关，而与同化时间窗的不同时刻有关。前文分析指出，起始时刻背景误差（εb）的谱响应特征由背景误差水平相关变换（U）决定，而任意时刻的背景误差pagenumber_ebook=118,pagenumber_book=991由切线性模式（M0,i）的动力演变和背景误差水平相关变换（U）共同决定，且误差水平相关特征应呈现出逐渐增大的特征。龚建东（2007）研究表明，次天气尺度到中尺度波的误差功率谱增大，表明同化时间窗的水平特征长度随时间演变是缩小的，这与Järvinen（2001）、Wee 等（2004）研究结论存在明显不一致，表明背景误差水平相关变换（U）给定不合理。

3 同化时间窗功率谱演变特征诊断
3.1 集合资料同化生成扰动样本
pagenumber_ebook=118,pagenumber_book=991
图1 风场（a，单位：m2/s2）与无量纲气压（b）分析增量功率谱在6 h 同化时间窗的演变（纵轴为功率谱方差，采用常用对数坐标，T0至T6分别对应同化时间窗起始（t0）至终止（tn）时刻）
Fig.1 Variance spectra of analysis increment of vector wind（a，unit：m2/s2）and non-dimensional pressure（b）evolved in the 6 h data assimilation time window（T0to T6denote times from the beginning（t0）to the end（tn）of the assimilation window）

背景误差水平相关特征需通过新息向量计算（Järvinen，2001）或预报样本来估计。传统上多采用在同一检验时刻不同预报时效的预报差异样本来估计6 h 背景误差协方差。Fisher（2003）研究表明，这种方法估计背景误差水平相关时，随着预报差异的时长增加，会出现水平相关特征长度明显高估的情况，特别是对1500 km 以内的水平相关高估更为明显。一种可行方法是产生6 h 预报的集合样本来估计背景误差。为准确估计6 h 背景误差水平相关特征，文中使用扰动观测和扰动海温的方法，建立集合资料同化系统（Ensemble of Data Assimilation，EDA）来生成统计样本。因所有样本都是6 h的预报，与实际背景场为6 h 的预报时效一致，模式的偏差特征也一致，可以比较容易去除系统性偏差，从而准确估计背景误差的水平相关特征（Fisher，2003）。Isaksen 等（2010）研究表明，基于扰动观测的集合资料同化方法与集合卡曼滤波方法等价，该方法生成的统计样本其分析与预报误差协方差与集合卡曼滤波的分析与预报误差协方差的公式表达一致。

扰动观测和海温分析是在已有GRAPES 全球4DVar 系统基础上通过扰动观测资料和扰动海温分析来表示观测及海洋下边界的不确定性对分析和预报的影响。利用不同扰动观测产生多组4DVar分析，在模式积分时扰动海温分析来考虑模式不确定性获得预报样本。扰动观测和扰动海温分析可以表示为分析步骤和预报步骤

pagenumber_ebook=119,pagenumber_book=992
pagenumber_ebook=119,pagenumber_book=992
3.2 集合资料同化在GRAPES 中的具体实现
GRAPES 全球4DVar 所使用的观测资料包括探空、小球测风、地面、船舶、飞机、云导风等常规资料，以及欧美微波温度垂直探测资料AMSU-A、红外高光谱垂直探测资料IASI 与AIRS、GPS 掩星折射率资料等（龚建东等，2016；Zhang，et al，2018）。对每种资料的每个观测值按照所给定观测误差叠加高斯型随机扰动，并假定观测之间误差不相关，也即不同观测之间的随机扰动不存在依赖关系。式（9）在具体实现时，选择质量控制后的观测资料，按照集合资料同化的样本数量，对每个观测叠加随机扰动值，生成观测数量相同、观测数值不同的多组观测，每个集合资料同化样本使用其中的一组观测资料。由于相对湿度观测有非负和不能超饱和限制，在叠加高斯型随机扰动时，当相对湿度扰动后的数值小于5%时设为5%，而当出现超饱和时设最大观测为100%。为节省计算量，EDA 使用低分辨率全球四维变分同化系统，其外循环分辨率与GRAPES 全球4DVar 的内循环分辨率一致。EDA外循环分辨率为1.0°，内循环分辨率为1.5°，垂直60 层。EDA 为逐6 h 循环分析系统，每次循环生成40 个集合预报样本，统计时总计使用了240 个样本。

在集合资料同化的分析预报循环中，考虑海温分析的不确定性对集合发散度的影响。GRAPES全球预报使用美国环境预报中心（NCEP）的每日海温分析（Reynolds，et al，2007；Maturi，et al，2017），其分析误差为0.1—0.7 K，其中在赤道东太平洋和赤道西太平洋、西北太平洋、赤道大西洋、大西洋西北部、印度洋西部等区域，海温分析误差为0.4—0.6 K，其他区域分析误差较小。Reynolds 等（2007）研究表明海温分析误差的水平相关尺度为100—200 km，考虑到海表温度（SST）的分析误差具有一定的空间分布特征和时间相关特征，对海温的扰动不是按照每个格点独立叠加随机扰动，而是计算不同日期海温分析的差异。本研究中使用2015 年以来的逐日NCEP SST 分析数据，计算每5 d 间隔的海温分析差异，通过质量控制剔除差异较大的数据，并将最大海温扰动幅度控制在0.9 K 以下，从而构成海温分析扰动样本数据集，按照不同月份归类数据集，按照序号排列每个扰动样本。式（10）在具体实现时，选择随机正整数生成算法，按照集合资料同化样本数量的一半生成样本序号子集，并扩展为数值符号相反的另一半子集，从而实现海温扰动随机抽取、正负成对叠加扰动。每个集合资料同化样本使用其中的一个样本序号对应的海温值。

3.3 误差水平相关的谱特征分析
由误差方差公式（7）与（8），水平误差水平相关可以表示为

pagenumber_ebook=120,pagenumber_book=993
背景误差使用集合资料同化获得的240 个6 h 预报样本pagenumber_ebook=120,pagenumber_book=993计算样本集合平均，并获得扰动样本pagenumber_ebook=120,pagenumber_book=993此外，背景误差水平相关变换（U）所给出的二阶自回归相关函数也可以在谱空间展开获得误差功率谱系数。这样可以比较背景误差水平相关变换（U）及背景误差pagenumber_ebook=120,pagenumber_book=993的功率谱的差异。

值得注意的是，由于在集合资料同化中对观测和海温分析叠加了扰动来表述观测和海温分析的不确定性，集合资料同化获得的样本pagenumber_ebook=120,pagenumber_book=993在短波部分随机噪声较大，具有较大误差。Raynaud 等（2009）、Bonavita 等（2011）在使用这些集合样本前进行了低通滤波，抑制短波部分噪音的影响。文中在使用集合样本时，采用基于三角函数的低通滤波器

pagenumber_ebook=120,pagenumber_book=993
式中，Ntrunc是最大谱截断波数，与不同变量和模式层次有关。集合样本的分辨率为1.0°，对应匹配的高斯网格分辨率为1.125°，截断波数可以达到106 个，本研究中最大谱截断波数也取为106，滤除随机噪声的同时，保留天气到中尺度主要扰动波数的扰动振幅。对滤波后的方差进行调整，使得ρ(0)=1。以500 hPa 高度的分析变量流函数、势函数、无量纲气压和比湿为例，其水平相关的功率谱的分布特征如图2 所示。可见，以滤波后的集合资料同化扰动样本pagenumber_ebook=120,pagenumber_book=993（简称EDAFLT）为参考，SOAR 在谱空间功率谱的主要特征是在0—10 波方差比较大，而在20 波以上的天气尺度波和短波部分方差明显较小。SOAR 代表了起始时刻背景误差水平相关功率谱特征，与代表6 h 背景误差水平相关功率谱特征的EDAFLT 相比，存在20 波以上的天气尺度波和短波部分方差低估问题。谱空间背景误差方差低估，将直接造成在对应尺度的分析增量减小，造成分析场不是最优估计。

为比较SOAR 定义中不同水平相关特征长度对功率谱的影响，分别将特征长度参数缩短至原值的0.75 倍、0.5 倍。以特征长度缩小至0.5 倍（简称为SOAR05）为例，结果如图2 所示，在30 波以上的短波部分功率谱明显增大，与SOAR 相比更接近EDAFLT 的结果，但15 波以下的长波部分功率谱显著较小，对流函数和势函数在0—3 波出现功率谱数值明显小于EDAFLT 的结果。GRAPES 全球4DVar 数值试验表明，0—3 波的合理给定对1—10 d预报效果非常关键，这与减缓模式背景场存在大尺度系统性偏差有关。因而需要发展合理给定的起始时刻背景误差水平相关功率谱模型，使得同化时间窗起始与终止时刻的背景误差水平相关特征一致。

4 谱空间误差融合对分析的影响
4.1 谱空间误差融合
上述分析指出，现有GRAPES 全球4DVar 的6 h同化时间窗的起始与终止时刻的背景误差水平相关特征不一致，一种优化的办法是将集合资料同化获得的样本pagenumber_ebook=120,pagenumber_book=993所计算的水平相关的功率谱的方差与背景误差水平相关变换（U）所给出的二阶自回归相关函数的功率谱的方差进行线性加权，融合二阶自回归相关函数在20 波以下方差特征和集合资料同化样本在20 波以上方差特征，使得起始时刻功率谱的方差分布更为合理，但对水汽场不做融合。在两者融合过程中，重点考虑在对流层集合资料同化计算的水平相关的功率谱，在实际实现时，比较了不同权重系数对分析预报的影响，最后采用在对流层10%的权重来自集合资料同化样本计算的功率谱，而在平流层这一比例降低到0.2%。对流层和平流层采用双曲正切曲线光滑过渡。融合后的谱空间方差进行调整，确保 ρ（0）=1。图2 中Merge为误差融合后的起始时刻水平相关功率谱的特征。对流函数、势函数以及无量纲气压，比较SOAR与Merge，可以发现在0—20 波，融合后的结果与SOAR 的功率谱分布特征一致，而在20 波以上功率谱数值要明显大于SOAR，但要略低于EDAFLT 功率谱。比较SOAR05 与Merge 可以发现，在40 波以上的短波部分两者接近。

4.2 误差融合对分析场的影响
pagenumber_ebook=121,pagenumber_book=994
图2 流函数（a）、势函数（b）、无量纲气压（c）与比湿（d）在500 hPa 高度的水平相关的功率谱方差特征（SOAR 为二阶自回归相关函数，SOAR05 为水平特征长度缩小50%，EDA 为集合资料同化样本估计值，EDAFLT 为EDA 滤波后的估计值，Merge 为融合后的数值）
Fig.2 Variance spectra of horizontal correlation for stream-function（a），velocity potential function（b），non-dimensional pressure（c）and specific humidity（d）at 500 hPa geopotential height（SOAR for second-order auto-regression function，SOAR05 for de-correlation length scale reducing 50%，EDA for ensemble member estimation，EDAFLT for EDA with spectral filter，Merge for blending variance spectra）

试验选取2018 年7 月9 日的资料来验证融合水平相关对分析增量功率谱的影响，结果如图3 所示。比较同化时间窗起始时刻与终止时刻，采用SOAR 水平相关模型的分析增量功率谱，对风场在起始时刻30 波以下的功率谱要明显大于终止时刻，而在50 波以上的功率谱，终止时刻要明显大于起始时刻，在同化时间窗内则存在明显的不一致。无量纲气压的分析增量，同样在40 波以上的功率谱终止时刻要明显大于起始时刻。采用融合水平相关分布特征后，风场在起始与终止时刻的分析增量在谱空间表现更为一致，特别是在0—80 波非常接近，无量纲气压0—80 波也是明显好于SOAR 水平相关模型。对比图2 可知，融合后的背景误差水平相关功率谱在天气尺度到短波部分数值增大，这会造成分析增量在该尺度的数值增大。

pagenumber_ebook=122,pagenumber_book=995
图3 风场（a，单位：m2/s2）与无量纲气压（b）的分析增量功率谱方差在6 h 同化时间窗起始时刻（t0）与终止时刻（tn）的变化
Fig.3 Variance spectra of analysis increment for vector wind（a，unit：m2/s2）and non-dimensional pressure（b）evolved at the beginning（t0）and ending（tn）of the 6 h time window

以上分析表明，起始时刻背景误差水平相关变换（U）的特征极为关键。但在目标函数极小化迭代计算中因为起始时刻背景误差水平相关变换（U）的结构与动力模式生成的特征不一致，由式（4）可知，目标函数对控制变量（v）的梯度受到水平相关变换转置矩阵（UT）的控制，在起始时刻由背景误差水平相关变换（U）而不是动力模式生成的特征生成分析增量。当背景误差水平相关变换定义不合理时，这一过程会损失分析增量在起始时刻的信息。由图1分析增量的功率谱特征可见，这种信息损失主要影响到6 h 同化时间窗的前3 h。

4.3 单点分析增量在格点空间特征
对于融合后的功率谱，可以分析其在格点空间的实际表现。对球面上的任意两点的水平相关，根据Boer 等（1983）的研究，相关函数仅仅依赖于点之间的距离角，它可以用勒让德级数展开到谱空间。具体做法是给出以北极点为中心的相关函数在全球高斯格点的分布，然后运用谱格变换得到谱空间的相关函数。同样，已知融合后的功率谱，通过谱格逆变换，可以获得以北极点为中心的水平相关函数的分布。由于高斯格点在北极点没有数值，需要对极点进行数值拟合计算，由计算结果（图4）可见，融合后的水平相关与SOAR 相比在距离北极点较近时体现出水平相关尺度更小的特征，而在远离北极点时水平相关系数下降变慢，与SOAR 相关系数的下降特征一致。这表明融合后的水平相关特征呈现出一种相关系数随距离下降由快到慢、相关特征由“窄”变“宽”的“长尾型”分布特征，这种变化特征与多尺度水平相关模型的特点一致（Purser，et al，2003）。

4.4 格点空间的分析增量
以纬向（u）风场和温度场为例考察格点空间分析增量的变化，结果（图5）可见，在同化时间窗的起始时刻，基于SOAR 水平相关模型的分析增量不但最大分析值的数值偏小，而且分析尺度大，分析增量的细节不足，反映出天气尺度到次天气尺度的分析特征不足。而采用融合水平相关模型的分析增量，一方面分析增量数值明显偏大，表明在靠近观测的位置分析场更靠近观测，观测更多地被分析场吸收，另一方面天气到次天气尺度的分析特征明显增强，分析呈现出更多的细节，这表明在起始时刻融合水平相关特征可以更好地吸纳观测信息。比较温度分析增量，可以发现融合水平相关模型虽然也增加了分析增量的细节，但与SOAR 水平相关模型的差异并没有风场显著。其原因在于风场计算是通过对流函数和势函数的水平求导求解的，由图4可见，融合水平相关模型的“长尾型”的分布特征与SOAR 的分布特征存在不同，两者的水平导数也存在差异（图略），造成两个模型的分析增量特征不同。温度分析增量是通过无量纲气压分析增量计算得到的，并不涉及水平求导计算，两者的分析增量特征较为接近。

pagenumber_ebook=123,pagenumber_book=996
图4 流函数（a）、势函数（b）、无量纲气压（c）与比湿（d）的水平相关系数分布
Fig.4 Horizontal correlation of stream-function（a），velocity potential-function（b），no-dimensional pressure（c）and specific humidity（d）

进一步分析目标函数随迭代步数下降的收敛情况（图6）可见，采用融合水平相关模型后，目标函数随迭代步数的下降明显改善，比较目标函数中的背景项（Jb）与观测项（Jo），可以发现主要变化来自观测项（Jo），而背景项（Jb）的变化较小，这表明采用融合水平相关后初始场的积分轨迹与观测的拟合更好。目标函数对分析的梯度，两种方法的区别不大。在引言中提到起始和终止时刻的背景或分析误差协方差的总体统计特征应当相近才能获得更好的分析结果，本研究表明，合理给定背景误差协方差，使得起始和终止时刻背景误差协方差的统计特征相近，将改善积分轨迹在同化时间窗与观测拟合，使得分析结果更合理。积分轨迹在同化时间窗与观测拟合以及起始和终止时刻背景误差协方差的统计特征相近，是4DVar 的两个重要方面，这一点在以往的研究工作中强调不够。

pagenumber_ebook=125,pagenumber_book=998
图5 起始时刻（t0）水平风场（a、b.单位：m/s）与温度（c、d.单位：K）的分析增量（a、c.SOAR 水平相关模型，b、d.融合水平相关）
Fig.5 Horizontal wind（a，b.unit：m/s）and temperature（c，d.unit：K）analysis increment at the beginning time t0（a，c.analysis with SOAR；b，d.analysis with blending variance spectra）

pagenumber_ebook=125,pagenumber_book=998
图6 归一化后的目标函数随迭代步数下降的收敛情况
Fig.6 Converge of normalized cost-function with iteration steps of minimization

5 融合水平相关对分析预报的影响
5.1 对分析与背景误差的影响
前面的分析指出，采用融合水平相关可以更好地吸纳观测信息，这一点可以从分析误差与6 h 预报的背景误差表现出来。分别采用SOAR 水平相关模型（对照试验）和融合水平相关模型（对比试验）进行2016 年6 月的比较试验。以欧洲中期天气预报中心（ECMWF）的再分析场为大气“真值”，比较4DVar 同化后分析误差与背景误差的变化，结果（图7）表明，对温度分析采用融合水平相关特征后，在700—150 hPa 的对流层分析误差和背景误差明显减小，在500 hPa 上下对比试验的背景误差与对照试验分析误差相当。同样对风场分析，在700—150 hPa 的对流层分析误差和背景误差也明显减小，在300 hPa 改进尤为明显，分析误差和背景误差显著改进。这表明引入融合水平相关模型，观测信息被更有效地吸收，分析误差和背景误差显著减小。

5.2 对预报的影响
进一步考察采用SOAR 水平相关模型和融合水平相关模型对预报的影响。对2016 年6 月每日12 时进行10 d 预报（图8），以南、北半球500 hPa位势高度、温度与风场为检验对象，结果表明融合相关模型的引入显著改善了北半球预报的距平相关系数，分析对预报的影响改善在前120 h 通过95%的置信度检验，对均方根误差的改进更为明显，对前192 h 的预报都有显著影响（图略）。对南半球也有改善作用，对距平相关系数改进有影响但不显著，对均方根误差的影响超过2 d。对北半球温度场也有改善，对前192 h 的预报有显著影响，南半球影响超过96 h（图略）。对风场的影响与温度场非常类似，对北半球影响更为显著，南半球相对弱一点。对赤道地区的影响总体比较弱（图略）。

6 总结与讨论
pagenumber_ebook=125,pagenumber_book=998
图7 温度（a）与东西风场（b）的6 h 背景场与分析相对于ECMWF 再分析的均方根误差垂直分布特征（实线为分析误差，虚线为背景误差；黑线为采用SOAR 模型的结果，红线为采用融合相关模型的结果）
Fig.7 Vertical distributions of 6 h background and analysis RMSEs of temperature（a）and zonal wind（b）against ECMWF ERA-Interim reanalysis（The solid-line indicates analysis error，the dashed-line is for background error，the black line denotes results from SOAR，and the red line is for results from blending variance spectra）

pagenumber_ebook=126,pagenumber_book=999
图8 500 hPa 10 d 预报效果（a、b.北半球与南半球距平相关系数，c、d.北半球温度（单位:K）与风场（单位：m/s）的均方根误差；4DVar1 为对照试验，4DVar2 为对比试验；两个方案评分差异落在红色方框外代表通过95%的信度检验）
Fig.8 500 hPa 10 d forecast（a，b.ACC of north-hemisphere and south-hemisphere；c，d.RMSE of temperature（unit：K）and vector wind（unit:m/s）of north-hemisphere；4Dvar1 for control experiment，4DVar2 for comparative experiment；differences outside the hallow bars are significant at the 95% confidence level）

相对于3DVar 不考虑观测与背景场时间对应的不足，4DVar 资料同化方法的一个显著特征是考虑了每个观测与背景场在时间上的多时刻对应关系，使得新息向量的计算更为准确。GRAPES 全球4DVar 发展过程中，对4DVar 分析效果的评估也重点放在这个方面，以目标函数收敛到极小值、分析误差减小作为度量4DVar 分析质量改善的重要依据。文中再次检视了4DVar 的式（4）—（6）定义，结合Lorenc 等（2005）研究所指出的4DVar 在同化时间窗可以隐式描述背景误差协方差随时间演变的特征，而背景误差水平相关特征是背景误差协方差的一个重要方面，本研究以此为切入点，重点关注起始和终止时刻背景误差水平相关功率谱的特征是否合理。

通过对GRAPES 全球4DVar 的分析增量在同化时间窗的不同时刻的诊断，发现对同化时间窗的不同时刻分析增量的功率谱差异很大，在起始时刻分析增量在短波部分明显偏小，而随着时间演变短波部分分析增量逐步增大，在终止时刻分析增量的短波部分明显增大，这种特征与目标函数极小化迭代计算步数和分析变量无关，而与同化时间窗的不同时刻有关，表明背景误差水平相关模型定义存在不合理之处。为准确估计6 h 背景误差水平相关特征，文中使用扰动观测和扰动海温的方法建立集合资料同化系统（EDA）来生成统计6 h 集合样本，因所有样本都是6 h 的预报，模式偏差特征一致，可以容易地去除系统性偏差，从而准确估计背景误差的水平相关特征。研究发现SOAR 在谱空间功率谱的主要特征是在0—10 波方差比较大，而在20 波以上的天气尺度波和短波部分方差明显偏小，存在低估情况。谱空间背景误差方差低估将直接造成在对应尺度的分析增量减小。针对该不足，文中提出将集合资料同化获得的样本所计算的水平相关的功率谱的方差与背景误差水平相关变换所给出的二阶自回归相关模型的功率表的方差进行融合，使得谱空间方差分布更加合理。在两者融合过程中，重点考虑在对流层集合资料同化计算的水平相关的功率谱，对流层10%的权重来自集合资料同化样本计算的功率谱，而在平流层这一比例降低到0.2%。融合后的功率谱显著增加了20 波以上的天气尺度及短波部分，而在0—20 波的功率谱基本保持不变。采用融合水平相关分布特征后，起始与终止时刻的分析增量在谱空间表现更为一致，特别是在0—80 波非常接近，明显好于SOAR 水平相关模型。融合后的水平相关特征在格点空间呈现出一种相关系数随距离下降由快到慢、相关特征由“窄”变“宽”的“长尾型”分布特征，这种变化特征与多尺度水平相关模型的特点一致。格点空间的分析增量也表明，采用融合水平相关的分析增量，一方面分析增量数值明显偏大，表现在靠近观测的位置分析场更靠近观测，观测信息更多地被分析场吸收，另一方面天气到次天气尺度的分析特征明显增强，分析呈现出更多的细节结构，这表明在起始时刻融合水平相关特征可以更好地吸纳观测信息。积分轨迹在同化时间窗与观测拟合以及起始和终止时刻背景误差协方差的统计特征相近，是4DVar 的两个重要方面，这一点在以往的研究工作中强调不够。融合水平相关对分析预报的影响表明，引入融合水平相关模型，观测信息被更有效地吸收，分析误差和背景误差显著减小，显著改善北半球预报的距平相关系数，对均方根误差的改进更为明显。对南半球也有改善作用，但不显著。北半球温度场、风场影响更为显著，南半球相对弱一点，对赤道地区的影响总体比较弱。

Järvinen（2001）研究表明，在观测空间的6 h 同化时间窗随着预报时效延长背景误差水平相关特征呈现出逐渐增大的特征，其原因与背景误差随切线模式的时间积分出现误差能量跨尺度的传播，从小尺度向更大尺度波传播有关。本研究通过融合水平相关模型改善了原有相关模型在50 波以上的功率谱在终止时刻要明显大于起始时刻，在同化时间窗内则存在明显不一致的不足，但没有发现从小尺度向更大尺度波传播的特征，这一点在未来的研究中还需要进一步诊断误差能量的传播特征来更合理地表述背景误差协方差的演变。

参考文献

龚建东.2007.资料同化中二维特征长度随模式分辨率变化的分析研究.大气科学，31（3）：459-467.Gong J D.2007.The analysis on variation of horizontal de-correlation length with model resolution in data assimilation system.Chinese J Atmos Sci，31（3）：459-467（in Chinese）

龚建东，王瑞春，郝民.2016.温湿统计平衡约束关系对GRAPES 全球湿度分析的作用.气象学报，74（3）：380-396.Gong J D，Wang R C，Hao M.2016.The impact of a balance constraint between temperature and humidity on the global humidity analysis in GRAPES.Acta Meteor Sinica，74（3）：380-396（in Chinese）

龚建东，刘永柱，张林.2019.面向四维变分资料同化的NSAS 积云深对流参数化方案的简化及线性化研究.气象学报，77（4）：595-616.Gong J D，Liu Y Z，Zhang L.2019.A study of simplification and linearization of the NSAS deep convection cumulus parameterization scheme for 4D-Var. Acta Meteor Sinica，77（4）：595-616（in Chinese）

刘艳，薛纪善.2019.GRAPES 的新初始化方案.气象学报，77（2）：165-179.Liu Y，Xue J S.2019.The new initialization scheme of the GRAPES. Acta Meteor Sinica，77（2）：165-179（in Chinese）

刘永柱，张林，金之雁.2017.GRAPES 全球切线性和伴随模式的调优.应用气象学报，28（1）：62-71.Liu Y Z，Zhang L，Jin Z Y.2017.The optimization of GRAPES global tangent linear model and adjoint model. J Appl Meteor Sci，28（1）：62-71（in Chinese）

王金成，庄照荣，韩威等.2014.GRAPES 全球变分同化背景误差协方差的改进及对分析预报的影响：背景误差协方差三维结构的估计.气象学报，72（1）：62-78.Wang J C，Zhuang Z R，Han W，et al.2014.An improvement of background error covariance in the global GRAPES variational data assimilation and its impact on the analysis and prediction：Statistics of the three-dimensional structure of background error covariance.Acta Meteor Sinica，72（1）：62-78（in Chinese）

王瑞春，龚建东，张林等.2015.热带风压场平衡特征及其对GRAPES 系统中同化预报的影响研究 II：动力与统计混合平衡约束方案的应用.大气科学，39（6）：1225-1236.Wang R C，Gong J D，Zhang L，et al.2015.Tropical balance characteristics between mass and wind fields and their impact on analyses and forecasts in GRAPES system.Part II：Application of linear balance equation-regression hybrid constraint scheme.Chinese J Atmos Sci，39（6）：1225-1236（in Chinese）

薛纪善.2009.气象卫星资料同化的科学问题与前景.气象学报，67（6）：903-911.Xue J S.2009.Scientific issues and perspective of assimilation of meteorological satellite data.Acta Meteor Sinica，67（6）：903-911（inChinese）

张华，薛纪善，庄世宇等.2004.GRAPES 三维变分同化系统的理想试验.气象学报，62（1）：31-41.Zhang H，Xue J S，Zhuang S X，et al.2004.Idea experiments of GRAPES three-dimensional variational data assimilation system.Acta Meteor Sinica，62（1）：31-41（in Chinese）

张林，刘永柱.2017.GRAPES 全球四维变分同化系统极小化算法预调节.应用气象学报，28（2）：168-176.Zhang L，Liu Y Z.2017.The preconditioning of minimization algorithm in GRAPES global fourdimensional variational data assimilation system.J Appl Meteor Sci，28（2）：168-176（in Chinese）

Bannister R N.2008a.A review of forecast error covariance statistics in atmospheric variational data assimilation.I：Characteristics and measurements of forecast error covariances.Quart J Roy Meteor Soc，134（637）： 1951-1970

Bannister R N.2008b.A review of forecast error covariance statistics in atmospheric variational data assimilation.II：Modelling the forecast error covariance statistics.Quart J Roy Meteor Soc，134（637）：1971-1996

Boer G J，Shepherd T G.1983.Large-scale two-dimensional turbulence in the atmosphere.J Atmos Sci，40（1）：164-184

Bonavita M，Raynaud L，Isaksen L.2011.Estimating background-error variances with the ECMWF Ensemble of Data Assimilations system：Some effects of ensemble size and day-to-day variability.Quart J Roy Meteor Soc，137（655）：423-434

Daley R.1991.Atmospheric Data Analysis.Cambridge： Cambridge University Press，457pp

Derber J，Bouttier F.1999.A reformulation of the background error covariance in the ECMWF global data assimilation system.Tellus A，51（2）：195-221

Fisher M.2003.Background error covariance modelling∥Proceedings of the ECMWF Seminar on Recent Developments in Data Assimilation for Atmosphere and Ocean.Shinfield Park，Reading：ECMWF

Gauthier P，Thépaut J N.2001.Impact of the digital filter as a weak constraint in the preoperational 4DVar assimilation system of Météo-France.Mon Wea Rev，129（8）：2089-2102

Gauthier P，Tanguay M，Laroche S，et al.2007.Extension of 3DVar to 4DVar： Implementation of 4DVar at the meteorological service of Canada.Mon Wea Rev，135（6）：2339-2354

Isaksen L，Bonavita M，Buizza R，et al.2010.Ensemble of data assimilations at ECMWF.ECMWF Technical Memorandum 636

Järvinen H.2001.Temporal evolution of innovation and residual statistics in the ECMWF variational data assimilation systems.Tellus A，53（3）：333- 347

Lopez P，Moreau E.2005.A convection scheme for data assimilation： Description and initial tests.Quart J Roy Meteor Soc，131（606）：409-436

Lorenc A C，Rawlins F.2005.Why does 4D-Var beat 3D-Var?.Quart J Roy Meteor Soc，131（613）：3247-3257

Maturi E，Harris A，Mittaz J，et al.2017.A new high-resolution sea surface temperature blended analysis.Bull Amer Meteor Soc，98（5）：1015-1026

Purser R J，Wu W S，Parrish D F，et al.2003.Numerical aspects of the application of recursive filters to variational statistical analysis.Part II：Spatially inhomogeneous and anisotropic general covariances.Mon Wea Rev，131（8）：1536-1548

Rabier F，Järvinen H，Klinker E，et al.2000.The ECMWF operational implementation of four-dimensional variational assimilation.I：Experimental results with simplified physics.Quart J Roy Meteor Soc，126（564）：1143-1170

Rawlins F，Ballard S P，Bovis K J，et al.2007.The Met Office global fourdimensional variational data assimilation scheme.Quart J Roy Meteor Soc，133（623）：347-362

Raynaud L，Berre L，Desroziers G.2009.Objective filtering of ensemblebased background-error variances.Quart J Roy Meteor Soc，135（642）： 1177-1199

Reynolds R W，Smith T M，Liu C Y，et al.2007.Daily high-resolutionblended analyses for sea surface temperature.J Climate，20（22）：5473- 5496

Talagrand O，Courtier P.1987.Variational assimilation of meteorological observations with the adjoint vorticity equation.I：Theory.Quart J Roy Meteor Soc，113（478）：1311-1328

Wee T K，Kuo Y H.2004.Impact of a digital filter as a weak constraint in MM5 4DVar：An observing system simulation experiment.Mon Wea Rev，132（2）：543-559

Zhang L，Liu Y Z，Liu Y，et al.2019.The operational global four-dimensional variational data assimilation system at the China Meteorological dministration.Quart J Roy Meteor Soc，145（722）：1882-1896

Zhang L H，Gong J D，Wang R C.2018.Diagnostic analysis of various observation impacts in the 3DVar assimilation system of global GRAPES.Mon Wea Rev，146（10）：3125-3142

GONG Jiandong ZHANG Lin WANG Jincheng

National Meteorological Center/Numerical Weather Prediction Center，Beijing100081，China

Abstract Temporal evolution characteristics of the GRAPES global 4DVar analysis increment in the spectral space is investigated,and the impact on the analysis and forecast is also analyzed when the background error horizontal correlation characteristics at the beginning and end of the data assimilation time window are obviously different.Three types of horizontal correlation model,i.e.,the second-order autoregressive model (SOAR),the statistical model from the samples generated by ensemble of data assimilation(EDA),and the blending results of SOAR and EDA,are compared.The results show that the information of synoptic scale in the analysis increment is obviously underestimated by the SOAR model.For the horizontal correlation power spectrum calculated by blending of the SOAR and the EDA,the results show a multi-scale horizontal correlation characteristics,which can better absorb observation information and significantly improve the forecast of geopotential height and temperature in the northern hemisphere.The wind forecast is also improved in the southern hemisphere,while a neutral impact is found in the tropical region.The above results indicate that the merging horizontal correlation scheme developed in this paper is reasonable and practical.

Key words GRAPES 4DVar，Background error horizontal correlation，Variance spectra blending，Numerical experiment

2020-02-03 收稿，2020-07-06 改回.

* 资助课题： 国家重点研发计划项目（2018YFC1506700、2018YFC1506705）。

作者简介： 龚建东，主要从事数值预报资料同化、集合预报研究。E-mail：gongjd@cma.gov.cn

龚建东，张林，王金成.2020.背景误差水平相关结构对四维变分资料同化的影响研究.气象学报，78（6）：988-1001

Gong Jiandong，Zhang Lin，Wang Jincheng. 2020. An impact study of background error horizontal correlation structure on 4DVar. Acta Meteorologica Sinica ，78（6）:988-1001

中图法分类号 P435，P456.7