基于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 個林場的撂荒地以碎小分布為主，大面積撂荒現象的發生率較低，界址鎮、烏逕鎮和鄧坊鎮的撂荒率最低；雄州街道的耕地撂荒現象明顯，耕地撂荒程度較高，需引起相關部門的重點關注。