基于WorldView-2影像估算黄土高原六道沟小流域退耕还林(草)后的切沟发育速率

黄婷婷, 杨 扬†, 史扬子, 曹 琦, 冯 博, 刘宝元, 刘瑛娜

(1. 北京师范大学地理科学学部 地表过程与资源生态国家重点实验室, 100875, 北京;2. 北京师范大学地理科学学部地理学院, 100875, 北京; 3. 中国科学院 水利部 水土保持研究所, 712100, 陕西杨凌; 4. 沈阳农业大学水利学院 辽宁省水土流失防控与生态修复重点实验室, 110866, 沈阳)

切沟是黄土高原主要的土壤侵蚀类型之一。在黄土丘陵沟壑区，切沟侵蚀的产沙量通常占小流域总产沙量的50%以上，部分地区甚至可以达到70%以上[1-2]。21世纪以来，GPS技术发展迅速，其中实时动态差分GPS(real-time kinematic GPS, RTK GPS)因其定位精度高、作业速度快等优势，成为切沟监测的重要手段之一[3]；然而，RTK GPS测量成本高、野外负担重，多用于较小时空尺度的监测研究。近年来，高分辨率遥感产品不断涌现，兼具监测时间长、覆盖范围大等特点，为研究较大时空尺度的切沟分布特征与发育速率提供便捷可靠的数据源[4-5]。利用高精度野外测量仪器如三维激光扫描仪[6]、RTK GPS[7]评估遥感影像的切沟解译精度，分析误差来源并总结解译经验，对于提高遥感影像解译精度具有重要意义。

切沟发育与集水区面积、降雨强度、土地利用和土壤物理性质等密切相关[8-9]。黄土高原切沟侵蚀严重，退耕还林(草)工程大规模实施后，土地利用格局改变，必然导致切沟侵蚀的相应变化。笔者以六道沟小流域为研究对象，利用RTK GPS进行切沟实地测量，评价WorldView-2影像的切沟解译精度；在此基础上，对比2012年和2018年WorldView-2影像的切沟解译结果，估算退耕后的切沟发育速率并分析其影响因素，以期为治沟治坡决策、水土保持措施配置提供科学依据。

1 研究区概况

六道沟小流域(E 110°21′～110°23′，N 38°46′～38°51′)位于陕西省神木市以西14 km(图1)，属黄土高原北部水蚀风蚀交错区，是典型的盖沙黄土丘陵沟壑地貌。小流域面积689.2 hm2，主沟道长4.2 km，自南向北流入窟野河的一级支流三道河。小流域属中温带半干旱气候，夏秋季节暴雨集中，冬春季节干旱多风。自2002年退耕还林(草)工程在神木市全面实施以来，小流域内大量坡耕地退耕；同时，林草生态恢复与人工造林种草等生物措施广泛开展，植被盖度增加[10]。

图1 六道沟小流域与RTK GPS所测切沟的地理位置Fig.1 Geographic locations of Liudaogou catchment and gullies measured by RTK GPS

经实地踏查，笔者在六道沟小流域共选取33条切沟开展RTK GPS测量(图1b)。2018年10月，在道路两侧及步行易达的山坡测量10条代表性切沟，切沟长度介于8～71 m之间，面积范围25～613 m2，主要土地利用类型为草地。2019年5月，在这10条切沟的基础上，根据沟缘植被、沟壁特征等另外选取23条切沟，长度范围扩大至7～200 m，面积20～2 500 m2。

2 材料与方法

2.1 WorldView-2影像的切沟解译验证

2.1.1 RTK GPS测量 采用水平测量精度<1 cm、垂直测量精度<2 cm的中海达iRTK2经典版差分GPS对所选切沟开展实地测量。按大致0.5 m的间距测量切沟沟缘线各点的地理坐标，同时在沟缘曲折变化剧烈处适当增加测量点。将所测各点坐标导入ArcGIS 10.2，分别提取各切沟的长度、最大表面宽度、周长和面积。

2.1.2 WorldView-2影像处理与解译 以2018年10月23日的WorldView-2(WV2)影像为数据源，对0.5 m分辨率的全色数据和2.4 m分辨率的多光谱数据进行正射校正后，利用Gram-schmidt方法融合校正后的全色和多光谱数据，得到0.5 m分辨率的真彩色影像。在此基础上，根据形状、尺寸、阴影、色调和纹理等对影像进行切沟目视解译。

2.1.3 WorldView-2影像解译精度评估 为评估WV2解译各切沟参数p的精度，计算其相对误差

(1)

式中：pwv和pgps分别为基于WV2和RTK GPS得到的切沟参数p。

在此基础上，采用配对制图误差指数E分析WV2和RTK GPS描绘切沟轮廓的差异[11]：

(2)

式中：Awv∪Agps为基于WV2和RTK GPS得到的切沟并集面积，hm2，Awv∩Agps为相应的切沟交集，hm2；E取值范围为0～100%，0表示二者轮廓完全吻合，100%表示完全不吻合。

2.2 小流域集水区划分及切沟解译

2.2.1 集水区划分及其属性提取 基于1∶1万地形图生成5 m分辨率的DEM(digital elevation model)，将六道沟小流域划分为若干集水区，并参照Strahler[12]河网分级原则进行集水区分级，即单条干沟控制的集水区为Ⅰ级集水区；2个Ⅰ级集水区交汇形成Ⅱ级集水区；不同级别的集水区交汇，等级为原来集水区中级别较高者。

利用DEM提取各集水区的平均坡度、坡长、坡向、剖面曲率和平面曲率作为背景因子。因2012年无可用Landsat数据(http://www.gscloud.cn)，以2013—2018年夏季30 m 分辨率的Landsat 8 OLI影像为数据源提取六道沟小流域的植被盖度,计算小流域多年平均植被盖度后，再根据各集水区范围得到其对应植被盖度。

参照小流域2010年和2017年的遥感影像和土地利用数据[10]，建立遥感解译标志，分别对2012年和2018年的WV2影像进行判读，获得相应的土地利用数据。由此计算各集水区的综合土地利用动态度[13]，以量化集水区d的土地利用在2012—2018年间(Δt)的总体变化速率：

(3)

式中：Kd为集水区d的综合土地利用动态度，%/a；Δt为研究期时间，a；ΔUij为集水区d在Δt内土地利用类型i转变为土地利用类型j的面积，hm2；n为集水区d的土地利用类型总数；Sd为集水区d的面积，hm2。

2.2.2 集水区切沟解译 以2012年11月23日与2018年10月23日的WV2遥感影像为数据源，分别解译各集水区内的切沟，影像处理及解译方法同2.1.2。在此基础上，统计各集水区各切沟的长度、最大表面宽度、周长和面积，分别加和得到对应集水区的切沟总长度、总最大表面宽度、总周长和总面积。

2012—2018年间(Δt)集水区d的参数p年均变化率

(4)

式中：Δd,p为集水区d在Δt内的参数p年均变化率，%/a；pd,wv2018和pd,wv2012分别为2018年和2012年得到的集水区d参数p的总和。

3 结果与分析

3.1 WorldView-2影像的切沟解译精度评价

3.1.1 WorldView-2与RTK GPS所得切沟参数对比 相比2018年10月RTK GPS测量的10条切沟数据，同时期WV2解译得到的切沟长度、最大表面宽度、周长和面积相对误差分别介于0.02%～3.2%、0.1%～9.2%、1.6%～11.7%和0.1%～9.3%之间，平均值分别为1.3%、3.8%、5.4%和3.5%。配对制图误差指数介于7.4%～24.4%之间，平均15.9%，小于Fiorucci等[14]以RTK GPS为基准得到的18%的WV2制图误差指数，表明同时相WV2影像与RTK GPS得到的切沟轮廓非常接近，可用于切沟解译。

对比2018年10月和2019年5月利用RTK GPS测得的10条切沟数据发现，2个时期切沟变化不大，切沟长度、最大表面宽度、周长和面积变化率分别介于0.004%～3.5%、0.1%～3.4%、0.7%～3.1%和0.5%～6.7%之间，平均值依次为1.3%、1.4%、1.9%和2.5%。尽管小流域位于黄土高原水蚀风蚀交错区，4—5月风蚀强盛，7—8月水蚀严重，但土壤侵蚀仍以水蚀为主，占全年土壤侵蚀总量的80%以上[15]。此外，切沟是水力侵蚀的结果，风蚀对坡面切沟形态的影响较小，主要表现为对切沟沟岸的吹蚀与磨蚀。因此，笔者认为2018年10月至2019年5月期间切沟发育不明显。由于2019年5月所测切沟数量多、规模代表性强，以该时期的RTK GPS结果为基准，验证2018年10月WV2影像的切沟解译精度。

根据配对t检验，2种方法所测切沟长度与最大表面宽度无显著差异，周长和面积则分别在0.01和0.05的置信水平上差异显著。统计切沟形态参数的相对误差发现，切沟长度、最大表面宽度、周长和面积的相对误差范围分别为0.1%～10.6%、0.2%～16.7%、0.3%～14.4%和0.2%～17.4%，平均1.9%、2.9%、6.0%和3.8%；且误差主要集中在5%以内，分别占所测切沟总条数的94%、79%、58%和76%。配对制图误差指数范围为5.3%～39.6%，平均17.3%，同样<18%[14]，说明WV2影像解译可用于六道沟小流域的切沟研究。

尽管WV2与RTK GPS得到的切沟周长与面积差异显著，但对应散点的分布与1∶1线十分接近(图2)。为准确描述这2个形态参数，以RTK GPS结果为基准，分别建立WV2所得切沟周长与面积的转换方程(图2)。结果表明，2个方程均显著，且决定系数都十分接近1，可用于对应参数的转换。

图2 WV2影像与RTK GPS提取的切沟周长和面积及对应转换方程Fig.2 Linear conversions of gully perimeter and area extracted by WV2 image and RTK GPS

3.1.2 WorldView-2切沟解译误差的影响因素 对比分析WV2与RTK GPS提取的各切沟轮廓，WV2的切沟解译误差可能主要与沟头活跃与否、沟缘植被类型及沟壁土体状况有关。切沟沟头可分为不活跃和活跃2种类型[16]。不活跃沟头具圆角沟缘，无陷穴，沟头边壁倾斜，沟壁及沟底有植被覆盖；活跃沟头表现为锐利沟缘，有陷穴，沟底有下切、沉积和流水痕迹。相比之下，具不活跃沟头的切沟沟缘在影像上更不易识别(图3a)。测量的33条切沟中，此类切沟共5条， WV2提取的各切沟参数平均误差均在5%左右，配对制图误差指数平均为22.7%，切沟轮廓存在较大差异。具活跃沟头的切沟共28条，WV2与RTK GPS得到的轮廓较吻合(图3b)，配对制图误差指数平均仅16.3%；且除面积平均误差略高于5%外，其余形态参数的误差皆<3.5%。

植被可有效减缓切沟发育，但当沟缘被树冠遮蔽时，会给影像解译带来极大的困难[17]。六道沟小流域中的切沟沟缘多有草本和灌木植被覆盖，沟底则为乔灌草混合分布。当沟缘被灌木覆盖时，从影像上难以准确识别(图3c)。对9条沟缘以灌木为主的切沟而言，WV2解译的最大表面宽度和周长都与RTK GPS的对应结果之间存在显著差异(P<0.05)，各形态参数的平均误差介于2.0%～6.4%之间。余下24条沟缘被草本植物覆盖的切沟的解译精度则更高(图3d)，各参数的平均误差与配对制图误差指数均更小。一般来说，沟底的草本、灌木并不会影响切沟沟缘的识别，但沟底乔木在影像上多近似圆形(图3d)，当靠近沟缘时易与植被覆盖的沟缘混淆，需仔细分辨。

实地调查发现，切沟沟壁因植被稀少或沟岸塌陷出现的土体裸露往往在影像上呈白色或灰白色(图3e和f)，沟缘较易区分。出现此类状况的19条切沟的长度、最大表面宽度和面积相对误差仅3%；周长误差略高，为5.5%。其余14条沟壁无裸露的切沟各参数误差均较高。

3.2 切沟发育速率

3.2.1 集水区基本情况 将六道沟小流域划分为50个集水区，包括37个Ⅰ级集水区、12个Ⅱ级集水区和1个Ⅲ级集水区(图4e)，面积平均23.4 hm2。其中，Ⅰ-17集水区面积最小，仅为2.7 hm2；Ⅱ-3面积最大，高达120.5 hm2。集水区平均坡度介于12.1°～25.1°之间，坡度较大的集水区主要为分布在小流域北部的Ⅰ级集水区，如Ⅰ-21和Ⅰ-27，坡度均超过22°。平均坡长范围12.8～33.0 m，小流域西侧的集水区坡长总体高于东侧。平均剖面曲率和平面曲率分别介于-0.4～0.3和-0.1～0.1之间。2013—2018年均植被盖度介于27.9%～78.2%之间，Ⅰ-24土地利用以独立工矿用地为主，植被盖度最低；Ⅰ-34草地、灌木林地和有林地面积广，植被盖度最高。综合土地利用动态度在1.6%～12.8%/a之间变化。土地利用转移趋势主要表现为草地转出为灌木林地。

图3 WV2影像与RTK GPS提取的典型切沟轮廓Fig.3 Typical gully boundaries extracted by WV2 image and RTK GPS

图4 六道沟小流域集水区划分及2个典型集水区2012与2018年解译切沟分布Fig.4 Division of drainage areas in Liudaogou catchment and interpreted gullies in two typical drainage areas in 2012 and 2018

3.2.2 切沟发育状况 根据2012和2018年的WV2解译结果，统计各集水区的切沟总长度、总最大表面宽度、总周长和总面积，并基于3.1.1建立的转换方程对WV2提取的切沟总周长和总面积进行转换，计算得到各参数对应的年均变化率分别介于-2.1%～1.2%、-1.7%～0.5%、-1.8%～1.4%和-0.2%～1.5%之间，平均-0.02%、0.01%、0.01%和0.2%。这与Li等[8]利用QuickBird影像解译得到的黄土高原蔡家川流域2003—2010年切沟周长与面积年均变化率十分相近，后者范围分别为-1.0%～1.2% 、0.1%～1.4%。

沟头溯源侵蚀与沟岸侧向侵蚀是切沟侵蚀在水平方向上的2种主要形式，其结果分别体现在切沟总长度和总最大表面宽度。统计各集水区切沟的沟头后退与沟岸扩张距离发现，2012—2018年各集水区的沟头后退总距离介于-47.4～29.1 m之间，沟岸扩张总距离介于-15.7～10.4 m。其中，9个集水区(18%)的切沟总长度和10个集水区(20%)的总最大表面宽度不增反减，尤以Ⅰ-22和Ⅰ-2集水区最为剧烈， 沟头分别后退-47.4和-36.5 m，沟岸分别扩张-5.3和-15.7 m，这分别与集水区内独立工矿用地的挖方扩张、径流小区的填土扩建有关。另有26个集水区(52%)的切沟总周长减少，这可能与相邻切沟的兼并有关[8]。切沟总长度、总最大表面宽度、总周长和总面积年均增长率主要集中在0～0.5%内，分别占集水区总数的78%、80%、40%和78%，其中Ⅰ-23切沟扩张最为迅速，沟头后退9.4 m，沟岸扩张2.5 m，总长度、总周长和总面积年均变化率皆达到1.2%以上(图4a和b)。

3.2.3 切沟发育的影响因素 小流域各集水区的切沟发育速率存在明显差异。由于小流域面积较小，气候和土壤的空间差异不大，笔者仅探讨地形、植被盖度和土地利用等属性对切沟发育的影响，建立各属性与切沟各参数变化率的相关系数矩阵(表1)。

表1 集水区切沟总参数变化率与集水区属性相关系数

切沟发育速率与集水区面积通常呈显著的正相关关系[8]。然而，笔者发现切沟各参数变化率与集水区面积的均无显著相关(表1)。原因可能是集水区面积仅粗略反映了流向沟头的径流规模，实际径流量还受集水区形状、地形、土地利用、土壤性质与降雨空间分布的综合影响[9]。

坡度和坡长对坡面汇水过程及径流量大小等有重要影响[18]。切沟总长度变化率与集水区平均坡度呈显著负相关关系(P<0.01)，切沟总面积变化率与集水区平均坡长显著正相关(P<0.05)，即切沟扩张速率随坡度减小、坡长增加总体呈增大趋势。当其他条件一致时，坡度可影响坡面的承雨面积，坡度越大，承雨面积越小，降雨量与产流量相应减小[19]。陡峭坡地不易受牲畜与人类扰动可能也是切沟扩张速率与坡度呈负相关的一个原因[20]。随着坡长增加，汇流面积增大，径流量增加[21]，往往会加剧切沟扩张。

在坡向方面，迎风坡(E、S、SE、SW)集水区的切沟扩张速率总体高于背风坡(W)(图5c)。在六道沟小流域，切沟主要发生在夏季强降雨期间，此时当地以东南风为主，冬春季则多为西北风。

曲率在很大程度上控制着水流和泥沙的运移方向[22]，尤其是剖面曲率，直接影响径流的加速与减速，进而作用于土壤的侵蚀和沉积。切沟总长度和总面积变化率与集水区平均剖面曲率呈显著正相关(表1)，即剖面曲率越大，切沟扩张越快。平面曲率影响径流的汇聚和分散，对切沟发育的影响更为复杂。对各集水区的平面曲率进行分级统计发现，切沟总面积变化率随其总体无明显变化(图5d)。

图5 切沟总面积变化率随主要集水区属性的变化Fig.5 Change rate of total gully area with main drainage area properties

植被可增加地表粗糙度，促进入渗，减少地表径流，从而减缓切沟侵蚀[23]。尽管集水区年均植被盖度与切沟各形态参数变化率之间均无显著相关关系(表1)。但按照SL190—2007《土壤侵蚀分类分级标准》将植被盖度划分为4级：<45%、45%～60%、60%～75%、>75%，统计不同植被盖度的切沟总面积变化率发现，后者随植被盖度的增加而减小(图5e)。当植被盖度<45%时，切沟总面积变化率平均高达0.7%/a；当植被盖度>45%时，切沟面积变化率均降至0.2%/a左右，表明盖度大于45%时，植被可显著抑制切沟发育，这与陈一先等[24]利用1976与1997年1∶1万地形图估算陕北黄土区切沟发育速率得到的结论一致。

除切沟总周长变化率外，综合土地利用动态度与切沟各形态参数变化率均无显著相关关系(表1，图5f)。2012—2018年间土地利用变化多表现为草地转出为灌木林地，但这2种土地利用类型的产流系数并无显著差异[25]，在降水特征无明显趋势性变化的背景下，径流特征和切沟侵蚀也不会发生太大的变化。

4 结论

1)WV2目视解译提取切沟形态参数具有较高的精度。与RTK GPS实测数据相比，WV2解译切沟长度、最大表面宽度、周长和面积的相对误差主要集中在5%以内，分别占切沟样本总数的94%、79%、58%和76%。 当切沟沟头活跃、沟缘以草本植物为主或沟壁存在土体裸露时，WV2对切沟的解译精度较高，误差较小。

2)2012—2018年间，小流域各集水区切沟总长度、总最大表面宽度、总周长和总面积的年均变化率分别介于-2.1%～1.2%、-1.7%～0.5%、-1.8%～1.4%和-0.2%～1.5%之间，平均-0.02%、0.01%、0.01%和0.2%。

3)小流域切沟扩张速率与集水区平均坡长和平均剖面曲率呈显著正相关，与平均坡度呈负相关，与坡向、植被盖度也存在一定的相关关系。当位于背风坡或植被盖度>45%时，切沟扩张速率较低。