水文响应单元划分对SWAT模型总氮模拟效果的影响

郑思远，王飞儿，2*，俞 洁，董秋仪

（1.浙江大学环境与资源学院，杭州 310058；2.浙江省水体污染控制与环境安全技术重点实验室，杭州 310058；3.浙江省环境监测中心，杭州 310012）

分布式模型具有一定的物理基础，充分考虑空间异质性，是目前非点源污染量化研究最直接有效的工具之一，为非点源污染控制和管理的定量化研究提供了有效的技术手段[1]。在众多的模型中，SWAT（Soil and Water Assessment Tool）模型因其在水文和非点源污染模拟上有较好的适用性，而得到广泛的应用[2-4]。SWAT模型主要通过将相同坡度、土地利用方式和土壤特征的地块划分成水文响应单元（Hydrological response units，HRU），独立计算出每个单元的水量、泥沙以及污染物负荷，进而得到子流域的输出总量[5]。SWAT模型的模拟受到如模型结构、输入数据和模型参数等诸多因素的影响[6-7]，模型的模拟结果也存在一定的不确定性。其中空间离散化是实现SWAT模型模拟的首要步骤和关键步骤，并且会影响后续径流、泥沙和污染物模拟的准确性。

SWAT模型的空间离散化主要是根据数字高程模型（Digital elevation model，DEM）进行子流域的划分，在子流域划分的基础上根据土地利用-土壤类型-坡度阈值组合生成更小的HRU，表示单一地面覆盖、单一土壤类型的研究单元。关于空间离散化的不确定性研究目前较多关注于子流域层面，包括DEM和集水面积阈值的选用[8-10]，研究发现过多或过少的子流域划分均会导致模拟的偏差，对产沙模拟的影响较大，而不同地区合适的子流域划分方案不一，均需符合实际情况[11-13]，但目前针对HRU划分层面的研究相对较少。而HRU作为SWAT模型的基本计算单位，对于模型的模拟结果有直接的影响。王琳等[14]对子流域和HRU划分的综合研究也发现，SWAT模型模拟结果对HRU的敏感性高于子流域划分，因此有必要单独针对HRU划分层面开展研究。已有的研究结果表明，土地利用、土壤类型和坡度阈值会影响HRU的划分，不同的HRU划分方案对SWAT模型的模拟结果产生影响，较小的土地利用和土壤类型的阈值组合能得到较好的模拟效果，而针对坡度阈值的研究相对较少[15-17]。但目前单独针对HRU划分层面，分别考虑土地利用、土壤类型和坡度阈值对总氮模拟效果影响规律的研究较少，且在苕溪流域开展的SWAT模型研究采用的土地利用-土壤类型阈值组合不一，有5%～10%、10%～5%和10%～15%等划分组合[18-20]，对坡度阈值的选择也不明确，因此不同方案的合理性和适用性仍有待实证。本文在相对合理的子流域划分条件下，针对HRU划分层面开展探讨，拟通过设置18种不同的HRU划分阈值范围组合，以东苕溪上游区域为实证，以模型可靠性为约束准则，较为全面地探讨坡度、土地利用和土壤类型阈值分别对总氮模拟效果影响的规律，并综合确定该流域总氮模拟相对合适的HRU划分方案。

1 材料与方法

1.1 研究区域概况

苕溪位于浙江省北部，属太湖水系，是太湖流域重要的入湖河道，有东、西苕溪两大支流，主流长度157.4 km，流域总面积4 576.4 km2。本文选取东苕溪瓶窑上游段流域为研究对象，包含苕溪上游南苕溪、中苕溪和北苕溪3段，流域面积1393 km2（119°31′～119°37′E、30°21′～30°23′N）（图 1）。该流域属于典型的亚热带季风气候，年平均气温15.5～15.8℃，多年平均降雨量为1460 mm。流域的地带性土壤为黄红壤，主要农作物为水稻。研究区内所占面积最大的土地利用类型是林地，达70%；其次是水稻田，占16%；城镇建设用地占总面积的8%；除此之外的水域、果园所占面积相对较小，两者共占6%左右。流域内非点源污染物主要来源于农业、畜禽养殖业、渔业和农村生活等。

1.2 SWAT建模数据

SWAT模型的输入数据包括空间数据和属性数据，空间数据包括数字高程（DEM）图（30 m分辨率）、土壤类型图（1∶50 000）、土地利用图（1∶10 000）等，均来自浙江大学环境与资源学院遥感所；属性数据包括2011—2014年的气象数据（桥东村、临安、横畈、余杭、莫干山等气象站点的逐日降雨量、日最高最低气温、日平均风速和相对湿度）、土壤数据（土壤理化性质）、社会经济数据（各县市区人口、农业种植施肥量、畜禽养殖量、水产养殖产量等）、点源污染数据（污水处理厂和工业点源的氮磷排放量）以及水文观测数据（瓶窑断面日均径流量和月均水质数据）等，其中土壤类型及其理化性质参考第二次全国土壤普查及《浙江土壤志》，另有水稻田管理数据根据实地调查得到。

1.3 划分方案设计

根据DEM图对研究区域进行子流域划分，并结合实际的水系进行修正，使子流域的划分相对符合实际情况，共得到26个子流域，具体见图1。在子流域划分的基础上，根据坡度、土地利用与土壤类型3个基本要素的阈值（分别记为A%、B%、C%）组合进行HRU的划分，具体步骤见图2[15]。当子流域中要素某一类型的面积占比小于设定的阈值，这个类型将被合并到超过该阈值的其他类型当中，从而将同一子流域内具有相同类型要素的HRU合并在一起[21]。

图1 研究区位置图Figure 1 Location of study area

为了探讨合理的HRU划分阈值组合模式，本研究分别设置了A、B、C 3组方案，每组方案分别对坡度、土地利用和土壤类型阈值变化的影响进行了探究（表1）。本研究设置的阈值范围为5%～30%，A组在土地利用和土壤类型阈值均为5%的条件下，以5%为梯度，对坡度分级要素设置了从5%～30%的变化范围，通过比较总氮模拟值与观测值的符合度探究坡度阈值对总氮模拟的影响并确定坡度分级的合理阈值；B组设定土壤类型阈值为5%，在A方案确定的坡度阈值条件下，同样以5%为梯度，对土地利用阈值从5%～30%变化，探究合理的土地利用阈值；C组在A、B方案确定的坡度和土地利用阈值条件下，探究土壤类型阈值变化对总氮模拟的影响。由此得到18种HRU的划分阈值组合情景，通过不同方案下总氮的模拟结果，分别探讨坡度、土地利用和土壤类型阈值对总氮模拟的影响，并综合确定本研究相对合理的HRU划分方案。

图2 SWAT模型HRU划分流程图Figure 2 Process of HRU delineation by SWAT model

1.4 方案评价方法

为了更好地反映研究区流域特征，以坡度、土地利用和土壤类型阈值均为0的方案作为基准情景进行模型参数率定。2012—2013年为率定期，2014年为验证期进行径流和总氮率定。采用SWAT-CUP的SUFI2率定方法，主要以决定系数（R2，式1）和纳什系数（NSE，式2）衡量，R2表示实测数据和模拟数据趋势的一致性，NSE用于衡量实测数据和模拟数据的拟合程度。

为控制变量，所有HRU划分方案的SWAT模型构建均采用相同的数据输入和一致的子流域划分，模型参数均采用基准情景率定得到的参数值，以此运行模型，得到不同HRU划分方案的总氮模拟值。采用纳什系数和均方根差与标准差的比值（RSR，式3）对各方案下模型总氮输出结果与观测值进行评估，进而探讨不同HRU划分方案对SWAT模型总氮模拟效果的影响。

表1 HRU划分阈值组合方案设计Table 1 Different threshold schemes for HRU delineation

式中:Qobs，i为观测值；Qsim，i为模拟值；Qobs为所有观测值的平均值；Qsim为所有模拟值的平均值；n为观测数据个数。NSE与R2越接近1，RSR越接近0，模拟的效果越好。一般来说，R2大于0.6，NSE大于0.5，RSR小于0.7，模型模拟效果相对较好[22]。

2 结果与分析

2.1 SWAT模型率定验证

在坡度、土地利用和土壤阈值均为0的基准情景下进行SWAT模型率定:径流率定期R2为0.81，NSE为0.52，RSR为0.69，验证期R2为0.89，NSE为0.53，RSR 为 0.68（图 3）；总氮率定期R2为 0.60，NSE 为0.60，RSR为0.64，验证期R2为0.65，NSE为0.54，RSR为0.68（图4）；模拟值与观测值一致性较好，但峰值模拟存在一定的偏差，但总体来说，满足模型模拟要求。具体参数取值见表2。

2.2 不同的HRU划分对总氮模拟的影响

18种HRU划分方案（除土地利用阈值为30%的B-6方案）总氮模拟的NSE值大于0.5，RSR值小于0.7，符合模型模拟可信度要求（图5）。A方案中，不同坡度情景下，NSE和RSR分别在0.51～0.57和0.65～0.70范围内变化，变化趋势较为明显（图5a）。随着坡度分级阈值的增加，NSE升高，RSR下降，表明较为粗放的坡度分级有利于缩小总氮模拟结果与实际的差值。在坡度分级阈值达到20%时，NSE最大，RSR最低，随后两者略有下降或上升，但变化相对趋于稳定。因此，A方案中较为粗放的坡度阈值（20%～30%）对于总氮模拟较为合适，且以20%相对最佳，故将20%作为B、C方案确定的坡度阈值。

图3 基准情景下径流率定结果Figure 3 Calibration result of runoff in basic scheme

图4 基准情景下总氮率定结果Figure 4 Calibration result of TN in basic scheme

表2 SWAT模型参数取值Table 2 Parameters selected for SWAT model

总氮对土地利用阈值变化相对敏感，B方案中NSE和RSR的变化幅度相对较大，分别在0.49～0.57和0.65～0.71范围内（图5b）。随着土地利用类型阈值的增大，NSE减小，RSR增大，总氮的模拟效果大幅度降低，表明在较小的阈值下，总氮模拟的效果较好，土地利用阈值增大会导致总氮模拟结果偏离实测值。当划分阈值上升到30%时，NSE骤降，出现了NSE＜0.5、RSR＞0.7的情况，模型模拟存在一定的失真现象。土地利用阈值在5%～10%时，NSE降幅不大，由此，B方案中较小的土地利用阈值范围（5%～10%）对总氮模拟相对合适，且以5%相对最佳，并以此作为C方案确定的土地利用类型阈值。

C方案中，NSE和RSR分别在0.55～0.57和0.65～0.67范围内变化，土壤阈值的影响相对较小（图5c）。随着划分阈值的增加，整体NSE呈现下降趋势，RSR呈现上升趋势，但在土壤阈值为15%～20%时略有波动。从模拟结果来看，较小的土壤类型阈值有利于提高总氮模拟效果。C方案中，较小的土壤类型阈值范围（5%～10%）对总氮模拟相对合适，且以5%相对最佳。但在实际的应用过程中，由于土壤阈值影响相对较小，考虑到模型的运行效率，土壤阈值可以适当放大。

图5 不同HRU划分方案对总氮模拟的影响Figure 5 Effect of HRU delineation on TN simulation

综合A、B、C 3组方案的总氮模拟情况，较小的土地利用阈值和土壤类型阈值，较为粗放的坡度阈值有利于减小总氮模拟值和观测值的偏差，提高总氮模拟效果，在18种方案中，坡度阈值20%，土地利用阈值和土壤阈值均为5%的方案相对最为合适。

3 讨论

SWAT模型的HRU划分是通过设定土地利用阈值，清除各子流域内的次要土地利用类，设定土壤类型阈值来清除特定土地利用面积中的次要土壤类，设定坡度阈值来清除特定土地利用和土壤类型面积中次要的坡类，进而得到单一的土地利用、土壤类型和坡度组合，各阈值的影响原因相似，均通过确定相对合适的划分阈值，使HRU的划分更符合实际情况，进而得到较好的模型模拟结果。从A、B、C 3组18种HRU划分方案来看，土地利用阈值对总氮模拟的影响最大，且在B方案中当土地利用阈值为30%时存在一定的模拟失真现象，因此，以土地利用阈值为例，对其影响的原因进一步分析和讨论。

在HRU划分中，通过设定土地利用阈值，当某类土地利用类型的面积小于设定阈值，这类土地利用类型将被合并到超过该阈值的其他土地利用类型当中。从表3中发现，随着土地利用阈值的增加，耕地、果园、城镇和水域部分面积会合并到林地这一土地利用类型中，导致耕地、果园、城镇和水域的总面积均下降，而林地的面积占比不断升高。通过相关性分析发现，总氮模拟NSE同耕地面积占比呈显著正相关，同林地面积占比呈极显著负相关（表4）。当土地利用类型划分阈值达到25%时，果园这一土地利用类型缺失，城镇和水域的面积占比低至近0，耕地面积占比也有所下降，而林地的面积占比不断升高，由于林地对总氮负荷的贡献相对较小，因此较大的土地利用类型划分阈值会导致模拟的总氮负荷偏小，整体模拟效果变差。由于在农作物的播种、生长、收成等不同时期需要进行施肥等管理措施，总氮模拟NSE与耕地这类土地利用类型的面积占比呈显著正相关，耕地面积占比在土地利用阈值为30%时骤降，降至3.11%，存在较大的地理信息缺失，也导致了该阈值条件下模型总氮模拟NSE降到0.5以下，存在一定的模拟失真现象。

表3 不同土地利用类型面积（占比）随划分阈值的变化Table 3 Each land use area（ratio）under different threshold schemes

表4 不同土地利用面积占比和总氮模拟NSE相关性Table 4 Correlation analysis between different landuse area ratio and NSE of TN simulation

4 结论

（1）在18种HRU划分方案中，坡度阈值20%、土地利用阈值和土壤类型阈值均为5%的方案对于本研究总氮模拟来说相对最佳。

（2）土地利用阈值对总氮模拟的影响相对最大，其通过对土地利用类型的筛选，使总氮产量高低不同的HRU面积发生变化来影响总氮负荷的模拟量，其中耕地由于施肥管理等措施，其面积占比是主要的影响因素。

（3）土壤类型阈值的影响相对较小，在实际的模型应用过程中，兼顾模型的模拟精度和运行效率，在合理范围内尽量选择小的土地利用阈值，坡度和土壤类型阈值可以适当放大。