基于Google Earth Engine與機(jī)器學(xué)習(xí)的黃土梯田動態(tài)監(jiān)測

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

（寧夏大學(xué) 農(nóng)學(xué)院，寧夏 銀川 750021）

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

1 研究區(qū)概況

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

2 數(shù)據(jù)源與研究方法

黃土梯田動態(tài)監(jiān)測的流程可分為4個主要功能模塊：遙感數(shù)據(jù)加載、數(shù)據(jù)預(yù)處理、分類算法優(yōu)選、序列優(yōu)化。各模塊從上到下，層層遞進(jìn)，最終實(shí)現(xiàn)黃土梯田動態(tài)監(jiān)測(圖1)。

圖1 黃土梯田動態(tài)監(jiān)測流程圖Figure 1 Flowchart of dynamic monitoring of loess terraces

2.1 數(shù)據(jù)源

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

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

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

2.2 研究方法

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

表1 光譜指數(shù)計(jì)算方法Table 1 Calculation methods of spectral index

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

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

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

2.2.5 精度驗(yàn)證方法 采用混淆矩陣的方法，以總體精度、Kappa系數(shù)、生產(chǎn)者精度和用戶精度等指標(biāo)作為識別精度評價依據(jù)。具體計(jì)算方法參考文獻(xiàn)[18]。

2.3 植被覆蓋度計(jì)算

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

3 結(jié)果與分析

3.1 不同機(jī)器學(xué)習(xí)算法的精度

表2為隨機(jī)抽取的1 051個樣點(diǎn)的驗(yàn)證結(jié)果。4種精度指標(biāo)均為隨機(jī)森林算法最高，決策樹算法次之，支持向量機(jī)算法最小。隨機(jī)森林算法基于樣點(diǎn)檢驗(yàn)的精度分別為：梯田的生產(chǎn)者精度94.46%、梯田的用戶精度89.03%、總體精度94.10%、Kappa系數(shù)為0.87，都遠(yuǎn)大于另外2種算法。因此，后文采用隨機(jī)森林機(jī)器學(xué)習(xí)算法進(jìn)行梯田遙感識別。

表2 不同機(jī)器學(xué)習(xí)算法識別結(jié)果的樣點(diǎn)驗(yàn)證精度Table 2 Sample points verification accuracy of the results of different machine learning algorithms

3.2 識別結(jié)果與人工勾繪斑塊的驗(yàn)證

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

3.3 LandTrendr擬合效果示例

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

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

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

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

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

3.5 梯田使用時間分析

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

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

4 討論

4.1 機(jī)器學(xué)習(xí)識別精度的影響與優(yōu)化

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

4.2 梯田面積變化分析

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

5 結(jié)論

基于GEE云平臺，使用隨機(jī)森林機(jī)器學(xué)習(xí)算法與LandTrendr算法，可以高效、準(zhǔn)確地實(shí)現(xiàn)長時間序列、大尺度的黃土梯田動態(tài)監(jiān)測。相比1988年，研究區(qū)2019年梯田面積減少45.90%，植被覆蓋度增長52.44%，說明近30 a梯田農(nóng)業(yè)比例逐漸降低，生態(tài)環(huán)境持續(xù)向好發(fā)展。