澜沧江中下游流域植被NDVI时空演变特征

袁泽申，陈晨晨，李伯根，王雨春

(1.河北工程大学水利水电学院，河北 邯郸 056004；2.中电建生态环境集团有限公司，广东 深圳 518102；3.云南省水文水资源局，云南 昆明 650000；4.中国水利水电科学研究院，北京 100038)

植被作为土地覆盖最重要的组成部分，在水土保持、调节大气、维持生态系统稳定等方面发挥着非常重要的作用[1]。大规模植被变化使流域下垫面格局发生剧烈变化，对区域水循环条件和水文过程产生深远影响[2]。植被变化还会影响流域水沙过程，进而影响水库调度运行过程，因此，研究植被变化对流域水电开发具有重要意义。植被动态研究是土地覆盖研究的重要内容，归一化差值植被指数Normalized Difference Vegetation Index(NDVI)作为植被生长和植被覆盖度的最佳指标[3]，能够反映时空尺度的植被覆盖度信息，被广泛应用于生态、环境和农业领域。澜沧江流域是湄公河流域的源头地区，而澜沧江流域的上游地区是高原地区，植被覆盖度极低，整个澜沧江流域的植被大量集中在中下游地区，所以澜沧江中下游流域的植被动态是整个湄公河流域生态环境质量的关键因素。在过去的20年时间里，中国兴建了大量的水利水电设施，其中极大部分的梯级水库坐落在澜沧江流域。梯级大坝的建设改变了河道的走向，而河道泥沙输运作为水利发电安全的重要影响因素，又被植被覆盖的变化程度影响。因此，研究澜沧江中下游流域的植被覆盖，对整个澜沧江流域乃至湄公河流域的生态环境保护和水利水电设施建设都有重要意义。

澜沧江被称作“东方多瑙河”，是中国发展与东南亚关系的重要纽带，同时也是湄公河的重要上游源头，其生态环境变化对亚洲地区生态及国际关系都有深刻影响。因此，对澜沧江中下游流域地区植被的分布规律及其动态演变经行研究，具有重大的生态价值和现实意义。OUYANG等[4]分析了湄公河上游植被覆盖的垂直差异变化，得出了湄公河上游流域植被覆盖NDVI在3 200 m以下和4 500 m以上的海拔高度对气候因子的变化较为敏感。LI等[5]分析了梯级水电站对澜沧江中下游河岸和高地植被的影响，发现梯级水电站对河岸和高地植被分布的影响比单个水电站更复杂，上游对植被分布的影响可能是与洪水相关的永久性变化，对植被分布的下游影响可能是与水文状况变化相关的动态相互作用。张景华等[6]使用2000—2010年MODIS NDVI数据和气象站台数据发现，澜沧江流域植被生长对气候响应表现出明显的滞后效应以及不同植被类型受气温和降水影响程度。

近年来，对澜沧江流域的植被覆盖研究多基于不同气候分区，较少有研究从整个研究区的尺度上对植被覆盖整体进行等级划分，进而分析其时间的变化和空间的转移，且研究对近几年的时间尺度涉及较少。本次研究基于1998—2018年NDVI数据使用Mann-Kendall模型耦合Sen’s Slope斜率估计分析了研究区域植被覆盖时空变化特征，分析其变化趋势和突变时间，对未来澜沧江中下游流域生物多样性保护、水循环过程研究和完善流域水资源管理措施提供数据支撑。

1 研究区域概况与方法

1.1 研究区域概况

澜沧江流域地处中国西南部，流经西藏、云南两省(自治区)，出中国国境后被称为湄公河，是东南亚最大的国际河流[7]。澜沧江的云南下半段，属于澜沧江的中下游流域，流域面积达78 000 km2，见图1。研究区域地形南低北高，地貌复杂多样，山高谷深，地势起伏较大[7-8]，加之西南季风的影响使得流域内降雨量呈现自下游向上游递减的趋势[5]。根据中国地理空间数据云网站2015年的地形和土地覆盖分布数据，研究区域的主要土地覆盖是森林(占总面积的47%)和草地(占总面积的44%)。主要植被类型为亚热带和热带山地针叶林、亚热带和热带常绿阔叶林、高山嵩草草甸和杂草草甸[4]。

1.2 研究方法

1.2.1Mann-Kendall检验

Mann-Kendall是非参数检验统计，用于评估时间序列变化的显著性。最初由Mann在1945年提出[9]，后由Kendall进一步完善[10]，其优点是不需要特定的样本分布，不受离群值的干扰，在环境科学领域中得到了广泛应用[11-12]。其统计检验统计量Z值计算公式如下：

(1)

测试趋势是在特定的显著性水平上进行的。当|Z|>Z1-a/2，拒绝原假设，时间序列存在显著趋势。Z1-a/2由标准正态分布表得到，本研究采用显著性水平α= 0.05，Z1-a/2=1.96。

(2)

其中，n为时间数据的个数，xi和xj分别为时间序列i和j(j>i)中的数据值，确定n(n-1)/2个sgn(xj-xi)插值函数。

(3)

(4)

其中，n为数据点的个数，P为绑定组的个数，ti为第P组的数据值个数。如果没有绑定组，则可以忽略此汇总过程[13]。

1.2.2突变检验

首先定义待检验时间序列为X(x1，x2，…,xn)，以序列中后一个数值大于其之前所有数值的个数计数值为新序列Pk，见式(5)、(6)：

(5)

式中，i=2,…,n；j=1,…,n-1。

(6)

然后将Pk进行累加求和得Sk，见式(7)：

(7)

计算Sk的均值和方差，见式(8)、(9)：

(8)

(9)

求统计量UFk，见式(10)：

(10)

按照待检验时间序列的逆序列，重复式(6)—(10)进行计算，所得结果进行倒序并求负，得到UBk统计量序列。本次突变研究给定显著性水平α=0.05，临界值U0.05=±1.96。将UFk和UBk统计量以及显著性统计区域临界线(±1.96)，绘制到一张图上。UFk和UBk的交点就是变异的开始时间，如果交点在临界值之间，则突变达到显著水平。在确定具体突变时间后，将突变检验公式载入到ArcGIS中，计算出突变年份NDVI的UFk和UBk统计量的栅格影像，利用栅格计算器筛选处于显著性统计区域临界线内UFk和UBk交点的栅格影像，来进一步确定研究区域NDVI发生突变年份时，发生突变的具体区域。

1.2.3Sen’s Slope斜率估计

使用Sen’s Slope 斜率估计[14]估算2000—2019年研究区域NDVI的年际变化。Sen’s Slope斜率估计是估算多年NDVI时间序列的所有连续数据组合的斜率的中位数，见式(11)：

(11)

式中，Slope为像元回归方程的斜率，Ai为第i年的NDVI的中位数，n为研究的时间长度。如果Slope>0，则NDVI呈上升趋势，否则，呈下降趋势。Sen’s Slope斜率估计作为一种非参数方法，已被证明在估计随时间变化时非常可靠[11]。

1.3 子流域划分

为了了解研究区域子流域的多年NDVI的时空分布变化，利用QGIS(Quantum GIS)软件载入SWAT(Soil and Water Assessment Tool)模型，根据研究区域地形和主要水系分布进行子流域划分，经过水文计算后共得10个子流域S1—S10，子流域S1—S4、S5—S10分别位于澜沧江中游和下游流域。

1.4 数据来源

以澜沧江中下游流域为研究对象，主要采集了该流域1998—2018年21年间的NDVI数据。NDVI数据来源于地理空间数据云中Landsat 4-8 TM卫星影像，下载范围按照澜沧江流域经纬度(北纬21°30′至32°40′，东经94°至101°50′)下载。该影像提供7个波段的影像图，取近红外和红波段影像来计算NDVI(NDVI=(IR-R)/(IR+R))。利用ArcGIS将所得卫星数据进行预处理，批量拼接、格式转换、重投影等，最后再采用最大值合成法MVC[15-16](Maximum Value Composites)得到1998—2018年逐年NDVI数据，与中国科学院资源环境科学与数据中心的NDVI数据进行对比矫正后，再将其定义投影为WGS_1984_UTM_Zone_47N，从而进行分析研究。

2 结果与分析

2.1 澜沧江中下游NDVI时间变化特征

研究区域1998—2018年NDVI时空分布特征见图2。利用ArcGIS软件将研究区域各年的平均NDVI值提取出来见表1。

年份NDVI年份NDVI19980.7720090.8119990.7620100.8020000.7520110.8020010.7620120.8220020.7620130.8320030.7820140.8620040.7920150.8620050.8020160.8620060.7920170.8720070.8020180.8620080.80——

研究区域1998—2018年平均NDVI线性回归分析结果见图3。研究表明，NDVI整体上呈微弱上升趋势，平均增长速度为5.8×10-3/a。而在1998—2018年中，NDVI最大的是2017年为0.87，最小的是2000年为0.75，总平均值为0.81，其中有8 a的年平均NDVI大于总平均水平，分别是2009、2012—2018年。研究区域1998—2018年平均NDVI分段趋势变化分析结果见图4，将21年数据平均分为5段，除1998—2002年段出现下降趋势外，其余4段均呈上升趋势。1998—2006年2个分段的年份NDVI无显著性，其余3段均满足5%的显著性水平。其中2010—2014年NDVI值出现明显提升，5年增加了0.066，达到了年平均增加值的2.7倍。以上结果表明：研究区域NDVI是在2000年后开始出现增长趋势，且在2011—2014年有明显提升。

通过Mann-Kendall趋势检验对1998—2018这21年研究区域的NDVI进行趋势分析，总体21年的NDVI的检验统计量U为0.15，小于临界判别统计量Uα/2=1.96(p=0.05)，所以澜沧江中下游流域的NDVI在1998—2018年呈不显著递增趋势。Mann-Kendall突变分析的结果见图5，其中在显著性水平p=0.05下的临界值为±1.96，UF与UB在显著性水平内出现交叉，这表明突变年份为2001年。综上所述，研究区域植被覆盖率呈上升趋势但递增不显著，且在2001年发生突变后开始上升，虽然在2006、2010、2011、2015、2018年有下降，但并不影响整体的上升趋势。

2.2 澜沧江中下游NDVI空间分布特征

采用子流域划分来研究NDVI的空间分布特征，更能体现研究区域NDVI的区域性特征。将研究区域划分为10个子流域，10个子流域的NDVI在1998—2018年内整体上均呈缓慢递增趋势，见图6，其中S8、S9和S10这3个子流域的NDVI年平均递增速度略低于整个研究区域的递增速度，它们的年平均递增速度分别是4.5×10-3、4.4×10-3、5.2×10-3/a；而子流域S6的NDVI年平均递增速度最高，其NDVI的递增速度为7.3×10-3/a。由上述数据可知：在整个研究区域内，中部地区的NDVI年递增速度较大，两头较低，其中以南部的S8和S9 2个子流域最为明显。

各年份的NDVI变化在空间上均呈现自北向南逐级递增的趋势，并有突出的区域特征性，以高程为列，结合图1可看出，高程在1 700 m以上的区域，S1、S2、S3子流域在同一时期NDVI相比较其他区域大多时候处于低水平状态；而高程在1 300 m又是另一个分界线，研究区域内的S8、S9、S10 3个子流域是高程1 300 m以下的区域，在同一时期NDVI相比较其他区域大多时候处于高水平状态。从图7中可以看出中下游子流域平均NDVI范围在0.73～0.89，大部分处于0.82～0.86，每个子流域从1998—2018年都呈现递增的趋势，但是研究区域上游子流域的变化趋势与研究区域整体的变化趋势更为接近。

将Sen’s Slope斜率估计与Mann-Kendall检验模型耦合分析研究区域NDVI的空间变化特征。由表2可知，严重退化区域占总体的0.11%，轻微退化区域占总体的1.39%，稳定不变区域占总体的0.77%，轻微改善区域占总体的89.9%，明显改善占总体的7.83%。从以上数据可得知，澜沧江中下游流域NDVI呈轻微改善趋势，小部分区域有明显改善的特征。

表2 研究区域NDVI空间变化特征

分析NDVI在像元尺度上的变化趋势，见图8。从面积上看，轻微改善区所占比例最大，其次是明显改善区，严重退化区、轻微退化区和稳定不变区所占比例极小，明显改善区域主要分布在研究区域的下半段，属于降水较多的低纬度地区，剩余部分零星分布在研究区中部地区，退化区域大部分是处于城市群附近。再对研究区进行空间Mann-Kendall突变检验，得出22个突变点，突变点均出现在退化区域，其中子流域S1中有7个突变点，S7中有2个突变点，S8中有4个突变点，S9中有5个突变点，S4、S6、S10均有1个突变点，其余子流域没有突变点。且大部分出现在高纬度地区的大理州和低纬度地区的西双版纳州，也与趋势分析中NDVI年增速相对应。以上结果可知：研究区域几乎没有植被退化情况，且总体的植被覆盖正在缓慢提升，而提升较为缓慢的是处于高海拔地区的北部区域的大理以及NDVI本就处于高水平的南部区域的西双版纳。

3 讨论

总体而言，研究区域多年(1998—2018年)平均NDVI呈现出明显的区域特征性，其变化梯度由北向南递增，与高程变化相反。研究区域高程范围在480～4 137 m，从图7可看出区域性特征的2条分界线分别在1 700、1 300 m，高程在1 700～4 137 m区域的NDVI在同一时期处于低水平状态，而高程在480～1 300 m区域的NDVI在同一时期处于高水平状态，3个区域的NDVI有明显的层次感，充分体现了研究区域NDVI的区域特征性。在低高程地区NDVI处于高水平状态可能是充足的降水使土壤中的含水量大部分时间处于饱水状态，从而为植被的生长提供大量养分[6,17]；而在高高程地区气温极低，降水少，生长环境较为恶劣，从而抑制了植被的生长，使其NDVI处于低水平状态。但是近20年的全球气温升高也让高高程地区的植被生长环境得到缓解，在高高程地区植被有所增长但是速度缓慢[18-19]。研究区域大部分地区的NDVI增加，表明1998—2018年的植被变化总体呈上升趋势。这一结果与之前研究的中国西部地区NDVI值的变化趋势相似[20]。时间序列的趋势分析以及Mann-Kendall突变检验得出结果，2000年后出现增长趋势，证明2000年的西部大开发战略以及2002年的退耕还林政策的成功实施，使得研究区域的生态环境不断改善。植被变化与气候变化、人类活动和地理因素的影响密不可分[21]。除了自然因素的影响，人类发展对NDVI的影响也是不可忽视的，近30 a来的城市化和工业化建设，人类破坏了大量的植被以获取经济上的建树，当城市规模较大时，辐射效应对NDVI的负面影响显著，城市对植被破坏的影响远大于直接占用对NDVI的影响。城市规模的扩大导致建设用地的增加，使地表覆盖发生明显改变，由植被水域等变化为人工不透水表面，例如柏油路面、水泥路面和建筑物顶等。城市辐射引起的植被退化问题应引起重视。由于时间限制，没能在NDVI、地区GDP以及人口密度上做出分析，用以深度研究人为因素对NDVI的影响。

4 结论

基于地理空间数据云中Landsat 4-5 TM卫星影像分析澜沧江中下游流域1998—2018年NDVI时空演变机制，得出结论如下：①用Mann-Kendall分析得出研究区域1998—2018年NDVI呈不显著递增，在2001年出现突变；②将研究区域划分为10个子流域进行区域分析得知研究区域NDVI均呈现自北向南逐级递增的趋势，有明显的区域特征性；③Mann-Kendall模型耦合Sen’Slope估计对研究区域的NDVI进行空间分布特征分析得出的结果表明，研究区域整体上无突出变化，一直呈缓慢上升趋势，大理和西双版纳出现小范围的退化聚集，突变点大多出现在城市附近，其突变产生可能与城市发展有关。