土地沙化是中国北方的主要生态环境问题，其引起的土地资源减少、土壤质量下降以及沙尘暴等问题，直接影响区域生态系统的可持续发展^［1-2］。第六次全国荒漠化和沙化监测结果显示，截至2019年，全国沙化土地占国土面积的17.58%，具有明显沙化趋势的土地占2.91%。沙区生态状况呈现整体好转、改善加速态势，但土地沙化形势依然严峻，部分沙区沙化程度仍在加剧。随着全球气温升高、降水量减少、极端天气频发以及人类过度垦荒等影响^［3］，中国北方农牧交错带土地沙化问题日趋严重，尤其是东部科尔沁沙地和浑善达克沙地，其位置相邻、土地沙化相对集中，是中国当前防沙治沙的主要阵地。科尔沁-浑善达克沙地生态系统脆弱，频繁的沙尘暴限制了区域社会经济的可持续发展。因此，深入研究区域土地沙化时空变化格局及其主要驱动机制，对制定更加科学和有效的防治措施具有重要意义。

目前，区域尺度上土地沙化遥感监测方法主要包括遥感分类后比较法和基于遥感指数的趋势分析法等^［4］。分类后比较法一般通过对比两个时期土地沙化程度的变化来揭示土地沙化过程^［5］。基于遥感指数的趋势分析法中，常用的监测指数包括植被覆盖度^［6］（Fractional Vegetation Cover， FVC）、荒漠化指数^［7］（Desertification Index， DI）、荒漠化差值指数^［8］（Desertification Difference Index， DDI）和温度植被干旱指数^［9-10］（Temperature Vegetation Drought Index， TVDI）等，其中，植被覆盖度由于指示性强、易于计算、长时间序列可获取并与生态系统服务相关等优势，在土地沙化程度研究中得到广泛应用^［7］。与遥感指数同时使用的趋势监测方法主要有一元线性回归趋势法^［11］和Sen+Mann-Kendall趋势分析^［6］，与传统的一元线性回归趋势法相比，通过Sen趋势度和Mann-Kendall趋势检验结合可提高分析的准确性。土地沙化过程驱动力一直是沙化研究的重点，传统的研究方法包括主成分分析法、相关分析法、残差法等^［7，12］，但这些方法在分析影响因素的空间差异性和量化驱动因子方面存在一定局限性。地理探测器模型可以通过研究空间分异性来定量识别各驱动因子以及它们之间的相互作用，从而对土地沙化进行解释，该模型已成功应用于中国干旱和半干旱地区土地沙化驱动力研究^［13］。

对于科尔沁和浑善达克单个沙地的研究已经得到了学者广泛的关注。研究发现科尔沁沙地1975—2013年经历了发展-逆转-发展-逆转的波动变化过程^［14］，2000年以来，植被得到了一定程度的恢复和改善，沙化趋势得到一定的遏制^［15］。浑善达克沙地植被对气候变化敏感，尤其是生长季早期，1982—2015年沙地植被总体呈恢复趋势，沙地植被改善面积大于退化面积^［16］。20世纪80年代末期到2000年是土地沙化快速发展期，2000年以后，政府加强了防沙治沙力度，实施了退耕还林还草、京津风沙源工程、禁牧等措施，对沙化逆转起到了一定的作用^［17］。已有研究多集中于局部尺度土地沙化的发展态势，但对于两沙地在整体上是否相连、土地沙化的主要驱动力等还缺乏明确的回答，而这对于当前中国东部沙地防沙治沙工作的科学开展十分重要。

因此，本研究以科尔沁-浑善达克沙地为研究区，以FVC作为衡量土地沙化程度的指标，分析2000—2020年研究区土地沙化时空变化格局及驱动力，明确两相邻沙地之间的扩张状况，为研究区生态环境保护与可持续高质量发展提供参考依据。

本研究以科尔沁沙地和浑善达克沙地涉及的所有县级行政单位为研究区。科尔沁沙地和浑善达克沙地位于中国北方农牧交错带，两大沙地横亘在内蒙古高原中东部的辽河两岸，总面积约28.09万km²，包括4个市27个县（旗）。区域地势西高东低，海拔180~1 400 m，地表富含疏松沙质沉积物。整体上属于温带大陆性半干旱季风气候，年均温1.3~7 ℃，年降水量140~500 mm，年蒸发量1 500~3 000 mm，年均风速3.2~4.5 m·s^-1。地带性土壤是栗钙土，但是受沙化影响，现分布最广、面积最大的是风沙土^［11］。

两沙地在主要植被类型上存在差异。科尔沁沙地位于蒙古高原与中国东北平原的过渡地带，是中国半干旱半湿润地区典型的农牧交错带，植被类型以草原植被为主，植被种类较为丰富。常见的灌木、半灌木植被有黄柳（Salix gordejevii）、沙杞柳（Salix kochiana）、东北木蓼（Atraphaxis manshurica）、冷蒿（Artemisia frigida）等；草本植被包括糙隐子草（Cleistogenes squarrosa）、芦苇（Phragmites australis）、狼尾草（Pennisetum alopecuroides）、牛鞭草（Hemarthria altissima）等；沙地森林植被主要有榆树疏林、旱柳和山杏等，各旗县均有成片分布。浑善达克沙地位于内蒙古高原东南部，处于中国北方典型草原区向旱作农业区的过渡地带，植被类型主要为沙地草原和荒漠植被。植被以大针茅（Stipa grandis）、克氏针茅（Stipa krylovii）、羊草（Leymus chinensis）、冰草（Agropyron cristatum）等地带性植被和榆树（Ulmus pumila）、小叶锦鸡儿（Caragana microphylla）、黄柳（Salix gordejevii）、差巴嘎蒿（Artemisia halodendron）等沙地植被为主。

本研究使用的数据主要包括NDVI、气象数据、DEM数据、土壤类型、植被类型、中国GDP空间分布千米网格数据集（GDP）和中国人口空间分布公里网格数据集（POP）7类。NDVI数据来源于美国国家航空航天局的MOD13Q1数据产品（https：//lpdaac.usgs.gov/），空间分辨率为250 m，时间分辨率为16 d。由于NDVI最大值更能反映植被状况，采用年最大值合成法（Maximum Value Composite， MVC）得到2000—2020年逐年NDVI数据。气象数据来源于中国气象数据网（http：//data.cma.cn），包括研究区内及周边154个气象站点的降水、平均气温和日照时数。DEM数据来源于中国科学院资源环境科学与数据中心（https：//www.resdc.cn/），空间分辨率为250 m。土壤类型、植被类型、GDP^［18］和POP^［19］均来源于中国科学院资源环境科学与数据中心（https：//www.resdc.cn/），空间分辨率均为1 km。以上所有栅格数据投影坐标系统一为Krasovsky_1940_Albers，空间分辨率重采样/插值为250 m。

植被覆盖度（FVC）是表征地表植被状况的重要参数，而NDVI是植被生长状态的最佳指示因子，并与植被覆盖度具有很好的相关性^［20］。目前，利用像元二分模型计算植被覆盖度的方法应用广泛且精度较高^［21-22］。本研究基于NDVI数据，采用像元二分模型计算获取2000—2020年逐年研究区植被覆盖度（FVC）。

式中：NDVIveg为植被全覆盖时的NDVI值；NDVIsoil为无植被覆盖时的NDVI值。根据前人研究^［23］，提取累计频率为5%的NDVI值作为NDVIsoil，累计频率为95%的NDVI值作为NDVIveg。

土地沙化程度分级以植被覆盖度为主要指标，土地沙化与植被覆盖度呈负相关关系，即植被覆盖度越高，土地沙化程度越低^［7］。本研究在参考已有土地沙化划分标准与方法的基础上，结合研究区生态环境特点^［24-25］，以植被覆盖度为指标将土地沙化等级划分为未沙化、轻度沙化、中度沙化、重度沙化4级（表1）。

本研究采用Theil-Sen Median趋势分析和Mann-Kendall显著性检验相结合的方法对研究区2000—2020年土地沙化的变化趋势进行分析。Theil-Sen Median是一种非参数统计的稳定趋势计算方法，该方法计算效率高，能够有效应对数据收集过程中的测量误差和不规则性，在长时间序列的趋势分析中得到了广泛应用^［26-27］。Theil-Sen Median的计算公式为：

式中：Median为取中值函数；β是通过Theil-Sen Median估计的斜率，当β>0时，时间序列呈上升趋势，反之呈下降趋势；i和j表示2000—2020年的不同年份；xi和xj为i年和j年植被覆盖度像元值。

Mann-Kendall显著性检验的计算公式为：

式中：Z为标准化统计量；n为时间序列数据的个数，本研究为21；sgn表示符号函数。当n≥8时统计量S近似正态分布，V（S）表示S统计量的方差。在α水平下，如果Z>Z1-α/2，表示无趋势的假设被拒绝，即趋势变化显著。Z1-α/2是标准正态分布函数在α置信水平下的分布表所对应的值。当|Z|>1.96时，表明趋势通过了95%的显著性检验。结合β值和|Z|值，本研究将土地沙化趋势分为明显逆转、轻微逆转、基本稳定、轻微加剧和明显加剧等5个等级（表2）。

地理探测器是探测空间分异性并揭示其背后驱动力的方法^［28］。Song等^［29］进一步解决地理探测器空间数据离散化和空间尺度效应问题，并开发了基于最优参数的地理探测器模型。本研究选取平均气温、日照时数、降水量、高程、GDP、人口密度、植被类型和土壤类型7个因子，利用参数最优地理探测器中的因子探测和交互探测，探究不同驱动因素对土地沙化的解释程度以及它们之间的相互作用。首先利用Arcmap创建渔网工具创建3 km×3 km渔网点，剔除空值点后采样各因子数据并输入地理探测器中，初设等间距、自然间断、分位数、几何间隔、标准差5种离散化方法，3~8类间断数量，进行因子探测与交互探测。其中，因子探测的计算公式为：

式中：q为驱动因子X对因变量Y的空间分布的解释力；h为驱动因子的分层；Nh为驱动因子X的层h，N为研究区的单元数；δh2和δ2分别是因变量Y在某一驱动因子的层h方差和全区域的方差。q值越大，表示自变量对因变量值的解释力越强，相反则越弱。

交互探测是量化驱动因子X1和X2的共同作用对因变量Y的解释力，识别双因子交互作用时是否存在增强或减弱解释力的现象。它的原理是分别计算因子X1和X2的q值（q（X1）和q（X2）），然后计算因子X1和X2交互作用的q值（q（X1∩X2）），再通过比较单一因子的q值和双因子交互的q值来判断双因子间交互作用特点，即q（X1）、q（X2）和q（X1∩X2）的关系，交互作用可分为非线性减弱（q（X1∩X2）<min［q（X1），q（X2）］）、单因子非线性减弱（min［q（X1），q（X2）］<q（X1∩X2）<max［q（X1），q（X2）］）、双因子增强（q（X1∩X2）>max［q（X1），q（X2）］）、独立（q（X1∩X2）=q（X1）+ q（X2））、非线性增强（q（X1∩X2）>q（X1）+q（X2））。

为了分析科尔沁沙地和浑善达克沙地土地沙化的时空变化特征以及探讨两相邻沙地是否呈现相连趋势，本研究在整个研究区和沙地核心区条带两个尺度上分析了不同等级沙化土地的空间分布、面积及占比。研究区为科尔沁沙地和浑善达克沙地覆盖的所有行政范围；沙地核心区条带，即基于沙地核心分布区划定一个宽200 km、长950 km的条带，以两沙地的中点克什克腾旗为中心，以50 km为单位，分别向东部科尔沁沙地和西部浑善达克沙地方向共划分19个单元（图1）。

2000—2020年研究区土地沙化程度呈现西重东轻的分布格局（图2）。2000年研究区土地沙化面积为22.43万km²，占研究区总面积的79.86%。重度、中度、轻度、未沙化面积占比分别为18.98%、28.64%、32.24%、20.14%（表3）。重度沙化主要分布在研究区西部与东南部，具体包括苏尼特右旗、二连浩特市、苏尼特左旗、阿巴嘎旗和锡林浩特市；中度沙化与轻度沙化在各个旗县内零散分布（图2A）。到2020年，研究区沙化土地面积减少至16.03万km²，占研究区总面积的57.08%。重度、中度、轻度、未沙化面积占比分别为6.33%、21.38%、29.37%、42.92%。重度沙化土地在空间上呈明显收缩态势（图2B）。

研究期间内，典型沙地条带上不同程度的沙化土地具有显著的空间分布规律，两期沙化土地面积占比均由最西端的100%降低至中心条带达到最低值，继而向东增加至200 km处，随后又呈递减趋势（图3）。中心条带位于克什克腾旗中东部，2000—2020年，沙化土地面积占比由37.70%降低为16.25%，一直存在沙化现象，主要为轻度沙化。重度沙化地区主要在中心条带向西300~450 km，中度沙化主要在中心条带向西150~300 km，中心条带向西150 km以东均以轻度沙化为主（图4）。

2000—2020年，研究区土地沙化表现出显著逆转趋势。平均植被覆盖度从2000年的36%增加到2020年的53%，年均增加0.55个百分点（P<0.01）；最大值出现在2018年，为57%，最小值出现在2009年，为33%。从两个沙地来看，科尔沁沙地植被覆盖度均值在45%~67%波动（图5），同期沙化逆转趋势更为显著，年均增加0.75个百分点（P<0.01）；浑善达克沙地植被覆盖度均值为25%~43%，年均增加0.39个百分点（P<0.1）。

研究时段内区域沙化逆转面积显著大于沙化加剧面积，沙化逆转（β≥0.0005）面积为20.67万km²，占研究区总面积的73.59%。沙化明显逆转的区域面积为8.76万km²（31.17%），主要分布在科尔沁沙地东部地区和多伦县；沙化轻微逆转的区域面积为11.92万km²（42.42%），广泛分布于浑善达克沙地中东部和科尔沁西北部；沙化基本稳定的区域面积为4.28万km²（15.24%），集中分布在浑善达克沙地的苏尼特右旗、二连浩特市和苏尼特左旗；沙化加剧的区域面积为3.14万km²（11.17%），零星分布在库伦旗、阿鲁科尔沁旗、巴林右旗、苏尼特左旗东南部和克什克腾旗西南部。从空间分布上看，东部科尔沁沙地的土地沙化动态变化相对于西部浑善达克沙地表现出向好的趋势（图6）。

从沙地条带上来看，中心条带以西的浑善达克沙地段，植被覆盖度均无显著增长趋势；中心条带土地沙化有显著逆转趋势，植被覆盖度年均增加0.35个百分点（P<0.05，图7）。中心条带以东的科尔沁沙地段，土地沙化逆转更为显著；中心条带以东150~450 km处植被覆盖度均呈显著上升趋势且增长速率较高，年均增加0.76~1.08个百分点（P<0.05）。

从单因子探测结果来看，各因子对土地沙化的解释程度均有显著影响（P<0.001，图8A），因子探测器的q值为降水量（66.87%）>人口密度（51.26%）>GDP（43.65%）>高程（39.08%）>日照时数（38.95%）>土壤类型（33.64%）>植被类型（21.35%）>气温（19.33%）。因此，2000—2020年，降水量是影响土地沙化的首要因素，人口密度对土地沙化也有重要影响，仅次于降水量。GDP、高程、日照时数、土壤类型对土地沙化的影响居中，植被类型和气温对土地沙化的解释力较小。

交互探测结果表明，所有因子的交互作用均为双变量增强，任何两个因子之间的交互作用对土地沙化的解释力（q值）都大于任何单个因子（图8B）。这说明影响土地沙化的驱动因子不独立，它们之间相互作用时具有较强的解释力。降水量与其他各因子的交互作用和人口密度与其他各因子的交互作用对土地沙化的影响最为明显。降水量与土壤类型（72.13%）的交互作用最大，其次为降水量与人口密度（71.19%）、降水量与气温（71.07%）、降水量与气温（70.93%）、降水与GDP（70.45%）、降水与高程（69.87%）、降水与植被类型（68.84%）、人口密度与土壤类型（64.43%）、人口密度与日照时数（62.53%）、人口密度与气温（61.83%）、人口密度与高程（60.12%）。虽然气温、日照时数、高程、GDP、植被类型、土壤类型之间的交互作用增强，但它们的解释力相对较低。

2000—2020年，研究区土地沙化整体呈显著逆转的趋势，西部浑善达克沙地土地沙化状况较东部科尔沁沙地更为严峻。浑善达克沙地是中国北方沙尘暴的主要策源地，由于早期过度开垦、放牧和乱砍滥伐导致区域土地严重沙化^［30］。同丽嘎等^［8］发现，1985—2017年浑善达克沙地土地沙化面积呈现先增长后下降的趋势，整体沙化情况发生好转，其中重度沙化土地面积明显减少。科尔沁沙地地处半干旱半湿润区，原生植被多为由稀疏的耐旱树种组成的疏林草原，其对于维持科尔沁沙地生态平衡有重要作用^［31］。由于过度放牧、开垦和樵采等，该区顶级植被群落类型——疏林草原植被受到破坏，诱发了不同类型以及不同程度的土地沙化^［32］，段翰晨等^［33］发现科尔沁沙地在2000—2010年经过治理后，土地沙化总体呈现逆转趋势，沙化面积也有所减少，但轻度沙化土地有向南移动的趋势。政府部门采取“禁牧休牧”和“草畜平衡”政策，实施诸如三北防护林、京津风沙源治理、退耕（牧）还林还草等重点林业工程建设，天然草原退化和沙化趋势得到有效缓解^［34］。王涛等^［35］研究科尔沁地区沙化过程驱动因素，认为自2000年以后国家实施生态环境保护政策，加强对区域生态环境的保护和重建力度，使得林地和草地面积增加，土地沙化面积明显减少，加之载畜量的降低，减轻了对草场的破坏，沙漠化呈现逆转趋势。本研究还发现科尔沁沙地和浑善达克沙地相接处中心条带地区经过21年的治理后，仍存在沙化现象，沈贝贝等^［6］对内蒙古草原植被覆盖度时空分布的研究结果也呈现出类似的空间分布特征。

本研究发现科尔沁-浑善达克沙地土地沙化过程受自然和人为因素共同作用，其中降水量和人口密度是影响区域2000—2020年土地沙化空间格局的主导因子。科尔沁沙地和浑善达克沙地处于干旱半干旱地区，降水是土壤水分补给的主要来源，是沙地植被正常生长和退化植被恢复重建的主要限制因子^［6］。随着牧区人口数量的不断增长，农牧民主要依靠农牧业生产获得最基本的生活资料，在沙地单位面积产量较低的条件下只能扩大开垦面积，导致土地压力增大，生产能力降低；牧区则盲目扩大牲畜量，造成严重超载的草场退化，使得沙地生态环境不断恶化^［8］。但另一方面，人口规模的增长对沙化逆转产生了积极的影响，人口密集和经济发达的地区沙化明显得到改善，这与对土地沙化防治的人力投入和退耕还林还草等政策有关^［36］。降水量与土壤类型的交互作用对土地沙化的影响最大，因为它们共同决定了土地的抗风蚀能力，不同土壤类型的结构、质地和理化性质存在较大差异，影响其肥力、含氧量高低以及其内部有机质含量与水分利用效率^［37］，降水量直接影响土壤的湿润程度和植被覆盖情况，在降水量较少、植被相对较少的地区，土地更容易受到沙化的影响^［38］。由于数据的限制和区域尺度的制约，本研究未能定量分析政策、重点林业工程建设等对土地沙化的影响，但人口密度、GDP、土壤类型对土地沙化的影响一定程度上也反映了政策、重点林业工程建设等的重要性。

为减轻风沙危害，打好打赢两大沙地歼灭战，统筹山水林田湖草沙一体化保护和修复，筑牢北方重要生态安全屏障，需要采取科学有效的防沙治沙措施，因地制宜、科学推广应用行之有效的治理模式。由于浑善达克沙地的脆弱性和可逆性，在治理的同时必须加大保护力度。针对草地为主的区域，可以采取人工造林种草、封禁保护和埋设沙障等综合措施，在牧区制定合理的放牧方式，采取禁牧或分区轮牧、限制超过核定的适宜载畜量等措施^［39］。而科尔沁沙地，在生态系统恢复过程中需要加强对土地利用的合理区划，对以耕地为主的区域应优化生态环境良性循环的种植结构、改进作物耕作方式、推广土壤改良技术、实施黑土地保护性耕作技术并采用高效节水设施灌溉。同时禁止发展高耗水、粗放经营的农牧林草种植业，推动生态与经济协调发展。优化农、林、牧业结构，根据当地实际情况发展沙地经济和林下经济，推广绿色防虫方法代替农田化学药剂使用。积极推动沙地治理、生态转化和科技应用的有机结合，推进山水林田湖草沙一体化保护和系统治理^［40］。同时，应当加强对科尔沁沙地和浑善达克沙地相接处的生态保护，控制沙化发展。在治理过程中，需重点考虑降水量和人口密度的影响，制定相应的政策措施，加强对资源环境承载力的科学评估。加强跨区域的合作与协调，实现区域间土地治理的协同发展，共同应对土地沙化问题。

科尔沁-浑善达克沙地西部土地沙化状况较东部更为严峻，土地沙化程度呈西重东轻的分布格局；2000—2020年区域沙化逆转趋势的土地面积占区域总面积的73.59%，主要分布在浑善达克沙地中东部和科尔沁沙地东南部；沙化加剧趋势的土地面积占比为11.17%，呈零星分布状态。科尔沁沙地和浑善达克沙地相接处仍存在轻度沙化现象，但呈显著逆转趋势。降水量和人口密度是影响科尔沁-浑善达克沙地2000—2020年土地沙化空间格局的主导因子；所有因子的交互作用均呈双变量增强，降水量与土壤类型的交互作用影响最大。建议在科学评估区域资源环境承载力的前提下进一步加大对浑善达克沙地的治理力度，优化科尔沁沙地产业结构，生态与经济协调发展，促进山水林田湖草沙一体化保护和系统治理。