基于集成學習和多時相遙感影像的枸杞種植區分類

史飛飛， 高小紅， 肖建設， 李宏達, 李潤祥， 張昊

(1.青海師范大學地理科學學院，西寧 810008； 2.青海省氣象科學研究所，西寧 810008； 3.青海省自然地理與環境過程重點實驗室，西寧 810008； 4.青藏高原地表過程與生態保育教育部重點實驗室，西寧 810008； 5.青海省防災減災重點實驗室，西寧 810008;6.高原科學與可探突發展研究院，西寧 810008)

0 引言

青海柴達木盆地作為我國枸杞的新興種植區，近年來因枸杞種植規模的無序擴大以及市場開拓不足致使供需失衡嚴重，已顯現出枸杞價格持續低迷和田間消極管理引發減產等問題，而實現枸杞種植面積監測與有序調控是保障枸杞市場平穩、持續發展的關鍵手段[1-2]。與傳統地面調查方法相比，遙感技術具備宏觀、經濟和高效等優勢，已廣泛應用于作物的分類、種植結構提取和種植面積監測等方面[3]。

作物遙感識別與分類是獲取農業土地利用信息的基礎，近年來采用遙感觀測資料進行作物分類的研究較多，其大致可分為2大類[4-5]： 一種為選用作物生長關鍵期內的單時相影像進行作物分類； 另一種為選用作物生長期內的多時相影像，經過構建時序的植被指數進行作物分類。受我國作物種植結構復雜性和作物光譜相似性影響，在采用單時相影像開展作物分類時，存在“最佳識別期”不易選取以及易出現“同譜異物，同物異譜”的問題，限制了該方法的分類精度[6]。而在同一區域內因不同作物的生長周期普遍存在差異，采用多時相影像能夠獲取光譜與物候特征信息，具備獲取更佳作物分類精度的潛力[7]。在該方法應用初期以采用具有高時間分辨率的MODIS數據為主，并對全球廣泛種植的大宗作物開展了遙感監測[8-10]，但因MODIS數據空間分辨率較低，且在種植結構復雜的中小規模種植區易出現大量混合像元，致使分類精度難以提升[11-12]。伴隨遙感技術的發展，近年來涌現出利用多時相的Landsat7/8[13-15]、GF-1[16-17]和Sentinel-2[18-21]等中高空間分辨率的影像對香梨、葵花、苜蓿等特色經濟作物進行提取的案例。通過對比上述研究發現，時序高分影像數據能獲得豐富的作物冠層光譜、紋理和物候等特征信息，提高了作物識別的準確性，但也面臨如單一數據源易存在關鍵物候期影像缺失，采用決策樹和神經網絡等分類器存在分類特征利用不充分等問題[22-23]。針對上述問題，多源數據應用能較好地解決關鍵時相數據的缺失，并且已有研究對Landsat8 OLI，Landsat7 ETM+和Sentinel-2A MSI等傳感器，進行了波段特性對比以及數據間差異消除方法的構建[24-25]，能為多源數據綜合提供方法參考，但同時嘗試引入具有更強學習與泛化能力的分類器可能是提升作物分類精度的關鍵[23,26]。任何單一的強分類器都有其自身的優勢與不足，集成學習(ensemble learning)方法能綜合運用各分類器的優勢，以提高分類精度，這也是當前研究的熱點方向[27]。集成學習依照個體學習器的形成方式可劃分為串行或并行的序列化，前者典型代表為提升法(boosting)，而后者的代表是自舉匯聚法(bootstrap aggregating,bagging)和隨機森林(random forest,RF)，其中RF對多維特征和數據噪聲具有更強的魯棒性，已廣泛應用于作物分類并取得較好效果[28-29]，而采用boosting算法進行作物分類鮮有嘗試。

本文以柴達木盆地典型枸杞種植基地諾木洪農場為研究區，選取Landsat8 OLI和GF-1 WFV影像構建作物生長期內時序植被指數數據，并采用輕量級梯度提升機(light gradient boosting machine,LightGBM)、梯度提升決策樹(gradient boosting decision tree, GBDT)、極端梯度提升(extreme gradient boosting,XGBoost)、RF、支持向量機(support vector machine,SVM)和多層感知分類器(multilayer perceptron classifier,MLPC)共6種分類器，旨在探討集成學習方法在開展枸杞種植區分類時的適應性。

1 研究區概況

諾木洪農場位于柴達木盆地都蘭縣境內，生態環境十分純凈，是全國規模較大的枸杞種植基地，被譽為“中國枸杞之鄉”。農場地理坐標E96°15′～96°35′，N36°20′～36°30′，海拔為2 745～2 790 m，總面積約200 km2，諾木洪農場及周邊枸杞種植區分布狀況見圖1。農場屬高原大陸性氣候，輻射強，晝夜溫差大，氣候干燥，年平均氣溫約5 ℃，降雨量約51 mm，而蒸發量在2 200 mm以上。農場地處沖洪積扇前緣的細土平原帶，土壤肥沃并被大量開墾進行枸杞種植，南部為山前荒漠戈壁，而北部為沼澤地，諾木洪河為柴達木流域內第六大河流，由南至北穿過農場為其生產生活提供充足水源[30]。

圖1 枸杞種植區分布

2 數據源與研究方法

2.1 數據采集與預處理

2.1.1 影像數據與預處理

為構建具備高時間分辨率的時序數據集，選取研究區內作物生長階段4—11月份內的6景Landsat8 OLI影像和3景GF-1 WFV影像，其中Landsat8 OLI影像在谷歌地球引擎(Google Earth Engine,GEE)平臺下載(https: //code.earthengine.google.com)，空間分辨率30 m，重訪周期16 d，GF-1 WFV數據在中國資源衛星應用中心下載(http: //www.cresda.com/)，空間分辨率為16 m，波段數為4，重訪周期為2 d。數據詳情如表1所示。

表1 影像數據列表

①OLI為陸地成像儀(多光譜)，WFV為多光譜寬幅相機。

GF-1 WFV數據利用ENVI5.5軟件進行預處理，以鄰近時期的Landsat8 OLI數據和30 m DEM數據為參考，選用自動采集控制點的正射校正工具(RPC orthorectification using reference image, RPC)進行校正，并將GF-1 WFV進行輻射定標和FLAASH(fast line-of-sight atmospheric analysis of spectral hypercubes)大氣校正。Landsat8 OLI數據在GEE平臺中可直接獲取地表反射率數據，但仍需在GEE平臺中通過編寫腳本使用多通道掩模算法(C function of mask, CFMask)對影像中的云、陰影和雪進行處理[31]。最后將預處理后的GF-1 WFV和Landsat8 OLI影像空間分辨率重采樣為30 m，并將投影統一轉換為Albers等面積投影。

2.1.2 實地野外調查數據

為獲取準確且均勻分布的地面樣點調查數據，用于分類器的訓練和精度驗證，按1 km×1 km格網將研究區劃分為264個單元，于2017年9月在諾木洪農場開展實地調查，通過隨機選取格網以及參考調查地點的可通達性，共對101個格網開展地面調查，期間利用GPS對主要地物類型進行定位并拍照，使用ASD FieldSpec 4型地物光譜儀獲取枸杞作物的冠層光譜。因研究區內分布大范圍的沼澤地和荒漠戈壁致使部分格網單元無法開展實地調查，因而利用GF-1影像通過目視判識獲取該格網內的地物類型信息。最終通過調查將研究區待分類對象劃分為幼齡期枸杞地、壯齡期枸杞地、行道樹、城鎮、沼澤地、戈壁灘和公路共7種類型，并獲得各類地物樣點共1 378個，樣點分布見圖2。

圖2 研究區樣點分布

2.1.3 物候期觀測數據

2016—2018年在諾木洪農場選取試驗田開展作物物候觀測，觀測地枸杞樹齡為5 a左右，圖3展示了枸杞樹的生長周期，一般枸杞樹從5月上旬開始進入萌芽期，5月下旬進入春梢生長期，6月下旬進入開花盛期，7月上旬—9月下旬進入果熟期，期間枸杞樹不斷開花結果，當地果農一般開展2～3茬摘果，10月初進入落葉始期，至11月初果樹開始進入休眠期。

圖3 枸杞作物物候期

2.2 植被指數選取

植被指數能增強光學遙感影像中的植被信息，并被廣泛用于描述生物物理特性的變化[32]。其中歸一化植被指數(normalized difference vegetation index, NDVI)是監測植被覆蓋度和生長狀況的最佳指數，但受計算方法的影響容易在高植被覆蓋區產生飽和，而增強植被指數(enhanced vegetation index，EVI)加入了土壤背景調節參數和大氣修正參數，能彌補NDVI指數的不足[33]。使用植被指數進行作物物候期監測時，如對地面植被類型及覆蓋度狀況不了解，一般應綜合選取NDVI和EVI指數互為補充[33]，2種植被指數的計算公式分別為：

(1)

(2)

式中:ρNIR，ρRed和ρBlue分別為近紅外、紅光波段和藍光波段的反射率值；L，C1和C2分別為土壤背景調節參數、大氣校正紅光參數和大氣校正藍光參數，其中L=1，C1=6，C2=7.5。

2.3 時序植被指數重構

受云霧和水汽因素影響，利用Landsat8 OLI和GF-1 WFV影像構建的NDVI/EVI時間序列數據會夾雜大的波動和噪聲，需進行平滑和去噪[34]，本文使用SG(Savitzky-Golay)濾波[18]，其計算公式為：

(3)

2.4 分類器構建與精度評價

在Visual Studio Code軟件中使用sklearn (Scikit-learn)，LightGBM和XGBoost等第三方庫，實現了LightGBM，GBDT，XGBoost，RF，SVM和MLPC分類器搭建。當前集成學習是機器學習領域的研究熱點，它具有準確性高、速度快、參數少等優點，已在葉綠素反演和森林地上生物量估算等方面開展了應用[35-36]，其中LightGBM，GBDT和XGBoost均是在Boosting集成學習方式基礎上的改進算法，但在基分類器、列采樣方式、決策樹生長策略、殘差擬合方法和并行策略等方面存在差異[37]； RF是以決策樹作為基分類器，并采用bagging方式進行集成的，能夠處理高維數據，并且無需大訓練樣本、反復的模型訓練和參數調整便能獲取較好的分類精度[37]； SVM和MLPC是典型的非參數化機器學習方法，相比傳統的參數化方法具有較強的學習及泛化能力，適用于土地利用、作物分類和參量反演等場景[38]。

使用sklearn庫中的樣本分割方法，按照2∶8將樣本集隨機劃分為測試集和訓練集。同時分類器訓練時超參數的設定對于模型分類效果影響十分關鍵，一般在參數集合較大時采用隨機格網搜索法能夠高效選定候選超參數，文中設定10折交叉驗證和10次重復來評估超參數值的所有可能組合。通過實驗LightGBM的boosting類型設定為GBDT，并對學習率、最大葉子節點數、樹深度和子模型樹設定值范圍并進行尋優； GBDT和RF均對學習率、樹節點分裂最小樣本數、樹終節點最少樣本數、樹深度和子模型樹參數進行尋優； XGBoost中設定類別數為7，對樹深度、學習率和損失函數進行尋優； SVM的核函數設定為徑向基函數，并對懲罰系數進行尋優； MLPC對隱含層樹及節點數、激活函數、權重優化器和學習率進行尋優。

精度評價主要從像素尺度出發通過構建混淆矩陣，計算總體分類精度(overall accuracy,OA)、制圖精度(producer accuracy,PA)、用戶精度(user accuracy,UA)和Kappa系數進行分類精度評價[39]。

2.5 分類特征重要性評估

利用sklearn中的元轉換器并與LightGBM，GBDT，XGBoost和RF分類器的特征重要性評分方法配合，用于評估不同時相NDVI/EVI特征在作物分類中的重要性。而SVM和MLPC分類器缺乏特征重要性評分能力，利用單變量特征選擇通過對單變量的統計測試來抉擇較好的分類特征。

3 結果與分析

3.1 時序植被指數曲線特征

采用SG方法逐像元重建了時序NDVI/EVI圖像，并利用地面調查樣點提取了研究區內7種地物類型的平均時序NDVI/EVI曲線(圖4)，經過分析重建后的NDVI/EVI曲線更加平滑，消除了“鋸齒”和“毛刺”等噪聲，表明SG方法具有較好的濾波和特征保持能力； 壯齡/幼齡期枸杞在不同生長階段，其NDVI/EVI曲線表現的物候特征變化趨勢較為接近，但增減幅度存在差異，如壯齡期與幼齡期枸杞在休眠期NDVI/EVI值均為最低，5月上旬—中旬枸杞由芽開放期轉至展葉盛期，該階段NDVI/EVI曲線上升趨勢遲緩，5月下旬—6月下旬，枸杞由春梢開始生長至春梢開花盛期，該階段NDVI/EVI曲線均呈快速增長，但壯齡期枸杞增幅較大，7月中旬—8月中旬為夏果形成與成熟期，該階段枸杞NDVI/EVI曲線均達到峰值并保持穩定，其中壯齡期枸杞峰值高于幼齡期。9月中旬—下旬枸杞由秋果成熟盛期逐漸過渡至落葉始期，枸杞NDVI/EVI曲線逐漸開始降低。而到11月上旬，枸杞進入落葉盛期，枸杞NDVI/EVI曲線持續快速下降，伴隨氣溫驟降以及降雪過程的來臨，枸杞進入休眠期。研究區內行道樹和沼澤地植被NDVI/EVI曲線也均表現出明顯的物候特征，其中在田間道路兩旁種植的行道樹主要以挺直高大的新疆楊為主，生長期內NDVI/EVI值均高于枸杞等其他植被類型，而沼澤地植被的NDVI/EVI值也略高于幼齡期枸杞； 研究區內公路、城鎮和戈壁灘非植被類型的NDVI/EVI值均低于0.2，并且無明顯峰谷特征。

(a) NDVI曲線 (b) EVI曲線

圖4 NDVI和EVI曲線

3.2 時序植被指數曲線可分性

對比分析時序NDVI/EVI曲線對于枸杞等不同植被類型間的區分程度發現，因NDVI和EVI指數的計算公式導致同一類型的時序NDVI曲線值較EVI值整體偏高[19]。7種類型的時序NDVI/EVI曲線特征均存在差異，但是壯齡期枸杞與行道樹在時序EVI曲線中的差異更加明顯，并主要表現在枸杞生長期后期(8月下旬—11月上旬)。為定量分析NDVI/EVI指數在枸杞作物不同生長期內表現的差異性，選取壯齡期枸杞的NDVI/EVI時序曲線，并計算枸杞下一生長期與當前時期NDVI/EVI指數的增加或減少幅度，結果如圖5所示，分析發現在枸杞芽開放期至開花盛期(5月11日—6月28日，生長初期)，EVI與NDVI指數增幅呈快速上升趨勢，而后在枸杞夏果形成至成熟盛期(7月14日—8月2日，生長中期)，EVI與NDVI指數增幅降低其數值達到生長期內最大，在枸杞秋果形成至落葉盛期(8月14日—11月3日，生長中后期)，EVI和NDVI指數降幅逐漸呈上升趨勢。總體而言，EVI指數曲線相較于NDVI指數曲線在枸杞物候期的前期和中期增幅均較大，兩者差異并不明顯，而在中后期EVI指數曲線的降低幅度明顯高于NDVI指數曲線，兩者差異逐漸顯現。

圖5 NDVI及EVI指數增加與減少幅度

3.3 分類結果與精度評價

3.3.1 分類結果

利用NDVI/EVI時序數據采用6種分類器進行枸杞種植區分類，結果如圖6及圖7所示，在12種分類結果中各地物類型的細節展現總體較好，其中幼齡期枸杞種植區連片分布在農場南部，多為新開墾的枸杞幼苗種植區，壯齡期枸杞分布在農場大部地區，并被行道樹分割呈現規則的棋盤狀分布，上述3種類型邊界提取完整且準確。沼澤地主要分布在農場的北部，南部為廣闊的戈壁灘，而城鎮主要散落分布在枸杞種植區內，上述4種類型在部分分類結果中出現了相互混淆。對農場南部國道G109的提取存在間斷并易出現碎斑，其原因為道路寬度較窄采用30 m空間分辨率影像對其識別能力存在不足。

(a) LightGBM分類結果(b) GBDT分類結果

(c) XGBoost分類結果(d) RF分類結果

(e) SVM分類結果(f) MLPC分類結果

圖6 基于NDVI時序數據的分類結果

(a) LightGBM分類結果(b) GBDT分類結果

(c) XGBoost分類結果(d) RF分類結果

(e) SVM分類結果(f) MLPC分類結果

圖7 基于EVI時序數據的分類結果

3.3.2 精度評價

對12種分類結果進行精度驗證與對比，結果如圖8及圖9所示，分析表明以LightGBM+EVI的OA最高達到91.67%，Kappa系數為0.90，其對幼齡期枸杞和壯齡期枸杞的PA分別為91.0%和98.0%，而MLPC+EVI的OA最低為83.33%，Kappa系數為0.81，但對幼齡期枸杞和壯齡期枸杞的PA也分別達到82.0%和89.0%，表明上述6種分類器與2種分類數據組合均能實現枸杞種植區的準確提取。文中分類結果按分類精度由高至低依次為： LightGBM+EVI，GBDT+EVI，XGBoost+EVI/GBDT+NDVI，LightGBM+NDVI，XGBoost+NDVI，RF+EVI，RF+NDVI，SVM+EVI，MLPC+NDVI，SVM+NDVI和MLPC+EVI。基于6種分類器對采用NDVI/EVI時序數據下的平均分類精度進行對比分析， LightGBM和GBDT的平均OA均為90.40%，Kappa系數均為0.89，而分類效果最差的MLPC，其平均OA也達到84.06%，Kappa系數為0.82，表明6種分類器對枸杞種植區及其他地物類型均有較強區分能力，按照分類器精度由高至低依次為： LightGBM/GBDT，XGBoost，RF，SVM和MLPC。對研究區不同地物類型的精度進行分析，在12種分類結果中幼齡期枸杞的平均PA為91.08%，壯齡期枸杞的平均PA為94.17%，表明選用的分類器均能實現枸杞種植區的精準提取，通過混淆矩陣分析也存在壯齡期枸杞種植區易于與沼澤地和行道樹類型發生混淆，而幼齡期枸杞易于與戈壁灘產生錯分。研究區內行道樹、公路、戈壁灘和沼澤地的平均PA也較高，分別達到了86.58%，88.25%，85.00%和85.42%，而6種分類器對城鎮的提取效果普遍不佳，平均PA僅為79%，農場城鎮規模較小且散落分布，并且其圖像特征也與公路、戈壁灘類型較為接近，當影像空間分辨率不足以及分類器區分能力較差時易發生上述類型邊界提取模糊和錯分現象。比照NDVI和EVI這2種時序數據的分類效果，在同一分類器下，使用EVI時序數據的OA將高出1.09～2.54個百分點，同時在對壯齡期枸杞、行道樹和沼澤地等類型進行提取時，采用EVI時序數據的PA較采用NDVI時序數據時的PA平均高1.70個百分點。

(a) 基于時序NDVI的分類器制圖精度(b) 基于時序EVI的分類器制圖精度

(c) 基于時序NDVI的分類器用戶精度(d) 基于時序EVI的分類器用戶精度

(a) 6種分類器總體分類精度(b) 6種分類器Kappa系數

3.4 分類時相優選與評價

通過模型特征重要性評分方法分析了不同時相NDVI/EVI數據對于分類的貢獻程度，并對優選后的部分時相特征進行分類和精度驗證，結果見圖10，其中NDVI-OA與NDVI-Kappa分別為分類器特征選取前模型的OA和Kappa系數，NDVI-OA-S與NDVI-Kappa-S分別為分類器特征選取后模型的OA和Kappa系數，分析發現GBDT+EVI，XGBoost+NDVI，SVM+NDVI，SVM+EVI和RF+EVI共5種分類器與數據組合方式下，優選5～6個關鍵時相的NDVI/EVI數據能獲得與采用全部時相NDVI/EVI數據相同的分類精度，并且GBDT+EVI和XGBoost+NDVI在選取6個特征后的OA均達到89.85%，Kappa系數為0.88，對幼齡期枸杞和壯齡期枸杞的PA分別為94.0%和96.0%，表明通過特征優選能在保證分類精度的同時進一步降低數據冗余。LightGBM通過特征重要性評估能夠選擇更少的時相特征，但對比特征優選前后的OA下降了3～8個百分點。

(a) 特征選取前后總體分類精度(b) 特征選取前后Kappa系數

ofmodelsbeforeandaftertheselectionoftemporalfeatures

對分類器選取的時相次數進行統計分析(圖11)，其中6月下旬和9月下旬時相被選取次數最多，該時段分別處于枸杞開花盛期和枸杞秋果成熟盛期，應為枸杞種植區分類的關鍵時相。同時通過對比NDVI和EVI這2種數據源的時相特征選取分布規律，針對EVI數據選取的特征時相集中分布在6—9月，而NDVI時序數據選取的時相特征集中分布在5—8月。

(a) NDVI(b) EVI

圖11 時序NDVI/EVI的時相選擇次數統計

4 討論

中高分辨率遙感數據為作物識別提供了豐富數據源，同時蓬勃發展的集成學習方法將推動作物識別精度的提升。首先，文中綜合GF-1 WFV和Landsat8 OLI數據進行互補以獲取時序較為完整的植被指數數據，前期宋軍偉等[40]已對上述2種衛星的反射率數據進行對比，表明其對應波段均存在顯著相關，并且對典型地物類型的可分離性也十分接近。本文在此研究基礎上進一步對其計算的植被指數進行差異分析，具體是采用相同目標區域的同日過境數據，并對比分析了影像預處理前后2種數據提取的NDVI/EVI相關程度(圖12)，結果表明經影像預處理后提取的GF-1 NDVI/EVI與Landsat8 NDVI/EVI間的相關性被大幅提高，其決定系數由0.73～0.74上升至0.97～0.98，表明上述2種傳感器的植被指數具有較高的一致性，并且在數據互補使用過程中開展輻射定標、大氣校正、幾何精糾正和異常值剔除等影像預處理能夠提高數據間的一致性，而在后續研究中可通過建立GF-1和Landsat8衛星對應波段間的轉換關系以更好地消除數據源間的差異； 同時，為獲取準確且靈敏的植被指數，選用同類研究中普遍使用的NDVI和EVI指數，通過對比發現EVI指數在枸杞生育期的中后期更為靈敏，并與行道樹等高覆蓋植被類型間的時序曲線差異更明顯，這與白燕英等[14]認為在生育期高峰期和高植被覆蓋下EVI指數更加靈敏的結論較為一致。并且為進一步分析NDVI和EVI這2種時序數據對作物分類精度的影響，采用6種分類器進行評估，結果表明在同一分類器下使用EVI時序數據的分類精度普遍高出1.09～2.54個百分點，EVI指數用于作物分類具有更好的區分性； 最后，選用LightGBM，GBDT，XGBoost，RF，SVM和MLPC 6種分類器和NDVI/EVI這2種時序數據開展枸杞種植區分類，各分類器的平均OA分別為90.4%，90.4%，89.31%，86.96%，85.14%和84.06%。表明4種集成學習分類器相較于SVM和MLPC分類器，可獲取更高的作物識別精度。目前，集成學習已在土地利用[41]、葉綠素反演[35]和地上生物量估算[36]方面進行應用并獲得較好精度，而本研究則表明LightGBM，GBDT，XGBoost和RF分類器在枸杞作物分類中也具有較好適用性。

(a) 預處理前NDVI點密度 (b) 預處理后NDVI點密度

(c) 預處理前EVI點密度 (d) 預處理后EVI點密度

5 結論

綜合應用Landsat8 OLI和GF-1 WFV影像，構建作物生育期內時序NDVI/EVI數據，并采用6種分類器對枸杞種植區進行分類，主要結論如下：

1)LightGBM，GBDT，XGBoost和RF集成學習分類器在枸杞種植區分類中能夠獲得更高精度，6種分類器精度由高至低依次為： LightGBM/GBDT，XGBoost，RF，SVM和MLPC，在分類結果中以LightGBM+EVI的OA值最高達到91.67%，Kappa值為0.90，其對幼齡期枸杞和壯齡期枸杞的PA值均達到91%以上。

2)時序NDVI/EVI數據均能較好區分枸杞作物以及種植區內其他地物類型，但EVI指數在枸杞作物生長中后期更為靈敏，并在相同分類器下使用EVI時序數據能夠獲得更好的枸杞制圖精度。

3)選用GBDT，XGBoost和RF分類器的特征重要性評分方法進行NDVI/EVI數據分類時相特征優選后，能夠保證在分類精度不損失的同時降低數據冗余提高分類效率。

4)文中僅選用單一年份的4—11月影像構建時序植被指數數據，而在后期研究中應考慮使用多年份的時序數據以進一步提高分類結果的可靠性。