1 引言

利用雷达开展降水估测的研究国内外已经开展了很多年^[1]，主要分为三类：概率配对法、Z-I关系法和基于Z-I关系的数学校正法，这三类方法都是根据雷达观测的反射率因子Z估算降水强度I。Z-I关系法^[2-6]使用雷达反射率因子和雨量计观测的降水量，通过对数线性化拟合法、最优关系法等方法，求得适用于当地雷达的Z-I关系进行降水估测，在业务系统中得到了非常广泛的应用。概率配对法由Miller^[7]提出，该方法假设任一点任一时刻出现雨强为I的可能性与出现雷达反射率因子为Z的可能性相等，根据两者概率相等的原则可配出一系列的Z_i和I_i对，针对每个配对统计确定Z与I的对应关系。郑媛媛等^[8]根据雷达不同距离段不同仰角PPI资料与雨强I建立概率配对序列，然后对所有Z和I的配对利用线性回归统计方法建立它们之间的对应关系。王艳兰等^[9]进一步简化该方法，不刻意寻找雷达反射率因子Z与雨强I之间的函数关系，而是将不同的Z和I进行简单数值对应，他们将回波反射率5~55 dBz以间距5 dBz划分为8个等级，根据统计确立了上述8个回波等级分别对应8个等级降水值。高晓荣等^[10]和梁建茵等^[11]通过订正技术对概率配对法降水估计进行订正，比单纯使用雨量计及与雷达拼图回波强度得到更好的降水估测精度。

汪瑛等^[12]指出一个好的定量降水估测算法，不一定非要建立在长期的观测资料基础上，即定量降水算法的统计代表性可能不在于遍历各种资料，而在于能否实时刻画降水系统的时间变化。基于此种考虑，建立了动态分级Z-I关系法，该方法不需要长期历史资料，将上一时刻每个等级回波反射率从10~75 dBz之间以5 dBz为间隔分为13个等级；根据Z=a_mI^bn进行60×38组的雷达QPE和CTF2，然后根据CTF2最小原则判据选取最优的一组分级Z-I关系用于雷达回波反射率反演降水。该方法只考虑降水强度和降水性质在很短时间内的延续性，因此能精准抓住一些极强的短时强降水的回波点。

本文针对动态分级Z-I关系法中的算法，根据数据的实际范围，提出更高效的算法，并从数学角度推导出动态分级Z-I关系法中使用CTF2做判别式得到结果偏小的幅度，为此提出使用偏差和方差作为判别式的快速动态分级法（Fast Dynamic Categorical metthod，FDC)，该方法得到的结果偏差与均方根误差较小，计算速度快，能较好地表现出降水的时空分布特征，快速动态分级法由于使用6 min的Z-I关系，对较弱回波的降水有一定估测能力。

从时间变化角度看，动态分级Z-I关系法、快速动态分级法和基于Z-I关系的数学校正法属于动态方法，概率配对法、Z-I关系法则属于静态方法。在使用相同时间长度雷达资料进行估测时，动态分级Z-I关系法由于使用的是小时Z-I关系，空估情况少，漏估多，快速动态分级法使用的是逐6 min Z-I关系，情况为空估多，漏估少；使用逐小时所有自动站雨量合计值进行评估时，快速动态分级法由于使用判别函数的差异，偏差、绝对误差和均方根误差等指标都比动态分级Z-I关系法更为优异。

2 资料及其质量控制

广西崇左市位于南宁雷达站西南侧，共有7个国家级自动站和146个区域自动站，各站点与雷达站的距离为40~200 km。本文使用从全国综合气象信息共享平台（CIMISS）获取的两种资料：南宁多普勒雷达（SA）的组合反射率产品和崇左辖区的自动站分钟数据，资料时间段为2016年8月。南宁多普勒雷达每次体扫时间约为6 min，雷达1小时的完整数据可能有9~11个时次。

分钟雨量资料与雷达资料的匹配方式是T-6 min到T-1 min累计雨量与雷达回波反射率反演的降水相匹配，将这一匹配得到的回波等级与降水关系用于T到T+5 min时间段的降水估测，并将各小时累计得到的小时雨量与实况小时雨量进行评估。为保证资料的完整性，只对小时内雷达资料达9个以上时次和该时次无缺测的自动站进行统计，共获得628个时次的样本数据。每小时累计的雨量估测值按下式计算：

(1)

其中QPE_1h为1 h累积的雨量估测值，QPE_i为第i次6 min估测的雨量值，m为估测的次数。

3 快速动态分级法 3.1 动态分级Z-I关系法的算法优化

动态分级Z-I关系法将回波反射率从10~75 dBz之间以5 dBz为间隔分为13个等级，然后以经典的Z=aI^b公式，使用变步长的方法计算出各级强度回波对应的雨强，根据判别函数CTF2找出一组最优的a、b关系，并以这组a、b关系应用到下一个时次的QPE计算。动态分级Z-I关系法思路上强调了雷达回波与降水的关系与最临近时次的相关度优于历史统计平均，这一思路从机器学习角度看待则是一种类推方法，用时空关系最相似的回波降水关系进行估测。CTF2定义如下：

(2)

其中H_i为Z-I关系反演降水，G_i为实测降水，i为自动雨量站序列，N为估测的点数。

由于动态分级法中Z-I关系法只与当前的数据有关，不涉及到历史统计，可以通过优化算法提高计算速度。动态分级Z-I关系法有这样的特点，每级的最优a、b一经确定，这级强度的回波反射率Z即对应的一个确定的雨强I，I的合理范围应当介于本时次观测最小值到最大值之间，即I∈[I_min，I_max]。我们注意到，气象业务中雨强的分度值为0.1 mm，CTF2计算结果与a、b无关，所以计算过程可以跳过a、b的求解过程，使用在合理的区间内逐级回波对应的雨强I作为循环变量，步长=0.1，求出最优Z-I值。这一算法优化将计算由幂运算简化成常数运算，对每一级回波的Z-I关系计算循环次数也从固定的m×n次(如汪瑛等^[4]方法为60×38次循环)减至动态变化的(I_max-I_min)/0.1+1次。分级法采用任何判别函数都可以使用这一优化算法。

3.2 动态分级Z-I关系法整体估测偏小

判别函数CTF2在许多文献里得到应用^[13-16]，一般认为，正是因为判别函数CTF2中包含的一次幂项起到的特殊作用，使求出的Z-I关系实现了最优化^[14]。由于高量级雨量比低量级的雨量少，所以在不分级的情况下，使用CTF2判别函数得到Z-I最优关系求出的估测结果，在大多数情况下可以降低低值区的高估，对高值区的影响幅度较小，降低统计回归效应(statistical regression effect)的影响。但是动态分级Z-I关系法由于各级回波对应的雨量区间幅度变窄，CTF2的作用受到影响。

验证优化后动态分级Z-I关系法的QPE结果，发现得到的QPE值整体偏小。在动态分级Z-I关系法中，同级Z反演降水的值相等，将反演值H_i简记为H，G为这一时间段这一级回波对应自动站观测数据的均值，将H与G的差值记作ΔH，H=G+ΔH。注意到此时H、G、ΔH均为常数，将其代入（2）, 整理可得，

由于，

所以

(3)

由于每一时次的逐级的∑(G-G_i)²不变，所以CTF2只受（ΔH²+ΔH)影响，易得到当ΔH=-0.5时CTF2最小，即只要H=G-0.5即可得到理论上最小的CTF2，当G_{i, min}-G < 0.5时，得到的“理论最优值”会落在观测值区间外。

3.3 快速动态分级法

为解决以上问题，本文提出一种思路源于动态分级Z-I关系法的改进算法，为便于区别，将其命名为快速动态分级法。

对系统性能评估的方法有很多，在放弃使用判别函数CTF2后，本文采用偏差(Bias)和方差(Variance)作为衡量指标。

偏差描述的是预测值的期望与真实值的差距，如式（4），

(4)

方差描述的是预测值的变化范围，离散程度, 如式（5），

(5)

一个好的系统，需要同时衡量偏差和方差，在快速动态分级法中，各级回波对应的雨量均有唯一的值可以同时满足偏差为0且方差最小，证明如下。

对于具体的某一时刻某一级的降水要满足Bias=0，必须有，

可得，

(6)

即对于每一级的估测值H等于该级自动站观测值的均值时，偏差为0。

再研究H与G的方差：

由式（7）可知，只有在H=G的时候方差最小，即H无论取任何值，方差都不可能比H=G的方差更小。综上所述，H=G时，同时满足了偏差等于零且方差最小，这就是快速动态分级法的最优解，所以快速动态分级法各级Z-I关系的计算简化到分级求均值，只需要一个循环即可得到所求时次的Z-I值。这与将每级回波与降水关系通过通常意义上的最小二乘法求得的解是一致的。

4 试验与结果分析 4.1 试验方法与评价指标

对于有可能某些等级范围内的回波在上一时刻没有出现，而下一时刻出现了，在动态分级统计中可能没有该回波等级的Z-I关系对应，动态Z-I关系法使用的是历史统计的Z-I关系结果代替，快速动态分级法则使用最近一次出现的该等级的Z-I关系代替，这一处理方法使得当系统运转开始，只要某一等级对应的雨强出现，即可获得对应的Z-I关系并随着系统运转一直动态更新，使实时估测从根本上脱离了对历史统计数据的使用。

为验证快速动态分级法的估测能力，开展试验进行验证，为描述误差和精度特征，采用了相关系数、偏差、相对误差、均方根误等常用指标，由于QPE还应用在雷达定量降水预报上，所以本文参考预报评价指标，增加TS(Threatscore，风险评分)、空估率、漏估率等指标。

相关系数：

偏差：

相对误差：，

均方根误差：，

TS：，

空估率：，

漏估率：，

式中R_a(i)为第i点上的分析值，R_g(i)为第i点上的观测值，R_a、R_g分别为分析值和观测值的均值，n为计算的点数，N_a、N_b、N_c、N_d定义见表 1。

4.2 试验1

试验1比较的是快速动态分级法与动态分级Z-I关系法，试验共设计三组，第一组使用FDC方法，逐6 min建立Z-I关系用于下一个6 min的雨量估测；第二组使用动态分级Z-I关系法，使用的是上一小时建立的QPE关系估测本小时雨量，记为DC-60 min；第三组使用动态分级Z-I关系法，但使用与FDC方法相同的时间段建立小时Z-I关系，记为DC-6 min。三组建立QPE关系的数据时间段与进行估测时间段的关系见图 1。本文在动态分级Z-I关系法的小时回波反射率采用张培昌等^[8]使用的平均值Z进行计算。

4.2.1 点对点检验分析

表 2给出以阈值=0 mm/h时几种方法对非弱回波（Z小时均值≥10 dBz）和较弱回波（Z小时均值∈[0 dBz, 10 dBz））的晴雨估测结果。

在非弱回波的晴雨降水估测上，TS评价的是估测技巧，得分越高越优，评价结果：FDC>（DC-6 min）>（DC-60 min）。漏估指的是估测降水值≤阈值，但是实况值>阈值，这项指标越小越优，PO评价结果：FDC>(DC-6 min)>(DC-60 min)。空估指的是估测出现大于阈值的降水，但是实况并未出现，这项指标越小越优，FAR评价结果：(DC-6 min)>(DC-60 min)>FDC。快速动态分级法漏估少于动态分级Z-I关系法，但空估较多。

快速动态分级法由于以6 min的回波与降水关系更新Z-I关系，能在一定程度上估测出弱回波造成的降水，但空估率偏高（FAR=0.651），动态分级Z-I关系法无弱回波降水估测能力。

使用偏差等指标进行检验分析，阈值=0 mm/h，结果见表 3。在不分级情况下FDC方法偏差最小，为-0.13 mm/h；相对误差DC-6min最小，与FDC方法差异不大。分级情况下，实况降水≥40 mm/h时为FDC方法较优，实况降水 < 5 mm/h则为DC-6 min的偏差最小；相对误差比较FDC与DC-6 min差别不大，均优于DC-60 min。均方根误差(DC-6 min)>FDC>(DC-60 min)。

DC-6 min所有指标均比DC-60 min好，表明动态分级Z-I法时间更接近的Z-I关系更优，说明汪瑛等^[12]提出的雷达回波与降水的关系与最临近时次的相关度优于历史统计平均思路的正确性。

图 2为1 h雷达定量估测降水与自动站观测降水绝对误差的频率分布，各方法小时QPE误差大部分介于[-2.5，2.5]，快速动态分级法偏于正值区，动态分级Z-I关系法偏于负值区。

由于各自动站的估测结果有差异，使用相关系数对各自动站的降水实况与估测值进行评估，最大值0.955, 最小值0.238，分布见图 3，只有两站的相关系数小于0.5，快速动态分级法估测得到的结果表明大多数自动站相关性较好。

4.2.2 逐小时整组检验分析

用点对点的检验分析只能说明估测方法对单点的估测能力，我们将试验时间段内参与试验的站点小时雨量合计进行分析，逐小时所有站点雨量的合计值在一定程度上能评价估测方法对天气系统在各时间段的整体估测效果。

从表 4看出，当使用整组进行检验分析时，快速动态分级法的偏差、绝对误差和均方根误差指标均为最优。各方法的逐小时整组合计值按时间序列排列估测值走势均围绕实况值，为便于观察，按实况值降序排列（图 4），可见FDC和DC-6 min的估测值与实况值非常接近，DC-60 min相对而言与实况值的离散度较大，DC-6 min基本都比实况值小，DC-6 min的离散度大于FDC。从图 4比较DC-6 min的稳定性和DC-60 min的波动幅度，可以看出随着时间的推移，同等级回波与降水的关系会发生变化，这也许是雷达回波外推进行降水预报的能力随时间增长而快速降低的原因。当将降水估测从时间段上整体考虑，快速动态分级法与实况值的拟合度最优。

4.3 试验2

为验证快速动态分级法对无自动站区域的QPE能力，本文设计一种兼顾交叉验证（Cross-validation）与留出法(hold-out)的测试方法，首先将降水相关系数最小的三个站点摒弃，用剩余的150个自动站采用快速动态分级法进行评估，每轮试验随机分成两组，A组100站，B组50站，分别进行以A、B组互为训练集与测试集更新Z-I关系，再将两组试验结果合并得到合并组数据进行评估，试验共进行十轮，共得到二十组的试验数据和十组合并数据，各站点重复分入A、B组情况见图 5。A、B组不平均分配是根据自动站分布比较稀疏的特点而设计，比较训练集和测试集比例差异较大时候的估测能力，由于测试集未使用训练集的雨量数据更新Z-I关系，多次随机分组试验可以体现快速动态分级法对无自动中区域的实际QPE能力。

表 5为十轮合并统计分析的单点的逐小时分级与不分级评估结果，在训练集与验证集不重复的情况下，偏差为-0.14 mm/h，均方根误差为3.61 mm/h。与试验1对比，可以看出快速动态分级法在使用训练集与验证集合不重合的情况下，得到的结果会稍微变差，但是幅度不大，这一试验表明快速动态分级法在对无自动站区域的降水估测能力良好。

表 6给出了A、B两组逐次分组与合并试验的整组分析结果。试验看出，每轮次试验同组中RAE、RMSE、r波动幅度不大，Bias受训练集合测试集的分组影响有波动。

比较A、B两组，可见当测试集与训练集的比例变大时，整组的估测性能稍差，由于现实情况中，无自动站区域比有自动站区域大得多，说明无自动站区域的估测能力可能有降低。

为与试验1比较，列出A、B合并组Bias差异最大的两组数据（图 6）, 可见快速动态分级法的QPE能力比较稳定。

5 结论

本文研究表明，在算法设计时，应当优先考虑到数据的实际意义，据此设计的算法能大幅度提高计算效率。通过数学推导的方式，给出了以偏差和方差为判别式的快速动态分级法的最优解，将计算强度降低到一次循环分级求均值，并且避免使用判别函数CTF2造成得到的QPE结果偏小，随着计算时间跨度的精细化，多次累积后总体偏小严重的情形。

不分级情况下点对点快速动态分级法得到的QPE值偏差为-0.13 mm/h，均方根误差为3.45 mm/h。快速动态分级法由于使用的是逐6 min雨量数据与雷达回波反射率因子建立Z-I关系，对Z-I关系的变化比较敏感，对弱回波降水具有一定估测能力，在晴雨估测时快速动态分级法漏估少于动态分级Z-I关系法，但是空估较多。将逐小时整组自动站的雨量合计值作为整体考量时，快速动态分级法降水估测效果表现最优。

设计兼顾交叉验证与留出法的测试方法的多次随机分组试验表明，用于建立Z-I关系的自动站都不重合的时候，估测结果比使用所有自动站的进行QPE得到的结果稍差，但幅度不大，这说明快速动态分级法对于无自动站区域的估测能力值得肯定。

由于计算速度快，快速动态分级法可以进一步应用到机器学习领域作为数据预处理的方法。