在城镇化进程中，土地利用的变化可能导致土地利用模式和景观格局的改变，故社会经济发展与生态环境保护的平衡日益受到关注[1]。景观生态风险是指在自然或人为因素影响下景观格局与生态过程相互作用过程中可能产生的不利生态后果，其依托景观生态学的生态过程与空间格局的耦合关联视角，更加注重生态风险的时空异质性和尺度效应，致力于实现多源生态风险的综合表征及其空间可视化[2]。
目前，相关学者针对流域[3]、城市[4]、典型地貌区[5]和生态脆弱区[6]等展开了较为深入的研究，为推进生态共建共治，优化土地利用格局提供了参考。关于景观生态风险防范与治理的研究内容在不断扩展，但关注典型城市的人地互动区因旅游业快速发展带来的景观生态风险响应机制的研究较少。旅游业的发展在一定程度上推动了经济社会的发展，但其本身具有资源依托和能源消耗特征，同时需不断提高基础设施建设水平以满足现代旅游需求，土地利用方式也随之改变[7]。学者针对以旅游驱动的土地利用变化对区域生态系统的影响展开了深入研究，如朱顺顺等[8]采用功能价值法，分析了黄山风景区生态系统服务价值的动态变化特征及其旅游业发展的增值效应；钟莉娜等[9]以武夷山市为例，搭建了生态系统服务视域下旅游导向型城镇化的研究框架；李志远等[10]测度了31个省的旅游生态韧性水平，凸显了旅游生态系统复合性及动态演化特征。综上，已有研究多关注旅游与生态系统服务、生态韧性等的关系，较少涉及与生态风险之间的关系。科学评价景观生态风险是保障地区生态安全的重要前提。生态旅游是生态系统价值转换的重要方式，因此，探究旅游地区开发建设与生态系统安全间的潜在联系，对于该地区土地资源可持续管理与生态旅游可持续发展具有重要意义。
鉴于此，本研究以典型旅游城市H为对象，通过构建景观生态风险评价指标体系，识别该地区景观生态风险时空分布格局，并探究其空间联结效应，借助地理探测器揭示景观生态风险与旅游城市自然社会因素间的潜在联系，为该地区及生态相似区域在旅游业驱动的快速城镇化进程中防范和治理生态风险提供参考。
1 材料与方法
1.1 研究区基本情况
研究区区域总面积9 807 km²，地形以山地为主，海拔在61～1 806 m，高差较大，主要山体呈东西走向。境内旅游资源丰富多样且品级较高，拥有2处世界遗产、1处世界地质公园、3处国家5A级景区和3处国家级风景名胜区，是典型的旅游城市。研究区旅游总收入占比较高。旅游业的快速发展促进了城镇化进程，该地景观格局受开发影响。
1.2 数据来源
本研究中土地利用数据来源于武汉大学CLCD数据集（https://essd.copernicus.org/articles/13/3907/2021/），空间分辨率30 m。高程与坡度数据来源于地理空间数据云（https://www.gscloud.cn/）ASTER GDEM数据集，空间分辨率30 m。GDP（Gross domestic product）、人口密度、年均温度和年均降水数据分别来源于中国科学院资源环境科学与数据平台（https://www.resdc.cn/）中国GDP空间分布公里网格数据集、中国人口空间分布公里网格数据集和中国气象要素年度空间插值数据集，空间分辨率1 km。归一化植被指数（Normalized vegetation index，NDVI）数据来源于地理资源数据云（www.gis5 g.com）30 m分辨率逐年最大NDVI数据集，空间分辨率30 m。行政区划数据来源于国家基础地理信息中心（https://www.ngcc.cn/）。
1.3 研究方法
1.3.1 土地利用变化特征分析
土地利用转移矩阵可定量描述一定时间跨度内研究区内各地类的转移情况，体现特定时间内土地利用变化的方向与流量，帮助了解研究区内地类转换趋势，其通过空间相交实现[11]，计算如式（1） 。
式中， 为研究初期i地类转为j地类的面积，转移矩阵中行为初期时第m种地类，列为末期第n种地类。
1.3.2 景观生态风险量化
景观生态风险不仅受到外部环境的影响，还取决于景观内部的脆弱性，为表征研究区景观生态风险，基于景观格局指数构建景观生态风险评价模型。景观格局指数是基于土地利用格局计算得出的用于表征研究区景观组成与空间配置特征的定量指标，通过景观格局指数构建景观生态风险评价模型，从空间上定量反映人类活动改造景观格局与生态风险的联系[12]。同时，为空间量化及合理分析研究区景观生态风险，以平均斑块面积2～5倍的原则综合考虑研究区面积大小及计算量，选取5 km×5 km格网对研究区进行划分，共得到475个评价单元。借鉴已有文献，从研究区景观结构出发，选取景观干扰度、景观损失度和景观脆弱度指数构建景观生态风险指数，并将各评价单元的景观生态风险指数作为该评价单元中心质点的景观生态风险值，在此基础上分析研究区景观生态风险的时空演变趋势。计算如式（2 ）～（7 ）。
式中， 用于表征景观类型受外界干扰程度的大小，其值越大越容易引发生态风险，其由景观破碎度 、景观分离度 和景观优势度 构成， 、 和 为对应景观格局指数的权重，需保证其和为1，根据已有文献分别赋值0.5、0.3和0.2； 反映景观受分割造成的碎片化程度， 为第i类景观的斑块个数， 为第i类景观的总面积； 表征景观在空间上的离散程度， 为景观总面积； 反应景观类型在研究区的重要程度， 为有i类景观的单元个数与单元总个数的比值， 为i类景观的斑块个数与斑块总数的比值， 为i类景观的面积与单元总面积的比值； 为景观损失度，表示外部因素致使景观受损的程度，其由景观脆弱度 和景观干扰度 构成， 表示景观生态系统内部的易损性，通过专家打分法得出，耕地、林地、草地、水域、建设用地和未利用地赋值归一化后分别为0.19、0.10、0.14、0.24、0.05和0.29； 是第k个评价单元内的景观生态风险指数， 是第k个评价单元内第i类景观的面积， 是第k个评价单元内的地类总面积。
1.3.3 空间自相关分析
为描述研究区景观生态风险指数在特定位置的属性值是否与相邻空间的数值存在显著的空间联结效应[13]，借助Geoda软件的空间自相关分析功能，通过计算三期莫兰指数（Moran’ I）表征景观生态风险值在整个研究区内的全局空间相关性， 指数的值域为[-1,1]，其值为负表示空间变量呈现空间负相关，其值为正表示空间变量呈空间正相关，其值为0则表示空间变量不存在空间自相关。通过绘制局部空间自相关LISA空间聚类图反映研究区景观生态风险值在局部的聚集情况，聚集类型分为高—高、高—低、低—高、低—低和不显著5种类型。全局自相关系数及局部自相关系数的计算如式（8 ）～（9 ）。
全局
局部
式中， 与 是i和j像元上的对应空间属性变量的观测值， 是区域空间变量的均值， 是i和j像元的空间邻接权重矩阵，由矩阵 构成，其中 为空间像元个数。
1.3.4 景观生态风险空间分异影响因素分析
地理探测器由因子探测器、交互作用探测器、风险探测器和生态探测器组成[14]，本文主要利用其中的因子探测和交互探测两个模块。该模型是一种可定量探测地理现象空间分异性并分析其驱动机理的统计学方法，结果采用q统计定量衡量影响生态系统服务的因子贡献值大小及其交互作用强弱，q值取值区间为[0，1]，q值越大表明解释力越强，其具体计算如式（10） 。
式中，h为变量或因子的分层 ， 和 分别为全区单元数和第h层的单元数， 和 分别为第h层和全区变量的方差，SSW为该层内的方差之和，SST为全区的总方差。本研究从自然和社会两个层面共选取数字高程模型（Digital elevation model，DEM）、GDP、NDVI、年均降水、坡度、人口密度和年均温度共7个驱动因子，对应以x 1～x 7表示，因地理探测器适用于类别数据，因此采用SPSS软件的K值聚类法将7个驱动因子分为6个层级，并分别赋予1～6的值，提取至对应评价单元中。
2 结果与分析
2.1 土地利用变化特征
研究区2000—2020年土地利用格局如图1所示。2000—2020年，研究区主要土地利用类型为林地，高森林覆盖率决定了其优异的自然环境禀赋，为游客提供了良好的旅游体验；研究区耕地主要位于各区的建设用地聚集区周围，主要分布于A、B、D、E和G地区；水域主要为H区北部的T湖，草地与未利用地面积占比较小。从时间维度来看，20年间研究区建设用地增加明显，由原来的各区已有建成区中心向外拓展，其中A、D、E和F地区建成区增加范围较大。
2000—2020年研究区土地利用方式的转变情况见表1。部分耕地主要转出地类为林地与建设用地，在前一个10年耕地着重用于补充生态用地，在后一个10年耕地主要用于开发建设；林地面积呈先增后降、总体呈增加趋势，由8 796.46 km²增加至8 986.83 km²，主要与耕地发生互换；水域面积有所增加。
表1 研究区2000—2020年土地利用转移矩阵单位：km2 |
| 年份 | 用地类型 | 地类面积 | ||||||
|---|---|---|---|---|---|---|---|---|
| 耕地 | 林地 | 草地 | 水域 | 建设用地 | 未利用地 | 总计 | ||
| 2000—2010年 | 耕地 | 824.73 | 240.48 | 0.44 | 8.95 | 33.72 | 0 | 1 108.32 |
| 林地 | 48.90 | 8 743.43 | 0.50 | 0.71 | 2.92 | 0 | 8 796.46 | |
| 草地 | 0.35 | 1.06 | 1.74 | 0.01 | 0.05 | 0 | 3.21 | |
| 水域 | 0.98 | 1.52 | 0 | 86.51 | 0.23 | 0 | 89.24 | |
| 建设用地 | 1.48 | 0.35 | 0 | 3.57 | 58.60 | 0 | 64.00 | |
| 未利用地 | 0 | 0 | 0 | 0 | 0 | 0.67 | 0.67 | |
| 总和 | 876.44 | 8 986.83 | 2.68 | 99.75 | 95.52 | 0.67 | ||
| 2010—2020年 | 耕地 | 766.00 | 57.08 | 0.26 | 2.91 | 50.20 | 0 | 876.45 |
| 林地 | 181.10 | 8 798.01 | 0.46 | 0.79 | 6.46 | 0 | 8 986.83 | |
| 草地 | 0.25 | 0.48 | 1.60 | 0.01 | 0.34 | 0 | 2.68 | |
| 水域 | 2.22 | 0.54 | 0 | 96.53 | 0.46 | 0 | 99.75 | |
| 建设用地 | 1.91 | 0.28 | 0 | 0.72 | 92.61 | 0 | 95.52 | |
| 未利用地 | 0 | 0 | 0 | 0 | 0 | 0.67 | 0.67 | |
| 总和 | 951.49 | 8 856.38 | 2.32 | 100.95 | 150.08 | 0.67 | ||
2.2 景观生态风险时空演变特征
运用克里金法将研究区评价单元的景观生态风险指数进行空间插值，得出2000—2020年研究区景观生态风险空间分级（图2），分级标准按2000年自然断点法统一划分，分为低风险区（0.024 68≤ERI）、较低风险区（0.024 68<ERI≤0.041 81）、中风险区（0.041 81<ERI≤0.067 92）、较高风险区（0.067 92<ERI≤0.116 06）和高风险区（ERI>0.116 06）。为便于比较分析研究区各时间节点的景观生态风险时空演变特征，通过面积统计法得出各等级风险的面积及占比和空间叠加法得出风险的空间变动情况。
从整体空间分布来看（图2），研究区景观生态风险总体维持在较低水平，但不同等级的风险区空间分异特征明显，其中低风险区与林地分布具有较高程度的重叠，除境内T湖流域、G地城区，B区、E区、F区和D地集中建设区及研究区东南边陲外，大部分地区生态系统较为稳定，景观生态风险等级较低，这些地区以生态保护优先为原则，严格限制建设活动，景观斑块整体性较强，生态安全得到保障。较高以上风险区面积占比较小，主要位于T湖地区及西南边陲地区，此类地区由于地类交错和水域内生性脆弱等致使易诱发生态风险问题。研究区内人类活动较为频繁的耕地和建设用地集中区风险等级为中度，源于建设用地对周边地类的干扰及耕地景观内生性中度脆弱度。从时间维度来看（表2），研究区生态低风险区面积在20年间持续上升，由68.47％上升至74.83％，较低风险区占比则由21.45％降至17.52％，中度及以上风险区面积明显减少，占比由10.09％下降至7.64％。综合来看，20年间研究区景观生态风险呈现逐步改善的趋势，以生态优先的发展模式进一步巩固了林地的景观整体性，减轻了对生态系统的侵扰。
表2 研究区2000—2020年各风险等级面积及占比 |
| 风险等级 | 2000年 | 2010年 | 2020年 | |||
|---|---|---|---|---|---|---|
| 面积/km² | 占比/% | 面积/km² | 占比/% | 面积/km² | 占比/% | |
| 低风险区 | 6 889.32 | 68.47 | 7 408.57 | 73.63 | 7 529.32 | 74.83 |
| 较低风险区 | 2 158.02 | 21.45 | 1 749.91 | 17.39 | 1 763.30 | 17.52 |
| 中风险区 | 908.41 | 9.03 | 823.06 | 8.18 | 684.47 | 6.80 |
| 较高风险区 | 87.76 | 0.87 | 62.28 | 0.62 | 67.70 | 0.67 |
| 高风险区 | 18.71 | 0.19 | 18.51 | 0.18 | 17.48 | 0.17 |
根据景观生态风险的空间变化来看（图3），2000—2020年，研究区景观生态风险等级呈现明显降低趋势，降低区域集中在东南部集中建成区、A区T湖外围、G地建成区周边和D地中部；风险升高区域位于A区南部及C区北部。按时间维度划分，0～10年间景观生态风险降低态势尤为明显，第10～20年间T湖中风险区向外略有扩张，A区南部因建设扰动和草地损失度的提高景观生态风险等级出现上升。总体看来，20年间景观生态风险的降低主要源于部分耕地被替换为景观生态稳定性更高的林地与建设用地，研究区合理的生态建设与治理工作使得其生态系统稳定性不断得到提高。
2.3 景观生态风险空间自相关
通过计算全局空间自相关和局部空间自相关得出研究区景观生态风险的空间依赖程度和空间聚集分布情况。由Moran’I可得出，研究区景观生态风险呈现较弱的空间正相关性，且逐年降低，表明景观生态风险在空间上的依赖程度下降，与景观生态风险的好转情况相对应。由局部空间自相关LISA图（图4）可看出，研究区景观生态风险以“低—低”聚集为主，低值聚集区主要位于C地中部、G地南部和G地北部，这部分区域位于山林腹地，区域景观类型单一，受人为干扰程度较小。随着时间的推移，C地南部低值聚集区数量减少，中部低值聚集区数量逐渐增加，原因可能是山林内人类活动不断增强，对林地产生了一定干扰，中部林地密度提高致使景观生态风险降低。“高—高”聚集区位于研究区东南边陲，且数量较为稳定，“低—高”聚集区零星分布于高值聚类区周边。
2.4 景观生态风险影响因素
采用地理探测器定量识别研究区景观生态风险空间分异的影响因素，并判别对景观生态风险起到关键影响的自然和社会因素。由表3可知，GDP、人口密度等社会因素的p值大于0.1，未通过显著性检验，表明社会因素对研究区景观生态风险未产生明显影响。其余5个自然因素均通过显著性检验，2000年因子驱动力排名q值为年均温度>坡度>NDVI>DEM>年均降水，研究区植被覆盖度高，植物种类丰富，受温度、地形影响较大。2000—2020年因子驱动力排名q值未发生变化，DEM的驱动力随时间推移逐渐提高，研究区为山区，海拔高，陡坡多，地势平坦地区受人类开发活动影响其景观生态风险变动较大；NDVI、年均降水量和坡度的解释力逐年递减，这可能源于研究区植被覆盖率高且较为稳定，人类活动对土地的干扰与投入不断加大导致影响植被覆盖的自然因素的影响逐渐降低。
表3 研究区2000—2020年景观生态风险因子探测结果 |
| 驱动因素 | 2000年 | 2010年 | 2020年 | |||
|---|---|---|---|---|---|---|
| q值 | p值 | q值 | p值 | q值 | p值 | |
| DEM | 0.492 | 0.002 6 | 0.518 | 0 | 0.546 | 0 |
| GDP | 0.210 | 0.562 0 | 0.239 | 0.466 0 | 0.123 | 0.754 7 |
| NDVI | 0.525 | 0.008 7 | 0.538 | 0.007 0 | 0.403 | 0.026 0 |
| 年均降水量 | 0.435 | 0 | 0.385 | 0.003 3 | 0.349 | 0.006 6 |
| 坡度 | 0.530 | 0 | 0.531 | 0 | 0.434 | 0.002 8 |
| 人口密度 | 0.310 | 0.224 6 | 0.319 | 0.204 8 | 0.182 | 0.472 9 |
| 年均温度 | 0.566 | 0 | 0.583 | 0 | 0.605 | 0 |
研究区景观生态风险空间分异的因子交互探测结果如图5所示，结果显示，各驱动因子交互时均呈现增强效应，同时具备非线性增强与双因子增强关系，表明驱动因子间的交互作用可明显提升对景观生态风险时空分异的解释能力，即说明研究区景观生态风险的空间分异格局是多因子协同作用的结果。其中，自然因子间的交互作用是导致研究区景观生态风险空间分异的关键诱因，2000年驱动力大于0.8的因子交互类型有x 1∩x 3、x 1∩x 4、x 1∩x 5、x 2∩x 7、x 3∩x 4、x 3∩x 6、x 3∩x 7、x 4∩x 7、x 5∩x 6、x 5∩x 7和x 6∩x 7，驱动力大于0.7的因子交互类型共有15种；2000—2020年因子间的交互作用在数量上未下降但在强度上有所减弱，两两因子交互解释强度均有不同程度的下降。综合来看，因子间的交互作用总体呈现减弱趋势，这与单因子探测中多个因子的解释力下滑相对应，自然因子间的交互作用解释力较强，但自然因子与社会因子的交互带来的解释力提升明显高于自然因子间的交互，表明社会因子虽不能独立对研究区景观生态风险产生明显影响，但其与自然因子之间的协同作用是导致景观生态风险分异的关键原因。这表明研究区地形复杂，植被覆盖度高，人类活动的增强与自然环境变化的交互影响促成了现有景观生态风险格局。
3 结论与讨论
旅游业的快速发展可能致使土地空间格局和生态安全格局在不断演化和重构，人类活动对景观格局的改造和区域生态系统造成了明显影响。通过评估景观生态风险，厘清景观生态风险的影响因素，有助于为区域旅游业高质量发展和土地资源合理利用提供理论支撑。基于2000、2010和2020年三期土地利用数据，借助景观生态风险指数、空间自相关分析和地理探测器等，刻画了研究区景观生态风险的时空演变格局，探究了景观生态风险的自然和社会驱动因素，得出如下结论。（1）2000—2020年，研究区林地和水域等面积得以增加。景观生态风险整体呈现较低水平，各等级风险的空间分异特征明显，山林腹地景观整体性高，发生景观生态风险的概率低，水域及人类活动强度较高的区域，因景观内生性脆弱和人类活动扰动等，景观生态风险处于中度以上水平。20年间研究区景观生态风险明显降低，生态系统稳定性不断提高。（2）研究区景观生态风险呈较弱的全局空间正相关性，且在20年间正相关依赖程度不断降低。景观生态风险的局部空间自相关以位于研究区西南部山林腹地的“低—低”聚集类型为主，且呈减少趋势，高值聚集区位于研究区东南边陲地带。（3）2000—2020年，研究区景观生态风险的时空分异格局主要由年均温度、坡度、NDVI、DEM和年均降水量等自然因素决定，且因子两两交互均呈现出双因子交互增强效应。研究区景观生态风险空间分异是多因子协同作用的结果，自然因子与社会因子交互对景观生态风险格局产生明显影响。
研究区依托其丰富的自然风光与历史人文底蕴成为典型的旅游城市，依托旅游业推进城市经济发展的本质是促进生态环境和社会经济的协调发展。旅游业的发展依托良好稳定的生态环境，需划定“三区三线”，协调土地利用结构和布局，权衡经济效益与生态环境间的关系，在建设用地集中的B、E地及各县的旅游片区，做好城市规划工作，维护城市周边地带的景观整体性，建立绿色生态廊道作为城区与主体生态功能区之间的缓冲。在T湖、H风景区等重要生态保护区加强生态风险监测，依据环境容量和资源承载力严格控制开发强度，建立可持续的生态友好的土地管理和景观维护机制，尽可能降低对生态地类的侵扰，有效提升生态系统的稳定性。