澜沧江下游景观破碎化时空动态及成因分析

刘世梁，安南南，侯笑云，董世魁，赵爽，许经纬

北京师范大学环境学院水环境模拟国家重点实验室，北京 100875

运用遥感影像获取地表植被覆盖现状及其变化信息，对于揭示地表空间分异规律、探讨景观格局动态具有重要的现实意义。目前，在景观格局研究中应用最为广泛的两种方法：一是景观格局指数方法，二是空间自相关分析法（张松林等，2007a，2007b；冯异星等，2010；吴莉等，2013；周天墨等，2013）。景观格局主要是指构成景观生态系统或土地利用/覆被类型的形状、比例和空间配置（胡巍巍等，2008），通常利用景观格局指数进行描述。目前，通过计算多时相不同景观类型的景观指数，分析时间序列景观格局动态的研究较多（王小东，2011；崔世杰，2013）。而利用移动窗口法可以将景观指数进行空间化展示，使景观的空间格局更加清楚，因此成为研究景观空间连续分布的重要方法。移动窗口法是通过设置窗口边长的大小，来计算窗口内所选的景观指数。目前移动窗口法已广泛应用在流域、城市、干旱河谷和梯级水电站等空间梯度格局研究中（张琳琳等，2010；刘琦等，2012；陶宇等，2013；张玲玲等，2014）。

在区域尺度上，通过遥感影像反演植被指数来表征地表覆盖信息的研究也越来越广泛，归一化差植被指数（Normalized Different Vegetation Index，NDVI）（崔林丽等，2010；宋富强等，2011；王晓利等，2013）是较为常用的指数。通过运用空间自相关法分析NDVI的空间聚集性及关联性，来理解景观异质性的研究较为广泛。空间自相关法用于检验具有空间位置的某一要素的观测值是否显著地与其相邻空间点上的观测值之间表现出相对位置的依赖关系（邱炳文等，2007；张松林等，2007a，2007b）。张雪艳等（2009）、马燕飞等（2010）、王晓利等（2013）采用空间自相关指数来分析区域NDVI的空间分布格局及聚集特征。

但是，运用两种方法对比分析景观异质性的研究相对较少。Chao et al.（2014）利用移动窗口法和空间自相关法分析城市景观破碎化程度，认为两种方法在反映景观格局和空间异质性方面有着良好的可比性，但未利用NDVI进行指示性分析，同时基于多源数据的空间分析仍需要加强。本研究基于遥感数据和景观指数，对比分析两种方法在提供景观结构特征信息方面的可靠性，以此来分析景观的空间异质性。

另外，NDVI的变化与土地利用/土地覆盖之间关系及其驱动因子分析也受到广泛关注（毛德华等，2011；梁尧钦等，2012）。大量研究表明，植被分布随地形梯度的变化可以揭示地形因子对景观演变的影响程度（张秋玲等，2007；陈利顶等，2008）。因此，本研究通过综合选择多个地形因子，分析NDVI变化与外部环境的相互关系，为了解景观格局的空间分布规律奠定基础。

2015年，景洪市林业用地面积占全市土地总面积的 68.3%（李江龙，2015），随着工业化、城镇化等人类活动的加剧，该市林地生态系统监控与管理问题更加突出。本研究结合多时段遥感影像资料，通过提取连续型数据（NDVI）、离散型数据（景观类型）和地形因子等信息，对比移动窗口法和空间自相关法定性分析区域景观分布格局和 NDVI的空间关联特征，了解植被景观异质性，探寻可稳定反映景观异质性的方法，以期从案例研究角度获得林业资源管理的科学依据。

1 数据来源及研究方法

1.1 研究区概况

景洪市位于云南省西南部，地处横断山系纵谷区南端，澜沧江大断裂带两侧，最高海拔约为2400 m，澜沧江和流沙河为流经该区的主要河流。该地区属于热带、亚热带季风气候（刘世梁等，2006），长夏无冬，年均气温为 18.6～21.9 ℃，年均降水量为1200～1700 mm，季节分配极不均匀，一年之内有明显的旱季和雨季。植被类型主要包括热带雨林、热带季雨林、南亚热带常绿阔叶林、南亚热带针阔叶混交林、竹木混交林、灌木林和草丛等。该地区是中国进入东南亚各国的主要通道和枢纽，交通建设发展迅速，人口密度增长快速。本研究区域为60 km×60 km的正方形，主要包括景洪市的中部和勐海县的东部地区，地跨 21°40′～22°20′N，100°30′～101°10′E。以重点旅游干线和重点集镇为骨架（耿红生，2011）。2000—2010年，研究区内林地资源被大量占用，建设用地快速扩张，景观格局发生了明显变化。

1.2 数据源

采用多时相Landsat遥感数据分类，获得所需要的景观类型信息，根据2000年、2005年和2010年3个时段的TM影像（影像分辨率为30 m），进行人工目视判读与监督分类，并结合2010年实地调研验证，获取3个时期研究区景观类型图，其影像的分类精度达 85%。本文参照 IGBP土地利用/土地覆被分类标准（张景华等，2011），将该区景观类型划分为常绿针叶林、常绿阔叶林、人工林、城市绿地、草地、水域、耕地、建筑用地和未利用地（图1）。

图1 研究区2000年、2005年和2010年景观类型图Fig. 1 The landscape type maps of study area in 2000, 2005 and 2010

利用2000—2010年的SPOT VEGETATION逐旬NDVI数据（空间分辨率为1 km），综合考虑研究区域数据源及研究数据的有效性，在研究区域进行重采样，生成间隔为1 km的点图层，叠加3期NDVI栅格图层提取每个点上的NDVI值。同时，基于DEM图层，提取海拔、坡度、坡向、地形起伏度及地表粗糙度等5种地形因子。

1.3 研究方法

1.3.1 移动窗口法

利用 Fragstats 4.2选取反映景观水平的景观指数，选择边长为300 m矩形窗口进行计算，通过对窗口内选中的景观特征进行计算，输出数据为对应所选景观指数的栅格图。由于景观指数类型繁杂，大部分指数在指示景观格局特征时，存在灵敏度不高、指示含义重叠等问题（何鹏等，2009）。根据研究的实际需要，选取斑块密度（PD）、边缘密度（ED）、最大斑块指数（LPI）、景观形状指数（LSI）、香农多样性指数（SHDI）、蔓延度指数（CONTAG），拟从破碎度、形状和多样性角度探讨研究区景观格局的空间变化特征（何鹏等，

2009；张立平等，2014）。

1.3.2 空间自相关法

空间自相关法通常采用全局和局部指标进行度量。由于全局指标在一定程度上会掩盖局部状态的不稳定性，本文采用局部指标来分析研究区域内NDVI在不同空间单元之间的自相关程度。用以表示空间自相关的指数和方法很多，其中 Moran’s I指数是在分析自相关性中应用最广泛的空间自相关判断指标。局部Moran’s I指数计算公式如下：

式中，xi和xj分别是变量x在配对空间单元i和j上的取值；x是变量x的平均值；n是空间单元总数；wij是相邻权重。本文空间权重矩阵采用的是K-nearest-neighbors，邻域数目为4。局部Moran’s I指数的取值范围为-1～1，正值表示该空间单元与邻近单元的属性值相似（“高—高”或“低—低”），负值表示该空间单元与邻近单元的属性值不相似（“高—低”或“低—高”）（李慧等，2011），该值为零时，该空间单元呈独立随机分布。对于空间是否存在自相关性，常采用z检验（Southworth et al.，2004）进行标准化。检验公式如下：

式中，E(I)为期望值；Var(I)为方差。

本文采用多时相SPOT-VEGETATION NDVI数据，通过对多个时段的NDVI与地形因子的空间自相关性进行比较，进一步解释NDVI变化与地形因子之间的生态意义，这使得NDVI空间自相关研究更有价值。局部 Moran’s I的计算通常在 Geoda 0.9.5-i软件中进行。

2 结果与分析

2.1 移动窗口法分析景观格局的空间变化

2000—2010年，各景观指数在变化不明显（图2）。其中，CONTAG和LPI变化明显的区域主要分布在研究区域的中南部，该地区分布着广泛的常绿阔叶林，在2000—2010年间，CONTAG和LPI值不断增大，说明该地区的各类型的景观斑块离散程度较低，连接度较好。

以 2010年 CONTAG指数为例，研究区中心CONTAG值较低，LPI值也较低，ED和LSI值都较高，说明研究区中心的景观破碎化程度较高，异质性高。主要是因为景洪市和景洪水电站距离研究区中心较近，工业化、城镇化速度加快，人类干扰活动较多，土地利用随之变化，斑块类型更为丰富。在研究区中心的四周和北侧中部区域，CONTAG值较高的区域，LPI值也较高，ED、LSI、SHDI及PD值则较低，说明该地区的斑块较大且完整，连通性较好，景观边缘较规则且景观多样性较高。其中，研究区中心的四周出现这种现象的原因主要是该地区主要的景观类型是人工林，在规范的管理模式下形成了较规则的斑块。而在研究区域北侧中部出现该现象主要是由于该地区景观类型大多以常绿阔叶林、常绿针叶林为主，这些林地主要是天然林，林地形成历史较长久，加上受到的人类活动干扰较少，因此景观破碎化程度较低。在研究区的西侧中部地区，CONTAG值较低，LPI值也较低，ED和LSI值都较高，说明该地区的景观斑块分散，连通性较差，其主要原因是该地区景观类型较为多样，主要包括灌丛、常绿阔叶林和常绿针叶林等，斑块类型较为丰富。

2.2 研究区不同景观类型的NDVI年际变化

利用 2000年、2005年和 2010年的年 NDVI最大值的平均值绘制时间变化柱状图（图3）。整体上，研究区不同景观类型历年NDVI值均在0.6以上。2000—2010年间，各景观类型的植被覆盖趋于上升态势。其中，未利用地的植被覆盖呈明显上升趋势。

2.3 研究区时间序列NDVI整体空间自相关性

图2 基于移动窗口法计算的2000年、2005年和2010年景观指数栅格图Fig. 2 Raster grid mapsof landscape indices based on moving window method in 2000, 2005 and 2010

分析研究区3期年最大NDVI的空间自相关指数（图 4）。2000年、2005年和 2010年的年最大NDVI的局部 Moran’s I指数分别为 0.781、0.791和0.856，说明NDVI均呈现出显著的空间正自相关，即该区域植被覆盖都呈现集聚状态，在空间上具有较强的关联性。研究期间，NDVI的局部Moran’s I指数呈逐年升高趋势，说明该区域NDVI的空间集聚格局随着时间推进变得更加显著。

进一步分析NDVI的局部空间自相关特征，以更好地探索NDVI的局部空间聚集模式和规律（图5）。其中，以2010年为例，结合景观类型图，NDVI呈“高－高”自相关的地区大多为常绿阔叶林分布区，这些地区植被覆盖总体较好，NDVI高值区集聚模式较明显；NDVI呈“低—低”自相关的地区大多为建筑用地、耕地分布广泛的地区，这些地区是人类活动聚集的地区，景洪市及景洪水电站均位于建筑用地区域，景洪市快速扩张及景洪水电站对周边植被的破坏，都导致这些地区生态脆弱而敏感，植被覆盖较少，而耕地呈现出“低—低”自相关模式的原因主要为该地区人类活动剧烈，农作物多为“一年两熟”型，导致该地区的植被覆盖不稳定，空间集聚性较差。NDVI呈“低—高”和“高—低”聚集模式的区域在研究区内分布较少。

图3 2000—2010年间不同景观类型NDVI变化Fig. 3 NDVI changes of different landscape types from 2000 to 2010

图4 NDVI的局部空间自相关Moran’s I散点图Fig. 4 Moran’s I scatterplot of local spatial autocorrelation of NDVI

图5 不同时期NDVI空间关联局部指标（LISA）分布图Fig. 5 LISA distribution map of NDVI in different years

2.4 移动窗口法和空间自相关对比

分别统计各景观类型的CONTAG和LSI值，以进一步对比分析局部Moran’s I指数和景观指数在反映景观破碎化时的可靠性。景观破碎化程度较高时，表现为低CONTAG值，高LSI值。由于草地在2005年后几乎不存在，本文忽略不计。以2010年为例（表 1，表 2），常绿阔叶林的局部Moran’s I指数平均值为0.53，表明常绿阔叶林空间分布较集聚，空间关联性较好。同时，常绿阔叶林的蔓延度指数CONTAG为51.96，景观形状指数 LSI为1.35，这也表明常绿阔叶林的斑块较少且完整，破碎化不明显。裸地和建筑用地的局部 Moran’s I指数平均值分别为-3.44和-3.54，说明未利用地斑块分布较为分散，空间聚集性差。同时，未利用地的蔓延度指数 CONTAG分别为24.15和39.24，景观形状指数LSI分别为1.62和1.41，在一定程度上表明未利用地的景观破碎化较为严重。

表1 2010年不同景观类型的局部Moran’s I指数描述性统计Table 1 Descriptive statistics of local Moran’s I index by landscape types in 2010

2.5 景观空间分布与地形因子之间的关系

地形在一定程度上决定着景观类型的空间分布格局，不同景观类型的空间聚集和异常现象都与地形密切相关。地形因子与土地利用聚集或异常特征的出现关系密切，进而影响景观类型空间自相关格局（谷建立等，2012）。提取了海拔、坡度、坡向、地形起伏度和地面粗糙度等地形因子（表3），分析了NDVI与地形因子之间的关系。海拔、坡度、坡向、地形起伏度和地面粗糙度与三期景观类型的 NDVI的 Moran’s I指数分别为0.262、0.271、-0.018、0.302和 0.186。该研究区的 NDVI与海拔、坡度、地形起伏度和地面粗糙度的Moran’s I指数达到了0.18以上。此外，NDVI与这4种地形因子呈正相关关系。由此可以说明，本研究区植被覆盖度一定程度上由海拔、坡度、形起伏度和地面粗糙度决定。

表3 不同时期地形因子与NDVI的Moran’s I统计值Table 3 Statistics of Moran’s I valuesbetween NDVI and topographical factors in different years

3 讨论

在景观类型和NDVI数据基础上，结合遥感和GIS等技术，利用移动窗口法和空间自相关法，研究时间序列景观格局的变化规律，并比较了两种方法在指示景观破碎化方面的可比性，为进一步探索指示景观破碎化的合适方法提供了依据。同时，通过对景观格局与地形因子的相关性进行分析进一步解释了环境影响景观格局的生态意义，使得景观异质性研究更有价值。

分析多个景观指数栅格图，发现常绿阔叶林和未利用地等景观类型，其局部自相关指数高的区域，CONTAG值也高，而LSI值较低；反之亦然。这种现象表明，移动窗口法计算的景观指数和空间自相关指数在指示景观破碎化时具有可比性。2000—2010年间，研究区域中心偏东南方向的人工林聚集区域CONTAG指数有逐年变大的趋势，而LSI、ED等指数则有逐年变小的趋势，这表明，在这期间该区域的植被景观受到良好的保护，这和当地政府采取的“退耕还林”的政策密切相关，退耕还林后该地区植被类型的组成和结构发生了变化。

表2 2010年不同景观类型的景观指数描述统计Table 2 Descriptive statistics of landscape indices based on landscape types in 2010

4 结论

（1）研究区面积和移动窗口的大小，将会影响景观破碎化的计算结果。在采用移动窗口法计算时分别尝试150、200、300、500、1000和1500 m的边长，通过对比分析发现，300 m的窗口边长能保留梯度特征又不至于使景观指数出现较大波动，可以通过景观指标的变动特征真实反映空间格局的变化。

（2）统计对比景观指数和空间自相关指数后发现，两种方法在指示景观破碎化上均具有良好的适用性。并且发现，区域内植被的分布一定程度上由海拔、坡度、地形起伏度和地面粗糙度所决定。但影响植被景观变化的因素比较复杂，在今后的工作中可以结合土壤、水文、生物等其他自然地理要素和人为因素，分析植被景观的多尺度变化特征，从而进行合理的景观布局。

（3）景观生态过程是动态的、发展延续的过程，将静态的格局分析赋予动态变化的属性，是把格局分析与生态过程研究结合起来的重要分析方法。今后，将结合时间序列的景观数据和社会经济数据等，建立过程影响因子动态变化图谱，从而构建新的景观动态变化格局，以分析其变化特征并实行空间格局的优化。