气候变化是当今时代面临的首要生态问题
[1]，其加剧了生态环境退化，提升了自然灾害的发生频率
[2]。与此同时，世界经济局势处于调整期，全球总人口仍在高速增长
[3]，气候、经济、人口的多重压力，使得粮食安全的重要性不断凸显
[4]。土地是粮食生产的基础，中国东北黑土区是世界四大黑土区之一
[5]，具有肥力高、分布性状好、适宜农作物生长的特征
[6]，其粮食产量和耕地面积均约占中国总量的25%，是中国的主要粮食产区和商品粮基地
[7-9]。漫川漫岗黑土区地跨辽宁、吉林和黑龙江三省，主要包含松嫩平原和辽河平原北部，区内年降水量约500 mm，径流易于集中。在自然原因和人类改造地表的双重作用下，土层遭到水土流失破坏，并使黑土功能在近几十年来发生严重退化
[10-11]，东北黑土区因水土流失而使粮食减产约10%
[12]，坡耕是漫川漫岗黑土区水土流失的重要原因，并以侵蚀沟为主
[13]，造成可耕作面积减少、机械化效率降低
[14]，严重威胁粮食和生态安全。
根据《东北黑土区侵蚀沟治理专项规划（2016—2030年）》
[15]，获取侵蚀沟的总体空间分布是开展治理的前提，利用遥感数据提取侵蚀沟的方法主要有：目视解译、基于阈值分析的半自动提取，以及基于机器学习的自动提取。其中，目视解译方法虽精度较高，但需耗费大量人力和时间；半自动提取方法一般为通过设置影像灰度值及地形特征的阈值，将侵蚀沟从背景中分离出来；自动提取方法一般为利用已知侵蚀沟样本训练得到侵蚀沟自动提取模型，并开展自动提取，主要包括：：像元分类法、面向对象的分类法、图像分类法（判断某个矩形图像内是否含有侵蚀沟）和语义分割法。目前关于侵蚀沟自动提取的研究较少，为探索适宜漫川漫岗黑土区的侵蚀沟自动提取方法流程，并验证自动提取效果，以哈尔滨市巴彦县北部为研究区，利用夏季、春秋季的高分辨率卫星多光谱影像、数字高程模型（digital elevation model，DEM）和第三次全国国土调查数据（三调），通过目视解译获得了侵蚀沟的空间分布，构建36个并优选得到27个特征变量，在研究区内选择某一区域作为预测区，将目视解译获得的侵蚀沟作为训练样本，采用面向对象的分类法，开展了侵蚀沟的自动提取研究。以期验证面向对象分类方法在区域性侵蚀沟自动提取方面的应用效果，并获取适于面向对象分类方法的特征变量组合，助力黑土保护及粮食安全。
1 研究区概况
巴彦县行政隶属于哈尔滨市，位于松花江北岸，松嫩典型黑土带内
[16]。将巴彦县北部作为研究区（
图1），包括黑山镇、万发镇、天增镇、洼兴镇、德祥乡、山后乡、红光乡、兴隆镇共8个乡镇，总面积约1 600 km
2。研究区具有中温带大陆性季风气候，夏季凉爽短暂、冬季寒冷漫长，年平均降水量约560 mm，6—9月为主要降雨期，11月—次年1月为主要的降雪期。区域地貌上位于张广才岭与松嫩平原的过渡带，高程位于150~470 m，东北高、西南低，地貌类型包括低山丘陵、平原及其过渡带内的漫川漫岗。
2 数据方法
2.1 数据来源
用于开展侵蚀沟目视解译的春秋、夏季高分辨率卫星多光谱影像来自于黑龙江测绘地理信息局，其空间分辨率为1 m，时相为2021年5月—2022年9月；用于构建驱动特征变量的ALOS DEM来自美国国家航空航天局，其空间分辨率为12.5 m，是目前空间分辨率最优的免费数据；用于辅助开展侵蚀沟目视解译及自动识别的“三调”成果数据来自中国地质调查局哈尔滨自然资源综合调查中心，如
表1所示。各类数据经投影变换和重采样处理，具有正确统一的坐标系统和空间分辨率。
2.2 研究方法
在收集数据的基础上，开展信息构建及自动提取，其中信息构建包括侵蚀沟目视解译和驱动特征变量、识别特征变量的构建，自动提取包括特征变量优选、图像分割、随机森林建模提取及提取后处理，如
图2所示。
2.2.1 信息构建
（1）侵蚀沟解译。利用春秋、夏季高分辨率卫星多光谱影像，采用目视方法判读和勾勒侵蚀沟轮廓。侵蚀沟为地表的汇水源头，随着其向下游的延伸和汇水面积的增大，会形成稳定的地表径流，容易与河流、沟渠发生混淆，因此需将“三调”用地类型作为辅助数据，只将在“三调”中河流、沟渠之外的能够目视判别的自然沟道作为侵蚀沟。侵蚀沟分类体系多样
[17-20]，根据发育程度可划分为细沟（边界不明显，深度小于犁底层，可被农业机械跨越并通过耕作消除）、浅沟（0.3~0.5 m宽，可被农业机械跨越但不能通过耕作消除）、切沟（大于0.5 m宽，呈条带状，横截面呈V形，沟底与坡面的剖面形态相近）、冲沟（大于0.5 m宽，呈条带状，横截面呈V形，沟底相对于坡面的剖面形态更加下凹，存在流水沉积物）、坳沟（大于0.5 m宽，呈宽条带或面状，横截面呈U形，沟底剖面相对平缓，存在大量流水沉积物）（
表2）。结合卫星多光谱影像的空间分辨率，此次将切沟、冲沟、坳沟作为目视解译对象。
（2）特征变量构建。利用DEM数据，通过通用土壤流失方程（universal soil loss equation，USLE）
[21-22]计算得到高差、坡度、坡度因子、坡长因子、汇水线密度、距汇水线距离6个驱动特征变量，如
表3所示。灰度共生矩阵（gray-level cooccurrence matrix，GLCM）包含影像某一滑动窗口内各个方向距离上的灰度值共生特征，像元处的GLCM统计值能够反映影像中的纹理特征
[23-24]，此次在春秋、夏季影像的蓝、绿、红、近红外（near-infrared，NIR）波段和计算得到的植被指数（normalized difference vegetation index，NDVI）
[25]的基础上，采用5*5像元的滑动窗口和2像元的步长，拟计算蓝、绿、红、近红外、NDVI的GLCM并统计0°、45°、90°、135°共4个方向的均值、均一性、非相似性、熵、角二阶矩均值，共得到30个识别特征变量，如
表3所示。
2.2.2 自动提取
（1）特征变量优选。将目视解译得到的5 000余条侵蚀沟以及随机拾取的5 000处半径为20 m的非侵蚀沟圆形作为原始训练集，侵蚀沟和非侵蚀沟样本分别赋值为1和0，采用5折交叉验证法将原始训练集划分为5组训练-训练验证集，采用随机森林算法，分别利用5组训练-训练验证集训练得到5个预测模型，并计算5个预测模型的平均绝对误差（mean absolute error，MAE）
[26]。将MAE均值作为精度参考，并采用递归特征消除算法
[27-28]开展特征变量优选，该方法属于迭代算法，能够降低特征变量数量并提高模型精度。在第1轮递归特征消除中，分别利用
n个特征变量以及去除每个特征变量后剩余的
n-1个特征变量开展侵蚀沟预测与验证，得到
n+1组参考精度，当去除某个特征变量后能够得到最优的参考精度时，则去除该特征变量，开展第2轮递归特征消除。在第2轮递归特征消除中，分别利用
n-1个特征变量以及去除每个特征变量后剩余的
n-2个特征变量开展侵蚀沟预测与验证，得到
n组参考精度，当去除某个特征变量后能够得到最优的参考精度时，则去除该特征变量，开展第3轮递归特征消除。依次类推，直到去除特征变量后并不能优化参考精度为止，则此时剩余的特征变量即为优选得到的特征变量。
（2）侵蚀沟自动提取。利用优选得到的特征变量组合，开展图像分割，经小区试验，采用光谱特征0.6、光滑度0.5的权重和10的分割尺度，获取与侵蚀沟实际边界较为一致的分割图斑[
图3（a）]。将研究区划分为13个区[
图3（b）]，根据侵蚀沟的空间分布特征，将1~5区和7~13区内与侵蚀沟分布较为一致的约300 000个分割图斑以及随机拾取的约300 000个非侵蚀沟分割图斑作为原始训练集，侵蚀沟图斑和非侵蚀沟图斑分别赋值为1和0，并利用随机森林算法开展模型训练。随机森林属集成学习算法，其应用广泛，具有数据无需归一化的优势
[29-30]，此次将决策树数量设置为1 000，节点选择的特征变量数量为9。将第6区内分割得到的图斑作为预测集，利用训练得到的随机森林模型，开展侵蚀沟自动提取，并将预测值大于0.90的图斑作为侵蚀沟原始提取结果。由于侵蚀沟的细长带状几何特征易与道路、河流、水渠、建筑等混淆，因此需参照“三调”用地类型，擦除易混淆地物并剔除散小图斑，得到侵蚀沟最终提取结果。将第6区内的侵蚀沟目视解译结果作为参考，构建混淆矩阵，并采用用户精度、生产者精度和Kappa系数验证侵蚀沟的自动提取效果。
3 结果
3.1 信息构建
目视解译结果显示，侵蚀沟主要分布于河谷平原和低山丘陵过渡带的漫川漫岗区，用地类型以坡耕地为主，随机生成的非侵蚀沟圆形分布于全区（
图4）。将侵蚀沟与各类特征变量进行空间分析，得到不同特征变量对侵蚀沟的区分能力。其中，驱动特征变量方面（
图5），侵蚀沟处的像元主要位于20~40 m高差区间、0°~5°坡度区间、0~1坡度因子区间、0~8坡长因子区间、3.7~4.7 km/km
2汇水线密度区间及0~50 m距汇水线距离区间。识别特征变量中的波段特征方面，经初步目视判断，夏季真彩色影像（
图5、
图6）中的侵蚀沟主要呈翠绿、灰黑色，耕地呈深绿色，道路呈灰白色，水体呈灰黑、墨绿色；春秋季真彩色影像中的侵蚀沟主要呈翠绿、灰黑色，耕地呈土黄色，道路呈灰白色，水体呈灰黑、墨绿色；夏季NDVI中的侵蚀沟处像元值低于周边的耕地，低于或近于水体周边的草地；春秋NDVI中的侵蚀沟处像元值高于或近于周边耕地，低于或近于水体周边的草地。识别特征变量中的GLCM特征方面，夏季绿GLCM均值、夏季NIR GLCM均值、春秋蓝GLCM熵、春秋绿GLCM熵、春秋红GLCM熵、春秋NIR GLCM熵、春秋NDVI GLCM熵共7个特征变量中的侵蚀沟处像元值具有相对较高的数值，夏季红GLCM均值、春秋NIR GLCM均值、春秋NIR GLCM均一性、春秋NDVI GLCM均一性共4个特征变量中的侵蚀沟处像元值具有相对中等的数值，夏季NDVI GLCM均值、春秋蓝GLCM角二阶矩、春秋绿GLCM角二阶矩、春秋红GLCM角二阶矩、春秋NIR GLCM非相似性、春秋NIR GLCM角二阶矩、春秋NDVI GLCM非相似性、春秋NDVI GLCM角二阶矩等8个特征变量中的侵蚀沟处像元值具有相对较低的数值，春秋NDVI GLCM均值中的侵蚀沟处像元值可能具有相对较低或较低的数值，如
图5、
图6、
表4所示。
3.2 自动提取
3.2.1 特征变量优选
以目视解译得到的侵蚀沟及随机拾取的半径为20 m的非侵蚀沟圆形（
图4）为训练样本，结合目视初步选择的特征变量（
表4），设置随机森林模型的决策树数量为1 000，节点处随机选取的特征变量数量为特征变量总数的1/3
[31-33]，利用随机森林模型开展训练并将其验证结果作为参考精度，采用递归特征消除方法开展特征变量优选。第1次迭代去除汇水线密度变量后，MAE达到最低，为0.169 2（
图7）；第2次迭代去除坡度因子变量后，MAE达到最低，为0.168 0；第3次迭代去除坡度变量后，MAE达到最低，为0.167 4；第4次迭代去除春秋NDVI非相似性变量后，MAE达到最低，为0.167 2；第5次迭代去除春秋NDVI均一性变量后，MAE达到最低，为0.166 4；第6次迭代去除春秋NDVI熵变量后，MAE达到最低，为0.165 4；第7次迭代去除春秋NDVI角二阶矩变量后，MAE达到最低，为0.165 4；第8次迭代去除夏季NDVI变量后，MAE达到最低，为0.164 9；第9次迭代去除春秋NIR均值变量后，MAE达到最低，为0.164 6；第10次迭代去除任一变量后，MAE均未低于0.164 6，迭代终止，最终优选得到27个特征变量（
表5）。
3.2.2 侵蚀沟自动提取
以第1~5区和7~13区内侵蚀沟、非侵蚀沟图斑为训练样本，设置随机森林模型的决策树数量为1 000，节点处随机选取的特征变量数量为特征变量总数的1/3，针对第6区内的分割图斑，利用随机森林模型开展训练并进行自动提取，得到侵蚀沟原始提取结果，如
图8（a）所示。侵蚀沟原始提取结果由侵蚀沟、耕地及大量的道路、河流、水渠、建筑物组成，在第6区内的总面积约16.46 km
2，利用“三调”用地类型数据，在擦除道路、河流、水渠、建筑物缓冲区并剔除散小图斑后，得到侵蚀沟最终提取结果，在第6区内的总面积约4.06 km
2，空间分布较为均匀，如
图8（b）所示。
以第6区内目视解译得到的侵蚀沟为参考，计算得到自动提取侵蚀沟与目视解译侵蚀沟的混淆矩阵，结果显示，侵蚀沟自动提取的用户精度为54.28%，生产者精度为72.16%，总体精度为99.62%。由于侵蚀沟面积占第6区的比例很小，会产生类别不均衡效应，使得总体精度参考意义不大，因此采用Kappa系数反映总体的自动识别效果，结果显示Kappa系数为0.61，通过该自动提取方法能够获得区域侵蚀沟的总体分布特征。
4 讨论
4.1 特征变量与侵蚀沟关系
4.1.1 驱动特征变量
高差为侵蚀强度评价的重要特征
[34]，研究区侵蚀沟处的高差与其他黑土区相近，主要位于10~40 m
[35-36]；坡度因子和坡长因子可由坡度、坡长计算得到，能够分别反映斜坡的坡度和长度特征，其为流水侵蚀提供了直接动力；汇水线为通过相应分辨率的DEM获取的地表径流位置，强降水引起的地表径流是驱动侵蚀沟产生的主要特征之一
[37]，而汇水线密度能够反映区域内地表径流的密集程度，且侵蚀沟的发育位置理论上应与某条汇水线重合，因此汇水线密度和距汇水线距离可用于区分侵蚀沟与周边地物。研究区内侵蚀沟具有在漫川漫岗区相对集中分布的特征，空间统计分析结果显示，侵蚀沟处的高差、坡度、坡度因子、坡长因子、汇水线密度、距汇水线距离均位于某一特定区间，而特征变量优选结果显示，坡度、坡度因子、汇水线密度的加入会降低侵蚀沟自动提取的精度，因此三者被从侵蚀沟自动提取特征变量组合中剔除，原因在于：①0°~5°的坡度区间、0~1的坡度因子区间属于区域背景特征，同时，相似坡度、坡度因子具有沿等高线展布的特征，与侵蚀沟具有横向交切的关系，从而使侵蚀沟与非侵蚀沟处的坡度、坡度因子均以0°~5°和0~1为主，导致侵蚀沟处的坡度、坡度因子虽具有统计特征，但不具有区分侵蚀沟、非侵蚀沟的作用；②距汇水线距离反映的是某一像元和某条一维的汇水线之间的空间关系，相近的距汇水线距离趋向于呈带状、树枝状分布，更符合侵蚀沟呈条带状延伸的实际，而汇水线密度更多的是反映某一像元和多条一维汇水线组成的二维面状区域的空间关系，其计算结果趋向于呈团状分布，不利于区分侵蚀沟与周边地物，同时其劣势进一步被DEM的空间分辨率误差放大。
4.1.2 识别特征变量
灰度共生矩阵（GLCM）能够反映影像的高级纹理特征
[23-24]，其中均值代表纹理的规律性，值越大则纹理规律性越强，熵代表纹理的复杂度，值越大则纹理越杂乱，角二阶矩代表纹理的规律性和频率，值越大则纹理越规律和细致，均一性代表灰度的一致性，值越大则灰度值越相近，非相似性能够反映线性对比度，值越大则线性对比度越高
[38-42]。通过GLCM中侵蚀沟与主要地物处的像元值相对大小可知（
表4），侵蚀沟在GLCM均值特征中具有中低-中高值，反映具一定规律性的纹理特征，在GLCM熵特征中具有中高值，反映较杂乱的纹理特征，在GLCM角二阶矩特征中具有低值，反映粗糙的纹理特征，在GLCM均一性特征中具有中值，反映较为相近的灰度值，在GLCM非相似性特征中具有低值，反映较差的线性纹理对比度。变量优选结果显示，春秋NIR均值、春秋NDVI均一性、春秋NDVI非相似性、春秋NDVI熵、春秋NDVI角二阶矩等5个GLCM特征变量因不利于侵蚀沟自动识别而被剔除，这5个变量反映的是春秋季侵蚀沟及周边地物在植被特征方面的纹理差异，虽然在春秋季时部分侵蚀沟内因生长有植被而易于通过目视与裸露的耕地相区分，但春秋季侵蚀沟与周边地物在植被长势纹理的一致性、对比度、杂乱度、规律性方面的差异并不容易通过统计特征相区分，不同物候时期的影像及由其构建的特征变量对侵蚀沟的区分效果并不一致。
4.2 侵蚀沟自动提取
可以用于侵蚀沟自动提取的方法包括像元分类法、面向对象分类法、图像分类法、语义分割法，其中像元分类法将每个像元作为研究客体，其训练预测过程中的数据量大、效率低、不包含几何信息；面向对象分类法的思路为先将图像进行分割，然后以每个对象为研究客体，降低了数据量、效率较高、引入了几何信息，但分类精度易受分割效果影响；图像分类法属于卷积神经网络算法（深度学习），能够充分提取侵蚀沟及背景的空间几何、纹理关系，但该方法只能判断某一场景图像区域内是否包含侵蚀沟，无法提取具体轮廓；语义分割法同样采用了卷积结构，其将场景图像作为输入，将像元级别的侵蚀沟轮廓作为输出，在地物分类研究中逐渐趋于主流。此次以第6区内目视解译的侵蚀沟为参考，验证了面向对象分类方法在区域侵蚀沟自动提取方面的应用效果，并针对该方法初步筛选得到适宜的地貌、水文、光谱、灰度共生矩阵特征组合。测试显示用户精度为54.28%，生产者精度为72.16%，Kappa系数为0.61，表明该方法和特征组合能够获取区域性侵蚀沟总体特征，但提取精度仍具有较大的提升空间，可在今后进一步应用语义分割方法并与之对比。此次自动提取误差主要来源于2方面：①为分割图斑的边界误差，表现为图斑边界呈锯齿状，分割参数未能适宜整个研究区而导致图斑边界与实际侵蚀沟边界具有偏移；②为特征变量组合无法完全区分侵蚀沟与周边地物，此次训练过程中使用的正样本均为侵蚀沟，但对比
图8中的侵蚀沟原始提取结果和最终提取结果，可见自动提取模型仍会将大量道路、河流、水渠、建筑误识为侵蚀沟，在识别特征变量纹理特征之外，应进一步增强驱动特征变量尤其是距汇水线距离数据的分辨率，从而更加有效的区分侵蚀沟与非侵蚀沟。侵蚀沟自动提取虽不能获取与目视解译同样精细的侵蚀沟轮廓，但能够快速获取侵蚀沟总体分布及空间位置，便于辅助制定治理防护计划，促进黑土保护与科学管理。
5 结论
以哈尔滨市巴彦县北部为研究区，通过构建特征变量、特征优选、图像分割及随机森林算法，开展了侵蚀沟的自动提取研究，得到以下结论。
（1）侵蚀沟所处地貌属低山丘陵与河谷平原过渡带的漫川漫岗区，主要分布于坡耕地内。
（2）高差、坡长因子、距汇水线距离能够较好地反映侵蚀沟与周边地物的地貌、水文特征差异；受沟内及沟边地表覆盖的影响，侵蚀沟表现为稍具规律性的、杂乱粗糙的纹理特征，能在一定程度上通过灰度共生矩阵的统计特征与周边地物相区分。
（3）通过图像分割及机器学习的方法流程，能够自动提取获得区域性侵蚀沟分布特征，但易将侵蚀沟与道路、河流、水渠、建筑等几何、纹理特征相近的地物相混淆，因此需充分利用高精度DEM构建地貌、水文类的特征变量，从侵蚀沟驱动要素方面入手提高自动提取的效果。
PDF
(42454KB)
