石漠化是指在潜在自然因素和人类不合理活动共同作用下，自然植被景观被严重破坏，土地生产力急剧下降，基岩大面积裸露，地表出现类似荒漠景观的极端土地退化现象^[1]。石漠化地区基本意味着生存环境的丧失，由于耕地遭受破坏、粮食产量下降、生态系统功能退化，往往陷入自然灾害频发与经济衰退的恶性循环^[2]。因此，精确获取石漠化信息对区域生态安全和经济社会持续发展具有重要意义。

遥感是以电磁波和地球表面物质相互作用为基础，探测、分析和研究地球资源与环境，揭示地表各要素空间分布特征与时空变化规律的一门科学技术^[3]，具有大面积同步观测、时效性强、受人为干扰小等优势，已成为石漠化信息提取的重要手段。目前，国内外学者针对石漠化问题，从提取方法^[4-7]、演变特征^[8-9]、驱动因素^[10-11]等方面开展了大量研究并取得诸多成果。但由于研究视角、方法、区域差异等原因，多数研究仍存在以下不足：1)高精度研究空间范围较小，对大范围石漠化信息提取研究不够全面；2)影像时相、指标选择单一，未构建长时序数据对石漠化演变特征进行分析；3)石漠化区域影像阴影清晰，光谱混合现象严重，仅通过植被指数等光谱特征具有很大的局限性^[12]。已有研究表明，引入纹理、地形等多特征进行提取分析，能够有效解决石漠化地区因地形复杂而导致提取精度不足的问题^[13]，同时Sentinel-2影像因其在时间、空间、光谱分辨上的优势，被广泛应用于农作物提取和地物分类中^[14]。当前，以谷歌地球引擎(Google earth engine, GEE)为代表的云平台的出现，打破了传统遥感数据获取和预处理的方式，其强大的并行运算能力和海量在线遥感数据集为大区域、长时序、高时空分辨率遥感研究提供了新的思路^[15]。

针对复杂地形石漠化监测中存在的人工调查周期长、费时费力、分类特征单一、区域尺度提取结果精度低等问题，借助GEE平台、以Sentinel-2影像作为主要数据源、结合DEM数据，提出一种基于多特征组合优选方案的石漠化信息提取方法，以期为石漠化防治提供理论参考和科学依据。

云南省昭通市(E 102°52′~105°19′, N 26°55′~ 28°36′)地处滇东北金沙江下游乌蒙山区腹地、四川盆地与云贵高原过渡地带(图 1)，辖9县1市1区，行政区划面积2.3万km²。境内地势西南高、东北低，最低海拔267 m，最高4 040 m。受地形和季风影响，兼有高原、温带、亚热带3种气候特征，年均气温介于11~21 ℃之间，年均降雨量介于660~1 230 mm之间。植被以亚热带针叶林、阔叶林、灌木林为主，区内土壤类型复杂多样，主要包括红壤、黄壤、棕壤、黄棕壤等。黄棕壤是主导土壤类型，广泛分布于昭通全境，黄壤主要分布在威信、镇雄一带，红壤集中分布在巧家、鲁甸县。典型的峰丛洼地地貌为石漠化形成和扩张提供了物质条件。石漠化较为严重的区域主要位于昭通市西南地区，其他区域均有分布。

图 1(Fig. 1)

图 1 研究区位置和高程 Fig. 1 Location and elevation of the study area

Sentinel-2多光谱(multi spectral instrument, MSI)影像来源于GEE平台，该平台融合谷歌服务器提供的云计算能力以实现对庞大数据集的在线调用计算。Sentinel-2 MSI影像产品级别为经过正射校正和大气校正的Level-2A级，空间分辨率为10、20和60 m。一颗卫星的重访周期为10 d，2颗互补，重访周期为5 d，其中QA60去云波段可有效去除山区多云、多雾天气造成的污染像元。研究区位于乌蒙山区腹地，局部地区海拔高度超过3 000 m，常年多云多雨，受云雨天气的影响，导致实际影像可用率低(含云量＜20%)。综合时间、空间、属性信息进行影像筛选，最终获得能够覆盖研究区全域的405景影像，时间范围集中于2020年1—5月之间，逐景去云后按月中值合成、拼接。

由NASA官网获取的ALOS DEM数据(https://search.earthdata.nasa.gov/search)，空间分辨率为12.5 m。为了与Sentinel-2影像空间分辨率保持一致，将其重采样至10 m后上传至GEE平台。

三套非同源土地覆被遥感监测数据产品包括欧空局和ESRI公司基于Sentinel-1和Sentinel-2影像制作的全球土地覆被数据产品；中科院空天院基于GEE平台，利用随机森林分类算法制作的先导专项土地覆被数据产品^[16]。以上3套产品用来确定随机森林分类算法的输入样本类型，各产品详细信息如表 1所示。

表 1(Tab. 1)

表 1 3套非同源土地覆被遥感监测数据产品 Tab. 1 Three non-homologous remote sensing monitoring land cover data products 产品名称 Product name 空间分辨率 Spatial resolution/m 类别数 Classification number 制作年份 Production year 制图精度 Mapping accuracy/% ESA_2020_10m 10 11 2020 74.4 Esri_2020_10m 10 10 2020 85.0 中科院先导专项数据产品 Leading special data products of Chinese Academy of Sciences 10 6 2020 > 90.0

表 1 3套非同源土地覆被遥感监测数据产品 Tab. 1 Three non-homologous remote sensing monitoring land cover data products

首先利用Sentinel-2反射率影像提取光谱指数特征，基于灰度共生矩阵(gray level co-occurrence matrix, GLCM)和主成分分析法(principal component analysis, PCA)提取纹理特征(texture feature, TF₁)，利用DEM数据提取地形特征，包括坡度(slope)和坡向(aspect)。其次，根据输入特征不同，构建以下5种方案：方案1包括波段、光谱指数特征；方案2在方案1的基础上加入植被覆盖度(fractional vegetation cover, FVC，C_fv)和基岩裸露率(bedrock exposure rate, BER，R_be)；方案3在方案2的基础上加入TF₁；方案4在方案3的基础上加入slope、aspect；方案5为以JM距离算法为基础的特征优选组合方案。同时，结合3套非同源土地覆被数据产品和Google Earth影像对输入样本进行优化。最后，利用随机森林分类算法完成石漠化土地提取并进行精度评价，研究方法流程如图 2所示。

图 2(Fig. 2)

GLCM: Gray level co-occurrence matrix. PCA: Principal component analysis 图 2 石漠化提取流程 Fig. 2 Flow chart of rocky desertification extraction

1) 波段特征。选择Sentinel-2可见光波段B2、B3、B4；红边波段B5、B6、B7、B8A；近红外波段B8；短波红外波段B11和短波红外波段B12共10个波段共同构成分类的波段特征。

2) 光谱指数。利用GEE平台计算以下光谱指数：归一化植被指数(normalized difference vegetation index, NDVI，I_ndv)、增强型植被指数(enhanced vegetation index, EVI, I_ev)、归一化水体指数(normalized difference water index, NDWI，I_ndw)、归一化岩石指数(normalized differential rocky index, NDRI，I_ndr)、绿色叶绿素指数(green chlorophyll vegetation index, GCVI，I_gcv)^[17]、地表水体指数(land surface water index, LSWI，I_lsw)^[18]、裸土指数(bare soil index, BSI，I_bs)以及反照率(Albedo，A_lbedo(α))^[19]。各光谱指数具体计算公式及说明如表 2所示。

表 2(Tab. 2)

表 2 光谱指数计算公式及特征说明 Tab. 2 Calculation formula and characteristic explanation of spectral indexes 光谱指数 Spectral index 计算公式 Calculation formula 特征说明 Feature description 归一化植被指数  NDVI Indv=(B8-B4)/(B8+B4) 反映植被生长状态  Reflecting vegetation growth status 增强型植被指数  EVI Iev=2.5×(B8-B4)/(B8+6×B4-7.5×B2+1) 反映植被冠层结构变化  Reflecting the changes of vegetation canopy structure 归一化水体指数  NDWI Indw=(B3-B8)/(B3+B8) 增加水体边界提取能力  Increasing the extraction capacity of water boundary 归一化岩石指数  NDRI Indr=(B11-B8)/(B11+B8) 提高裸露岩石识别能力  Improving the recognition ability of exposed rock 绿色叶绿素指数  GCVI Igcv=B8/B3-1 反映植被光合作用活性  Reflecting the photosynthetic activity of vegetation 地表水体指数  LSWI Ilsw=(B8-B11)/(B8+B11) 反映植被冠层及土壤表层含水量  Reflecting the water content of vegetation canopy and soil surface 裸土指数  BSI Ibs=((B11+B4)-(B8+B2))/((B11+B4)+(B8+B2)) 突出农业用地和非农业用地之间的区别，提高识别裸露土壤和休耕地的能力  Highlighting the distinction between agricultural and non-agricultural land and improving the ability to identify bare soil and fallow land 反照率  Albedo Albedo(α)=((0.356×B2)+(0.130×B4)+(0.373×B8)+(0.085×B11)+(0.072×B12)-0.018))/1.016 通过地表吸收的辐射能量反映土壤水分和暴露状况  Reflecting the soil moisture and exposure by surface absorption of radiation energy 注：B2为蓝光波段；B3为绿光波段；B4为红光波段；B8为近红外波段；B11为短波红外1；B12为短波红外2。Notes: B2 is the blue band; B3 is the green band; B4 is the red band; B8 is near-infrared band; B11 is shortwave infrared 1; and B12 is shortwave infrared 2.

表 2 光谱指数计算公式及特征说明 Tab. 2 Calculation formula and characteristic explanation of spectral indexes

3) 植被覆盖度和基岩裸露率。FVC和BER均采用像元二分模型，以NDVI和NDRI作为中间指数进行计算。FVC和BER计算式^[20]分别为：

式中：I_ndvmax和I_ndvmin分别为研究区内最大和最小的NDVI值；I_ndrmax和I_ndrmin分别为研究区内最大和最小的NDRI值。

4) 地形特征。利用GEE平台在线编辑代码，由DEM计算坡度、坡向，作为波段共同参与分类。

5) 纹理特征。基于GEE提供的GLCM计算出研究区影像的14个纹理特征统计量。选择角二阶矩(angular second moment, ASM)、对比度(contrast, CONTRAST)、相关性(correlation, CORR)、熵(entropy, ENT)、方差(variance, VAR)、逆差矩(inverse difference moment, IDM)、总和平均(sum average, SAVG) 7个纹理特征统计量，采用主成分分析，选择荷载比最高的TF₁作为最终参与分类的纹理特征。

根据研究区实际情况，将土地利用/覆被类型分为耕地、林地、草地、裸地、水体(河流、湖泊、水库等)、人造地表(建设用地、道路等)和石漠化土地7类。其中针对前6种地类，结合上述3套土地覆被数据集进行样本选取，将3种数据产品分类一致的纯像元作为选择范围，在范围内随机、均匀选取；石漠化土地样本借助Google Earth影像，通过目视解译选取。石漠化土地样本点在选择过程中采用多限制条件组合的方式依据石漠化等级划分标准^[21]，按照植被覆盖率(＜40%)、基岩裸露率(>60%)以及坡度(>18°)3个指标特征的综合信息来确保样本选择的准确性。最终，共选取样本数量为5 009个，其中水体209个、耕地817个，林地617个、草地718个、裸地659个、人造地表489个、石漠化土地1 500个。所有满足条件的样本随机分为70%和30%，分别作为随机森林分类算法的训练样本和验证样本数据。

在样本满足准确性、随机性和均匀性的前提下，采用JM(Jeffries-Matusita)距离对不同输入特征的可区分性进行评价并以此来确定最佳特征组合。JM距离J_m的具体计算式^[22]为

式中B为某一特征的巴氏距离(Bhattacharyya distance)。J_m可表示特征间的可分离程度，其值在0~2，值越大表示可分性越好。在特征满足正态分布的前提下，不同类别特征的巴氏距离B^[23]为

式中：e_k为k类特征的均值; δ_k²为k类特征的方差。

随机森林是一种机器学习算法，具有预测几千个解释变量的能力。它以决策树为单元，将多个决策树聚合起来进行分类，可以对数量庞大的较高维数据进行分类^[16]。与其他的分类算法相比，随机森林具有运算、训练速度快且不易过拟合的优点，并且算法内部暗含对每个指标分类效果的判别权重以突出利于分类的特征^[24]。

基于验证样本数据(图 3)，采用混淆矩阵，计算总体精度(overall accuracy，OA)、生产者精度(producer accuracy，PA)、用户精度(user accuracy，UA)和Kappa系数4个评价指标对石漠化土地提取结果进行检验。其中OA和Kappa系数用来评价总体分精度，PA和UA用来评价各土地覆被类型之间的错分、漏分误差。

图 3(Fig. 3)

图 3 验证样本空间分布图 Fig. 3 Spatial distribution of validated samples

为避免特征冗余对分类结果造成影响，采用JM距离算法对26个输入特征进行可分离性评价。计算训练数据集中每类样本点的均值和方差，得到不同输入特征下每种土地利用/覆被类型与石漠化土地之间的JM距离(表 3)，最终依据JM距离均值大小优选14个特征构建方案5。

表 3(Tab. 3)

表 3 不同输入特征下各地类与石漠化土地之间的JM距离 Tab. 3 JM distance between different input characteristics and rocky desertification 序号 Order No. 输入特征 Input feature 水体 Water bodies 林地 Forest 耕地 Cultivated land 草地 Grassland 人造地表 Artificial surfaces 裸地 Bare land 平均JM距离 Average JM distance 1 BSI 1.930 1.985 1.942 1.918 1.946 1.938 1.940 2 Albedo 1.923 1.961 1.938 1.953 1.947 1.927 1.940 3 TF1 1.880 1.900 1.920 1.930 1.940 1.870 1.900 4 Slope 1.900 1.840 1.800 1.880 1.840 1.820 1.840 5 B12 1.704 1.977 1.833 1.794 1.863 1.720 1.820 6 NDRI 1.792 1.810 1.824 1.756 1.883 1.810 1.813 7 B11 1.684 1.966 1.801 1.779 1.820 1.753 1.800 8 FVC 1.750 1.834 1.780 1.725 1.831 1.824 1.800 9 NDVI 1.749 1.811 1.798 1.746 1.820 1.807 1.789 10 B2 1.740 1.720 1.770 1.830 1.830 1.780 1.780 11 B4 1.690 1.740 1.860 1.760 1.820 1.720 1.770 12 BER 1.749 1.811 1.795 1.724 1.819 1.741 1.770 13 B6 1.710 1.870 1.700 1.780 1.770 1.720 1.760 14 B8 1.710 1.900 1.700 1.770 1.760 1.710 1.760

表 3 不同输入特征下各地类与石漠化土地之间的JM距离 Tab. 3 JM distance between different input characteristics and rocky desertification

由表 3中可见，在14个输入特征中，平均JM距离最大对应的输入特征是BSI、Albedo，其次是TF₁和slope，平均JM距离最小时对应的输入特征是B6和B8。石漠化土地与水体、林地、耕地、人造地表以及裸地分离度最高时对应的输入特征是BSI；与草地分离度最高和最低时对应的输入特征分别是Albedo和BER；综合来看，石漠化土地与其他土地覆被类型JM距离在1.9以上的输入特征(slope、TF₁、BSI和Albedo)对分类精度提升效果较显著，其他输入特征对分类精度提升也能够起到积极作用。

5种方案分类精度如表 4所示，可以看出不同特征组合得到的分类精度存在一定差异。方案1中石漠化土地的PA和UA均未达到85%，说明仅靠光谱特征与指数特征并不能够有效提取石漠化土地。方案3在方案2的基础上加入基于GLCM和PCA提取的纹理特征后，其总体分类精度和Kappa系数相比于方案1分别提高4.3%和0.05，其中石漠化土地的PA和UA分别提高1.7%和5.8%，说明纹理特征能够有效减少石漠化土地与其他地类之间的混分现象。方案4在方案3的基础加入地形特征后，总体精度和Kappa系数分别提高4%和0.04，其中石漠化土地的PA和UA相比于方案3分别提高8.5%和5.9%，说明地形特征能够显著提高石漠化土地的分类精度。方案4虽然能够得到较高的分类精度，但是参与分类的输入特征包括波段、光谱指数、纹理、地形等其他特征，存在特征冗余和运算量大的问题。方案5在方案4的基础上通过JM距离算法优选输入特征，不仅解决了特征冗余问题，且该方案的总体精度和Kappa系数分别达88.0%和0.85。

表 4(Tab. 4)

表 4 5种不同分类方案精度统计 Tab. 4 Accuracy statistics of 5 different classification schemes % 方案Scheme 精度指标Accurancy index 石漠化土地Rocky desertification land 水体Water bodies 耕地Cultivated land 草地Grass land 林地Forest 人造地表Artificial surfaces 裸地Bare land 总体精度Overall accuracy Kappa系数Kappa coefficients 1 PA 81.1 100.0 81.2 70.7 70.6 85.3 82.2 79.4 0.75 UA 72.1 97.9 78.5 81.7 84.1 80.2 87.6 2 PA 80.0 96.1 82.5 70.0 75.8 76.3 91.4 80.8 0.76 UA 76.8 100.0 78.5 78.0 80.0 81.3 84.6 3 PA 82.8 96.1 83.1 71.0 77.8 75.9 91.8 83.7 0.80 UA 77.9 99.0 81.0 80.4 79.2 84.0 84.1 4 PA 91.3 95.7 78.7 90.4 76.4 89.0 93.3 87.7 0.84 UA 83.8 100.0 85.7 91.6 87.9 89.4 90.0 5 PA 91.2 95.7 80.0 89.6 76.9 89.5 92.7 88.0 0.85 UA 83.8 100.0 91.7 87.1 87.1 87.9 90.5

表 4 5种不同分类方案精度统计 Tab. 4 Accuracy statistics of 5 different classification schemes

随机森林分类算法中最佳决策树的获取不仅能够提高分类精度而且可以减少不必要的计算量。针对方案5将决策树起始数量设置为10，终止数量设置为200，以10为步长执行循环迭代，最后得到分类精度与决策树数量点线图(图 4)，可以看出，决策树在10~40棵时，分类精度不断升高，40~80棵时分类精度产生波动，80棵时分类精度达到最大值，因此决策树最优值设置为80棵。

图 4(Fig. 4)

图 4 决策树数量与分类精度关系 Fig. 4 Relationship between number of decision trees and classification accuracy

对方案5中参与随森林分类的所有输入特征进行重要性分析(图 5)，可以看出，坡度特征对分类精度贡献最大，其次是纹理特征TF₁、NDVI和BSI，光谱特征中B4和B6波段贡献相对较小。

图 5(Fig. 5)

图 5 输入特征对分类精度的重要性排序 Fig. 5 Importance ranking of input feature to classification accuracy

表 5为多特征优选组合的方案5精度评价混淆矩阵，可以看出：1)草地的分类精度相对较低，同时草地、林地和裸地与石漠化土地之间存在混分现象，主要是因为植被覆盖度在同一年不同季节都会表现出明显差异，筛选到的可用影像时间范围集中在1—5月的旱季，这一时期研究区植被生长状态较差，易将草地误认为石漠化土地，导致样本点误选，造成混分；2)除草地外，其他土地覆被类型的PA值均在80%以上，其中石漠化土地、水体、裸地的PA能够达到90%以上。水体、林地、耕地、草地、人造地表以及裸地的UA都在87%以上，能够取得较好的分类结果；3)对于方案5，总体分类精度为88.0%，Kappa系数为0.85，石漠化土地的PA和UA均能够达到80%以上，满足分类精度要求。

表 5(Tab. 5)

表 5 方案5精度评价混淆矩阵 Tab. 5 Precision evaluation confusion matrix of scheme 5 类别 石漠化土地  Rocky desertification land 水体  Water bodies 林地  Forest 耕地  Cultivated land 草地  Grass land 人造地表  Artificial surfaces 裸地  Bare land 石漠化土地  Rocky desertification land 690 0 20 0 31 0 14 水体  Water bodies 0 90 0 0 0 3 1 林地  Forest 54 0 233 7 2 0 0 耕地  Cultivated land 10 0 16 349 3 7 1 草地  Grass land 64 0 3 16 268 0 0 人造地表  Artificial surfaces 0 0 0 7 0 203 18 裸地  Bare land 5 0 0 2 1 14 307 生产者精度  Producer accuracy/% 91.2 95.7 80.0 89.6 76.9 89.5 92.7 用户精度  User accuracy/% 83.8 100.0 91.7 87.1 87.1 87.9 90.5 总体精度  Overall accuracy/% 88.0 Kappa系数  Kappa coefficient 0.85

表 5 方案5精度评价混淆矩阵 Tab. 5 Precision evaluation confusion matrix of scheme 5

图 6a显示由方案5提取的昭通市石漠化总体空间分布情况。空间上，石漠化土地广泛分布于除昭通东北部以外的地区，且在巧家和鲁甸一带较为集中，全区石漠化土地面积为2 820 km²。分类结果图 6b显示，林地、耕地和草地面积分别占所有土地利用/覆被类型面积的43.98%、31.09%和23.07%，石漠化土地面积占11.11%，位于第四位，反映出研究区石漠化严重程度。从图 6c中可以看出，巧家县杨家坪村、龙潭村附近的石漠化土地均能够被准确识别，但也存在少量细碎的混分情况(包括“椒盐”现象)。总体上，分类结果比较符合实际情况。

图 6(Fig. 6)

图中红色代表石漠化土地提取结果，绿色实线代表石漠化土地边界范围。 In the figure, red indicates the extraction result of rocky desertification land, and green solid line indicates the boundary range of rocky desertification land 图 6 基于JM距离优选特征组合分类结果及石漠化典型区域对比 Fig. 6 Classification results of optimal feature combination based on JM distance and comparison of typical areas of rocky desertification

目前，利用遥感技术进行石漠化信息提取的相关研究对特征数量选择存在一定差异，例如李丽等^[25]以云南省隆林县为研究区，NDVI为提取特征，采用像元线性分解模型估算植被覆盖率的方法得到石漠化空间分布情况，但仅选用单一指标易造成石漠化土地与其他地类之间的混分；凌成星等^[7]以贵州省务川县为研究区，综合EVI、NDWI、碳酸盐岩裸露面积、地形坡度、植被覆盖率等多种特征对石漠化进行提取，最终提取精度达到88.2%。

本研究主要探讨利用多特征组合优选进行石漠化信息提取的方法。首先根据输入特征不同，构建5种分类方案，后一方案在前一方案的基础上逐步增加特征并进行精度对比；其次依据JM距离算法优选输入特征，确定最终分类方案5；最后采用最优特征组合完成研究区石漠化信息自动提取。与方案4对比，方案5总体精度和Kappa系数均有所提高，说明最优特征组合能够有效提高石漠化信息提取的精度。

笔者提出的研究方法兼顾光谱、纹理、地形等特征，在一定程度上克服了使用单一分类特征的缺陷，石漠化信息提取精度显著提高，但是仍有以下不足：1)受多云多雾天气的影响，研究区长时间序列的高质量光学遥感数据难以获得，可能会造成最终分类精度的不稳定性；2)随机森林分类算法得到的结果是基于像元提取的，难以避免“椒盐”现象发生；3)样本选择过程中结合的3套土地覆被数据产品存在影像时间与分类精度的差异，对最终分类精度有一定影响。结合合成孔径雷达不受云雾天气影响、穿透性和抗干扰性强的特点进行区域尺度的石漠化信息提取是下一步工作的重点。

1) 该方法综合考虑传统石漠化提取中存在的不确定性以及交叉信息对分类结果的影响，在区域尺度能够取得较高的分类精度，是石漠化信息提取的有效手段。

2) 结合土地利用/覆被数据与GEE高分辨率底图对样本选择进行优化，可以有效改善高原山地石漠化分布破碎区易出现错分、漏分、误差大的现象。JM距离算法在特征优化方面作用显著，不仅能够有效保留不同地类重要特征，还能解决特征冗余造成运算效率低的问题。