基于聯合變化檢測的耕地撂荒信息提取與驅動因素分析

楊 通 郭旭東 岳德鵬 汪曉帆 韓圣其

(1.北京林業大學精準林業北京市重點實驗室， 北京 100083；2.中國土地勘測規劃院自然資源部土地利用重點實驗室， 北京 100035)

0 引言

現階段撂荒地提取方法主要分為：①基于遙感圖像的分類方法，如目視解譯[1-2]、監督分類[3-8]、面向對象分類[9-10]等。②基于變化檢測的方法，如直接變化檢測[5]、分類后變化檢測[11-13]、植被指數變化檢測[14]、非遙感數據(主要為GIS數據)的變化檢測等[15]。

分類方法的前提是基于撂荒地與其他地類存在顯著差別，但撂荒地的地表覆被復雜，樣本可分離度低，容易與裸地、草地、灌木等地類混分。變化檢測方法的前提是基于耕地在監測周期內發生撂荒，但對于監測周期之前已經發生的撂荒地無能為力，且容易受非耕地變化噪聲影響。此外，遙感數據本身也存在一定限制，中低分數據空間分辨率低，不足以提取撂荒地；高分數據重訪周期長、幅寬小，難以保證大區域的時空覆蓋，且易受地形、云層干擾。再加上撂荒地存在多種類型(完全撂荒、半撂荒、過渡撂荒)[1]，不同區域的作物特征和種植制度也有差別，因此，用傳統方法提取撂荒地困難。

本文提出基于多源數據的聯合變化檢測方法提取撂荒地。首先利用多源高分遙感數據的觀測頻度優勢，保證研究區的時空覆蓋和晴空數據比重；其次利用地理信息調查數據圈定耕地范圍和退耕還林范圍，以減少非耕地變化噪聲干擾；最后，針對不同類型的撂荒地分別提出年內檢測法和年際檢測法，并耦合為聯合變化檢測法。利用該方法提取撂荒地的距離特征、高差特征、鄰域特征及灌溉特征，并對上述特征進行統計分析及顯著性分析，以期得到區域撂荒主導因素，為定向提升撂荒地管理提供依據。

1 數據與方法

1.1 研究區概況

圖1 研究區地理位置圖Fig.1 Geographical location map of research area

1.2 數據來源

遙感數據源包括高分一號衛星(GF-1)8、16 m多光譜影像；高分二號衛星(GF-2)4 m多光譜影像；資源三號衛星(ZY-3)6 m多光譜影像；哨兵二號衛星(Sentinel-2A)10 m多光譜影像。獲取2014—2017年每年春(5月)、夏(7—8月)、秋(9月下旬—10月上旬)多源遙感影像106景，依據“無云數據優先、高分數據優先、高可靠性數據優先”原則篩選可用影像63景，經處理得到覆蓋全縣的12期高分辨率遙感影像數據，其中，國產衛星數據來源于中國土地勘測規劃院，哨兵數據來源于歐空局(ESA)官網，有效載荷技術指標見表1。

表1 衛星有效載荷技術指標Tab.1 Satellite payload technical index

遙感數據源共計4種衛星、5種傳感器，各傳感器均具備可見光到近紅外波段，且各波段中心波長一致，具備數據融合的基礎。地理信息數據包括全國土地利用變更調查數據、退耕還林數據。輔助數據集包括GDEMDEM的30 m DEM數據、設施占地數據、統計數據及問卷調查數據。其中，土地變更調查數據來源于中國土地勘測規劃院，退耕還林數據來源于和林格爾縣林業局，設施占地數據來源于地方網站，DEM數據來源于地理空間數據云。數據源時相分布見表2。

1.3 原理與方法

變化檢測方法是利用地物在時間維度的變化情況，回避了易混地類樣本可分離度低的問題，且該方法受地形影響較小，相對適宜提取撂荒地。但是該方法不能提取監測周期以前發生的撂荒地，也易受非耕地變化因素干擾。為解決這兩個問題，首先需要將ROI定位在耕地范圍，然后在ROI內提取正在種植耕地，兩者相減可得未種植耕地，即撂荒地。該方法有賴于對農作物的精確提取，需對研究區的主要農作物類型、特征、種植制度予以統計分析，在此基礎上確定數據時相和變化閾值。此外，利用多源數據的高時空分辨率和高晴空數據比，可有效提高監測精度。最后，對不同類型的撂荒地需采取不同的提取策略。參考歐盟環境政策研究所關于撂荒地的定義[1]，并結合研究區的實際情況，本研究將連續2年及2年以上未種植的耕地歸入完全撂荒；將一年未耕種的耕地歸入半撂荒；將退耕還林、設施占地造成的短期撂荒歸入過渡撂荒。用年內檢測法提取完全撂荒，用年際檢測法提取半撂荒，用退耕還林、設施占地數據提取過渡撂荒，對上述方法進行耦合分析，聯合變化檢測撂荒地提取方法技術體系如圖2所示。

表2 數據源時相分布Tab.2 Data source phase distribution

圖2 技術路線圖Fig.2 Technical roadmap

1.3.1多源數據預處理

所選4種高分衛星在所需波段上的中心波長一致，具備數據融合的基礎。通過野外踏勘建立聯合平差控制點，用七參數法將所有坐標系統轉換為WGS84坐標系統。根據傳感器的定標方程和定標系數，將其記錄的量化DN值轉換成對應視場的表觀輻射亮度，再經大氣校正獲得地表反射率。使用ATCOR軟件實現GF-1、GF-2、ZY-3等國產數據的大氣校正；使用Sen2Cor軟件實現Sentinel-2A的大氣校正。利用ENVI5.3軟件對所有GF-1、GF-2、ZY-3等國產數據進行幾何校正。Sentinel-2A的L1C級產品為幾何精校級，可作為參考影像校正其他數據。

多源數據融合包括高精度幾何相對校正與相對輻射歸一化。其中高精度幾何相對校正以Sentinel-2A數據為基準自動選擇控制點以優化國產數據的RPC參數，并結合DEM實現高精度幾何相對校正。相對輻射歸一化:在輻射定標和大氣校正基礎上，將所有數據重采樣到10 m空間分辨率，升采樣采用二次線性插值模型，降采樣采用像素聚合模型。依據重疊區域直方圖匹配方法將同組數據進行拼接、勻色和羽化，增強同組數據的輻射一致性。

1.3.2主要農作物生長特征統計

撂荒地提取精度依賴于種植耕地的提取精度。樣本縣主要農作物為玉米、馬鈴薯、大豆，還有少量谷子，其余零星播種作物不納入統計(表3)。主要農作物生長期的NDVI特征統計如圖3所示。由圖3可知，玉米和谷子的NDVI較高，馬鈴薯次之，大豆最低。5月所有作物NDVI略有下降，可能是播種期除草所致；6—8月NDVI持續增長，于8—9月間達到峰值；9—10月收割期NDVI陡降；11月天氣轉冷，雜草枯死，NDVI繼續降低。撂荒地的NDVI在8月以前持續上升，于8月達到峰值，8月后逐漸下降。撂荒地與農作物在8月峰值的區別度較大，在5—8月、8—10月的變化梯度差異更為顯著。總體而言，研究區的主要農作物與撂荒地的NDVI峰值和變化梯度具有較高的區別度，以此作為閾值分割的基礎和聯合變化檢測的前提。

表3 主要農作物特征Tab.3 Characteristics of main crops

圖3 農作物生長期NDVI曲線Fig.3 NDVI curves of crops during growth period

1.3.3閾值分割算法實驗

圖4 NDVI閾值分割算法對比Fig.4 Comparison of NDVI threshold segmentation algorithms

閾值分割是變化檢測的核心，分割精度直接決定方法總體精度[16]。目前較為成熟的自動閾值分割算法包括基于直方圖形狀的Otsu算法[17]、基于力矩的Tsai算法[18]、基于信息熵的Kapur算法[19]及基于直方圖雙峰的Kittler算法[20]。本研究采用基于樣本統計的閾值分割方法。以年內檢測法為例，5種算法進行對比如圖4所示。由圖4可知，Kapur算法與Kittler算法振幅較大，極不穩定。Otsu算法與Tsai算法振幅較小，且具有較好的一致性。樣本統計分割閾值經2017年實地驗證精度為97.6%(2.2節)，故以2017年的閾值分割及地物信息為基準，控制2014—2016年的閾值分割?？傮w而言，研究區的閾值分布區間為0.3～0.4，最接近此區間為Tsai算法。

1.3.4聯合變化檢測

聯合變化檢測由年內檢測法和年際檢測法組成，其原理是：農作物由于耕作措施，在成熟期的NDVI最大值通常高于自然植被(雜草、灌木)。且農作物收割后，其NDVI在秋季會發生陡降，而自然植被不收割，故農作物在生長期的NDVI變化梯度顯著高于自然植被(圖3)。基于此，分別做3次變化檢測，按數據時相分為2次年內檢測和1次年際檢測，對檢測結果做交集、并集處理以提取撂荒地。

年內檢測法：以土地變更調查數據中的耕地地類為總耕地，通過遙感影像提取正在種植耕地，用總耕地減去正在種植耕地，間接提取未種植耕地，即撂荒地。首先，以耕地矢量圖層對遙感圖像進行掩膜，去除耕地范圍外的變化地物干擾，提取總耕地。然后，利用農作物NDVI的季節變化分別進行春-夏、夏-秋兩次變化檢測，對兩次檢測結果取交集，提取“正在種植耕地”。由圖3可知，春、夏、秋NDVI變化梯度最大的時相分別是5、8、10月，宜選用相應月份的遙感數據。變化檢測采用圖像差值法，平滑核為3，最小聚類值為20，閾值分割以農作物統計結果和算法實驗為參照。最后，以總耕地減去正在種植耕地得到未種植耕地，即撂荒地。

年際檢測法：撂荒地在撂荒前和撂荒后的地表覆被不同，NDVI變化明顯。經實地調研，研究區的撂荒地地表覆被主要為草地，新撂荒的地塊雜草較為稀疏，長期撂荒的地塊雜草較為茂密?；诖?，以耕地矢量圖層對遙感圖像掩膜，將ROI限定在耕地范圍內，并對前后兩年的夏季遙感影像ROI進行變化檢測，將NDVI明顯小于前一年的耕地斑塊提取出來，即當年新增加的撂荒地。

聯合變化檢測方法提取全類型撂荒地的計算式為

A=AC+AH+AT

(1)

其中

AC=A1∪A2-AT-AH

(2)

AH=A2-A2∩AT

(3)

AT=T+B

(4)

式中A——總撂荒AC——完全撂荒

AH——半撂荒AT——過渡撂荒

A1——年內檢測法提取撂荒斑塊集合

A2——年際檢測法提取撂荒斑塊集合

T——退耕還林B——設施占地

2 結果與分析

2.1 聯合變化檢測提取結果

利用多源數據、聯合變化檢測方法對研究區的撂荒地進行全類型提取。以2014年為基準年，全類型撂荒地提取結果見圖5(2014年沒有退耕還林工程，2017年尚未完成小班驗證)。由聯合變化檢測結果可知，種植耕地主要分布于和林格爾縣西北部土默川平原，及縣域內兩條河流沿岸區域。西北部平原因地形平整、易于灌溉和機械化種植，撂荒較少。同時，河岸附近耕地多為沖洪積地貌，地形相對平整，方便耕種且農作物產量較高，幾乎不存在撂荒現象。從河岸向兩側山區輻射，因機械化種植的難度加大，農作物產量降低，導致撂荒地逐漸增多，可見地形是影響和林格爾縣山區丘陵耕地撂荒的重要自然因素。

圖5 聯合變化檢測方法提取結果Fig.5 Extraction result maps of joint change detection method

2.2 精度評價

圖6 驗證點位Fig.6 Verifying point bitmap

結合野外實地驗證和目視解譯兩種方法評價精度。在提取的2017年撂荒地圖斑上，用ArcGIS 10.5軟件生成均勻分布的450個隨機點(圖6)，其中年際檢測法提取230個點，年內檢測法提取220個點，最小點距為150 m。對450個驗證點中的120個點進行野外實地考察驗證，其余330個驗證點由Sentinel-2A(分辨率10 m)數據進行目視解譯驗證。經驗證，聯合變化檢測方法于2017年提取撂荒地總體精度為97.6%，說明撂荒地提取結果真實可信，方法可行。此外，基于圖5b分類結果生成200個隨機點以評估分類精度。參照數據為同年GF-2數據(1 m全色、4 m多光譜)，評估方式為目視驗證，評估結果見表4。

基于表4計算2015年的總體分類精度為94.5%，Kappa系數為0.88。保障方法有效性及觀測精度需注意：①研究前期進行實地踏勘，建立農作物及撂荒地樣本知識庫。②多源數據的組合原則：無云數據優先、高時空數據優先、高質量數據優先。③以高質量數據(Sentinel-2A)為基準，盡量削弱多源數據的系統差異，平滑NDVI突變。④閾值分割以各地類樣本統計為基礎，綜合考慮時相跨度和年際降雨量變化。

表4 混淆矩陣Tab.4 Confusion matrix

2.3 區域撂荒綜合分析

完全撂荒主要分布于西南及東南部的黃土丘陵區以及東北部山區。值得注意的是，完全撂荒在一些區域(西南、東南)已呈集中連片之勢，并完全取代了種植耕地，說明此區域人口大量遷出、村莊衰落，農業完全廢棄。

半撂荒在全縣零星分布，分布于平原區的半撂荒多為休耕地，待地力恢復后繼續耕種；分布于山區丘陵的半撂荒則有很大概率在第2年繼續撂荒，最終轉為完全撂荒。

過渡撂荒在研究區有兩種形式：退耕還林和設施占地。研究區每年退耕還林666.7 hm2耕地，各個鄉鎮平均分配，山區鄉鎮較多。設施占地方面較大型的占地項目為“盛樂國際機場項目”，計劃征收2 266.7 hm2土地，主要涉及西北4個行政村，截至目前尚未征收完畢。對已有退耕還林及設施占地數據進行數字化、坐標轉換、配準等處理，并與其他數據耦合分析以提取過渡撂荒。對種植耕地和全類型撂荒地進行統計如圖7所示，其中，農作物播種面積數據來源于和林格爾縣2014—2017年的國民經濟和社會發展統計公報。

圖7 種植耕地及全類型撂荒地面積Fig.7 Area of cultivated land and all types of abandoned land

與統計公報中的全年農作物播種面積相比，所述聯合變化檢測法提取的種植耕地面積與之相仿，4年組間平均偏差為4.9%，表明所述方法的提取結果較為可信。研究區撂荒類型以完全撂荒為主，半撂荒和過渡撂荒較少；撂荒地多分布于山區丘陵，平原分布較少。該縣2014年撂荒率為37.8%，2015年為32.4%，2016年為36.7%，2017年為39%，4年平均撂荒率為36.5%，撂荒情況較為嚴重。總體撂荒規模相對穩定，組內偏差最大的2015年為11.1%，其余3年偏差較小。結合農戶問卷和野外踏勘發現，造成研究區高撂荒率的原因是多方面的。內源性因素為：①和林格爾縣山區、黃土丘陵區地貌接近80%，耕作難度大、人力成本高。②研究區的耕地質量較差，絕大部分為11～15等耕地，土質貧瘠，產出較低。③研究區為典型的農牧交錯區，農牧文化并存，一些農民更愿意棄耕還牧。④和林格爾縣距自治區首府呼和浩特僅60 km，根據人口遷移的“引力理論”[21]，較高的人口規模及GDP、較短的距離，會對樣本縣人口產生較強的吸附力。外源性因素為：①隨著城市化進程的加快，勞動力老齡化逐年加劇。②較高的種植成本及較低的機械化程度導致山區耕地流轉困難。平原區的半撂荒多為主動撂荒，具備復耕潛力；山區的完全撂荒則多為被動撂荒。高撂荒率會進一步加速村莊邊際化，直至村莊完全廢棄。

2.4 空間特征提取及統計分析

以聯合變化檢測法提取的撂荒地斑塊為基礎，結合30 m DEM數據，以及圖6的隨機驗證點和二調數據的旱地、水澆地、建筑用地、村莊邊界等信息，可以提取4種空間特征指標，分別是：撂荒地塊距村莊中心距離、撂荒地塊距村莊中心高差、撂荒地塊的灌溉類型(旱地/水澆地)、撂荒地塊的鄰接關系。考慮到坡度可由距離和高差求得，導致多重共線性，故不予提取?？臻g特征統計信息見圖8。

由圖8a可知，大多數撂荒地塊分布在所屬村莊的2 km以內。究其原因，平原區人口稠密，村莊之間距離很近；而山區的部分村莊完全廢棄，全部撂荒。計算距離平均值發現，研究區的撂荒地距所屬村莊中心的平均距離為1.17 km，可見距離對撂荒的影響并不顯著。

以村莊中心高程為基準高程，基于30 m DEM數據對撂荒地距所屬村莊的高差進行統計，得出研究區撂荒地平均高于村莊基準面40.8 m，且近40%的撂荒地塊高于村莊基準面60 m以上(圖8b)，撂荒地受高差因素影響顯著。結合實地調研分析，高差對撂荒的驅動作用力限制了農用機械的使用，增加了勞動力成本。

以二調數據為本底，對研究區耕地、撂荒地的旱地/水澆地屬性類別進行統計(圖8c)，得出全縣耕地的旱地比例為72.3%，全縣撂荒地的旱地比例為89.6%。撂荒地旱地比例遠高于耕地旱地比例，說明灌溉條件是決定耕地撂荒與否的重要因素之一。

圖8 空間特征統計四維雷達圖Fig.8 Four-dimensional radar maps of statistical spatial characteristics

由圖8d可得，樣本縣81.9%的半撂荒地塊緊鄰完全撂荒地塊，83.1%的完全撂荒地塊彼此相鄰。一方面農戶決策易受相鄰其他農戶決策的影響；另一方面，局部區域的耕地內源性因素相似，撂荒地會向相似條件的地塊蔓延。

3 結論

(1)提出了耦合年內檢測法和年際檢測法的聯合變化檢測方法，為基于遙感的區域撂荒地識別提供了快速有效的技術方法與體系。該方法對撂荒地的提取精度為97.6%，全類型分類精度為94.5%，Kappa系數為0.88。

(2)研究區撂荒地距村莊中心的平均距離為1.17 km，平均高差為40.8 m，旱地比例為89.6%(耕地為72.3%)，鄰接比例超過81%，結果表明研究區的撂荒驅動因素以高差因素、鄰域因素和灌溉因素最為顯著，距離因素不顯著。針對多源數據的特征提取及顯著性分析，有助于判別區域撂荒主導因素，為撂荒驅動力研究、定向提升撂荒地管理方法提供依據。

(3)研究區2014—2017年撂荒率均超過32%，南部黃土丘陵區部分村莊完全廢棄，全部撂荒。