论文快递｜黄强等：基于可解释机器学习的黄河源区径流分期组合预报

黄河源区年径流量占全流域的34.1%，是黄河重要的产流区和水源涵养区，素有“黄河水塔”之称［1－2］。黄河流域属于资源性缺水地区，源区来水的变化对全流域水资源开发、利用和管理影响巨大，加强黄河源区径流预报研究至关重要。此外，全球升温对黄河源区在内的冰冻圈、大气圈、水圈、生物圈等多个圈层产生了显著影响［3］，绿色低碳发展已经成为实现温控目标的全球共识。黄河源区跨越我国一、二级阶梯，地势落差大，水能资源蕴藏量达2 000 万kW。在水能资源富集区规划建设水风光一体化清洁能源基地，是黄河等我国主要流域绿色低碳转型的重要举措［4］。流域水电站是改善风、光出力随机波动性的调节中枢，径流则是决定其调节能力的主控因素，准确预报径流情势可为流域清洁能源高效利用、绿色低碳转型提供关键技术支撑。因此，开展黄河源区径流预报研究对于流域水资源科学调配和绿色低碳发展具有重要意义。

近年来，基于机器学习的数据驱动模型在径流预报领域快速发展，并得到广泛应用［5］。Lu 等［6］将XGBoost（Extreme Gradient Boosting） 模型、RF（Random Forest）模型、AdaBoost（Adaptive Boosting）模型应用于富春江水库的日出库流量预测，结果表明XGBoost 模型预测精度最高且预报能力比较稳定。Han 等［7］将深度学习方法作为预测模型误差后处理器，校正了俄罗斯河径流预报结果，有效提升了预测精度，并有助于揭示多种不确定因素对模型预测结果的影响。

虽然机器学习善于模拟高维非线性复杂系统，在径流预报研究中表现良好，但是其结构复杂且具有“黑箱”属性，无法真正理解和刻画水文过程及其物理规律，导致预测结果难以被充分信任［8］。机器学习的可解释性分析旨在揭示模型输入因子对预报对象的驱动机制，提高机器学习模型决策过程的透明度。Lundberg 等［9］将博弈论与局部解释性关联起来，提出了机器学习模型的SHAP 可解释性分析框架，阐明了各输入变量对机器学习决策的影响机理。Liu 等［10］构建了一个月预见期的XGBoost 径流预测模型，并使用SHAP 可解释性分析方法计算了SHAP 总效应值、主效应值、交互作用值和损失值，量化了XGBoost 径流预报模型中气象要素、历史径流、全球气候因子等输入变量对径流预报结果的贡献率。Jing 等［11］将LSTM 模型应用于汉江上游日径流预测中，利用IG（Integrated Gradient）可解释性方法研究了径流预测的潜在决策机制。

基于上述研究，笔者以黄河源区唐乃亥和玛曲水文站为研究对象，考虑季节性融雪变化影响，建立中长期径流分期组合机器学习预报模型，构建机器学习预报模型的可解释性分析框架，定量识别不同预报因子的贡献率及多因子的复杂交互作用。

黄河发源于青藏高原的巴颜喀拉山北麓，全长5 464 km，流域面积79.5 万km2。其中，唐乃亥站以上为黄河源区（见图1）。源区河道总长1 553 km，面积达12.19 万km2，占黄河流域总面积的16.2%。唐乃亥站多年平均来水量200 亿m3，年径流量占全流域的34.1%［2］，源区多年平均径流量分割统计显示，降雨、地下水补给（基流）、冰雪及冻土融水分别占年总径流量的63.15%、26.18%和9.17%［12］。黄河源区湖泊数量众多，水资源丰沛。河段落差达3 745 m，水能资源丰富［13－14］，是我国九大清洁能源基地之一。

图１　黄河源区

 

研究数据涵盖了水文、气象、积雪覆盖、大尺度环流因子等多源数据，时段跨度为1960—2020 年。具体的数据类型及其来源见表1。

表１　主要数据及其来源

 

基于积雪覆盖变化和弹性系数精细划分径流的年内预报时段，构建基于XGBoost 的分期组合预报模型，并使用SHAP 可解释性分析方法研究不同预报因子对径流的驱动机制。技术路线见图2。

图２　技术路线

 

根据文献［15］提出的方法，推导降水、融雪水当量弹性系数，定量表征径流对降水、融雪水当量变化的敏感程度。 计算公式如下：

式中：φ 为干旱指数；ET0 为潜在蒸散发量；P 为降水量；εP为降水弹性系数；εSW为融雪水当量弹性系数；r为融雪径流占总径流的比例；n 为下垫面参数，主要与流域下垫面情况有关，参考以往研究［16］，n 取2.5。

黄河源区年内不同时期的融雪径流占比差异较大，为提升中长期径流预报的精度，本研究以融雪水当量弹性系数和流域积雪覆盖率变化为判断依据，将年内的径流预报时段划分为融雪影响期和非融雪主导期，分别构建径流预报模型。 融雪影响期的划分标准为：1）融雪水当量弹性系数大于0；2）时段内流域积雪覆盖率呈下降趋势；3）流域内积雪覆盖率不小于年内最大积雪覆盖率（其中，唐乃亥以上流域为32%、玛曲以上流域为38%）的10%；4）融雪径流占比达到10%。参考文献［17］，采用下式估算融雪径流比。

式中：SWr为融雪径流占比，SW 为融雪水当量。

采用分期组合预报策略，分别对融雪影响期和非融雪主导期开展径流预报，再将两个时期的预报结果组合，得到年内所有月份的预报结果。 分期组合预报策略见图3。

图３　分期组合预报策略示意

 

备选的预报因子分为3 类，包括水文气象因子（Ⅰ类，历史月平均流量、月累计降水量等8 种）、大尺度环流因子（Ⅱ类，8 种大尺度环流因子与太阳黑子指数）和融雪水当量（Ⅲ类）。 针对融雪影响期和非融雪主导期，采用点间互信息法和相关系数法依次筛选预报因子的最佳组合及最优滞时。 需要注意的是，非融雪主导期的融雪径流占比低，融雪水当量可不作为该期的主要预报因子。

从“预报因子类型”和“是否进行预报因子筛选”两个角度出发，设置多种组合方案，在融雪影响期和非融雪主导期分别优选出性能最好的预报模型，具体的方案设置见表2。

表２　融雪影响期和非融雪主导期预报模型的优选方案

 

基于优选出的融雪影响期和非融雪主导期预报模型，建立中长期径流分期组合预报模型。 此外，为评估分期组合预报策略在研究区的适用性，对比分期组合模型和传统不分期模型的预测性能，对比方案设置见表3。

表３　分期组合模型和传统不分期模型的预报性能对比方案

 

上述方案的径流预报模型均基于XGBoost 模型构建，通过粒子群算法迭代寻优XGBoost 模型的超参数。模型的率定期设定为1960—1999 年，验证期为2000—2014 年，测试期为2015—2020 年。

为了深入挖掘预报模型的决策机制并提升可解释性，基于XGBoost 模型构建径流分期组合预报模型，并采用SHAP 方法研究其可解释性分析框架。

通过计算SHAP 值来衡量单变量预报因子、多因子组合对预报结果的贡献程度，由此可开展预报结果的归因分析［18］。 对于一个已经训练好的机器学习预报模型，各预报因子的SHAP 值可由下式计算：

式中：βi为预报因子i 的SHAP 值，M 为所有预报因子的集合，为所有预报因子的个数，S 为不包含因子i 的任意预报因子子集，为子集中预报因子的个数，S ∪｛i｝ 为包含因子i 的任意预报因子子集，fS∪｛i｝（xS∪｛i｝）为由预报因子子集S∪｛i｝得到的预报结果，fS（xS）为由预报因子子集S 得到的预报结果。

此外，SHAP 机器学习可解释性分析框架还可以分解单一预报因子的SHAP 主效应值和两种预报因子交互作用的SHAP 交互作用值。

式中：βi，j为预报因子i 和j 的SHAP 交互作用值，T 为不包含因子i 和j 的任意预报因子子集，T∪｛i，j｝为包含因子i 和j 的任意预报因子子集，fT∪｛i，j｝（xT∪｛i，j｝）为由预报因子子集T ∪｛i，j｝得到的模型预报结果，fT（xT）为由预报因子子集T 得到的模型预报结果，βi，i为预报因子i 的SHAP 主效应值。

根据SHAP 值、SHAP 主效应值和SHAP 交互作用值，从各预报因子对预报结果的贡献程度、预报结果对预报因子的依赖关系、预报因子间的交互作用等角度对建立的径流预报模型进行可解释性研究。

基于水量平衡原理，计算了唐乃亥站和玛曲站以上流域不同月份降水弹性系数和融雪水当量弹性系数（见表4）。结果显示两者均为正值，说明降水、融雪水当量与径流正相关，弹性系数计算结果具有合理性。融雪水当量弹性系数在唐乃亥站以上流域和玛曲站以上流域的年内变化过程相似，于4—5 月达最大值，此时径流对融雪变化最敏感。如图4 所示，融雪径流占比从3 月起显著上升，在4—5 月达到最大（40%～50%）。由此可见，融雪水是黄河源区地表径流的重要组成部分［19－20］。

表4　唐乃亥站和玛曲站以上控制流域降水和融雪水当量弹性系数

 

图４　融雪径流占比

 

如图5 所示，根据2.1 节所述的年内预报时段划分方法，识别出融雪影响期为3—6 月，非融雪主导期为7 月—次年2 月，唐乃亥站和玛曲站以上流域的预报时段划分结果一致。

图５　唐乃亥站和玛曲站以上流域的预报时段划分

 

3.2.1 融雪影响期和非融雪主导期预报模型优选

根据表2 设定的预报模型方案集，分别在融雪影响期和非融雪主导期优选出预报性能最佳的模型方案。

首先，利用点间互信息法初筛各方案（见表2）的径流预报因子。图6 显示了各预报因子的点间互信息（PMI，Pointwise Mutual Information）值。在融雪影响期和非融雪主导期，PMI 值较大的前一个月流量（Qt－1）、气温、蒸发等气象要素与唐乃亥站、玛曲站流量表现出显著相关性。此外，在融雪影响期，融雪水当量的PMI 值均超过0.4（排名第6），也是重要的预报因子。

图６　不同类型预报因子的ＰＭＩ 值

 

分别取两时期内PMI 值较大的10 种因子作为对应径流预报时期的预报因子。采用相关系数法确定各方案初筛得到的预报因子的最优滞时，不同滞时预报因子与径流的相关系数如图7 所示。由图7 可见，在融雪影响期和非融雪主导期，滞时为一个月的前期流量、前期降水与当月流量之间的相关系数均为1～12 个月中的最大值（大于0.68），关系密切。在融雪影响期，径流与前一个月融雪水当量的相关系数达0.75。根据不同滞时预报因子与流量的相关系数大小，确定了不同模型方案预报因子的最优滞时，见表5，其中，Qt－1表示滞时为1 个月的前期流量，其余因子采用类似的表达方式。

表５　不同方案下预报因子组合及其最优滞时

 

图７　不同滞时（前期１～ １２ 个月）预报因子与径流的相关系数

 

根据表5 中的不同方案预报因子集合，采用粒子群优化算法优化XGBoost 的超参数，分别率定融雪影响期和非融雪主导期的径流预报模型。 并基于测试期的纳什效率系数（NSE）、确定系数（R2）、均方根误差与标准偏差之比（RSR）和均方根误差（RMSE），分别筛选融雪影响期和非融雪主导期最优预测模型，结果见表6。

表６　各方案下融雪影响期和非融雪主导期模型的预报性能（测试期）

 

由表6 可以发现：在融雪影响期，以水文气象要素和融雪水当量作为预报因子时（TNH－Sw－3－S 和MQSw－3－S 方案），模型预报能力均达到最优；在非融雪主导期，将水文气象要素作为预报因子时（TNH－Pr－1－S和MQ－Pr－1－S 方案），其模型预报能力在两站均为最佳。

通过筛选预报因子并以最优滞时作为输入（对应表6 中加星号的方案），模型预报精度可得到提升。例如，在唐乃亥站，经预报因子筛选后的模型方案TNH－Sw－3－S 和TNH－Pr－1－S 与未筛选方案TNH－Sw－3 和TNH－Pr－1 相比，NSE 分别提高约17%和6%、RSR 分别降低约44%和17%、RMSE 分别减少约44%和17%。

由表6 还可知：TNH－Sw－3－S、TNH－Pr－1－S 和MQ－Sw－3－S模型方案表现较佳；在玛曲站非融雪主导期内，虽然MQ－Pr－1－S 的R2略低于MQ－Pr－2，但其余3 个指标均更优，可视为最优模型。据此，得到研究区中长期分期组合预报模型方案如下：由方案TNH－Sw－3－S 和TNH－Pr－1－S 可得到唐乃亥站性能最优的分期组合模型（记为TNH－F－0－S）；由方案MQ－Sw－3－S和MQ－Pr－1－S 得到玛曲站性能最佳的分期组合模型（记为MQ－F－0－S）。

3.2.2 分期组合预报模型性能分析

将得到的最优分期组合模型与传统不分期模型（建模方案见表3）进行对比，评估分期组合预报策略在研究区的适用性及其预报性能提升效果。

在构建不分期模型时，运用点间互信息和相关系数方法筛选各模型方案（表3 中Z－1－S、Z－2－S、Z－3－S 和Z－4－S 方案）中预报因子类型及其最优滞时，进而采用PSO 优化的XGBoost 机器学习算法构建预测模型，4 种方案不分期模型性能见表7。由表7 可见，与不分期模型相比，分期预报模型在唐乃亥站和玛曲站的NSE 分别提高3%和10%、RSR 分别降低10%和17%、RMSE 分别减少10%和17%，两站径流预报的准确率得到提高。综上可知，本研究构建的分期组合预报模型（即表7 中TNH－F－0－S 和MQ－F－0－S）的预报精度高，整体性能较传统不分期模型明显提升，分期组合预报策略在研究区具有良好适用性。

表７　分期组合模型与不分期模型的性能对比（测试期）

 

为进一步减小径流预报误差，提升分期组合预报模型的精度，采用分位数映射方法校正径流预报误差，结果见表8。根据《水文情报预报规范》（GB／T 22482—2008）对预报精度的分级标准，在校正后的TNH－F－0－S分期组合模型的确定系数R2为0.926，达到甲等精度；MQ－F－0－S 分期组合模型的确定系数R2为0.850，接近甲等精度。

表８　误差校正前后最优预报模型在测试期的性能指标

 

分期组合模型预报功能的实现， 需要利用XGBoost 机器学习方法拟合预报因子与径流的复杂响应关系。为了深入理解XGBoost 在径流预报中的决策逻辑和工作机制，提升分期组合预报模型的可理解性和透明度，基于上文优选的最优组合预报模型及其预测结果，计算各预报因子的SHAP 值、SHAP 主效应值和SHAP 交互效应值，依次解析组合预报模型的单因子贡献度、因子依赖关系和多因子间交互作用。

3.3.1 单因子贡献度

SHAP 值可以衡量某一模型预报因子对预报结果的贡献程度，值越大表示其对预报结果的贡献越大。计算各预报因子的SHAP 均值，并按其绝对值大小进行排序，如图8 所示。由图8 可知，在不同水文站和不同预报时段，尽管预报因子对预报结果的贡献度略有差异，但排名前列的预报因子类型相同。综合SHAP 均值来看，预测因子按贡献由大到小依次为降水、前一个月流量、蒸发、气温、相对湿度、融雪水当量、气压和风速。

图８　不同预报时期预报因子的ＳＨＡＰ 值

 

采用SHAP 方法进行预报模型可解释性分析的优势还在于能够可视化模型的决策路径。图9 中每条折线代表预报模型对样本的决策路径，颜色表示预报流量的大小（蓝色表示小流量预报值，红色表示大流量预报值）。线条方向的变化指示预报因子的影响效果，向左（SHAP 值＜0） 表示会减小预报值，向右（SHAP 值＞0）则增大预报值。结果表明，在小流量模型中，预报因子的影响较一致，归纳其典型决策路径为：Q（减小）、PRE（减小）、EVP（减小）、TEM（增大）、RHU（增大）、WIN（增大）、PRS（减小）、SW（减小）；但随着预报流量的增大，其预报因子的决策路径变得更加复杂。

图９　模型样本的决策路径

 

3.3.2 预报值对因子的依赖关系

绘制SHAP 部分依赖图，可以识别预报结果与预报因子间线性或非线性、单调或非单调的复杂依赖关系，如图10 所示。图10 中，横、纵轴分别代表预报因子数值与SHAP 主效应值。可以发现：流量预报值与前期降水整体成单调递增的依赖关系，与蒸发则成单调递减的依赖关系，且在较大范围内以线性依赖关系为主。除单调线性关系外，预报径流与风速还存在非单调、非线性的复杂依赖关系。此外，预报流量值与前期流量、相对湿度整体同样成单调递增的关系，与气温、气压成非单调、非线性的复杂依赖关系。受篇幅限制，此处不再赘述。

图１０　ＳＨＡＰ 部分依赖图

 

对于成单调依赖关系的预报因子，其依赖关系在各径流预报时段基本相似，但变化阈值存在一定差异。以降水为例，构建的预报模型在不同时期均显示出了单调递增的线性依赖关系，在同一预报时期内预报因子阈值接近，推测与研究区高寒山区下垫面对径流调节机制的影响有关［21］。例如，在融雪影响期，气温较低且冰雪、冻土未完全消融，土壤下渗能力弱，少量降水即可产生径流，在局部依赖图中表现为SHAP 主效应值开始增长的阈值较小但之后的响应斜率较大；而在非融雪主导期，气温和蒸发强度较大，土壤下渗能力增强，需要更加充足的降水才能产生地表径流，因此相应的阈值也较大。

3.3.3 两因子的交互作用

计算预报因子两两组合的SHAP 交互作用值，揭示多因子交互作用对预报流量的影响。图11 为降水与其他预报因子的SHAP 交互作用散点图（以唐乃亥站非融雪主导期TNH－Pr－1－S 模型、玛曲站非融雪主导期MQ－Pr－1－S 模型为例）。图中，纵轴表示SHAP交互作用值，散点颜色表示预报因子（降水）值，横轴是与降水组合的预报因子。SHAP 交互作用值的绝对值越大，交互作用越强；交互作用值为正表示有增大径流量的效果，负值则相反。

图１１　降水与其他预报因子的ＳＨＡＰ 交互作用

 

以SHAP 交互作用绝对值判断交互作用强度，发现降水与气温、降水与蒸发、降水与前期流量的交互作用较强。其中，在唐乃亥站非融雪主导期TNH－Pr－1－S 模型中，前期降水与气温交互作用值达－100，表明在低气温下的降水以积雪为主，导致径流预报值容易出现减小的情况。在玛曲站非融雪主导期MQ－Pr－1－S模型中，前期降水与蒸发交互作用值约为70，表明当降水量较大、蒸发较少时，会导致径流量增大。

同时还发现SHAP 交互作用散点分布具有拖尾式或阶跃式的特征。拖尾式的特点为在某一区间内交互作用变化不显著（SHAP 交互作用值基本保持不变），区间外连续变化显著（ 散点连续分布）， 如TNH－Pr－1－S模型中降水与蒸发、降水与风速的组合。阶跃式的特点为在某一区间内散点分层（分区）现象明显，反映出当给定某一降水时交互作用变化对横坐标预报因子的变化不敏感，此时交互作用强度主要由降水量决定，如MQ－Pr－1－S 模型中降水与流量、降水与气压组合等。

以黄河源区唐乃亥站和玛曲站为研究对象，提出了分期组合预报策略，构建了中长期径流分期组合预报机器学习模型及其可解释性分析框架，主要结论如下。

1）黄河源区存在显著的季节性积融雪过程，基于水量平衡原理并结合积雪覆盖率、融雪水当量变化，定量划分年内的径流预报时段，其中融雪影响期为3 月至6 月，非融雪主导期为7 月至次年2 月。

2）与传统模型相比，采用分期组合预报模型能提高黄河上游唐乃亥站和玛曲站径流预报精度。与不分期模型相比，分期组合预报模型在唐乃亥站和玛曲站的NSE 分别提高3%和10%、RSR 分别降低10%和17%、RMSE 分别降低10%和17%，经误差校正后唐乃亥站预报模型R2 为0.926，玛曲站预报模型R2 为0.850。

3）基于SHAP 机器学习可解释性分析框架，识别出预报因子对径流预报结果的贡献程度从大到小依次为降水、前一个月流量、蒸发、气温、相对湿度、融雪水当量等。径流对不同预报因子的依赖关系复杂，其中降水、前期流量和相对湿度以单调递增的依赖关系为主，蒸发在一定范围内成单调递减的线性依赖关系，不同预报因子之间交互作用具有拖尾式或阶跃式的特征。