土地利用与土地覆盖变化深刻影响着全球陆地生态系统结构与功能^［1-2］，尤其在干旱半干旱区，人类活动与自然胁迫的叠加效应加剧了景观破碎化与生态风险^［3-4］。因此，景观生态风险评价成为衡量土地利用类型、空间分布对生态环境的影响程度、潜在风险的重要手段。

景观生态风险（ERI）是指景观要素镶嵌的异质性、格局演变的规律性，以及在外部风险源干扰下响应生态环境的综合生态过程^［5］。随着自然因素和人类活动对景观格局影响不断加剧，生态风险评价成为近年来景观格局-生态过程机制研究的热点。目前，生态风险评价方法主要围绕分析风险源、风险汇和景观格局展开^［6-8］，学者们使用不同方法衡量土地利用的生态风险。建立生态风险指数，根据土地利用结构和生态风险强度参数评价生态风险^［9］，该方法得到广泛认可与应用。根据土地利用变化（LUCC）生态风险综合指数，Fan等^［10］认为黑河中游地区高风险区主要分布在北部，均值为6.60。此外，景观结构指数已被应用于反映生态系统结构衍生的土地利用生态风险水平^［11-12］，景观生态风险由景观干扰度指数（E_i ）、景观损失度指数（R_i ）等构成，包括景观破碎度指数（C_i ）、景观分离度指数（N_i ）和景观优势度指数（D_i ）^［13-14］。现有研究多聚焦城市或流域尺度^［14-16］，针对生态脆弱沙地系统的风险动态及驱动机制研究较少。

毛乌素沙地是黄河“几字弯”风沙防治攻坚战的主战场，地处干旱半干旱区，流沙分布广泛，属于典型生态脆弱区，1953—1976年由于人类过度放牧、乱砍滥伐，毛乌素沙地发生沙漠扩张^［17］，进一步加剧了毛乌素沙地荒漠化程度；1978年，国家启动“三北防护林体系建设工程”，毛乌素沙地作为重点治理区域^［18］；20世纪80年代，治理模式进一步创新；2000年以来毛乌素沙地治理被纳入京津风沙源治理、退耕还林还草等国家级生态工程^［19］；目前仍存在沙地地下水位严重下降、人工固沙林衰退风险较大、沙区可持续发展压力大等问题。因此本文通过研究毛乌素沙地景观生态风险时空演变及驱动力，通过构建高分辨率（3 km×3 km）生态风险网格体系，并融合地理探测器与景观生态风险开展交互因子识别，以沙地风险评价-过程解析-政策验证为研究框架，为风沙防治与生态安全屏障建设提供理论支撑与实践范式。

毛乌素沙地（图1）横跨内蒙古鄂尔多斯市、陕西榆林和宁夏盐池县，地处37°20′—39°23′N，107°20′—111°30′E。作为中国四大沙地之一，其总面积为4.22万 km²，呈现干旱半干旱大陆性气候特征，年降水量为250~440 mm，季节分布不均；年平均气温为6.0~9.0 ℃，有效积温（>10 ℃）为3 146.2 ℃，年蒸发量1 800~2 500 mm；年平均风速为2.1~3.3 m·s^-1，以西北风为主。地貌由固定、半固定、流动沙丘构成，总体地势自西北向东南倾低。由于土壤养分含量低，抗侵蚀能力差，毛乌素沙地长期处于植被稀疏、沙漠化严重的状态，为典型的荒漠化地区。植被主要包括半灌木［黑沙蒿（Artemisia ordosica）、沙柳（Salix cheilophila）等］、多年生草本［沙鞭（Psammochloa villosa）、赖草（Leymus secalinus）等］和一年生植物［雾冰藜（Grubovia dasyphylla）、角蒿（Incarvillea sinensis）等］。

本研究所采用的数据包括1990、1995、2000、2005、2010、2015、2020、2023年土地利用数据，数据来源于中国科学院资源环境数据平台（https：//www.resdc.cn），空间分辨率为30 m×30 m；参照《土地利用现状分类》（GB/T21010—2017）的分类方法，并结合毛乌素沙地景观特点，将土地利用类型分为耕地、草地、林地、水域、建设用地和未利用地；高程和坡度数据则基于地理空间数据云（https：//www.gscloud.cn）获取；气温和降水量数据来源中国气象服务中心（http：//data.cma.cn）；归一化植被（NDVI）和夜间灯光数据来自中国科学院资源环境数据平台（https：//www.resdc.cn）；人口密度数据来源于Worldpop （https：//www.worldpop.org/）；GDP密度数据来自《中国县域统计年鉴》。

本文采用网格法对毛乌素沙地景观生态风险进行评价，风险评价小区样本面积应为研究区平均景观斑块面积的2~5倍^［20-21］较为适宜，本文结合实际情况采用3 km×3 km划分网格将毛乌素沙地划分为4 483个风险评价小区。

采用景观干扰度、破碎度、分离度、优势度、脆弱度、损失度指数构建沙地景观生态风险评价指标体系^［13-14］。其中，景观生态风险指数计算公式为：

式中：ERI为研究区景观生态风险评价指数；A_ij 为第j个风险评价小区第i类景观的面积；A_i 为第i个风险评价小区的面积；R_i 为第i类景观的景观损失度指数。

景观干扰度指数（E_i ）表示不同景观所代表的生态系统受到外部干扰的程度。以景观格局分析为基础，通过景观破碎度、分离度和优势度指数叠加来反映不同景观所代表的生态系统受到干扰的程度。

式中：C_i 为景观破碎度指数；N_i 为景观分离度指数；D_i 为景观优势度指数。a、b、c分别为C_i 、N_i 和D_i 的权重，且a+b+c=1，依据相关参考文献［22］，对a、b和c分别赋以0.5、0.3和0.2的权值。

其中，景观破碎度指数（C_i ）公式为：

式中：n_i 为景观类型i的斑块数；A_i 为景观类型i的总面积。

景观分离度指数（N_i ）公式为：

式中：n_i 为景观类型i的斑块数；A_i 为景观类型i的总面积；A为景观总面积。

景观优势度指数（D_i ）公式为：

式中：Q_i =景观类型i的斑块数/斑块点数；M_i =景观类型i所占样方数/斑块总数；L_i =景观类型i的面积/总面积。

景观脆弱度指数（F_i ）：表示不同景观所代表生态系统内部结构的易损性。针对毛乌素沙地实际情况，在他人研究成果^［23］的基础上，将研究区脆弱度分为6个等级，从低到高依次为：建设用地、林地、草地、耕地、水域和未利用地，为消除不同量纲对数据的影响，对其进行了归一化处理（表1）。

景观损失度指数（R_i ）：将各指数叠加，用来表示在受到自然和人为双重干扰时各景观类型所代表的生态系统自然属性损失程度，其公式为：

式中：E_i 为景观干扰度指数；F_i 为景观脆弱度指数。

地理探测器优势在于解释性，特别是针对分析空间分异性和量化驱动因素方面，利用单因子空间分异性，可有效识别多因子之间的关系^［24］，通过叠加数据，检验数据空间分布的耦合性，探测变量数据之间的关系，本研究利用地理探测器揭示毛乌素沙地景观生态风险空间分异性的驱动过程。

单因子探测模型。本文应用单因子探测模型评估影响因子对景观生态风险的影响力。具体公式如下：

式中：SSW是层内方差之和；SST是全区总方差；N是样本总数；L是分层数。q值在（0，1］，q值大小与影响因子对景观生态风险的空间分异解释力相同，q值越大，表明该因子对沙地景观生态风险差异的影响越显著。

因子交互探测模型。主要用于探测不同影响因子之间的交互作用，在一定程度上能较好地量化分析不同因子之间交互作用对研究现象空间分布格局的影响。

结合毛乌素沙地地形地貌、自然条件和经济发展情况^［25］，本文选取的自然因素有高程、坡度、年降水量、年平均气温、归一化植被指数；根据沙地城镇化水平和经济发展情况^［25］，选取GDP密度、夜间灯光指数、人口密度和人为干扰度作为社会经济因素。

利用创建渔网形成的4 483个网格中心点作为采样点，提取每个采样点对应的空间X和Y坐标属性值；利用重分类工具中的自然断点法，将影响因子X_i 按照其大类标准分为9类。

人为干扰度是反映人类活动对生态环境的影响强度，通常基于与人类活动相关的指标，如土地利用变化、建筑物密度和工业活动等^［26］。本文基于景观类型面积计算人为干扰度，具体公式如下：

式中：H_i 为景观类型i的人为干扰度；E_i 为景观类型i的干扰度值；A_i 为景观类型i的面积；A为景观生态风险管理小区的总面积。

1990—2023年，毛乌素沙地土地利用发生显著变化，表现为草地由东向西显著增加，未利用地呈由东向西显著减少（图2）。其中，草地面积变化呈减少—增加—减少变化，2023年研究区内神木市、榆阳区、横山区、靖边县、定边县、伊金霍洛旗、乌审旗东部和北部的草地面积较1990年显著增加，毛乌素沙地近30年来草地面积增加主要来源于未利用地（图3），而未利用地面积则呈现出增加—减少的变化，在榆阳区、靖边县、乌审旗和伊金霍洛旗显著减少，这一时期的变化与“三北”防护林、退耕还林等一系列生态建设项目密切相关。耕地面积表现为减少—增加；林地和水域用地整体面积占比较少且变化不明显；建设用地面积整体占比偏少，但2023年较1990年建设用地在榆阳区和神木市的面积占比有较显著增加。

沙地景观格局指数时间序列图结果表明（图4），1990—2023年，耕地、林地、草地和建设用地的各景观指数整体上变化一致，景观破碎度、分离度、干扰度和损失度整体呈下降趋势，表明这4种景观类型地域分布趋于集中，抵御外部干扰力增强，生态风险变小，景观优势度整体上表现为上升趋势，景观优势逐渐增加；水域和未利用地各景观指数整体上变化趋势相同，景观破碎度、分离度、干扰度和损失度整体呈上升趋势，这表明两种景观类型的地域分布趋于分散化，随机散落分布现象不断加剧。同时，景观优势度整体表现为下降趋势。景观脆弱度指数来源于不同等级景观类型的归一化处理，故数值多年保持一致。

本研究ERI值1990年为0.0350~0.2190，1995年为0.0283~0.2053，2000年为0.0261~0.2113，2005年为0.0255~0.2046，2010年为0.0242~0.2239，2015年为0.0240~0.2127，2020年为0.0240~0.1848，2023年为0.0239~0.2168。根据沙地各风险小区ERI所处的范围，利用ArcGIS 10.8提供的自然断点法将生态风险划分为5个等级：0.02~0.04为低风险区，0.04~0.05为较低风险区，0.05~0.07为中等风险区，0.07~0.11为较高风险区，0.11~0.22为高风险区。从时间上看（图5），整体上沙地ERI有显著改善，高风险区面积减少不明显，但较高和中等风险区面积均有显著减少，而较低和低风险区面积显著扩大。分时段看，1990—1995年变化最为明显，较低风险区面积增加10 026 km²，中等、较高和高风险区面积分别缩小5 391、7 767、441 km²；1995—2005年，中等风险区面积变化最大，分别减少3 690、2 655 km²；2005—2020年，低风险区面积变化最大，分别增加4 284、4 158、2 475 km²；2020—2023年，各等级景观风险区域变化不大。

从空间上看（图6），1990年毛乌素沙地未利用地广布，生态风险程度高；较高、中等风险区分布变化最为显著，由沙地整体集中分布变为以该区西部为主的零散化分布；随着“三北”防护林等一系列生态措施实施，2023年该地域以草地分布为主，生态风险程度低，脆弱度低，生态系统稳定性显著增加，人为修复措施显著；未利用地景观以沙地为主，现存面积较大，沙地主要分布在鄂托克旗和鄂托克前旗，这与2023年高风险区域分布吻合。

为明确各等级生态风险转换关系，本文采用生态风险转移矩阵分析各生态风险等级转移情况（图7、表2），除各等级风险区自身不变外，生态风险转化类型有：低-较低、低-较高、较低-低、较低-较高、中-低、中-较低、较高-低、较高-较低、较高-中、较高-高、高-中、高-较高12种类型。研究期内沙地生态风险等级由中等向低转化最为显著，转化总面积为13.88×10³ km²，占沙地总面积的34.16%，由较高、高向较低、低转化总面积为9.4×10³ km²，占沙地总面积的23.14%。由中等、较高和高转向低、较低等级沙地总面积为57.30%，说明近年来毛乌素沙地实施的林草地封育保护、补植补种和防风固沙等措施，使沙地半数以上的生态环境得到改善，生态风险显著降低，但由于沙地原本生态环境脆弱性导致小部分生态保护效应的滞后性，由低、较低转向较高、高等级总面积为0.19×10³ km²，占总面积的0.48%。

生态风险等级面积净变化整体上低等级面积变化最显著，增加面积为20.39×10³ km²，较低等级面积增加9.61×10³ km²，中等级风险面积减少最多为19.59×10³ km²，较高、高等级面积减少面积和为10.41×10³ km²。可看出生态风险呈由高向低转化趋势，生态风险表现出逐步降低趋势。

景观生态风险测算结果表明毛乌素沙地生态风险演变主要受自然因素和人为因素的影响，为进一步明确沙地景观生态风险变化的影响因素，结合沙地实际情况，选择高程、坡度、年降水量、年平均气温、归一化植被指数、GDP密度、夜间灯光指数、人口密度和人为干扰度等9个指标进行定量分析。

利用地理探测器分析不同因素对沙地景观生态风险的影响程度及各因素之间的相互作用关系，得到毛乌素沙地景观生态风险驱动因素单因子探测结果及贡献度（表3）。1990—2023年，在自然因素中，归一化植被指数对毛乌素沙地景观生态风险因素贡献率最显著，均在30.00%左右，第二、三主导驱动因素是年降水量、年平均气温，坡度对沙地景观生态风险驱动因素贡献度呈下降趋势，说明“三北”防护林、退耕还林还草、天然林资源保护等相关生态保护措施工程在毛乌素沙地取得了显著成果；在社会经济因素中，人口密度对沙地景观生态风险驱动因素单因子贡献度呈上升趋势，而GDP密度贡献度呈下降趋势，说明人为活动对沙地生态景观风险影响显著，但是随着相关生态修复措施的实施，降低了人类活动对沙地生态系统的风险。

整体上看，研究期内自然因素均比社会经济因素解释力更强，归一化植被指数、高程和年降水量这些景观分布格局的主要影响因子具有较强解释力，说明自然景观分布对沙地景观生态风险分布具有较强影响力，其对沙地景观生态风险的变化起主导作用。

选取9个指标进行单因子探测，结果表明1990—2023年不同指标对毛乌素沙地景观生态风险影响存在显著差异，为进一步探究在自然因素和社会经济因素演进过程中是否存在多元因素共同作用，对9个指标进行交互因子探测（图8）。相较于单因子探测结果，各驱动因素交互作用对景观生态风险解释力呈增强态势，表明毛乌素沙地景观生态风险发展受多因素综合影响。各因子间交互作用均表现为双因子增强和非线性增强两种类型。图8表明：1990年，X₁∩X₃（0.280）、X₁∩X₄（0.280）、X₁∩X₆（0.259）、X₁∩X₇（0.259）交互作用值相对较高，对景观生态风险变化有较强解释力，说明自然和人类活动对景观分布有重要影响，X₁∩X₃、X₁∩X₄交互作用的q值最高；1995年，解释力排位前3的是X₁∩X₃（0.279）、X₁∩X₄（0.279）、X₁∩X₆（0.259），与1990年相比，整体上看大部分交互值均有提升；2000年，解释力排位前3的是X₀∩X₄（0.419）、X₀∩X₆（0.393）、X₀∩X₇（0.404），与1995年相比，人类活动对景观生态风险分布解释力减弱；2005年，解释力排位前3的是X₀∩X₁（0.394）、X₀∩X₆（0.390）、X₀∩X₇（0.394），该时期坡度、GDP密度、人口密度和夜间灯光交互作用解释力有所减弱；2010年和2015年因子交互影响q值前3的均是X₀∩X₁（0.416、0.377）、X₀∩X₄（0.425、0.390）、X₀∩X₆（0.423、0.369），其中2010年是研究期间交互因子解释力q值最高年份；2020年和2023年因子交互影响q值前3的均为X₀∩X₁（0.376、0.397）、X₀∩X₃（0.367、0.389）、X₀∩X₄（0.369、0.385），与1990年相比，人为活动对沙地景观生态风险分布交互因子解释力有明显提升。

1990—2023年整体上，NDVI（X₀）与其他因子交互作用解释力最为显著，其中NDVI与人为干扰度、NDVI与人口密度解释力均呈先增加后减少的趋势，且均在2010年为最高值X₀∩X₆（0.423）、X₀∩X₈（0.352），该时期是“三北”防护林体系四期工程的收官年份，进一步验证了政策干预对脆弱生态环境修复起到的重要作用。

1990—2023年，毛乌素沙地呈现显著景观转型，未利用地面积缩减，主要转化为草地类型，辅以耕地扩张与建设用地增长，这与修晓敏等^［27］研究结果一致。空间分异特征显示，草地、耕地和未利用地是主要的景观类型，耕地主要分布在沙地东部和南部，未利用地主要分布在沙地西部，主要是由于西部位于鄂托克旗、鄂托克前旗，在两旗有多处连绵的高大沙丘；草地面积占沙地总面积最大，分布呈分散状，Gao等^［28］研究指出，通过“三北”保护林工程、新一轮“退耕还林还草”工程、京津风沙源治理和沙化土地禁养地试点工程，沙地生态修复方面取得了重大成果。总体来看，景观类型的变化反映了自然环境和政策干预对景观格局的影响，Mu等^［29］研究指出政策对内蒙古景观格局变化有显著影响，而景观格局变化又是导致景观生态风险变化的重要原因，表明制度性干预对脆弱生态系统的重塑效能。

自1978年实施各项生态建设以来，经过多年防沙治沙工作，目前毛乌素沙地植被覆盖率达到80%^［30］，动植物生物链和种群得到有效恢复，生态环境改善也带动了沙产业建设，设有沙地生态修复示范区，已成为国家沙地生态修复示范的典型^［31］。

从时间上看，草地分离度、干扰度和损失度各时期均为最高，均呈上升趋势，未利用地各指数在各时期变化都较显著；研究期间沙地内景观生态风险有显著改善，较高风险区域面积整体上呈缩小态势，较低风险区域则呈扩大态势；从空间上看，中等风险区面积显著缩小，低风险区显著增大。高风险区分布在地势相较陡峭区域，低风险区主要分布在沙地东部和南部，地类景观以草地为主。随时间变化，高风险区主要由陕西的张家畔镇、东坑镇、堆子梁镇、石洞沟乡和安边镇、鄂尔多斯市的尔林兔镇转移到鄂托克旗、鄂托克前旗，呈现出由沙地南部、北部向沙地西部转移趋势，高风险和景观生态风险等级变化活跃区域集中分布在人类活动较强、自然地势陡峭的区域，Zhao等^［32］研究指出高风险区域分布受地貌影响，广泛分布在荒漠区和无植被覆盖的区域，这与毛乌素沙地高风险区分布特征一致；低风险、较低风险区域自然生态环境较好，这些区域以草地和林地为主，Wang等^［33］研究表明低风险区域主要分布在草地、林地，表明风险降低区域主要发生在植被覆盖度相对较高、自然条件良好、生态系统稳定且完整性高的区域。

从单因子探测结果看，自然因素是影响沙地景观生态风险变化的主导因素，1990—2000年，影响沙地的主导因素是年降水量和年平均气温，2000—2023年则是植被覆盖度，这与Yan等^［34］和Li等^［35］的研究结果一致，表明年降水量、年平均气温和植被覆盖度是影响景观生态风险变化的重要因素。社会经济因素中，人口密度对沙地景观生态风险变化影响较显著，Yao等^［36］研究结果佐证本条结论，该地区多种植经济作物、从事畜牧养殖活动和参与沙地治理^［37］，景观类型随人类生命周期过程而变化，进而影响景观生态风险程度的变化。

从交互因子探测结果看，与单因子作用相比，因子交互作用效果显著，尤其是归一化植被指数与其他影响因素的相互作用，Chang等^［38］认为归一化植被指数同其他影响因子间的交互作用解释力较强。与单一的自然因素作用相比，自然因素、社会经济因素两者交互作用解释力更强，这与Xia等^［39］研究结论一致，荣维康等^［40］指出植被恢复对景观生态风险具有重要作用。国家实施相关工程修复毛乌素沙地的生态环境（图9），使得植被覆盖大幅提升、水土保持与沙地稳定和生物多样性恢复等成效显著，促进了当地经济发展与产业转型，使得农牧民增收、增加就业机会，因此植被与人为干扰交互作用对景观生态风险的影响显著。

根据上述结果，应根据不同风险区域的特点（表4），结合当地环境条件和社会发展制定景观生态风险预防策略，对于高风险区，即沙地西部，主要受自然条件主导，通过建立自然保护区、采取防沙治沙措施等来实施沙地的生态安全屏障保护工作，可缓解气候变化造成的生态退化，并降低该区景观生态风险等级。在沙地东部的低风险区域，景观生态风险与人类活动相关的趋势呈增加态势，要确定合理的生态红线，合理调配人、生产和生态空间，协调减少人类对生态环境的影响。

草地是毛乌素沙地的优势景观类型，呈分散状分布，2023年较1990年面积增加5.70×10³ km²。沙地景观类型以低风险等级为主，2023年较1990年低风险区面积增加2.04×10⁴ km²，整体上景观生态风险呈降低态势，其中1990—1995年较高风险区在乌审旗、神木市和榆阳区变化较为显著，面积减少7.77×10³ km²，但从1995年到2010年较高风险区在乌审旗西部、鄂托克旗和鄂托克前旗呈上升趋势，从2010年开始呈下降趋势；高风险区呈下降趋势，2023年集中分布在鄂托克旗和鄂托克前旗，2023年较1990年面积减少0.05×10³ km²，较低风险区呈散状分布于毛乌素沙地。自然因子中，归一化植被指数（NDVI）贡献率最高（29.8%~33.3%），年降水量次之（5.3%~15.2%）；社会经济因子中，人口密度贡献率由4.5%升至12.9%，GDP密度影响减弱；自然-人为因子交互作用显著增强风险解释力，表明植被恢复与地形-气候协同效应主导生态风险演变。