土地利用空间格局和类型组成不仅对坡面水文过程具有直接的控制作用^[1]，而且也是影响流域内水分循环和水体质量的重要因素^[2]。因此水土流失与非点源污染负荷不仅与流域土地面积结构有关，同时与土地利用格局密切相关^[3]。土地利用中植被空间分布格局能够改变径流和泥沙运移路径的连通性，在流域的不同坡度设立植被条带，控制水土流失的效率有所不同^[4]。如林地、耕地等自然与半自然斑块破碎度越低、面积比越大，对流域内非点源污染的缓冲与抑制作用越强^{[5, 6]}。

水文模型模拟是研究水土流失和非点源污染的重要手段，SWAT模型是美国农业部(United States Department of Agriculture，USDA)农业研究署(Agricultural Research Service，ARS)开发的基于GIS平台的流域分布式水文模型，能够模拟较长时段内复杂大面积流域的水文过程，广泛应用于分析和预测因土地利用变化等引起的径流、泥沙及非点源污染的变化。Mahesh Sahu等通过对Walnut Creek流域设置情景模拟，假设研究区某坡度带土地利用类型发生改变，利用SWAT模型分析和预测景观格局变化下的水文变化过程^[7]。在国内外相关研究中，往往默认SWAT模型能够模拟不同土地利用空间分布变化对水文过程的影响^{[8, 9, 10, 11, 12]}，在此基础上分析不同坡度种植植被条带对水土流失和非点源污染的控制效率，从而为流域土地利用管理决策提供重要支持。但是，由于SWAT模型在一定空间尺度上概化了土地利用空间差异^[13]，因此有必要了解模型对土地利用景观格局变化的敏感性，从而判断模型模拟各种情景假设的精度。

本文以丹江口库区的老灌河流域为例，模拟两幅土地利用图，利用Fragstats 3.3计算和分析两者在景观格局指数方面的差异；并分别将其代入 SWAT模型，通过对比结果和相关关系分析模型对不同景观空间配置的敏感性；对比基于真实与模拟土地利用图的SWAT模型校准过程与参数变化，分析模型精度，针对模拟土地利用情景假设提出建议，增加模型对不同坡度土地利用景观空间变化的敏感性，从而提高模型精度。

老灌河流域位于东经110°49′ — 111°50′，北纬33°18′— 33°59′之间，总面积2612 km²(图 1)，是河南境内一支重要的丹江口水库源头；河流落差大，水力资源丰富。年平均气温12.1 ℃，年日照2103 h，年均降水量862.8 mm，无霜期198 d。该区地形复杂，北部是海拔高、坡度大的中低山地，南部是灌河谷地，两侧为低山丘岭区。

图 1 老灌河流域示意图 Fig.1 Map of Lao Guanhe watershed

图选项

SWAT作为分布式模型，通过设定阈值将研究区划分为多个具有一定拓扑关系的子流域，每个子流域包含至少一个以上的水文响应单元(HRU)^[14],每个水文响应单元包含唯一的土地利用类型、土壤类型和坡度等级且假定有统一的水文过程^[15]。因此，景观格局变化在SWAT模型中主要表现为HRU在形状、面积、分布等方面的变化；通过对每个HRU的负荷进行计算，汇总在一起确定子流域的负荷。

SWAT模型中主要包含有水文过程子模型和土壤侵蚀子模型。其中对于水文循环的陆地阶段，模型中采用水量平衡表达式进行模拟^[14]：

点击浏览公式

式中,SW_t为土壤最终含水量(mm)；SW₀为土壤初期含水量(mm)；t为时间步长(d)；R_day为第i天的降水量(mm)；Q_surf为第i天地表径流量(mm)，与流域各土地利用类型面积和平均坡度直接相关；E_a为第i天蒸散发量(mm)；W_seep为第i天透过土壤层的渗透量和侧流量(mm)；Q_gw为第i天的地下水含量，mm。

SWAT模型采用修正的通用土壤流失方程来分别计算每个HRU的土壤侵蚀量^[14]，公式如下:

点击浏览公式

式中,m_sed为单个HRU的总土壤侵蚀量(t)；Q为地表径流(mm/h)；q_feng为洪峰径流(m³/s)；A_HRU为HRU面积(hm²)；K为土壤可蚀性因子；C为植被覆盖与管理因子；P为水土保持措施因子；LS为地形因子，与HRU平均坡度直接相关；CFRG为粗糙度因子。

景观格局指数是度量景观格局的指标，是对景观格局信息的高度浓缩，主要包括类型组成和空间配置两个方面。景观类型组成不涉及空间因素，主要包括土地利用类型的面积与比例；景观的空间配置则主要对土地利用类型的分布情况、空间排列以及各斑块的形状和相对位置等进行刻画。本文选取土地利用类型组成指数包括流域面积(TA)和最大斑块面积(LP)，空间配置指数包括斑块密度(PD)、景观丰度(PR)、边缘密度(ED)、总边缘长度(TE)和最大形状指数(LSI)^{[16, 17]}。

SWAT模型的输入数据包括数字高程图(DEM)，土壤数据，土地利用数据和气象数据。

DEM是对地形地貌的一种离散的数学表达，数据来源于中国科学院计算机网络信息中心国际科学数据服务平台，精度为30 m(图 2)。

图 2 研究区DEM Fig.2 DEM of study area

图选项

土壤数据是生成水文响应单元HRU的重要基础，通过栾川、西峡、卢氏3县的1 ∶ 5万土壤图的数字化得到78种土壤类型，最后根据河南省土壤分类标准合并生成65种土壤类型(图 3)，其中编号42的中砾质薄层淡岩黄砂石土面积占流域23.56%，为研究区主导土壤类型。各类土壤的物理属性数据通过查阅上述3县的土壤普查报告得到。

图 3 土壤分布图 Fig.3 The distribution of soil map

图选项

实际土地利用数据(Land Use/ Land Cover data,LUC)通过对2000年Landsat-TM影像监督分类(数据来源于中国科学院计算机网络信息中心国际科学数据服务平台，精度30 m，日期：2000年4月12日，条带号：125，行编号：37)，并结合国土资源部门土地利用现状图生成，得到10种土地利用类型(图 4)。

土地利用类型面积和平均坡度直接影响模型各子流域的产流、产沙结果，因此分析SWAT模型对不同景观空间配置的敏感性应首先排除上述因素对模型结果的影响。研究利用ArcGIS 9.3的Create Random Raster(整型随机分布)功能随机生成两幅土地利用图LUC-A和LUC-B，设置LUC-A和LUC-B各子流域各土地利用类型面积与2000年实际土地利用一致，并保持两幅模拟图之间各子流域中每种土地利用类型的平均坡度对应相同。LUC-A和LUC-B最小象元分别为90 m和30 m，保证两幅模拟图在景观空间配置上存在较大差别。LUC-A和LUC-B在研究区相同区域的局部放大图如图 4所示，分辨率均为30 m，可以直观清晰地看出二者在景观格局上存在较大差别。

图 4 实际土地利用图和模拟土地利用局部放大图 Fig.4 Land use map and locally enlarged simulated land use map LUC: Land Use/ Land Cover data

图选项

气象数据来自栾川、西峡、卢氏三县气象站，时段为1971年—2000年，主要包括：逐日降雨量、最高气温、最低气温、太阳辐射强度、相对湿度和风速数据，逐日太阳辐射强度通过日照时间计算^[18]得到。水文数据来自河南省西峡站1993—2000年逐日观测数据。

表 1为LUC-A和LUC-B各子流域景观格局指数。从表中看出，在两幅模拟生成的土地利用图中，景观类型组成指数的流域面积(TA)和最大斑块面积(LP)在各子流域中都是相同的；空间配置指数的斑块密度(PD)，景观丰度(PR)，边缘密度(ED)，总边缘长度(TE)和最大形状指数(LSI)在各子流域中差别明显，符合研究对两幅模拟土地利用图的要求。其中，在相同子流域中，LUC-B的PR值与PD值更大，反映出模拟的土地利用破碎度更高，说明在该模拟图中各土地利用类型分布更加破碎；LUC-B的ED值较LUC-A更大，表明异质景观斑块之间单位面积的边缘长度更长；由于除13号子流域，各子流域占主导地位的景观类型均为林地，更容易截留泥沙与氮磷，从而降低污染物汇入河道的风险。LUC-B中较大的景观类型指数LSI值不仅表明景观的斑块分布更加破碎，同时景观的类型分布也较为破碎。由于在景观设计中剔除了坡度影响，根据以往研究^{[1, 2, 3, 4, 5, 6, 19]}，在其他因素相同的情况下，基于不同的土地利用景观类型，流域产流、产沙和非点源污染的结果应该有所差别。

由于真实的土地利用图与模拟图相比，各景观类型的平均坡度不同，SWAT模型结果受其他参数影响较大，不能直接对比，因此首先对比模型未校准前模拟土地利用图的结果。

利用土地利用模拟图LUC-A和LUC-B分别建立SWAT模型。根据数字高程模型，研究区被划分为13个子流域(图 1)。然后通过在各子流域内部分别叠加模拟土地利用图(LUC-A、LUC-B)和土壤类型图生成HRU，最后输入气象数据对每个HRU分别模拟，经河网汇集得到整个流域2000年的产流、产沙和非点源污染(有机氮)结果。为保证模拟土地利用图的每个斑块都参与运算，设置该阈值为0%，同时选择主导土壤类型代表整个流域，以排除土壤因子对HRU划分的影响。采用平均精度和Nash-Stucliffe(E_ns)确定性系数对两个模型模拟结果的吻合度进行评价，平均精度和E_ns 越接近1，两个模型模拟结果越接近，当E_ns >0.75,可以认为模型模拟效果好。模型运行结果分别如下(表 2)。

从表 2发现，未校准前，基于两幅模拟土地利用方式的SWAT模型模拟结果差别很小，对各子流域径流、泥沙及非点源污染的平均精度分别为1、0.97和0.98，E_ns均为1。两个模型模拟的结果几乎是相同的。两个模型各子流域间产沙、有机氮的微小相对误差与研究区主要土地利用类型林地的平均坡度差别呈现显著的正相关性。因此可以发现，当输入的土地利用方式只有景观空间配置差异时，SWAT模型模拟结果几乎不发生改变。

分别将两个模型各子流域模拟结果与对应的子流域景观格局指数进行spearman相关性分析，结果如表 3。

可以看出，同一景观格局指数与基于不同土地利用模拟图的SWAT模拟结果间的关系值几乎是相同的。景观类型组成的子流域面积TA和子流域最主要景观所占面积LP与SWAT模拟结果之间有着较为显著的相关性，但在景观水平下，景观空间配置的主要参数与SWAT模拟结果之间没有显著的相关关系。SWAT模型不能很好地反映景观空间变化对产流、产沙和非点源污染的影响。

进一步对比LUC-A、LUC-B模型HRU划分可以看出，各子流域基于不同景观格局类型生成的HRU个数都是相等的(表 2)。虽然土地利用类型都是由一个或多个空间上不连续的斑块组成，但每个斑块在SWAT模型中并没有形成独立的HRU，模型通过聚集相同的土地利用类型和土壤类型构成唯一的水文响应单元，从而模糊了土地利用类型下各斑块的形状、地理位置等空间配置。以利用SWAT模型模拟土地利用情景模拟为例，当增加某坡度带如河流岸边带林地面积时，HRU数量不会改变，模型更多是反映代表林地的HRU因面积增加，平均坡度降低而减少非点源污染的结果，忽视了林地不同的景观空间位置对非点源污染的调节作用，从而使模拟的结果很难真实反映不同坡度带林地的实际作用。

因此，当研究者基于不同坡度假设土地利用情景模拟时，应利用SWAT模型HRU划分中Mutiple Slope选项，将流域划分为多个坡度类，从而增加子流域HRU的数量，缩小每个土地利用类型斑块的坡度与所属坡度带平均坡度差值，当增大或减小不同坡度的土地利用面积，能够显著地修改HRU平均坡度，增加了模型对不同坡度土地利用的敏感性，从而有利于对岸边带等特殊坡度的模拟。

利用SWAT模型，分别对LUC和LUC-A进行模拟，在生成HRU时，定义所有阈值为0，保证所有的土地利用类型和土壤类型参与运算。经模型敏感性分析，SWAT校准需要调整6个与径流模拟有关的参数，分别为CN₂(径流曲线数)，SOL_AWC(土壤有效水容量)，ESCO(土壤蒸发补偿系数)，GW_REVAP(地下水revap系数)，REVAPMN(发生revap的浅层含水层水位阈值)和GWQMN(发生回归流所需的浅层含水层水位阈值)。利用西峡水文站1993—1996年实测径流数据校准两个SWAT模型，1997—2000年实测径流数据用于模型验证。结果表明，当两个模型的6个校准参数都相等时，基于LUC和LUC-A的SWAT模型径流模拟在校准期分别为0.88和0.86，验证期分别为0.90和0.87，2个模型模拟的径流结果与观测值的拟合度都很高(图 5)。即便是不符合真实流域景观特征的LUC-A，SWAT模型也能够在验证期内较为准确地模拟老灌河流域的径流状况。

图 5 西峡站点月平均模拟值与实测值的比较 (校准期：1993—1996，验证期：1997—2000) Fig.5 Simulated and observed runoff on a monthly step at Xixia station (Calibration stage: 1993—1996,validation stage: 1997—2000)

图选项

由于缺少实测泥沙数据，研究设置基于实际土地利用(LUC)的SWAT模型产沙数据为校准数据，对基于LUC-A的SWAT模型进行校准，经敏感性分析确认，需要对以下参数进行调整，结果如表 4。

从表 4可以看出，在基于LUC-A的模型中，河道的单位降雨侵蚀力造成的土壤侵蚀量(Ch_erod)和河道水流冲蚀之比(Ch_cov)以及USLE方程中的覆盖和管理因子(USLE_C)和河道挟沙能力计算的待定幂指数相比目标模型都有不同程度的降低。即，模型通过其他方面的弥补来缩小因景观格局指数变化造成的产沙差别，从而造成上述因子可能与实际的进一步背离。

同时还可以发现，基于LUC-A的模型对LUC的产沙结果拟合度很高：E_ns达到0.93。如果后者能够通过调整泥沙参数反映真实的产沙情况，那么前者也能够经参数率定使模型反映该流域同时期的产沙状况；但两个模型在河道挟沙能力和水土保持因子等参数上的差别，会影响模型对其他指标的计算，从而产生误差。

本研究表明，由于模型采用每个HRU负荷单独计算，汇总确定流域总负荷的方式，使SWAT模型能够反映每种景观类型因平均坡度改变，引起产流、产沙等指标的变化，但无法反映除坡度和面积变化外，各类型斑块之间因景观空间格局改变对流域非点源污染的影响。因此，在利用SWAT模型分析模拟景观变化时，不应默认模型能够模拟景观空间格局改变对流域水文过程的影响。研究者可以通过划分坡度带，增加子流域HRU的数量，缩小每个土地利用类型斑块的坡度与所属坡度带平均坡度差值，提高模型对不同坡度土地利用的敏感性。

通过进一步研究发现，在保持流域景观类型组成比例不变的情况下，基于模拟土地利用类型与真实土地利用类型的SWAT模型，经过校准都能够很好地反映该流域的产流情况；同时，模型通过对其他参数的修改，使两个模型的产沙数据也能够高度拟合。但二者在泥沙参数上的差异，反映出当模型对一种数据的分析存在不足与误差时，会通过其他参数的调整来达到新的平衡。因此，一个能够反映当地近年产流产沙数据的SWAT模型，不一定是建立在当地各种真实情况的基础上产生的，而是各个参数之间平衡的结果。