基于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的运用提供了更广阔的前景。