CMA-TYM混合En3DVar方案的设计和初步试验*

瞿安祥 麻素红 张 进

1 中国气象局地球系统数值预报中心，北京 100081 2 灾害天气国家重点实验室，北京 100081

提 要： 基于CMA-TYM的3DVar系统,发展了利用扩展控制变量引入流依赖背景误差协方差(集合扰动成员统计表达)的混合En3DVar同化方案。测试显示,单点台风中心气压数据同化会引起风场非对称增量形成，以及导致传统3DVar方案认为不相关的湿度增量出现。台风个例试验表明,混合En3DVar方案能提取台风内零散观测资料信息，并按实际台风动力特征和分布区域向四周传播，从而影响分析场中台风涡旋强度、结构。同时，与3DVar方案相比，混合En3DVar方案对提高台风路径、强度预报效果明显。

引 言

2012年国家气象中心将基于中尺度数值模式CMA-MESO开发的区域台风模式系统CMA-TYM投入业务运行。近些年来检验表明，该系统120小时内路径强度预报性能接近国际水平，在业务应用中具有重要参考价值。但与此同时，CMA-TYM在进一步提高其预报能力方面遇到了诸多技术难题，一个最主要的困难是如何提高初始台风涡旋质量。虽然存在丰富的、时空密集的卫星探测数据，但受云和降水污染，这些资料并不能被分析系统有效吸收和同化。同时，受限于传统变分同化技术局限性，部分台风探测数据同化反而会对模式预报产生负面影响(Aberson,2008)。

进入21世纪，一种引入集合预报扰动表达流依赖背景误差协方差的混合变分同化方案被发展起来(Hamill and Snyder，2000；Lorenc，2003；Etherton and Bishop，2004；Wang et al，2007；Wang，2010)，相比于传统方案使用气候统计(静态)背景误差参数，混合方案将集合预报扰动引入变分同化系统，显著改善背景误差协方差对实时天气系统描述能力，因而对观测数据信息的提取和传播，更符合实际天气系统分布特征。科学研究(Wang et al，2008a;2008b；Buehner et al，2010a;2010b； Zhang and Zhang，2012;Zhang et al,2013)和业务应用(Clayton et al,2013；Bonavita et al,2012;Wang,2011;Wang et al,2013；Kleist and Ide，2015)都表明，混合变分同化方案会明显提高分析场质量及模式预报水平。随着技术成熟，混合变分同化方案也开始应用于台风模拟研究(Wang,2011；Hamill et al，2011；Li et al,2012)取得不错预报效果。

因此，在CMA-TYM 3DVar基础上，本文设计发展了通过三维Alpha控制变量引入集合预报扰动(隐式流依赖背景误差协方差表达)的混合En3DVar同化方案。然后利用单点台风中心气压数据进行测试，验证En3DVar方案对观测信息的提取和传播，是否与集合预报扰动统计的背景误差协方差特征相符合，以及形成的分析增量是否与实际台风环流分布及动力学属性相匹配。最后基于新建立CMA-TYM En3DVar系统进行台风个例试验，来分析其改进台风初始结构方面的能力，以及提高模式路径和强度预报方面的效果。

1 CMA-TYM 3DVar同化系统

CMA-TYM衍生于中国气象局自主研发的新一代全球与区域同化预报系统GRAPES(Global/Regional Assimilation and Prediction System)，是一个全可压非静力有限差分格点模式(薛纪善等，2008；陈德辉等，2012)。水平方向为经纬网格点球面坐标，垂直方向为高度地形追随坐标。主要特点包括：模式预报变量在水平方向采用Arakawa-C 跳点格式，垂直方向采用非均匀Charney-Philips跳层格式，时间上采用半隐式-半拉格朗日差分方案。CMA-TYM是一个覆盖西北太平洋、南海和部分东亚大陆的区域台风数值模式系统(张进等，2017；麻素红等，2018;麻素红，2019；麻素红和陈德辉，2018)。水平分辨率为0.12°，垂直方向为50层，模式层顶约为33 000 m，参考大气采用水平平均廓线。模式物理过程包括(郭云云等，2015；聂皓浩等，2016；郑晓辉等，2016)：包含水汽、雨、雪、云水、云冰、霰的WSM6显式微物理方案、YSU边界层方案、Monin-Obukhov近地面层方案、Goddard短波辐射、RRTM长波辐射，陆面过程采用SLAB热量扩散方案。另外，考虑0.12°分辨率不能完全描述台风对流特征，Meso-SAS积云对流参数化方案(Pan et al,2014)也考虑在内。

针对CMA-MESO发展的CMA-TYM 3Dvar系统(马旭林等，2009；张华等，2004；庄世宇等，2005；庄照荣等，2006)，其变量水平、垂直分布与模式完全一致。虽然CMA-TYM为非静力模式，预报变量包括水平和垂直风、温度、比湿、无量纲气压(Exner函数)、云水等。但考虑到观测资料属性和变分同化可操作性，3DVar选取水平风(u,v)、位温θ、比湿q、无量纲气压π作为同化分析变量，其中质量场变量：无量纲气压π或位温θ只能二选一(本文选取了无量纲气压π)。

2 混合En3DVar同化方案

以一种简单形式，用权重平均(集合预报扰动统计的和气候统计的背景误差协方差)后的背景误差协方差替换3DVar目标函数中相应项，就可以获得混合En3DVar同化方案直接表达公式：

(1)

其中

x′=x-xb

β1+β2=1

(2)

(3)

满足梯度为零目标函数极小化条件：

HTR-1[H(x′+xb)-yo]=0

(4)

Hx′=H(x′+xb)-H(xb)

(5)

分析增量场最优解为：

HTR-1[yo-H(xb)]

(6)

通过引入扩展控制变量α=(α1,…,αm)T，Lorenc(2003)以一种巧妙方式获得了混合En3DVar同化方案间接表达式：

(7)

可以看出，只要选择合适局地化矩阵模拟方法，基于扩展控制变量技术的混合En3DVar方案很容易在传统3DVar方案上融入实现(具体实现细节见下节)。

3 混合En3DVar同化方案实现

R-1[H(Uυ+xb)-yo]

(8)

接下来为消除分析状态变量(π,u,v,q)之间相关，通过物理算子Up将其变换为互不相关分析变量：流函数ψ、势函数χu(与ψ不平衡部分)、πu(与ψ不平衡部分)、比湿q。然后假设ψ，χu，πu，q自身空间相关结构在水平和垂直方向上具有可分离属性，利用垂直模态映射算子Uv与水平滤波算子Uh的Kronecker积来共同模拟实现。经过上述系列变换，目标函数及其梯度计算就转为对新控制变量υ的迭代求解，分析增量x′的计算表达式为：x′=Uυ=UpUvUhυ。

(9)

x′的计算表达式为：

(10)

目标函数及其梯度的计算就转换为对新控制变量(υ,w)的迭代求解。

3)读取各种类型观测资料yo并设定相应观测误差；

4)计算新息向量：d=yo-H(xb)；

5)设置分析控制变量(υ,w)迭代初值，其中w由m个向量场w1,w2,…,wm组成；

6)计算基于控制变量(υ,w)表达的目标函数值；

7)计算目标函数关于分析控制变量(υ,w)梯度值；

8)基于获得的梯度值，计算梯度下降方向和最优步长，并更新控制变量(υ,w)迭代值。后重复计算目标函数和目标函数梯度值(步骤6和步骤7)，直至达到收敛标准；

9)计算最终分析值：xa=xb+x′。

计算资源方面，混合En3DVar方案利用的是三维α控制变量，在传统3DVar方案分析状态变量总维数为n前提下，会额外增加n×m维分析状态变量空间，且随着集合预报成员数目m增加，空间还在增大，这对计算资源硬件需求是一个比较大的挑战。不过，在采用并行计算环境下，这个需求会随着并行核增加而变得相应容易满足。

4 集合预报扰动数据来源

本试验所需集合预报扰动场数据，源自96个成员的全球(谱模式)集合预报系统T574-EnKF(1)T574-EnKF集合系统是NCEP基于采用EnSRF技术建立的全球集合预报-同化循环滚动业务系统。(Whitaker et al,2008；Wang et al,2013；Kleist and Ide，2015)的，通过CMA-TYM积分获得。具体方案为，首先基于T574-EnKF系统输出的分析场集合，利用降尺度技术获得CMA-TYM初值集合。然后，考虑En3DVar分析更多关注台风内中小尺度信息提取和传播，而降尺度形成的初值仅包含全球模式可分辨的(较粗分辨率)天气尺度信息，因此需通过CMA-TYM模式短时间积分，来获得具有中小尺度信息的集合预报场(Pu et al,2016;Lu et al,2017)。

图1显示的是1306号台风温比亚(Rumbia)形成初期，基于96个集合预报扰动场统计表达的：背景场状态变量π,u,v,q在850 hPa层上误差(集合预报离散度)分布信息。从图中可以看出，各个状态变量误差水平分布范围，具有明显台风环流特征：无量纲气压 (图1a)和风场u(图1b)、v(图1c)都显示越靠近台风中心，误差值越大，并且呈现闭合环流形状。很明显，这些状态变量误差分布与模式对台风系统描述能力吻合，台风内部模拟准确程度要远低于外围环流准确程度。同时，图1d显示的比湿q也符合台风降雨云系分布的螺旋形状(螺旋雨带降水模拟也是常见台风模拟误差来源)。

从图2中可以看出，在垂直方向上，台风中心气压误差与中低层无量纲气压π误差存在正相关，与高层存在负相关，这与台风环流中低层气旋辐合、高层反气旋辐散结构相吻合。同时台风中心气压误差与风速u,v分量误差在不同象限区域显示着有正有负的相关，并且正负值所在区域与台风风压动力学平衡关系相匹配。

图1 由集合预报扰动样本统计的背景场状态变量 (a)π,(b)u,(c)v,(d)q在850 hPa层上的误差分布 (：台风中心位置，下同)Fig.1 Background error distribution of (a) π, (b) u, (c) v, (d) q at 850 hPa level calculated by the ensemble forecast perturbation (：center of typhoon, same as below)

图2 由集合预报扰动样本统计的背景台风中心气压与状态变量(a)π，(b)u，(c)v，(d)q 在模式垂直层上(沿台风中心经、纬向剖面)的误差相关性分布Fig.2 Distribution of model vertical levels (meridional and zonal profile along typhoon center) of background error correlation between typhoon central pressure and (a) π, (b) u, (c) v, (d) q calculated by the ensemble forecast perturbation

5 台风中心气压的单点同化分析

依据前文算法步骤，本文发展实现了混合En-3DVar同化方案。为了验证其合理性和正确性，选取单点台风中心气压数据进行同化测试，并将分析结果和传统3DVar方案进行了对比分析。

图3 同化台风中心气压资料后，3DVar方案形成的(a)π,(c)u增量与En3DVar方案形成的 (b)π,(d)u增量在地面层上的水平分布Fig.3 Surface increment of (a) π and (c) u analyzed by 3DVar and (b) π and (d) u analyzed by En3DVar with typhoon central sea level pressure assimilation

图4 同图3，但为在模式垂直层上经向沿台风中心的分布Fig.4 Same as Fig.3, but for meridional distribution along typhoon center at the model vertical levels

6 实际观测资料同化效果分析

基于新建立的CMA-TYM En3DVar系统，本文进行了实际观测资料同化试验。试验模式区域范围为5°～49°N、93°～156°E，覆盖部分大陆和太平洋、南海的东南亚地区，水平格点数为531×451。试验时间是2013年6月28日12：00 UTC(1306号台风温比亚生成初期)。试验冷启动数据来自T1148-GSI(2)为本文试验提供冷启动数据来源的NCEP全球EnKF集合预报系统和全球确定性预报系统，是一套采用高低双分辨率运行的双向耦合系统，出于业务计算资源和运行时效考虑，EnKF系统采用较低分辨率的T574模式运行，而确定性预报采用高分辨率的T1534模式运行。本文出于试验资源的考虑，确定性预报选择了T1148模式配置运行。全球数值预报模式系统(提供初始场和侧边界)和T574-EnKF全球集合预报系统(提供集合预报场)。同化所需观测资料来自实时数值预报业务常规观测资料集，包括探空、地面站、洋面浮标和船舶资料、飞机报和云导风、洋面风等。除此之外，由预报员实时分析的台风中心海平面气压也作为常规资料进入了分析同化系统。

γ(k1,k2)=exp(-d2/L2)

(11)

式中：d为模式层k1和k2的几何距离，L为垂直局地化特征长度。为获得最优取值，本文对100～1 000 km 范围内的水平局地化特征长度进行了测试，最终发现在260 km情况下会显示非常合理的分析增量。L也采取了类似的测试确定取值为36。

图5、图6显示的是En3DVar系统同化观测资料后，状态变量π,u,v,q在水平(地面层)和垂直方向上(沿台风中心经纬向剖面)的增量分布。从π来看，水平方向上(图5a)在台风中心周围形成一个气压增量负值区；垂直方向上(图6a)，气压增量呈现中低层有效加深，高层反气旋弱增强的分布形式，匹配台风低层辐合高层辐散的质量场属性。从u，v分析增量来看，图5b和5c显示在水平方向上，伴随着气压场加深，围绕台风中心形成明显气旋性环流。实际数据显示增量分布并不均匀，东部、北部区域更大一些，呈现明显非对称加深结构。图6b和6c显示在垂直方向上，u，v分析增量紧密围绕台风中心，且遍布整个中低层，同样非对称特征也非常明显。从湿度q分析增量来看,水平方向上(图5d)其分布区域与台风环流螺旋状云系特征紧密匹配,垂直方向上(图6d)，在台风中心附近产生的分析增量遍布眼区及周围环流区域(看似“杂乱”的增量分布实际和台风中小尺度对流区域十分吻合)。

7 台风预报试验

为评估混合En3DVar方案对台风路径强度预报影响，基于CMA-TYM模式系统对1306号台风温比亚进行了预报试验。试验分两种方案进行，一种为3DVar方案，另一种为En3DVar方案，即分别应用3DVar、En3DVar形成的分析场进行模式积分预报。运行方案为，在“温比亚”生命史期间，每天选取00 UTC和12 UTC时刻进行同化分析和预报试验。预报所需数据同样来自全球模式系统T1148-GSI和全球集合预报系统T574-EnKF。

图5 同化常规观测资料后，En3DVar方案形成的(a)π，(b)u，(c)v，(d)q增量在地面层上的水平分布Fig.5 Surface increment of (a) π, (b) u, (c) v, (d) q analyzed by En3DVar with conventional data assimilation

图6 同图5，但为在模式垂直层上沿台风中心的(a,b)经向、(c,d)纬向分布Fig.6 Same as Fig.5, but for (a, b) meridional and (c, d) zonal distribution along typhoon center at the model levels

图7显示在“温比亚”生命史期间，CMA-TYM模式分别应用3DVar系统和混合En3DVar系统产生的预报路径和实际观测路径情况。从图中可以看出，两种方案对36 h内台风路径预报效果差别不大，但是对长时效的48～96 h预报，混合En3DVar方案表现较好，特别在“温比亚”生命史早期几个预报时刻，相比于3DVar方案，混合En3DVar方案产生的预报路径更趋向于实际观测，有效纠正其右偏趋势。而在“温比亚”生命史后期，尽管混合En3DVar方案优势表现的并不如前期那么明显，但是实际检验数据显示，其路径预报还是好于3DVar方案。

从强度预报检验来看(图8)，无论是中心气压还是近地面最大风速，混合En3DVar方案产生的预报结果都好于3DVar方案。特别在“温比亚”生成初期，混合En3DVar方案就提前72 h准确预报出其在近海快速增强爆发过程，并且这种表现一直延续贯穿后期所有预报时刻，直至台风登陆消亡。反观3DVar方案，尽管也预报描述出台风后期增强过程，但是这个趋势缓慢且不明显，而且在数据量级(气压与风速大小)上也和实际观测相差甚远。

图7 CMA-TYM系统分别基于En3DVar和3DVar方案预报的 1306号台风温比亚移动和实况路径对比Fig.7 Tracks predicted by CMA-TYM model with En3DVar and 3DVar scheme for No.1306 Typhoon Rumbia compared with the observed track

图8 2013年6月28日至7月1日CMA-TYM系统分别基于En3DVar和3DVar方案预报的1306号台风温比亚 (a)中心气压,(b)最大地面风速和实况对比Fig.8 Central pressure (a) and maximum surface wind (b) predicted by CMA-TYM model with En3DVar and 3DVar scheme for No.1306 Typhoon Rumbia compared with the observed data from 28 June to 1 July 2013

8 结论与讨论

基于CMA-TYM 3DVar同化系统，本文设计发展了通过三维Alpha控制变量引入集合预报扰动(隐式流依赖背景误差协方差表达)的混合En3DVar同化方案。在匹配现有3DVar方案算法流程的基础上，具有计算高效、易于实现的特点。

单点台风中心气压数据同化测试表明，相比于3DVar方案，混合En3DVar方案形成的气压场分析增量明显分布于台风环流区域内，风场分析增量结构上更加匹配台风动力学特征：数值量级上更大且非对称结构明显。更重要的是，气压分析增量会导致湿度场(比湿)增量的出现(传统3DVar方案假设两者是不相关的)。由此可见，在集合预报扰动统计表达的背景误差协方差能准确表征背景台风误差条件下,台风内部资料的同化就能够获得与台风环流范围相匹配的分析增量,并且这些分析增量内部之间的相关平衡符合实时台风动力学特征。

应用常规观测资料的同化试验显示,混合En3DVar方案能有效提取台风区域内零散观测信息，并且这部分信息会按实时台风特征向周围传播出去，从而影响分析场台风结构、强度变化。这使得通过有限观测资料信息同化改进初始台风涡旋质量成为可能。同时，台风个例预报试验表明，无论是路径还是强度，混合En3DVar产生的预报结果都要明显好于3DVar预报结果。

需要指出的是，混合变分同化方案形成的分析增量实质是多个集合成员预报扰动的线性组合(被限制在集合成员扰动所张的子空间)，其同化效果与集合预报扰动统计表达的背景误差协方差的准确程度密切相关。只要有限集合预报扰动成员能够抓住实时天气系统在模式格点尺度体现的主要误差特征，那么基于其统计所得的误差协方差矩阵就具有实时天气系统的流依赖型表达能力。