干扰是自然界中无时无处不在的一种现象, 是在不同时空尺度偶然发生且不可预知的事件, 直接影响生态系统的结构和功能演替^[1]。自20世纪50年代Elton^[2]和MacArthur^[3]提出生态系统稳定性理论后, 景观稳定性的评价研究成为生态学领域备受关注的热点问题之一。等级斑块动态范式理论认为生态系统是由不同性质、不同大小且相互兼容的斑块镶嵌体所组成的等级系统, 从状态上强调生态系统是一个不断变化的动态过程, 并由斑块镶嵌体所组成的等级系统的总体变化体现出来, 景观结构特征则在很大程度上描述了景观发生变化的可能性, 斑块在各尺度上的动态变化反映了景观系统中总体的动态变化特征, 进而影响景观稳定性水平^[4]。而以往景观稳定性研究多是借用生态系统稳定性的一些概念, 如抗性、持久性、振幅、韧性、弹性、脆弱性、惰性等, 但由于在景观平衡中各种干扰因素下景观恢复能力在时空尺度上具有较强的复杂性, 景观稳定性程度往往取决于各类型景观斑块的稳定性程度及其组成结构能否保持景观生态过程的连续性与功能稳定性^[5-8]。

近年来, 国内外学者从景观结构和格局角度, 在景观格局变化对景观稳定性的影响开展了大量的研究。Hermosillaa等^[9]对加拿大森林景观稳定性研究发现, 景观破碎化增加的区域表现出森林覆盖和平均森林斑块面积减少、森林斑块数量和森林与非森林的连接数量增加的状态, 进而导致景观稳定性降低；Sharma等^[10]在吉美多吉国家公园森林景观稳定性的评价中也证实了森林斑块数量增多、边缘对比度高且整体呈支离破碎的森林区域, 资源压力加大, 生态系统和山区生物多样性受到影响, 景观稳定性水平降低；同时, 大量在城市、流域、森林、湿地、荒漠、绿洲等景观稳定性评价研究也证实景观破碎化程度的加剧, 会导致景观稳定性水平的降低^[11-17]。然而, 不同研究对景观稳定性的研究尺度、计算方法或表征方式多样, 景观稳定性评价模型所采用的指标、因子及各类指数都有所不同^[18-20]。李余等^[21]以聚集指数、最大斑块指数、斑块密度和形状指数为基础, 采用自组织地图(SOM)与k均值聚类相结合的方法识别森林空间格局, 并从空间格局变化结果得出空间格局相对完整的林区具有最大稳定性；张岳等^[22]利用斑块密度、最大斑块指数、景观形状指数和邻接度指数识别森林空间格局类型, 采用频繁序列模式挖掘方法对其时间序列进行检测, 汾河流域森林景观在相对完整的景观单元具有较高的稳定性, 其演变过程多表现为聚集且空间位置相对接近的状态。而在政策与建设开发影响较大的城市区域, 城市土地利用的变化导致斑块数量与小斑块数量增多, 加剧城市景观破碎化程度, 是导致景观稳定性变化重要原因^{[11-12, 23-24]}。并且大规模支离破碎的城市景观相较于连通的相对稳定的景观系统而言, 其抵御能力与恢复能力更弱^[25-26]。高强度干扰下形成的景观破碎化很大程度会导致生物多样性丧失, 降低生态系统服务频率和能力, 威胁生态系统的演化进程^[27-29]。总体上, 从景观破碎化与景观稳定性的关系出发, 筛选具有代表性的景观格局指数, 构建景观稳定性评价模型, 近年来在基于景观尺度进行稳定性评价方面得到普遍共识和广泛应用。

黔中喀斯特山地城市贵阳市地处我国西南喀斯特片区腹地, 地貌类型以山地、丘陵为主, 城市景观表现为“城中有山、山中有城”的“山—城”镶嵌特征。喀斯特山地景观破碎化和复杂性程度高是影响喀斯特山地城市景观生态系统稳定性水平的典型因素^[30]。尤其是随着近年来城市化进程不断加快, 黔中喀斯特山地城市景观格局发生剧烈变化, 景观的破碎化程度进一步加剧, 对脆弱的喀斯特生态环境造成巨大影响, 进而对城市景观稳定性造成极大威胁^[31-32]。如何在快速城市化导致喀斯特山地城市景观破碎化加剧的过程中, 守住生态和发展两条底线, 保障喀斯特山地城市景观稳定性, 已成为现阶段喀斯特山地城市亟待解决的重要科学问题。因此, 非常有必要从景观生态学的角度, 系统评价城市化发展和生态保护耦合背景下喀斯特山地城市景观稳定性及其时空异质性特征, 以科学指导喀斯特山地城市景观空间格局的优化与生态安全格局的构建, 为喀斯特山地城市景观生态安全和社会经济可持续发展提供科学依据。

基于此, 本研究以典型的黔中喀斯特山地城市贵阳市为对象, 以生态文明城市建设和快速城市化发展为节点, 依据景观生态学、GIS空间信息技术和地统计学等理论和方法, 采用2008年、2013年和2017年Landsat遥感图像, 系统分析贵阳市快速城市化和生态文明示范城市建设双重影响下景观稳定性时空变化特征与演变规律, 以期为喀斯特山地城市景观生态修复、国土空间规划与生态文明城市建设提供参考。

1 材料与方法 1.1 研究区概况

贵阳市(106°07′—107°17′E, 26°11′—27°22′N)是贵州省省会城市, 位于云贵高原黔中山原丘陵中部, 平均海拔约1100m。气候属于亚热带湿润温和型气候, 年平均气温15.3℃, 夏季平均气温24℃, 年平均相对湿度77%。贵阳市下辖6区3县和1个县级市, 国土总面积约8034km², 2020年森林覆盖率达55%。2000年至今, 贵阳市城镇化率不断上升, 在岩溶地貌城市强烈发育的过程中, 其城镇化过程中的环境问题发生显著的变化^[33]。

1.2 数据源与数据预处理

研究采用贵阳市2008年的Landsat-5 TM、2013年和2017年的Landsat-8 OLI遥感图像为主要数据源, 图像来源于地理空间数据云网站(http://www.gscloud.cn/), 行列号为127/041和127/042, 不同时期影像获取季相一致, 且含云量小于5%, 图像质量良好, 空间分辨率为30m×30m, 所有图像坐标系统统一为高斯-克吕格投影CGCS2000_3_Degree_GK_Zone_36。

基于ENVI5.3软件平台, 对各时期遥感影像进行辐射定标、大气校正、图像裁剪等预处理, 依据《土地利用现状分类》(GB/T21010—2017)标准, 结合研究区实际, 运用支持向量机监督分类方法, 将景观类型划分为耕地、林地、草地、水域、建设用地和未利用地等6类, 得到研究区2008—2017年景观类型空间分布图(图 1), 运用随机抽样法在各时期景观类型分布图中分别选取150个验证点(共450个验证点), 通过野外调查和Google地球高分辨率图像, 采用混淆矩阵和Kappa系数精度验证的方法对分类结果精度进行评价, 精度验证结果表明3期图像分类结果总体精度均达到90%以上, Kappa系数大于0.86, 分类结果精度满足研究需要。

1.3 研究方法 1.3.1 景观格局指数计算

景观指数高度浓缩景观格局信息, 能够有效地反应景观结构组成与空间配置^[34], 可以通过相应的景观格局指数定量反映不同水平上的景观格局演变趋势^[35-36]。综合考虑研究区实际情况, 从景观结构、形状与异质性三个方面, 选取反映土地利用景观格局特征斑块密度(PD)、总边缘对比度指数(TECI)、蔓延度指数(CONTAG)、景观形状指数(LSI)、Shannon多样性指数(SHDI)和Shannon均匀度指数(SHEI)等6个指标。利用ArcGIS 10.8软件的Fishnet工具, 为保障景观样本大小、数量足够反映研究区景观信息的前提下, 依据景观样本面积应为斑块平均面积2—5倍的原则^[37-38], 按4km×4km大小划分评价单元, 剔除边缘碎小格网, 共得到587个规则格网单元。进而, 基于Fragstats 4.2景观格局分析软件, 分别计算2008年、2013年和2017年景观格局指数, 各指数的计算公式及生态学意义详见文献^[39-42]。

1.3.2 景观稳定性评价模型构建

在利用景观指数建立评价模型时, 单个指数并不能较全面的反映系统景观稳定性, 各指数应当具有相互独立性且单个指数能够较好的描述景观格局特征^{[40, 43-44]}。本研究采用Spearman秩相关分析方法, 选取相关性较小的景观指数构建景观稳定性评价模型, 结果表明各指数间的相关性均通过双侧显著性检验(P=0.000＜0.05), 并且CONTAG和TECI、PD和TECL的相关系数在2008年、2013年和2017年均为最小(图 2)。因此, 选取CONTAG、TECI和PD指数, 并依据等级斑块动态理论, 景观镶嵌体的蔓延度越高、斑块密度和总边缘对比度越小, 景观系统越不容易因干扰或扰动而发生改变, 景观稳定性水平越高^[11], 构建景观稳定性评价模型公式如下：

式中, S代表景观稳定性(Landscape Stability), C代表蔓延度指数(CONTAG), P代表斑块密度(PD), T代表总边缘对比度指数(TECI)。S值越高越趋于稳定状态, 反之则趋于不稳定状态。

1.3.3 景观稳定性等级划分

借鉴罗攀^[45]对景观稳定性等级的划分方法, 将2008年、2013年和2017年研究区587个评价单元的景观稳定性指数, 运用样条插值法插值分析并归一化处理, 得到贵阳市2008—2017年景观稳定性空间分布图(图 3), 采用自然断点法, 将景观稳定性分为不稳定、较不稳定、较稳定、稳定、极稳定5个等级(图 4)。

1.3.4 数据统计分析

采用SPSS 25.0多元统计分析软件, 不同年份各景观指数进行Shapiro-Wilke test正态性检验, 结果发现各年份的PD、TECI和CONTAG指数呈非正态分布(P=0.000＜0.05), 而Kruskal-Wallis H非参数检验方法的分析结果表明拒绝零假设, 三组数据之间至少存在一组及以上有显著差异(P=0.000＜0.05), 因此, 对2008、2013和2017年3组独立样本进行Kruskal-Wallis H多重比较中的成对比较, 分析不同年份间各指数变化特征。同时, 本研究基于ArcGIS 10.8软件平台, 采用探索性空间数据分析的全局空间自相关分析(Moran′s I指数)、局部自相关分析(LISA指标)和空间关联指数Getis-Ord G_i^*分析^[46-47], 揭示景观稳定性的空间差异和空间关联程度, 识别不同景观稳定性等级的冷热点时空聚集分布特征。

2 结果与分析 2.1 景观格局指数变化特征

结合贵阳市2008—2017年景观类型空间分布图(图 1), 统计分析不同时期各景观格局指数的变化(图 5)。由图 1可以知：10年间, 贵阳市景观结构发生明显变化, 表现为耕地面积的大量减少(1543.85km²)和林地、建设用地的持续增加(分别增加787.29km²和619.20km²), 优势景观类型由2008年的耕地和林地转变为2017年的林地主导型景观格局。

从各景观指数的变化来看(图 5), 2008—2017年, PD指数值呈小幅减小(P＞0.05)和显著增大的趋势(P＜0.05), 且2017年PD指数离散化程度增加, 数据波动大, 景观破碎化程度加剧；TECL指数持续增大, 且2017年显著高于2008年和2013年(P=0.00＜0.05), 说明研究区人工化的高对比度斑块边界增加且更加分散, 但其数值分布呈集中化趋势, 变异程度有所减缓；CONTAG指数呈先增后降的变化特征, 且2013—2017年变异系数有明显增大的趋势(P=0.00＜0.05), 说明这一时期景观破碎化程度有减弱的趋势。

2.2 景观稳定性总体变化特征

统计分析2008—2017年景观稳定性时空变化特征(表 1)。结果表明：2008—2017年研究区景观稳定性等级主要以不稳定和较不稳定为主, 二者占研究区总面积的86%—97%, 不同时期景观稳定性水平变化趋势不同, 呈现出“先降低后升高”的总体变化趋势, 其中：2008—2013年以景观稳定性减弱为主要特征, 表现为较不稳定及以上等级面积均不同程度减小, 共减小801.46km²；2013—2017年景观稳定性水平得到极大提高, 较稳定及以上等级面积大幅增加, 并基本恢复至研究初期2008年的景观稳定性水平, 其中较稳定和极稳定等级增幅最大, 分别增加111.57%和123.98%。总体上, 2017年较稳定及以上等级的面积比例由2008年的9.65%增加至2017年的13.89%, 增幅达43.92%, 景观稳定性水平呈现逐渐好转的趋势。

从各等级景观稳定性空间分布图可知(图 4)：研究区北部和中东地区, 该时期稳定和极稳定等级面积最大, 约为10年间稳定和极稳定等级面积变化的47.03%—56.11%。2013年较稳定以上等级同比增长-73.48%, 零散分布于研究区南部的乌当区、东南部的花溪区和中西部的观山湖区；2017年较稳定及以上等级的分布范围剧增, 同比2013年增幅达442.59%, 且主要集中分布于贵阳市东部和南部各区县。

进一步统计分析不同时期各景观稳定性等级的面积转移特征, 可以看出(图 6)：2008—2013年以不稳定等级的大量转入和其他稳定性等级的大量转出为主要特征, 其中, 极稳定等级全部转为不稳定等级, 较不稳定、较稳定和稳定等级转入不稳定等级的面积达1858.63km²；2013—2017年以不稳定和较不稳定景观的转出, 其他各景观稳定性等级不同程度的转入为主, 其中：不稳定和较不稳定等级分别转出面积约1354.12km²、1603.14km², 约占这一时期总转出面积的98%；此外, 除较稳定等级向不稳定和较不稳定等级的少量转移外, 其它各景观稳定性等级均表现出低稳定性等级向高稳定性等级转化的趋势。

2.3 景观稳定性时空异质性特征 2.3.1 空间自相关分析

采用ArcGIS 10.8空间统计工具Moran′s I对贵阳市2008—2017年景观稳定性空间分布进行空间自相关分析, 结果显示：2008—2017年三期的全局自相关Moran′s I指数分别为0.4005、0.2008、0.2499, Z值分别为13.5350、7.2531、8.3219, 且通过5%的显著性检验, 表明贵阳市近10年景观稳定性具有一定的空间正相关性, 并存在明显的集聚性。但10年间Moran′s I指数呈先降后升、总体下降的变化, 2008—2013年景观稳定性的空间聚集性下降, 空间自相关程度减小, 空间趋同性减弱, 2013—2017年的变化则与之相反, 景观稳定性聚集状态逐渐增强。

为进一步探测贵阳市不同等级景观稳定性聚集区的位置和空间演化规律, 采用ArcGIS10.8软件的Anselin Local Moran′s I工具得到2008年、2013年和2017年的局部自相关LISA聚集图(图 7)。由图 7可知：2008—2017年贵阳市景观稳定性空间聚集状态以无明显聚集为主, 均占研究区总面积的82%以上, 2008年景观稳定性水平的“高-高”聚集区主要分布在息烽县北部、乌当区中东部和开阳县东北部区域, 面积约392.34km², 占研究区总面积的约4.88%；2013年“高-高”聚集区面积大幅减小, 面积约259.17km², 较2008年减少133.17km², 且分布区域发生明显变化, 息烽县北部和开阳县东北部的“高-高”聚集转为“低-低”聚集, 乌当区“高-高”聚集区分布范围缩小, 而花溪区中东部和西南部出现“高-高”聚集；2017年的“高-高”聚集区在2013年的基础上有所扩大, 面积增至302.19km², 并在开阳县东南部新增零散分布的“高-高”聚集区, 进一步说明这一时期贵阳市景观稳定性水平有所回升。“低-低”聚集区的分布在2008—2017年呈现由零散分布向相对集中分布的变化趋势, 2008年“低-低”聚集区面积约516.88km², 零星散布于研究区西部和南部各区县内；2013年“低-低”聚集区相对聚集分布于清镇市西部、开阳县东北部、息烽县北部和中心城区大部, 面积增至1009.39km²；至2017年“低-低”聚集区面积约1087.72km², 分布范围再次发生变化, 息烽县北部分布扩大, 清镇市西部“低-低”聚集区消失, 中心城区的“低-低”聚集状态呈零散分布, 但在清镇市中部、东南部和花溪区西部出现相对集中的“低-低”聚集。

2.3.2 稳定性热点分布与变化分析

利用ArcGIS10.8软件的关联指数Getis-Ord G_i^*分析2008年、2013年和2017年贵阳市景观稳定性的热点空间分布, 通过90%、95%和99%的景观稳定性热点分布置信水平可知(图 8)：2008年景观稳定性的热点区面积为520.78km², 占研究区总面积的6.5%, 是贵阳市10年间景观稳定性热点区聚集面积分布最大的时期, 主要分布于息烽县北部和乌当区中部, 部分散布于开阳县东北方向；2013年热点区聚集面积为227.58km², 相对2008年减少293.20km², 主要集中在乌当区和散布于花溪区和开阳县的边缘；2017年热点区聚集面积明显增大, 占研究区总面积的6.23%, 约是2013年热点区聚集面积的1.2倍；主要分布于开阳县、乌当区, 零星散布于清镇市和花溪区边缘区域。近10年的贵阳市景观稳定性无冷点聚集, 不显著的区域占研究区总面积的94%以上。

3 结论与讨论 3.1 讨论

景观是由许多干扰形成的不断演替变化的斑块集合, 除自然形成的干扰外, 人类活动对景观演替的方向具有很大影响^[43]。在城市化建设中, 不同土地利用的决策引导景观类型的发展方向, 使得景观格局的演替途径在不同时期发生变化, 进而导致景观稳定性的时空异质性^[30]。贵阳市在生态文明示范城市建设和快速城市化发展双重影响下, 景观类型发生明显变化, 表现出城市化和生态化趋势, 这与王志杰等^[47]、韩会庆等^[48]、任红玉等^[32]人的研究结果一致, 反映了贵阳市快速城市化和生态保护与修复双重影响的结果。相关研究表明, 景观斑块相对集中, 斑块完整性较好, 连通性较强, 有助于提高总体景观稳定性^[49-50]。而景观结构的变化, 加之喀斯特山地城市本底景观的破碎性, 10年间贵阳市景观破碎化程度加剧, 景观斑块密度增大, 斑块蔓延度和聚集度降低, 人工景观边界愈发明显, 导致研究区景观稳定性整体处于较低水平。

等级斑块动态理论认为, 景观破碎化程度越高, 景观稳定性越差, 而景观斑块相对完整的区域, 景观稳定性越高^[11-12]。本研究中, 景观稳定性热点区主要集中在城镇建设用地或林地斑块相对集中连片、完整的区域, 这与京津冀地区景观稳定性研究的结果基本一致^[11]。另外, 受不同时期城市景观时空格局及其驱动因素的差异影响, 不同时期景观稳定性冷热点分布也随着城市化的强度及生态保护工程的实施在空间上发生变化。2008—2013年间, 贵阳市处于快速城市化发展阶段和生态保护修复相关措施的实施阶段, 导致景观类型之间的转移频繁, 尤其是城镇区域建设用地面积快速增加, 但受喀斯特山地城市开发利用空间不足的限制, 建设用地以地势相对平坦的山间坝地或谷地以点缀式或跳跃式扩展^[51], 空间分布相对分散, 加剧了景观破碎化程度, 并且在退耕还林等生态保护与修复区域(如清镇市西部和开阳县东北部区域), 在相关工程措施的实施期, 耕地和林地两种景观类型交织, 分散分布的耕地和尚未连通连片的林地导致这些区域景观破碎化加剧, 使该区域景观稳定性水平明显降低；2013年以后, 随着城市化进程不断深入, 相关生态保护工程相继完成, 城镇化水平较高的区域建设用地相对集中连片, 林地斑块趋于完整和连通, 景观破碎化趋势有所减缓, 景观稳定性水平有所提高, 但同一时期, 城市化扩展强度较大的息烽县北部、清镇市中部和东南部、以及花溪区西南部, 城市建设用地斑块增加且镶嵌于自然景观基质中, 导致景观破碎化程度加剧, 稳定性水平降低(图 1, 图 3)。因此, 基于喀斯特山地城市自然景观背景的破碎性, 在快速城市化和生态文明建设的耦合时期, 为提高喀斯特山地城市景观稳定性水平, 增强景观抗干扰能力, 在今后的城市化建设、生态保护与修复过程中, 应着重加强景观斑块的相对完整性, 集约化建设用地和耕地, 修补破碎的斑块, 强化生态用地斑块连通性, 提高景观蔓延度或聚集度, 减缓景观破碎化程度, 提升景观稳定性水平, 为构建喀斯特山地城市景观生态安全和生态网络格局奠定稳定的景观基础, 促进喀斯特山地地区生态与发展的协调与平衡。

景观稳定性是景观生态学领域复杂又重要的命题, 科学的评价模型构建是景观稳定性研究的关键内容^{[11, 52]}。虽然景观格局指数在景观生态学研究的应用广泛, 但准确解释景观格局指数的生态学涵义, 建立景观格局指数与生态过程的联系, 依然是景观生态学研究亟待解决的难题^{[10, 45]}。由于景观的时空异质性特征, 不同区域反映景观稳定性的指标会有所差异。景观生态学理论认为景观格局的破碎化和复杂性是衡量景观稳定性的评价标准之一, 一般可以利用景观破碎度、聚集度等指标来刻画^[53]。本研究依据等级斑块动态理论的相关原理, 采用秩相关分析方法, 筛选出景观蔓延度、总边缘对比度和斑块密度等3个指数构建景观稳定性评价模型, 其中景观蔓延度指数主要反映景观组分的空间配置特征和某种优势景观类型的连续性, 是衡量景观斑块类型空间聚集特征的重要指标；而斑块密度和总边缘对比度指数直接反映景观被分割的破碎化程度和景观空间异质性特征, 一定程度上反映了人为活动对景观的干扰程度^[53]。因此, 综合景观蔓延度、总边缘对比度和斑块密度3个指数构建的景观稳定性评价模型, 对于揭示景观破碎化及其变化与景观稳定性的关系具有一定的代表意义。但时空异质性是景观稳定性的典型特征, 并且景观指数对不同异质性景观的表征意义不同, 本研究针对典型的黔中喀斯特山地城市所筛选出的景观格局指数, 虽能较好地表征喀斯特破碎景观的稳定性, 但其普适性及其科学内涵尚需在不同区域景观稳定性评价的研究中进一步验证。同时, 景观稳定性具有抽象性的特征, 景观稳定与否是其内部斑块之间物质流动、能量流动和信息交换综合作用的结果, 因此, 需结合研究目标和研究对象的实际情况, 辅以多种研究方法, 揭示景观格局特征与生态过程的关系, 基于景观格局角度构建科学合理的景观稳定性评价模型亟需进一步深入研究。

3.2 结论

2008—2017年贵阳市景观格局发生剧烈变化, 主导型优势景观类型由2008年的耕地和林地转变为2017年的林地为主导类型, 且建设用地面积持续增加。在研究时段内, 贵阳市景观稳定性整体水平较低, 以不稳定和较不稳定等级为主, 表现出“先快速降低、后缓慢提高, 总体略有降低”的变化趋势。此外, 贵阳市景观稳定性的空间分布与演变具有明显的空间聚集效应, 全局Moran′s I指数介于0.2008—0.4005之间, 但聚集趋势呈“总体减弱, 小幅提高”的特征。