基于Google Earth Engine的典型峰丛洼地石漠化时空演变与驱动因子分析——以西畴县为例

徐红枫, 王 妍, 刘云根, 肖羽芯, 黎 舟

(1.西南林业大学, 昆明 650224; 2.云南省山地农村生态环境演变与污染治理重点实验室, 昆明 650224)

我国西南岩溶地区是石漠化分布最广的区域[1]，其主要地貌特征是地表起伏大、多山地[2]、多碳酸盐岩地区，土壤层浅薄，土地贫瘠[3]。根据《岩溶地区石漠化综合治理规划大纲(2006—2015)》区划，峰丛洼地石漠化治理区是我国石漠化综合治理八大类型区之一，该区出露碳酸盐岩古老、坚硬、层厚质纯，且连片分布，水文系统具有典型的二元结构，地表水发达而地下水发育，因此缺水、少土，耕地资源匮乏，人地矛盾的冲突尤为严重[4]。西畴县地处云南省文山壮族苗族自治州中部，位于云贵高原向广西盆地过渡的斜坡地带，是典型的峰丛洼地地区，石漠化问题一直制约着西畴县的经济、社会、环境可持续发展[5]。30年来，西畴县持续进行石漠化治理工作，人们通过退耕还林、坡改梯、经济林营造等一列石漠化治理工程[6]，探索出了“六子登科”、“五法治水”等石漠化综合治理模式，不仅取得良好的生态和社会效益，西畴人民在治理石漠化所产生的西畴精神也被列为云南脱贫攻坚精神的典型代表。

近些年来已经有不少学者针对石漠化信息提取及驱动因子进行了相关研究，在石漠化信息提取方面苏旺德等以ETM，DEM和地质岩性数据构建了南汀河石漠化[7];罗旭玲等认为基于NPP，NDVI、地表反照率和坡度可以较为精准地反演石漠化[8];为了揭示石漠化的驱动因子，李阳兵等基于人机交互解译，对黔中高原后寨河地区的石漠化演变驱动因子进行了分析，得出较多坡耕地的存在是石漠化的驱动因素[9]；王正雄等基于地理探测器对岩溶槽谷地区的石漠化驱动因素进行分析，得出岩性和土地利用、岩性和坡度的组合共同驱动岩溶槽谷地区的石漠化形成[10]。以往的研究多使用传统的方法进行数据处理，处理效率不高，在大区域和长时间跨度的研究上局限性较大，基于云计算的可用性，本文借助Google Earth Engine云计算平台，以典型峰丛洼地石漠化地区西畴县作为研究区，旨在揭示其1990—2020年30年间的石漠化演变规律和驱动因子，为后续的石漠化演变及更大尺度生态恢复研究提供一定的参考和技术指导。

1 研究区和研究方法

1.1 研究区概况

西畴县位于东经104°22′—104°58′，北纬23°06′—23°37′，其下辖有2个镇、7个乡：西洒镇、兴街镇、蚌谷乡、莲花塘乡、新马街乡、柏林乡、法斗乡、董马乡、鸡街乡，总面积1 506 km2。其中石漠山区面积达99.9%，岩溶面积达75.4%，境内岩溶峰丛溶蚀洼地，一般连片出现，山峦起伏，上层基岩裸露，地形封闭破碎，数百个水淹洼塘和大小不等的山间小盆地镶于丛山峦障中。且洼地多为椭圆或条形，宽数十米至数百米；溶洞、竖并、天窗等岩溶现象较为普遍。

1.2 数据源与研究方法

1.2.1 Google Earth Engine简介 Google Earth Engine(GEE)是谷歌旗下的云计算平台，可以快速、批量处理大量数据[11]。GEE上提供了BP数量级的公开地球观测数据，有超过500万张影像，200多个公共数据集可供免费使用[12]。GEE的API接口支持JavaScript和python，可以通过编程调用数值、数组、矩阵、栅格数据和图表等相关的多种模块和函数，用于实现影像预处理、图像分类、影像时间序列分析等操作。网址(https:∥earthengine.google.com/)。

1.2.2 数据源与研究方法 Landsat陆地观测卫星具有30 m良好空间分辨率，且具有长的时间覆盖，所以基于Google Earth Engine调用长时间序列的Landast影像作为数据源，对西畴县1990—2020年近30年的石漠化演变进行分析。其中植被覆盖度和岩石裸露率使用了Landsat 5/8 Surface Reflectance Tier 1数据集，该数据集已经使用LaSRC进行了大气校正，道路和建筑用地的提取使用了Landsat 5/8 Collection 1 Tier 1 TOA Reflectance数据集。坡度数据在GEE上基于SRTM Digital Elevation 30 m数据计算得到，年均降水量数据单位为0.1 mm，人口密度数据单位为人/km2，GDP数据单位为万元/km2，以上栅格数据均下载自中科院资源环境科学与数据中心(网址为：http:∥www.resdc.cn)。

基于GEE的超强计算能力，以及其提供的BP级免费公开遥感数据，通过JavaScript API编程在线调用Landsat 5影像309景和Landsat8影像174景，以云分量20为阈值进行影像筛选、然后去云、进行影像集的自动拼接和裁剪、对水体进行掩膜，并基于像元二分模型计算得到研究区植被覆盖度、岩石裸露率，通过计算指数和纹理特征实现道路和建筑用地的提取，最后将结果导出，下载到本地使用ENVI 5.3和ArcGIS 10.5进行后续处理，通过叠加制图得到了西畴县1990—2020年石漠化等级分布图，整体流程见图1。

图1 石漠化制图流程图

2 石漠化指标

2.1 植被覆盖度

(1) 归一化植被指数(NDVI)。石漠化过程往往伴随着植被的退化[13]。归一化植被指数(NDVI)可以较好地反映地表植被状况[14]，公式为：

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

(1)

式中:NDVI为像元归一化植被指数值；NIR为近红外波段；RED为红光波段。

(2) 运用像元二分模型[15]。在年度合成NDVI频率累积表上取频率5%对应的NDVI值作为NDVImin，取累积频率为95%的NDVI值作为NDVImax[16]，公式为：

FVC=(NDVI-NDVImin)/(NDVImax-NDVImin)

(2)

2.2 岩石裸露率

(1) 归一化岩石指数(NDRI)。基于石漠化地区的地物特性，张晓伦提出用Landsat数据的短波红外(SWIR)与近红外(NIR)波段，计算归一化岩石指数(NDRI)提取石漠化分布结果精度高[17]。公式为:

NDRI=(SWIR-NIR)/(SWIR+NIR)

(3)

(2) 岩石裸露率(Fr)。在年度合成NDRI频率累积表上取频率卫5%对应的NDRI值作为NDRImin，取累积频率为95%的NDRI值作为NDRImax，计算公式为：

Fr=(NDRI-NDRImin)/(NDRImax-NDRImin)

(4)

2.3 建筑用地与道路提取

基于Google Earth Engine，以Landsat 5 Collection 1 Tier 1 TOA Reflectance(1990年、2000年和2010年)和Landsat 8 Collection 1 Tier 1 TOA Reflectance为数据源(2020年)提取西畴县不透水面(建筑用地和道路)，计算归一化差异水体指数(MNDWI)，来提取水体[18]，计算归一化建筑指数(NDBI)来提取部分建筑特征[19]，并运用灰度共生矩阵(GLCM)计算其他纹理特征[20]，然后选取样本点，取其中70%的样本点作为训练样本，使用GEE上的随机森林分类器进行训练，另外的30%样本点作为验证样本，生成混淆矩阵，并计算Kappa系数。样本点分成两类：第一类为不透水面，主要包括道路和建筑用地；第2类为透水面，主要包括植被、裸地、农用地水体等。样本点的拾取都在GEE平台上完成，2020年样本的拾取基于Google高精度在线地图，1990年、2000年、2010年样本点的拾取则是基于Landsat5。

Kappa系数通常用于一致性检验，也可以用于衡量分类的精度，Kappa系数越大，分类精度越高。公式如下：

(5)

式中:po是每一类正确分类的样本数量之和除以总样本数,也就是总体分类精度；pe代表偶然一致性,通常Kappa系数在0～0.200代表极低的一致性,0.210～0.400代表一般的一致性,0.400～0.600代表中等的一致性,0.610～0.800代表较高的一致性,0.810～1代表几乎完全一致。

不透水面提取结果如下：1990—2020年四期不透水面提取的总体精度分别为0.816，0.860，0.865，0.982，Kappa系数分别为0.627，0.707，0.730，0.963。1990年、2000年、2010年三期提取精度在0.62～0.80之间，2020年的提取精度在0.963四期提取结果一致性较高，基本可以满足研究需要。

2.4 水体掩膜

归一化植被指数范围在-1至1之间，通常地表云、雪以及水体覆盖时，近红外波段的反射率低于红波段的反射率，所以NDVI值为负数，所以在GEE上计算植被覆盖度和岩石裸露率时对NDVI小于0的部分进行掩膜，从而达到剔除水体的目的。

在参考前人标准[21]的基础上，将石漠化划分为无、潜在、轻度、中度、重度和极重度石漠化6个等级，并剔除了水体、道路和建设用地，得到了1990年、2000年、2010年、2020年总共四期的石漠化等级分布图(图2)。在Google Earth Pro上验证解译精度，然后制作混淆矩阵，并计算Kappa系数，2020年使用的验证影像来源为Google Earth高精度影像，1990年、2000年、2010年使用的影像数据源为Landsat 5。四期石漠化等级分布图的总体精度分别为0.715，0.735，0.778，0.895，Kappa系数分别为0.681，0.671，0.745，0.721。1990—2020年四期石漠化等级分布图精度均在0.610～0.800，一致性较高。

图2 1990-2020年石漠化等级分布

3 结果与分析

3.1 石漠化时间变化特征

由图3可知，第一阶段1990—2000年的10年间各级石漠化面积都有一定减少，无石漠化土地增加明显，增加了19.349%；在第二阶段2000—2010年，2010年石漠化分布情况相较2000年又有所恶化，无石漠化土地面积减少明显，减少了22.874%，重度石漠化面积增加了8.940%；在第三阶段，2010—2020年期间，西畴县的无石漠化地面积上升明显，较2010年面积增加了46.015%，潜在、轻度、中度、重度石漠化土地面积较2010年显著减少。

3.2 石漠化转移矩阵

转移概率矩阵是指客观事物从一个状态转移到另一个状态，可以清晰表示不同时期同一区域石漠化类型相互转换的数量化关系[22]。为分析不同等级石漠化相互之间的转化关系，利用ArcGIS 10.5构建西畴县石漠化转移矩阵，在Origin 2021作图得到石漠化转移图(图4)。

图3 1990-2020年各级石漠化面积比重

从图4可以看出：(1) 1990—2000年西畴县的潜在石漠化转向无石漠化转变最显著，面积160.809 km2；其次，变化较大的是轻度石漠化转换成无石漠化面积84.767 km2；其余类型之间的转换面积相对较小，总体来说石漠化都在呈好转趋势。(2) 2000—2010年，无石漠化转变为潜在石漠化变化最显著，面积144.746 km2；无石漠化转变为轻度石漠化的面积次之，有115.608 km2；除此之外，潜在和轻度石漠化向中度石漠化转化面积较大分别为72.322 km2，75.303 km2，同时也有较大面积的中度石漠化向重度石漠化转化，有79.321 km2。(3) 2010—2020年主要趋势为潜在、轻度和中度石漠化转变为无石漠化石，面积分别为186.49 km2，195.916 km2，186.985 km2；同时也有66.168 km2重度漠化转化为无石漠化；(4) 1990—2020年潜在、轻度、中度、重度、极重度石漠化均有向无石漠化转变，无石漠化面积增加总计648.476 km2，截至2020年西畴县无石漠化土地面积共计940.854 km2；修复难度较大的极重度石漠化也得到了很好的修复，向重度、中度、轻度及潜在石漠化转变，其中有20.208 km2转变为无石漠化。

3.3 石漠化空间演变分析

为探究西畴县近30年石漠化空间演变格局，将石漠化等级降低的区域视为好转，石漠化等级提升的区域视为恶化，石漠化等级未发生变化的区域视为不变(图5)。总的来看，1990—2020年西畴县的石漠化有51.715 km2发生恶化，有850.183 km2石漠化土地发生好转，在空间上整体都呈现好转，也有少部分土地呈现恶化，主要分布在在莲花塘乡西北部、新马街乡的东北部，鸡街乡和西洒镇的交界处也存在部分零星的恶化土地。其中1990—2000年西畴县的东北部石漠化改善较大，包括鸡街乡、董马乡和法斗乡，西南部存在恶化，尤其是兴街镇和新马街乡，还有蚌谷乡和莲花塘乡交界处；2000—2010年北部和西部大面积石漠化土地都存在不同程度的恶化，尤其是董马乡、西洒镇、蚌谷乡以及法斗乡北部；2010—2020年西畴县石漠化总体是发生改善，尤其是西洒镇、蚌谷乡和兴街镇，蚌谷长箐片区、兴街江龙片区、罗汉冲片区、三光片区等石漠化综合治理示范区的设立治理效果十分显著。

图5 1990-2020年石漠化变化

3.4 贡献因子获取与分析

地理探测器(GeoDetector)是由王劲峰团队开发的地理探测器软件[23]。其理论核心是通过空间异质性来探测因变量与自变量之空间分布格局的一致性，据此度量因变量对自变量的解释度即q值。q的值域为[0，1]，值越大，说明这个因子对Y值的解释力越强。其模型：

(6)

为分析西畴县石漠化驱动因子，分别选取了西畴县2010年坡度、年均降水量、人口密度、GDP作为自变量[10]。在ArcGIS 10.5上将4个变量栅格采用自然间断法重分类成6类，然后用500 m网格点对石漠化类型分布栅格与变量栅格进行采样提取。

从单因子探测结果来看GDP(0.626)>年均降水(0.577)>人口密度(0.056)>坡度(0.022)，多因子交互探测如表1所示，坡度和GDP的交互作用下，对石漠化的影响最大(q=0.645)，GDP和人口密度对石漠化形成的影响则次之(q=0.639)。

表1 多因子交互探测结果

4 讨论与结论

4.1 讨 论

Google Earth Engine(GEE)存储有长时间序列的多源遥感数据，且计算能力强大，随着近几年云计算的兴起，越来越多的学者开始使用GEE云计算平台开展相关研究。张滔等[11]基于GEE对京津冀地区2001—2015年植被变化与城镇扩张的关系进行研究实现了数据的快速处理；裴杰等[24]利用GEE计算长时间序列的Landsat-NDVI对南洞地下河流域1988—2016年29年的植被覆盖度变化进行了监测研究；陈黔等[25]基于GEE与机器学习对毛乌素沙地灌木覆盖度进行估算，有效、快速地区分了灌木并对其植被覆盖度进行估算。在岩溶地区的大尺度石漠化监测方面，GEE可以解决数据收集困难、数据预处理复杂和效率低下的问题，与传统的石漠化遥感信息提取相比，利用GEE云计算平台提取石漠化表征因子，只需要通过代码的编写，就能进行大尺度遥感数据的批量处理与计算，使得在相同数量级数据处理下效率得到大幅提升，GEE势必成为实现大尺度石漠化监测的有效工具。

经过近30年的石漠化治理，2020年西畴县的无石漠化土地面积已经达到总面积的62.474%，石漠化治理取得了很大成效。1990—2000年的10年间，石漠化得到了一定的治理但是成效并不明显，可能的原因是石漠化治理虽然已经开展，但是由于缺乏统一的规划和资金，所以成效甚微。2000—2010年，石漠化情况有所恶化，此时西畴县进入了经济发展的新阶段，峰丛洼地大量耕地转化成建筑用地[26]，人们为谋求更多的生存空间，任意开采开荒，破坏林草植被，加剧了石漠化形成；后10年时间(2010—2020年)，西畴石漠化治理成效显著，这和政府的规划有很大的关系，西畴县一期石漠化综合治理工程(2011—2014年)，对磨合、三光、长箐3个小流域实施了封山育林、植树造林、坡改梯等措施，2016年制定《西畴县全面深化生态文明体制改革实施方案》，并相继实施了蚌谷长箐片区、兴街江龙、罗汉冲、三光片区等石漠化综合治理工程，探索出了“六子登科”、“五法治水”的石漠化综合治理模式。鉴于植物措施是巩固土壤、减少水土流失的有效措施[27]，后续应当继续推进石漠化治理工程，并加大对植树造林与经济林营造的扶持力度，在治理石漠化的同时增加人民的收入。

在驱动因子方面GDP(q=0.626)和年均降水(q=0.577)对西畴石漠化的形成贡献较大，GDP的贡献度最高，说明经济的发展对自然资源的影响较大，如果不能处理好经济发展与生态保护的关系，极易加剧石漠化的进程；其次是降水，其决定了植被生长的条件，尤其是西畴县属于典型的峰丛洼地，极度缺水，降水对植被的生长尤为重要。坡度对西畴县石漠化形成的贡献并不大(q=0.022)，但坡度与GDP交互作用下对石漠化的影响最大(q=0.645)，通常情况下，在适合人类活动的坡度范围内，坡度越陡，人类的活动越容易造成水土流失[28]；人口密度和GDP的交互作用对石漠化的形成影响则次之(q=0.639)，一般人口密度越大、GDP越高，经济活动越频繁，越容易造成植被破坏与水土流失，形成石漠化。

4.2 结 论

(1) 从时间上来看1990—2020年西畴县石漠化演变的总体趋势是先小幅减少，再扩张，然后大幅减少。1990—2000年石漠化转移的主要方向是潜在石漠化和轻度石漠化向无石漠化转移；2000—2010年，石漠转移方向是无石漠化向潜在和轻度石漠化转移；2010—2020年潜在、轻度和中度石漠化向无石漠化石转移最显著。从空间上来看西畴县石漠化分布北重南轻，在石漠化治理成效上，经过近30年的治理西畴石漠化整体都呈现好转，石漠化改善面积为850.183 km2，尤其是西洒镇、蚌谷乡和兴街镇，蚌谷长箐片区、兴街江龙片区、罗汉冲片区、三光片区等石漠化综合治理示范区治理效果十分显著。在莲花塘乡西北部、新马街乡的东北部、鸡街乡和西洒镇的交界处也存在零星恶化，恶化面积为51.715 km2。

(2) 截至2020年，西畴县无石漠化土地面积已经达940.854 km2，占西畴县总面积的62.474%，无石漠化面积增加总计648.476 km2。说明西畴县的石漠化治理工作较为成功，西畴县人民因地制宜探索出的“六子登科”、“五法治水”等石漠化综合治理模式效果明显。

(3) 地理探测器的结果显示年均降水量(q=0.626)和GDP(q=0.577)是驱动西畴县石漠化形成的主要驱动因子，多因子交互探测结果显示坡度和GDP交互作用(q=0.645)，人口密度和GDP交互作用(q=0.639)这两对组合共同驱动石漠化形成。

(4) 鉴于Google Earth Engine(GEE)云计算平台提供的BP级海量遥感数据和强大的计算能力，本研究利用GEE进行石漠化表征因子的快速提取，分析西畴县1990—2020年石漠化的动态变化及驱动因子，旨在利用GEE快速地实现大尺度地学问题分析，为后续大尺度的石漠化遥感制图提供有力的技术支持。