CMA-MESO三维变分同化系统2 m相对湿度资料同化研究*

徐枝芳 吴 洋 龚建东 蔡 怡

1.中国气象局地球系统数值预报中心，北京，100081 2.中国气象科学研究院灾害天气国家重点实验室，北京，100081 3.国家气象中心，北京，100081 4.天津市津南区气象局，天津，300350 5.成都信息工程大学，成都，610225

1 引言

暴雨灾害在中国频发，具有突发性、强度大等特点，常造成人员伤亡和巨大的经济损失，能否准确预报暴雨是社会关注的热点之一。暴雨系统多是在多种水平尺度系统相互作用条件下产生的，如导致暴雨产生的各类天气尺度系统中，大尺度系统为暴雨提供了有利的环境条件，而中、小尺度系统是造成暴雨的直接系统（陆汉城等，2004）。当前，数值预报模式是成功预报暴雨系统的重要技术支撑。提高数值预报模式对暴雨等强对流系统的预报能力可从两方面着手，一方面是提高模式分辨率，开发适合高分辨率模式的物理过程参数化方案和动力框架，增强模式预报精度；另一方面是改善模式初始条件，提高模式对中小尺度系统的刻画能力（Vendrasco，et al，2015）。资料同化利用各种观测资料能为模式提供一个尽可能准确的初值，是提高数值模式预报技巧的一种有效手段。因此，随着数值模式的不断发展完善，资料同化为数值模式提供一个更加精准的初值场十分重要（陈东升等，2004）。

地面观测要素均为模式变量，且相对探空等非常规观测具有高时、空分辨率，可满足数值预报的快速同化发展需求，因此地面观测是快速循环同化系统的重要资料来源，在数值模式中所起的作用也愈加显著。目前中国的2800多个国家级地面自动气象站均为全要素观测站，其资料经过严格质量控制可进入数值模式。由于现有数值预报模式还不能准确地描述大气地表层的过程，且对起伏地形在模式中的表述（模式采用平均地形）与实际地形并不匹配，模式地形高度与观测站高度存在明显的差异，该差异是地面资料同化技术难点之一。中国地形状况复杂，在地面资料同化中不解决这个问题，资料同化有可能带来负效果（徐枝芳等，2007a，2007b，2009）。近年来，一些研究者在寻找解决这两种地形高度差异在地面资料同化中产生负影响的方案，如Lazarus等（2002） 在ADAS系统 （the ARPS（Advanced Regional Prediction System） Data Analysis System）客观分析复杂地形下地面观测和探空观测资料时，考虑了高度差异，并在进行分析时所取的客观分析权重系数中不仅包含了水平距离权重系数，还增加了模式分析层高度与观测站地形高度差权重系数以及模式地形与实际观测地形高度差权重系数，有效地降低了高度差异在资料客观分析过程可能带来的负影响；Devenyi（2003）和Benjamin等（2004）采取气象要素局地递减率的方法，将地面观测要素值由实际地形高度订正到模式地形高度上，解决地面资料同化中地形差异问题；徐枝芳等（2007a，2007b，2009）在Guo等（2002）方法的基础上，将地形高度的代表性误差加入到观测误差中，解决郭永润方法中地形差异问题。张鑫宇等（2021）针对模式与观测地形高度差异问题，对WRF_DA中的地面资料同化订正方案进行优化，提高地面观测气温及风速订正值的合理性，从而改进预报。2 m气温资料相对而言具有较好的连续性，而2 m相对湿度资料变化不连续性增加了同化难度，针对2 m相对湿度同化单独开展的研究工作非常少。

近年来，中国数值预报自主创新发展取得了显著进展，建立了以中国气象局全球同化预报系统（CMA-GFS）（原GRAPES全球同化预报系统，GRAPES_GFS）和中尺度天气数值预报系统（CMA-MESO）（原GRAPES区域中尺度数值预报系统，GRAPES_MESO）为核心的数值预报业务体系（沈学顺等，2020）。CMA-MESO自2006年业务化以来，不断升级改进（徐枝芳等，2013；黄丽萍等，2017；许晨璐等，2017），逐渐取代引进的数值模式系统（李泽椿等，2014；沈学顺等，2020），在日常天气预报中发挥着越来越重要的作用。到2020年，中国区域范围CMA-MESO水平分辨率为3 km，垂直51层，具备了同化雷达、卫星等多种非常规观测资料的能力，而地面资料中仅地面气压（连治华等，2010）、部分2 m相对湿度（观测站地形高于模式地形部分）参与变分同化，地面降水则是用nudging方法进行同化（吴亚丽等，2015）。通过对中国东部区域范围CMA-MESO三维变分同化系统2 m相对湿度同化试验分析发现，超过一半2 m相对湿度资料因为观测站高度低于模式面而不能使用。由于湿度资料是不连续变化，因此如何把观测站高度低于模式面的2 m相对湿度资料同化入模式是一大难点。

对边界层内的气象要素研究（盛裴轩等，2003）发现，在对流边界层的主体部分各种气象要素的梯度都很小，混合层中温度、湿度和风随高度接近均匀分布。模式系统中相对湿度与高度和稳定度是否存在相关？文中采用T639分析场（T639L60全球中期数值预报系统，0.28125°×0.28125°）对2015年夏季（6—8月）华北、长江中下游和华南地区模式最低层和2 m相对湿度差与理查森数（判定大气稳定度）的关系进行统计分析发现，模式最低层和2 m相对湿度的差与理查森数有很好的相关，因此，在该统计结果基础上建立CMA-MESO 3DVar系统2 m相对湿度资料同化方案，并采用该系统进行1个月（2018年7月）连续集成试验分析，评估同化方案的效果。

2 理查森数与近地层相对湿度差关系分析

在白天小风少云的天气下，太阳对下垫面的加热将导致感热通量向上的输送，逐渐形成不稳定层结的边界层，在不稳定层结中大气湍流较强。因此，在对流边界层主体各气象要素的梯度很小，在中等以上不稳定时大气湿度随高度变化较小 （盛裴轩等，2003）。文中利用理查森数作为判定大气稳定度的变量，理查森数（Ri）为（程麟生，1994）

式中，g为重力加速度，Δz为模式最低层到地面的高度差，θvm为模式最低层的虚位温，θm为模式最低层位温，θvo为地面的虚位温，Um、Vm为模式最低层的风场分量。当Ri＜0时，为不稳定层结；当Ri=0时，为中性层结；当Ri＞0时，为稳定层结。

文中采用T639（T639L60全球中期数值预报系统，0.28125°×0.28125°）分析场中等压面分析最低层（1000 hPa）和2 m相对湿度数据。考虑到中国地面观测主要分布在中东部地区，因此统计分析2015年夏季（6、7和8月）华北地区（35°—45°N，110°—120°E）、长江中下游（28°—35°N，110°—120°E）和华南地区（22°—28°N，110°—120°E）模式最低层相对湿度和2 m相对湿度差与理查森数的关系。表1为不同时间不同地区大气不稳定度情况的占比，分为2 m相对湿度不低于90%和低于90%两种情况。4个时次（02、08、14、20时，北京时）分析发现，白天大气更不稳定，08和14时占比较大，20和02时占比较小。

从表1可见，白天大气不稳定，08和14时2 m相对湿度低于90%占比超过70%，晚上大气相对稳定，不稳定情况占比明显较低，华北和长江中下游地区占比低于10%。2 m相对湿度低于90%时的大气不稳定占比比不低于90%时高。大气不稳定情况下模式最低层相对湿度和2 m相对湿度的关系如何？下面以华南地区为例（其他地区结果大体一致）进行分析。

表1 不同时间大气不稳定度所占比例Table 1 Proportions of atmospheric instability at different times

图1a、c为华南地区08、14时2 m相对湿度低于90%时模式最低层相对湿度和2 m相对湿度差与理查森数的散点，此时大气不稳定占比为73.03%（表1）。08和14时2 m相对湿度大于等于90%时（图1b、d），大气不稳定点数明显比低于90%时少（占比38.7%，表1），而20和02时大气稳定情况占比较大（图1e、f）。当理查森数小于0时，相对湿度差集中在±20%，模式最低层相对湿度和2 m相对湿度比较接近，这和对流边界层主体各气象要素的梯度变化很小，在中等以上不稳定时大气湿度随高度变化较小 （盛裴轩等，2003）结论大体一致；理查森数大于0时，模式最低层相对湿度和2 m相对湿度差异则相对比较大。

图1 华南地区模式最低层相对湿度和2 m相对湿度差与理查森数散点（a、b.08时，c、d.14时，e、f.20时，g、h.02时；a、c、e、g.RH＜90%，b、d、f、h.RH≥90%）Fig.1 Differences between relative humidity at the lowest model level and at 2 m height as a function of the Richardson number （a，b.08：00 BT；c，d.14：00 BT；e，f.20：00 BT；g ，h.02：00 BT；a，c，e，g.RH＜90%；b，d，f，h.RH≥90%）

下文依据理查森数小于0时模式最低层相对湿度与2 m相对湿度差为±20%（同化系统给的低层观测误差）这一关系，尝试将观测站高度低于模式地形高度的2 m相对湿度观测资料由观测站高度直接插值到模式最低层高度，以解决2 m相对湿度资料由于模式地形与观测站地形高度差异而导致不能使用问题。

3 CMA-MESO三维变分同化系统2 m相对湿度同化方案设计

文中设计2组试验方案：（1）CTL试验：采用系统中原有2 m相对湿度同化方案。当观测站地形高度大于模式地形高度，将地面观测资料作为高空资料进入同化系统；观测站地形高度低于模式地形高度时将该观测站资料进行剔除，不进入同化系统。（2）RH试验：当观测站地形高度大于模式地形高度，将地面观测资料作为高空资料进入同化系统；当观测站地形高度低于模式地形高度时，依据前文公式计算得到背景低层理查森数，将理查森数小于0的地面观测资料高度由观测站高度处理到模式最低层，然后按照探空湿度同化方式进行同化，理查森数大于0的地面观测资料则和CTL试验一样做不使用处理。CTL试验和RH试验区别在于观测站地形高度低于模式地形高度且背景低层理查森数小于0的地面观测站资料同化。文中试验采用的是CMA-MESO 4.4版，模式系统水平分辨率为3 km，垂直不等间距51层，由于计算资源问题，试验范围为中国中东部（17°—50°N，102°—135°E）。采用的三维同化系统控制变量为U、V风场，地面气压（ps），气温（T）和假相对湿度（RH*）。同化分析资料为探空报（U，V，T，RH（相对湿度）），地面报（ps，RH），船舶报（ps，RH），飞机报（U，V，T），云导风，地基掩星反演可降水量（GPS/PW），雷达速度方位显示反演风廓线（VAD风），雷达径向风，风廓线雷达风，基于云分析系统（朱立娟等，2017）的雷达和卫星资料，地面自动气象站降水nudging同化。试验为每3 h同化分析24 h预报，00和12时（世界时）包含一次冷启动分析预报和一次暖启动分析预报。采用美国国家环境预报中心（NCEP） 0.5°×0.5°FNL 6 h预报场做背景场资料，同化分析预报试验时段为2018年7月1—31日，地面观测资料为国家气象信息中心收集的国家级地面观测站逐3 h气象数据。

4 试验结果分析

4.1 个例试验结果分析

为了分析2 m相对湿度资料同化方案改进前、后的差异，选取2018年7月14日的降水进行个例分析。7月14日00时—15日00时，华北南部、辽宁中东部、吉林中东部、四川盆地西部和北部、华南南部沿海以及海南岛等地的部分地区降大雨到暴雨、局地大暴雨。试验区域地面观测站共计1630个，CTL试验2 m相对湿度资料的利用率为43.2%，RH试验资料利用率为88.6%，约为CTL试验的2倍。从图2可见，RH试验增加的观测资料主要分布在中国地形相对复杂区域，且增加的观测站资料同化后的新息向量（背景减观测）与周围CTL试验原有同化观测的新息向量接近，新息向量的正负基本保持一致。

结合图2、3可知，2 m相对湿度同化对模式最低层的相对湿度分析影响很明显，2 m相对湿度大的区域，相对湿度分析增量为正值，2 m相对湿度小的区域，相对湿度分析增量为负值。低层相对湿度分析增量变化基本在±10%内，比2 m相对湿度新息向量的差值略小。RH试验增加了很多资料，但低层相对湿度分析增量变化不是很大，这可能与图2中新增观测资料的新息向量与CTL试验周围新息向量的正负基本保持一致、差异不是很大有关。

从前12 h降水预报情况（图4）来看，2组试验预报的6 h降水变化很小。对增加观测资料较多的区域分析发现，RH试验预报的6 h降水与实况更加接近，虚假降水强度明显减小（图4a、d紫色圈内），这和相对湿度分析增量减小表现一致，第2节分析表明，理查森数小于0时，2 m相对湿度小于90%占比更大，模式最低层相对湿度和2 m相对湿度比较接近，因而新增资料主要是相对湿度小于90%这部分。相对而言，模式低层相对湿度相对高些（图2中新息向量大于0的观测站），因此新增同化的2 m相对湿度资料降低了图4a、d中紫色圈附近模式低层相对湿度分析，使得虚假降水减少或消失。

图2 2018年7月14日同化2 m相对湿度的新息向量分布 （a.CTL试验，b.RH试验；单位：%）Fig.2 Distributions of 2 m relative humidity innovation （background minus observation） for 14 July 2018（a.CTL experiment，b.RH experiment；unit：%）

图3 模式最低层相对湿度分析增量 （a.CTL试验，b.RH试验与CTL试验差；单位：%）Fig.3 Relative humidity analysis increments at the lowest model level（a.CTL experiment，b.difference between RH experiment and CTL experiment；unit：%）

图4 2018年7月14日6 h累计降水 （单位：mm） 分布 （00—06时：a.CTL，b.RH，c.实况； 06—12时：d.CTL，e.RH，f.实况）Fig.4 6 h accumulated rainfall （unit：mm） from 00：00 to 06：00 UTC 14 July，2018 （a.CTL，b.RH，c.observation） and from 06：00 to 12：00 UTC 14 July 2018 （d.CTL，e.RH，f.observation）

降水常采用TS（Threat Score）、ETS（Equitable Threat Score）和Bias（Bias score）分析（韦青等，2020）。ETS评分是由TS评分改进而来，ETS和TS值越大，Bias值越接近1，表示对降水的预报越准确。根据24 h降水量划分为小雨[0.1 mm，10 mm）、中雨[10 mm，25 mm）、大雨[25 mm，50 mm）、暴雨[50 mm，100 mm）和大暴雨[100 mm，∞）5个量级，6 h降水量划分为小雨[0.1 mm，4 mm）、中雨[4 mm，13 mm）、大雨[13 mm，25 mm）、暴雨[25 mm，60 mm）和大暴雨[60 mm，∞）5个量级。

从24 h降水检验结果（图5）来看，RH试验降水预报明显改善，除大雨外其他量级降水TS评分都比CTL试验高，所有量级降水预报的Bias结果都优于CTL试验。

图5 2018年7月14日00时预报24 h累计降水检验结果 （a.TS，b.Bias）Fig.5 TS and Bias values of 24 h accumulated rainfall simulated in the experimental region at 00：00 UTC 14 July 2018（a.TS，b.Bias）

4.2 连续试验结果分析

在个例分析基础上，开展连续试验进行效果评估。由于多次循环同化会导致00和12时暖启动同化分析的背景场不一致，因此同化分析资料数量变化和分析场差异结果取自冷启动分析试验。从2018年7月1—31日00和12时2 m相对湿度资料参与同化分析资料数量（图6）可见，RH试验在不同高度上（依据观测站气压进行分类统计，大体代表不同地形高度）均有所增多，且低层高于高层，00时高于12时，这与第2节采用T639分析场分析理查森数小于0时，08时（北京时）资料占比高于20时（北京时）结果是一致的。

图6 2018年7月1—31日三维变分同化系统中不同高度同化分析2 m相对湿度资料数量随时间变化（a.CTL试验，b.RH试验）Fig.6 Time-dependent amount of 2 m relative humidity data used in 3Dvar system at different heights 1 to 31 July 2018（a.CTL experiment，b.RH experiment）

从2 m相对湿度新息向量以及分析残差（分析减观测）的垂直剖面偏差和均方根误差（RMSE）（图7）可见，2 m相对湿度新息向量在高层偏差为负，低层偏差为正，分析残差表现大体与新息向量一致，但偏差绝对值明显变小，且均方根误差也减小。RH试验新息向量和分析残差偏差及均方根误差明显改善，正、负偏差的大值区绝对值减小，同时均方根误差大值区减小。

图7 2 m相对湿度新息向量 （a—d） 和分析残差 （e—h） 垂直剖面的偏差 （a、b、e、f） 和均方根误差 （c、d、g、h）（a、c、e、g.CTL，b、d、f、h.RH；单位：%）Fig.7 Vertical cross sections of Bias （a，b，e，f） and RMSE （c，d，g，h） of 2 m relative humidity innovation（a—d） and residual analysis（e—h） （a，c，e，g.CTL；b，d，f，h.RH；unit：%）

由以上分析可知，RH试验通过增加2 m相对湿度同化资料数量，分析残差的偏差和均方根误差大值区的值明显减小，有效改进了数值模式初始分析场。图8和9是冷启动和暖启动两种情况下中国中东部区域24 h降水检验结果。与CTL试验相比，RH试验00时冷启动除大雨外，其他量级降水ETS评分略有提高，偏差减小并接近1。00时暖启动除小雨ETS评分略有下降外，中雨、大雨、暴雨ETS评分均有所增大，偏差在小雨和中雨相当，大雨以上量级则减小。由此可见，改进的2 m相对湿度资料同化方案能提高降水预报技巧。

图8 2018年7月1—31日00时 （世界时） 冷启动24 h累计降水检验结果（a.ETS，b.Bias）Fig.8 ETS （a） and Bias （b） values of 24 h accumulated rainfall simulated in the experimental region for the period from 00：00 UTC 1 to 00：00 UTC 31 July 2018 with cold start

从00时冷、暖启动6 h降水检验ETS评分（图10、11）来看，冷启动时RH试验对前6 h降水预报影响较大，各量级降水ETS评分均有明显提升，随着预报时效延长，对降水的影响逐渐减弱，小雨和中雨基本处于弱的正影响，大雨（13 mm）以上量级的ETS则有些下降，与前面分析的个例结果一致，即RH试验增加的资料主要是相对湿度低于90%部分，对大雨量级以上的降水正影响小。RH试验每3 h同化都有新的2 m相对湿度资料增加，不断改进模式的分析和预报效果（冷启动结果显示前6 h基本为正贡献），因而RH试验暖启动试验逐6 h降水ETS评分检验基本表现为弱的正贡献（大部分量级降水ETS评分在各时段都有所提高）。6 h降水检验偏差（图略）表现与24 h降水检验表现相似，冷启动时2组试验偏差差异明显，RH试验偏差更接近1，而暖启动时2组试验无明显差异。

续图7Fig.7 Continued

图9 同图8，但为2018年7月2—31日00时 （世界时） 暖启动结果Fig.9 Same as Fig.8 but for the period from 00：00 UTC 2 to 00：00 UTC 31 July 2018 with warm start

图10 2018年7月1—31日00时 （世界时） 冷启动逐6 h累计降水ETS评分检验结果（a.00—06时，b.06—12时，c.12—18时，d.18—00时）Fig.10 ETS values of 6 h accumulated rainfall simulated in the experimental region for the period from 00：00 UTC 1 to 00：00 UTC 31 July 2018 with cold start（a.00：00—06：00，b.06：00—12：00，c.12：00—18：00，d.18：00—00：00）

图11 同图10，但为2018年7月2—31日00时 （世界时） 暖启动结果Fig.11 Same as Fig.10 but for the period from 00：00 UTC 2 to 00：00 UTC 31 July 2018 with warm start

5 结论与讨论

中国地面观测资料的站点分布密集、观测频率高、资料实时性强，能很好地反映边界层内大气的动力和热力变化，是高分辨率模式同化分析的重要资料之一，但由于模式地形与观测站地形高度差异的影响，地面观测资料在数值模式中的应用还不够充分。文中分析了2015年6—8月T639（0.28125°×0.28125°）分析场中低层相对湿度和2 m相对湿度差与理查森数的关系，依据分析结果在CMA-MESO三维变分同化系统中建立解决模式与观测站地形高度差异负影响的2 m相对湿度同化方案，并进行1个月（2018年7月）连续数值试验评估新2 m相对湿度同化方案效果。得到以下主要结论：

（1）当理查森数小于0时，相对湿度随高度变化较小，模式低层相对湿度与2 m相对湿度的差在模式低层相对湿度同化观测误差（±20%）范围内。08和14时（北京时）理查森数小于0的占比高于20和02时（北京时）。2 m相对湿度小于90%时的大气不稳定占比比大于等于90%时高。

（2）个例和批量试验同化分析结果表明，新2 m相对湿度同化方案应用后同化分析2 m相对湿度资料数量明显增多，低层多于高层，08时高于20时（北京时），新息向量和分析残差偏差绝对值和均方根误差均减小。

（3）个例和批量试验降水结果分析显示，新2 m相对湿度同化方案方法对降水预报技巧有改善，ETS（TS）评分提高，Bias值减小，冷启动预报对前6 h预报降水改善更为明显，暖启动预报由于不断增加观测资料，新的2 m相对湿度同化对24 h时段内每6 h降水预报都有改善，影响时效明显长于冷启动预报。

致 谢：感谢国家气象中心王瑞春博士、陈炯博士在同化方案实现过程中给予的帮助。