SWAT模型子流域划分引起的土地利用变化对径流和输沙模拟结果的影响

钟小敏, 钟科元, 高怡婷, 李煜连, 涂 浙, 黄 莉, 帅 红,2

(1.湖南师范大学 地理科学学院, 长沙 410081; 2.湖南师范大学 地理空间大数据挖掘与应用湖南省重点实验室, 长沙 410081; 3.赣南师范大学 地理与环境工程学院, 江西 赣州 341000; 4.中国科学院大学 西北生态环境资源研究院, 兰州 730000)

在人类活动与气候变化的共同影响下，土地利用/覆被变化及其水文效应成为研究热点[1-2]。合理评价土地利用变化对水文过程的影响对流域水土资源合理利用具有重要意义[3]。分布式水文模型因其较强的物理机理，通过计算单元划分反映流域下垫面及水文要素的空间异质性[4]，能模拟流域不同下垫面特征下的水沙过程，近年来已被广泛应用于土地利用变化的水文效应研究[5-6]。分布式水文模型根据流域内地形、土壤和土地利用等空间分布将流域划分为一定数量的计算单元，计算单元内地形、土壤和土地利用具有一致性，通过这些单元求解模型方程，以此来预测流域响应。然而地形、土壤和土地利用等输入数据的分辨率往往小于这些计算单元的空间分辨率[7]，以致每个计算单元的参数值计算会存在某种程度的集总，模型通过对输入数据进行空间集总来实现流域离散化，而输入数据集总程度对模拟结果具有一定的影响[8]。因此，分析水文模型计算单元划分导致的流域离散化对模拟结果的影响，对减少模拟不确定性、构建高精度的水文模型具有重要意义。

在分布式水文模型发展初期，一些研究已经注意到了计算单元划分引起的流域离散化对模拟结果的影响[9-10]。部分学者对分布式水文模型计算单元划分引起的下垫面空间聚合对流域产流产沙过程的影响开展研究[11-12]，普遍认为子流域划分对径流的影响比较小[13-14]，对泥沙的影响比较大[15-16]。张雪松等[17]发现子流域划分通过改变CN值的空间集总程度，从而影响流域产流量，同时通过改变水文响应单元数量和面积以及通用土壤流失方程中各因子参数值的统计值，对流域产沙量造成影响。Chen等[18]从模型机理角度分析了SWAT模型子流域划分对模拟结果的影响，认为径流量和输沙量受输入参数空间集总与模型结构两个因素的共同影响。钟科元等[19]在桃溪流域进行AnnAGNPS模拟分析发现土壤/土地利用参数空间聚合主要通过改变土壤侵蚀因子和地表覆被与作物栽培管理因子等参数的取值与空间分布，从而导致年平均输沙量减少。

综上所述，子流域划分引起的空间离散化主要分为地形、土壤和土地利用要素的空间离散化，子流域划分方案通过影响这3种下垫面要素，进而影响模拟结果。模型对下垫面属性刻画的准确性是影响模型不确定性的重要因素，已有研究主要关注子流域划分引起的下垫面空间离散化对模拟结果的综合影响，缺乏探讨土地利用、土壤和地形空间离散化分别对水文过程的影响。同时，研究表明东江源是广州和香港地区3 000多万居民的主要饮水源，然而大规模的矿山开采导致植被受到大面积破坏[20]，鉴于此，本文以东江源流域水背水文站以上流域为研究区，构建流域日尺度产流产沙SWAT模型，通过设置11种情景，输入单一地形和土壤数据，以揭示SWAT模型子流域划分引起的土地利用空间离散化对流域产流产沙过程的影响，研究结果为减少模型的不确定性、确定最优计算单元提供科学参考。

1 研究区概况与数据来源

1.1 研究区概况

本文以江西省境内东江源头水背水文站(24°47′56″N，115°40′40″E)以上流域(以下称东江源流域)为研究区(图1)，流域发源于江西省寻乌县桠髻钵山[21]，流域控制面积为904.4 km2。东江源是深圳、广州以及香港等地3 000多万人的主要饮用水来源，该水源的安全与东江流域、珠三角流域经济的繁荣稳定密切相关[20]。研究区地形以山地和丘陵为主，地势北高南低，海拔在179～1 488 m。土地利用类型主要有林地、草地、耕地、建设用地和水域等，其中以林地最为典型，占总面积的82.52%，土壤类型主要有泥红壤、黄壤、石质土、渗育质水稻土等，其中以泥红壤为主，占总面积的80.41%。研究区属亚热带季风气候，降水年内分配不均，主要集中在5—6月份，水资源较为丰富，1979—2017年多年平均气温为20℃，年平均降雨量为1 526.3 mm。该地区人类活动广泛，土地利用变化显著，对其径流和泥沙产生不同程度的影响。

图1 研究区站点分布

1.2 数据来源

空间数据主要包括DEM、土壤、土地利用和气象水文数据。流域DEM数据来源于地理空间数据云(http:∥www.gscloud.cn/)，其空间分辨率为30 m；土壤和土地利用数据来源于江西省水土保持研究院。本文所用的水文气象数据为水背水文站上游12个雨量站点2008—2016年逐日降雨数据(源于赣州市水文局)、水背水文站1980—2016年日径流量和输沙量数据(源于赣州市水文局)、寻乌国家气象站1980—2016年逐日最高气温、最低气温、太阳辐射和风速等数据(源于国家气象局)。

2 研究方法

2.1 SWAT模型

1994年美国农业部农业研究中心开发出建立在地理信息系统基础之上的SWAT模型[22]。模型根据流域地形、河流等进行子流域划分，再通过子流域关键自然要素的输入建立自然条件下的水文模型，并以不同的时间尺度模拟各子流域的径流、输沙等。模型主要采用Green&Ampt入渗方法和径流曲线数(curve number)法计算地表径流量，其中径流曲线数法应用最为广泛[14,23]；采用MUSLE计算地表输沙量[24]。

2.2 模拟方案

SWAT模型基于以下两个方面实现流域离散化：一是基于设定的阈值划分子流域，阈值越小，子流域划分数量越多；二是在每个子流域内进一步划分水文响应单元(HRUs)，HRUs具有相同的地形、土壤和土地利用属性和参数[8]。为进一步分析SWAT模型子流域划分引起的土地利用空间离散化对模拟结果的影响，本文在前人的研究基础上[15,25]，将土壤固定为泥红壤，坡度设置为单一坡度，设置11种情景(表1)以分析土地利用空间离散化对模拟结果的影响。11种方案的子流域集水面积阈值分别为4 000，3 200，2 000，1 500，1 100，800，400，200，100，50，25 hm2，对应的子流域数量分别为11，17，29，39，47，69，137，267，524，1 007，1 955个，同时设定土地利用面积阈值为5%来确定子流域中HRUs的数量。

表1 各土地利用空间离散化方案的特征参数

2.3 土地利用空间离散化对SWAT模型模拟的影响评价

SWAT模型参数变化对径流量和输沙量的模拟结果会产生影响[12]，为避免参数变化对模型输出结果的影响，本文使用模型默认的参数，以实现不同方案下参数的一致性。本文采用11种子流域划分方案，模拟2008—2016年各方案下的产流产沙过程。同时，以年均径流深度和年均输沙量为年指标，年最大1 d、连续最大5 d、连续最大7 d和洪水过程的径流与输沙量为洪水过程指标，采用Manner-Kendall趋势检验法分析径流量和输沙量的变化趋势。采用相对误差(RE)来分析SWAT模型子流域划分引起的土地利用空间离散化对径流和输沙模拟的影响。相对误差RE的计算公式如下：

(1)

式中:Vmin为模拟的最小汇水面积阈值的模拟值;Vi为不同子流域划分方案下的模拟值。

3 结果与分析

3.1 土地利用类型的离散效应

在各划分方案下，流域内各土地利用面积表现出不同变化趋势(图2)。当子流域数量由11个增加至1 955个时，林地面积呈先增加后减少趋势，11种方案下模拟的林地面积均高于实际面积；随着子流域数量的增加，耕地面积呈显著减少趋势(p<0.01)，其面积由128.49 km2减少至123.08 km2，减幅为4.21%；而面积较小的地类(草地、建设用地和水域)，随着子流域数量的增加面积快速增加(p<0.01)，草地、建设用地和水域的增幅分别为16.29%，145.86%，110.55%，并逐渐接近实际地类面积。其中，草地面积模拟值在子流域数量为1 955个时与实际草地面积接近，建设用地和水域在子流域数量分别为29，137个时开始出现，在各方案下的面积均低于实际地类面积。

注：虚线为实际地类面积。图2 子流域划分对土地利用类型的面积变化影响

由于东江源流域地处山区，林地面积大且分布集中，概化程度小；耕地面积较破碎，随着离散程度增加，概化为林地或其他地类的几率逐渐增大，因此耕地面积呈显著减少趋势(p<0.01)；面积较小的地类(草地、建设用地和水域)面积随着离散程度的增加，呈显著增加趋势(p<0.01)，并与实际地类面积越来越接近。

3.2 土地利用空间离散化对径流量的影响

东江源流域不同子流域划分方案下，年均径流深度呈波动下降趋势(图3)，由1 023.61 mm降至1 011.62 mm，在子流域数量为267个后趋于稳定，研究方案中最大RE为1.85%，土地利用空间离散化导致年径流量减少，且随着离散程度的增加，其影响越大。

图3 不同土地利用离散化水平下的年均径流深度

本文选取2016年间的一次暴雨洪水事件(2016年3月21日)，以进一步分析土地利用空间离散化对洪水过程的影响(图4)。从图中可以看出，随着土地利用空间离散化水平的提高，径流过程基本重合，土地利用空间离散化对径流过程的影响不显著。

图4 洪水事件日径流量

SWAT模型采用径流曲线数法模拟地表径流，其中径流曲线数(CN)是影响径流量的主要敏感性参数，降水截留能力越大，CN值越小，不同土地利用类型对降水截留能力存在较大差异，从大到小依次为林地、草地、耕地和建设用地[26]。土地利用空间离散化导致土地利用类型分布在空间和数量上存在差异，改变了CN值的空间集总程度，使其统计特征值发生变化，进而影响流域产流量[17]。通过图5可知，随着土地利用空间离散程度的增加，CN值呈不显著减少趋势，总体变化趋势与年均径流深变化基本一致(图3)，子流域数量在11，17，29个时年均径流深度和CN值呈相似的波动变化。由3.1节分析可知，研究区林地类型占比最大，为总面积的82.52%，草地、建设用地和水域面积小且分布分散，该区域土地利用类型的分布使得各方案下CN值在62.50～63.00变动，总体变化不大(<0.13%)，呈不显著减少趋势，该变化导致其对径流和洪峰影响不大(<1.85%)。这与Lin等[27]研究结果相似。

图5 土地利用空间离散化对CN值的影响

3.3 土地利用空间离散化对输沙量的影响

随着土地利用离散程度的增加，年均输沙量呈显著下降趋势(p<0.01)(图6)，输沙总量由1.5×106t减少至9.04×105t；子流域数量为11，17，29，39，47，69，137，267，524，1 007个的年均输沙量RE分别为65.84%，54.40%，45.05%，37.17%，33.52%，29.17%，19.10%，15.08%，8.99%，5.46%，其中最大RE为65.84%。土地利用空间离散化导致年均输沙量减少，并且随着离散程度的增加该影响越显著。

图6 不同组合方案下的年均输沙量

年最大1 d、连续最大5 d和连续最大7 d输沙量随着土地利用空间离散化水平的提高呈显著下降趋势(p<0.01，图7)，与子流域数量为11时相比，子流域数量为17，29，39，47，69，137，267，524，1 007，1 955个的年最大1 d输沙量分别减少了11.01%，16.59%，23.31%，24.52%，23.49%，28.49%，36.49%，39.13%，41.12%，41.45%(图7)；连续最大5 d的输沙量分别减少了9.04%，15.90%，21.87%，23.85%，24.84%，31.52%，36.40%，40.19%，42.41%，45.30%(图7)；连续最大7 d的输沙量分别减少了7.28%，11.04%，17.47%，19.22%，21.19%，28.38%，33.83%，30.89%，33.30%，40.58%(图7)。随着土地利用空间离散化水平的提高，年最大1 d、连续最大5 d和连续最大7 d输沙量显著下降(41.45%，45.30%，40.58%)，且年最大1 d输沙量下降趋势最显著(p<0.01)。

图7 年最大1 d、连续最大5 d和连续最大7 d输沙量

在2016年3月21日洪水事件中，当子流域数量由11增加到1 955个时，输沙量显著减少(图8)，且洪水事件达到峰值时(2016年3月21日)最显著(图9)，随着土地利用空间离散化水平的提高，输沙量由3.13×105t减少到1.66×105t，减幅达47.16%(p<0.01)，土地利用空间离散化使输沙峰值显著下降(p<0.01)。

图8 2016年3月21日洪水事件日输沙模拟量

图9 土地利用空间离散化对洪水事件日输沙洪峰的影响

SWAT模型通过土壤流失方程对每一个HRUs中由降水和径流产生的土壤侵蚀量和输沙量进行估算，并通过河道汇流到流域出口[22]。USLE_C，USLE_K，USLE_Ls和USLE_P是影响输沙量的主要参数。本文通过输入单一的土壤数据和地形数据从而固定了USLE_K，USLE_Ls和USLE_P参数，因此土地利用空间离散化主要通过USLE_C影响流域输沙过程。由图10可以看出，随着土地利用空间离散化水平的提高，USLE_C呈显著下降趋势(p<0.01)，土地利用空间离散化通过影响USLE_C从而引起输沙量的变化。这一结论与其他学者的研究[10,12]相似。子流域划分引起的土地利用空间离散化，通过改变USLE_C，进而影响输沙量模拟结果。

图10 USLE_C随子流域变化情况

模型将流域划分为一定数量的子流域，径流、输沙模拟结果受到不同子流域划分方案的直接影响，准确确定子流域划分方案的标准尚未统一。部分学者根据流域总面积中集水面积阈值所占比重确定最优子流域划分方案[13,28]，还有部分学者认为可通过模拟结果随子流域数量变化的“拐点”，来确定最优子流域划分方案[19]，如郝芳华[9]、胡连伍[11]、李曼曼[14]等发现可根据“拐点”位置确定最优子流域的划分水平，即当子流域数量不断增加时，输沙模拟值的变化存在一个“拐点”，随着子流域数量的继续增加，模拟结果变化不显著。本文结果发现，在东江源流域，径流量在子流域数量不低于267个时，随着土地利用空间离散化水平的提高，年均径流深度逐渐趋于稳定，而输沙量变化幅度较大，主要是由于东江源流域地处山区，土地利用复杂且破碎，空间差异大，需要更详细的子流域划分方案，以此充分体现该研究区的空间差异，但子流域划分过细容易产生虚假水系，难以准确模拟流域产流产沙过程，子流域数量划分过少则对流域描述不够充分，导致模拟结果不稳定，难以达到理想的预测精度[29]。子流域数量为1 955个时，草地、建设用地和水域面积模拟值已逐渐接近实际面积值，但是并未发现输沙量“拐点”的存在。由于受计算机运行能力的限制，子流域继续增加则模型无法运行，因此无法确定研究区最优集水面积阈值。最优计算单元划分受流域面积、下垫面复杂程度等因素的综合影响，因此，为减少分布式水文模拟的不确定性及确定最优子流域划分方案，不同研究区需要具体分析。

4 结 论

(1) 土地利用离散化使得流域内各类土地利用的面积表现出不同变化趋势。随着土地利用离散程度的增加，林地面积呈先增加后减少趋势，耕地面积呈显著下降趋势(p<0.01)，草地、建设用地和水域面积均呈显著增加趋势(p<0.01)，并随离散程度增加越接近实际面积。

(2) 年均径流量和输沙量模拟值对土地利用空间离散化的敏感性较强(p<0.01)，且输沙量对土地利用空间离散化的敏感程度大于径流量，径流量模拟值最大RE为1.85%，输沙量模拟值最大RE为65.84%。

(3) 土地利用空间离散化对年最大1 d、连续最大5 d和连续最大7 d输沙量的影响随着离散程度的增加呈显著下降趋势(p<0.01)，其中对年最大1 d的影响最为显著。土地利用空间离散化导致泥沙峰值显著下降，而对洪峰流量影响不显著。

(4) 本研究结果表明，东江源流域子流域数目为1 955个时，未发现“拐点”的存在，受计算机运行能力的影响，无法确定研究区最优集水面积阈值，因此不同研究区需要具体分析，以减少分布式水文模拟的不确定性，确定最优子流域划分方案。