土地沙漠化是制约区域高质量发展的重要因素^［1］。植被是陆地生态系统的主体，是其中最为活跃和敏感的因子^［2-3］。土地沙漠化的发展与治理一定程度上取决于植被的破坏和恢复，长时序的植被变化在一定程度上可以代表区域土地沙漠化的总体变化^［4-5］。在全球变化的背景下，由于植被稀疏，风沙灾害严重，风沙区的植被覆盖变化对气候变化更为敏感。在大规模的人类活动（如生态工程建设、土地复垦、过度放牧等）干扰下，气候变化和人类活动叠加使得风沙区的植被变化及驱动机制更加复杂^［6-8］。因此，开展风沙区的植被变化动态监测及影响因素评估具有重要的意义。

随着对地观测技术的快速发展，研究长时序的植被指数（如归一化植被指数NDVI）变化是监测区域植被变化的主要途径^［9-12］。受气候变化和大规模的人类活动影响，区域植被变化的驱动机制成为当前研究的热点。当前主要采用相关分析法、残差法等来研究影响植被变化的驱动因子^［13-14］。邵艳莹^［15］利用GIMMS NDVI3g数据，基于相关分析法研究了1982—2013年中国北方四大沙地的植被变化及影响因子。Yu等^［16］利用残差法分析了京津风沙源区的植被变化及驱动机制，发现降水是限制区域植被生长的主要因子，植被恢复是自然因素和人类活动共同作用的结果。然而，这些方法在对影响因素的空间异质性方面和驱动因子量化方面挖掘仍存在一定的不足。地理探测器模型可以通过研究地理现象的空间分异来分析其背后的驱动力^［17］。地理探测器共分为4部分，分别为因子探测、生态探测、因子交互作用探测和风险探测^［18］。以植被覆盖变化为例，因子探测用来分析每个影响因子对植被覆盖空间分异的单独影响^［19］；生态探测用来分析两个因子对植被覆盖空间分布的影响是否存在显著差异；因子交互作用探测用来分析两两因子之间的交互作用；风险探测用来分析每个因子的子区域之间是否对植被覆盖空间分布存在显著差异，用于探测植被覆盖最好的区域^［20］。因此，地理探测器模型成为目前研究植被变化影响因素的主要手段之一。如庞静^［21］利用地理探测器模型研究了新疆地区植被变化的定量归因，发现自然因素与人类活动因子交互作用也以双协同作用和非线性协同作用为主，二者的共同作用导致了区域植被覆盖变化。

浑善达克沙地是中国北方四大沙地之一，地处中国北方荒漠化区中部，区域生态环境脆弱，土沙转换频繁，是离北京最近的沙源^［22］。浑善达克沙地区域气候变化显著，近30年出现了气温上升、降水量波动下降的趋势^［23］。为改善区域生态环境，国家在该区实施了“退耕还林工程”和“京津风沙源治理工程”等一系列生态修复工程。大规模生态工程建设措施主要有造林工程、围栏封育、生态移民等^［24］，这些措施的实施必然会对区域的植被覆盖产生重要影响。因此，浑善达克沙地大规模的生态工程建设和气候变化为研究生态脆弱区的植被变化及影响因素提供了良好的试验区。尽管前人对浑善达克沙地的植被变化及驱动机制进行了研究^［25-28］，当前研究仍存在一定的不足：①近20年，研究区的植被是如何变化，研究的时间尺度有待进一步扩充；②自然因素和人类活动的交互作用如何影响区域植被变化，有待进一步加强。因此，基于2000—2018年的MODIS NDVI数据，结合气象数据和社会经济数据，利用地理探测器模型，选取7个因子来研究浑善达克沙地植被覆盖变化及其驱动力。其中，自然因子选择相对湿度、降水量和平均气温；地形因子选择坡度因子；在人类活动因子中，选择累计造林面积密度作为造林工程的因子，选择年末牲畜头数密度来间接作为围栏封育、禁牧轮牧政策的因子，选择人口密度间接作为生态移民政策的因子。研究结果可为浑善达克沙地的生态修复和区域可持续发展提供一定的科学依据。

浑善达克沙地位于内蒙古高原东部，东临大兴安岭，南依阴山山脉，向西一直延伸到二连盆地中部，北部为平坦的干草原。基于中国1∶200万沙漠分布图，参考杨小平等^［29］的研究，划定了浑善达克沙地的边界（图1）。在行政区划上，沙地分布在内蒙古锡林郭勒盟苏尼特右旗、苏尼特左旗、镶黄旗、正镶白旗、阿巴嘎旗、正蓝旗、锡林浩特市、多伦县及赤峰市的克什克腾旗和河北省围场满族蒙古族自治县。浑善达克沙地地处东亚夏季风北缘，属于温带半干旱气候，年降水量为150—450 mm，年均温为0—5 ℃，降水稀少且气候灾害频发，大风天气频繁。年蒸发量为2 000—2 700 mm，降水与蒸发比例严重失调。干燥度为1.2—2.0，全年无霜期为104—120 d^［30］。浑善达克沙地的地势由东南向西北降低，区域内的主要生境类型有流动沙丘、半固定沙丘、固定沙丘、丘间低地和湿地。该地区主要地带性土壤为栗钙土，其次为棕钙土，非地带性土壤为风沙土。沙地西部植被较少，东部植被较茂盛，沙地的植物种类随生境的不同有较大的差异，但仍以草原植物为主，针、阔叶乔木，榆树疏林等超地带性植物明显^［31-32］。

NDVI数据为2000—2018年16 d的MODIS NDVI数据集，空间分辨率为250 m。该数据集来自美国USGS数据中心发布的MOD13Q1数据产品（https：//lpdaac.usgs.gov/）。选择MRT（MODIS Reprojection TOOL）工具完成MODIS NDVI数据的拼接与投影。在此基础上，运用最大值合成法（Maximum Value Composite，MVC）合成2000—2018年每月的NDVI数据。由于浑善达克沙地冬季植被较少，且受积雪的影响，冬季NDVI产品无法真实反映区域的植被状况。因此，参考前人的研究，选择4—10月生长季的NDVI均值进行研究^［33］。

2000—2018年地面气候月值数据来自中国气象数据共享网（http：//data.cma.cn），气象因子（表1）包括相对湿度（X₁）、降水量（X₂）和平均气温（X₃）。由于原始数据为站点数据，采用反距离权重插值法（Inverse Discharge Weight，IDW）将气象数据插值为分辨率为1 km的数据。DEM数据来自美国国家航空航天局（NASA）和国防部国家测绘局（NIMA）共同测量的SRTM（Shuttle Radar Topography Mission）90 m数据产品（http：//dwtkns.com/srtm/）。利用ArcGIS软件中3D分析数据得到了研究区的坡度数据（X₄，表1）。

人类活动因子（表1）选择人口密度（X₅）、累计造林面积密度（X₆）和年末牲畜头数密度（X₇）。人口密度、年末牲畜头数和各县面积来自2000—2018年内蒙古自治区和河北省统计年鉴。造林面积来自2002—2018年的中国林业统计年鉴。由于区域从2001年开始实施京津风沙源治理工程，累计造林面积密度（hm²·km^-2）为2002—2018年各县累计造林面积除以各县面积得到。年末牲畜头数密度（头·km^-2）为2000—2018年各县年末牲畜头数的均值除以各县面积得到。

地理探测器是探测空间分异性，以及揭示其背后驱动力的一组统计学方法^［17-18］。

（1）因子探测将植被覆盖作为因变量，各因子作为自变量，则影响区域植被覆盖分异决定值q为：

式中：q是度量空间分异性的指标；h=1，2，···，L，L为分类数目；n_h和n分别为层h和全区的样本单元数；σ_h²和σ²分别为层h和全区的方差。q的取值范围为［0，1］，q值越大，说明影响因子对植被覆盖变化的解释力越强。

（2）生态探测主要以F统计量来衡量两个因子之间是否存在显著差异：

式中：N_x1及N_x2分别表示两个因子的样本数量；L₁和L₂分别表示变量x₁和x₂分层数目。

（3）交互作用探测用于识别因子之间的交互作用，即评估自然因子和人类活动共同作用（增加或减弱）或者相互独立作用对植被覆盖空间分布的解释力。通过计算两种因子的q值之和与影响因子两两的联合q值，可识别不同影响因子之间对区域植被变化的交互作用，来判断两个因子的交互作用是增加了对植被覆盖变化的影响还是减弱了对植被覆盖变化的影响，或者两个因子是独立起作用的。“交互关系（∩）”代表两个影响因子x和y之间的交互关系，Min（q（x），q（y）），Max（q（x），q（y））和q（x）+q（y）将坐标轴分成4个区间，由q（x∩y）在4个区间中的位置确定交互关系（表2）。

（4）风险探测用于判断两个因子子区域间的属性均值是否有显著的差别，用于搜索NDVI好的区域。风险探测用t统计量来检验：

式中：Y¯h表示子区域h内的属性均值；n_h为子区域h内的样本数量；Var表示方差。

时间上，2000—2018年浑善达克沙地的植被覆盖呈波动上升的趋势，增加趋势为0.01/10a（R²=0.11，P=0.169）。最小值出现在2001年，最大值出现在2012年（图2A）。空间上，2000—2018年浑善达克沙地的植被覆盖变化呈现一定的异质性（图2B）。苏尼特右旗的东部、苏尼特左旗的西部和克什克腾旗的部分区域植被有明显退化的趋势，而阿巴嘎旗、锡林浩特市、正蓝旗、正镶白旗、多伦县和围场满族蒙古族自治县的大部分区域植被则呈增加的趋势。从图2C可以看出，浑善达克沙地显著退化的区域集中在研究区东部（即克什克腾旗），显著改善的区域集中在研究区的中部和南部（即多伦县、阿巴嘎旗、正蓝旗和正镶白旗部分区域）。研究区68.83%的区域植被覆盖在增加，其中21.19%的区域呈显著增加的趋势。研究区31.17%的区域植被覆盖在减少，其中仅3.32%的区域呈显著减少的趋势（图2C）。该研究结果与前人研究结果一致。例如，元志辉等^［26］分析了2000—2014年浑善达克沙地植被覆盖的变化，发现无论是生长季还是每月的均值，浑善达克沙地的NDVI均呈增加的趋势。Ma等^［27］也发现2000—2015年浑善达克沙地的NDVI呈增加的趋势，植被覆盖增加的区域也集中在研究区的中部和南部。

影响浑善达克沙地植被覆盖的主要影响因子分别为相对湿度（0.566）、降水量（0.495）、累计造林面积密度（0.491）、年末牲畜头数密度（0.464）、人口密度（0.235）、平均气温（0.082）和坡度（0.056）。在自然因素中，水分是控制浑善达克沙地植被生长的主要因子。相对湿度和降水量变化可解释区域植被变化的56.6%和49.5%。相对湿度、降水是驱动浑善达克沙地植被覆盖年际波动的主要因素^{［25，27］}。水分参与植被的光合和蒸腾作用等生理生化过程，水分减少会导致植被光合速率下降、有机物产量降低^［32］。在干旱半干旱地区，降水是影响植被活动的主要限制因子^［33］。在人类活动因素中，累计造林面积密度和年末牲畜头数密度可解释区域植被变化的49.1%和46.4%。2001年后，中国政府在浑善达克沙地实施了退耕还林、人工造林、围栏封育、禁牧、轮牧等政策，并实施了草原生态补奖政策来控制当地的牲畜头数，其目的是通过给农牧民发放补贴，鼓励他们减少牲畜保护草原生态系统^［34］。因此，退耕还林工程和禁牧、轮牧等措施对区域的植被恢复起到了一定的积极作用。大规模的生态工程建设对浑善达克沙地的植被恢复产生了重要的影响。人类活动对浑善达克沙地26.02%的区域植被恢复起到了积极作用^［27］。生态工程建设措施是导致京津风沙源区（浑善达克沙地是其中一部分）植被覆盖增加主要因素^［16］。张思源等^［34］利用地理探测器模型研究发现2000—2015年家畜密度对内蒙古自治区NDVI的变化的解释力有小幅度的上升。本文进一步利用生态探测器分析各个因子对区域植被覆盖变化的显著性差异。影响植被变化的前4个主要因子（相对湿度X₁、降水量X₂、累计造林面积密度X₆和年末牲畜头数密度X₇）之间没有显著差异，而两个主要人类活动因子（累计造林面积密度X₆和年末牲畜头数密度X₇）与其他3个因子（平均气温X₃、坡度X₄和人口密度X₅）之间存在显著差异（图3）。因此，相对湿度、降水量、累计造林面积密度和年末牲畜头数密度是影响区域植被覆盖变化的主要因子。

为了说明不同驱动因子之间是会增强、减弱或单独影响区域的植被变化，利用因子交互探测器识别了不同驱动因子对浑善达克沙地植被变化的交互作用。从图4可以看出，不同因子之间的交互作用q值均大于每个单独因子的q值，表明驱动因子的两两交互作用均会增加对区域植被变化的解释力。例如，降水量（X₂）与累计造林面积密度（X₆）交互作用解释力达到了54.5%，高于降水量单独影响解释力的49.5%和累计造林面积密度单独影响解释力的49.1%。

从表3可以看出，自然因素和人类活动因子之间的交互作用以双协同作用和非线性协同作用为主。其中，相对湿度（X₁）与其他因子的交互作用均以双协同作用为主。除相对湿度与坡度的交互作用外，相对湿度（X₁）与其余5个因子（X₂，X₃，X₅，X₆，X₇）的交互作用的解释力均达到了60%以上。降水量（X₂）与年末牲畜头数密度（X₇）交互作用解释力达到了62.2%，累计造林面积密度（X₆）与其余因子交互作用的解释力均达到了50%以上。因此，自然因素和人类活动对浑善达克沙地植被变化的影响既不是单方面起作用，也不是二者的简单叠加，而是相互增强或非线性增强效应。该研究结论与前人的研究类似^［34-36］。例如，张思源等^［34］利用地理探测模型分析了2000—2015年内蒙古自治区的NDVI变化及驱动力，发现不同因子之间的交互作用呈现双因子增强或非线性增强。

本文利用风险探测器分析不同驱动因子对浑善达克沙地植被生长的适宜范围，用于搜索植被覆盖较优的区域。风险探测器结果给出了某一因子不同等级之间NDVI均值及不同等级之间是否存在显著性的差异（图5）。以降水量为例，降水量被分成11个等级，数值越大表示降水量越大。从图5B可以看出，第11等级（降水量为416—435 mm）时，NDVI的值最大。随着降水量的增加，植被覆盖度呈增加的趋势。统计检验也表明，第11等级降水量的NDVI均值与其他等级的NDVI均值之间存在显著差异，进一步证明降水量为416—435 mm时，浑善达克沙地的植被覆盖最好（表4）。其他因子也做类似的分析，从而找到利于浑善达克沙地植被生长适宜的范围（表5）。从图5可以看出，随着相对湿度和降水量的增加，植被覆盖均呈增加的趋势。植被生长最适宜相对湿度范围为53.7%—58.1%，此时NDVI最大值为0.345。降水量适宜范围为416—435 mm，此时NDVI最大值为0.442。随着累计造林面积的增加，植被覆盖呈增加的趋势，当累计造林面积密度范围为31.1—38.9 hm²·km^-2时，NDVI值最大，为0.440。而随着年末牲畜头数的增加，植被覆盖则呈波动的变化，当年末牲畜头数密度为58—60头·km^-2，NDVI值最大，为0.384（图5）。随着其他因子的变化，植被覆盖呈现一定的波动变化，其中当平均气温为5.5—5.7 ℃、坡度为5.7—34.7、人口密度为50—57人·km^-2时，区域NDVI的值达到最大（表5）。

地理探测器虽有独特优势，但也存在一定不足，例如不同的离散化分类方法会对结果产生一定的影响^［30］。因此，未来需利用不同离散化分类方法分析影响区域植被覆盖变化的驱动因子。此外，不同时期年末牲畜头数受自然因素和国家政策等因素共同影响，在考虑其对植被覆盖的影响时，应分为不同时期进一步研究。同时，人类活动除生态工程建设外，还涉及其他的一些因子，未来需要加入更多的人类活动因子来分析影响植被覆盖变化的驱动因子。

2000—2018年浑善达克沙地的植被覆盖呈增加的趋势，增加趋势为0.01/10a。其中，21.19%的区域呈显著增加的趋势，集中在研究区的阿巴嘎旗、锡林浩特市、正蓝旗、正镶白旗、多伦县和围场满族蒙古族自治县的大部分区域。而研究区苏尼特右旗的东部、苏尼特左旗的西部和克什克腾旗的部分区域植被有退化的趋势。

自然因子方面，相对湿度和降水量是影响区域植被覆盖的主要因素，解释力分别达到了56.6%和49.5%。人类活动因子方面，累计造林面积密度和年末牲畜头数密度是影响区域植被覆盖的主要因素，解释力分别达到了49.1%和46.4%，表明生态工程的实施对区域植被覆盖产生了重要的影响。

自然因素和人类活动对浑善达克沙地植被变化的交互作用以双协同和非线性协同为主。其中，相对湿度和3种人类活动因子（累计造林面积密度、年末牲畜头数密度和人口密度）、降水量和年末牲畜头数密度均以双协同为主，解释力均超过了60%。而平均气温和年末牲畜头数密度以非线性协同为主，解释力也达到了60%。