基于Budyko互补关系法的嘉陵江径流变化归因研究

胡淼, 王润瑄, 姬广兴

(1.教育部人文社科重点研究基地,辽宁师范大学 海洋经济与可持续发展研究中心,辽宁 大连 116029;2.辽宁省“海洋经济高质量发展”高校协同创新中心,辽宁 大连 116029;3.河南农业大学 资源与环境学院,河南 郑州 450046)

目前,径流的非平稳性在全球范围内普遍存在[1]。在全球变化背景下,蒸发、降水以及径流的时空分布不断发生变化[2],同时人类通过修建拦水坝和引水渠及城镇化建设等活动改变了流域的下垫面结构,从而对河川的径流过程产生了深刻影响[3]。因此,在全球变化背景下对影响流域径流变化的因素进行定量计算,对于人类认识自然、保护自然以及采取有效措施实现区域可持续发展目标具有重大意义。

当前对流域径流变化的原因分析已成为全球的研究热点,主要分析方法有水文模型法[4]、流域对比试验法[5]、多元统计分析法[6-7]和经验模型法[8-9]等。水文模型法的模型主要有SWAT模型、HBV模型、VIC模型等,这些模型具有不错的物理机制并且适用于不同范围、不同复杂度的流域水文模拟,但是这些方法均需要较多的数据且对数据的精度要求较高,并且模型参数与结构存在较大的不确定性[10]。流域对比试验法的物理意义明确,但不适用于大流域范围。多元统计分析法简单,但需要长期大量的观测数据作为样本以支撑统计学分析,并且该方法缺乏对流域径流量变化物理意义的解读。对于经验模型法,当前具有代表性的是基于Budyko假设的水热耦合平衡方程法,该方法内含参数较少、计算过程简单且具有明确的物理意义,目前已成为长时间序列下研究流域径流变化原因的常用方法[11-14]。该方法可以定量评估气候变化和人类活动对径流变化的敏感性和贡献率,目前应用范围广泛。但该方法在计算偏导数部分时所利用的方法不唯一(主流的方法有全微分法和外推法),导致计算结果与实际观测结果出现一些难以解释的残差。针对此问题,文中借鉴ZHOU S等[15]的相关成果,在Budyko假设的基础上建立了一种互补关系模型,在该模型中将影响径流变化的因素划分为降水、潜在蒸散发和下垫面3类,以便在理论计算上消除误差,提高径流变化归因分析的准确度。

嘉陵江是长江支流中流量仅次于岷江且流域面积最大的河流[16]。由于人类活动和全球气候变化的影响,嘉陵江流域人地关系矛盾日益突出。本文基于Budyko假设的互补关系模型量化气候因素和下垫面因素对径流变化的影响,并对下垫面因素中的植被与土地利用状况进行分析讨论,该研究结果可为嘉陵江流域水资源可持续利用与规划管理提供借鉴。

1 流域与数据

1.1 嘉陵江流域简介

嘉陵江干流全长1 345 km,流经陕、甘、川、渝4省(市),流域集水面积约16万km2,平均年径流量682亿m3。嘉陵江流域上游地势以高山为主,平均海拔超过3 000 m,中下游地势趋于平缓,海拔多在200～400 m,水流亦趋于平缓。嘉陵江流域气候属于亚热带季风气候类型,该地区降水丰沛,雨热同期,季节分配极不均匀,山地上植被较浓密,丘陵地区上的植被覆盖率相对较低。文中的水文站选用的是位于嘉陵江干流下游的北碚水文站(经度106°396′E、纬度29°806′N),如图1所示。

图1 嘉陵江流域及站点分布示意图

1.2 数据来源

1961—2020年嘉陵江北碚水文站的径流数据来源于长江水利网(http://www.cjw.gov.cn)和长江水文水资源公报。基础气象数据采用流域范围和周边的37个气象站点在研究时段内所测得的气象数据,包括2 m高处风速、最低和最高温度、净辐射量、实际水气压和日降水量等,这些数据来源于中国气象局网站(http://www.cma.gov.cn)。嘉陵江流域的DEM数据与1980年和2020年土地利用数据分别来自于中国科学院科学数据中心(https://www.casdc.cn)与中国科学院资源环境科学数据中心(https://www.resdc.cn),数据分辨率都为1 km。嘉陵江流域1982—2019年的归一化植被指数(normalized difference vegetation index,NDVI)数据集来源于网站https://www.ncdc.noaa.gov/cdr/terrestrial/normalized-difference-vegetation-index。

2 研究方法

2.1 趋势和突变分析

潜在蒸散发数据通过Penman-Monteith(彭曼公式)计算得到,计算方法参考文献[17]:通过Kriging(克里金)插值得到流域年降水量和年潜在蒸散发数据;数据的趋势分析通过线性倾向估计法和滑动平均法进行,滑动平均法可以削弱部分极值所造成的波动影响,从而能更清晰看出数据变化趋势[18]。该方法的原理简单,在此不再赘述。

年径流序列的突变时间点通过Pettitt突变检验法确定。Pettitt突变检验法是一种基于Mann-Whitney非参数秩检验法[19],其是对原始序列数据所获得的秩进行统计分析来确定突变点。其中Mann-Whitney非参数统计量Ut,N的计算式为:

(1)

用Ut,N来检测两个样本x1、x2、…、xt和xt+1、xt+2、…、xN是否来自同一序列。根据式(2)和式(3)分别算得统计量Kt,N和检验两个样本累计分布函数的最大差值所求的累计概率p,即:

Kt,N=max|Ut,N|, 1≤t≤N;

(2)

(3)

2.2 气候和下垫面对径流影响的估算

2.2.1 下垫面参数n的计算

采用CHOUDHURY B[20]和YANG H B等[21]在Budyko假设基础上所提出的单参数Budyko方程以及流域水量平衡方程来计算下垫面参数n的值,具体推算过程参见文献[13],其表达式为:

(4)

式中:R为径流深,mm;P为降水量,mm;ET0为流域内指定时段潜在蒸散发量,mm;n为表征流域下垫面状况的参数。n值的大小与地形因素、土地表面渗水能力以及植被状况等相关[22]。在式(4)中,若已知R、P和ET0,则可求出参数n。

2.2.2 径流对气候和下垫面的弹性系数

(5)

(6)

(7)

2.2.3 气候和下垫面对径流变化的贡献率

ZHOU S等[15]通过假设P和ET0相互独立,提出径流对降水和潜在蒸散发的弹性系数具有互补关系的方程,文中利用该方程计算贡献率,该互补关系的方程为:

(8)

由于气候条件和流域条件对径流改变的贡献不能分开独立分析,原因是径流发生改变的路径不能被确切知道。因此,引入加权因子α来表征这种变化,具体如下:

(9)

式中α的取值范围为[0,1]。基于ZHOU S等[15]的实际测算,α取0.5时的模拟效果最好。

降水量的变化量(ΔRP)、潜在蒸散发量的变化量(ΔRET0)、下垫面特征值的变化量(ΔRn)的计算公式分别如下:

(10)

(11)

(12)

式中下角标1、2分别代表基准期和变化期。

贡献率Cx的计算式为:

(13)

2.2.4 土地利用转移矩阵

根据1980年与2020年嘉陵江流域土地利用数据,借助ArcGIS 10.5软件中的融合与相交工具,计算得出嘉陵江流域土地利用类型在这两个时间段内的变化情况。土地利用转移矩阵可以用于分析某一区域土地利用类型的转入与转出情况以及统计其面积变化[24-25]。

3 结果与讨论

3.1 水文气象要素趋势分析与突变点检验

嘉陵江流域1961—2020年径流深变化趋势如图2(a)所示。由图2(a)可看出:嘉陵江流域年径流深R呈逐年递减趋势,递减速率为11.378 mm/(10年);从径流深5年滑动平均可以看出,年径流深表现出规律性的丰-枯水年交替特征;年径流深最大值为1983年的683.02 mm,最小值为1997年的199.79 mm;对年径流深进行非参数M-K趋势分析,统计量Z值为-1.32,绝对值小于1.96,未达到95%的置信水平,表明嘉陵江流域年径流深的下降趋势不显著。

图2 年径流量变化趋势和Pettitt突变点检验

对嘉陵江流域1961—2020年径流深进行Pettitt突变检验,结果如图2(b)所示。图2(b)中显示:嘉陵江径流深突变发生的时间为1985年,并且达到了α=0.01的显著水平。

嘉陵江流域水文气象要素特征见表1。由表1可知:1961—2020年嘉陵江流域的平均径流深为427.52 mm,极值比为2.42;流域降水以及蒸发等气象因素直接影响河川径流,对嘉陵江流域年降水量和年潜在蒸散发量时间序列进行M-K趋势检验,其Z值分别为-0.26和0.62,均未呈现显著的变化趋势。

表1 嘉陵江流域水文气候特征

嘉陵江流域的年降水量和年潜在蒸散发量的趋势分析如图3所示。由图3可知:嘉陵江流域的年降水量变化速率为-0.348 mm/(10年),最大值为1 132.22 mm,出现在1981年,最小值为625.14 mm,出现在1997年;潜在蒸散发量的变化速率为2.869 mm/(10年),最大值为960.22 mm,出现在2013年,最小值为783.43 mm,出现在1989年。

图3 年降水量与年潜在蒸散发量变化趋势

3.2 径流的敏感性分析

根据Pettitt突变检验的结果,将1961—1985年确定为基准期,将1986—2020年确定为变化期,对影响流域径流的三大因素进行弹性系数计算,结果见表2。由表2可以看出:变化期相对于基准期来说,嘉陵江流域其潜在蒸散发量、径流深以及降水量均减少,分别减少0.13、68.01、33.22 mm,下垫面特征值增加了0.16;从弹性系数值可以看出,嘉陵江流域径流变化量与降水量(P)呈正相关、与潜在蒸散发量(ET0)和下垫面参数(n)值呈负相关,说明n值增大将致使流域径流深减小;变化期与基准期相比,降水量的弹性系数εp从1.49增大到1.60,表明当降水量同样增加10.0%时,1985年以前会使嘉陵江流域径流深增加14.9%,1985年之后则会使其增加16.0%,而潜在蒸散发量的弹性系数εET0和下垫面的弹性系数εn均减小,分别从-0.50减小到-0.63、从-0.66减小到-0.76,说明当变化期的潜在蒸散发量和下垫面特征参数值增加10.0%,将会导致嘉陵江流域径流深分别减少6.3%和7.6%。由于径流对各个影响因子敏感程度的强弱可由弹性系数绝对值的大小反映出来,因此,从表2中可看出,近60年来,嘉陵江流域径流深对降水量、潜在蒸散发量和下垫面的敏感性均增强。对比三者弹性系数的绝对值大小,其中εp的绝对值最大,其次为εn的,εET0的最小,说明流域径流对降水补给最为敏感,降水量的多少直接影响着径流深的大小,径流对潜在蒸散发的敏感性最弱。

表2 嘉陵江流域各影响因子弹性系数计算结果

图4(a)—(c)中反映了嘉陵江流域径流深对潜在蒸散发量、下垫面特征值和降水量的弹性系数随时间的变化特征。图4(a)—(c)中εET0和εn的值都小于0,εp的值都大于0,εET0值的范围为-0.84～-0.32,εn值的范围为-1.12～-0.47,εp值的范围为-1.84～-1.32;潜在蒸散发的速率为-0.028/(10年),下垫面与降水弹性系数的变化速率分别为-0.011/(10年)和0.028/(10年),表明嘉陵江流域径流深对潜在蒸散发量和降水量的敏感性增速最快。三者的弹性系数的绝对值大小随时间都呈现出增大趋势,下垫面弹性系数的变化趋势不显著,径流深对降水量和潜在蒸散发量的弹性系数呈现显著增加(Z>2.58)的趋势。在研究时段内径流深对潜在蒸散发量、下垫面和降水量的敏感性增强的原因为:嘉陵江流域的生态稳定性在减弱,随着人类活动强度的增大,对流域资源的过度利用使得区域生态系统平衡失调[26]。与之前相比,同样微小的变化会导致流域更大的风险或波动。

图4 各因素弹性系数与干旱径流指数年际变化

干旱指数与径流系数变化情况如图4(d)所示。干旱指数是ET0与P的比值,可以反映气候干旱的程度,径流系数是R与P的比值,反映了流域降水转化为径流的比率。由图4(d)可看出:干旱指数在1961—2020年呈现上升趋势,说明嘉陵江流域气候类型向暖干化方向变化,与全球变暖趋势一致,干旱指数的上升速率为0.003/(10年);径流系数呈现下降趋势,下降的速率为-0.012/(10年),径流系数值减小说明嘉陵江流域降水转变为径流的量减少,综合反映出了流域的降水量更多地消耗于下垫面、蒸发等自然地理要素和人类活动当中。

3.3 径流变化归因分析

通过式(10)—(12)对嘉陵江流域径流深、降水量、潜在蒸散发量和下垫面特征值进行计算,所得结果见表3。由表3可知:降水和下垫面所引起的北碚站径流深分别减少24.57 mm与43.48 mm,潜在蒸散发引起的北碚水文站径流深增加了0.036 mm。与基准期相比,在变化期,嘉陵江流域径流深减少了14.53%,模拟的径流变化量(ΔR)为-68.01 mm,与实际观测的径流深变化量(ΔR′)相同,说明此种计算方法所得到的结果相比于传统的微分法计算结果更为准确。对于嘉陵江流域径流深减少68.01 mm来说,各影响因子的贡献率分别为:降水量的贡献率为36.12%,下垫面变化的贡献率为63.93%,潜在蒸散发的贡献率为-0.05%,说明潜在蒸散发对径流的减少起到了负贡献的作用。

表3 嘉陵江流域各影响因素对径流量变化的贡献率

综上表明,引起嘉陵江流域径流深减少的主要因素为流域下垫面状况的改变,下垫面通过对流域降水的拦截、填洼下渗等影响流域的径流汇水过程,这可能与嘉陵江流域近60年来植被覆盖率增加、土地利用类型改变和人类活动对水资源的直接提取利用有关。次要因素为流域降水量的减少,降水的多少对径流的影响是直接的。对流域径流变化影响最小的是潜在蒸散发。

3.4 讨论与建议

3.4.1 归一化植被指数与径流的关系

嘉陵江流域径流减少的主要原因为流域下垫面状况的改变,下垫面参数n与归一化植被指数(NDVI)的年际变化以及两者关系如图5所示。

图5 下垫面参数n与NDVI的年际变化和拟合结果及嘉陵江流域1980—2020年土地利用类型变化情况

从图5(a)中可以看出:嘉陵江流域1961—2020年下垫面参数n值的变化表现为递增趋势,递增速率为0.039/(10年),通过M-K趋势检验,其Z值为3.44,大于2.58,说明其递增趋势显著。文献[23]表明,流域下垫面参数n值的大小主要与植被覆盖程度、土壤类型、地形坡度以及地表渗水透水能力等因素有关。一般认为在几十年的时间变化序列中,较大流域研究范围内其地形坡度以及土壤的类型等因素不会发生显著的变化,流域下垫面参数n值的变化主要与植被覆盖程度有关[27],且植被覆盖度与NDVI具有极强的一致性。

图5(b)反映出嘉陵江流域1982—2019年年均NDVI值的变化趋势,其递增速率为0.031/(10年),Z值为8.01,大于2.58,也呈现显著的递增趋势。NDVI值递增表明嘉陵江流域植被覆盖程度整体呈增加趋势,这可能与“退耕还林、荒山造林”与“长治”工程等政策在嘉陵江流域取得较好的成效有关。21世纪,嘉陵江流域植被覆盖度为0.70～0.81,覆盖程度良好[28-29]。植被覆盖度的增加可以增大植物对降水的截留量,延缓地表径流的时间,起到对水分的涵养作用,增大水分的蒸散发量,减小降水对径流的作用,导致径流减少。下垫面参数n与年均NDVI的拟合结果如图5(c)所示。由图5(c)可看出,回归系数a与截距b分别为1.090 8和0.683 7,决定系数R2为0.074 1,说明下垫面参数n与年均NDVI并不具备较好的正相关关系,因此不能进行相关的定量计算。宁怡楠等[30]根据WANG S等[31]拟合得到的黄河流域河龙区间植被覆盖度与径流系数的函数关系,定量分析了植被覆盖度对径流减少的贡献率。但该函数关系的建立仅适用于人类活动影响微弱的区域,对于流域范围大、人类活动影响强的区域具有局限性。图5(c)可以说明,嘉陵江流域下垫面参数n的值并不主要由植被覆盖度决定,也可能由人类活动所导致的土地利用类型改变所决定。

3.4.2 土地利用类型与径流的关系

嘉陵江流域1980—2020年土地利用类型变化状况如图5(d)和图5(e)所示。

从图5(d)中可以看出:嘉陵江流域主要发生的是耕地、林地与草地三者之间的土地利用类型变化,耕地与林地的面积少量增加,草地面积略微减少。耕地向林地的转出面积最多,其次是林地向耕地的转出面积,然后是草地向耕地的转出面积;1980年的耕、林、草地面积占比分别为44.52%、31.13%和22.58%,2020年的耕、林、草地面积占比分别为43.24%、31.78%和21.88%。从图5(e)中可以看出:嘉陵江流域主要土地利用类型为林地、耕地和草地。近40年来,嘉陵江河道两岸大片耕地转变为建设用地,产生该现象的原因主要是嘉陵江流域自20世纪80年代以来经济的快速发展推动并加快了城镇化建设,并且我国在近几十年相继在嘉陵江流域进行了水能开发;土地利用结构和开发强度的变化,都对流域的产汇流机制造成了影响;农业、工业以及人类生活用水量的增加也使得河川径流减少。总的来说,由于经济发展和人口增长所导致的嘉陵江流域水资源的开发利用是径流减少的重要原因。

4 结论与建议

文中通过运用Budyko水热耦合平衡方程的互补关系模型从气候因素和下垫面因素方面定量揭示1961—2020年嘉陵江流域径流减少的原因,主要结论如下:

1)1961—2020年,嘉陵江流域径流深和降水量都呈现减少的趋势,潜在蒸散发量呈增大趋势。

2)1961—2020年,嘉陵江流域径流深发生突变的时间点为1985年,嘉陵江流域径流对降水敏感性最强,其次为下垫面状况,敏感性最弱的是潜在蒸散发量。

3)影响嘉陵江流域径流减少的三大因素中,下垫面因素的贡献率为63.93%、降水因素的贡献率为36.12%,潜在蒸散发因素的贡献率为-0.05%。说明嘉陵江流域径流减少主要受下垫面作用影响。嘉陵江流域径流减少的原因有植被覆盖度增加、土地利用类型空间结构上的转变和人类工农业活动上的取用水。

嘉陵江流域内经济发展水平相对滞后,在当前区域经济发展需求强烈的背景下,需要多保障敏感河段内的生态流量,对嘉陵江流域河湖湿地多加保护,防止生态稳定性持续下降。另外,要统筹优化嘉陵江流域水资源的开发保护,统筹全流域水资源配置,合理分配生产、生活和生态用水,优化水资源开采和利用形式,提高水资源的利用效率,防止粗放利用。文中研究的不足之处在于对流域径流变化的归因分析中未考虑人工修建的水利工程设施对水域流量的拦截作用,希望未来能够找出水利设施以及土地利用类型与径流变化的函数关系,从而更加精确地定量分析径流变化的原因。