基于Google Earth Engine與多源遙感數(shù)據(jù)的 土地覆蓋變化監(jiān)測

熊元康，張鴻輝*，梁宇哲，羅偉玲，洪 良

（1.廣東國地規(guī)劃科技股份有限公司，廣東 廣州 510650；2.廣東省土地調(diào)查規(guī)劃院，廣東 廣州 510075）

土地利用/覆被變化（LUCC）記錄了人類在地球表面的空間格局活動，是導(dǎo)致生物多樣性減少、氣候變化、生態(tài)環(huán)境演變、生物化學(xué)循環(huán)乃至全球變化的主要因素[1-2]?？焖?、準(zhǔn)確、全方位地獲取一個地區(qū)的土地覆蓋變化信息，可為該地區(qū)的社會經(jīng)濟(jì)發(fā)展、生態(tài)環(huán)境建設(shè)、國土空間規(guī)劃等提供重要的基礎(chǔ)數(shù)據(jù)。近年來，遙感技術(shù)因其具有監(jiān)測范圍廣、實時性強以及使用成本低等優(yōu)點已被廣泛應(yīng)用于土地覆蓋變化監(jiān)測中[3]。

目前，根據(jù)變化監(jiān)測對象的不同，基于遙感影像的土地覆蓋變化監(jiān)測方法可分為3大類：①基于像素級的土地覆蓋變化監(jiān)測方法，直接利用像元的光譜特征值，通過多時相的遙感數(shù)據(jù)進(jìn)行影像間的差值、比值等代數(shù)運算，進(jìn)而獲取差異影像，再利用經(jīng)驗或通過樣本集獲得的閾值進(jìn)行土地覆蓋變化監(jiān)測，如于冰[4]等通過對GlobeLand30數(shù)據(jù)進(jìn)行重編碼等操作，提出了一種基于像元轉(zhuǎn)換的土地覆蓋變化監(jiān)測方法；②基 于特征級的土地覆蓋變化監(jiān)測方法，通過不同地類所表現(xiàn)出來的光譜特征、紋理特征等進(jìn)行土地覆蓋變化監(jiān)測，與基于像素級的土地覆蓋變化監(jiān)測方法相比，該方法不易受遙感影像時相變化的影響，差異信息更突出、監(jiān)測精度更高，是目前土地覆蓋變化監(jiān)測的主流方法之一，如Jose D A[5]等結(jié)合Landsat數(shù)據(jù)以及L波段SAR數(shù)據(jù)計算得到各類指數(shù)（NDVI、EVI等），再利用隨機(jī)森林分類器成功提取了緬甸南部Tanintharyi地區(qū)的土地覆蓋動態(tài)變化信息；③基于對象級的土地覆蓋變化監(jiān)測方法，主要針對高分辨率的遙感影像，利用多尺度分割生成不同的基元對象，再結(jié)合基元對象的特征進(jìn)行土地覆蓋變化監(jiān)測，與前 兩類方法相比，該方法能獲取更加豐富的特征信息，便于提高分類精度，如吳田軍[6]等基于SPOT 4的多時相數(shù)據(jù)，利用對象級的土地覆蓋變化監(jiān)測方法成功提取了廣東省東莞市東北區(qū)域2005－2008年的土地覆蓋變化信息。

目前，利用中高分辨率遙感影像在我國南方地區(qū)進(jìn)行大范圍土地覆蓋變化動態(tài)監(jiān)測的研究較少（四季多雨，無云數(shù)據(jù)更加難以獲取），且多集中于利用統(tǒng)計數(shù)據(jù)與土地利用數(shù)據(jù)相結(jié)合的方法進(jìn)行土地覆蓋變化動態(tài)監(jiān)測[7]。為了解決大范圍土地覆蓋變化監(jiān)測所面臨的低分辨率遙感影像混合像元嚴(yán)重、中高分辨率遙感影像監(jiān)測范圍小和重返周期長等問題，本文以廣東省為研究區(qū)，基于Google Earth Engine （GEE）云平臺，結(jié)合中等分辨率的Sentinel-1/2和 Landsat 7/8等多源遙感數(shù)據(jù)，通過構(gòu)建多元時間序列影像的方式進(jìn)行土地覆蓋變化動態(tài)監(jiān)測。GEE云平臺是一個能在大尺度范圍下進(jìn)行遙感數(shù)據(jù)處理和分析的云平臺，提供了一個完備的集成環(huán)境。目前，國內(nèi)外多個研究組已基于該平臺開展了各種對地觀測研究，如水稻遙感制圖[8]、耕地遙感提取[9]、農(nóng)作物種植結(jié)構(gòu) 提取[10]等。

1 研究區(qū)概況與數(shù)據(jù)來源

1.1 研究區(qū)概況

廣東省地處我國大陸最南部，東鄰福建，北接江西、湖南，南臨南海，西連廣西；屬于熱帶或亞熱帶季風(fēng)氣候區(qū)。廣東省以廣州市為中心，東起潮州市、西至雷州市，橫跨21個地級市（圖1）；陸地面積為17.97萬km2；地勢表現(xiàn)為北高南低；地貌類型復(fù)雜，北部、東北部和西部均有較高山脈，中部和南部沿海地區(qū)多為低丘、臺地或平原[11]。隨著廣東省各地區(qū)城市化和工業(yè)化的快速推進(jìn)，耕地和林地不斷被蠶食，土地可墾率降低，土地后備資源嚴(yán)重不足，導(dǎo)致人地矛盾日益突出[12]。因此，快速、準(zhǔn)確、全方位地獲取廣東省的土地覆蓋變化信息，可為其社會經(jīng)濟(jì)發(fā)展、生態(tài)環(huán)境建設(shè)、國土空間規(guī)劃等提供重要的基礎(chǔ)數(shù)據(jù)。

圖1 研究區(qū)位置與驗證樣本點的空間分布 （審圖號：GS(2016)2929）

1.2 數(shù)據(jù)來源

1.2.1 遙感數(shù)據(jù)

由于廣東省四季多雨，無云數(shù)據(jù)難以獲取，且中高分辨率遙感影像監(jiān)測范圍小、重返周期長，采用單一數(shù)據(jù)源的中高分辨率影像不能滿足大范圍土地覆蓋變化監(jiān)測的需求，因此本文結(jié)合中高分辨率的 Sentinel-1/2和Landsat7/8等多源遙感數(shù)據(jù)進(jìn)行研究區(qū)內(nèi)的土地覆蓋變化監(jiān)測。遙感影像數(shù)量如表1所示。

表1 本文采用的遙感影像數(shù)量

1.2.2 地面參考數(shù)據(jù)

本文結(jié)合野外實地考察、航拍影像（0.2 m）、Google高清影像、Planet高清影像（5 m）以及土地利用數(shù)據(jù)，將研究區(qū)內(nèi)的土地覆蓋類型分為6大類，如表2所示，并根據(jù)上述數(shù)據(jù)集獲得訓(xùn)練數(shù)據(jù)集。

表2 土地覆蓋類型分類表

1.3 數(shù)據(jù)預(yù)處理

本文利用Sentinel-2、Landsat 8、Landsat 7自身的質(zhì)量評價波段（QA）進(jìn)行相應(yīng)的去云處理。同時，由于不同傳感器在設(shè)計上的差異，導(dǎo)致不同傳感器具有不同的光譜響應(yīng)函數(shù)（SRF），如Sentinel-2和Landsat 7數(shù)據(jù)在紅光波段的均方根誤差（RMSE）超過了8%[13]，因此為了減小不同光譜響應(yīng)函數(shù)帶來的誤差，本文以Landsat 8數(shù)據(jù)為基礎(chǔ)，對Sentinel-2和Landsat 7數(shù)據(jù)相應(yīng)波段進(jìn)行相關(guān)的線性轉(zhuǎn)換[14-15]，并將所有數(shù)據(jù)的分辨率重采樣至30 m。

2 研究方法

2.1 多元時序特征選擇

土地覆蓋變化監(jiān)測與分析主要的依據(jù)是遙感影像時間序列的變化特征信息，因此在構(gòu)建時間序列時，時序特征的選擇顯得尤為重要。通過對研究區(qū)內(nèi)土地覆蓋變化類型進(jìn)行分析發(fā)現(xiàn)，研究區(qū)內(nèi)主要的土地覆蓋變化過程為：耕地或林地轉(zhuǎn)化為建設(shè)用地、耕地轉(zhuǎn)化為魚塘或魚塘轉(zhuǎn)化為耕地、耕地轉(zhuǎn)化為園林地等。鑒于此，本文選擇的時序特征為：

1）歸一化植被指數(shù)（NDVI）。NDVI常被作為特征參數(shù)來評估地表植被的生長狀況[16]，能反映紅光波段（植物吸收強烈）與近紅外波段（植被反射強烈）之間的關(guān)系，因此能良好地區(qū)分植被與建筑用地等。其計算公式為[17]：

式中，ρNIR為近紅外波段的反射率；ρRed為紅光波段的反射率。

2）歸一化建筑指數(shù)（NDBI）。NDBI常被用于城鎮(zhèn)建設(shè)用地提取，能反映城鎮(zhèn)建設(shè)用地在近紅外和中紅外波段具有高反射率值，而其他地物具有低反射率值的趨勢。其計算公式為[18]：

式中，ρMIR為中紅外波段的反射率；ρNIR為近紅外波段的反射率。

3）改進(jìn)的歸一化差異水體指數(shù)（MNDWI）。MNDWI是徐涵秋[19]為了改進(jìn)NDWI[20]提取水體信息的局限性（影像中有建筑陰影的水體提取效果不佳）而提出的，能明顯增強水體與建筑物的反差，減少背景誤差，因此能良好地區(qū)分水體與建筑等。其計算公式為：

式中，ρGreen為綠光波段的反射率；ρMIR為中紅外波段的反射率。

4）SAR影像時序特征。針對研究區(qū)內(nèi)長時間序列無云數(shù)據(jù)難以獲取等問題，結(jié)合合成孔徑雷達(dá)（SAR）不受云雨霧等自然條件影響以及全天候等特性，利用不同地類的不同散射機(jī)理[21]來區(qū)分水體、植被以及建筑等土地覆蓋類型。本文將Sentinel-1的VV極化數(shù)據(jù)作為研究區(qū)土地覆蓋變化監(jiān)測的雷達(dá)數(shù)據(jù)。為了將VV極化數(shù)據(jù)與NDVI、NDBI以及MNDWI數(shù)據(jù)更好地結(jié)合起來，本文利用式（4）對其進(jìn)行歸一化處理。

式中，maxρVV為研究區(qū)內(nèi)SAR長時序數(shù)據(jù)上的最大VV極化散射系數(shù)；minρVV為研究區(qū)內(nèi)SAR長時序數(shù)據(jù)上的最小VV極化散射系數(shù)；ρVV為研究區(qū)內(nèi)長時序上某一時期的VV極化散射系數(shù)。

2.2 分類特征光譜響應(yīng)分析

根據(jù)選取的時序特征，通過Google的高清影像，本文選取研究區(qū)內(nèi)典型的土地覆蓋變化類型進(jìn)行時序光譜響應(yīng)分析，結(jié)果如圖2所示，可以看出，在耕地轉(zhuǎn)化為建設(shè)用地的過程中，NDVI和VV后向散射系數(shù)在不斷下降，而MNDWI和NDBI在地類表現(xiàn)為耕地時，隨時間的變化出現(xiàn)了較大的浮動，在地類轉(zhuǎn)化為建設(shè)用地時，隨時間的變化表現(xiàn)較為平穩(wěn)；在林地轉(zhuǎn)化為建設(shè)用地的過程中，NDVI和MNDWI隨時間的變化呈下降趨勢，而NDBI和VV后向散射系數(shù)隨時間變化呈逐步上升趨勢；在林地轉(zhuǎn)化為裸地的過程中，NDVI、MNDWI和VV后向散射系數(shù)隨時間變化呈逐步下降趨勢，而NDBI呈逐步上升趨勢；在水體轉(zhuǎn)化為裸地的過程中，NDVI、NDBI和VV后向散射系數(shù)隨時間變化呈上升趨勢，而MNDWI呈下降 趨勢。

圖 2 各種土地覆蓋類型的時間序列

2.3 土地覆蓋類型提取方法

本文分別利用隨機(jī)森林分類（RF）、支持向量機(jī)（SVM）以及回歸分類決策樹（CART）進(jìn)行土地覆蓋類型提取，其中RF和SVM已被廣泛應(yīng)用于土地覆蓋類型分類中，并取得了良好的分類結(jié)果，而CART易于表達(dá)且能很好地解釋某些特定的規(guī)則[22]。

RF是一種集成學(xué)習(xí)方法，采用Bootstrap抽樣技術(shù)從原始數(shù)據(jù)集中抽取訓(xùn)練集，并通過訓(xùn)練集構(gòu)建CART決策樹。影響RF分類器性能和效率的主要參數(shù)為決策樹數(shù)量、候選特征子集以及葉節(jié)點最小樣本數(shù)。本文將這3個參數(shù)設(shè)置為：RF分類器的決策樹數(shù)量=100；候選特征子集=4；葉節(jié)點最小樣本數(shù)=1。SVM是一種非參數(shù)機(jī)器學(xué)習(xí)算法，核心是找到一個最優(yōu)的超平面作為高維空間中的決策函數(shù)，進(jìn)而將輸入向量分成為不同的類別。影響SVM性能和效率的主要參數(shù)包括核函數(shù)的選擇和cost參數(shù)。本文選擇線性函數(shù)作為SVM的核函數(shù)，并將cost參數(shù)設(shè)置為100。CART是一個在二進(jìn)制遞歸分區(qū)過程中成長的決策樹，通過子集內(nèi)變量的最大方差和最小方差將訓(xùn)練數(shù)據(jù)集分成不同的類別。樹的最大深度參數(shù)決定了CART模型的復(fù)雜性，大的深度可能具有更高的精度，但也會增加過度擬合的風(fēng)險。綜合以往的研究和實際需求，本文將樹的深度參數(shù)設(shè)置 為10。

2.4 精度評價

對遙感影像分類結(jié)果進(jìn)行精度評價是一個十分重要的工作[10]。為了保證精度驗證的準(zhǔn)確性，本文選擇總體精度（OA）和F-Score作為土地覆蓋變化監(jiān)測結(jié)果的精度評價指標(biāo)。這兩個指標(biāo)均來源于混淆矩陣，其中OA用以評價整體算法的有效性，而F-Score用以評價每一類的分類精度。

3 研究結(jié)果與分析

3.1 精度分析

為了對土地覆蓋類型變化監(jiān)測結(jié)果進(jìn)行精度評價，本文通過分層隨機(jī)采樣在研究區(qū)內(nèi)生成了1 000個驗證樣本點，如圖1所示，并通過人工目視解譯判讀其是否為土地覆蓋類型發(fā)生變化的區(qū)域，進(jìn)而獲得587個非變化樣本點和413個變化樣本點。由精度評價結(jié)果可知（表3），RF方法變化監(jiān)測的OA為83.5%，F(xiàn)-Score為0.82；SVM方法變化監(jiān)測的OA為73.5%，F(xiàn)-Score為0.58；CART方法變化監(jiān)測的OA為70.9%，F(xiàn)-Score為0.51。因此，RF表現(xiàn)最優(yōu)，SVM次之，而CART結(jié)果最不理想。

表3 分類精度統(tǒng)計

3.2 土地覆蓋類型提取結(jié)果

根據(jù)多元分類特征，在GEE云平臺上分別融合生成2017年、2018年研究區(qū)內(nèi)的多元分類特征影像（49個特征波段、12個NDVI時間序列波段、12個MNDWI時間序列波段、12個NDBI時間序列波段、12個VV后向散射系數(shù)時間序列波段，1個Slope波段），并同時利用收集到的訓(xùn)練數(shù)據(jù)集進(jìn)行RF、SVM以及CART訓(xùn)練，以建立研究區(qū)內(nèi)不同土地覆蓋類型的監(jiān)測模型。在此基礎(chǔ)上，利用訓(xùn)練好的RF、SVM以及CART分類模型進(jìn)行2017年、2018年研究區(qū)內(nèi)的土地覆蓋類型監(jiān)測，結(jié)果如圖3所示。

圖3 廣東省2017－2018年土地覆蓋類型監(jiān)測結(jié)果

3.3 土地覆蓋變化提取結(jié)果

根據(jù)2017年、2018年廣東省土地覆蓋類型分類結(jié)果，本文采用基于空間疊置分析和統(tǒng)計原理分析的分類后變化檢測方法進(jìn)行以a為單位的土地覆蓋變化監(jiān)測，結(jié)果如圖4、5所示。在此基礎(chǔ)上，本文對研究區(qū)內(nèi)的土地覆蓋變化區(qū)域的空間分布狀況進(jìn)行統(tǒng)計發(fā)現(xiàn)，研究區(qū)內(nèi)土地覆蓋類型發(fā)生變化的區(qū)域主要集中在珠三角地區(qū)（廣州市等9個地區(qū)），其次為粵東地區(qū)（潮州市等4個地區(qū)），最后為粵西地區(qū)（湛江市等3個地區(qū)）和粵北山區(qū)（清遠(yuǎn)市等5個地區(qū)）。由于珠三角地區(qū)人口的集聚和工業(yè)的發(fā)展，其土地覆蓋變化類型主要為耕地或林地轉(zhuǎn)化為建設(shè)用地或工業(yè)用地；由于粵東部分地區(qū)人口眾多，其土地覆蓋變化類型主要為耕地或林地轉(zhuǎn)化為建設(shè)用地；由于粵西地區(qū)和粵北山區(qū)是重要的農(nóng)業(yè)生產(chǎn)空間，其土地覆蓋變化類型主要為耕地轉(zhuǎn)化為魚塘或魚塘轉(zhuǎn)化為耕地，以及耕地轉(zhuǎn)化為園林地。

圖4 廣東省2017－2018年土地覆蓋變化監(jiān)測結(jié)果

圖5 珠江口2017－2018年土地覆蓋變化監(jiān)測結(jié)果

為了發(fā)掘影響土地覆蓋變化的潛在因素，結(jié)合2004－2018年廣東省的各類統(tǒng)計數(shù)據(jù)，本文對廣東省的土地覆蓋變化進(jìn)行分析發(fā)現(xiàn)，影響廣東省土地覆蓋變化空間分布的主要原因為人口增長、區(qū)域社會經(jīng)濟(jì)發(fā)展以及政策變化，如圖6所示。在人口增長方面（圖6a），珠三角地區(qū)人口年均增長99.2萬人，粵東地區(qū)人口年均增長4.6萬人，粵西地區(qū)人口年均增長 7.8萬人，粵北山區(qū)人口年均增長8.8萬人；在區(qū)域社會經(jīng)濟(jì)發(fā)展方面（圖6b），珠三角地區(qū)的年均生產(chǎn)總值增長為5 031億元，粵東地區(qū)的年均生產(chǎn)總值增長為386億元，粵西地區(qū)的年均生產(chǎn)總值增長為 438億元，粵北山區(qū)的年均生產(chǎn)總值增長為281億元； 在政策變化方面，主要體現(xiàn)為政府通過制定發(fā)布相關(guān)的政策來干預(yù)、調(diào)整土地的使用（如加強基礎(chǔ)設(shè)施建設(shè)等），由研究區(qū)內(nèi)各地區(qū)的固定資產(chǎn)投資（基本建設(shè)投資、更新改造投資、房地產(chǎn)開發(fā)投資）情況可知（圖6c），珠三角地區(qū)的年均固定資產(chǎn)投資增長為 1 833億元，粵東地區(qū)的年均固定資產(chǎn)投資增長為 397億元，粵西地區(qū)的年均固定資產(chǎn)投資增長為290億元， 粵北山區(qū)的年均固定資產(chǎn)投資增長為242億元，進(jìn)而導(dǎo)致研究區(qū)內(nèi)建設(shè)用地與耕地面積比值的不斷增加，即建設(shè)用地面積不斷增加，而耕地面積不斷減 少（圖6d）。

圖6 影響廣東省土地覆蓋變化的潛在因素

4 結(jié) 語

本文以廣東省為研究區(qū)，通過分析研究區(qū)內(nèi)的主要土地覆蓋變化類型，選取NDVI、MNDWI、NDBI、VV后向散射系數(shù)的時序影像以及DSM影像作為研究區(qū)土地覆蓋類型監(jiān)測的特征波段影像，同時結(jié)合訓(xùn)練數(shù)據(jù)集、RF、SVM以及CART對研究區(qū)內(nèi) 2017年、2018年的土地覆蓋類型進(jìn)行了監(jiān)測，并在此基礎(chǔ)上進(jìn)行了研究區(qū)內(nèi)的土地覆蓋變化監(jiān)測。通過對土地覆蓋變化監(jiān)測結(jié)果進(jìn)行分析，本文得出的主要結(jié) 論為：

1）構(gòu)建的49個分類特征波段以及RF分類器能較好地適應(yīng)于研究區(qū)內(nèi)的土地覆蓋類型監(jiān)測。較SVM和CART而言，RF具有最高的土地覆蓋變化監(jiān)測精度，其土地覆蓋變化監(jiān)測OA為83.5%，F(xiàn)-Score為 0.82。

2）對研究區(qū)內(nèi)土地覆蓋變化監(jiān)測結(jié)果進(jìn)行分析發(fā)現(xiàn)，研究區(qū)內(nèi)土地覆蓋發(fā)生變化的區(qū)域主要集中在珠三角地區(qū)，其土地覆蓋變化類型主要為耕地或林地轉(zhuǎn)化為建設(shè)用地或工業(yè)用地；其次為粵東地區(qū)，其土地覆蓋變化類型主要為耕地或林地轉(zhuǎn)化為建設(shè)用地；最后為粵西地區(qū)和粵北山區(qū)，其土地覆蓋變化類型主要為耕地轉(zhuǎn)化為魚塘或魚塘轉(zhuǎn)化為耕地，以及耕地轉(zhuǎn)化為園林地。

3）結(jié)合統(tǒng)計年鑒等相關(guān)數(shù)據(jù)，對研究區(qū)內(nèi)的土地覆蓋變化進(jìn)行分析發(fā)現(xiàn)，影響研究區(qū)內(nèi)土地覆蓋變化空間分布的主要原因為人口增長、區(qū)域社會經(jīng)濟(jì)發(fā)展和政策變化。