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）

$\begin{split} J({\text δ}{{X}}^{\mathrm{a}})=&\frac{\,1\,}{\,2\,}{{\text δ}{{X}}^{\mathrm{aT}}}{{B}}_{0}^{-1}{\text δ}{{X}}^{\mathrm{a}}+\frac{\,1\,}{\,2\,}\sum\limits_{{i}=0}^{{n}}\Big[{({{d}}_{{i}}\!-\!{{H}}_{{i}}{{M}}_{0,{i}}{\text δ}{{X}}^{\mathrm{a}})}^{\mathrm{T}}{\text •}\\&{{R}}^{-1}({{d}}_{{i}}\!-\!{{H}}_{{i}}{{M}}_{0,{i}}{\text δ}{{X}}^{\mathrm{a}})\Big]+\!{J}_{\mathrm{c}} \end{split}$

(1)

式（1）右端第1项为背景项（ $ {J}_{\mathrm{b}} $ ），第2项为观测项（ $ {J}_{\mathrm{o}} $ ）， $ {J}_{\mathrm{c}} $ 是基于数字滤波的弱约束项。式中 ${i} $ 为时间 $ {t}_{{i}} $ 的下标， $ {t}_{0} $ 与 $ {t}_{{n}} $ 为6 h同化时间窗的起始与终止时刻，取为00、06、12、18时（世界时，下同）的前后各3 h。 $ {\text δ}{{X}}^{\mathrm{a}}={{X}}_{0}^{\mathrm{a}}-{X}_{0}^{\mathrm{b}} $ 为起始时（ $ {t}_{0} $ ）的分析增量， $ {{X}}_{0}^{\mathrm{b}} $ 为背景场， $ {{X}}_{0}^{\mathrm{a}} $ 为分析场。 $ {{B}}_{0} $ 是背景误差协方差，由统计模型给出（王金成等，2014）。 $ {{M}}_{0,{i}} $ 是非线性预报模式 $ \mathcal{M} $ 的切线性模式，表示 $ {t}_{0} $ 到 $ {t}_{{i}} $ 的时间积分，H是非线性观测算子 $ \mathcal{H} $ 的切线算子，R是观测误差协方差。新息向量 ${{d}}_{{i}}={{y}}_{{i}}^{\mathrm{o}}-{\mathcal{H}}_{{i}}{\mathcal{M}}_{0,{i}}{({X}}_{0}^{\mathrm{b}})$ ，其中 $ {{y}}_{{i}}^{\mathrm{o}} $ 是 $ {t}_{{i}} $ 时的观测值，由于4DVar是时间上的连续分析， $ {t}_{{i}} $ 时刻背景场与分析场分别表示为 ${{X}}_{{i}}^{\mathrm{b}}={\mathcal{M}}_{0,{i}}{({X}}_{0}^{\mathrm{b}})$ ，以及 $ {{X}}_{i}^{\mathrm{a}}= {{X}}_{{i}}^{\mathrm{b}}+{{M}}_{0,{i}}{\text δ}{{X}}^{\mathrm{a}} $ 。

4DVar中采用分析变量的变量变换方式来避免显式求解背景误差协方差矩阵（ $ {{B}}_{0} $ ）的逆矩阵（张华等，2004；王瑞春等，2015）。利用根号矩阵表达 $ {{B}}_{0}={{B}}_{0}^{1/2}{{B}}_{0}^{\mathrm{T}/2} $ ，以及变量变换 $ {\text δ}{{X}}^{\mathrm{a}}={{B}}_{0}^{1/2}{v} $ ， $ {v} $ 为控制变量，是归一化的无量纲向量，在GRAPES全球4DVar中， ${v}=({\psi },{{\chi }}_{\mathrm{u}},{{\pi }}_{\mathrm{u}},{q})$ ，分别为流函数、非平衡势函数、非平衡无量纲气压以及比湿场。根号矩阵（ $ {{B}}_{0}^{1/2} $ ）用误差模型近似表达为 $ {{B}}_{0}^{1/2}={P}{{\varSigma }}_{\mathrm{u}}{U} $ ，其中 $ {P} $ 为物理变换，主要实现非平衡物理变量向物理变量的转换以及从流函数与势函数到风场的物理变换， $ {{\varSigma }}_{\mathrm{u}} $ 为非平衡部分的背景误差方差， $ {U}={{U}}_{\mathrm{h}}{{U}}_{\mathrm{v}} $ 分别表示水平与垂直相关变换，为表述简单起见，文中 $ {U} $ 仅表示水平相关变换。于是 $ {t}_{{i}} $ 时分析增量表示为

$ {\text δ}{{X}}_{{i}}^{\mathrm{a}}={{M}}_{0,{i}}{\text δ }{{X}}^{\mathrm{a}}={{M}}_{0,{i}}{P}{{\varSigma }}_{{u}}{U}{v} $

(2)

略去式（1）中的 $ {J}_{\mathrm{c}} $ 项，于是目标函数表述为

$ \begin{split} J({v})=&\frac{\,1\,}{\,2\,}{{v}}^{\mathrm{T}}{v}+\frac{\,1\,}{\,2\,}\sum\limits_{{i}=0}^{{n}}\Big[{({{d}}_{{i}}-{{H}}_{{i}}{{M}}_{0,{i}}{P}{{\varSigma }}_{{u}}{U}{v})}^{\mathrm{T}}{\text •}\\&{{R}}^{-1}({{d}}_{{i}}-{{H}}_{{i}}{{M}}_{0,{i}}{P}{{\varSigma }}_{{u}}{U}{v})\Big] \end{split}$

(3)

目标函数相对于 $ {v} $ 的梯度是

$ {\nabla }_{\mathrm{v}}J={v}+\sum\limits_{{i}=0}^{{n}}{{U}}^{\mathrm{T}}{{\varSigma }}_{{u}}{{P}}^{\mathrm{T}}{{M}}_{{i},0}^{\mathrm{T}}{{{H}}_{{i}}^{\mathrm{T}}{R}}^{-1}({{d}}_{{i}}-{{H}}_{{i}}{{M}}_{0,{i}}{P}{{\varSigma }}_{{u}}{U}{v}) $

(4)

由式（2），同样 $ {t}_{{i}} $ 时背景误差 $ {{\varepsilon }}_{{i}}^{\mathrm{b}} $ 可以表示为

$ {{\varepsilon }}_{{i}}^{\mathrm{b}}={{M}}_{0,{i}}{{\varepsilon }}^{\mathrm{b}}={{M}}_{0,{i}}{P}{{\varSigma }}_{{u}}{U}{{\varepsilon }}^{\mathrm{v}} $

(5)

式中， $ {{\varepsilon }}^{\mathrm{b}} $ 为 $ {t}_{0} $ 时背景误差， $ {{\varepsilon }}^{\mathrm{v}} $ 为 $ {t}_{0} $ 时归一化无量纲向量 $ {v} $ 的随机扰动， $ {E}=\left\langle {{\varepsilon }}^{\mathrm{v}},{{{\varepsilon }}^{\mathrm{v}}}^{\mathrm{T}} \right\rangle $ 是单位矩阵。在4DVar极小化循环迭代过程中，切线与伴随模式积分在同化时间窗隐式演变背景误差协方差，于是表示为

$ \begin{split} {{B}}_{{i}}&=\left\langle {{\varepsilon }}_{{i}}^{\mathrm{b}},{{{\varepsilon }}_{{i}}^{\mathrm{b}}}^{\mathrm{T}} \right\rangle ={{M}}_{0,{i}}{{{B}}_{0}{M}}_{0,{i}}^{\mathrm{T}} \\& {={M}}_{0,{i}}{P}{{\varSigma }}_{{u}}{U}{E}{{U}}^{\mathrm{T}}{{{\varSigma }}_{{u}}^{\mathrm{T}}}{{P}}^{\mathrm{T}}{{M}}_{0,{i}}^{\mathrm{T}} \end{split}$

(6)

由式（6）可见，在 $ {t}_{0} $ 时通过统计模型给出背景误差协方差 $ {{B}}_{0} $ ，而在同化时间窗的其他时刻，背景误差协方差 $ {{B}}_{{i}} $ 由切线模式的动力约束来控制。

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

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

2.2 分析增量的功率谱分布特征

4DVar同化时间窗计算不同时刻分析增量 $ {\text δ}{{X}}_{{i}}^{\mathrm{a}} $ 在球谐谱空间的功率谱特征。按照球谐谱空间的误差协相关表达，球面任意两点相对于球心夹角为 $ \theta $ 的误差协相关可以定义为

$ f\left(\theta \right)=\sum\limits_{{n}}\sum\limits_{{m}=-{n}}^{{n}}{{\delta }_{{n}}^{{m}}}^{2}{Y}_{{n}}^{{m}}=\sum\limits_{{n}}\sqrt{2{n}+1}{\delta }_{{n}}^{2}{Y}_{{n}}^{0} $

(7)

式中， $ {m} $ 为纬向波数， $ {n} $ 为二维总波数，而 $ {Y}_{{n}}^{{m}} $ 是球谐谱空间的正交基底， $ {\delta }_{{n}}^{2} $ 为二维总波数为 $ {n} $ 的方差， $ {Y}_{{n}}^{0} $ 为纬向平均的球谐谱函数，相当于0阶勒让德函数（ $ {P}_{{n}}^{0} $ ），误差方差可以导出为

$f(0) = \mathop \sum \limits_{{n}} \sqrt {2{{n}} + 1} \delta _{{n}}^2P_{{n}}^0(1) = \mathop \sum \limits_{{n}} \delta _{{n}}^2({2{{n}} + 1} )$

(8)

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

3 同化时间窗功率谱演变特征诊断 3.1 集合资料同化生成扰动样本

背景误差水平相关特征需通过新息向量计算（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分析，在模式积分时扰动海温分析来考虑模式不确定性获得预报样本。扰动观测和扰动海温分析可以表示为分析步骤和预报步骤

$ \begin{split} J\left({{\text δ}{{\tilde{ X}}^{\rm{a}}}} \right) =& \frac{\,1\,}{\,2\,}{\text δ}{\tilde{ X}^{\rm{aT}}}{{B}}_0^{ - 1}{ \text δ}{\tilde{ X}^{\rm{a}}} + \frac{\,1\,}{\,2\,}\mathop \sum \limits_{{{i}} = 0}^{{n}} {\Big[\left({{{\tilde{ d}}_{{i}}} - {{H}_{{i}}}{{M}_{0,{{i}}}}{\text δ}{{\tilde{ X}}^{\rm{a}}}} \right)^{\rm{T}}}{\text •}\\&{{R}^{ - 1}}\left({{{\tilde{ d}}_{{i}}} - {{H}_{{i}}}{{M}_{0,{{i}}}}{\text δ}{{\tilde{ X}}^{\rm{a}}}} \right)\Big]\\[-12pt]\end{split}$

(9)

$ {{\tilde{{X}}}}^{\mathrm{b}}={\mathcal{M}}_{0,{n}}\left({\tilde{{X}}}_{0}^{\mathrm{a}}\right)+\tilde{{q}} $

(10)

式中， ${\tilde{{d}}}_{{i}}={{y}}_{{i}}^{\mathrm{o}}+{{\eta }}_{{i}}-{\mathcal{H}}_{{i}}{\mathcal{M}}_{0,{i}}$ （ ${\tilde{{X}}}_{0}^{\mathrm{b}} $ ）为叠加扰动观测后的新息向量，其中 $ {{y}}_{{i}}^{\rm{o}} $ 为 $ {t}_{{i}} $ 时的观测， $ {{\eta }}_{{i}} $ 为叠加的观测误差扰动，服从期望值为0、方差为 $ {R} $ 的观测误差分布， $ {{y}}_{{i}}^{\rm{o}}+{{\eta }}_{{i}} $ 为叠加扰动后的观测。 $ \tilde{{q}} $ 为叠加的海温分析随机扰动，服从期望值为0、方差为海温分析误差的随机扰动。（ $ \tilde{\bullet } $ ）表示对一组样本中的任意一个进行分析或预报，形成分析或预报场。定义δ $ \tilde{{X}}^{\mathrm{b}} = $ ${\tilde{{X}}}^{\mathrm{b}}- {\overline{{\tilde{{X}}}^{\mathrm{b}}}} $ ，为去除样本平均值（ $ \overline{{\tilde{{X}}}^{\mathrm{b}}} $ ）后的预报扰动样本。由式（9）和（10）可见，通过扰动观测和扰动海温分析方式获得一组预报扰动样本，可以用来分析6 h背景误差水平相关的特征。

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），水平误差水平相关可以表示为

$\rho \left(\theta \right) = \frac{{f\left(\theta \right)}}{{f\left(0 \right)}} = \frac{1}{{f\left(0 \right)}}{ \sum \limits_{{n}}}_{}\;{ \sqrt {2{{n}} + 1} \delta _{{n}}^2Y_{{n}}^0}$

(11)

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

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

$ f\left(\mathrm{n}\right)={\mathrm{c}\mathrm{o}\mathrm{s}}^{2}\left(\frac{\ {\text{π} }\ }{2}\mathrm{min}\left({n},{{N}}_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{n}\mathrm{c}}\right)/{{N}}_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{n}\mathrm{c}}\right) $

(12)

式中， $ {{N}}_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{n}\mathrm{c}} $ 是最大谱截断波数，与不同变量和模式层次有关。集合样本的分辨率为1.0°，对应匹配的高斯网格分辨率为1.125°，截断波数可以达到106个，本研究中最大谱截断波数也取为106，滤除随机噪声的同时，保留天气到中尺度主要扰动波数的扰动振幅。对滤波后的方差进行调整，使得 $ \rho \left(0\right)=1 $ 。以500 hPa高度的分析变量流函数、势函数、无量纲气压和比湿为例，其水平相关的功率谱的分布特征如图2所示。可见，以滤波后的集合资料同化扰动样本δ $ {\tilde{{X}}}^{\mathrm{b}} $ （简称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同化时间窗的起始与终止时刻的背景误差水平相关特征不一致，一种优化的办法是将集合资料同化获得的样本（ $ {\tilde{{X}}}_{0}^{\mathrm{b}} $ ）所计算的水平相关的功率谱的方差与背景误差水平相关变换（U）所给出的二阶自回归相关函数的功率谱的方差进行线性加权，融合二阶自回归相关函数在20波以下方差特征和集合资料同化样本在20波以上方差特征，使得起始时刻功率谱的方差分布更为合理，但对水汽场不做融合。在两者融合过程中，重点考虑在对流层集合资料同化计算的水平相关的功率谱，在实际实现时，比较了不同权重系数对分析预报的影响，最后采用在对流层10%的权重来自集合资料同化样本计算的功率谱，而在平流层这一比例降低到0.2%。对流层和平流层采用双曲正切曲线光滑过渡。融合后的谱空间方差进行调整，确保 $\ \rho$ （0）=1。图2中Merge为误差融合后的起始时刻水平相关功率谱的特征。对流函数、势函数以及无量纲气压，比较SOAR与Merge，可以发现在0—20波，融合后的结果与SOAR的功率谱分布特征一致，而在20波以上功率谱数值要明显大于SOAR，但要略低于EDAFLT功率谱。比较SOAR05与Merge可以发现，在40波以上的短波部分两者接近。

4.2 误差融合对分析场的影响

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

以上分析表明，起始时刻背景误差水平相关变换（U）的特征极为关键。但在目标函数极小化迭代计算中因为起始时刻背景误差水平相关变换（U）的结构与动力模式生成的特征不一致，由式（4）可知，目标函数对控制变量（v）的梯度受到水平相关变换转置矩阵（U^T）的控制，在起始时刻由背景误差水平相关变换（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的分布特征存在不同，两者的水平导数也存在差异（图略），造成两个模型的分析增量特征不同。温度分析增量是通过无量纲气压分析增量计算得到的，并不涉及水平求导计算，两者的分析增量特征较为接近。

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

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 总结与讨论

相对于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波以上的功率谱在终止时刻要明显大于起始时刻，在同化时间窗内则存在明显不一致的不足，但没有发现从小尺度向更大尺度波传播的特征，这一点在未来的研究中还需要进一步诊断误差能量的传播特征来更合理地表述背景误差协方差的演变。