山东省植被覆盖度长时序动态监测与分析

刘晨曦,陈关州,王 铜,张晓东,*,朱 坤,王 庆

1.长江大学地球科学学院,湖北 武汉 430100

2.武汉大学测绘遥感信息工程国家重点实验室,湖北 武汉 430079

植被覆盖度（Fractional Vegetation Coverage，FVC）指植被（包括叶、茎、枝）在地面的垂直投影面积占统计区总面积的百分比，是目前描述植被群落及生态系统较为简单和有效的生物物理参数[1]。近年来，由于人类过度开垦、城镇扩张等因素[2]，我国部分地区面临着日益严峻的生态环境问题[3]。长时间、大范围监测植被覆盖度，对加强生态文明建设，形成合理的生态布局，帮助我国早日实现碳达峰、碳中和具有重要意义[4]。

传统的植被覆盖度监测方法一般采用目测法、采样法以及仪器法[5,6]，但这些方法不仅耗时耗力还会受到天气的影响，同时精度也难以保证，更无法实现监测结果时空的连续性。遥感数据因其时效性强、覆盖范围广等特点，已被广泛用于连续时间序列、多区域大尺度植被覆盖的相关研究[7]。在空间尺度方面，经典的遥感影像植被覆盖度监测方法大多基于县域、市域[8,9]展开，对大范围区域研究较少；在时间尺度方面，目前的研究大部分为基于间隔年份、或相对较短时间年份[10,11]。同时，遥感影像常用的单次单景的下载方式，也造成了数据处理步骤繁琐、影像处理周期过长、人力成本较高等问题。

GEE（Google Earth Engine）、笛卡尔实验室Descartes Labs 以及中国的地球大数据挖掘分析系统Earth Data Miner 和PIE-Engine 等遥感云平台提供了全球尺度、多源遥感影像数据，其强大的计算能力、日更新速度及大范围云计算资源[12,13]能辅助进行快速、大尺度的植被覆盖度监测。相较于经典的遥感影像处理模式，全新的遥感大数据处理平台解决了以往大尺度地学分析应用研究中收集数据难、解译效率低等弊端，可不受时间和空间限制，快速、批量处理时序遥感数据[14-16]。

Vermeulen LM 等基于GEE 利用有限的野外数据和开源遥感图像提供了一种进行FVC 估计的稳健方法[15]。Michelle Ern Ang Li 等基于GEE 平台首次绘制了具有粗主题分辨率的高水平土地覆盖图[17]。Daojun Zhang 等利用Moran 指数分析了空间限制下，植被间对水分和养分竞争的现象，并利用Hurst 指数分析了未来植被覆盖的发展趋势[18]。王华、刘垚燚等利用GEE 平台及Landsat 遥感影像，在短时间内分别分析了所研究地区约30 年间的植被覆盖度的时空变化特征与驱动力[19,20]。上述方法表明，利用GEE 为代表的遥感大数据平台进行植被覆盖度长时序动态监测是可行的、高效的。

本文以山东省为实验区域，基于GEE 云平台获取所需遥感影像，结合像元二分模型，分析植被覆盖时空变化特征，通过分析变异系数与Hurst 指数，探究近20 年来该省植被覆盖的变化趋势。本文对山东省植被覆盖变化实现了长时序动态监测，可以为山东省合理处理经济发展与生态建设的关系提供科学依据和辅助决策支持。

1 区域概况与数据来源

1.1 研究区域概况

山东省（北纬34°22.9′～38°24.01′，东经114°47.5′～122°42.3′，面积约为157900 km²），是位于我国华东地区的一个沿海省份，也是我国13 个重要的粮食主产区之一。据《山东省第一次全国地理国情普查公报》报告，山东省植被覆盖比例达77.54%，且陆地的生态系统复杂多样。

1.2 影像数据源

本文基于GEE 云平台，在原始数据集中筛选出云量相对较少的Landsat 系列、MODIS、Sentinel 2 影像集，获取了2000 年1 月至2020 年12 月山东省内共27614 景影像，其中包括Landsat 5、Landsat 7、Landsat 8 影像2704 景、5553 景、2616 景，以及Sentinel-2 卫星影像9504 景和MODIS 影像7237景。经过去云、拼接等处理，并针对不同分辨率的影像，统一采样到250 m 分辨率，构建成了山东省长时序的多光谱影像数据序列。

2 研究方法

本文针对长时序大范围的海量遥感影像处理，基于GEE 平台收集山东省2000-2020 年间的MODIS、Landsat 系列、Sentinel-2 等多源遥感影像数据并进行去云、镶嵌、融合等数据预处理；提取山东省21 年的长时序植被覆盖度，并与已有土地覆盖数据集GEOV1 FVC 进行精度评定；在此基础上对山东省21 年间的植被覆盖度进行上述监测分析（方法包括：最小二乘法、空间自相关、变异系数、Hurst 指数）。

2.1 植被覆盖度计算

本文通过计算每幅影像中的每个像元的归一化植被指数（Normalized Difference Vegetation Index，NDVI），在计算NDVI 时，先计算各个数据源的NDVI 再直接聚合得到同一分辨率的NDVI。采用NDVI 最大值合成方法获取研究范围内年度NDVI 值，再利用像元二分模型估算植被覆盖度，最后计算不同像元对应的年度FVC 合成数据。具体流程如图1 所示。

图1 植被覆盖度计算流程图Fig.1 FVC calculation flow chart

2.1.1 NDVI 计算NDVI 能有效消除因扫描时传感器自身角度和大气辐射等因素的干扰，适用于区域性尺度的植被动态监测和对FVC 的反演及其变化，其波段计算表达式为：

式中，NIR表示近红外波段的反射率值，R表示红波段的反射率值[21]。

2.1.2 基于最大值合成的多时相NDVI 融合 为消除云、大气、太阳高度角、卫星的观测角度、不同来源传感器的波段响应函数不一[22,23]对NDVI 的“衰减”作用[24]，最大限度满足植被信息提取需求，本文对每年的多时相NDVI 数据进行逐像元最大值合成，为接下来的特征信息提供数据基础。

2.1.3 基于像元二分模型的FVC 估算 本文依据反演模型自身所具有的普遍性和特性，选择精度较好的混合像元二分模型。像元二分模型计算公式为：

式中，S为地表光谱信息，Sv为植被地表光谱信息，Ss为非植被地表光谱信息。

对于仅包含非植被与植被两种成分构成的混合像元，植被信息可表示为纯植被像元贡献的信息Sveg和FVC的乘积，非植被部分的信息即为纯非植被光谱信息Snveg和（1-FVC）的乘积，即:

其中，Sveg包含了植被的各方面信息，如植被结构、植被类型等；Snveg包含非植被区域内的各类信息，如土壤类型、湿度、亮度、颜色等。

2.1.4 基于NDVI 的FVC 估算完整流程NDVI 像元二分模型估算植被覆盖度的相关研究中己被广泛采用，将NDVI 计算公式代入（4）可近似得到:

式中，NDVI为混合像元内的植被指数值；NDVIveg为纯植被像元内的植被指数值，NDVInveg表示纯非植被像元（例如土壤像元）的植被指数值。理论上，NDVI的最小值为NDVInveg，最大值为NDVIveg，但由于噪音点和影像上NDVI 的灰度干扰产生异常值，须选定合适的置信区间来确定NDVIveg和NDVInveg[20]。通过总结前人研究经验和研究区实际情况[25,26]，本文对利用像元计算得到的NDVI值进行排序，根据计算结果选取前5%对应的NDVI 值为NDVInveg，后95%所对应的NDVI 值为NDVIveg。为了分析研究区FVC 在空间上的分布特征、时间上的变化特征，根据本研究区内植被生长特征以及植被覆盖实际情况，对山东省植被覆盖度划分为5 个等级[27,28]。植被覆盖度分级表如表1 所示。

表1 研究区植被覆盖分级表Table 1 FVC classification standard in the study area

2.2 植被覆盖度精度评定方法

本文使用GEOV1 FVC 数据（2013 年内与2017 年）与本文计算得到的FVC 数据进行对比验证，评定指标为相关系数[29]（http://land.copernicus.eu/global/products/fcover）。

2.3 植被覆盖度变化趋势分析

对于获取的FVC 数据，本文采用最小二乘方法、可视化等方式，分析研究区域内FVC 变化趋势；而后通过空间自相关分析，评估植被空间分布模式；最后结合变异系数波动程度、Hurst 指数的长程依赖性对植被覆盖度的稳定性与持续性进行分析。流程如图2 所示：

图2 FVC 变化趋势流程图Fig.2 Flow chart of FVC change trend

2.3.1 基于最小二乘法的相关分析 本文采用线性回归形式，即最小二乘线性拟合来分析每个像元FVC 的变化趋势[30]，其计算公式：

式中，Slope 为回归方程的斜率，表示研究时序内植被的变化趋势；n为监测年份总数，i为监测年序号，FVCi表示第i年的植被覆盖度。Slope＞0 表示研究区植被覆盖度为增长态势；Slope＜0 意味着植被覆盖度呈下降态势；Slope 接近0，植被覆盖度变化不明显。

2.3.2 空间自相关分析 为评估植被覆盖度的分布模式，本文使用空间自相关分析的方法，通过计算Moran 指数值，使用z得分和p值来对该指数的显著性进行评估。空间自相关的Moran 指数统计可表示为：

其中，zi是要素i的属性与其平均值的偏差，wi,j是要素i和j之间的空间权重，n等于要素总数，So是所有空间权重的聚合：

统计的zI得分按以下形式计算：

其中：

当Moran 指数＞0 及z得分越大时，表明要素的分布呈现聚集的模式分布，反之则说明其更加趋向于在空间上离散分布。当Moran 指数=0 时，空间分布呈随机性[31,32]。

2.3.3 基于CV 指数的相对波动程度分析 变异系数指数（Coefficient of Variation Index，CV 指数）可用于揭示栅格像元值的相对波动程度[33]。CV 指数值越大，则表明植被生长季FVC 像元值波动性越大，FVC 值越不稳定；反之，则说明时间序列数据分布较为集中，FVC 越稳定。计算公式如下：

式中，CV指数为栅格像元值的变异系数；FVCi表示第i年植被覆盖度的最大值；为研究区每个栅格单元的植被覆盖度在2000-2020 年内的平均值。

2.3.4 基于Hurst 指数的长程依赖性分析Hurst 指数（H 值）是定量描述时间序列信息长程依赖性的有效方法之一[34]。本文中的Hurst 指数采用“极差/标准差”（R/S）分析法进行计算[35]，用来表示研究区FVC 时间序列的变化持续性。具体计算流程如下：

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

式中：H为Hurst 指数，可在双对数坐标系中用最小二乘法拟合求得，H值介于0～1 间，当H=0.5时，认为序列没有自相关性，H＞0.5 时，序列有很强的正相关性，当H＜0.5 时，序列有很强的负相关性。

3 结果与分析

3.1 FVC 提取精度评定

本文将所计算的FVC 数据与GEOV1 FVC 数据进行验证。计算得出的FVC 数据与GEOV1 FVC数据的相关系数在2013 年为0.62，在2017 年为0.64，基本符合本文的精度要求。上述结果表明了本文提取FVC 数据的有效性，后文将在此基础上进行时序动态监测分析。

3.2 植被覆盖度时间变化特征

对2000-2020 年中每年12 个月份的NDVI 影像做最大值合成处理，取最大值作为每年的NDVI值，基于改进的像元二分模型计算植被覆盖度，对植被覆盖度变化趋势情况进行分析，计算得到的全省年均植被覆盖度和年际变化率如图3 所示。

图3 2000 年-2020 年山东省FVC 变化趋势Fig.3 Change tendency of the FVC in the Shandong province during 2000-2020

从图3 中可以看到，山东省2000-2020 年FVC 整体年均增长速率约为0.4%。2020 年达到植被覆盖度的最高值（71.5%），2008 年是植被覆盖的最低值（60.3%），平均值为65.6%。简而言之，山东省近21 年的FVC 呈现出波动增长的趋势。为了进一步分析山东省内部植被覆盖度的详细变化，针对每个市的年均植被覆盖度和年际变化率分析，结果如图4 所示。

图4 2000 年-2020 年各市年均FVC 及其变化趋势Fig.4 The average annual FVC of each city and its changing trend from 2000 to 2020

3.3 植被覆盖度时空分布特征

图5 为依据表1 制作的山东省21 年各年度的最大植被覆盖度分级图。

图5 2000 年-2020 年植被覆盖度分级图Fig.5 Classification map of FVC from 2000 to 2020

根据植被覆盖度划分等级，统计每个年度各等级植被覆盖度的面积及比例，结果如图6 所示。

图6 2000-2020 年不同等级植被覆盖度的区域占比情况Fig.6 Percentages of FVC under different levels during 2000-2020

由图5、图6 可以看出，根据植被覆盖度空间分布差异，可以将2000 年至2020 年植被覆盖度的变化分为3 个阶段。第一阶段为2000-2006 年，山东省低植被覆盖分布范围为中部地区（主要为济南、淄博、潍坊地区）、沿海地区（主要为盐场，包括烟台、滨州、东营）。整个研究区域，高植被覆盖区占比平均值为32.0%；低植被覆盖区占比平均值为10.4%。第二阶段为2007-2012 年，此时山东省中部地区存在部分低植被覆盖区域并且出现低植被覆盖区域往南部地区扩展现象，沿海地区植被覆盖度已经相对升高，高植被覆盖区仍分布最广，占比为30.5%。第三阶段为2013-2020 年，此阶段总体植被覆盖度较高，高植被覆盖区占比为43.7%。低植被覆盖区、中低植被覆盖区、中度植被覆盖区、中高植被覆盖区的平均占比依次为8.0%、8.2%、15.1%、25.0%。

总体而言，低、中低植被覆盖范围随着经济发展与生态保护政策实施下，先略微上升后下降。到18、19、20 年山东省植被覆盖整体呈现以高、中高覆盖为主，其他覆盖区零星分散分布现象。

3.4 植被覆盖度时空变化趋势

本文采用最小二乘法对每个像元2000 年至2020 年间的植被覆盖度变化趋势进行直线拟合，最终获得每个像元植被覆盖度随时间的变化趋势。在具体实验中，首先每三年进行一次拟合，最后对这21 年的所有数据进行拟合。根据检测结果将Slope 分为5 个等级，如表4 所示。

表4 研究区变化趋势分级表Table 4 Slope classification standard in the study area

在此基础上，本文绘制了如图7 所示的山东省植被覆盖度变化趋势的分级分布图，以直观反映时空变化趋势。

图7 植被覆盖度变化趋势的分级分布Fig.7 Hierarchical distribution of the changing trend of FVC

为了进行定量分析，对每个变化等级所占的面积比例进行统计，得到结果如图8 所示。

图8 2000-2020 年植被覆盖变化趋势面积比例统计Fig.8 Statistics of changing area ratio for FVC during 2000-2020

由图7、图8 可以看出，2000-2020 年，山东省植被覆盖度变化呈现基本不变和轻微改善的特征，其中基本不变面积占比76.3%，轻微改善面积占比15.4%，轻微退化面积占比8.2%，退化面积占比0.1%，改善面积可忽略不计。

由图7、图8 可以发现，在以三年为间隔的变化分析中，2000 到2020 每三年植被覆盖度变化趋势呈现相对较大的变化，不同区域呈现不同范围、程度的增长和减少，而在全部21 年的变化分析中，无显著变化区域的面积占75%以上，这是因为在不同时期植被的退化和增加区域有所不同，在21年的长时序尺度上互相抵消，这也与4.5 节中Hurst 指数反应的变化区域相互验证，表明该结论的内在自洽性。

3.5 植被覆盖度的空间分布模式

针对植被覆盖度的空间分布结构和模式，所计算的Moran 指数及其z得分如图9 所示。由图9 可看出山东省2000-2020 年的Moran 指数范围在0.15～0.45、z得分在10～40，这表明山东省内植被空间分布呈现聚集模式分布。其中2000、2008 聚集度相对较高，2018、2019、2020 年聚集度相对较低。该结果与图5 综合分析表明，山东省的植被覆盖度呈现全局增长的趋势，而非某个地区的特定贡献。

图9 2000-2020 年山东省Moran 指数与z 得分Fig.9 Moran index and z-score in Shandong Province from 2000 to 2020

3.6 植被覆盖度的稳定性与持续性

根据实验中CV 指数的分布规律，依据表5 对其划分了五个等级。

表5 研究区CV 指数分级表Table 5 CV Index classification table of the study area

根据上表所计算得到的CV 指数等级空间分布图如图10 所示。

研究区FVC 的各稳定性面积占总面积比例分别为高稳定性（33.5%）、中高稳定性（43.6%）、中度稳定性（15.5%）、不稳定性（5.7%）、极不稳定性（1.7%）。从空间分布看，高稳定性主要分布在山东省西部；中高稳定性主要位于山东省南部；中度稳定性主要位于东营市；不稳定性与极不稳定性主要位于济南市、淄博市、潍坊市、青岛市与烟台市以上城市的市中心及周边地区。通过图10 可以看出，变异系数变化与距离城市中心距离长度呈现负相关，本文推测这是由于人类活动给FVC 的稳定性带来相应的影响。

图10 2000～2020 年山东省CV 指数等级分布图Fig.10 Distribution map of CV index in Shandong Province from 2000 to 2020

本文将Hurst 指数分为0～0.4，0.4～0.6，0.6～1 这三个区间，所得到的Hurst 指数等级分布如图11 所示。

图11 山东地区植被覆盖度Hurst 指数Fig.11 Hurst index of FVC in Shandong Province

由图11 可以看出，山东省大部分区域FVC 变化的Hurst 指数值在0.6～1 之间，占总研究区面积的85.4%，几乎覆盖整个山东省，这表明大部分地区未来变化趋势与过去变化趋势一致、区域持续性较明显；Hurst 指数值在0.4～0.6 之间的面积占比为14.3%，主要分布在各市的市中心及其周边地区，这部分地区与过去变化趋势不存在较强的相关性，持续性与反持续性不明显；Hurst 指数值在0～0.4 之间的面积占比为0.3%，主要包括济南市市中区、历下区，昌乐县宝城街道，济宁市市区、新兖镇、南屯煤矿，临沂市罗庄区，区域反持续性明显。Hurst 指数值的大小与城市区域距离呈现正相关，表明人类活动会影像植被覆盖度的变化，这与CV 指数的分析结果一致。同时，仍然存在许多Hurst 指数较大的城市，这说明山东省部分城市在发展过程中兼顾了生态环境的保护。

4 结论

本文对山东省2000-2020 年间的植被覆盖度进行了计算，并基于土地覆盖数据集进行了精度评定，在精度满足要求的基础上，利用最小二乘法、变异系数、Hurst 指数等对其时空变化特征进行了分析，得出以下结论：

（1）山东省2000-2020 年间总体FVC 变化趋势呈上升趋势，但内部每个市的FVC 变化趋势有所不同；

（2）2000-2020 年间，山东省植被覆盖度变化范围为60.3%～71.5%，因此山东省FVC 属于中高覆盖度等级。其中中部地区与沿海的低FVC 占比较多，并存在低FVC 区域由中部向南部扩展的趋势。最终的低植被覆盖区域呈现零星分布状态；

（3）2000-2020 年间，山东省FVC 的变化趋势呈现基本不变或轻微改善的特征，其中基本不变面积占比约为76%，轻微改善面积占比约为15%。呈现向好趋势发展；

（4）2000-2020 年间，山东省FVC 分布的Moran 指数由0.45 降至0.15，表明其分布模式由聚集性分布向随机分布靠近，并结合FVC 的时序结果得到结论：山东省的植被覆盖度呈现全局增长的趋势，而非某个地区的特定贡献；

（5）2000-2020 年间，山东省FVC 变化稳定性以高稳定性和中高稳定性为主，占比分别为33.5%、43.6%。山东省大部分地区FVC 变化趋势有显著的持续性；各市周边地区呈现较多的随机性，表明其难以预测；各市中心等区域的FVC 变化趋势存在显著的反持续性。

本文对山东省2000-2020 年间的植被覆盖度进行了详细分析，并通过诸多时空分析方法对其内部发展趋势与分布模式进行评价。分析结果表明山东省大部分区域植被覆盖度无论是整体还是局部均在稳定增长，仅少部分城区存在着城市发展与生态平衡之间的矛盾。综上所述，通过GEE 平台，能够高效率、低成本的实现大范围、长时间周期的FVC 的提取与时空监测分析。