1990—2019年阿拉尔垦区植被覆盖时空变化特征分析

宋 奇，高 琪，马自强，王 楠，王明玥，彭 杰*

(1.塔里木大学植物科学学院，新疆 阿拉尔 843300；2.北京大学地球与空间科学学院，遥感与地理信息系统研究所，北京 100871；3.浙江大学环境与资源学院，浙江 杭州 310058)

植被作为全球陆地生态系统的重要组成部分，其通过光合作用、呼吸作用与生物圈其他自然要素形成紧密联系，在生态系统的能量流动和物质循环中扮演着重要角色[1-2]。作为连接土壤、水、生物和大气等要素的重要纽带，植被不仅具有维持区域气候稳定和生态平衡的功能，而且能够反映区域绿洲和荒漠2大生态系统的相互作用过程，为区域生态研究提供重要信息[3-4]。此外，植被还具有防风固沙、涵养水源、调节区域小气候、治理大气污染和美化环境等众多生态功能[5-6]。明确植被覆盖时空变化规律，探明气候变化和人类活动对植被覆盖变化的影响，对区域环境质量评价和生态系统维护具有重要意义[7]。

目前，研究植被覆盖变化的主要方法有地面调查法和遥感监测法[8]。传统地面调查法需投入大量人力、物力和财力且不适合于大面积研究。随着遥感技术的不断发展，遥感监测法因其具有信息更新快、覆盖范围广、获取便捷、成本低、利用价值高等优势，成为监测植被覆盖变化的重要手段[9-10]。利用遥感手段分析植被覆盖空间分布特征已得到广泛运用[11]。在植被覆盖变化分析中，归一化植被指数(Normalized difference vegetation index,NDVI)是植被覆盖度遥感监测方法中最常用的植被指数，综合反映了区域植被覆盖密度和植被生长状况，是描述植被空间分布密度的最佳指标。此外，植被覆盖度(Vegetation Fractional Coverage,VFC)也是对区域植被长势进行量化的重要指标[12]，它是指单位面积内植被树冠在地面形成的垂直投影面积占单位面积的百分比。通过像元二分模型[13]计算归一化植被指数并对植被覆盖度进行提取的方法被广泛运用于植被覆盖动态监测研究[14]。

多数植被覆盖变化研究的处理方法是对单一的遥感影像进行NDVI提取展开相关研究[15-17]。然而，这样的研究未能捕捉连续性信息，容易造成关键拐点的遗漏，而且由于云覆盖、植被物候差异、耕地轮作休耕和病虫害的影响，单一遥感影像无法反映当年植被覆盖最全面的情景，难以准确分析研究区的植被覆盖情况。因此，基于长时序多时相合成数据进行的植被覆盖变化研究仍需进一步加强。

阿拉尔垦区地处我国天山山脉和塔克拉玛干沙漠的过渡带，属于西北干旱区典型的人工绿洲。阿拉尔垦区不仅是新疆主要的棉花(Gossypiumspp.)及特色林果生产基地，更是防止塔克拉玛干沙漠北移的重要生态屏障，在维持整个新疆地区生态系统平衡中发挥着重要作用。近些年，垦区人口不断增加，机械化作业水平逐步提高，人们为满足生活所需，盲目开垦土地资源，人工植被面积(农作物和经济林等)不断增加，自然植被面积(天然林地和草地等)不断减少，生态危机日益严重。

本研究收集了阿拉尔垦区1990—2019年共405景Landsat影像，并对其进行年内NDVI最大值合成，利用像元二分模型计算逐年的植被覆盖度，运用多年平均、图像差值法和面积加权重心转移模型分析垦区植被覆盖的空间分布特征、空间变化特征和重心转移情况，并探讨了气候变化和人类活动对垦区植被覆盖变化的影响。旨在探明人工绿洲区域的植被覆盖变化过程，为实现垦区水土资源合理利用、植被保护、生态安全和可持续发展提供理论支撑。

1 材料与方法

1.1 研究区概况

阿拉尔垦区(80°30′0″～81°58′0″ E，40°22′0″～40°57′0″ N)，地处天山南麓，塔克拉玛干沙漠北缘，内含胜利、上游和多浪3大水库。垦区面积为4 105.91 km2，主要土地利用类型为草地和林地，分别占垦区面积的25.67%和14.83%。垦区主体部分以塔里木河为中轴线，次级河流和灌溉渠道为纽带，构成密集网状水体廊道，形成了较为成熟的人工农田绿洲景观。垦区气候具有典型暖温带极端大陆性干旱荒漠气候特征，温暖干燥、光热资源丰富，日照率58%～69%、年均日照时数2 556.3～2 991.8 h、年均蒸发量1 876.6～2 558.9 mm、年均降水量40.1～82.5 mm。垦区多受扬尘和沙暴等恶劣天气影响。

1.2 数据来源与预处理

从美国地质勘探局USGS网站(http://giovis.usgs.gov)中获取阿拉尔垦区1990—2019年间所有Landsat遥感影像(影像轨道行/列号：146/32，分辨率：30 m)，从中筛选出云覆盖度低于100%的可用影像共计405景，图1为各月份的遥感影像云覆盖及每年的影像数量情况，其中1990—1998年选用Landsat 5 TM影像，1999—2012年选用Landsat 7 ETM+影像，2013—2019年选用Landsat 8 OLI影像。利用ENVI 5.3对影像逐一进行辐射定标、大气校正、几何校正、裁剪和图像增强等预处理操作。1990—2019年阿拉尔垦区(区站号51730)气象数据来自《地面气象记录月报表》，土地利用现状遥感监测数据是我国精度最高的土地利用遥感监测产品，此类遥感监测数据来自中国科学院资源环境科学数据中心的共享平台。

图1 遥感影像的云覆盖度及影像数量情况

1.3 数据处理与分析

1.3.1NDVI最大值合成 对筛选出的405景可用影像进行预处理后，逐景提取NDVI[18]，计算公式为：

(1)

其中，NIR近红外波段反照率(Near infrared spectrum,NIR)，R为红波段反照率(Spectrum,R)。

进而将每年各月份影像进行NDVI周年最大值合成，得到1990—2019年间30景NDVI合成后影像(图2)。由图2所示，合成后影像的植被覆盖度更高，更能全面反映当年植被覆盖最全面的情景，有利于植被覆盖度的计算，从而得到更精细的解译结果。

图2 NDVI最大值合成流程图

1.3.2植被覆盖度计算 基于NDVI和植被覆盖度之间存在的极显著线性相关关系，运用像元二分法模型对植被覆盖度进行计算[19]。其公式为：

(2)

式中，NDVIsoil为无植被覆盖的裸土像元NDVI值，NDVIveg为全植被覆盖像元NDVI值。取累计直方图的 5%和 95%为置信度区间，确定研究区的NDVIsoil和NDVIveg值[20-21]。

参考中华人民共和国水利部2008年发布的《土壤侵蚀分类分级标准》，结合本研究区实际情况和相关研究[22-24]，将植被覆盖度阈值按照等间距法划分为5个等级：裸地或极低植被覆盖(0～0.15)，等级为1；低植被覆盖(0.15～0.3)，等级为2；中植被覆盖(0.3～0.45)，等级为3；高植被覆盖(0.45～0.6)，等级为4；极高植被覆盖(0.6～1)，等级为5。

1.3.3图像差值法 为反映研究区植被覆盖空间变化特征，利用图像差值法计算不同影像之间的植被覆盖变化ΔFg，差值大于0表示植被覆盖增加，小于0表示植被覆盖减少[25]。其公式为：

ΔFg=Fyears2-Fyear1

(3)

式中，Fyear1为前一时期植被覆盖度等级，Fyear2为后一时期植被覆盖度等级。

为突出1990—2019年植被覆盖空间变化情况，对1990年和2019年植被覆盖度提取结果进行差值量化分析：其中，ΔFg=3为极度改善；ΔFg=2为中度改善；ΔFg=1为轻微改善；ΔFg=-3为极度退化；ΔFg=-2为中度退化；ΔFg=-1为轻微退化；ΔFg=0为稳定。

1.3.4植被覆盖面积加权重心转移分析 重心是地理学中描述地理要素或空间对象转变的重要指标。重心的动态转移变化反映了地理要素空间分布的转移轨迹，加权重心则是通过对重心坐标赋予权重来表示地理现象分布的不均匀性，选用面积加权重心模型计算不同等级植被覆盖度的重心，并反映研究时段内垦区植被空间演变过程[26]。其公式为：

(4)

(5)

式中，X，Y分别为垦区植被分布重心的经纬度坐标；Ci表示第i个垦区植被分布斑块的面积；Xi，Yi分别表示第i个垦区植被分布斑块分布重心的经纬度坐标。

2 结果与分析

2.1 NDVI年内合成最大值验证和植被覆盖精度检验

考虑到每个年份使用的遥感影像数量有限，有可能造成年内NDVI被低估，因此对本文的NDVI年内合成最大值进行验证。将本文的NDVI年内最大值与中分辨率成像光谱仪(Moderate resolution imaging spectroradiometer,MODIS)的NDVI产品计算的年内最大值进行对比分析，确保本文计算的年内NDVI最大值未被低估。获取了阿拉尔垦区2010年间所有MODIS的NDVI产品进行最大值合成，随后与Landsat影像的NDVI最大值进行对比，如图3所示，分别将两产品的NDVI最大值合成后结果进行两处细节放大，从中可以明显看出Landsat产品的NDVI最大值合成后结果更清晰，NDVI值更高，效果最佳。并且Landsat和MODIS产品都是16 d一期，两产品的空间分辨率不同，Landsat是30 m，MODIS是250 m，明显Landsat影像的成像效果更佳，因此本文选择Landsat影像进行NDVI最大值合成。

图3 中分辨率成像光谱仪和遥感影像产品的归一化植被指数最大值合成后比较

为验证所得植被覆盖情况的可靠性，根据阿拉尔垦区的土地利用/覆被现状及相关参考文献[27-28]，经综合对比分析，建立了一套适合于阿拉尔垦区的解译标志。以2010年为例，各植被覆盖类型的地表特征、影像特征和实地调研情况如表1所示。

表1 阿拉尔垦区解译标志

根据该解译标志建立各植被覆盖类型所对应的验证样本，样本点个数及分布情况如图4所示，各样本点的分布遵循均匀、全方位覆盖原则，最终各样本所占的像素点个数分别为：极高覆盖100 603个、高覆盖60 387个、中覆盖10 394个、低覆盖53 570个、极低覆盖240 571个，同理将其余年份逐一建立验证样本进行精度验证。

图4 每类样本点个数和分布图

进而对其建立混淆矩阵进行各植被覆盖类型的精度评价。由表2可知，阿拉尔垦区的总体精度(Overall accuracy,OA)为88.50 %，Kappa系数为0.83。相关文献表明[29-30]，总体精度OA值和Kappa系数均大于0.8时，所得的结果可用于相关研究，因此本研究的结果可信。

表2 阿拉尔垦区2010年植被覆盖精度评价

3.2 植被覆盖空间分布特征及变化趋势

基于NDVI像元二分模型，计算出垦区1990—2019年不同等级的植被覆盖度，以此为依据进行等级划分，得到阿拉尔垦区1990—2019年的植被覆盖等级分布情况(图5)。由图5可知，阿拉尔垦区的植被空间分布总体以塔里木河为轴线，从内向外扩展。极高和高植被覆盖区主要集中于垦区中部，呈现出大片块状集中分布的趋势；中植被覆盖区主要集中于高植被覆盖区外围，植被类型多以人工防护林为主；低和极低植被覆盖区集中于垦区南、北两端，主要为稀疏的天然植被。

图5 阿拉尔垦区不同时期植被覆盖等级分布

统计得到的阿拉尔垦区平均植被覆盖度变化和不同等级植被覆盖面积变化情况(图6)。1990—2019年，阿拉尔垦区平均植被覆盖度呈波动式增加的趋势，拟合优度R2为0.9042，平均植被覆盖度在0.1和0.5之间变化，但存在年际差异：1991最低，约为0.1286；2018年达到峰值，约为0.4656。

极低植被覆盖区面积呈线性减少趋势明显(R2=0.9639)，面积由1990年的3 050 km2减至2019年的1 226 km2，降幅59.8 %。低植被覆盖区面积呈增加趋势显著(R2=0.9729)，增幅450%。30年间中植被覆盖区面积的变化波动较大，整体呈减少趋势(R2=0.1475)。高植被覆盖区面积的变化呈明显增加趋势(R2=0.9858)，增幅162%。极高植被覆盖区面积的年际间变化波动性减少(R2=0.9514)，面积变化呈增加趋势，增幅1 880 %(图6)。对比可知，在不同覆盖等级的面积变化中，极低植被覆盖区的面积减少最显著，高植被覆盖区的面积增加最显著。

图6 阿拉尔垦区平均植被覆盖度和不同等级植被覆盖面积变化

2.3 植被覆盖空间变化特征

由于每年的植被覆盖波动较大，因此用多年平均的方式来表征不同阶段阿拉尔垦区的指标覆盖变化和不同等级的植被覆盖面积情况。在ArcGIS软件中对阿拉尔垦区植被覆盖进行多年平均，结果如图7所示。

图7 阿拉尔垦区植被覆盖的多年平均变化情况

统计多年平均中不同等级的面积变化(表3)可知，1990—1995年，以极低覆盖为主(2 431.12 km2)，主要分布在垦区南北两端的偏远地区；1995—2000年，极低覆盖面积减少，面积为1 918.70 km2；2000—2005年，极低覆盖面积在所有覆盖等级中依然为最大值(1 694.10 km2)；2005—2010年，垦区西北部开始出现植被，以中覆盖类型为主；2010—2015年，极高覆盖面积高达772.73 km2；2015—2019年，极高覆盖面积成为所有覆盖等级中的最大值(1 113.94 km2)；1990—2019年，整体来讲，垦区还是以极低覆盖为主，而其他密度较高的覆盖类型主要分布在塔里木河区域。

表3 阿拉尔垦区植被覆盖多年平均中不同等级的面积变化情况 单位：km2

综上可知，阿拉尔垦区植被变化不仅存在时间阶段性差异，同时存在区域性差异。为进一步分析阿拉尔垦区植被覆盖面积的变化情况，运用转移矩阵对1990年和2019年不同等级的植被覆盖面积进行计算。1990—2019年，阿拉尔垦区植被覆盖主要以极低覆盖向高覆盖转移和高覆盖向极高覆盖转移为主，具体转移量为：极低植被覆盖转向高植被覆盖的面积为1 119.48 km2，高植被覆盖转向极高植被覆盖的面积为792.49 km2。从最终变化量来看，极低植被覆盖的面积减少量最大，减少了1 861.54 km2，极高和高植被覆盖面积的增加量较大，分别增加了1 541.54 km2和1 376.37 km2(表4)，垦区的植被覆盖情况整体上表现出极低植被覆盖区的面积逐渐减少，高和极高植被覆盖区的面积逐渐增加趋势。

表4 1990—2019年阿拉尔垦区植被覆盖度等级转移矩阵

通过差值量化提取的阿拉尔垦区植被覆盖面积变化看出(图8)，阿拉尔垦区中心区域周围成为垦区植被覆盖面积退化的热点区域。同时垦区中心区域因经济相对发达，周边乡镇和公路干线集中，所以植被覆盖面积的退化程度更为显著。相对而言，离塔里木河流域较远的区域成为垦区植被覆盖面积改善的热点区域。

图8 1990-2019年阿拉尔垦区植被覆盖变化等级量化分布

2.4 植被覆盖重心转移

重心分布变化能够从空间上反映垦区植被分布的时空演变特征[31]。阿拉尔垦区植被覆盖重心变化结果显示(图9)，近30年，阿拉尔垦区植被覆盖重心不断朝东北方向转移，但在2005年后开始往相反的西南方向来回移动。具体来说，1990—2005年向东北方向转移；2005—2010年向西南方向转移；2010—2015再次向东北方向转移；2015—2019再次向西南方向转移。在各时段中，1990—1995年重心转移最明显，直线转移距离为4.12 km。这主要是因为1994年在阿拉尔垦区的东北区域大面积开垦耕地，致使垦区植被覆盖重心向东北方向移动，而在2006年国家在垦区西南区域大面积增设团场，使得垦区西南方向的覆被面积迅速增加，重心向西南方向转移，在2005年后垦区重心在东北和西南方向之间来回转移。

图9 1990-2019年阿拉尔垦区植被覆盖重心转移变化

3 讨论

3.1 气候变化影响

阿拉尔垦区人为种植作物多以高和极高植被覆盖的形式呈现，天然疏林灌丛和草地多以低植被覆盖为主，而自然气候的变化主要影响天然植被的生长[32]。因此，阿拉尔垦区年均降水量的增加对垦区天然疏林灌丛和草地产生了积极作用，这也正是垦区低植被覆盖向中植被覆盖转移的原因。

气候变化对植被覆盖的影响，还间接表现在对流域径流的影响[33]。阿拉尔垦区年均气温在1990—2019年呈现出增加趋势，特别是在2004年后，平均升高了1.09℃[34]。塔里木河水源主要来自天山山脉的冰雪融化，温度的升高将导致流域径流量的增加[35]，从新疆塔里木河流域管理局提供的径流量数据显示，2005年塔里木河年均径流量比2004年增加了27.6×108m3，为垦区植被扩张提供了一定的水源保障，有利于植被覆盖面积的增加，同时也是垦区近30年平均植被覆盖度增加的主要原因。

3.2 人类活动影响

人类活动对生态环境造成的影响与土地利用有关[36]。从中国科学院资源环境科学数据中心获取了阿拉尔垦区1990，2000，2010和2018年的土地利用数据(表5)。1990—2018年阿拉尔垦区耕地和城乡工矿居民用地的面积增加幅度最大，其中耕地面积从1990年的61.59 km2(1.5%)增至2018年的1 180.04 km2(28.74 %)，净增加1 151.3 km2，说明阿拉尔垦区植被覆盖面积的增加主要得益于耕地的扩张，扩张最明显的区域位于垦区西北部(图8)。从表5中还可看出，垦区未利用土地面积减少。据1990—2018年阿拉尔垦区土地利用变化转移矩阵计算，阿拉尔垦区1 750 km2的自然草地和2 667 km2的未利用土地被转移成为耕地。

表5 阿拉尔垦区不同时期土地利用结构

人类活动除通过开垦农田、种植作物影响垦区植被覆盖面积变化外，城乡基建对垦区植被覆盖面积变化的影响也十分显著[37-38]。耕地和城乡基建的扩张促进了垦区经济发展，也带来了巨大的用水压力，造成垦区生态用水不足，在垦区人工植被面积不断增加的同时天然植被面积不断减少。这种以牺牲天然植被(草地)换取人工植被(耕地)的方式，破坏了垦区的正常生态平衡。干旱区天然植被作为维持垦区生态稳定的重要屏障，其面积的减少势必加快了垦区荒漠化进程，威胁到垦区的长期发展。同时，塔里木河的过度引流造成下游输水减少，严重威胁下游的生态安全和可持续发展。

4 结论

本研究对阿拉尔垦区植被覆盖时空变化特征进行分析，主要结论为：1990—2019年，阿拉尔垦区平均植被覆盖度呈增加趋势。极低植被覆盖区的面积呈明显的减少趋势。植被覆盖空间变化存在时序性和区域性差异。整体上，垦区植被覆盖重心整体向东北转移。气候变化对阿拉尔垦区植被覆盖有一定的促进作用，但人类活动的影响最为直接。大面积耕地开垦是垦区植被覆盖面积增加的主要原因，这种以牺牲自然草地换取耕地的方式破坏了垦区的生态稳定和可持续发展。