基于遥感的撂荒耕地提取方法比较研究

刘鑫鹏，李栋梁，赵天成，2，黄 晴*，赵 童，王思悦，颜 翔，杨竣程

(1.南京晓庄学院 环境科学学院，江苏 南京 211171；2.南京师范大学 地理科学学院，江苏 南京 210023)

0 引言

全球气候持续变暖以及由此引起的生态环境问题，成为21世纪人类社会面临最为严重的挑战之一.已有研究表明，全球气候持续变暖将引起小麦、玉米和水稻等世界主要粮食作物大幅减产[1， 2].研究表明，温度升高会使昆虫繁殖和代谢加速，如果全球平均气温升高2℃，害虫导致的全球小麦、玉米和水稻产量损失将分别增加46%、31%和19%[3].粮食安全问题已经成为人类社会可持续发展的关键问题，特别是新冠肺炎疫情全球蔓延以及地缘政治冲突更加剧了人们对于粮食安全问题的担忧[4].

近年来，随着城乡一体化发展和社会经济转型的快速推进，在农业生产效益低下与农村劳动力向外转移等多重因素的影响下，农民耕种的积极性不断减退，大量农村土地闲置，耕地撂荒现象不断加剧，在经济欠发达，农业生产条件落后，地形复杂的地方更为严重.中国以仅占世界7%的耕地，养育了约占世界22%的人口.我国非常重视粮食安全问题，并始终将农业放在国民经济的基础地位.耕地撂荒现象的不断加剧，将对我国粮食供给产生较大影响[5].因此，开展对撂荒耕地遥感监测研究，有利于大范围及时掌握撂荒耕地的空间分布，进而为农业部门科学制定应对措施提供依据，为农业相关规划以及政策制定提供参考.对我国坚守18亿亩耕地的红线[6]，维护粮食安全具有十分重要的意义.

耕地撂荒是在经济发展以及城市化推进过程中出现的一种情况，在全球各地都普遍存在，在许多国家和地区，如欧洲东部地区[7]、地中海地区[8]、拉丁美洲[9]以及美国、日本和澳大利亚[10]等，都有耕地撂荒情况出现.许多学者也围绕耕地撂荒的监测开展了深入研究，一些学者采用实地调查法，如雷锟等利用农户评估法对12个典型村庄耕地撂荒情况开展农户调查，并分析山区撂荒耕地形成的主要原因[11]，尽管实地调查法具有较高的精度，但难以大面积开展.遥感具有大面积、多时相观测的特点，特别是多源遥感数据为大范围开展耕地撂荒监测提供了重要支撑.基于遥感的监测方法成为当前耕地撂荒监测的主要途径.如程维芳等采用MODIS NDVI数据产品和ALOS遥感数据，并根据不同土地利用类型植被生长周期特征和NDVI特征，建立了撂荒地识别方法[12]；Alcantara等以罗马尼西亚以南地区作为研究区，使用Landsat影像绘制土地利用图，对1990—2005年度耕地的撂荒程度进行了分析，得出该研究区的撂荒率为21.1%[7]；牛继强等基于支持向量机的方法，通过对4个不同日期遥感影像分类结果的比较来提取季节性及常年性撂荒地[13]；Estel等基于MODIS时间序列数据，使用随机森林分类的方法提取撂荒耕地，发现2001～2012年欧洲耕地撂率达到24.4%[14]；肖国峰等人采用CART决策树分类方法，基于Landsat数据和HJ1A数据提取撂荒地空间分布并进一步获得了持续撂荒时间分布和撂荒地复垦区域[15]；朱小文等以甘肃省平凉市崆峒区为研究区，通过解译Landsat TM/OLI 遥感影像分析了该地区近20年耕地撂荒的规模和分布特征[16]；王玲玉等利用不同分辨率数据提取了石漠化地区撂荒耕地的时空分布，撂荒地综合识别精度达90.7%[17].总体上，多源遥感数据已被广泛应用于不同尺度的撂荒耕地遥感监测，但监测方法仍然存在较大差异，不同方法的精度和适应性仍不明确.

西南地区是我国典型的丘陵山地地带，经济发展水平相对较低，耕地撂荒现象较为严重.本文以清镇市为典型研究区，通过分析研究区耕地主要作为的种植特征，结合不同季节的哨兵2号遥感影像数据，通过构建单时相波段法、植被指数变化检测法和CART 决策树分类方法对研究区撂荒耕地进行提取，分析了不同的提取精度及适应性，研究成果将为其他相似区域撂荒耕地监测提供实践参考.

1 研究区和数据来源

1.1 研究区概况

本文选取贵州省贵阳市下辖县级市清镇市为研究区(图1).清镇市位于贵州省中部，地处黔中腹地，东经106°07′至106°33′，北纬26°21′至26°59′，总面积1386.6 km2.研究区地势南高北低，大部分海拔在1180～1450 m之间，自然地貌属于以山地、丘陵为主的丘陵盆地；气候属于亚热带季风湿润气候，年平均气温14.1 ℃，季风明显；由于地形相对崎岖，不适合大面积种植，农业机械化程度较低，经济水平相对落后，外出务工人数较多，区域内耕地撂荒现象较为明显.

图1 研究区地理位置图

1.2 数据与处理

本研究利用数据包括遥感影像数据、土地利用数据以及农作物种植数据.遥感数据采用的是Sentinel-2A卫星影像，数据来源于欧空局(https://scihub.copernicus.eu)，影像共4期，成像时间分别为2019年4月9日和8月12日，2019年10月01日和12月15日，能够较好地反映研究区一年中春夏秋冬 4个季节撂荒耕地和正常耕地的季节变化，并获取不同季节的光谱特征差异[18].为消除大气和光照等因素对地物反射的影响，得到地物真实反射率[19]，本研究使用的Sentinel-2A卫星数据均经过大气校正处理，并利用影像的红光波段和近红外波段计算NDVI值.

土地利用数据来源于中国科学院资源环境数据中心(https://www.resdc.cn/)30 m分辨率的遥感解译数据，同时结合Google Earth影像和经过波段融合的2 m分辨率高分1号遥感影像作为样地采集数据的基础数据.本研究利用土地利用数据对影像耕地范围进行掩膜，并选取了73个撂荒地样地和27个耕地样地进行光谱分析及结果验证.

农作物种植数据包括研究区农作物种类，种植熟制以及物候特征，数据来源于清镇市政府网站公开数据及文献总结[20，21]，并作为研究区耕地数据的辅助数据，详细情况见表1.

表1 农作物物候特征统计表

2 研究方法

2.1 单时相波段特征检测法

单时相波段特征检测法是以单个季节的遥感影像作为基础，利用地物光谱特征获取光谱反射率，形成光谱曲线来形象化区分不同地物.由于不同地物在不同光谱波段中的反射率存在差异，可通过观察每个地物在某个波段的反射率差值来构建识别模型.本研究在掩膜后的耕地范围内对当地农作物典型地物的特征进行统计分析，根据光谱特征曲线差异来构建识别模型.

2.2 植被指数变化检测法

植被指数变化检测法是根据同一植被在不同时间的观测，来确定其变化的处理过程.通过检测可以找出感兴趣的变化信息，滤除作为干扰因素的不相干的变化信息.由于需要处理同一地点不同时间的多幅影像，通常处理的数据量更多、数据异质性更强、地物情况更复杂[22].归一化植被指数(NDVI)是应用最为广泛的植被指数之一，能够较为准确的反映植被的变化.由于不同类型和不同季节农作物生长状况不同，在影像上的光谱特征也存在较大差异，在植被生长季夏季，种植农作物的区域与长满野草的撂荒耕地，其NDVI差异显著，作物区NDVI明显大于撂荒地NDVI.而在播种季春季以及收获季秋季，作物区NDVI明显减小，因此可利用农作物NDVI的季节变化对撂荒耕地进行识别[23].本研究通过分析撂荒耕地和正常耕作耕地在春季和夏季、夏季和秋季连续2个及以上种植季末NDVI的差异来有效的提取撂荒耕地.

2.3 CART决策树分类法

CART(Classification and regression tree)决策树分类法通过建立二叉树，对分类特征进行二分，迭代生成决策树.决策树分类方法在一定程度上可以解决遥感影像中“同物异谱”、“异物同谱”的问题.CART 决策树由于模型简单，且具有自动根据训练样本确定分类阈值，自动建立决策树，不受其他因素影响等优势，分类精度相对较高[24].

3 结果分析

3.1 单时相波段特征提取分析

利用4个季节不同时间点的哨兵2号遥感影像，分别提取了撂荒耕地、水稻、玉米以及蔬菜等地物的光谱特征曲线，如图2所示.夏季影像中撂荒地光谱特征与其他地物差异明显，其中波长为0.740～0.865 μm范围内光谱反射率高于其他地物，并且红光(b4)波段反射率最低，低于其他地物.春季、秋季和冬季的影像的地物光谱反射率差异较小，但总体上各获取的地物光谱曲线趋势相似，很好的反应了这三个季节农作物之间的地物光谱反射率存在着相似性.

图2 典型地物光谱曲线图

根据夏季影像中撂荒地的光谱特征构建出识别模型来提取试验区的撂荒地，将红边波段与红光波段的光谱反射率差值大于0.2的地方定义为撂荒地.提取结果如图3所示.基于单时相波段特征提取的撂荒耕地栅格118382900个，分布范围广，提取的斑块非常破碎，撂荒耕地多分布于清镇市的西北部的新店镇和中部地区的王庄布依族苗族乡.通过对比高分辨率影像，发现所提取的大部分地区都没有撂荒的迹象，主要是一些其他作物种植的耕地与撂荒耕地表现为相似的光谱特征，误判率较高.

3.2 植被指数变化检测法的提取分析

对耕地与撂荒地不同季节的NDVI变化情况进行统计，如表2所示.夏季是作物的主要生长季，利用夏季和春季以及夏季和秋季的NDVI值差值，能够有效识别耕地和撂荒地.统计结果表明，耕地和撂荒地不同季节NDVI差异的最小值和均值均表现出较明显的不同，夏季与春季的NDVI差值最小值为0.29，小于此值的区域不是耕种的土地，而是撂荒地，同时，夏季与秋季的NDVI差值最小值为0.30，小于此值的区域也不是耕种的土地.从最小值看，两个不同季节NDVI与夏季NDVI差值的最小值较为相近，通过对不同类别不同样点的数值进行了多次的实验，本研究将阈值设定为0.36.即NDVI差值小于0.36的区域为撂荒地.

表2 耕地与撂荒地不同时相NDVI差值的统计分布情况

根据上述提取模型，对研究区撂荒耕地进行了提取，结果如图4所示.基于植被指数变化检测法提取的撂荒耕地栅格51203800个，显著小于单时相波段特征提取结果.撂荒耕地大部分分布于清镇市的西北地区，并且呈现出破碎的斑块分布.高分辨率影像进行对比表明，清镇市的东南地区为地势相对较低的县中心地区，坡度也较小，耕地的质量和耕种率相对较高，而东部为森林覆盖区，出现撂荒的现象较少.撂荒地较为集中的新店镇其地势中部和东南部较高，坡度大，自然环境条件恶劣是撂荒的主要原因.

图3 基于单时相波段特征提取的撂荒地分布图

图4 基于植被指数变化检测法提取的撂荒地分布图

3.3 CART决策树分类法提取分析

根据单时相波段撂荒耕地的光谱特征分析结果，夏季光谱值区分度最强，同时结合变化监测NDVI变化分析结果，本研究利用哨兵2号夏季影像的7个波段以及夏季NDVI构建决策树模型，并对试验区遥感影像构建ROI(感兴趣区).CART决策树自动建立规则，详见图5.

图5 决策树自动构建规则

基于CART决策树分类法提取的撂荒耕地分布如图6所示，共提取撂荒耕地栅格60010700个.总体上，提取的撂荒耕地分布主要分布于西北部的乡镇，分布较为破碎，空间分布上与基于变化监测方法提取的结果相似.

图6 基于决策树模型构建的撂荒地分布图

3.4 撂荒耕地提取方法对比分析

3.4.1 提取精度分析

本文采用基于混淆矩阵的精度验证方法对提取结果进行精度验证.评价指标包括总体分类精度和制图精度.总体分类精度是指对每一个随机样本，所分类的结果与检验数据类型相一致的概率，是具有概率意义的一个统计量.制图精度是表示相对于检验数据中的任意一个随机样本，分类图上同一地点的分类结果与其相一致的条件概率.

利用本文选取的73处撂荒地和27处耕地的样地数据，建立精度验证混淆矩阵如表3所示.单时相波段法提取的撂荒地制图精度仅为58.90%，耕地的制图精度为85.19%，总体精度也只有67%；植被指数变化检测法提取的撂荒地制图精度为79.45%，耕地的制图精度为88.89%，总体精度为82%，相比单时相波段法，提取精度有明显的提升；CART决策树分类法提取的撂荒地制图精度只69.86%，而耕地的制图精度为81.48%，分类总体精度73%.从分类精度看，基于作物季节变化特征的植被指数变化检测法具有较高的提取精度.

表3 精度验证混淆矩阵

3.4.2 不同方法对比

单时相波段法仅从不同地物的光谱特征值差异角度进行提取，容易受到其他光谱值相近的地物的影响，因此分类精度较低，总体精度仅为67%，显著低于变化监测法和CART决策树分类方法；植被指数变化检测法充分利用了作物生长季和非生长季NDVI值的差异，一定程度上减少了相近地物对于提取结果的影响，提取精度较好，达到82%，CART决策树分类法综合考虑了光谱值以及NDVI的差异，提取精度虽然高于单时相波段法，但低于植被指数变化检测法.因此，考虑植被生长季节变化的植被指数变化检测法更适用于撂荒耕地的提取.

图7 共同提取的撂荒耕地分布图

本研究将三种不用方法提取的结果进行叠加来分析不同方法共同提取的撂荒耕地的空间分布.单时相波段特征法与植被指数变化检测法提取的撂荒地重叠区面积为69920个栅格点，单时相波段特征法与CART决策树分类法提取的撂荒地重叠区面积为82003个栅格点，而植被指数变化检测法提取的撂荒地与CART决策树分类法提取的撂荒地的共同撂荒面积为600107个栅格点.植被指数变化检测法与CART决策树分类法提取的撂荒地叠加区面积显著大于与单时相波段法叠加的结果.将三种方法同时提取的结果叠加，得到研究区撂荒耕地空间分布如图7所示，由图可知，研究区撂荒耕地多分布于西北方向的新店镇、流长苗族乡等地势相对较高，坡度较大的地区.自然地理环境欠佳、交通不便以及农业生产效益低，大量农村人口外出务工成为这些地区耕地撂荒的主要原因.

4 结论

本文以我国西南地区的贵州省清镇市为例，基于一年中不同季节哨兵2号遥感数据，综合地物光谱特征、植被归一化指数(NDVI)特征，构建了基于单时相波段特征检测法、植被指数变化检测法和CART 决策树分类方法的撂荒耕地提取方法，并对提取结果进行了对比分析.形成以下结论：

(1) 撂荒耕地与正常种植耕地的光谱特征具有明显的季节性差异.撂荒耕地由于季节性作物种植和生长的差异，其光谱特征差异明显，波长为0.740～0.864 μm范围内的红边波段(b6)到近红外波段(b8)的光谱反射率高于其他地物，并且红光(b4)波段反射率最低，低于其他地物.

(2) 基于植被生长季相变化的提取方法能够很好的识别出撂荒耕地.基于不同地物光谱特征值差异的单时相波段法容易受到其他光谱值相近的地物的影响，撂荒耕地的提取精度显著低于基于植被生长季相变化的植被指数变化检测法以及综合光谱差异和NDVI变化的CART决策树分类法，且植被指数变化检测法提取精度要高于CART决策树分类法.