2000—2020年广州城市绿地生态系统时空变化及其对地表蒸散的影响

2023-04-18刘昊东杨俏敏臧传富

热带地理 2023年3期

刘昊东，杨俏敏，臧传富

（华南师范大学地理科学学院，广州 510631）

蒸散是地表水和生态水吸热相变向大气扩散的水通量，作为地表与大气之间的联系机制（Mueller et al., 2013），是水循环中最重要的通量变量之一。每年约有60%的陆地降水以蒸散的形式进入大气（Jung et al., 2010），参与全球水循环；而在干旱气候下，该值可能超过90%（Oki et al., 2006）。此外，蒸散也是调节全球天气和气候的关键因素。大范围的地表蒸散变化会对降水、地表温度、大气环流产生显著影响（Shukla et al., 1982）；并且与干旱、热浪（Seneviratne et al., 2006; Vautard et al.,2007）等极端水文气象事件都有较为密切的联系（朴世龙等，2019）；区域地表蒸散的研究对水资源的优化配置和管理具有一定的参考价值（Liu et al.,2019）。

土地利用和土地覆盖是影响蒸散和其他水文过程的重要因素（Mlba et al., 2019）。在大范围尺度下，王麒麟等基于1982—2016年连续遥感观测，认为区域地表蒸散与该地森林覆盖率的变化之间存在正相关关系（Wang et al., 2021）；Teuling等（2019）在对欧洲平均蒸散发和径流量的影响因素的研究中指出，欧洲大部分区域蒸散量的增加是地区规模造林、森林覆盖率提高和气候变化的综合结果。随着城市化进程的推进，城市绿地的时空变化成为城市尺度下研究地区土地利用和土地覆盖的关键问题。城市绿地包括人工植被、半自然以及自然植被（臧传富等，2020），是城市整体的重要有机组成部分，也是支持生物多样性的关键栖息地（Lepczyk et al.,2017），同时为人类的健康生活提供了不可或缺的生态服务（Semeraro et al., 2021）。此外，城市绿地蒸散作为地区水文过程和生态过程相互作用机理的重要环节（臧传富等，2020），对区域各水文循环过程的再分配和再调节起一定促进作用。在高强度人类活动的介入下，城市土地覆盖、城市绿地在一定时期的类型转化过程更加多样且复杂，其对区域内地表蒸散等生态水文过程的总体影响具有更高的不确定性。

广州是珠三角城市群的核心城市。2000—2020年，广州制定了“南拓、北优、东进、西联”的空间发展战略，相继将原白云、番禺、从化3个县级市纳入市区范围，进入多核心快速发展的城市化阶段，城市土地利用和土地覆盖发生剧烈变化，城市化进程较为典型，该过程对区域地表蒸散的影响也更具研究意义。因此，本研究基于Landsat、MODIS和气象站监测数据对广州城市绿地时空变化及地表蒸散状况进行分析，探讨该时期广州复杂地表覆盖转化过程对区域蒸散的综合影响。以期为城市化进程下城市绿地生态建设和海绵城市布局提供科学依据。

1 研究区概况、数据与方法

1.1 研究区概况

按照现行行政区划广州共有荔湾、越秀、海珠、天河、白云、黄埔、番禺、花都、南沙、从化、增城11个市辖区，面积约7 434.4 km2（图1）。2020 年广州常住人口约1 874 万，地区生产总值约为25 019 亿元人民币（广州市统计局等，2021），是广东的政治中心和经济中心。改革开放近40 年来，广州成为城市化进程推进的热点区域，人地矛盾较为凸出，人类活动对地表景观的改造较为剧烈。

图1 研究区域及气象站点分布Fig.1 The study area and the distribution of the meteorological stations in and around Guangzhou

广州的地形地貌以平原、丘陵为主，最高海拔为1 210 m。地处亚热带季风性气候区，雨热同期；多年平均气温为22.5℃，多年平均降水量为1 968 mm，降水主要集中在5—9月；多年平均日照时数为1 643 h。

1.2 数据来源

采用的遥感数据为下载自美国地质调查局Earth Explorer 数据中心的Landsat Collection 2 Level-2 的地表反射率数据集，条带号为122043、122044。该数据通过LEDAPS或LaSRC算法进行大气校正且经过高精度的几何校正，可直接用于遥感解译。此外，部分高分二号遥感影像亦被用于解译标志体系建立、训练样本及验证样本选取等遥感解译辅助工作。

蒸散数据选用下载自美国国家航空航天局Earth Data数据中心的500 m空间分辨率MOD16A2 GF 8 d数据集；遥感卫星轨道号为h28v06，获取时段为2000—2020 年。MOD16 数据是在Penman Monteith 公式的基础上改进获得的蒸散数据产品（Mu et al., 2011）。

MOD16A2GF 是MOD16 数据的改进型，蒸散模拟效果较前代产品更优。自发布以来，MODIS数据产品经过充分验证并在世界范围内广泛运用（Li et al., 2017），姜艳阳等（2017）亦在中国范围内分流域对MOD16 数据集进行质量评估。将MOD16A2GF/ET 原始数据投影至WGS84/UTM 49N 投影坐标系，导入Arcmap 根据广州行政边界矢量文件按掩膜提取，在去除填充值后使用栅格计算器运算得到研究时段内历年广州地表蒸散的空间分布数据；通过统计得出历年广州地表年平均蒸散量以及5期遥感解译结果中各地表覆盖类型的年平均蒸散量和多年平均蒸散量。

蒸散过程不仅受下垫面类型、植被覆盖情况影响，还受供水条件、能量供给条件、水汽输送条件3方面因素控制。选用降水量、气压、温度、风速、日照时数、相对湿度6个主要气象气候因子，以年为单位进行统计，探究其与地表蒸散的相关性。数据主要来自中国气象局国家气象信息中心发布的中国国家级地面气象站基本气象要素日值数据集（V3.0）。由于2020 年4 月之后该数据集不再生产，2020 年4—12 月的数据来自中国气象局全国2 000多个站点的内部网络共享版本数据。

1.3 遥感解译分类体系

在土地城市化和产业城市化进程等社会经济条件的影响下，广州城市绿地景观破碎化程度较高；加之地处亚热带季风气候区，雨热条件充沛，生态演替速度较快，在中分辨率地表资源卫星遥感影像下容易出现不同绿地类型交错混杂分布的小面积斑块；而且，根据广东统计信息网中《关于广东省第二次全国土地调查主要数据成果的公报》（广东省国土资源厅等，2014）数据，在2007—2009 年的调查中广州市内草地面积为49 km2，占研究区域总面积的0.65%；此外，中国科学院土地资源分类系统中所定义的大规模草地在研究区分布较少。因此，结合研究区实际情况，在遥感解译时将广州的地表覆盖分为林地、非林地植被、建设用地或未利用地、水体4种类型，各类型具体定义见表1。

表1 遥感解译分类体系Table 1 Remote sensing interpretation and classification system

1.4 遥感影像解译与统计

城市内部绿地类型破碎化较高且具有高度的空间异质性（林冰钰等，2021），TM、ETM+、OLI图像所能提供的绿地斑块内部纹理信息比较有限。一般情况下，在标准假彩色合成的影像中，林地像元的色调往往较非林地植被像元偏暗。然而，一些处于特定生长阶段的湿地、耕地像元色调可能与林地类似（Huang et al., 2008），易导致遥感解译训练样本选取的偏差。为消除上述误差的影响，采用多时相复合法（Singh, 1989），每期选择2个研究区无云时段的遥感影像进行叠加处理以用于监督分类（表2）。

表2 遥感影像获取时段Table 2 Dates acquired of the remote sensing images

在遥感内业解译工作开始前，于2021年3—11月在研究区以Landsat8 2020-02-18 所获取OLI 影像及高分二号卫星2019-12-12所获取的影像为参考进行外业采样，共采集有效样地107处。在区域层面，样地分布覆盖广州的10个市辖区（花都区和白云区的地表分布类型和模式较为相似，且距离市区较远，遂前期外业调查过程中未前往花都区）；在地面覆盖类型层面，样地覆盖建筑用地、公园林地、公园草地、山体林地、果园、耕地、撂荒地、湿地等研究区常见地表类型。遥感内业解译工作结束后，参考已出版的土地利用图（Fan et al., 2008）及MODIS Vegetation Continuous Field（VCF）全球植被分布数据对研究区的历史时期情况进行大致了解，在Google Earth 对应年份的高分辨率卫星影像上随机选择一定数量的各地类验证样本，对每期分类结果进行基于地面真实情况的分类后精度验证。5 期分类总体精度分别为94.42%、87.92%、84.39%、93.11%、91.55%，kappa系数分别为0.90、0.83、0.76、0.90、0.88。由2020年解译所得的林地面积与遥感作业总面积算得当年广州森林覆盖率为41.75%，与2020年广州市林业和园林局有关绿化统计数据（广州市林业和园林局，2022）公布的41.6%的相对误差仅为0.36%，证明本文遥感影像解译精度可满足研究需求。

1.5 气象数据处理与统计

为了提高气象数据分析统计的精度，选用广州及周边地区共38个气象站点的数据（见图1），采用逐月平均值插补法对各气象站数据进行缺测数据插值；并对各站点历年的年降水量、年平均气压、年平均气温、年平均风速、年日照时数、年平均相对湿度进行统计。基于统计结果，采用普通克里金空间插值法得到广州区域内21年间各年份的空间气象数据，计算其与年蒸散量的皮尔逊相关系数（Pearson, 1920），并进行显著性检验，以探究蒸散量和各气象因子的相关关系。皮尔逊相关系数计算公式为：

式中：Rxy表示变量x和y之间的相关系数；xi表示第i年的气象数据；yi表示第i年的蒸散量；xˉ、yˉ分别表示21 年间各气象数据和蒸散量的平均值；i表示年变量。相关系数Rxy>0 表示两变量存在正相关关系，相关系数Rxy<0 表示两变量存在负相关关系。P< 0.05 表示两变量存在显著相关性，P< 0.01 表示两变量存在极其显著相关性。

2 结果分析

2.1 广州市城市绿地时空变化特征

图2显示，近20年来广州城市绿地总面积呈减少趋势，共减少984.09 km2。就城市绿地系统内部而言，非林地植被覆盖面积总体有所减少，林地覆盖面积有所增加，林地的占比由2000 年的43%提升至2020年的62%。

图2 2000—2020年广州市地表覆盖类型及城市绿地变化Fig.2 The changes of land cover and urban green space in Guangzhou during 2000-2020

以2000、2020年作为始末节点，依据遥感解译分类体系的4种地表覆盖类型转移情况建立转移矩阵，可知：近20年来中4类地表之间互有转入转出（表3），主要表现为城市化进程下大量非林地植被转出为以城乡、工矿、居民用地为主的建设用地或未利用地，一部分转化为林地。

表3 2000—2020年广州市地表覆盖类型转移矩阵Table 3 The transition matrix of land cover of Guangzhou during 2000-2020 km2

2.2 广州市地表蒸散的时空变化特征及其对城市绿地变化的响应

从时间上看，广州年均蒸散量总体呈波动上升趋势。期间最大值出现在2016 年，为774.49 mm；最小值出现在2005年，为600.52 mm（图3）。从空间上看，研究区年蒸散量具有较高的空间分异性，基本呈北高南低的分布趋势。高值区主要分布在从化、增城、黄埔、花都、白云林地覆盖的丘陵地区，低值区分布在以耕地为主的农业生产区和以建设用地或未利用地覆盖为主的城市建成区（图4）。

图3 2000—2020年广州市地表蒸散及气象因子年际变化Fig.3 The inter-annual variations of evapotranspiration and meteorological variables in Guangzhou during 2000-2020

对各非水体地表覆盖类型的年均蒸散量（图5）、多年平均蒸散量进行统计，结果表明不同地表覆盖类型的蒸散能力存在明显差异，多年平均蒸散量由高到低分别为林地（919.10 mm）、非林地植被（612.32 mm）、建设用地或未利用地（303.91 mm），与已有研究（张猛等，2018）所得的大致规律相符。一方面，林木的根系具有更强的吸水能力，可以供给植物进行更强的蒸腾作用（黄葵等，2019）。另一方面，研究区的林地生态群落主要为亚热带常绿阔叶林、亚热带针叶林；在群落的垂直结构上，乔木层以下有林下灌木层、草本层多层复合，其单位面积植物生物量高于单一的草本、灌木群落。综上，本研究区域内林地的蒸散量较其他地类更高。

图5 2000—2020年各地表类型年均蒸散量变化Fig.5 The changes of different land covers' annual average evapotranspiration during 2000-2020

在探究蒸散对城市绿地时空变化的响应时，为了降低气象因子年际变化对分析的影响，将研究时段分为2000—2003、2004—2007、2008—2011、2012—2015、2016—2020 年5 个阶段，各阶段广州地表多年平均蒸散量分别为636.80、624.58、648.63、674.62、719.41 mm。以各期数据中研究区域内林地面积占陆地总面积比例为变量x、对应阶段广州市多年平均蒸散量作为变量y计算皮尔逊相关系数，算得相关系数为0.64，显示广州市地表蒸散量与广州市城市绿地结构的变化存在一定相关性。

以2000、2020年作为始末节点对广州市地表覆盖类型的转移变化情况进行分析，近20年来区域内共有2 555.96 km2的土地发生过地表覆盖变化；其中非林地植被转林地、林地转非林地植被、建设用地或未利用地转非林地植被、林地转建设用地或未利用地、非林地植被转建设用地或未利用地5类地表覆盖变化总计面积为2 285.82 km2，占所有发生地表覆盖变化区域面积的89%。近20年来广州地表蒸散变化量的空间分布与不同蒸散能力地表覆盖类型之间的转化具有较高一致性（图6）。增城区北部、白云区东部及从化区东南部、西部等林地扩张区域地表蒸散量有较大程度的增加。花都区中部、白云区西部、黄埔区和增城区南部、番禺区、南沙区等非林地植被集中转化为建设用地或未利用地的区域，其地表蒸散量出现较大幅度的减少。

2.3 气象因子对广州市地表蒸散的影响

蒸散量与各气象因子的皮尔逊相关分析结果如表4所示。2000—2020年研究区域与蒸散量呈正相关的气象因子为降水量、相对湿度和气温，与蒸散量呈负相关的气象因子为气压、日照时数、风速。

由广州周边气象站点记录分析可知，近20年来区域气象要素都大致向利于地表蒸散量增加的方向变化。但值得注意的是，在此过程中蒸散量与各气象要素间显示显著及以上相关的区域占研究区域总面积的比例偏低。占比最高的要素为相对湿度，达到48.70%；占比最低要素为气温，仅为0.78%；各要素间平均值仅为21.33%。

Li等（2017）在对2001—2013年中国地表蒸散研究中指出，该时段气候变化对蒸散的影响较土地利用和土地覆盖变化更为显著。但本文结果（表4）显示：虽然气象因子对地表蒸散过程具有不可忽视的重要影响，但在城市尺度下气象因子波动并不能被解释为近20年来广州市大部分区域蒸散量变化的主要驱动因素。由此可见，2000—2020年广州城市尺度下地表覆盖变化对地表蒸散的影响或更为显著。

表4 2000—2020年广州地表蒸散量与各气象因子相关性面积占比Table 4 The area percentage of different correlation between ET and meteorological variables during 2000-2020 in Guangzhou %

3 讨论

数据结果显示，2000—2020年广州林地规模呈增长趋势，从2000 年的2 534 km2增长至2020 年的3 031 km2。这与近20年来广州市相关规划与林业生态建设工作密切相关。期间广州市逐步推进包括迹地更新、低效林改造、荒山造林、中幼林抚育、大田育林等造林绿化工作。2003 年广州市宣布启动“青山绿地”工程，重点工作为城市林带、城市林区、园林的建设和采石场整治（广州年鉴编纂委员会，2004）；2007年第二期“青山绿地”工程启动，规划范围扩大至全市建成区以外的全部范围，建设内容拓展至沿海防护林体系建设、新农村绿色家园建设、城市森林修复、花卉苗木产业等（广州年鉴编纂委员会，2007）。通过统计《广州市年鉴》2004—2011 年林业生态建设板块相关数据得出，2003—2010 年两期“青山绿地”工程三年计划期间，广州共完成造林绿化规模约228.73 km2；为协调城市整体生态文明发展，该工程重点推进流溪河水源涵养林、增江河水源涵养林的建设、林分改造及保护工作；其在空间上与本研究监测到的林地增加热点区域相吻合。

随着“绿色亚运”“花园城市”、创建国家森林城市等工作的推进，2000—2020 年广州市公益林、碳汇林、城市生态景观林带、森林公园建设任务也逐步落实。此外，加强森林资源管护、推进森林病虫害检疫防治、完善林业法治、加强森林公安执法队伍建设、完善公益生态林补偿办法、建设完善森林生态监测体系等举措也是广州森林生态建设、绿化造林工作推进的重要保证。

值得注意的是，2000—2005年广州市林地增加规模为438 km2，占21 年间林地面积增长总量的88%。但2004—2007年广州多年平均蒸散量并没有出现明显上升，反而较前一阶段有所下降；这一现象可能与森林的生长周期有关。李根等通过研究2001—2013 年中国土地利用/覆被变化对蒸散和气候变化的影响发现，森林转化为其他土地利用类型引起的蒸散量化远大于由其他土地利用类型转化为森林引起的蒸散量变化（Li et al., 2017）。苏联科学家krestovsky 通过对森林采伐后恢复期的年均蒸散量和径流量变化统计得出krestovsky 生态水文函数（Shiklomanov et al., 1988），该函数揭示了在森林恢复初期，区域内蒸散量并不会有显著的增加，反而在最初5—10 年可能有所减少，此后才逐渐上升，并在50—60年后达到干扰前的峰值，其后保持稳定或有所下降。而对于不同林分、不同物种组成的森林生态系统，其生态水文函数曲线的具体形态和特征也存在差异。Molchanov 等（1963）的研究也显示，苏联的橡树、云杉的蒸散量峰值出现在60岁，梣木和白杨树的峰值出现在20～40岁。随着林木生长时间的增加，其后3个阶段广州多年平均蒸散量也呈现递增趋势。

叶面积指数（LAI）可以定量描述植物群落水平上的叶片密度，是研究植物群落的重要参数（王希群等，2005），其空间分布格局变化在一定程度上反映区域植被覆盖度的变化和植被恢复状况。为印证研究区域森林生长周期相关问题，笔者在增城区北部划出一面积为381 km2的验证区域。区域内2000年林地面积为167 km2，占区域总面积的43%；2005 年林地面积为252 km2，占比提高至66%；2006—2020年由于区域开发、设施建设等原因林地面积稍有回落，在225 km2上下波动。

对验证区域内2000、2005、2010、2015、2020年5期MOD15A2H叶面积指数求年平均值进行统计分析。图7-a显示，近20年增城北部集中造林区内LAI的空间分布格局较为稳定，并无发生剧烈变化。图7-b 显示，区域内各像元LAI 值的相对频率分布呈双峰形态，高峰出现在值为1附近的低值区，次高峰出现在值为3～4.5 附近的高值区；图中出现的高值区次高峰对应2000年区域内植被覆盖度较高的像元，低值区的高峰及其移动则反映2000—2005年新增林地的恢复情况。2000、2005、2010年对应相对频率曲线的低值区高峰具有较高重叠度，该时段新增林地的LAI 并无明显增加（图7-a）。而2015、2020年对应频率曲线低值区高峰向高值方向发生整体偏移；由此可知，规模造林后，确实存在一定时间阶段区域内群落叶片面积密度并无明显的增加，即植被尚未进入快速恢复的阶段。

图7 2000—2020年增城北部集中造林区域LAI空间分布及相对频率曲线Fig.7 LAI's spatial distribution and relative frequency curve of the concentrated afforestation area in the north of Zengcheng during 2000-2020

气象条件的波动也可能是导致第二阶段广州多年平均蒸散量偏低的原因之一。从下垫面供水条件的角度看，2004—2007年广州的年降水量和年均相对湿度都持续处于相对低位，可见该时期由大气提供的水汽补给相对较少；与此同时，该时期年平均风速较21年间大部分年份相对偏高（见图3）。上述因素可能导致当年的土壤水分不足，从而对地表植被蒸散作用起限制作用。从植物的能量获取角度看，2004—2007年广州年均温度较前一阶段有所降低，这也可能导致第二阶段广州市多年平均蒸散量较第一阶段有所下降。

此外，由于其在空间分辨率上的历史局限，选用长时间序列中分辨率遥感数据无法在城市尺度下对地表环境进行细致的解译；加之大尺度蒸散数据在不同小尺度区域研究中存在无法避免的误差，本研究更多偏重于在宏观尺度上对相关现象、趋势的阐述，在定量上存在不足。未来可结合实地测量数据对现有的蒸散模型进行优化，确保其在小尺度区域上仍能有较高的精度；并辅以更高分辨率的遥感数据及能表现城市绿地多样性的解译结果建立模型，以更精确量化城市尺度下地表覆被变化对蒸散等生态水文过程的影响。