黄河为中国第二条大河,被誉为母亲河,受多变的水资源条件和剧烈的人为活动影响,自20世纪90年代以来,黄河干流实际来水量不断减少,且下游出现连年断流现象,最严重的是1997年黄河入海口断流时间达到了226天[1]。与此同时,黄河在流经黄土高原时会携带大量泥沙,导致河道淤塞,洪涝灾害事件频繁发生[2]。为了防治水土流失,1999年我国投资了约870亿人民币在全国范围内实施了退耕还林还草工程,生态环境脆弱的黄土高原成为首批试点区域之一[3,4,5]。随着工程的实施,1999-2010年黄土高原累计造林面积已达1.89×107hm2,且该区域植被生长季NDVI以0.15/年的速率增加[6,7]。得益于植被的水土保持效益,截止到2015年黄河河道径流量减少约26%,输沙量减少约为21%,黄土高原地区脆弱的生态环境得到修复[8]。随着植被覆盖面积持续增加导致径流量持续下降,已有报道[4]发现,该地区的植被恢复已接近阈值,该区域将无法实现生态环境用水和人类生活用水之间的平衡发展。因此,定性并定量解析某时间段内气候变化和人类活动对径流变化的影响程度是实现流域综合治理的关键问题。

气候及土地利用等自然-人为因素对流域径流产生显著影响[9,10],且因其存在时空变化,使径流的变化也存在时空差异[11]。Liu等[12]以全国流域为研究对象,定量解析气候变化与土地利用对径流量的影响,结果表明:在1980—2000年期间,气候变化对径流变化起主导作用;2001—2014年,人类活动对径流变化的影响日渐显著,气候变化和人类活动对径流变化的贡献率分别为53.5%和46.5%,但在黄河流域发现人类活动对径流的贡献率(69.7%)显著大于气候变化(30.3%)。李志等[13]以黄土高塬沟壑区的黑河流域(泾河支流)为研究对象,发现人类活动是驱动径流量减少的主要因素(贡献率为76%)。师忱等[14]以滦河流域为研究对象,结果表明,人类活动是影响径流减少的主要因素,其贡献率为57.26%,远大于气候变化的影响。还有相关学者[15]以黄土高原清水河流域为研究对象,发现土地利用变化对年径流减少的贡献率(53.21%)高于气候变化(46.79%),对流域径流变化原因深入分析表明,土地利用变化主要是因乔木林地面积增加和潜在蒸发散增加造成的。然而,对黄土高原吕二沟流域的研究[16]结果显示,人类活动对生态水文影响作用逐渐减少,而气候变化对其影响逐渐加剧。对黄土高原径流量变化归因分析发现,植被恢复对径流量的减少起着关键性的作用[14,15]。植被层表面通过将降水分配为林内降雨量、树干流量和冠层截留等过程,地下根系则通过改变土壤理化性质,增强土壤的入渗能力等影响产流[17]。综上所述,在流域尺度、气候变化和人类活动对径流量变化的比重并不相同,需采用定性与定量相结合的方法开展径流量的变化成因研究[18]。因此,本文基于黄土高原昕水河流域多年水文气象资料,采用趋势分析、双累积曲线法以及生态水文分析法,围绕流域尺度水量平衡,研究气候变化和土地利用对黄土高原中部昕水河流域径流量的影响,分析径流变化趋势并定量解析气候变化和土地利用对径流变化的贡献率,为合理开展流域综合治理、生态恢复以及水土资源管理提供科学依据。

1、材料与方法

1.1研究区域概况

昕水河为黄河的一级支流,发源于山西省临汾市蒲县摩天岭,地理位置为110°30′—111°27′E,36°05′—36°55′N,干流全长约为178km,面积为3992km2,高程为699～2010m(图1)。

图1昕水河流域地理位置、地形、气象站及水文站

昕水河流域属于黄土高原丘陵沟壑区,地形破碎,沟壑纵横,呈明显的条带型,主要土壤类型为褐土,水土流失严重。该流域为暖温带半湿润性大陆性季风气候,年均气温在7～10℃,年均降水量约为533mm,主要集中在夏季。流域植被类型以松栎林、灌木草本和经济林为主,植被包括栎树,油松,刺槐等。本研究数据资料来源见表1,其中水文数据来自黄河水利委员会的1958—2015年逐日径流数据,气象数据来自中国气象科学数据共享服务网1958—2015年日序列数据,DEM数据源自地理空间数据云。其中,流域面降雨量采用泰森多边形进行计算[19],潜在蒸散发选择应用最广的彭曼公式进行计算[14]。

表1昕水河流域数据资料及其来源

1.2研究方法

1.2.1趋势分析及突变点

采用Mann-Kendall(M-K)非参数检验方法分析昕水河流域1958—2015年径流量、降水量和潜在蒸发散的变化趋势,并初步确定径流的突变点。M-K非参数突变检验作为常规非参数检验方法,被广泛应用于气象水文时间序列的趋势分析。其中,

式中:n为时间序列个数;S为检验统计值;Z为M-K统计值。当Z>0时,则时间序列呈上升趋势,Z<0时,则时间序列呈下降趋势;当Z的绝对值大于等于1.28,1.64,2.32时,即通过置信区间90%,95%,99%的显著性检验,具体计算公式见文献[20]。M-K非参数突变检验是检验水文时间序列突变点的方法,该方法首先计算时间序列的正序列的统计值UF和逆序列的统计值UB,当统计值UF和UB序列曲线出现交点时,且交点在临界线之间时,此交点所在的年份即为可能发生突变的时间[20]。

双累积曲线法是检验2个参数间关系是否一致性及其变化的常用方法[21]。其原理是2个变量按同一时间长度逐步累加,一个变量的累积值为横坐标,另一个变量的累积值作为纵坐标,可以生成一条直线,其斜率为2个变量对应点的比例常数。若双累积曲线的斜率发生突变则意味着2个变量之间发生了改变,那么斜率发生突变点所对应的年份就是2个变量累积关系出现突变的时间[21]。本文采用修订后的双累积曲线法确定流域径流量的突变点。

1.2.2土地利用

本文选取了1986年、1995年、2006年的TM遥感影像和2015年Landsat8遥感影像,运用ENVI5.1和ARCGIS10.2等软件进行土地利用信息的获取,结合昕水河流域土地利用现状,将流域划分为林地、灌木林地、草地、耕地、建筑用地及居民区和水体6种土地利用类型。本文采用土地利用变化率(K)来研究某一时段内土地利用类型面积的变化规律[22],计算公式为:

K=100×(Ub−Ua)Ua×100%(2)

式中:Ua和Ub分别为研究期初和研究期末某一土地利用类型的面积(km2)。

1.2.3气候变化和土地利用对径流变化的定性分析

生态水文分析方法简称TS法[23],是一种概念性模型,主要通过利用流域内水-热收支的耦合关系鉴定流域内植被与气候变化对流域水循环的影响。计算公式为:

Pex=(P-ET)/P(3)

Eex=(PET-ET)/PET(4)

式中:P为降水量(mm);PET为潜在蒸发散量(mm);ET为实际蒸发散量(mm);Pex和Eex代表流域中剩余或是未利用的水和能源。对于正在经历植物群落演变的原生流域生态系统,TS法的原理是假定一种在环境中该群落尽量减少未使用的水和能源的趋势。因此,在该生态系统中任何人为活动的连续的变化和气候变化都会使生态水文的性质发生变化,具体以Pex、Eex的形式体现。然而,植被变化或者管理措施的实施通常改变的是实际蒸发散而不是降雨和潜在蒸发散,因此生态水文对土地利用的变化表现为:(1)Pex、Eex同时增加,与森林砍伐有关;(2)Pex、Eex同时减少,与植被恢复与重建有关。而对于气候变化:(1)Pex增加、Eex减少,由于降雨与潜在蒸发散的比率随时间上升;(2)Pex减少、Eex增加,由于降雨与潜在蒸发散的比率随时间下降见图2。由于气候变化和土地利用变化使Pex和Eex产生不同程度的改变,从而使土地利用变化和气候变化的影响得以区分[23]。

图2生态水文法中气候和土地利用变化下的生态水文变化的概念模型[23]

1.2.4气候变化和土地利用对径流变化的定量分析

采用水量平衡原理方法定量评估气候变化和土地利用变化对径流变化的影响[16]。具体公式为:

ΔQTOT=ΔQC+ΔQU(5)

式中:ΔQC-是气候原因引起的年均径流量的变化(mm);ΔQU-是土地利用变化原因引起的年均径流量的变化(mm);且ΔQTOT对于一个流域可认为是由径流突变点前后2个阶段的年均径流之差[15](mm),计算公式为:

ΔQTOT=QAC−QBC(6)

式中:QAC-为突变点后年均径流量(mm);QBC-¯为突变点前年均径流量(mm);ΔQC-为假设土地利用不发生变化由气候变化引起的年均径流量变化量[24](mm),计算公式为:

ΔQC=βΔP+γΔPET(7)

式中:ΔP和ΔPET是降雨和潜在蒸发散的变化量(mm);β、γ为径流响应于降雨和潜在蒸发散变化量的灵敏系数。其计算公式为:

β=1+2x+3wx2(1+x+wx2)2      (8)    

γ=−1+2wx(1+x+wx2)2        (9)

式中:x为干旱指数,x=PET/P;w为与下垫面条件有关的一个综合参数,与地形、植被和土壤等的因素有关,考虑到在研究区的58年时间内,地形和土壤都不会发生大的变化,故植被变化是影响w的主要因素,且w是各地类wi的面积加权和[15,16]。

2、结果与分析

2.1土地利用变化

由表2可知,1986年流域土地利用类型以耕地、草地、灌木林地和林地为主,四者的比例分别为29.30%,40.82%,22.66%,7.04%。到1995年,建筑用地和居民区、灌木林地和林地的面积较1986年分别增加16.12%,15.6%,6.28%,而耕地和草地的面积分别减少3.40%,7.39%。较之1995年,2006年建筑用地及居民区、草地、灌木林地和林地分别增长11.77%,6.20%,12.93%,22.32%,而耕地呈现大幅度下降(26.17%)。到2015年,建筑用地及居民区、草地、灌木林地、林地较2006年分别增加32.11%,2.19%,2.00%,3.01%,耕地减少8.66%。相比于1985年,2015年灌木林地、林地、草地、建筑用地及居民区分别增加33.15%,33.91%,0.51%,71.46%,而耕地减少34.86%。总体来说,昕水河流域林地、灌木林地和建筑用地及居民区显著增加,耕地明显减少,下垫面条件得到很大的改善。

表2昕水河流域不同时期土地利用对比

2.2突变检验

M-K非参数检验结果表明,1958—2015年,流域径流量呈显著下降的变化趋势(Z=-5.84,p<0.0001),且平均每年减少0.69mm。但降水量(Z=-0.72,p=0.31)和潜在蒸发散(Z=-0.50,p=0.88)没有显著趋势变化(图3)。突变检验分析结果显示:UB和UF有且只有1个交点(1980年左右),但在临界线范围外(图4),因此需要采用双累积曲线法进一步确定径流量的突变点。降水-径流双累积曲线分析结果表明(图5),1974年和2000年的直线斜率发生明显偏折,由此判断1974年和2000年为径流变化突变点。所以将昕水河流域径流研究时序划分为:1958-1974年、1975—2000年、2001—2015年3个阶段。其中,第1阶段1958—1974年可作为昕水河流域径流序列的基准期。

2.3气候变化和土地利用对流域径流影响的定性分析

在生态水文分析法中,气候变化通过降水量和潜在蒸发散来影响径流量,而土地利用变化是以实际蒸发散的形式来影响径流量。若水文过程中的剩余降水量和剩余能量值发生明显改变时对实际蒸发散有显著影响,其中,由图6a可知,流域1958—2015年逐年剩余的水能量数据点分布状况。图6b中3个点分别代表基准期(1958—1974)、第2阶段(1975—2000)和第3阶段(2001—2015)Pex和Eex的平均值,可以得出从基准期到第2阶段的走向中Pex呈下降趋势,Eex呈上升趋势,这段时间内剩余的水减少,剩余的能量增加,因此气候变化为第2阶段年间径流变化量的主导因素。从基准期到第3阶段的走向为Pex和Eex都呈下降趋势,第3阶段剩余的水和剩余的能量都呈减少的趋势,因此土地利用变化为第3阶段径流变化的主导因素。

图3年径流深、降水量和潜在蒸发散变化及其线性趋势

2.4气候变化和土地利用对流域径流影响的定量分析

基于昕水河流域年径流量突变点,将整个时间序列分成3个阶段,并比较这3个阶段内气候因素和径流量的平均值、标准差和变异系数(表3)。由表3可知,第2阶段年均降水量、标准差、变异系数分别为493.45mm,99.87,0.20,比基准期(第1阶段)分别减少10.74%,28.80%,16.67%;第3阶段年均降水量、标准差、变异系数分别为553.64mm,108.84,0.20,比基准期分别增加0.15%,减少22.41%,16.67%。第2阶段年均潜在蒸发散1041.6mm,50.68,0.05,比基准期分别减少0.99%,20.29%,16.67%;第3阶段年均潜在蒸发散1040.47mm,41.34,0.04,比基准期分别减少0.43%,34.98%,33.33%。

图4年径流量M-K突变检验

注:UF为其序列正序M-K检验的统计值变化曲线;UB为其序列逆序M-K检验的统计值变化曲线;虚线为临界线。

依据下垫面土地利用类型的改变求得研究区内w参数结果为0.88。基于公式(8)和公式(9)求得昕水河流域径流对降水的敏感性系数β为0.35,对潜在蒸发散的敏感性系数γ为-0.10。敏感性结果表明,降水每增加1mm,流域径流将增加0.35mm;潜在蒸发散每增加1mm,流域径流将减少0.10mm。

图5昕水河流域年有效降雨-径流双累积曲线

图61958-2015年Pex和Eex的散点图

基于水量平衡原理计算气候变化和土地利用对径流的影响。由表4可知,昕水河流域第2阶段和第3阶段年均径流量分别为26.94,16.83mm,与基准期相比分别减少26.48,36.59mm。在第2阶段,气候和土地利用变化分别使径流量减少19.37,7.11mm,气候因素对流域径流量减少的贡献率为73.14%,其中降水和潜在蒸发散对径流量减少的贡献率分别为76.91%和-3.78%,而土地利用变化对流域径流量减少的贡献率为26.86%。在第3阶段,气候变化使径流量增加了1.4mm(其中降水使径流量增加0.29mm,潜在蒸发散使径流量增加1.11mm)和土地利用变化使径流量减少37.98mm,土地利用变化对流域径流量减少的贡献率为103.81%,而气候因素对流域径流量减少的贡献率仅占-3.81%。因此,可以认为第2阶段气候变化是径流减少的主导因素,第3阶段土地利用变化是径流减少的主导因素。

表4昕水河流域气候和土地利用变化对径流量的影响

3、讨论

昕水河流域1958—2015年径流量呈显著减少的趋势,对其归因分析,发现第2阶段(1975—2000年)气候变化对径流量减少的贡献率为73.14%,降水量的减少和潜在蒸发散的增加是径流变化的主要原因。张淑兰等[25]同样发现在泾河流域20世纪80,90年代,气候变化是径流减少的主要原因,其贡献率为80.9%。陈利群等[26]通过研究黄河源区气候和土地覆被变化对径流的影响结果表明,20世纪90年代以前,气候变化对径流减少的贡献率在95%以上。另外,Niedda等[27]对半干旱地中海生态系统意大利巴拉茨湖流域研究结果发现,气候变化因素对年均径流量减少的贡献率达到72%,这也与本研究结果相同。降水是昕水河流域径流的主要来源,当降水每减少1mm,流域径流将减少0.35mm。第2阶段(1975—2000年)流域降水量呈大幅度减少趋势(59.36mm),因此昕水河径流的减少受气候因素影响较为显著。

2001—2015年间,昕水河流域径流量变化主要受人类活动的影响,在黄土高原清水河流域[15]、泾河流域[28]、河西走廊的石羊河流域[29]所开展的径流对气候和土地利用变化的响应均得到同样的结论。影响流域径流变化的人类活动主要包括水利工程建设、土地利用变化以及水资源的开发利用等[30,31]。昕水河流域自1980年以来土地利用变化较大,特别是1999年实施退耕还林还草工程后,流域内林地、灌木林地、草地面积分别增加76.55,156.51,126.87km2,而耕地减少362.6km2(表2)。林地面积增加,会导致植被耗水量变大以及径流量的显著减少[25]亦或导致潜在蒸发散增加及径流的减少[15]。本研究中,2000年以后,尽管降水量略有增加,潜在蒸发散呈下降趋势,但是径流较之基准期下降幅度为68.49%,气候变化对径流减少的贡献率为-3.81%,即气候变化导致了径流增加。而人类活动的影响如植被恢复引起耗水量增加是导致径流减少的绝对主要因素,其贡献率为103.81%。同时,随着经济发展和人口密度增大,流域用水量逐渐增加,也导致径流量的减少。

王国庆等[32]采用有序聚类分析法对昕水河流域的研究结果显示,大宁站实测径流变化以1965年和1979年为分界研究径流变化及其对环境要素的响应。本文通过采用修订后的双累积曲线法确认的突变年份与王国庆等[32]关于昕水河流域的研究结果不一致,这主要是由两方面原因造成的;一方面是由于数据系列原因,本文的研究时间序列为1958—2015年,王国庆等[32]研究时间序列为1955—2008年;另一方面是检验方法的差异,王国庆等[32]研究是采用有序聚类分析方法推估水文序列的可能显著干扰点,该方法只针对于径流序列推求其最优分割点,未考虑降水等气候因素对径流的影响。本文采用修订的双累积曲线法来计算突变点,该方法不仅考虑了有效降雨和径流之间的关系,而且将降雨和径流数据依时间进程累加处理,对随机过程起到滤波作用从而消弱随进噪声[21]。郑培龙等[30]采用双累积曲线法确定了籍河流域径流的突变点,将研究时序划分为基准期和措施期;郭爱军等[28]采用双累积曲线法与滑动偏相关系数法的突变结果进行对比验证,提高了诊断结果的准确度。另外,考虑到降水量和潜在蒸发散是气候变化中影响径流的最重要因素,本文采用定性和定量相结合的方法将气候因子及人类活动(包括植被变化或者管理措施等)与径流变化紧密联系起来,使得分析径流对气候变化和人类活动的响应结果更为合理,但本研究中并未考虑流域自身的产流机制[31]。因此,在未来的研究中,还应结合水文模型对昕水河流域气候变化和土地利用对径流的影响进一步研究。

4、结论

(1)昕水河流域4期土地利用变化结果表明,该流域林地、灌木林地、草地、居民区及建筑用地面积呈增加趋势,耕地面积呈下降趋势,退耕还林还草工程的大规模开展,改变了流域下垫面特征。

(2)昕水河流域年径流量呈显著下降的趋势(Z=-5.84,p<0.0001),突变年份分别为1974年和2000年,而年降水(Z=-0.72,p=0.31)和潜在蒸发散(Z=-0.5,p=0.88)变化趋势均不显著。昕水河流域径流对气候因素的敏感性结果表明,降水每增加1mm,流域径流将增加0.35mm;潜在蒸发散每增加1mm,流域径流将减少0.10mm。

(3)人类活动和气候变化分别是不同阶段驱动流域径流变化的关键因素。1975—2000年,气候因素是流域径流减少的主导因素,其贡献率为73.14%,其中降水量、潜在蒸发散的贡献率分别为76.91%,-3.78%。2001—2015年,人类活动是径流减少的主导因素,其贡献率为103.81%。虽然降水量有所增加,但径流仍然显著下降,显示气候变化已不再是径流减少的主要原因,而土地利用变化已成为驱动径流变化的主导因素。

参考文献：

[2]刘昌明,张学成.黄河干流实际来水量不断减少的成因分析[J].地理学报,2004,59(3):323-330.

[5]王雅舒,李小雁,石芳忠,等.退耕还林还草工程加剧黄土高原退耕区蒸散发[J].科学通报,2019,64(增刊1):588-599.

[6]易浪,任志远,张翀,等.黄土高原植被覆盖变化与气候和人类活动的关系[J].资源科学,2014,36(1):166-174.

[7]赵安周,张安兵,刘海新,等.退耕还林(草)工程实施前后黄土高原植被覆盖时空变化分析[J].自然资源学报,2017,32(3):449-460.

[13]李志,刘文兆,郑粉莉,等.黄土塬区气候变化和人类活动对径流的影响[J].生态学报,2010,30(9):2379-2386.

[14]师忱,袁士保,史常青,等.滦河流域气候变化与人类活动对径流的影响[J].水土保持学报,2018,32(2):264-269.

[15]唐丽霞,张志强,王新杰,等.晋西黄土高原丘陵沟壑区清水河流域径流对土地利用与气候变化的响应[J].植物生态学报,2010,34(7):800-810.

[16]邓文平,李海光,余新晓,等.黄土高原吕二沟流域土地利用/覆被和气候变化对径流泥沙的影响[J].水土保持研究,2011,18(4):226-231.

[17]余新晓.森林生态水文研究进展与发展趋势[J].应用基础与工程科学学报,2013,21(3):391-402.

[19]芮孝芳.水文学原理[M].北京:高等教育出版社,2013.

[20]张盛霖,邓高燕,黄勇奇.Mann-Kendall检验法在Excel中的实现与应用[EB/OL].北京:中国科技论文在线,2014.http://www.paper.edu.cn/releasepaper/content/201406-448.

[21]穆兴民,张秀勤,高鹏,等.双累积曲线方法理论及在水文气象领域应用中应注意的问题[J].水文,2010,30(4):47-51.

[22]谢余初,巩杰,赵彩霞,等.干旱区绿洲土地利用变化的生态系统服务价值响应:以甘肃省金塔县为例[J].水土保持研究,2012,19(2):165-170.

[25]张淑兰,王彦辉,于澎涛,等.泾河流域近50年来的径流时空变化与驱动力分析[J].地理科学,2011,31(6):721-727.

[26]陈利群,刘昌明.黄河源区气候和土地覆被变化对径流的影响[J].中国环境科学,2007(4):559-565.

[28]郭爱军,畅建霞,王义民,等.近50年泾河流域降雨-径流关系变化及驱动因素定量分析[J].农业工程学报,2015,31(14):165-171.

[29]周俊菊,雷莉,石培基,等.石羊河流域河川径流对气候与土地利用变化的响应[J].生态学报,2015,35(11):3788-3796.

[30]郑培龙,李云霞,寇馨月,等.黄土高原藉河流域径流对气候和土地利用变化的响应[J].水土保持通报,2016,36(2):250-253.

[31]郭晓英,陈兴伟,陈莹,等.气候变化与人类活动对闽江流域径流变化的影响[J].中国水土保持科学,2016,14(2):88-93.

[32]王国庆,张建云,李雪梅,等.黄土高原昕水河流域径流变化及其对环境要素的响应[J].水土保持研究,2014,21(3):192-196.

冷曼曼,张志强,于洋,査同刚,亢小语.昕水河流域径流变化及其对气候和人类活动的响应[J].水土保持学报,2020,34(03):113-119+128.

基金:国家科技支撑重点项目(2015BAD07B030303).