基于MODIS NDVI的长江中游区域植被动态及与气候因子的关系

植被可定义为陆地地表上所覆盖的所有植物群落的总称, 是陆地生态系统中的主体, 其生长受到自然和人为的共同影响, 对环境变化和气候变化十分敏感[1-2]。植被长期处在动态演变之中, 并在生物圈的物质循环、能量流动、生态系统稳定性上发挥着重要的作用, 气候变化是植被活动的重要影像因素, 研究植被的动态变化和气候因子的关系, 对区域生态平衡维护有重要的现实意义[3-4]。归一化植被指数(Normalized Difference Vegetation Index, NDVI)是反应和监测植被生长状况和动态变化的重要指示因子, 对地表植被的动态变化非常敏感(包括植被覆盖度、生物量和叶面积指数等), 在时间序列的植被监测方面有着广泛的运用, 已成为全球变化研究的重要方向之一[5-7]。

国内外学者利用NDVI数据, 从不同的时间和空间尺度上对各个区域的植被动态变化规律进行了研究, 探讨了各个不同研究区的NDVI变化特征、演变趋势和对气候变化的响应机制[8-10]。陶帅等通过分析2000—2015年长江上游植被变化特征, 指出该区NDVI未来可能出现波动的区域占一半以上, 海拔和年均温是影响NDVI变化的主要因素[11]。袁喆等探讨了2000—2015年长江流域NDVI与气象因子的关系, 结果表明NDVI呈逐年增加的趋势, 长江流域水量丰沛, 热量条件是影响植被生长的限制性因子[12]。Mao等和Stenhman等研究表明全球的植被变化受人为活动影响较大, 特别是一些大型的植树造林活动[13-14]。John和Ichii等分析全球范围内NDVI与气候变量之间的关系, 结果显示区域北部中高纬度地区, 植被NDVI与气温呈显著正相关, 而在北部和南部半干旱地区, 植被NDVI与降水呈显著负相关[15]。Fensholt等利用NDVI数据分析了半干旱地区植被的季节变化趋势和气候驱动因素, 表明全球平均半干旱地区的植被覆盖呈增加趋势, 降水是该地区植物生产的主要限制因素[16]。王桂钢等通过分析天山地区气候参数对草地季节变化影响的滞后性特征, 指出当旬降水和前一个月降水累积量对于草地长势的影响最大[17]。

不同区域尺度上的NDVI响应特征存在差别, 已有研究对长江流域植被变化进行了相关分析, 然而针对长江中游地区植被演变趋势和响应机制还未做充分探讨。本研究选取长江中游地区(湖北省、湖南省和江西省), 基于1999—2015年MODIS NDVI数据和气象资料, 分析长江中游地区NDVI的变化特征和演变趋势, 辨析了气象因子与NDVI在年际尺度和月际尺度上的关系, 从而全面揭示了长江中游地区植被动态变化规律, 探讨了植被对气候变化的响应机制, 以期为长江流域应对气候变化影响、区域生态环境建设与保护措施的制定, 提供科学参考依据。

长江中游区域包含湖北省、湖南省和江西省(24°25′N—33°16′N, 108°24′E—118°23′E), 总面积为56.46万km2(图 1)。区域气候类型属于亚热带季风气候, 降水量达到1000—1600 mm/a, 年平均气温在16—18℃[18]。长江中游区域地形主要以山地丘陵与平原为主, 海拔最高处位于鄂西山地一带, 最低处位于中部平原地区, 其中山地丘陵占50%以上, 地势西高东低, 平均海拔约为1497 m[19]。2018年地区总人口为17438万人, 占全国总人口12.7%, 生产总值为9.8万亿元, 城市化水平达50%以上, 人类经济活动频繁, 对区域扰动较大[20]。

长江中游地区林木资源丰富, 2015年湖北省、湖南省和江西省的森林覆盖率分别为39.61%、59.82%和64.00%[21]。主要植被资源有水杉(Metasequoia glyptostroboides)、冬青(Ilex chinensis Sims)、银杉(Cathaya argyrophylla)、珙桐(Davidia involucrata Baill.)、银杏(Ginkgo biloba L.)、杜仲(Eucommia ulmoides)、山茶(Camellia japonica)、南方铁杉(Tsuga chinensis var. tchekiangensis)竹柏(Podocarpus nagi)、江西槭(Acer kiangsiense)等。

本研究所用的NDVI数据产品为1999—2015年MODIS NDVI数据, 来源于美国国家航空航天局(NASA)的EOS/MODIS植被指数产品数据MOD13A2(https://wist.echo.nasa.gov/ap)。数据空间分辨率为1 km×1 km, 时间分辨率为16d, 获取的数据利用Envi 5.3对原始数据进行预处理。气象数据来源于国家气象科学数据中心(http://data.cma.cn), 研究区附近的57个基准气象台站的逐月数据(月平均气温、月降水及月辐射等), 时间序列为1999—2015年。

Theil-Sen Median趋势分析计算时间序列数据的所有n×(n-1)/2成对组合之间的斜率的中位数。其计算公式为[22]：

式中, ß是植被变化趋势, NDVIj是第j年的NDVI值, NDVIi是第i年的NDVI值。如果ß＞0, 则NDVI有增加的趋势, 表明在时间段内植被一直在改善或恢复。如果ß＜0, 则NDVI呈下降趋势, 表明在时间段内植被呈现退化的趋势。

Mann-Kenddall检验可以判断趋势显著性, 其样本无特定的分布要求, 其结果也更为准确, 不受异常值的干扰[23]。计算公式如下：

设定NDVIj, j=1999, 2000, …, 2015, 定义Z统计量为：

式中, NDVIi和NDVIj表示像元i年和j年的NDVI值；n是时间序列的长度；sgn是一个符号函数。统计量Z的值在(-∞, +∞)范围内。在给定显著性水平α下, 当|Z|>Z1-α/2时, 表示时间序列在α水平上存在显著的变化。在本研究中, 取α =0.05, 判断在0.05置信水平上即|Z|＞1.96时, 1999—2015年区域NDVI变化趋势的显著性。

变异系数通过计算标准差与平均值的比值而来, 又称为离散系数, 是量度各观测值与单位均值离散程度的数学指标[24]。

式中, CV是指NDVI值的变异系数, NDVIi是第i年的NDVI, 是区域1999年至2015年的年际平均NDVI。当CV值较大时, 数据更加分散且植被年际变化较大；当CV值很小时, 表明数据更紧凑植被更稳定。

基于重标极差法(R/S)的Hurst指数是定量描述时间序列长期依赖性的有效方法[25], 其计算公式是：

定义时间序列NDVI(t), t=1, 2, …, n., 对于任意正整数τ, 定义其均值序列：

计算累积离差：

计算极差：

计算标准差：

计算Hurst指数：

式中, 通过使用最小二乘法拟合等式log(R/S)n=a+H×log (n) 获得H值, 其中H为赫斯特指数。根据Hurst[25]和Mandelbrot[26]的观点, Hurst指数值的范围为(0, 1), 有3种类型：当H=0.5时, 表明NDVI时间序列是随机的序列, 没有可持续性, 表明了此变化趋势未来的时间序列与研究的时间序列无关；当H＞0.5时, 反映NDVI时间序列具有可持续性, 表明未来趋势变化与当前NDVI时间序列趋势相同；当H＜0.5时, 此类未来变化趋势有可能保持基本稳定或者发生逆转, 本研究中在空间变化上将其归为变化趋势不确定[26]。

利用1999—2015年研究区及周边57个国家级气象站逐日记录的气象数据(气温、降水、相对湿度、日照时数)进行最大值合成计算出月最大值, 基于月最大值求取年均值, 形成站点逐年数据集, 利用ArcGIS软件的克里金插值法得出气象因子空间分布。

地球系统是由多要素所构成的复杂系统, 系统中任何一个要素的变化都会引起其它要素也产生变化, 偏相关分析便可有效解决上述问题, 单独研究两要素之间的相关性时, 偏相关将其它要素的影响视为常数[27]。计算公式为

式中, rxy为变量x与y之间的相关系数, i为样本数, 、分别为NDVI多年均值和各气象因子均值。

式中, rxy、rxz、ryz分别为要素x与y、x与z、y与z之间的相关系数；rxy·z为将要素z固定后要素x与y的偏相关系数, 通常采用t检验法对偏相关系数进行显著性检验。

1999—2015年NDVI数据的年平均值变化成体呈现增加的趋势, 最低值为2001年的0.72, 最高值为2015年的0.80, 整个研究时段内的平均值为0.76。研究区植被覆盖向好的趋势发展, 可大致分为3个阶段：第一阶段为缓慢下降期(1999—2001年)、第二阶段为波动上升期(2002—2011年)、第三阶段为平稳上升期(2011—2015年)。区域NDVI在2000—2005年、2006—2010年、2011—2015年的增长率分别为4.86%、2.23%、4.81%。在2000—2005年与2011—2015年时间段内, 研究区植被得到了显著的改善(图 2)。

利用1999—2015年的年平均NDVI数据, 计算16年间平均值得到长江中游地区NDVI空间分布图(图 3)。研究区植被覆盖格局整体呈现“西高东低、北高南低”的分布特征, 西部地区海拔较高, 以山区为主, 天然林、灌木和草地等分布广泛, 植被生长良好, NDVI值较高。北部地区以平原为主, 是中国重要的商品粮基地, NDVI值同样较高。NDVI低值区主要分布在三生省会城市(湖南省长沙市、湖北省武汉市和江西省南昌市)及其附近城市地区, 此部分区域多数经济发展迅速, 人工扰动较大。16年间, NDVI值在0.1—0.5段占比1.40%, NDVI值在0.5—0.7段的区域分别占11.45%, NDVI值从0.7—0.8的区域占55.06%, NDVI值大于0.8的区域占32.09%。

长江中游地区1999—2015年NDVI空间格局变化特征表现为, 植被明显改善区域占最大比例(67.39%), 主要分布在研究区西部及东南部；稳定不变及轻微改善区域占区域面积的四分之一(25.79%), 主要以江汉平原、洞庭湖平原、鄱阳湖平原为中心分布；退化区域占区域面积的6.82%, 零星散布在各大城市的中心区位置, 并以中心区向外扩散分布；严重退化位置以三省省会武汉市、长沙市和南昌市为中心向外辐射, 表明城市规模的扩大对NDVI带来显著影响(表 1和图 4)。

长江中游区域1999—2015年NDVI稳定性整体表现为“高低波动并存, 低波动占比高, 区域性明显”。空间格局变化特征表现为：低波动区(0.014＜CVNDVI＜0.047)与较低波动区(0.047＜CVNDVI＜0.054)占全区面积的59.32%, 植被变化较小；中度波动区(0.054＜CVNDVI＜0.087)占全区32.59%, 环簇于较高波动区边界并向外延伸；较高波动区(0.087＜CVNDVI＜0.250)占7.00%, 在区域中部形成连接三省省会的带状分布, 其余呈斑块状分布在各城市中心, 较高波动区NDVI表现为轻微退化及严重退化(表 2和图 5)；高波动区(CVNDVI＜0.250)占比仅为1.09%。

长江中游地区1999年至2015年间, 植被Hurst小于0.5的区域占总面积的60.54%, 大于0.5的区域占39.46%。将NDVI变化趋势与Hurst指数结果耦合叠加, 得到研究区NDVI变化趋势与持续性的耦合空间图(图 6), 将结果分为6类(表 3)：持续性严重退化和持续性轻微退化, 其比例为0.79%和2.47%, 主要以大斑块分布在长沙、南昌、武汉省会发达城市及其周边；持续性稳定不变, 其比例达1.42%；持续性轻微改善和持续性明显改善, 其比例为9.78%和25%, 主要分布在低山丘陵区；未来变化趋势不确定, 此类占比高达60.54%, 均匀分布在各个省市。

对长江中游地区1999—2015年NDVI与各气象因子在年际变化尺度上, 分别做了相关、偏相关分析(表 4)。结果表明, 区域整体NDVI与气温、降水呈正相关, 偏相关系数分别为0.488、0.277, 与相对湿度、日照时数呈负相关, 偏相关系数分别为-0.646、-0.532。NDVI与气温的相关性大于与降水的相关性。NDVI与相对湿度的相关性大于日照时数的相关性, NDVI与各气候因子的相关系数及偏相关系数均未达到显著水平(P＞0.05)。表明在年际变化尺度上, 长江中游地区整体区域NDVI对同期气候变化的响应不是特别强烈。

分析每月气候因子与NDVI之间的相关性, 可以更好地揭示气候因子是如何在月际尺度上影响NDVI变化的。对1999—2015年12个月各月份的NDVI年际变化与当月及前1—3月的气温、降水、相对湿度、日照时数进行双变量相关性分析, 得出对应的相关系数表(表 5)。结果显示, 1—4月份NDVI生长受降水和日照影响较为显著, 并且具有明显的滞后性；6月NDVI与当月相对湿度具有显著的相关性；7月NDVI受5月降水的滞后性影响；8月与当月的气温、相对湿度、日照时数均存在显著的相关性, 受气候因子的影响最大；9月受当月降水和相对湿度的共同影响；10月受9月的降水、相对湿度和日照的影响；11月与各个气象因子均无显著相关关系；12月受当月相对湿度与日照时数的影响。降水、相对湿度和日照时数是影响长江中游区域NDVI季节变化的主要气候因子, 降水和日照时数的时滞性较为明显。

从空间相关性来看, 研究区绝大部分地区NDVI年际变化与气候因子的相关性并不显著, NDVI与气温、降水呈正相关的面积比例大于负相关的面积比例(表 6和图 7)。NDVI与气温和降水呈正相关的面积比例占71.82%和61.56%, 其中, 在P＜0.01(气温：1.85%；降水：0.63%)和P＜0.05(气温：5.39%；降水：2.82%)水平上达到显著。NDVI与相对湿度和日照时数呈负相关的面积比例占66.38%和60.51%, 其中, 在P＜0.01(相对湿度：1.90%；日照：0.32%)和P＜0.05(相对湿度：4.66%；日照：2.28%)水平上达到显著。总体来看, 长江中游地区NDVI与各气候因子达到显著的相关性比例分别为(图 7和图 8)：7.53%(气温)、7.31%(相对湿度)、4.12%(降水)、3.76%(日照时数)。表明在研究区大面积范围内气温升高、降水增多有利于植被的动态生长。

1999—2015年, 长江中游地区NDVI呈总体上升趋势, 这与湖北省[28]、湖南省[29]、江西省[30]、长江流域[12]等区域的研究果相一致。长三角、珠三角等沿海城市群的NDVI年均值多在0.3—0.5之间, 与之相比, 长江中游地区NDVI相对较高, 多年均值为0.76。从NDVI的变化趋势来看, NDVI稳定不变及轻微改善区域占区域面积的四分之一(25.79%), 主要以江汉平原、洞庭湖平原、鄱阳湖平原为中心展布, 这三大平原是我国主要商品粮基地, 符合农作物轻微变化特征。退化区域占区域面积的6.82%, 以三省省会武汉市、长沙市、南昌市为中心向外辐射, 表明城市规模的扩大对植被覆盖带来显著影响。

从NDVI的波动性来看, NDVI较高和高波动区占7.00%, 在区域中部形成连接三省省会的带状分布(长沙市-岳阳市-咸宁市-武汉市-黄冈市-九江市-上饶市-南昌市), 其余呈斑块状分布在各城市中心, 对应轻微退化及严重退化区域。中度波动区占全区32.59%, 环簇于较高波动区边界并向外延伸, 低波动区与较低波动区占全区面积的59.32%, 植被变化较小, 与NDVI稳定不变及轻微改善相吻合。从NDVI的未来变化趋势来看, 未来变化趋势不确定的区域占研究区一半以上, 这与前人在长江流域的结果一致[11]。持续性明显改善区域占研究区四分之一, 其斑块分布在各个省市内, 主要分布在区域西部山地区, 十堰市、神农架林区等地。持续性严重退化区域(0.79%)则主要以大斑块分布在长沙、南昌、武汉省会发达城市, 在襄樊市、宜昌市、常德市等人类活动频繁的较发达城市呈小斑块分布。

Nemani等和Liu等认为水热气候条件是影响陆地植被覆盖空间格局的驱动因素[31-32]。本文研究表明, 在年际变化上气温是影响研究区植被生长的主要气候因素, 这与崔林丽等[33]、李本纲等[34]和陈云浩等[35]对我国华东及周边地区植被NDVI对气候响应的研究结果基本一致, 其原因可能是研究区降水较丰富, 降水足够满足植被生长需求, 而热量的差异便是造成NDVI差异的主要驱动因素。本文研究表明, 降水对该区植被生长的影响为滞后2—3个月, 日照时数为滞后2个月, 相对湿度为滞后3个月, 表明植被NDVI变化也受土壤水分和养分有效性的影响。这与Bao等[36]从全球、区域等多尺度研究证实结果一致, 植被覆盖对气候变化的反馈存在一定的滞后效应。从空间相关性来看, 研究区与各气候因子达到显著相关的比例为：气温＞相对湿度＞降水＞日照时数。

(1) 在年际尺度上, 研究区近16年NDVI呈上升趋势平均值为0.76, 植被状况较好。在空间尺度上, 分布呈现“西高东低、北高南低“的分布特征。

(2) 从变化趋势来看, 研究区NDVI改善区域占比一半以上, 稳定性表现为“高低波动并存, 低波动占比高, 区域性明显“。区域NDVI未来变化趋势不确定的区域占比60.54%, 需要引起持续关注和重视。

(3) 从植被与气候因子的相关性来看, 年际尺度上, 研究区NDVI与各气象因子关系均不显著；月际尺度上, 降水、相对湿度和日照时数是影响NDVI变化的主要因子, 降水和日照时数有明显的时滞性；空间相关性上, 植被变化与气温、降水呈正相关的面积大于负相关, 与相对湿度、日照时数呈负相关的面积大于正相关。