黑土小流域沟道分布遥感监测及主控因素研究

温艳茹，余强毅，杨扬，张斌，吴文斌*

（1.中国农业科学院农业资源与农业区划研究所/农业农村部农业遥感重点实验室，北京 100081；2.北京师范大学地表过程与资源生态国家重点实验室，北京 100875）

我国黑土区生产全国四分之一的粮食和三分之一的商品粮，对保障粮食安全至关重要。我国黑土区现有切沟近30万条，且61%出现在耕地，损毁耕地占0.5%～3%[1]。切沟侵蚀使耕层变薄、地力下降，直接威胁黑土可持续发展[2-3]，同时，侵蚀土壤在下游堆积会加剧洪涝灾害和环境污染，且侵蚀也会导致碳排放增加，不利于国家“碳中和”目标的实现[4-5]。然而，黑土区地形漫川漫岗、夏季降雨集中、冬春冻融交替、起垄耕作、土壤黏重，加之长期高强度利用，导致切沟分布复杂，主控因素多变，其他区域（如黄土高原）沟蚀理论和防治措施无法在此合理应用[6-7]。因此，迫切需要系统开展黑土区切沟分布及其主控因素研究，为明晰黑土沟蚀发育过程和针对性阻控沟蚀提供科学支撑。

切沟侵蚀是在水力、风力、冻融等作用下，地表径流集中冲刷形成的永久性沟道[8-9]，与浅沟相比，切沟不能被常规耕作逾越。水力主导下沟头沿着坡面汇流通路向坡上溯源前进增加切沟的长度，沟壁崩塌向两侧扩张增加切沟的宽度，沟底水力下切增加切沟的深度，使得切沟从坡上向坡下垂直于等高线分布[10]。切沟分布特征反映了降水侵蚀力和汇水面的影响，随着新的汇流通路形成，产生二级甚至三级切沟，形成复杂的切沟系统[11]。除了水力外，东北黑土切沟形成还受冻融、特殊的多面向地形以及土壤耕作等因素的影响，切沟形态呈复杂的“S”形态[12]。随着3S 技术的发展，切沟的监测和评价手段也由传统的卷尺、侵蚀基准针和断面测定仪等，陆续发展到GPS、高分遥感、摄影测量、机载/地面激光扫描等新技术[8-14]。其中，遥感影像监测，尤其是无人机摄影监测因具有测量快、时空分辨率高和精度高等显著优势，广泛用于监测切沟形态、时空分布以及构建三维数字高程模型。

地形是影响切沟发生发展的重要因素，可通过汇水面积、坡长、坡度等因子影响下垫面形态，从而影响径流冲刷过程。88%的黑土区为小于5°的长缓坡，坡长可达几百米甚至上千米，汇流面积较大，在降雨作用下地表易形成集中汇流，冲刷发育成切沟[8-10]。有研究表明切沟密度随着汇水面积和坡长的增加呈现先增加后减少的趋势，临界面积约29 hm2[15]。而坡度主要影响地表土壤与雨滴的接触角、地表径流流速及径流对土壤的剪切力，侵蚀速率随坡度的增加先增加后减少。WANG 等[16]的研究指出，黑土坡面坡度为3°～5°时，坡面沟壑密度最大。此外，切沟密度也表现出随坡长增加呈先增加后降低的趋势，汇水区地形不同，临界坡长有一定的差异。许晓鸿等[17]在吉林梅河口的研究发现切沟密度在200～299 m 最大，而顾广贺等[18]指出沟密度的临界坡长为900 m。一方面，黑土区小流域切沟侵蚀发育与地形因子间的定量关系表现出一定的空间变异，随着地形条件不同，定量结果并不太一致；另一方面，不同汇水条件下地形因子，如汇水面积、坡长、坡度，对切沟发育的相对重要性尚不明确。因此有必要综合分析流域尺度上切沟侵蚀的主控因素，这对于开展小流域侵蚀沟综合治理尤其是有针对性地阻控切沟发育十分关键。

本研究选取黑龙江省海伦市光荣村小流域264条切沟为研究对象，采用遥感影像和地形图目视解译人工判读，结合野外详查、地貌学分析和统计分析的方法，量化小流域尺度上的汇水面积、坡长、坡度高差等变量，探讨机耕道影响的地表径流对切沟分布的贡献，揭示黑土切沟发育的主控因素，为深入理解沟蚀过程、科学监测和有效阻控黑土耕地切沟发育提供关键科学支撑。

1 材料与方法

1.1 研究区概况

本研究选取典型黑土区中心区域——黑龙江省海伦市光荣村小流域（47°23′N，126°51′E）（图1），面积约为24 km2。研究区属北温带大陆性季风气候，夏季高温多雨，冬季寒冷干燥。中国科学院海伦农田生态系统国家野外科学观测研究站1965—2010 年气象数据显示，年均降水量为543 mm，降雨集中于6—9月。年均气温为1.5 ℃，平均最低气温（-21.8 ℃）和最高气温（21.7 ℃）分别在1月和7月。研究区为黑土区典型漫川漫岗地貌，坡缓坡长（坡度＜7°，坡长200～1 000 m）。土壤类型为典型黑土，土壤质地为黏壤土，其表层（0～20 cm）黏粒含量为35%～45%，有机质含量为42.1 g·kg-1[19]，母质为黄土状亚黏土，泥页岩为第四纪河湖相沉积物[20]。

图1 研究区地理位置Figure 1 The geographical location of the study area

研究区域原生植被主要为草地和零散分布的落叶林。19 世纪90 年代后期坡上开始零散开垦耕地，为更好获取光照，采取南北向顺坡垄，这在一定程度上加速了土壤侵蚀。新中国成立后，快速增长的人口对耕地需求增加，至1955 年，研究区绝大部分坡上已被开垦为耕地，到1982 年前后所有原始林地被开荒殆尽[21]。研究表明，黑土侵蚀沟主要形成于大规模开垦耕地的20世纪五六十年代，1968—2018年，研究区侵蚀沟密度增加1 倍（1.2～2.3 km·km-2），侵蚀沟规模增大且更加复杂[7]。

1.2 数据与处理

本研究基于2018 年4 月卫星影像（地面分辨率0.3 m）和1982 年等高线地形图（比例尺1∶10 000，等高距1 m）。等高线地形图数据来自黑龙江省测绘地理信息局，在前期使用时已完成拼接和裁剪[21]。卫星影像数据来源于谷歌Digital Globe 公司，已完成正射校正。通过高精度控制点对两期影像进行几何纠正和地理配准，在ArcGIS 环境中转换成WGS_1984_UTM 投影。根据野外实测结果更新检验解译结果。

1.3 研究方法

本研究采用遥感目视解译、地貌学野外详查结合统计分析的方法，量化黑土区小流域地形地貌特征、高差、坡度、地表径流和土地利用等对沟道侵蚀发育的影响。

1.3.1 目视解译

本研究基于2018 年遥感影像，通过目视解译获取2018 年沟道、机耕道等线状分布特征。主要沿沟底中心线绘制，复杂沟道系统中单条沟道独立计数，并单独测量长度。切沟的辨认标准为遥感影像上有线状斑驳、颜色表现为略白，伴随沟底处有沉积区、颜色表现为偏深或黑色、形状多扇形。研究区土地利用主要以耕地、林地和村级居民点为主，其中耕地约占研究区总面积85%以上。本研究中土地利用类型主要提取林地，基于2018 年遥感影像目视解译、人工判读和绘制林地轮廓，在ArcGIS 中计算各汇水区的林地面积。

1.3.2 地貌学分析

基于1982年等高线地形图在ArcGIS中划分基本汇水区范围，沿等高线弯曲方向，绘制出各个汇水区之间的分水线。通过遥感影像解译、地形图分析、野外详查，获取小流域特征、沟道特征和驱动因素，量化流域尺度贡献沟道发育的主控因素。每个汇水区的变量组包括，①侵蚀沟道形态组：总沟长、主沟长、沟面积和沟密度；②汇水区特征组：汇水区面积、周长和最长距离；③高差和坡度组（地貌特征）：汇水区高差、主沟高差、汇水区平均坡度和主沟平均坡度；④土地利用组：汇水区林地面积，其中切沟特征为y变量，数据收集和处理方法见表1。

表1 研究区内汇水区特征指标及其收集方法Table 1 Sub-catchment properties defined based on digital elevation model and ArcGIS

本研究结合目视解译和等高线图地形地貌分析，获取侵蚀沟道分布、发育过程及其驱动因素的关系，把等高线与侵蚀沟道交汇情况归类（图2）。①等高线向高值方向凸出，地形上体现为凹形坡，集中汇聚水流对应深度5 m 以上的侵蚀沟；②等高线向低值方向凸出，地形上体现为凸形坡，一般是高地，水流分散。同时，根据野外详查经验，将影响沟道发育因素分为6 类：①管流，主要地下径流汇聚；②地表径流主要驱动，明显汇聚地形下的侵蚀沟，归为地表径流；③地表径流主要驱动，结合机耕道影响，如机耕道位于沟道汇水线附近，本研究认为其有可能对地表径流驱动的沟道发育有影响，则归为地表径流（机耕道）；④机耕道明显引起沟道发育，沟头连接机耕道，经过野外确认，归为机耕道；⑤机耕道作为主要原因引起地表径流，如平行坡情况，沟道发育不受汇聚水流影响，但明显受到机耕道的影响，或者沟头连接机耕道，归为机耕道（地表径流）；⑥不确定因素，如不是以上5种情况，归为不确定的其他因素。最后，将特征信息加入ArcGIS属性表中进行下一步分析。

图2 侵蚀沟道发育方向与等高线的曲面情况Figure 2 Classification of gully sections and main gully using the contour line

1.3.3 野外详查

结合2017 年6 月、8 月的典型侵蚀沟道野外调查结果，于2018年4月10—18日对研究区264条侵蚀沟道、机耕道、林地等进行实地确认和野外详查（图3）。因调查时积雪已融化，春播尚未开始，无农作物和积雪覆盖影响，因此详查以地形图解析的小流域单元展开，主要内容包括沟道是否存在、位置、实际数量、长度、活跃类型等。具体流程：①携带彩色打印的卫星影像、侵蚀沟道分布图、谷歌底图、地形图等逐条走访，完成“侵蚀沟道调查数据表”；②携带手持GPS 从沟头向沟尾测量整条沟长，并在沟头、沟尾处定点测量，完善数据表；③区分活跃期和非活跃期侵蚀沟道，依据沟头、沟壁和沟底植被覆盖程度，沟头是否有裂隙和崩塌，沟头上方是否有汇水痕迹，及沟底是否有短期流水（积雪融水，降水）痕迹，满足四条中的两条，即可判断为活跃沟道[22]；④重点调查林地内侵蚀沟道是否存在、活跃程度和长度等。

图3 黑土典型切沟侵蚀Figure 3 Typical gullies in the Mollisol region

1.3.4 数据统计分析

采用Excel 2010进行数据处理，用SPSS 25.0进行Pearson 相关性分析、多元线性回归分析和方差分解分析。用Origin 2020b绘制侵蚀沟道与汇水区特征相关性矩阵图、等高线与切沟驱动因素的相关性热图。利用R studio 和Inkscape 绘制切沟与汇水区特征方差分析图。

2 结果与分析

2.1 汇水区特征对切沟分布的影响

研究区域总面积23.73 km2，2018 年林地面积217.1 hm2，占研究区面积的9.1%。区域共有切沟264条，切沟总长度为52.1 km，其中活跃切沟的长度为41.8 km，占总长度的80.2%。林地中活跃沟长约19.3 km，占活跃切沟总长度的46.2%。研究区分为19 个汇水区，具体侵蚀沟道形态和汇水区特征见表2。不同汇水区特征差异较大，平均面积为1.25 km2，变化范围为0.06～4.33 km2；平均最大坡长为1.43 km，平均坡度为3.6°，平均高差为44.0 m。汇水区内的切沟平均总长度2.74 km，平均主沟长0.96 km，平均沟面积3.06 hm2，平均沟密度2.69 km·km-2，变化范围为0.99～6.68 km·km-2。

表2 研究区19个汇水区及切沟特征Table 2 Properties of the sub-catchment and the total gullies in the study area

对19 个汇水区的切沟特征与汇水区地形、高差、坡度等进行相关性分析（图4），结果表明，切沟总沟长、主沟长、沟面积与汇水区面积、周长、最长距离均呈显著正相关，即汇水面积、周长和坡长越大，切沟长度和面积越大。切沟特征与汇水区高差无显著相关性。总沟长、主沟长、沟面积与汇水区坡度、林地面积占比均呈显著负相关。汇水区坡度在1.3°～7.9°之间变化，坡度越小，切沟长度和面积越大，缓坡促进汇水，使得切沟规模较大。此外，沟密度与主沟高差、坡度呈显著正相关，切沟主沟高差在9.4～44.9 m 之间变化，主沟平均坡度在0.7°～11.2°之间变化，主沟高差和坡度越大，切沟越密集。

图4 切沟侵蚀特征与汇水区特征的相关性矩阵Figure 4 Correlation matrix between the gully networks and sub-catchment properties

2.2 汇水区地形特征对切沟特征的影响

通过多元线性回归分析方法，获取汇水区内的切沟总沟长、主沟长、沟面积、沟密度与汇水区相关显著变量之间的回归方程（表3），探究影响切沟特征的主要原因。结果表明，总沟长、主沟长、沟面积与汇水区特征密切相关（R2分别为0.928、0.965 和0.929，P＜0.001），沟密度主要与汇水区高差、坡度相关（R2为0.443，P＜0.05）。其中，汇水区最长距离对总沟长、主沟长和沟面积具有显著的正向影响。最长距离为分水岭最远点和汇水区的出水口点的距离，即最大坡长，平均值为1.43 km，变化范围为0.30～3.19 km。

表3 多元线性回归方程Table 3 Multiple linear regression equation

通过对汇水区特征进行方差分解分析，获取各变量组对切沟特征的贡献率（图5）。汇水区特征变量组分别解释了总沟长、主沟长和沟面积的88.2%、95.0%和88.3%，模拟效果较好。其中，汇水区特征组（面积、周长和最长距离）对总沟长和主沟长的贡献率分别为32.9%和42.0%，即汇水面积、坡长对切沟长度的影响显著。土地利用（林地面积占比）对总沟长、主沟长和沟面积的贡献率较小，分别为0.2%、2.0%和0.6%。

图5 汇水区变量组对总沟长（a）、主沟长（b）、沟面积（c）的贡献率Figure 5 Contribution rates of variable category of sub-catchment on total length of gullies（a），length of main gully（b）and gully area（c）

2.3 地表径流、机耕道对切沟分布的影响

结合切沟分布与等高线交汇划分的3 类情况（汇聚的凹形坡、平行坡和分散的凸形坡）与地表径流的6 种情况，获得切沟侵蚀驱动因素热图（图6），再结合研究区机耕道和切沟的空间分布图（图7），结果表明，发育在汇聚凹形坡的切沟长度占总长度的68.1%（深度＜5 m 的坡谷30.7%，深度＞5 m 的谷地37.4%），平行坡占30.3%，分散地形占1.6%。地表径流驱动是导致研究区内切沟发育的最主要因素（长度占总长的63.0%）。除地表径流外，机耕道为目前亟需关注的驱动因素，黑土切沟多伴生于农田道路一侧。研究区与机耕道有关的切沟长度占总长度的23.3%，其中位于农田道路旁切沟有11.8%，可能受机耕道影响的切沟有11.5%。在凹形坡谷（深度＜5 m），地表径流是最显著的影响因素（20.4%），机耕道只占5.0%；在凹形谷地（深度＞5 m），除地表径流（27.5%）外，机耕道占9.8%。此外，在平行坡面上，地表径流的影响占15.1%，机耕道的影响占7.4%。

图6 等高线的曲面情况与切沟侵蚀驱动因素热图Figure 6 Heatmap of orientation,curvature types and driving factors in total gully length

图7 2018年研究区切沟与道路的空间分布图Figure 7 Gully and tractor road distribution in the study area in 2018

3 讨论

目前东北黑土区切沟侵蚀唯一权威数据为全国第一次水利普查的侵蚀沟专项调查公告（2013 年），基于2.5 m 分辨率遥感影像和野外核查的结果显示，东北黑土区有大于100 m 的侵蚀沟29 万余条。有研究指出实际数量比公告数量高1 倍，即实际侵蚀沟数量约60 万条[2]。究其原因，一是遥感解译难以将林下或林边沟道提取出来；二是2.5 m分辨率的遥感影像难以将伴随道路的切沟全部解译[2]。本研究目视解译结合野外详查结果表明，林下活跃切沟占切沟总长的46.2%，机耕道旁切沟占总沟长的10%左右。

研究区林地面积占总面积9.1%，林地中活跃切沟长度占活跃切沟总长的46.2%，同时，林地面积与沟长度、沟密度和沟面积呈显著负相关，但贡献率小于2%。与研究区切沟长期发育特征相比，1968—2018年切沟密度增加了1倍，林地面积增加7%[7]。然而，光荣村农田道路变化是切沟密度提高的最主要原因。大量研究表明，退耕还林还草通过增加入渗、减少地表径流、增加土壤抗蚀能力等作用防治土壤侵蚀，且在一定程度上稳定了活跃切沟[23]。然而，增加造林还草面积并不能从根本上解决切沟侵蚀问题。一方面，研究区林地主要分布在沿海伦河的坡度陡降区域，该区域在20 世纪60 年代和90 年代修建农业梯田，退耕还林还草。已有研究表明海伦市光荣村可见废置的农业梯田或林地下有活跃切沟发育[12]。农业梯田废置后疏于定期检查与维护，导致梯田面有倾斜坡度或梯田梗缺失，进而改变汇水线或增大汇水面积，这是林地下活跃切沟发育的主要原因。另一方面，小流域地形起伏，加之切沟形态复杂，造林时很少依据流域景观格局进行科学设计，使植树初期地表覆盖率低，不能很好地截留水源。黑土侵蚀沟综合治理仍然缺乏针对切沟的长期定位研究数据。

本研究的汇水区平均面积为1.25 km2，平均最大坡长为1.43 km，平均坡度为3.6°，属于典型黑土漫川漫岗坡耕地地貌。切沟长度主要受到汇水面积和坡长的影响，其中坡长影响更大，而切沟密度主要受到主沟高差和坡度的影响（图8）。一般来说，汇水区面积和坡度越长，区内切沟的总长度越大。以往研究表明流域/汇水区面积是影响切沟规模和发育规律的重要因素[8]。地形地貌特征是地表径流过程的主要影响因素，切沟发生发育主要由坡面的水文动力过程控制[9-12]。本研究发现，对于3 hm2左右的汇水区，汇流面积和长缓坡仍然是未来侵蚀沟治理需要考虑的主要因素。因此，小流域综合治理侵蚀沟必须考虑拦截长坡，减少长坡导致的地表径流汇聚，特别是配套排水系统的设计和建设。

图8 黑土区小流域沟道形成和发育示意图Figure 8 Gully initiation and development in the catchment scale of Mollisol region

研究区内在汇聚凹形坡下受地表径流主要影响的沟长度占总长的63.0%，受农田机耕道影响的切沟长度占23.3%。丁超等[24]对东北典型黑土区1 049 条切沟进行遥感解译和实地详查，结果显示耕地切沟占55%，路边切沟占34%，林地切沟占10%。研究区的农田机耕道多位于地块边缘，多为土路。农用拖拉机等碾压后形成车辙，成为暴雨后的水流通道。同时，农田道路地表被机械压实，研究表明重力机械的压实土壤容重可达1.4 Mg·m-3，孔隙度大幅降低[25-26]。加之黑土黏重，黏粒含量可达40%，导致入渗速率低至2.5×10-5m·s-1，直接阻碍春季融雪径流和夏季降雨径流入渗[17，27-28]。农田道路硬化后汇集地表径流，加速沟头形成并溯源上行，沟道向下深切，沟壁崩塌增宽。

本研究以东北黑土区海伦市光荣村内264 条切沟为研究对象，基于高分辨率遥感影像和数字高程模型，结合实地详查，分析切沟分布、长度、密度等与地形特征、机耕道等的关系。黑土区切沟长度约100～1 500 m，宽度约10～20 m，深度约5～15 m，受云雾、种植作物以及树木的影响，基于高分辨率遥感影像绘制切沟存在一定的误差。此外，本研究基于遥感目视解译结合地貌学方法，将切沟分布与等高线的关系划分为3 类：汇聚的凹形坡（切沟深度＞5 m 和＜5 m）、平行坡和分散的凸形坡。同时，结合地形和土地利用类型，人工判读影响切沟发育的6 类因素，包括汇聚地形下地表径流汇集以及机耕道等。基于前人研究结果，为保证结果的准确性，切沟的判读和绘制需要熟悉研究区域和野外切沟调查经验丰富的人员进行，同时，必须结合野外详细的逐条切沟测量，进而确认切沟位置和长度等特征[7，12，29]。

本研究结合以往研究成果，提出以下建议：

（1）将农田道路周边大型活跃切沟防控作为黑土保护的重中之重。黑土地区坡耕地土壤侵蚀和大型活跃切沟必须协同治理，才能更好地发挥作用。切沟预防的关键措施在于建设农田道路配套排水管道，为此，需要加快制定黑土区侵蚀沟治理技术标准，特别是黑土坡耕地农田道路合理布局和强制性配套排水系统的建设规范。

（2）增加林木种植密度、控制坡面坝埂截水面积可增强水分就地入渗，减少因地表局部汇流而形成切沟。造林保持水土应该考虑山水林田湖草综合因素，结合小流域设计进行综合治理，合理布局林草地、耕地和机耕道。同时，建议增修垂直于梯田埂的截留土坝，以拦截径流，增加入渗；对于已经还林的农业梯田，应该给予更多关注，重点是要定期检查和维护。

（3）把黑土区切沟发生高风险区域作为黑土地保护工程重点区域。根据地形和地质条件确定切沟发生高风险区域，优先纳入黑土地保护工程实施范围。在黑土区切沟发生高风险区域加大农田道路排水系统配套建设和高标准侵蚀沟整治试验示范力度，为大面积推广侵蚀沟治理技术打好基础。

4 结论

（1）黑土区海伦市小流域汇聚地形下，水力主导切沟沟头沿着坡面汇流向坡上溯源增加切沟长度，地表径流直接影响的切沟长度占总长的63.0%。汇水区面积和坡长是影响切沟长度的最主要因素，其中坡长影响最大，而主沟坡度和高差是影响沟密度的显著因素。

（2）黑土区林地下活跃切沟不容忽视，2018 年研究区林地面积217.1 hm2，占研究区总面积的9.1%，林地活跃沟长达19.3 km，占活跃切沟总长的46.2%。对于坡度较陡的汇水区，林木多分布在切沟两侧，密度不足，垄作导流，使得林地内也会出现地表径流汇集，形成切沟侵蚀的退化地貌。

（3）目前，机耕道是黑土切沟不可忽视的影响因素。黑土切沟多伴生于农田道路一侧，研究区农田道路旁切沟长度占其总长的23.3%。农田道路的地表硬化后汇集地表径流，加速沟头形成并溯源上行，沟道向下深切，沟壁崩塌增宽。切沟侵蚀导致耕地破碎，农田道路切断或被迫移位，不仅损毁耕地，而且阻碍现代农业的发展。侵蚀沟综合治理必须考虑小流域景观布局和山水林田湖草综合治理，合理布设农田道路。