淮河流域1999—2018年NDVI时空演变特征研究

高 峥,刘赛艳,秦 璇,徐柳昕

(扬州大学水利科学与工程学院,江苏 扬州 225009)

植被是全球碳循环的主要参与者之一,在陆面水分运移和能量传输中发挥着巨大作用。植被覆盖变化情况是区域生态环境变化的重要指示特征,能直接反映出区域生态环境对变化环境的响应[1]。目前,对植被覆盖变化的监测方法主要有地表实测和遥感测量[2]。其中,地表实测对人力和物力的消耗较大,且只能提供小尺度范围内的植被结构和分布状况,因此具有明显的局限性[3]。相比之下,遥感测量则多是以归一化差异植被指数(Normalized Difference Vegetation Index,NDVI)作为地表植被的特征参量,直观地揭示地表植被覆盖的变化趋势。

作为反映区域植被覆盖综合状况的一种指标,NDVI数据不仅获取简便,而且能够实时地反映植被分布特征及变化,在大尺度地表植被监测和研究方面有着巨大的优势[4]。目前,国内外学者基于NDVI时间序列,开展的研究主要集中在3个方面:①不同时间尺度(如月、季、旬、年等)NDVI的变化规律[5-6];②不同气象因子(如降水、气温等)对NDVI的影响[7-8];③不同人类活动(如土地利用类型、农作物种植类型等)对NDVI的影响[9-10]。这些研究表明,受当地气象因子和人类活动的影响,NDVI变化具有典型的区域特征,因此从区域尺度上研究植被覆盖的变化特征更有现实意义。尽管如此,以往研究多注重NDVI的时间变化特征,尤其是整体上升或下降的变化趋势及其空间分布,而对NDVI空间特征向量或者说空间模态变化规律的相关研究却不多。事实上,NDVI空间模态变化规律对于了解流域植被覆盖状况、生态水文变化过程、开展区域生态环境保护与治理具有重要的意义[1,10]。经验正交函数分解法(Empirical Orthogonal Function,EOF)能够刻画要素(如降水)场的空间分布特点,反映不同区域之间的空间变化关系,而且分解的空间结构多具有明确的物理意义[11],已被广泛用于气象水文分析领域[12]。因此,本研究拟采用EOF方法研究区域NDVI变化情况,揭示区域NDVI空间模态变化规律。

淮河流域地处中国南北气候过渡带,自然植被分布具有明显的地带性特点,研究该流域植被覆盖的变化特征,对于当地生态环境的保护与改善、自然资源的合理利用具有重要的意义。纵观现有研究成果,对淮河流域植被覆盖的研究多集中在流域整体尺度NDVI的变化规律及其影响因子上[13],而且植被数据较为陈旧,对流域不同子流域NDVI空间模态变化特征的研究也较少。因此,本文选取淮河流域为典型研究对象,拟在淮河流域1999—2018年的NDVI时间序列的基础上,采用Mann-Kendall(MK)趋势检验法和EOF法,研究该流域不同时间尺度上植被覆盖的变化趋势和空间模态演变特征,以期为流域植被分布特征、生态环境保护等提供参考。

1 研究区概况

淮河流域位于中国东部,地处111°55′～121°20′E,30°55′～36°20′N,西起桐柏山、伏牛山,东临黄海,南部与长江流域接壤,北部与黄河流域相邻,流域总面积为27万km2。流域总体属于暖温带半湿润季风气候区,春冬干旱少雨,夏秋闷热多雨,冷暖、旱涝转变剧烈。流域上中游主要的植被类型为亚热带落叶阔叶林,流域农作物多为双季水稻、油菜、豌豆、小麦等。

本文研究区域为淮河全流域,淮河流域2013年土地利用类型见图1,数据来源于中国科学院资源环境科学与数据中心网站(https://www.resdc.cn/)。同时,为了明晰淮河流域各子流域的植被覆盖变化特征,将淮河流域划分为4个子流域:上游、中游、下游和沂沭泗水地区。

图1 淮河流域土地利用类型及其4个子流域

2 数据与研究方法

2.1 数据来源与处理

NDVI数据主要来源于中国科学院资源环境科学与数据中心网站(https://www.resdc.cn/)提供的月植被空间分布数据集(1999—2018年)。数据空间分辨率为1 km,是基于连续时间序列的SPOT/VEGETATION NDVI卫星遥感数据,采用最大值合成法生成的月度(1—12月)植被数据集[14]。该数据集质量良好,应用广泛。基于ArcGIS软件平台和Python编程语言等对数据集进行了批量格式转换、裁剪等处理工作。

2.2 研究方法

MK法是一种应用范围很广的非参数趋势检验方法,具体可见文献[15],此处不再赘述。采用MK趋势检验法分析淮河流域NDVI的变化趋势时,置信水平取95%。此外,采用Sen趋势估计法[16]计算了淮河流域NDVI序列的线性倾向率。

对某一物理量X进行采样研究,采样的空间点数量为m,时间点数量为n。建立该物理量的样本矩阵Xm×n=[X1,X2,…,Xn],其中X的第j个真实空间场为Xj=[x1j,x2j,…,xmj]T,xij(i=1～m,j=1～n)表示在第i空间点、第j时间点对X的采样值。结合物理量X,阐述EOF法[14]的步骤如下。

a)对矩阵Xm×n进行距平化处理,得到样本距平矩阵A:

(1)

其中aij为xij的距平化值,即:

(2)

式中xik——在第i空间点、第k时间点对X的采样值。

b)计算矩阵A与其转置矩阵AT的交叉积,得到协方差矩阵C:

(3)

c)计算矩阵C的特征根λj(j=1～m,λ1>λ2>…>λm)和对应特征向量Vj(即X的第j个空间模态)所构成的矩阵V=[V1,V2,…,Vm],其中Vj=[v1i,v2i,…,vmi]T,vij(i=1～m,j=1～m)表示Vj的第i个元素。矩阵C与V满足:

C×V=V×Λ

(4)

其中Λ为对角矩阵,即:

(5)

d)计算空间模态Vj的方差贡献率δj,确定物理量X的主要空间模态:

(6)

由式(6)可知:δj越大,空间模态Vj的重要性越突出。

e)检验主要空间模态之间差异的显著性。为了判断空间模态是随机的噪声还是有物理意义的信号,需要进行显著性检验。利用North法[16]对特征值的误差范围进行显著性检验,在95%置信度水平下的特征根误差为:

(7)

式中 Δλi——特征值λi的误差范围,i=1～m。

由式(7)计算各特征值的误差范围,若某一特征值的误差范围与其前后特征值的误差范围没有交集,则对应的空间模态通过差异显著性检验,即该空间模态是有物理意义的。

本研究采用EOF法研究淮河流域NDVI空间变化特征,选取方差贡献率大于20%的空间模态作为淮河流域NDVI的主要模态,在95%置信水平下检验各主要模态差异的显著性。

3 结果与讨论

3.1 淮河流域NDVI的时间变化特征

淮河全流域及各子流域1999—2018年月NDVI的多年平均值见图2。由图2可知,淮河全流域及各子流域月NDVI的多年平均值变化均呈M形。就淮河全流域而言,月NDVI的峰值分别出现在4月(0.65)和8月(0.74),双峰之间的谷值出现在6月份(0.46)。具体而言,淮河流域2—4月春季来临,植被生长迅速,月NDVI值快速增加;5—6月第一茬小麦和水稻开始收割,月NDVI值明显减少,并在6月下降至谷值;7—8月第二茬水稻、夏玉米等粮食作物生长,月NDVI值开始恢复并进一步增加,在8月达到第2个峰值;自9月步入秋冬季节后,随着第二茬粮食作物的收割,全流域月NDVI值再次下降。由此可见,淮河全流域月NDVI的变化受该流域作物熟制的重要影响,这与严四英等[17]研究的结果一致。

图2 淮河流域1999—2018年月NDVI的多年平均值

由图2还可以得出,月NDVI均值由大到小依次为上游、中游、下游和沂沭泗水地区,这种均值差异跟各子流域的植被类型特征密切相关。具体而言,淮河流域上、中游地区林地较多,而沂沭泗水地区的丘陵地貌较多、岩体占比大而植被少[18]。此外,从双峰之间的谷值来看,上游月NDVI变化幅度最小,中、下游变化幅度最大。月NDVI变化幅度的差异主要是因为上游农田占比相对较少,受作物熟制影响较小,而中下游地区绝大多数的土地覆盖类型为农田,受作物熟制影响较大[18]。

淮河全流域及各子流域春季(3—5月)、夏季(6—8月)、秋季(9—11月)、冬季(12月至次年2月)和年度NDVI的变化过程见图3。由图3可以看出,淮河全流域及各子流域春季、夏季、秋季和冬季NDVI的数值分别位于0.46～0.68、0.54～0.75、0.40～0.60和0.20～0.48范围内,夏季NDVI最高,冬季NDVI最小。此外,淮河流域年度NDVI的数值位于0.42～0.58范围内,其中上、中游NDVI最大,沂沭泗水地区NDVI最小。

图3 淮河流域各季节和年度NDVI过程

淮河全流域及各子流域春、夏、秋、冬季和年度NDVI的变化趋势见表1。由表1可知,除植被生长缓慢的冬季外,淮河流域四季NDVI总体呈显著上升趋势。在4个子流域中,仅下游夏、冬季NDVI呈不显著下降趋势,而其他子流域四季NDVI均呈上升趋势,尤其是春秋两季,多为显著上升趋势。淮河下游与其他子流域存在差异的主要原因是:下游耕地面积占比较高,且多为水旱两熟连作作物,夏冬两季NDVI变化规律易作物熟制影响[18]。此外,淮河流域(除下游以外)年度NDVI均呈显著上升趋势,年度NDVI线性倾向率由大到小依次为上游(0.05/10a)、中游(0.04/10a)、沂沭泗水地区(0.02/10a)和下游(0.01/10a)。

表1 淮河流域各季节和年度NDVI的变化趋势

总体上,淮河流域NDVI时间序列在1999—2018年呈持续增加态势,这和当前“全球变绿(Global Greening)”[1,17]的研究结果是一致的。尽管有研究表明,全球气候变暖中大气中的二氧化碳增加能够让植物生长得更加旺盛,但是人类活动尤其是植树造林、退耕还林还草的贡献也是不可忽略的[1]。结合淮河流域近20 a退耕还林以及其他水土保持措施的实施情况,可知国家相关政策的实施在淮河流域植被恢复与生态保护工作中取得了较好的成效[19]。

3.2 淮河流域NDVI空间变化特征

将淮河4个子流域1999—2018年NDVI和四季NDVI集,按照(子流域×年份)的形式组成5个4×20数据矩阵,然后对这5个矩阵进行EOF分解,获得年和季节NDVI的主要空间模态。

淮河流域年NDVI的空间模态见表2。由表2可知,淮河流域年NDVI第一个特征向量的方差贡献率达到72.1%,第二个特征向量的方差贡献率达到22.8%,两者累计方差贡献率达到了94.9%,且两者的误差范围没有交叠。因此,第一和第二模态代表了淮河流域年NDVI的主要空间变化特征。

表2 淮河流域年NDVI的空间模态

淮河流域年NDVI的第一、二模态见图4。由图4a可知,第一模态特征向量值均为正值,表明淮河子流域近20 a植被覆盖变化一致,呈现同增同减的分布模式。特征向量高值中心位于上游,次高值中心位于沂沭泗水地区,中游地区为低值中心,说明上游与沂沭泗水地区NDVI变化幅度较大,中游和下游地区的NDVI变化幅度较小。由图4b可知,第二模态特征向量值大致以上中游交界处为界,东部为负值区,西部为正值区。其中,正值中心出现在淮河流域上游区域,负值中心出现在下游区域,特征向量值从内陆向沿海不断减小,表明淮河流域植被覆盖状况是由内陆向沿海递减,这和实际情况是相符的。

a)模态一空间分布

b)模态二空间分布图4 淮河流域年NDVI第一、二模态

淮河流域四季NDVI的主要模态见表3。由表3可知,淮河流域春、夏、秋季NDVI只有第一个特征向量的方差贡献率都超过了20%,因此第一模态代表了3个季节NDVI的空间分布特征。淮河流域冬季NDVI有2个特征向量的方差贡献率超过20%,累积方差贡献率达到91.2%,而且两者误差范围无重合区间。因此,第一、第二模态代表了淮河流域冬季NDVI的空间分布特征。

表3 淮河流域季节NDVI变化的主要空间模态

淮河流域四季NDVI的第一模态见图5。总体而言,淮河流域四季NDVI变化一致,呈现同增同减的空间分布模式。高值中心一般都位于上游,低值中心多位于中下游地区,特征向量值由西向东呈递减趋势,说明上游地区NDVI变化幅度较大,中下游地区的NDVI变化幅度较小,植被覆盖度由内陆向沿海变化逐渐降低。

a)春季

b)夏季图5 淮河流域四季NDVI的第一模态

c)秋季

d)冬季续图5 淮河流域四季NDVI的第一模态

淮河流域冬季NDVI第二模态见图6。淮河流域冬季NDVI第二模态的特征向量值有正有负,西部和东北部为负值区,中部和东部为正值区。其中,负值中心出现在淮河流域上游区域,正值中心出现在中游区域,特征向量值从内陆向沿海依次减小,反映出淮河流域冬季NDVI变化是由内陆向沿海递减的。

图6 淮河流域冬季NDVI的第二模态

4 结论

基于1999—2018年NDVI时间序列,采用MK趋势检验法和EOF法分析了淮河流域植被覆盖时空演变规律。

a)受流域作物熟制影响,淮河流域月NDVI的多年平均值整体呈M状。由于植被类型不同,淮河子流域NDVI的均值和变化幅度都存在一定的差异。

b)除植被生长缓慢的冬季外,淮河流域四季NDVI总体呈显著上升趋势。在农田面积占比较多的下游,夏、冬两季NDVI呈不显著下降趋势,而在淮河其他子流域,四季NDVI均呈上升趋势,尤其是春、秋两季,多为显著上升趋势。

c)淮河流域年NDVI总体呈显著上升趋势,和全球变绿的趋势是一致的,同时也表明近20年来退耕还林及其他水土保持措施,在淮河流域植被恢复与生态保护工作取得较好的成效。

d)淮河流域NDVI空间模态分布整体上同增同减,东部与西部差异较大的形式。且上游NDVI变化幅度较大,中、下游NDVI变化幅度较小,植被覆盖度由内陆向沿海逐渐降低。