基于GEE的桂林市主城區(qū)熱環(huán)境變化定量遙感分析

婁佩卿,付波霖,何宏昌,高二濤,范冬林,唐廷元,林星辰,閉 璐

(桂林理工大學(xué) 測(cè)繪地理信息學(xué)院，廣西 桂林 541006)

0 引 言

城市化進(jìn)程所伴隨的土地利用類(lèi)型變化對(duì)城市熱環(huán)境動(dòng)態(tài)變化的影響成為城市生態(tài)環(huán)境保護(hù)的重要問(wèn)題之一。城市熱環(huán)境是城市生態(tài)環(huán)境的一個(gè)重要參數(shù),近年來(lái)城鎮(zhèn)化進(jìn)程不斷改變下墊面形態(tài)及土地覆蓋類(lèi)型,對(duì)城市熱環(huán)境產(chǎn)生了深刻影響[1-3]。衛(wèi)星遙感可在短時(shí)間內(nèi)提供對(duì)地觀測(cè)數(shù)據(jù),是獲取地表溫度的重要手段[4-7]。采用衛(wèi)星遙感對(duì)城市熱環(huán)境及其影響因素進(jìn)行動(dòng)態(tài)分析具有重要意義。

當(dāng)前城市化進(jìn)程導(dǎo)致的土地利用變化帶來(lái)一系列生態(tài)環(huán)境問(wèn)題,如植被覆蓋面積大幅縮減[8]、河流湖泊大量消失[9]、熱島效應(yīng)持續(xù)增強(qiáng)[10]等。針對(duì)這些問(wèn)題,近年來(lái)國(guó)內(nèi)外學(xué)者開(kāi)展了關(guān)于土地利用變化對(duì)城區(qū)熱環(huán)境變化影響的研究。MODIS及AVHRR地表溫度產(chǎn)品多被用于地表溫度反演,但其空間分辨率較低,僅適用于大區(qū)域研究[11-15]。2013年發(fā)射成功的Landsat 8衛(wèi)星搭載熱紅外傳感器(TIRS),提供大量較高空間分辨率數(shù)據(jù),更適宜于城區(qū)范圍的地表溫度反演研究?；贚andsat 數(shù)據(jù)的地表溫度反演算法很多,如輻射傳輸方程法、單窗算法、劈窗算法及多通道多角度算法等。Jiménez-Muoz等基于Landsat 數(shù)據(jù)對(duì)單窗算法和劈窗算法的反演精度進(jìn)行對(duì)比,發(fā)現(xiàn)隨著大氣水汽含量增加,劈窗算法精度要略高于單通道算法[16];Yu等對(duì)輻射傳輸方程法、劈窗算法及單通道算法進(jìn)行對(duì)比分析,結(jié)果表明輻射傳輸方程法精度最高[17]。美國(guó)地質(zhì)調(diào)查局(USGS)指出TIRS 11波段存在定標(biāo)不穩(wěn)定性,不建議采用劈窗算法反演地表溫度等[18]。雖然王修信等[19]和梁保平等[20]已對(duì)桂林市地表溫度進(jìn)行了反演研究,但未對(duì)土地利用變化與植被覆蓋度及地表溫度之間的影響機(jī)制進(jìn)行分析,且亟需一套可實(shí)時(shí)快速進(jìn)行地表溫度反演的研究技術(shù)。

谷歌地球引擎(GEE)是用于衛(wèi)星影像及其他空間數(shù)據(jù)解算的開(kāi)源智能云平臺(tái)[21]。GEE提供全球尺度的Landsat TM/OLI、Sentinel-1/2、MODIS和DMSP/OLS燈光數(shù)據(jù)等多空間尺度、多源的遙感數(shù)據(jù),容量達(dá)到PB級(jí)別,且每天都在增加。GEE云平臺(tái)打破了傳統(tǒng)遙感軟件下載數(shù)據(jù)、預(yù)處理、信息提取、分析應(yīng)用從而獲取專(zhuān)題信息的定式[22]。通過(guò)GEE進(jìn)行JavaScript及Python語(yǔ)言編程,可不受時(shí)間和空間限制,快速、批量處理數(shù)據(jù),為遙感工作者提供極大便利。

本文基于GEE采用多時(shí)相Landsat遙感數(shù)據(jù)對(duì)桂林市主城區(qū)熱環(huán)境進(jìn)行了動(dòng)態(tài)分析,并探索土地利用類(lèi)型、植被覆蓋度及地表溫度等主要影響因素之間的關(guān)系,以此作為城市生態(tài)環(huán)境優(yōu)化的重要參考。研究結(jié)果對(duì)桂林市城鎮(zhèn)化建設(shè)及綠地建設(shè)提供技術(shù)支撐,并為減緩城市熱島效應(yīng),建設(shè)宜居城市提出合理建議。

1 數(shù)據(jù)源和處理方法

1.1 研究區(qū)域概況

桂林市位于廣西壯族自治區(qū)東北部,東經(jīng)109°36′50″—111°29′30″,北緯24°15′23″—26°23′30″,市中心海拔154 m,具有典型的喀斯特(巖溶)地貌, 是世界著名的旅游城市, 東盟自由貿(mào)易區(qū)門(mén)戶(hù)城市, 還是“一帶一路”和黔粵湘桂交界的重要連接點(diǎn)[23]。 本文以桂林市主城區(qū)(秀峰區(qū)、 疊彩區(qū)、 象山區(qū)、 七星區(qū))為研究區(qū), 總面積約為286 km2(圖1a)。 主城區(qū)地處低緯, 屬亞熱帶季風(fēng)氣候, 年平均氣溫近19.1 ℃, 7、 8月最熱, 平均氣溫約為28 ℃, 1、 2月最冷, 平均氣溫約為9 ℃。

圖1 研究區(qū)位置(a)與2018年11月17日Landsat 8全色影像(b)Fig.1 Location of the study area(a) and the full-color image of Landsat 8(b) on November 17, 2018

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

(1)Landsat影像:由GEE平臺(tái)提供的2010年Landsat TM影像、2014與2018年Landsat OLI影像?；贕EE API編程,篩選成像區(qū)間為當(dāng)年11月1日 —12月1日且云量小于8%的Landsat影像,最終選取的Landsat影像成像時(shí)間分別為2010年11月14日、2014年11月14日和2018年11月17日。

(2)其他數(shù)據(jù):SRTM DEM數(shù)據(jù)空間分辨率為30 m,由NASA提供,在GEE中可以直接調(diào)用;MODIS L1B Calibrated Reflectance產(chǎn)品同樣可以在GEE中直接調(diào)用,其主要用于估算研究區(qū)內(nèi)大氣水汽含量,成像時(shí)間分別為2010年11月14日、2014年11月14日和2018年11月17日,與Landsat影像保持對(duì)應(yīng)一致;訓(xùn)練樣本數(shù)據(jù)及精度驗(yàn)證數(shù)據(jù)自Google Earth Pro中目視解譯選取;地表溫度數(shù)據(jù)源于中國(guó)氣象數(shù)據(jù)網(wǎng)(http://data.cma.cn/)提供的國(guó)家級(jí)氣象站點(diǎn)地面氣象資料日值數(shù)據(jù)集,獲取時(shí)間與Landsat 一致,主要用于單窗算法的地表溫度反演。

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

基于GEE平臺(tái)的API編程篩選出已進(jìn)行過(guò)輻射定標(biāo)及大氣校正的Landsat大氣表觀反射率(TOA)影像, 調(diào)用GEE平臺(tái)中的SRTM DEM數(shù)據(jù)并計(jì)算歸一化植被指數(shù)(NDVI)及歸一化水體指數(shù)(NDWI)作為進(jìn)行隨機(jī)森林分類(lèi)時(shí)的特征變量以提高分類(lèi)精度。

于Google Earth Pro中對(duì)2010、2014及2018年研究區(qū)高分辨率影像進(jìn)行目視解譯,進(jìn)而選取樣本。將桂林市全城土地利用類(lèi)型劃分為建設(shè)用地、植被、未利用土地及水體4類(lèi),按照各種土地利用類(lèi)型在研究區(qū)內(nèi)的面積占比選取適量樣本點(diǎn)(表1)。

表1 土地利用類(lèi)型及樣本數(shù)Table 1 Types of land use and numbers

將樣本點(diǎn)以*.kml格式從Google Earth Pro中導(dǎo)出, 在ArcMap 10.2中轉(zhuǎn)化為*.shp格式, 導(dǎo)入GEE API中進(jìn)行調(diào)用, 將其中70%作為訓(xùn)練樣本, 30%作為驗(yàn)證樣本。

2 研究方法

2.1 基于GEE的隨機(jī)森林土地利用分類(lèi)

隨機(jī)森林算法是由多棵決策樹(shù)構(gòu)成的一種靈活且使用方便的機(jī)器學(xué)習(xí)算法, 即使未進(jìn)行參數(shù)調(diào)優(yōu), 也可獲得較好的分類(lèi)精度[23]。 其基本原理為: ①采用Bootstrap方法自原始數(shù)據(jù)有放回的抽取N組訓(xùn)練集, 抽取的每組訓(xùn)練集約為整體的2/3; ②利用這N組訓(xùn)練集構(gòu)建決策樹(shù),進(jìn)而產(chǎn)生由N棵決策樹(shù)構(gòu)成的森林,每棵決策樹(shù)生長(zhǎng)過(guò)程中均會(huì)從全部的M個(gè)特征變量中選取m個(gè)(m

在GEE中調(diào)用隨機(jī)森林分類(lèi)器,以土地利用類(lèi)型為目標(biāo)變量,并以Landsat影像的各波段加上SRTM DEM數(shù)據(jù)和NDVI及NDWI作為訓(xùn)練過(guò)程中輸入的特征變量,將決策樹(shù)的數(shù)量(ntree)設(shè)置為500,變量默認(rèn)數(shù)量(mtry)為特征變量總數(shù)的平方根[24],進(jìn)而對(duì)桂林市主城區(qū)進(jìn)行土地利用分類(lèi)。

2.2 基于GEE的地表溫度反演

2.2.1 植被覆蓋度及地表比輻射率 在GEE中計(jì)算得到3個(gè)時(shí)期研究區(qū)內(nèi)的NDVI,再根據(jù)NDVI推算出研究區(qū)內(nèi)植被覆蓋度PV,最后采用由Sobrino提出的NDVI閾值法[25-26]計(jì)算地表比輻射率ε

(1)

式中:NIR、R為的近紅外、紅光波段的反射值；NDVISoil為裸土或無(wú)植被覆蓋區(qū)域的NDVI值;NDVIVeg為植被完全覆蓋的像元的NDVI值, 即純植被像元的NDVI值。 選取經(jīng)驗(yàn)值即NDVISoil=0.05、NDVIVeg=0.70(當(dāng)不足0.70時(shí)即為最大像元值)。 當(dāng)某個(gè)像元NDVI>0.70時(shí),PV取值為1; 而當(dāng)NDVI<0.05時(shí),PV取值為0。

2.2.2 地表溫度 地表溫度指大氣與地表交界處的溫度,除入射太陽(yáng)輻射強(qiáng)度外,還與地物的發(fā)射輻射和反射輻射有關(guān),但這三者都受植被覆蓋度的影響,因而都與植被覆蓋度密切相關(guān)。本文選用精度較高的由傳統(tǒng)輻射方程改進(jìn)的單窗算法(mono-window algorithm)于GEE 中進(jìn)行對(duì)桂林市主城區(qū)地表溫度(land surface temperature,LST)進(jìn)行反演[27-28],USGS指出TIRS 11波段存在定標(biāo)不穩(wěn)定性[18],故本文將原始算法中所用的Landsat TM/ETM+第6波段改進(jìn)為L(zhǎng)andsat 8影像的第10波段,并對(duì)參數(shù)a、b進(jìn)行修正[29]，則實(shí)際地表溫度T為

T={a(1-C-D)+[b(1-C-D)+C+D]T10-DTa}/C,

(2)

式中:a=-62.735 657,b=0.434 036[28];C=ετ, 其中,τ為大氣透射率;D=(1-τ)[1+(1-τ)τ];T10為L(zhǎng)andsat 8的熱紅外波段10的像元亮度溫度, K。冬季(11月—次年1月)桂林市主城區(qū)所在位置(中緯度)的大氣平均作用溫度(Ta)與接近地面處(2 m左右)氣溫(T0)存在線性關(guān)系[29]

Ta=19.270 4+0.911 18T0,

(4)

T10=K2/ln(1+K1/Lλ)。

(5)

式中:K1、K2均可從Landsat 8數(shù)據(jù)的MTL頭文件中獲取(K1=774.89,K2=1 321.08);Lλ為波段10的像元亮溫。

綜上, 采用單窗算法進(jìn)行地表溫度反演需要大氣平均作用溫度、大氣透過(guò)率及地表比輻射率3個(gè)參數(shù)。大氣平均作用溫度可通過(guò)式(4)求得，大氣透過(guò)率可通過(guò)大氣水汽含量計(jì)算,而Landsat 8數(shù)據(jù)則可在NASA公布的網(wǎng)站(http://atmcorr.gsfc.nasa.gov)查詢(xún)，地表輻射率可通過(guò)式(1)求得。

3 結(jié)果與分析

3.1 土地利用動(dòng)態(tài)變化特征

為初步驗(yàn)證在GEE平臺(tái)中基于隨機(jī)森林算法的分類(lèi)結(jié)果,導(dǎo)出誤差矩陣精度評(píng)價(jià)結(jié)果,精度評(píng)價(jià)標(biāo)準(zhǔn)分為總體精度與Kappa系數(shù)。3個(gè)時(shí)期的總體精度分別為86.27%、90.78%及91.45%,Kappa系數(shù)分別為0.79、0.88及0.89。為進(jìn)一步驗(yàn)證分類(lèi)結(jié)果準(zhǔn)確性,在Google Earth Pro中的高分辨率影像選取驗(yàn)證點(diǎn)(每個(gè)時(shí)期選擇60個(gè)驗(yàn)證點(diǎn),每類(lèi)別15個(gè)),目視解譯判斷驗(yàn)證點(diǎn)地物類(lèi)別,導(dǎo)入GEE 中開(kāi)展分類(lèi)精度驗(yàn)證及空間一致性分析。結(jié)果表明，2010、2014、2018年的驗(yàn)證精度分別為83.3%、81.6%和88.3%，分類(lèi)精度均高于80%,滿(mǎn)足本文需求,混淆矩陣見(jiàn)表2。

表2 2010—2018年桂林市主城區(qū)土地利用精度 驗(yàn)證混淆矩陣Table 2 Land use accuracy verification confusion matrix of the main urban districts in Guilin from 2010 to 2018

在對(duì)桂林市主城區(qū)3個(gè)時(shí)期的Landsat影像進(jìn)行土地利用分類(lèi)時(shí)發(fā)現(xiàn)，桂林市某機(jī)場(chǎng)的地類(lèi)屬性并不一致,主要有兩個(gè)原因:①未利用土地與建設(shè)用地易發(fā)生混淆;②該機(jī)場(chǎng)在2010—2018年期間不斷建設(shè),水泥地面必然逐漸發(fā)生變化。

由圖2、 3可知, 桂林市主城區(qū)2010年植被覆蓋面積為184.76 km2, 2018年為129.47 km2, 較2010年縮減55.29 km2, 面積占比變化為-19.74%, 出現(xiàn)負(fù)增長(zhǎng); 2010年建設(shè)用地面積為50.98 km2, 2018年為135.49 km2, 較2010年增長(zhǎng)84.51 km2, 面積占比變化為30.19%, 面積正增長(zhǎng); 2010年未利用土地面積為21.23 km2, 2018年為6.38 km2, 較2010年減少14.85 km2, 面積占比變化為-4.49%, 呈現(xiàn)負(fù)增長(zhǎng); 2010年水體面積為23.03 km2, 2018年為6.38 km2, 較2010年縮減16.65 km2, 面積占比變化為-5.94%, 出現(xiàn)負(fù)增長(zhǎng)。 隨著城市發(fā)展, 人類(lèi)伐林用以造房及造路(建設(shè)用地及未利用土地面積增長(zhǎng)迅速), 占用河流、 池塘進(jìn)行圍湖造田(水體面積大大縮減)等活動(dòng)嚴(yán)重破壞了桂林市主城區(qū)的生態(tài)環(huán)境。

圖2 桂林市主城區(qū)3個(gè)時(shí)期土地利用類(lèi)型Fig.2 Land use types of three periods in the main districts of Guilin

3.2 植被覆蓋動(dòng)態(tài)變化特征

將2010、 2014及2018年的植被覆蓋采用密度分割劃分為高覆蓋度、 較高覆蓋度、 中等覆蓋度、 較低覆蓋度及低覆蓋度5個(gè)等級(jí)。 由圖4、 5可知, 2010年桂林市主城區(qū)以高植被覆蓋度為主, 面積為140.64 km2; 較高覆蓋度45.34 km2、 中等覆蓋度39.52 km2及較低覆蓋度35.64 km2相差不足5 km2; 而所占面積最少的為低覆蓋度, 面積18.86 km2。

圖3 土地利用類(lèi)型占地面積(a)和面積占比變化(b)Fig.3 Coverage by land use type(a) and land change in area percentage(b)

圖4 桂林市主城區(qū)3個(gè)時(shí)期植被覆蓋度分級(jí)Fig.4 Vegetation coverage classification in three periods in the main districts of Guilin

2018年桂林市主城區(qū)以中等植被覆蓋度69.47 km2及較高覆蓋度69.38 km2為主, 相比2010年中等植被覆蓋度面積擴(kuò)張10.52%, 較高植被覆蓋度面積擴(kuò)張8.48%; 較低覆蓋度58.34 km2及高覆蓋度56.97 km2相差不足5 km2, 相比2010年較低覆蓋度面積擴(kuò)張8.10%, 高覆蓋度縮減29.89%; 而所占面積最少的為低覆蓋度, 面積31.80 km2, 相比2010年擴(kuò)張4.59%。 隨著城市化進(jìn)程的發(fā)展,高植被覆蓋度面積大幅縮減, 其余4個(gè)植被覆蓋度分級(jí)均有所增長(zhǎng), 植被覆蓋狀況持續(xù)惡化。

3.3 熱環(huán)境動(dòng)態(tài)變化特征

將2010、2014及2018年的地表溫度采用密度分割劃分為高溫區(qū)、較高溫區(qū)、中等溫區(qū)、較低溫區(qū)及低溫區(qū)5個(gè)等級(jí),其中低溫區(qū)及較低溫區(qū)主要分布于植被覆蓋度較高的地區(qū)及水體周?chē)?而中等溫區(qū)、較高溫區(qū)及高溫區(qū)主要分布于建筑區(qū)及未利用土地等植被覆蓋度低的地區(qū)。由圖6、7可知,2010年桂林市主城區(qū)以較低溫區(qū)為主,面積為96.48 km2,其次為中等溫區(qū)(80.39 km2)、 低溫區(qū)(48.74 km2)及較高溫區(qū)(45.12 km2), 而所占面積最少的為高溫區(qū), 面積10.27 km2。 2018年桂林市主城區(qū)以較低溫區(qū)為主, 面積為86.20 km2, 其次為中等溫區(qū)(72.84 km2)、 較高溫區(qū)(66.70 km2)及低溫區(qū)(35.92 km2), 而所占面積最少的為高溫區(qū), 面積24.37 km2, 相比2010年，較低溫區(qū)(縮減3.67%)、 中等溫區(qū)(縮減2.70%)及低溫區(qū)(縮減4.58%)面積縮減, 而較高溫區(qū)(擴(kuò)張7.71%)及高溫區(qū)(擴(kuò)張5.04%)面積擴(kuò)張。 隨著城市化進(jìn)程的發(fā)展及植被覆蓋度狀況的嚴(yán)重惡化, 桂林市主城區(qū)地表溫度變化顯著, 平均溫度由2010年的29.24 ℃及2014年的29.67 ℃升高為2018年的30.53 ℃。

圖5 各級(jí)植被覆蓋度占地面積和面積變化率Fig.5 Coverage by land use type and cand ratio change of different vegetation coverage at different levels

圖6 桂林市主城區(qū)3個(gè)時(shí)期地表溫度分級(jí)Fig.6 Surface temperature classification in the main districts of Guilin in 2010, 2014 and 2018

4 結(jié) 論

本文采用基于GEE云平臺(tái)的單窗算法計(jì)算得到桂林市主城區(qū)的植被覆蓋度及地表溫度,并與基于GEE云平臺(tái)的隨機(jī)森林算法土地利用分類(lèi)結(jié)果進(jìn)行對(duì)比分析,得出如下結(jié)論:

圖7 各級(jí)溫區(qū)占地面積和面積變化率Fig.7 Coverage by land use type and district ratio change of different temperature zones at different levels

(1)2010—2018年桂林市主城區(qū)平均溫度呈上升趨勢(shì),增加1.29 ℃,且各級(jí)別溫區(qū)由低溫區(qū)、較低溫區(qū)及中等溫區(qū)轉(zhuǎn)化為較高溫區(qū)及高溫區(qū)。這與桂林市主城區(qū)的城鎮(zhèn)化進(jìn)程導(dǎo)致的植被覆蓋度下降有直接關(guān)系。

(2)較低溫區(qū)及低溫區(qū)主要分布于植被及水體覆蓋區(qū)域,而中等溫區(qū)、較高溫區(qū)及高溫區(qū)主要分布于建設(shè)用地及未利用土地覆蓋區(qū)域,故增加植被覆蓋度,涵養(yǎng)水源及減少未利用土地面積可有效緩解桂林市主城區(qū)熱環(huán)境的惡化。

(3)2014—2018年,高植被覆蓋度面積大幅縮減與城鎮(zhèn)化進(jìn)程進(jìn)一步加快有直接關(guān)系,其中高植被覆蓋度面積占比縮減31.34%,建設(shè)用地面積增長(zhǎng)30.19%。主城區(qū)未利用土地面積縮減4.49%進(jìn)一步說(shuō)明了城鎮(zhèn)化建設(shè)的加快。應(yīng)合理利用土地,減少土地浪費(fèi),通過(guò)植樹(shù)造林增加植被覆蓋度。

(4)基于GEE云平臺(tái)的隨機(jī)森林算法不僅具有運(yùn)算快速的有點(diǎn)且具有較好的分類(lèi)精度,2010、2014及2018年的驗(yàn)證精度均高于80%,是一種行之有效的土地利用分類(lèi)技術(shù)。

綜上所述,為改變桂林市主城區(qū)逐漸惡化的熱環(huán)境,應(yīng)制定更為合理的土地利用政策,優(yōu)化建設(shè)用地空間利用情況,減少未利用土地面積,通過(guò)植樹(shù)造林增加植被覆蓋度,進(jìn)而影響地表徑流來(lái)涵養(yǎng)水源,改善桂林市熱環(huán)境逐漸惡化的問(wèn)題。