1 引　言

背景误差协方差简称背景误差，它的准确估计直接影响资料同化的分析性能。背景误差一般使用气候历史样本进行统计，并考虑到误差协方差结构的复杂性，采用高度模型化的误差自相关、协相关等模型和各向同性均质误差假设来简化协方差结构，使得背景误差可以显式表达并用于资料同化系统中（Lorenc，1986；Bannister，2008）。这样处理的不足是在变分资料同化的每个分析循环中，虽然数值预报背景场在不断更新，误差特征也在不断变化，但同化框架中背景误差的特征却固定不变。导致背景误差在不同时刻不断变化的因素较多，但其明显受到所在地天气尺度系统的影响，天气系统的斜压不稳定成为背景误差增长的主导因素，呈现出随不同天气流型变化的“流依赖”特征（Kalnay，2003，第5.6节）。如果能在同化框架中引入具备流依赖特征的背景误差，将会帮助同化框架更加有效地确定背景场和观测的相对权重，也能更加合理地将观测信息在空间上传递出去，从而提高分析质量和预报技巧。

随着集合卡尔曼滤波（EnKF，ensemble Kalman filtering）等集合资料同化技术的快速发展与逐步成熟（Evensen，1994；Anderson，2001；Houtekamer，et al，2005），利用集合样本的短期预报来估计集合预报误差协方差（简称集合预报误差）作为背景误差更为合理。从使用集合预报误差角度分析，在现有资料同化方法中，三维变分资料同化（3DVar）使用气候背景误差协方差。4DVar在同化时间窗起始时刻使用气候背景误差协方差，并利用切线性和伴随模式在同化时间窗内实现背景误差的时间演变。基于4DVar，通过扰动观测资料和海温、叠加随机物理过程等来表示观测资料、海洋下边界条件及模式的不确定性，产生多组分析和预报的集合样本来估计集合预报误差，这种方法被称为集合方式的四维变分资料同化方法，简称集合资料同化（EDA，ensemble of data assimilations）方法。在4DVar同化时间窗的起始时刻，同时使用气候背景误差和EnKF或EDA方法生成集合样本估计的集合预报误差，使得在同化时间窗的起始时刻就具有随天气形势演变特征的流依赖动态背景误差，这种方法被称为集合四维变分同化（En4DVar，ensemble four dimensional variational data assimilation），其与4DVar的主要差别在于是否使用集合预报误差。而在同化时间窗内，使用集合预报在时间序列上的多组集合样本来直接估计整个同化时间窗内不同时刻的背景误差，这种方法被称为四维集合变分同化（4DEnVar，four dimensional ensemble variational data assimilation），其与4DVar和En4DVar的主要区别是不再引入切线性和伴随模式。熊春晖等（2013）、Fairbairn等（2014）对这些方法的差异做了详细介绍。

目前，国际上领先的业务变分资料同化系统大多已实现从静态的、表示气候特征的背景误差向基于集合样本估计的、能表示随天气形势演变特征的流依赖动态背景误差的升级（熊春晖等，2013）。英国气象局采用在目标函数中增加扩展控制变量的方式将集合转置卡尔曼滤波（ETKF）估计的集合预报误差与气候背景误差结合，发展了En4DVar系统（Clayton，et al，2013），加拿大气象局也发展了EnKF与4DVar结合的En4DVar同化系统（Buehner，et al，2010a，2010b）。业务数值预报中心的试验结果表明，考虑集合预报误差改善了4DVar的分析质量，特别是在南半球、热带地区，并对台风路径预报也有显著改进。欧洲中期天气预报中心与法国气象局使用基于4DVar的集合资料同化方法来生成集合样本，采用谱滤波来滤除噪声以保留有用的天气系统变化的集合预报误差并用于4DVar，显著地改善了分析和模式预报质量，且正贡献在全球范围10 d预报中都有效（Bonavita，et al，2012）。此外，英国气象局（Lorenc，et al，2015；Bowler，et al，2017）、美国环境预报中心（NCEP）还进一步采用扩展控制变量方式将集合预报误差用于全球4DEnVar（Wang，et al，2013；Kleist，et al，2015a，2015b），加拿大气象局也发展并业务应用了4DEnVar系统（Buehner，et al，2013）。

应用集合预报误差虽然可以通过不同方式实现，但事实证明这些方式之间是等价的（Buehner，2005；Wang，et al，2007）。可以看到，在上述业务实践中，Lorenc（2003a，2003b）提出的增加扩展控制变量方法是引入集合预报误差的常用方案。增加扩展控制变量可以使得分析增量有更大的空间自由度，从而帮助更好地引入观测信息。使用集合样本来估计集合预报误差有优势，但其不足也很明显。Lorenc（2003a）指出，当集合样本数较少时预报误差容易被低估，且存在虚假遥相关，需要进行误差协方差的局地化，而局地化容易破坏集合样本中所包含的动力平衡关系（Kepert，2009）。

GRAPES于2018年6月实现了全球4DVar业务运行（Zhang，et al，2019），中国成为世界上为数不多的采用四维变分同化的全球业务预报中心之一。4DVar方法的基本原理表明，虽然其能在同化时间窗内通过切线性与伴随模式积分来隐式计算流依赖背景误差，但其初始时刻的背景误差一直是气候态的。要克服这一缺陷，需要在4DVar中引入集合预报误差，发展GRAPES全球En4DVar技术，从而实现在同化时间窗起始时刻就应用流依赖背景误差。

在GRAPES全球4DVar基础上进行拓展，有效使用集合样本估计的集合预报误差来改进全球分析质量是本研究的重点。文中主要研究这几个方面的问题：首先，GRAPES全球4DVar采用增量分析方案（Courtier，et al，1994），在具体实施上有自己的特殊性（薛纪善等，2008），有必要研究在该系统框架中增加扩展控制变量的科学方案。虽然已有一些研究基于GRAPES区域3DVar开展了使用集合预报误差的工作（Chen，et al，2015），但全球4DVar中尚未开展相关研究。其次，通过增加集合样本数来降低抽样误差，并对集合预报误差给予较大权重是各业务中心同化技术发展的趋势。考虑到扩展控制变量的空间维数与集合样本数有关，集合样本数越大对应的扩展控制变量的空间维数越大。在大于100个集合样本数的条件下，如何在GRAPES全球4DVar中既能有效使用集合预报误差，又使得扩展控制变量的维数增加受到限制，以减少4DVar系统的计算与存储消耗。最后在GRAPES全球4DVar引入扩展控制变量并对其进行局地化时，需要研究局地化对分析场平衡造成的影响，并发展减缓局地化影响的办法。

2 GRAPES全球4DVar背景误差公式表达

GRAPES全球模式为水平经纬度格点模式。GRAPES全球4DVar的气候背景误差矩阵（B）可以表示为 $ {{B}}={{{{B}}}^{1/2}{B}^{{\rm{T}}/2}\approx {{U}}{{U}}}^{{\rm{T}}} $ 。B^1/2是B的平方根矩阵，在具体实现时通过一系列误差结构函数（或变量变换矩阵）来完成。从极小化求解的控制变量（ $ {{v}} $ ）到分析变量（ $ {{\delta }}{{X}}_{1} $ ）的变换可以表示为

$ {\rm{\delta}}{{X}}_{1}={{U}}{{v}}={{{U}}}_{\rm{P}}{{{U}}}_{\rm{K}}{{{\varepsilon }}_{\rm{b}}{{{U}}}_{\rm{h}}{{U}}}_{\rm{v}}{{v}} $

(1)

式中， $ {\rm{\delta}}{{X}}_{1}=({\rm{\delta}}u,{\rm{\delta}}v,{\rm{\delta}}\pi,{\rm{\delta}}q{)}^{{\rm{T}}} $ 分别是纬向风与经向风、无量纲气压及比湿的分析增量，为 $ n $ 维空间向量。 $ {{v}}=({v}_{\rm{\psi }},{v}_{{\rm{\chi }}_{\rm{u}}},{v}_{{\rm{\pi }}_{\rm{u}}},{v}_{\rm{q}}{)}^{{\rm{T}}} $ 为控制变量，包括流函数、非平衡速度势函数、非平衡无量纲气压及比湿。 $ {{{U}}}_{\rm{h}} $ 是水平变量变换，在全球模式中采用水平谱滤波形式来模拟背景误差水平相关的各向同性分布特征，实现分析增量信息由观测点向周围模式格点的传播。水平谱滤波需要采用水平高斯网格，与水平经纬度网格不同，因而在GRAPES全球4DVar中包含水平经纬度网格与水平高斯网格两套网格。 $ {{{U}}}_{\rm{v}} $ 为垂直变量变换，采用经验正交函数分解（EOF）方式表达背景误差垂直相关结构。 $ {{\varepsilon }}_{\rm{b}} $ 为对角矩阵，表示非平衡部分的背景误差均方根误差。式（1）中水平谱滤波和经验正交函数分解的非平衡分析变量的具体实现过程可由下式表示

$ {\rm{\delta}}{{X}}_{\rm{u}}={{{\varepsilon }}_{\rm{b}}{{{U}}}_{\rm{h}}{{U}}}_{\rm{v}}{{v}}={{{\varepsilon }}_{\rm{b}}{{G}}_{{\rm{s}}\to {\rm{d}}}^{1}{{S}}}^{-1}\sqrt{{{{B}}}_{{\rm{s}}}}{{{E}}{{\varLambda }}}^{\frac{\,1\,}{\,2\,}}{{S}}{{W}}_{\rm{g}}^{-\frac{\,1\,}{\,2\,}}{{v}} $

(2)

式中， $ {\rm{\delta}}{{X}}_{\rm{u}} $ 是非平衡分析变量， $ {{E}}$ 、 $ {{\varLambda }} $ 是背景误差垂直相关矩阵（ $ {{R}}_{\rm{v}} $ ）通过经验正交函数分解求得的特征向量与特征值，其表达与模式垂直层次有关，而与格点空间或谱空间无关，可作用在每个谱系数上。 $ {{S}}$ 、 $ {{{S}}}^{-1} $ 分别是从格点空间到谱空间及从谱空间到格点空间的变换。 $ {{{B}}}_{{\rm{s}}} $ 是谱空间的水平相关系数矩阵。 $ {{W}}_{\rm{g}} $ 是勒让德变换权重系数。式（2）表明， $ {{v}} $ 是垂直特征模态上的投影系数。 $ {{G}}_{{\rm{s}}\to {\rm{d}}}^{1} $ 是从控制变量（ $ v $ ）所在的水平高斯格点到水平经纬度格点的插值算子。

式（1）中 ${U}_{\rm{K}}$ 表示背景误差中变量间的协相关部分，也即平衡算子，可表示为

$ {\rm{\delta}}\pi ={\rm{\delta}}{\pi }_{\rm{b}}+{\rm{\delta}}{\pi }_{\rm{u}}={{N}}{\rm{\delta}}\psi +{\rm{\delta}}{\pi }_{\rm{u}} $

(3)

$ {\rm{\delta}}\chi ={\rm{\delta}}{\chi }_{\rm{b}}+{\rm{\delta}}{\chi }_{\rm{u}}={{M}}{\rm{\delta}}\psi +{\rm{\delta}}{\chi }_{\rm{u}} $

(4)

式中， $ {{N}} $ 是流函数对无量纲气压解释部分的平衡算子，通过风压线性平衡方程或平衡方程与统计回归求解两种方法结合的方式给出（王瑞春等，2015）。 $ {{M}} $ 是流函数对势函数解释部分的平衡算子，通过统计方式给出。通过 $ {{{U}}}_{\rm{K}} $ 变量变换将相互独立非平衡分析变量 $ ({\rm{\delta}}\psi,{\rm{\delta}}{\chi }_{\rm{u}},{\rm{\delta}}{\pi }_{\rm{u}},{\rm{\delta}}q{)}^{{\rm{T}}} $ 变换为分析变量全量 $ ({\rm{\delta}}\psi,{\rm{\delta}}\chi,{\rm{\delta}}\pi,{\rm{\delta}}q{)}^{{\rm{T}}} $ 。考虑到实际观测及模式预报变量是风场，因而还需进行物理变换 $ {{{U}}}_{\rm{P}} $ ，从流函数、速度势函数变换为纬向风与经向风。此外， $ {{{U}}}_{\rm{P}} $ 还包括由静力平衡关系从无量纲气压（ $ {\rm{\delta}}\pi $ ）诊断得出位温（ ${\rm{\delta}} \theta $ ）

$ {\rm{\delta}}{\theta }_{\rm{v}}=\frac{{{c}}_{p}{\theta }_{{\rm{v}}{\rm{b}}}^{2}}{{g}}\left(\frac{\partial \left({\rm{\delta}}\pi \right)}{\partial {{z}}}\right) $

(5)

式中， ${{c}}_{p}$ 是摩尔定压热容，g是重力加速度。

式（2）中经验正交函数分解和水平谱滤波的参数由背景误差垂直相关模型和水平相关模型决定。误差垂直相关矩阵（ $ {{R}}_{\rm{v}} $ ）的构造由垂直相关模型（ $ {\rho }_{{\rm{l}},{\rm{k}}} $ ）决定，与垂直层次所在的高度有关，表示为

$ {\rho }_{{{l}},{{k}}}=\frac{1}{1.0+{{K}}_{{{z}}}{\left({{{z}}}_{{{l}}}-{{{z}}}_{{{k}}}\right)}^{2}} $

(6)

$ {{K}}_{{{z}}}={{{K}}\left(\frac{{{RT}}}{{g}}\right)}^{2} $

(7)

式中，z_l、z_k是任意模式层l、k的平均高度。K是经验系数，当K=0时，所有垂直层次间的相关均为1，当其数值增大时，垂直相关结构逐步变得更加局地。垂直相关结构也可以通过误差样本统计给出（王金成等，2014）。龚建东等（2020）在GRAPES全球4DVar基础上发展了EDA方法生成多组集合样本，也可以用于统计垂直相关结构。

定义新息向量 ${{d}}_{{n}}={{y}}_{{n}}^{\rm{o}}-{\cal{H}}_{{n}}{\cal{M}}{({{X}}}^{\rm{b}})$ ，其中 $ {{y}}_{{n}}^{\rm{o}} $ 是 $ {t}_{{n}} $ 时的观测值，在高分辨率外循环进行计算。定义 $ {\cal{H}} $ 为从模式状态变量到观测空间的观测算子，其线性化算子表示为H。定义 $ {\cal{M}} $ 为非线性预报模式，其切线性模式表示为 $ {{M}} $ 。于是GRAPES全球4DVar的目标函数可以表示为

$ J\left({{v}}\right)\!=\!\frac{\,1\,}{\,2\,}{{{v}}}^{{\rm{T}}}{{v}}\!+\!\frac{\,1\,}{\,2\,}\sum \limits_{{{n}}=0}^{{L}}{\left({{d}}_{{n}}\!-\!{{{H}}}_{{n}}{{M}}{{U}}{{v}}\right)}^{{\rm{T}}}{{R}}^{-1}\left({{d}}_{{n}}\!-\!{{{H}}}_{{n}}{{M}}{{U}}{{v}}\right)\!+\!{J}_{\rm{c}} $

(8)

式中，R是观测误差协方差。 $ {t}_{0} $ 与 $ {t}_{{L}} $ 分别为同化时间窗起始与终止时刻。通过极小化求解式（8）可以求得 $ {t}_{0} $ 时刻分析增量 $ {\rm{\delta}}{{X}}_{1} $ 。 $ {J}_{\rm{c}} $ 是目标函数约束项，采用数字滤波方案约束重力波引发的不平衡（刘艳等，2019）。后文为公式表述简洁而略去 $ {J}_{\rm{c}} $ 项。

由式（1）与（8）可知，在4DVar的 $ {t}_{0} $ 时刻，分析增量（ $ {\rm{\delta}}{{X}}_{1} $ ）利用气候背景误差（B）或变量变化矩阵（ $ {{U}} $ ）只能获得反映气候背景误差的分析增量，而在 $ {t}_{{n}} $ 时刻的分析增量表示为 $ {\rm{\delta}}{{X}}_{{n}}={{M}}{\rm{\delta}}{{X}}_{1} $ ，通过模式动力约束可以部分反映流依赖特征。

3 GRAPES全球4DVar引入集合预报误差 3.1 扩展控制变量的引入

定义预报变量 $ {{Z}}^{\rm{f}}=\dfrac{1}{\sqrt{{{ N}}-1}}\left\{{\rm{\delta}}{{X}}_{1}^{\rm{f}},\;\;\cdots,\;{\rm{\delta}}{{X}}_{{{N}}}^{\rm{f}}\right\} $ ， $ {{N}} $ 为集合样本总数， $ {\rm{\delta}}{{X}}_{{i}}^{\rm{f}} $ 为第 $ {{i}} $ 个集合样本相对于集合平均的偏差，则集合预报误差可以表示为 $ {{P}}_{\rm{e}}^{\rm{f}}={{Z}}^{\rm{f}}{{{Z}}^{\rm{f}}}^{{\rm{T}}} $ 。定义分析增量 ${\rm{\delta}} {{X}}_{2} $ ，它表示由 $ \left\{{\rm{\delta}}{{X}}_{1}^{\rm{f}},\; \cdots,\;{\rm{\delta}}{{X}}_{{{N}}}^{\rm{f}}\right\} $ 为基底所张成 $ {{N}} $ 维空间上的分析增量， $ {\rm{\delta}}{{{X}}_{2}={{Z}}}^{\rm{f}}{{\alpha }} $ 。向量 $ {{\alpha }}=({\rm{\alpha }}_{1},\cdots,{\rm{\alpha }}_{{{N}}}{)}^{{\rm{T}}} $ 表示分析增量 ${\rm{\delta}} {{X}}_{2} $ 在 $ {{N}} $ 维基底上的投影系数。考虑到集合预报误差 $ {{P}}_{\rm{e}}^{\rm{f}} $ 是 $ {{n}}\times {{n}} $ 维矩阵，但其矩阵的秩只有 $ {{N}} $ 维，直接使用 $ {{P}}_{\rm{e}}^{\rm{f}} $ 会造成虚假的遥相关，在距观测位置较远的位置上产生虚假的分析增量。引入“Schur”乘积（元素对元素的乘积，表示为“ $ \circ $ ”），抑制虚假的遥相关，缓解 $ {{P}}_{\rm{e}}^{\rm{f}} $ 的不满秩问题（Lorenc, 2003a，2003b）

$ {\rm{\delta}}{{X}}_{2}=\frac{1}{\sqrt{{{N}}-1}}\sum \limits_{{{i}}=1}^{{{N}}}{{\rm{\delta}}{{X }}}_{{i}}^{\rm{f}}\circ {{{\alpha }}}_{{i}} $

(9)

这时 $ {{{\alpha }}}_{{i}} $ 是 $ {\rm{\delta}}{{X}}_{{i}}^{\rm{f}} $ 对应的局地化分析变量。受限于高性能计算机的存储能力，引入局地化分析变量后需要考虑对其维数进行控制。GRAPES全球模式水平采用Arakawa C网格，垂直方向采用Chaney-Philips模式分层，不同分析变量 $ ({\rm{\delta}}u$ 、 $ {\rm{\delta}}v$ 、 $ {\rm{\delta}}\pi$ 、 $ {\rm{\delta}}q) $ 所在的模式网格空间位置不同，如果对同一集合样本 $ {\rm{\delta}}{{X}}_{{i}}^{\rm{f}} $ 的不同分析变量都取同一个局地化分析变量 $ {{{\alpha }}}_{{i}} $ ，并视质量点δ $\pi$ 所在位置为共同位置，可以降低 $ {\alpha }_{{i}} $ 的长度。

定义扩展控制变量 $ {{{v}}}_{{i}}^{\rm{\alpha }} $ ，与局地化分析变量 $ {{{\alpha }}}_{{i}} $ 的关系是

$ {{{\alpha }}}_{{i}}={{C}}^{1/2}{{{v}}}_{{i}}^{\rm{\alpha }}={{{U}}}^{\rm{\alpha }}{{{v}}}_{{i}}^{\rm{\alpha }} $

(10)

式中， $ {{C}} $ 为局地化矩阵， $ {{C}}^{1/2} $ 为其平方根矩阵，可通过变量变换矩阵 $ {{{U}}}^{\rm{\alpha }} $ 构造。 $ {{{U}}}^{\rm{\alpha }} $ 可进一步表示为垂直局地化和水平局地化矩阵 $ {{{U}}}^{\rm{\alpha }}={{{U}}}_{\rm{h}}^{\rm{\alpha }}{{{U}}}_{\rm{v}}^{\rm{\alpha }} $ 。 $ {{{U}}}_{\rm{h}}^{\rm{\alpha }}{{{U}}}_{\rm{v}}^{\rm{\alpha }} $ 的作用与 $ {{{U}}}_{\rm{h}}{{{U}}}_{\rm{v}} $ 的作用一致，同样可以采用水平谱滤波和垂直经验正交函数分解实现。于是式（9）表示为

$ {\rm{\delta}}{{X}}_{2}=\frac{1}{\sqrt{{{N}}-1}}\sum\limits _{{{i}}=1}^{{{N}}}({\rm{\delta}}u,{\rm{\delta}}v,{\rm{\delta}}\pi,{\rm{\delta}}q{{)}_{{i}}^{{\rm{T}}}}^{\rm{f}}\circ {{{U}}}_{\rm{h}}^{\rm{\alpha }}{{{U}}}_{\rm{v}}^{\rm{\alpha }}{{{v}}}_{{i}}^{\rm{\alpha }} $

(11)

结合式（11），以 $ {{{v}}}_{{i}}^{\rm{\alpha }} $ 为扩展控制变量的目标函数定义为

$\begin{split} J\left({{{v}}}_{1}^{\rm{\alpha }}{,\cdots,{{v}}}_{{{N}}}^{\rm{\alpha }}\right)=&\frac{\,1\,}{\,2\,}{\sum\limits _{{{i}}=1}^{{{N}}}{({{v}}}_{{i}}^{\rm{\alpha }})^{{\rm{T}}}}{{{v}}}_{{i}}^{\rm{\alpha }}+\\&\frac{\,1\,}{\,2\,}\sum\limits _{{{n}}=0}^{{L}}{\left({{d}}_{{n}}-{{{H}}}_{{n}}{{M}}{\rm{\delta}}{{X}}_{2}\right)}^{{\rm{T}}}{{R}}^{-1}\left({{d}}_{{n}}-{{{H}}}_{{n}}{{M}}{\rm{\delta}}{{X}}_{2}\right)\end{split} $

(12)

式（12）相较于式（8）的定义表明，对以 $ {{{v}}}_{{i}}^{\rm{\alpha }} $ 为控制变量的变分同化系统，其计算量大小与集合样本数有关。需要构造 $ {{{v}}}_{{i}}^{\rm{\alpha }} $ ，使得其维数远小于 $ {{{\alpha }}}_{{i}} $ 的维数，以减少存储和计算的负担。从局地化分析变量 $ {{{\alpha }}}_{{i}} $ 到扩展控制变量 $ {{{v}}}_{{i}}^{\rm{\alpha }} $ 的变量变换（称为 $ {{T}}^{\rm{\alpha }} $ 变换）可以表示为

$ {{{{v}}}_{{i}}^{\rm{\alpha }}={{T}}^{\rm{\alpha }}{\bf{\alpha }}}_{{i}}{={{{T}}_{\rm{v}}^{\rm{\alpha }}{{T}}}_{\rm{h}}^{\rm{\alpha }}{{\alpha }}}_{{i}} $

(13)

由于水平与垂直的局地化的尺度都要明显大于气候背景误差相关的尺度，因而水平局地化可以只取前若干波数，垂直局地化只取主导特征向量，这样 $ {{T}}^{\rm{\alpha }} $ 变换是一个降低维数的过程， $ {{{v}}}_{{i}}^{\rm{\alpha }} $ 的空间维数可以远小于 $ {{{\alpha }}}_{{i}} $ 的空间维数。 $ {{{T}}_{\rm{v}}^{\rm{\alpha }}{{T}}}_{\rm{h}}^{\rm{\alpha }} $ 的逆变换是 $ {{{U}}}_{\rm{h}}^{\rm{\alpha }}{{{U}}}_{\rm{v}}^{\rm{\alpha }} $ ，用在式（11）中。

3.2 局地化分析变量实现方式

在GRAPES全球4DVar中实现局地化分析变量有两种方式，一是只进行水平局地化而不进行垂直局地化。此时，扩展控制变量（ $ {{{v}}}_{{i}}^{\rm{\alpha }} $ ）蜕变为水平二维变量， $ {{{\alpha }}}_{{i}}={{{U}}}_{\rm{h}}^{\rm{\alpha }}{{{v}}}_{{i}}^{\rm{\alpha }} $ ， $ {{{v}}}_{{i}}^{\rm{\alpha }} $ 的维数显著减少，因而可以不考虑水平局地化的降维处理。也即，水平局地化所用的谱滤波的波数与4DVar内循环水平谱滤波的波数相同， $ {{{v}}}_{{i}}^{\rm{\alpha }} $ 定义在控制变量（ $ {{v}} $ ）所在的水平高斯网格点上， $ {{{\alpha }}}_{{i}} $ 定义在内循环低分辨率水平经纬度网格点上。这在GRAPES全球4DVar中表示为

$ {{{{\alpha }}}_{{i}}{={{G}}}_{{\rm{s}}\to {\rm{\alpha }}}^{1}{{S}}}^{-1}\sqrt{{{{B}}}_{{\rm{s}}{\rm{\alpha }}}}{{S}}{{W}}_{\rm{g}}^{-\frac{1}{2}}{{{v}}}_{{i}}^{\rm{\alpha }} $

(14)

式中， $ {{{B}}}_{{\rm{s}}{\rm{\alpha }}} $ 表示水平局地化变量在谱空间的水平相关系数矩阵。

另一种办法是进行三维局地化—同时考虑水平和垂直方向的局地化。此时，通过在式（13）的 $ {{T}}_{\rm{h}}^{\rm{\alpha }} $ 水平谱滤波中取更少波数， $ {{T}}_{\rm{v}}^{\rm{\alpha }} $ 垂直经验正交函数分解中取前几个主导特征向量来大幅度缩减控制变量（ $ {{{v}}}_{\rm{i}}^{\rm{\alpha }} $ ）的维数。在GRAPES全球4DVar中表示为

$ {{{\alpha }}}_{\rm{i}}={{{G}}_{{\rm{s}}\to {\rm{\alpha }}}^{2}{{{{S}}}_{\rm{\alpha }}^{-1}}}\sqrt{{{{B}}}_{{\rm{s}}{\rm{\alpha }}}}{{{{E}}}_{\rm{c}}{{{\varLambda }}}_{\rm{c}}^{\frac{1}{2}}}{{{S}}}_{\rm{\alpha }}{{W}}_{\rm{g}}^{-\frac{1}{2}}{{{v}}}_{\rm{i}}^{\rm{\alpha }} $

(15)

式中， $ {{{v}}}_{\rm{i}}^{\rm{\alpha }} $ 定义在比内循环水平高斯网格更低分辨率的水平高斯网格点上，而 $ {{{\alpha }}}_{\rm{i}} $ 定义在内循环低分辨率水平经纬度网格点上。在水平方向采用低分辨率谱空间，则变量 $ {{{v}}}_{\rm{i}}^{\rm{\alpha }} $ 的水平格点数也相应减少。下标 ${\rm{c}}$ 表示对垂直局地化相关矩阵 $ {{C}}_{{\rm{s}}}={{E}}{{{\varLambda }}}^{\frac{1}{2}}{{{E}}}^{{\rm{T}}} $ 进行截断，仅取前若干个主导特征向量与特征值，对应垂直方向相关“更宽”的模态。算子 $ {{G}}_{{\rm{s}}\to {\rm{\alpha }}}^{2} $ 为从高斯格点到低分辨率经纬度格点的双线性插值算子。

在实际应用时，当局地化所采用网格的水平分辨率与内循环的分辨率不一致时，需要重新定义一组网格用于局地化变量，这会使得GRAPES的网格定义更为复杂。一种替代的办法是在水平局地化网格分辨率不变的情况下，仅在垂直局地化时使用较少的主导特征向量。Clayton等（2013）的研究表明，当经验正交函数分解前若干个特征值对总误差的解释比例达到95%时可以获得满意的结果。目前的GRAPES全球4DVar中，水平局地化与内循环的网格一致，垂直方向按照解释比达到90%来截取主导特征向量。试验表明前8个特征值占比已经超过90%，后续试验中取前8个特征模态进行垂直局地化。

3.3 GRAPES全球En4DVar公式表达

采用增加扩展控制变量，是在具有气候背景误差特征的分析增量上增加扩展控制变量给出的具备随流型演变的分析增量，最终获得的分析增量是这两部分增量的加权之和（Lorenc，2003b；Wang，et al，2007）

$ {\rm{\delta}}{{X}}_{1}=\sqrt{{\rm{\beta }}_{\rm{c}}}{{U}}{{v}} $

(16)

$ {\rm{\delta}}{{X}}_{2}=\sqrt{{\rm{\beta }}_{\rm{e}}}\frac{1}{\sqrt{{{N}}-1}}\sum\limits _{{{i}}=1}^{{{N}}}{{\rm{\delta}}{{X }}}_{{i}}^{\rm{f}}\circ {{{U}}}^{\rm{\alpha }}{{{v}}}_{{i}}^{\rm{\alpha }} $

(17)

其中，气候背景误差（ $ {{B}} $ ）与集合预报误差（ $ {{P}}_{\rm{e}}^{\rm{f}} $ ）所占的比重分别为 $ {\rm{\beta }}_{\rm{c}}$ 、 $ {\rm{\beta }}_{\rm{e}} $ ，满足关系 $ {\rm{\beta }}_{\rm{c}}+{\rm{\beta }}_{\rm{e}}=1 $ 。

由新息向量定义，分析值在观测空间投影（ $ y_{{n}} $ ）表示为

$ {{y}}_{{n}}={\cal{H}}_{{n}}{\cal{M}}{({{X }}}^{\rm{b}}+{\rm{\delta}}{{X}}_{1}+{\rm{\delta}}{{X}}_{2}) $

(18)

于是 $ y_{{n}} $ 与观测（ $ {y}_{{n}}^{\rm{o}} $ ）的偏差表示为

$ {{y}}_{{n}}^{\rm{o}}-{{y}}_{{n}}={{d}}_{{n}}-{{{H}}}_{{n}}{{M}}({\rm{\delta}}{{X}}_{1}+{\rm{\delta}}{{X}}_{2}) $

(19)

考虑了扩展控制变量的目标函数表示为

$\begin{split} J\left({{v}},{{{v}}}_{1}^{\rm{\alpha }}{,\cdots,{{v}}}_{{{N}}}^{\rm{\alpha }}\right)=&\frac{\,1\,}{\,2\,}{{{v}}}^{{\rm{T}}}{{v}}+\frac{1}{2}{\sum\limits _{{{i}}=1}^{{{N}}}{({{v}}}_{{i}}^{\rm{\alpha }})^{{\rm{T}}}}{{{v}}}_{{i}}^{\rm{\alpha }}+\\&\frac{\,1\,}{\,2\,}\sum\limits _{{{n}}=0}^{{L}}{\left({{y}}_{{n}}^{\rm{o}}-{{y}}_{{n}}\right)^{{\rm{T}}}}{{R}}^{-1}\left({{y}}_{{n}}^{\rm{o}}-{{y}}_{{n}}\right)\end{split} $

(20)

由式（20）及式（16）—（17），对控制变量（ $ {{v}} $ ）和扩展控制变量（ $ {{{v}}}_{{i}}^{\rm{\alpha }} $ ）的梯度表示为

$ {\nabla }_{\bf{v}}J={{v}}+\sqrt{{\rm{\beta }}_{\rm{c}}}{{{U}}}^{{\rm{T}}}\sum\limits _{{{n}}=0}^{{L}}{{{M}}}^{{\rm{T}}}{{{H}}}_{{n}}^{{\rm{T}}}{{R}}^{-1}\left({{y}}_{{n}}^{\rm{o}}-{{y}}_{{n}}\right) $

(21)

$ {\nabla }_{{{{v}}}_{{i}}^{\rm{\alpha }}}J={{{v}}}_{{i}}^{\rm{\alpha }}+\sqrt{{\rm{\beta }}_{\rm{e}}}{\frac{1}{\sqrt{{{N}}-1}}\left({{\rm{\delta}}{{X }}}_{{i}}^{\rm{f}}\circ {{{U}}}^{\rm{\alpha }}\right)}^{{\rm{T}}}\sum\limits _{{{n}}=0}^{{L}}{{{M}}}^{{\rm{T}}}{{{H}}}_{{n}}^{{\rm{T}}}{{R}}^{-1}\left({{y}}_{{n}}^{\rm{o}}-{{y}}_{{n}}\right) $

(22)

对式（21）和（22）再次求偏导数，计算目标泛函相对于控制变量的二阶偏导数（公式略），构成海森矩阵可以发现，将权重系数 $ {\rm{\beta }}_{\rm{c}}$ 、 $ {\rm{\beta }}_{\rm{e}} $ 放在分析增量的计算中（式（16）、（17）），而不是放在目标函数（式（20））背景误差与集合预报误差项上，主要是改进海森矩阵的条件数，加快变分同化系统收敛。试验表明，在GRAPES全球4DVar中若将权重放在背景误差与集合预报误差项上，当 $ {\rm{\beta }}_{\rm{c}}$ 、 $ {\rm{\beta }}_{\rm{e}} $ 取值差异大时将显著影响极小化迭代收敛。

最后，在GRAPES全球4DVar中具体实现的步骤是：首先由外循环计算新息向量 $ {{d}}_{{n}} $ ，并按照下列步骤计算

（1）给定 $ {{v}} $ 和 $ {{{v}}}_{1}^{\rm{\alpha }}{,\cdots,{{v}}}_{{{N}}}^{\rm{\alpha }} $ 初值为0；

（2）由控制变量（ $ {{v}} $ ）通过式（16）计算 $ {\rm{\delta}}{{X}}_{1} $ ，由扩展控制变量（ $ {{{v}}}_{1}^{\rm{\alpha }}{,\cdots,{{v}}}_{{{N}}}^{\rm{\alpha }} $ ）通过式（17）计算 $ {\rm{\delta}}{{X}}_{2} $ ；

（3）由式（19）计算分析值在观测空间投影（ $ {{y}} $ ）与观测（ $ {{y}}^{\rm{o}} $ ）的偏差（ $ {{y}}_{{n}}^{\rm{o}}-{{y}}_{{n}} $ ）；

（4）由式（20）计算目标函数J（ $ {{v}},{{{v}}}_{1}^{\rm{\alpha }}{,\cdots,{{v}}}_{{{N}}}^{\rm{\alpha }} $ ）；

（5）由式（21）计算 $ {\nabla }_{{{v}}}J $ 并由式（22）计算 $ {\nabla }_{{{{v}}}_{{i}}^{\rm{\alpha }}}J $ ；从而获得目标函数对所有控制变量的梯度（ $ {\nabla }_{{{v}},{{{v}}}_{1}^{\rm{\alpha }}{,\dots,{{v}}}_{{{N}}}^{\rm{\alpha }}}J $ ）；

（6）利用目标函数（ ${{J}} $ （ $ {{v}},{{{v}}}_{1}^{\rm{\alpha }}{,\cdots,{{v}}}_{{{N}}}^{\rm{\alpha }} $ ））及梯度（ $ {\nabla }_{{{v}},{{{v}}}_{1}^{\rm{\alpha }}{,\cdots,{{v}}}_{{{N}}}^{\rm{\alpha }}}J $ ），由最优下降算法求解 $ {{v}} $ 和 $ {{{v}}}_{1}^{\rm{\alpha }}{,\cdots,{{v}}}_{{{N}}}^{\rm{\alpha }} $ ；

（7）回到第（2）步，进入下一次循环。如精度满足需求，则退出循环。

由扩展控制变量（ $ {{{v}}}_{1}^{\rm{\alpha }}{,\cdots,{{v}}}_{{{N}}}^{\rm{\alpha }} $ ）通过式（17）计算 ${\rm{\delta}} {{X}}_{2} $ 时，实际计算是先由式（14）或式（15）计算并存储 $ {{{\alpha }}}_{{i}} $ 值，再由式（9）通过“Schur”乘积计算 $ {\rm{\delta}}{{X}}_{2} $ 。当 $ {{{\alpha }}}_{\rm{i}} $ 是二维空间变量时可以将所有的 $ {{{\alpha }}}_{{i}} $ 叠放形成三维数组，这时数组垂直方向表示不同集合样本对应的局地化分析变量。而当 $ {{{\alpha }}}_{{i}} $ 是三维空间变量时，在大规模高性能计算机上存储 $ {{{\alpha }}}_{{i}} $ 将变得非常困难。在实际计算时对每个集合样本利用式（15）由 $ {{{v}}}_{{i}}^{\rm{\alpha }} $ 计算 $ {{{\alpha }}}_{{i}} $ 后，利用式（9）计算该集合样本对应的分析增量（ $ {\rm{\delta}}{{{X}}_{2}^{i}} $ ），每次对集合样本循环计算时，每增加集合样本后获得的分析增量由下式求得

$ \sum {{\rm{\delta}}{{X}}_{2}^{i}}=\sum {{\rm{\delta}}{{X}}_{2}^{{i}-1}}+{{\rm{\delta}}{{X}}_{2}^{i}}\qquad{{i}}=1,{{N}} $

(23)

采用以上的计算策略时，通过循环集合样本，全部样本只需要存储一个 $ {{{\alpha }}}_{{i}} $ 即可，大幅度减少了对变量存储的需求。然而这种处理需要采用循环计算的方式反复计算 $ {{{{U}}}_{\rm{h}}^{\rm{\alpha }}{{U}}}_{\rm{v}}^{\rm{\alpha }} $ 并依次求得每个 $ {{{\alpha }}}_{{i}} $ ，在程序设计上更为复杂，这与二维空间变量的处理方式有明显不同。

4 水平局地化对地转平衡约束关系的影响 4.1 水平局地化尺度影响

GRAPES全球4DVar的气候背景误差的水平相关模型采用融合水平相关模型（龚建东等，2020），水平局地化模型采用高斯型相关函数。借鉴Xu等（2008a，2008b）提出的在分析时刻前后利用集合样本在时间序列上的扩展来增加集合样本数的方法，试验利用20个NCEP全球集合预报的6与9 h预报的集合样本，并采用扣除对应时间集合平均的扰动样本来减缓因预报时长不一致带来的大尺度场存在的相位误差，总计40个样本。任选2018年7月9日太平洋北部地区一次锋面天气过程进行分析。由集合样本估计出来的集合预报误差如图1所示，在300 hPa高度水平面上（图1a、b，约为模式35层）集合样本估计的气压与风场误差最大的位置位于冷锋槽及槽后位置，在沿54°N的垂直剖面上（图1c、d），气压误差随高度升高而递减，槽后明显大于槽前，纬向风场的误差主要位于风速急流区附近。

将单个观测资料放置在冷锋锋面后的西北平直气流上，位于（54°N，173°W）的模式35层（约250 hPa）位置（图1）。取气候背景误差与集合预报误差的权重分别为0.1与0.9，这时分析增量主要体现集合预报误差的作用，并考虑不同水平局地化尺度参数来研究其对分析增量的影响，试验结果如图2所示。水平局地化采用高斯相关分布特征，在水平局地化尺度的4倍距离的位置，水平相关系数基本降为0。当水平相关尺度参数取值较大时（如28000 km），地球上最远大圆距离（20000 km）处的分析增量受局地化作用将会衰减到原值的70%，水平局地化基本不起作用，可以作为不进行水平局地化的试验结果。由于集合样本数有限，直接由集合样本估计 $ {{P}}_{\rm{e}}^{\rm{f}} $ 会产生虚假的遥相关，在距观测位置较远的距离上会产生虚假的分析增量（图略）。在单点试验中分析增量不仅分布在观测点周围，在远离观测点的位置也有较大的分析增量（图2a）。试验发现当局地化水平相关尺度在1500 km以下时距离观测点较远处的虚假分析增量明显受到抑制。例如：图2b给出的水平相关尺度在850 km时的分析增量，由于所选单点在锋面后且集合预报误差的权重占主导，因而分析增量呈现出随天气流型演变分布的特征—即分析增量沿气流方向延展较长，而在垂直气流方向分析增量范围受到压缩。当不考虑集合预报误差，仅有气候背景误差时在同化时间窗起始（ $ {t}_{0} $ ）分析增量的特点如图9a所示。

4.2 非平衡变量空间局地化

Mitchell等（2002）、Clayton等（2013）研究表明，直接对水平风场及气压进行局地化时，由于集合样本中的水平风与气压变量之间已经包含大尺度准地转平衡约束关系，进行水平局地化容易破坏这种关系。以纬向风（ ${\rm{\delta}} u $ ）为例，其地转分量（ $ {\rm{\delta}}{u}_{\rm{g}} $ ）与气压（ $ {\rm{\delta}}p $ ）的关系可以表示为

$ {\rm{\delta}}{u}_{\rm{g}}=-\frac{1}{{{f}}{\rm{\rho}}}\frac{\partial{\rm{\delta}} p}{\partial {\rm{y}}} $

(24)

式中，f为科里奥利参数， $ {\rm{\rho}} $ 为空气密度。当对无量纲气压和风场分量进行局地化时，等式右侧包含了对局地化函数的偏导数项而左侧没有，从而破坏对准地转平衡关系的表述。Clayton等（2013）研究表明，如果对相互独立的分析变量或非平衡变量进行局地化时，将减缓局地化造成的影响。在GRAPES全球4DVar中，先对 $ {{\rm{\delta}}{{X }}}_{\rm{i}}^{\rm{f}} $ 进行物理逆变换，再将其中所包含的准地转平衡部分抽取出来，仅保留变量间不相关的非平衡部分。

$ {{\rm{\delta}}{{X}}_{{u}}^{{i}}={({\rm{\delta}}\psi,{\rm{\delta}}{\chi }_{\rm{u}},{\rm{\delta}}{\pi }_{\rm{u}},{\rm{\delta}}q{)}_{{i}}^{{\rm{T}}}={{T}}}_{\rm{K}}{{T}}}_{\rm{P}}{{\rm{\delta}}{{X }}}_{{i}}^{\rm{f}} $

(25)

$ {{T}}_{\rm{P}} $ 主要实现从分析变量 $ ({\rm{\delta}}u,{\rm{\delta}}v,{\rm{\delta}}\pi,{\rm{\delta}}q{)}^{{\rm{T}}} $ 到分析变量 $ ({\rm{\delta}}\psi,{\rm{\delta}}\chi,{\rm{\delta}}\pi,{\rm{\delta}}q{)}^{{\rm{T}}} $ 的物理逆变换。 $ {{T}}_{\rm{K}} $ 变换是抽取 $ ({\rm{\delta}}\psi,{\rm{\delta}}\chi,{\rm{\delta}}\pi,{\rm{\delta}}q{)}^{{\rm{T}}} $ 中流函数与势函数平衡关系，以及流函数与无量纲气压平衡关系，给出相互独立的非平衡分析变量 $ ({\rm{\delta}}\psi,{\rm{\delta}}{\chi }_{u},{\rm{\delta}}{\pi }_{u},{\rm{\delta}}q{)}^{\rm T} $ 。此时对非平衡分析变量进行水平局地化，由于平衡部分已经被抽取出来，就不会影响平衡特征。在水平局地化后，重新使用 $ {{{U}}}_{\rm{K}} $ 变量变换，计算流函数对应的气压平衡部分。于是式（11）可以表示为

$ {\rm{\delta}}{{X}}_{2}={{{U}}}_{\rm{P}}{{{U}}}_{\rm{K}}\left\{\frac{1}{\sqrt{{{N}}-1}}\sum\limits _{{{i}}=1}^{{{N}}}\left({{{{T}}_{\rm{K}}{{T}}}_{\rm{P}}{\rm{\delta}}{{X }}}_{{i}}^{\rm{f}}\right)\circ {{{U}}}_{\rm{h}}^{\rm{\alpha }}{{{U}}}_{\rm{v}}^{\rm{\alpha }}{{{v}}}_{{i}}^{\rm{\alpha }}\right\} $

(26)

图3给出了分别在分析变量空间与非平衡分析变量空间进行局地化时，对分析平衡影响的差异。为方便比较起见，同时给出水平局地化相关尺度28000 km的结果，作为不进行水平局地化的对比试验。由图3c、f可见，当不进行局地化时，同化时间窗起始时刻（ $ {t}_{0} $ ）的地面气压分析增量与Z=35层（约250 hPa）的分析增量对应，但随着切线性模式积分，在同化时间窗终止时刻（ $ {t}_{{L}} $ ）分析增量大值区向东北移动，同时在锋面前端出现接近−0.6 hPa的气压增量。由图1c可知，集合样本估计的集合预报误差在模式低层的数值远大于35层，单点分析增量在模式低层出现比观测位置更强的分析增量。而当水平局地化尺度减小时，以水平局地化相关尺度850 km的结果为例，在同化时间窗起始时刻（ $ {t}_{0} $ ），模式第一层同样出现与高层分析增量对应的分析增量大值区，但在同化时间窗终止时刻（ $ {t}_{{L}} $ ），在锋面前端出现−3 hPa的气压增量，并在远离观测位置区域出现大片区域的气压增量。该气压增量围绕观测位置以环形分布，其与观测位置的距离随水平局地化尺度变化。进一步分析发现其位置与高斯水平相关模型二阶导数的最大值处对应。高斯相关模型的二阶导数是表征气压与风场的准地转耦合关系，分析图3d表明，在分析变量空间进行局地化将产生虚假的分析增量。而当在非平衡分析变量空间进行局地化时（图3b、e），可以发现其结果与不进行水平局地化的结果相似，分析增量正值区向东北方向移动，同样在锋面前出现−0.6 hPa的分析增量。上述结果表明将集合预报误差引入GRAPES全球4DVar时，采用在非平衡变量空间进行局地化来保证气压与风场的准地转耦合关系非常重要。

5 三维局地化对分析平衡的影响 5.1 垂直局地化特征尺度

GRAPES全球4DVar采用静力平衡分析方案，根据静力平衡关系通过式（5）由无量纲气压增量计算虚位温增量。垂直局地化的相关结构可以采用式（6）和（7）所给出的垂直相关模型来实现，也可以直接利用模式预报的统计样本来估计。王金成等（2014）比较了垂直相关模型与同一检验时刻的不同预报时效统计样本估计出的垂直相关结构的差异，结果表明统计样本估计的垂直相关结构更加合理。文中采用EDA方法生成的集合样本来估计垂直相关结构（龚建东等，2020）。由集合样本估计的垂直相关结构（ $ {{R}}_{\rm{v}}^{{\rm{s}}} $ ）存在因样本数有限造成的虚假遥相关，文中利用相关范围较宽的垂直相关模型（ $ {{R}}_{\rm{v}} $ ）对距离远的统计噪音进行衰减限制，表示为 $ {{{R}}_{\rm{v}}\circ {{R}}}_{\rm{v}}^{{\rm{s}}} $ ，“ $ \circ $ ”表示“Schur”乘积，结果如图4b所示。采用垂直相关模型给出的垂直相关结构（ $ {{R}}_{\rm{v}} $ ）在不同层次变化较为平滑，而集合样本统计获得的垂直相关结构（ $ {{R}}_{\rm{v}}^{{\rm{s}}} $ ）在边界层相关较强，在边界层之上对流层中层垂直相关迅速减弱，到对流层中高层垂直相关又进一步变强。总体上，统计的垂直相关结构明显窄于相关模型。借鉴Clayton等（2013）的方法，使用流函数的垂直相关结构进行所有变量的垂直局地化，垂直局地化结构表示为 ${{C}}_{{\rm{s}}}=({{{R}}_{\rm{v}}\circ {{R}}}_{\rm{v}}^{{\rm{s}}}{)}^{{P}}$ ，这里使用参数P控制垂直局地化的幅度，当参数P＜1时，垂直相关变宽，试验中取为0.5。这时，局地化垂直相关结构如图4c所示。比较垂直相关模型（图4a）给出的垂直局地化结构，可以发现两者在边界层比较类似，而在对流层中层与高层差异较大。在实际应用垂直局地化矩阵（ $ {{C}}_{{\rm{s}}} $ ）时，为了避免模式层顶附近的层次存在较大噪音的干扰，不考虑模式层顶最高的5层。

进行垂直局地化时，要对无量纲气压增量（ ${\rm{\delta}} \pi $ ）进行局地化，由于式（5）右端涉及到 $ {\rm{\delta}}\pi $ 在垂直方向的梯度，对 $ {\rm{\delta}}\pi $ 进行局地化会显著改变梯度的计算，进而影响到虚位温的计算。但由于虚位温是通过静力平衡关系计算得到的，因而垂直局地化后的分析场仍然满足静力平衡关系。使用单点试验来分析垂直局地化及局地化相关尺度对分析的影响。由图5可见，如果垂直局地化相关尺度较大（式（7）中经验系数K取为1）时，单点观测产生的气压分析增量分布在整个模式层次。对比图1c可以发现，由于集合预报气压误差从观测位置向地面呈现出逐步增大的趋势，垂直局地化对气压分析增量的影响非常大。对比图1d可以发现，纬向风场误差大值区在垂直方向位于模式30—40层（400—200 hPa），相应的纬向风场分析增量也在30—40层，随垂直相关尺度的影响比较小。当垂直局地化相关尺度适当时（K值为7—20），分析增量被约束在观测点附近，气压分析增量与风场分析增量对应较好。采用集合样本估计的垂直局地化结构对气压分析增量在垂直方向的限制作用更强。垂直局地化的影响表明，影响观测信息在空间传播的重要因素来自集合预报误差本身，垂直局地化主要起到约束虚假遥相关的作用。这也表明当考虑集合预报误差在GRAPES全球4DVar中应用时，集合样本质量及其估计的预报误差的合理性，对集合预报误差的应用起到关键作用。

进一步考察垂直局地化对位温分析增量的影响，结果如图6所示。由式（5）可知，位温分析增量一般位于无量纲气压增量（ $ {\rm{\delta}}\pi $ ）的垂直梯度最大的位置。当相关尺度取不同数值时， $ {\rm{\delta}}\pi $ 的垂直梯度发生了明显变化，由静力平衡关系导出的位温分析增量同样也发生明显变化。当垂直局地化相关尺度较大时位温分析增量最大在1 K左右，而当垂直局地化相关尺度较小时位温分析增量可以达到1.4 K左右。采用集合样本估计的垂直相关模型，其对无量纲气压与虚位温分析增量的影响，与垂直相关模型的影响较为一致。这表明垂直局地化结构对分析增量的影响较小，影响较大的是垂直相关尺度，需要合理给定。

5.2 分析场平衡特征分析

4DVar中低分辨率切线性模式积分时存在由诸多因素引起的高频重力振荡，当集合预报误差引入4DVar后，与之相随的水平与垂直局地化也是造成分析不平衡、激发高频重力振荡的原因之一。刘艳等（2019）研究表明，通过在4DVar中引入基于数字滤波方案的初始化过程，可以很好地抑制高频振荡过程，使得分析场在动力上更为平衡。由式（8）， $ {J}_{\rm{c}} $ 是目标函数约束项，其数值反映了对重力波的控制幅度。以单点试验为例比较4DVar与En4DVar，分析不同局地化方案对 $ {J}_{\rm{c}} $ 项收敛的影响，结果如图7所示。4DVar在单点观测的目标函数迭代下降27步后， $ {J}_{\rm{c}} $ 项数值减小到0.14。而在En4DVar中， $ {J}_{\rm{c}} $ 项收敛过程要慢于4DVar，最终收敛数值在0.2左右，也要明显大于4DVar的结果。采用不同局地化相关尺度的结果类似，不过在非平衡分析变量空间进行局地化的结果要略好于在分析变量空间进行局地化。这表明采用集合预报误差后，切线性模式积分产生的高频重力振荡要强于4DVar。不过由于目标函数中 $ {J}_{c} $ 项的存在，分析场仍能保证动力平衡约束关系。

这里将4DVar退化为3DVar，不考虑切线性模式积分的影响，以进一步评估水平与垂直局地化对分析场平衡特征的影响。利用引入集合预报误差的3DVar同化实际常规资料与卫星资料，并用分析场驱动GRAPES全球模式，通过模式启动初期地面气压的时间倾向来分析局地化对分析场平衡特征的影响（图8）。由于背景场中未引入新的同化增量，由其驱动的模式轨迹平衡特征较好，这里作为对照试验。从图8中可以看到，当在分析变量空间进行水平局地化时，模式积分初始时刻的地面气压倾向明显大于背景场启动的结果，且在5 h积分后仍然高于背景场结果。而在非平衡分析变量空间进行水平局地化时，平衡特征明显改善。这进一步说明，在非平衡分析变量空间进行水平局地化更有优势。当在水平局地化基础上进一步考虑垂直局地化时，地面气压倾向增加十分显著。这表明相较于水平局地化，垂直局地化是造成分析不平衡的主要原因。对不同大小的垂直局地化相关尺度，随局地化特征尺度变宽，不平衡程度会略有所减缓，但改善不大。这表明垂直局地化造成的不平衡与垂直局地化的相关尺度有关，但其影响较小。此外，采用集合样本统计给出的垂直局地化结构对分析平衡的改善也有限（图略）。

6 集合预报误差对分析的影响 6.1 集合预报误差权重影响

气候背景误差与集合预报误差所起的作用由其权重系数 $ {\rm{\beta }}_{\rm{c}}$ 、 $ {\rm{\beta }}_{\rm{e}} $ 来决定。当 $ {\rm{\beta }}_{\rm{e}} $ 逐步增大时，集合预报误差逐步占主导。试验利用40个集合样本，当 $ {\rm{\beta }}_{\rm{c}}=1 $ 而 $ {\rm{\beta }}_{\rm{e}}=0 $ 时，同化时间窗起始时刻的分析增量呈现各向同性（图9a）。随着 $ {\rm{\beta }}_{\rm{e}} $ 的值逐步增大，分析增量明显沿着风速较大的方向进行延伸，而与风速垂直的方向分析增量受到抑制（图9c）。起始时刻分析增量呈现出随集合背景误差变化的特征，也即背景误差与背景场的动力演变机制有关，具备随天气系统流型而变的特征。对于仅采用气候背景误差的分析增量，由于切线性模式动力演变机制的存在，在同化时间窗终止时刻分析增量也呈现出随流型变化的特点（图9b），而对于集合预报误差占主导作用的分析增量（图9c、d），终止时刻的分析增量随流型变化的特点更加明显。对比图9c和d，终止时刻的分析增量与起始时刻的分析增量在特征上更加一致，流依赖分析增量的中心位置随气流向下游传播，这是使用集合预报误差对4DVar改进的一个主要特征。

6.2 计算效率影响

对N个集合样本的抽样统计，则对应的蒙特卡洛抽样误差是 $ \dfrac{1}{\sqrt{N}\,\,} $ ，当集合样本数越小时，抽样误差越大。当集合样本数大于100时，抽样误差低于10%（图10）。文中研究的重点是在大于100个集合样本情况下引入扩展控制变量后对计算效率的影响。由于在4DVar目标函数迭代循环的总计算时长中，切线性模式和伴随模式的计算占主导。首先关闭切线性模式和伴随模式将4DVar退化为3DVar，利用实际常规与卫星观测资料评估增加扩展控制变量及不同规模集合样本数对并行计算效率的影响。此时在内循环水平分辨率设为1.0°，气候控制变量（ $ {{v}} $ ）定义在水平高斯网格上，垂直方向取58个特征模态。迭代收敛判据取为目标函数梯度下降到初始梯度的 $ {10}^{-4} $ ，采用256个CPU的情形下，完成一次分析需迭代65步，内循环墙钟时间194 s。在此基础上引入三维扩展控制变量（ $ {{{v}}}_{{i}}^{\rm{\alpha }} $ ），水平方向采用与气候控制变量（ $ {{v}} $ ）一致的高斯网格，垂直方向采用前8个主导特征模态，此时每个集合样本对应的三维扩展控制变量的数目为水平高斯网格数的8倍。当集合样本取20、40、60、120、180个时，气候控制变量与扩展控制变量之和相对于气候控制变量分别增加至1.7、2.4、3.1、5.1、7.1倍，墙钟时间分别增加至1.4、1.8、2.1、2.9、4.1倍。因为垂直局地化所选择的垂直模态数明显小于气候控制变量（ $ {{v}} $ ）的垂直相关模态数，因而扩展控制变量的格点数增加有限。而当采用二维局地化扩展控制变量时，同样20—180个集合样本，气候控制变量与扩展控制变量之和相对于气候控制变量增加至1.1、1.2、1.3、1.5、1.8倍，墙钟时间增加至1.1、1.1、1.2、1.4、1.6倍。比较三维与二维扩展控制变量，在集合样本数较大时二维扩展控制变量的计算效率更高。然而正如前文分析，采用三维扩展控制变量时分析增量在垂直方向的分布更加合理。这表明文中所采用的扩展控制变量引入方式在有限增加扩展控制变量数量的同时，较好地解决了计算量和计算效率的问题，具有实用性。对于采用气候控制变量的4DVar，由于切线性模式和伴随模式的计算占主导，相对于3DVar的194 s，内循环墙钟时间增加到1340 s，计算时间增加至6.9倍。比较采用扩展控制变量的4DVar和仅采用气候控制变量的4DVar，对60个集合样本，引入扩展控制变量后墙钟时间从1340 s增加到1751 s，墙钟时间增加约30%。

7 总结与讨论

GRAPES模式是中国自主发展的非静力格点模式，全球模式采用三/四维变分资料同化作为分析方法。文中研究在现有四维变分资料同化的基础上，有效使用集合样本估计预报误差扩展变分资料同化的分析能力。围绕GRAPES全球4DVar在具体实施上的特殊性，探讨了在引入扩展控制变量过程中，如何在大于100个集合样本数条件下，增加扩展控制变量的同时能显著降低扩展控制变量的维数增加，并减少对计算机存储和计算的压力，提高算法的实用性。此外，重点研究了水平局地化与垂直局地化对分析场平衡性的影响，通过单点试验、常规与卫星观测资料试验，分析了局地化对地面气压倾向的影响。

研究表明，可以充分利用GRAPES全球4DVar气候背景误差结构模型在水平方向采用谱滤波和在垂直方向采用经验正交分解的特点，利用谱滤波实现水平局地化，利用经验正交分解实现垂直局地化。考虑到水平局地化所采用的相关模型更宽，对谱滤波可以采用更少的波数、更粗的高斯网格来降低扩展控制变量的空间维数。考虑到需要独立定义扩展控制变量所依托的网格，在实施上较为复杂，目前GRAPES变分同化系统只用了较宽的水平特征长度模型，没有再定义独立的高斯网格。对垂直局地化，可以通过仅采用前8个主导特征模态来降低空间的维数，此时主导特征模态方差解释占比超过90%。在引入20—180个集合样本的情形下，对二维局地化，控制变量增加至1.1—1.8倍，墙钟时间增加至1.1—1.6倍。对三维局地化，仅考虑垂直局地化降低空间维数，扩展控制变量增加至1.7—7.1倍，墙钟时间增加至1.4—4.1倍。对100—200个集合样本，引入扩展控制变量的计算机墙钟时间增加可控，文中所发展的引入扩展控制变量的方式具备实用性。

为了减少局地化操作对分析平衡的负面影响，发展水平局地化方案时，先将风压准地转平衡关系抽取出来，然后在非平衡分析变量上进行局地化操作，最后再将风压准地转平衡关系加回到非平衡分析变量上。试验结果表明，相比直接在分析变量空间进行水平局地化，本研究发展的在非平衡分析变量空间进行局地化的方案能更好保持分析场的平衡特征。研究还发现，引入垂直局地化会对分析场平衡特征产生较大影响，但依靠GRAPES全球4DVar目标函数中重力波弱约束控制，分析场仍能保证动力平衡约束关系。

气候背景误差与集合预报误差所起的作用由其权重系数 $ {\rm{\beta }}_{\rm{c}}$ 、 $ {\rm{\beta }}_{\rm{e}} $ 来决定。随着 $ {\rm{\beta }}_{\rm{e}} $ 的值逐步增大，单点试验表明分析增量明显沿着风速较大的方向延伸，而在垂直风速方向分析增量受到抑制。分析增量呈现出随天气形势变化的特征。将权重系数 ${\rm{\beta }}_{\rm{c}}{\text{、}} {\rm{\beta }}_{\rm{e}}$ 放在分析增量的计算中，而不是放在目标函数的背景误差与集合预报误差项上，主要是改进海森矩阵的条件数，加快变分同化系统收敛。

致　谢：感谢数值预报中心苏勇、孙健、王金成博士对本研究给予的支持。