引言

近年来，利用新一代多普勒天气雷达反演降水技术已成为雷达资料应用研究的重点，也是发展雷达技术的主要研究方向之一(李柏等，2013)。通过雷达反演降水可及时获取大面积高分辨率降水信息，弥补降水观测站点稀疏的不足。雷达反演降水是精细化气象服务的基础资料，对天气气候的监测和分析、强降雨灾害的评估、数值预报模式的检验和校准以及水文分析等具有很高的应用价值(姚燕飞等，2007；白先达等，2011；万玉发等，2013; 彭涛等，2010；庄薇等，2013)。

雷达定量估测降水是利用雷达回波资料建立雷达回波与降水的Z-I关系，间接获得降水资料。雷达回波是对空中雨滴大小及密度的一种电磁学间接测量，由于雨滴粒径分布及雨滴下落末速度的变化性很大，不同性质的降水适用不同的Z-I关系，即使相同的回波强度所对应的降水强度也可能不同(仲凌志等，2007)。雨滴谱的自然可变性是雷达测量降水不确定性的一个重要原因(万玉发等，2013；李峰等，2014；施红蓉等，2014)。2006年全国自动站资料开始实时上传(任芝花等，2010)，降水观测的时间和空间密度加大，加密降水资料与雷达回波时空分布有很好的对应关系(刘雨佳等，2014；Luo et al，2013)，使雷达估测降水的精度不断提高并稳定，尤其有利于提高积云降水过程降水反演的精度(东高红等，2012b)。目前建立Z-I关系的方法主要分为两种，一种是统计学方法，如采用最优化法(袁晓清等，2010；东高红等，2012a)、概率配对法(郑媛媛等，2004；王丽荣等，2008)等统计多个降水过程资料建立经验Z-I关系。此外，一些研究还根据降水天气形势、不同季节、雷达探测距离和不同降水类型(仲凌志等，2007；陈小华等，2013)对雷达回波进行分类，建立本地化的Z-I关系。为进一步提高估测降水精度，研究者还利用平均校准法(张亚萍等，2013)、变分同化法(马慧等，2008)及卡尔曼滤波校准法(邵月红等，2009；王春青等，2010)进行订正，这些均可在一定程度上改善估测降水的精度，但校正方法一般计算复杂，消耗时间长。另一种方法是动态分级Z-I关系法(汪瑛等，2011)，该方法建立在快速更新资料的基础上，实时频繁变动Z-I关系中的A、b值，不依赖于长期历史资料，只考虑很短时间(6 min)内降水强度和降水性质的延续性。动态分级Z-I关系方法将上一时次雷达回波反射率因子按强度分级，根据最优化方法选取最优的一组A、b值，将该最优分级Z-I关系用于下一时刻估测降水中，该方法能改进对短时强降水的估测，尤其是强降水低估严重的问题(胡胜等，2012)。

本文在以前研究成果的基础上，提出了利用日常业务应用的雷达二次产品回波反射率因子和加密自动站降水资料，确定随时间和空间变化的动态Z-I关系，通过自动站分钟降水数据和雷达回波数据建立动态Z-I关系，来定量估测无自动站观测格点的降水量。与汪瑛等(2011)不同，本文不仅考虑了Z-I关系时间上变化，还考虑了其空间上的动态变化，利用最优化方法得到某时次自动站位置上的A、b值，同时依照与雷达回波强度最接近及与自动站距离最接近的原则，将得到的A、b值延伸应用到该时次无自动站观测的雷达回波格点上，从而实时得到具有业务应用价值的高分辨率定量估测降水资料，为精细化预报服务及强降雨积涝风险评估提供支撑。

1 资料和方法 1.1 资料

本文以石家庄市为研究区域，包括石家庄市区和周边17个县(市)，石家庄位于太行山东麓，地形比较复杂，西部为太行山区，东部为平原。降水资料为石家庄市191个自动站观测的逐分钟雨量数据，雷达数据为石家庄新一代天气雷达CINRAD-SA逐6 min分辨率的雷达二次产品数据，包括组合反射率因子和低层5个仰角(0.5°、1.5°、2.4°、3.4°和4.3°)反射率因子。自动站雨量观测数据和雷达二次产品数据分别经过检验订正和剔除杂波等质量控制，为日常业务应用数据。自动站和雷达站分布见图 1，图中字母为17个县(市)的首字母简称，雷达站位于石家庄新乐市。

选择2009—2012年4年共29个降水过程为研究对象，其中21个降水过程出现在主汛期7—8月，8个过程出现在4、5、6、9月。根据降水性质划分，21个降水过程为对流降水，8个降水过程为层云降水。

1.2 方法 1.2.1 数据预处理

在建立Z-I关系之前，首先将雷达回波与自动站降水数据在时间和空间上建立一一对应关系。第一步，将原始极坐标雷达数据转化到经纬度坐标上，空间分辨率为0.01°×0.01°。第二步，空间上，取自动站周围半径为0.015°范围(9个格点)内雷达回波反射率因子的平均值与区域自动站降水建立数据对(以下简称雷达-自动站数据对)。这与袁晓清等(2010)和梁建茵等(2011)选取方法一致，他们的研究发现自动站上空9个格点的平均雷达回波强度与地面降水强度关系最为密切；时间上，参照梁建茵等(2011)方法，认为雷达资料与该分钟观测到的前6分钟降水量相对应，即雷达资料扫描过程中的降水量，从而得到一一对应的逐6 min雷达-自动站数据对。

1.2.2 动态Z-I关系算法

动态Z-I关系算法是以Z=A×I^b关系式为基础，利用逐6 min的雷达-自动站数据对资料调节A、b系数，使估测降水与区域自动站观测降水的误差最小，从而确定该雷达回波-自动站降水数据对的系数A、b值，并将该系数延伸应用到适用范围内。该算法主要分为以下两个步骤：

(1) 基于Z-I关系，将每组逐6 min雷达-自动站数据对资料应用到该关系式中，调节动态变化的A、b值，选取雷达回波反演降水与自动站降水误差最小时的A、b值为适用于该时刻该组雷达回波-自动站降水数据对Z-I关系的系数。为确保A、b变动范围合理性，限定A、b的调节范围分别为[150, 400]和[0.8, 2.4]，调节间隔分别取0.1和0.01。

(2) 每组雷达-自动站数据对都能得到一组A、b值，每6 min可得到191组A、b值，将A、b值延伸应用到该时刻各个雷达回波格点上，延伸的原则为与雷达回波反射率因子值最接近的合理距离范围内，当该格点内有自动观测站点时，该格点降水直接使用自动站观测降水。

为减小雷达回波从极坐标转换为经纬度坐标可能造成的误差，用雷达格点周围半径为0.015°范围内(9个格点)雷达回波反射率因子的平均值代替该格点值。考虑到中小尺度降水系统的空间尺度，在雷达格点周围半径为0.5°的范围，找到与该雷达回波反射率因子平均值最接近，且差值不大于5 dBz的雷达-自动站数据对，取对应的A、b值，即考虑雷达回波强度的分级；如果0.5°范围内没找到，继续扩展到整个研究区域，同样找与该雷达回波反射率因子平均值最接近，且差值不能大于5 dBz的雷达回波-自动站降水数据对，取对应的A、b值；如果整个区域都没有找到对应自动站，重新在限定区域范围内寻找匹配格点来确定对应A、b值，代入Z-I关系进行降水反演，得到整个研究区域内的0.01°×0.01°分辨率的降水资料。

1.2.3 误差检验方法

交叉检验是检验数据精度的常用评估方法，本文使用交叉检验方法(Chen et al, 2008)对雷达估测降水的精度进行对比检验。首先从所有自动站中随机抽取10%的站点做为检验站点，再从剩下的自动站中抽取10%，重复9次，形成10小组自动站，即每组站点数均为全部站点的10%(本文前9组19个自动站，最后一组20个)，直到所有的自动站均被抽取过一次。考虑地形因素的影响，要求随机抽取的10组自动站空间分布均匀。利用动态Z-I关系方法，利用每个小组自动站抽取后剩余的90%自动站降水与雷达回波数据，生成每个时刻所有格点上的反演降水量，利用检验站点自动观测降水与雷达估测降水进行对比检验，直到所有站点全部完成检验。选择以自动站点滑动1小时累积降水量为单位，将以自动站为中心半径0.015°范围内的雷达估测降水平均值与自动站观测降水进行对比，分析其误差。

2 结果分析 2.1 雷达回波定量估测降水的空间分布

为检验动态Z-I关系估测降水的精度，对比不同仰角雷达估测降水和自动站观测降水的空间分布。图 2和图 3分别给出了2012年7月4日对流性和9月25日稳定性两个降水过程中1小时雷达反演降水和自动站降水空间分布图。2012年7月4日降水空间分布不均匀，局部地区降水强度大，持续时间短，14：36—15：36的1小时降水量主要分布特征为石家庄中心位置有一条南北向，跨越0.5°经度，且降水量大于10 mm·h^-1的降水带；9月25日降水强度较均匀，全区均有降水，但降水持续时间长，10：42—11：42的1小时降水量大部分地区在5 mm·h^-1以下，市南北中间位置有一条东北—西南向降水量大于10 mm·h^-1窄降水带。通过对比图 2和图 3中(b)与(c)~(g)表明，组合反射率因子及5个仰角反射率因子估测降水均能基本反映出地面降水的分布特征，反演出强降水中心，明显改善反演降水偏小的问题，说明动态Z-I关系方法估测降水在空间分布上能抓住地面降水的主要特征，但不同仰角强降水中心的强度和降水范围有所不同。其中，0.5°和1.5°仰角估测降水分布特征与实况最接近，组合反射率因子反演降水显示的强降水范围和强度均较观测实况偏大，2.4°~4.3°仰角则随着仰角的增大，反演降水的强度逐渐减少，降水范围也逐步缩小并集中，已不能显示距离雷达站较远地区[如赞皇(ZH)]的降水分布特征。

从空间分布上看，表明雷达低层0.5°和1.5°仰角反射率估测降水与自动站观测降水最接近，低层仰角的雷达数据对天气监测的应用价值更大，这与陈垚森(2011)的结果一致。以下选取0.5°和1.5°仰角雷达估测降水资料进行对比检验。

2.2 雷达定量估测降水误差分析

利用交叉检验方法综合分析29个降水过程，分析滑动1小时1.5°仰角雷达回波定量估测降水与自动站观测降水的误差，图 4为绝对误差的频率分布图。由图可见，绝对误差在±0.1 mm·h^-1之间的发生频率为41%，绝对误差在±2.5 mm·h^-1之间的发生频率为87%，仅2%的绝对误差超过10 mm·h^-1。从误差分析来看，动态Z-I关系雷达定量估测降水方法得到的降水与地面观测降水比较接近，精度较高。

除综合分析所有降水过程的误差外，根据降水过程的性质，区分对流降水和层云降水，检验动态Z-I方法对不同类型降水的估测能力。计算各过程1.5°仰角反射率反演降水与实况降水的误差和相关系数，结果表明，对流降水和层云降水绝对误差分别为-0.54和-0.65 mm·h^-1，平均为-0.60 mm·h^-1，对流降水误差小于层云降水，估测降水与实况的相关系数达到了0.64和0.66，说明该方法对对流降水的估测能力较层状云稳定性降水的精度偏高。

2.3 不同大小等级降水的估测能力

以前雷达估测降水方法往往无法估测出强降水区域，大雨估测偏小，而强降水对于灾害的预防和评估至关重要(汪瑛等，2011；张亚萍等，2013)。本文对1小时降水按大小分级，检验该方法对不同等级降水量的估测能力。参照(寿绍文等，2006)1小时降水等级表，将自动站实测的滑动1小时降水分为中雨(2~5 mm·h^-1)、大雨(5~10 mm·h^-1)、暴雨(10~20 mm·h^-1)和大暴雨(>20 mm·h^-1)，在一次降雨过程中按照不同降水量级的降水随机选取100个样本，检验雷达回波反射率对不同等级降水的估测能力(图 5)。从图中可见，雷达估测降水与实测降水的变化趋势基本一致，经过相关性t检验，相关系数通过了0.01的显著性水平检验。说明动态估测降水能反映降水实况的时空分布，较好反演出不同等级降水量，显示强降雨中心，说明该方法具有一定的实用性，弥补了强降水明显低估的不足。

表 1给出两种性质降水中不同等级降水的绝对误差分析，可以看出，两种性质降水误差均在一个较小范围，在层云和对流降水过程中一致，说明动态Z-I关系反演降水能抓住降水的强度变化，对10 mm·h^-1以上的降水反演结果有明显的改善，误差明显减少，尤其能表现出20 mm·h^-1的降水大值中心，且对流云降水相对误差小于层云降水，说明动态Z-I关系对对流降水估测能力高于层云降水，这主要表现为对10 mm·h^-1以上大值降水估测能力的提高。

3 个例分析

选取石家庄2011年8月9日强降水过程，利用0.5°和1.5°仰角反射率估测降水，从空间分布和时间序列上与自动站实测降水对比检验。由2011年8月9日15：42—16：42的1小时降水量空间分布图(图 6)可以看到，该小时降水的主要特点是有一条南北向偏西倾斜且穿过石家庄市区，降水量大于10 mm·h^-1的降水带(图 6b)，这与雷达回波分布图(图 6a)上的强回波区域一致。其中，石家庄市区及南部元氏县(YS)存在两个强降水中心(>40 mm·h^-1)，雷达反演降水显示的不同等级降水量面积，形态和落区基本与自动站降水相同，尤其强降水中心的位置和强度也一致，该方法明显改善了反演降水偏小的问题。地面观测降水(图 6b)无井陉(JX)西南角的强降水中心，可能由于该地区为山区，地面无观测数据所致。

除了检验雷达估测降水的空间分布外，选取不同的自动站点，检验了雷达反演降水与自动站观测降水的时间分布序列(图 7)，除B1010、B1014和B1031站为降水双峰极值外，大部分站点观测降水均为明显的单峰极值，反演降水曲线与自动站观测降水形状相近，能显示该站点的降水强度变化，从开始出现降水到小时降水量增大为50~70 mm·h^-1，两者都有较好的拟合。

4 结论与讨论

(1) 本文利用日常业务应用的雷达二次产品组合反射率因子和低层5个仰角反射率因子资料与加密自动站降水资料，建立了随时间和空间变化的动态Z-I关系，利用最优化方法得到某一时次自动站位置上的A、b值，依据与雷达回波强度最接近及距离自动站最接近的原则，将A、b值延伸应用到该时次无自动站观测的雷达回波格点上，从而获取实时高分辨率定量估测降水资料，提高了雷达二次产品资料在业务应用中的价值。

(2) 通过对比检验雷达估测降水的空间分布特征、误差分析及对不同等级降水量的表现能力，表明动态Z-I关系估测降水能抓住降水的时空分布特征，其空间形态和落区基本与自动站观测降水相同，尤其强降水中心的位置和强度两者具有一致性；对不稳定性对流降水和稳定性层云降水的估测精度较高，明显提高了对10 mm·h^-1以上降水的反演能力，能够反演40 mm·h^-1以上的强降水，弥补了一般雷达估测降水方法明显低估强降水的不足，可在业务中实时监测降水过程，这为强降雨积涝及山洪地质灾害等风险评估提供了数据支撑。

(3) 虽然动态Z-I关系估测降水能力有明显的改进，但部分自动站点上的误差仍然存在，这可能由于两者的观测方式、时空分布存在的差异、自动站布局的非均匀性和地形等因素影响所致，另外，雷达资料坐标转换也可能对反演降水产生误差，需要在业务应用中不断优化。