科尔沁左翼中旗土地沙化时空动态变化及驱动力分析

韩阿茹汗 包斯琴 马雅婷 那日嘎拉 张佳华

(内蒙古农业大学沙漠治理学院，内蒙古 呼和浩特 010018)

土地沙化是指因风力、重力等外营力作用，导致地表组成物质中的细颗粒损失，或者原有地表因为外来沙砾覆盖等原因而导致土地退化的过程[1，2]。科左中旗境内有着中国4大沙地之一的科尔沁沙地，是防风治沙的重要区域，曾是水草丰沛的科尔沁大草原，进入21世纪以来，受人口剧增、滥垦滥牧等多种因素的影响，导致土地退化严重、水土流失加剧，生态环境急剧恶化[3]。科尔沁沙地歼灭战是习近平总书记亲自部署的“三北”攻坚战3大标志性战役之一，2023年7月31日，内蒙古自治区建设我国北方重要生态安全屏障促进条例正式发布，明确表示了土地沙化荒漠化综合防治应因地制宜，突出重点治理，全面提升生态系统的质量和稳定性，因此分析科左中旗土地沙化时空动态变化，了解土地沙化的空间分布具有现实意义。

随着遥感技术的快速发展，卫星图像成为监测和分析土地利用时空变化的重要数据支撑，谷歌地球引擎遥感大数据平台因其获取数据便捷，处理计算数据准确，已广泛成为遥感监测的重要工具[4]。目前对于土地沙化的监测方法，野外实地调查与遥感目视解译结合，这种方法耗时耗力，不适用于长时间序列的监测；多指标综合监测法，单一的地表参量指标不能够全面反映区域沙化，因此想要监测结果更加合理，应选取多个地表参量进行共同监测。如，曾永年等[5]利用归一化植被指数(Normalized Difference Vegetation Index，NDVI)和地表反照率指数(Surface Albedo)构建了特征空间，选取沙漠化遥感监测差值指数(Desertification Difference Index，DDI)作为监测指标。刘晓茜等[6]通过建立Albedo-NDVI特征空间沙化监测模型法对腾格里沙漠东南缘格林滩绿洲20年间沙化土地进行动态监测及分析。该模型充分利用了多维遥感信息，可以自动提取土地沙化信息，可以避免单一地表参量的误差，有利于土地沙化的监测与分析。

科尔沁左翼中旗作为科尔沁沙地的重要部分，亟需对土地沙化监测与驱动力分析。基于此，本文以科尔沁左翼中旗为研究区，通过提取归一化植被指数与地表反照率指数分析2000—2021年土地沙化的时空演变特征及转移特征，并利用地理空间探测器模型探究自然驱动因子和社会经济驱动因子对土地沙化的影响，进而对土地沙化提供新的治理思路，为防沙治沙提供科学依据。

1 材料与方法

1.1 研究区概况

科尔沁左翼中旗位于内蒙古自治区通辽市(E121°08′～123°32′，N43°32′～44°32′)，地处大兴安岭东南边缘，西辽河北岸，地势西北高东南低。旗内气候属中温带大陆性季风气候，气候干燥多风沙。土壤类型主要包括草甸土、风沙土、栗钙土、盐土、沼泽土，植被类型有沙生植被、草甸植被、盐生植被等。

图1 研究区地理位置

1.2 数据来源

1.2.1 遥感数据

科尔沁左翼中旗2000年、2005年、2010年、2015年、2021年的Landsat影像数据均下载自GEE(Google Earth Engine)云平台。其中，2000年、2005年、2010年为TM影像，2015年和2021年为OLI影像，分辨率均为30m。

1.2.2 驱动因子数据

本文选择自然驱动因素和社会经济驱动因素的10个驱动因子栅格数据进行土地沙化驱动研究。其中，年平均气温和年降水量数据通过国家气象科学中心(http：//www.geodata.cn)途径获得，人口空间分布、GDP、土壤类型、土壤侵蚀、平均风速、日照时数均来源于中国科学院资源环境科学与数据中心(https：//www.resdc.cn)，坡度数据使用Arcgis 10.8坡度功能计算，DEM数据来源于地理空间数据云(http：//www.gscloud.cn)。

1.3 研究方法

1.3.1 Albedo-NDVI特征空间

地表反照率(Albedo)是重要的遥感反演参数，反映了地球表面对太阳短波辐射反射特效。Albedo受到地表粗糙度、地表含水量因素影响而发生变化。参照梁顺林[7]建立的反演模型来估算地表反照率，计算公式：

Albedo=0.356B+0.130R+0.373NIR+0.085SWIR1+0.072SWIR2-0.0018

(1)

式中，B为蓝光波段的反射值；R为红光波段的反射值；SWIR为短波红外反射值，利用Arcgis栅格计算器计算出各年份的地表反照率。

归一化植被指数(NDVI)是一种使用较为广泛的植被指数，是现阶段充分体现植被具体生长情况的一种指数，依据该指数可知植被的覆盖量、植物量以及生长状况等数据[8]。可通过遥感影像近红外波段和红光波段的反射值计算，计算公式：

NDVI=(NIR-R)/(NIR+R)

(2)

式中，NIR为近红外波段的反射值；R为红光波段的反射值。

为了后续数据处理的方便，对二者进行归一化处理：

N=[(NDVI-NDVImin)/(NDVImax-NDVImin)]×100%

(3)

A=[(Albedo-Albedomin)/(Albedomax-Albedomin)]×100%

(4)

式中，N为归一化后的NDVI；A为归一化后的Albedo；NDVImin、Albedomin分别为NDVI和Albedo的最小值；NDVImax、Albedomax分别为NDVI和Albedo的最大值。

对归一化后的Albedo和NDVI数据，利用Origin软件进行线性回归分析并绘制特征空间散点图，可以得出，在Albedo-NDVI特征空间中，Albedo和NDVI有着明显的负线性关系，计算得到拟合方程Albedo=a×NDVI+b，a为方程斜率，b为方程截距，即随着土地沙化程度的加剧，NDVI的值不断变小，而Albedo的值逐渐变大。高中低等级的植被覆盖区域可以明显的划分，反映了沙化土地的地表信息和水热组合状况。其线性拟合方程：

Albedo2000=-0.6243×NDVI+0.8697

(5)

Albedo2005=-0.6490×NDVI+0.8552

(6)

Albedo2010=-0.7785×NDVI+0.9322

(7)

Albedo2015=-0.6413×NDVI+0.8742

(8)

Albedo2021=-0.4038×NDVI+0.6445

(9)

图2 Albedo-NDVI特征空间散点图

1.3.2 沙漠化遥感监测差值指数计算

在以上构建的Albedo-NDVI特征空间中不难发现，两者的关系为线性负相关，若将在代表土地沙化变化趋势的垂直方向划分，即在线性拟合的趋势线上划分二维空间，由此反演出沙漠化差值指数模型(DDI)[9]，通过此模型可以区分不同等级的沙化土地，进而可以对不同程度的沙化土地进行监测和分析，计算公式：

DDI=(-1/a)×NDVI-Albedo

(10)

式中，DDI为沙漠化差值指数；a为特征空间线性拟合的斜率。计算结果表达式：

DDI2000=1.6015×NDVI-Albedo

(11)

DDI2005=1.5408×NDVI-Albedo

(12)

DDI2010=1.2844×NDVI-Albedo

(13)

DDI2015=1.5593×NDVI-Albedo

(14)

DDI2021=2.4765×NDVI-Albedo

(15)

结合《第五次全国土地荒漠化和沙化监测技术规定》和Jenk自然间断点法[10]，将DDI分为极重度沙化、重度沙化、中度沙化、轻度沙化和非沙化。自然间断点法具有分级的良好界限，可以把同一类别中的数据值差异值差异最小化，把类别之间的数据值差异最大化。

1.3.3 Theil-Sen Median趋势分析和Mann-Kendall检验分析法

Theil-Sen median斜率估计被广泛运用在分析长时间序列数据的趋势分析中，是一种使用较普遍计算非参数统计的方法。作为一种检验非参数统计方法，Mann-Kendall优点在于对趋势是否线性不做要求，也无需测量值服从正态分布[11]。Sen趋势分析和M-K显著性检验已成为长时间序列数据趋势判断的重要方法，与线性回归分析相比，具有很强的抵抗数据误差的能力，不易受异常值干扰，在一定程度上提高了准确性。对于土地沙化的变化趋势类型，根据实际情况并参考郭晓娜研究方法[12]划分为严重退化、轻微退化、稳定不变、轻微改善、明显改善。

1.3.4 土地沙化转移矩阵

土地沙化转移矩阵可知其在不同阶段，存在不同程度沙化类型相互转化的特征，反映不同程度沙化之间的转移速率和土地沙化的时空格局变化，其表达式：

(16)

式中，i与j分别为研究始末的土地沙化类型；Pij为土地沙化类型i转为土地沙化类型j的数量。

1.3.5 地理探测器

地理探测器是王劲峰等[13]通过研究得出一种全新的统计学方法，用以诠释背后驱动因子，其核心理念是探讨若某个自变量影响某个因变量，在空间分布上面，自变量与因变量存在相似性特征，可基于空间异质性原理分析各驱动因素对土地沙化的影响。根据计算得到的q值判断各影响因子的影响力大小，依据q值可知，自变量影响因变量的程度，两者呈现正相关关系；p值则代表统计学中的置信度检验标准，p值<0.05则表示为差异显著。地理探测器由因子、交互、风险以及生态等4种探测手段组成，本文使用因子探测进行分析，其计算公式：

(17)

2 结果与分析

2.1 科尔沁左翼中旗土地沙化空间分布特征

利用Arcgis 10.8软件，通过沙化差值指数法计算得到研究区2000年、2005年、2010年、2015年和2021年5期土地沙化遥感影像，空间分布图如图3所示。

图3 科尔沁左翼中旗土地沙化空间分布图

由图4可知，2000—2021年土地分布情况出现较大变化，有一部分土地沙化程度好转，但也有一部分土地沙化程度出现恶化。在西部和东北部出现极重度和重度土地沙化；西北部和东部存在中度沙化土地，其分布较为分散；中部和西南部主要分布的是轻度沙化土地；西南部则主要分布的是非沙化土地。总体而言，研究区土地沙化程度类型主要以中度沙化和轻度沙化土地分布为主。

图4 2000—2021年科尔沁左翼中旗土地沙化等级转移矩阵

2.2 科左中旗土地沙化面积变化特征

由表1可以看出，在研究区中，不同等级沙化土地面积均出现不同变化。2000—2021年，研究区极重度沙化土地面积出现先减后增加的情况，相较于2000年，2021年极重度沙化土地面积减少了0.02%；重度沙化土地面积整体上呈增长趋势，在2015年下降至19.72%，随后在2021年又上涨到21.38%；中度沙化土地面积整体呈减少的趋势，2021年占比26.15%；轻度沙化土地面积整体呈下降趋势，在2015年增加到27.94%，随后在2021年又下降到20.70%；未沙化土地面积在2021年增加到23.53%。

表1 2000—2021年科尔沁左翼中旗沙化土地面积及比例

2.3 土地沙化程度转移分析

由图5可知，2000—2005年，中度转重度和轻度转中度较为显著，转化面积分别为901.71km2和813.19km2，同时也有699.88km2的中度转为轻度；2005—2010年，轻度转中度和重度转中度较为明显，转化面积分别为804.82km2和761.56km2；2010—2015年，面积转化较大的是中度转轻度和重度转中度，转化面积分别为888.78km2和751.10km2，说明此阶段有较多面积的沙化土地好转；2015—2021年主要是轻度和中度之间的转化；2000—2021年，整体上也是轻度和中度之间的转化，有704.61km2的中度转成了轻度，说明总体上土地沙化有好转趋势。

图5 2000—2021年科尔沁左翼中旗土地沙化趋势空间分布图

2.4 土地沙化变化趋势分析

通过Theil-Sen Median趋势分析和Mann-Kendall检验结果对2000—2021年研究区土地沙化22年的变化趋势进行分析，从表2可以看出，变化趋势以稳定不变为主，占研究区总面积的46.71%。呈改善趋势的土地占比为30.75%，主要分布于研究区中部；但也有22.52%的土地呈退化趋势，主要分布于研究区西部和东北部。其中，明显改善土地和严重退化土地均占0.01%。总的来看，沙化状况趋于稳定发展，说明应在加大治理沙化土地的同时，坚持边治理边巩固治理成果的原则。

表2 2000—2021年科尔沁左翼中旗土地沙化趋势分析

表3 科尔沁左翼中旗土地沙化驱动因素解释力统计

2.5 土地沙化驱动因素分析

利用地理探测器模型对选取的驱动因子进行离散化处理和分类发现，研究区土地沙化受到多种驱动因素的影响。由图3可知，各个驱动因子显著性影响土地沙化的程度，且得出其影响程度的顺序为土壤类型>DEM>平均气温>平均风速>日照时数>年降水量>GDP>人口空间分布>土壤侵蚀>坡度。其中，土壤类型、DEM和年平均气温3种因子影响程度均超过10%，是主要影响因子。社会经济因子中GDP和人口空间分布q值排序靠后，这说明自然因素比社会经济因素对土地沙化影响程度更大。

3 讨论

本文基于GEE大数据平台，通过Landsat遥感影像提取了研究区NDVI和Albedo数据，从土地沙化时空演变、转移矩阵、地理探测器模型等方面对科尔沁左翼中旗土地沙化进行监测与驱动因素分析。综合分析，2000—2021年研究区土地沙化主要分布在中部和东北部，整体呈稳定趋势发展，改善土地面积大于退化土地面积，这与众多研究科尔沁沙地沙化的研究结果较相似[14]。土壤类型、DEM和年平均气温3种因子是土地沙化的主要影响因子。科尔沁左翼中旗处于我国北方农牧交错带，生态环境脆弱。从研究方法来看，选用归一化植被指数和地表反照率2种地表参量相结合的方法，综合了对植被敏感的可见光和近红外波段的反射光谱信息，二者拟合的线性函数具有明显的负相关关系，此研究结果与曾永年等[5]研究结果一致，曾永年等指出，随着植被指数的增加，地表反照率减小，土地沙化过程得到了明显的反映。李嘉豪等[15]利用该指数模型研究得到小滦河流域土地沙化动态变化，本研究也认为，沙漠化差值指数可作为土地沙化的监测指标。地理探测器模型现已成为主流地理分析模型，利用此模型可定量化地揭示不同驱动因子对土地沙化的影响大小。本文选取了地理探测器模型中的因子探测模块，后续应结合其他3种模块做进一步分析。土地沙化是我国当前最为严重的生态环境问题之一，近年来，随着退耕还林还草、实施“山水林田湖草沙一体化保护和修复工程”等一系列治理措施，科尔沁沙地治理成效显著，但仍有小部分区域改善不明显。因此，应加大治沙防沙治理力度，因地制宜，阻止土地沙化进一步发展，优化生态环境。

4 结论

2000—2021年，研究区沙化土地主要是轻中度沙化且集中在东北部和中部。对比2015年和2021年数据可知，虽然增加了2.80%的重度和极重度的沙化土地，但同时增加了8.52%的非沙化土地，说明沙化土地有好转的迹象；2000—2021年，研究区沙化土地变化趋势以稳定不变为主，面积占比46.71%，轻微改善面积占比30.75%，轻微退化面积占比22.52%；通过地理探测器模型探测，在自然因素和社会经济因素中，研究区土地沙化主要受到自然因素影响，其中土壤类型是最主要的驱动因子，q值为0.1279，其次是DEM和年平均气温，q值分别为0.1244、0.1055。