土地沙化严重影响生态环境，可导致干旱、半干旱和半湿润地区土地的物理、化学和生物潜力降低^［1］。中国荒漠化土地和沙地面积分别占国土面积的25%以上和16.7%以上^［2］。土地沙化不仅侵占当地居民生存土地，也会威胁当地环境安全、经济增长，增加贫困风险。干旱半干旱地区植被稀疏、水资源稀缺、生态环境脆弱，土地沙化现象频发^［3］。土地利用类型能够反映人类活动与自然环境的交互作用。以往自然资源与环境情况决定土地利用方式，但现代人类科技进步引起改造自然的能力上升，不合理的改造方式引起资源环境的不良变化，在土地覆盖方面表征为土地沙化^［4］。从土地利用角度出发探究土地沙化空间分布特征，确定土地沙化空间分布及转换影响因子，是区域土地沙化防治的基础。

学界对于土地沙化的原因争论不止，主要的观点集中在自然主导、人为主导、自然和人类活动共同作用这3种观点上^［5-8］，多数学者认为土地沙化是气候变化和人类活动共同作用的结果。目前，常见的荒漠化驱动力研究有主成分分析法^［9-12］ 、残差分析法^［13］、层次分析法^［14］、灰色关联度分析法^［15］等。董玉祥^［12］使用主成分分析法对中国半干旱地区6个典型沙漠化区域进行研究，发现中国半干旱区的现代沙漠化是自然与人为因素共同作用的过程。Xu等^［13］使用残差分析法选取NPP为指标，利用NPP及潜在NPP与实际NPP值之差衡量了多尺度下鄂尔多斯高原沙漠化过程中人类活动与气候变化作用程度，发现1980—2000年沙漠化恢复的主因是气候因素。周华坤等^［14］利用层次分析法以江河源区高寒草地为研究对象，研究了高寒草地的退化原因，发现超载过牧和暖干化气候是导致高寒草地退化的主导因子。章予舒等^［15］使用灰色关联度分析方法以沙漠化土地面积变化为母序列计算各驱动因子关联系数，发现疏勒河流域沙漠化最主要的驱动因素是牲畜头数。使用主成分分析法可以定量研究驱动因素的贡献率，通过主成分分析数据压缩技术将多个驱动因子变量转换为少数的综合指标，分析过程需要大样本数据，分析精度受数据完整性影响^［16］。残差分析法则通过对比预测值与实际值的残差来判定人为因素的贡献度，数据以栅格单元储存，不受行政单位界限牵制，但需要调整预测参数，不同的评判基准也会引起误差^［17］。层次分析法需要专家知识判断权重，权重的设定会影响分析结果。灰色关联度分析法是以关联系数确定驱动因子影响程度的分析方法，分析过程所需样本量小，但需要现行确定各项指标的最优值，主观性过强，且部分指标最优值难以确定。对荒漠化驱动因子的研究多是单驱动因子贡献分析，对于驱动因子的交互作用研究相对较少。地理探测器是一种从空间分异性角度来探测其背后驱动力的统计方法，最初用于探究公共健康问题^［18］，目前也被广泛使用在土地利用^［19］、地质^［20］、植被覆盖^［21］、气象^［22］等方面。引入地理探测器可以从空间分层异质性定量探测影响沙地空间分布的驱动因子及各因子之间的交互作用，从而能够定量分析多因子交互作用程度。

阿拉善盟境内有巴丹吉林、腾格里、乌兰布和三大沙漠，是中国主要的沙尘源区，确定该区域的沙地时空分布格局及土地沙化驱动因子，对系统地厘清土地沙化过程及其驱动机制，对该区域科学合理进行沙漠化防治具有重要的现实意义^［23］。本文将借助地理探测器对阿拉善盟1980—2015年的沙地时空格局演变进行探究，探求阿拉善盟沙化驱动因子。

阿拉善盟地处内蒙古自治区最西部，下辖阿拉善左旗、阿拉善右旗、额济纳旗，北与蒙古国毗邻，东南临宁夏回族自治区，西南临甘肃省，地理位置37°24′—42°47′N、97°10′—106°53′E（图1），总面积23.69万km²，海拔900~1 400 m，西南高，东北低。气候干旱少雨，夏热冬寒，风大沙多，蒸发强烈，属于典型大陆性气候。年均气温6.8~8.8 ℃，年降水量40~200 mm，雨季集中在7—9月，该段时间降水占全年降水量的59%~75%。境内主要河流为黄河与额济纳河（黑河下游），土地类型主要为荒漠、戈壁，分别占土地面积的40.6%、27.4%^{［11，24］}。

阿拉善盟土地利用数据、DEM及矢量数据来源于资源环境数据云平台（www.resdc.cn），其中土地利用数据采用了1980、1990、2000、2010、2015年数据，并将原来的6个大类25个小类重分类为耕地、林地、草地、水域、建设用地、沙地、未利用地。

选取的驱动因子为9个自然因子（年均气温、年降水量、年平均风速、年日照时间、距水域的距离、土壤类型、海拔、NDVI、土壤湿度）和1个人为因子（人口密度）。其中人口密度来源于中国人口密度数据集；土壤类型利用《1∶100万中华人民共和国土壤图》数字产品；年均气温、年降水量、年平均风速、年日照时间由中国气象数据网日值数据集插值获取；NDVI由美国国家航空航天局NOAA-NDVI3G数据集^［25］获取，因采集年份受限，1980年NDVI数据使用1982年的代替；土壤湿度由欧洲空间局气象变化项目ESA-CCI土壤湿度数据集^［26-28］获取；距水域的距离利用ArcGIS10.2空间分析模块中的欧氏距离处理相关矢量数据得到。所有数据均依据阿拉善盟行政区矢量边界进行裁剪，以3 km间距采样点采样。

地理探测器对于驱动因子的输入要求为类型量，需要将所选驱动因子连续变量离散化为类型量，离散化方法如表1所示。

地理探测器是由王劲峰等^［29］提出的探测空间分异性以及解释其背后驱动力的一组统计学方法，除了对单因子贡献率识别，也可以用于识别多因子两两之间的交互作用，主要包括分异及因子探测器、交互作用探测器、风险区探测器、生态探测器4个探测器。本研究使用前3个探测器。

因子探测器。用于探测各驱动因子X对因变量Y的空间分异性的解释力度，驱动因子的影响力由q_D，U表示，q_D，U值的物理含义是驱动因子X解释了100%的因变量Y。计算公式如下：

式中：q_D，U为影响因子探测力指标；N_D，i为次一级区域样本数；N为整个区域样本数；m为变量Y或因子X的次级区域个数；δu2为研究区所有样本值的方差；δu2Di为次一级区域内样本的方差。q_D，U的取值范围为［0，1］，q_D，U值越大，驱动因子的影响力越大。

交互作用探测器。识别不同驱动因子之间的交互作用，评估两因子共同作用时对因变量Y的作用是增强或减弱。通过判断q（X₁）、q（X₂）与q（X₁∩X₂）的关系，确定影响增强和减弱的类型。

风险区探测器。判断两个子区域间的属性均值是否有显著的差别，使用t统计检验分区间因变量Y的均值是否有显著性，也可以搜索在各驱动因子分区间因变量Y的平均值，通常可以寻找均值的高低值区域。

土地利用转移矩阵是将土地利用变化类型转移面积按矩阵形式列出，可以反映不同用地类型之间的变化关系^［30］。

式中：Si,j为研究期内第i类土地向第j类土地转化的面积；n为土地利用的类型数量；i为研究期初土地利用类型；j为研究期末的土地利用类型。

研究步骤包括：①对土地利用图层进行空间分析，利用土地利用转移矩阵确定出1980—2000年与2000—2015年沙地的变化趋势，判断沙地转移方向；②将土地利用类型划分为沙地与非沙地做二值化处理，得到沙地空间分布图；计算2000—2015年沙地变化区域变化率，得到沙地变化图；③将所选定的驱动因子进行离散化分类制图；④利用地理探测器，判断各驱动因子对1980—2015年沙地空间分布影响程度、各驱动因子对2000—2015年沙地变化的影响程度。利用交互作用探测器计算两两驱动因子之间的交互作用（表2）。利用风险探测器寻找平均值超过0.5的区域，该区域沙化风险概率超过50%。

1980年沙地面积占阿拉善盟总面积的40.72%，2000年沙地面积占40.93%，2015年占40.79%，表明该区域沙地在1980—2000年略有增加，2000年以后略有减少（表3）。20世纪70—90年代，阿拉善盟人类活动增加，农林畜牧业发展，是造成1980—2000年沙地变化的主要原因^［31］。2000年后，在腾格里、巴丹吉林、乌兰布和等三大沙漠周边陆续建设了“三北”防护林工程、天然林保护工程、退耕还林工程并设立了森林生态效益补偿制度及造林补贴试点项目、沙地封禁保护补助试点等保护和建设项目^［32］，与2000年以后沙地面积减少相对应。

1980—2000年，新转入的沙地集中在阿拉善盟东南部腾格里沙漠内部水域周围，转出的沙地面积较少，无明显集聚性。1980—2000年沙地变化活跃区域集中在阿拉善左旗腾格里沙漠内部（图2）。

根据土地利用转移矩阵，1980—2000年，按照面积排序，沙地的转出土地利用类型依次为未利用地、耕地、草地，转入面积依次为草地、未利用地、水域（表4）。

2000－2015年，沙地转入转出活跃，且有一定聚集性，集中在阿拉善左旗与阿拉善右旗（图3）。该时期新转入的沙地集中在腾格里沙漠内部水域周围，与1980—2000年转入范围相似。转出的沙地集中在巴丹吉林、腾格里、乌兰布和沙漠边缘地区。阿拉善盟2010年后实施沙漠锁边工程^［32］，集中在腾格里沙漠东南部边缘、乌兰布和西南部边缘，这与沙地转移至未利用地的范围一致。沙地的转入面积依次为未利用地、草地、水域、耕地、林地、建设用地，转出面积依次为未利用地、草地、林地、水域、耕地、建设用地（表5）。

各因子对2015年沙地分布影响程度的贡献程度NDVI（18.09%）>年日照时间（17.87%）>土壤湿度（17.13%）>年均气温（11.17%）>土壤类型（9.99%）>年降水量（8.04%）>距水域的距离（7.26%）>海拔（5.7%）>风速（4.69%）>人口密度（0.49%）。NDVI、年日照时间、土壤湿度是影响沙地空间分布的主要决定因子，其次是年均气温、土壤类型、距水域的距离、年降水量、海拔、风速，人口密度对沙地空间分布影响较弱，不足1%（表6）。

1980—2000年各驱动因子q_D，U值较为相似，2010—2015年各驱动因子q_D，U值较为相似。除2015年土壤湿度的影响力略微下降，其余各年土壤湿度的影响力均居第一，土壤湿度作为地表径流、地下水量分布情况的反映，是影响沙地空间分布的主要影响因子。各年人为因子影响力小于自然因子，影响沙地空间分布的主要影响因子是自然因子（图4）。

选定的驱动因子之间的交互作用q_D，U值均大于单因子作用下的q_D，U值，因子交互结果呈现双因子增强和非线性增强（表7）。沙地空间分布影响因子交互作用解释力度大于30%的次序为X1⋂X9>X2⋂X4>X4⋂X5>X1⋂X4，表明自然因子在交互作用下增加了对沙地空间分布的影响力度。人文因子与自然因子的交互作用也增加了对沙地空间分布的影响程度，交互作用下的影响力均明显大于人文因子单独作用下的影响程度。两因子交互作用均呈现双因子增强与非线性增强效应，表明对于沙地分布影响不单单是各驱动因子独立作用结果，更多是各驱动因子相互作用影响。

将1980－2015年沙化概率大于50%的分布区域统计为沙化风险区域。从自然条件来看，沙化风险区域1980—2000年年均气温为8.25~9.28 ℃，2000年以后温度出现了上升，未来防治需要重点注意统计出的2000—2015年9.05~10.18 ℃沙化区域。年平均风速2.3~3.45 m·s^-1属于沙化风险区域，年日照时间大于3 063 h的区域沙化风险概率增大，土壤类型为初育土、水成土、人为土、养殖场的土地属于沙化风险区，容易发生沙化。海拔733~1 127 m也易成为沙化区。植被稀少（0.079<NDVI<0.109）地区更易发生沙化，土壤湿度0.1~0.23的地区也易发生沙化。从人类活动来看，1980—2000年沙化风险区分布在人烟稀少地区，2000年以后沙化区人口出现明显增长，人类活动对沙地分布产生了一定影响（表8）。

2000－2015年，影响阿拉善盟沙地变化的驱动因子影响力从大到小依次排列为人口密度（2.9%）、土壤类型（1.95%）、土壤湿度（1.66%）、NDVI（1.26%）、年降水量（1.04%）、年均气温（0.65%）、海拔（0.65%）、距水域的距离（0.56%）、年平均风速（0.46%）、年日照时间（0.38%）。研究时段沙地变化是人为因子和自然因子共同作用的结果，人口密度对沙地变化具有较强的影响力，土壤类型、土壤湿度、NDVI、年降水量影响力略低于人口密度，海拔、距水域的距离、年平均风速、年日照时间等对沙地变化的影响力微弱（表9）。

从交互探测的结果来看，影响2000—2015年沙地变化的驱动因子之间的交互作用也均大于各因子独立作用（表10）。人为因子与自然因子交互作用对沙地变化的影响力明显增强。其中人口密度与年降水量、土壤湿度、NDVI、土壤类型呈双因子增强，年均气温、年日照时间、年平均风速、海拔、距水域的距离的交互作用呈非线性增强。综上，2000—2015年沙地的变化是人为与自然因子共同作用的结果。

1980—2000年，阿拉善盟沙地以转入为主，沙地增量大于减量。2000—2015年，阿拉善盟沙地转入转出活动活跃，但减量大于增量，沙化趋势得到了遏制。转移区域集中在各大沙漠边缘，结合万炜等^［11］和周建秀等^［31］的研究结果，人类活动如建设锁边林、退耕还林、禁止放牧等是影响变化的主因。2000年后影响沙地转移的主要驱动因子是人为因子（人口密度），自然因子（土壤类型、土壤湿度、NDVI、年降水量等）的影响略低于人文因子，沙地的变化是人为与自然因子共同作用的结果。

2000年以前土壤湿度为阿拉善盟沙地空间分布的首要影响因子，NDVI、土壤类型、年平均风速、年均气温、年日照时间、年降水量为次级影响因子。2000年后除年日照时间的影响力增加、年平均风速的影响力略微下降外，其余因子与2000年前相似。沙地空间分布驱动单因子影响力中自然因子影响较大，但自然因子、人为因子的独立作用弱于彼此间交互作用的影响力，说明该区域沙地的空间分布是自然因子与人为因子协同作用下形成的。

风险探测器探测出了沙地空间分布概率大于50%的区域，后续的沙地防治区域划定与防护政策实施应该优先考虑沙化风险区域。

本文讨论了阿拉善盟1980—2015年沙地变化情况，使用地理探测器分析了影响该地区沙地空间分布、沙地变化的驱动因子。将地理探测器引入到沙地的驱动机制研究中，实现了定量分析各驱动因子影响、厘清沙地驱动机制的可能，其双因子探测器也提供了各因子交互作用下的影响力，提供了一种新的研究思路。受数据收集难度，本文选用的驱动因子缺少了放牧因子，阿拉善盟放牧对沙地产生的影响未能考虑。后期研究可以考虑纳入放牧因子，完善人为因子的选择。