水文变异条件下潦河生态流量计算研究
2021-03-25刘成建陈明昊
王 强, 夏 瑞, 邹 磊, 陈 焰, 张 远, 刘成建,4, 陈明昊,5
1.中国环境科学研究院水生态保护修复研究室, 北京 100012 2.武汉大学, 水资源与水电工程科学国家重点实验室, 湖北 武汉 430072 3.中国科学院地理科学与资源研究所, 陆地水循环及地表过程重点实验室, 北京 100101 4.西北大学城市与环境学院, 陕西省地表系统与环境承载力重点实验室, 陕西 西安 710127 5.中国科学院重庆绿色智能技术研究院三峡生态环境研究所, 重庆 400714
气候变化和高强度人类活动加速了流域水文循环过程,进而改变了降水和地表水资源的时空分布,导致径流序列的一致性遭到了破坏[1-3]. 地表水文循环过程的显著变化,改变了流域产汇流过程的原有属性,导致径流序列出现变异,给流域生态系统带来严重影响. 此外,随着社会经济的发展和人类日益增长的用水需求,高强度的人类活动(如跨流域调水、水利工程建设、城市化建设等)对流域水文过程的影响将进一步加剧,给变化环境下水资源的适应性管理带来了严峻挑战[4-6]. 面对水资源量时空分布不均、水质污染严重、水生态环境日益恶化等问题,科学、合理估算河流生态流量将为保障河流生态系统健康良性发展、促进流域生态建设和水资源合理配置提供重要依据[7-9].
自20世纪七八十年代以来,生态流量的相关研究开始蓬勃发展. 据不完全统计,全球生态流量的计算方法超过200种,大致可分为4类,包括利用历史流量资料的水文学方法(如Tennant法[10]、Texas法、NGPRP法等)、根据河道断面水力特性的水力学法(如湿周法[11]、R2CROSS法等)、以指示物种所需的水力条件为依据的栖息地法(如WUW法、IFIM法[12]等)以及综合考虑研究区生态环境整体的综合法(如BBM法[13]等). 虽然各种方法都有各自的生态学理论基础,但均有一定的适用条件和不足,目前针对河道生态流量的计算还并没有统一的标准,各种方法仍在不断发展中,不同方法的适用性也存在较大的争议.
由于水文资料易于获取、不需要开展现场测定数据,能够利用长序列径流资料分析统计特性,且具有较强的通用性,因此,水文学方法在全世界范围内得到了广泛的推广和应用. 如朱才荣等[14]采用9种水文学方法估算了汉江中下游基本生态流量,提出了表征生态需水动态特征的生径比概念,对不同水文学方法生态流量进行生径比分析;李紫妍等[15]采用变化范围法(RVA)、多年日流量资料排频法等分析了汉江子午河生态流量过程,研究该河道维持多种生态环境目标所需的生态流量;龙凡等[16]通过建立年月丰枯遭遇的Copula联合分布,求得不同典型年条件下各月丰平枯的条件概率,利用改进的FDC法计算了丰平枯典型年的生态流量过程;王俊钗等[17]以生径比为标准,运用流速法、月最小径流法、逐月频率设定法选取动态生态流量过程;张涛等[18]采用多项式各时段的流量频率曲线,建立流量与频率的函数关系,以某枯水期频率计算该时段河道的最小生态流量. 但上述研究均未充分考虑水文序列的一致性条件是否发生变化. 随后,国内外部分学者开展水文变异条件下的生态需水研究,如肖才荣等[19]根据秩和检验法分析东江水文变异,并以变异点前各月平均流量序列进行分析推求生态流量;李剑锋等[20]采用滑动秩和检验法分析黄河干流水文变异,并对河道内生态流量进行研究;刘剑宇等[21]分析了鄱阳湖流域的水文变异及其成因,并估算了鄱阳湖流域五河的生态流量;孟钰等[22]分析了淮河干流的水文变异,以典型鱼类长吻鮠为保护对象,采用改进的FLOWS法对水文变异前后长吻鮠的生态环境进行对比分析.
由于受到气候变化和人类活动的影响,流域自然水文气象序列往往容易发生突变,从而导致水文序列的时空变异,显著影响当地生态系统完整性与河流健康状态. 然而,目前针对水文变异条件下的生态流量研究大多是针对变异点前水文序列进行分析,但往往变异点前水文序列数据有限,限制了以长序列资料为数据基础的水文学方法应用. 鄱阳湖是我国第一大淡水湖,与洞庭湖并称“长江双肾”,对长江中游水环境和水生态保护具有十分重要的作用. 在气候变化和强人类活动作用下,鄱阳湖流域的管理和保护愈发重要,是当前长江保护修复的关键环节之一. 基于此,该研究以江西省鄱阳湖西北部的潦河为例,通过实测水文数据分析水文序列的变异性,提出了基于水文变异条件下河道内生态流量的计算方法,旨在为探索潦河生态需水标准的制定提供数据参考和技术支撑.
1 研究区域与数据来源
潦河流域位于江西省鄱阳湖区西北部,为修水的最大支流,又称为奉新江. 潦河流域以万家埠水文站位出口点,地理位置为114°87′E~115°38′E、28°38′N~29°04′N,流域面积 3 548 km2. 潦河流域地势西高东低,中游为高山丛林地貌,植被良好,下游为丘陵岗地,植被稀少,水土流失严重. 流域内降水量较为充沛,属于半湿润地区,当地产生暴雨的天气系统主要是台风、高空槽、低压冷锋等,特大暴雨往往由多种天气系统配合形成,多年平均年降水量为 1 600 mm,但受地形影响,降雨分布极不均匀,上游山区最大年降雨量记录超过 2 000 mm. 潦河流域属于暴雨多发区,为江西省四大暴雨中心之一,加上河床坡降较大,洪水汇流速度快,因其特殊的地理位置、气候因素和洪水特点,洪灾频繁发生.
该研究收集了潦河流域及其周边11个气象站1963—2014年的逐日气象观测数据,站点分布较为均匀,能够基本反映流域气象水文变量的时空变化,部分缺测数据通过与相邻气象站点相对应的气象变量序列建立回归关系进行插补. 径流数据来自于万家埠水文站1953—2014年的水文要素摘录表,包括日均流量和月均流量,用于对模型进行参数率定和模型验证. 水文模型采用的数字高程模型(DEM)、土地覆盖和土壤类型数据均来自中科院资源环境科学数据中心(http:www.resdc.cn),流域位置及站点分布如图1所示.
图1 潦河流域位置及站点分布Fig.1 The location of Liaohe River Basin and stations
2 研究方法
2.1 水文变异诊断方法
径流变异点检测方法众多,常见的有基于贝叶斯理论的里海哈林法、Mann-Whitney非参数检验法、有序聚类法、RS分析法、F检验法、T检验法、滑动秩和检验法和Mann-Kendall非参数统计检验法、Pettitt突变点检验法、SNHT(Standard Normal Homogeneity Test)方法、Buishand检验法和Buishand秩和突变检验法等. 单一检测方法在径流突变点检验中存在局限性,因此该研究采用Mann-Kendall非参数统计检验法、Pettitt突变点检验法、SNHT方法、Buishand检验法和Buishand秩和突变检验法对突变点进行综合评估. 其中,Mann-Kendall非参数统计检验法所需样本无需遵从预先假定的分布,且不受少数异常值的干扰,计算简便实用,广泛应用于水文气象时间序列的趋势和突变检验. 该方法具体计算步骤:对于具有n个样本的水文气象时间序列X,其观察值序列为{xi,i=1,2,3,…,n},假设该时间序列相互独立且有相同连续分布,则按照时间序列顺序和逆序构建统计量UF和UB.
(1)
(2)
其中,
(3)
式中,xi为第i个水文气象变量,n为序列长度,Sk为秩序列,E(Sk)为方差,VarSk为期望,UF和UB为统计量. 通过进一步分析统计量序列UFk和UBk可以获取时间序列的趋势变化情况,可以明确突变时间.
2.2 分布式时变增益模型
由于气候变化和高强度人类活动的影响,目前水文站实测径流过程对于表征天然径流情况存在较大局限,因此需要对测站以上受人类活动影响的水量进行还原,即为天然径流还原. 目前,采用天然径流还原方法有多种,其中应用最为广泛的是分项调查法、蒸发差值法、降雨径流关系法和水文模型法等[23-24]. 其中,水文模型是对自然界中复杂降雨径流水文循环过程的概化和近似描述,以流域产汇流理论为基础,可以综合考虑气象因素、地形、土地利用变化、土壤类型等的空间不均匀性,具有较为完善的物理基础,且随着模型软件的不断成熟,其获取方式更为便捷,应用更加广泛.
该研究采用分布式时变增益模型(Distributed Time-Variant Gain Model, DTVGM)[25-27]在径流突变检验的基础上进行径流还原. DTVGM模型是将XIA等[28-30]提出的水文非线性时变增益模型(TVGM)通过DEMGIS平台进行拓展,进而推广到流域水文变量时空模拟的分布式水文模型. 该模型通过建立不同下垫面状况(土地利用、土壤类型等)与水文非线性产流参数之间的关系,一方面具有分布式模型概念性模拟的特点,另一方面具有水文非线性系统分析适应能力强的优点,该模型已经在国内外诸多流域水文模拟方面进行了应用和验证. 针对分布式模型水文模拟,首先根据潦河流域数字高程模型(Digital Elevation Model,DEM)对流域进行子流域划分,进而选取两种评价指标——Nash-Sutchliffe效率系数(NSE)[31]和相关系数(r)作为评估模型性能的评价指标对模型进行参数率定和验证. 其中,NSE反映了观测值和对应模拟值之间的拟合程度,可用式(4)计算:
(4)
(5)
2.3 生态需水计算方法
该研究主要采用逐月频率计算法估算生态流量,逐月频率计算法即对多年各月的流量序列进行排频和适线分析,一般取90%的保证率作为河道内生态流量. 该方法一方面可以考虑河流生态系统对水量的要求,另一方面还能够根据不同月份径流的变化特征,考虑不同月份生态环境需水的不同要求. Green等[32]研究指出,任何一种概率分布函数只可能对某种分布或子样本容量的检验效果比较好,并不存在一种占绝对优势的概率分布函数. 因此,该研究采用两参数正态分布(ND)、伽玛分布(GD),三参数广义极值分布(GEV)、皮尔逊Ⅲ型分布(P-Ⅲ)、对数正态分布(LN3)函数对经模型还原后的各月径流序列进行分析,各分布详细介绍如表1所示. 为了寻求最适宜于描述月径流过程的概率分布线型,该研究采用Kolmogorov-Smirnov(K-S)、Anderson Darling(A-D)和概率点据相关系数(PPCC)3种检验法[33]进行拟合优度检验分析. 其中,PPCC检验法是假设待检验样本服从某种分布,排序后的观测值与假设的理论分布值的概率点据呈近似线型相关,相关系数r即可用来比较不同假设分布函数线型的优劣.
表1 该研究采用的概率分布函数
3 结果与讨论
3.1 突变检验与成因分析
采用Mann-Kendall非参数统计检验法在显著性水平α=0.05下对潦河万家埠站1953—2014年年径流时间序列进行趋势突变分析,结果(见图2)表明,万家埠站径流在1972年发生了突变. 突变前,万家埠站径流在1953—1972年间呈显著下降趋势,1972年后径流量迅速增加. 突变点前后多年平均径流深分别为910和1 020 mm,相比突变点前,突变点后的流域年径流量增加了12%. 采用Pettitt突变点检验法、SNHT方法、Buishand检验法和Buishand秩和突变检验法对万家埠站水文序列进行突变点诊断,均检测出1972年为突变点,进一步增加了突变点存在的可信度,该结果与肖丽英等[34]得到的变异检验结果一致.
图2 潦河流域年径流的变化趋势及突变点分析Fig.2 Trends analysis and break-point detection of annual runoff in Liaohe River Basin
潦河位于江西鄱阳湖流域,该区域自1950s以来,气候呈现出一定程度的波动性变化特征. 潦河流域降水量呈增加趋势,降水量在20世纪60年代初期和90年代初期发生突变,而蒸发量呈现下降趋势,夏季尤为明显[35]. 鄱阳湖流域20世纪60年代气温处于下降趋势,70年代初期气温有所回升,尤其在90年代后呈现显著增温时期. 潦河降水量的突变增加和蒸发量的突变减少以及气温的增加,导致径流在20世纪70年代初期发生突变,径流呈现增加趋势[36-37].
3.2 水文模拟与径流还原
根据潦河万家埠站水文序列突变检验结果,选择突变点前的1953—1971年为基准期,用于率定和验证DTVGM模型参数[38]. 考虑到部分雨量站数据缺测,选择1963—1968年(6年)作为模型率定期,1969—1971年(3年)作为模型验证期. 采用SCE-UA算法[39],以NSE及r作为目标函数,对所构建水文模型进行率定. 表2列出了模型最优参数,表3列出了模型在率定期和验证期的模拟精度评估,图3给出了模型的径流模拟结果.
表2 DTVGM模型参数意义及最优参数结果
表3 潦河流域模型率定期和验证期结果
图4 潦河流域1963—2014年月径流模拟过程Fig.4 Simulated monthly hydrograph during 1963-2014 in Liaohe River Basin
图3 潦河流域率定期和验证期月径流模拟过程Fig.3 Simulation of monthly runoff during calibration and validation period in Liaohe River Basin
从表3可以看出:针对日径流过程,模型率定期和验证期NSE分别为0.78和0.85,r分别为0.88和0.92;针对月径流过程,模型率定期和验证期NSE分别为0.85和0.91,r分别为0.92和0.96. 可见,率定期和验证期模拟径流和实测径流拟合程度较高,从而验证了DTVGM模型. 能够模拟潦河流域水资源的天然分布状况[40],为进一步开展生态流量研究奠定了坚实的基础. 固定DTVGM模型参数,输入1960—2014年流域降水量,将1960—1962年设置为模型预热期,即可获取基于水文模型的1963—2014年潦河流域的月径流过程(见图4).
3.3 生态流量计算结果
该研究采用GD、ND、GEV、P-Ⅲ及LN3五种分布对潦河流域的月径流过程进行适线分析,采用线性矩法获取各分布参数[41]. 由于不同分布针对同一样本系列(如7月径流序列)的拟合存在显著差异,因此需分析不同分布的拟合优度,对比不同分布的拟合效果,进而选择最优分布进行生态流量估算. 以多年7月径流序列为例,将5种分布及样本序列绘制在一起(见图5).
图5 潦河流域万家埠站7月径流序列5种概率分布函数拟合结果Fig.5 The fitting diagrams of 5 probability distribution functions of runoff in July at Wanjiabu Station, Liaohe River Basin
由于生态流量主要取90%保证率条件下的月径流过程,因此该研究重点聚焦于频率为0.5~1.0之间的适线效果. 从图5(b)可以明显看出,5种分布中,GEV、P-Ⅲ与LN3分布在低径流中拟合效果较好,GD分布虽常用于月降水的拟合[42],但拟合效果明显劣于GEV和P-Ⅲ分布. 5种分布中,ND分布拟合效果最差,这也从侧面证明了水文序列为偏态分布,与已有研究结果[43]相一致. 为了定量分析不同分布对月径流过程的拟合效果,该研究采用K-S检验、A-D检验及PPCC检验法对5种分布的适线性进行检验. 以潦河万家埠站7月径流序列为例,选用3种检验方法对5种分布进行拟合优度检验,结果如表4所示. 综合3种拟合优度检验方法的分析结果表明,GEV分布为7月径流序列的最优概率分布函数. 同时,绘制5种分布下理论值与观测值之间的相关关系(见图6),结果表明,P-Ⅲ、GEV与LN3分布对样本拟合效果最好,GD分布与ND分布拟合效果较差,这与表4拟合优度分析结果相一致.
表4 潦河流域万家埠站7月径流序列概率 分布函数拟合优度检验结果
选用5种分布对潦河流域逐月径流序列进行分析,采用3种检验方法对频率分布函数进行拟合优度综合检验. 检验分析表明,针对不同月份、不同检验指标其最优分布函数存在差异,如对万家埠站7月径流序列分布函数进行检验,其中K-S、A-D检验下GEV分布为最优,而PPCC检验法下LN3分布为最优(见表4). 因此,需综合分析多种方法对概率分析函数的拟合优度进行检验更为合理,且每个月份的水文序列相互独立,也需要分别进行概率分析. 采用相同方法对经还原的潦河万家埠站逐月径流序列进行拟合优度检验分析,即可获取相应水文序列的最优分布函数,结果如表5所示. 由表5可见,不同月份径流的最优分布不尽相同,其中GEV分布在6个月中表现为最优分布,GD、P-Ⅲ和LN3均在2个月中为最优分布. 选取逐月最优分布频率曲线,以90%保证率对应的流量作为生态基流,获取如表5所示的生态流量结果.
将最优分布的逐月频率计算法与Tennant法、最枯月平均流量法和7Q10法等常用方法进行对比,各方法所求生态流量及水量如表5所示. 其中Tennant法是较为传统且经典的水文学方法,该方法将一年分为两个阶段,根据鄱阳湖流域主要鱼类的繁殖规律,选定10月—翌年3月为一般用水期,4—9月为鱼类产卵期. 根据Tennant法对河流流量的描述,选定“非常好”对应的生态流量作为参照标准,即选用多年平均流量的50%作为鱼类产卵期的生态流量,一般用水期采用多年平均流量的30%作为河道内生态流量[44]. 最枯月平均流量法以河流最小月平均实测径流量的多年平均值作为河流的基本生态环境需水量. 7Q10法根据各年连续7日最低径流流量绘制理论频率关系分布曲线,但由于该标准对国内河流要求比较高,因此采用近10年最枯月平均流量作为生态流量. 对比表5中各种方法得到的生态流量可知,该研究采用最优分布的逐月频率计算法推求的年生态流量最大,为1 548.5×106m3;处于“非常好”水平的Tennant法计算的年生态流量次之,为1 282.7×106m3;7Q10法和最枯月平均流量法获取的年生态流量最低,分别为706.8×106和620.6×106m3. 从生态需水的季节性特征分析而言,最枯月平均流量法和7Q10法所给出的每月生态流量均相同,无法反映出不同月份生态流量之间的差异,其中最枯月平均流量法获取生态流量的Tennant法等级位于“差”的水平,而7Q10法获取生态流量的Tennant法等级位于“中等”水平. 基于“非常好”水平的Tennant法推求的生态需水过程,将一年分为汛期和枯水期,虽然考虑了生态需水在汛期和枯水期的差异,但仍未考虑水文变量的季节性变化特征,因此对于水资源量年内分布极度不均的流域仍需进一步分析[45]. 该研究采用最优分布的逐月频率计算法进行推求生态流量,根据逐月90%保证率条件下径流作为月生态流量,所获取的生态流量能够考虑生态需水的动态变化,该方法所求得年生态流量与“非常好”等级下Tennant法计算所得年生态流量相近,且远大于其他两种水文学方法得到的年生态流量,表明该生态流量能够满足水生生物的正常需求,计算结果具有确定性和合理性[18,20-21].
图6 潦河流域万家埠站7月径流序列5种频率曲线理论值与观测值散点图Fig.6 The scatter plot of observed and theoretical values of 5 probability distribution functions for runoff in July at Wanjiabu Station, Liaohe River Basin
表5 不同方法下潦河流域生态流量计算结果比较
4 结论
a) 采用多种水文序列突变检验方法对潦河万家埠站进行检验,发现潦河万家埠站径流在1972年发生变异,主要是由于降水量的突变增加和蒸发量的突变减少,导致径流在20世纪70年代初期发生突变,径流量呈现增加趋势.
b) 基于突变点前基准期对DTVGM模型进行参数率定和验证,率定期和验证期模拟径流和实测径流拟合程度较高,从而验证了DTVGM模型能够较好地模拟潦河流域水资源的天然分布状况.
c) 基于多种拟合优度指标对分布函数进行分析,选用最优分布函数开展生态流量计算,生态流量能够满足水生生物的正常需求. 与Tennant法、最枯月平均流量法、7Q10法等方法比较,基于最优分布函数的生态流量结果更具确定性与合理性,可为变化环境下水资源的适应性管理提供重要依据.