融合Sentinel-2与ASTER GDEMV2数据的坡耕地遥感提取方法研究——以子洲县为例

刘翼遥，吴太夏

（河海大学，南京 211100）

黄土高原是黄河流域的重要组成部分，同时也是我国水土流失问题最严重的地区。根据资料显示，黄土高原总面积为64万km2，但是截至2017年，其水土流失面积高达47.2 km2，侵蚀模数大于8 000 t/km2·a的极强度以上的水蚀面积占全国同类面积的64.95%[1]。黄土高原的土壤流失问题不仅与土质松软、植被覆盖度过低有关[2-4]，区域内存在大量坡耕地也是导致水土流失的重要原因[5]。

由于坡耕地的存在，使得黄土高原跑水、跑土、跑肥现象时有发生[5]。坡耕地与其他耕地类型本质的区别是其具有一定坡度，坡度属性使得坡耕地水土流失和面源污染问题十分严重。截至2012年，黄土高原坡耕地总面积占耕地总面积的2/3，平均土壤侵蚀模数2.5万t/km2·a[6]。治理坡耕地是治理水土流失最重要的步骤之一。1999年我国实行退耕还林，将易水土流失的坡耕地停耕恢复植被，黄土高原水土流失得到有效改善，效果显著[7]。但是2017年黄土高原无定河流域发生了“7·26”特大暴雨，造成陕西省子洲县等多个地区严重内涝和土壤侵蚀。子洲县坡耕地水土流失，形成了许多细沟和小切沟，严重的土壤侵蚀强度可达110 000 t/km2，坡耕地多的支沟的侵蚀强度至少为46 000 t/km2，坡耕地少的支沟为12 000 t/km2，前者是后者的3.8倍[8]。本次特大暴雨事件也反映了子洲县坡耕地面积占比依旧较大，并且分布在陡峭的沟坡，沟蚀严重，造成大量水土流失。因此，准确掌握坡耕地以及在每个坡度级的面积大小、分布位置，并在相应区域做出防护措施，就能防止悲剧再次发生，对减少土壤侵蚀和水土流失也具有重大意义。

然而坡耕地具有较为分散、不够集中、地形破碎的特点[9]，并且坡耕地普遍分布在较为陡峭的区域，这无疑给人工调查增加了难度。而遥感技术具有大范围观测地表的能力，具有全天候、多时相、涵盖信息多且连续观测的特点，并且在耕地等地物信息分类提取等行业内有大量应用，已经取得良好成果[10-12]。王祯等[6]使用Landsat 4/5结合Landsat 8数据对延安市坡耕地进行提取。陈正发等[9]使用资源环境数据云平台的土地利用数据对云南坡耕地进行识别，并对其质量进行评价。时亚坤等[13]使用GFSAD数据提取坡耕地，研究退耕还林对粮食安全的影响。这些坡耕地提取结果分辨率大多都为30 m，对于那些破碎、分散的坡耕地的提取结果不高。杨萌等[14]使用无人机航拍获取高分辨率影像提取岔巴沟和清水沟小流域的坡耕地。使用无人机航拍虽然分辨率和精度较高，但是只针对于面积较小的区域，无法对大区域做到快速提取。从2015年起，欧空局发射哨兵（Sentinel）系列卫星，其最高分辨率可达10 m，重访周期为5 d。哨兵系列卫星普遍应用于地物的提取，取得了不错的成果。

本研究融合Sentinel-2遥感影像和ASTER GDEMV2数据在子洲县建立一种大范围的坡耕地快速提取模型，对坡耕地进行大范围的识别。得到坡耕地的空间分布情况，完成研究区坡耕地时空特征的分析，为黄土高原依然存在的坡耕地的退耕还林和水土流失防治工作提供参考依据。

1 研究区概况与数据源

1.1 研究区概况

子洲县位于陕西省榆林市南部地区，全县总面积达2 024 km2，位于北纬37°15′30″~37°50′00″，东经109°29′08″~110°07′30″之间，其地理位置如图1所示。子洲县属于典型的陕北黄土高原丘陵沟壑区，地貌主要以梁、峁、沟和川为主，研究区内山区占总面积的95%，剩下的5%为川区[15]。子洲县整体地势较高，海拔最高为1 339 m，最低为886 m。总体来说，子洲县地貌特征较为统一，便于进行大范围的整体研究。

图1 研究区地理位置

子洲县是农业县，农业人口占总人口的93%，县内黄土性土壤广泛分布，占总土地面积的89.97%，宜于发展农业。子洲县主要作物为秋粮收获的玉米、大豆、谷子等，占比可达全年粮食产量的90.82%。根据《子洲县土地利用总体规划（2006—2020年）》文件显示[16]，全县基本农田约为551 km2。并且大部分耕地地理位置所处坡度较高，受到地形地势影响，耕地较为分散，不易组成统一、集中的区域，所以需要从县域角度进行坡耕地研究。

研究区位于温带半干旱区，属于大陆性季风气候，气候干燥，日照充足，冷热温差较大。子洲县年均气温9.1℃，年平均降雨量为430.8 mm，降水分配极不均匀，5—9月的降雨量占全年总降雨量的90%以上，多集中于几场高强度、短历时的暴雨中。子洲县几乎每年都会有干旱、霜冻、暴雨、大风和沙尘暴灾害性天气发生，土壤侵蚀严重，加之研究区的坡耕地较多，极易造成水土流失，使得子洲县成为榆林市水土流失的重点区域。

1.2 遥感数据来源

为满足子洲县坡耕地大范围信息提取的空间分辨率要求（10 m）和研究坡耕地变化的时间分辨率要求，本文选择Sentinel-2数据进行研究。Sentinel-2卫星是欧空局发射的双星卫星，Sentinel-2A和Sentinel-2B两颗卫星互为补充，重访周期可达到5 d。Sentinel-2卫星可覆盖从可见光波段到短波红外的13个光谱波段，不同波段的空间分辨率分别为10 m、20 m和60 m。哨兵2号卫星为总作物的提取和制图提供了较为充足且质量较好的数据源，其具体参数见表1。

表1 Sentinel-2卫星具体参数

由于研究区诸如玉米、大豆和谷子等主要作物于每年的9月至10月开始收割，故本文选择2018年至2021年每年8月前后的四幅影像进行耕地的提取。研究区范围包含了2张哨兵2号影像，在保证影像质量的前提下选择了2幅S2A影像和2幅S2B影像进行实验，数据获取的信息见表2。

表2 Sentinel-2数据获取信息

本文所采用的Sentinel-2数据下载于欧洲航天局（https：//scihub.copernicus.eu/）。接着使用欧空局官方提供的Sen2cor插件对数据进行大气校正，最后使用SNAP对校正好的数据进行拼接、按掩膜裁剪以及格式转换等操作，以供后续使用。

1.3 DEM数据

本文所使用的DEM数据是先进星载热发射和反射辐射仪全球数字高程模型（Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation）。ASTER GDEM数据由日本METI和美国NASA联合研制并免费面向公众分发，是目前唯一覆盖全球陆地表面的高分辨率高程影像数据，并在全球对地观测研究中取得了广泛的应用。V2版则采用了一种先进的算法对V1版GDEM影像进行了改进，对V1版中存在的错误做了很好的矫正，提高了数据的空间分辨率精度和高程精度。本次研究所使用的30 m分辨率的ASTER GDEMV2数据下载于地理空间数据云网站（http：//www.gscloud.cn/home）。

1.4 样本与验证数据

样本数据的好坏直接影响到提取结果的精度，应选择那些具有代表性的纯净像元进行实验[17]。为验证子洲县耕地遥感提取的精度，我们实地采集了耕地及非耕地（包括水体、不透水体、林地和其他）的样本，之后，再结合Google Earth高清卫星地图分辨率产品共选择了206个耕地样本和428个非耕地样本，各类样本均匀分布于整个研究区，它们的位置如图1所示。以上两部分数据用于耕地提取模型的建立和精度验证。其中，Google Earth高清卫星地图分辨率产品为17级，其空间分辨率为2.39 m。

2 研究方法

2.1 地物分类体系

本次研究参照GB/T 21010—2017《土地利用现状分类》的12个一级分类、73个二级分类规定对研究区地物进行分类。由于使用的卫星产品分辨率有所不同，在真彩色图像下，同一处耕地的情况存在着区别，展示了不同的耕地特征。因此在研究中，将耕地具体分为未种植作物、种植作物两种。研究区的地物类型划分出种植作物耕地、未种植作物耕地、不透水体、水体、林地及其他几类，各种地物类型的卫星影像图如图2所示。

2.2 耕地提取方法

本文采用基于规则的面向对象提取方法，建立耕地提取决策树。面向对象的决策树遥感分类方法首先对影像进行分割，合并，将很多对象合并为一个整体；然后针对每一个对象构建分类规则，建立分类决策树，从而实现对地物的分类提取[18]。

本次研究以分割阈值45.0和合并阈值85.0进行实验，再结合由光谱特征、纹理特征、几何特征、归一化差异植被指数及水体指数分析得到的划分标准，进行适宜的耕地提取规则设定，构建子洲县耕地提取决策树，接着再根据DEM影像进行坡度分级，达到提取研究区坡耕地的结果。归一化差异水体指数（NDWI）可以实现水体与陆地的分离[19]。与其他水体指数相比，NDWI能有效抑制植被信号，增强水分信号，其公式如公式（1）所示。归一化差异植被指数（NDVI）常常被用于对农作物进行分类，以及对半干旱地区降水量的预测研究，占有相当重要的位置[20-21]，其公式如公式（2）所示。

式中：NIR为影像的近红外波段；R为影像红外波段；G为绿波段。

利用上述指数区别出水体和有植被地，再通过纹理特征和光谱特征，可以进一步区分出种植作物耕地和未种植作物耕地，决策树如图3所示。得到耕地所在位置后，根据DEM影像进行坡度分级，与提取出来的耕地进行相交处理，进而得到坡耕地的分布位置。

图3 耕地提取决策树

图3所示阈值中，NDWI为归一化水体指数，NDVI为归一化植被指数，Spectral Mean（Bi）为第i波段的光谱均值，Bi_ME为第i波段的纹理均值，Retangular Fit为矩形化拟合程度。

2.3 精度评价

在地物提取工作结束后，必须对得到的提取结果进行精度评价，验证其是否有令人满意的准确性。学者常用混淆矩阵和Kappa系数进行精度评价工作。矩阵行与列分别表示分类得到结果与实际地物情况，对角线上数据表示得到正确分类的数量，利用对角线数据、行总值、列总值等数据，可以得到总体精度、Kappa系数等评价指标，混淆矩阵见表3，总体精度、Kappa系数计算公式如公式（3）、公式（4）所示。

表3 混淆矩阵

式中：N为样本的总数，n为一共有多少种分类，Xii为对角线数据，Xi+与X+i分别为某种分类的列/行总值。Kappa系数的范围为0~1，得到系数越大，表明分类一致性越高，精度相较而言更好。当Kappa系数小于0.2，表示该分类结果较差；当Kappa系数大于0.2小于0.6时，分类一致性一般；当Kappa系数大于0.6小于0.8时，表示分类精度较高；当Kappa系数大于0.8时，该情况下可以近似于分类结果与实际趋于一致。

3 结果与分析

3.1 精度评价结果

按上述方法对子洲县2018年至2021年耕地进行提取，并计算提取结果的精度。精度评价结果如图4和表4所示。

表4 耕地提取结果精度评价

图4 子洲县2018年至2021年耕地提取结果精度评价

从图表中的结果可知，2018年至2021年耕地提取精度都较好，提取结果的总体精度都大于80%，并且每年的Kappa系数都大于0.6，表明提取的耕地与实际较为一致，结果可结合DEM数据进行坡度分级，提取子洲县2018年至2021年坡耕地。

3.2 坡耕地提取结果

以1984年《土地利用现状调查技术规程》为技术依据，可将耕地以坡度为标准，划分为5种等级，分别为小于等于2°；2~6°；6~15°；15~25°；大于25°进行坡度分级。目前多少坡度范围的耕地属于坡耕地，不同研究人员给出的标准也不尽相同[22-24]。本次研究根据罗光杰等[25]的研究，将坡度大于6°的耕地认为坡耕地，进行坡耕地的提取，其结果如图5所示。

图5 2018年至2021年子洲县坡耕地提取结果

使用ArcGIS的相交工具计算4年中每个坡度级上坡耕地的面积和占比，其结果如表5、图6所示。

从表5和图6结果可知，子洲县在2018年至2020年间，耕地分布情况变化并不明显，基本都维持在520~530 km2左右，仍然维持原来的高程、坡度等级进行耕种活动。2020年至2021年，坡耕地面积从519.02 km2增加到569.48 km2，变化较为明显，总体呈增加趋势。但从每个坡度级的坡耕地面积占当年坡耕地总面积来说，变化较不明显，4年都基本维持在同一数量级上。从每个坡度级占比来看，4年中6~15°坡度级的占比最大，都大于50%，最高可达57.28%，其次是15~25°和大于25°坡度级。该种现象与研究区沟壑起伏变化大等情况有关，该因素也直接影响子洲县耕地分布总体较为零散，以未来发展角度，不宜进行更高效率的大机器耕作方法，对于农业县来说，对其发展并不友善。根据陈正发[24]的研究，大于25°的坡耕地极易发生水土流失，应当严格实行退耕还林。但是这4年中子洲县大于25°的坡耕地依然存在，并且面积都大于30 km2，在2021年甚至达到了44 km2，需要当地政府在平衡耕地红线的情况下，更好更快落实退耕还林、还草的相关政策。

表5 2018年至2021年每个坡度级中坡耕地面积占比

图6 2018年至2021年每个坡度级中耕地面积

4 结论

本研究以陕西省榆林市子洲县作为研究区域，利用Sentinel-2和ASTER GDEMV2作为数据来源，基于规则的面向对象分类方法，确定影像分割、影像合并阈值大小，以及各种地物类型的提取规则，并最终建立关于种植耕地、未种植耕地及非耕地的提取决策树，基于实地勘察数据和googleearth影像，对提取结果进行验证，得到Kappa系数高于0.65的分类结果，可认为一致性较高，精度较好。接着对2018至2021年，子洲县同期耕地时空变化情况进行分析，发现子洲县在2018年至2020年间，耕地分布情况变化并不明显。2021年，坡耕地面积增加到569.48 km2，变化较为明显。4年中6~15°坡度级的占比最大，都大于50%，最高可达57.28%，其次是15~25°和大于25°坡度级。并且这4年中子洲县大于25°的坡耕地依然存在，并且面积都大于30 km2，在2021年甚至达到了44 km2。