基于贝叶斯方案的雷达反射率反演水汽及其同化试验

张诚忠，薛纪善，冯业荣，黄燕燕，戴光丰

(1.中国气象局广州热带海洋气象研究所/广东省区域数值天气预报重点实验室，广东广州510640；2.中国气象科学研究院，北京100081)

1 引 言

水汽是数值天气预报模式模拟各种物理过程中主要的变量与预报对象，比如在云微物理、辐射传输以及对流等过程中水存在的三相变化及其伴随的动力热力效应会进一步影响大气的动力、热力场分布；水汽初值的偏差直接影响云的分布和降水量级的大小，包括水汽要素在内的模式初始场正确性对于预报有决定性影响[1-2]。

雷达资料可以提供高分辨率的降水信息，如果能够很好地将其同化到模式中，对提高中小尺度天气预报能力将起到十分重要的作用。由于观测的雷达反射率不是三维变分系统的分析变量，在同化过程中需要引入变量的物理变换。这种变换可以作为观测算子直接引入同化系统，进行雷达反射率的直接同化；该方案需要涉及到复杂的湿物理过程，过程包含对流启动开、关选项，具有很强的非线性与不连续性，这给伴随模式的编写和实现带来了诸多的不确定性。尽管如此，一些研究[3-5]在伴随模式的处理中，采用了将湿物理过程中的切换时间（switching time）与模式的前向积分时间保持一致的技术处理及雷达反射率观测算子的简化等，利用4Dvar和3Dvar直接同化一部或多部雷达反射率，对比模拟试验显示同化雷达资料后，对对流系统和降水预报起到中性到正作用。另一种同化方案为将雷达反射率反演成大气水汽廓线，再进一步将其引入三维变分同化系统的间接同化；该方案同化过程相对简单，占用计算资源不大，因此一些数值预报业务机构采纳后一种（间接同化）方法[6-7]。

在雷达间接同化过程中关于雷达反射率的水凝物反演研究并不多，Brewster等[8]曾提出一种简单云分析方案，用于由雷达反射率反演云水混合比和雨水混合比。近年来，贝叶斯理论逐渐应用到雷达反射率水汽反演研究上，其思路主要体现为在给定观测资料概率密度函数、背景场概率密度函数（先验密度函数）情况下，求出前两种条件下的条件概率密度函数（后验密度函数）。该理论提供了一套基于概率统计方法来融合同化各种不同来源信息的资料同化方法[9]。Caumont等[10]基于贝叶斯方法建立了雷达反演水汽的一维变分同化系统，通过致洪暴雨个例的模拟试验与分析，认为同化一维反演的水汽后有效提高短期降水预报水平。Wattrelot等[7]则在前者的研究基础上，对反演的水汽廓线进行稀疏化、前期质量控制与后期的质量控制及采用3Dvar同化方法，通过反算一个月的模拟对比试验，结果反映同化雷达水汽后改善了分析场的质量，对前12小时的降水预报有正影响。为充分利用欧盟的所有雷达资料，Martin等[11]统一处理欧洲范围内不同格式的雷达资料，建立一套资料处理与质量控制、同化流程，试验结果与前人的研究结果类似，Martin的工作使得在欧盟范围内从不同国家、不同格式雷达资料投入业务成为可能。

国内有少量类似的反演研究工作，文献[12]考虑在反射率大于30 dBz情况下根据水汽与饱和水汽关系为水汽的正演观测算子进行反演；文献[13]则提出将雷达反射率转化成降水率，通过建立包括大尺度凝结方案和简化对流参数化的线性组合为观测算子反演大气水汽廓线；这些研究成果已逐步应用到基于GRAPES的逐时同化分析与模式预报循环（TheCycleofHourly Assimilation and Forecast，简称“CHAF”）业务系统中[14]。目前国内尚无开展基于贝叶斯理论的雷达水汽反演研究工作，本文拟通过Caumont等[10]提出的贝叶斯方法将观测雷达反射率反演为湿度廓线，经3Dvar将其引入模式，对一个月的华南地区降水预报进行模拟试验，探讨该反演湿度方案在处理华南地区雷达资料的合理性，以及考察该区域雷达反射率资料同化对短临降水预报的影响，为下一步确定获取更合理的水汽信息和短临降水预报效果改进提供经验和技术参考。

2 雷达资料及处理

华南区域精细预报模式范围所包括的雷达如图1所示，其中构成雷达观测网的大部分雷达为Sa或Sb波段，以及少量的Sc波段等共35部地基多普勒雷达，覆盖了模式三分之二的区域范围。雷达资料每隔6分钟收集一次。

雷达资料预处理：对每个体扫的原始资料在PPI面上进行初步的质量控制；根据最近相邻的体扫资料、相邻的几个格点反射率因子(以后简称反射率)判断并剔除反射率奇异点或孤立点，地物回波，空间平滑，之后再将其插值到水平分辨率为1 km的笛卡尔直角坐标的格点上，垂直方向分辨率为模式产品输出格式的标准层共11层(1 000 hPa，850 hPa，……，100 hPa)，以便于下一步的同化反演。

图1 模式范围雷达站点分布图

3 雷达反射率的间接同化方法

雷达反射率的间接同化包括采用贝叶斯方法由雷达反射率反演出水汽（相对湿度）廓线（1D）、并通过GRAPES_3Dvar系统[15]同化生成模式初值场两部分，即1D+3Dvar方法。

3.1 贝叶斯理论

该法则就是对给定未知参数的先验信息(边缘概率密度分布函数P(x)）与观测信息(概率密度分布函数P(y|x))以及边缘概率密度分布函数P(y)，求出后验信息（概率密度分布函数）P(x|y)，表达式如下[16]：

3.2 基于贝叶斯理论的雷达水汽一维反演

文中设矢量Χ代表要估算的模式状态变量，Χture代表真实大气的变量，y0代表一组观测变量。给定一组观测数据y0，模式变量Χ最佳估计可表示：

按照贝叶斯法则，上式可以写成：

在假设模式误差、观测误差遵循高斯分布和非相关的情况下：

其中，R为误差协方差矩阵,h为观测算子。

方程（3）进一步可写成:

根据公式（5）和（6），基于雷达反射率反演的相对湿度可表示如下：

根据公式（7），可以理解为每个模式格点上反演的相对湿度yhu为其周围各个格点（模式预报）相对湿度Xihu的权重线性组合，各格点的权重由实况雷达反射率与模式模拟反射率之间差值决定，反射率差值越小，权重越大。其中，yhu和yz分别为反演出的大气相对湿度、观测反射率；i，j均代表分析范围内所包含模式格点，h为反射率的观测算子，即模式模拟反射率的正演算子。为了更好地模拟反射率，正演计算中考虑了雨水、雪、冰、霰等粒子的反射率因子[17]，误差协方差R在本文中设为0.2 dBz。

3.2.1 一维反演的反演步骤与方法、一致性分析

一维湿度反演前，将观测雷达反射率插值到与模式水平分辨率一致（3 km）的各标准层上；对某一特定的层面上（如700 hPa），依次反演出在以雷达站为中心、半径为150 km的范围内所包括的模式网格点大气相对湿度。针对符合上述范围的模式格点，以该点为中心的18 km×18 km范围内按式（7）反演出该点大气相对湿度。在具体的一维反演中，对如下三种情况分别处理：模式模拟与观测均有回波的情况下，与观测反射率最接近的背景场对应格点的相对湿度所赋予的权重最大，同时抑制与观测反射率相差大的格点对反演湿度的贡献率；对观测有回波、模式模拟无回波情况，将其相对湿度调整为100%，这些调整只局限于高度高于模式抬升凝结高度的点，以下则不进行调整，以确保云底以下降水能在不饱和空气中蒸发；对于两种均无回波的情况下，相当于18 km×18 km范围内的背景场的湿度平均，引进和同化模式和观测点均没有回波点的湿度的主要目的是抑制相对湿度的正增量外扩范围。暂不考虑模式有回波，观测无回波的情况。

一般而言，一维反演方法获取的要素如湿度或反射率与其对应的观测变量一致性、协调性分析是非常必要的[7]。图2分别为在2016年8月1日15时（世界时，下同）的观测场减背景场（O-B）、观测场减分析场（O-A）的反射率差异分布图，其中分析的反射率场由公式（7）计算。显然，雷达观测与背景场模拟的回波差值没有呈高斯分布，存在正偏差，相关系数低（0.65）。经过一维反演生成的雷达分析场，（O-A）差异分布遵循高斯分布，两者相关系数提升至0.85，反映该一维反演方案是合理的，提高了反演变量的一致性。

3.2.2 简单质量控制

由于雷达反射率与相对湿度是非线性关系，即相同的相对湿度对应的降水级别不尽相同（甚至出现有雨、无雨），也就会出现雷达反射率增量为正（分析的雷达-模拟雷达大于0）而相对湿度增量（反演湿度-背景场湿度小于0）为负（符号相反）的点，这些也不是我们所希望同化的资料；为确保进入3Dvar的反演湿度具有更明确的物理意义，对进入3Dvar系统的资料进行简单的质控：只选雷达反射率增量与相对湿度增量符号一致的点；最后对进入三维变分系统的资料进行稀疏化（每隔18 km取一点），减少它们之间的水平相关，并设定相对湿度误差为20%。

图2 2016年8月1日15时雷达观测减背景场（O-B，a）与分析场减观测差异分布图（A-O，b）。

4 个例试验与批量试验

4.1 模式概况与试验设计

文中所用模式为基于GRAPES的华南精细预报模式，其垂直层次55层，水平分辨率3 km。模式垂直坐标采用高度地形追随坐标，水平网格变量配置格式为Arakawa-c格点分布，垂直方向为Charney-Philips跳层设置，计算方案采用半隐式半拉格朗日的时间差分方案。物理过程包括Duhia短波辐射和RRTM长波辐射以及SLAB陆面过程，MRF边界层方案，ARAKAWA对流参数化方案和WSM6微物理过程。

表1 试验设计

为减缓模式冷启动的Spin-up问题，上述的试验均提前3小时进行模式冷启动，将其3小时的预报场为模式初值场，在此基础上进行有无同化雷达资料的对比试验。模式背景场与边界条件均为水平分辨率为0.125°的ECMWF资料。试验中有、无同化资料的试验分别以BYE、CTRL代表。降水预报评估采用威胁指数(TS)，TS=命中/（命中+空报+漏报），预报命中越高则TS评分越高。预报偏差为（命中+空报）/（命中+漏报），偏差接近1为最好。

4.2 台风个例试验结果分析

4.2.1 湿度增量分析

图3为台风“妮妲”在2016年8月1日15时的3 km高度观测的雷达反射率（图3a）与模式预报的反射率（图3b，8月1日12时起报的3小时预报回波，没有同化任何资料）、相对湿度（图3c）。图3d为一维反演的相对湿度场与预报相对湿度（背景场）差值。与实况相比，模式预报的雷达回波没有报出粤东与闽西交界处的回波，粤西的回波范围、强度和落区与实况偏差不小。通过一维反演后，湿度增量主要分布在这两区域，使得该区域的湿度有5%～40%不等的增湿。台风眼西侧、南侧模拟的雷达反射率比实况偏强，反演后该区的湿度减少5%～20%，范围小。以上反映一维反演后背景场反射率（雨带）偏强地区湿度变干（减湿），而反射率偏弱或没有地区的背景场湿度变得更湿（增湿），符合我们预期的结果。

反演的相对湿度通过3Dvar系统分析后，其水汽增量如图4所示。正的水汽增量基本分布在实况回波有而模式模拟没有的区域，以及台风眼附近模式预报回波偏弱的地区；负的湿度增量出现在台风眼西侧南侧附近。显然，在这个个例中1D+3Dvar方法不仅对实况观测有回波而背景场无回波的区域有着不同程度的增湿订正，也订正了背景场湿度偏高的地区。

图3 2016年8月1日15时3 km高度的回波与相对湿度 a.实况回波；b.模式预报的回波(单位：dBz）；c.背景场相对湿度(%)；d.反演同化前后相对湿度增量（单位：%）。

图4 2016年8月1日15时进行1D+3Dvar分析后3 km高度水汽增量 单位：g/kg。

4.2.2 对短临降水预报影响

图5为台风“妮妲”登陆深圳（8月1日20时登陆）前的4小时的雷达回波分布图。期间台风中心沿着西北方向移动，逐渐靠近深圳海岸。观测回波除了出现在台风眼附近并呈环状圆形之外，在粤闽交界、粤西等分别有较强的回波（图5a1～5a4）。其中粤西的强回波中心随时间反时针地向西南方向移动，到19时已移至茂名到海上一带，导致该区域的强降水（图略）。控制试验（图5b1～5b4）中基本模拟出台风中心及其移动方向，但其登陆时间略比实况早，没有模拟出粤西的回波带，模拟的台风眼附近回波带结构比实况紧凑，强度偏强等。经同化雷达反演水汽资料后（图5c1～5c4），粤西的回波在模式运行1小时后能激起对流发展，并在接下来的3小时后回波演变基本与实况吻合，强度范围略小。可见模式初值时刻水汽同化后，有能力模拟出粤闽交界与粤西的强回波带。6小时累积降水预报（图略）显示在同化水汽信息后，明显加强茂名区域的降水预报，闽粤交界的降水也有十几毫米的增加，更接近实况雨带的分布。

逐小时降水TS进一步支持上述分析结果。如图6所示，在大于1 mm降水预报中，同化水汽试验在前9小时的评分TS均高于控制试验TS；这种优势在大于5 mm降水预报上更突出，前12小时的同化试验逐时评分TS分别高于控制试验的逐时TS。

图5 2016年8月1日16—19时的逐小时3 km高度雷达回波 a1～a4.实况；b1～b4.控制试验；c1～c4.同化雷达资料试验（单位：dBz）。

图6 台风“妮妲”登陆期间2016年8月1日15时起报的逐小时降水TS BYE代表同化雷达反演资料，CTRL代表不同化资料，图中纵坐标为TS，横坐标为预报的时间。

4.3 1D+3Dvar同化效果评估-批量试验

4.3.1 1D+3Dvar同化对分析场、预报场相对湿度的影响

图7为有无同化雷达反射率的分析场、3小时、12小时相对湿度预报场平均差异廓线图。同化雷达资料后，模式初值场湿度差异平均廓线增湿（正值）最明显出现在850～400 hPa之间，使得大气对流层中低层湿度变得更大（正的影响）；到3小时和12小时预报，其预报场平均湿度差异正值仍出现在700～400 hPa之间，说明初始场对流层中低层湿度增湿的影响可以延续到12小时预报、到24小时预报，其中低层平均差值变为很小或者负值（数值小，图略）。

图7 2016年5月1—31日期间有无同化试验的00时分析场，3小时、12小时相对湿度预报场平均差异廓线图CTRL为控制试验，BYE为同化试验。单位：%。

4.3.2 对逐时降水预报和2 m温度、10 m风预报影响

图 8为 2016年 5月 1—31日期间的有（BYE）和无（控制试验CTRL）同化雷达反射率1～12小时逐时降水预报平均TS。在大于1 mm的降水预报中，同化试验的降水TS在第1到第12小时均高于控制试验（图略）。在大于5 mm的降水预报TS，前10小时的TS基本上高于控制试验。对大于1 mm的降水预报，同化试验空报偏多；而对大于5 mm的预报，同化试验的预报偏差比控制试验更接近于1，尤其在前10小时的预报（图9）。以上分析可认为同化试验对前12小时的降水预报有正的效果。

相对湿度的变化对2 m温度、10 m风预报的影响有限或很小（图10）。从1个月的1～12小时平均2 m温度RMSE看，同化和控制两试验的RMSE均在1.6～2.4℃之间，同化试验从第2～11小时略比控制试验小，但它们之间相差0.5°以内，优势不明显。同样地，对一个月的10 m风的1～12小时预报评估发现两试验之间的10 m风RMSE相差甚小（图略），可忽略不计。

图8 2016年5月1—31日期间00时起报的1～12小时之平均逐小时降水大于1 mm、大于5 mm的TS CTRL为控制试验，BYE为同化试验，纵坐标为TS，横坐标为预报的时间。

图9 2016年5月1—31日期间00时起报的1～12小时之平均逐小时降水大于5 mm的预报偏差CTRL为控制试验，BYE为同化试验。图中纵坐标为预报偏差值，横坐标为预报的时间。

图10 2016年5月1—31日期间00时起报的1～12小时之平均逐小时2 m温度RMSE CTRL为控制试验，BYE为同化试验。纵坐标为温度RMSE（单位：°），横坐标为预报的时间（单位：小时）。

基于1个月的统计分析认为，一维反演湿度资料同化对华南地区短临降水预报在降水预报偏差保持不变情况下对降水TS有积极影响作用，但对于近地层的风速、2 m温度预报效果甚微或不明显。

5 小结与讨论

雷达资料是目前为数不多有能力为高分辨率预报模式提供高分辨率信息资料之一。为充分利用该资料所包含的丰富中小尺度信息，文中基于雷达反射率与贝叶斯方法反演出大气相对湿度，并将其引入3Dvar系统进行同化分析，为高分辨率模式提供初值场。

通过台风“妮妲”个例试验的分析，得出的初步结果。（1）贝叶斯方法有效地订正了实况有回波而模式预报弱（无）回波区域的大气湿度趋于合理，增加背景场的湿度，减小模拟回波比观测偏强的区域的大气湿度。（2）同化大气湿度后（1D+3Dvar），模式在前6小时预报的台风外围回波分布、演变更合理，改进了降水雨带的分布与强度。

1个月批量试验结果如下。（1）经雷达资料同化后，大气对流层中低层（850～400 hPa）增湿明显，其增湿影响程度可延续12小时以上；（2）其降水预报（大于1 mm和大于5 mm）在前12小时的逐小时降水TS均比控制试验的高，而小于1 mm的预报偏差则变大，空报偏多，大于5 mm更接近1或少变；（3）同化后对地面10 m风场和2 m温度的预报影响不大。

如上所述，1D+3Dvar同化方案对小级别的降水预报空报偏多，为抑制这些空报，文中亦考虑非降水区的湿度同化，即假设雷达反射率小于等于0 dBz为无降水区，这样有可能将一些微弱降水区误判为无降水区，抑制效果仍然不佳；另外，对雷达有明显回波、模式无回波时强制将相应点的大气相对湿度调整为100%，并无动力、热力场的配合，这些均导致同化分析后模式预报空报偏多。为此，在下一步模拟试验中，将同时考虑雷达反演湿度、雷达反演潜热以及风场等信息的同化，以期抑制虚假对流发展，以减少降水空报问题。