基于M-TWDTW 模型的粤北山区耕地撂荒遥感监测研究

林妙萍，杨颖频，吴志峰，徐国良，黄永芳，许志明

（1.自然资源部华南热带亚热带自然资源监测重点实验室，广东 广州 510670；2.广东省国土资源测绘院，广东 广州 510670；3.广州大学 地理科学与遥感学院，广东 广州 510006）

现阶段利用遥感技术手段对撂荒地进行监测识别的研究成果众多[1-3]，常用的方法包括基于遥感图像分类方法、基于变化监测方法和基于时间序列分析方法三大类，其中，基于遥感指数时间序列分析方法能够表征植被覆盖随时间变化的规律，有效识别农作物与撂荒地自然植被的差异，具备较强的植被生理学机理[4-5]。

时间加权动态时间规整（TWDTW）模型是近年来时间序列分析中具有代表性的方法之一[6]，可以比较两条不同的时间序列曲线变化轨迹的相似性，适用于大尺度范围，还能在一定程度上规避物候差异的影响，且样本量依赖度低，被广泛应用于作物分类的研究中[7-8]。但在撂荒地监测方面，因撂荒时长不一，植被覆盖度不同，常规TWDTW方法极易造成漏分情况。

为了解决此类局限性，本研究以耕地NDVI 全年均值为重要特征，对现有TWDTW 距离计算模型进行改进，提出基于均值校正的TWDTW 模型的撂荒地遥感监测方法。为验证该方法在撂荒地的监测识别精度和适用性，本文选取广东省韶关市南雄市作为研究区开展撂荒地监测识别提取实验，采用Sentinel-2 时间序列数据，利用大量样本对撂荒地的提取结果进行精度验证，并进行研究区的撂荒地空间制图，分析撂荒地的空间分布特征及其数量情况，在一定程度上有利于国家相关自然资源部门开展耕地信息把控、耕地资源可持续利用和耕地保护监测监管工作。

1 研究区概况与数据源

1.1 研究区概况

南雄市地处广东省北部（113°55′30″~114°44′38″E，24°56′59″~25°25′20″N），大庾岭南麓，毗邻江西。地势上大致呈南北高中间低，南北群山环绕，中部以丘陵平原为主，丘陵沿浈江延伸分布，形成狭长的“南雄盆地”。南雄市位于北回归线稍偏北的低纬地区，属于亚热带季风气候区，冬无严寒，夏无酷暑，年平均气温在17~26℃。南雄市现有耕地总面积37 454.24 hm2[9]，主要农作物有水稻、花生、烟叶，是国家和省双料“产粮大县”[10]。由于经济结构转型和农村劳动力析出，南雄市耕地撂荒现象日趋严重，因此对该地区开展耕地撂荒的遥感监测识别研究具有重要的现实意义。

1.2 遥感时间序列数据

撂荒耕地监测提取使用的遥感数据源来自高分辨率多光谱成像卫星Sentinel-2，采用遥感大数据云计算平台（GEE）获取。为了充分利用碎片化的有效数据，选取了所有整体云量小于80%的影像，共262景，涵盖了51个观测时相，能够反映撂荒地与非撂荒地在全年的光谱特征，对该数据集进行NDVI 指数计算，形成了涵盖耕地农作物生长周期的NDVI 时间序列数据。考虑到云覆盖的影响，需要对Sentinel-2 时序数据开展云掩膜、异常值剔除、NDVI 时序重建等预处理工作，具体方法见2.1小节。

1.3 无人机影像及样本数据

通过实地调研结合无人机影像目视解译的方式，在研究区共选取了945 个样本点，其中撂荒地样本445 个，非撂荒地样本500 个。主要分布在古市镇、湖口镇、油山镇、乌迳镇。按照6∶4比例，将样本集划分为训练集和测试集。训练集用于寻找撂荒地与非撂荒地的最佳距离分割阈值，测试集用于验证撂荒地识别结果精度。

1.4 土地利用数据

本研究收集了来自广东省国土资源测绘院的第三次全国土地调查数据和土地利用变更调查统计数据作为本底使用数据，获取南雄市的耕地空间分布范围。

2 研究方法

本文研究方法包含以下几个步骤：①利用遥感云平台GEE直接获取Sentinel-2 NDVI时间序列数据并进行预处理；②结合无人机获取的撂荒地样本数据集，构建撂荒地参考NDVI 时序曲线；③采用M-TWDTW方法，对撂荒地NDVI 参考曲线进行校正，并计算未知像元NDVI 时序曲线与校正后参考时序曲线的距离（相似性程度）；④通过决策树二分类方法获取最佳分割阈值，提取撂荒地，进行精度验证和撂荒地空间分布制图分析。

2.1 NDVI时序数据预处理

本研究基于GEE 直接获取Sentinel-2 数据，根据研究区范围进行裁剪，将同1 d 的影像镶嵌为一景影像，对该数据集进行NDVI 指数计算。NDVI 具体的计算公式为：NDVI=（NIR-R）/（NIR+R），其中NIR 为近红外波段，R 为红光波段。由于Sentinel-2数据在成像过程中会受到云噪声的干扰，因此基于QA 波段对NDVI 数据进行掩膜处理，利用多时相NDVI 数据构建NDVI 时序曲线后，再参照Ma[11]等判定异常值点的方法，进一步剔除时序曲线上的异常值，并通过16 dNDVI 数据最大值合成处理方法对去云和剔除后的时相点漏洞进行填补，以保持数据完整，形成涵盖耕地农作物生长周期的NDVI 时间序列数据。

2.2 撂荒地NDVI参考时序曲线构建方法

基于上述预处理后的NDVI 时序数据的曲线特征，利用训练样本集中的撂荒地样本，对每个观测时相NDVI 数据进行箱型图分析，选取其中的上下四分位数（25%~75%）区间的数据，获取各个时相的撂荒地NDVI 中位数，最终形成撂荒地NDVI 参考时序曲线。

2.3 基于M-TWDTW模型的时序曲线距离计算方法

常规TWDTW 算法的基本原理可归纳为：根据已构建的NDVI时间序列曲线数据，逐像元计算NDVI时序曲线与撂荒地NDVI 参考时序曲线的距离，并通过阈值判定像元的类型。然而，常规TWDTW 方法无法克服因撂荒时长不同，所造成的撂荒样本NDVI 时序曲线绝对值差异，极易造成漏分情况。

为了避免此类不足，本研究基于均值校正策略，对现有TWDTW 距离计算模型进行改进，提出基于M-TWDTW模型的撂荒地遥感监测方法。实现方法为：计算撂荒地NDVI参考时序曲线（TSref）的全年均值，记作meanref，计算待分类像元时序曲线（TSpixel）的全年均值，记作meanPixel，获取两者的差值meanPixel-meanref，记作△diff。以待分类像元的均值为标准，对撂荒地NDVI 参考时序曲线沿纵轴进行平移，平移距离为△diff，使得校正后的新参考曲线与待分类像元的全年均值一致，获得新的撂荒地参考曲线，记作TSref_new。计算校正后参考曲线TSref_new与TSpixel之间的TWDTW 距离。改进过程和修正后的新参考曲线如图1 所示，具体的算法过程如下：

图1 基于M-TWDTW模型的撂荒地参考时序曲线修正过程

定义An为参考时序曲线，Bn为未知像元时序曲线：

式中，ak和bk分别为An和Bn的第k个节点值，n、m为时间序列的长度，此处n与m相等，An和Bn的距离矩阵为：

式中，Sij为两时间序列的基距离，一般采用欧氏距离的平方和时间权重的乘积作为TWDTW 的基距离，即：

式中，wij为权重值；g为增益因子，值越大则对匹配点间隔差异的惩罚越大；c=|i-j|为距离因子，mc一般为时间序列的中间节点，本文设为100。

设时间序列A和B 在(i,j)处的累积距离为Q(i,j)，公式如下：

在寻找最小累积距离路径L时，应满足：

1） max{m,n}≤L≤m+n-1。

2）若Ql=Sij,Ql+1=Si'j'则0 ≤i'-i≤1,0 ≤j'-j≤1

本文基于M-TWDTW 模型，逐像元计算待分类像元NDVI时序曲线与校正后撂荒地NDVI参考时序曲线的距离值，获得的距离值越小，则相似性程度越高，代表该耕地像元为撂荒地的可能性越大。

2.4 决策树二分类模型

基于M-TWDTW模型，计算所有撂荒地样本与非撂荒地样本NDVI时序与撂荒地NDVI参考时序曲线的距离，通过CART 决策树方法不断进行二分类划分，建立有效的分类规则，获取最佳分割阈值（Dist）。构建的撂荒地监测识别规则：若距离大于阈值Dist，则判定为非撂荒地；若距离小于等于阈值Dist，则判定为撂荒地。

2.5 精度评价

分别获取基于TWDTW 模型和M-TWDTW 模型方法的分类混淆矩阵，从混淆矩阵中提取总体精度（OA）、Kappa系数、制图精度和用户精度，对模型精度进行评价。具体计算公式如下：

式中，N为像素总数；TP为真正类，表示样本的真实类和模型预测结果都是正类；TN为真负类，表示样本的真实类为正类，模型预测结果为负类；FP为假正类，表示样本的真实类为负类，模型预测结果为正类；FN为假负类，表示样本的真实类和模型预测结果都是负类。OA为真实分类像素占区域总像素的比例。Kappa 系数用于真实值和预测值的一致性检验，衡量分类精度，当Kappa 系数处于0~0.2 区间时，一致性极低；处于0.21~0.4 区间时，一致性一般；处于0.41~0.6 区间时，一致性中等；处于0.61~0.8 区间时，表示高度的一致性；处于0.81~1 区间时，表示真实值和预测值几乎完全一致。

制图精度是指分类器将整个影像的像元正确分为撂荒地类的像元数（对角线值）与撂荒地类真实参考总数（混淆矩阵中撂荒地类列的总和）的比率；用户精度是指正确分到撂荒地类的像元总数（对角线值）与分类器将整个影像的像元分为撂荒地类的像元总数（混淆矩阵中撂荒地类行的总和）比率。

3 结果与讨论

3.1 撂荒/典型作物的NDVI时序曲线

图2 展示了典型撂荒地样本NDVI 时序曲线及相应的无人机影像，图3 为不同轮作制度下农作物类型的NDVI 时序曲线图。由图2 可知，撂荒地的NDVI 时序曲线在全年时相上的波动变化较平缓，未呈现出剧烈波动的作物生长物候特征。图3a 中撂荒地植被稀疏，可以看出其撂荒时间较短，NDVI 最大值约在0.5 左右，时序曲线振动幅度较小；图3b 中撂荒地植被密集，NDVI 最大值已接近0.8，属于撂荒时间较长的撂荒地，较图3a 的时序曲线振动幅度稍大。

图2 典型撂荒地样本的NDVI时序曲线及无人机影像

图3 不同作物类型的NDVI时序曲线图

从图3的农作物类型的NDVI时序曲线图来看，双季稻的NDVI时序曲线特征明显呈现出2个波峰，分别在4~7月和7~10月；红薯的生长发育周期有发根缓苗期、分枝结薯期、茎叶盛长期和茎叶衰退薯块膨大期等四个阶段，一般在6月进行发根缓苗，10月茎叶逐渐衰退，薯块迅速膨大，进入收获期；烟叶从移栽到成熟收获一般历时5~6个月，从2月底开始移栽，到6月初收获，其余时间种植其他作物；对于单季作物的花生而言，春季的4 月是最佳种植时间，生长至8 月则到采收期。非撂荒地的NDVI 时序曲线变化趋势与作物生长物候阶段相对同步，随着作物生长发育、成熟衰老呈现出先上升后下降的变化特征。总体而言，撂荒地NDVI时序曲线与非撂荒地NDVI时序曲线无论是在波峰数量上还是在曲线振动幅度上都有着较大的差异。

3.2 撂荒地NDVI参考时序曲线

如图4 所示，本研究基于撂荒地样本构建了撂荒地全年的NDVI 参考时序曲线。从整体来看，撂荒地NDVI 参考时序曲线无明显的作物生长物候特征，曲线在全年的波动变化、起伏相对稳定。撂荒地参考时序曲线的具体变化可归纳为4 个阶段：第一阶段的NDVI 值从1 月开端到2 月中下旬大约50 d内呈下降变化趋势，且在2 月中下旬达到最低值0.2，说明在冬季时期撂荒耕地中的荒草可能出现枯萎现象，因此绿色植被减少导致植被指数降低；在第二阶段，NDVI 数值紧接着约在一个月内（3 月）迅速上升至0.5，早春来临，使得绿草迅速茂密生长，直接反映在NDVI 数值的迅速上升中；第三阶段NDVI 数值保持小幅度波动至11 月初，撂荒耕地被杂草覆盖，没有农作物的明显物候特征，NDVI 时序曲线变化平缓；最后是第四阶段的冬季，NDVI 数值逐渐下降。

图4 撂荒地NDVI参考时序曲线

3.3 撂荒地识别结果与精度对比

本研究采用撂荒地和非撂荒地样本验证M-TWDTW模型对撂荒地的识别精度，并与常规TWDTW模型方法进行对比。计算样本集中各个样本与撂荒地参考时序曲线的距离，统计撂荒地样本与非撂荒地样本的距离直方图（图5）。

图5 撂荒地与非撂荒地样本NDVI时序与撂荒地NDVI参考时序的距离直方图

图5a 基于TWDTW 方法的距离中，撂荒地样本的距离计算结果分散分布在多个区间，与非撂荒地的距离计算结果直方图存在较多重叠区域，不利于阈值的分割，易导致较大的精度误差；相比之下，图5b 基于M-TWDTW 方法的距离集中分布在连续的区间，样本数量更接近正态分布，撂荒地与非撂荒地的距离直方图重叠区域较少，有利于选取最佳分割阈值。

基于以上直方图，采用决策树二分类方法，获取撂荒地与非撂荒地的时序距离最佳阈值，结果显示：基于TWDTW模型方法的距离分割阈值为1.948 1；经过均值校正后，利用M-TWDTW模型的距离分割阈值为1.769 9。本文利用实地调研和无人机影像目视解译结合方式获取的撂荒地与非撂荒地样本，分别对TWDTW 模型和M-TWDTW 的识别精度进行验证，结果分别如表2、3所示。结果表明，当采用M-TWDTW模型对撂荒地进行监测识别时，其总体精度（87.04%）和Kappa 系数（0.74）比采用TWDTW 模型识别撂荒地的总体精度（76.72%） 和Kappa 系数（0.53）更优，撂荒地的制图精度为90.45%，用户精度为83.42%，非撂荒地的制图精度为84.00%，用户精度为90.81%。深入分析可知，TWDTW 模型仅能比较不同的植被时序曲线变化轨迹绝对值的距离相似性，未考虑撂荒时长不同的撂荒样本NDVI 时序曲线绝对值的差异，在一定程度上出现漏分和误分，因而TWDTW 模型方法识别撂荒地的精度偏低，对于撂荒时长不一的地区的适用性不强； 而基于M-TWDTW 模型方法，通过对撂荒地NDVI 参考时序曲线均值的校正，重点捕捉撂荒地的NDVI 时序变化形态特征，弱化了撂荒时长不一引起的NDVI 绝对值差异的影响，避免因撂荒时长不一引起的漏分或误分等误差，因此识别精度较高，说明在不同撂荒时长的耕地撂荒监测识别中，M-TWDTW 模型方法适用性更强。

表2 基于TWDTW模型的精度验证混淆矩阵表

表3 基于M-TWDTW模型的精度验证混淆矩阵表

3.4 撂荒地空间分布制图

基于上述方法进行南雄市耕地撂荒遥感监测，获得南雄市撂荒地空间分布制图结果（图6）。根据实验结果数据，以图斑为单元进行统计，撂荒耕地平均图斑面积为2 200.83 m2，非撂荒耕地平均图斑面积为2 456.01 m2。在空间分布特征上，南雄市撂荒地集中分布在中部区域，存在大面积撂荒现象。南雄市各镇/街道中，雄州街道、湖口镇和黄坑镇的撂荒地平均图斑面积较大，大面积撂荒现象的发生率较高；由于3 个林场（山门林场、帽子峰林场和泷头林场）的耕地面积较少，撂荒地以碎小分布特征为主，大面积撂荒现象的发生率较低，特别是山门林场，撂荒地平均面积最小。结合图6 撂荒地与非撂荒地空间分布和耕地面积统计结果得知，乌迳镇、湖口镇和油山镇的耕地面积最大，其撂荒地总面积也最大。从撂荒率来看，泷头林场、帽子峰林和山门林场的撂荒率最高，界址镇的撂荒率最低，其次是乌迳镇和邓坊镇。值得注意的是，雄州街道不仅撂荒平均面积最大，撂荒发生率也高，这说明雄州街道的耕地撂荒现象明显，耕地撂荒程度较高，相关部门需对其重点关注。

图6 南雄市2021年撂荒地空间分布图

4 结 论

本研究基于均值校正策略，提出基于M-TWDTW模型的撂荒地遥感监测方法，该方法通过重点捕捉撂荒地的NDVI 时序变化形态特征，减少了因撂荒时长不一造成的漏分和误分，在耕地撂荒监测识别方面具有较强的鲁棒性。综合采用Sentinel-2 时间序列数据，结合实地调研与无人机影像目视解译的方式，利用样本数据开展实验和精度验证工作，最后进行研究区的撂荒地空间制图，分析撂荒地的空间分布特征及其数量情况。以下为本研究的主要结论：

1）撂荒地的NDVI时序曲线在全年时相上的波动变化较平缓，无剧烈波动的作物生长物候特征；非撂荒地的NDVI 时序曲线变化趋势与作物生长物候阶段相对同步；撂荒地与非撂荒地NDVI 时序曲线在波峰数量和曲线振动幅度上都有着较大的差异。

2）撂荒地NDVI 参考时序曲线变化在一年中大致分为4 个阶段：第一阶段因冬季荒草枯萎，NDVI值呈下降趋势；在第二阶段早春时期绿草迅速生长，NDVI 数值迅速上升；第三阶段撂荒耕地被杂草覆盖，NDVI 数值保持小幅度波动；第四阶段NDVI数值逐渐下降。

3）与TWDTW 方法相比，基于M-TWDTW 方法的撂荒地样本在距离直方图中的数量分布更接近正态分布，与非撂荒地的重叠区域更少，更利于选取最佳分割阈值；M-TWDTW 方法弱化了撂荒时长不一引起的NDVI 绝对值差异的影响，重点捕捉NDVI时序曲线的形态特征，对耕地撂荒监测识别精度更高；基于改进的M-TWDTW 模型的撂荒地制图精度为90.45%，远高于基于常规TWDTW 模型的撂荒地制图精度66.85%。

4）2021 年南雄市撂荒地空间分布集中在中部区域，存在大面积撂荒现象；3 个林场的撂荒地以碎小分布为主，大面积撂荒现象的发生率较低，界址镇、乌迳镇和邓坊镇的撂荒率最低；雄州街道的耕地撂荒现象明显，耕地撂荒程度较高，需引起相关部门的重点关注。