监测植被覆盖变化是评价景观条件和区域生态环境质量的有效途径, 也是气候变化和地表碳循环等诸多领域研究的基础^[1]。植被覆盖度定义为单位面积内植被冠层垂直投影所占的百分比^[2]。它指示了植被的茂密程度和植物进行光合作用面积的大小, 是量化植被群落生长态势和人类对自然胁迫效应的综合指标, 在全球变化与陆地生态系统响应(GCTE)和国际地圈生物圈计划(IGBP)等研究中具有重要地位^[3-4]。天地一体化的遥感监测和高光谱、高分辨率影像应用的日益成熟, 使得植被覆盖变化信息的获取更加准确、灵敏^[5]。NDVI(归一化植被指数)作为当前主流的植被状况反映指标, 在尺度、时效以及表征植被的生物物理特征等方面具有明显优势^[6]。国内外学者利用NDVI监测植被覆盖变化已取得丰硕成果：Tucker等^[7]研究1981—1999年NDVI序列发现, 由于气候变暖的原因, 北半球中纬度植被覆盖增长趋势明显。Park等^[8]分析东亚北部1982—2006年植被生长季NDVI发现以1990s中期为分界点, 植被活动由强变为弱。

《京津冀协同发展规划纲要》指出生态环境保护作为京津冀协同发展三个率先突破之一, 以“环”(三大区域构建国家公园环)、“带”(京津冀湿地群将实现连接)、“廊”(绿色生态廊道互通京津冀)、“养”(修复西北部生态涵养区)等方式构建京津冀城市群生态体系。晏利斌等^[9]在京津冀地区, 发现植被覆盖经历衰退期(1999—2002年)后, 在2003—2006年间表现出抬头趋势。2006年后京津冀地区植被覆盖能否延续这一趋势, 未来又将保持怎样的发展趋势, 需做深入的研究。已有研究多是针对京津冀区域NDVI均值的时序分析, 其缺陷是不能详实地描述植被覆盖时空格局的动态变化, 并且NDVI影像的异常值会引起变化趋势分析的偏差^[10], 同时相关的研究缺少对植被覆盖变化趋势的预测。鉴于上述情况, 本文基于MODIS NDVI遥感数据, 采用线性回归分析、稳定性分析、R/S分析等数理方法, 反演2005—2015年京津冀地区植被覆盖时空演变趋势, 评价了区域植被覆盖的稳定性, 并在此基础上对植被未来变化趋势进行评估和预测, 为京津冀城市群地区形成生态完整性网络、构建区域生态安全格局提供科学依据。

1 数据与方法 1.1 研究区概况及数据来源

京津冀城市群地区包括北京、天津和河北省, 位于我国华北地区(36°03′—42°40′N, 113°27′—119°50′E), 东临渤海, 西倚太行, 北接燕山, 南面华北平原。包括北京市、天津市以及河北省11个地级市, 土地面积21.8×10⁴ km², 常驻人口1.1×10⁸人, 是北方经济规模最大、最具有活力的地区；境内地貌复杂, 高原、山地、丘陵、盆地、平原等类型齐全, 地势自西北向东南呈阶梯状倾斜；气候属于典型大陆性季风气候, 寒暑悬殊, 雨量集中。孟丹等^[11]研究该地区植被覆盖变化驱动因素, 得出结论, 区域内非气候因子驱动占有绝对优势, 表明人类活动对植被变化贡献巨大。截至目前, 该地区共建设18个国家级自然保护区, 且北部地区从属三北防护林华北段, 生态安全战略性意义重大。

遥感数据来源于国际科学数据镜像网站(http://www.gscloud.cn)的MODND1M数据产品, 该产品是将空间分辨率为500 m的MODIS(Terra星)NDVI每天数据, 采用国际通用的最大合成法(MVC)排除云和大气等干扰得到月NDVI产品, 最后通过京津冀地区行政区掩膜, 裁剪出2005—2015年逐月NDVI的栅格影像。其中, 2015年9月份的遥感影像质量不佳, 以近3年同月影像均值代替。为了避免雪盖等因素的干扰, 更加准确反映植被覆盖状况, 在充分了解当地自然和物候特征的基础上, 选取植被生长季(4—10月份)作为研究时段。

1.2 植被覆盖度

植被覆盖度和NDVI之间存在极显著的线性相关关系, 建立二者之间的转换关系, 可接提取植被覆盖度信息^[12]。像元二分法模型^[13-14]对影像的辐射校正影响不敏感, 且计算简便, 是基于像元线性分解计算植被覆盖度的一种有效方法。它的原理是将遥感传感器所获取的地表光谱信息分解成两部分, 完全由绿色植被所贡献的信息NDVI_v和无植被覆盖裸地所提供的信息NDVI_s。像元的NDVI信息则由这两部分合成得到, 即：

(1)

单个像元中, 假设植被覆盖面积比例(植被覆盖度)为C_i, 则裸地覆盖面积比例为1-C_i

再设一个全部由植被覆盖的纯像元的信息为S_veg, 则对一个由植被与裸地组成的混合像元, 植被所贡献的信息NDVI_v就等于S_veg与C_i的乘积, 即

(2)

同理设一个全部由裸地覆盖的纯像元信息为S_soil, 则对由植被与裸地组成的混合像元, 裸地所贡献的信息NDVI_s就等于S_soil与1-C_i的乘积, 即

(3)

将(2), (3)带入(1)式变换则得到C_i的计算公式：

(4)

S_veg和S_soil的取值是像元二分模型应用的关键。由于大气云层、地表湿度和光照等条件的影响, S_soil不再是一个接近于0的定值, 其变化范围通常在-0.1—0.2之间。对于纯植被像元来说, 植被类型及其构成、植被的空间分布和植被生长的季相变化都会造成S_veg值的时空变异。不同的研究对S_veg和S_soil的取值方法有较大的差异, 本文采用一种广泛应用的方法, 根据整幅影像上NDVI的灰度分布, 以0.5%置信度截取NDVI的上下限阈值分别近似代表S_veg和S_soil。

1.3 植被覆盖度年际变化趋势

基于像元的一元线性回归分析可以模拟区域时空格局变化趋势^[15]。本文以年植被覆盖度和时间序列建立一元线性方程, 模拟2005—2015各年份植被生长季覆盖度变化趋势。计算公式如下：

(5)

式中, θ_slope是回归方程的斜率, 若值为正, 表示植被覆盖度有增加趋势；反之, 植被指数有下降趋势。n为监测年数；C_i表示第i年的植被覆盖度。趋势的显著性采用F检验, 显著性代表变化趋势的可信程度。统计量计算公式为：

(6)

式中, 称为误差平方和, 称为回归平方和, y_i为第i年覆盖度的像元值, ŷ_i为其回归值, y为监测时段覆盖度的平均值, n为监测年数。根据检验结果将变化趋势分为6个等级：极显著减少(θ_slope＜0, P＜0.01)；显著减少(θ_slope＜0, 0.01＜P＜0.05)；不显著改善(θ_slope>0, P>0.05)；不显著退化(θ_slope＜0, P>0.05)显著增加(θ_slope>0, 0.01＜P＜0.05)；极显著增加(θ_slope>0, P＜0.01)。

1.4 植被覆盖度稳定性

变异系数是量化一组数据离散程度最常用的一种方法, 值的大小反映了数据序列的波动幅度。本文基于单个像元在研究时段内植被覆盖度的变异系数, 模拟京津冀地区植被覆盖度的稳定性。计算公式如下：

(7)

式中, σ为标准差, 计算公式为, C为研究时段多年平均植被覆盖度。采用聚类分析Natural Breaks(Jenks)将变异系数(C)分为低(C≤0.045)、较低(0.045＜C≤0.063)、中(0.063＜C≤0.083)、较高(0.083＜C≤0.133)、高(C≥0.133)5个等级。

1.5 植被覆盖度演变趋势预测

R/S分析法(重新标度极差分析法)由Hurst在研究水文观测资料时提出, 是定量描述时间序列信息依赖性的一种分型理论, 目前在经济学、水文学、地质学气候学等领域有着广泛应用^[16-17]。基本原理如下：

对于一个时间序列{ξ(t)}, t=1, 2…, 等于任意正整数τ≥1, 定义均值序列：

(8)

(9)

(10)

(11)

引入无量纲的比值R/S, 若存在H使得R/S=(cτ)^H成立, 则说明{ξ(t)}存在Hurst现象, H为Hurst指数。在双对数坐标系中lnτ, lnR/S使用最小二乘法拟合, 求得像元Hurst指数。Hurst指数取值范围：

(1) 若0＜H＜0.5, 表明植被覆盖度时间序列具有反持续性, 过去变量与未来趋势呈负相关, 序列有突变跳跃特性, H值越接近于0, 反持续性越强；

(2) 若H=0.5, 表明植被覆盖度时间序列为互相独立的随机序列；

(3) 若0.5＜H＜1, 表明植被覆盖度时间序列具有长期相关特性, 过程具有持续性, H值越接近于1, 持续性越强。

2 结果与分析 2.1 植被覆盖度时空变化特征

2005—2015年京津冀地区植被覆盖度呈上升趋势(0.065/10年), 植被覆盖度最低出现在2006年(0.63), 最高出现2014年(0.71), 平均值为0.672。增长过程出现明显的波动(图 1), 可大致分为三个阶段：第一次飞跃期(2006—2008)、第二次飞跃期(2009—2011)和平稳波动阶段(2011—2015), 其中两次飞跃期增长速率均为0.05。在研究时段内, 2009年出现较为明显的下降是由于该地区遭受接近历史极限的干旱灾害^[18], 对植被造成了严重的影响。经历两次飞跃期后, 植被基数已较高, 2011年后保持良好的平稳波动姿态。

2005—2015年京津冀地区平均植被覆盖度区域差异性比较明显(图 2)。北部滦河流域水热条件充足, 植被覆盖度整体情况较好, 植被覆盖度较小的区域除了分布在张家口市西北和环渤海滩涂地区外, 其余分布于各大城市建成区及其周边区域, 如北京、天津、石家庄、唐山等, 说明了城市的发展对植被覆盖造成了较大的影响。以市行政区为单位, 分区统计植被覆情况, 结果表明：衡水以0.75的覆盖度居于首位, 张家口最低为0.55；北京、石家庄、保定等发达城市植被覆盖度均在0.7以上, 而环渤海城市植被覆盖度相对较低。

2.2 植被覆盖度变化趋势和稳定性评价

随着“三北防护林”、“退耕还林还草”、“绿色河北”等政策不断的实施, 京津冀地区植被有较明显的改善(图 3, 表 1)。空间格局演变特征主要表现为：(1)极显著改善区域比重最大(47.45%), 以“滦河流域—大马群山—太行山脉”东西连横, 南北相接, 环绕于京津冀北部和西部；东南部以衡水、邢台东部为轴形成一条明显绿带；其余分布于北京、石家庄、天津城市中心区。(2)不显著变化区域占全区面积的三分之一(33.9%), 主要以华北平原为中心展布, 该区域是历史悠久的农耕区, 不显著性符合农作物耕种特征。(3)退化区域比重为6.8%, 零星散布于各大城市中心区周边, 佐证了城市规模扩大对植被覆盖的显著影响。

京津冀地区2005—2015年稳定性整体表现为“高低波动并存, 低波动居多, 地域性明显”。空间格局(图 4)表现为：(1)低波动区(深绿色)和较低波动区(绿色)占全区面积的66.3%, 植被相对稳定, 基本吻合不显著改善区。(2)中度波动区(白色)和较高波动区(黄色)呈斑块状分布：较高波动区所占比重较小为7.4%, 主要分布在北京中心区、张家口南部、沧州东北一带, 植被表现为极显著改善；中度波动区比重为26%, 环簇于较高波动区边缘向外延伸。(3)高波动区(红色)比重仅占0.3%, 环渤海滩涂零星散布。

北部滦河流域以及西北部大马群山, 在全球变暖的大背景下, 加之20世纪80年代“三北防护林”工程建设的持续推进, 使得植被覆盖度显著增加, 人为护林养林使得大幅度波动变化可能性较小, 稳定性较高；衡水和邢台东部属于河北平原海河水系冲击平原中部, 在构造上属于冀中凹陷中部^[19], 年均气温明显高于周边^[20], 植被改善显著性优于周边区域, 形成一条明显的绿带；东南部华北平原是我国主要的农业耕种区, 农作物是该区域的主要绿色来源, 所以覆盖度变化不显著；北京、石家庄等城市建成区植被覆盖度有较显著的改善, 波动幅度较高, 说明近年来城市绿化对自然环境的改造取得了一定功效。

2.3 植被覆盖度演变趋势预测

Hurst指数(图 5)均值为0.55, 持续性序列比重占63.4%, 反持续性序列比重占36.6%。参考李双双等^[21]对同纬度陕甘宁地区植被覆盖变化的相关研究, 将京津冀地区Hurst指数进行等级划分, 弱持续性序列(0.35—0.65)占研究区总面积的72.4%, 植被的恢复若依靠单一因素(自然或人为干涉修复, 例如自然维持的原始森林、人类维持的农田), 其恢复序列表现出较强的持续性, 而京津冀地区的这种弱持续性证明了植被覆盖变化是在自然、人为等多种因素共同驱动下形成的。从空间分布来看(图 6), 持续性序列(绿色)占主体, 与反持续性序列(紫色)呈斑块状分布, 加之2005—2015年植被改善区面积高达93.2%, 说明该地区植被恢复、生态环境持续改善已成为主旋律, 进而预测京津冀大部分地区植被覆盖在未来将持续改善。强反持续性序列主要分布在张家口, 该区域属于坝上高原, 水热条件较差, 虽通过长年的环境建设使植被覆盖度有了较大的提高, 但相较于其他地区仍属于低覆盖区, 未来变化趋势还需要密切关注。

为了进一步了解京津冀地区植被覆盖变化趋势的可持续性, 本文将趋势分析和R/S分析结果重分类后, 进行叠加耦合分析, 得到如下结果(图 5, 表 2)：未来植被覆盖持续改善区域比重超过一半以上(58.8%), 反持续性改善比重为34.4%, 持续退化比重为4.8%, 整体情况较为乐观。从空间分布来看, 反持续改善(蓝色)主要分布在张家口、沧州以及保定东南地区；由于城市化、人口发展、经济结构调整等因素, 持续退化(红色)主要分布在天津、廊坊、沧州。以此为基础, 推断在2016—2020年期间, 如果没有较大的气候波动, 伴随着生态体系的建设, 京津冀地区将迎来第三次植被恢复飞跃期, 但增长幅度相较于前两次略低且逐步趋于平稳。

3 结论与讨论

本文基于MODIS NDVI遥感数据, 采用线性回归分析、稳定性分析、R/S分析等数理方法, 反演了2005—2015年京津冀地区植被覆盖时空演变趋势, 评价区域植被覆盖的稳定性, 并在此基础上对植被未来变化趋势进行评估和预测, 得到以下结论：

(1) 近11年来京津冀地区植被覆盖度呈上升趋势(0.065/10年), 增长过程分为三个阶段：第一次飞跃期(2006—2008)、第二次飞跃期(2009—2011)和平稳波动阶段(2011—2015)。北部滦河流域水热条件充足植被覆盖度较好, 受地带性规律影响的张家口地区植被覆盖度整体较低, 环渤海滩涂最差, 城市发展对植被覆盖影响显著。

(2) 2005—2015年京津冀地区植被覆盖极显著改善区域比重最大(47.45%), 不显著变化区域占全区面积的三分之一(33.9%), 退化区域比重为6.8%, 零星散布于各大城市中心区周边。稳定性整体表现为“高低波动并存, 低波动居多, 地域性明显”。低波动区和较低波动区占全区面积的66.3%；较高波动区所占比重较小为7.4%, 主要分布在北京市中心区、张家口南部、沧州东北一带；中度波动区比重为26%, 环绕于较高波动区边缘向外延伸；高波动区比重仅占0.3%, 环渤海滩涂零星散布。

(3) 京津冀地区植被以持续恢复为主, 持续性序列比重占63.4%, 反持续性序列比重占36.6%。植被覆盖持续改善区域比重超过一半以上(58.8%), 反持续性改善比重为34.4%, 持续退化比重为4.8%。反持续改善主要分布在张家口、沧州以及保定东南地区；持续退化主要分布在天津、廊坊、沧州。

京津冀地区植被持续恢复是本文主要结论之一, 本文虽然使用了最新的遥感影像, 但对于预测变化趋势来说, 11年的时间序列相对较短, 变化趋势的持续性是否能被充分反映有待于探讨和检验。再者, 影响植被覆盖变化的因素主要分为自然因素和人为因素^[22], 尤其是人为因素突出了变化趋势的不确定性, 人类活动在植被覆盖时空演化过程中相对贡献的量化评定还有待于进一步的探讨。近年来, 京津冀城市群地区在气候和城市扩张的双重压迫下, 水土资源短缺和环境污染已成为该地区典型生态问题, 植被恢复相应得到了较多的关注。气候的变化将会增加陆地生态系统发生大范围、严重以及不可逆转影响的概率^[23], 虽然京津冀地区经过多年的人工植被恢复减缓了气候变化的副作用, 但人工生态恢复也会存在不适应性和波动性的缺陷。未来, 以生态学视角综合分析植被恢复对气候变化和人类活动的响应, 量化其相对贡献, 解释相互作用机制, 在研究生态环境恢复的可持续性方面将具有更重要意义。