子流域划分在时间尺度上对年径流、泥沙及营养物的影响——以漳河上游流域为例

武海霞，李清雪，孙玉壮，冀雅珍

(1.河北省水资源高效利用工程技术研究中心，河北工程大学，河北 邯郸 056038；2.河北省煤炭资源综合开发与利用协同创新中心, 河北省资源勘探研究重点实验室, 河北工程大学，河北 邯郸 056038；3.山西省水利水电科学研究院，太原 030002)

0 引 言

目前，SWAT分布式水文模型在流域水文、生态和环境等领域得到广泛应用。然而，模型变量和参数具有很大的随机性和时空变异性，且这些变量和参数的变化会造成模拟结果较大的差异[1]。流域离散是分布式水文模型的基础，子流域的数量和面积对模拟结果影响差异极为显著，子流域离散方案产生的模拟偏差，可能会超过模型参数率定所产生的偏差，而且在地形复杂、坡度较大的流域，不同子流域数量使模拟精度存在不确定性[2]。关于子流域划分对流域模拟结果的影响，国内外学者进行了相关研究。一些学者认为子流域划分对流域产流模拟影响较小,而对产沙模拟影响较大[3,4]。郝芳华等和张雪松等的研究表明，当子流域划分数量达到一定水平时，增加子流域划分对模型输出结果的影响较小[5,6]。部分学者通过进一步分析，认为在不同集水面积阈值条件下其模拟结果存在阈值适用范围的上下限。Tripathi M和Jha M认为对流量、泥沙、NO3-N 和无机磷模拟的最适子流域面积占流域总面积的比例分别为4%、3%、2%和5%[7,8]。同时，基于SWAT模型的子流域划分方案对模拟结果的影响和研究尺度有一定的关系。Soni等的研究表明SWAT模型的产流预测对子流域和水文响应单元的数量具有一定响应，当研究的尺度太大时，流域参数的集总会影响径流曲线CN(Curve Number)值的分布，从而导致预测的不准确[9]。虽然国内外学者针对分布式水文模型子流域划分方案进行了相关研究，但由于子流域的数量和面积对模拟结果影响差异极为显著，模拟结果的不确定因素比较复杂，对于如何准确而有效地确定子流域划分数量并没有统一的结论和标准，而且在不同流域，结论并不完全一致。本次以大尺度复杂山区漳河上游流域为研究对象，分析了在不同水平年流域划分方案对SWAT模型模拟结果的影响，对准确模拟该地区的水沙情况奠定一定基础，可为该地区水资源水环境的进一步研究提供技术支撑。

1 流域概况及模型模拟

1.1 流域概况

本次研究选择漳河上游流域作为研究对象,漳河上游位于海河流域西南部,具体位置见图1。流域面积约17 730 km2，占海河流域总面积的6%。流域属温带大陆性季风型气候，四季分明。全年夏短冬长，多年平均气温7.4～10.3 ℃，无霜期170 d左右。流域多年平均降水量为576.9 mm，降水量总的分布趋势是南部大于北部，迎风坡多于背风坡和盆地中心。流域内支流众多，水系呈扇形分布，分东西两大支流。东支为清漳，该地区为石质山区，山高谷深，岩石裸露，坡陡流急，含沙量小；西支为浊漳，该区域为山丘盆地区，盆地内黄土覆盖较厚，植被较差，水土流失严重，洪水挟带泥沙较多。

漳河上游流域内以山区地貌为主，地形落差高达2 000多米，东邻太行山脉，西邻太岳山脉，中间是长治盆地。河流横穿太行山脉，大地构造处于太行山新华夏系第三隆起带，由一系列北东向平缓复式群皱组成，岩性地层出露主要有震旦系石英砂岩、泥岩，寒武系紫红色页岩，奥陶系石灰岩，河谷两岩阶地由第四系积物组成。

图1 漳河上游流域位置

1.2 模型数据和参数

模型所需的空间数据主要包括数字高程模型(DEM)、土壤数据、土地利用类型数据和气象数据。DEM使用国际科学数据服务网站提供的分辨率为30 m的数据，依照研究区经纬度范围下载，将多景的数据进行融合与拼接，使用数字化的研究区边界裁切出边界精准的漳河上游流域的DEM图。土壤数据来源于世界粮农组织网站的HWSD(Harmonized World Soil Database)数据。数据格式为栅格数据，分辨率为1 000 m，投影为WGS1984，土壤分类系统是FAO90。由于模型采用的是美国土壤分类系统，需要重新建立研究区的土壤属性数据库。土地利用类型数据由海河流域南部 1∶100 万土地利用图片数字化得到，并根据土地利用的一级分类，转化为 SWAT 模型能识别的模型代码，最终划分为 5种土地利用类型：耕地、林地、草地、水域和其他用地。河网水系图来源于国家基础地理中心提供的分辨率为30 m的数据。降雨、温度、相对湿度、风速及太阳辐射等气象数据来源于SWAT 官方网站及中国气象科学数据共享服务网。模型所需基础数据的详细描述见表1。利用Arcgis软件对DEM、土壤类型图和土地利用图作Albers等积圆锥投影统一地理坐标系WGS1984，并统一生成相同的栅格大小。

表1 模型所需基础数据

根据郝芳华[5]等的研究，由于现有的确定模型参数的最优化方法对实测资料的依赖性很大，只能反映模拟值与实测值的拟合程度，而不能揭示参数的物理意义，因此本次直接采用模型默认的参数值进行模拟。

1.3 子流域划分方案

本次研究基于最小集水面积( Critical Source Area，CSA)阈值划分出不同数量的子流域。最小面积阈值(CSA)是形成河流的最小集水面积，也可以表示为其与流域总面积的百分比，CSA越大子流域划分数量越少。

为深入分析流域离散化程度对SWAT模型模拟结果的影响，本次研究在前人成果的基础上，结合研究区特点选择6种子流域划分方法，即CSA分别为0.4%、0.5%、0.7%、1%、2%和10%，并通过确定土地利用和土壤类型面积阈值分别为 5% 和10%来确定子流域中 HRU 的数量。表 2 列出了6种流域离散化方案的特征参数。

表2 子流域划分及流域参数

根据CSA值对漳河上游流域进行子流域划分,获得6种不同的子流域划分结果见图2。

图2 流域不同子流域划分层次图

定义相对偏差为：

Re=(Vi-Vmin)/Vmin×100%

式中：Vi为不同子流域划分模拟得到的模拟值；Vmin为子流域数目最小时得到的模拟结果。

本次研究选择的实测资料来自于课题组其他成员的分析成果。王喆[10]利用漳河上游近40年的气象资料进行面降雨量的趋势性分析，得出1979-1995年为枯水年，1996年为丰水年，1997-2013年为平水年。本次研究在时间尺度上选择1995、1996年和1997年的气象资料作为枯水年，丰水年和平水年的输入数据进行模拟结果的分析。

子流域划分方案下年均径流、泥沙及营养物在平水年的相对偏差(Re)见表3。

表3 不同子流域划分下年径流、泥沙及营养物在平水年的相对偏差(Re)

注：Re1、Re2、Re3、Re4分别为年均径流、泥沙、有机氮和有机磷相对偏差值。

2 结果与讨论

2.1 子流域划分方案对地表径流的影响

由图2可看出，不同子流域划分结果最直接的表现是各子流域内面积大小的变化及河网详实程度的描述。由表2可看出，子流域划分方案对径流模拟结果有一定的影响，最大偏差值为11.23%，但和泥沙负荷相比而言，子流域划分方案对径流的影响并不大，尤其在子流域个数为58～111范围内，影响很小，这与前人的研究结论一致。这主要是因为本文所选用的产流方程为径流曲线数法，其中决定产流量的敏感性参数为径流曲线数CN，而不同的划分方法对流域面积加权CN值的影响比较小。

图3 不同水平年径流量随子流域数量变化曲线

图3显示了在不同水平年径流量随子流域数量变化的趋势，总体表现为先增大，然后趋于稳定，最后呈减小趋势。具体表现为当子流域数量从5增加到25时，年径流量随之增大，丰水年增加的幅度最大，而枯水年呈略微减小趋势，表明这两者之间的子流域划分并不是最佳的子流域划分方案。当子流域划分数量进一步增加，子流域个数大于25，小于111，CSA值为流域总面积的0.5%～2%时，年径流量趋于稳定，且不同水平年趋于稳定的范围一致。当子流域划分数量超过111个时，径流量明显减少，尤其在丰水年减小更明显，而枯水年呈略微的增加趋势。这说明在不同水平年径流量随子流域划分数量变化的趋势并不完全一致。造成这种情况的原因可能是过多的子流域划分，形成了较多的过于狭长的子流域，使得模型出现了对水文过程的不真实模拟[11]。可见，在漳河上游流域，合理的子流域划分水平对径流模拟的影响并不大，但存在一个使径流模拟趋于稳定的子流域划分的上下限，当低于下限，高于上限时，模拟结果受划分数量的影响较大。

2.2 子流域划分方案对泥沙的影响

由表3可看出，泥沙负荷对子流域划分方案的敏感性强于径流量，最大偏差值高达30.31%。由图4可看出，在不同水平年泥沙负荷随子流域划分数量的增加逐步呈增加趋势，这和年均径流量的变化趋势并不完全一致。当子流域数量由5变到25时,泥沙负荷有增加趋势，当子流域个数在25～111范围时，泥沙负荷趋于稳定，这和年径流量趋于稳定的范围一致。但泥沙和年均径流在不同水平年随子流域数量变化的趋势并不完全相同，当子流域个数超过模拟结果趋于稳定的上限值111个时，年均径流量明显减少，而泥沙量明显增加，尤其在丰水年表现的更明显。这可能是因为子流域个数改变引起流域内高程的空间变化，从而影响到流域的坡面坡长和坡度，而坡度和坡长参数(LS 因子)对于泥沙量是敏感因子，能显著影响泥沙产量的预测。泥沙演算中的泥沙沉积和降解过程会影响产沙预测的结果，随着子流域划分数目的减小，流域描述简化，河网密度减小，从而影响泥沙演算，并使预测的准确性降低[12]。

图4 不同水平年泥沙随子流域数量变化曲线

2.3 子流域划分方案对营养物的影响

营养物氮、磷的变化原因比较复杂，除了受径流、沉积负荷的影响外，还受到自身复杂计算算法的影响，而有机态养分的迁移主要是通过吸附在泥沙等固体上进行的。

随着子流域划分数目的变化，营养物有机氮和有机磷负荷发生明显变化。由表3可看出，最粗略的子流域划分和最精细划分之间的相对误差，最大偏差值分别为26.53%和45.30%，而且有机磷的相对偏差值总体都比有机氮的大。 由此可见，子流域划分对营养物负荷的影响较大，这和前人的研究结果一致。

由图5、6可看出，在不同水平年，随着子流域划分数量的增加，流域出口营养物有机氮、磷总体上呈现增加的趋势，和泥沙负荷随子流域划分数量的变化趋势类似，营养物负荷随子流域数目变化趋于稳定的范围和泥沙负荷的范围一致，都在25～111。

图5 不同水平年有机氮随子流域数量变化曲线

图6 不同水平年有机磷随子流域数量变化曲线

有机氮、磷负荷与泥沙负荷是成比例的，但是由于有机氮磷的河道传输与泥沙的河道传输算法的不同，所以出口有机氮磷与泥沙并不一定表现出共同的趋势。径流量和泥沙负荷都随子流域划分方案的不同而产生一定变化，且各自具有不同的变化趋势和变化幅度(图 3、图 4)，因此营养物负荷量的变化原因很难归结为某一简单的原因。

2.4 最佳子流域划分方案的确定

对于子流域划分方案对SWAT模拟结果的影响，多数学者认为对径流量的影响较小，对泥沙和营养物的影响较大，但是如何准确有效地确定子流域划分数目至今并没有统一的标准。一般来讲，地形复杂的山区，由于地形破碎，空间差异大，需要更详细的子流域划分，以此充分体现该地区局部空间的差异。但CSA值过小，划分过细时会产生较多过于狭长或过小的虚假子流域，导致不真实的水文模拟，而当亚流域划分数量较少时，对流域描述不够充分导致模拟输出结果不稳定，难以达到理想的预测精度。本次研究在不同子流域划分方案下的模拟结果分析表明，在漳河上游流域，无论是径流、泥沙、还是营养物有机氮磷，对子流域划分方案趋于稳定的范围一致，都在25～111之间。当子流域的数量从5个增加到25个时，年均径流、泥沙和有机氮磷相对来说变化幅度较大，且在不同水平年表现出的变化趋势并不一致，这表明此区间并不是最佳的子流域划分方案。随着子流域划分数量的进一步增加，年均径流、泥沙和有机氮磷逐渐趋于稳定，当超过某一子流域划分方案时，划分的HRU数量急剧增加，计算单元划分过于精细，可能导致计算结果失真，同时花费更长的计算时间，降低计算效率。因此，在漳河上游流域，存在上下限阈值111、25，当子流域个数低于下限阈值，对空间的描述不够，难以满足流域规划与管理的要求，高于上限阈值，模拟结果出现失真。通过综合权衡下垫面特征、与其他数据的协调及模型运行能力等因素的基础上，根据模型模拟结果，最后确定漳河上游流域最佳子流域划分个数为79个，CSA取值为12 411 hm2，仅占流域总面积的0.7%，这说明漳河上游流域内复杂的地理空间环境需要更加详细的描述，这与实际需求相符。本次研究结果中子流域平均划分面积占整个流域面积的比值与以往的研究结论并不一致，说明利用子流域面积阈值设定不同的子流域划分方案存在地域差异，需要针对不同流域的实际地形、地貌、海拔等空间属性对模型进行初步的径流、泥沙和营养物质的敏感性分析，合理确定CSA的取值及子流域的划分方案，从而减少分布式水文模型模拟的不确定性，为进一步更准确地模拟水文过程奠定基础。

3 结 语

本文以中国海河流域西南部的漳河上游流域为研究对象，建立流域SWAT年尺度产流产沙模型，进一步分析不同CSA取值下，子流域划分数量对径流、泥沙及营养物模拟的影响。主要结论如下：

(1)从总体来看, 子流域划分层次对模拟结果的影响存在上限下限两个阈值25、111，介于两个阈值之间模拟结果相对稳定，超出上限阈值111，模拟结果出现失真；低于下限阈值25，对空间的描述不够，难以满足流域规划与管理的要求。子流域划分数量对年均径流量影响较小，对泥沙和营养物影响较大。

(2)在不同水平年，模拟结果对子流域划分数量的敏感性不一致，在丰水年最敏感，在枯水年最不敏感，建议在径流量较大的地区，要合理选择子流域划分个数。

(3)综合考虑模拟精度和模拟效率，漳河上游流域 SWAT 模拟的最合理子流域划分水平为79个，CSA值为12 411 hm2，占整个流域面积的0.7%。

□

[1] Manillapalli S. Effect of spatial variability on basin scale modeling[C]∥Third International Conference Workshop on Integrating GIS and Environmental Modeling. Proceedings, Santa Fe, New Mexico, 1996.

[2] Kim J G, Park Y, Yoo D, et al. Development of a SWAT patch for better estimation of sediment yield in steep sloping watersheds [J]. Journal of the American Water Resources Association，2009，45(4):963-972.

[3] FitzHugh TW,Mackay D S. Impacts of input parameter spatial aggregation on an agricultural nonpoint source pollution model[J]. Journal of Hydrology, 2000,(236):35- 53.

[4] 卢文喜，伊燕平，张 蕾，等． 不同亚流域划分数量对SWAT 模型模拟结果的影响[J]．水电能源科学，2010，28(10):23-25.

[5] 郝芳华，张雪松，程红光，等．分布式水文模型亚流域合理划分水平刍议[J]．水土保持学报，2003，12(4)：75-77．

[6] 张雪松,郝芳华,程红光,等.亚流域划分对分布式水文模型模拟结果的影响 [J]. 水利学报,2004,(7):119-120.

[7] Tripathi M, Raghuwanshi N, Rao G. Effect of watershed subdivision on simulation of water balance components[J].Hydrological Processes,2005,20(5):1 137-1 151.

[8] Jha M, Gassman P W, Secchi S, et al. Effect of watershed subdivision on swat flow, sediment, and nutrient predictions[J].Journal of the American Water Resources Association,2007,40( 3):811-825.

[9] Soni M. Pradhanang, Russell D. B. Effects of critical source area on sediment yield and streamflow [J].Water and Environment Journal, 2014,28:222-232.

[10] 王 喆．气候变化和人类活动对漳河上游地表水资源的影响研究[D]．河北邯郸：河北工程大学城市建设学院，2015．

[11] 耿润哲，李明涛，王晓燕，等. 基于SWAT模型的流域土地利用格局变化对面源污染的影响[J]．农业工程学报，2015,31(16):241-250.

[12] 张召喜，罗春燕，张敬锁，等．子流域划分对农业面源污染模拟结果的影响[J]．农业环境科学学报，2012，31(10)：1 986-1 993．