基于时序遥感的撂荒地监测及空间格局特征分析

肖文菊，杨颖频,2，吴志峰,3，郑少兰

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

随着城市化进程不断加快伴随着经济结构的调整，农村人口向城市转移，耕地撂荒现象日益严峻，给国家粮食生产安全造成严重威胁（Li et al.,2012; Yusoff and Muharam, 2015; Ustaoglu and Collier, 2015; Goga et al., 2019）。准确高效获取撂荒地的空间分布现状是防止耕地撂荒趋势进一步扩散的前提和基础（牛继强 等，2017；杨通，2020；Meijninger et al., 2022）。传统地面调查的方式效率低、周期长，遥感技术具有高效率、低成本对地观测的优势，为撂荒空间分布范围及面积监测提供了更为便 捷、有 效 的 手 段（Witmer, 2008; Watanabe and Saiga, 2009; Alcantara et al., 2012; Estel et al., 2015；陈欣怡 等，2018）。

当前，基于遥感技术的撂荒地监测方法可归纳为3类：1）监督分类方法，该方法主要基于遥感影像的空间、时间、光谱等多维特征，利用支持向量机（Li et al., 2012；牛继强 等，2017）、随机森林（Wu et al., 2020；张昊 等，2022）等机器学习算法，进行撂荒地识别。如张昊等（2022）基于地物的空间和光谱特征，采用随机森林分类方法提取青海省民和县的撂荒耕地分布信息，得到2018、2019、2020 年的撂荒地提取精度分别为86.9%、87.36%、85.92%。2）多时相范围叠加法，该方法通过对比2个或多个时相的土地利用数据提取耕地撂荒信息，多用于耕地撂荒情况的年际变化检测，且提取结果的可靠性主要依赖土地利用数据的精度（郑财贵 等，2010；史铁丑 等，2016）。如史铁丑等（2016）将重庆市2 期耕地地块空间范围叠加，在剔除退耕还林及森林工程图斑后，将剩余部分作为撂荒耕地。3）遥感时序分析方法，该方法通过分析撂荒地和非撂荒地在时间序列上的变化特征实现撂荒地的识别（Yusoff and Muharam, 2015; Dara et al., 2018; Xu et al., 2018；王红岩 等，2020；Wei et al., 2021；宋 宪 强 等，2021）。如 宋 宪 强 等（2021）根据春、夏、秋3个时相的耕地NDVI差值设置撂荒和非撂荒的分割阈值提取撂荒地，在四川凉山开展耕地撂荒监测；王玲玉等（2020）基于NDVI 时序数据，通过设定NDVI 时序峰值的取值范围提取撂荒地，并在贵州省息烽县开展撂荒监测。基于时间序列分析的撂荒监测方法，与一段时间内耕地作物与撂荒地中自然植被的物候差异紧密相关，具有较强的植被生理学含义。

中国华南地区水热条件丰富，作物种植类型多样、熟制结构复杂，耕地撂荒后自然演替生长的植被长势状态良好，进一步提高了种植耕地与撂荒地的区分难度。一方面，华南地区作物种植结构复杂，每种作物的物候期存在差别，冬季气候温和，在冬季种植瓜果、蔬菜等作物的情况十分普遍，因此，仅利用春、夏、秋季的某几个观测时相进行撂荒识别，很可能会造成误判，有必要将观测时间窗口扩宽到全年；另一方面，不同地块撂荒时间不尽相同，植被覆盖度存在差异，NDVI数值有所差别，如撂荒时间较短的地块杂草稀疏，NDVI水平较低，而撂荒时间较长的地块杂草生长繁茂，NDVI 呈现与作物生长旺季时相当的水平，因此，基于NDVI峰值的撂荒地识别方法在该地区也难以适用。而NDVI 时序振幅特征，即NDVI 时序最大值与最小值的差值，体现一段时间内耕地内部植被生长状态变化、生长发育速率等特征与植被生理变化状态紧密相关。

在遥感数据源方面，当前撂荒监测研究大多基于Landsat、MODIS 等中低分辨率影像（Witmer,2008; Alcantara et al., 2012; Estel et al., 2015；牛继强 等，2017；Dara et al., 2018; Wei et al., 2021），中低分辨率遥感数据能提供高时间分辨率的监测能力，捕捉耕地的季相变化特征，但受限于空间分辨率，在地表异质性高的地区难以适用（杨通 等，2019）。Sentinel-2数据遥感兼具较高的时空分辨率，可见光、近红外波段空间分辨率为10 m，显著提升了复杂地表区域的观测能力，双星观测时间分辨率为5 d，在耕地季相变化的动态观测方面具有较大的优势。

撂荒地空间分布格局能反映撂荒耕地在空间位置、空间形态等方面的特征，是探究耕地撂荒原因的重要手段。当前，相关研究多基于空间分析方法提取撂荒地空间分布格局特征，通过缓冲分析、密度分析、相关分析等方法提取撂荒地的空间位置、空间形态及空间关系特征。如牛继强等（2017）对罗山县子路镇的撂荒地进行了缓冲区分析，发现灌溉条件是影响主要因素；刘智丽（2020）对晋中祁县不同年份撂荒地空间分布格局变化进行对比分析，发现撂荒地在空间分布上呈现面积减小、斑块密度不断增加的趋势；董世杰等（2023）对中国撂荒梯田的空间格局特征的分异性进行了探究，发现撂荒梯田呈现“北低南高”的特征，南方地区山地丘陵地区撂荒严重；唐瑞等（2022）分析了阆中市撂荒地空间格局分异规律，发现低山区耕地撂荒率明显高于丘陵地区。计算撂荒地景观格局指数可为撂荒地空间分布特征分析提供基础。通过景观格局指数大小可分析撂荒地和非撂荒地的空间结构差异，为进一步探究耕地撂荒驱动因素提供数据基础。

鉴于此，本研究将撂荒时间≥1 a的耕地定义为撂荒地，提出了一种基于NDVI时序振幅的撂荒地识别方法，通过统计撂荒地与非撂荒耕地的NDVI振幅取值分布，划定NDVI 振幅的最佳分割阈值，从而构建撂荒地的识别规则，并以广东省湛江市坡头区为例，验证该方法的适用性。并在撂荒地空间分布制图的基础上，进一步探究撂荒地空间分布格局特征。以期通过扩充监测撂荒地的时间信息，捕捉植被在监测窗口内NDVI最大变化，获得较高的撂荒监测精度。

1 研究区域和数据源

1.1 研究区域

湛江市作为广东省农业生产主要基地，农业经济占有重要地位。其中坡头区是湛江市重要的农产品来源地。坡头区位于广东省西南部地区，雷州半岛东北部，湛江海湾东部，由一个半岛和一个海岛组成，（图1-a）。该区地处北回归线以南的低纬地区，属于热带季风气候，年均温在22.7～23.5℃。年均雨量1 395.5～1 723.1 mm，年均日照时数1 714.8～2 038.2 h，半岛地势平坦，整体耕种条件较好。该区总面积424 km2，其中耕地面积为145.18 km2，图1-b 为坡头区耕地分布情况，其种植类型以水稻、花生、红薯为主。随着经济结构的转型和劳动人口流失，坡头区耕地撂荒现象加剧，对该地区进行耕地撂荒监测具有重要的现实意义。

图1 研究区地理位置Fig.1 Geographical location of study area

1.2 数据源

1.2.1 遥感影像及耕地数据 基于地理计算云平台Google Earth Engine（GEE）共获取Sentinel-2 数据共12景（该数据已经过辐射定标、大气校正等预处理过程），计算NDVI指数，获取耕地NDVI时间序列数据。影像成像时间如表1所示，监测时间窗口涵盖作物生长的完整周期。

表1 Sentinel-2时间序列数据Table1 Time series data of Sentinel-2

以第三次全国土地调查数据中的耕地范围为基础，对坡头区大片农田中不同种植类型的耕地地块进行人工划分，得到耕地地块数据。

1.2.2 样本数据 撂荒地及非撂荒地样本主要来源于野外地面调查和人工目视解译，其中，基于实地调查获得撂荒地样本共200个，非撂荒地样本共200个。结合实地调查数据，观察撂荒地和非撂荒地在多时相Google Earth高分辨率影像上的空间特征，补充样本数量并构建撂荒地和非撂荒地样本集。如图2-a所示，因受到强烈的农业生产活动干预，种植耕地地块形态规整，内部纹理规则，因作物生长物候变化，地块内覆被变化明显。而撂荒地因无人工干预，发生撂荒后自然植被演替，纹理杂乱、边界不清晰，在时间维度上无明显物候特征（图2-b）。构建的样本集中共包含撂荒地和非撂荒地样本各300 个，按照2∶1 的比例随机分为2部分，分别用于构建撂荒地识别规则和验证精度，图3为精度验证样本的空间分布。

图2 不同时相撂荒地与非撂荒地影像对比Fig.2 Image comparison of abandoned land and non-abandoned land in different phases

图3 研究区地面验证点空间分布Fig.3 Ground survey points in research area

2 撂荒地识别及空间格局分析

本研究框架为：1）耕地NDVI 时序数据预处理；2）耕地NDVI时序振幅特征提取；3）撂荒地识别规则构建；4） 撂荒地空间格局特征分析（图4）。

图4 撂荒地监测技术流程Fig.4 Monitoring of abandoned land technical flowchart

2.1 耕地NDVI时序数据预处理

利用多时相Sentinel-2数据构建耕地的NDVI时序曲线。由于光学影像在成像过程中受到云的干扰，时间序列曲线通常包含一定“噪声”，为避免其对振幅统计结果的干扰，需将NDVI时序曲线上的异常值剔除。异常值剔除方法参照Ma和Veroustraete（2006）的方法，将短时间内急剧下降再上升的点判定为异常值点。结合试验区采集的样本点，分析其NDVI曲线变化趋势，经过多次实验，设计异常值点的判定方案：设NDVI 时序曲线上相邻3点的时间分别为Ti-1、Ti、Ti+1，其对应的NDVI 值分别为Vi-1、Vi、Vi+1。当符合以下2类情况之一：1）若Ti-Ti-1≤14 d，且Vi-1-Vi与Vi+1-Vi均＞0.2；2）若Ti-Ti-1≤21 d，且Vi-1-Vi与Vi+1-Vi均＞0.5，认定Ti为异常值点，将其从曲线上剔除，并利用前后邻近观测点的值对Ti处进行线性插补。基于预处理后的NDVI 时序曲线提取振幅特征，即利用NDVI时序曲线上的最大值减去最小值，用以表征耕地内全年的植被覆盖变化。

2.2 撂荒地识别规则构建

基于撂荒地与非撂荒地样本的振幅特征提取结果，划定用于撂荒地识别的最佳振幅阈值。阈值设定方法主要包括2步：1）分割阈值初始化：将所有撂荒地样本的最大振幅设定为初始阈值，记为Threshold0，表示在该阈值下所有撂荒地样本均能被正确识别，但可能存在部分非撂荒地被错误地判别为撂荒地。2）分割阈值最优化：在Threshold0基础上，以0.01 为步长，不断减小分割阈值，依据F1指数（全国科学技术名词审定委员会 等，2002）评价阈值可靠性，确定最佳分割阈值。F1指数是评价二分类模型中分类准确性的重要指标，可看作是模型精准率（precision）和召回率（recall）的一种加权平均，兼顾撂荒地与非撂荒地的识别精度，F1最大值为1，最小值为0。当F1越大时，分类效果越好，对应的阈值记作Thresholdoptimum。

式中：precision表示被识别为撂荒地的样本中，实际为撂荒地的概率，即模型精准率；recall表示在所有撂荒地样本中，被正确识别为撂荒地的概率，即召回率。

在选取最优的分割阈值后，构建的撂荒地识别规则为：当耕地NDVI 时序振幅≥Thresholdoptimum时，将其判定为非撂荒地；当耕地NDVI 时序振幅＜Thresholdoptimum时，将其判定为撂荒地。

2.3 撂荒地空间格局特征分析

计算撂荒地景观格局指数，分析坡头区撂荒地的空间格局特征，计算撂荒地景观格局指数，包括图斑总数（NP）、图斑总面积（TA）、平均图斑面积（MPS）、平均图斑形状指数（MSI）、平均图斑分维指数（MPFD）、聚集度指数（AI），分析撂荒地图斑在空间分布上的空间形态特征及空间集散程度。各指数计算方法为：

1）平均图斑面积（MPS）：

式中：ai为第i个图斑的面积；NP 为图斑个数，反映图斑的破碎化程度。MPS 值越小，说明斑块越破碎。

2）平均图斑形状指数（MSI）：

式中：Pi为图斑i的周长；MSI反映景观要素图斑的规则程度。正方形MSI 取值为1，MSI 值越接近1，说明图斑形状越规则。

3）平均图斑分维指数（MPFD）：

MPFD 可以度量图斑边界的复杂程度，MPFD值越大，说明图斑形态越不规则，反之图斑形状越规则。

4）聚集度指数（AI）：

式中：gii为相应景观类型的相似邻接斑块数量。AI可反映某一类景观斑块之间的连通度。AI值越大，说明该类型景观分布越密集，反之则越分散。

3 结果与讨论

3.1 NDVI时序曲线特征对比

坡头区主要作物类型包括水稻、花生、蔬菜和薯类。图5 展示了坡头区撂荒地与非撂荒地NDVI时序曲线。撂荒地在全年内NDVI 变化相对平稳，无明显作物生长物候特征。非撂荒地在作物生长窗口期表现出明显的波动变化，从图5所示的案例曲线看，水稻为两季作物，早稻生长季为4—7月，晚稻生长季为8—11 月，其余时间种植薯类。花生为单季作物，生长季为4—8 月，9 月至次年3 月种植蔬菜。红薯一般为6 月扎根缓苗，10 月上旬收获，11月直至次年5月以蔬菜种植为主。可见，撂荒地与非撂荒地在NDVI振幅特征方面存在较大差异。

图5 2020年不同轮作制度下作物NDVI时序曲线Fig.5 Crop NDVI time series curves under different crop rotation systems in 2020

3.2 振幅阈值确定

利用200 个撂荒地与200 个非撂荒地样本的NDVI 时序数据计算振幅特征，绘制振幅分布（图6）。

图6 NDVI时序振幅分布Fig.6 Histogram of amplitude of NDVI time series

撂荒地样本的NDVI振幅整体偏低，取值范围在0.16～0.545，非撂荒地样本的NDVI 振幅取值范围整体在0.29～0.744。如图6 所示，撂荒地与非撂荒地的NDVI振幅在取值范围上在整体上呈对称分布特征，撂荒地的NDVI振幅主要分布在0.1～0.4之间，非撂荒地的NDVI振幅主要分布在0.4～0.7，但二者的NDVI 振幅分布区间存在一定的交叉重叠。为了实现撂荒地与非撂荒地整体识别精度的最大化，依据F1指数动态设定NDVI振幅阈值，以撂荒地样本的最大振幅为起始阈值，非撂荒地样本的最小振幅为终止阈值，计算不同阈值下对应的F1指数。不同NDVI 振幅取值的F1指数如图7 所示。当NDVI 振幅为0.42 时，F1指数最高，F1为0.91，代表撂荒地与非撂荒地的二分类精度最高。此时，撂荒地样本的识别精度为91.83%，非撂荒地样本的识别精度为90.20%。

图7 不同阈值下F1计算结果Fig.7 F1 calculation results under different thresholds

3.3 撂荒地提取结果及精度验证

基于上述方法进行坡头区耕地撂荒遥感监测，获得坡头区撂荒地空间分布制图结果（图8）。经统计，2020 年坡头区撂荒斑块总数为1 514 个，撂荒耕地面积达14.65 km2，占总体耕地面积的10.1%。

图8 湛江市坡头区2020年撂荒地空间分布Fig.8 Spatial distribution of abandoned land of the Potou County in 2020

利用100 个撂荒地和100 个非撂荒地样本对遥感监测结果进行精度验证，混淆矩阵计算结果如表2 所示。识别的总体精度为91%，Kappa 系数为0.82。撂荒地的生产者精度为91.83%，非撂荒地生产者精度为90.20%。这说明基于NDVI振幅阈值分割方法能实现较高精度的撂荒地识别。

表2 基于NDVI振幅特征的精度验证混淆矩阵Table 2 Precision verification confusion matrix table based on NDVI amplitude characteristics

3.4 对比实验结果及精度验证

基于实验样本数据，分别采用多时相NDVI差值（宋宪强 等，2021）和NDVI峰值方法（王玲玉等，2020）进行撂荒识别。其中，基于多时相NDVI差值方法中，通过200个撂荒地实验样本与对应时相的NDVI进行叠加与相减，得出不同时相点撂荒地和非撂荒地的NDVI 差值，通过对NDVI 差值进行统计分析，得到春、夏2 季NDVI 差值最小为0.39，夏、秋2季NDVI差值最小为0.38。为保证所有撂荒样本全部落入撂荒区域，确定最佳分割阈值为0.39，当夏-春、夏-秋的NDVI差值均＜0.39时，即认定其为撂荒地，反之则为非撂荒地。基于全年NDVI 峰值提取实验中，对200 个撂荒地实验样本的NDVI 峰值进行统计分析，得到最大NDVI 峰值为0.68，即当植被生长期的NDVI 最大峰值＜0.68时为撂荒地，＞0.68时为非撂荒地。

对上述2种方法进行精度验证，采用与3.3小节中相同的验证样本，分别验证基于多时相NDVI差值方法和基于全年NDVI峰值方法的识别精度。结果发现，当利用多时相NDVI差值进行撂荒识别时，总体精度为84.00%，Kappa系数为0.68，撂荒地的生产者精度为86.17%，非撂荒地的生产者精度为82.08%（表3）。基于全年NDVI峰值特征识别的总体精度为74.00%，Kappa系数为0.48，其中撂荒地的生产者精度为77.91%，非撂荒地的生产者精度为71.05%。

表3 基于多时相NDVI差值和全年NDVI峰值的精度验证混淆矩阵Table 3 Precision verification confusion matrix based on multi temporal NDVI difference and NDVI peak value

从精度验证结果可看出，基于NDVI峰值方法识别撂荒地的精度最低，基于多时相NDVI差值提取方法精度有所提高，基于NDVI振幅特征提取方法精度最高。NDVI 峰值特征提取方法仅利用单个特征值作为判别依据，而不同地块撂荒时间不同，植被覆盖度存在差异，NDVI 峰值也各不相同，通过该方法无法准确识别。如一部分撂荒时间长的耕地其NDVI峰值较高，借助该方法易产生漏判。此外，一部分耕地种植类型为多年生树苗或果苗，当树苗稀疏、土壤裸露程度较高时，其NDVI水平较低，容易将种植耕地误判为撂荒地。基于多时相NDVI 监测方法利用春-夏-秋3 个时相点的植被信息，虽在一定程度上扩充了时间维度信息，但对于种植结构复杂、物候特征多样的种植区，该方法存在错判、漏判的可能性较高。如许多蔬菜为季节性作物，其种植和收获时间存在差异，仅通过春-夏-秋3个时相点无法准确捕捉其生长信息，易将部分蔬菜种植区误判为撂荒地。相较基于NDVI峰值和多时相NDVI 差值方法，基于NDVI 振幅特征的提取方法能捕捉植被全年生长变化情况，在一定程度上扩充了时间维度信息，同时结合了NDVI的最大值和最小值，可获得较高提取精度。

3.5 撂荒地空间分布格局分析

研究区的撂荒地空间格局指数结果如表4所示。对比撂荒地与非撂荒地的空间格局特征，撂荒地斑块为1 514个，约占耕地总图斑数量的19.8%；撂荒耕地平均图斑面积为0.97 hm2，非撂荒耕地平均图斑面积为2.13 hm2，撂荒地的平均图斑面积远小于非撂荒地，可看出坡头区撂荒地分布破碎，大面积撂荒现象较少。从平均图斑形状指数（MSI）看，撂荒地的形状偏离正方形的程度更高，说明坡头区撂荒耕地地块的形状较不规则。从平均图斑分维指数（MPFD）看，撂荒地的MPFD值高于非撂荒地，说明坡头区撂荒的耕地斑块多为形状、边界不规则的地块。从聚集度指数（AI）看，撂荒地斑块的分布聚集度值更小，相较于非撂荒地，撂荒地的空间分布更加零散。

表4 景观格局指数计算结果Table 4 Calculation results of landscape pattern indexes

4 结论

本文提出了基于光学时序振幅特征的耕地撂荒识别方法，根据撂荒地与非撂荒地NDVI时序振幅特征的差异，开展耕地撂荒监测，该方法具有较强的植被生理学含义。通过Sentinel-2 数据集构建耕地时间序列曲线，Sentinel-2 兼具较高的时空分辨率，在地表异质性较强的区域具有较好的优势。主要结论如下：

1）基于NDVI时序振幅特征能很好地体现撂荒地与非撂荒地的植被生长变化差异，为耕地撂荒监测提供有效特征，非撂荒地由于作物生长物候特征信息，其光学曲线形态起伏明显，振幅较大。而撂荒地在无人工干预条件下，在作物生长窗口期无农作物生长信息，其光学曲线形态平缓，振幅较小。

2）在撂荒地识别规则构建方面，采用F1指数衡量撂荒地与非撂荒地二分类精度，通过迭代方法动态设定的振幅阈值，以获取最佳分割阈值，该方法能有效识别撂荒地与非撂荒地，在广东省湛江市坡头区开展精度验证试验，撂荒地识别精度达91.83%，非撂荒地识别精度为90.20%。

3）将本研究方法应用于坡头区撂荒地分布制图，分析撂荒地的景观格局特征。坡头区撂荒耕地面积占总耕地面积的10.1%，撂荒地平均斑块面积为0.97 hm2，撂荒地块空间形态普遍不规则，空间分布零散、少有大面积聚集，呈现破碎化的空间分布格局。

本研究方法需要以获得较高频次的光学遥感数据为基础，可应用于种植结构复杂，轮作模式多样的区域。未来还可从以下方面进一步探究：①受限于云雨等天气因素影响，难以获得更加密集的光学时序影像，给振幅的提取带来一定的不确定性，未来可进一步结合多源光学遥感数据，协同光学和SAR时序数据，通过时空融合方法构建更高频次的时序观测数据集，为准确识别撂荒地提供数据基础；②本文主要利用了耕地的NDVI时序变化幅度特征来区分撂荒地和非撂荒地，未来可进一步挖掘二者在空间特征、光谱特征中的差异，增加分类特征维度，进一步提高撂荒地识别精度；③在判别耕地是否撂荒的基础上，可进一步探究季度撂荒、单年撂荒和多年撂荒的识别方法；④本文仅对撂荒地的空间分布特征进行景观格局分析，未来可进一步扩大研究区范围，探究耕地撂荒的驱动因素。