月径流预报建模方法对比分析——以嘉陵江北碚站为例

陈 雪 怡，陈 元 芳，王 文 鹏，邱 鹏

(1.河海大学 水文水资源学院，江苏 南京 210098； 2.重庆市水文监测总站，重庆 401147)

0 引 言

准确的中长期水文预报对流域水旱灾害防御、水资源综合利用、河湖水环境保护等水利业务工作具有重要意义[1]。根据预报变量的不同，径流预报可以划分为两类[2]：一类是基于单变量的水文预报，即发掘水文要素自身演变规律，常见的方法有时间序列法[3-4]、周期叠加法等；另一类是基于多变量的水文预报，即发掘大气环流、海陆状况等成因变量与水文要素之间的联系，常见的方法有随机森林[5]、人工神经网络、支持向量机[6]等。近年来，这两类方法在中长期径流预报中都有广泛的应用[7]。Fan等[8]建立自回归求和滑动平均模型预测了塔里木河干流径流量，并分析了温度与降水对其的影响。Wang等[9]将经验模态分解法与时间序列法相结合，对黄河上游进行长期径流预测。林剑艺等[10]将支持向量机应用于中长期预报中，在建模过程中引入径向基核函数简化非线性问题的求解，并应用SCE-UA算法辨识支持向量机的参数。赵文秀等[11]采用随机森林模型筛选预报因子，将筛选后的因子作为RBF神经网络的输入层进行长期径流预报。综上，国内外学者对两类预报方法均进行了大量探索，但何种方法更优并没有统一的结论。实际上，在具体流域的预报工作中，仍有必要综合考虑预报要素的水文特性、预报因子的可获得性、预报方法的稳健性等因素，通过定量对比分析选取合适的方法建立预报模型，实现精度可靠的中长期水文预报。

因此，本文以嘉陵江流域北碚水文站为例，建立基于时间序列模型与成因驱动的径流预报模型；同时应用时变权重方法，集合单一模型进行组合预报，通过对比单变量与多变量、单一与组合模型预报结果，分析引入多变量与组合模型是否有利于提高径流预报容错能力、稳定性与精度，研究成果可为径流预报的建模策略提供参考。

1 研究方法

本文选择时间序列法作为单变量预报方法，随机森林、BP神经网络与多元线性回归作为多变量预报方法，并采用组合预报方法对单一模型的权重进行实时调整从而对月径流预报结果进行优化。方法流程如图1所示。

1.1 水文时间序列预报方法

时间序列分析方法最早由耶尔于1927年提出。水文时间序列预报就是依据水文过程在相邻时间截口之间时序相依的统计特征，建立由序列前期值驱动的水文状态量预报模型。对于平稳的时间序列，可以构建自回归模型(Autoregressive Model，AR)、滑动平均模型(Moving Average Model，MA)与自回归滑动平均模型(Auto-Regressive Moving Average Model，ARMA)[12]。自回归滑动平均模型，记为ARMA(p，q)，具有以下结构：

(1)

式中：xt为平稳时间序列；εt为白噪声序列；φp，θq为模型参数，常用尤尔-沃尔克法估计；对于中心化序列xt，φ0=0。

对于非平稳的水文序列，通过差分运算将原序列转化为平稳序列，再通过白噪声检验后与ARMA模型组合，得到自回归滑动平均求和模型(Autoregressive Integrated Moving Average Model，ARIMA)，记为ARIMA(p，d，q)。

1.2 成因驱动的多变量预报方法

基于成因驱动的多变量预报模型是从径流成因的角度构建能够反映驱动因子与预报要素之间关系的多元相关预报模型，应用此类方法时驱动因子的筛选与预报模型的选择至关重要。首先根据径流的物理成因选择88项大气环流指数、26项海温指数、16项其他气象指数和月降水量作为初选影响因子，计算前一水文年各因子与径流量的相关系数，剔除无物理背景、不显著相关的影响因子并进行排序。对排序前50的影响因子构建随机森林模型，计算未被袋装法选中的样本(OOB数据)打乱前后均方误差的增加量(IncMSE)来定量评价因子的重要性[13]，取IncMSE排位靠前的因子作为预报模型输入变量。本文对比分析的模型为随机森林、BP神经网络与多元线性回归。

随机森林(Random Forest，RF)最早由Breiman于21世纪初提出[14]。随机森林由多棵决策树构成，其基本思想为：从原始样本中随机多次抽取样本，对每个重抽样本都建立决策树模型，组合多棵决策树对应变量的估计结果，当应变量的综合估计结果在样本训练集和验证集表现俱佳时，就构成了随机森林模型[15-16]。

BP神经网络模型具有3层结构：输入层、隐含层和输出层。输入层中每个神经元对应一项驱动因子，隐含层神经元个数通过“试错法”调整，输出层神经元个数为预报对象个数。将正向误差反向传播至网络层，调整连接节点权重，不断重复直至输出结果误差函数达到最小[17-18]。

多元线性回归就是建立多个自变量与因变量之间的回归方程，一般以误差平方和为目标使用最小二乘法估计参数，数学表达式为

y=α0+α1x1+α2x2+…+αnxn+ε

(2)

式中：αn为回归系数；xn为驱动因子值；ε为模型误差。

1.3 时变权重组合预测

权重组合预测将“好而不同”的单一模型预报结果加权求和，以获得更准确的预报效果。根据权重是否随时间改变，可以分为定常权重与时变权重两类。定常权重是指保持各单一模型的权重不变进行组合，该类方法的研究较成熟，但在某些时段组合预测精度较单一模型更差。时变权重则随着不同模型性能的变化而进行权重调整，通过动态预报有效规避定常权重组合预测常介于最优与最劣模型间的问题[19]。本文选择时变权重方法中的时变误差平方和(TSSE)方法，根据预报前序时段的误差来确定局部预报时段的动态权重，公式如下：

(3)

(4)

1.4 预测性能评价指标

采用合格率(QR，预报值与实际值相对误差在±20%之间为合格)、平均绝对百分比误差(MAPE)、平均相对误差(MRE)、均方根误差(RMSE)和决定系数(R2)5项指标作为评价模型预报精度的指标。采用均方根误差变化量(ΔRMSE)作为评价模型稳定性的指标。具体计算公式见文献[19]。

2 实例应用

2.1 研究区域与数据

本文预报对象为嘉陵江北碚站汛期月径流量，其月径流量数据来源于《中华人民共和国水文年鉴》。计算武都、略阳、广元、万源、阆中、巴中、达县、遂宁、高坪区和沙坪坝10个气象站(见图2)的算术平均降水量作为代表降水量，数据来源于中国气象数据网。收集全球逐月百项气象指数，数据来源于国家气候中心。本文以1979～2006年为训练期，以2007～2019年为测试期。

由于嘉陵江中下游地区的大规模水土流失综合治理、水库水电站等大型水利设施的修建等人类活动对下垫面条件的影响，嘉陵江的实测径流序列发生了突变，不能代表天然径流状况。因此考虑方差修正的预置白方法和Mann-Kendall检验对北碚站训练期降水序列与径流序列进行分析与一致性修正[20-22]，发现北碚站汛期(6～8月)降水量没有明显突变，而径流量在1986年发生了突变，如图3所示。

绘制训练期径流量-降水量双累计曲线图，根据突变前双累计曲线回归方程对突变后径流进行修正，如图4所示。按照修正前后总径流量的比例对各月径流量同倍比缩放，生成还原后的天然月径流量序列。对修正后的月径流量序列再进行突变检验，未发现新的突变点，满足一致性条件。

2.2 单变量与多变量模型预报效果

以经过还原得到的北碚站6～8月天然径流量作为目标变量，构建预报模型。对于单变量预报模型，参数p、q根据d阶差分序列的自相关图和偏相关图初选，再结合模型AIC准则进行确定，分别为ARIMA(2，2，0)、ARIMA(6，2，0)和ARIMA(1，1，0)。对于多变量预报模型，对初选因子按照IncMSE指标进行排序，选择指标值明显大于其他因子的因子集作为预报因子，6月和8月选定前8项因子，7月选定前7项因子，作为模型的输入变量，其具体名称如表1所列。经过反复测试，增加剩余因子不足以明显改善解释能力，反而增加模型复杂度，降低模型的泛化能力。构建随机森林模型时，对随机森林模型的两个待选参数进行优化，即子决策树数量Ntree与回归树节点划分数Mtry。挑选使OOB数据泛化误差最小时的Mtry值作为最优参数，确定Mtry后代入使模型误差基本稳定的Ntree值作为最优参数。本文中各月随机森林模型Mtry分别取4，4和3，Ntree均取1 000。构建BP神经网络模型时，隐含层神经元个数用试错法比选，使得训练误差和测试误差最小，网络结构分别为8-4-1，7-5-1和8-5-1。

表1 北碚站月径流预报因子Tab.1 Monthly runoff forecast factors of Beibei Station

单变量与多变量预测模型的散点图如图5所示，4种模型的预测值都均匀分布在1∶1线的两侧，且在不同月份内同一模型都表现出相似的准确性与适应性。对于不同变量的预报模型，基于多变量的预报模型精度优于基于单变量的预报模型。相对而言，ARIMA预报径流量的点据分布更为分散，说明其拟合程度相比多变量预报模型较差。基于两种模式下的径流预报的MAPE箱形图如图6所示。由图可知，单变量ARIMA模型的MAPE值中位数约为40%，而3种多变量模型的MAPE值中位数均小于20%，说明多变量预测模型具有更好的预报精度。因此，从径流成因角度尽可能地挖掘潜在的气候驱动因子，增加预报模型的信息量，对提高径流预报精度是有利的。

综合对比3种多变量模型的预报性能，如表2所列。从预测精度上看，RF在训练期的合格率最高、MRE与RMSE最低，在不同时段下合格率达到了82%以上、MRE在13%以下；BP在测试期的合格率最高、MRE与RMSE最低，在不同时段下合格率达到了76%以上、MRE在18%以下。总体上RF与BP模型表现均优于多元线性回归模型。从数据离散程度上看，RF的R2最接近1，说明其拟合程度最优。从模型稳定性上看，BP的ΔRMSE最小，7月与8月测试期的RMSE小于训练期，指示该模型具有较强的泛化能力。总之，在多变量预测模型中RF和BP神经网络模型的预报性能优于传统的多元线性回归模型。

表2 不同模型预报精度Tab.2 Forecast accuracy of different models

2.3 组合模型预报效果

目前能够公开获取且满足时效性要求的径流成因要素还比较有限，相关预报模型的建模依据仍以小容量样本为主，模型泛化容易出现过拟合情况。根据统计学原理，集成单一模型的预报结果，可以减少模型泛化的方差，有效降低预报不确定性。因此，选择对序列趋势模拟较准确的3种多变量模型进行组合预报，根据前文权重计算方法得到权重取值如图7所示。总体而言，RF模型的预报误差较小，组合权重占优，但随着预报时刻推移，各模型的组合权重会沿程发生变化。多元线性回归模型的权重呈现出平稳的特征，而两种机器学习模型的权重在各月都表现出非平稳特征，且在训练期和测试期都有明显变化趋势。RF模型的权重呈现先增大后减小的趋势，说明其预报误差在训练期呈现减小趋势，到了测试期开始呈现增大趋势，而BP模型则与之相反，这也验证了RF模型泛化能力差而BP模型稳定性强的特点。

图8比较了3种多变量模型和组合模型的预报性能。经组合预报后，训练期合格率达到了89%以上、MRE在13%以下；测试期合格率达到了84%以上、MRE在15%以下。合格率在训练期与测试期均得到了明显提高，高于或接近于最优的单一模型；定量误差在7月测试期时劣于BP神经网络模型，在其他月份各时期都有明显下降趋势，低于或接近于单一最优模型。因此，时变权重组合预报方法能有效集合单一预测方法提供的信息，结合各模型的优点，进一步提高基于多变量的预测模型的精度。

3 结 论

本文以北碚站月径流量预报为例，对比不同建模方法的预报性能，得到以下结论。

(1) 从单变量与多变量模型选择上看：多变量预报法优于单变量预报法。多变量预报模型的MAPE中位数均小于20%，说明从径流成因角度挖掘潜在预报因子、增加样本信息量，是提高月径流预报精度的有效途径之一。

(2) 从不同多变量模型表现上看：RF与BP模型在定量误差与模型稳定性方面均表现较好，说明针对成因机制复杂，样本空间非线性可分的预报对象，非线性映射模型常优于线性回归模型。

(3) 从单一与组合模型上看：时变权重能有效集合各单一模型的优势，使组合后的预报精度在训练期与测试期都高于或接近于最优的单一模型，在有效减小预报误差的同时提高了模型的稳定性。

本文探讨了天然径流量预报的建模方法，但在实际预报业务中，仍需进一步从水保措施治理面积、大中型水利工程建设数量、用水量、水库蓄水量等因素，定量分析人类活动对月径流量的影响及其与时间的关系，将天然径流量的预报结果与人类活动影响导致的变化量相结合，得到月径流过程的实际变化，以提高预报模型的实用性。