摘 要：
【目的】为科学探究流域水文情势变化驱动因子，识别受损河流生态退化驱动机制，【方法】以长江一级支流大宁河为研究区域，基于巫溪水文站1972—2016年长序列水文、土地利用、水利工程建设等数据，利用过程线、一维动态豪斯道夫分维数、Mann-Kendall检验、有序聚类、滑动T检验、滑动秩和检验等数理统计方法进行变异点分析及确定，结合实地调研、查阅当地县志等资料，通过对降雨径流关系、土地利用类型时空变化、小水电建设特征要素的分析论证。【结果】明确了区域水文一致性在1983年、2005年发生变异，区域水文变异驱动因子为小水电的兴建及运行，2005年后新建的小水电可占区域小水电数量的68.23%;基于IHA软件定量分析变异前后环境流量组分变化，发现大洪水组分变化最为剧烈，现状(变异后)大洪水出现次数为0,月低流量组分变化较小，改变度均为超过25,【结论】从而为后续变化环境下水文-生态响应关系、以河湖生态复苏为导向的水资源调配与管理的相关研究提供理论支撑。
关键词：
水文变异;驱动因子;流量组分;大宁河;小水电;土地利用;
作者简介：
彭娟莹(1981—),女，高级工程师，硕士，主要研究环境污染物检测与治理。
*于子铖(1993—),男，讲师，博士，主要从事于河湖生态保护与修复的研究。
基金：
湖南省科技重大专项(2018SK1010);
国家重点研发计划(2016YFC0401606,2016YFC0401408);
引用：
彭娟莹, 于子铖, 付意成. 大宁河巫溪段水文变异归因及其环境流量组分变化分析[J]. 水利水电技术(中英文), 2023, 54(9): 86- 98.
PENG Juanying, YU Zicheng, FU Yicheng. Attribution of hydrological variability in the Wuxi section of the Daning River and analysis of its environmental flow component changes[J]. Water Resources and Hydropower Engineering, 2023, 54(9): 86- 98.
0 引 言
我国幅员辽阔，区域水文水资源问题复杂多样。在变化环境下，水文系统的变化一般是由自然(气候、自然地理条件)变化、人类扰动共同作用而导致，近些年随着社会经济的快速发展，人类扰动通常被认为是引起变异的主要驱动因子。天然的水文循环受到人类扰动，原有的自然规律被破坏，从而致使水文系统发生改变，其中一个最大的表现就是天然水文时间序列的一致性被破坏。水文时间序列在某个时间节点前后，统计规律会发生显著的变异，这一变化时刻称为水文过程变异点。关于变异点诊断的方法很多，如Mann-Kendall检验法、有序聚类法、滑动秩和检验法、Brown-Forsythe检验法、游程检验法、滑动T检验法、累积距平法、滑动F检验法、Pettitt法、Lee-Heghinan法、Bayesian变点分析法、R/S分析法与Hurst指数等。为避免单一方法的局限性，谢平对水文时间序列展开了初步、详细、综合的变异分析，初步分析是对水文序列的随机性进行分析，详细分析为采取多种方法推求可能存在的变异点，最后通过综合对分析进行确定;孔波等采用Mann-Kendall变异诊断法和有序聚类法，诊断了北洛河流域枯汛交替期、汛期、汛枯交替期、枯期等分期径流系列变异点并划分了变异等级；陈隆吉等借用信息熵构建洪水过程均匀度模型，然后采用水文变异诊断系统分析洪水过程均匀度序列的变异规律；徐淑琴等利用水文变异诊断分析水文时间序列，确定人类不同扰动下的时间节点，并考虑东北地区具有冰封期的特点，推求了混合时间尺度下的天然、次天然及现状生态径流过程；为增加变异检验的可信程度，刘剑宇等采用有序聚类、Mann-Kendall检验、累积距平等8种方法对水文变异点进行分析。总体上，关于水文变异诊断方法的研究呈现多元化与多样化的趋势。
水文过程在人类、生态系统之间起着至关重要的连载作用，人类扰动水文过程，进而影响生态系统。反之，生态系统的变异也会影响水文情势的变化，进而使水文规律发生变异，偏离天然水文状态，从而影响人类社会。要实现人水和谐以及平衡三者之间的互馈关系，需了解其作用机理，明确其水文变异原因。引起水文变异的因子一般包括气象条件、水利工程、社会总耗水量、水土保持、土地利用等，不同区域需具体分析。流量组分是联系水文过程和生态过程的重要概念，描述了水文情势的不同部分与生态系统的相关性。常俊超等基于IHA/RVA方法，探究变化环境下黄河中游环境流量组分与水文情势变化情况；管新建等利用IHA/RVA方法进行了环境流量组分与水文生态特征分析，揭示了在水文变异后流量年平均值有所增加，低流量和高脉冲环境流量成分占比较大；张远等利用FLOWS法，以鱼类为保护目标，构建了流量组分与鱼类生态需求关系模型，计算了包含基流、脉冲流、平滩流和漫滩流4种组分在内的河流环境流量。
本文以大宁河巫溪水文站控制流域为研究区域，以长序列水文数据为基础，依据初步、详细、综合分析的路径，利用过程线法、Mann-Kendall检验法、有序聚类、滑动秩和检验等数理统计方法，科学确定水文变异点，结合区域降雨、土地利用、水利工程建设等数据，统筹考虑降雨-径流关系、土地利用时空分布变化、小水电数量与累积装机容量变化趋势，利用定性与定量相结合的手段探究水文变异的主要驱动因子，基于IHA软件定量分析变异前后环境流量组分变化。以期为区域水资源规划及管理、水生态系统恢复、水环境持续改善、水文预测等提供相关依据。
1 研究区域概况
大宁河地处东经108°44′—109°59′,北纬30°14′—31°44′之间，位于重庆市巫溪、巫山两县境内，全长为162.0 km, 流域面积为4 181.0 km2。本次研究区域为大宁河巫溪段(见图1)。巫溪县位于长江上游，属重庆市管辖，在其东北部。据统计，从每年3月的中下旬开始，降雨量会逐渐增加，径流也随之呈增大的趋势，4月是汛前过渡期，在5—9月会进入主汛期，径流量会出现大幅增加。本文依据向当地相关部门收集到的巫溪水文站1972—2016年水文径流、降雨数据，诊断分析大宁河巫溪段水文变异状况。巫溪水文站于1972年设站运行，位于巫溪县的城厢镇、大宁河干流，是大宁河中上游来水的重要控制站，控制区域面积约为2 440.15 km2,水系包括东溪河、西溪河、后溪河等。
图1 研究区域
2 研究方法
2.1 水文变异分析方法
2.1.1 初步分析
初步分析采用过程线法与一维动态豪斯道夫分维数法。过程线法依据目估，从而定性观察水文序列随时间变化的波动状况，以此初步判断序列的随机性；一维动态豪斯道夫分维数法源自分形理论，分形理论主要研究自然界的不规则现象及其内在规律，水文过程具有非线性与随机性，随机的径流流量过程表现出非线性的分形理论特征。
2.1.2 详细分析
推求水文变异的方法有很多，考虑单一检测方法存在局限性，本文基于对研究数据的初步分析，结合水文时间序列趋势与突变分析系统,筛选出Mann-Kendall检验法、有序聚类法、滑动T检验法、滑动秩和检验法进行巫溪水文站上游水文变异点进行详细分析。其中Mann-Kendall检验法为非参数检验方法，不需要样本遵从一定的分布，也不受少数异常值的干扰；有序聚类法旨在寻求左右分割点，需进一步对分割样本进行检验；滑动T检验法计算简单，考察两组样本均值的差异显著性；滑动秩和检验法不需要推求总体参数，对总体分布也没有特别要求。所选方法判断核心均为跳跃检验，对水文序列要求不高，具有较强的适用性。
Mann-Kendall(以下简称M-K法)检验法是一种非参数统计检验方法，对于具有n个样本量的时间序列x,构造一组秩序列,秩序列Sk表示的是第i时刻的数值大于j时刻的累计数。在时间序列随机独立的假定下，定义统计量UFk。UFi代表标准正态分布，在按照时间序列x的顺序x1,x2,…,xn下，推求统计量序列，结合显著性水平α下，查阅正态分布表，如果UFi>Uα,则表明序列存在显著变化。按x逆序xn,xn-1,…,x1,再次重复上述过程，与此同时令UBk=-UFk(k=n,n-1,…,1),UB1=0。如果UFk或UBk的值均大于0,那么可以说明序列呈上升的趋势，反之则说明呈现下降的趋势，超过临界直线(如α=0.5,则临界值为u0.05)时，说明变化趋势明显。超过临界曲线的范围确定为出现突变的时间范围，UFk和UBk两曲线的交点(在临界值内)为序列突变点；有序聚类法的实质为寻求序列的最优分割点，原则是保证同类型之间的离差平方和较小;滑动T检验法通过查验两组样本均值的差异是否显著来判定突变点，设x为具有n样本量的时间序列，基于实际分析，人为设置一个基准点，将x分为具有样本n1的x1、n2的x2,从而定义统计量T;滑动秩和检验法是水文序列变异诊断的主要方法,首先设定一个分割点τ0,建立τ0前后两个序列的分布函数F1(x)、F2(x),然后从两个序列中依次抽取样本，样本容量为n1、n2,令检验的原假设F1(x)=F2(x),再次确定n1、n2样本中的秩，最后，令n1、n2的秩和为T,然后对统计量T进行检验，当n1、n2-n1超出查表范围，可以用近似U检验。
2.2 变异归因分析
大宁河位于三峡库区，但大宁河巫溪段，尤其是巫溪水文站上游流域，水文情势受到三峡电站调节运行的影响可以忽略，不作为该流域水文变异考虑因素。结合2018—2021年对研究区域的实地调研、询问当地相关人员、查询当地县志等相关资料，研究区域气温、社会总耗水量、水土保持近些年均为发生较大改变，因而以下从降雨、水利工程、土地利用3方面进行水文变异驱动因子分析。
径流量受流域降雨与下垫面的共同影响，通常在无人为影响的阶段，径流量与降雨的变化趋势一致，具有较好的相关性。因此，首先对降雨数据进行突变检验；其次，建立区域降雨径流关系，明晰是否存在人为影响；最后，分析土地利用、水利工程建设变化是否与水文变异点一致，从而确定区域水文变异主要驱动因子。
2.3 环境流量组分变化分析
IHA将环境流组分分为34个指标(EFC),其释义与生态学意义如表1所列。利用IHA软对所有日流量数据按升序排频，取前50%为低流量事件，最小的20%被划分为月枯水流量，75%以后的流量为高流量，结合我国对洪水量级的划分，按不同重现期分别划分为小洪水、大洪水。结合不同指标对应的生态学意义，对不同阶段的环境流量组分进行分析，进而定量研究变异驱动因子对其的影响。
3 结果分析与讨论
3.1 水文变异点分析
3.1.1 水文变异初步分析
基于1972—2016年长序列水文数据，形成平均流量变化趋势图(见图2)。平均流量随时间变化整体呈减少趋势，在1983年取得该序列的最大值，为129.45 m3/s, 在2006年取得最小值，为34.36 m3/s, 从而初步判断1983年、2006年可能为水文序列发生变异的年份。
图2 平均流量变化趋势
依据基础数据结合相关公式，生成年径流量的一维动态豪斯道夫分维数趋势变化(见图3),起伏变化的趋势初步可以判断水文情势发生了突变。自1997年后Dft出现多出高值，初步分析可能在1997年前后发生了突变。
图3 年径流量的一维动态豪斯道夫分维数Dft趋势变化
3.1.2 水文变异详细分析
基于1972 —2016年径流量数据，生成巫溪水文站控制流域M-K曲线(见图4)。由UF曲线可见，自1985年前后，年径流量呈现整体下降的趋势，曲线大都在临界线之间，表明下降趋势不太显著。依据UFk和UBk两曲线的交点，结合UF曲线整体趋势，确定1977年、1997年、1998年、2001年、2003年、2005年、2007年、2009年为可能变异点。
图4 巫溪水文站上游年径流量M-K统计量曲线
依据有序聚类法计算公式得到巫溪水文站年径流量有序聚类结果变化趋势(见图5)。由图5知，1983年离差平方和最小，故而判定1983年为水文可能变异点。
图5 有序聚类结果趋势
依据滑动T检验法的计算公式，得到统计量T的变化趋势，最大值在2012年取得(见图6)。此外，在1984年、1994年、2006年存在局部最大值，因为判定1984年、1994年、2006年、2012年为可能变异点。
图6 滑动T检验方法结果趋势
依据滑动秩和检验的相关公式对长序列各点进行检验，结合巫溪水文站1972—2016年径流数据，得到表2,确定1983年、2005年为可能变异点。
3.1.3 水文变异点确定
将详细分析结果进行汇总(见表3),综合确定巫溪水文站上游水文变异点。由表3可以直观看出，2005年出现次数最多，可以判断为水文变异点；1983年、1984年出现次数相同，结合水文变异初步分析结果、累积径流量变化趋势，确定1983年为水文变异点，因而将水文序列分为1972—1982年、1983—2004年、2005—2016年三个序列。
3.2 变异驱动因子分析
3.2.1 基于降雨径流关系的驱动因子分析
以巫溪县1972—2016年巫溪气象站(测站编码57345)逐月降雨数据作为依据，分析自然变化对水文情势变异的影响，测站来源于中国气象数据网。结合图7,发现随着时间的变化，降雨整体呈现下降趋势，在1983年取得最大值，2001年取得最小值。
图7 年降雨变化趋势
基于降雨资料，结合M-K法相关计算公式，得到UF、UB统计量变化趋势(见图8)。由图8知，两统计量曲线存在交点，但其均未超过临界曲线，因而认为降雨未发生显著突变。结合多年径流量数据，得到多年径流深。在降雨无显著突变的背景下，建立降雨-径流双累积曲线，分别为y=1.154 7x+94.108 0(R2=0.998 5)、y=1.147 5x+1 793.800 0(R2=0.999 5)、y=1.005 1x+6 901.300 0(R2=0.996 6),具体如图9所示。由图9知，在1983年、2005年前后，曲线两次发生了偏离，以1972—1982年作为基准期，时间节点与上述所求变异点一致，从而也说明了两个阶段的水文变异均受到人类活动的影响。结合图8知，降雨未发生显著突变，由此可以判断两个阶段的水文变异驱动因子主要源自人为影响。依据研究区域实测1972—1982年、1983—2004年、2005—2016年径流资料，得到三个时间段年径流深均值R¯1、R¯2、R¯3为1 497.2 mm、1 443.1 mm、1 224.3 mm,由此可以算出与第一阶段径流深的差值ΔR2-1、ΔR3-1为-54.1 mm、-272.9 mm,2005—2016年径流变化明显强于前一个序列。
图8 M-K降雨突变检验
图9 降雨-径流累积双曲线
3.2.2 土地利用时空分析
结合研究区域概况，人类活动影响主要从土地利用和水利工程两个方面进行研究与分析。大宁河流域巫溪段1980年(20世纪70年代末)、1990年(20世纪80年代末期)、2000年(20世纪90年代末期)、2010年(21世纪初)、2015年的5期土地利用数据(30 m)从中国科学院资源环境科学数据中心(http: //www.resdc.cn/)获得。利用Arcgis进行巫溪县土地利用类型处理与分析。
结合图10可以看出：(1)1980—2015年土地利用/植被覆盖整体变化不太明显；(2)发现林地占比最大，主要集中在东北部，其次为耕地，主要集中在南部；(3)草地变化较为明显，西北部分区域由最初林地、耕地转化为草地；(4)城乡、工矿、居民用地主要集中在巫溪县城附近；(5)水域面积有所提升，主要表现在西溪河中上游区域。
图10 2015年大宁河流域巫溪段土地利用类型分布
由表4和图11可以直观看出：(1)5个时期土地利用面积占比大小一致，依次为林地、耕地、草地、水域、建设用地(城乡、工矿、居民用地)、未利用土地，林地占比均在58.00%以上，未利用土地较少，占比均在0.01%以下；(2)耕地面积逐年下降，由最初的27.90%下降至26.53%,下降了1.37%;(3)林地整体变化较小，呈现先下降后上升的趋势，变化了0.79%;(4)水域在2010年后有所增加，增加了0.05%;(5)随着社会经济的快速发展，城乡、工矿、居民用地2000—2010年期间变化较大，变为原来2倍；(6)1980—1990年，各土地利用类型变化均不大，由此可以判断，1983年的变异点与土地利用类型变化无相关关系；(7)以2000年土地利用类型占比为基础，2010年林地、耕地、草地等6个类型的变化率依次为3.43%、1.25%、0.65%、5.78%、92.06%、71.43%,占比较大的耕地、林地(两者相加占85%以上)变化较小，城乡、工矿、居民用地变化较大主要表现在巫溪县城附近用地的变化，且柏杨河汇入口在巫溪水文站下游，故认为此阶段土地利用类型变化可能会引起局部区域水文的变异，但对于确定的巫溪水文站上游2005年水文变异点无显著贡献。
图11 1980—2005年大宁河流域巫溪段土地利用类型对比
3.2.3 水利工程建设
巫溪县地形以山地为主，山地占93%,属于典型的中深切割中山地形，最低海拔139.4 m, 最高海拔2 796.8 m。起伏变化的地貌使其在水电站建设上具有得天独厚的优势，据调研统计，截止到2020年，巫溪县有130座小水电，1座中型水电站。为研究小水电对水文情势变化的影响，梳理统计了水文站控制流域内小水电(均为引水式)空间位置、建站时间、装机容量等相关信息。
针对小水电空间位置而言(见图12),大宁河(两河口—巫溪水文站)共建有11座，挡水坝均在支流上；西溪河流域共建有30座(其中中梁一级电站为中型电站),挡水坝8座在干流，其余均在支流；东溪河流域共建有26座，挡水坝9座在干流，其余均在支流；后溪河流域建有19座，挡水坝8座在干流，其余均在支流，巫溪水文站控制流域内共计小水电85座，中型水电站1座，其开发方式均为引水式。
图12 小水电空间分布示意
结合图13与图14,发现随着时间的变化，巫溪水文站控制流域内小水电数量、累积装机容量呈现一个明显递增趋势，在1977年、1983年、2004年存在拐点，尤其是2004年前后，有较大变动。1972—1977年，修建小水电6座，累积装机容量为6 505 kW;1978—1983年，修建小水电11座，累积装机容量新增19 575 kW;1984—2004年，修建小水电10座，累积装机容量新增13 810 kW;2005—2016年，修建小水电58座，占区域小水电总和的68.23%,累积装机容量新增304 065 kW,占区域累积装机容量的88.40%,该阶段小水电得到快速的发展。随时间变化小水电数量、累积装机容量与确定的水文变异点吻合，结合2018年、2019年、2021年在研究区域的实地勘察与梳理当地相关资料，可确定引起该流域水文变异的主要驱动因素是小水电的兴建及运行。从而将1972—1982年、1983—2004年、2005—2016年三阶段分为天然状态、弱干扰阶段、强干扰阶段。
图13 随时间小水电个数变化
图14 随时间小水电累积装机容量
3.3 环境流量组分变化分析
利用IHA软件进行计算分析，具体结果如表5所列。表中所提改变度为(弱干扰/强干扰阶段数值—天然状态数值)/天然状态数值，结合表5,发现大洪水组变化最为剧烈，月低流量组变化较小。
对比天然状态(1972—1982年)、弱干扰阶段(1983—2004年)环境流量组分差异，由图15知，变化程度前五的分别为E34(大洪水下降率)、E33(大洪水上升率)、E29(大洪水极大值)、E32(大洪水频次)、E16(特枯流量频次),主要影响组别为大洪水、特枯流量，对河岸区域的物种平衡、河漫滩的物理生境、河道的自然形态结构以及自然栖息地条件形成一定的影响。
图15 天然-弱干扰环境流量组分改变度排序
利用IHA计算天然状态(1972—1982年)、强干扰阶段(2005—2016年)环境流量组分差异，具体结果如表5所列、如图16所示。对比天然状态，强干扰阶段不存在大洪水，猜测原因是小水电的修建及调蓄。不考虑大洪水的情况下，改变程度前五的指标分别为E16(特枯流量频次)、E15(特枯流量发生时间)、E27(小洪水上升率)、E14(特枯流量持续时间)、E24(小洪水持续时间),分别属于特枯流量、小洪水组别，对河岸物种平衡、河漫滩物理生境、河流自然演变、水生动植物生活生存区域、河道自然形态结构、动植物群落分布均造成一定影响。
图16 天然-强干扰环境流量组分改变度排序
4 结 论
(1)以大宁河巫溪段为研究区域，基于长序列水文监测数据，利用多种数理统计方法进行水文变异点的初步、详细分析，避免了单一方法的局限性和偶然性，确定了研究区水文变异点为1983年、2005年。
(2)建立区域的降雨—径流关系曲线，明确了水文变异驱动因子主要源自人为影响。结合确定结果与现场踏勘观测，详细分析了土地利用、小水电建设特征随时间的变化，将变化的突变点与水文变异点进行对比，确定了小水电密集建设是区域水文变异的主要驱动因素，定量分析了变异前后环境流量组分的变化。
(3)以多方法综合分析水文变异点—降雨径流关系—对比驱动因子时间变化突变点为路径，为科学诊断区域水文变异原因提供了技术支撑。基于确定的水文变异点，可以明晰变异前天然水文序列的节律变化，了解河流的自然动态机制，可为以河湖生态复苏为导向的水资源调配与管理提供决策性支持。
水利水电技术（中英文）
水利部《水利水电技术（中英文）》杂志是中国水利水电行业的综合性技术期刊（月刊），为全国中文核心期刊，面向国内外公开发行。本刊以介绍我国水资源的开发、利用、治理、配置、节约和保护，以及水利水电工程的勘测、设计、施工、运行管理和科学研究等方面的技术经验为主，同时也报道国外的先进技术。期刊主要栏目有:水文水资源、水工建筑、工程施工、工程基础、水力学、机电技术、泥沙研究、水环境与水生态、运行管理、试验研究、工程地质、金属结构、水利经济、水利规划、防汛抗旱、建设管理、新能源、城市水利、农村水利、水土保持、水库移民、水利现代化、国际水利等。