汉江流域安康段降雨径流特征分析及预测

刘易文，李家科*，丁 强，郝改瑞

(1.西安理工大学省部共建西北旱区生态水利国家重点实验室，陕西 西安 710048；2.陕西省环境监测中心站，陕西 西安 710054)

降雨被认为是水文循环的主要组成部分，是河川径流的主要来源，径流对水资源规划、管理、开发、模型构建非常重要[1-3]。而识别降雨径流特征变化以及精确模拟河川径流是水资源开发利用、水文模型构建的基础工作之一。由于降雨径流受气候、人为等诸多不确定因素的影响，常常伴随着非线性变化和较强的随机性，属于高度的非线性系统[4-6]。传统的径流预测方法只能反映出线性时间序列和简单的非线性时间序列，而实际径流往往受多因素影响，径流序列与各影响因素之间的非线性特征较强。因此，对复杂系统而言，仅凭传统单一模型，难以精确模拟整个径流序列。基于深度学习和机器学习的组合模型在很大程度上缓解了径流模拟精度低的问题，成为现阶段广受关注的热点之一[7-8]。组合模型更能体现时间序列的非线性特征，而且具有很强的非线性映射能力，适用于随机性较强的水文过程模拟[9]。此外，基于降雨径流演变规律和趋势特征，组合模型模拟的径流序列更为准确、可靠[10]。

近年来，随着深度学习和机器学习的推广，众多学者利用人工神经网络所衍生出来的组合模型进行水文序列模拟。王少丽等[11]利用BP神经网络对年径流系数进行了模拟，并与线性主成分回归模拟进行比对，研究发现BP神经网络的模拟结果更好，精度更高。梁川等[12]利用偏最小二乘回归及神经网络对高寒区域河川径流进行预测，结果表明偏最小二乘回归的模拟结果合理，但神经网络的模拟精度优于偏最小二乘回归。李秀峰等[13]将径流序列进行小波降噪处理后，利用偏最小二乘回归对径流进行了预测，结果表明该方法精度高且稳定性较好，可作为径流预测的有效方法。

总体看，径流预测的研究热点仍为基于深度学习、机器学习的组合模型。但少有在突变成分附近或者在径流局部对模拟效果进行优化的组合模型。基于此，本文依据汉江流域安康水文站控制断面以上1956—2017年的实测降雨、气温及径流数据，采用Mann-Kendall突变检验法识别降雨径流突变成分，并结合R/S分析法对降雨径流进行趋势分析，利用Mrolet小波提取降雨径流变化主周期。选用能较好识别时间序列的偏最小二乘回归(PLSR)及BP神经网络——偏最小二乘回归(BP-PLSR)对径流进行模拟，根据已识别的突变成分，对比2种方法在峰值处和突变成分附近的模拟精度，以期为流域水资源管理、水资源利用提供依据以及支撑。

1 研究区概况

汉江属于长江一级支流，全长1 532 km，发源地位于陕西省宁强县冢山，是南水北调工程的重要水源。汉江流域涉及湖北、陕西、甘肃、四川、重庆5个省(市)，面积达到15.9万km2。汉江流域整体地势西北高，东南低，属于亚热带季风气候区，平均年降雨量达873 mm，降雨充沛，汛期集中在5—10月，汛期径流量可占全年75%，年际变化大。安康水文站(109°E，33.67°N)位于陕西省安康市汉滨区(图1)。

图1 汉江流域安康段

2 材料和方法

2.1 数据来源

本文主要收集了汉江流域安康水文站控制断面以上1956—2017年多年实测降雨、气温及径流数据资料。其中降雨、气温数据来自于中国科学院资源环境科学数据中心(http://www.resdc.cn/)，径流数据来自于长江流域水文年鉴。

2.2 突变识别及趋势分析

2.2.1Mann-Kendall突变检验法

M-K突变检验是一种非参数方法，不需要测试样本遵循特定的分布[14]。同时，也不会受到少数样本数值的影响，适用于多种类型且计算较为简单[15]。

2.2.2R/S分析法

在M-K突变检验的基础上，采用R/S法分析汉江流域安康水文站控制断面以上降雨径流的趋势特征。R/S分析法的优势在于可以定量描述时间序列的趋势性，该趋势特征由Hurst指数表征，具有良好的稳定性[16-17]。

2.3 降雨径流周期分析

Mrolet小波是基于傅里叶变换的基小波函数与高斯函数推导而来，很好地克服了传统谱分析方法的缺点，在实际应用中比实数形式小波更具优势。Mrolet小波能够提高降雨径流序列周期变化检验水平，如今被广泛运用于信号处理、模式识别、数据分析等领域。

2.4 径流预测

2.4.1偏最小二乘回归(PLSR)

PLSR有机结合了多种主流分析方法，包括多元线性回归分析、主成分分析、相关分析，该方法优于最小二乘回归[18]，且适用于时间序列的模拟。首先构建自变量和因变量的数据表，分别为X和Y；对二者实施线性组合，得到PLSR第一成分t1与u1，进行X对t1的回归以及Y对u1的回归，当第一成分对应的回归方程满足精度要求时终止计算；若误差较大，利用自变量数据表与第一成分回归后的残余信息进行第二次成分提取，重复上述步骤至误差最小[19-21]。并根据所提取的成分信息，形成回归方程。

2.4.2BP神经网络-偏最小二乘回归(BP-PLSR)

BP神经网络是多层前馈神经网络，通常由输入、隐藏、输出层构成基础网络结构。该方法的核心在于反向传递(Back Propagation)，即依靠误差反向传播不断优化神经网络的权重和阈值，以期减小误差。尤其对于非线性较强的径流序列，BP神经网络能更好地识别局部误差，避免结果过拟合。BP神经网络构建如下。

a)对数据进行归一化处理：

(1)

式中xmin、xmax——原始数据的最小值及最大值。

b)本例中选用了3层神经网络结构，选用S型传递函数(Log-Sigmoid)作为激活函数，激活正向传递过程，通过反传误差函数不断调节网络的权值以及阈值，使误差函数E最小。

(2)

yi=f(Sj)

(3)

(4)

(5)

(6)

式中Sj——正向传递函数；Wi,j——结点i与节点j之间的权值；Wj——第j个节点输出值；f——节点激活函数；dj——第j个输出结果；b——阈值；N——隐藏层节点数；m——输入层节点数；n——输出层节点数；a——常数(范围取1～10)。

c)PLSR能较好地表征时间序列，但该方法在径流预测时容易陷入局部最优，造成序列过拟合。考虑到BP神经网络Levenberg-Marquardt算法(L-M算法)具有自适应性噪声消除功能，借此可以较好地消除目标序列噪声，防止PLSR的输出结果过拟合。借助Matlab中的Neural Net工具箱，将降雨、径流及气温值作为输入序列，选择Levenberg-Marquardt算法对原始数据进行处理，此时输出值为降噪后的径流序列。处理后的径流序列并不等同于径流模拟，将其标准化后作为PLSR的目标序列运行模型。标准化公式如下：

(7)

(8)

式中F0——Y的标准化矩阵；E0——X的标准化矩阵；E(y)——Y的均值；E(xi)——X的均值；Sy、Sxi——Y、X的均方差；n——样本数。

3 结果与讨论

3.1 降雨径流过程分析

汉江流域安康水文站控制断面以上1956—2017年多年平均降雨为815.61 mm，流量为587.25 m3/s；其中年均最大降雨量出现在2010年(1 231.9 mm)，年均最大径流量出现在1983年(1 294.42 m3/s)，年均最小降雨量及径流量均出现在1999年，分别为525.8 mm、247.75 m3/s，二者极差分别为706.1 mm、1 046.67 m3/s。多年降雨、径流过程见图2。

根据图2a，1956—2017年降雨量出现了4个变化阶段。经历了增加—减小—增加—减小的变化过程。其中，1974—1983年及1999—2001年降雨的上升趋势最为明显，2个阶段平均值分别为898.9、844.8 mm。通过62年整体变化趋势，判断降雨呈现出整体增长趋势。图2b中，62年径流变化同样存在4个阶段，与降雨序列一致，呈现出相同的变化过程。90年代前，平均径流量大于62年整体均值的年份明显多于90年代后。说明90年代后，流量整体减少趋势明显。2011—2016年，累计减少255.3 m3/s。径流序列整体呈现下降趋势，降雨径流变化过程见表1。

表1 降雨径流变化特征

3.2 降雨径流突变成分识别

汉江流域安康水文站控制断面以上降雨径流突变识别见图3。结果显示：M-K检验结果中降雨突变成分过多，存在伪突变点，结合累计距平法剔除伪突变点(图4)。对比2种方法分析结果，得到降雨突变点出现在1973年、1984年、2002年。1973年为先减小后增加的转折点，从1974年开始，降雨呈现持续增加现象，直到1984年，降雨增加17.4%。1984年突变前后降雨相对变化为30%。随后降雨开始突变减小，直到第三次突变，降雨减小23.2%。

径流序列突变点出现在1977年、1979年、1985年。分析发现1977—1979年径流量呈现小幅增加趋势，结合累计距平检验结果，认为1979年为伪突变点，固将其剔除。突变成分主要体现为先增加后减小，直到1985年，径流突变减小37.9%。

结合了M-K突变检验与累计距平检验，发现累计距平检验结果与M-K检验结果呈包含关系，结合分析降雨径流数据，对伪突变点进行剔除。根据M-K突变检验，发现降雨径流2个序列的统计量UF与UB的交点均在显著性水平线之间，所以突变检验通过显著性分析。

降雨径流突变检验结果与安康市近50年降雨量变化趋势相吻合[22-23]。安康市60年代初与80年代初降雨量为极大时期，同时60年代至80年代为偏丰时期，而在80年代以后降雨偏少。其中1974年突变点经历了降水的先减少后增大，从1984年突变点开始，降水从极大值突变减小。经检验，径流突变点滞后于降雨突变，对比其他学者的径流突变点分析结果[24]，由于径流突变检验的时间段选择与方法不同，安康径流突变仅出现在1984年，与本文结果相差1年，且径流突变均滞后于降雨突变，分析结果相似。

3.3 趋势分析

选用R/S分析法，设置lgm为横坐标，lg(R/S) 为纵坐标，并结合最小二乘法绘制散点图。根据散点拟合直线，直线斜率即为Hurst指数。Hurst指数不同，对应的趋势变化则不同，其不同取值对应的趋势见表2。R/S计算结果见图5。

表2 Hurst指数取值

a)降雨

a)降雨

a)降雨

a)降雨

a)降雨

b)径流

对R/S结果进行分析可得：降雨Hurst指数接近0，说明与未来趋势相反。降雨突变减小的趋势减弱，降雨序列在研究时段内将呈现增长趋势；径流序列Hurst指数为0.025 1，H指数介于0～0.5之间，说明该序列未来趋势将与过去相反，且H指数接近0，表明这种反持续状态较为强烈(表3)。根据径流量突变分析可以看出径流整体呈现减少趋势，且第二突变成分之后再未出现径流量跳跃增加或跳跃减小，同时验证了这种突变趋势正在减弱，该结果与张珏等[16]的分析结果一致。

表3 Hurst计算结果

3.4 小波分析

利用Mrolet小波对降雨径流进行周期分析，获得降雨径流的周期性变化特征(图6)。

b)径流

1956—2017年，可看出降雨量有6、10 a 2个周期。在10 a以下降雨量周期变化规律较为明显，2个周期的波动控制着安康降雨在整个时间域内的变化特征，震荡周期附近都存在着明显的丰、枯交替变换，其中偏丰期、偏枯期交替突变的点出现在1958年、2013年。对于径流周期变化，62年内安康径流量存在周期14 a，此处震荡最强，为主周期。5 a处震荡较弱，为副周期。16 a以下径流周期变化规律较为明显，其中偏丰期、偏枯期交替突变的点出现在1979年、1983年。

3.5 径流预测

首先对气温序列、降雨及径流序列进行相关性检验，分析得到径流与气温以及降雨与径流之间均存在相关性(P<0.005)，通过了显著性检验(表4)，本例中所选用的自变量和因变量之间存在较为明显的相关性。自变量之间存在相关性会造成共线性，进而影响回归结果，而偏最小二乘回归会减小自变量共线性对结果造成的影响。

表4 相关性检验(新增降雨及气温之间的相关性检验)

基于1956—2017年62 a的实测年降雨及气温作为基础输入数据，分别利用PLSR及BP-PLSR 2种方法对径流量进行预测。将1956—2004年设置为训练期，2005—2017年设置为验证期。预测精度指标选择均方根误差(Root-Mean-Square Error,RMSE)及纳什效率系数(Nash-Sutcliffe Efficiency Coefficient ,NSE)，精度指标计算公式如下：

(9)

(10)

利用PLSR对径流进行模拟。首先根据平均值与均方差对原始数据进行标准化处理，提取出2个成分，得到PLSR回归方程：

y*=1207.9+0.9x1-83x2

(11)

式中y*——预测年径流量；x1——实测年降雨量；x2——实测年气温值。

PLSR径流模拟结果见图7。在模型训练期，RMSE为153.093，NSE为0.489；验证期内，RMSE为70.275，NSE为0.259；在整个模拟时期内，PLSR模拟结果精度较低，RMSE为152.182，NSE为0.456，纳什效率系数(NSE)低于0.5且R2为0.500 3。模拟精度低的主要原因为：62 a的时间周期内，径流量出现突变，在径流峰值前后出现了过拟合现象，且在突变点附近出现局部最优问题。1983年为径流最大年份，径流峰值处模拟值为944.58 m3/s，与实测值的误差为349.84 m3/s；在突变点附近，1977、1979、1985年的模拟值分别为678.30、781.00、635.19 m3/s，与实测值的误差分别为288.6、297.8、254.9 m3/s。

a)不同年份下径流量实测、模拟值

b)训练期实测-模拟

c)验证期实测-模拟

d)实测-模拟

采用BP-PLSR对径流进行模拟，模拟结果见图8。所生成的回归方程为：

y*=1333.7+0.6x1-78.7x2

(12)

式中y*——预测年径流量；x1——实测年降雨量；x2——实测年气温值。

在训练期RMSE为24.371，相比PLSR在验证期的结果减少了84.4%。NSE为0.784，明显大于PLSR在训练期的效率系数；在验证期，RMSE为12.343，比PLSR在验证期的结果减少了82.5%。NSE为0.891，BP-PLSR模拟结果精度明显优于PLSR。在整个模拟时期，RMSE为92.863，误差下降了40%。NSE明显增加，达到0.797，模拟结果良好。1983年峰值处，实测值为1 294.42 m3/s，模拟值为1 114.60 m3/s。相比PLSR，误差下降了48.6%。在突变点附近，1977年、1979年、1985年模拟值分别为613.94、498.77、593.92 m3/s。相比PLSR，误差分别降低了22.30%、94.76%、91.40%。

未对原始数据进行预处理，用PLSR进行模拟后，发现模拟结果精度较低，特别是在突变成分附近及峰值处的模拟结果精度低。对原始数据进行BP神经网络预处理，可以解决数据在峰值处的过拟合问题，能较好地避免模型在突变成分附近陷入局部最优。BP-PLSR模拟结果与实测值最为接近，模拟结果可信度高。

BP-PLSR模拟效果优于传统PLSR，模拟精度有所提高。但在径流局部区域难免出现过拟合现象，只有在1977年处，2种方法的模拟结果相对误差均大于50%(表5、图9)。在训练期，BP-PLSR对过拟合及局部最优问题控制较好，未出现明显泛化误差(表6)。由于PLSR容易在径流峰值处及突变点附近出现明显泛化误差，而BP-PLSR可以优化上述问题，能明显提高模拟精度。

表5 突变成分附近及峰值处的对比

a)不同年份下的径流量实测、模拟值

b)训练期实测-模拟

c)验证期实测-模拟

d)实测-模拟

图9 局部模拟结果

表6 不同模拟方法对应的径流模拟精度指标

将突变成分附近及峰值附近实测径流序列作为目标序列，选择对应时段降雨及气温值进行模拟，将结果回代到全序列以判断BP-PLSR的适用性(图10、11)，得到的回归方程为：

y*=1452.3+0.9x1-101.5x2

(13)

所选时段为1964—2002年，NSE系数为0.776，决定系数R2为0.809 4。

在突变点1977年、1979年、1984年处，模拟值与实测值误差为-224.185、-15.570、-21.310 m3/s，1977年模拟值误差较大。峰值处模拟与实测值误差为179.77 m3/s。整体模拟效果较好。将结果回代到全序列中时，得到NSE值为0.798，决定系数R2为0.829 5。整体模拟效果理想，但是仍然存在个别突变点及峰值处模拟误差较大的现象。通过2种模拟方式，即利用全序列体现局部效果及依靠局部回代全序列，都体现了较好的模拟结果，NSE系数均接近0.8，R2超过0.8，模拟结果良好，说明了该组合模拟适用性较好。

a)1956—2017年径流量实测、模拟值

b)实测-模拟

a)1964—2002年径流量实测、模拟值

b)实测-模拟

4 结论

a)62 a汉江流域安康水文站控制断面以上多年平均降雨为815.61 mm，1974—1983年、1999—2001年降雨上升趋势最为明显，2个阶段平均值分别为898.9、844.8 mm，降雨呈现出整体增长趋势。62 a径流均值为587.25 m3/s ，90年代前有19 a的径流量大于整体平均值587.25 m3/s，而90年代后只有5 a的径流大于587.25 m3/s。径流序列多个阶段均存在减少，径流整体呈现下降趋势。

b)降雨突变成分出现在1973年、1984年、2002年，径流序列突变成分出现在1977年、1985年。选择显著性水平为5%，降雨径流突变识别通过显著性检验。

c)利用R/S分析，降雨径流序列Hurst指数均介于0～0.5之间，说明未来将存在相反趋势；H指数接近0，表明这种反持续状态较为强烈。降雨径流趋势分析相关系数分别为0.941 5、0.894 3，趋势分析效果良好。

d)PLSR的径流模拟结果精度较低，RMSE为152.182、NSE为0.456。NSE值低于0.5，误差较大，模拟结果受到了过拟合及局部最优影响。对原始数据进行BP处理后，能有效提高模拟精度，RMSE为92.863、NSE为0.797；RMSE降低40%且NSE大于0.5。在峰值处，模拟精度提升48.6%；在突变点附近，精度提高了69.5%。模拟结果可信度增加，有效解决和避免了过拟合及局部最优问题。

5 结语

基于汉江流域安康水文站1956—2017年62 a实测降雨、气温及径流数据，利用M-K突变检验进行突变点识别，选用R/S分析法以及小波分析对安康水文站控制断面以上降雨径流序列进行趋势及周期分析。采用PLSR、BP-PLSR对径流进行模拟，并对两者进行了对比。由于水文序列随机性强，不确定因素较多，受气候、下垫面以及人类活动的影响大，传统方法难以较好模拟非线性序列，建议充分利用深度学习和机器学习，优化现有模拟方法，提高模拟精度。