植被覆盖度( fractiona of vegetation coverage, FVC)通常定义为植被(包括叶、茎、枝)在地面的垂直投影面积占统计区总面积的百分比^[1]。FVC不仅是描述地表植被覆盖的重要参数, 也是反映地表植被长势的重要参量, 同时还是草原退化监测最重要、最敏感的指标。我国是世界上草原资源最丰富的国家之一, 草原总面积近4×10⁸ hm², 占全国土地总面积的 40%。很多研究表明, 近年来我国草原大面积出现明显的退化。呼伦贝尔草原是我国保存较为完好的草原之一, 近年来也出现不同程度的退化。草原的退化往往表现为植被覆盖度的降低^[2,3,4,5], 基于多种平台搭载传感器的遥感观测技术很容易对大面积的草原进行实时、同步监测。

植被覆盖度的遥感监测方法主要有物理模型法和统计模型法。物理模型法是基于辐射传输内部物理机制, 过程复杂, 变量多, 难以在实际中应用。统计模型法中较常用的是回归模型法、植被指数法及像元分解法。回归模型法精度高, 但仅适用于小面积区域, 有一定的局限性^[3]。植被指数法直接利用植被指数, 近似地估算植被覆盖度, 优点是不需要建立回归模型, 缺点是精度相对较低。像元分解法利用混合像元分解的方法反演植被覆盖度。陈晋等^[6]和陈云浩等^[7]根据 TM 像元为非均一混合像元的特点, 提出基于土地覆盖分类的植被覆盖度估算亚像元模型, 并对北京市海淀区的植被覆盖度进行估算, 为大面积植被覆盖度的估算提供了一种有效途径。在像元分解方法中, 最具代表性的是像元二分模型。该模型基于线性的像元分解方法, 形式相对简单, 在很大程度上削弱了大气、土壤背景和植被类型的影响, 普适性较好^[7,8,9]。

对草原植被覆盖度进行遥感监测, 尤其是长期动态观测, 如果仅仅从某一时空来探讨草原植被与气候因子的内在关系, 则容易忽略时间序列包含的一些重要信息。国内一些学者利用 NOAA/AVHRR数据, 结合地面气象数据, 研究近20年来中国东北地区呼伦贝尔草原归一化植被覆盖指数(normalized difference vegetation index, NDVI)与温度、降水的相互关系。呼伦贝尔草原对气候变化的响应特征表明, NDVI与气温(尤其在春季)相关性显著^{[10,11,12,13,14,15]}。国外一些学者探讨植被信息与气象因子的内在关系, 论证了气候对植被的影响^{[10,11,12,13,14,15]}。上述研究大多是对某一时段内植被的NDVI数据进行间断性动态分析, 以此间接分析植被的生长状况。然而, NDVI对低覆盖植被的变化较敏感, 且饱和值较低^[8]。如果可以直接对草原植被覆盖度时间序列进行分析, 就可以在一定程度上规避这些不足, 更加清晰地反映草原退化的程度与过程。因此, 本文基于2000—2014 年呼伦贝尔草原核心区的新巴尔虎右旗、新巴尔虎左旗、陈巴尔虎旗和鄂温克族自治旗(简称“牧业四旗”)的 MODIS 产品 MOD13A2, 利用像元二分模型, 反演得到 15 年草原植被覆盖度的时间序列, 然后结合研究区的气象数据, 采用像元趋势分析方法和植被覆盖度异常变化检测算法, 研究15年间植被覆盖度的变化过程, 以期对研究区内草原的保护与发展提供依据。

1 数据与方法

1.1 研究区概况

呼伦贝尔草原位于内蒙古自治区呼伦贝尔市西南部, 地理坐标为: 115°31′—121°09′E, 46°10′—50°12′N, 西部与蒙古共和国接壤。呼伦贝尔草原地形总体上西高东低, 地势由西到东缓慢过渡。中部是波状起伏的海拉尔台地高平原, 其中包含大沙带和零星沙丘堆积。呼伦贝尔草原具有典型的温带大陆性季风气候, 光照充足, 雨热同季。土壤类型多为暗棕壤、黑土和暗色草甸土, 有机质丰富, 土层深厚。年均温为 3.1℃, 15 年间平均年降雨量为 275.5 mm, 夏季温暖多雨, 冬季严寒干燥, 导致草原生长周期短, 生长旺季集中在温暖多雨的夏季。本文选择天然牧场所在的呼伦贝尔草原核心区域“牧业四旗”为研究区(图1)。

显示原图| 下载原图ZIP| 生成PPT

图1   研究区边界示意图

Fig. 1   Geographic outline of study area

1.2 数据来源及其预处理

1.2.1 遥感数据与预处理

本研究采用的影像数据为研究区 2000—2014 年基于最大合成法(Maximum Value Compositing, MVC)的 MOD13A2 产品(https://ladsweb.nascom. nasa.gov/), 时间分辨率为 16天, 空间分辨率为 1 km。其中, 2000 年可用数据开始于 2 月。MOD13A2 产品已经过辐射校正和几何粗校正、每日每轨图像几何精校正以及除坏线等处理, 包括 4 个轨道数据: h25v03, h25v04, h26v03 和 h26v04。同时, 采用与MOD13A2 产品配套的研究区 2000—2014 年土地利用分类数据产品 MCD12Q1, 时间分辨率为 1 年, 空间分辨率为 1 km。运用MODIS Reprojection Tools (MRT)软件对下载的数据产品进行数据格式和投影转换, 完成影像数据的拼接。采用ERDAS 9.2软件, 对拼接后的影像进行裁剪等处理, 裁剪时采用的底图为覆盖该研究区的边界矢量图。

1.2.2 气象数据

本文使用研究区内各气象站点 2000—2014 年旬平均气温、最高气温、最低气温以及旬平均降水量。主要的气象站点包括海拉尔站点, 新巴尔虎左旗站点, 新巴尔虎右旗站点, 鄂温克族自治旗站点, 陈巴尔虎旗站点。数据来源于中国气象数据服务 平台。

1.3 研究方法

1.3.1 植被覆盖度遥感反演方法

首先, 采用像元二分模型^[6]反演草原植被覆盖度, 公式为

(1)

其中, NDVI_veg为全植被覆盖像元NDVI值, NDVI_soil为完全裸土的像元 NDVI 值。像元二分模型中, 两个必需参数的取值目前主要依靠对NDVI的数据统计。基于该研究区的土地利用数据, 计算每个单元NDVI 数据的频率累积值。经过反复实验对比, 土壤单元内取累积频率为 5%的 NDVI 作为 NDVI_soil, 土地利用单元内取累积频率为 95%的NDVI 值作为 NDVI_veg。最后, 依据式(1), 得到研究区 2000—2014年共计343幅植被覆盖度分布图。

利用MOD15 LAI产品和平均孔隙率模型, 对像元二分模型进行初步间接验证。

对于连续、均匀的冠层, 假定其叶片随机分布, 那么当一束光线穿过冠层时, 可以用 Poisson 模型计算孔隙率^[8]。这里的孔隙率指像元内土壤所占(像元面积)的百分比(假设该混合像元主要由植被和土壤两种组分随机构成)。孔隙率计算公式如下:

(2)

G_i函数表示的物理意义是: 在高度 Z 处, 叶面积体密度取 1 时, 向光子入射方向垂直面的平均投影值(假定叶子球形分布, 其取值为 0.5); μ_i 为太阳天顶角(可以通过查阅对应的元数据文档获得)的余弦值; LAI 为植被的叶面积指数。对于水平方向均匀、垂直方向分层的连续植被, 视线或射线方向上的孔隙率只与天顶角有关^[8]。

最终, 采用如下公式反演植被覆盖度:

(3)

为了验证像元二分模型的准确性和精度, 下载2014年7月覆盖研究区的MOD15 (LAI)数据产品, 随机选取对应的特征草地覆盖区域, 与像元二分模型得到的结果进行对比验证, 得到两者间的植被覆盖度散点图(图 2)。可以看出, 两种模型得到的结果相关性较强(R² 为 0.9847)。因此, 应用像元二分模型反演得到的植被覆盖度, 具有一定的准确性和可靠性。

显示原图| 下载原图ZIP| 生成PPT

图2   两种植被覆盖度反演结果的散点图

Fig. 2   Scattering points of two results estimating FVC

1.3.2 植被覆盖度的趋势变化

采用以下公式^[22]对连续多年观测的草原最高植被覆盖度进行动态趋势变化分析:

(4)

其中, n为监测的年数, C_i为第i年的草原植被覆盖度。利用式(4), 可以逐像元有效地求取多幅覆盖度影像数据中每个(特定)像元的灰度值(表征覆盖度值)在多幅连续影像时间序列数据中的变化趋势。θ_slope 表征草原植被覆盖度随时间的变化程度, 其值可正可负: 正值表示草原植被覆盖度随时间呈现增加的趋势, 其值越大, 表征增加趋势越明显; 负值表示草原植被覆盖度随时间呈现减少的趋势, 其绝对值越大, 表征减少趋势越显著。

1.3.3 草原植被覆盖度时间序列

将历年草原植被覆盖度结果与土地利用覆盖数据相结合, 选取草地类别的特征像元, 分区域对草原植被覆盖度进行时间序列提取。对提取出的原始时间序列进行 S-G 滤波处理, 得到不同区域草原植被覆盖度的时间序列。将提取出的各区域历年平均草原植被覆盖度时间序列进行平均, 得到整个研究区历年平均草原植被覆盖度时间序列。利用时间序列趋势项提取算法, 进行时间序列的趋势项提取。

1) S-G滤波原理。

Savitzky 等^[23]提出 S-G 滤波器, 又称最小二乘法或数据平滑多项式滤波器。S-G 滤波的设计思想是, 找到合适的滤波系数 P_i , 以保护高阶距来实现滑动窗内的最小二乘拟合。基本原理是: 通过取点X_i附近固定个数的点, 拟合一个多项式, 多项式在X_i的值给出它的光滑数值 g_i。基于 S-G 滤波原理, 覆盖度时间序列数据的S-G滤波公式如下:

(5)

式中, 为合成序列数据, Y_j+i 表示原始序列数据, P_i 为滤波系数, N 为滑动窗口所包括的数据点数量(2m+1)。S-G 滤波可以有效地剔除时序中的背景噪声, 保留原始时间序列中真实的植被生长信息。因此, 为了最大程度地保证趋势项提取算法和时序异常检测算法的精度, 我们对原始时间序列做 S-G 滤波处理。

2) 时间序列的趋势项提取算法。

趋势项主要反映时间序列的长期变化。Assi-makopoulos^[24]提出 Theta 算法, 通过 Theta 变换, 可将原始时间序列转换为反映长期趋势的新时间序列, 如式(6)所示:

(6)

其中, x_t为输入时间序列, y_t为输出时间序列, 为时间序列在t时刻的瞬时趋势, , N为样本空间。

Theta 算法每进行一次, 都会选出一个符合条件的序列值, 用其相邻序列值的平均值代替, 直到 为止。需要特别指出的是, 草原植被 覆盖度 15 年时间序列的长期变化不是线性的, 而是不规则曲线。经过反复试验, 我们选择 时停止运算, 此时的新序列即草原植被覆盖度时间序列的趋势项。

草原植被覆盖度受诸多自然和人为因素的影响, 在不同的时间段覆盖度会出现不同程度的剧烈 变化。因此, 本文采用基于最小二乘拟合的植被覆盖度异常变化检测算法, 对草原植被覆盖度时间序列进行检测, 提取时间序列中的覆盖度急剧变化点, 并在此基础上做进一步分析。

3) 基于最小二乘拟合的植被时序异常变化检测算法^[25]。

对已经得到的时间序列数据, 采用基于时间步长的后向检测方式。对每一个时间点 i (2≤ i ≤N-1), 遵循下列判定原则:

其中, ΔFVC_i代表i时刻的草原植被覆盖度。

基于式(7)~(9), 对每一个时刻点 i, 建立如下峰/谷点判定函数, 并进行峰谷点检测:

时间序列中所有的峰/谷点检测完毕后, 为了寻找时间序列中的转折点, 基于贝叶斯信息准则(Ba-yesian Information Criterion, BIC), 采用如下公式检测并提取覆盖度异常变化点:

最后, 将所有可能的覆盖度异常变化点一一检测出来, 并与影响草原植被覆盖度的气象因素进行相关性分析。

1.3.4 相关性分析

首先计算相关系数, 然后计算偏相关系数。植被覆盖度与温度或降雨量之间相关系数的计算公式如下:

其中, R_xy为两个变量的相关系数, x_i为第i年/月的植被覆盖度, y_i为第 i 年/月的温度或降雨量, 为多年/月植被覆盖度的平均值,为多年/月温度或降水的平均值, n 为样本数。采用偏相关分析方法, 控制其中一个变量, 分析其余两个量的相关性。基于降雨量的植被覆盖度与温度的偏相关系数以及基于温度的植被覆盖度与降雨量的偏相关系数计算公式如下:

其中, R₁₂₃为固定变量 3 后变量 1 和变量 2 之间的偏相关系数, R₁₂, R₁₃, R₂₃ 分别表示变量 1 与变量2、变量 1 与变量 3、变量 2 与变量 3 的相关系数。根据已有的植被覆盖度、气温和降雨量数据, 分析平均植被覆盖度与温度和降水的相关性。

显示原图| 下载原图ZIP| 生成PPT

图3   2014年研究区域草原植被覆盖度分布

Fig. 3   Spatial distribution of FVC of study area in 2014

2 结果与分析

2.1 草原植被覆盖度空间分布现状

本文对植被覆盖度的级别做如下定义: 0~20%为低覆盖度, 20%~40%为较低覆盖度, 40%~65%为较高覆盖度, 65%以上为高覆盖度。时间序列分析表明, 每年7—8月为草原的生长旺季, 植被覆盖度最高, 故求取 7—8 月平均值作为年平均覆盖度值。首先, 得出2014 年研究区草原平均植被覆盖度的空间分布状况(图 3)。然后, 对 2014 年研究区内不同等级植被覆盖度占草原总面积的比例进行统计(表 1)。可以看到, 研究区内以中等覆盖度为主, 其次是低覆盖度, 高覆盖度占的比例最小。低覆盖度和较低覆盖度主要分布在研究区中西部, 中等覆盖度和高覆盖度集中分布在研究区东部, 总体上呈现东高西低的格局, 特别是研究区西南部草原植被覆盖度相对于中东部偏低, 反映出西南部草原退化较严重。

2.2 15年来植被覆盖度变化趋势

研究区每年 7—8 月雨量充沛, 气温适宜, 草原植被生长最为旺盛; 因此, 7—8 月植被覆盖度的平均值代表当年植被覆盖的最高值。为了分析研究区15年来草原最高植被覆盖度的年际变化趋势, 逐像元求得每年 7—8 月的平均草原植被覆盖度, 得到2000—2014 年草原植被覆盖度变化趋势的空间分布(图 4)。可以看出, 研究区内草原植被覆盖度变化在空间分布上呈现东高西低的格局; 覆盖度降低较快的区域主要分布在研究区西部, 占总面积的18.71%; 覆盖度增加较快的区域主要分布在研究区的东部, 占总面积的 15.95%, 特别是东南部和东北部, 草原植被覆盖度年际增速较为明显, 反映该区域草原植被的生长状况良好。值得关注的是, 研究区西部草原植被覆盖度的减速较为显著, 表明 15年来草原退化趋势明显。

显示原图| 下载原图ZIP| 生成PPT

图4   2000—2014年植被覆盖度变化趋势的空间分布

Fig. 4   Distribution of the changing trend of FVC in 2000‒2014

进一步分析表明, 2000—2009 年(前 10 年)和2010—2014 年(近 5 年)两个时段, 研究区草原植被覆盖度的变化趋势呈现不同的规律(图 5)。

显示原图| 下载原图ZIP| 生成PPT

图5   15年平均植被覆盖度时间序列的趋势项提取结果

Fig. 5   Result of tendency term extracted from average FVC over 15 years

从图 6(a)可以看到, 前 10 年, 研究区东部草原植被覆盖度增加最显著; 中部和西部的个别区域也表现出显著增加的趋势, 但空间分布较为离散。研究区中西部, 特别是西部, 草原植被覆盖度降低的趋势最明显, 约占总面积的 17.15%。与 15 年草原植被覆盖度变化趋势分布(图 4)相比, 前 10 年覆盖度呈减少趋势的区域所占的比例明显高出 8.9%; 在研究区的中部, 前 10 年覆盖度降低的程度大于整个 15年时段。

从近 5 年草原植被覆盖度空间变化趋势分布(图 6(b))可以看到, 覆盖度增加较显著的区域主要分布在研究区的东北部、东南部以及中北部, 特别是东南部和东北部增加最为显著。草原植被覆盖度降低较明显的区域主要分布在研究区的西北部以及南部。与 15 年草原植被覆盖度变化趋势分布(图 4)相比, 近 5 年草原植被覆盖度呈减少趋势所占的面积高出 5.62%; 然而在研究区中部和东部, 植被覆盖度呈增加趋势的区域面积并无明显变化。

显示原图| 下载原图ZIP| 生成PPT

图6   2000—2009年(a)和2010—2014年(b)植被覆盖度变化趋势的空间分布

Fig. 6   Distribution of the changing trend of FVC in 2000‒2009 (a) and 2010-2014 (b)

对比前 10 年(图 6(a))和近 5 年(图 6(b))草原植被覆盖度变化趋势的空间分布, 可以看到, 在 15 年的时间跨度里, 近 5 年草原植被覆盖度呈增加趋势的区域集中在研究区中东部, 所占面积大于前 10年, 并呈现区域分块化特点; 近 5 年植被覆盖度降低较明显的区域集中在研究区西北部, 而在前 10年里, 覆盖度降低较显著的区域集中在中西部。同时发现, 呼伦湖以南及以东部分区域的植被覆盖度在前 10 年表现出减少趋势, 而在近 5 年里呈明显增加趋势; 植被覆盖度在前 10 年表现出轻度降低趋势的研究区南部, 在近 5 年里降低的幅度有所增加。进一步分析发现, 这两种截然不同的发展趋势并没有引起植被覆盖度总体水平的波动。

2.3 草原植被覆盖度时间序列

2.3.1 植被覆盖度时间序列总体规律

选取研究区内各个分区的草地特征像元, 分别提取其植被覆盖度, 得出相应的覆盖度时间序列。采用趋势项提取算法, 得到研究区 2000—2014 年平均草原植被覆盖度序列的趋势项(图 5)。结果表 明, 研究区 15 年间的草原植被覆盖度呈现先下降、后缓慢抬升的趋势。

图7 给出各分区的植被覆盖度时间序列。研究区东部植被覆盖度的峰值普遍高于西部, 且集中在7—8 月。年均草原植被覆盖度与年均气温和降水的平均相关系数分别为0.1362和0.6075 (表2), 其中植被覆盖度与降水呈正相关的面积占总面积的68.25%, 初步表明植被覆盖度变化主要与降水有关, 这与马玉玲等^[9]以及张戈丽等^[10]的研究结果基本上一致。

显示原图| 下载原图ZIP| 生成PPT

图7   研究区草原植被覆盖度时间序列

Fig. 7   Time-series of FVC of study area

2.3.2 草原植被覆盖度的异常变化

由于受其他外界因素的影响, 植被覆盖度在短期内可能会突然升高或降低。本文基于植被覆盖度时间序列, 引入植被时序变化检测算法^[25], 提取出覆盖度异常变化的转折点, 并对其分布规律做进一步分析。

显示原图| 下载原图ZIP| 生成PPT

图8   从时间序列中提取出的覆盖度异常变化点

Fig. 8   Abnormal changed points extracted from the FVC time-series

从图8和表3可以看出, 返青期内覆盖度的异常变化主要出现在 4—5 月, 与气温有明显相关性的占样本数的 83.33%; 生长旺季内覆盖度的异常变化主要出现在 7—8 月, 与降水量有明显相关性的占样本数的 66.67%; 枯萎期内覆盖度的异常变化主要发生在 9 月或 10 月, 与气温有明显相关性的占样本数的 50%, 与降水量有关的占样本数的25%。

3 结论与讨论

本文通过对 15 年的植被覆盖度进行分析, 得出呼伦贝尔草原植被覆盖度的空间分布特征: 从东到西逐渐降低, 草原西南部植被覆盖度最低。15 年平均植被覆盖度时间序列的趋势项提取结果显示, 呼伦贝尔草原植被覆盖度在年际尺度上呈现先降低、后缓慢升高的变化趋势。15 年来研究区草原植被覆盖度的年际变化状况分析结果表明: 植被覆盖度增加较明显的区域主要分布在草原东部, 植被覆盖度降低较明显的区域集中在草原西部。植被覆盖度异常变化点的检测结果显示: 返青期和枯萎期覆盖度的变化受温度的影响较大, 生长旺季覆盖度的剧烈变化主要与降水量有关。需要指出的是, 研究区西部草原退化明显, 加大植被保护力度对草原牧区的生态恢复和发展至关重要。

本文基于NDVI数据, 通过像元二分模型反演构建植被覆盖度时间序列, 并对 15 年来草原植被覆盖度时间序列进行时空变化分析和植被覆盖度异常变化点检测。本文初步在像元尺度上尝试应用草原植被覆盖度时间序列进行植被覆盖度异常点检测, 效果较好。另外, 时间序列趋势项的提取可以有效地检测出整个覆盖度时间序列中的时空变化转折点, 将原始时间序列按照转折点分成若干个不同时间尺度的子时间序列, 从而挖掘出更多的覆盖度时空变化信息。

本文也存在一些问题和不足: 1) 仅对植被覆盖度异常变化点与气象因素间的相关性做了分析, 初步探讨了植被覆盖度异常变化的时间分布规律及其影响因素, 下一步应对气候异常的影响进行更深入的分析; 2) 反演植被覆盖度采用的 MODIS 数据的时间分辨率为 16 天, 基本上反映出植被覆盖度年内和年际的主要变化特征, 但对更详细研究而言, 数据的时间分辨率偏低, 后续工作中需要改进; 3) 在分析影响草原植被覆盖度的因素时, 以气象因子为主, 暂时没有考虑对草原植被覆盖度影响机制比较复杂的人为因素和其他社会经济因素, 因而对2000—2009 年草原植被覆盖度总体上呈下降趋势, 2010—2014 年表现出上升趋势这一现象没能给出很好的解释。下一步研究将综合讨论气候、人为因素和社会经济因素对草原植被覆盖度的影响。

The authors have declared that no competing interests exist.

参考文献

文献选项

• 原文顺序
• 文献年度倒序
• 文中引用次数倒序
• 被引期刊影响因子