基于谷歌地球引擎的開封城區(qū)2010—2019年水體分布變化研究

周 珂, 柳 樂, 程承旗, 苗 茹, 楊 陽

(1.河南大學計算機與信息工程學院, 開封 475004； 2.北京大學工學院, 北京 100871; 3.河南大學, 河南省大數(shù)據(jù)分析與處理重點實驗室, 開封 475004; 4.河南大學， 空間信息處理工程實驗室, 開封 475004)

水是生命之源,任何城市的發(fā)展都離不開水域的變化。縱觀歷史,開封市的發(fā)展與水系工程建設息息相關(guān)。開封市作為八大古都之一、中國31個重點防洪城市之一、中原經(jīng)濟區(qū)以及鄭州大都市圈核心城市之一,城市發(fā)展受到文物保護以及生態(tài)保護等多種要求。目前隨著遙感技術(shù)的發(fā)展,遙感技術(shù)以其不受地理空間限制、覆蓋面廣、獲得數(shù)據(jù)速度快等優(yōu)點[1]已成為監(jiān)測水體、檢測水質(zhì)以及研究水體變化的重要手段[2-4],通過遙感影像去分析開封市城區(qū)水體分布的年際變化趨勢,對開封未來的城市發(fā)展以及生態(tài)環(huán)境的保護具有重要意義。

在遙感技術(shù)中常用的提取城市間水體的方法有單波段閾值法、多波段譜間關(guān)系法、水體指數(shù)法等方法[5]。隨著遙感技術(shù)的不斷發(fā)展,機器學習分類方法以其自動化程度高,分類精度較準確以及效率高等優(yōu)點在遙感科學領域取得廣泛應用。目前在大數(shù)據(jù)背景下海量遙感數(shù)據(jù)集的快速處理變得尤為重要,所以利用高性能云計算平臺使用機器學習分類方法對城市水體信息的提取研究很有必要。

谷歌地球引擎(google earth engine,GEE)是谷歌旗下的1款產(chǎn)品,是1個行星尺度的地球科學數(shù)據(jù)與分析平臺,具有強大的計算能力,將數(shù)PB(數(shù)據(jù)存儲單位)級的衛(wèi)星圖像和地理空間數(shù)據(jù)集與行星尺度分析功能結(jié)合在一起,使科學家、研究人員和開發(fā)人員可以使用它檢測地球表面的變化、繪制趨勢圖并量化差異[6]。在數(shù)據(jù)方面,GEE云平臺在云端存儲大量多源影像數(shù)據(jù)供用戶使用,包括Landsat系列、MODIS系列以及Sentinel等常用的影像數(shù)據(jù)集,且用戶可以上傳自己的數(shù)據(jù)在GEE進行計算分析等操作。在語言方面,GEE提供了JavaScript API和Python API, 只需了解JavaScript和python語言的基本語法和數(shù)據(jù)類型即可上手操作。目前,有許多科學工作者運用GEE云平臺進行了大量的地理科學研究工作。Luo等[7]研究了GEE在考古和文化遺產(chǎn)方面的應用；Liu等[8]使用GEE結(jié)合Landsat系列影像數(shù)據(jù)研究了越南玉顯縣1993—2017年的紅樹林年際變化檢測,得出結(jié)論:紅樹林和紅樹林-蝦塘面積呈現(xiàn)出負相關(guān)關(guān)系；裴杰等[9]使用GEE利用近2 000景30 m分辨率Landsat-NDVI長時間序列數(shù)據(jù)研究植被覆蓋度變化長時間序列遙感監(jiān)測。以上研究使用GEE云平臺都取得了較好的效果。

結(jié)合以上分析,利用GEE高性能云計算平臺,結(jié)合2010—2019年的Landsat系列影像數(shù)據(jù)，利用機器學習的方法對開封市城區(qū)的水體分布變化進行研究。

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

1.1 研究區(qū)域概況

選取河南省開封市城區(qū)作為研究區(qū)(圖1)。開封市位于113°5′15″～115°15′42″ E、34°11′45″～35°01′20″N,地處河南省東部。開封市與水有著千絲萬縷的聯(lián)系,境內(nèi)河流湖泊眾多,水資源豐富,分屬黃河和淮河2大水系,雖屬北城,卻擁有“一城宋韻半城水”的贊譽。

圖1 研究區(qū)Fig.1 Research area

1.2 數(shù)據(jù)來源及數(shù)據(jù)處理

1.2.1 Landsat系列影像數(shù)據(jù)

開封地處華北平原,屬溫帶季風氣候,四季分明,4—8月份的影像數(shù)據(jù)含云量較少,能夠較清晰地顯示各地物的地貌特征。為了更好地研究開封市城區(qū)的水體分布變化,選取的影像為4—8月份云量極少(均小于等于10%)的影像數(shù)據(jù)。根據(jù)研究對象的空間特性,研究所使用的影像數(shù)據(jù)來自Landsat系列的Landsat 5、Landsat 7以及Landsat 8影像數(shù)據(jù)。在2010—2011年選取Landsat 5的TM影像作為數(shù)據(jù)源；在2012年選取Landsat 7 ETM+影像作為數(shù)據(jù)源；在2013—2019年選取Landsat 8 OLI影像作為數(shù)據(jù)源。所有操作都在GEE上通過編寫代碼實現(xiàn),所獲得數(shù)據(jù)結(jié)果存儲在Google云盤中。Landsat 5 TM影像和Landsat 8 OLI影像各波段信息如表1、表2所示。

表1 Landsat 5 TM 各波段信息Table 1 Information of each band of Landsat 5 TM

表2 Landsat 8 OLI 各波段信息Table 2 Information of each band of Landsat 8 OLI

1.2.2 樣本數(shù)據(jù)

樣本數(shù)據(jù)均從Google Earth上選取。Google Earth由谷歌公司開發(fā)的一款虛擬地球平臺,可以為用戶提供免費高度清晰的遙感影像[10],用戶可以使用Google Earth查看全球范圍內(nèi)近10年來的高分辨率影像,很適合用來選取長時間序列地物分類的樣本。根據(jù)各地物在衛(wèi)星影像上的可見光反射率不同以及紋理信息以目譯的方法進行選取樣本,同時結(jié)合外業(yè)工作,在2010—2019年共選取1 227份樣本數(shù)據(jù),分布均勻,涵蓋了研究區(qū)的所有地物類型,用于訓練數(shù)據(jù)以及驗證分類結(jié)果。選取的所有樣本保存為kml格式，通過Arcgis軟件轉(zhuǎn)換成矢量文件上傳到GEE云平臺。樣本選取依據(jù)如表3所示。

表3 樣本選取依據(jù)Table 3 Sample selection basis

1.2.3 數(shù)據(jù)預處理

(1)影像數(shù)據(jù)處理。從GEE中下載的Landsat 影像數(shù)據(jù)均為一級校正的TOA遙感影像數(shù)據(jù),在GEE中選擇影像的過程如圖2所示。選取的2010—2019年影像信息如表4所示,選取的影像數(shù)據(jù)均為4—8月份的最佳影像數(shù)據(jù)。其中2012年遙感影像使用的是Landsat 7 ETM+影像數(shù)據(jù),據(jù)USGS官方公布,自2003年6月以來,由于掃描線矯正器(SCL)出故障導致Landsat 7影像數(shù)據(jù)出現(xiàn)空白區(qū)域,下載的影像缺失區(qū)域呈條帶狀，如圖3(a)所示,因此借助ENVI的Landsat_gapfill插件使用三角插值方法對2012年的影像數(shù)據(jù)進行填充修復,修復后的影像能夠?qū)θ笔У臈l帶空白地區(qū)進行填充,清晰展示各地物的地貌特征,如圖3(b)所示。

表4 2010—2019年影像信息Table 4 2010—2019 Image information

圖2 GEE中選取影像流程Fig.2 Select image flow chart in GEE

圖3 2012年影像修復前后對比Fig.3 Comparison before and after image restoration in 2012

(2)樣本數(shù)據(jù)處理。對上傳到GEE的樣本數(shù)據(jù)添加土地覆蓋標簽landcover屬性:水體為0,不透水面為1,植被為2,裸地為3。

2 研究方法

2.1 工作流程

工作流程如圖4所示。第1步,在GEE云平臺中獲取Landsat系列最佳時相數(shù)據(jù),并進行數(shù)據(jù)預處理。第2步,在Google Earth中制作訓練樣本以及驗證樣本。第3步,進行光譜特征構(gòu)建。第4步,使用分類回歸樹(classification and regression tree,CART)、隨機森林(random forest)、支持向量機(Support vector machine,SVM)3種分類器對樣本數(shù)據(jù)集進行訓練并對訓練后的結(jié)果進行精度驗證。最后,根據(jù)分類結(jié)果對開封市城區(qū)10年間水體分布變化情況進行研究。

圖4 工作流程Fig.4 Work flow chart

2.2 光譜特征構(gòu)建

在GEE云平臺計算經(jīng)過數(shù)據(jù)預處理后得到的研究區(qū)2010—2019年每景影像的歸一化水指數(shù)(normalized difference water index,NDWI)[11]、改進的歸一化水指數(shù)(modified normalized difference water index,MNDWI)[12]、歸一化植被指數(shù)(normalized vegetation index,NDVI)、歸一化建筑指數(shù)(normalized difference building index, NDBI)。將計算出的特征指數(shù)作為獨立的光譜波段添加進原始影像中作為分類特征。4種特征指數(shù)的計算方法如表5所示,RRed為傳感器紅光波段反射率；RGreen為綠光波段反射率；RNIR為近紅外波段反射率；RSWIR為短波紅外波段反射率。

表5 各指數(shù)計算公式Table 5 Calculation formula of each index

2.3 水體提取方法

利用機器學習分類的方法對開封市城區(qū)的水體進行提取,其中CART分類、SVM分類、RF分類在遙感領域具有廣泛的應用[13-15]。因此使用CART、SVM、RF 3種分類方法在GEE中以樣本數(shù)據(jù)的landcover屬性為標簽在圖像上疊加點以進行訓練, 使用與訓練相同的光譜特征對影像進行地物分類進而進行水體信息的定量提取,然后對3種分類器的提取水體的精度進行評價,選擇最優(yōu)的分類結(jié)果統(tǒng)計開封市水體面積及其他地物的面積，并對水體的提取結(jié)果結(jié)合空間分布進行分析。

2.3.1 CART分類

分類與回歸樹(CART)由Breiman等[16]提出。CART是應用廣泛的決策樹學習方法,由特征選擇、樹的生成和剪枝組成,既可以用來分類也可以用來回歸。CART分類時使用基尼系數(shù)(Gini)來選擇最優(yōu)的數(shù)據(jù)分割閾值劃分決策樹。假設樣本數(shù)據(jù)共有n個類別，Gini的計算公式為

(1)

式(1)中：ck表示在數(shù)據(jù)集D中第k類的個數(shù)。

基于GEE云平臺使用CART分類對樣本數(shù)據(jù)根據(jù)選擇的光譜特征進行訓練自動選擇合適的閾值構(gòu)建決策樹對地物進行分類,提取水體信息。

2.3.2 支持向量機(SVM)

支持向量機(SVM)由Cortes等[17]提出。支持向量機算法優(yōu)化的目標是尋找最大化的邊界[18-19]。邊界定義為分離超平面與其最近的訓練樣本之間的距離。在GEE中對SVM分類選用徑向基函數(shù)核(RBF kernel)[20]以及對應的其他參數(shù)使用樣本數(shù)據(jù)對影像進行分類,實現(xiàn)對開封市城區(qū)水體分布信息的定量提取,徑向基函數(shù)核定義為

(2)

2.3.3 隨機森林(RF)

RF算法由Breiman[21]提出，是1種組合了多棵決策樹的技術(shù)。因為該方法具有極好的準確率，能夠評估各個特征在分類問題上的重要性以及能提供快速、可靠的分類結(jié)果,所以在影像分類上具有廣泛的應用。在進行RF分類時參數(shù)的設置非常重要,幾個重要的參數(shù)有樹節(jié)點數(shù)n、最大葉節(jié)點數(shù)m、分割節(jié)點的特征數(shù)量v。在GEE中為了更好地定量提取水體信息,在進行RF分類時經(jīng)過多次實驗,將n設為50,將m設為無限制,將v設為0(輸入特征總數(shù)的平方根)。

3 結(jié)果與分析

3.1 數(shù)據(jù)處理效率分析

與傳統(tǒng)的處理遙感影像工具ENVI、ArcGIS相比,GEE集成多種遙感影像數(shù)據(jù),對遙感影像進行分析只需編寫少量代碼就能實現(xiàn),且GEE基本不占用本地資源。以2018年的影像處理為例,在GEE上加載2018年覆蓋研究區(qū)域的影像進行裁剪到可視化需要5 s左右的時間,從影像獲取到影像分類處理完畢并計算出結(jié)果的整個流程所需時間在15 s左右；而在本地環(huán)境(CPU:i5-4200H,內(nèi)存:8 GB,帶寬:100 MB)下,所用的2018年4月8日單景影像數(shù)據(jù)LC81230362018098LGN00,大小為914 MB,理論情況下,在百兆帶寬的網(wǎng)絡環(huán)境下下載時間為71.41 s左右,使用ENVI 5.3對影像進行裁剪、輻射校正、大氣校正花費的時間在96 s左右(排除掉手動操作的時間),再進行后續(xù)分類操作需要的時間遠比使用GEE處理的時間要長。使用GEE能夠極大地提高運行效率以及計算能力,效率改善達到10個量級以上。

3.2 精度評價與分析

根據(jù)研究區(qū)域范圍以及樣本數(shù)量結(jié)合算法本身特征將樣本數(shù)據(jù)按照7∶3的比例進行分配,從樣本數(shù)據(jù)中隨機選取70%的樣本數(shù)據(jù)作為訓練數(shù)據(jù)集進行分類的訓練,剩下30%的樣本數(shù)據(jù)作為驗證數(shù)據(jù)集用于驗證分類結(jié)果的精確度。由于用戶精度、總體精度以及Kappa系數(shù)能夠很好地反映地物的分類精度,所以使用基于混淆矩陣計算所得的水體提取精度、總體精度以及Kappa系數(shù)作為評價指標。

3種分類方法對開封市城區(qū)水體分類精度如圖5所示。使用CART分類器提取水體的平均精度為96.4%,最小精度為93.3%；使用RF分類器提取水體的平均精度為97.1%,最小精度為94.6%；使用SVM分類提取水體的平均精度為95.9%,最小精度為92.1%。3種分類方法對開封市地物分類的總體精度(OA)和Kappa系數(shù)如圖6、圖7所示,使用CART分類方法的平均總體精度為96.6%,平均Kappa系數(shù)為0.950；使用RF分類方法的平均總體精度為96.9%,平均Kappa系數(shù)為0.954；使用SVM分類方法的平均總體精度為95.1%,平均Kappa系數(shù)為0.930。經(jīng)過對比使用RF分類能夠更好地提取開封市城區(qū)的水體信息。

圖5 水體分類精度Fig.5 Classification accuracy of water body

圖6 3種分類器的總體分類精度Fig.6 Overall classification accuracy of three classifiers

圖7 3種分類器kappa系數(shù)Fig.7 Kappa coefficients of three classifiers

3.3 水體提取結(jié)果與分析

使用水體提取精度最好的RF分類方法分析開封市城區(qū)2010—2019年的水體分布變化情況。

如圖8所示,近10年研究區(qū)的景觀格局發(fā)生了較大變化,開封市城區(qū)的水體主要分布在北部以及中部,開封市西北高東南低的地勢造就了河流自西北向東南流的格局。通過表6的面積變化和圖9的折線走勢得出，開封市城區(qū)水體覆蓋率先減少后升高,總體上是增加的趨勢；植被與裸地總體上趨于減少；不透水面自2016年起較之前年份增長較多,總體上不透水面區(qū)域在不斷增多。

圖8 2010—2019研究區(qū)區(qū)各類土地分類結(jié)果Fig.8 Classification results of various types of land in the study area from 2010 to 2019

圖9 2010—2019年研究區(qū)各類土地覆蓋率Fig.9 Various land cover rates in the study area from 2010 to 2019

表6 2010—2019研究區(qū)各類土地面積Table 6 Area of various types of land in the study area from 2010 to 2019

2010—2011年水體和裸地面積在減少,植被和城鎮(zhèn)的面積在上升,北部的裸地部分被植被覆蓋,中部的裸地得到利用變成了不透水面,西北部的黃河流域由于水位下降部分轉(zhuǎn)變成裸地,水體面積從28.17 km2下降到23.75 km2；2011—2012年黃河水位恢復,中部和東部水體區(qū)域增多,總體水體覆蓋率基本恢復到2010年的水平,為26.65 km2;2012—2013年植被面積在增加,其他3類地物的面積均在減少,結(jié)合李肅秋等[22]的研究數(shù)據(jù)2013年平均降水量達到歷年的最低點335.3 mm，分析此年受降水量影響導致水體面積降至最低,為17.75 km2；2014年降水量增加水體面積有所恢復；2014—2016年,開封市政府為了更好地將水系連貫起來提升城市形象以及改善居民生活修建完成開封西湖,中部的部分裸地變成了水體區(qū)域,水體面積增加。2016—2018年,受黃河水體面積影響以及城市在往西擴建發(fā)展使得西部地區(qū)的水體區(qū)域在不斷減少。2018—2019年，城市在不斷發(fā)展不透水面的覆蓋率在提高，由44.93%上升到47.13%,城市的快速發(fā)展同時，環(huán)境以及人文生活都在發(fā)生改變。為了改善環(huán)境以及能更好地體現(xiàn)開封的“宋文化”,開封市通過“一渠六河”“十湖連通”等工程對城區(qū)的水系進行修整,水體面積增加從18.60 km2上升到31.16 km2。

總體上,隨著城市的發(fā)展,城市水體面積也在不斷發(fā)生變化,影響水體面積有2個主要原因,其一是由于每年的降水量不同導致研究區(qū)北部的黃河流域的水量每年在不斷變化,其二開封市古城保護和城市雙修建設,水體面積也在逐步增加。整體上,除掉黃河的影響,開封市的水體變化呈現(xiàn)好的發(fā)展態(tài)勢,這是鄭汴一體化期間開封市政府湖泊保護政策取得的顯著成果。

4 結(jié)論

研究了開封市城區(qū)的年際水體動態(tài)變化,探索GEE云平臺的年際水體變化監(jiān)測能力,基于GEE云平臺選取2010—2019年4—8月份的極少云(云量小于10%)Landsat系列影像數(shù)據(jù),使用CART、SVM以及隨機森林3種分類器提取開封市城區(qū)的水體分布信息,結(jié)論如下。

(1)GEE在云端存儲海量數(shù)據(jù),在進行科學研究時不需要在本地下載以及處理,只需將云端處理后的結(jié)果下載到本地即可,且其使用谷歌公司強大的計算能力,能夠?qū)﹂_封市城區(qū)的年際水體信息的提取提供有力保障。

(2)使用CART、RF、SVM 3種分類器都能夠?qū)﹂_封市城區(qū)年際水體信息進行有效提取,經(jīng)過對比使用RF分類提取的精度最高,其提取的水體平均精度為97.1%,提取水體的最小精度為94.6%,平均總體分類精度為96.9%,平均Kappa系數(shù)為0.954。

(3)開封市城區(qū)水體的面積及空間處于不斷變化中,水體面積總體趨勢在增長,水體區(qū)域主要分布在北部和中西部。開封市近年來大力推進“鄭汴一體化”的步伐,向西與鄭州相向發(fā)展,水體區(qū)域在不斷規(guī)劃完善。總體上導致水體面積不斷變化，主要因素既有北部黃河的水含量在不斷變化,也有開封城市化腳步的加快使城區(qū)的景觀格局發(fā)生改變。

由于使用的影像數(shù)據(jù)為中高分辨率影像數(shù)據(jù),不可避免地出現(xiàn)了混合像元的現(xiàn)象致地物分類受到影響,且影響細微水體區(qū)域的提取。下一步研究方向?qū)⑹褂酶叻直媛市l(wèi)星的影像數(shù)據(jù)進行水體的提取，將結(jié)合深度學習方法研究更加高效的水體提取算法,深入探索城市水體分布變化與城市生態(tài)環(huán)境變化之間的關(guān)系。