基于1990—2019年多时相影像的干旱区绿洲景观格局分析

宋 奇， 史 舟， 冯春晖， 马自强， 纪文君，彭 杰， 高 琪， 蒋学玮

（1.塔里木大学园艺与林学学院，新疆 阿拉尔 843300；2.浙江大学环境与资源学院，浙江 杭州 310058；3.塔里木大学农学院，新疆 阿拉尔 843300；4.北京大学地球与空间科学学院遥感与地理信息系统研究所，北京 100871；5.中国农业大学土地科学与技术学院，北京 100193；6.新疆昌吉州地质环境监测站，新疆 昌吉 831100）

全球变暖、冰川融化、生物多样性降低等生态环境问题日益突出，加之频繁的人类活动，打破了正常的生态平衡，区域生态环境逐渐恶化［1-3］，对世界景观格局造成不小的影响［4-6］。为探明全球变暖环境下生态环境所受的影响，有必要对区域景观格局进行分析，探索人类活动对景观生态系统的影响［7-9］。对景观格局信息的提取和分析不仅能够表征空间的复杂性而且也决定着生态系统中资源环境的分布情况。对景观格局进行研究，能在无序的景观斑块镶嵌中发现潜在的生态学意义和景观规律［10-11］。其中，区域景观分类是分析景观格局变化的基础［12］，传统地面调查方法不仅费时费力而且适用范围较小。随着遥感技术的应用，遥感影像数据以其低成本、获取方式便捷、覆盖范围广、具有长时间序列等诸多优势，成为景观格局变化分析中首选的数据源［13-15］。

已有相关研究［16-20］表明，目前景观格局研究的遥感数据源多为典型的时间断面数据，以基于像元的分类方法对长时序的区域景观格局变化进行分析。然而，景观格局的变化很多情况下并不是在关键年份发生的，典型的时间断面数据很容易遗漏关键节点信息，难以把握准确的景观格局变化情况，对连续性信息的把握不足［21］。此外，基于像元的分类方法（例如支持向量机分类方法）容易出现碎斑，分类后处理相对繁琐，容易导致分类后精度偏低，最终影响区域景观格局的分析。

阿拉尔垦区从地理方位来看，南抵塔克拉玛干沙漠，北靠天山山脉，具有不可替代的景观生态屏障功能。垦区受人类活动影响较大，区域内的植被以人工植被为主，自然植被很少，景观结构单一，生态环境较为脆弱，景观格局变化较大。因此，本研究选择西北地区典型人工绿洲阿拉尔垦区为例，对长时序区域景观格局变化信息进行分析，能为西北干旱地区相关研究提供一定的参考价值。

本研究获取阿拉尔垦区近30 a 所有可用的Landsat 遥感影像，共226 景。对其逐年进行NDVI最大值合成，由此得到连续30 a 长时间序列的数据源，并对比多时相面向对象、多时相基于像元和单时相面向对象3 种分类方法的精度，选取本区域精度最高的分类方法进行景观分类，进而得到1990—2019年间阿拉尔垦区的景观分类结果，以此分类结果分别从斑块类型、景观水平和空间分布3 方面对研究区景观格局信息进行深入分析，并通过累积距平法探明区域景观格局的突变信息，并检验了其突变点的真伪。本文以连续长时间序列和NDVI最大值合成为创新点进行区域景观格局分析，能够从数据源开始就提高景观格局信息的精确度，以此保证后续景观格局分析的可靠性。

1 研究区域与数据来源

1.1 研究区概况

阿拉尔垦区位于阿克苏河、和田河及叶尔羌河的三河交汇处，地理位置为40°22′～40°57′N，80°30′～81°58′E（图1），总面积4105.92 km2，地势呈西北高东南低走向，平均海拔1012 m，生态系统属于荒漠-绿洲型，土壤类型主要为砂壤土，水资源分布主要以塔里木河及胜利、多浪和上游三大平原水库为主。研究区地处欧亚大陆腹地，塔克拉玛干沙漠北缘，属暖温带极端大陆性干旱荒漠气候，气象多变，常有沙暴、扬尘和霜冻等自然灾害，全年降水少而蒸发大。阿拉尔垦区以棉花和红枣种植为主，是中国最大的长绒棉种植基地［22］。此外，有小面积种植香梨、苹果、葡萄、玉米、小麦、高粱等作物。

图1 阿拉尔垦区位置示意图Fig.1 Location map of Alar Reclamation Area

1.2 数据来源

获取1990—2019年南疆阿拉尔垦区（轨道号为146—32）空间分辨率为30 m 的Landsat 5、7 和8 所有遥感影像，为提高后期遥感影像分类的判读精度，从所有影像中筛选出226 景云量低于40%的遥感影像作为本文基础数据，其中各年遥感影像的云覆盖率情况及影像数如图2所示。

图2 遥感数据Fig.2 Remote sensing data

2 研究方法

2.1 NDVI最大值合成

首先获取1990—2019 年所有Landsat 遥感影像，从中筛选出226 景可用影像，在ENVI 软件中将本文获取的所有遥感影像（共226景）逐一进行辐射定标、大气校正、图像配准、裁剪等影像预处理操作，进而逐一提取归一化植被指数（Normalized Difference Vegetation Index,NDVI），其计算公式为：

式中：NIR为Landsat近红外波段反射率；R为红光波段反射率。

相关合成的流程如图3 所示，得到30 a 间连续时间序列的NDVI合成图（图4）。从NDVI合成图中可看出，多时相遥感影像合成图不仅消除了云量的影响，而且图像中的植被覆盖度更高，有利于后期景观分类。

图3 技术流程图Fig.3 Technique flow chart

图4 NDVI最大值合成图Fig.4 Maximum NDVI composite

将1990—2019年各年份的NDVI年际变化信息（图5）进行汇总，可以明显发现近30 a 间，阿拉尔垦区的NDVI 处于增加的趋势，从1990 年的0.1661增加到2019 年的0.4319，植被覆盖度有明显的增加趋势。

图5 1990—2019年阿拉尔垦区NDVI年际变化Fig.5 Annual variation of NDVI in Alar Reclamation Area from 1990 to 2019

2.2 分类方法

参考《新疆生产建设兵团统计年鉴》（1990—2019）中的土地利用情况、研究区地表分布特点和新疆已有案例［23-24］中的分类方法，选用面向对象分类方法进行景观分类。该分类方法是综合考虑各景观类型的光谱、形状和灰度值等特征，利用eCognition软件对多时相遥感影像进行多尺度分割，经过多次测试，通过ESP（Estimation of Scale Parameter）评价工具找到最佳分割尺度，选用面向对象的随机森林（Random Forest）分类算法对其进行信息提取。本区域的耕地、林草地、园地和非植被的NDVI值不同，由此对其进行分类；此外，水体、建设用地和未利用地的几何和纹理特征区别较大，以此作为这3类用地的分类依据，构建合理的影像分类规则，进行景观分类。

2.3 景观格局变化分析方法

参考已有案例［25］，分别从斑块类型水平和景观水平对阿拉尔垦区景观格局变化特征进行定量分析。选用常用的景观指数［26］进行分析，分别从景观破碎度、优势度、复杂度、聚合度、连通性和多样性方面体现研究区的景观格局变化情况。由于本研究的时间序列较长，在长时序的景观格局变化中很可能会出现突变现象，我们选用累积距平法［27］探明研究区的景观格局突变情况，计算公式如下：

式中：LRi为第i年的距平累积值；Ri为第i年的景观指数值；为景观指数序列的多年平均值。

为了验证所得突变点的可靠性，笔者对累积距平法所得的突变点进行验证［28］，保证所得结果的可靠性；最后对景观指数进行可视化分析，选用不同大小的移动窗口，所得结果不同，因此在对比不同的移动窗口后，筛选出最佳窗口来反映景观空间格局的空间变化特征。最终从景观破碎度、优势度、复杂度、聚合度、连通性、多样性、突变点和空间分布等方面综合分析阿拉尔垦区景观格局变化特征。

3 结果与分析

3.1 不同分类方法的分类结果及其精度评价

以2019年为例，参考高分遥感影像并对研究区进行实地调查后，在基于像元的分类方法中共选取723951 个分类样本，其中耕地有173940 个，林草地有131591 个，园地有135138 个，水体有139291 个，建设用地有17488 个，未利用地有126503 个，进而得到相对应的分类后结果，最后对3 种分类方法所得结果进行对比（图6）。从图中可以明显看出，多时相面向对象分类方法的分类结果最佳，而多时相基于像元和单时相面向对象的分类结果较为破碎，容易出现“错分”和“漏分”现象。

图6 3种分类方法的分类后结果比较Fig.6 The results of three classification methods were compared

随机选取阿拉尔垦区351173个验证样本，选用常用的混淆矩阵方法对分类后结果进行精度评价。从精度评价结果中得出多时相面向对象分类方法的总体精度和Kappa 系数最高，分别为96.57%和0.95，说明利用多时相面向对象分类方法进行景观分类能够有效区分不同的景观类型，分类结果可靠。此外，最优分类方法多时相面向对象分类方法的总体精度相较于多时相基于像元分类方法提高了7%，Kappa系数相较于多时相基于像元分类方法提高了0.11，各景观类型的生产者精度均得到不同程度的提高；相较于单时相面向对象分类方法提高了12.22%，Kappa系数提高了0.25。

多时相面向对象方法的详细精度评价情况如表1 所示。从表中可以看出，该最优分类方法的总体精度OA达到了98.79%。同时各类型的用户精度UA 和制图精度PA 均在90%以上，说明该分类方法在本研究区是可行的。但水体和未利用地之间容易混淆，这是因为阿拉尔垦区水体附近存在小面积湿地，很容易被错分成为未利用地，但对整个研究区而言，被错分的湿地面积较小，不影响整体分类结果。同理，其他历史年份均按此方法重新建立样本和规则，逐年进行分类，即1990—2019 年共建立30 次分类样本进行景观分类。参考新疆现有案例［29-30］中的分类情况对本文的分类结果进行精度对比，对每年的分类结果重新建立精度评价混淆矩阵，保证每年分类的用户精度、制图精度和总体分类精度均在85%以上。

表1 阿拉尔垦区用地类型精度评价Tab.1 Evaluation of land type accuracy in Alar Reclamation Area

3.2 斑块类型水平上的景观格局分析

图7为30 a间阿拉尔垦区在类型水平上的景观格局指数变化情况。（1）斑块密度（PD）：耕地、林草地和园地的斑块密度较大。30 a 间，耕地类型的斑块密度数值呈现出增加的趋势，说明耕地类型的景观趋于破碎化；园地类型的斑块密度数值在2005年发生突变，数值突然增大，这很可能是与2005 年的政策导向有关，红枣单价上涨，致使部分耕地转变为园地［31-32］，而水体、建设用地和未利用地的斑块密度数值变化较小，30 a 间没有明显的变化趋势。（2）最大斑块指数（LPI）：从最大斑块指数的数值变化情况来看，未利用地类型的最大斑块指数明显大于其他景观类型，但整体呈现出下降的趋势；此外，园地类型的数值整体呈现出上升的趋势。30 a 间，林草地类型在1994年出现了有史以来的最大值，而2017年的数值最小，出现最大值和最小值的可能原因是垦区最开始不断开荒扩地，后面未利用地的开垦速度降低，耕地的开源主要是林草地的转化。相比之下，水体和建设用地类型的最大斑块指数基本没有变化。（3）面积加权平均分维数（FRAC_AM）：林草地和园地的分维数明显高于其他类型，表明林草地和园地的景观斑块形状最复杂。30 a间，耕地、建设用地和未利用地的分维数都有所增加，而水体的分维数呈现出减少的趋势。（4）聚集度指数（COHESION）：园地类型在1990—1999 年的数值呈现出增加的趋势，说明此期间园地的聚合性增强。

图7 阿拉尔垦区在斑块类型水平上的景观格局指数变化Fig.7 Index change of landscape pattern at patch type level in Alar Reclamation Area

3.3 景观水平上的景观格局分析

30 a 间，研究区在景观水平上的指数变化情况如图8 所示。斑块个数（NP）呈现出波动式上升的趋势，说明区域内的景观碎斑有所增加；景观形状指数（LSI）同斑块个数呈现出相似的变化趋势，说明研究区景观斑块呈现出一定的复杂性；聚集度指数（CONTAG）呈现出明显的减小趋势，说明研究区近年来出现了较多的小斑块，而大型斑块的数量在减少，整体空间分布较为破碎；香农多样性指数（SHDI）整体上处于增加的趋势，说明研究区异质性有所增强。

图8 阿拉尔垦区在景观水平上的景观格局指数变化Fig.8 Index change of landscape pattern at landscape level in Alar Reclamation Area

图9为利用累积距平法进行突变点检验所得的结果，从图9a、图9b、图9d中可以看出，以2005年为界，累积距平呈现出先降后增的趋势；由此说明阿拉尔垦区在1990—2019年间的景观破碎度、复杂度和多样性呈现出了先降低后增加的趋势；图9c的曲线变化情况刚好相反，同样以2005 年为界，在2005年之前曲线呈现出一定的增加趋势，在2005年之后呈现出下降的趋势；由此可以看出，2005 年为阿拉尔垦区的突变年份，区域景观连通性呈先增后降的趋势。研究区从2004 年下半年开始大面积荒地被开垦为农田，2005 年开始种植农作物，农作物种植面积急剧增加，从而植被覆盖度增加，致使阿拉尔垦区的景观破碎度、复杂度和多样性显著增加，最终导致阿拉尔垦区的景观格局发生了明显的突变现象。本研究使用累积距平法对该突变点进行检验，验证其真伪，以保证研究的可靠性。我们将1990—2019年的连续时间序列数据以2005年为界，分为两个子时间序列，结合SPSS软件中的K-S两样本检验法对两子时间序列进行差异性检验，所得结果如表2所示，从表中可以看出，各指数的双尾显著性检验（sig）均低于0.05，说明两个子时间序列的差异性显著，所得结果可靠。

图9 阿拉尔垦区在景观水平上的突变分析Fig.9 Analysis of abrupt change at landscape level in Alar Reclamation Area

表2 双样本K-S检验结果Tab.2 Result of K-S test for two-independent samples

3.4 基于移动窗口的景观格局分析

在对景观指数进行可视化分析中，移动窗口大小的选择很重要，选用不同的移动窗口所得的景观格局图差异显著，经过不断的调试，对比在边长为300 m、600 m、900 m、1200 m、1500 m 和1800 m，共计6 种正方形移动窗下的景观格局结果，最终选用1200 m进行可视化分析，所得结果如图10所示。

图10 基于移动窗口的景观格局分布Fig.10 Landscape pattern distribution based on moving windows

过度的人类活动造成斑块破碎化，斑块密度能够反映出研究区景观破碎化程度，对斑块密度进行分析能够有效反映人类活动强度和土地利用程度。由于研究区中部有塔里木河流，使得靠近水源的区域农业发展更好，景观结构更为稳定，同时人类活动也更剧烈，从图10a可以看出，中间塔里木河区域的斑块密度更高，这很可能和近年来人类活动不断增强，工业用地和农田面积增加，然而阿拉尔垦区的工业废水和农业排水多是往偏远地区排放，导致偏远地区生长出零星的植物组团，植物多样性增加，这些植被并不是连片分布，而是呈现出离散的空间分布趋势，最终偏远地方的斑块数量增加。

同理对香农多样性指数进行空间可视化分析，所得结果如图10b 所示，在中部塔里木河区域香农多样性指数的值会更高一些，整体呈现出以塔里木河为中心向四周不断减弱的趋势。这和中部人类活动更频繁，土地利用结构更稳固有关，周边区域多为未利用地，景观结构单一，开发空间充足。

4 讨论

大多数研究使用的遥感影像数据往往只有3～6景，这种分析只能体现景观格局的大致演变趋势，不能在长时间尺度下得到详细的景观格局演变特征，容易遗漏连续长时间序列中的关键“拐点”信息［33-36］。本文对研究区近30 a 的所有Landsat 遥感影像数据进行了无差异下载，共使用226 景遥感数据参与研究区景观格局演变分析，并将同一年中各月份影像进行NDVI最大值合成，得到周年NDVI最大值合成的多时相影像。这种基于连续长时间序列的多时相遥感影像进行的景观格局分析，不仅能有效的弥补了以往基于3～6 景影像数据分析的不足，充分反映出研究区近30 a来景观格局的细部变化特征，而且能够避免由于不同作物生长旺季不同造成的“错分”现象［37-42］。

此外，在人为干扰影响较大的年份，如2005年，阿拉尔垦区景观格局发生了明显的突变现象，如果仅用少量的遥感影像数据进行景观格局分析，则不能有效捕捉研究区景观格局的细部变化特征。因此，在进行景观格局变化分析时，影像数据的时相选择至关重要，直接关系到研究结果的准确性和可靠性。

在景观格局变化的驱动机制方面：（1）滴灌技术的普及。2005 年以后，滴灌技术逐渐普及，不仅使得原先无法垦殖、地势较高的土地可以垦殖了，而且地势较低的、盐碱重的区域也能垦殖，致使绿洲耕地面积不断增加，未利用地面积不断减少，景观丰富度和异质性逐渐增加。（2）种植业结构的调整。如棉田向林果业调整对绿洲景观产生了一定的影响，使得园地的斑块密度、最大斑块指数、分维数、聚集度指数增加，致使绿洲景观的空间格局分布不均匀，景观聚合性增强。（3）农业技术进步对景观格局的影响。技术的进步（滴灌普及）使得原本少雨、高温、缺水的干旱区域大大提高了地下水的利用率，进而影响到土壤和植被，致使研究区景观丰富度和异质性增加，对绿洲景观的景观格局造成了一定的影响。

景观指数分析的客观性：结合研究区的实际情况以及研究重点，例如，本文对景观斑块数及其相关的斑块密度进行分析时，着重对研究区所占面积比例较大的耕地、林草地和园地进行了景观格局分析，并对其三者之间的相互转化进行了解释说明，而研究区中面积较小的水体和建设用地进行了简要介绍。因此，景观指数分析存在一定的客观性，不同的研究者分析，其结果差异显著。

5 结论

本研究对阿拉尔垦区近30 a间的景观格局变化进行分析研究，主要结论如下：

（1）对多时相面向对象、多时相基于像元和单时相面向对象3 种分类方法的比较分析表明，多时相面向对象方法的分类精度最高，总体精度为96.57%，Kappa 系数为0.95，较后两者分别提高了7%、0.11和12.22%、0.25，有效避免了“椒盐”现象以及植被物候差异造成的“漏分”和“错分”。

（2）从斑块类型水平上看，1990—2019 年阿拉尔垦区耕地景观趋于破碎化，未利用地优势度降低，园地和林草地的景观形状趋于复杂化，水体和建设用地的景观趋于均衡化；从景观水平上看，研究区景观形状趋于复杂化和破碎化，景观连通性降低，景观丰富度和异质性增加；垦区景观格局于2005 年发生突变，景观破碎度、复杂度和多样性显著增加，景观连通性降低。

（3）从空间分布方面看，塔里木河沿岸区域景观指数空间分布较为均匀，人类活动频繁，景观结构完整，景观格局呈现出以塔里木河沿岸区域为核心向偏远区域蔓延的趋势，其中向东北和西南方向的扩张趋势更为明显。