城市景观空间格局是包括干扰在内的各种生态过程在不同尺度上作用的结果，其作为一种较为复杂的镶嵌体结构，与人类的生存环境密切相关。景观格局动态变化分析是当前研究城市化剧烈发展对生物生境影响的主要方法^[1-3]，景观格局指数以其能够量化反映景观各要素的结构和空间配置得到广泛的应用^[4-6]。然而，在刻画具体城市景观变化及其生态过程时，需要选择适宜的分析尺度，即一定幅度下如何选择与生态过程关联的景观粒度，景观指数的选择须对粒度的变化具有很好的敏感性和指示性^[7-8]。

景观格局指数虽高度浓缩了景观格局信息，但不能深入反映景观要素间的内在联系及其生态过程，有学者指出连续变量用景观指数方法分析是目前景观格局研究中的缺陷，但此方面的研究对于植被的退化与恢复有重要的意义^[9-10]。特别是近20年来城市快速扩张对景观格局连通性产生影响所涉及连续性变量的分析，逐渐引起学者们的关注。连续变量的空间自相关分析可以有效地反映特定景观要素在功能或生态过程上的关联，是衡量和解释生态过程的重要指标^[11-13]。

目前，已有研究发现，梯度或连续性的模型方法能够准确地描述研究对象的景观异质性^{[2, 14]}。然而对于连续性的空间自相关分析方法和离散的斑块镶嵌方法进行关联性的研究还不多见。FAN et al^[9]提出仅基于粗略的土地覆被分类再利用传统的景观格局指数计算可能存在一定的精度问题，因此利用局部空间自相关指数和景观指数相结合的方法测定城市景观破碎化程度。KOWE et al^[2]采用多期美国陆地卫星(LANDSAT)影像数据，利用景观指数与表征局部空间自相关的关联性指标(local indicators of spatial association, LISA)图，对津巴布韦哈拉雷大都市区的植被斑块破碎化模式的时空变异特征进行研究。LI et al^[14]基于长时间序列的LANDSAT影像，应用空间自相关指数Moran′s I与景观指数结合的方法对我国官厅水库表面水体的缩胀进行研究。多数研究还发现自然植被在城市中为动植物提供自然栖息地，具有多种复合生态效益，是城市景观的重要组成部分^[2]。而归一化植被指数(normalized difference vegetation index，NDVI)是植物生长状态及植被空间分布密度的最佳指示因子^[15]，可作为典型的空间连续分布变量^[16]进行自相关分析。

因此，本文基于两期遥感数据，在最佳粒度分析的基础上，运用网格移动窗口法结合景观指数、空间自相关分析对常州市2010和2015年的景观格局和植被空间连通度特征的动态变化进行研究，从景观指数和空间自相关两个视角全面探究快速城市化发展中典型城市的景观格局空间变化特征，以及探讨两种方法在表达特定类型景观定量化之间的关联性与互补性。

常州市位于江苏省南部(119°08′~120°12′E，31°09′~32°04′N)，地处长江三角洲中心地带，行政辖区面积4 375 km²。市域南部为天目山余脉，西为茅山山脉，北为宁镇山尾部，中部和东部为广阔的平原、圩区。境内地势东北略低，西南略高，属高沙平原地貌，山丘平圩兼有^[17]。地属北亚热带季风气候区，气候温和湿润(年平均气温16.3 ℃，年平均降水量1 068.9 mm)，夏季炎热多雨，春末夏初时多有梅雨发生，冬季阴冷潮湿。2013年以来相继实施生态景观廊道建设、森林生态屏障建设等重点工程，2016年，被评为国家森林城市，2019年，绿化覆盖率达43.25%^[18]。常州市目前常住人口约570万，近年来城市土地利用开发强度较大，景观格局变化强烈。

选用2010年8月19日美国陆地卫星专题制图仪(thematic mapper，TM)所获取的多波段影像和2015年8月11日的LANDSAT-8 OLI影像，轨道号为119/038、120/038，共计4景影像作为数据源，来源于中国科学院计算机网络信息中心地理空间数据云平台(http://www.gscloud.cn/)。利用ENVI 5.0软件对数据进行波段合成、几何校正和影像拼接和掩膜裁剪等预处理，生成研究区域原始遥感影像图。采用最大似然法对影像进行监督分类，提取研究区土地利用类型特征，参照我国土地利用分类体系^[19]将研究区划分为耕地、交通用地、居民点及工矿用地、林地、水域及未利用地6类，在对地物分类错误进行修正后，经实地验证与比对，两期影像的解译整体精度达85%以上，均满足本研究分析要求。两期景观类型如图 1所示。

	图 1 2010和2015年常州市土地利用景观格局 Fig. 1 Changes in the landscape pattern of land use in Changzhou City in 2010 and 2015
图选项

本研究所使用的NDVI数据来源于中国科学院资源环境云平台提供的中国年度植被指数(http://www.resdc.cn/DataSearch.aspx)NDVI时序数据，空间分辨率为30 m，是基于SPOT卫星VEGETATION传感器获取的长时间序列NDVI影像，采用最大值合成法后裁剪生成2010和2015年研究区的NDVI数据。

为获取研究区的最佳分析粒度，采用重采样的方法，将2010和2015年的景观类型图分别进行像元的聚合，栅格单元的赋值采用中心属性值法。借鉴以往文献中关于相当规模区域粒度的设置^[20-22]，以20~400 m为粒度转换的考察区间范围，其中120 m以下每10 m为粒度增幅，以上为50 m为增幅，共生成16个不同粒度等级的栅格图，从而分析各斑块类型景观指数随粒度的变化规律。

基于FRAGSTATS软件，参考相关文献^{[8-9, 23]}并考虑常州市景观基底和廊道等特征的实际情况，从景观斑块的分布特征、结构与形状复杂性和邻近性等方面选取下列景观格局指数：包括斑块密度(patch density, PD)、景观形状指数(landscape shape index, LSI)、平均斑块分维数(mean patch fractal dimension index, FRAC_MN)、最大斑块指数(largest patch index, LPI)、景观类型比例(percentage of landscape, PLAND)和聚合指数(aggregation index, AI)等，统计分析后得到不同时期和空间粒度下景观格局指数的变化规律，分析获得各指数的粒度拐点，根据尺度域确定最佳分析粒度。同时，考虑到矢量数据在栅格转化过程中可能存在的精度损失，本文基于两个时期景观斑块的栅格面积变化，评估不同粒度水平下的整体精度损失^[24]，公式如下所示。

式中: L_i为第i类斑块类型面积损失率；A_i表示景观中第i类斑块类型的栅格面积；A_bi表示第i类斑块类型所基于的矢量面积。S_i表示区域内各斑块类型面积变化的差异指数，该值越大，空间粒度变化引起的信息损失就越大；n为研究区内斑块类型的数量。

在最佳分析粒度基础上，为从景观指数与空间自相关两个角度对比分析研究区的景观空间变化特征，以植被这一连续变量为例，基于fishnet网格移动窗口法，用两种方法提取植被类型的空间分布特征：(1)景观格局指数法测定移动网格窗口中植被覆盖类型的景观指数信息；(2)空间自相关方法计算反映连续植被信息的NDVI局部空间自相关指数^[9]。

参考前人研究与城市面积大小^{[3, 25]}，以5 km×5 km网格为基本研究单元，将研究区共划分为221个网格，基于移动窗口法利用FRAGSTATS 4.2软件批量计算不同时期各网格的景观格局指数^[26-27]，插值得到植被景观指数的栅格图像；应用Arcgis 10.0中ArcToolbox工具箱中的空间分析(spatial analysis)工具中的分区统计(zonal statistics as table)计算每个网格空间的NDVI指数平均值，两种方法的投影系统需一致。利用Geoda软件计算NDVI值的局部空间自相关指数(Moran′s I)。计算公式^[28]如下。

式中: x_i和x_j分别为景观要素x在空间单元i和j中的观测值，x为x的平均值，W_ij为相应权重，n为空间单元个数。I的值域为[-1, 1]，-1表示空间负相关，1表明空间正相关，而0表示随机分布，无空间自相关性。

运用SPSS 18.0对数据进行Kolmogorov-Simonov(K-S)正态分布检验和Pearson统计分析。

2010—2015年常州市在斑块类型水平上的土地利用景观格局指数在20~400 m空间粒度变化下具有显著的粒度效应，且不同斑块类型的粒度效应呈现出差异(图 2)。

	图 2 2010和2015年各斑块类型的景观指数对粒度变化的响应 Fig. 2 Response of landscape indices to grain size change at patch type level in 2010 and 2015
图选项

各斑块类型的斑块密度(PD)均随着粒度的增大明显减小，20~30 m范围内各斑块类型保持一个较稳定的高值，30 m左右出现尺度的转折点，粒度超过150 m以后，各斑块类型的PD值差异明显减小。PD值表征了斑块的破碎度与空间异质性程度，图 2(a)和图 2(b)显示居民点及工矿用地的PD值明显高于其他斑块类型，而林地、水域的PD值相对较低，说明建设用地的破碎度相对较大，而生态用地则相对团聚；但2010—2015年建设用地的面积增加，但PD值减少，可能与“两横一纵”城镇聚合发展等政策相关；而林地、耕地的PD值增大，则与旅游开发、建设用地扩张等影响，生态用地被分割，使得景观破碎化。

景观形状指数(LSI)是结合景观面积对景观总边缘长度或边缘密度的标准化度量，同时也是斑块聚散程度的量度，LSI值越大斑块越离散。图 2(c)和图 2(d)可见LSI随着粒度增大而总体呈指数型下降趋势。根据粒度效应曲线，景观形状指数在20~30 m范围内出现波峰值，30~150 m范围内，LSI值平缓下降，说明在粒度增加过程中各景观类型的边缘由锯齿状逐渐变为平滑；相较于斑块几何形状简单且聚集分布的林地和水域等类型，居民地及工矿用地由于农村居民点及部分工矿用地的离散分布，小斑块较多，斑块破碎化程度较高，形状复杂，其LSI值相对较大。从时间上看，LSI随粒度变化的平均值，居民及工矿用地(113.4~124.6)、耕地(107.8~116.0)、林地(43.1~55.7)均有增大趋势，这与城市化背景下建设用地的快速增加对生态系统的扰动有关。

平均斑块分维数(FRAC_MN)是衡量景观要素镶嵌结构复杂性的指标，分维数的理论值在1~2，值越大斑块结构越复杂。图 2(e)和图 2(f)分维数整体随着粒度变化呈非线性下降趋势，可分为3个阶段，在20~30 m范围内保持短暂的波峰值；20~90 m范围内急剧下降；90~400 m趋于平缓，整体呈下降态势，景观类型的斑块形状趋于规则。其中，交通用地的FRAC_MN在20~90 m范围内明显高于其他用地类型，且下降幅度最大(0.13)，这是由于交通用地呈带状分布且贯穿于整个城市，当采样粒度达到并超过斑块的宽度时(约30 m分界点)，带状斑块在重采样过程中断裂，最终在粒度的增大过程中被其他斑块吞并。90~400 m范围内交通用地的分维数逐渐低于其他用地类型，其原因可能是随着粒度增加，稍大交通用地斑块形状趋于规则与平滑，而较小斑块被吞并现象加剧所致。

最大斑块指数(LPI)是反映景观的优势类型和丰度等特征的指标，其值的变化可改变干扰的强度和频率，反映人类活动的方向和强弱。随着粒度增大，研究区各景观斑块类型LPI值变化较为平缓，但其中耕地类型随粒度变化呈阶梯式的增长方式较为明显，转折的粒度主要是30、70~90 m附近，这与粒度变化中大斑块对小斑块的兼并有关。从时间上看，2010—2015年耕地整体面积明显减少，由于受到“拓展南北，提升中心”等城镇发展政策和“三纵四横”等交通布局影响，城市发展中干扰了部分耕地的连片格局，使其最大斑块效应有所减弱。

总体看来，PD、LSI、FRAC_MN和LPI均表现出明显的粒度效应，前3个指数的粒度效应整体表现为曲线单调下降；最大斑块指数随粒度变化对优势类型(耕地)影响较为明显，具体影响表现为阶梯状的上升趋势。从粒度变化的不同阶段来看，20~30 m范围内表现出明显的梯台效应，30~90 m指数具有明显的上升或下降态势；90~400 m范围内各斑块类型的指数逐渐稳定，受粒度变化的影响逐渐减小。

研究景观的最佳粒度，需要在综合分析各个景观指数的粒度拐点的基础上划分最佳粒度尺度域。景观指数对粒度有明显响应时，应根据其第一尺度域确定最佳分析粒度，且选取第一尺度域内中等偏大的粒度，一方面能够很好地反映原图像的特征信息，另一方面可以减少计算量^[20]。因此，综合所选景观指数随粒度变化趋势，分析得到第一尺度域为20~30 m。常州市的面积和人口属于中等城市规模，主干道宽度为30~40 m，居民住宅和工矿用地平均长度30~50 m，市内的一些小型游憩林、生态公益林、行道树的宽度等也接近这个阈值范围，在粒度增加到与主干道宽度一致(约30 m)时，斑块易断裂，引起某些景观格局指数的突变；当粒度超过此阈值，则会被合并，其他类型相近大小的斑块也是如此。结合研究区的实际情况，以及在矢量数据转成栅格数据的过程中的精度损失变化情况(图 3)，最终选择30 m作为2010和2015年常州市景观格局的最佳分析粒度，为后面从景观指数和空间自相关两个角度深入分析景观格局空间分布特征作铺垫。

	图 3 2010和2015年不同栅格尺寸下景观面积变化精度 Fig. 3 Changes in landscape area accuracy under different grid sizes in 2010 and 2015
图选项

在最佳粒度的基础上，从景观指数和空间自相关两个视角对研究区的植被斑块进行深入分析，并探讨景观的连通性空间分布特征。

基于网络移动窗口法对6个具有代表性的景观指数(PLAND、PD、LPI、LSI、FRAC_MN、AI)进行分析，考虑到篇幅限制及指数之间含义的冗余性，图 4中展现了植被斑块类型4个景观指数的空间异质性分布特征，两期影像中植被斑块类型景观指数的分布规律主要表现为：(1)斑块密度(PD)和景观形状指数(LSI)在市域有3个热点高值区，分别为金坛区与武进区的交界处、武进区东南部、溧阳市中部区域，这些地区以林地和耕地为主，同时镶嵌着城镇居民用地，城镇化的发展阻断了自然植被的连通性，导致景观类型多样化、破碎度增大。(2)最大斑块指数(LPI)反映了植被类型斑块在溧阳市南部出现高值，这是由于天目山余脉的森林植被覆盖度较高，林地呈团聚状成片分布；聚合指数(AI)表现为南部高值聚合分布更加明显。

	图 4 2010和2015年植被斑块类型景观格局指数的空间分布 Fig. 4 Spatial distribution of landscape pattern indices of vegetation patches in 2010 and 2015
图选项

从时间上看，2010—2015年武进区西太湖区域，因中心城市建设用地急剧扩张，该区域出现了PD和LSI由低值区向高值热点区发展，表现出破碎度加剧的形势。而在常州市森林城市建设背景下，南部天目湖国家森林公园等自然资源得到保护，使得植被类型LPI和AI高值范围扩大稍往北延伸。

基于网格移动窗口提取2010和2015年的NDVI值，计算局部空间自相关指数Moran′s I，结果如图 5和图 6所示。两年局部空间Moran′s I值分别为0.74和0.67，从图 5中可以看出研究区NDVI在两个时期均表现出显著的空间正相关，即整体上植被覆盖呈聚集格局，但2015年的空间自相关程度明显低于2010年，这与人类活动逐步增强，在经济驱动下建设用地的快速扩张有关；住宅用地、产业用地及交通用地等逐步侵占耕地、林地和水域边缘，加剧了生境破碎化，使得植被等生态用地的连通性降低。

	图 5 2010和2015年常州市NDVI局部空间自相关Moran′s I散点图 Fig. 5 Moran′s I scatter plot of local spatial autocorrelation of NDVI in 2010 and 2015
图选项

	图 6 2010和2015年常州市NDVI的LISA集聚图 Fig. 6 LISA cluster map of NDVI in 2010 and 2015
图选项

两个时期NDVI指数的Moran′s I空间自相关特征如图 6所示，图 6(a)显示高-高聚集区主要分布在溧阳市南部，以林地覆盖为主，另外在金坛区的北部也零星分布着较强的空间正相关性区域；而新北区的南部、天宁区的西部和钟楼区以及武进区的中部表现出低值区集聚特征。图 6(b)与图 6(a)的区别在于金坛北部耕地区域高-高集聚效应消失；低-低自相关聚集区向武进滆湖区域蔓延。NDVI的高-高集聚区表明“十二五”时期《常州市重要生态功能保护区保护方案》和“一核八园”等政策的实施，以及在森林城市建设背景下生态功能区的保护措施，使得天目山余脉区域植被保护与恢复相对较好，NDVI空间区域之间彼此正向影响，呈高值团聚分布；而低-低自相关结果表明这类区域受到人类活动的影响，造成植被高度分散，连通性低，导致相互之间的关联性较小，形成低-低集聚区。两期莫兰指数的空间分布差异主要由于经济快速发展，城市化过程中西北部茅山风景区旅游业的发展带动交通系统日渐发达，导致城镇周边地区的耕地、林地被占用开发为建设用地，2010—2015年耕地缩减面积达7.9%，使得金坛区北部高-高集聚效应消失；而武进区滆湖等水域边缘的林地与绿地也随着中心城区建设用地规模扩张有一定程度的减少趋势，出现新的低-低集聚区。

因连续变量NDVI指数的高低值空间分布与林地植被相近，因此将林地植被类型的景观指数与NDVI的Moran′s I指数进行对比分析，基于同样的网格区域，分别计算2010和2015年两个时期景观指数与自相关指数，相关性结果如表 1所示。两年的植被Moran′s I与PLAND、LPI和AI均呈显著的正相关关系，2015年的正向相关性较2010年更为显著(P < 0.01)。而植被Moran′s I值与PD和LSI呈负相关关系，其中2015年植被Moran′s I指数增加时，PD、LSI指数极显著降低(P < 0.01)。植被空间自相关指数与表征斑块破碎度的指标呈显著负相关，而与表征斑块聚合与连通度的指标呈显著正相关关系，与实际情况相符。

不同之处在于，景观指数空间分布仅是镶嵌斑块的景观指数值在空间上的插值表现，空间自相关则对植被连续变量的空间关联性进行了定量化的表达，两者具有一定的关联性且互为补充。综合分析表明，景观指数与空间自相关方法相结合为综合全面分析城市景观格局的斑块镶嵌特征与空间自相关特征提供了有效的分析模式。

不同时期的景观类型随粒度变化分析可揭示景观格局的时空变化趋势。本文选取的各类景观指数，如PD、FRAC_MN和LPI(表征景观结构、配置和破碎度)，LSI(表征形态特征)等，对于不同的斑块类型具有明显的尺度效应^[29-30]。但各类景观指数随粒度变化响应的敏感程度与变化趋势有所差别，其中居民点及工矿用地属于建设用地，其斑块密度较大，斑块复杂度与破碎化程度相对其它用地类型较高，而交通用地为带状“廊道景观”，对粒度的响应较为敏感；耕地具有一定的景观“基底效应”，随粒度增大呈阶梯状上升，成片的耕地斑块逐渐兼并小斑块，使得粒度曲线出现明显的阶梯式转折点。2010—2015年，林地和耕地的景观形状指数LSI明显上升，说明受到城市发展过程中旅游开发和建设用地扩张等因素的影响，此类生态用地被分割，导致林地、耕地等斑块形状趋于复杂和破碎化，导致植被景观的连通度降低。

粒度效应研究表明，空间粒度为30 m时是景观指数第一尺度域的转折点，是确定粒度大小，进行景观格局分析的适宜尺度。该粒度分界点的面积转换精度损失最小，几乎保留了原有景观完整的格局信息^[24]。而格局与过程之间的相互影响具有尺度依赖性，只有当测量尺度和本征尺度相符时，格局与过程才能被可靠地揭示^[7]。因此，本文在对景观指数的粒度效应研究中，从时间和空间两个维度分析了各地类斑块的粒度效应，得出的适宜尺度为景观格局指数的空间分析奠定了基础^[8]。

与已有的研究相比，本文侧重于从景观的数量特征和空间自相关特征两个角度探讨快速城市化发展中研究区的景观动态发展规律。景观格局指数作为传统指标分析方法的代表，长期以来一直是理解景观格局和生态过程的主要方法^[31]，但景观指数无法描述土地覆盖中连续变量的关联性特征，分析过程可能会导致重要信息的丢失^[9]。因此，除了景观格局指数分析方法外，引入空间自相关方法来综合处理景观类型中的地类变量，其中植被斑块在城市景观中作为典型的连续变量具有空间自相关性和异质性特征，利用景观指数与空间自相关相结合的方法对其进行重点分析^{[2, 12]}，得出研究区生态用地的连通性变化规律，主要表现为高-高集聚区分布在市域南部森林植被集中区域；而低-低集聚区分布在东北部市区范围，连通度降低且破碎度较高等结果，需引起重视。

随着空间自相关Moran′s I增加，斑块密度PD和景观形状指数LSI显著降低，而斑块类型占景观面积比例PLAND、最大斑块指数LPI和聚合指数AI则显著增加(P < 0.05)，说明Moran′s I指数较大时，景观的破碎化程度与异质性程度较低，植被的连通性较好；反之，景观异质性高、斑块形状复杂，破碎化程度高，植被的连通度较差。此结论与董玉红等^[11]、FAN et al^[9]研究结果吻合。当然，相关系数的大小还与所采用的移动窗口的大小、遥感影像的精度与NDVI的空间分辨率等因素有关，因此，在最佳粒度情形下如何选用适当的幅度尺度来对景观格局进行研究，还有待进一步讨论与分析。

空间自相关分析是检验相邻空间点上观测值关联性的有效方法，能够对各类景观的结构关系进行定量化，在一定程度上弥补了景观指数分析的不足^{[9, 28]}。但研究中仅用了两期遥感数据，今后的研究中将选用多源、多时相的数据，通过对不同分辨率(尺度)、多时期的景观进行空间自相关分析，并进一步利用各种环境因子解释这一变化，深层次探讨引起景观空间自相关性变化的原因会更有价值。

2010—2015年各地类的景观指数随粒度的变化趋势揭示了研究区整体景观格局的时空动态发展规律。从空间上看，耕地随粒度体现出阶梯状的基底效应；从时间上看，建设用地扩张的同时斑块密度减小，体现了集中连片的“聚合”趋势，林地和耕地等生态用地斑块形状趋于复杂，局部连通度降低。空间粒度30 m是研究区景观格局的适宜分析尺度。

破碎化与连通度的空间分布特征通过基于移动窗口的景观指数空间插值和自相关分析获得，两年植被斑块的Moran′s I指数分别为0.74和0.67，说明其空间自相关性程度降低；主要表现为城市扩张使得金坛区北部耕地区域高-高集聚效应消失；低-低自相关聚集区向武进滆湖区域蔓延，该区东南部的景观形状指数LSI明显增加，用地类别趋于复杂，需要预防景观破碎化的发生；低-低集聚区(人类活动密集的市区)建议通过保护与增设市区行道树、街头游园绿地和古运河廊道景观绿化等方式，改善植被高度分散的现状，提升市区植被的连通度，保障城市生态流的畅通，进而实现城市的可持续发展。从整体规划的角度，各景观类型的破碎化程度加剧，应控制建设用地总量，提高集约用地水平，严格保护林地、优质耕地，提高植被斑块的连通度与系统整体的稳定性。

此外，景观指数与空间自相关方法相结合，从两个角度为综合全面分析城市景观格局斑块镶嵌特征与空间自相关规律提供了有效的分析模式。研究结果有助于加深城市复杂的异质性景观格局下对于连续生态过程的认识，对于预防城市化进程中的土地退化和破碎化，提出植被恢复策略，提升人居环境具有一定的指导意义。