2005—2030 年无定河流域土地利用变化对径流的影响

2023-03-10赵雪岩

人民黄河 2023年3期

赵雪岩，张鑫，孙媛

（西北农林科技大学水利与建筑工程学院，陕西杨凌 712100）

随着社会发展，土地利用变化通过调节产汇流和降水蒸发过程来影响水文循环和水资源量的变化［1－3］。自1960 年以来，无定河流域受人类活动（如退耕还林还草、防风治沙、开荒种田等）的影响，其水文循环发生了显著变化［4－7］。前人大多聚焦于无定河流域植被恢复前后水文通量的大幅度变化进行研究，如王计平等［4］认为1995 年前后无定河流域土地流转频繁；任宗萍等［5］认为水土治理措施使无定河支流大理河径流量在1996 年出现拐点；张守红等［6］认为水保措施是无定河流域径流量减少的主要原因。然而，水保措施实施后径流、生态的大幅变化趋势是不可持续的，国家山水林田湖草沙等自然资源共同发展演化的趋势是从大修大补的退耕还林（草）转向小微型、质量型发展［8］。

SWAT（Soil and Water Assessment Tool）模型是一个基于物理机制的、复杂的半分布式水文模型［9－16］，该模型能够模拟流域水文循环过程和土地利用变化过程等地表自然物理化学过程，并在国内外得到了广泛应用［17－20］。Yang 等［9］认为滦河流域退耕还林使夏季流量增大，其他季节流量减小；傅春等［10］研究表明抚河流域林地增加、耕地减少使径流量减少；杨倩等［11］研究表明汉江上游退耕还林使年径流量减少；Ma 等［12］研究表明岷江流域耕地增加、草地减少使产水量减少；李烁阳等［13－16］在渭河流域和汾河上游均开展不同角度的研究。

本文以无定河流域为研究对象，运用SWAT 模型，利用五期土地利用数据（2005 年、2015 年、2020年、2025 年、2030 年，其中2025 年和2030 年为CAMarkov 模型预测结果）驱动模型，分析4 种土地利用变化及其对径流的影响。

1 研究区概况与数据来源

1.1 研究区概况

无定河流域位于北纬37°2′—39°1′、东经107°47′—110°34′（见图1），多年平均径流量为11.51 亿m3。流域处在东部季风气候和西北干旱气候的过渡带，多年平均气温为7.8～9.6 ℃，多年平均降水量为300～500 mm，降水量由东南向西北递减，年际变化大，降水多集中在夏秋两季，年内分配不均［4－7］。

图1 无定河流域高程及水文站和气象站分布

1.2 数据来源与处理

1.2.1 空间数据处理

空间数据介绍详细见表1。其中，2005 年、2015年、2020 年三期土地利用数据用于搭建SWAT 模型的土地利用数据库，原始数据中的二级分类需要重分类为草地、林地、耕地、未利用地、居民用地和水域6 种类型（见表2，其中土地利用分类属性前的数字为代号）；2010 年土地利用数据用于构建CA－Markov 模型中的土地转移概率矩阵。HWSD 用于搭建土壤属性数据库，同时还需输入默认参数值，如上层盐度默认0.1等；依托HWSD，利用SPAW 软件进行计算后获取其他参数。另外，石灰性雏形土和石灰性冲积土根据属性的不同进一步细分为1 类和2 类。土壤类型和面积统计详见表3 和图2（空间参考坐标系统一设置为WGS 1984－utm－48n）。

表1 空间数据介绍

表2 土地利用数据重分类

表3 不同类型土壤的面积统计

图2 无定河流域土壤分类

1.2.2 气象水文数据处理

气象数据包括逐日降水量、最高／最低气温、太阳辐射、风速、相对湿度等，来源于中国气象数据网，气象站分布见图1。水文数据为逐日平均流量，来源于黄河流域水文年鉴（中游区上段），水文站分布见图1。两者时间段均为2006—2018 年。另外，各类气象数据需要按照站编号、站名、经纬度、高程做成索引表（txt文件），缺测数据使用天气发生器模拟进行插补。

2 研究方法

2.1 土地利用数据模拟

CA－Markov 模型结合了CA 模型模拟复杂系统空间变化的能力和Markov 模型的定量预测优势，可作为土地利用模拟工具，并在IDRISI 软件中实现，其实质是建立土地转移概率矩阵并利用卡帕系数Kappa评估概率矩阵模拟出的土地利用数据和实测数据的差异，从而建立模型用于输出预测数据［21－23］。

式中：PO为恰当被模拟的部分；PC为由偶然因素导致的被希望纠正的部分；PP为观测和模拟完全匹配的理想比例，取为1；Lij为i类土地利用类型转变为j类土地利用类型的概率；n为土地利用类型总数。

如果Kappa ＝1，说明二者是一致的；如果0.75 ≤Kappa ＜1，说明高度一致；如果0.5 ≤Kappa ＜0.75，说明中度一致；如果Kappa＜0.5，说明不一致［21］。

本文输出2025 年和2030 年的无定河流域土地利用数据的主要步骤［21－22］：①利用Markov 模型建立时间间隔为5 a（2005—2010 年）的过渡概率矩阵；②设定预测的起始年份为2010 年，并利用CA 模型设置循环迭代次数为5 次，分别用来模拟2015 年土地利用数据。③模拟2015 年土地利用数据与实际2015 年土地利用数据进行比较，以Kappa来评估模拟效果。

2.2 SWAT 模型

在SWAT 模型中，一个流域被划分为由河流网络相互连接的许多子流域，每个子流域根据土地利用类型、土壤类型和坡度进一步划分为水文响应单元（HRUs），模拟水的流动轨迹从HRUs 到子流域再到流域出口断面。水文过程基于如下水量平衡方程［18－20］：

式中：SWt为最终土壤含水量，mm；SW为初始土壤含水量，mm；R、Qsurf、ET、P、QR分别为渗漏量、径流量、蒸散发量、日降雨量、地下水含量，mm。

常用Nash－Suttcliffe 系数（Ens）、相关系数（r）［18］和百分比偏差（PBIAS）［24］来评估SWAT 模型的模拟效果。

式中：Qoi为i时的观测径流量；Qsi为i时的模拟径流量；为观测径流量平均值；为模拟径流量平均值；n为观测数据总数［24］。

当r ＝1 时，结果非常吻合；r≥0.5 时，模拟结果可以接受。当Ens＝1 时，模拟结果和实测值完全吻合；0.5＜Ens＜1 时，模拟结果可以接受。当PBIAS ＝0 时，模拟结果非常好；PBIAS ＜0 时，模拟结果偏大；PBIAS ＞0 时，模拟结果偏小［24］。

年内分配完全调节系数Cr可表示径流量年内分配不均匀程度，值越大不均匀程度越大［25］。

式中：t为月份；R（t）为t月的径流量；为多年平均径流量。

3 结果与分析

3.1 无定河流域土地利用变化分析

CA－Markov 模型运行结果表明：模拟的2015 年土地利用数据和实测数据的Kappa系数为0.930 7，说明二者高度一致。所以用2005—2010 年土地转移概率矩阵和2010 年土地利用数据分别迭代15、20 次来模拟2025 年、2030 年土地利用数据。

本文用于土地利用分析的数据分辨率为1 km。由表4 可知：①无定河流域的土地利用类型主要是草地、耕地和未利用地，其中草地占流域总面积的40%左右，耕地占30%左右，未利用地占22%左右，其余为林地、水域和居民用地，共占8%左右；②草地面积先增加后减小，耕地先增加后减小而后保持维定，未利用地先减少后增加，居民用地持续扩张且偶有波动；③水域面积变化很小，说明流域25 a 内的水库、河流、滩地基本保持不变。由图3 可知：①草地主要集中在流域中上游地区，并在下游的耕地中有少量零星草地；②耕地主要集中在流域中下游；③未利用地大多分布在中上游地区，呈点状分布，并且2030 年未利用地在东北部大块状出现，应引起关注；④居民用地增加主要集中在河道附近和流域东北部。

表4 不同时期土地类型面积统计

图3 各土地利用类型分布情况

本文涉及4 种土地利用变化方式，见表5（记为C1、C2、C3、C4）。为探究C1～C4 的规律，建立4 个土地转移矩阵（见表6～表9）。综合图4、表6～表9 可知：①居民用地虽然在流域内占比小，但是增长幅度最大，C1～C4 分别增长110%、220%、120%、210%；②C1～C4的转移矩阵中，互相转化的土地利用类型的面积远小于没有发生转化的面积，说明2005—2030 年土地流转仅在小范围内进行；③C1，未利用地面积的1.2%和草地面积的0.6%共同转向耕地，使耕地增加1.3%；④C2，4.7%未利用地和4.3%耕地转向草地，使草地增加3%，同时1.27%的草地转出为居民用地；⑤C3，未利用地减少10.4%，转出为耕地、居民用地和草地；⑥C4，1.1%草地和2.1%未利用地转向耕地，使耕地增加2.9%，与C1 相比变化幅度更大，同时4.3%的草地转为耕地、未利用地和居民用地。

表5 不同土地利用变化方式

表6 C1 土地利用转移矩阵 km2

表7 C2 土地利用转移矩阵 km2

表8 C3 土地利用转移矩阵 km2

表9 C4 土地利用转移矩阵 km2

3.2 SWAT 模型校准和验证

根据本文研究需要，并结合前人研究经验［26－30］，采用1 km 的土壤数据集和土地利用数据集构建模型，共划分18 个子流域（见图5）和372 个HRUs。把2006—2013 年设定为模型率定期，把2014—2018 年设定为模型验证期。2005 年土地利用数据用于模型初始构建，利用白家川水文站还原后的逐月平均流量（以下称为还原值）进行校准。同时使用SWAT－CUP自动校准和验证参数，见表10。由图6 可知：①率定期和验证期的多年月平均流量模拟值和还原值的变化趋势有较好的一致性，拟合的线率分别为0.913、0.853 3；②率定期r、Ens、PBIAS分别为0.85、0.84、5.4%，验证期r、Ens、PBIAS分别为0.89、0.88、5.0%。综上，在无定河流域建立的SWAT 模型有良好适用性，可以用来研究C1～C4 对径流的影响。

图5 子流域分布

表10 参数率定结果

图6 模型率定期和验证期流量的模拟值和还原值

3.3 分析土地利用变化对径流影响的试验设计

气象数据不变，分析C1～C4 对径流的影响，具体试验设计见表11。其中，定义Q1～Q5 为情景1～5 下得到的径流模拟值的时间序列，R1～R5 为由流量时间序列得到的径流深，ΔR1～ΔR4 为C1～C4 引起的径流深变化。

表11 试验设计

3.4 不同土地利用变化方式对径流的影响分析

由图7 可知：①ΔR1～ΔR4 的变化趋势基本一致，且C1～C4 都会导致流域径流量增加，径流深多年总增加量从大到小排序为ΔR2＞ΔR4＞ΔR3＞ΔR1，说明C1～C4 均对流域径流起促进作用，虽然C1～C4 中没有进行流转的土地占各自土地利用类型面积的88%以上（见图4、表6～表9），且草地、耕地和未利用地面积各有增减，对径流有不同程度的促进或抑制作用，但是居民用地面积均增加，说明居民用地面积增加对径流的促进作用占主导地位，城市化程度高的流域，不透水地面的面积大，下渗量和下渗强度降低，属于超渗产流，且居民用地主要分布在河道两侧，从而使汇流时间变短，径流量增大，这与杨满根［31］在淮河流域得出的结论一致。②ΔR2 多年变化总量达9.54 mm，共计2.74 亿m3，占流域内现有水量的23.81%，说明C2 对径流的促进作用最大。变化量最小的是ΔR1，多年变化总量达3.97 mm，C1 中居民用地面积增加1.1 倍，说明C1 对径流的促进作用最小。③ΔR3 整体上略大于ΔR1，变化总量为4.40 mm；ΔR4 整体上略小于ΔR2，变化总量为8.40 mm，针对C4 使ΔR4 增加较大的情况，应引起关注。

图7 土地利用变化引起的径流深变化

为进一步探究ΔR的年内变化情况，本文从多年月平均径流深和年内分配完全调节系数（Cr）两方面入手。对R1～R5 中每年相同月份求平均值，而后分别和情景1 做差比较，得ΔR的年内变化，由图8 可知：①ΔR1～ΔR3 中，多年月平均径流深的年内增加趋势一致，呈先减小后增大而后减小的趋势。多年月平均径流深的增大量排序为ΔR2＞ΔR3＞ΔR1＞ΔR4，ΔR2 的年内变化总量达0.92 mm，说明C1～C3 对ΔR的各月流量均起促进作用，且C2 对各月流量的促进作用最显著，原因是C2 草地增加最多，由耕地和未利用地转入，使地表粗糙度增大，这与孙铁军等［17］的结论一致。且经过多年水保措施改造后土壤含水量高［17］，雨水落在草地上形成超渗产流；居民用地面积增加最多，形成不透水面积最大［31］，两方面原因共同促使各月径流深增加量最多。②ΔR4 中，1—6 月径流深减小总量共计0.26 mm，7—12 月径流深增大总量共计0.68 mm，说明C4 对上半年径流起抑制作用，对下半年径流起促进作用，且促进程度比抑制程度更高，原因可能是C4 未利用地和草地向耕地转移，使耕地面积大量增加，农作物在春天需要大量使用河水灌溉［7］，河水从河道进入土壤和地下，从而起到抑制径流的作用；下半年降水量增加，灌溉用水量少，从而起到促进径流的作用，针对草地和未利用地大量向耕地转移，导致上半年各月径流减少的情况应引起关注。③从年内分布的整体来看，ΔR1～ΔR4 的最大月份均在8 月，最枯月份均在3 月；C1～C4 均使7—12 月ΔR远大于1—6 月的，降水是本流域水资源的主要补给方式，且多集中在夏季和秋季，河流进入汛期，C1～C4 对汛期径流的影响更大。

图8 C1～C4 引起的径流深变化（ΔR）的年内分布

依据Q1～Q5 中每月Q值按式（7）和式（8）计算每年的Cr值，而后分别和情景1 的Cr值做差得ΔCr，绘制2006—2018 年ΔCr分布图。由图9 可知：①2006—2018 年的ΔCr1～ΔCr4 均大于0，变化幅度在0.14 之间，说明一定程度上C1～C4 均使径流深的年内分配不均匀，原因是C1～C4 各种土地类型面积的12%左右发生流转，造成下垫面在小范围内变化。②2006—2018 年的ΔCr1、ΔCr3 值均小于ΔCr2、ΔCr4，说明C1、C3 使径流深年内分配相对来说更均匀些，虽然一年内大部分降水集中在夏秋两季，使河水上涨，但是C1、C3 使耕地面积不同程度增加，无定河流域地处陕北，多种玉米和马铃薯［7］，一年内均需消耗不同程度的灌溉用水量，因此使径流深年内分配不均匀程度较小。③ΔCr1～ΔCr4 中，2014—2018 年的差值整体均呈现上升趋势，ΔCr2 的上升趋势最为显著，说明C2 对2014—2018 年的径流深年内分配不均匀程度影响最大。