基于GEE的東南亞主產區(qū)橡膠林分布遙感提取

李偉光 張京紅 劉少軍 陳小敏 鄒海平 白蕤 呂潤

摘 ?要：為遙感監(jiān)測全球橡膠主產區(qū)橡膠的長勢，選取橡膠產量最大的泰國、馬來西亞、印度尼西亞三國，開展基于Google Earth Engine（GEE）的橡膠林分布遙感提取研究。通過目視解譯選擇典型樣本區(qū)，根據Landsat 7多波段光譜特征、MODIS NDVI反映的植被物候特征建立分類回歸樹CART分類模型提取橡膠林分布。精度評價顯示模型總體分類精度為95.8%，Kappa系數(shù)為0.94，生產者精度達到94.8%，用戶精度為88.2%，達到較高水平。模型提取結果顯示：橡膠林在泰國中部、南部半島，馬來半島的東部和南部地區(qū)，蘇門答臘島分布較為集中，而泰國北部、加里曼丹島及其他島嶼橡膠林相對稀疏。

關鍵詞：東南亞;橡膠;分布;Google Earth Engine

中圖分類號：S29 ? ? ?文獻標識碼：A

Abstract： In order to monitor the growth of rubber forest in the main rubber foresttree producing areas of the world by remote sensing， we carried out the remote sensing extraction of rubber forest distribution based on Google Earth engine （GEE） in Thailand， Malaysia and Indonesia， which are the three largest rubber producing countries in the world. The typical sample areas were selected through visual interpretation， and the classification regression tree CART classifica-tion model was established based on the multi-band spectral albedo of Landsat 7 and the vegetation phenological cha-racteristics reflected by MODIS NDVI. The overall classification accuracy of the model was 95.8%， the Kappa coeffi-cient was 0.94， the producer accuracy was 94.8%， and the user accuracy was 88.2%， reaching a high level. The extrac-tion results show that rubber plantations are concentrated in central and southern peninsula of Thailand， the eastern and southern parts of Malay Peninsula and Sumatra Island， while rubber plantations are relatively sparse in northern Thailand， Kalimantan Island and other islands.

Keywords： Southeast Asia; rubber; distribution; Google Earth Engine

DOI： 10.3969/j.issn.1000-2561.2021.11.037

天然橡膠不僅是重要的工業(yè)原料，更是國防和工業(yè)建設不可或缺的戰(zhàn)略能源。橡膠樹作為重要的經濟作物在熱帶地區(qū)農民增收和經濟發(fā)展中發(fā)揮著越來越重要的作用。東南亞地區(qū)產出的橡膠占全球總產量90%以上;其中泰國、印度尼西亞、馬來西亞三國占全球割膠面積的72%[1]。掌握該地區(qū)橡膠種植、長勢等情況，對我國政府部門及時制定或調整產業(yè)發(fā)展政策，保障天然橡膠資源供應安全，維護橡膠產業(yè)的健康發(fā)展具有重要意義。利用遙感技術提取天然橡膠林種植的空間分布信息，獲取橡膠林種植的時空變化特征，是開展橡膠林長勢、災害、產量遙感監(jiān)測的重要前置條件[2-4]。

在遙感提取國內橡膠林種植分布方面，張京紅等[5]利用Landsat-TM衛(wèi)星影像采用監(jiān)督分類的方法提取了海南島2008年天然橡膠種植面積信息。田光輝等[6]利用MODIS EVI數(shù)據構建橡膠樹物候特征參數(shù)提取了海南天然橡膠林分布情況。楊紅衛(wèi)等[7]利用高分辨率遙感影像紋理和多光譜特征提取了海南島某農場的精細化分布。在我國另一橡膠主產區(qū)云南西雙版納地區(qū)，廖諶婳等[8]、余凌翔等[9]、Senf等[10]也分別通過不同方法開展了橡膠林分布提取，并分析了橡膠林的擴張情況[11]。針對我國境外橡膠林分布情況，李陽陽等[12]利用MODIS數(shù)據及橡膠林的物候特征提取了老撾北部地區(qū)橡膠林分布及擴張情況。李宇宸等[13]應用決策樹方法對Landsat OLI多時相遙感影像數(shù)據表征的橡膠樹物候特征提取了中老緬交界區(qū)橡膠林分布情況。從提取采用的數(shù)據和分類依據來看，多數(shù)研究采用MODIS中分辨率影像多時相植被指數(shù)作為區(qū)分橡膠林與天然林的差異特征[12-17]，這種方法適用于大范圍橡膠林分布提取;還有采用雷達衛(wèi)星[18]、高分辨率遙感影像紋理[4， 8-9]作為分類依據，提取小范圍、精細化分布情況的報道。以上研究主要針對我國海南、云南和周邊國家橡膠種植區(qū)，而對橡膠種植最多的泰國、印度尼西亞、馬來西亞地區(qū)分布提取的研究較少。

東南亞橡膠主產區(qū)（泰國、印度尼西亞、馬來西亞）空間范圍大、遙感影像產品數(shù)據量巨大，傳統(tǒng)的數(shù)據下載、單機運算成為提取該地區(qū)橡膠林分布情況的瓶頸。因此本研究利用Google Earth Engine（GEE）遙感數(shù)據處理云平臺[19]，以Landsat OLI及MODIS NDVI數(shù)據作為數(shù)據源，融合光譜、物候2種特征作為分類依據，通過云計算技術解決大尺度、長時序的海量遙感數(shù)據處理問題，快速、準確提取東南亞地區(qū)橡膠林分布。

1 ?材料與方法

1.1 ?研究區(qū)概況

研究區(qū)包括橡膠種植面積最廣、產量最高的泰國、印度尼西亞、馬來西亞三國，位于北緯20～南緯10°、東經96～140°之間。泰國地處中南半島，為熱帶季風氣候，年均氣溫24～30 ℃，常年溫度不低于18 ℃，平均年降水量約1000 mm;11月至次年2月受較涼的東北季風影響比較干燥，3—5月氣溫最高，可達40～42 ℃;7—10月受西南季風影響，是雨季;農作物一般在雨季播種，旱季收獲。印度尼西亞和馬來西亞地處馬來群島，屬熱帶雨林氣候，終年高溫多雨，年平均溫度25～ 27 ℃，無四季分別，北部受北半球季風影響，7～9月降水量豐富，南部受南半球季風影響，12月、1月、2月降水量豐富，年降水量1600～2200 mm。

1.2 ?數(shù)據源及預處理

研究區(qū)地處熱帶，云量較多，遙感影像數(shù)據質量較差。因此選用GEE平臺提供的Landsat 7 2012—2014年大氣層頂影像產品融合數(shù)據集[20]。該數(shù)據集為GEE對經NASA輻射定標、幾何校正、云雪陰影掩膜處理的大氣層頂反射率產品（TOA）用Simple compose方法融合而成。

植被指數(shù)選用MODIS MOD13Q1數(shù)據集（2011—2016年）。該數(shù)據集是由美國國家海洋和大氣管理局先進的中分辨率成像光譜儀（MODIS）生成的歸一化植被指數(shù)（NDVI），空間分辨率為250 m，時間分辨率為16 d[21]。

純凈、典型的樣本是高精度分類的關鍵。為保證具有足夠數(shù)量和質量的樣本點用于分類以及精度驗證，本文利用Google Earth高分辨率影像數(shù)據目視解譯具有代表性、典型性的純凈像元作為樣本點。

1.3 ?方法

2 ?結果與分析

2.1 ?典型橡膠林影像

在高分辨率的遙感影像下，可以清晰發(fā)現(xiàn)橡膠林、大型農作物農田、天然森林具有顯著的差異（圖2）。天然森林植被密度高，亮區(qū)連片、陰影呈散落斑點狀;橡膠林及大型作物農田在種植時被人工排列成整齊的行列，具有明顯的行列紋理特征。橡膠林與大型農作物相比，都成行成列，但大型農作物的農田一般行列間距更大;在冠層形態(tài)上也有顯著不同，橡膠林在一行中更密集。若對橡膠林與大型農作物存在不確定，查閱Google Earth歷史影像即可發(fā)現(xiàn)不同：橡膠林歷史影像變化不大，均為行列種植的大型樹木;農田歷史影像變化差異相對更大，這是由于農作物生長更快、種植的作物更換更加頻繁。其他水體、城鎮(zhèn)用地（裸地）在影像上具有更加清晰明顯的差異。通過目視解譯識別這些差異可以選定典型樣本區(qū)，共選取橡膠樣本像素點107個，天然森林113個，水體30個，人工建筑用地（裸地）240個，農田223個。

2.2 ?不同地表覆蓋類型的光譜信息對比

遙感影像分類過程中首先考慮研究區(qū)目標地表覆蓋類型，本文將地表覆蓋類型劃分為農田類、天然森林類、橡膠林類、水體類及城鎮(zhèn)用地類（圖3）。通過統(tǒng)計典型地表覆蓋純凈像元的光譜信息不同波段的反射率可以發(fā)現(xiàn)：水體和城鎮(zhèn)用地在波段1～3顯著高于植被覆蓋的區(qū)域，而在波段4又低于植被覆蓋區(qū)域。這個差異可以準確區(qū)分植被覆蓋區(qū)和水體、城鎮(zhèn)用地。波段6為亮溫數(shù)據，幾種地表覆蓋類型差距較小，普遍在300 K附近。農田在波段4的反射率明顯小于橡膠林、天然森林，易與天然森林和橡膠林區(qū)分。橡膠林與天然森林二者曲線走勢及數(shù)值都非常相近，單純依靠不同波段的光譜特征區(qū)分橡膠林與天然森林較為困難。因此，需進一步選擇能反映物候差異的NDVI時間序列才能較好識別橡膠。

2.3 ?森林與橡膠林NDVI時間序列分析

不同的植被在不同季節(jié)或生育期表現(xiàn)出不同的生理特征，比如長葉、落葉，這些變化能夠通過多時相植被指數(shù)時間序列的變化曲線來表示。本文提取了上半年12期16 d合成的NDVI數(shù)據中間值（圖4），對比發(fā)現(xiàn)農田在所有時相的植被指數(shù)均小于天然森林及橡膠林;天然森林NDVI值在不同時間普遍穩(wěn)定在0.7～0.8附近;橡膠林在第4～6期有一個低值時段。橡膠林在春季一般有落葉、第1蓬葉抽發(fā)的物候現(xiàn)象。在中國區(qū)域橡膠林一般12月份開始落葉，至2月落葉過程完成，3—4月第1蓬葉抽發(fā)完成。橡膠林落葉、新葉抽發(fā)時間主要受氣溫影響，氣溫越高落葉時間越晚，新葉抽發(fā)速率越快。分析東南亞地區(qū)典型橡膠林NDVI曲線可以發(fā)現(xiàn)，該地區(qū)的橡膠林與中國區(qū)域具有相似的物候特征，在2月下旬至3月間完成落葉、新葉抽發(fā)過程。因此可以通過NDVI時間序列反映出的橡膠林與天然森林不同物候特征區(qū)分二者。

2.4 ?分類精度

基于目視解譯的典型樣本數(shù)據，采用分類回歸樹CART方法對典型樣本進行分類，并根據分類結果建立混淆矩陣（表1）。依據精度評價計算公式（3）、公式（4），該分類結果的總體分類精度為95.8%，Kappa系數(shù)為0.94。從分類精度評價指標來看，這個分類結果精度能夠滿足空間分析與實際應用需求。從橡膠林樣本生產者精度來看，達到93.8%;用戶精度達到88.2%，也達到較高的水平。從橡膠林分類誤差的來源看，主要發(fā)生在農田、森林和橡膠林之間，特別是森林與橡膠林之間?？傮w而言，這個分類模型的精度滿足大范圍提取橡膠林分布的要求。

2.5 ?東南亞橡膠林分布狀況

利用GEE云計算平臺根據上述建立的CART分類模型提取的東南亞地區(qū)（泰國、馬來西亞、印度尼西亞）橡膠林空間分布見圖4。從圖中可以發(fā)現(xiàn)泰國中部、南部半島橡膠林分布較為集中，另外東部小范圍地區(qū)較為集中，其余地區(qū)零星分布。馬來西亞的橡膠林主要分布在馬來半島的東部和南部，而加里曼丹島北部地區(qū)橡膠林相對稀疏。印度尼西亞橡膠林在蘇門答臘島分布較為集中，特別是蘇門答臘島的南部東側，加里曼丹島分布相對稀疏。遙感提取東南亞三國橡膠林分布與三國橡膠產區(qū)的文字描述[23-24]基本一致，該分布圖可以作為開展橡膠林長勢、災害等遙感監(jiān)測的基礎數(shù)據。

3 ?討論

本研究利用GEE云計算平臺，通過目視解譯高分辨率遙感影像選取典型樣本區(qū)，分析Landsat影像和MODIS NDVI時間序列差異，建立CART分類回歸樹分類模型，提取橡膠林分布信息。目視解譯發(fā)現(xiàn)在高分辨率遙感影像下，橡膠林具有獨特的行列紋理特征，以區(qū)別于其他植被覆蓋區(qū)。天然森林與橡膠林在Landsat多波段光譜曲線較為相似，與其他地物特征顯著不同。MODIS NDVI時間序列反映的植被物候特征表明，東南亞地區(qū)的橡膠林與我國境內的類似，在2月下旬至3月間完成落葉、新葉抽發(fā)，相應時段的NDVI有一低值時段。利用以上影像特征建立的CART分類回歸樹模型，分類精度達95.8%。模型提取的橡膠林在泰國中部、南部半島，馬來半島的東部和南部地區(qū)，蘇門答臘島分布較為集中，而泰國北部、加里曼丹島及其他島嶼橡膠林相對稀疏。提取的橡膠林分布信息與文獻[23-24]中橡膠林分布及相關文字描述相吻合。

研究區(qū)域地處熱帶，全年云量大、大氣中水汽含量高，嚴重影響光學遙感影像及產品的質量，很難獲取特定時間段內高質量的大范圍無云影像。為克服缺少光學影像問題，嘗試利用雷達等主動遙感影像來研究土地分類[25-26]。雷達遙感對土壤水分變化較為敏感，用于大范圍分類時，會因土壤水分差異而影響精度，常用于小范圍、高精度提取。橡膠作為一種多年生常綠植物，在空間分布上相對穩(wěn)定。本研究所采用的Landsat 7 3年大氣層頂影像產品融合數(shù)據集和MOD13Q1 NDVI時間序列數(shù)據集，均為多時次遙感影像通過simple compose方法合成的數(shù)據集，反映的是地物光譜特征的中間值，具有一定的穩(wěn)定性，在本分類中具有較好表現(xiàn)。另一個影響橡膠林分布提取精度的原因是農田、天然森林、橡膠林之間的混淆[26]。當森林上空懸浮的薄云未達到云識別閾值時，會造成NDVI值的下降，可能會誤識為橡膠林。東南亞三國森林覆蓋率高，在55%～75%以上，即使有小比例森林辨識為橡膠林，也會造成用戶精度較大下降。在這3個國家中橡膠林占國土面積的比例較小，而且種植的集中程度低于中國的海南及西雙版納，這也給橡膠林分布的提取帶來了困難。未來需要進一步實地調研獲取第一手典型樣方，結合地形、橡膠林的年齡等信息來構建更高精度的分類模型，進一步對東南亞橡膠主產區(qū)開展長勢和產量遙感監(jiān)測，以滿足政策制定、貿易判斷等決策需求。

參考文獻

[1] 莫業(yè)勇， 楊 ?琳. 2019年國內外天然橡膠產銷形勢[J]. 中國熱帶農業(yè)， 2020， 93（2）： 8-12.

[2] 陳小敏， 陳匯林， 李偉光， 等. 海南島天然橡膠林春季物候期的遙感監(jiān)測[J]. 中國農業(yè)氣象， 2016， 37（1）： 111-116.

[3] 劉少軍， 張京紅， 蔡大鑫， 等. Landsat 8在橡膠林臺風災害監(jiān)測中的應用[J]. 自然災害學報， 2016， 25（2）： 53-58.

[4] 張京紅， 張明潔， 劉少軍， 等. 風云三號氣象衛(wèi)星在海南橡膠林遙感監(jiān)測中的應用[J]. 熱帶作物學報， 2014， 35（10）： 2059-2065.

[5] 張京紅， 陶忠良， 劉少軍， 等. 基于TM影像的海南島橡膠種植面積信息提取[J]. 熱帶作物學報， 2010， 31（4）： 661-665.

[6] 田光輝， 李海亮， 陳匯林. 基于物候特征參數(shù)的橡膠樹種植信息遙感提取研究[J]. 中國農學通報， 2013， 29（28）： 46-52.

[7] 楊紅衛(wèi)， 童小華. 高分辨率影像的橡膠林分布信息提取[J]. 武漢大學學報（信息科學版）， 2014， 39（4）： 411-416， 421.

[8] 廖諶婳， 李 ?鵬， 封志明， 等. 西雙版納橡膠林面積遙感監(jiān)測和時空變化[J]. 農業(yè)工程學報， 2014， 30（22）： 170-180.

[9] 余凌翔， 朱 ?勇， 魯韋坤， 等. 基于HJ-1CCD遙感影像的西雙版納橡膠種植區(qū)提取[J]. 中國農業(yè)氣象， 2013， 34（4）： 493-497.

[10] Senf C， Pflugmacher D， Van D L， et al. Mapping rubber plantations and natural forests in Xishuangbanna （Southwest China） using multi-spectral phenological metrics from MODIS time series[J]. Remote Sensing， 2013， 5（6）： 2795- 2812.

[11] Kou W， Dong J， Xiao X， et al. Expansion dynamics of deciduous rubber plantations in Xishuangbanna， China during 2000—2010[J]. GIScience & Remote Sensing， 2018， 55（6）： 905-925.

[12] 李陽陽， 張 ?軍， 劉陳立， 等. 老撾北部5省橡膠林提取及時空擴張研究[J]. 林業(yè)科學研究， 2017， 30（5）： 709-717.

[13] 李宇宸， 張 ?軍， 薛宇飛， 等. 基于Google Earth Engine的中老緬交界區(qū)橡膠林分布遙感提取[J]. 農業(yè)工程學報， 2020， 36（8）： 174-181.

[14] 劉曉娜， 封志明， 姜魯光. 基于決策樹分類的橡膠林地遙感識別[J]. 農業(yè)工程學報， 2013， 29（24）： 163-172， 365.

[15] 陳匯林， 陳小敏， 陳珍麗， 等. 基于MODIS遙感數(shù)據提取海南橡膠信息初步研究[J]. 熱帶作物學報， 2010， 31（7）： 1181-1185.

[16] Kou W， Liang C， Wei L， et al. Phenology-based method for mapping tropical evergreen forests by integrating of MODIS and landsat imagery[J]. Forests， 2017， 8（2）： 34.

[17] Zhai D L， Yu H， Chen S C， et al. Responses of rubber leaf phenology to climatic variations in Southwest China[J]. In-ternational Journal of Biometeorology， 2019， 63（5）： 607-616.

[18] Dong J， Xiao X， Sheldon S， et al. Mapping tropical forests and rubber plantations in complex landscapes by in-tegrating PALSAR and MODIS imagery[J]. ISPRS Journal of Photogrammetry and Remote Sensing， 2012， 74（1）： 20-33.

[19] 韓冰冰， 陳圣波. 基于GEE時間序列遙感影像分類方法研究[J]. 世界地質， 2020， 39（3）： 706-713.

[20] Landsat 7 3-year TOA percentile composites[EB/OL]. （2021-1-1） .https：//developers.google.com/earth-engine/datasets/catalog/LANDSAT_LE7_TOA_3YEAR.

[21] 李偉光， 田光輝， 鄒海平， 等. 海南島典型植被區(qū)EVI特征及其對氣象因子的響應[J]. 中國農學通報， 2014， 30（35）： 190-194.

[22] 于莉莉， 孫立雙， 張丹華， 等. 基于Google Earth Engine的環(huán)渤海地區(qū)土地覆蓋分類[J]. 應用生態(tài)學報， 2020， 31（12）： 4091-4098.

[23] Food and Agriculture Organization of the United Nations. National Agro-Economic Zoning for Major Crops in Thail-and[M]. Rome， 2017.

[24] Kaur A. Plantation systems， labour regimes and the state in Malaysia， 1900-2012[J]. Journal of Agrarian Change， 2014， 14（2）： 190-213.

[25] 陳幫乾， 李香萍， 肖向明， 等. 基于PALSAR雷達數(shù)據與多時相TM/ETM+影像的海南島土地利用分類研究[J]. 熱帶作物學報， 2015， 36（12）： 2230-2237.

[26] Razali S M， Marin A， Nuruddin A A， et al. Capability of integrated MODIS imagery and ALOS for oil palm， rubber and forest areas mapping in tropical forest regions[J]. Sen-sors， 2014， 14（5）： 8259-8282.

責任編輯：謝龍蓮