2000–2020 年珠三角30 m 空间分辨率植被覆盖度时空变化数据集

沈明潭，谭炳香*，侯瑞霞，黄逸飞

1.中国林业科学研究院资源信息研究所，北京 100091

引 言

珠三角在国家整体发展中具有重要的战略地位，自改革开放以来经历了快速的经济增长、人口增长和城市化扩张，导致当地植被急剧退化和破坏，从而严重影响了生态系统的服务和功能，因此动态监测珠三角植被的时空变化和分布格局对当地的生态建设具有重要的参考意义。在已公开发布包含珠三角地区的植被覆盖度数据集产品中，高吉喜等基于空间分辨率为250 m 的MODI3Q1 NDVI（Normalized differential vegetation index）数据，采用NDVI 最大合成法的方式进行中国区域的植被覆盖反演[1]。与该数据集相比，本研究选用空间分辨率更高的30 m Landsat5 TM 和Landsat8 OLI 数据，利用中值合成法获取研究区的影像数据，更能精确估测出珠三角地区的植被覆盖度，更科学有效地反映出珠三角地区的植被覆盖状况。

植被覆盖度（Fraction Vegetation Cover, FVC）是指地面植被（包括树叶、茎和树枝）的垂直投影面积占统计区域总面积的百分比，可以反映地表植被的生长状况，是监测生态环境的重要基础数据，在评估地表覆盖变化、土地利用方式、景观和生态系统功能等方面有重要指示作用[2]，被广泛用于监测植被变化、区域环境质量评估、气候预测、土壤侵蚀评估等。遥感技术因其获取速度快、覆盖范围广、能够提取大量信息而广泛应用于研究植被动态监测、生态环境评估和土地覆盖分类等领域[3]，是监测区域和全球FVC 最高效的方法[4]，可以动态监测植被覆盖变化及其影响因素。通过长时间序列和大规模的FVC 研究来评估特定区域的生态环境条件和质量，在揭示区域环境条件的演变过程和预测未来变化趋势等方面具有重要意义。

1 数据来源和处理方法

本研究采用以下主要步骤对珠三角植被覆盖度进行分析：（1）在GEE（Google Earth Engine）平台调用Landsat 数据，基于NDVI 像元二分法估算珠三角2000、2005、2010、2015 和2020 年植被覆盖度；（2）利用变异系数分析珠三角5 个时期植被覆盖度变化的波动性；（3）利用Slope趋势分析来探究植被覆盖度变化的时空变化特征；（4）采用Hurst 指数预测研究区植被覆盖度未来的发展趋势。

1.1 数据来源

Landsat5 TM 和Landsat8 OLI 数据具有较高的空间分辨率（30 m），但由于研究区地处亚热带区域，受到长时间大面积高云量的影响，难以获取研究区完整清晰的影像。因此本研究基于GEE 平台，利用一年中云量低于20%的数据，按照成像时间进行排列，选取位于中间时段的像元，合成出质量较高的研究区影像（见表1）。调用“LANDSAT/LT05/C01/T1_SR”数据集，分别利用2000-01-01 至2000-12-31 期间的141 景Landsat5 TM 影像合成2000 年珠三角年度中值影像，2005-01-01 至2005-12-31 期间的119 景Landsat5 TM 影像合成2005 年珠三角年度中值影像，2010 年的数据量较少，需要扩大时间段才能得出质量较好的影像数据，因此利用2009-11-01 至2010-12-31 期间的101 景Landsat5 TM 影像合成2010 年珠三角年度中值影像。调用“LANDSAT/LC08/C01/T1_SR”数据集，分别利用2015-01-01 至2015-12-31 期间的162 景Landsat8 OLI 影像合成2015 年珠三角年度中值影像，2020-01-01 至2020-12-31 期间的152 景Landsat8 OLI 影像合成2020 年珠三角年度中值影像。由于按照时间段内云量低于20%的数据进行筛选，每个时间段内符合的景数不同，且能完整合成出研究区清晰无云的影像数据需要的景数也不同，因而每个期间所调用的影像景数会有差别。所用的Landsat 数据集都已经过辐射定标和大气校正，分别裁剪出珠三角5 个时期的合成影像（见图1）。

图1 珠三角Landsat8 标准假彩色影像Figure 1 Standard false colour images of Landsat8 in the Pearl River Delta

表1 影像合成信息Table 1 Image synthesis information

1.2 数据处理方法

1.2.1 植被覆盖度估测

本研究基于GEE 平台生成珠三角年度中值植被覆盖度影像[5-6]。主要的计算过程为：（1）输入珠三角范围；（2）构建Landsat 遥感影像的时间序列数据集。分别调用LANDSAT/LT05/C01/T1_SR、LANDSAT/LC08/C01/T1_SR 数据集，根据时间段进行筛选数据，然后进行去云处理；（3）生成珠三角年度中值遥感影像数据集。对构建后的影像数据集按不同的年份进行分组，逐像元计算生成每个时期的年度中值影像；（4）计算NDVI。调用normalized Difference 函数计算NDVI（公式1）；（5）去除水体。为了避免水体对植被覆盖度估测的影响，调用mask 函数将水体进行掩膜处理；（6）计算珠三角FVC。选用基于NDVI 的像元二分法（公式2），分别计算每个时期的植被覆盖度；（7）导出计算结果。

1.2.2 归一化差值植被指数

Kriegler 等人[7]提出了一种简单的波段变换：用近红外光反射率减去红光反射率，再除以近红外光反射率加红光反射率，称为归一化差异植被指数（NDVI）。NDVI 对植被的生长状态和空间分布特征非常敏感[8]，可以在很大程度上消除地形和大气的干扰，并与植被覆盖度有显著的线性相关[9]，在大规模植被动态监测等方面具有明显优势，NDVI 值越高，植被覆盖度越高，是目前应用最广泛的植被指数[10]。计算公式为：

式中，NDVI是归一化差异植被指数，NIR和R分别为传感器在近红外光和红光区域记录的反射率。

1.2.3 像元二分模型

植被覆盖度与NDVI 之间具有显著的线性相关[11]，利用植被指数与像元分解模型相结合的方法，具有一定的物理意义、反演算法简单[12]，而且在很大程度上弱化了大气、土壤背景等因素的影响，被广泛应用于区域和全球范围内植被覆盖度提取[13-14]。因此本研究利用基于NDVI 的像元二分模型[15]对研究区植被覆盖度进行估算：

式中FVC为植被覆盖度，本研究选取置信区间为[5%, 95%]的NDVI 值作为NDVIsoil和NDVIveg[16]。NDVIsoil表示无植被覆盖像元NDVI 值，NDVIveg表示纯植被覆盖像元NDVI 值。

1.2.4 变异系数

对研究区5 个时期的植被覆盖度数据逐像元计算变异系数，探究研究区植被覆盖度的稳定性[17]。CV 值越大，植被变化的波动性越大，将波动性划分为4 种类型，见表2。计算公式如下所示：

表2 植被覆盖度波动变化类型划分表Table 2 Classification of FVC fluctuation types

式中CV植被覆盖度的变异系数，n 是研究期数，FVCi第i 年的植被覆盖度值，F̅̅̅VC̅̅̅是研究期内植被覆盖度均值。

1.2.5 Slope 趋势分析

采用一元线性回归分析法中的斜率公式计算要素（y）的变化趋势（Slope）并进行显著性检验，根据Slope值和P 值对植被覆盖度的变化趋势进行划分，将变化趋势分为7 种类型[18]，见表3。其中Slope的计算公式如下：

表3 植被覆盖度变化趋势类型划分表Table 3 Classification of the trend types of FVC changes

式中，Slope表示要素的年际变化趋势；i为时期变量；n是研究期数，本研究为5 个时期；yi为第i年的要素。Slope>0 和Slope<0 分别表示要素序列随着时间增加和减少，Slope的绝对值越大表明要素变化越大。

采用F 检验对变化趋势进行显著性检验，当检验结果P<0.05 为显著变化，0.050.1 为不显著变化。统计量计算公式如下：

式中，U为回归平方和，Q为残差平方和，yi为第i年的要素，ŷi为要素的回归值，y̅为要素的平均值，n表示期数。

1.2.6 Hurst 指数

基于R/S 分析法的Hurst 指数[19]，能够用于定量描述时间序列信息的长期相关性，目前已广泛运用于分析区域植被覆盖未来变化趋势，Hurst 指数取值范围存在以下3 种情况：（1）若0.5

时间序列{FVC(t)},t=1, 2, 3,…,n，对于任意一个正整数τ=1，定义均值系列：

累计离差：

极差：

标准差：

考虑比值R(τ)/S(τ)≅R/S，若存在R/S∝τH，则说明时间序列{φ(t)}, t=1, 2, 3,…,n存在Hurst 现象，其中H 值称为Hurst 指数。参考刘明霞等[21]将植被覆盖度未来变化趋势分为7 类，见表4。

表4 植被覆盖度未来发展趋势类型划分表Table 4 Classification of future development trends of FVC

2 数据样本描述

2.1 数据组成

本数据集包括：（ 1 ） PRD_FVC_2000.tif 、 PRD_FVC_2005.tif 、 PRD_FVC_2010.tif 、PRD_FVC_2015.tif、PRD_FVC_2020.tif，分别为2000 年、2005 年、2010 年、2015 年和2020 年5 个时期的年度植被覆盖度数据；（2）PRD_FVC_CV.tif 为植被覆盖度变异系数数据；（3）PRD_FVC_Slope_F.tif 为植被覆盖度变化趋势数据；（4）PRD_FVC_Hurst.tif 为植被覆盖度未来变化趋势数据。上述数据空间分辨率为 30 m，地理坐标系为 WGS1984，投影坐标系为WGS_1984_UTM_Zone_49N，数据保存为GeoTIFF 格式，整个数据集解压后大小约为2.75 GB。

2.2 样本描述

基于GEE 平台，利用公式1 和公式2 进行珠三角植被覆盖度反演，得到2000 年、2005 年、2010年、2015 年和2020 年5 个时期的年度植被覆盖度影像（见图2）。基于珠三角5 个时期的植被覆盖度反演结果，利用公式3 对珠三角植被覆盖度变化的波动性进行分析，根据表2 将变异系数划分为4 种等级，结果见图2(a)；利用公式4 计算出珠三角5 个时期的植被覆盖度的时空变化趋势，采用F检验（公式5、6 和7）对变化趋势进行显著性检验，根据表3 的将变化趋势划分为7 个等级，结果见图2(b)；利用公式4 计算出珠三角植被覆盖度的变化趋势，再利用Hurst 指数（公式8、9、10 和11）来反映植被覆盖度变化趋势的持续性，从而揭示珠三角植被覆盖度未来的发展趋势，根据表4 将珠三角发展趋势划分为7 种类型，结果见图2(c)。

图2 珠三角5 期植被覆盖度（2000、2005、2010、2015、2020）、变异系数(a)、变化趋势(b)、Hurst 指数(c）空间分布Figure 2 Spatial distribution of vegetation coverage in five phases (2000, 2005, 2010, 2015, 2020), coefficient of variation (a), trend (b) and Hurst index (c) in the Pearl River Delta

3 数据质量控制和评估

本数据集参考植被覆盖度的相关研究成果，在GEE 平台利用中值合成法得到研究区清晰无云的影像，提高反演数据源的质量和精度。植被覆盖度的实际采样点观测值与像元二分法的估测值有较好的一致性，利用像元二分法对植被覆盖度进行反演是科学可靠的方法[22-23]。因此选择基于NDVI的像元二分法进行珠三角地区植被覆盖度反演，根据研究区的气候条件和植被覆盖状况，在5%的置信水平下，采用累计频率在5%～95%的NDVI 值作为上下阈值，在一定程度上减少遥感影像的噪点影响[24]。基于NDVI 的二分像元模型对图像辐射校正的影响不敏感，具有较高的精度和易操作性，不需要实际的植被覆盖度数据建模，仅需简单计算，易于解释[18]，目前广泛应用于区域、省、流域、国家或全球范围等不同尺度的植被覆盖变化研究。参考Song 等学者的研究[25]，利用相关系数R2和均方根误差RMSE 进行产品的精度验证和一致性检验。

3.1 精度验证

本研究在珠三角范围内随机生成100 个采样点，通过world imagery 2020-07-22 高分辨率影像进行精度验证，将相同坐标位置的2020 年植被覆盖度遥感反演值和高分辨率影像目测值进行回归分析，两者相关系数R2=0.90，均方根误差RMSE=0.03，表明本数据集的精度基本满足区域性植被覆盖度研究的需要。

3.2 一致性检验

关于产品的一致性检验，目前公开发布的包含珠三角地区的植被覆盖度产品都为粗分辨率的，比如GEOV2 的2000 m 植被覆盖度产品，GLASS 的500 m 植被覆盖度产品，高吉喜团队的中国区域250 m 植被覆盖度数据集（2000–2022）。本数据集以高吉喜团队2020 年植被覆盖度数据为对比产品，在珠三角范围内随机生成1000 个点，提取本研究2020 年植被覆盖度数据和对比产品相同位置的值进行一致性检验，相关系数R2=0.69，均方根误差RMSE=0.09，表明两个产品整体上有较好的一致性，但根据world imagery 2020-07-22 影像可以看出，本研究的FVC 数据相较于对比产品表现出了更多的细节（见图3）。一致性不高的原因有：（1）虽然都采用了像元二分法，但NDVI 的取值方式不一样，对比数据集采用的是NDVI 最大值合成法，而本研究利用中值法合成影像，然后再计算NDVI；（2）两个数据集的空间分辨率不一样；（3）数据源不一样。

图3 world image 影像（a）、本研究2020 年FVC（b）、对比产品2020 年7 月FVC（c）Figure 3 world image data (a), 2020 FVC (b) of this study, July 2020 FVC (c) of comparative products

4 数据价值

本数据集提供了珠三角2000、2005、2010、2015 和2020 年5 个时期植被覆盖度的反演结果、变异系数、变化趋势和Hurst 指数，可用于定量研究珠三角植被覆盖度的时空变化特征和预测未来变化趋势，探究珠三角植被覆盖度变化的机理，更好地了解珠三角植被变化的动态过程，对植被资源的管理及生态环境的保护具有重要的指导意义。在珠三角快速的城市化进程中，本数据集具有重要的实际应用价值。在城市规划方面，可以评估城市化对植被生态系统的影响，有助于指导城市的规划和发展。在生态保护建设方面，植被覆盖度研究有助于科学地认识生态系统及生态过程，可以指导生态环境的保护和恢复工作。在气候变化研究方面，监测植被变化可以深入了解气候变化对生态系统的影响，为可持续发展提供关键信息。

5 数据使用方法和建议

本数据集以GeoTIFF 格式进行存储，在常用的GIS 软件（ArcGIS, QGIS 等）及其他数据分析软件（R, MATLAB 等）都可进行查看和再分析，可以实现2000、2005、2010、2015 和2020 年珠三角植被覆盖度的定量评价和时空变化特征分析，可以应用于珠三角植被覆盖度动态变化与影响因子之间的相关性探究，以及定量分析自然因子和人为因子对珠三角植被覆盖度变化的相对影响程度等方面的研究。植被覆盖度的取值范围为[0,1]，为了监测地表植被覆盖细微的时空变化，同时参考相关的植被覆盖度数据集[1]，本研究FVC 值保留6 位有效数字，用户在使用过程中可根据实际需要自行决定保留几位有效数字。

数据作者分工职责

沈明潭，男，云南文山人，硕士研究生，研究方向为林业遥感技术与应用。主要承担工作：数据处理，数据论文撰写。

谭炳香，女，北京市人，博士，研究员，研究方向为森林资源遥感监测与森林信息遥感提取等。主要承担工作：数据论文修改。

侯瑞霞，女，山西翼城人，博士，副研究员，研究方向为GIS 与数据分析。主要承担工作：数据论文撰写建议。

黄逸飞，女，湖南邵阳人，硕士研究生，研究方向为林业遥感技术与应用。主要承担工作：基础数据收集。