基于Google Earth Engine與機器學習的黃土梯田動態監測

李萬源，田 佳，馬 琴，金學娟，楊澤康，楊鵬輝

（寧夏大學 農學院，寧夏 銀川 750021）

寧夏回族自治區固原市位于黃土高塬溝壑區[1]。長期的過度放牧、不合理耕作，導致該地區植被稀疏、水土流失加劇[2]，嚴重影響了當地社會經濟發展和生態安全。梯田有效緩解了農業生產帶來的水土流失問題[3]，從20世紀80年代開始，固原市實施了大面積的坡改梯工程[4]。加之2000年開始實施的國家退耕還林還草工程[5]，該地區的水土流失問題有所緩解，生態環境持續向好[6]。隨著遙感技術的快速發展，如何從遙感影像中高效、準確、大尺度地獲取梯田時空分布信息，對于指導農業生產、水土保持監測和防治水土流失具有重要的意義。傳統的梯田遙感識別主要采用目視解譯[7]，該方法精度較高，但存在耗時耗力、成本高、方法復用性差等問題，目前更多用來采集機器學習(machine learning)的樣本[8]。近年來，大部分學者采用面向對象或基于像元的監督識別技術，利用決策樹(CART)、隨機森林(RF)、支持向量機(SVM)、深度學習(DL)等[9?11]機器學習算法，先學習采集的樣本，然后利用學習好的模型識別新的樣本。面向對象技術較基于像元識別技術，不僅依靠地物的光譜特征，還利用像元和像元之間的關系提高識別精度，識別過程更加復雜，影像分辨率要求更高[7]。但是，無論采用哪種方法進行梯田遙感識別，基本上都是基于單機處理，普遍存在遙感數據獲取困難、預處理復雜、性能限制等問題[9]，難以開展大尺度的遙感識別研究。為了解決這些問題，Google公司借助其強大的計算資源與海量數據存儲，推出了遙感云平臺Google Earth Engine(GEE)[12]。借助該平臺，研究人員可以極大擴展自身原有研究的覆蓋范圍，提供國家乃至全球尺度的研究成果[13]。目前，GEE在大尺度森林變化監測、土地利用類型分類、人類居住地動態監測等[14?16]方面應用廣泛，但大尺度梯田遙感識別未見相關報道。為此，本研究在GEE平臺支持下，利用Landsat時間序列數據和SRTM數字高程模型(digital elevation model，DEM)，建立每年時間序列影像的百分位數特征。對比3種機器學習算法的分類精度大小，選擇分類精度最高的識別結果，應用LandTrendr時序算法逐像元擬合修正時間序列，實現固原市1988?2019年度梯田動態監測的目的。研究結果可為黃土丘陵地區梯田的高效、準確識別和水土保持監測、評價提供參考。

1 研究區概況

固原市 (35°14′～36°31′N，105°19′～106°57′E)位于寧夏回族自治區南部的六盤山地區，轄原州區、西吉縣、隆德縣、涇源縣、彭陽縣，國土面積1.05 萬km2。屬大陸暖溫帶半干旱氣候，年均氣溫6.3 ℃，年均降水量493.5 mm，降水量由東南向西北遞減，年均蒸發量1 472.9 mm，年均無霜期152.0 d。域內地形南高北低，溝壑縱橫，黃土丘陵面積達67.9%。地帶性土壤以黑壚土為主，但嚴重的土壤侵蝕導致土壤母質層出露，黃綿土廣布。植被總體上由東南的半濕潤森林草原區向西北的干旱半干旱草原區過渡[4]。

2 數據源與研究方法

黃土梯田動態監測的流程可分為4個主要功能模塊：遙感數據加載、數據預處理、分類算法優選、序列優化。各模塊從上到下，層層遞進，最終實現黃土梯田動態監測(圖1)。

圖1 黃土梯田動態監測流程圖Figure 1 Flowchart of dynamic monitoring of loess terraces

2.1 數據源

2.1.1 Landsat影像 使用 T1級別 (質量最高)的 Landsat地表反射率數據 (surface reflectance, SR)。該數據產品已經過幾何校正、輻射校正和大氣校正，空間分辨率30 m，時間分辨率16 d。由于Landsat 5/7/8衛星的服務年限不同，1988?2011年使用Landsat 5影像，2012年使用 Landsat 7影像，2013?2019年使用Landsat 8影像，共使用1 690景影像。

2.1.2 高程數據 采用 30 m 空間分辨率的數字高程模型，具體編號為SRTMGL1_003。

2.1.3 樣本數據 地類僅分為梯田和其他 2類。通過Google Earth Pro提供的高清歷史影像，利用目視解譯法采集樣本數據。樣本數據包括樣點數據和斑塊數據。樣點數據按時間分為2010?2014年地類屬性相同和2000年的樣點，以滿足Landsat 5/7/8不同衛星分別進行機器學習樣本訓練的需求。樣點采集遵循以下原則：①在研究區生成5 km方形格網，以使樣點分布均勻；②保持樣點100 m以內屬性相同。樣點數據共2 673個，梯田樣點1 040個，其他樣點1 633個。斑塊數據為6個隨機分布的5 km×5 km 正方形區域，參考 Google Earth Pro中 2019年厘米級高清遙感影像人工勾繪以及實地驗證。

2.2 研究方法

2.2.1 合成影像 選擇 Landsat對應衛星影像的紅波段 (Br)、綠波段 (Bg)、藍波段(Bb)、近紅外 (Bnir)、短波紅外1(Bswir1)、短波紅外2(Bswir2)6個光譜波段；再經裁邊(壞像元)、光譜指數計算(計算方法如表1)、去云后，針對黃土梯田全年季相變化特點[17]，統計每年度內時間序列影像百分位數特征融合影像[18]，即逐像元對某一波段1 a內所有觀測值取其10%、25%、50%、75%、90%百分位數，獲得該像元位置該波段對應的5個指標波段；再與6個地形特征波段組合，即由數字高程計算得到的海拔、坡度、坡向，以及 3個 3×3、7×7、11×11像元窗口內地形起伏度波段。共計61個特征波段。

表1 光譜指數計算方法Table 1 Calculation methods of spectral index

2.2.2 機器學習 3種機器學習算法為隨機森林、決策樹、支持向量機，GEE均有內建，可直接調用。另外，針對不同衛星分別進行機器學習，把樣點數據分年度映射到對應合成影像并匯總(如Landsat 5包括2000、2010和2011年的樣本)，再按9∶1劃分樣本，90%的樣本用于分類器訓練，10%的樣本用于精度驗證。

2.2.3 LandTrendr算法 LandTrendr算法將以年時間序列的值進行分割、逐段擬合、平滑[19]，獲取單個像元在整個研究時間段內的整體變化特征。具體介紹參考文獻[19]。

2.2.4 識別結果優化 應用前文分類精度最高的機器學習算法，對研究區1988?2019年逐年進行梯田遙感識別。為減少極端天氣和人類活動導致識別錯誤，利用地類在時間序列上連續、穩定的特征，使用LandTrendr算法[19]對識別結果的時間序列(概率為0～1的浮點)擬合平滑處理。參考中國水土保持措施分類[20]，提取坡度＞2°和坡度＜25°區域的梯田，以減少溝壑地及塬地的誤分。

2.2.5 精度驗證方法 采用混淆矩陣的方法，以總體精度、Kappa系數、生產者精度和用戶精度等指標作為識別精度評價依據。具體計算方法參考文獻[18]。

2.3 植被覆蓋度計算

植被覆蓋度(fractional vegetation cover, FVC)采用歸一化植被指數和像元二分模型計算。具體計算方法參考文獻[21]。

3 結果與分析

3.1 不同機器學習算法的精度

表2為隨機抽取的1 051個樣點的驗證結果。4種精度指標均為隨機森林算法最高，決策樹算法次之，支持向量機算法最小。隨機森林算法基于樣點檢驗的精度分別為：梯田的生產者精度94.46%、梯田的用戶精度89.03%、總體精度94.10%、Kappa系數為0.87，都遠大于另外2種算法。因此，后文采用隨機森林機器學習算法進行梯田遙感識別。

表2 不同機器學習算法識別結果的樣點驗證精度Table 2 Sample points verification accuracy of the results of different machine learning algorithms

3.2 識別結果與人工勾繪斑塊的驗證

表3顯示：去除交界100 m緩沖區后的驗證精度高于未去除時(0 m)的驗證精度。另外，經LandTrendr處理后梯田的生產者精度、梯田的用戶精度、總體精度和Kappa系數分別為：81.75%、85.97%、93.33%、0.80，均大于LandTrendr處理前的驗證精度。

3.3 LandTrendr擬合效果示例

選擇3個不同位置來展示LandTrendr算法擬合效果(圖2)，位置A原始識別結果在1994、2002、2004年被錯誤識別為其他類型，位置B原始識別結果在1997年被錯誤識別為其他類型，在2015年被錯誤識別為梯田類型。經LandTrendr算法處理后，這些錯誤類型均被校正。位置C原始識別結果與經LandTrendr算法處理后的結果均為其他類型，識別類型沒有變化。

圖2 3個不同位置的原始識別結果及使用LandTrendr算法處理后的概率Figure 2 Classification probability of the original results and the results of using LandTrendr algorithm at 3 different positions

3.4 近30年梯田面積和植被覆蓋度

經LandTrendr算法處理后的研究區梯田面積(圖 3)變化趨勢更穩定，從 1988年 5 816.59 km2減少到 2019年 3 146.72 km2，年均減少 90.85 km2·a?1。1988?2019年，研究區植被覆蓋度則呈現不斷增加的趨勢，與梯田面積變化趨勢相反。另外，處理前、處理后的梯田面積與植被覆蓋度極顯著(P＜0.001)相關，其相關系數分別為?0.50和?0.75。

圖3 1988?2019年研究區梯田面積與植被覆蓋度變化Figure 3 Variations of annual terraces area and annual mean fractional vegetation cover in the research area from 1988?2019

3.5 梯田使用時間分析

圖4顯示了研究區1988?2019年梯田使用時間長短的分布。從整體上來看，梯田主要分布在六盤山山脈兩側，且西部的梯田使用時間較東部更長。從局部來看，南部的涇源縣區域，梯田零星分布，使用時間相對較短；西部西吉縣的溝谷條帶、中部的六盤山山脈、北部原州區清水河的河谷沖積平原(紅色部分)能明顯區分出來。

圖4 1988?2019年研究區梯田使用時間分布示意圖Figure 4 Distribution of time to use terraces in the research area from 1988?2019

4 討論

4.1 機器學習識別精度的影響與優化

已有的梯田遙感監測研究[3?4]受限于單機處理性能和準確的歷史樣本采集，其研究內容往往時間短、區域小，限制了長時間序列、大尺度遙感監測的應用與發展。本研究使用模型遷移法，針對每一個傳感器獨立訓練機器學習分類模型，減少了樣本采集的難度，得以實現黃土梯田動態監測。然而，機器學習的識別精度主要受樣本量、特征、機器學習算法的影響[7]。本研究利用多年采樣法增加樣本量，選取最優機器學習算法，得到較高的識別精度。另外，關于特征選取，我們前期使用了最大值、最小值、眾數、中位數、平均數等多種特征融合方法，但識別精度均低于本研究的百分位數特征融合。而對于深度學習，我們在本地電腦使用相同樣本集，多次構建深度學習模型并訓練，然而識別精度也低于本研究的隨機森林。最后引入LandTrendr算法逐像元擬合時間序列軌跡，有效校正了如圖2中的異常值，提高了識別精度。而且，在斑塊驗證數據與樣點采集時同樣保留100 m空間誤差時，消除2種利用類型相鄰區域地理配準誤差帶來的系統錯誤后，基于2019年斑塊檢驗總體精度93.33%，與樣點驗證總體精度94.10%相當，說明訓練好的模型隨時間遷移應用，識別性能不會降低。

4.2 梯田面積變化分析

整體來看，研究區1988?2019年梯田面積呈減少趨勢，植被覆蓋度則呈現逐步增長趨勢，梯田面積與植被覆蓋度極顯著相關(P＜0.01)，說明梯田面積減少有助于生態環境向好發展。局部來看，研究區在1988?1996年梯田面積年均減少69.02 km2·a?1，與該時期寧南山區逐步退耕還林還草時間一致；1997?2000年梯田面積年均增長 91.60 km2·a?1，與該時期耕地面積增長趨勢相同[22]；在 2001?2005年梯田面積下降較快，梯田面積年均減少 250.51 km2·a?1，遠高于 1988?2019年年均減少速率 90.85 km2·a?1，且2001?2005年植被覆蓋度年均增長速率是1988?2019年植被覆蓋度年均增長速度的4倍，與寧夏“退耕還林工程”生態政策大力實施的時間節點相符；從2007年開始，為鞏固退耕還林工作，持續推進生態文明建設，研究區梯田面積下降減緩[6]。另外，研究區西部的梯田使用時間較東部更長，這可能與東部年降水量達650 mm，而西部年降雨量不到450 mm[4]，在東部進行梯田退耕后有利于提高植被成活率有關。

5 結論

基于GEE云平臺，使用隨機森林機器學習算法與LandTrendr算法，可以高效、準確地實現長時間序列、大尺度的黃土梯田動態監測。相比1988年，研究區2019年梯田面積減少45.90%，植被覆蓋度增長52.44%，說明近30 a梯田農業比例逐漸降低，生態環境持續向好發展。