基于時序遙感的撂荒地監測及空間格局特征分析

肖文菊，楊穎頻,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時序變化幅度特征來區分撂荒地和非撂荒地，未來可進一步挖掘二者在空間特征、光譜特征中的差異，增加分類特征維度，進一步提高撂荒地識別精度；③在判別耕地是否撂荒的基礎上，可進一步探究季度撂荒、單年撂荒和多年撂荒的識別方法；④本文僅對撂荒地的空間分布特征進行景觀格局分析，未來可進一步擴大研究區范圍，探究耕地撂荒的驅動因素。