全球经济发展引发的生态环境问题已成为国际社会关注的核心议题^[1]，有研究表明，由人类活动引起的景观类型变化是影响陆地生态系统结构与功能的重要途径之一^[2，3]，其通过改变景观格局、干扰生态过程，直接或间接驱动生态风险的形成^[3]. 因而，量化人类活动和自然灾害对生态系统结构功能的复合胁迫，精准识别环境系统潜在危害与空间脆弱性，可以为国土空间管控优化与区域生态韧性建设提供科学决策支撑^[4，5]. 在生态风险评估体系中，景观生态风险（landscape ecological risk，LER）突破传统风险评估对单一要素或局部区域的局限^[6]，通过系统解析自然-社会复合干扰下对景观格局动态与生态过程的耦合机制，实现了生态风险空间规律的量化、区域生态敏感性的识别及综合性评价体系的构建^[7，8].

现阶段关于LER研究主要围绕两方面开展：一是探究不同景观类型与LER响应关系，二是分析LER时空演变规律及驱动机制^[9，10]. 已有研究表明，在黄土高原地区，受生态工程的积极推动，林草地生态系统的LER持续下降，而城市和农田区域LER因人类活动持续保持在高风险状态^[11]. 与此同时，LER与各驱动因子之间的关系直接或间接影响区域景观格局优化方案制定和生态环境建设实践的开展. 在旱地生态系统中，基于地理探测器的研究进一步揭示，地形和气候条件共同影响LER空间分布特征^[12]；此外，国外有学者亦通过建立“自然-景观格局-人类社会”的多维度生态风险评价体系，并通过空间聚类和主成分分析识别其驱动因素^[13]. 因此，深入研究LER时空演变规律，解析其驱动机制，对保障区域生态安全与可持续发展具有重要意义. 当前，用于探究驱动因子研究方法主要包括地理探测器（Geodetector）^[14]、偏最小二乘法结构方程（PLS-SEM）^[15]、普通最小二乘法（OLS）^[16]、随机森林回归（RFR）^[17]、地理加权回归（GWR）^[18]和多尺度地理加权回归模型（MGWR）等^[19]. 有研究通过RFR评估在人为和自然双重因素驱动下LER演变趋势，但缺乏从空间角度对驱动因子具体作用方向探究^[20]；相比之下，地理探测器可以量化多因子间的交互作用，进而揭示驱动因子的协同效应^[21]. GWR和MGWR模型基于局部回归框架，能够有效揭示各驱动因子对LER影响的空间非平稳性，为解决空间异质性问题提供有效途径^[22]，而MGWR允许不同驱动因子采用各自独立的空间平滑尺度，从而解决了GWR单一空间尺度的局限，因此其结果相对更为合理.

固原市位于西安、兰州和银川这三省会城市构成的三角地带中心，是宁夏副中心城市及宁南区域中心城市，其经济位于宁夏5个地级市生产总值末位，但其增速最快，是全省最具发展潜力的城市之一^[23]. 通过多年的生态建设，固原市生态环境有一定程度的改善，但水土流失面积占土地总面积的80%，生态系统脆弱性突出、稳定性差^[24]. 当前针对宁南山区城市区域景观风险研究较少，尤其在驱动因子多尺度调控上还尚未深入探讨. 因此本研究以固原市为例，基于2000年、2005年、2010年、2015年和2020年这5期景观类型数据，评估区域LER演变特征；采用地理探测器和MGWR筛选出多重共线性问题较小的驱动因子，并分析其空间影响力；从自然和社会经济双重视角揭示固原市LER驱动机制，以期为优化固原市景观布局、维护生态健康提供理论支持.

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

固原市位于宁夏最南部（105°19′~106°57′E，35°14′~36°31′N，图 1），地处黄土高原沟壑区，地势南高北低、沟壑纵横. 总面积达1.05万km²，下辖原州区、西吉县、泾源县、隆德县和彭阳县这5个区县. 研究区年平均降水为300~500 mm^[25]，主要集中在6~9月. 固原市常住人口约为114.21万人，地区生产总值463.73亿元，其中原州区人口总数和GDP均是全市最高^[26].

1.2 数据来源

本研究所用数据主要包括固原市景观类型数据、气象数据和社会经济数据等，所有数据均涵盖年份2000年、2005年、2010年、2015年和2020年. 景观类型数据来源于中国科学院资源环境科学与数据中心（https://www.resdc.cn/），归一化植被指数（NDVI）来源于国家生态科学数据中心（https://nesdc.org.cn/）^[27]. 数据高程模型（DEM）来源于地理空间数据云（https://www.gscloud.cn/），坡度数据来源于对DEM数据提取^[28]. 年平均降水数据来源于国家青藏高原科学数据中心（https://data.tpdc.ac.cn/）. 人口密度和国内生产总值（GDP）数据来源于地理遥感生态网（https://gisrs.cn/）. 夜间灯光数据来源于美国航空航天局（NASA）和美国国家海洋和大气管理局（NOAA）（https://www.noaa.gov）^[29]. 所用数据均以30 m× 30 m分辨率重新采样，投影坐标系统一定义为WGS_1984_UTM_Zone_48N^[30].

1.3 景观类型变化分析

单一景观类型动态度是用于定量描述某一特定区域的某一景观类型在一段时间的变化速度和幅度^[31]，计算公式如下：

K=Ub-UaUa×1T×100%

(1)

式中，K为单一景观类型动态度；U_a和U_b分别为初期和末期的某景观类型的面积；T为时段.

综合景观类型动态度^[32]是衡量某一特定区域景观类型综合变化程度，反映区域内各景观类型之间相互转化的综合速度，其值越大，则区域内景观类型变化越强烈. 计算公式如下：

LC=∑i=1nΔLUi-j2∑i=1nLUi×1T×100%

(2)

式中，LC为综合景观类型动态度；LU_i为初期i类景观类型的面积；ΔLU_i-j为研究时段内i类用地类型转为其他的面积；T为时段.

景观类型转移矩阵是定量描述特定时间段内某一区域各景观类型之间相互转化关系^[33]. 计算公式如下：

Sij=S11…S1n⋮⋮Sn1…Snn

(3)

式中，S_ij为研究初期年第i类景观类型转移为末期年第j类景观类型的面积；i和j分别表示初期和末期的景观类型（i=1，2，…，n；j=1，2，…，n）；n为景观类型数量.

1.4 景观生态风险评价模型 1.4.1 景观生态风险评价模型构建

由于LER具有空间异质性，利用ArcGIS 10.7软件将研究小区按照5 km× 5 km^[34]的网格划分为508个评价单元. 根据固原市2000年、2005年、2010年、2015年和2020年固原市景观类型数据，通过Fragstats 4.2计算各景观类型的景观破碎度、景观分离度和景观分维度指数^[35]. 景观生态风险指数具体计算公式如下：

Ri=Fi×Si

(4)

Si=aCi+bNi+cDi

(5)

Ci=ni/Ai

(6)

Ni=A2AiniA

(7)

Di=2lnPi4/lnAi

(8)

式中，R_i为景观损失度. F_i为景观脆弱度，通过借鉴前人研究并结合研究区实际情况将建设用地赋值为1，林地为2，草地为3，耕地为4，未利用地为5，水域为6，并进行归一化处理. S_i表示景观干扰度，其中a、b和c分别表示景观破碎度、景观分离度和景观分维数指数的权重且a+b+c=1，根据前人研究将a、b和c分别设置为0.6、0.3和0.1^[36]. C_i、N_i、D_i、n_i、A_i和P_i分别为地类i的景观破碎度、分离度、分维数、斑块数量、斑块面积和斑块周长. A为生态小区的总面积.

根据景观指数构建景观风险模型（LER），并采用自然断点法将研究区划分为5个风险等级（低风险：LER≤0.018 0，较低风险：0.018 0 < LER≤0.030 1，中风险：0.030 1 < LER≤0.042 2，较高风险：0.042 2 < LER≤0.054 3，高风险：LER > 0.054 3）^[37]. 景观风险模型计算公式如下：

LER=∑i=1nAkiAkRi

(9)

式中，A_ki为第k个风险小区中地类i的面积；R_i为地类i的景观损失度指数.

1.4.2 空间自相关分析

空间自相关性通过量化空间单元之间属性值的相似性或相关性而来判断这些属性值在空间上是否存在聚集或分散. 本研究基于GeoDa软件计算全局莫兰指数（global Moran's I）指数与局部莫兰指数（local Moran's I）并绘制高低聚类图^[38]，计算公式如下：

$ \text { global Moran's } I=\frac{n \sum\limits_{i=1}^n \sum\limits_{j=1}^n w_{i j}\left(x_i-\bar{x}\right)\left(x_j-\bar{x}\right)}{\sum\limits_{i=1}^n \sum\limits_{j=1}^n w_{i j} \sum\limits_{i=1}^n\left(x_i-x\right)^2} $

(10)

$ \text { local Moran's } I=z_i \sum\limits_{j=1}^n w_{i j} z_j \quad(i \neq j) $

(11)

式中，n为风险评价单元数量；x_i和x_j分别为风险评价单元i和j的景观生态风险值，x为景观生态风险评价均值；z_i和z_j分别为风险评价单元i和j的标准化景观生态风险值；w_ij为空间权重矩阵.

1.5 驱动因素分析 1.5.1 地理探测器

区域景观生态风险受多重因素的影响，传统的线性回归对数据质量要求较高，地理探测器可以探测因变量空间异质性，定量分析各因素对景观风险的影响，计算公式如下^[39]：

q=1-1Nσ2∑h=1LNhσh2

(12)

式中，q为驱动因子解释力，值域为0~1，值越大表明该因子解释力越强；L为因子分区数量；N_h为分区h内样本数，N为样本总数；σ²和σ²_h分别为流域LER方差和分区h内的LER方差.

1.5.2 多尺度地理加权回归（MGWR）

本研究选取MGWR模型探究各因素对景观生态风险的空间分异性规律^[40]，初步筛选出6个多重共线性问题较小因子，包括自然因素（归一化植被指数、降水和坡度）和社会经济因素（人口密度、GDP和夜间灯光指数），探究这些因子对LER响应，模型构建如下：

yi=β0ui,vi+∑k=1pβbwkui,vixik+εii=1, 2,…,n

(13)

式中，yi 为ui,vi 第i个样本地理坐标；βbwkui,vi 为第i个样本第k个参数的回归系数，其绝对值大小代表驱动因子对LER的影响力，正负代表影响方向^[41]；εi 为第i个样本空间单元的随机误差项；p 为样本总数；xik 为第k个样本的第i个解释变量.

本文采用MGWR 2.2软件构建驱动因子与LER直接关系，核函数和带宽选择标准采用为最常用的二次函数和赤池信息量准则AICc^[42]（AICc用来评估模型拟合优度，其值越小拟合优度越好）.

2 结果与分析 2.1 景观类型动态变化特征

2000~2020年固原市景观类型以耕地和草地为主（图 2），主要分布在西吉县和原州区，耕地面积呈现下降趋势，草地总体变化幅度不大，分别占全市总面积的43.78%和45.76%. 林地主要分布在泾源县和原州区，其面积于2000~2010年呈现增长趋势，2010~2020年有略微下降但总体稳定. 本研究期内水域和建设用地主要分布在西吉县和原州区，面积占比最小但变化速度较大，其单一景观类型动态度绝对值均大于1. 未利用地总体变化程度最大，2015~2020年单一景观类型动态度达23.28. 研究期内景观类型综合变化速度较快，2015~2020年期间景观类型变化程度最高，2000~2005年、2005~2010年和2010~2015年综合景观类型动态度均分别为0.22、0.15和0.17（表 1）.

2000~2020年景观类型以耕地转出为主，转出面积为474.02 km²，转移为林地、草地、水域、建设用地和未利用地分别占比6.19%、75.24%、16.08%、2.81%和0.44%. 林地转移面积最小，89.60%面积未发生改变；建设用地扩增明显，54.93%来源于耕地，38.83%来源于草地；水域主要来源于耕地转化.（表 2）

4个研究期内景观类型转移主要特征是耕地转移为草地（图 3），2005~2010年和2015~2020年间耕地转移为草地面积较多，分别转移了227.69 km²和317.25 km². 2000~2005年间林地增幅最大，主要来源于耕地和草地，2005~2010年间建设用地变化最明显，增长面积占总转移面积的30.39%，2010~2015年未利用地最不稳定，转移率达3.01%. 2015~2020年林地基本不变，草地增长最快.

2.2 景观生态风险时空特征分析 2.2.1 景观生态风险时空分布特征

2000~2020年研究区LER指数由0.045 3降至0.032 2，生态风险程度有所改善. 高风险和较高风险集中在西吉县和原州区北部，面积由68.76%降至58.38%，中风险由固原市中部向西南方向扩增，低风险和较低风险主要分布在彭阳县北部和泾源县六盘山一带（图 4），面积随年份不断增加. 2020年5个风险区面积为：中风险区域（33.64%） > 较低风险区域（24.74%） > 低风险区域（18.91%） > 较高风险区域（16.96%） > 高风险区域（5.75%）.

2.2.2 景观生态风险时空演变特征

2000~2020年LER呈现高等级向低等级转移趋势（图 5），高风险和较高风险转移面积随年份逐渐降低；中风险转入面积最多，占总转入面积36.75%. 2000~2005年和2010~2015年LER改善面积最大，分别占总面积30.80%和30.54%. 前期LER降低区域主要分布在泾源县和彭阳县，后期则以西吉县为主要风险改善区域（图 6）.

2.2.3 景观生态风险空间自相关分析

本研究期内Moran's I散点图主要集中在第一、三象限，说明研究区的LER指数在空间上呈正相关，主要呈现“高-高”和“低-低”聚类模式，存在显著的空间集聚性和空间分异性. 2000~2020年的Moran's I分别为0.488、0.532、0.571、0.620和0.685（图 7）. 随年份的增加，空间聚集性不断增强. 局部自相关LISA图表明研究区“高-高”聚类模式主要分布在西吉县和原州区，“低-低”聚类模式主要分布在彭阳县和泾源县（图 8）. 2000~2020年“低-低”聚集区占比相差不大，（平均值为28.66%），“高-高”占比由28.06%下降至22.73%. “低-高”和“高-低”聚集区占比较低，均小于5%，表明研究区的高风险和低风险区比较集中.

2.3 景观生态风险驱动因子分析 2.3.1 多元线性结果分析

各驱动因素多元线性回归结果表明，2000~2020年各驱动因素膨胀系数（VIF）均小于7.5，显著性均小于0.05，说明这6个驱动因子之间不存在多重共线性且对LER影响具有显著性（表 3）.

2.3.2 地理探测因子分析

基于地理探测器因子分析，利用q值定量评估驱动LER空间格局各因素的影响力及其随时间变化特征. 结果显示，所有因子均通过0.01显著性水平检验（P < 0.01）. 单因子结果表明（表 4），本研究期内NDVI和降水解释力始终高于其他因子，是LER的首要驱动因子，人口密度和GDP的解释力随时间不断增加，坡度对景观生态风险解释力呈波动趋势，夜间灯光对景观生态风险解释力不大，2020年各因子q值从大到小为：NDVI（0.589） > 降水（0.371） > 人口密度（0.308） > GDP（0.251） > 坡度（0.171） > 夜间灯光（0.106）.

双因子交互作用显示，研究期间随机两个因子的交互作用均大于单一因子作用（图 9），说明固原市LER受多个因子的综合作用，NDVI与各个因子交互均大于其他双因子交互作用，其中NDVI与人口密度和降水交互作用较大，NDVI与GDP和夜间灯光交互作用呈上升趋势，LER受自然因子和社会经济因子共同影响. 变化最明显的是降水与人口密度交互作用，解释力大小随年份先增加后降低，2015年q值达到最高为0.497，2020年q值0.385.

2.3.3 多尺度地理加权回归分析

对比OLS、GWR和MGWR这3种模型模拟效果，MGWR的AICc相较于GWR和OLS降低了76和186，因而MGWR模拟精度更高. 其次，OLS模型中各驱动因素均为全局尺度，共用一个带宽；GWR模型虽可体现空间异质性，但各驱动因素均用一个最优带宽480；MGWR模型各驱动因素均有专属带宽，因而各因素空间作用尺度更精确（表 5）.

OLS不涉及局部回归且整体R²仅为0.153，其解释能力较弱. MGWR的R²达到0.818，比GWR高0.18；比较MGWR和GWR局部拟合效果，MGWR局部区域细节表现效果更好，局部R²超过0.5的达到80.06%. 因此，筛选出来的驱动因子对LER具有较强的解释力（图 10）. 综上所述，MGWR拟合效果更好，相较于OLS和GWR更能反映各驱动因素与LER直接关系.

自然因素方面，NDVI总体上与LER呈负相关，降水和坡度呈正相关（图 11）. NDVI在空间上具有很强的负向作用，且随着年份增长，负向作用愈强烈. 坡度负向驱动作用较广，正向驱动作用主要集中在地势较高的六盘山一带，研究期内其影响面积不断扩大. 降水量以正向驱动力为主，主要集中在西北部，该区域降水不稳定，易发生极端暴雨事件.

社会经济因素对景观生态风险影响空间差异显著.人口密度和GDP对北部原州区和西部西吉县主要呈正向作用，对于人口密度和GDP较低区域呈负向作用，说明在原州区和西吉县等经济水平处于较高水平的区域，增加人口活动会加剧生态风险；而对于经济相对滞后区域，适当提高人口密度有助于降低生态风险. 夜间灯光指数在空间上对景观生态风险起正向作用，其作用强度变化较大.

3 讨论 3.1 景观生态风险格局演变分析

景观生态风险演变聚焦于特定区域在自然与人为因素作用下，LER随时间推移所呈现出的动态变化^[43]. 研究区高风险和较高风险区域主要分布在原州区和西吉县且在空间上有明显聚集性，这主要是由于二者景观类型以耕地为主，受人为影响大，加之过度耕地开垦导致植被覆盖度降低，容易发生水土流失^[44]. 原州区作为固原市中心城市，其建设用地占全市的42.46%，有相关研究发现建设用地会严重挤占自然系统生存空间，使得景观破碎化和分离程度增加，物种栖息地丧失及生态廊道受阻，导致区域稳定性降低^[45].

2020年固原市低风险和较低风险区共占比43.65%，相较于2000年增加了29.20%，且主要集中在彭阳县和泾源县，这些区域植被覆盖度较高，其中泾源县森林覆盖率达50.10%，景观破碎度和分离度较低，受人为干扰小，生态系统较稳定，因而LER较低^[46]. 彭阳县耕地面积占全县41.18%，其“88542”隔坡梯田工程通过地形重塑和景观格局优化，提高植被覆盖率，避免单一耕地导致的生态退化^[47].

空间自相关结果表明，“高-高”聚集区主要集中在原州区和西吉县且空间聚集性不断增强，因此在生态建设时要有限针对重点区域，合理规划景观格局，此外对“不显著”区域建立合理调控机制，防止生态风险反弹.

3.2 景观生态风险驱动因素分析

为了确定各驱动因子影响力，采用了地理探测器模型进行初步归因；利用MGWR模型聚焦空间异质性，分析不同地理单元内驱动因素对景观生态风险的影响程度，发现各驱动因素对景观生态风险的作用方向呈现正向与负向双重属性.

研究结果表明，固原市景观生态风险受自然因子和社会经济因子共同影响，其中NDVI和降水是主要影响因子，这与谢兴龙等^[48]研究结果一致. 空间上NDVI与LER呈负效应（回归系数 < 0.06），则区域植被覆盖度越高，生态风险越低；说明植被覆盖度可以有效地稳定景观格局进而降低区域生态风险^[49]. 降水在空间上主要与LER呈正效应，尤其在南部六盘山湿润区（年降水量约600~800 mm），该区域山地多于平原且极端暴雨事件频发，而在北方干旱区域（年降水量不超过400 mm），降水呈现一定负效应，因此适当降水有助于植被恢复，从而提高生态系统稳定性. 坡度整体上对LER呈正效应，尤其是西北坡度陡峭区域，在暴雨作用下容易发生沟道侵蚀，侵蚀沟的增生会分割景观单元，降低生态功能，增加自然灾害频率^[50].

社会经济因素对LER影响通常与人类活动密切相关，地理探测器单因子结果表明社会经济因子对LER的影响较弱，其q值均小于0.4，但与自然因子交互作用后，其影响力显著增强. 人口密度与NDVI交互作用最强，q值达0.728. 有研究表明，高强度人类活动对生态系统持续压力会加剧城市风险区的扩增^[51]. 局部空间上看，社会经济因子对LER主要起正向驱动效应且主要集中在原州区，而对于人类活动较弱的南部区域，社会经济因子对LER起抑制作用. 这是由于社会因子与城市建设密不可分，相关研究表明，社会经济因素对植被净初级生产力贡献率达41%，因而，人类可以通过生态工程提高区域植被恢复进而降低区域生态风险.

本研究基于固原市景观类型和LER，参考前人相关研究，针对区域可持续发展提出以下建议：从区域景观类型布局来看，耕地脆弱性和敏感性最高，最易发生侵蚀，针对高风险区域，应该合理规划耕地结构和调整区域景观布局，在满足农业生产同时最大限度降低耕地损耗. 中风险区域应强化生态保护与修复力度，通过景观破碎化治理，推进生态廊道与基础设施网络构建，同步实现生态资源的可持续利用. 低风险区域应以生态监测为主，维持现有的良好生态状况同时避免区域过度开发，防止潜在风险出现. 协调生态建设与社会经济可持续发展，对研究区进行功能划区，强化政府部门在风险管理中主导监督作用，建立健全生态补偿机制.

4 结论

（1）2000~2020年间景观类型以耕地为主，变化特征主要是耕地减少林地增加，水域、建设用地和未利用地面积占比较小但其变化速度较快. 景观类型之间转移剧烈，耕地转出面积最多，耕地转移为林地、草地、水域、建设用地和未利用地分别占比6.19%、75.18%、16.08%、2.81%和0.44%.

（2）2000~2020年间高风险和较高风险集中在西吉县和原州区北部，中风险由中部向西南方向扩增，低风险和较低风险集中在彭阳县北部和泾源县六盘山一带. 研究区的LER随时间而有明显改善，空间上表现明显聚集性，均以“高-高”和“低-低”聚集性为主.

（3）地理探测器结果表明，影响景观生态风险主要因子是NDVI、降水和人口密度，双因子的交互作用均大于单因子，NDVI和人口密度的交互作用最强. MGWR模型拟合效果优于OLS和GWR模型，NDVI对LER起抑制作用，极端暴雨和不利地形会导致LER增加，社会经济因素对LER主要起正向驱动作用.