黄土高原水土流失严重, 新中国成立以来, 先后开展了淤地坝建设、小流域综合治理、退耕还林(草)工程等一系列的水土流失治理措施^[1-2]。研究表明黄土高原一系列的水土流失治理措施显著减少了流域的侵蚀产沙量, 同时也在一定程度上使得黄河流域年径流量有所下降^[3-4]。黄河潼关水文站径流量由多年平均426亿m³锐减到近年来的231亿m^3[3]。黄河上游和中游年径流量分别从1919—1985年的251亿m³和168亿m³减少为1985—2010年的159亿m³和78亿m^3[4]。径流变化引起流域水文水资源系统的变化, 对流域水资源的开发利用产生了深远影响^[5-7]。因此, 研究流域径流演变规律及其驱动力因子有助于深刻认识水循环过程及变化特征, 同时也为流域水资源利用和生态环境建设提供科学依据。

无定河是黄河的一级支流, 流域位于毛乌素沙地南缘和黄土高原北部地区, 是黄河中游典型的风水两相复合侵蚀区^[8], 也是黄土高原水土保持措施实施的重点区域。全流域按地貌和水土流失类型可分为河源涧地区、风沙区和丘陵沟壑区(图 1)。其中西部河源涧地区(占流域面积12.2%)和东南部丘陵沟壑区(占流域面积34.3%)均为厚层黄土覆盖, 以水蚀过程为主。东南丘陵区侵蚀强烈, 年输沙量超过流域总输沙量的70%；西北部风沙区(占流域面积54.3%)以风力侵蚀为主, 侵蚀模数较小, 年输沙量不足流域总输沙量的6% ^[9]。目前, 国内学者就无定河径流变化、周期特征及影响因素等虽然已有相关研究^[10-13]。然而, 针对流域不同地貌类型区径流变化及形成原因的分区研究十分有限^[14-15]。由于下垫面条件的不同, 西北部风沙区和东南部丘陵沟壑区的流域具有明显不同的径流来源和产流机制^[14]。因此, 对不同地貌类型区径流变化趋势及其影响因素的探讨对于揭示全流域水沙演变规律及其形成原因具有重要意义^[15]。本文通过收集无定河及其典型支流海流兔河(风沙区)和大理河(丘陵区)出口控制水文站白家川、韩家峁和绥德站1960—2012年的降水量、蒸散量和径流量等气象水文资料, 利用Mann-Kendall(MK)非参数检验和启发式分割算法分析无定河干流及其不同地貌区支流径流量变化特征及其差异, 采用弹性系数法探讨影响流域径流量变化的主要因素及其贡献率, 并利用Hurst指数法预测流域干流及不同地貌区支流径流量变化趋势, 以期为无定河流域不同地貌区水资源管理和生态环境建设提供科学依据。

1 研究区域概况

无定河流域发源于陕西省定边县白于山北麓, 干流全长491.0 km, 流域面积30261 km²。流域出口控制站为白家川水文站, 控制流域面积为29662 km²。无定河流域属于温带大陆性干旱半干旱季风气候类型, 多年平均降水量为387.8 mm, 平均年径流深36.3 mm。流域水土流失面积23137 km², 平均侵蚀模数6090 t km^-2 a^-1, 是黄河粗泥沙的主要来源区之一^[15]。无定河流域自1950年开展水土流失治理工作以来, 治理范围日益广泛。20世纪50—60年代为起步阶段, 此阶段治理范围不大, 初步完成治理面积2153 km²；20世纪70—80年代为初具规模的治理阶段, 此阶段主要进行了沟道治理, 修建了大量淤地坝, 共建成淤地坝5929座, 占现有淤地坝总数的51.1%。80年代开始, 无定河成为国家重点治理区, 开始了大范围的治理工作, 截至1996年, 区内共退耕还林还草6734 km², 修建梯田966 km², 修成淤地坝11710座, 建成库容100万m³以上水库74座, 总库容14.9亿m³, 总治理面积8364 km², 占全流域水土流失面积的36.4%^[15]。

海流兔河和大理河分别为无定河位于西北部风沙区和东南部黄土丘陵区的一级支流(图 1)。其中海流兔河流域主要土壤类型为风沙土, 干流长约80 km, 流域面积2600 km², 出口控制站韩家峁水文站控制流域面积2452 km², 多年平均降水量367 mm, 年均潜在蒸散量1233 mm, 年径流量0.85亿m³；大理河流域主要土壤类型为黄绵土, 干流长170 km, 流域面积3906 km², 出口控制站绥德水文站控制流域面积3893 km², 多年平均降水量443 mm, 年均潜在蒸散量1199 mm, 年径流量1.4亿m³。截止2010年风沙区海流兔河流域林草面积约为1800km²；大理河流域林草面积为2300 km², 梯田面积为185km², 淤地坝约270多座^[15]。

2 研究方法及数据来源 2.1 研究方法 2.1.1 Mann-Kendall非参数检验

Mann-Kendall(MK)检验被广泛应用于降水、径流和水质等水文气象序列的趋势变化分析^[16-18], 基于MK检验的统计值Z值和P值来反应变化趋势。当Z大于0表示该序列呈现增加趋势, Z小于0表示序列呈现减少的趋势, 而P大于0.05表示该序列增加或者减少不显著, 而P小于0.05表示增加或者减少趋势显著。本文基于MK检验分析年降水、潜在蒸散发和年径流量的变化趋势。

2.1.2 启发式分割算法

某一时间序列X(t)由N个点组成, 从左到右分别计算每分个点左边和右边部分的平均值分别为μ₁(i)和μ₂(i)及标准差分别是s₁(i)和s₂(i), 则i点合并偏差S_D(i)可表示为：

(1)

式中, N₁、N₂分别表示i点左边与右边部分的点数。用t检验的统计值T(i)量化i点左右均值的差异：

(2)

式中, T值越大, 则说明该点左右两边两子序列的差异越明显。计算T中最大值T_max所对应的统计显著性P(T_max), 其计算公式如下：

(3)

由蒙特卡洛模拟可得：n=4.19lnN-11.54且δ=0.40。其中, N表示序列的长度, v=N-2, I_x(a, b)是不完全β函数。预先设定一个临界值P₀(P₀可取0.05—0.95), 当P(T_max)≥P₀, 则在该点处将此序列分割成左右两个均值差异较大的子序列, 否则不进行分割。对新的两个子序列不断进行迭代并重复以上操作, 直到子序列的长度小于l₀(l₀的取值则不应小于25)时便停止对其分割^[19]。

2.1.3 径流变化归因分析

假定流域多年蓄水量变化忽略不计, 基于Choudhury-Yang公式和流域水热耦合平衡方程^[20-22], 利用弹性系数法将流域径流变化分解为降水量影响, 潜在蒸散发影响以及下垫面变化等人类活动的影响。

2.1.4 重标度极差分析法

Hurst指数是一种时间序列的统计方法^[23], 可以用来定量表征序列的持续性。基本原理为：对一个时间序列X_t, t=1, 2, …, n。对于任意正整数τ≥1, 构造一个均值序列 < X>_τ, 其累积离差Y(t, τ)的具体算法为：

(4)

(5)

Hurst发现用标准差S除极差R建立一个无量纲比率, 满足关系式:

(6)

式中, α为常数, H为Hurst指数。对于不同的H, 意味着序列有不同的趋势变化：当H=0.5时, 表明序列是完全独立的, 即序列是一个随机过程；当0≤H < 0.5时, 意味着未来的变化状况与过去相反, 即反持续性, H值越接近0, 反持续性越强；反之, 1≥H>0.5时, 意味着未来的变化状况与过去一致, 即持续性, H越接近1, 接续性越强。Hurst指数能很好的揭示出时间序列中持续性或者反持续性的大小, 由此总结了Hurst指数分级表持续性和反持续性都分为5级(表 1)^[24]。

2.2 数据来源

1960—2012年白家川、绥德水文站、韩家峁的年径流量均来自黄河流域水文年鉴。无定河流域面降水量由16个雨量站(李家河、丁家沟、曹家岔、青阳岔、赵石窑、孟家湾、曹坪、韩家峁、补浪河、羊羔山、靖边、赵石畔、横山、涧峪岔、绥德、白家川)算术平均法求得；大理河流域面降水量由绥德、曹坪、李家河、青阳岔4个雨量站算术平均法求得；海流兔河流域由韩家峁、补浪河2个雨量站算术平均法求得。

气象数据为横山站、榆林站、绥德站1960—2012年的逐日气象资料包括最高气温, 最低气温、平均气温、相对湿度、风速、日照时数, 来自于中国气象科学数据共享服务网。潜在蒸散发量通过Penman-Monteith公式计算得到, 无定河流域潜在蒸散发量为三站的算术平均求得, 大理河流域为绥德站, 海流兔河流域为横山站。由于各流域面积不同, 径流量差异较大, 为便于比较, 本文采用径流深(径流量/流域面积)反映流域径流量变化趋势。

3 结果分析 3.1 降水量、蒸散量和径流量变化特征分析

无定河及其支流大理河和海流兔河流域1960—2012平均年降水量分别为387.7, 442.6 mm和367.4 mm, 位于黄土丘陵区的大理河流域较风沙区的海流兔河流域年均降水量高75.2 mm。三条流域年均降水量均在1960s年代和2000s年代高于其他时段(图 2), 但其降水量变化的MK检验无显著趋势(P>0.05)(表 2)。

无定河、大理河和海流兔河三条流域1960—2012平均年蒸散量分别为1195.5, 1199.2 mm和1233.3 mm, 大理河流域年均蒸散量较海流兔河低34.1 mm。无定河和大理河流域年蒸散量呈增加趋势, 这种增加在1990s年代和2000s年代增强(图 2)。大理河流域年蒸散量MK检验结果在P < 0.01水平显著, 无定河流域年蒸散量MK检验结果在P < 0.1水平显著(表 2)。相比之下, 海流兔河流域年平均蒸散量在1970s—1990s年代呈增加趋势, 但2000s年代流域年平均蒸散量有所下降。总体上, 海流兔河1960—2012年平均蒸散量的MK检验不显著(表 2)。

无定河、大理河和海流兔河三条流域1960—2012平均年径流深分别为36.3, 35.5 mm和3.4 mm, 大理河流域年均径流深较海流兔河高31.1 mm。3条流域年径流深变化的MK检验均呈极显著趋势(P < 0.01)(表 2)。其中, 无定河流域从1960s年代到2000s年代年径流深呈连续下降趋势；大理河年径流深在1970s和1980s持续下降, 但在1990s年代有所回升, 随后在2000s年代再次下降；海流兔河年径流深在1970s至1990s年代持续下降, 但在2000s年代后有所回升。此外, 位于黄土丘陵区的大理河流域径流深年际变化较风沙区的海流兔河大。

3.2 径流突变点及变化归因分析

无定河流域径流变化突变点在1979年和1996年, 海流兔河流域径流变化突变点在1971年和1990年, 而大理河流域径流变化突变点在1971年(图 3)。根据突变点将无定河流域径流序列分为3个时期, 1960—1979年为基准期, 1980—1996年和1997—2012年为研究期。为方便对比, 将大理河和海流兔河流域径流序列也分为3个时期, 1960—1970年为基准期, 1971—1990年和1991—2012年为研究期。与基准期相比, 无定河流域1980—1996年研究期径流深减少了30.5%；1997—2012年径流深减少了46.4%。相比基准期, 在1972—1990年大理河流域径流深的减少量为29.2%, 大于海流兔河流域的20.5%；而在1991—2012年, 大理河流域的径流深变化量为31.6%, 小于海流兔河流域的38.6%(表 3)。

采用弹性系数法对三条流域两个研究期的径流变化进行归因分析(表 3)。无定河流域1980—1996年的径流深减少了14.5 mm, 人类活动影响占主要作用, 贡献率为66.8%；其次为降水量, 占到减少总量的30.3%；蒸散量变化仅占减少总量的2.8%；1997—2012年径流深减少了22.2 mm, 其中人类活动导致径流量减少98.2%, 降水量增加使得径流总量增加1.7 mm, 蒸散发导致的径流减少量占9.5%。大理河流域径流深在1972—1990年减少了13.6 mm, 在1991—2012年减少了14.7 mm, 两个时期减少量相差不大；而且两个时期人类活动影响造成的径流量的减少占到50%左右, 也相差不大；但是气候变化对径流量的影响中, 1972—1990年, 降水的减少对径流下降的贡献占到47.9%, 蒸发变化贡献不足1%；在1991—2012年, 降水和蒸发变化分别贡献了21.8%和29.0%。海流兔河径流深在1972—1990年减少了0.9 mm, 在1991—2012年减少了1.7 mm, 是1972—1990年减少的一倍。分析两个时期径流减少归因发现, 1972—1990年, 降水量占到了33.3%, 蒸散量占到22.2%, 而人类活动占到了44.4%；而1991—2012年, 人类活动的影响是造成径流减少的主要原因, 占到82.4%, 气候变化中蒸散量贡献为17.6%, 由于该时期平均降水量与基准期一样, 降水对径流减少的贡献为0。

3.3 基于Hurst指数的径流量未来趋势分析

无定河、大理河和海流兔河流域年降水量、蒸散量、径流量序列的长期相关性特征为持续性, H值大于0.5, 未来的变化与过去的变化趋势一致(表 4)。3个流域年降水变化趋势的持续性强度都较弱；无定河和海流兔河流域蒸散量增加趋势的持续性强度也较弱, 但是大理河流域的蒸散量增加趋势的持续性较强；无定河流域径流量减少趋势持续性最强, 海流兔河流域径流量减少趋势的持续性较强, 而大理河流域径流量减少趋势持续性较弱(表 4)。

4 讨论

黄土丘陵区的河流是无定河流域侵蚀产沙的主要源区, 也是水土保持措施实施的重点区域。人类活动对黄土丘陵区的大理河流域1972—1990年和1991—2012年两个时期径流减少的贡献率分别为48.6%和51.7%(表 3), 研究表明梯田、林地、草地和坝地建设等水土保持治理措施是流域侵蚀产沙减少的主要原因, 也是导致流域径流减少的重要原因。大理河流域的水土流失治理始于1960s年代, 流域水土保持措施面积在1970s和1980s年代快速增加, 尤其是坝地面积在1970s年代较1960s年代增加了2倍多(表 5)^[25]。对大理河流域3条典型支流不同年代淤地坝数量统计发现, 流域85%以上淤地坝修建于1970s年代(表 6)。截止到2002年, 大理河流域共有淤地坝3100余座, 控制70%的流域面积。与1960s年代相比, 仅坝库拦蓄作用导致大理河年径流量下降800万m³, 约占流域径流量的10%^[26]。高海东等^[27]阐述了淤地坝对径流过程的影响机理, 指出随着坝地的淤积, 沟道形状由原来的V型沟道逐渐演变为U型沟道, 其比降降低, 过流断面面积变大；淤地坝作为障碍性节点, 显著削弱了流域的洪峰流量。因此, 1970s年代大规模水土保持治理措施, 尤其是淤地坝建设显著减少了大理河流域侵蚀产沙, 也在一定程度上减少了流域径流量, 成为大理河径流量在1971年出现拐点的重要原因。另外, 由于1970s—1980s年代建设的淤地坝拦蓄寿命大部分小于20年, 随着淤地坝的淤满, 其拦蓄作用下降, 导致大理河年径流量在1990s年代有所回升。在相同干旱指数情况下, 林草覆盖率越大, 流域产水越少^[28]。植被恢复对黄河径流减少的贡献达26%^[29]。2000s年后, 国家退耕还林还草政策实施, 大理河流域水土保持建设达到新的高潮, 使得流域年径流量在2000s年代再次出现下降。此外, 大理河流域20世纪90年代后期蒸散量开始显著增加(图 2), 也加剧了流域径流量的下降。

与大理河流域相比, 位于风沙区的海流兔河流域水土流失微弱, 水土保持措施很少, 流域内并无淤地坝分布(图 1)。然而, 海流兔河流域第一个径流突变点也发生在1971年, 人类活动对海流兔河1972—1990年和1991—2012年两个时期径流减少的贡献率分别为44.4%和82.4%。据统计, 海流兔河流域目前共有10座较大型水利工程, 主要用于农田灌溉。其中, 在1970s年代和1990s年代先后修建4座, 而在1980s年代和2000s年代各有1座(表 7)^[30]。坝库的修建年代与海流兔河径流变化的两个拐点1971和1990基本一致, 表明流域内过度的农田灌溉用水的增加可能是引起径流减少的主要原因。此外, 海流兔河径流量在2000s年代有所回升, 其原因一方面与流域2000s年代降水量增加, 蒸散量下降有关(图 2)；另一方面也与流域退耕还林还草后农地面积下降, 灌溉需水减少有关。

无定河干流及其支流大理河和海流兔河径流未来仍然呈减少趋势, 流域水资源利用形势日趋严峻。归因分析表明, 人类活动是引起无定河流域20年来径流减少的主要原因。未来迫切需要从流域社会、经济和生态环境可持续发展的角度制定水资源利用规划, 开展节水型流域建设, 加强雨水、工业用水、农业用水和生活用水资源的循环高效利用。其中, 在大理河流域要优化现有植被建设布局, 利用乡土树、草种逐步替代耗水高、生长差的植被类型, 减少流域蒸散发, 减缓径流下降。在海流兔河流域要适当控制农田灌溉面积, 提高农田灌溉用水效率, 在必要情况下增加退耕还林(草)面积, 减少灌溉用水的同时提高流域水源涵养能力。

5 主要结论

本研究发现1960—2012年无定河及其不同地貌区典型支流海流兔河和大理河年径流量均发生显著下降趋势, 但其在不同年代变化略有差异。其中, 位于黄土丘陵区的大理河年径流量在1960s—1980s持续下降, 但在1990s年代有所回升, 2000s年代再次下降；而风沙区的海流兔河年径流量下降一直持续到1990s年代, 在2000s年代有所回升。尽管大理河和海流兔河径流突变点均出现在1971年, 但不同地貌区流域径流下降的原因有所不同。在黄土丘陵区, 大规模的水土保持治理措施, 尤其是淤地坝建设等人类活动是大理河径流量在1971年出现拐点的重要原因；而在风沙区, 过度的农田灌溉用水增加是引起海流兔河径流变化的主要原因。Hurst指数表明, 未来无定河流域及其不同地貌区支流海流兔河和大理河年径流量均持续下降。