1. 引言
在我国,水土流失面积达294.9万平方千米,其中东北黑土区作为我国重要的粮食生产基地,因长期坡耕地开垦、降水集中等因素,黑土流失速率远超成土速率,已成为制约区域粮食安全与生态保护的核心瓶颈[1]。径流侵蚀能量作为驱动土壤侵蚀的关键动力因子,直接决定了径流对土壤颗粒的剥离能力与搬运效率,其大小与空间分布特征是揭示水土流失机制、制定精准水土保持措施的核心依据。不同于“降水侵蚀力”聚焦于降水对地表的冲击作用,径流侵蚀能量更侧重于地表径流形成后,通过水力剪切、拖拽等物理过程对土壤产生的侵蚀作用,尤其在流域尺度下,受地形、土壤、土地利用及水文过程的综合影响,径流侵蚀能量呈现显著的空间异质性。李占斌[2]等学者提出的径流侵蚀功率模型,通过采纳洪峰流量模数这一核心参数,实现了对流域下垫面状况与水文动力过程的全面融合,有效描绘了水力侵蚀的动力学原理。鲁克新[3] [4]等人的深入研究显示,与传统的降雨侵蚀力指标相比,径流侵蚀功率在流域尺度次暴雨产沙模型的侵蚀动力因子的应用上显现出更高的适宜性。基于此理论框架构建的次暴雨输沙模型进一步验证,径流侵蚀功率与流域输沙模数之间呈现显著的幂函数相关性,即径流侵蚀强度的增加与流域输沙模数的增长呈现出一致性趋势。因此,精准刻画流域尺度径流侵蚀能量的空间分布规律,是理解“水文过程–侵蚀动力–土壤流失”耦合关系的关键环节,也是提升水土流失模拟与预测精度的前提。
流域作为水土流失发生与调控的基本单元,其水文循环与侵蚀过程的复杂性要求研究需兼顾“过程模拟”与“空间解析”。传统径流侵蚀能量研究多基于小区观测或经验公式(如基于径流深、流速的能量计算模型),虽能获取局部区域的能量特征,但难以反映流域内下垫面异质性(如坡度、土壤质地、土地利用类型)对径流侵蚀能量的影响,且空间扩展性差,无法满足大尺度流域水土保持规划的需求。随着分布式水文模型的发展,以SWAT (Soil and Water Assessment Tool)为代表的模型凭借其“分布式结构–物理过程驱动–多尺度适配”的优势,能够整合气象、地形、土壤、土地利用等多源数据,实现对流域水文过程(产流、汇流)与土壤侵蚀过程的动态模拟,为流域尺度径流侵蚀能量的定量计算与空间分布解析提供了有效工具[5]。国内外学者多使用SWAT模型进行流域水沙研究。龚珺夫研究团队[6] [7]运用SWAT水文模型对无定河与延河流域水沙运移规律及年径流侵蚀功率的时空变异特性进行了系统模拟,研究结果表明:研究区域内年径流侵蚀功率的空间分布格局呈现出明显的差异性特征。Fetene等[8]识别埃塞俄比亚里布流域的侵蚀热点区域以确定保护措施优先级,并评估减少土壤侵蚀的最佳管理措施。
本研究基于SWAT模型对拉林河流域的模拟结果,通过将单场次暴雨事件的径流侵蚀功率参数进行年尺度扩展,量化了各子流域的年度径流侵蚀功率值,系统解析了其空间分异规律与尺度效应规律,旨在揭示不同空间格局与尺度条件下水力侵蚀过程的形成机制,进而为拉林河流域水土资源的差异化管控、水土保持工程的科学布局及生态恢复策略的制定提供理论支撑。
2. 材料及方法
2.1. 研究区概况
拉林河(图1)属松花江流域二级支流,发源于吉林省东部张广才岭西麓(主峰老爷岭,海拔1696 m),自东南向西北流经吉林省吉林市(舒兰市、永吉县)、黑龙江省哈尔滨市(五常市、双城区、阿城区),最终在黑龙江省肇东市五站镇附近汇入松花江干流。流域地理坐标介于东经126˚00′~128˚30′、北纬44˚30′~46˚15′之间,总流域面积约1.98万平方千米(其中吉林省境内占42%,黑龙江省境内占58%),干流全长约448 km,河道平均比降为0.25% [9]。
注:该图基于自然资源部标准底图服务网站下载的审图号为GS(2019)4345号的标准地图制作,底图无修改。
Figure 1. Location of the Lalin River Basin and distribution map of meteorological and hydrological stations
图1. 拉林河流域位置及气象水文站点分布图
年均气温2.5~4.0℃,年内温差大。1月均温−18℃至−22℃ (极端最低温−40.9℃),7月均温22℃至24℃ (极端最高温37.5℃);大于等于10℃积温为2600~2900℃,无霜期130~150天,满足一年一熟的农业生产需求(主要作物为玉米、水稻、大豆)。年均降水量500~700 mm,空间上呈“上游多、下游少”的分布规律(上游山区年均降水650~700 mm,下游平原500~550 mm);时间分配极不均匀,6~8月汛期降水占全年的60%~70%,且多以短时强降雨(日降水量 ≥ 50 mm)形式出现,单次暴雨量可达100~150 mm,易形成高强度地表径流,是驱动土壤侵蚀的主要气候因子。年均蒸发量(20 cm蒸发皿)为1100~1300 mm,春旱(4~5月)期间蒸发量可达降水量的3~5倍,易导致土壤表层干燥疏松,抗蚀性下降;年均风速2.5~3.0 m/s,春季大风(风速 ≥ 10 m/s)日数可达15~20天,偶发风蚀,但流域侵蚀仍以水力侵蚀为主。流域内设有多个水文与泥沙观测站,其中蔡家沟站(位于黑龙江省五常市,控制流域面积1.56万平方千米)是干流核心控制站,拥有1956年至今的连续月径流、泥沙观测数据,为SWAT模型的校准与验证提供了关键数据支撑[10]。
2.2. 数据来源与处理
拉林河流域SWAT模型的构建依赖于多源基础数据的整合与预处理,所需数据集包括数字高程模型(DEM)、土地利用类型、土壤属性、气象观测资料以及水文监测记录等,各类数据的具体信息及来源详见表1。考虑到模型对空间数据一致性的严格要求,本研究将所有栅格数据统一转换为WGS_1984_UTM_Zone_52N投影坐标系,并采用GCS_WGS_1984地理坐标系统作为基准。拉林河流域的高程信息、土地利用空间分布及土壤类型划分情况参见图2~4。
Table 1. Research Data and Sources
表1. 研究数据及来源
输入数据 |
精度 |
数据来源 |
数字高程模型 |
30 m |
欧洲航天局 |
土地利用数据 |
30 m |
中国资源与环境中心 |
土壤数据集 |
1000 m |
世界土壤数据库 |
降雨量、风速、气温、太阳辐射、相对湿度等 |
逐日 |
中国大气同化驱动数据集(CMADS V1.0) |
径流量 |
逐月 |
蔡家沟水文站 |
注:该图基于自然资源部标准底图服务网站下载的审图号为GS(2019)4345号的标准地图制作,底图无修改。
Figure 2. DEM of the Lalin River Basin
图2. 拉林河流域DEM
注:该图基于自然资源部标准底图服务网站下载的审图号为GS(2019)4345号的标准地图制作,底图无修改。
Figure 3. Land use map of the Lalin River Basin
图3. 拉林河流域土地利用图
注:该图基于自然资源部标准底图服务网站下载的审图号为GS(2019)4345号的标准地图制作,底图无修改。
Figure 4. Soil distribution map of the Lalin River Basin
图4. 拉林河流域土壤分布图
2.3. 模型构建
本研究选取ARCSWAT2012软件为研究平台,在运用DEM数据进行河网提取过程中,将最小集水面积阈值设定为36981 hm2,由此将研究区域拉林河划分为19个子单元。在该单元划分过程中,系统自动生成包含河网信息的Reach图层及子流域信息的Watershed图层。基于上述两个图层数据,可直接获取各子流域出口断面上游河道长度及其控制面积。此外,通过河长与控制面积的比值计算,可获得河网密度这一水文特征参数。在将流域坡度分级为五个等级(0~5%、5%~10%、10%~15%、15%~20%、20%~99.99%)的基础上,融合土壤类型图层与土地利用(LUCC)图层信息,采用[11]优势地面覆盖与优势土壤类型组合法构建HRU (水文响应单元)。设定土地利用和土壤类型的阈值为5%,以此标准筛选并剔除子流域内面积百分比低于阈值的次要土地利用及土壤类型,从而确保仅保留覆盖超过阈值的土地单元。随后,对剩余的土地利用和土壤类型面积进行重新分配,以确保模拟覆盖子流域的全域面积[12] [13],随后,19个子流域被细分为326个水文响应单元(HRU)。进而,将自建的土壤数据库、气象数据库及天气发生器等相关属性数据导入SWAT模型,以执行模拟分析。在经过模型适用性验证并确认其适用于研究区域之后,本研究基于模型模拟结果展开区域层面的研究。模型输出的数据包括模拟期内各子流域的月径流量、泥沙负荷及污染物浓度等参数,本研究主要选取月径流量数据进行分析。
为保障SWAT水文模型的模拟精度,本研究选取拉林河蔡家沟水文站2008年作为模型的预热时段,并采用2009~2013年各月实测平均径流数据对模型参数进行优化率定,随后利用2014~2016年同期实测月均流量数据对模型模拟效能进行独立验证。
2.4. 敏感性分析
为保障模型预测精度并优化计算效能,本研究通过参数敏感性分析对关键影响因素进行了筛选与识别。研究中采用了LH-OAT方法[14]开展参数敏感性评估,该方法有效整合了拉丁超立方抽样(LH, Latin Hypercube)与单因素逐次变化(OAT, One-factor-At-a-Time)分析法的双重优势,其技术路径为对LH抽样获取的每个样本点执行单因素扰动分析,并通过各部分敏感性指标的平均值表征整体敏感程度。基于p值和t检验值构建参数敏感性评价体系,其中p值趋近于0且t绝对值增大标志着参数敏感性的增强。研究最终筛选出对模型模拟结果具有较高敏感性的参数,用于后续的模型参数率定与验证。
2.5. 模型率定与验证
在模型运行完毕后,需要对模拟输出的数据进行校准,并对模型的可靠性进行检验。本研究选取蔡家沟水文站2009至2013年的月径流数据作为模型校准的依据,而2014至2016年的月径流数据则用于模型的验证。在识别出研究区域内的关键敏感参数后,本研究运用SWAT-CUP软件中的SUFI-2算法对这些参数进行了校准优化,该参数率定工作主要在模型的率定阶段完成。在获取最佳参数校正值后,将其代入验证期数据进行验证分析。本研究通过决定系数R2和纳什效率系数NSE (Nash-Sutcliffe efficiency Coefficient)两项指标对模型在拉林河流域的应用效能进行了科学评估,其中R2表征变量间的关联程度,其数值越趋近于1,表明关联性越强;而NSE则用于量化实测值与模拟值的一致性程度,其数值越接近1,说明模型拟合精度越高。依据相关文献[15]-[17]的判定标准,当R2 > 0.6且NSE > 0.5时,可判定模拟结果具有较高的可靠性与精确度。当模型在率定期与验证期均表现出良好的模拟性能时,则表明所确定的参数校正值可有效适用于该研究区域的模拟分析。
2.6. 流域径流侵蚀功率
本论文采纳了李占斌等人[2]所提出的次暴雨径流侵蚀功率理论,其表达式如下:
(1)
(2)
(3)
式中E为次暴雨径流侵蚀功率,单位为m4/(s∙km2),具有功率的量纲;
为洪峰流量模数,单位为m3/(s∙km2),其数值等于次暴雨洪水洪峰流量Qm (m3/s)与流域面积A (km2)的比值;为次暴雨平均径流深,单位为m;Q为次暴雨平均流量,单位为m3/s;
为时段长度,此处为1秒。
径流深H和洪峰流量Qm是常用的反映次暴雨洪水过程的关键参数,但它们均不足以全面反映降雨过程与地表条件交互作用的影响。次暴雨径流侵蚀功率则结合了两者的优势,能够更准确地反映水力侵蚀的作用。本研究将次暴雨降雨侵蚀功率的概念扩展至年度尺度,旨在探究侵蚀能量在流域尺度上的空间分布。通过SWAT模型对2008年至2016年间19个子流域的月径流量进行模拟,并将每个子流域每年的径流过程简化为单次月尺度暴雨事件,以年平均径流深作为
、以年内流量最大月的径流量作为
,将二者作为参数代入公式(1)至(3),从而获得改进的年尺度径流侵蚀功率计算公式(4)。
(4)
(5)
(6)
式中
为年径流侵蚀功率,单位为m4/(s∙km2);
为年内月最大流量模数,单位为m3/(s∙km2),其数值等于年内最大月流量
与流域控制面积
的比值;
为年平均径流深,单位为m;
为年平均流量,单位为m3/s,
的计算方法为年内各月径流量之和除以12;为时段长,在本研究中以月为单位,按平均每月30天计算,
= 2592 × 103 s。
3. 结果与分析
3.1. 模型率定验证
基于拉林河干流蔡家沟水文站2008至2016年系列观测资料,本研究构建了拉林河流域SWAT水文模型。选取2009~2013年月均流量序列作为模型参数率定时段,运用LH-OAT方法开展参数敏感性识别,最终筛选出10个敏感性参数参与模型率定,具体参数取值详见表2所示。
Table 2. Sensitivity ranking of SWAT model parameters in the Lalin River Basin
表2. 拉林河流域SWAT模型参数敏感性排名
参数名称 |
参数含义 |
敏感排名 |
调参方法 |
参数范围 |
最佳校准值 |
ALPHA_BF |
基流消退系数 |
1 |
V |
−0.17~0.81 |
0.40 |
CN2 |
SCS径流曲线数 |
2 |
R |
0.24~0.47 |
0.34 |
EPCO |
植物蒸腾补偿系数 |
3 |
V |
0.31~1.21 |
0.73 |
HRU_SLP |
水文响应单元的平均坡度 |
4 |
R |
0.47~1.11 |
0.99 |
SFTMP |
降雪气温 |
5 |
V |
−3.83~1.74 |
−0.78 |
GWQMN |
浅层地下水径流系数 |
6 |
R |
−1.94~-1.18 |
−1.31 |
CANMX |
最大冠层蓄水量 |
7 |
R |
−0.57~0.39 |
0.26 |
USLE_P |
通用土壤流失方程水土保持措施因子 |
8 |
R |
−1.10~-0.41 |
−0.81 |
SPCON |
河道中泥沙输移的线性参数 |
9 |
V |
−0.003~0.004 |
−0.0028 |
SOL_BD |
表层土壤容重 |
10 |
R |
−0.68~0.047 |
−0.34 |
表2数据显示,在拉林河流域径流模拟过程中,基流消退系数对模拟结果表现出最为显著的影响,该参数主要表征地下水对河道的补给能力并控制基流退水速率。SCS径流曲线数CN2位居第二,该参数与流域径流量呈现正相关关系,即CN2值增大会导致地表径流量上升。植物蒸腾补偿系数EPCO则主要调控植物在水分胁迫条件下从深层土壤获取水分以满足蒸腾需求的能力。此外,水文响应单元的平均坡度HRU_SLP、降雪气温阈值SFTMP以及浅层地下水径流系数GWQMN等参数对模拟结果亦表现出较强的敏感性[18] [19]。
本研究将蔡家沟水文站所获取的月均流量观测数据划分为预热阶段(2008年)、参数率定阶段(2009年至2013年)以及验证阶段(2014年至2016年)三个不同时期。采用SWAT-CUP软件,以参数率定阶段的数据为基础,对所选取的十个关键参数进行优化调整,通过迭代优化,确保模型的决定系数R2超过0.6且纳什效率系数NSE超过0.5,以此作为模型模拟精度达到可接受水平的标准。继而对调整后的参数集实施应用,执行SWAT模型的验证流程,并将所得月度流量模拟数据与实际观测值进行对照分析,复算R²与NSE指标。在验证阶段,若R2和NSE指标均符合既定阈值要求,则可断定SWAT模型输出结果的可靠性得到保障,进而依据这些指标确定的最优参数配置适用于后续月份流量模拟的计算过程。
通过对蔡家沟水文站2009~2013年模型率定期与2014~2016年验证期的月均流量实测值与模拟值进行对比分析(图5),结果显示:率定期决定系数R2为0.65,纳什效率系数NSE为0.64;验证期R2达0.75,NSE为0.7。根据所列评价指标,SWAT模型在模拟拉林河流域月流量序列方面表现良好,显示出该模型在研究区域内具有较高的适用性。通过该模型的分析,成功收集了2009至2016年间拉林河流域入江口以及19个子流域出口断面的月径流序列及年径流深度等关键水文数据。
Figure 5. Comparison chart of simulated and measured average monthly runoff at Caijiagou Station
图5. 蔡家沟站月均径流量模拟与实测对比图
3.2. 径流侵蚀功率空间分布
SWAT模型模拟生成了拉林河流域19个子流域在2009~2016年期间逐月的径流数据序列。基于该输出结果,可识别各子流域年内最大月径流量
,并通过公式(5)进一步计算得出月最大流量模数
;采用公式(6)可获得各子流域的年均径流深
;随后,借助改进后的年径流侵蚀功率计算模型(公式4),确定了研究区内各子流域逐年径流侵蚀功率值。为揭示区域水文侵蚀特征的一般性规律,对2009~2016年共8年的计算结果进行算术平均处理,最终获得各子流域的多年平均径流侵蚀功率(如图6所示)。
注:该图基于自然资源部标准底图服务网站下载的审图号为GS(2019)4345号的标准地图制作,底图无修改。
Figure 6. Map of annual average runoff erosion power in the watershed
图6. 流域年均径流侵蚀功率分布图
可以明显看出,拉林河流域径流侵蚀功率的空间分异规律呈现出显著的梯度特征:干流区域高于支流区域,下游地区高于上游地区,西北部高于东南部;具体表现为干流沿线区域自河源向下游方向,径流侵蚀强度呈递增态势;源头地带的侵蚀功率值明显低于干流主体区域;从空间格局来看,流域西北部各支流子单元的侵蚀动能普遍高于东南部子流域。
3.3. 径流侵蚀功率的空间尺度效应
然而,对于上面得到的径流侵蚀功率空间分布,2号和8号子流域却呈现出“支流大,干流小”的反常现象。因此本研究为揭示可能存在的空间尺度效应对水文过程的影响,选取子流域出口断面集水区面积作为空间量化参数[20] [21],建立其与径流侵蚀功率的耦合关系模型。基于Watershed矢量数据,系统解算各子流域出口断面上游集水范围,将所得集水区面积作为空间尺度特征值,与各单元多年平均径流侵蚀功率进行回归建模,具体数值关系如图7所示。
Figure 7. Fitting analysis chart of sub-basin area and multi-year average runoff erosion power
图7. 子流域面积与多年平均径流侵蚀功率拟合分析图
可见,径流侵蚀功率与控制面积呈现出幂函数关系:
(7)
在模型分析中,R2指数达到0.6202,表明随着控制区域的扩大,径流侵蚀能量呈现下降趋势,且下降速度逐渐减缓,最终趋于稳定,这暗示着在汇流面积达到某一临界点后,进一步增加面积对径流侵蚀能量的影响趋于微弱,其敏感性显著减弱。为精确界定这一临界面积,本研究对拟合曲线进行了导数求解,得到导数方程:
(8)
若
,则表明子流域出口断面的多年平均径流侵蚀功率变化速率逐渐趋于恒定状态。经计算分析,当子流域的集水面积超过150平方公里时,
将小于0.01,此时子流域出口断面的多年平均径流侵蚀功率不再随子流域面积的扩大而显著降低。同时本研究又通过Python使用变化点检测(Change-Point Analysis)模型来识别幂函数关系中的斜率变化点加以验证。从而可以认为拉林河流域中,影响多年平均径流侵蚀功率的空间阈值位于子流域集水面积150平方公里处。
基于上述研究表明,拉林河流域径流侵蚀功率表现出明显的尺度依赖性特征,其年均值随子流域面积的扩张呈现出先急剧下降后趋于平稳的变化模式,具体表现为:当集水面积小于150 km2时,径流侵蚀功率随面积增加而显著降低;而当面积超过150 km2这一临界阈值后,其递减速率明显放缓。这一发现与鲁克新等[3]建立的流域次暴雨输沙模型理论相吻合,该模型揭示了流域次暴雨径流侵蚀功率与输沙模数之间存在显著的幂函数相关性,即前者数值增大将导致后者同步提升。因此,在实施拉林河流域生态修复工程时,建议优先对出口控制断面集水面积小于150 km2的子流域进行重点治理,从而实现水土保持效益的最大化。
4. 讨论
4.1. 年径流侵蚀功率空间分布规律
在本项研究以及对龚珺夫等人[6] [7]针对无定河和延河流域的径流侵蚀功率空间分布所进行的探讨中,揭示了不同流域年径流侵蚀功率在地域上的差异之处。研究表明,河流区域年径流冲刷能力的地域分布受到地表特征与环境条件的制约,从而导致各流域呈现出其特有的空间分布模式。
4.2. 年径流侵蚀功率的空间尺度效应
为更全面地考察子流域出口位置处的径流侵蚀功率与其集水规模之间的内在联系,本研究相较于龚珺夫所采用的研究方法,转而采纳李占斌[2]与程圣东[22]提出的年度径流侵蚀功率评估体系,并将年度径流深度作为反映径流侵蚀动能的核心变量,以期对流域侵蚀动力学特征进行更为精确的量化描述。结果表明,拉林河子流域的年均径流侵蚀功率与其出口控制断面集水范围之间呈现出明显的幂函数关系特征。值得注意的是,不同流域的径流侵蚀功率对集水面积的响应阈值存在明显差异。在拉林河流域这一特定研究区域内,当子流域集水面积低于临界阈值时,其空间分布特征多表现为各支流汇入干流的出口位置,对应于上述图中径流侵蚀功率高的子流域,这一发现对拉林河流域生态治理实践提供了重要的参考价值。
4.3. 不足与展望
本调研运用SWAT模型对拉林河区域年内径流侵蚀强度的空间分布特性及规模影响开展了模拟研究,然而,对于流域多年平均径流侵蚀功率空间分布机制并未进行深入剖析。特别是,关于坡度比、土地利用及土壤类型等变量对侵蚀功率分布的影响尚未得到充分探究。此外,本研究所采用的气象数据源自中国大气同化驱动数据集,如能替换为精度更高的中国气象网站实测数据,模型模拟精度有望提升。同时,本研究未能构建拉林河流域产沙输沙过程的模拟模型,从而限制了径流侵蚀功率指标在流域输沙模数计算中的应用价值。为弥补此研究局限,后续工作应着重拓展相关理论体系,以期提升流域侵蚀产沙预测精度并为区域生态治理实践提供更为完善的理论依据。
5. 结论
(1) 本研究基于拉林河流域构建的SWAT水文模型,在率定期(2009~2013年)中,决定系数R²和纳什效率系数NSE分别达到0.65与0.64,而在验证期(2014~2016年)则分别为0.75和0.70。研究结果表明,所采用的SWAT模型在本研究所涉及的地域内展现了良好的契合度,其输出数据与拉林河流域的实测径流特性呈现出较好的一致性。
(2) 拉林河流域的年均径流侵蚀功率表现出显著的空间差异性,具体表现为:下游区域侵蚀功率高于上游,干流区域强于支流,而西北地区则明显大于东南地区。
(3) 拉林河流域各子流域单元的径流侵蚀功率均值随其集水面积扩大呈幂函数衰减趋势,但当子流域面积达到150 km2这一阈值时,侵蚀功率对子流域集水面积的敏感程度显著降低并逐渐平稳。研究表明,子流域径流侵蚀功率与集水面积间存在显著的幂函数相关性。基于此,建议优先对集水面积小于150 km2的中上游子流域实施坡耕地改造、植被恢复等水土保持措施与生态治理工程,以实现对水土流失的有效控制。
(4) “侵蚀功率”不完全等同于“实际侵蚀量”,150 km2的阈值定义为一个“指示性”或“潜在”的临界范围,而非一个精确的、普适的管理红线,需要考虑空间异质性,因地制宜,但仍可以为类似于拉林河流域自然地貌以及气象环境的流域治理提供一定的参考价值。