基于Google Earth Engine云計算的城市群生態質量長時序動態監測——以粵港澳大灣區為例

王 淵,趙宇豪,2,吳健生,2,*

1 北京大學城市規劃與設計學院, 城市人居環境科學與技術重點實驗室, 深圳 518055

2 北京大學城市與環境學院, 地表過程分析與模擬教育部重點實驗室, 北京 100871

生態環境是指影響人類生產生活和生態系統發展的各種生態因素的總和,與人類生存和發展所處的環境以及社會可持續發展息息相關[1]。地球進入人類世后[2],人類活動對全球生態環境的影響愈發強烈[3-5],其主要是通過改變地表狀況從而影響生態環境,而城市化是造成地表變化的重要原因之一[6-7]。中國是目前城市化速度最快的國家之一,1978 到2019年間,城市化率增長了 42.68%[8],然而傳統城市化過度側重于發展速度,導致生態環境被破壞。隨著新時代中國經濟和城鎮化開始由高速增長轉向高質量發展[9],國家日益重視生態環境評估與保護。其中生態環境部根據生物豐度指數,植被覆蓋指數,水網密度指數,土地退化指數和環境質量指數構建的生態指數(Ecological Index, EI),在區域生態環境質量(Ecological Environmental Quality, EEQ) 評估方面得到了廣泛應用[10]。在實際應用中,學者們雖然可以根據研究區域對指標和權重進行了不同的調整,但普遍面臨評價指標提取困難、數據空間精度較低和數據更新慢等問題[11]。而及時動態地監測生態環境狀況,明晰生態環境變化的特點和趨勢,對生態環境管理和社會可持續發展具有重大意義。

衛星遙感具有大面積、快速和周期性的重復觀測的優點,已被廣泛應用于生態學研究領域,使得生態環境質量評價工作得到了改善[12-13]。目前,已創建了多種遙感指數來量化生態狀況,如歸一化植被指數(NDVI)、增強植被指數(EVI)、永久植被分數(PVF)和干旱條件指數(SDCI)等,但大多都是以某個特定的生態相關主題為導向,不能全面的評估區域的綜合生態狀況[14]。而基于遙感信息并結合了綠度、濕度、干度、熱度的遙感生態指數(remote sensing based ecological index, RSEI)[15]可以比較好的解決上述問題。遙感生態指數指標容易獲取且計算簡便,無需人為設定權重和閾值,是一種客觀、快速、簡便的城市生態質量的監測和評價技術[16-17]。通過該指數來評估生態質量的時空分異,已在不同尺度得到了廣泛的應用[12,18-19]。然而在遙感生態指數的應用中還存在一些問題。首先,遙感影像普遍面臨云遮擋問題,去云較為困難,直接去云則會造成云遮擋區域數據缺失；其次,不同景影像的獲取時間有所差異,拼接有一定困難,可比性不足。為規避上述問題,部分研究通常選取云量較少的小塊區域某幾個時間點的數據開展研究[20-21],而云量較多區域的大范圍長時序研究則相對較少[22]。近年來,谷歌云計算的Google Earth Engine(GEE)平臺的快速發展為遙感數據提供了強大處理平臺[23],基于GEE的圖像處理可以較好的改善遙感數據缺失、多云、色差、時間不一致等問題。在大范圍長時間序列的遙感應用研究,GEE更具優勢,極大的縮短了圖像處理時間,提高了工作效率[24]。

粵港澳大灣區(Guangdong-Hong Kong-Macao Greater Bay Area, GBA)的發展歷程可以看做中國城市化的縮影,對其進行長時序的生態質量時空變化評估,不僅可以助力該區域的建設和綠色發展,對全國其他地區也有借鑒意義。同時,粵港澳大灣區是典型的南方多云地區,基于GEE平臺的影像數據處理和遙感生態指數計算可以為其他類似地區的提供參考。因此,本文借助GEE平臺,通過逐年遍歷1988—2018年粵港澳大灣區3530景Landsat遙感影像,對其去云、計算指標、疊加提取中值、鑲嵌等處理,提取年度遙感生態指數,評價了近三十年來該區域的生態質量的動態變化,分析了生態質量的時空異質性并比較城市間發展態勢的差異,以期在快速城市化背景下,為土地管理和生態保護提供理論依據。

1 研究區概況與數據來源

1.1 研究區概況

粵港澳大灣區(21°32′—24°26′N,111°20′—115°24′E)位于中國華南地區(圖1),屬亞熱帶季風氣候,年平均氣溫為22°C,雨季集中在4月至9月[25]。其由香港、澳門兩個特別行政區和廣州、深圳、珠海、佛山、惠州、東莞、中山、江門、肇慶九個珠三角城市組成,總面積5.6萬km2,2019年末常住7264.92萬人,GDP總量高達11.59萬億元人民幣,是中國開放程度最高、經濟活力最強的區域之一,在“一帶一路”建設和國家發展大局中具有重要地位[26]。

圖1 研究區Fig.1 Study area

1.2 數據來源與預處理

采用的數據及來源如表1所示。其中,Landsat數據來源于美國地質調查局(United States Geological Survey, USGS),在GEE平臺集成,空間分辨率為30m,主要包括原始圖像(Digital Number, DN),經輻射校正的大氣層頂表觀反射率圖像(Top of Atmosphere Reflectance, TOA Reflectance)和經輻射校正和大氣校正的地表反射率圖像(Surface Reflectance, SR),均已完成幾何精校正。由于Landsat7 衛星的ETM+機載掃描行校正器(ScanLinesCorrector, SLC)自2003年故障,導致數據條帶的部分丟失,為了盡量避免傳感器之間差異的影響,研究中使用1988—2011年Landsat5衛星的TM數據和2013—2018年Landsat8衛星的OLS和TIRS數據,不含2012年。

表1 數據集說明Table 1 Dataset descriptions

預處理包括去云和水體掩膜兩個部分。首先對Landsat影像利用CFMASK算法生成的質量評估波段QA進行去云處理,它通過標示哪個像素可能受儀器或云層影響,從而提高了科學研究的完整性(https://www.usgs.gov/land-resources/nli/landsat)。具體過程為：選擇出有云陰影覆蓋、有云并且云層置信度為中等的像元,將其像元值設置為0。其次,由于水體會影響主成分載荷,根據徐涵秋提出的改進的歸一化水體指數(Modified Normalized Difference Water Index, MNDWI)對水體掩膜[27]。

2 研究方法

遙感生態指數(RSEI)由歸一化植被指數(NDVI)、濕度分量(WET)、地表溫度(LST)和干度指數(NDBSI)構成,分別反映與人類生存息息相關的綠度、濕度、熱度和干度四種生態要素。其所選的指標完全基于遙感信息,容易獲得,且計算過程無需人工干預,因此結果客觀可靠、可比性強[15]。

2.1 城市綠度

城市綠度是指城市范圍內為植被覆蓋的區域,且具有一定的生態服務效益,對生態質量有積極的影響[28]。歸一化植被指數(Normalized Vegetation Index,NDVI)據植被葉面在紅光波段的吸收和近紅外波段的反射特性構建,能反映植物生物量、葉面積指數以及植被覆蓋度,因此,以NDVI 表征城市綠度指標,計算方法見公式(1)。

(1)

式中ρred、ρNIR分別表示TM 影像和OLI影像影像所對應的紅波段和近紅外波段的反射率。

2.2 地表濕度

基于纓帽變換的濕度分量(WET)可反映地表水體條件,特別是土壤的濕度狀態,基于 TM 和OLI數據進行 WET 提取[29-30],計算方法見公式(2)(3)。在提取之前,使用MNDWI對水體進行掩膜,綜合前人研究[31]和目視解譯,設置閾值為0.15,計算方法見公式(4),使WET反映真實的陸地地表濕度狀況。

WET(TM)=0.0315ρblue+0.2021ρgreen+0.3102ρred+0.1594ρNIR-0.6806ρSWIR1-0.6109ρSWIR2

(2)

WET(OLI)=0.1511ρblue+0.1972ρgreen+0.3283ρred+0.3407ρNIR-0.7117ρSWIR1-0.4559ρSWIR2

(3)

(4)

式中ρblue、ρgreen、ρred、ρNIR、ρSWIR1、ρSWIR2分別表示TM 影像和OLI影像所對應的藍波段、綠波段、紅波段、近紅外波段、短波紅外 1 波段和短波紅外 2 波段的地物反射率。

2.3 城市干度

建筑物是城市人工生態系統的重要組成部分,建筑不透水面取代原有自然生態系統,導致了地表的“干化”,因此利用建筑裸土指數代表“干度”。應用建筑指數(index-based built-up index, IBI)[16]和 裸土指數(soil index,SI)[32]合成干度指標,記為干度指數(normalized difference built-up and soil index, NDBSI),計算方式如下所示

NDBSI=SI+IBI/2

(5)

(6)

(7)

式中ρblue、ρgreen、ρred、ρNIR、ρSWIR1、ρSWIR2分別表示TM 影像和OLI影像所對應的藍波段、綠波段、紅波段、近紅外波段、短波紅外 1 波段和短波紅外 2 波段的地物反射率。

2.4 地表溫度

地表溫度(Land Surface Temperature, LST)是地球能量收支的重要組成部分。是反映地表環境的重要參數,本研究以經過反演的地表溫度表征熱度指標。由于Landsat8衛星的TIRS傳感器有兩個熱紅外波段,因為TIRS band 10波段相對于TIRS band 11波段而言位于一個較低的大氣吸收區域,具有更高精度的大氣透過率,因此選取Landsat5的第6波段和Landsat8的第10波段作為地表溫度反演通道。LST采用統計單窗模型(statistical mono-window model,SWM)進行反演[33]。采用植被覆蓋度計算地物比輻射率[34-35],植被覆蓋率由NDVI計算得到。SWM具體計算方法如下所示。

(8)

(9)

Fveg=NDVI-NDVIbare/NDVIveg-NDVIbare2

(10)

式中Tb表示TIRS通道的大氣頂層反射值,ε是地表比輻射率,系數Ai、Bi、Ci由線性回歸確定,是對10類TCWV(i=1,...,10)進行的輻射傳輸模擬[36]。NDVIbare為裸露土壤或建筑表面的NDVI值,NDVIveg為全植被覆蓋區的NDVI值,當NDVI大于0.7時,取植被覆蓋度為1。當NDVI小于0時,取植被覆蓋度為0。

2.5 綜合生態指數構建

對去云之后的遙感影像分別計算各遙感指標,通過影像疊加、提取中值(式11),能夠消除異常值和改善色差影響。

(11)

式中,n表示該像素點處某一年衛星訪問的次數,即包含的影像總數,pix表示第i張影像的像素值,pix_f表示取中值后最終的像素值。

最后,對單一指標采用極差標準化去除量綱,依據主成分變換(式12)構建遙感生態指數。為方便比較,將RSEI標準化處理至[0,1]之間,以0.2為間隔,將RESI劃分為差(0—0.2)、較差(0.2—0.4)、中等(0.4—0.6)、良(0.6—0.8)和優(0.8—1.0)五個等級[37]。

RSEI0=PCAfNDVI,WET,NDBSI,LST

(12)

(13)

3 結果分析

3.1 模型檢驗3.1.1 RSEI模型表現

將標準化后的各年度綠度、濕度、干度、熱度指標進行波段合成,PCA變換后得到主成分分析結果,等間隔抽取樣本進行檢驗。由表2可知,1988、1998、2008、2018四個年份的第一主成分貢獻率分別為72.40%、65.75%、77.72%和85.34%,表明主成分集中了大部分特征,且每個指標在第一主成分上的載荷分布均勻,因此基于第一主成分提取的信息來表征RSEI是合理的。

表2 主成分分析結果Table 2 The result of principal component analysis

3.1.2相關性檢驗

采用平均相關度模型檢驗 RSEI的適用性,相關系數越接近1,表明模型的綜合代表程度越高,適宜性越強[15]。計算方法見式(14)。

(14)

表3 各指標相關性矩陣Table 3 Correlation matrix among RSEI and 4 factors

計算各指標之間的Pearson相關系數,在99%置信水平P值均為0.000,均通過顯著性檢驗。平均相關度計算結果表明,RSEI的平均相關度最大,多年的平均值為0.761,其中1988、1998、2008和2018年分別為的0.744、0.681、0.780和0.840。其次是NDBSI,平均相關度為0.604,WET、LST、NDVI的相關度較低,分別為0.478、0.407、0.340,這表明綜合多個指標的RSEI比單一指標更適合用來評價生態質量。RSEI與各分指標的相關系數中,NDVI和WET與RSEI呈正相關,NDBSI和LST與RSEI呈負相關,表明RSEI值越大態質量越好。NDBSI的相關系數最大,相關系數的絕對值均大于0.9,其主要反應裸土、建筑分布情況,進表明城市擴張會造成生態質量的惡化。

3.2 生態質量時空變化

1988—2018年粵港澳大灣區RSEI的統計值分布如圖2所示,箱線圖中展示了RSEI的下四分位數、中位數、上四分位數、上限和下限,限制線以外的異常值不顯示。圖中折線表示RSEI的均值變化,多年RESI均值為0.632,表明其總體生態質量處于良好水平。1988—1991年間,生態質量呈現上升趨勢,1991—2000年生態質量逐漸下降,2000—2002年呈上升趨勢,到2003年RSEI達到近30年最大值0.802,達到優秀水平。此后直至2018年表現為波動下降的趨勢。

圖2 1988—2018年粵港澳大灣區遙感生態指數RSEI的統計值分布Fig.2 Statistical distribution of RSEI in the GBA from 1988 to 2018

為確保階段劃分的科學性,使用基于梯度下降法的最優迭代算法對RSEI值進行自動分段擬合。兩種分段擬合曲線和精度檢驗結果如圖3所示。圖3(a1)顯示通過直接觀察圖2初步劃分的斷點位置,即1991年、2000年和2003年。圖3(b1)顯示迭代后的自動擬合結果,斷點為1991年、2000年和2002年。圖3(a2)和(b2)比較了兩種方法的擬合精度,從分段擬合曲線的R2值來看,自動線性擬合的精度更高。由于2002至2018年時間跨度較大,新增2009年作為時間節點進行后續分析。

圖3 1988—2018年粵港澳大灣區的遙感生態指數RSEI分布Fig.3 Distribution of RSEI in the GBA from 1988 to 2018

RSEI值的空間分布如圖4所示,圖中生態質量較低的紅色區域在前期分布相對較為分散,與各城市建成區所處地理位置保持高度一致,呈現出明顯的多中心分布。而后期主要集中分布在佛山、中山、廣州南、東莞、深圳和惠州的市中心和惠陽區,表現為明顯的空間集聚分布特征。同時從圖4可以看出,研究區域內的生態質量空間異質性高,且城市區域的生態質量明顯低于非城市區域。

圖4 1988—2018年粵港澳大灣區生態質量變化Fig.4 Remote sensing ecological index (RSEI) changes of the GBA from 1988 to 2018

3.3 生態質量變化檢測

為了進一步分析生態質量變化的時空分布狀況將RSEI的 五個等級,采用差值法對多時相的生態質量進行變化檢測(表4)。

表4 生態質量等級轉移矩陣Table 4 Transition matrix of ecological quality levels

圖5和圖6表示變化檢測結果,紅色、黃色和綠色分別表示生態質量退化、不變和提高。1988—1991年,有43%的區域生態質量提高,生態質量下降的區域占6%,空間上分布較為零散。1991—2000年,生態質量提高的區域占10%,生態質量下降的區域達45%,分布在廣州、佛山、東莞、深圳的大部分區域。2000—2002年,有8%的區域生態質量退化,主要分布在肇慶市東南部、惠州市東北部以及深圳市南部沿海區域,區域生態質量提高的區域占43%,不同程度地分布在每個城市。2002—2018年生態質量總體下降,其中2002—2009年,生態質量下降區域占35%,2009—2018年,生態質量下降區域占36%,空間分布較為均勻,每個城市有不同程度的下降。從1988—2018年的變化檢測結果來看,三十年來共有36%的區域生態質量下降,廣州、東莞、佛山、珠海生態質量退化明顯,肇慶西北部、深圳、惠州等區域也有不同程度的退化。

圖5 1988—2018年生態質量變化檢測Fig.5 Change detection of ecological quality in the GBA from 1988 to 2018

圖6 1988—2018年粵港澳大灣區生態質量變化Fig.6 Changes of ecological quality in the GBA from 1988 to 2018

3.4 生態質量分區統計

統計各縣(區)尺度上三十年來生態質量變化檢測結果的均值,將其分為五個級別,分別為總體改善、基本不變、輕度退化、中度退化和重度退化。由圖7可知,生態質量存在明顯空間分異。其中佛山、中山、東莞、廣州南部以及珠海市金灣區、江門市江海區為重度退化,空間上集聚分布在粵港澳大灣區的中心,這些區域海拔較低,城市擴張范圍大,綠地保有量較少,城市發展過程中生態質量遭到破壞。澳門、深圳市的寶安區、龍華區、光明區、珠海市蓬江區、肇慶市端州區下降幅度次之,表現為中度退化。江門市大部分地區、佛山西部的高明區以及與其毗鄰的肇慶市高要區、廣州北部的從化區、惠州市博羅縣等生態質量整體改善,這些區域相對而言山地居多,以自然景觀為主,城市開發強度較弱,此外,深圳市南山區和福田區的生態質量也得以改善。研究區北部和香港生態質量整體基本不變,這可能是因為香港本身處于高度城市化水平,開發空間有限,且香港對生態用地保護嚴格,土地開發阻力和難度較大,而研究區北部多山區,基本沒有開發利用。

圖7 1988—2018年不同區域生態質量平均變化結果Fig.7 Results of average variation of ecological quality in different regions from 1988 to 2018

4 討論4.1 研究亮點

以往在運用RSEI生態指數對區域生態質量時空變化進行監測時,經常會受到遙感圖像質量的制約,為消除圖像質量的影響,研究多采取在年度影像中挑選質量較好且季節相近的影像數據[20,38],或者僅進行小范圍的研究[39-40],而對于云量較多的區域往往會造成數據缺失,通常采用鄰近年份代替的方法[21],這限制了RSEI指數的準確性、可比性以及研究的區域范圍和時間序列長度。而本研究在計算粵港澳大灣區1988—2018三十年來的RSEI時,基于Google Earth Engine(GEE)提供的強大云計算能力可以較好的避免上述問題。RSEI指數中的綠度、濕度、干度、熱度等指標均是基于全年遙感數據計算的結果,通過疊加全年所有圖像,提取中值的方法,改善了去云、圖像色差等常見問題。同時,該處理方法同樣可以作為季度、月度等周期性時間段的數據處理的參考,增加了研究結果在時間序列上的可比性。GEE的使用可以極大的提高影像處理的效率、在大范圍長時間序列的遙感應用研究中極具優勢,為RSEI的運用提供了更廣闊的前景。

4.2 模型指標關系

生態質量與四個分指標密切相關。根據相關性分析結果,NDVI和WET對RSEI起正向作用,即在指標的數量區間內,隨著NDVI和WET增大,RSEI相應增加,生態質量越好。相反,隨著NDBSI和LST的增加,RSEI隨之減小,生態質量下降。從標準化指標的系數來看,其中NDBSI系數的絕對值最大,對RSEI的影響最大。在一定程度上證實了城市擴張帶來的城市景觀的變換會對生態環境產生負面影響,因此城市規劃和設計者應該兼顧生態環境,權衡城市發展和生態環境質量的關系,促進可持續的城市發展。需要指出的是,植被質量對生態質量的影響很大,因此NDVI合成通常根據植被物候期提取特定季節的數據。而粵港澳大灣區處于亞熱帶季風氣候區,終年溫暖濕潤,植被物候特征不明顯[41],同時,其位于南方多云地區,遙感成像易受天氣條件影響,數據缺失問題更加突出,僅選擇生長季數據,會導致研究區部分年份出現較多空值。因此本文在NDVI指標提取過程中對全年所有影像數據進行合成。相反,在北方地區,植被物候明顯,且云量較少,可以根據研究區的植被生長季節的起訖日期,進而確定相應的NDVI數值[42]。

4.3 不足與展望

已有的評價指標體系仍有待完善之處。首先,四項指標是否可以全面表征區域生態質量值得商榷,因此下一步可以考慮增加指標數量[43]以提升指標對區域生態質量的代表性。其次,四項指標本身的評價方式也可以進一步提升。如溫度和濕度指標對生態質量的影響應該具有條件限制,地表溫度對人類健康風險的影響存在閾值,過高或過低均會導致死亡率增加[44],可以考慮根據研究區狀況制定合適的閾值之后再進行歸一化[45]。此外,城市生態質量的綠度、濕度、干度、熱度等指標是從二維結構方面的考慮,可以考慮引入城市植被三維信息等空間度量進行優化。

在數據處理上,當研究區域過大、像素點過多時,基于GEE的大范圍像素級的處理可能會遇到內存溢出、處理速度慢等問題,可以采用圖像分割并行計算等方法進行改進。文中主要側重于使用GEE平臺實現長時間大范圍的生態質量監測,僅對1988—2018年來大灣區生態質量變化趨勢和原因進行了簡單論述,但并沒有對驅動力做具體的分析和討論,這是下一步將要開展的工作。

5 結論

本文基于Google Earth Engine平臺,對1988—2018年來粵港澳大灣區共3530景Landsat遙感影像進行批量去云處理,逐年提取綠度、濕度、干度和熱度等遙感指標,構建生態指數RSEI,并基于梯度下降法的最優迭代算法確定生態質量變化拐點,評價了粵港澳大灣區近30年間四個不同時段內的區域生態質量時空變化。主要結論如下：

(1)從時間上看,粵港澳大灣區三十年來生態環境質量呈現先上升、再下降、后上升、最后下降的總體趨勢。空間上,生態質量具有明顯的空間異質性,主要呈現出西北部和東北部的生態質量高,中部城市建成區生態質量低的狀態,粵港澳大灣區應加強區域協調發展,形成優勢互補的可持續發展模式。

(2)1988—2018年期間佛山東南部、中山、東莞、廣州南部以及珠海市金灣區、江門市江海區的生態質量表現為重度退化, 澳門、深圳市的寶安區、龍華區、光明區、珠海市蓬江區和肇慶市端州區表現為中度退化,肇慶市西北部、深圳市東部、珠海市香洲區、江門市新會區等區域表現為輕度退化,江門市大部分地區、佛山西部的高明區以及肇慶市高要區、廣州北部的從化區、惠州市博羅縣等生態質量整體改善,香港生態質量基本不變。

(3)1988—1991年,生態質量上升,生態質量提高和下降區域分別占43%和6%,1991—2000年,生態質量下降,生態質量提高和下降區域分別占10%和45%,下降區域分布在廣州、佛山、東莞、深圳的大部分區域,2000—2002年生態質量上升,生態質量提高和下降區域分別占43%和8%,下區域主要分布在肇慶市東南部、惠州市東北部以及深圳市南部沿海區域,2002—2018年生態質量總體下降,其中2002—2009年,生態質量下降區域占35%,2002—2018年,生態質量下降區域占36%。

(4)Google Earth Engine平臺可以較好的改善遙感數據缺失、多云、色差、時間不一致等問題,極大的提高影像處理的效率、在大范圍長時間序列的遙感應用研究中極具優勢,為RSEI的運用提供了更廣闊的前景。