土壤侵蚀是全球公认的最严重的环境问题之一，其破坏水土资源，加剧旱涝等自然灾害，威胁人类的生存和发展。嘉陵江流域地质地貌条件复杂、降水丰富、土壤可蚀性强，且不合理生产活动多，使得流域内水土流失面积和土壤侵蚀总量曾常居长江流域前列，水土流失治理刻不容缓^[1]。

当前有许多量化土壤侵蚀对流域影响的模型，如通用土壤流失方程、修正土壤流失方程和中国坡面土壤流失方程等^[2]。降雨侵蚀力反映降雨及其产生的径流剥离和携带土壤颗粒的能力，其作为水蚀动力指标被广泛应用于土壤侵蚀模型及土壤侵蚀分析中^[3-4]。但是降雨侵蚀力的计算通常需要高精度且不易获取的长序列场次降雨数据，资料处理较繁琐，且其仅通过雨滴的击溅效应表征土壤侵蚀作用，并未体现水蚀动力过程中的径流侵蚀和径流输沙作用^[5-6]，具有一定的局限性。与降雨侵蚀力相比，径流侵蚀功率能更好的反映水蚀动力特性，对于侵蚀动力机制的反应更敏感，数据要求更低^[7]。目前已有众多学者基于径流侵蚀功率理论研究黄土高原地区的土壤侵蚀空间分布特征及尺度效应，取得较好成果^[8-9]。长江流域的降水特征及下垫面条件与黄河流域不同，其侵蚀产沙和主要驱动机制也会有所差异^[10]，而此类研究在长江流域相对较少。因此，利用径流侵蚀功率理论研究长江流域地区的侵蚀情况对揭示其水沙响应机理具有重要意义。

笔者以嘉陵江子流域渠江流域为例，基于SWAT(Soil and Water Assessment Tool)模型模拟流域径流，计算并分析径流侵蚀功率的时空分布特征、空间尺度效应及主要影响因子，以期为流域有效识别重点产沙区、合理开发利用水土资源等做出贡献。

渠江流域面积为3.92万km²(E 106°28′~109°00′，N 30°00′~32°48′)，属亚热带湿润季风气候。其雨季集中于7—9月，多年平均降雨量为1 078 mm，输沙量主要来自汛期，多年平均输沙模数为347 t/km²。流域地势东北高西南低，源头地势坡度较大，至下游浅丘区比降逐渐减小，土壤类型以棕壤、黄棕壤和紫色土为主。渠江流域水系及水文站点分布见图 1。

图 1(Fig. 1)

图 1 渠江流域水系及水文站点分布图 Fig. 1 Distribution of river system and hydrological stations in Qujiang River Basin

笔者采用的数据主要包括摘自《长江流域水文年鉴》的2008—2018年渠江流域上、中、下游巴中、七里沱和罗渡溪水文站实测日流量资料，由中国大气同化驱动数据集整理的逐日降水、风速、温度等气象数据，以及中科院地理空间数据云、中科院数据中心遥感影像、世界土壤数据库提供的30 m分辨率地形数据、1∶10万土地覆盖数据和1∶100万土壤类型等空间数据。

笔者使用SWAT模型^[11]模拟渠江流域2009—2018年月径流，进而分析计算流域径流侵蚀功率。先依据集水面积、坡度等级、土地利用与土壤类型将渠江流域划分为77个子流域和1 377个水文响应单元。然后，选取2008年作为预热期进行参数预调，2009—2013年作为率定期，2014—2018年作为验证期。在SWAT-CUP中选择径流相关参数进行参数敏感性分析和参数率定，并运用决定系数R² (R-Square)、纳什效率系数NS(Nash-Sutcliffe efficiency coefficient)、偏差比例PBIAS(percent bias) 作为评价指标检验径流模拟精度。一般认为R²>0.6、NS>0.5、PBIAS≤±25%时，结果较好。满足精度要求时，输出渠江各子流域出口的模拟月径流过程。

径流侵蚀功率属于经验模型指标中的侵蚀动力因子，其计算原理及应用详见文献[7]。依此推出的季径流侵蚀功率和年径流侵蚀功率计算公式分别见式(1)和式(2)。

式中：E_s为季径流侵蚀功率，m⁴/(km²·s)；Q_s为季最大月流量模数，m³/(km²·s)；H_s为季径流深，m。

式中：E_y为年径流侵蚀功率，m⁴/(km²·s)；Q_y为年最大月流量模数，m³/(km²·s)；H_y为年径流深，m。

选用系统聚类方法^[12]和斯皮尔曼相关系数^[13]对径流侵蚀功率影响因子进行合并归类与相关性分析。斯皮尔曼相关系数的绝对值越大，指标之间相关性越高，其表达式见式(3)。

式中：ρ为斯皮尔曼相关系数；n为系列长度；d为按升序或降序排列后，同次序指标间的差值，量纲均为1。

率定期、验证期巴中、七里沱和罗渡溪水文站模拟与实测月径流量值对比见图 2~4。分析可知，巴中水文站率定期R²=0.96、NS=0.92、PBIAS=16.2%，验证期R²=0.70、NS=0.69、PBIAS=16.43%，率定结果较好，验证结果评级为可信；七里沱水文站率定期R²=0.95、NS=0.93、PBIAS=-0.3%，验证期R²=0.82、NS=0.80、PBIAS=-14.81%，率定结果极好，验证结果评级为较好；罗渡溪水文站率定期R²=0.93、NS=0.92、PBIAS=-3.5%，验证期R²=0.87、NS=0.84、PBIAS=-16.72%，率定结果极好，验证结果评级为较好。综上所述，在率定期和验证期内，3个水文站模拟所得月径流过程与实测月径流过程吻合度较高，表明所构建的SWAT模型能较真实地反映渠江流域水文情势。

图 2(Fig. 2)

图 2 巴中水文站月径流量模拟值与实测值对比图 Fig. 2 Comparison of simulated and measured monthly runoff at Bazhong station

图 3(Fig. 3)

图 3 七里沱水文站月径流量模拟值与实测值对比图 Fig. 3 Comparison of simulated and measured monthly runoff at Qilituo station

图 4(Fig. 4)

图 4 罗渡溪水文站月径流量模拟值与实测值对比图 Fig. 4 Comparison of simulated and measured monthly runoff at Luoduxi station

基于SWAT模型得到2009—2018年渠江77个子流域季尺度和年尺度的多年平均径流侵蚀功率，其空间分布见图 5和6。从时间角度看，全流域多年平均年径流侵蚀功率为0.049 m⁴/(km²·s)，总体上大于各季度多年平均径流侵蚀功率；在季尺度，第3季度多年平均径流侵蚀功率最大(0.028 m⁴/(km²·s))，其次为第2季度、第4季度和第1季度，其变化周期与水文周期一致。显然，第3季度的径流侵蚀功率对全年贡献最大，应当重视该季度的水土流失治理。在空间尺度上，流域全年和第3季度的多年平均径流侵蚀功率均呈现出“北部大、南部小；西部大、东部小；上游大、下游小”的空间分布特征。第4季度的空间分布虽在南北和上下游部分与全年相似，但在东西部表现出“东部大、西部小”的特征。第1季度和第2季度的空间分布与全年相反，表现为“南部大、北部小；东部大、西部小；下游大、上游小”，但因其径流侵蚀功率较小，对年尺度空间分布影响也较小。

图 5(Fig. 5)

4个季度的季径流侵蚀功率量级不同，值越大侵蚀越严重，其总体范围为0.000 01≤Es≤0.080 00。下同。 Magnitude of seasonal runoff erosion power varies among four periods, with larger values indicating more severe erosion. The overall range is 0.000 01≤Es≤0.080 00. The same below. 图 5 渠江流域多年平均季径流侵蚀功率空间分布图 Fig. 5 Spatial distribution of multi-year average runoff erosion power at season scale in Qujiang River Basin

图 6(Fig. 6)

图 6 渠江流域多年平均年径流侵蚀功率空间分布图 Fig. 6 Spatial distribution of multi-year average runoff erosion power at year scale in Qujiang River Basin

渠江流域的径流侵蚀功率空间差异明显，可能存在空间尺度效应，即径流侵蚀功率随流域控制面积变化可能呈现出一定规律。首先利用SWAT模型计算2009—2018年渠江77个子流域出口断面及以上控制面积的全年和第3季度的多年平均径流侵蚀功率，结果见图 7。

图 7(Fig. 7)

图 7 渠江各子流域出口断面及以上控制面积的多年平均径流侵蚀功率空间分布图 Fig. 7 Spatial distribution of multi-year average runoff erosion power for control area at and above the outlet cross-section of each Qujiang River subbasin

由图 7可知，对全年和第3季度，渠江各子流域出口断面及以上控制区域的多年平均径流侵蚀功率总体呈现出“上游大，下游小；干流大，支流小”的分布特点。为量化径流侵蚀功率与流域控制面积的相关关系，选取侵蚀较严重的渠江干流和支流大通江为研究对象，将流域控制面积分别与全年和第3季度的多年平均径流侵蚀功率进行拟合分析。干流的拟合结果分别见式(4)、式(5)和图 8，大通江的拟合结果分别见式(6)、式(7)和图 9。

图 8(Fig. 8)

图 8 渠江干流多年平均年和第3季度径流侵蚀功率与流域控制面积关系图 Fig. 8 Fitting results of multi-year average runoff erosion power and watershed control area at year and season scale(season Ⅲ) in the main stream of Qujiang River

图 9(Fig. 9)

图 9 渠江支流大通江多年平均年和第3季度径流侵蚀功率与流域控制面积关系图 Fig. 9 Fitting results of multi-year average runoff erosion power and watershed control area at year and season scale(season Ⅲ) in Datong River

式中：E_my、E_my1分别为干流和大通江多年平均年径流侵蚀功率，m⁴/(km²·s)；E_ms、E_ms1分别为干流和大通江多年平均季径流侵蚀功率，m⁴/(km²·s)；A为流域控制面积，10³ km²。

由图 8和9可知，渠江干流和支流大通江存在较显著的空间尺度效应，其多年平均径流侵蚀功率与流域控制面积之间均呈幂指数关系，且干流相关关系优于支流大通江，全年相关系数与第3季度接近。随着流域控制面积的增大，多年平均径流侵蚀功率逐渐减小并趋于平缓，减小速率由大变小，说明存在空间阈值。为确定阈值，对式(4)~(7)求导，导数方程见式(8)~(11)。

经计算可知，渠江干流在全年和第3季度的空间阈值分别为8 549.4和4 834.9 km²，大通江在全年和第3季度的空间阈值分别为8 504.4和6 223.5 km²。即当流域控制面积小于空间阈值时，|E|≥0.001，多年平均径流侵蚀功率随流域控制面积增加而迅速减小，反之则变化平缓，并逐渐减小至某一稳定值。其中干流的全年和第3季度径流侵蚀功率稳定值分别为0.04和0.02 m⁴/(km²·s)，大通江的全年和第3季度径流侵蚀功率稳定值分别为0.056和0.035 m⁴/(km²·s)。

为探究渠江流域径流侵蚀功率时空分布及空间尺度效应的主要驱动因素，选用SPSS(statistical product and service solutions)软件对多年平均的年径流侵蚀功率进行聚类分析。本文共选取地形(75°以上坡度面积比例)、流域形态(圆度)、气象(各子流域上的多年平均降水量)、侵蚀动力(多年平均径流侵蚀功率)、土地利用(林草地面积比例)和土壤(易受侵蚀土壤面积比例)6类聚类因子，经分析将渠江77个子流域划分为3种聚类类型，其空间分布见图 10。可知，第1聚类与第2聚类主要分布于上游地区，其面积所例分别为26.0%和19.3%，多年平均径流侵蚀功率分别为0.058和0.052 m⁴/(km²·s)，均大于多年平均径流侵蚀功率中值0.049 m⁴/(km²·s)，该地区径流侵蚀能量较大，受侵蚀程度较剧烈。第3聚类是渠江流域的主导聚类，主要分布在中下游地区，其面积比例为54.7%，流域内多年平均径流侵蚀功率为0.046 m⁴/(km²·s)，小于多年平均径流侵蚀功率中值，该地区植被类型多样，下垫面条件较好，受侵蚀程度较小。

图 10(Fig. 10)

图 10 渠江流域聚类空间分布图 Fig. 10 Clustering spatial distribution of Qujiang River Basin

统计侵蚀动力因子与其余聚类因子的斯皮尔曼相关系数见表 1，结果表明气象因子、地形因子和流域形态因子为渠江流域径流侵蚀功率排名前3的主要影响因子。且在显著水平为0.05时，第1聚类中的气象因子与径流侵蚀功率呈显著相关，说明降雨的空间分布不均性对该地区径流侵蚀能量影响较大；第2聚类中气象因子、地形因子和流域形态因子与径流侵蚀功率呈显著相关，说明该地区75°以上坡度面积比例很大，流域内地势陡峭，降雨空间分布差异较显著，流域形态偏圆形，产汇流过程较快，能显著影响径流侵蚀功率的空间分布和规律；第3聚类中流域形态因子与气象因子与径流侵蚀功率呈一般相关，其对该聚类地区径流侵蚀功率的影响远不如第2聚类地区，说明该地区整体地势较平缓，各子流域多年平均降雨量差距不大，流域形态较狭长，产汇流过程相对较慢。

表 1(Tab. 1)

表 1 斯皮尔曼相关系数计算结果表 Tab. 1 Spearman correlation coefficient calculated results table 聚类名 Cluster name 地形因子 Topographical factor 流域形态因子 Watershed morphology 气象因子 Meteorology 土地利用因子 Land use 土壤因子 Soil 第1聚类Cluster Ⅰ -0.451 -0.484 0.621* -0.280 -0.022 第2聚类Cluster Ⅱ 0.636* 0.622* 0.678* 0.441 -0.336 第3聚类Cluster Ⅲ 0.068 0.349* 0.336* -0.261 0.072 注：*表示在0.05水平上相关性显著。Notes: * indicates significant correlation at the 0.05 level.

表 1 斯皮尔曼相关系数计算结果表 Tab. 1 Spearman correlation coefficient calculated results table

1) 所构建的SWAT模型径流模拟精度较高，在渠江流域适用性好。巴中、七里沱、罗渡溪水文站率定期R²、NS系数均在0.92及以上，PBIAS系数均＜16.20%；验证期R²、NS系数均在0.69及以上，PBIAS系数均在16.72%以下。这说明该模型能较真实地反映渠江流域水文情势。

2) 渠江流域多年平均年径流侵蚀功率总体上大于季径流侵蚀功率，其中第3季度>第2季度>第4季度>第1季度。流域全年和第3季度的多年平均径流侵蚀功率空间分布均呈现出“北部大、南部小；西部大、东部小；上游大、下游小”的特征；第4季度的空间分布在东西部呈现出“东部大、西部小”的特征，其余部分与全年相似；第1、第2季度的空间分布与全年分布规律相反。

3) 渠江干流及其支流大通江的径流侵蚀功率具有较为显著的空间尺度效应，全年和第3季度的多年平均径流侵蚀功率与流域控制面积之间呈幂指数均关系。当流域控制面积分别＞8 549.4和4 834.9 km²时，干流全年和第3季度的多年平均径流侵蚀功率随着面积增加变化幅度极小，并逐渐稳定于0.04和0.02 m⁴/(km²·s)；当流域控制面积分别大于8 504.4和6 223.5 km²时，大通江全年和第3季度的多年平均径流侵蚀功率随着面积增加变化幅度极小，并逐渐稳定于0.056和0.035 m⁴/(km²·s)。

4) 气象因子、地形因子和流域形态因子对渠江流域径流侵蚀功率的分布有主要影响作用。流域上游区因地势陡峭、降水分布不均、形态易于产汇流而表现出较大的径流侵蚀功率；下游地区因地势平坦，产汇流过程缓慢受侵蚀情况较轻。