清华大学 水沙科学与水利水电工程国家重点实验室, 北京 100084
收稿日期：2019-03-18
基金项目：国家重点研发计划项目（2016YFC0402406）
作者简介：李二辉(1987-), 男, 博士后
通信作者：傅旭东, 教授, E-mail:xdfu@tsinghua.edu.cn
摘要：针对SWAT（soil and water assessment tool）模型逐日产沙模拟中含沙量-流速关系的尺度匹配问题，给出了日尺度产沙模拟的最大含沙量封闭关系。该文采用黄河中游砒砂岩区9个水文站长系列实测水沙资料，确定日输沙率-日流量的幂函数型关系曲线，提供了SWAT模型日尺度最大含沙量的封闭关系。并将封闭的SWAT模型用于皇甫川流域水沙产输过程模拟，采用空间嵌套的皇甫站和沙圪堵站实测资料进行模型验证。结果表明：皇甫站率定期和验证期日输沙过程的Nash-Sutcliffe效率系数（Nash-Sutcliffe efficiency coefficient，NSE）分别为0.70和0.66，决定系数R2分别为0.74和0.68；沙圪堵站率定期和验证期NSE分别为0.78和0.72，R2分别为0.78和0.74；未进行封闭模型的沙圪堵站率定期和验证期NSE分别为0.36和0.26，R2分别为0.39和0.26。表明封闭后的SWAT模型显著提高了模拟效率，能够用于黄土高原砒砂岩地区泥沙逐日过程模拟。
关键词：SWAT模型输沙率-流量关系曲线封闭关系日尺度皇甫川流域
Closure and validation of a SWAT model for daily scale sediment simulations
LI Erhui, WANG Bingjie, FU Xudong
State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing 100084, China
Abstract: A closed relationship for daily maximum sediment concentration was developed to solve the scale matching problem between the sediment concentration and the flow velocity in the SWAT model. A power function was used to fit the daily flow rate to the sedimentation flux at 9 gauging stations located in the soft sandstone region of the middle Yellow River to provide a closed relationship for the daily sedimentation rate. The closed SWAT model was then used to simulate the sediment transport in the Huangfuchuan watershed of the Yellow River. The model was validated by comparing the observed and simulated sediment loads at the Huangfu and Shagedu stations. The results showed that the Nash-Sutcliffe efficiency coefficient (NSE) was 0.70 and the determination coefficient (R2) was 0.74 during the calibration period and were equal to 0.66 and 0.68 during the validation period for the daily scale sediment load at the Huangfu station. The closed SWAT model predictions at the Shagedu station correlated with the observed sediment load with NSE=0.78 and R2=0.78 for the calibration period and NSE=0.72 and R2=0.74 for the validation period. The unclosed SWAT model had NSE=0.36 and R2=0.39 for the calibration period and NSE=0.26 and R2=0.26 for the validation period at the Shagedu station. The results indicate that the closed SWAT model improves the simulation accuracy for sediment transport within the basin to improve applications of the SWAT model in the soft sandstone region for sediment simulation.
Key words: SWAT modelsediment rating curveclosed relationshipdaily scaleHuangfuchuan basin
泥沙的侵蚀、搬运和沉积会对流域地貌发育过程产生影响。开展流域产沙过程模拟，有助于认识泥沙历史变化过程规律和未来演变趋势，对流域治理和河流管理具有重要意义。
黄河是中国第二大河，以含沙量高而著称。黄河流经的黄土高原地区土质疏松，地形破碎，植被覆盖度低，暴雨集中，使得流域洪水涨洪快，洪峰流量高，历时短，并携带大量泥沙。据冉大川等[1]研究，黄河中游泾河流域1966年7月24日一次暴雨平均次降水量为150.8 mm，占全年降水量的26.7%，而次洪量占年径流量的46.1%，产沙量占年输沙量的70.3%。由于黄土高原场次降雨、洪水泥沙持续时间短暂，开展流域产沙模拟应尽可能采用小的时间步长进行研究。而现在的流域治理规划不仅要考虑一次降水过程的水文响应，并且还要揭示较长系列的水文周期及水沙情势的演变趋势，其涉及的区域往往较大，需要考虑产流产沙模型的时空尺度。同时，较小的时间尺度虽能将流域的降水、产流产沙过程描述清楚，但数据量往往过大，计算效率低，成本较高。现在以日尺度为时间步长的水文实测资料较容易获取，并且能够反映长时间的水沙演变趋势。因此，开展日尺度泥沙输移过程模拟能够兼顾实测资料限制和长时间过程的产沙预测的需求。
近些年兴起的分布式流域模型，考虑了参数与过程的空间异质性，根据流域地形、土地利用/覆被、土壤等将流域离散成较多小单元，利用数学物理方程描述物质循环的各个过程，比较接近客观世界[2]。美国农业部开发的SWAT(soil and water assessment tool)分布式水文模型，其模型结构中包含了地表产汇流、侵蚀、泥沙输移等过程，能够模拟不同气候、土地利用/覆被、土壤条件下的水文、泥沙过程，加之源码公开，便于模型修改和开发，因此应用十分广泛[3-6]。SWAT模型以日为步长进行计算，应用SWAT模型模拟泥沙输移的研究中[7-9]常采用简化的Bagnold方程[10]计算河道输沙能力，其定义为河段上可输送的最大泥沙量是河道最大流速的函数：
| ${{\mathop{\rm conc}\nolimits} _{{\rm{sed}}, {\rm{ch}}, {\rm{mx}}}} = \alpha \cdot v_{{\rm{ch}}, {\rm{pk}}}^\beta .$ | (1) |
已有大量研究幂函数形式的输沙率-流量关系(Qs=aQb，其中Qs为输沙率，Q为流量)中系数a和指数b的影响因素。Yang等[15]研究长江流域不同水文站点输沙率-流量关系参数，认为水沙关系参数与所选河道地形相关。Hu等[16]研究长江流域输沙率-流量关系时空分布特征，结果表明：输沙率-流量幂函数方程系数和指数具有时空依赖性，人类活动会对方程中系数和指数产生较大影响。Bussi等[17]研究泰晤士河含沙量与流量关系指出，含沙量-流量幂函数方程系数与洪水发生次序相关。Zhang等[18]研究黄河中游输沙率-流量关系指出，输沙率-流量关系曲线呈现区域性特征，幂函数方程中的系数和指数取值与地貌类型相关，其中系数反映流域可侵蚀物质量，随流域面积的变化而变化。由此可知，在SWAT模型框架下，Bagnold方程中的系数和指数在不同河段取相同值会导致河道冲淤计算不准确。本研究通过实测的日尺度输沙率-流量关系转化为SWAT模型中Bagnold方程关系，考虑参数的空间依赖性，并选择黄河中游皇甫川流域对模型进行封闭，研究模型在泥沙模拟中的适用性。
1 SWAT模型泥沙输移模块封闭为实现Bagnold方程中的系数和指数按不同河段分别赋值，需研究系数和指数随流域特性的关系。Zhang等[18]研究指出在不同地貌类型区，输沙率-流量关系幂函数方程中的系数取值随流域面积的变化而变化，而指数与流域面积相关性不大。因此，本研究选择流域面积作为影响因子，建立系数和指数与面积的关系。
选择位于黄河中游砒砂岩地区的9个水文站的实测逐日悬移质输沙率-流量资料，通过幂函数拟合确定了输沙率-流量关系的系数a和指数b，拟合的幂函数方程决定系数R2的变化范围为0.53~0.86，说明拟合结果可以接受。在此基础上，建立了水文站控制面积和系数a之间的回归方程，如图 1a所示，表明系数a的取值与空间尺度有关。据此，可以采用系数a与流域面积的回归方程，结合流域内不同河段出口的控制面积，计算不同河段出口处的系数a。图 1b为指数b与水文站控制面积之间的相关关系，由图可知，指数b与控制面积之间无显著相关性，因此，本研究中指数b取平均值2.22。
 |
| 图 1 输沙率-流量幂函数方程系数a和指数b与控制面积关系 |
| 图选项 |
以上确定的系数a和指数b是基于输沙率-流量Qs=aQb关系得到，而SWAT模型中需要的是含沙量与流速concsed, ch, mx=α·vch, pkβ关系，因此需要将输沙率-流量关系转换为含沙量-流速关系。利用断面河相关系式流速与流量关系v=kQm，及输沙率-流量关系Qs=aQb的系数和指数可得出含沙量-流速关系的系数和指数，即系数α=a·k-(b-1)/m，指数β=(b-1)/m。对于河相关系式流速-流量关系中的系数k和指数m的计算，采用数字高程模型(digital elevation model, DEM)计算的不同河段几何参数确定。通过以上过程，可确定不同河段的系数值。通过对SWAT模型源代码的修改，实现对不同河段系数和指数的赋值。
2 研究区与数据2.1 研究区概况皇甫川流域位于黄河中游上部(见图 2)，集水面积3 246 km2。区域降水量350~450 mm，其中80%以上集中在6—9月，且以暴雨形式出现，导致洪水、产沙显著，是黄河的粗泥沙主要来源区之一。流域多年平均(1954—2015年)径流量1.255亿m3，输沙量0.387亿t，径流量和输沙量年际差异较大，年内分配也不均匀。径流量和输沙量均集中在6—9月份，分别占全年的81.9%和98.6%。皇甫川流域地貌类型复杂，分为砒砂岩丘陵沟壑区、黄土丘陵沟壑区和沙化黄土丘陵沟壑区。流域内设有2个水文站，分别为皇甫站(控制面积3 175 km2)和沙圪堵站(控制面积1 351 km2)。
 |
| 图 2 (网络版彩图)研究区地理位置和雨量站、水文站及气象站分布图 |
| 图选项 |
2.2 数据来源SWAT模型建模所需的数据涉及空间数据和属性数据，包括水文气象数据、DEM数据、土壤类型数据、土地利用/覆被数据等。皇甫、沙圪堵水文站1978—1989年径流、输沙资料来自黄河流域历年水文年鉴，气象数据采用河曲站逐日气温、风速、相对湿度及日照时数等，数据来源于中国气象局气象数据中心。DEM数据(30 m)来自中国科学院计算机网络信息中心地理空间数据云平台。土地利用/覆被图采用1978年由Landsat影像通过非监督分类法解译的流域土地类型，分别为草地、耕地、林地、居民地、沙地、水体和未利用地7类。NDVI数据采用1983年AVHRR数据产品并采用像元二分模型计算植被覆盖度[19]，与土地利用图叠加，获得不同植被覆盖度类型的土地利用。最终使用的土地利用/覆被图类型分为9类，分别为低覆盖度草地(面积占比74.8%)、中覆盖度草地(4.8%)、低覆盖度林地(0.8%)、中覆盖度林地(0.03%)、耕地(3.2%)、居民地(0.27%)、沙地(4.9%)、水体(3.5%)和未利用地(7.7%)。土壤数据来自黄土高原生态环境数据库，土壤类型共分为7类，分别为黑垆土(面积占比0.63%)、栗钙土(3.31%)、黄绵土(19.31%)、新积土(29.8%)、风沙土(17.99%)、粗骨土(24.4%)和石质土(4.56%)。
3 模型检验3.1 模型参数确定基于物理过程的SWAT分布式水文模型考虑了气候、地形地貌、土壤植被等的空间差异性，通过数学物理方程描述坡面产流、产沙和河道汇流输沙等物理过程，封闭这些物理过程模型的关键参数或物理量，依赖于试验或观测资料率定，或者参考前人研究成果给出。
SWAT模型中坡面产沙计算采用修正通用土壤流失方程(modified universal soil loss equation, MUSLE)：
| $\begin{array}{l}{\rm{sed}} = 11.8 \cdot {\left( {{Q_{{\rm{surf }}}} \cdot {q_{{\rm{peck }}}} \cdot {\rm{ area}}{{\rm{ }}_{{\rm{hru}}}}} \right)^{0.56}} \cdot \\K \cdot C \cdot P \cdot {\rm{LS}} \cdot {\rm{CFRG}}{\rm{.}}\end{array}$ | (2) |
与坡面产沙相关的因子中，K与土壤本身性质相关，LS与地形相关，CFRG与土壤颗粒组成相关，这些参数理论上根据土壤条件确定，其值应该保持不变。而作物覆盖因子C和耕作因子P与土地利用/覆被有关，应随流域土地利用/覆被的改变而改变。前人研究已经给出了不同土地利用/覆被类型下的C和P值。江忠善等[20]根据径流小区观测资料，分别建立了草地和林地C与植被覆盖因子veg之间的关系。
草地：
| ${C = {{\rm{e}}^{ - 0.0418({\rm{veg}} - 5)}}.}$ | (3) |
| ${C = {{\rm{e}}^{ - 0.0085{{({\rm{veg}} - 5)}^{1.5}}}}.}$ | (4) |
通过坡面泥沙参数的选取，可控制坡面产沙精度；采用实测资料对不同河段Bagnold方程中系数和指数进行确定，空间分布式赋值可控制流域内部冲淤过程精度；在此基础上对流域出口输沙过程进行模拟，从而保证流域产沙、输移、冲淤过程的准确性。
在Li等[14]的研究中，采用7年(1978—1984年)校准期和5年(1975—1989年)验证期来评估SWAT在皇甫川流域的适用性。模型模拟结果显示，月径流的NSE在率定期和验证期分别为0.83和0.81，月输沙的NSE在率定期和验证期分别为0.8和0.83。该研究中河道输沙Bagnold方程的系数和指数取值在各河段相同，其值分别为0.007和1.49。在该研究基础上，本研究对已构建皇甫川SWAT模型中关于河道输沙Bagnold方程的系数和指数进行空间分布赋值，并对影响产沙输沙的主要因子进行率定，包括修正土壤流失方程(MUSLE)中的耕作因子P、作物覆盖因子C，以及河道的侵蚀和覆盖因子CH_COV1和CH_COV2。采用皇甫站1978—1984年实测输沙量资料作为校准期来优化土壤侵蚀和输移的参数，1985—1989年作为验证期来评估模型适用性。
模型参数率定完成后，采用决定系数R2和NSE评价模型的模拟精度。模型评价指标计算如下：
| ${R^2} = \frac{{{{\left[ {\sum\limits_i {\left( {{W_{{\rm{si}}}} - \overline {{W_{\rm{s}}}} } \right)} \left( {W_{{\rm{si}}}^\prime - \bar W_{\rm{s}}^\prime } \right)} \right]}^2}}}{{\sum\limits_i {{{\left( {{W_{{\rm{si}}}} - \overline {{W_{\rm{s}}}} } \right)}^2}} \sum\limits_i {{{\left( {W_{{\rm{si}}}^\prime - \bar W_{\rm{s}}^T} \right)}^2}} }}, $ | (5) |
| ${\rm{NSE}} = 1 - \frac{{\sum\limits_i {{{\left( {{W_{{\rm{si}}}} - W_{{\rm{si}}}^\prime } \right)}^2}} }}{{\sum\limits_i {{{\left( {{W_{{\rm{si}}}} - {{\bar W}_{\rm{s}}}} \right)}^2}} }}.$ | (6) |
3.2 模型率定与验证基于SWAT模型，采用封闭后的泥沙模拟算法进行泥沙模拟计算。封闭后的SWAT模型皇甫站泥沙模拟结果见图 3。率定期(1978—1984年)和验证期(1985—1989年)实测值与模拟值趋势变化表现出了较好的一致性，基本反映了皇甫川流域输沙量的变化趋势。率定期日输沙量的NSE为0.70，决定系数R2为0.74；验证期日输沙量的NSE为0.66，决定系数R2为0.68。表明封闭后的SWAT模型对皇甫川流域日输沙量的模拟效果较好，满足精度要求。
 |
| 图 3 (网络版彩图)皇甫站日输沙量模拟值与实测值对比 |
| 图选项 |
为比较流域内部河段模拟效果，选择流域内部沙圪堵站的1978—1989年输沙量实测值与模拟值进行对比(见图 4)。由图可知，封闭后的模型输沙量模拟值和实测值的趋势变化较一致，率定期日输沙量的NSE为0.78，R2为0.78；验证期日输沙量的NSE为0.72，R2为0.74。表明封闭后的SWAT模型能够模拟流域内部输沙过程。
 |
| 图 4 (网络版彩图)沙圪堵站日输沙量模拟值与实测值对比 |
| 图选项 |
同时比较了Li等[14]前期研究中采用未封闭SWAT模型对沙圪堵站输沙量的模拟结果(见图 5)。由图可知，未封闭的SWAT模型率定期日输沙量的NSE仅为0.36，R2为0.39；验证期NSE和R2均为0.26。与封闭后的模型日输沙量模拟效果对比可知，封闭后的模型日输沙量模拟精度较未封闭模型有所提高，表明封闭后的SWAT模型提高了模拟流域内部输沙过程的精度，使模型输沙模拟的适用性有所提升。本研究的意义在于，当流域只有一个水文站时，只能通过该水文站水沙资料进行模型参数率定，问题是虽然模型可以通过参数率定将该水文站的水沙过程进行还原，但关注该水文站以上流域内部水沙输移过程时，往往误差较大，本研究可为解决这一问题提供方法参考。同时，在进行输沙量模拟时，需要选择与时间尺度匹配的输沙率-流量关系进行计算。
 |
| 图 5 (网络版彩图)沙圪堵站未封闭模型输沙量模拟值与实测值对比 |
| 图选项 |
4 结论本文将SWAT模型用于河道日均输沙率过程模拟，采用黄河中游砒砂岩区实测的日均输沙率-日均流量关系封闭日均最大含沙量，建立了模型参数随流域面积变化关系，体现了输沙率-流量关系的空间不均匀性、日内流量过程的非恒定性。利用黄河中游皇甫川流域嵌套水文站点的实测资料检验表明，封闭后的SWAT模型显著提高了模拟效率，具有良好适用性。
参考文献
| [1] | 冉大川, 刘斌, 罗全华, 等. 泾河流域水沙变化水文分析[J]. 人民黄河, 2001, 23(2): 9-11. RAN D C, LIU B, LUO Q H, et al. Hydrological analysis of water and sediment evolution in Jing River basin[J]. Yellow River, 2001, 23(2): 9-11. DOI:10.3969/j.issn.1000-1379.2001.02.004 (in Chinese) |
| [2] | 徐宗学, 程磊. 分布式水文模型研究与应用进展[J]. 水利学报, 2010, 41(9): 1009-1017. XU Z X, CHENG L. Progress on studies and applications of the distributed hydrological models[J]. Journal of Hydraulic Engineering, 2010, 41(9): 1009-1017. (in Chinese) |
| [3] | PATI A, SEN S, PERUMAL M. Modified channel-routing scheme for SWAT model[J]. Journal of Hydrologic Engineering, 2018, 23(6): 04018019. DOI:10.1061/(ASCE)HE.1943-5584.0001657 |
| [4] | FURL C, SHARIF H, JEONG J. Analysis and simulation of large erosion events at central Texas unit source watersheds[J]. Journal of Hydrology, 2015, 527: 494-504. DOI:10.1016/j.jhydrol.2015.05.014 |
| [5] | ADHIKARI U, NEJADHASHEMI A P. Impacts of climate change on water resources in Malawi[J]. Journal of Hydrologic Engineering, 2016, 21(11): 05016026. DOI:10.1061/(ASCE)HE.1943-5584.0001436 |
| [6] | MARHAENTO H, BOOIJ M J, RIENTJES T H M, et al. Attribution of changes in the water balance of a tropical catchment to land use change using the SWAT model[J]. Hydrological Processes, 2017, 31(11): 2029-2040. DOI:10.1002/hyp.11167 |
| [7] | DURU U, ARABI M, WOHL E E. Modeling stream flow and sediment yield using the SWAT model: A case study of Ankara River basin, Turkey[J]. Physical Geography, 2018, 39(3): 264-289. DOI:10.1080/02723646.2017.1342199 |
| [8] | VERMA S, BHATTARAI R, BOSCH N S, et al. Climate change impacts on flow, sediment and nutrient export in a great lakes watershed using SWAT[J]. Clean-Soil Air Water, 2015, 43(11): 1464-1474. DOI:10.1002/clen.201400724 |
| [9] | YESUF H M, ASSEN M, ALAMIREW T, et al. Modeling of sediment yield in Maybar gauged watershed using SWAT, northeast Ethiopia[J]. Catena, 2015, 127: 191-205. DOI:10.1016/j.catena.2014.12.032 |
| [10] | WILLIAMS J R. SPNM, a model for predicting sediment, phosphorus, and nitrogen yields from agricultural basins[J]. Water Resources Bulletin, 1980, 16(5): 843-848. DOI:10.1111/j.1752-1688.1980.tb02497.x |
| [11] | QIU L J, ZHENG F L, YIN R S. SWAT-based runoff and sediment simulation in a small watershed, the loessial hilly-gullied region of China: Capabilities and challenges[J]. International Journal of Sediment Research, 2012, 27(2): 226-234. DOI:10.1016/S1001-6279(12)60030-4 |
| [12] | DUTTA S, SEN D. Application of SWAT model for predicting soil erosion and sediment yield[J]. Sustainable Water Resources Management, 2018, 4(3): 447-468. DOI:10.1007/s40899-017-0127-2 |
| [13] | YU X, XIE X H, MENG S S. Modeling the responses of water and sediment discharge to climate change in the upper Yellow River basin, China[J]. Journal of Hydrologic Engineering, 2017, 22(12): 05017026. DOI:10.1061/(ASCE)HE.1943-5584.0001590 |
| [14] | LI E H, MU X M, ZHAO G J, et al. Effects of check dams on runoff and sediment load in a semi-arid river basin of the Yellow River[J]. Stochastic Environmental Research and Risk Assessment, 2017, 31(7): 1791-1803. DOI:10.1007/s00477-016-1333-4 |
| [15] | YANG G F, CHEN Z Y, YU F L, et al. Sediment rating parameters and their implications: Yangtze River, China[J]. Geomorphology, 2007, 85(3-4): 166-175. DOI:10.1016/j.geomorph.2006.03.016 |
| [16] | HU B Q, WANG H J, YANG Z S, et al. Temporal and spatial variations of sediment rating curves in the Changjiang (Yangtze River) basin and their implications[J]. Quaternary International, 2011, 230(1-2): 34-43. DOI:10.1016/j.quaint.2009.08.018 |
| [17] | BUSSI G, DADSON S J, BOWES M J, et al. Seasonal and interannual changes in sediment transport identified through sediment rating curves[J]. Journal of Hydrologic Engineering, 2017, 22(2): 06016016. DOI:10.1061/(ASCE)HE.1943-5584.0001466 |
| [18] | ZHANG S Y, CHEN D, LI F X, et al. Evaluating spatial variation of suspended sediment rating curves in the middle Yellow River basin, China[J]. Hydrological Processes, 2018, 32(11): 1616-1624. DOI:10.1002/hyp.11514 |
| [19] | 高健健, 穆兴民, 孙文义. 1981—2012年黄土高原植被覆盖度时空变化特征[J]. 中国水土保持, 2016(7): 52-56. GAO J J, MU X M, SUN W Y. Temporal and spatial variation characteristics of vegetation coverage in the loess plateau from 1981 to 2012[J]. Soil and Water Conservation in China, 2016(7): 52-56. DOI:10.3969/j.issn.1000-0941.2016.07.022 (in Chinese) |
| [20] | 江忠善, 王志强, 刘志. 黄土丘陵区小流域土壤侵蚀空间变化定量研究[J]. 土壤侵蚀与水土保持学报, 1996, 2(1): 1-9. JIANG Z S, WANG Z Q, LIU Z. Quantitative study on spatial variation of soil erosion in a small watershed in the loess Hilly Region[J]. Journal of Soil and Water Conservation, 1996, 2(1): 1-9. (in Chinese) |
| [21] | 李斌兵, 郑粉莉, 龙栋材, 等. 基于GIS纸坊沟小流域土壤侵蚀强度空间分布[J]. 地理科学, 2009, 29(1): 105-110. LI B B, ZHENG F L, LONG D C, et al. Spatial distribution of soil erosion intensityin Zhifanggou small watershed based on GIS[J]. Scientia Geographica Sinica, 2009, 29(1): 105-110. DOI:10.3969/j.issn.1000-0690.2009.01.017 (in Chinese) |
| [22] | 郑捷, 李光永, 韩振中, 等. 改进的SWAT模型在平原灌区的应用[J]. 水利学报, 2011, 42(1): 88-97. ZHENG J, LI G Y, HAN Z Z, et al. Application of modified SWAT model in plain irrigation district[J]. Journal of Hydraulic Engineering, 2011, 42(1): 88-97. (in Chinese) |
