中国南方喀斯特地区存在着非常严重的石漠化问题^[1-3], 在此地区的自然灾害和贫穷发生率极高, 严重制约着该地区的发展。目前, 研究其空间演变规律已成为治理石漠化的重要前提, 以往学者利用不同方法分析了石漠化的空间演变规律：熊康宁等基于2000年Landsat TM影像和实际调查数据, 将贵州省石漠化分为6级并探讨了贵州省石漠化的现状及其空间演变趋势^[4]；胡宝清等基于2000年TM影像提取了广西省石漠化信息, 并利用数理统计等方法探讨了石漠化的分布演变规律, 得出了人类活动是石漠化产生的主要原因^[5]；白晓永等基于1986、1995、2000年3期贵州省石漠化数据, 利用转移矩阵等方法探讨了贵州省石漠化的时空演变规律^[6]。可以看出, 以往的研究多基于实际数据和数理统计的方法综合分析, 所取得的成果对早期的石漠化治理提供了巨大的科技支撑, 然而研究方向多集中于省域级大尺度低空间分辨率上, 且研究年限多在20世纪80年代—21世纪初, 所取得的研究成果和经验方法对已治理了20余年的石漠化指导意义不足。

此外, 如何快速高效地获取石漠化各因子的贡献率一直是众多学者的研究重点, 在以往研究中学者们通过数学分析的手段总结出了很多经典的方法。蓝安军等基于主成分分析、相关性分析, 并建立了有关石漠化的动力指数公式, 对贵州省的石漠化众多影响因子分析得出了包括植被覆盖度、喀斯特面积、河网密度、耕地、未利用地、草地、海拔、土地开垦率和25°以上的土地面积等9个因子是导致石漠化产生的驱动因素^[7]；胡宝清等基于主成分分析, 主要分析了岩性、地貌、水热条件、土壤和NDVI等5个指标, 计算得到石漠化地质背景区的环境生态指数回归模型^[8]；皇甫江云等基于专家打分的方法, 对各等级石漠化对应的影响因子归类打分, 得出草地石漠化中各因子的贡献率^[9]。可以看出, 以往的经典方法受人为影响较大, 且由于喀斯特环境的复杂性, 石漠化影响因子回归模型无法推广。

因此, 针对以往研究中存在的两大问题, 本文选取贵州省中部典型喀斯特地区普定县为研究区, 基于Landsat、DEM等30 m分辨率数据, 提取了植被覆盖度、岩石裸露率和坡度等3个最具代表性的指标来精确反演石漠化^[10-12], 旨在解决以下3个问题：(1)寻找更加高效快速的石漠化模拟反演方法。(2)引入地理探测器, 采用更科学的方法确定各因子的贡献率。(3)研究各等级石漠化分布规律, 分析新形势下贵州典型喀斯特地区的石漠化演变速率, 揭示石漠化的时空演变特征(技术路线见图 4)。研究结果中获得的各因子贡献率为反演石漠化提供了参考, 发现的石漠化分布演变规律为石漠化治理提供了方向。

1 研究区与数据源 1.1 研究区概况

位于贵州省中西部地区的普定县(图 1), 全县国土面积约1081.8 km², 喀斯特面积为869.64 km², 年均降水量约1200 mm, 年均温15.3℃, 平均海拔为1100 m, 地势南北高, 中间低。截止2013年, 普定县人口密度达到440人/km², 远远高于我国西部其他地区。由于岩溶地貌的广泛发育, 石灰土比重占到了63%, 连年的高强度降水和人类活动, 使得普定县石漠化极为严重, 2008年, 普定被确定为全国100个、全省55个石漠化综合治理试点县, 众多科研单位在此地区开展石漠化治理工作, 普定县成为了石漠化防治研究的典型地区。

1.2 数据源

本研究使用的实验数据主要包括：普定县岩性图、普定县石漠化分布图、Landsat TM/ETM+、DEM等30 m空间分辨率数据, 具体见表 1。

2 研究方法 2.1 石漠化指标 2.1.1 岩石裸露率

岩石裸露作为石漠化地区最显著的特征, 在石漠化信息的提取中有着至关重要的作用, 本文拟采用NDRI来获取研究区岩石裸露率, NDRI(Normalized Difference Rock Index)是指归一化岩石指数, 基于喀斯特碳酸盐岩的光谱信息在短波红外波段具有较大差异^[13-14], 基于NDVI思想, 在ENVI 5.2软件下进行波段运算, 根据公式(1)得到研究区的NDRI：

(1)

根据像元二分模型原理, 假设像元由岩石裸露和非裸露两部分组成, 根据NDRI, 计算得到基岩裸露率(Fr)：

(2)

式中, NDRIr为全由岩石裸露得到的像元值, NDRIo为全无岩石裸露得到的像元值。在统计得到的所有像元值中, 取频率1%的像元值为NDRIo, 取累计频率为99%的像元值为NDRIr。

2.1.2 FVC(Fractional Vegetation Cover)

FVC植被覆盖度指的是植被冠层在地面的垂直投影面积与所处地总面积之比^[15-16], 同样根据像元二分模型原理, 本文拟采用NDVI来获取研究区植被覆盖度, 公式如下：

(3)

在计算得到像元的NDVI值后, 统计其百分比, 取1%的像元值为NDVIo, 作为没有植被存在的像元值, 取99%的像元值为NDVIr, 作为完全由植被覆盖的像元值。

2.1.3 Slope

本文选择的DEM数据为GDEMDEM 30 m分辨率数字高程数据, 数据标识为ASTGTM_N26E105, 在ArcGIS下对其进行投影转换为ALBERS等面积投影, 坐标系为Krasovsky_1940_Albers, 掩膜提取后制作得到坡度图。

2.2 贡献率获取

地理探测器是探测和利用空间分异性的工具^[17], 在空间一致的情况下, 通过采集样本点上各自变量和因变量的值, 可以利用公式计算自变量对因变量的解释能力大小, 公式表达式为：

(4)

式中, x为自变量, 即NDRI、FVC、slope等因子等级, y为石漠化等级。N、N_i分别为因变量和自变量的样本点总数, δ²、δ_i²为分别为y和x的栅格值方差, i和h分别为第i类因子等级及其等级总数, Q值取值在0—1之间, 值越大说明自变量x对因变量y的解释能力越强。为快速准确地获取贡献率, 采用自然间断法将自变量分为10级, 在ArcGIS 10.2中以500 m为间距均匀布点, 并且布点位置仅在已发生石漠化区域, 石漠化影像裁切于喀斯特科学数据中心2000年贵州省石漠化分布图(图 2), 各因子等级分布图由2001年Landsat遥感影像获得(图 3)。

2.3 石漠化分类定级标准

基于前人工作建立石漠化分级分类标准^{[10, 18-21]}(表 2), 通过地理探测器得到各因子贡献率之后, 在ArcGIS 10.2软件平台下按照下表进行重分类, 并叠置分析, 基于普定县岩性图和土地分类结果剔除非喀斯特地区和水体、城乡建筑等不可能发生石漠化的地方, 得到石漠化空间分布图。

2.4 石漠化空间动态变化分析

石漠化动态变化无论是在时间上还是在空间上都是一个复杂的过程, 传统的数字分析表达能力欠佳, 因此本文在探讨石漠化时空演变时将数学模型与制图综合在一起, 旨在探讨以下科学问题：第一是明确石漠化因子的变化特征, 以探讨因变量对自变量的影响大小；第二是明确石漠化的转移方向；第三是明确各等级石漠化的转移快慢。因此, 有必要对石漠化的因子变化和时空转变分别构建模型讨论。

2.4.1 石漠化因子变化

本文中石漠化的影响因子主要由植被覆盖度、岩石裸露率、坡度等组成, 而坡度由于变化较小无法捕捉, 故不在此讨论。在石漠化演变过程中, 像元中的植被覆盖度和岩石裸露率等级必然也在发生着变化, 因此, 两者的变化必然与石漠化演变具有密切关系, 基于此, 构建下述公式：

(5)

式中, Q指两个不同年份的因子在栅格计算器下相减后的值, q₁指一段时间内该像元变化前的等级, q₂指变化后的等级。

2.4.2 石漠化转移矩阵

在某一时间段内, 石漠化的演变方向与演变规模大多较复杂, 且时常因地而异, 因此需要定量化地讨论以明确其大致的演变方向和趋势, 构建转移变化矩阵。

(6)

式中, S为面积, 单位为km²；m, n为一段时期内前期和末期石漠化的等级；S_{k, k}为在该时期内石漠化的转移面积, 单位为km²。

2.4.3 石漠化转移变化速率

仅仅探讨石漠化的转移方向和面积对于治理石漠化的指导意义是远远不够的, 石漠化的演变快慢反映了一种类型石漠化在该地区的变化发展模式, 在横向比较上更具有科学价值, 也更能说明石漠化的演变情况。因此需要建立模型以探讨石漠化的演变速率。

(7)

式中, P_i为某一类型的石漠化演变速率, 单位为km²/a；S₂为一定时间段内末期的面积, S₁为前期的面积, 单位为km；T为该时间段, 单位为a。

2.5 技术路线 3 结果与分析 3.1 权重分析

在地理探测器中根据石漠化和各因子等级计算得到各自所占的权重, 结果表明, 三个指标中植被覆盖度和岩石裸露率主要控制了石漠化的分布：FVC(0.47)>NDRI(0.45)>Slope(0.15), 由于三者相加大于1, 故各因子权重所占的百分比即为本文中得出的石漠化贡献率：FVC(44%)、NDRI(42%)、Slope(14%)(表 3)。

3.2 石漠化分布及各因子等级图

按照表 2已划分的石漠化等级标准, 在ENVI 5.2、ArcGIS 10.2软件下制图分类, 得到普定县1990—2015年FVC、NDRI、Slope等级图, 并在地理探测器中获得其各自的贡献率, 在ArcGIS下叠置分析, 并结合Landsat影像, 在ENVI下得到土地利用分类图, 剔除了水体、城乡建筑地等不可能发生石漠化的地区, 得到各年份石漠化空间分布图。

3.2.1 NDVI提取植被覆盖度

首先在ENVI 5.2下波段运算得到的1990—2015不同年份的NDVI, 然后再根据公式(3)计算得到植被覆盖度, 并在ArcGIS10.2软件下利用Reclassify工具进行重分类, 按照表 2将其分为＞71%、51%—71%、31%—50%、21%—30%、11%—20%、＜10%共六类并分别赋值为1、2、3、4、5、6, 得到1990—2015年植被覆盖度分布图(图 5)。由图中可以看出, 越靠近河流和城镇地区, 植被覆盖度的等级越低。

3.2.2 NDRI提取岩石裸露率

首先根据公式(1), 在ENVI 5.2下波段运算得到的1990—2015不同年份的NDRI, 再利用公式(2)计算得到岩石裸露率, 并在ArcGIS下重分类, 将其按照表 2分为＜20%、21%—30%、31%—50%、51%—70%、71%—90%、＞90%共六类并分别赋值为1、2、3、4、5、6, 得到1990—2015年岩石裸露率分布图(图 6)。由图中可以看出, 越靠近河流和城镇地区, 岩石裸露率的等级越高, 且1997年的岩石裸露面积明显最小。

3.2.3 DEM提取坡度

在投影变换后, 利用ArcGIS软件中的Slope工具根据DEM数据提取坡度, 并按照表 2的标准, 重分类得到普定县坡度空间分布图(图 7), 由图中可以看出, 坡度较陡的区域多沿着普定县SW—NE方向分布, 总体呈现出南北坡度大, 中部较平坦。

3.2.4 石漠化分布图

在得到各年段的植被覆盖度、岩石裸露率和坡度分布后, 利用ArcGIS 10.2中的叠加分析工具, 按照地理探测器计算得出的各因子贡献率叠加制图, 并剔除非喀斯特地区和水体、城乡建筑等不可能发生石漠化的区域, 得到1990、1997、2006和2015年的石漠化空间分布图(图 8)。

3.2.5 精度验证

为检验反演效果, 本文选取了喀斯特科学数据中心提供的2005年普定县石漠化实际分布图作为标准, 将制作得到的2006年数据进行了对比, 由图 9和表 4可以看出, 石漠化在总体上反演精度较高, 达到了80%, 该石漠化反演方法可信度较高。

3.3 时空演变分析 3.3.1 石漠化面积统计分析

根据石漠化空间分布图, 利用栅格计算器计算得到普定县不同年份的各等级石漠化面积, 统计如下表 5, 可以看出, 无石漠化面积占比最大最稳定, 而极重度石漠化由于面积小, 变化面积不明显, 而潜在石漠化、轻度石漠化和中度以上石漠化的面积在不同年段发生着剧烈地变化。1990—1997年间, 石漠化面积减少了70.95 km²；1997—2006年间, 石漠化面积增加了216.45 km²；2006—2015年间, 石漠化面积又减少了111.05 km²；1990—2015年间, 石漠化面积增加了34.50 km²。综上, 普定县轻度以下石漠化分布面积在喀斯特地区均有大量分布, 而中度以上多分布在三岔河流域沿岸附近；在1990—2015年间, 普定县石漠化演变复杂剧烈, 经历了由好转到恶化再到好转的阶段, 演变总体趋势在好转, 但局部在恶化, 说明治理石漠化不可掉以轻心, 还需要加大治理力度。

3.3.2 石漠化因子时空演变分析

根据图 5以及公式(5), 在ArcGIS 10.2软件下运用栅格计算器进行减法运算, 将各年份植被覆盖度相减, 得到1990—1997、1997—2006、2006—2015、1990—2015年间的植被覆盖度变化图(图 10), 其中, 正值表示该地区的植被覆盖度等级在升高, 负值则表示在下降。可以看出, 1990—1997年植被覆盖度大都在减少, 而在2006—2015年却大都在增加, 在1997—2006年间, 植被覆盖度的变化则具有明显的地理位置特征。

再次运用栅格计算器工具, 将图 6中的各年份岩石裸露率相减, 得到1990—1997、1997—2006、2006—2015、1990—2015年间的变化图(图 11), 同样, 正值表示该地区的岩石裸露率等级在升高, 负值则表示在下降。可以看出, 1997—2006年岩石裸露率等级大都在升高, 而在其他年段则大都在降低, 结合普定县行政区分析, 升高的地区多位于普定县中东部城关镇一带。

3.3.3 石漠化转移变化时空演变分析

根据公式(6), 在栅格计算器下通过栅格叠加得到某一年份对应另一年份的石漠化面积, 统计结果见图 12(NKRD—无石漠化；PKRD—潜在石漠化；LKRD—轻度石漠化；MKRD—中度石漠化；SKRD—重度石漠化；ESKRD—极重度石漠化；Non—karst：非喀斯特区域)和图 13, 可以看出：1990—1997年, 石漠化好转的比例(25.06%)＞恶化(16.76%)；1997—2006年, 石漠化好转(9.73%)＜恶化(37.83%)；2006—2015年, 石漠化好转(28.26%)＞恶化(14.20%)；1990—2015年, 石漠化好转(24.65%)＜恶化(25.47%)；此外, 在图中可以看出由无、潜在石漠化转变为重度以上石漠化的面积虽然不大, 但占比极高, 结合影像分析(图 14), 发现此种演变方式集中于三岔河流域沿岸附近, 由于河道拓宽, 导致沿岸一带植被减少、基岩裸露率增大, 且此地区坡度均较大, 是石漠化发生急剧变化的主要原因。而重度以上石漠化转变为无、潜在石漠化面积较少, 结合影像分析(图 14), 发现此种演变方式多发生于城关镇一带, 此处靠近普定县城, 且海拔较低, 有利于土壤积累, 极少部分较为严重的石漠化地区突变为了无石漠化或潜在石漠化, 笔者认为是居民的开垦荒地导致附近的石漠化土地演变为了耕地或林地, 因此可以看到人类活动的两面性。

综上, 在石漠化演变过程中, 轻度以下石漠化演变最为剧烈, 尤其是潜在石漠化和轻度石漠化, 演变面积大, 方向难测；石漠化的演变方式因地而异, 在靠近河流沿岸地区, 石漠化常发生突变而恶化, 而在城乡居民地附近, 常发生重度石漠化突变为无石漠化。

3.3.4 石漠化转移变化速率

再次运用栅格计算器计算得到普定县石漠化转移速率变化图(图 15), 并统计其转移矩阵如图 16, 可以看出重度以上石漠化演变速率极低, 几乎等于0, 这也验证了前人的工作基础^{[6, 22]}。在1990—1997年间, 轻度-潜在石漠化变化速率最大, 为11.2 km²/a；1997—2006年间, 潜在-轻度石漠化变化速率最大, 为16.5 km²/a, 2006—2015年间, 轻度-潜在石漠化变化速率最大, 为10.4 km²/a, 1990—2015年间, 潜在-轻度石漠化变化速率最大, 为3.22 km²/a。计算各等级石漠化平均变化速率得出(图 17), 潜在石漠化恶化的速率最大, 为2.75 km²/a, 其中, 在1997—2006年间中速率达到了16.5 km²/a, 演变十分迅速。此外, 在各年段上, 中度以上石漠化整体均在好转, 而轻度以下石漠化则反复无常, 尤其是潜在石漠化, 好转的速率极小, 而恶化则十分迅速, 因此需要加强对轻度以下石漠化的监察与防护。

4 讨论 4.1 石漠化时空演变分析

1) 1990年至1997年间, 石漠化面积减少了71.0 km², 实际上是在好转, 而在1997年至2006年间, 石漠化面积又快速增加了217.0 km², 尤其是潜在向轻度石漠化的转变面积最大, 速度最快, 这段时间正处于经济发展的黄金期, 人类不合理活动强度的增大是导致恶化的主要原因, 直到2005年后, 政府开始实施退耕还林政策, 大量石漠化土地被禁止开发, 石漠化面积才又开始逐渐减少^[23-24], 从2006年至2015年, 石漠化面积又减少了110.0 km²。整体来看, 从1990年到2015年, 石漠化面积增加了34.5 km², 普定县石漠化经历了由好转到恶化再到好转的演变, 总体趋势虽然在好转, 但局部地区仍在恶化, 石漠化的治理不可掉以轻心。

2) 石漠化转移变化平均速率最快的是潜在石漠化, 平均速率为2.75 km²/a, 其中在1997—2006年间, 速率达到了16.5 km²/a；变化速率最慢的是极重度石漠化, 平均速率为0.0002 km²/a, 这也验证了前人的工作成果, 石漠化变化速率与等级高低成反比。此外, 好转速率最快的是中度石漠化, 平均速率为0.80 km²/a, 恶化速率最快的是无石漠化, 平均速率为3.10 km²/a, 由图 17看出, 轻度以下石漠化变化方向不易确定, 而中度以上石漠化则几乎都在好转, 因此, 在治理中度以上石漠化的同时, 更要注重轻度以下尤其是潜在石漠化的防治工作。

3) 结合石漠化因子转移变化图与石漠化空间变化图可以看出, 植被覆盖度与岩石裸露率的变化成相反关系, 在植被覆盖度增加高的区域, 石漠化演变多以好转为主, 多集中在普定县中部、东南一带；而在岩石裸露率增加高的区域, 石漠化演变多以恶化为主, 多集中在三岔河一带。

4) 由图 8分析得出, 在石漠化空间分布上, 本研究发现普定县中度以上石漠化多集中分布在三岔河沿岸附近, 而轻度以下石漠化则分布较散乱, 无太大规律。

5) 由图 13、图 14看出, 在三岔河附近, 常发生着无或潜在石漠化突变为重度石漠化的现象, 这往往是由河道拓宽, 雨水冲刷, 导致无或潜在石漠化地区植被剧减, 岩石裸露率暴增导致。此外, 在城关镇靠近县城一带, 也会出现重度石漠化突变为无石漠化的现象, 笔者认为由于此处地势较低, 易于土层堆积, 人类的开垦荒地等城镇化活动使得重度或极重度石漠化也变为了无或潜在石漠化, 从而这也体现了人类活动对石漠化演变的积极一面。

4.2 展望与不足

本文使用的数据由于依赖于遥感影像, 因此难免受到卫星过境时间和天气云量等的影响, 加之获得的植被覆盖度、岩石裸露率等常受到季节变化的干扰, 很难去计算全年的平均值, 此为最大的局限性。本文采取的方法在总体上反演石漠化精度较高, 因子权重较为合理, 在喀斯特中小尺度地区内反演石漠化具有一定的优势。

5 结论

本文基于Landsat、DEM等数据提取得到的普定县1990—2015年4期石漠化影像展开分析, 重点讨论了普定县的石漠化空间演变规律, 主要认识如下：1)在1990—2015年间, 普定县石漠化演变复杂剧烈, 经历了由好转到恶化再到好转的阶段, 总体演变趋势虽然好转, 但局部地区仍旧在恶化, 石漠化治理力度需要加大。2)普定县中度以上石漠化集中分布在三岔河沿岸附近, 轻度以下石漠化则无明显地理分布特征。3)潜在和轻度石漠化演变方向复杂, 潜在石漠化变化速率最大(2.75 km²/a), 最大时在1997—2006年间达到了16.5 km²/a。4)植被覆盖度与岩石裸露率的变化呈负相关, 前者主导了石漠化的好转, 后者控制了石漠化的恶化。5)石漠化演变存在突变, 在三岔河沿岸附近, 存在无/潜在石漠化有突变为重度石漠化的现象；在靠近城乡居民地附近, 也存在重度以上石漠化突变为无石漠化的现象。6)植被覆盖度、岩石裸露率、坡度等能较准确反演石漠化, 贡献率分别为44%、42%、14%。