引言

随着社会的发展，各行各业对气象服务的需求越来越多，要求也越来越高。同时，随着高时空分辨率的数值预报模式在预报业务中的使用，人们希望气象部门能提供更为准确的天气预报产品。然而，受模式误差、输出误差、模式稳定性以及预报员的主观分析等原因影响，如果直接使用数值预报模式中的输出结果进行预报，其准确率相对较低。为了得到客观化和定量化的预报效果，目前，对数值模式产品进行客观订正的释用是解决该问题的主要方法之一。实践证明：通过对数值预报产品的释用，其预报能力比模式的直接输出预报有明显提高^[1]。同时，由于暴雨是我国最重要气象灾害之一，开展客观定量的数值预报产品暴雨预报释用方法研究具有重要意义。

到目前为止，很多学者开展了利用数值预报产品进行客观定量降水预报(quantitative precipitation forecast，QPF)的释用方法研究^[2-6]，包括强降水的客观定量预报释用方法研究，取得了较多的研究成果，并在实际预报业务中得到应用。如David等^[7-8]利用分级技术从NCEP多模式短期集合预报系统中输出较可靠的3 h概率定量降水预报(PQPF)，并对此进行分析发现分级技术对降水事件能提供更多预报技巧和可靠性；赵声蓉等^[9]以多个模式的降水量预报产品作为神经网络的预报因子，利用神经网络方法建立一种数值预报产品释用的定量降水预报模型，预报效果提高明显；唐晓文等^[10]根据深厚湿对流系统长时间维持将产生强降水这一配料法的预报原理以及中国不同区域的天气特点，选取对强降水(不小于25 mm·d^-1)有显著影响的4类因子(水汽因子、动力因子、不稳定因子及热力因子)，在一定的物理条件约束下，利用经验和统计相结合的方法建立配料综合指数与强降水之间的关系；Jerome等^[11-12]利用美国国家气象局气象发展研究室开发的高分辨率全球预报系统(Global Forecasting System，GFS)得到美国本土的4 km格距的高分辨率MOS定量降水预报；孙靖等^[13]利用分等级消除偏差法，并采用混合训练期和60 d滑动训练期方案，对2012年欧洲中期天气预报中心(E-CMWF)数值模式夏季1—5 d的降水预报进行订正试验，试验结果表明：该方法对25 mm·d^-1以上降水预报的ETS评分提高明显。

从上述国内外研究现状可以看到，基于数值预报产品释用的客观定量降水预报方法是目前短期客观定量预报的重要发展方向。然而，目前国内外还没有一种成熟、有效并得到预报业务技术人员普遍认同的数值预报产品释用方法。为此，本文尝试提出一种基于最大相关最小冗余度(maximum relevance minimum redundancy，MRMR)算法和随机森林回归算法(random forest regression，RFR)相结合的数值预报产品释用方法。该方法以广西89个站点的降水作为预报对象，采用最大相关最小冗余度算法从ECMWF集合预报的51个成员筛选出若干个成员，这些成员与预报对象具有相关性最大、成员间冗余度最小的优点。将提取出的预报成员作为随机森林回归算法的输入因子进行释用。

1.   资料与方法

1.1   资料

本文所用资料是ECMWF集合预报模式的逐日08:00(北京时，下同)和20:00的24 h和48 h降水量预报产品(格距为0.5°×0.5°，水平范围为21°~26.25°N，104.75°~112°E)；选取试验样本时段为2012年4月—2016年9月，共2176个样本(去除数据缺失样本)，其中2012年4月—2015年12月为建模样本，2016年1—9月(共526个样本)为独立样本进行预报试验。

1.2   方法

统计预报中，有两方面因素影响其预报效果：一是预报因子的选取，即所选取的模型预报因子与预报对象之间的相关性强弱以及预报因子间的复共线性大小问题；二是预报模型本身的问题，即所选取模型的拟合能力和泛化性能强弱。为此，本文首先采用最大相关最小冗余度(MRMR)算法，从多个预报因子中选取与预报对象相关性大而它们之间的预报信息重叠少(复共线性小)的若干个预报因子作为预报模型的输入。其次，在预报模型的构建上，考虑采用具有泛化能力和抗噪声能力强、训练时间短且不容易陷入过拟合的随机森林回归算法。

本文提出的基于最大相关最小冗余度算法和随机森林回归算法相结合的数值模式预报释用方法主要步骤如下：

① 采用三次多项式插值方法将每个ECMWF集合预报成员插值到站点上。

② 利用ECMWF集合预报成员插值的平均值，对预报对象的历史样本序列进行分类。

③ 将分类后的预报对象样本对应的因子矩阵采用最大相关和最小冗余度算法，对每个站点的ECMWF集合预报51个成员进行筛选。

④ 将步骤③选出的预报成员作为随机森林算法的模型输入进行预报建模，并输出预报结果。

1.2.1   最大相关最小冗余度算法

预报因子是统计预报模型的重要组成部分，预报因子的选取直接影响预报效果。本文将ECMWF集合预报产品51个成员的输出值作为预报模型的输入因子，成员之间的相关性和冗余性会影响预报模型的预报能力，为此本文尝试采用MRMR算法对集合预报的成员进行筛选。

MRMR算法是信息论中典型的特征降维算法^[14]，其主要思想是从特征空间中寻找与目标类别相关性最大而相互之间冗余性却最少的m个特征^[15]，他们之间的相关性和冗余性利用互信息^[16]衡量。互信息用于衡量两个随机变量之间的相互约束程度。对于给定S={x_i|i=1, …, m}为特征集合和目标类别Y，则特征集合S中的特征与目标类别Y的最大相关度以及度量特征集合S中各特征间的互相关度分别为式(1)和式(2)。

Peng等^[15]定义算子(D, R)结合相关性因子D和冗余性因子R，即将式(1)减式(2)进行组合，从而得到最大相关性最小冗余度准则(MRMR)：

1.2.2   随机森林回归算法

统计预报模型构建方法的选取，对释用效果也有直接影响，本文采用多决策树组合而成的随机森林回归算法(random forest regression，RFR)进行构建预报模型，RFR算法是由Breiman于2001年提出的一种非线性统计方法^[17]，该算法具有抗噪声强、预报结果稳定的优点。RFR是利用自举法从原始样本中抽取多个训练样本子集，对每个样本子集分别进行决策树建模。进一步通过组合多棵决策树进行预测，并通过取平均值得到最终预测结果^[18]，其本质与气象上的集合预报思想相近，是将多棵决策树建模得到的预测结果进行集成。

随机森林回归算法同时还具有计算速度快、泛化性能好以及参数少等优点，且不容易出现人工神经网络的过拟合现象。目前该方法在农业、水文和医学等众多领域得到广泛应用^[19-21]。然而，该方法在降水预测中的应用报道较为少见。为此，本文尝试采用该方法进行数值模式的降水预报产品进行释用。

RFR算法是通过自举法抽样技术，由随机向量θ_i(即回归决策树)生长形成{h(X, θ_i), i=1, …, k}的组合模型。预测变量是一个数值型变量，与其分类模型不同，其预测值是通过k棵回归决策树的预测结果取平均值得到的。RFR算法实现见图 1。

图  1  RFR算法流程图

Figure  1.  RFR algorithm flow chart

下载: 全尺寸图片 幻灯片

2.   建模试验

本文以广西89个气象站未来24 h(每日2个预报时次为08:00和20:00)降水量作为预报对象，重点研究暴雨以上量级的降水预报。研究区域和站点分布如图 2所示。

图  2  研究区域和站点分布

Figure  2.  The target area and station distribution

下载: 全尺寸图片 幻灯片

2.1   单站暴雨集合预报建模样本及预报因子处理

经过统计分析，2012年4月—2016年9月共2176个样本中(去除资料缺失样本)，广西89个基本气象站24 h(统计时段包含当日08:00—次日08:00和前日20:00—当日20:00)内出现1个站以上、降水量达到暴雨以上的有4775个站次。本文的研究对象24 h降水量达到暴雨以上量级的情况比较频繁。为提高暴雨预报的准确率，在进行单站暴雨集合预报建模时，对模型的建模样本采取分类处理方式，具体建模步骤如下：

① 对ECMWF集合预报的每个成员(共51个成员)，利用前一日48 h累积降水量预报场减去24 h累积降水量预报场，从而获得该成员当日的未来24 h降水量预报场R₂₄。

② 采用多项式插值方法，将降水量预报场R₂₄插值到89个气象站。此时，每个预报对象(气象站)都对应有51个预报因子(51个ECMWF集合预报成员的插值)，记为F₅₁。

③ 对于第k(k=1, …, 89)个预报对象Y(气象站), 如果该站点的F₅₁平均值大于α mm, 则在预报对象Y的历史样本序列(其样本序号记为Ω)中，选出降水量大于β mm的样本序号集Ω₁(Ω₁⊂Ω)。若该站点的F₅₁平均值小于α mm，则令Ω₁=Ω。

④ 利用步骤③得到的预报对象Y的建模样本序号集Ω₁，可得到与之对应的由F₅₁组成的因子矩阵S'和预报量序列Y'。

⑤ 将S'和Y'带入MRMR算法，求出这51个因子中(集合预报成员)与预报量Y'相关性最大和冗余度最小的n个因子(成员)。

⑥ 以步骤⑤选出的S'和Y'为基础，采用随机森林回归模型算法建立集合预报模型。

⑦ 将步骤⑤选出的因子(成员)对应的预报样本代入步骤⑥训练好的预报模型中，得到对该站未来24 h的降水量预报。

2.2   结果分析

在进行建模试验时，为了能够更好地了解本文提出方法的释用能力，在此采取交叉检验法，即将2012年4月—2015年12月的样本作为独立样本进行与实际预报相同的独立预报试验，对2016年1—9月进行业务预报试验。在进行试验计算前，要先确定出2.1节计算步骤中的几个阈值参数(n，α，β)。参数n用于选取若干个集合预报成员作为模型的预报因子参与最后的训练计算，根据文献[22]的研究结果，参与集合预报计算的成员，其数量控制在8~10个比较合适。本文先设定参数n为10，分析另外两个阈值参数α和β不同取值时对预报模型的影响情况。同时，为了能够将预报结果与数值预报产品进行比较，本文采用常用的多项式插值方法将ECMWF集合预报成员的格点资料插值到站点上，并将所有成员插值到站点后取其平均值的预报方法称为AVI(average value of the interpolation)。表 1为本文提出的新方法(称为MRMR-RFR)预报结果与相应预报时段的ECMWF集合预报51个成员插值到站点后取其平均值(AVI方法)的TS评分，表 1中TS评分越高，表示对应方法预报的准确性越高。考虑到ETS评分方法可有效去除随机降水概率对评分的影响，同时也可以避免使用气候概率的情况，本文也对试验预报结果采用ETS评分进行对比。

表  1  2012—2015年暴雨以上量级降水交叉独立预报TS评分

Table  1.  TS of cross independent sample test forecast of rainstorm from 2012 to 2015

年份	MRMR-RFR									AVI
	(n=10, α=15, β=10)			(n=10, α=20, β=15)			(n=10, α=25, β=15)
	TS	ETS		TS	ETS		TS	ETS		TS	ETS
2012	0.12	0.11		0.15	0.12		0.14	0.11		0.07	0.06
2013	0.11	0.10		0.13	0.12		0.13	0.12		0.09	0.08
2014	0.11	0.11		0.17	0.14		0.13	0.12		0.08	0.07
2015	0.12	0.11		0.14	0.12		0.14	0.10		0.06	0.06

下载: 导出CSV 

| 显示表格

对数值模式产品的解释应用，是在承认其具有一定预报能力的基础上进行的。因此，假设当集合预报的成员插值平均达到15 mm以上(阈值α取15 mm以上)时，实况降水量才可能出现暴雨以上降水的情况。同时，由于模式预报结果与实际降水量差异存在，一般设定阈值参数β＜α。从表 1的统计结果中可以看到：两种评分方法的统计结果中，ETS评分整体上略低于TS评分，这是因为ETS评分方法对空报和漏报都有惩罚。相比AVI方法的ETS评分较TS评分偏低幅度，MRMR-RFR释用方法ETS评分偏低幅度更大一些，说明MRMR-RFR释用方法空报的次数多于AVI方法。②本文选取的3组阈值参数，其预报结果在暴雨以上量级降水预报(降水量不小于50 mm)的TS和ETS评分相差不大，说明这两个阈值参数对其预报能力的敏感性不明显，只要取值在一定范围之内，MRMR-RFR释用方法TS和ETS评分均高于AVI方法的预报结果，这其中提高幅度最大的是α取20 mm，β取15 mm时，即当集合预报51个成员的插值平均预报降水量达到20 mm以上时，选取实况降水量大于15 mm的样本进行建模，该组参数的预报结果在2012年4月—2015年12月的交叉独立预报中，其暴雨TS评分比AVI方法分别提高了0.08，0.04，0.09和0.08，而相应的ETS评分方法也比AVI方法的ETS评分分别提高了0.06，0.04，0.07和0.06。

由上述的分析可知，参数α和β在一定范围内进行取值，MRMR-RFR方法的预报效果稳定。为了进一步考察参数n取值对预报模型稳定性的影响，首先固定参数α和β的取值(取TS，ETS评分最高的参数组合：α=20 mm，β=15 mm)选取不同的n值进行上述同样本的试验分析，其中n的取值分别尝试取9个和11个(8~10个的附近^[22])，其统计结果见表 2。

表  2  不同参数n下2012—2015年暴雨交叉独立预报TS评分

Table  2.  TS of cross independent sample test forecast of rainstorm under different n from 2012 to 2015

年份	MRMR-RFR
	(n=9, α=20, β=15)			(n=10, α=20, β=15)			(n=11, α=20, β=15)
	TS	ETS		TS	ETS		TS	ETS
2012	0.15	0.13		0.15	0.12		0.14	0.12
2013	0.13	0.11		0.13	0.12		0.14	0.13
2014	0.14	0.12		0.17	0.14		0.15	0.13
2015	0.12	0.11		0.14	0.12		0.16	0.14

下载: 导出CSV 

| 显示表格

分析表 2可知，MRMR-RFR方法在3个不同n取值情况下，该释用方法交叉独立样本的试验结果中，3组参数组合的TS，ETS评分互有高低，统计结果相对比较稳定。由此可知，参数n在一定范围的取值，MRMR-RFR释用方法的预报结果也保持稳定。

为了进一步检验模型的预报性能，采用2012年4月—2015年12月试验结果中TS，ETS评分相对较高的参数组合(n=10, α=20, β=15)，对2016年1—9月进行业务预报试验。预报结果见表 3。

表  3  2016年1—9月单站暴雨以上量级降水业务预报TS，ETS评分

Table  3.  TS and ETS of single-station forecast of rainstorm using different methods from Jan 2016 to Sep 2016

月份	MRMR-RFR (n=10，α=20, β=15)			AVI		降水量不小于50 mm站次
	TS	ETS		TS	ETS
1	0.24	0.18		0.00	0.00	60
2	0.50	0.45		0.00	0.00	1
3	0.07	0.05		0.06	0.05	11
4	0.08	0.05		0.03	0.02	105
5	0.17	0.15		0.13	0.11	203
6	0.17	0.14		0.06	0.03	231
7	0.13	0.10		0.01	0.00	114
8	0.24	0.22		0.18	0.15	206
9	0.07	0.04		0.00	0.00	36
1—9	0.15	0.11		0.08	0.06	967

下载: 导出CSV 

| 显示表格

由表 3可知，本文提出的MRMR-RFR释用方法在的逐次独立样本预报检验中，9个月TS，ETS评分全部超过0分，并且在出现暴雨站次达到100个以上的4—8月，该释用方法的TS，ETS评分均明显高于AVI方法。统计可知，MRMR-RFR释用方法在5—8月TS，ETS评分均高于0.10。该释用方法1—9月平均的TS，ETS评分比AVI方法评分分别提高0.07和0.05。

本文选取2016年8月两次暴雨预报情况进行分析。一次过程选取受2016年第4号台风妮妲影响而造成的大范围暴雨，具体时段为2016年8月2日20:00—3日20:00，共有30个气象站降水量达到50 mm以上，图 3为两种方法预报对比。另一次过程选取2016年8月12日20:00—13日20:00的非台风类的一般性暴雨强降水过程，共有27个气象站降水量达到50 mm以上，图 4为该次过程两种方法预报对比。

图  3  2016年8月2日20：00—3日20:00 24 h降水实况与预报

(a)实况，(b)MRMR-RFR释用方法预报，(c)AVI方法预报

Figure  3.  Observation and prediction of the case from 2000 BT 2 Aug to 2000 BT 3 Aug in 2016

(a)observation, (b)prediction of MRMR-RFR, (c)prediction of AVI

下载: 全尺寸图片 幻灯片

图  4  2016年8月12日20:00—13日20:00 24 h降水实况和预报

(a)实况，(b)MRMR-RFR释用方法预报，(c)AVI方法预报

Figure  4.  Observation and prediction of the case from 2000 BT 12 Aug to 2000 BT 13 Aug in 2016

(a)observation, (b)prediction of MRMR-RFR, (c)prediction of AVI

下载: 全尺寸图片 幻灯片

由图 3可以看到，对于8月2日20:00—3日20:00预报中，MRMR-RFR释用方法和AVI方法预报降水量大于50 mm的落区基本重合，覆盖了全区的大部分地区，然而，MRMR-RFR释用方法比AVI方法在广西西南部和东南部这两个区域更接近实况。

由图 4可知，对于8月12日20:00—13日20:00 AVI方法预报出该降水过程的雨带，但所有站降水量均未超过50 mm。而实况有27个站降水量超过50 mm。MRMR-RFR释用方法预报出13个站降水量超过50 mm，空报16个站，漏报14个站，TS评分为0.30(ETS评分为0.16)。即MRMR-RFR释用方法具有正预报技巧。

综上所述，MRMR-RFR释用方法对ECMWF集合预报产品释用，其预报技巧在大部分情况下为正技巧，可以一定程度上提高ECMWF集合预报的暴雨预报能力。这与本文采用的因子选取方法、建模样本分类以及预报模型建立的方法有关，通过这些处理，既可以使模型更加专注于强的降水过程，同时又能在预报强降水的建模样本中消除部分降水量很小或较小的样本，减少建模样本中的噪声(干扰)。

3.   小结

1) 采用最大相关最小冗余度的互信息技术进行因子选取，可以提高入选因子的预报信息，同时入选因子与预报量具有最大的相关性。构建模型时采用可调参数极少、计算速度快且有很好的非线性拟合能力和泛化性能的随机森林算法。

2) 独立样本预报试验表明：本文提出的MRMR-RFR释用方法相对于ECMWF集合预报产品插值方法(AVI方法)，能更好地判断暴雨的落区及落区范围大小。

本文提出的对数值预报产品的MRMR-RFR释用方法较为简单，便于预报人员在业务中使用。