喀斯特石漠化是只在亚热带脆弱的喀斯特背景下的一种土地退化过程^[1-2], 表现为自然因素和人类活动综合作用下的地表植被破坏, 水土物质流失^[3], 土地生产力下降, 岩石大面积裸露类似荒漠化景观^[4]。石漠化带来的危害性严重^[5], 已成为我国西南地区生态可持续发展的屏障之一^[6-7]。因此, 精确找准石漠化空间分布, 合理划分石漠化等级是治理石漠化的关键^[8]。

为了能更好为石漠化治理服务提供理论基础, 众多学者对石漠化信息提取进行了研究, 马士彬等^[9]基于Logistic-CA-Markov模型对六枝特区进行了石漠化空间变化规律研究；陈飞等^[10]对典型喀斯特槽谷区石漠化时空演变及未来情景模拟；陈希等^[11]基于多源遥感数据的石漠化信息提取研究, 认为多源数据支持下的石漠化信息提取更准确；李丽等^[12]利用中巴光学遥感数据提出了一种植被基于植被覆盖度的石漠化信息提取方法。但是, 以上方法更多的是局限于小区域, 若按照同样的方法对大尺度的区域进行石漠化解译, 数据的完整获取存在一定的困难。

所以, 为了探究石漠化信息提取方法是否在大尺度研究区域上也能得出较好的结果, 本文选择西南喀斯特地区为研究区, 利用归一化差分植被指数(NDVI)，植被净初级生产力(NPP), 坡度和地表反照率反演2000—2015年四期西南石漠化图, 分析其时空变化规律, 依据不同土地利用、坡度和生态保护区来对石漠化分布进行分析。以期在上述分析中能得出石漠化分布规律, 为研究区生态环境修复提供一定的依据。

1 材料与方法 1.1 研究区概况

西南喀斯特地区包括贵州、云南、广西、湖南、湖北、重庆、四川、广东八省, 位于102—111°E, 23—32°N, 北起秦岭山脉南麓, 南至广西盆地, 西至横断山脉, 东抵罗霄山脉西侧^[13], 面积约54万km^2[14], 然而其地表喀斯特出露面积就达51万km^2[15], 有24个生态功能保护区(图 1)。石漠化集中分布于黔、桂、云三个地区, 以黔地区石漠化面积高达32427 km²最为严重^[16]。年降水量在1000 mm以上, 属于亚热带季风湿润气候^[17]。该区域生态位置很重要, 其不仅是珠江源头以及长江水源的而重要补给区, 还是南北水调水源区和三峡库区^[18]。西南喀斯特环境的脆弱性和易伤性, 以及不合理的人为活动, 导致生态环境遭到严重破坏, 引起一系列生态问题, 例如, 水土流失、森林生态系统退化、土地石漠化等^[19]。

1.2 数据来源及预处理

1:50000岩性图来自喀斯特科学数据中心(http://www.karstdata.cn/), 目的是为了裁剪出非喀斯特区域, 减少石漠化发生区域的误判。

数值高程模型(DEM)和地表反照率数据获取自地理空间数据云(http://www.karstdata.cn/)。DEM空间分辨率为90 m, 为了与其他数据拥有相同的分辨率, 对其重采样为1000 m, 进行坡度提取, 并把坡度分为0°— 2°、2°— 6°、6°— 15°、15°— 25°、>25, 共5个等级。地表反照率空间分辨率为1 km×1 km。

MODIS-NDVI数据源于美国航空航天局(https://www.nasa.gov/)的地球观测系统数据集中的MOD13Q1 L3级数据产品, 其空间分辨率为250 m × 250 m, 时间分辨率为16d。应用MODIS Reprojection Tools(MRT)软件对原始数据的hdf格式转换为tif格式。

土地利用数据来自美国国家航空航天局(https://www.nasa.gov/), 空间分辨率为300 m × 300 m, 为了减少数据量, 在进行石漠化图生成前就已经利用湿地, 水域和建设用地裁剪掉不发生石漠化的区域, 所以, 在分析土地利用中石漠化变化时就没有列出湿地、水域和建设用地三种地类。

NPP数据来源于美国NASA EOS/MODIS的2000—2015年产品(http://files.ntsg.umt.edu/data/NTSG_Products/MOD17/), 空间分辨率为1 km×1 km, 以年为时间分辨率, 单位为kgC/m², 对下载好的NPP乘以0.0001得到实际值。

西南生态功能保护区数据来自中国科学院资源环境科学数据中心(http://www.resdc.cn/)。

以上数据均转换为相同的地理坐标(WGS_1984)和投影坐标(Albers Equal-Area Conic), 并且空间分辨率统一采样为1 km×1 km。

1.3 方法 1.3.1 石漠化解译路线

本文石漠化解译方法基于地理探测器, 具体流程为：(1)准备好反演石漠化需要的基础数据, 包括NDVI、NPP、地表反照率、数值高程模型(DEM)和土地利用/土地覆盖变化(LUCC), 利用岩性图裁剪出喀斯特。由于认为LUCC中的建设用地、湿地和水域不发生石漠化, 所以, 用建设用地、湿地和水域进一步裁剪, 这样既减少石漠化发生在非喀斯特和不发生石漠化地区的误判, 也减少数据量。(2)在ArcGIS平台根据已有的贵州省石漠化图进行样本点信息的提取, 利用已有的石漠化样本点信息对裁剪好的NDVI、NPP、坡度和地表反照率一一进行样本提取。对提取好的样本信息导出ArcGIS另存为表格。打开地理探测器进行权重分析, 自变量选择石漠化信息, 因变量选择四个变量, 在0.01置信度内得出NDVI、NPP、Slope和Surface Albedo的Q分别为46.2、36.3、11和16.5, 根据Q值计算出相对贡献率, 分别为0.42、0.33、0.1和0.15(表 1)。(3)依据得出的因子权重反演出西南喀斯特地区2000—2015年间的四期石漠化空间图。基于王世杰^[20]等石漠化分级体系, 将石漠化分为无石漠化、潜在石漠化、轻度石漠化、中度石漠化、重度石漠化和极重度石漠化, 分别对应序号1、2、3、4、5和6。

1.3.2 因子贡献率计算

地理探测器是探测和利用空间分异性的工具, 由分异及因子、交互作用、风险区和生态探测四个部分组成^[21]。其分异及因子探测某因子X多大程度上解释了属性Y的空间分异。因此可以利用公式来计算反演石漠化的因子对其解释能力大小^[22]。表达式如下：

(1)

式中a为自变量, b为石漠化图, N为因变量的样本点总数, N_i为自变量的样本点总数, δ²、δ_i²分别为a和b的栅格值方差, i和h分别为第i类因子等级及其等级总数。

1.3.3 石漠化转移矩阵

石漠化转移矩阵能够清晰明了地反映研究区域石漠化的的具体流向与各种石漠化等级的转入面积与转出面积^[10], 进而分析石漠化变化的数量特征。具体公式如式；

(2)

式中：P_ij为石漠化转移矩阵；n为石漠化的类型数；i、j分别为研究期初与研究期末的石漠化等级序号。

1.3.4 石漠化变化量和变化率

用变化量和变化率两个指标对石漠化变化进行定量分析^[23]。公式如下：

(3)

(4)

式中：S_a和S_b分别表示前期和后期土地石漠化的面积, ΔS和K分别表示2000—2015年间石漠化的变化量和变化率。如果变化量和变化率为负表示增加, 为正则表示减少。

2 结果与分析 2.1 石漠化时空分布 2.1.1 石漠化随时间的变化

由表 2可知, 西南石漠化以轻度和中度为主, 其中两种等级石漠化面积之和占西南八省总面积的比例在4个时期分别为76.76%、76.34%和79.83%、84.98%, 由比例可知, 轻度和中度石漠化面积在2000—2005年间减少, 2005—2015年间增加, 总体呈现上升趋势。按照比例来看, 轻度先减少后增加, 总体呈上升趋势, 轻度比例依次为41.80%、37.01%、42.58%、49.57%, 中度呈现增加的趋势, 比例依次为34.96%、39.33%、37.25%、35.41%。

由表 3可知, 2005、2010、2015年间石漠化土地分别为13.22万km²、12.03万km²和10.70万km², 国家公报2005、2011、2016年间石漠化土地分别为12.94万km²、11.98万km²、10.06万km²。利用本研究减去公报石漠化土地面积得到石漠化的变化, 分别为0.28万km²、0.05万km²、0.64万km², 用变化面积除以本研究的面积得到变化率分别为2.09%、0.38%、5.98%。

3.1.2 石漠化空间变化

由表 4可知, 2000—2015年间贵州石漠化占比最大, 分别为31.46%、27.00%、25.10%和27.66%；云南次之28.01%、21.79%、23.61%和26.36%；广西分别为22.96%、18.00%、16.04%和15.70%。从占比来看, 整体呈现增加趋势的有湖南、湖北、重庆和四川, 其中以2010年为节点, 表现出先增加后减少；广州和其他地区石漠化占比呈现减少的现象, 但是云南和广东以2005年转点, 表现出先增加后减少, 贵州也表现出减少的趋势, 但是在2010以后也有增加。

以贵州为中心向四周辐射的西南八省石漠化空间分布(图 2)显示, 石漠化发生地区最明显的为贵州, 云南和广西三个省份。石漠化发生可以分为两个时间段, 即2000—2005年间是石漠化的增长阶段, 2005—2015年为石漠化减少阶段。研究时段内, 贵州, 云南, 广西石漠化逐渐呈现减少趋势, 其他区域呈现先增加后减少的趋势。

3.2 不同土地利用中石漠化变化

不同的土地利用类型中人类活动强度有很大的差异, 由表 5可知, 2000—2015年间无石漠化面积最多, 潜在石漠化次之。石漠化主要发生在耕地和林地两种土地利用类型上, 其他用地上石漠化发生面积最少, 耕地和林地上无石漠化呈现增加趋势, 草地和其他用地上的潜在石漠化情况则相反。

2000—2015年间已发生石漠化面积变化, 耕地上轻度石漠化先从2000年的31989.49 km²减少到2005年30653.93 km², 然后又增加到2015年31841.68 km²。中度石漠化从2000年25636.75 km²增加到2005年26935.07 km², 然后又减少到2015年22179.91 km²。重度石漠化从2000年1999.33 km²增加到2005年13663.93 km², 然后又减少到2015年7462.70 km²。极重度石漠化呈现减少趋势。由表 5可知, 林地和草地大致与耕地上已发生面积变化一样, 呈现出轻度先减少后增加, 中度和重度先增加后减少, 极重度呈现减少趋势。其他用地中轻度、中度、重度和极重度石漠化呈增加趋势。发生石漠化的影响因素众多, 但是人类活动的强度也促进石漠化的加速形成, 所以, 耕地中石漠化面积分布较多。

由图 3可知, 土地利用类型不同, 石漠化的分布具有很大的差异。从各石漠化等级面积比例来看, 不同的地类中比例最大的是无石漠化, 潜在石漠化次之。已发生从土地利用类型上的石漠化分布比例来看, 耕地上的石漠化发生比例最大, 草地次之, 然后是林地, 其他用地最少。因为耕地是人类活动较强的地类, 且很多耕地分布在一定坡度的地表上, 坡度越大, 地表覆盖的土层就越薄, 冲刷越强, 水土流失快, 所以石漠化发生面积较大。其他用地中的极重度石漠化比大, 其他用地一般为岩石裸露率大的难利用土地, 地表植被覆盖少, 土层薄, 极重度石漠化发生概率大。

3.3 不同坡度等级内石漠化分布

由表 6可以看出, 不同坡度石漠化类型在不同时间段内的变化差异很大, 2000—2015年间无石漠化在6°—15°坡度上发生面积最大, 分别有53580.7 km²、57962.97 km²、60238.36 km²和59526.82 km²和35146.43 km²。在>25°坡度上最少, 分别为4231.59 km²、4311.38 km²、4404.40 km²和4373.18 km²。在0°—6°坡度上的土地由于接近平地, 所以, 无石漠化面积也很大, 潜在石漠化分布与无石漠化分布具有相似的规律。已发生石漠化的土地在坡度上的面积分布从大到小依次为6°—15°、2°—6、15°—25°、< 2°、>25°, 轻度石漠化、中度石漠化、重度石漠化和极重度发生规律相似, 都在6°—15°坡度上达到最大值。但是2000—2015年间不同等级已发生石漠化在6°—15°上具有不同的演变规律, 其中轻度石漠化以2005年为节点先减少后增加, 由2000年16954.40 km²减少到2005年13854.05 km², 减少了3100.35 km², 后又增加到2015年16056.45 km², 增加了12956.10 km²。中度石漠化逐渐减少, 由2000年15776.06 km²减少到2015年11340.23 km², 减少了4435.83 km²。重点石漠化在2000—2005年间增加了2437.51 km², 2005至2015年间减少了5030.03 km²。极重度也表现出逐渐减少的现象, 2000—2015年减少了649.36 km²。

由图 4不同坡度下石漠化占比可知, 不同坡度内以无石漠化和潜在石漠化为主。已发生石漠化土地上, 轻度石漠化主要分布在 < 25°坡度内, 极重度石漠化随着坡度的变大逐渐增加。2000年中度石漠化在>2°坡度土地上都有分布, 重度石漠化在不同坡度都有发生。2005年中度石漠化随着坡度的变陡逐渐增加, 在>25坡度上, 中度石漠化大于轻度石漠化, 重度石漠化主要分布在坡度 < 25°的土地上, 在>25°坡度上的比例最小, 坡度6—15范围内达到最大。2010—2015年间中度石漠化随着坡度变大呈现增加的趋势, 2010年中度石漠化在随着坡度变大而增加, 重度石漠化在坡度>25°的土地上占比最少。2015年中度石漠化随坡度变大逐渐增加, 其中>25°坡度的中度石漠化占比最大, 重度石漠化主要分布在坡度>2°土地上。综上分析, 石漠化主要发生在较为平坦的地方, 喀斯特区域的山间坝子较少, 平坦的地方土地利用率较高, 人类活动更加频繁, 随着坡度的增加, 地表起伏变大, 人类活动较少, 石漠化发生的概率变小。

2.4 西南生态保护区石漠化变化

由表 6可知, 2000年和2015年生态保护区石漠化面积最大的是西南喀斯特区, 面积分别为27481.86 km²和21738.65 km²。从石漠化变化量来看, 呈现增加趋势的生态保护区有云南东川、三峡库区、岷山-邛崃山、南水北调中线工程区、大别山地区、武陵山地区。其中三峡库区增加最多, 南水北调中线工程区次之, 分别增加1641.22 km²、1466.29 km², 增加最少的是大别山, 为0.51 km²。石漠化减少最多的的生态保护区是西南喀斯特区5743.22 km², 南岭山地次之, 减少393.34 km², 减少面积最少的是秦岭山地生态保护区, 减少2.11 km²。从石漠化增加的变化率看, 变化率最大的是岷山-邛崃山生态功能保护区, 变化率为770.47%, 大别山变化率为439.15%, 南水北调中线工程区变化率为121.45%, 变化率最少的是武陵山山地区, 变化率为6.81%。从石漠化减少的变化率来看, 减少最多是横断山南部生态保护区, 变化率为82.58%, 桂西南石灰岩区变化率次之, 为56.62%, 变化率最少的是秦岭山地生态保护区, 为8.09%。综上所述可知, 石漠化减少最多或是增加最多的生态功能保护区不一定就是变化率最大, 这说明石漠化生态修复呈现出边治理边恶化的情况, 此外, 石漠化生态修复不能只注意已经发生石漠化的区域, 还应注意未发生石漠化区, 要做到石漠化恶化地区得到改善, 未发生石漠化地区得到预防。

由图 5可知, 不同的保护功能区内石漠化转移变化明显的是水土保持生态功能保护区。综上分析, 2000—2015年间西南生态保护功能区是发生石漠化减少最大的地方, 这是由于研究期间国家和政府开展一系列的生态工程, 如2003年退耕还草工程、2005年集体林权制度改革、2008年石漠化综合治理、2014年全国生态保护与建设规划(2013—2020年)、2015年生态文明体制改革总体方案的等, 这些生态修复工程对西南喀斯特石漠化的治理起着很大的作用。三峡库区生态功能保护区是石漠化增加最多的地方, 石漠化公报显示2016年山峡库区石漠化土地面积为5630km², 比2011年少了1040km², 这与本研究的2015年三峡库区石漠化3436.79 km²相比, 多2193.21 km²。

3 讨论 3.1 土地利用对石漠化空间分布的影响

西南喀斯特地貌分为中高山、断陷盆地、岩溶高原、岩溶峡谷、峰丛洼地、岩溶槽谷、峰林平原、溶丘洼地以及局部分布的石林等^[24-25]。奇特的地形地貌决定了喀斯特地区平原坝子会被城市、村镇居民点以及大型工厂等利用, 耕种类型以山地旱地为主。土地利用类型以林地、耕地和草地为主^[26]。本研究得出林地与耕地中石漠化发生面积最多, 这与陈欢等研究万山石漠化与土地利用后, 指出发生石漠化的土地类型主要为旱地和林地(有林地、灌木林地)相符^[27]；李阳兵等研究喀斯特地区土地利用与石漠化关系后, 指出石漠化土地中以灌木林地、低覆盖度草地、未利用地、山区旱地、石质旱地等土地利用类型占较大比例结果相近^[28]。耕地对石漠化的影响除了与人口密度有关外, 还与农村聚居的密集程度、单位面积上的壮劳力有关^[29]。但是, 其根本的原因是土地承载力低下, 人口压力超过了土地承载力, 洼地耕地资源不足, 迫使当地耕种坡耕地^[30]。因为人类在坡耕地的耕作活动不仅破坏了山地旱地的土壤肥力、抗侵蚀性、颗粒组成、水稳定团聚体、地表植被等^[25], 还促进了土壤地表、地下流失和生物流失, 间接驱动坡地石质化, 犁耕运移土壤和农耕促进土壤流失, 直接和间接的驱动坡地石质化^[31]。林地对石漠化的影响主要在于土地利用类型的转化。当植被从有林地向疏林地、灌木林地和草地转化时, 植被覆盖度慢慢降低, 景观类型变得单一, 生态系统稳定性差^[32], 土壤侵蚀发生, 引发石漠化景观产生。

3.2 坡度对石漠化空间分布的影响

石漠化发生较严重的地区主要分布在贵州西南部及云南东部一带, 因为该地区石灰岩比重较大、地表土层较薄、成土速率慢、降水丰沛、山地面积较广、耕地比重大。地表坡度大小决定了土壤侵蚀强度, 但是随着坡度的越来越大, 土壤侵蚀强度增加会使得地表土壤变薄, 基岩裸露, 石漠化发生。本研究得出潜在石漠化主要发生在 < 2°和>25°坡度上, 已发生石漠化主要发生在6°—15°之间, 重度石漠化和极重度在>25°坡度上变化幅度较小。王明明等研究后寨河流域喀斯特石漠化演变特征及其关键表征因子与驱动因素后, 指出石漠化区域主要集中在2°— 22°坡度上^[33]。李阳兵等研究王家寨小流域石漠化与坡度的关系后, 指出石漠化集中分布在15°— 25°坡度范围内, 认为簸箩地区轻度、中度、强度集中分布在 < 5°、15°—25°的坡度上^[34]。除此之外, 李阳兵等在研究了后寨河流域石漠化演变特征后, 指出15°—25°、>25°坡度范围内的强度石漠化比例基本保持不变^[35]。周忠发等研究喀斯特地区石漠化与地形坡度关系后, 指出地形坡度陡, 人类活动频繁的地区石漠化面积大, 石漠化程度高^[36]。国家农业坡度上限为25°, 也就是说>25°以上的耕地实施退耕还林还草政策。但是, 综合上述分析, 石漠化面积发生最大的坡度范围是6°—15°, 在15°—25°范围内石漠化发生也严重。因此, 石漠化严重区域可以适当的把退耕还林还草范围降得再低一点, 将喀斯特地区降低退耕还林坡耕地坡度标准为15°^[37]。

3.3 西南石漠化动态变化的影响因素分析

石漠化的发生除了存在复杂的自然因素, 还有人类活动的影响。研究时段内石漠化演变存在2000—2005年的增加、2005—2015年减少的两个阶段, 但是石漠化整体情况呈现改善趋势。由图 6可知, 石漠化加剧阶段, 研究区农村人口众多, 对土地依附加强、土地需求较大, 不合理的开垦方式, 对森林的乱砍乱伐、林地减少。其次, 洪水自然灾害的爆发, 引起土壤大面积侵蚀发生, 导致土地破坏、加剧洪涝灾害淤塞湖泊、影响开发利用, 最大的后果就是引起石漠化的发生, 导致土地石漠化加剧。石漠化的减少阶段, 由图 6可知, 1999年底“退耕还林还草”和西部地区“山川秀美”工程等生态修复政策得到实施^[38], 2003年退耕还草工程、2005年集体林权制度改革和“绿水青山就是金山银山”的科学概论被提出、2008年石漠化综合治理、2014年全国生态保护与建设规划(2013—2020年)、2015年生态文明体制改革总体方案、2015年国务院关于全国水土保持规划(2015—2030年)的批复等生态修复政策的实施。随着城镇化的发展, 大量农村人口外出务工、务农人员减少使得耕地逐渐撂荒；沼气、太阳能和液化气等普及, 改变了传统的能源获取方式, 提高了人们能对现有能源的利用率, 降低了人们对木柴资源的依附, 减少对林地的砍伐和毁坏, 林地得到增加, 石漠化逐渐好转。

3.4 石漠化解译方法的优缺点

由表 5可知, 贵州省石漠化反演结果中, 轻度石漠化为22019.65 km², 中度石漠化为3730.00 km², 重度石漠化为3001.07 km², 极重度石漠化700.57 km²。与陈起伟^[39]等人的结果对比, 分别相差3131.35 km², 1403.60 km², 728.93 km², 87.43 km²。2000年贵州石漠化反演结果中轻度石漠化为21924.83 km²、中度石漠化12061.78 km²、重度石漠化2525.29 km²、极重度石漠化592.82 km², 依据李建存^[40]等人文献石漠化显示, 其中重度包含了本文的重度和极重度石漠化, 所以, 石漠化对比分别为轻度、中度和重度, 分别相差1075.17 km²、1061.78 km²、118.10 km²。综合以上分析, 本文石漠化结果误差在可控范围内, 因此基于地理探测器得出反演石漠化因子权重进行石漠化解译具有一定的优点, 使用方法简单, 解译快速, 结果覆盖范围广, 数据易获取。缺陷在于石漠化分布结果图存在不确定性, 因为采用的数据空间分辨率较低, 喀斯特地区地表复杂, 精度较低的石漠化图不能在小区域上得到精确地体现。在今后的研究中, 我们将进一步加强石漠化解译的精度, 争取反演解译出更精准的石漠化图, 为生态治理提供一定的支撑。

5 结论

(1) 西南石漠化等级分布以轻度和中度为主。其发生比例变化先减少后增加, 总体呈上升趋势。从空间分布来看, 石漠化发生分布面积最广的是贵州, 其次为广西和云南。研究时段内的石漠化发生具有两个时间段, 即2000—2005年的增加时段, 2005—2015年的减少时段。研究时段内, 西南石漠化面积整体呈现改善趋势。

(2) 从不同土地利用来看, 2000—2015年间无石漠化面积最多, 潜在石漠化次之。石漠化主要发生在耕地和林地两种土地类型上, 耕地因为人类活动很强, 石漠化的发生比例最大。其他用地上石漠化发生面积最少, 但是极重度石漠化在其他用地上的发生比例很大, 平均在11%左右。

(3) 从不同坡度来看, 石漠化的发生并不是随坡度的增加而呈现在正相关。各类型石漠化在坡度等级分布中的变化如下：< 15°坡度范围内, 轻度石漠化以2005年最少, 中度和极重度以2015年石漠化分布面积最少, 重度石漠化以2010年最少；15°—25°坡度范围内的强度石漠化在2000—2015年也有发生, 说明此坡度范围也受到垦殖等强烈人为干扰。>25°各等级石漠化面积变化不大。总之, 在坡度6°—25°之间石漠化发生面积最大。

(4) 从生态保护区来看, 2000和2015年西南喀斯特生态保护区是石漠化面积分布最多的区域, 分别为27481.86 km²和21738.65 km², 最少的是大别山山地生态功能保护区。从变化量来看, 增加最多的是三峡库区, 增加1641.22 km², 减少最多的是西南喀斯特生态功能保护区, 减少5743.22 km²。

(5) 本文基于地理探测器得出的石漠化反演因子贡献率, NPP为0.33、NDVI为0.15、地表反照率为0.42、坡度为0.1。