李亚娇， 沈昞昕， 李家科， 李 莹， 刘易文， 周 翔
(1. 西安科技大学 建筑与土木工程学院， 陕西 西安 710054；
2. 西安理工大学 西北旱区生态水利国家重点实验室， 陕西 西安 710048)
景观格局是不同大小、形状和类型的景观斑块的空间排列。它是各种人类活动与自然因素在不同尺度上耦合的结果[1]。景观格局的变化影响着生态系统的水循环过程，主要体现在土地利用情况的改变以及在它的影响下形成的水文过程，如表面蒸散、土壤水分渗透和地下水循环等，从而影响径流的形成和转化，并进一步影响水文循环过程[2]。21世纪以来，随着经济的快速发展和人类活动的不断增强，大规模的垦荒、城市建设等措施改变了各流域内的景观格局。这对流域水资源的形成、转化和时空分布产生了重大影响，导致水土流失、干旱、洪涝灾害等极端水文事件频发[3]。因此，研究水文特征对景观格局变化的响应关系逐渐成为近年来全球生态环境研究的热点和前沿问题，对流域水资源管理具有重要的现实意义和科学价值[4]。
目前，关于人类活动的水文效应研究主要集中在各要素对径流的影响上[5-6]，其方法包括水文特征参数时间序列法[7]、实验流域法[8]和流域水文模型模拟方法。其中，流域水文模型模拟方法考虑了地理因素的空间差异，能有效模拟自然界复杂的水循环[9]。因此，通过构建流域水文模型，特别是基于物理基础的分布式水文模型，可以有效分析景观格局变化的径流效应。与其他水文模型相比，SWAT模型具有较强的物理机制，已被广泛用于模拟土地利用和气候变化背景下的流域水文特征响应分析[10-11]。汉江是长江最大的支流，其上游的丹江口水库是我国南水北调中线工程的水源地，兼顾重大民生工程，具有重要的战略意义。近年来，随着人类活动影响的加剧，汉江流域出现了植被退化、土壤侵蚀等一系列生态环境问题。因此，研究该地区景观格局变化对径流、泥沙的影响，对保障该地区水生态安全、促进水土资源的合理开发利用和保护具有理论和实践意义。此前，叶加俊[12]和陈欢[13]等利用水文模型，通过设定不同气候、不同景观配置，对汉江上游各水文要素的变化特征进行了研究。然而，该流域径流及泥沙对景观格局变化的响应机制尚不明确。因此，基于前人研究，利用SWAT模型对汉江上游流域1980—2017年月径流及月泥沙变化过程进行模拟。基于该时间段内景观格局及水沙变化情况，探讨该流域水沙变化对景观格局的响应机制，并找出影响水沙变化的关键因素，从而为该流域水土资源的合理开发利用和水环境保护政策的制定提供科学依据。
汉江发源于陕西省汉中市宁强县，在陕西境内自西向东穿行于北部的秦岭山地与南部米仓山、大巴山之间，地形起伏大，峡谷盆地交错分布。在陕西境内属亚热带季风性气候。汉江支流众多，大多属于山溪性河流，河床狭窄，水流湍急。主要汇入支流有沮水、玉带河、漾水河、褒河、湑水河等。当地降水主要集中在夏秋两季，暴雨、洪涝灾害频发。汉江不仅是我国南水北调中线工程的重要水源地，也是引汉济渭工程的引水水源地，是我国的重要水源保护区之一。随着长江经济带城镇化进程的推进和南水北调工程的实施，汉江上游流域的生态环境受到严重威胁，植被退化频现，泥石流等地质灾害频发，造成了严重的水土流失及荒漠化现象。本研究主要关注洋县断面以上的汉江流域上游地区(图1)，该地区是丹江口水库的水源地，生态与经济发展的战略地位显著[14]。
图1 汉江流域上游位置示意图Fig.1 Schematic diagram for upstream location of Hanjiang River Basin
2.1 数据来源与处理
本研究采用的汉江上游流域1980年、1990年、2000年、2010年、2017年的土地利用数据来源于中科院资源环境科学与数据中心 (http://www.resdc.cn)。该土地利用数据包含耕地、林地、草地、水域、建设用地和未利用地6个一级类别以及水田、旱地等25个二级土地利用类别。统一采用WGS_1984 地理坐标系，经检验，各期数据精度均达到85%以上。利用ArcGIS软件将所有二级地类重分类合并成6个一级地类进行评价。土壤数据取自寒区旱区科学数据中心(http://data.casnw.net)世界土壤数据库(HWSD)的中国土壤数据集(v1.1)；
数字高程数据(DEM)采用地理空间数据云(https://www.gscloud.cn)中的地理高程数据，分辨率为30m；
气象数据取自中国气象数据网(http://data.cma.cn)的中国地面气象日值数据集(V3.0)，时间尺度为1980—2017年；
水文数据(包括径流量及输沙量数据)来自汉中市水文局。
2.2 分析方法
2.2.1SWAT模型构建
利用所需数据构建SWAT模型后，使用SWAT-CUP 软件，以月径流量、月输沙量为目标，对模型进行率定和验证，一般选取决定系数(R2) 和Nash-Sutcliffe效率(ENS)来评价模型精度[15]。使用汉中市水文局提供的武侯镇站和洋县站的径流量、输沙量实测数据进行对比计算。其中，径流量数据时间尺度为2003—2017年，输沙量数据时间尺度为2003—2015年。分别将1980年、1990年、2000年、2010年的土地利用数据代入模型，模拟汉江上游流域1980—2017年径流量及输沙量特征。对比观测值与模拟值，参考其他学者研究所得出的最佳模拟参数值和推荐参数范围值进行重复率定，直至模型精度合理。
2.2.2景观格局特征分析
子流域包含最全面的景观特征，因此是分析包括各景观类型占比及景观格局指数在内的景观格局与水文特征变化的最适尺度[16-17]。
因此，基于DEM数据，同时结合当地地理特征，使用SWAT模型的水系提取功能进行子流域划分，流域出口设定为洋县水文站，最终将研究区划分为25个子流域，如图2所示。
图2 汉江流域上游子流域划分Fig.2 Sub watershed division in the upper reaches of Hanjiang River Basin
参考国内外关于水沙变化对景观格局响应的研究成果[18-20]，在景观组成上选取耕地、林地、草地、水域、建设用地、未利用地6种土地利用类型；
在景观配置内选择表征斑块破碎程度的斑块密度(PD)、表征景观斑块间连通性的聚合度指数(AI)、表征景观斑块复杂程度的景观形状指数(LSI)、表征景观斑块聚集程度和连通性的蔓延度指数(CONTAG)和结合度指数(COHESION)以及表征景观异质性的香农多样度指数(SHDI)。
这些景观指标简单且通用，能较为全面地反映各子流域的景观特征。各景观格局指数的生态学意义参见文献[21]。
2.2.3冗余分析
冗余分析(RDA)目前被广泛用于确定不同指标之间的相互关系。
本文以25个子流域的多年径流量、输沙量及多年景观格局指数为样本，对各景观格局指数与流域水沙变化间的相关关系进行冗余分析。
冗余分析图中，箭头线的长度代表景观格局对水沙变化的影响程度，箭头线越长，影响越大；
景观格局箭头线方向与水沙箭头线方向的夹角代表两者的相关程度，夹角大于90°，呈负相关关系，小于90°，呈正相关关系，等于 90°，则不存在关系[22-23]。
3.1 流域产水产沙时空分布格局
3.1.1模型率定及验证
将1980—2017年的日降水、温度、相对湿度和太阳辐射等气象数据以及土壤、土地利用数据等导入SWAT模型。利用SWAT-CUP 2012的SUFI-2方法对模型进行率定。
首先，采用全局灵敏度分析方法选择模型参数进行率定；
然后，利用实测数据对径流量和输沙量进行标定。通过决定系数(R2)、Nash-Sutcliffe效率(ENS)来评估率定和验证的有效性。
在率定和验证后，基于4期土地利用情景运行模型，如基于1980年土地利用情况计算1980—1990年径流量、输沙量数据，以此类推，得到了1980—2017年径流量及输沙量的模拟结果[24]。
在模型率定和验证阶段，将模拟结果与两个站点实测的月径流量、月输沙量值进行对比。图3为武侯镇站及洋县站径流量及输沙量的拟合过程；
表1为评价结果。各站点径流及泥沙的R2>0.7，ENS>0.7，因此，模型适用于模拟汉江上游流域的径流泥沙状况。
图3 武侯镇及洋县站模拟值与实测值对比Fig.3 Comparison between simulated and measured values at Wuhou town and Yangxian stations
表1 SWAT模型校准与验证结果Tab.1 Calibration and inspection result by the SWAT model
3.1.2流域水沙时间分布特征
汉江上游受季风型大陆性气候影响，呈现四季分明的特征，其北部为山地暖温带温和湿润气候区，南部为北亚热带温热湿润气候区；
雨量充沛，具有降水分配不均、年际变化大的特点[25]。根据流域出水口洋县站径流量及输沙量数据，分析流域水沙时间分布特征可知，7月～9月为该流域产沙期，河川输沙量占全年总量的73.9%，径流量占全年总量的64%。经计算，该流域径流的年内分配不均匀系数(Cv)为0.83，说明径流年内分配较不均匀；
泥沙的年内分配不均匀系数为1.06，说明泥沙年内分配极不均匀。图4为流域水沙及降雨量变化趋势图。总体而言，径流量和输沙量总体呈增加趋势，降雨量呈减少趋势；
同时，径流量与输沙量时序保持一致。经spearman相关分析检验，径流量与输沙量相关性为0.875，降雨量与径流量相关性为0.561，均通过了p<0.01的显著性检验，呈显著性相关。体现了研究区降雨径流的高度一致性及地表河道内多水多沙、少水少沙的特点[26]。
图4 汉江流域上游水沙变化趋势Fig.4 Variation trend of runoff and sediment in the upper reaches of Hanjiang River Basin
3.1.3流域水沙空间分布特征
为了更好地判断径流量及输沙量对不同景观特征的响应，并消除子流域大小对区域径流量及输沙量的影响，基于SWAT模型的模拟数据，分别计算25个子流域1980—2017年的多年平均径流模数、输沙模数和降雨量，如图5所示。结果表明，该流域多年平均径流模数为154.18 L·s-1/(km2·a)，多年平均输沙模数为61.27 t/(km2·a)，由《土壤侵蚀分类分级标准》(SL 190—2007)[27]可知，汉江流域上游的土壤侵蚀强度属于微度。此外，径流模数、输沙模数及降雨量分布不均，呈现明显的“南北高、东西低”的分布特征。根据李朝月等[28]的研究，地形因子(平均坡度、坡长、高程等)与径流模数及输沙模数存在显著正相关，其原因在于随着坡度增大，降雨入渗减小，地表径流增大，进而导致径流冲刷能力增强，土壤泥沙被侵蚀、输移的强度增大。
图5 多年平均径流模数和输沙模数的空间分布Fig.5 Spatial distribution of annual average runoff modulus and sediment transport modulus
同时，大量研究表明，影响径流量的最直接因素为降水，而降水与径流的演变规律保持高度一致，因而子流域的径流量及输沙量受降雨量的影响显著。并且，该流域内北侧为秦岭南麓，南侧为米仓山北麓，土地利用类型以林地为主，其中针叶林居多，阔叶林、灌木林较少。针叶林林冠小，降水对地表更容易形成直接冲刷，所以流域的输沙量有明显增大[29-30]。
对各子流域径流模数、输沙模数进行统计后发现，流域内各子流域的径流量及输沙量分配不均，具有较强的空间分布特征，这与研究区的地形因子、气候特征及土地利用方式等因素的分布方式密切相关。
3.2 土地利用特征分析
3.2.1土地利用变化分析
汉江上游流域总面积约为15 940 km2，景观类型以耕地、林地、草地为主，三者面积之和占比95%以上；
未利用地面积过小，在0.05%以下(图6)。以2017年土地利用数据为例，耕地和建设用地主要分布在流域的中部和南部地区(子流域10、18、21和22)；
林地主要分布在流域北部和南部地区(子流域1、2、3、5和6等)，主要为各支流水源区；
草地面积在各子流域均占比20%以上(图7)。在该流域两侧的秦岭山脉南麓和大巴山北麓，林地草原密布，典型树种包括侧柏、马尾松、云杉，还有柑橘等农业经济作物，草地则以稀疏覆盖草原地类为主。耕地集中分布于广袤的汉中盆地内，以水田为主，当地农业生产水平较高，城镇化水平低，从而对自然景观的破坏程度小；
但近十几年来，城镇化的快速发展也导致建设用地面积显著增长，且在空间上呈现集聚现象，主要集中分布于汉江沿岸。水域占地面积逐渐增长，但直到2017年，占地仍不足1%。尽管当地河网密布，但经实地调查，多为村镇内的人工沟渠及诸多河道宽度较小的支流。这些支流与村镇内的人工沟渠共同构成了汉江上游流域河网。
图6 1980—2017年土地利用变化情况Fig.6 Land use change from 1980 to 2017
图7 各子流域景观类型占比Fig.7 Proportion of landscape types in each sub watershed
据资料显示，汉中市为省内主要矿产资源聚集地之一，目前已发现矿产92种，引来众多中小企业从事矿产冶炼开采。尽管环保部门多次督促其进行地表整治，但植被覆盖率仍然较低。经开发后形成的裸地零星分布于山地间，与河滩采砂的滩涂地一同构成未利用地。
3.2.2土地利用转换分析
表2为土地利用类型转移矩阵。景观格局组成结构的动态变化最直接的体现就是土地利用类型的变化转移。最显著的转移为各地类向建设用地的转移，其中，耕地向建设用地转移的面积占转移总面积的44.15%，未利用地向建设用地转移的面积占28.09%。除此之外，土地利用类型的转移主要发生在耕地与草地之间，林地与草地之间；
其中，耕地与草地间的相互转化面积均超200 km2，林地与草地间的相互转化面积均超100 km2。
表2 1980—2017年土地利用类型转移矩阵Tab.2 Land use type transfer matrix from 1980 to 2017 单位：km2
结合土地利用类型面积及实际情况分析可知，当地主要依赖第三产业，但随着21世纪的到来，相当一部分人口流入城镇，沿河耕地也在此时被开发为建设用地；
同时，部分草地被开垦为耕地，以维持当地农业经济发展。
3.3 景观格局特征分析
对流域1980—2017年的景观格局指数进行计算，结果如表3所示。从景观水平上看，1980—2017年，PD值总体减少，其值经历了先增加再减少随后再增加的过程；
20世纪人口增长及社会经济发展所引起的城镇化以及21世纪后积极实施退耕还林等生态保护措施是导致斑块密度先增加后降低的主要原因。从景观形状上看，LSI值逐渐增加，说明景观形状趋于复杂；
频繁的人类活动导致了耕地、建设用地等对其他土地利用类型的侵占，这使得景观斑块形状愈加复杂和不规则。从景观集聚程度上看，COHESION、AI、CONTAG值逐渐减少，表明研究区内景观不同斑块分布的聚集程度前期基本保持稳定，在21世纪初则逐渐减弱，离散程度随之增大；
主要原因在于城镇化使农村居民点聚集的同时，令其他建设用地诸如工矿、采石场等的分布也更加分散；
同时，植树造林使得林地、草地斑块入侵耕地、建设用地斑块，使得聚集程度有所下降。SHDI值也表现出相似的变化趋势，同样印证了景观形状越来越复杂，但变化并不明显。总体而言，1980—2017年间，研究区内景观特征表现出斑块数减少、破碎化程度增加、景观分割度升高、聚集度下降和景观复杂度上升的趋势[31]。
表3 流域景观指数变化特征Tab.3 Variation characteristics of watershed landscape index
3.4 对景观格局的响应分析
3.4.1景观格局与流域水沙间的冗余分析
利用Canoco 5软件对景观组成因子、景观配置因子、土地利用类型与流域水沙间的相关关系进行冗余分析。统计结果表明，所选的景观格局指标对流域水沙变化的总解释率较高，达到了98.03%。
由图8可知，COHESION、CONTAG、AI均与径流量、输沙量呈正相关，其中COHESION与径流量及输沙量的相关关系显著；
这可能与当地以水田为主有关，相同斑块类型间良好的连通性及水田与林地等不同斑块类型间的紧密连接，使得空间距离较小，有助于调节水文过程[32]。LSI与径流量呈正相关关系，与输沙量呈不明显的正相关关系；
表明随着城市化进程的推进，对自然景观的改造力度加大使得流域内自然景观的水文调蓄能力(如降水截留等)减弱，景观格局形状逐渐复杂化，更易造成水土流失[33]。SHDI与径流量呈正相关，与输沙量呈负相关；
据文献[28]、[34]的研究结果，景观格局多样性的提高能削弱土壤侵蚀的形成和运输能力。PD与径流量相关性不强，而与输沙量呈显著负相关；
主要原因在于该流域内主要发生的是林地、草地被耕地及建设用地入侵造成的景观斑块破碎，而林地、草地景观斑块的破碎化不利于斑块景观内部的物种繁衍及景观整体的生态维护，进而增加了水土流失风险。总之，流域内景观格局的斑块间越连通、团聚性越好并且斑块形状越复杂，径流量及输沙量的增幅越大，越易造成水土流失；
当景观格局多样性提高时，径流量随之增大且输沙量减少，对水土保持具有显著的积极意义。
注：FARM为耕地占比；
GRASS为草地占比；
FOREST为林地占比；
CON为建设用地占比；
UNUSED为未利用地占比。图8 流域水沙变化与景观格局的冗余分析(RDA)Fig.8 Redundancy analysis of watershed runoff and sediment change and landscape pattern index (RDA)
耕地面积与径流量呈正相关，与输沙量呈不显著正相关；
这是由于人类活动范围扩大，开垦了新的耕地，使得原本的土壤结构被破坏，对径流下渗率、土壤含水量和地下水补给量等水文要素造成了影响，从而导致地表径流增加；
同时，耕地的植被覆盖度普遍较低，易发生水土流失[35]。林地面积与径流量呈不显著负相关，与输沙量呈负相关，说明林地面积的增加对集水固沙有积极作用；
这是由于林地植被根系发达，能够扎入深层土壤，且拥有强大的蒸腾能力，其下渗量与蒸发量均较大，能汲取较多土壤中的水分，从而调节地表径流；
土壤含水量的增加也会对泥沙的产生与输移造成阻滞，同时，树冠对降雨存在拦截作用，雨滴降落后打击在地层表面时动能已经减弱，降低了地表径流流速，同时还可防止土壤溅蚀，提高了土壤对径流侵蚀的抵抗力[36]。草地面积与径流量呈不显著正相关，与输沙量呈负相关；
虽然草地根系不发达，但在土壤和根系的双重作用下，也可在一定程度上增加降雨入渗，提高土壤含水量和地下水补给量，从而增加产流；
同时，存在于土壤表层的根系也提高了土壤对径流侵蚀产沙的抵抗力，从而减少了产沙[37]。建设用地面积与径流量呈不显著正相关，与输沙量呈负相关；
原因在于随着当地经济的持续发展及城市化进程的加快，建设用地不断扩张，各类用水使得水资源开采量加大，一定程度上减少了径流量；
同时，大量农村人口聚集在城镇，城镇化使得不透水面面积大幅增加，清除了土壤侵蚀的来源；
此外，相关工业设施及居民市政配套的建设需要大量河沙，这也使得流域水沙迅速减少[38-39]。未利用地面积与输沙量呈正相关，与径流量呈不显著正相关关系；
流域内的未利用地多为矿山开采后的荒地，荒地普遍植被覆盖度较低，地表结皮较多，造成地表径流难以入渗，其径流量和产沙量最大[40]。
3.4.2景观类型和格局指数与水沙变化的逐步回归分析
基于子流域进行径流量、输沙量与景观组成、景观配置之间的冗余分析，通过逐步回归分析(表4)，在各景观格局指数中找到关键指数。同时，对各参数进行共线性检验[41]，若共线性严重，则表示回归模型中的解释变量间存在高度相关关系，这将使模型难以准确估计。通常利用容差及方差膨胀系数(VIF)来衡量回归模型中的共线性程度，当容差≤0.1或VIF≥10时，代表变量间共线性严重。结果表明，各参数的容差介于0.12～0.85之间，VIF介于1.0～9.6之间，显然，各变量间不存在严重的共线性。进入径流量回归模型的景观配置因子有LSI和CONTAG，景观组成因子有耕地占比(FARM)与建设用地占比(CON)；
其中，LSI和CONTAG达到显著性水平(p<0.05)。进入输沙量回归模型的景观配置因子有LSI和AI，景观组成因子为草地占比(GRASS)；
且均达到显著性水平(p<0.05)。在径流量回归模型中，LSI、CONTAG和建设用地为正效应，耕地为负效应；
在输沙量回归模型中，LSI、AI为正效应，草地为负效应；
其中LSI都是最主要因子。逐步回归分析结果表明，景观复杂程度是促使水沙变化最重要的景观格局指数，其值越大，景观形状越复杂。人类活动越频繁、流域内景观越破碎、景观类型越复杂，水土流失越严重。因此，应避免人类活动对大范围林地、草地等景观的破坏，并提高不同景观间的连通性，降低景观破碎化程度，从而有效防止水土流失。
表4 景观格局与水沙变化的逐步回归分析Tab.4 Stepwise regression analysis of landscape pattern and water sediment changes
本文从子流域尺度分析了汉江流域上游景观格局与水沙变化的关系，也从一定程度上揭示了水土流失发生的原因。在水土流失综合治理的过程中，调整土地利用结构要遵循治理与开发相结合的基本原则；
同时，要积极实行退耕还林等生态保护措施，还要注意当地林地、草地的形状、位置特征等景观格局要素的变化。虽然子流域尺度上的景观格局与水沙变化分析能够进一步细化流域景观格局改变的具体空间区位，但其他学者[42]的研究表明，当地的人类活动，如对河道径流的拦截等，也会影响流域内的水循环规律，进而对径流泥沙的输运产生影响。因此，应尝试从多个尺度进行景观配置，从而为水资源统筹规划与可持续发展提供参考。
本文利用 SWAT 模型还原了汉江流域上游1980—2017年径流量及输沙量的变化过程，并分析了38年来该流域水沙变化对景观格局的响应关系。
1) 利用SWAT模型模拟汉江流域上游径流量及输沙量的变化，率定结果满足要求，说明该模型在汉江流域上游的适用性较好。1980—2017年，汉江流域上游径流量及输沙量时空分布不均，具体体现在径流及泥沙年内分配不均，存在变异性；
在年际变化上，年径流量及输沙量都呈缓慢上升趋势。各子流域间的径流量及输沙量分配不均，具有较强的空间分布特征，这与研究区的地形因子、气候特征及土地利用方式等因素的分布方式密切相关。
2) 汉江流域上游的主要土地利用类型为耕地、林地和草地。21世纪以前，由于人类活动范围的扩大，大量林地及草地转化为耕地及建设用地。21世纪以来，当地采取了退耕还林等生态保护措施，效果显著。因此，土地利用类型的转换也主要发生在这三者之间。同时，流域内各景观格局指数基本稳定，景观总体结构未发生重大变化。但也表现出斑块数减少、破碎化程度增加、景观分割度升高、聚集度下降和景观复杂度上升的特征。
3) 景观组成因子及景观配置因子与流域水沙间的相关关系各不相同。从整体上来看，植被覆盖率提高可有效减少径流量及输沙量，从而防止水土流失现象的发生。景观格局的斑块间越连通、团聚性越好并且斑块形状越复杂，径流量及输沙量的增幅越大，越易造成水土流失；
当景观格局多样性提高时，径流量随之增大且输沙量减少，对水土保持具有显著的积极意义。经逐步回归分析发现，景观复杂程度(LSI)是影响水沙变化最重要的景观格局指数。人类活动越频繁、流域内景观越破碎、景观类型越复杂，水土流失越严重。因此，应避免人类活动对大范围林地、草地等景观的破坏，降低景观破碎化程度，提高景观间的连通性，从而有效防止水土流失。