基于数据驱动的陕北风沙区秃尾河流域径流变化归因

余煌浩，李彬权

(河海大学水文水资源学院，江苏 南京 210098)

受气候变化和高强度人类活动影响，中国北方河流江河径流量总体上呈现显著性减少趋势，在黄河流域，人类活动是径流变化的主要驱动因素[1]；
同时在空间分布上，人类活动的影响贡献自上游往下游逐渐增大[2]。在黄河中游黄土高原地区，大部分支流流域径流均表现为锐减趋势，人类活动影响为主导因素。鲍振鑫等[3]基于VIC模型对窟野河流域径流进行归因分析发现，影响期与天然期相比，人类活动对径流减少的贡献为61%～76%。Li等[4]运用数理统计方法对无定河和芦河流域径流归因发现，人类活动对两流域径流减少的贡献率分别为85%～90%和83%～86%。在渭河流域，Gao等[5]采用双累积曲线法计算得到人类活动对径流减少的贡献率为83%。然而在下垫面植被条件较好的延河流域，Wu等[6]基于SWAT模型计算得到气候变化的贡献率为56%，为径流减少的主要影响因素。薛帆等[7]基于Budyko假设和分形理论对北洛河流域径流进行归因识别，研究得出流域上游和下游径流减少主要受人类活动影响，而中游主要受气候变化影响。目前，径流变化归因分析基本是建立在水文模拟基础上的，很大程度上受模拟精度的限制。然而，黄土高原地区水文模型精度普遍不高[8-9]，特别是风沙区等特殊地貌条件下的水文模型适用性不强。本文从数据驱动途径改进风沙区典型流域(以陕北秃尾河为例)径流模拟方法，选择适用性较强的径流模拟方法对径流变化进行归因分析。

秃尾河发源于陕西省神木县瑶镇西北的公泊海子，全长140 km，流域面积3 294 km2(图1)，属于干旱半干旱风沙区，多年平均面雨量为404 mm，雨量在年内分配不均，汛期(6—9月)占比达75%[10]。研究数据包括90 m分辨率的数字高程数据DEM(来源于地理空间数据云，https://www.gscloud.cn/)、2个气象站的关键气象要素数据(来源于中国气象数据网，http://data.cma.cn/)以及9个雨量站的日降水和2个水文站的日径流数据。日潜在蒸散发量根据Penman-Monteith方程[11]计算。各站点数据根据泰森多边形法处理得到流域面平均系列。

1.1 趋势分析及突变检测

采用Mann-Kendall检验法(M-K法)[12]对降水、潜在蒸散发和径流系列进行趋势分析。根据M-K检验统计值Z大于0(或小于0)，确定系列呈上升(或下降)趋势，显著性水平α=0.1。

突变检测选择滑动t检验、有序聚类法(OC法)、Pettitt法、标准正态均一性检验(SNHT法)这4种方法[4,12]，并采用水文变异综合诊断的专家评分法[13]确定最终的突变点，具体步骤为：①每种方法评分总分为1，对处在可信度较低区间的点给予较低的评分，对在可信度较高区间的点给予较高的评分；
②若检验方法对样本检验结果所处区间无要求，则对检验的变异点给予相同的评分；
③对没有可信度的检测结果评分为0，然后采取加权平均和归一化处理计算各可能变异点的综合权重，综合权重最大的点即为最可能突变点。

1.2 天然期径流模拟方法

采用降水蒸发因子法、成因分割法、赵文林法、张经之法这4种方法[4,10]对天然期径流系列进行拟合，同时引入潜在蒸散发因子改进后3种方法，使其更为合理。

1.2.1降水蒸发因子法

直接建立径流量与降水量、潜在蒸散发量之间的统计关系，见式(1)：

Q=a×P+b×E+c

(1)

式中Q、P、E——年径流深、年降水量、年潜在蒸散发量；
a、b、c——天然期径流系列的回归参数。

1.2.2改进成因分割法

按径流的不同形成过程，成因分割法将径流分解为地下径流(基流)和地表径流(表流)，分别建立降水与地下径流、地表径流的相关关系，见式(2)、(3)：

Q=K1×Pm1+K2×P+C1

(2)

在式(2)中加入潜在蒸散发因子，则可改写为：

Q=K1×Pm1+K2×P+K3×E+C1

(3)

式中Q、P、E——年径流深、年降水量、年潜在蒸散发量；
K1、K2、K3、m1、C1——天然期径流系列的回归参数。

1.2.3改进赵文林法

赵文林法将径流系列划分为汛期径流和非汛期径流，分别建立汛期径流与汛期降水、上一期非汛期径流的相关关系，以及非汛期径流与非汛期降水、前一个汛期雨量与前一个汛期径流深之差的相关关系，见式(4)—(8)：

QX=K4×PXm2+K5×QFpast

(4)

QF=K7×PFm3×HAm4

(5)

引入潜在蒸散发因子，则可改写为

QX=K4×PXm2+K5×QFpast+K6×EX+C2

(6)

QF=K7×PFm3×HAm4+K8×EF+C3

(7)

Q=QX+QF

(8)

式中 QX、QF——汛期径流深和非汛期径流深；
PX、PF——汛期降水量和非汛期降水量；
EX、EF——汛期潜在蒸散发量和非汛期潜在蒸散发量；
QFpast——前一时段非汛期径流深；
HA——前一个汛期降水量与前一个汛期径流深之差；
K4、K5、K6、K7、K8、m2、m3、m4、C2、C3——天然期径流系列的回归参数。

1.2.4改进张经之法

张经之法根据汛期、非汛期降水量及汛期最大1日降水量，建立年径流量，见式(9)、(10)：

Q=K9(PX×fm5+PFm6)+C4

(9)

引入潜在蒸散发因子，则可改写为：

Q=K9(PX×fm5+PFm6)+K10×E+C4

(10)

式中Q、E——年径流深和年潜在蒸散发量；
PX、PF——汛期降水量和非汛期降水量；
f——汛期最大一日降水量与汛期降水量的比值；
K9、K10、m5、m6、C4——天然期径流系列的回归参数。

最优回归参数根据非线性回归拟合得到。输入模型表达式和参数初始值，经过多次迭代直到残差平方和最小时停止迭代，得到最优回归参数。

1.3 径流变化的归因分析

对比天然期、人类活动影响期的实测和模拟径流系列，建立径流变化的归因分析方法如下[3-5]：

ΔQT=QHR-QB

(11)

ΔQH=QHR-QHN

(12)

ΔQC=ΔQT-ΔQH

(13)

(14)

(15)

式中 ΔQT——径流变化总量；
ΔQH——人类活动对径流的影响量；
ΔQC——气候变化对径流的影响量；
QB——天然期实测径流量；
QHR——人类活动影响期的实测径流量；
QHN——人类活动影响期模拟的径流量；
ηH、ηC——为人类活动和气候变化对径流影响的贡献比例。

2.1 趋势分析与突变检测

1961—2015年流域年尺度降水和潜在蒸散发呈不显著下降趋势(统计特征量Z分别为-0.19和-0.16)，相应的年径流深为显著下降趋势(Z=-8.23)，见图2。4种方法检测得到秃尾河高家川站1961—2015年年径流系列的突变点主要有1979、1983、1996年(表1)，它们的综合权重分别为0.575、0.250、0.175，据此可确定主突变点为1979年、次突变点为1996年(1983年与主突变点位置接近，暂不考虑)。将整个径流系列划分为天然期(1961—1979年)、影响期I(1980—1996年)和影响期II(1997—2015年)3个阶段。支童等[14]在对秃尾河流域的降水、潜在蒸散发和径流系列进行趋势分析和突变检测时，也得出同样的结论。据此可判断上述具有一定可靠性。

突变检测方法可能突变年份评分滑动t检验19790.519960.5OC法19790.819960.2Pettitt法19831.0SNHT法19791.0

2.2 天然期径流系列拟合结果

采用上述7种数据驱动的统计方法对天然期年径流系列进行拟合，选用Nash-Sutcliffe效率系数(NSE)、Pearson相关系数(PCC)、均方根误差(RMSE)和平均绝对误差(MAE)进行拟合精度评定；
利用非线性回归拟合得到各方法的最优参数值见表2，拟合精度统计结果见表3。结果表明，改进后3种数据驱动方法的NSE、PCC、RMSE和MAE 4种精度指标较原方法均有大幅度提升；
总体上，改进赵文林法的拟合精度最高(NSE=0.77、PCC=0.88、RMSE=9.10、MAE=7.08)，可用于后续径流变化归因研究。

2.3 径流减少的归因分析

以实测降水资料为数据驱动，利用改进赵文林法重构2个影响期的径流系列，见图3，计算得到气候变化和人类活动对径流减少的贡献见表4。在影响期Ⅰ，年径流模拟值略高于实测值，表明该时期人类活动的影响不显著(对径流减少的贡献为19.1%)，而气候变化影响占主导地位(贡献率为80.9%)。相较于影响期I，影响期Ⅱ的多年平均年径流量比天然期减少幅度更大(-47.1%)；
年径流模拟值远高于实测值，表明该时期人类活动的影响较大，对径流减少的贡献提高至58.8%，占主导地位。

从1980—2015年逐年气候变化和人类活动影响贡献柱状图来看，影响期Ⅰ中多数年份气候变化贡献率都大于人类活动；
影响期Ⅱ中多数年份的人类活动贡献率均大于气候变化的贡献，且呈增长趋势，表明人类活动对径流减少的影响程度呈增长趋势。对比天然期与整个影响期径流变化，结果表明多年平均年径流减少37%，气候变化和人类活动影响的贡献率分别为54.2%、45.8%。

秃尾河流域在影响期Ⅰ中径流减少的主要驱动因素为气候变化，表现为多年平均降水量(339 mm)远低于天然期的数值(420 mm)；
但自1999年黄土高原退耕还林还草政策实施后，秃尾河流域植被覆盖面积显著上升，从而导致影响期Ⅱ中径流减少的主要驱动因素转变为人类活动[14]。

表2 7种数据驱动方法的参数取值

表3 天然期年径流拟合精度统计

注：表格中括号内容为未改进前拟合精度。

时期实测值/mm模拟值/mm减少比例/%人类活动的贡献绝对值/mm比例/%气候变化的贡献绝对值/mm比例/%天然期125.10124.80影响期Ⅰ92.8899.0525.86.1719.126.0580.9影响期Ⅱ66.17100.8447.134.6758.824.2641.2整个影响期78.78100.0037.021.2245.825.1054.2

利用M-K趋势分析方法发现，秃尾河流域1961—2015年年尺度降水、潜在蒸散发系列呈不显著下降趋势，而流域出口水文站高家川站年径流量呈显著减少趋势；
利用水文变异综合诊断方法得到年径流系列的突变点为1979年和1996年。引入潜在蒸散发因子，改进成因分割法、赵文林法和张经之法，增强了数据驱动方法在陕北风沙区秃尾河流域径流模拟中的适用性，其中改进赵文林法精度最高(NSE=0.77)。影响期(1980—2015年)与天然期(1961—1979年)相比，多年平均年径流量减少比例为37%，其中气候变化的影响贡献占主导地位(54.2%)；
但在影响期II(1997—2015年)，人类活动的影响更大，对径流减少的贡献比例为58.8%。