年径流量序列还现方法研究

薛树文,曹升乐,王利朵,刘 阳

(山东大学土建与水利学院,山东济南250061)

年径流量序列还现方法研究

薛树文,曹升乐,王利朵,刘 阳

(山东大学土建与水利学院,山东济南250061)

受人类活动的影响,不同时期年径流量的影响因素不同。特别是由于大型工程的影响,实测年径流量序列已不能真实反映河川径流的天然特性。基于年径流量序列的变化成因分析,通过对暂态成分的排除,得到代表性、一致性较好的新序列,用于后期水资源评价及计算。以黄河高村站实测径流序列为例,分别对暂态成分进行线性拟合及分阶段线性拟合并排除,结果显示,分阶段线性拟合方式较第一种更加合理。计算结果表明还现后的年径流量序列均值与波动范围均与实际情况相符。

水文序列；一致性；暂态成分；年径流量；黄河高村站

0 引 言

人类活动[1]破坏了实测径流的一致性,在以往进行水资源调查评价时,需要对实测径流资料进行还原处理。这种还原法的正确性毋庸置疑,但是实际工作中发现存在计算精度不太可靠、应还原的项目不能全部考虑及成果使用不方便等问题[2]。《全国水资源综合规划细则》中提到的还原(向前还原)方法和后期提出的还现(向后还原)均是将水文序列恢复到过去和现在下垫面条件下的天然径流；而实际工作中,大部分径流计算均是基于现有工程和取引水的前提下进行的,这种还原和还现方法不能很好地适应环境变化的需求及实际工作的需要。为此,国内学者进行了大量的研究,主要有突变点前后系列与同参数关系分析法[2-3]、水文序列分解合成法[4-7]、水文模型法[8-12]。本文基于水文序列分解合成法,提出了一种根据水文序列变化成因分析的非一致性径流序列还现方法,可以较好地适应环境变化和后期工作需求。

1 研究方法

水文现象随时间变化的过程称为水文序列,一定时段内的水文序列同时受自然气候、地理地质、人类活动等共同影响,资料的形成和变化受多种因素的影响。再复杂的水文现象都可以分解为两部分,即确定性成分和随机性成分[13]。水文序列的确定性成分包括周期成分和非周期成分,随机性成分包含平稳成分和非平稳成分,其常用线性叠加的形式表示。即

Xt=Nt+Pt+St

(1)

式中,Nt为确定性的非周期成分(包括趋势、跳跃、突变)；Pt为确定性的周期成分(包括简单周期、复合周期、近似周期)；St为纯随机成分。

一般来讲,水文序列的随机成分主要受气候、地质等因素的影响,相对稳定,变化周期漫长。因此,其随机性成分的一致性是相对较好的；而其确定性成分主要受气候周期和人类活动的影响。本文主要研究年径流量变化过程,周期成分不明显,因此确定性成分多为非周期成分,水文学中又将非周期成分称为暂态成分。暂态成分的变化规律可以在较短的工程年代内发生变化,可能是由于植树造林、水土保持、工农业引水等活动引起的缓慢的渐变,也可能是河道内水工建筑物的修建引起的剧烈的突变,在水文序列分析中常表现为趋势、突变、跳跃等,因此水文序列中暂态成分的一致性往往较差[4]。基于以上分析,本文假设水文序列由一致性较好的随机成分和一致性较差的暂态成分组成[5- 6]。

若序列中呈现趋势或跳跃成分,亦即暂态成分；则意味着序列的一致性受到破坏。假设水文序列由暂态成分和随机成分组成,即

Xt=Nt+St

(2)

则当暂态成分Nt仅为跳跃成分时,Nt为一确定常数；当暂态成分Nt仅为趋势成分时,Nt为时间t的函数；当同时出现跳跃和周期成分时,Nt为实际序列。基于最小二乘法拟合得到的函数[6],暂态成分Nt可用多项式来描述。即

Nt=a+b1t+b2t2+…bptp

(3)

式中,a为常数；b1，b2，…，bp为回归系数。

由式(3)可见,实际工作中Nt可以是常数、线性函数,也可以是非线性函数。

本文采用线性趋势的相关系数检验法[6]、Kendall秩次相关检验法[14-16]、滑动平均法、有序聚类分析法[17]、秩和检验法[13]对序列暂态成分中的趋势成分及跳跃成分进行识别和显著性检验；然后,分别对暂态成分进行线性拟合、分阶段线性拟合并排除,将水文序列还现到现在水平。

由于人类活动常常是一个渐变的过程,受人类活动影响的水文要素也会呈现一个缓变过程。因此,对于实测水文序列,可采用线性方程来拟合。即

(4)

不同时期人类活动的强度可能存在明显差异,对水文要素的影响也截然不同。若用一条直线去拟合一个长序列,可能不能很好反映人类活动的影响,拟合误差可能较大。因此,可根据社会发展的不同时期,如中国的改革开放前后可分为两个时期,或者大型工程建成前后也分为两个时期,将序列分为几段来进行拟合。即

(5)

ak-1·nk-1+bk-1=ak+bk,k=1,2,…,m

(6)

式中,m为所分段数；nk为第k段序列长度；ak,bk分别为第k段序列拟合参数。

按式(5)、(6)式求解拟合方程,可避免每段拟合存在节点处不连续的问题。根据其目标函数及约束条件,进行非线性规划的求解[18],即可得到最优的参数值。对于序列较长的非线性规划的求解可借助软件完成。实际工作中,根据各阶段的实际情况,还可继续增加约束条件,以满足实际工作要求。

2 实例分析

黄河发源于青藏高原,流经9省(自治区),自山东东营注入渤海。截至2000年,黄河流域已建成大中小型水库及塘堰坝等蓄水工程近20 000座,引水工程9 860处,提水工程23 600处,保证了沿岸50多座大中型城市和420个县(旗)城镇及733.3万hm2耕地的用水[19]。众多的引提水工程及各种水利枢纽的修建导致黄河干流水文站尤其是下游水文站实测资料的一致性较差。为此,本文选取黄河下游高村站1951年～2014年共64年实测径流资料进行分析。

2.1 暂态成分的识别及显著性检验

暂态成分,又称确定性非周期成分,主要包括趋势、跳跃和突变(跳跃的一种特殊情况)。因此,本节主要对趋势成分及跳跃成分进行识别及显著性检验。

2.1.1 趋势成分的识别及显著性检验

对黄河高村站1951年～2014年实测年径流资料使用滑动平均法进行趋势识别,并用Kendall秩次相关检验验证趋势成分的显著性,黄河高村站年径流局部趋势变化曲线见图1。

图1 黄河高村年径流局部趋势变化曲线

2.1.2 跳跃成分的识别及显著性检验

对于同时存在趋势、跳跃成分的序列,根据序列实际值,基于最小二乘法进行拟合。

2.2 暂态成分的拟合

根据最小二乘法估计回归系数,可得黄河高村站年径流量回归方程为：Nt′=-5.234 4t+518.92；t=1，2,…,64。

2.3 暂态成分的排除

若某水文序列一致性较好,则该水文序列各值会围绕均值上下波动。本文将序列还现到2014年水平,假设新序列均值过直线最后一点(t=64)的一条水平线,其值为N64′=183.92。它反映该序列现状平均水平,因此该年径流量序列的暂态成分即为Nt′-N64′。亦即,Nt=-5.234 4t+335.001 6;t=1,2,…,64。

2.4 新序列的生成

如前文所述,将水文序列划分为一致性的随机成分和非一致性的暂态成分,根据需要,排除水文序列Xt中的暂态成分,即可得到一致性较好的新序列,新序列；Xt新=Xt-Nt=Xt+5.234 4t-335.001 6；t=1,2,…,64。原始序列与新序列对比见图2。

图2 黄河高村站年径流量原始序列与新序列趋势变化(1)

3 分析与改进

3.1 结果分析

由图2可以看出,采取以上方法计算得到的新序列一致性较好,但是存在以下问题：

(1)线性拟合不能很好地代表年径流量的变化趋势,新序列均值186.51亿m3/a,与2011年～2014年实测均值286.38亿m3/a相差较大。

(2)线性拟合导致新序列的某些结果与实际不符,如1997年出现14.43亿m3/a的极低流量,需修正。

3.2 方法改进

基于以上问题,本文提出了一种根据原始径流特征,分阶段进行线性拟合的改进方法：

根据黄河高村站年径流量变化特征,将1951年～2014年年径流量资料以突变点(1985年)和小浪底水库全面建成年份(2001年)为分界点,分三个阶段进行线性拟合。

1951年～1985年,年径流量波动加大,进行Kendall秩次相关检验得U=-1.46,取显著性水平α=5%,可得该时段黄河年径流量存在减少趋势但不明显。这种减少趋势主要是上游植树造林、水土保持减水及取引水的缓慢增加所致[20]；同理可得,1985年～2001年径流量由于经济的快速发展,地表水的大量取引及地下水的过量开采,年径流量呈现明显的下降趋势；2002年～2014年径流量因小浪底水库的调蓄作用,而使年径流量呈现明显的上升趋势。对这三个阶段的暂态成分分别进行线性拟合,从而形成新的序列。

具体的实施方法如下：首先将1951年～1985年径流量资料还现到1985年水平；再将还现得到的1951年～1985年的径流量与1986年～2001年的径流量还现到2001年水平；最终将1951年～2000年的径流量与2001年～2014年的径流量还现到2014年水平。之后，根据各阶段趋势不同的特点,分三次进行暂态成分的拟合及排除,最终将1951年～2014年的径流量资料的还现到2014年水平,形成新的序列。

1951年～1985年Nt1=-3.77t+504.05；t=1,2,…35

1985年～2001年Nt2=-14.06t+386.14；t=1,2,…17

2001年～2014年Nt3=12.75t+134.43；t=1,2,…14

根据各阶段线性回归方程,进行暂态成分的排除和新序列的生成,新生成的序列与原始序列趋势变化见图3。

图3 黄河高村站年径流量原始序列与新序列趋势变化(2)

3.3 对比分析

对比两种方法,发现改进方法可以较好解决第一种方法遇到的问题。首先采用改进方法生成新序列均值为314.37亿m3/a,年径流量变化范围为[152.21,734.79]亿m3/a,最值分别出现在1960年和1964年。与2011年～2014年实测均值286.38亿m3/a相比,误差为9.77%。改进方法计算得到的年径流量变化范围较原始序列的[103.41,873.17]缩减,缩减率为24.32%。可以看出,采用改进方法计算得到的新序列代表性和一致性较好,可以很好地适应工程计算的需求。

4 结 语

本文基于水文序列组成成分分析,假设非一致性的水文序列由一致性较好的随机序列和一致性较差的确定性趋势成分组成,通过对趋势成分的线性、分阶段线性拟合及趋势成分的排除,进行非一致性水文序列的还现计算,从而得到如下几点结论：

(1)两种方法计算得到的新序列均有较好的一致性,但是由于趋势拟合的合理性问题,第一种方法计算得到的新序列与实际情况相差较大,改进方法计算结果较为理想；

(2)采用改进方法生成新序列均值为314.37亿m3/a,与最近几年(2011年～2014年)实测均值286.38亿m3/a较为接近。生成的新序列波动振幅减小,与小浪底水库修建运行的实际情况相符。

(3)分阶段拟合法可根据未来径流量的变化趋势,根据水文序列成因分析,继续增删阶段,以更好的适应未来水资源规划及计算的需求。

[1]水利部水利水电规划设计总院. 全国水资源综合规划技术细则[R]. 北京： 水利部水利水电规划设计总院, 2008: 148- 154.

[2]陆中央. 关于年径流量系列的还原计算问题[J]. 水文, 2000, 20(6)： 9- 12.

[3]韩瑞光, 丁志宏, 冯平. 人类活动对海河流域地表径流量影响的研究[J]. 水利水电技术, 2009, 40(3)： 4- 7．

[4]谢平, 陈广才, 夏军. 基于趋势分析的非一致性年径流量序列频率计算方法[M]. 郑州： 黄河水利出版社, 2005: 80- 85．

[5]谢平, 陈广才, 雷红富. 变化环境下基于趋势分析的水资源评价方法[J]. 水力发电学报, 2009, 28(2)： 14- 19．

[6]谢平, 陈广才, 夏军. 变化环境下非一致性年径流序列水文频率计算原理[J]. 武汉大学学报： 工学版, 2005, 38(6)： 6- 10．

[7]胡明义, 梁忠民. 基于跳跃分析的非一致性洪量系列的频率计算[J]. 东北水利水电, 2011(7): 38- 40．

[8]王国庆, 张建云, 刘九夫, 等. 气候变化和人类活动对河川径流影响的定量分析[J]. 中国水利, 2008(2)： 55- 58．

[9]王国庆. 气候变化对黄河中游水文水资源影响的关键问题研究[D]. 南京： 河海大学, 2006．

[10]张爱静. 东北地区流域径流对气候变化与人类活动的响应特征研究[D]. 大连： 大连理工大学, 2013．

[11]王亮, 高瑞忠, 刘玉才, 等. 气候变化和人类活动对滦河流域内蒙段河川径流的影响分析[J]. 水文, 2014(3)： 70- 79．

[12]梁忠民, 胡义明, 王军. 非一致性水文频率分析的研究进展[J]. 水科学进展, 2011, 22(6)： 864- 871．

[13]王文圣, 丁晶, 金菊良. 随机水文学[M]. 2版. 北京： 中国水利水电出版社, 2008．

[14]于延胜, 陈兴伟. 基于Mann-Kendall法的水文序列趋势成分比重研究[J]. 自然资源学报, 2011.09, 26(9)： 1585- 1591．

[15]HAMED K H. Trend detection in hydrologic data: The Mann-Kendall trend test under the scaling hypothesis[J]. Journal of Hydrology, 2008, 349(3- 4)： 350- 363．

[16]BURN D H, HAG ELNUR M A. Detection of hydrologic trends and variability[J]. Journal of Hydrology, 2005, 255(1- 4)： 107- 122．

[17]张敬平, 黄强, 赵雪花. 漳泽水库水文序列突变分析方法比较[J]. 应用基础与工程科学学报, 2013, 21(5)： 837- 841．

[18]尚松浩. 水资源系统分析方法及应用[M]. 北京： 清华大学出版社, 2006.

[19]张学成, 潘启民. 黄河流域水资源调查评价[M]. 郑州： 黄河水利出版社, 2006．

[20]王乐平. 基于小波变换的黄河下游水沙变化特征及其成因分析[D]. 太原理工大学, 2015．

[21]田垅, 刘宗田. 最小二乘法分段直线拟合[J]. 计算机科学, 2012, 39(6)： 482- 484．

[22]薛丽红. 基于最小二乘法的分段直线拟合算法[J]. 贵阳学院学报： 自然科学版, 2015, 10(4)： 9-10．

[23]王继强. 基于LINGO的最小二乘拟合的数学规划解法[J]. 信息技术, 2009(8)： 74- 79．

(责任编辑 陈 萍)

Research on Forward Restore Method of Annual Runoff

XUE Shuwen, CAO Shengle, WANG Liduo, LIU Yang
(School of Civil Engineering, Shandong University, Jinan 250061, Shandong, China)

Affected by human activities, the influence factors of annual runoff is different during different periods. Especially because of the influence of large engineering, the composition of annual runoff series can not reflect the natural characteristics of runoff. Based on the analysis on change reason of runoff series, a consistency and representative new sequence can be got by eliminating transient components for later evaluation and calculation of water resources. Taking the runoff series of Gaocun Hydrometric Station as example, different ways are used to fit and eliminate the transient components, and the results show that the linear fitting method in stages is more reasonable than linear fitting. The calculation results also show that the average and change range of current annual runoff is consistent with actual situation．

hydrological series; consistency; transient component; annual runoff; Gaocun Hydrometric Station of Yellow River

2016- 06- 02

水利部公益性行业科研专项资助项目(201501054)

薛树文(1991—),男,山东聊城人,硕士研究生,研究方向为水文学及水资源；曹升乐(通讯作者)．

TV121.4

A

0559- 9342(2017)05- 0021- 04