2000-2015年中国干旱半干旱地区沙漠化进程驱动力研究

王晓峰 ，张明明，尹礼唱，黄鹏程，勒斯木初，周潮伟

1. 长安大学地球科学与资源学院，陕西 西安 710054；2. 陕西省土地工程重点实验室，陕西 西安 710054；3. 兰州交通大学测绘与地理信息学院，甘肃 兰州 730070

中国是全世界受沙漠化危害严重的国家之一（王涛等，2001），2015年国家林业局公布的《中国荒漠化和沙化状况公报》显示，中国沙漠化防治形势依然十分严峻。沙漠化进程具有复杂性，并且影响着人地耦合系统的可持续发展，其驱动力研究一直是国内外学者关注热点。研究表明，埃及El-Dakhla绿洲沙漠化主要受气候因素的影响（Ismael，2015）；墨西哥中部半干旱高地沙漠化进程受人类活动主导（Becerril et al.，2015）；而蒙古HognoKhaan、印度 Dharmapuri地区沙漠化受气候因素和人类活动共同影响（Lamchin et al.，2017；Shoba et al.，2016）。中国干旱半干旱地区是连接中国与丝绸之路经济带沿线国家的重要枢纽，风蚀沙化、土地盐碱化等问题突出（李春梅等，2002），在该区域及周边地区已开展了大量沙漠化驱动力研究。结果表明，科尔沁沙地、呼伦贝尔沙地沙漠化进程主要受气候变化的影响（王永芳等，2016；李姣等，2017），亚洲干旱区东部和西部、三北防护林工程区沙漠化也对气候变化的响应敏锐（Wang et al.，2017a；黄麟等，2018）；黑河流域、青海刚察县沙漠化受人为因素主导（刘蔚等，2008；田丽慧等，2013），近年来玛曲高原沙漠化逆转主要也是人为因素所致（胡梦珺等，2016）；而黄河源地区、毛乌素沙地沙漠化受自然因素和人类活动共同影响（胡光印等，2011；Liang et al.，2016）。沙漠化驱动力有明显的时空差异，目前在中小尺度已有较好的突破，但在更大尺度的定量评估相对较少，并且已有定量研究主要采用传统的相关性（Wang et al.，2017a）和主成分分析（胡梦珺等，2016）等方法，其驱动机制有待进一步探析。

地理探测器是空间分异性探测及驱动力分析的重要方法（Wang et al.，2016a）。该方法无过多假设条件制约，克服了传统统计方法处理类别变量的不足（湛东升等，2015）。相比传统的相关性和主成分分析，地理探测器充分挖掘了地理要素同其驱动因素之间的关联信息，它不仅可以分析各个驱动因素对地理要素的解释力，而且能够揭示驱动因素两两之间的交互作用对该地理要素的影响。地理探测器早期被应用于疾病影响因子评估方面（Wang et al.，2010），近年逐渐被应用于社会、生态环境等领域（Xu et al.，2018；Shrestha et al.，2017）。Liang et al.（2016）和 Du et al.（2016）利用地理探测器研究了毛乌素沙地和晋北生态脆弱山区沙漠化的驱动力，参考其研究，在前期用MODIS合成产品提取中国北方2000年、2015年两期沙漠化土地的工作基础上（Wang et al.，2017b），以气象数据、人口、GDP、DEM、土壤类型为基本资料，利用地理探测器方法，对 2000-2015年中国干旱半干旱地区沙漠化进程驱动因素进行定量分析，揭示大尺度上沙漠化进程的驱动力，为丝绸之路经济带沿线地区防沙治沙、建设丝绸之路经济带生态屏障提供重要参考。

1 资料与方法

1.1 研究区概况

中国干旱半干旱地区位于 36°44'-49°57'N，73°26'-123°55'E，面积约243万平方千米，主要包括大兴安岭以西，昆仑山-阿尔金山-祁连山和长城一线以北的广大地区（李飞等，2011）。区内大部分地区属于温带大陆性气候，降水量稀少，植被覆盖低，沙漠、裸地、盐碱地等广布。中国八大沙漠和四大沙地皆分布于此，沙漠化问题突出。为了保证行政区划的完整性，本研究选取中国干旱半干旱地区所覆盖的西北六省区为研究区，包括新疆、青海、甘肃、陕西、宁夏、内蒙古（图1）。

1.2 数据来源及数据处理

本研究沙漠化数据来源于前期成果（Wang et al.，2017b），该数据为基于 MODIS合成产品、利用决策树提取的2000年和2015年的中国北方沙漠化土地，其Kappa系数为0.856，分类精度良好，结果准确可靠，满足本研究应用精度。参考前人有关中国北方地区沙漠化驱动力的研究成果（表1），本研究选取气象数据、人口、GDP、坡度、坡向和土壤类型为沙漠化进程的驱动因素。

本研究气象数据来源于中国气象科学数据共享服务网（http://data.cma.cn/），包括 2000-2015年研究区222个站点的日照时数、平均最高气温、平均最低气温、平均气温、平均风速、平均相对湿度和年降水量。为保证数据的完整性，个别缺测数据选用相邻两年的平均值替代。人口、GDP和土壤数据来源于中国科学院资源环境科学数据中心（http://www.resdc.cn/），包括2000-2015年的人口公里网格数据，2000年、2005年、2010年的GDP公里网格数据和土壤类型，其分辨率均为1 km。坡度和坡向利用1 km分辨率DEM提取。

表1 沙漠化进程驱动因素选取参考Table 1 Reference of driving factors in desertification process

图1 研究区示意图Fig. 1 Map of study area

为研究气候因素对沙漠化进程的影响，选用专用于气象数据空间插值的 ANUSPLIN软件对气象数据进行空间插值处理（刘志红等，2008）。自然因素和人类活动对沙漠化地区生态系统的影响具有滞后性（陈秋红等，2008），长期数据能更好揭示地区的沙漠化驱动机制。因此，本研究计算了2000-2015年的气象插值数据、人口的平均值作为2000-2015年影响沙漠化进程的气候、人口因素。2015年GDP公里网格数据缺失，参考官方公布的中国历年GDP增长率，利用2010年公里网格数据推算出2015年GDP，并以2000年、2005年、2010年和2015年GDP的平均值作为2000-2015年影响沙漠化进程的GDP因素。

1.3 研究方法

1.3.1 沙漠化变化率

图谱分析法已成为研究沙漠化的一种重要方法（康文平等，2016）。本文首先对解译的沙漠化程度数据进行编码（表2）。

表2 沙漠化程度编码Table 2 Coding of degree of desertification

然后，计算终止年份与起始年份的沙漠化程度编码差值。计算公式如下：

式中，Sij表示i-j年份的沙漠化程度编码差值；Si、Sj分别表示i、j年份沙漠化程度编码。

为定量研究沙漠化进程的驱动力，需对沙漠化动态变化进行数量化描述。根据相关研究（康文平等，2016；Liang et al.，2016），本研究以沙漠化变化率（Desertification Change Ratio）对沙漠化动态变化进行量化表示。由于沙漠化数据分辨率为1 km，最终确定用10 km×10 km方格对沙漠化程度编码差值进行运算，得到每个方格的沙漠化变化率。每个方格为一个样本，其沙漠化变化率计算公式如下。

式中，DCR表示该方格的沙漠化变化率；s表示该方格内的沙漠化程度编码差值；Ns表示该方格沙漠化程度编码差值为 s的象元个数；smax为沙漠化程度编码差值的最大值，其为定值 4；Nsum为该方格的象元总个数，其为定值100。

1.3.2 地理探测器

地理探测器包括风险探测、因子探测、生态探测和交互探测4个模块（王劲峰等，2017），其模型定义为：

式中，qD,DCR表示驱动因素D对沙漠化变化率（DCR）的解释力，即驱动因素D对沙漠化进程的影响程度；N和σ2分别表示研究区样本总数（10 km×10 km方格数）和沙漠化变化率的方差，本研究N为21559；i表示驱动因素D划分的类/组；NDi和σ2分别表示驱动因素D中i类的样本数和沙漠i化变化率方差；m是驱动因素 D的分类/组总数。qD,DCR∈[0, 1]，qD,DCR越大说明驱动因素D对沙漠化进程的影响程度越大，反之则越小；qD,DCR=1表示驱动因素D完全决定了沙漠化进程，qD,DCR=0表示驱动因素D和沙漠化进程没有关联。

风险探测器运用t统计量检验子区域间沙漠化变化率的均值是否存在明显差异；因子探测器用于分析驱动因素对沙漠化变化率的解释力；基于F统计量检验的生态探测器用于比较D1和D2驱动因素对沙漠化进程的影响是否存在显著差异；交互探测器用于分析驱动因素之间的交互作用，即两个驱动因素共同作用时是否会减弱或增强其对沙漠化进程的解释力，或两个驱动因素对沙漠化进程的影响相互独立，判据公式为：

非线性减弱：

单因子非线性减弱：

双因子增强：

式中，∩是D1和D2驱动因素的叠加运算。

地理探测器适用于自变量（驱动因素）为类型量、因变量为数值量的情况（王劲峰等，2017）。本研究中的沙漠化数据为类型量，而驱动因素中气象数据、人口、GDP、坡度为数值量，故需将沙漠化数据转化为能反映沙漠化动态变化的数值量，而将气象数据、人口、GDP、坡度因素离散化为类型量。基于两期沙漠化土地，利用图谱分析法计算2000-2015年的沙漠化程度编码差值，进而利用沙漠化变化率计算得到反映沙漠化动态变化的数值量。对连续型的气象插值数据、人口、GDP和坡度运用非监督分类方法将其离散化为类型量，并根据相关标准对坡向和土壤类型进行分类/组处理。

结合相关研究（Cao et al.，2013；李俊刚等，2016），本研究选用等距离法（EI）、几何间隔法（GI）、自然断点法（NB）、百分位数值法（QU）、标准差法（SD）对气象数据、人口、GDP和坡度进行离散化，由于样本总量充足，故选用前4种方法将驱动因素分为4-20组，SD方法使用1、1/2、1/3、1/4倍标准差进行分组，并运用地理探测器计算分类的驱动因素对沙漠化变化率的解释力，选取解释力最大的分组作为驱动因素的最佳离散化法。随着驱动因素分类数目的增多，驱动因素对沙漠化进程的解释力逐渐趋于稳定（图2）。当驱动因素的分类/组数目增加到20以上时，驱动因素对沙漠化变化率的解释力可能进一步提高，即分析结果会进一步改善，但解释力的改进程度极小，故目前分类/组已能满足分析精度。对于日照时数、平均最低气温、平均风速，解释力最大的是GI，分20类；对于人口、GDP，解释力最大的是GI，分4类；对于平均相对湿度、坡度，解释力最大的是QU，分18类；对于平均最高气温、平均气温、降水量，解释力最大的分别是 NB、SD、GI，分别分 17类、18类（1/3倍标准差）、19类。

根据土纲将土壤类型分为 12组，分别是淋溶土、半淋溶土、钙层土、干旱土、漠土、初育土、半水成土、水成土、盐碱土、人为土、高山土和其他类型土。其中，土纲代码在 20以上的区域范围小，故将其合为其他类型土，该土壤主要为西北盐壳。坡向按照 ArcGIS坡向分类原理，分为平地、北、东北、东、东南、南、西南、西、西北9类。分组的气象数据、人口、GDP、坡度、坡向、土壤类型如图3所示。

2 结果与分析

2.1 2000-2015年沙漠化变化动态

基于中国北方两期沙漠化土地，利用图谱分析法计算沙漠化程度编码差值，进而计算得到沙漠化变化率（图4）。2000年，沙漠化面积为157.17×104km2，其中轻度、中度、重度、极重度沙漠化面积分别为 37.78×104、14.85×104、11.66×104、92.88×104km2；2015年，沙漠化面积为 147.40×104km2，其中轻度、中度、重度、极重度沙漠化面积分别为31.58×104、17.00×104、19.00×104、79.82×104km2。

2000-2015年，沙漠化面积减少了 9.77×104km2，其中轻度、极重度沙漠化面积分别减少了6.20×104和 13.06×104km2，而中度和重度沙漠化面积分别增加了2.15×104km2和7.34×104km2。整体上，2000-2015年，沙漠化动态变化面积为 80.81×104km2。沙漠化逆转面积为38.63×104km2，其中沙漠化明显逆转（沙漠化程度编码差值为-4和-3）面积为23.90×104km2，不明显逆转（编码差值为-2和-1）面积为14.73×104km2。沙漠化发展面积为42.18×104km2，其中 14.22×104km2区域沙漠化明显发展（编码差值为3和4），27.96×104km2区域沙漠化不明显发展（编码差值为1和2）。

图2 不同分类情况下驱动因素对沙漠化变化率的解释力QFig. 2 Q value of driving factors classified in different discretization method toDCR

图3 驱动因素最佳离散化分类/组结果Fig. 3 Driving factors classified bythe best discretization method

图4 沙漠化空间分布及其变化Fig. 4 Spatial distribution and change of desertification

2000-2015年，塔克拉玛干沙漠南部和北部边缘、柴达木盆地东部、巴丹吉林沙漠、阴山北麓-浑善达克沙地生态功能保护区西南部等地区沙漠化明显发展（沙漠化变化率大于0.5），毛乌素沙地西南和东北部、浑善达克沙地东部及东北部等地区沙漠化明显逆转（变化率小于-0.5）。

2.2 2000-2015年沙漠化进程驱动力分析

图5 风险探测器分析结果（平均沙漠化变化率）Fig. 5 Result of the risk detector（average DCR）

风险探测器分析结果（图5）表明：2000-2015年，在中国干旱半干旱地区，随着日照时数的增大，沙漠化有所发展，但过低的日照时数也会加速沙漠化，这可能与过低的日照影响沙区植被恢复与生长有关（曹永强等，2018）。随着气温升高，沙漠化明显发展。研究区平均风速小的地区反而沙漠化明显发展，这可能与大风日数有关（董玉祥，2001），或与过小的风速不利于植被生长有关（Wang et al.，2017a），或与该地区沙漠化受气温、降水等因素主导有关（Wang et al.，2016b）。随着降水量和平均相对湿度的增大，沙漠化程度大幅度逆转，但过高的降水量也会加速沙漠化，这可能与过高的降水量会加剧地区土壤侵蚀有关。人口多和 GDP高的地区沙漠化明显改善，这与沙漠化防治的人力与经济投入有关。坡度过低或过高的地区沙漠化加速，这与相关研究结果具有一致性（段翰晨等，2018）。缓坡地区沙漠化的发展可能与缓坡有利于风沙侵蚀物的搬运有关，而陡坡地区沙漠化的发展可能与风蚀物不易沉降有关。坡向对沙漠化进程的影响极小，平地、西北、正北坡向的地区呈沙漠化发展，其他地区沙漠化有不同程度改善。平地上沙漠化发展可能与平地更易于风沙侵蚀物的搬运有关，而西北、正北坡向沙漠化发展可能与研究区长期盛行西北风有关。盐碱土、漠土、其他类型土（主要为西北盐壳）、高山土、水成土地区呈沙漠化发展，其中盐碱土、漠土、其他类型土地区沙漠化发展明显，其平均沙漠化变化率分别为 0.1061、0.0612、0.0522。人为土、钙层土、半淋溶土地区沙漠化明显逆转，其平均沙漠化变化率分别为-0.1118、-0.1061、-0.1051，其中人为土地区沙漠化明显逆转可能与退耕还林还草政策有关。

因子探测器的q值（表3）为：降水量（0.1676）＞平均相对湿度（0.1347）＞日照时数（0.0855）＞土壤类型（0.0817）＞人口（0.0801）＞平均最高气温（0.0797）＞平均气温（0.0681）＞GDP（0.0659）＞平均风速（0.0555）＞平均最低气温（0.0475）＞坡度（0.0207）＞坡向（0.0041）。因此，2000-2015年，降水量是影响沙漠化进程的首要因素，其对沙漠化变化率的解释力为16.76%。同时，平均相对湿度和日照时数对沙漠化进程也有重要影响，仅次于降水量。土壤类型对沙漠化变化率的解释力为8.17%，其对沙漠化进程的影响也不容忽视。社会经济因素中，人口对沙漠化进程的影响大于GDP。平均最高气温对沙漠化变化率的解释力略微大于平均气温，说明极端高温对沙漠化进程的影响更明显。平均最高气温、平均气温对沙漠化进程的影响略大于平均风速、平均最低气温，而平均风速对沙漠化进程的影响略大于平均最低气温。地形因素中，坡向对沙漠化进程的影响极小。

表3 因子探测器分析结果Table 3 Result of the factor detector

2000-2015年，平均相对湿度、降水量和其他气候因素对沙漠化进程的影响有明显差异（表4）。任意两个驱动因素的交互作用都是增强型（图6），即任意两个因素的交互作用对沙漠化进程的影响都大于单个因素。其中，日照时数和人口、日照时数和GDP、平均最高气温和平均气温、平均相对湿度和降水量、平均相对湿度和人口、平均相对湿度和GDP、平均相对湿度和土壤类型、降水量和人口、降水量和GDP、降水量和土壤类型、人口和GDP、人口和坡度、人口和土壤类型、GDP和坡度、GDP和土壤类型、坡度和土壤类型呈双因子增强的交互作用，其他驱动因素均呈非线性增强的交互作用。

表4 生态探测器分析结果Table 4 Result of the ecological detector

气候因素之间的交互作用对沙漠化进程的影响显著，其中平均相对湿度、降水量、日照时数与其他气候因素的交互作用对沙漠化进程的影响更明显，即平均相对湿度、降水量、日照时数和其他气候因素共同作用时，其共同作用对沙漠化变化率的解释力大幅度提高。

3 讨论与结论

3.1 讨论

图6 交互探测器分析结果Fig. 6 Result of the interactive detector

沙漠化驱动力具有明显的时空差异性。在小尺度上，2001-2015年，呼伦贝尔沙地沙漠化主要驱动力为降水量的年际变化（李姣等，2017）；2000-2013年，科尔沁沙地沙漠化进程主要受水热因素的影响（王永芳等，2016）；2000-2015年，锡林郭勒盟地区沙漠化逆转主要驱动力为降水量的增多与国家政策的支持，而 2010年以后沙漠化的发展为人类活动所主导（王艳琦等，2018）；2000-2012年，山西北部生态脆弱山区风蚀沙漠化主要驱动因素为土壤类型、降水量和风速，地质、气候因素对沙漠化进程的影响比人类活动大（Du et al.，2016）；2000-2015年，鄂尔多斯市沙漠化逆转驱动力以自然因素为主，而人为因素为辅（冯坤等，2018）；毛乌素沙地东部沙漠化为人类活动所致，而其西部地区沙漠化为自然因素所主导（Liang et al.，2016）。在中等尺度上，1964-2014年，玛曲高原沙漠化扩张是自然和人为因素共同作用的结果，而 1998-2014年沙漠化逆转主要受人为因素影响（胡梦珺等，2016）；1975-2015年，阿拉善高原影响沙漠化进程最主要的因素为人类活动，同时气候变化也是另一重要因素（万炜等，2018）。在大尺度上，2000-2015年，“三北”防护林等生态工程实施导致风速减弱等气候因素变化从而减少土壤风蚀量的贡献约 85%（黄麟等，2018）；亚洲干旱区东部和西部沙漠化进程对气候变化的响应较为敏锐（Wang et al.，2017a）。

本研究认为，2000-2015年，中国干旱半干旱地区沙漠化进程受自然和人为因素共同作用，但自然因素为沙漠化主要驱动力，这与相关研究具有一致性（冯坤等，2018；Du et al.，2016）。中国干旱半干旱地区沙漠化对气候响应敏感，这与 Wang et al.（2017a）研究结果相吻合，其中降水量影响最大，这与相关研究具有一致性（李姣等，2017；王艳琦等，2018），说明在分析大空间尺度驱动力方面地理探测器是一种科学有效的新方法。作为一种分析驱动因子的新统计学方法，地理探测器不仅能探测驱动因素的解释力，还能分析驱动因素的交互作用（王劲峰等，2017）。中国干旱半干旱地区沙漠化进程受到平均相对湿度、降水量、日照时数与其他气候因素的交互作用更明显。前期相关学者运用该方法分别研究了毛乌素沙地、晋北生态脆弱山区的沙漠化驱动力（Liang et al.，2016；Du et al.，2016），参考其研究，本研究将该方法应用于中国干旱半干旱地区，分析了大空间尺度上沙漠化进程驱动力，可为沙漠化及其他地理要素的驱动力分析提供参考。此外，参考相关研究（康文平等，2016；Liang et al.，2016），本文提出了沙漠化变化率指数，其对沙漠化研究有一定的参考意义。

政策、科技投入也会影响沙漠化进程（张靖等，2013；宁宝英等，2018）。由于数据的限制和区域尺度的制约，本研究未能定量分析政策、科技等对沙漠化进程的影响，但人口、GDP、土壤类型对沙漠化进程的影响一定程度上也反映了政策、科技等的重要影响。此外，沙漠化对气候变化的响应具有空间差异性（王永芳等，2016）。由于方法的局限，本研究未能分析沙漠化对气候因素响应的空间差异。本研究结果可为丝绸之路经济带沿线地区防沙治沙、建设丝绸之路经济带生态屏障提供重要参考。

3.2 结论

（1）2000-2015年，沙漠化面积减少了9.77×104km2，其中轻度、极重度沙漠化面积分别减少了6.20×104、13.06×104km2，而中度和重度沙漠化面积分别增加了2.15×104km2和7.34×104km2。整体上，沙漠化逆转面积为38.63×104km2，其中沙漠化明显逆转面积为23.90×104km2，不明显逆转面积为14.73×104km2，而沙漠化发展面积为 42.18×104km2，其中14.22×104km2区域沙漠化明显发展。

（2）2000-2015年，沙漠化进程是自然和人为因素共同作用的结果，但以自然因素为主要驱动力，其中沙漠化对气候因素响应尤为敏感，而降水量对沙漠化进程影响最大，平均相对湿度的影响次之。整体上，影响土地沙漠化的驱动力依次为气候因素、土壤类型、社会经济因素、地形因素。

（3）任意驱动因素的交互作用对沙漠化进程的影响均大于单一驱动因素。气候因素之间的交互作用对沙漠化进程影响显著，而平均相对湿度、降水量、日照时数与其他气候因素的交互作用更明显。