基于MODIS影像喀斯特石漠化状况研究——以贵州省为例

吴林霖，官云兰，李嘉伟，袁晨鑫，李 睿

(1.东华理工大学测绘工程学院，南昌 330013；2.流域生态与地理环境监测国家测绘地理信息局重点实验室，南昌 330013；3.广东南方数码科技股份有限公司，广州 510665)

0 引言

喀斯特地貌在我国又称为岩溶地貌，一般指碳酸盐岩分布地区或存在流经石灰岩的地下水所形成的特有地貌现象[1]。贵州省是我国喀斯特地貌分布最广的省份，比“喀斯特”名称起源地的前南斯拉夫的喀斯特地貌分布面积占比还大，其石漠化问题也最为突出[2]。由袁道先院士首先提出了石漠化的概念，即含有植被、土地覆盖的喀斯特地区经各类因素使得岩石裸露的过程，并指出了石漠化是中国南方亚热带喀斯特地区严峻的生态问题[3]。我国对石漠化重视较早，取得了一定的治理进展，研究方法也由最初的人工调查转变为基于3S集成技术的石漠化时空演化与动态变化分析，以及石漠化信息提取和目标识别。李水明等[4]基于K2结构学习算法结合样本分类信息获取石漠化特征子集，并以此提取石漠化信息；蓝安军等[5]分析了植被覆盖率与石漠化的负相关关系并引入石漠化动力指数来表征石漠化的分级和现状；李松等[6]基于1992年、2001年2个时期的TM/ETM+影像，在进行数据预处理后采用图像差值法对石漠化进行了变化检测，并提取了变化信息。当前对喀斯特石漠化的研究主要集中在石漠化分布状况、动态监测、预警分析、演化机制和生态治理等方面[7]，但大多是在小区域尺度范围内进行的，仅依据遥感提取的石漠化信息做单因素的石漠化特征反演，获得的成果有一定的局限性；并且由于选取的指标各不相同，研究结果难以进行相互间的分析、比较。

为此，本文利用MODIS影像大空间尺度的特点，综合考虑人地因素对石漠化的影响，引入4类石漠化评价因子(植被覆盖度(vegetation fractional cover，VFC)、基岩裸露度、坡度和人口密度(population density，PD))建立石漠化评价模型，研究贵州省近10 a间的石漠化分布状况，提取石漠化信息并对石漠化进行变化检测，进而分析石漠化变化特征与时空演化趋势。快速、高效地了解大范围的石漠化现状，监测石漠化变化情况是提供石漠化治理方向和改善生态环境的基础，在石漠化研究中具有重大的意义。

1 研究区及数据源

1.1 研究区概况

贵州省地处长江和珠江两大水系的分水岭地区，云贵高原东部，气候温暖湿润，属亚热带湿润季风气候[8]。全省总面积为176 167 km2，地势西部最高，中部稍低，并分别向北部四川盆地、东部湖南低山丘陵和南部广西盆地逐步过渡，东西呈三级阶梯，南北有两大斜坡，平均海拔为1 100 m，最大高差达2 763 m，形成了一个总体地势较高，内部分异较大，深受河流切割的亚热带喀斯特高原山区[9]。贵州省是一个碳酸盐岩广泛分布、喀斯特强烈发育的区域，碳酸盐岩总厚度约17 000 m，占全省沉积岩总厚度的70%以上，其中纯碳酸盐岩厚度约为碳酸盐岩总厚度的62%[10-11]。

1.2 数据源

本文使用的数据包括遥感影像、数字高程模型(digital elevation model，DEM)数据、人口数据及地理矢量数据等。遥感影像采用EOS系列卫星上搭载的MODIS中分辨率成像光谱仪传感器获取的陆地标准产品MOD09A1 500M 地表反射率8 d合成产品，提供了1—7波段的500 m空间分辨率影像数据，参数如表1所示。选取2007年、2010年和2016年获取时间相近的3期影像进行后续处理。

表1 MODIS 500 m地表反射率8 d合成产品参数Tab.1 MODIS 500 m surface reflectance 8-days composition product

DEM数据为SRTM(shuttle Radar topography mission)90 m空间分辨率原始高程数据，由美国国家航空航天局和国防部国家图像测绘局联合测量。DEM数据与MODIS数据均来源于中国科学院计算机网络信息中心地理空间数据云平台(http://www.gscloud.cn)。人口数据来源于贵州省统计局《贵州省统计年鉴》获取的2007年、2010年和2016年各市(州)年末常住人口。地理矢量数据来源于国家基础地理数据库，包括省界和市界。

2 评价因子

目前并没有确定的石漠化评价体系，而遥感石漠化监测需要考虑遥感技术的特点和石漠化评价的应用需求[12]。基岩裸露度和VFC是石漠化最直观的地面表现特征，故为石漠化评价的关键指标[13]。坡度会影响土壤侵蚀，是诱发石漠化的原因之一[14]。石漠化是自然因素和人为因素共同作用的结果，人类活动对自然环境有着不可磨灭的影响，人口剧增给土地带来的压力以及不合理的人类活动，会加速石漠化的进程[15]。因此，本文选取VFC、基岩裸露度、坡度和PD作为评价因子，根据层次分析法(analytic hierarchy process,AHP)建立石漠化评价模型，提取石漠化信息。

2.1 VFC

VFC是指植被(包括叶、茎、枝)在地面的垂直投影面积占统计区总面积的百分比[16]。岳跃民等[17]已验证植被指数对于提取石漠化信息的有效性，故本文选用归一化差值植被指数(normalized difference vegetation index，NDVI)估算VFC，并以李苗苗[18]在像元二分法基础上研究的VFC提取方法建立模型。由于岩溶地区土被分布浅薄且不连续，在提取VFC时可以暂不考虑土壤的影响[19]。本文假设研究区可近似取VFCmax=100%，VFCmin=0，则VFC的计算公式为

(1)

式中NDVImax和NDVImin分别表示研究区NDVI最大和最小值。由于噪声的影响，NDVImax和NDVImin一般会取一定置信度范围内的最大和最小值。根据置信区间5%～95%，使用ENVI软件查看不同时期NDVI直方图，获取不同时期NDVI最大和最小值如表2所示。

表2 NDVI最大和最小值Tab.2 Maximum and minimum of NDVI

将表中的NDVImax和NDVImin代入式(1)，可估算出研究区2007年、2010年和2016年的VFC。

2.2 基岩裸露度

基岩裸露度采用Rikimaru[20]提出的裸土指数(bare soil index，BSI)估算。由于裸土的光谱特征在中红外波段反射率最高，BSI通过高反射率波段减去低反射率波段从而突出裸土信息，并且加入了蓝光波段进行归一化，计算结果中裸土的亮度(DN值)最高，其次为建筑用地，易于区分出裸土、植被与水体。BSI公式为

(2)

式中r1，r2，r3和r6分别为MODIS数据红光、近红外、蓝光和短波红外波段的反射率值。

2.3 坡度

坡度是地面陡缓的程度，是衡量石漠化的重要因子[14]。采用地理空间数据云平台DEM数据切割模型对DEM裁剪，再使用ENVI软件对DEM拼接，最后生成坡度图。坡度计算公式为

(3)

式中：Slope为坡度；fx和fy分别为水平方向和垂直方向高程变化率。

2.4 PD

PD是指单位面积土地上居住的人口数[15]。本研究的人口信息从贵州省统计年鉴获得(2007年、2010年、2016年各市(州)年末常住人口)，通过ArcGIS软件将信息添加入县市级矢量数据属性表中，再经由核密度分析计算研究区各市级单位的人口分布情况，获取研究区的PD。

2.5 评价因子归一化

不同的评价因子数值大小不在同一量级，量纲不统一，直接对石漠化程度进行评价较为困难。因此，需要对各评价因子进行归一化，使其具有相同的量纲进而参与石漠化评价模型计算。将编码值设为1—10，统一各评价因子的量值，以数值大小表示对石漠化的影响程度，影响越大，编码值越大；影响越小，编码值则越小。坡度、基岩裸露度和PD对石漠化起正向影响，数值从低到高分为1—10的等级。VFC值越高表示植被越茂密，对石漠化起反向影响；VFC值越高，石漠化程度越小，反之则越大。因此，将VFC值从高到低分为1—10的等级。

3 石漠化评价模型建立

石漠化评价模型是基于各类归一化评价因子，通过比较评价因子重要性来综合计算各评价因子权重而建立的，以Saaty[21]提出的AHP来确定评价因子权值。该方法是一种定性和定量相结合的目标决策分析方法。首先，根据相对重要性的比例标度评判各因子的重要性建立判断矩阵，比例标度如表3所示；其次，经过一系列计算逐步剔除主观元素，使计算的权值尽可能客观公正；最后，所求出的权值还需要对其进行一致性检验，计算出一致性检验系数CR，当CR<0.1时，判断矩阵通过一致性检验，否则不满足一致性。

表3 比例标度Tab.3 Proportion criteria

①8,6,4,2,1/2,1/4,1/6,1/8为以上重要性标度的中间值。

具体计算步骤如下：

1)根据表3所示的比例标度结合研究区实际情况将评价因子两两进行比较确定相对重要性数值来建立判断矩阵A，如表4所示。

表4 判断矩阵及权值Tab.4 Judgment matrix and weight

2)计算判断矩阵各行的乘积，即

(4)

式中：i和j分别为判断矩阵的行和列号；n取值为4，因为选择了4个评价因子建立模型。

3)计算Mi的n次方根，即

(5)

(6)

式中wi为第i行对应的评价指标权值。

5)将上述结果代入公式(7)求出最大特征根，即

(7)

式中(Aw)i为矩阵A第i行所有元素与向量w第i行对应权值的乘积。

6)计算一致性系数CR，检验一致性。其表达式为

(8)

(9)

式中RI通过查表可得，当n=4时，RI=0.89。当CR<0.1时，通过一致性检验，否则修改判断矩阵比例标度，直到通过一致性检验。根据上述步骤计算出CR=0.060 726 704，通过一致性检验，说明判断矩阵比例标度合理，权值可行。根据计算的权值建立石漠化评价模型D，即

D=w1VFC+w2BSI+w3Slope+w4PD，

(10)

式中：w1=0.55，为VFC所对应的权值；w2=0.28，为BSI所对应的权值；w3=0.12，为坡度所对应的权值；w4=0.05，为PD所对应的权值。

4 石漠化状况分析与评价

4.1 石漠化分级标准

对石漠化程度分级，基于不同的影像数据和评价因子有不同的分级标准。刘芳等[22]在Landsat8热红外影像基础上，确定不同石漠化程度的亮温差，并依据亮温差建立了评价标准；李丽等[23]基于NDVI估算VFC，并依据研究区的VFC建立了石漠化评价标准。本文沿用现行的分级方法[24]，将石漠化程度划分为4个等级，即无石漠化、轻度石漠化、中度石漠化和重度石漠化。根据归一化编码值(1—10)分为4类，如表5所示。

表5 石漠化分级标准Tab.5 Hierarchy standard of rocky desertification extent

4.2 3个时期石漠化状况分析

根据上述方法及分级标准，得到2007年、2010年及2016年3期贵州省不同等级的石漠化空间分布状况，如图1所示。

(a)2007年 (b)2010年 (c)2016年

图1 石漠化分级结果

Fig.1Classificationofrockydesertificationextent

1)2007年石漠化分布。2007年全省大部分区域为轻度石漠化，其次为中度石漠化。轻度石漠化占全省总面积的50.80%；中度石漠化与重度石漠化占全省总面积的比例分别为42.16%和6.99%，其中凯里市和遵义市南部存在较少的石漠化现象，中度石漠化主要分布在遵义市东北部、毕节市和六盘水市，兴义市、贵阳市也有少量分布；重度石漠化主要分布在毕节市，少数分布于安顺市与六盘水市的交界区、遵义市和贵阳市。

2)2010年石漠化分布。相较于2007年而言，全省范围的重度石漠化区域有所减少，约占全省总面积的3.69%；中度石漠化约占全省面积的40.53%；而轻度石漠化则占全省面积55.67%。其中毕节市依旧表现为大面积的中度石漠化，重度石漠化相较2007年有部分减少；六盘水市和兴义市西部中度、重度石漠化分布范围增加；都匀市与遵义市主要为中度石漠化和极小片的重度石漠化；安顺市与凯里市基本为轻度石漠化和无石漠化。

3)2016年石漠化分布。2016年的石漠化情况明显改善了许多，轻度石漠化占全省面积的79.71%，中度石漠化和重度石漠化分别减少到了17.35%和2.87%。其中毕节市由重度石漠化大范围分布转变为轻度石漠化，使轻度石漠化所占比重明显增加，重度石漠化仍主要分布在西部地区；遵义市、铜仁市和都匀市的石漠化情况较为相似，有少许中度石漠化，而大范围分布轻度石漠化；贵阳市重度石漠化分布没有明显变化，说明贵阳市重度石漠化程度较稳定；兴义市西北部中度石漠化分布较多；凯里市主要为轻度石漠化，中度石漠化有零散分布。

4.3 石漠化总体变化

通过图像差值法，将对应像元的灰度值相减，获取不同时相的差异以达到对2景不同时期影像进行石漠化变化检测的目的。对2007年、2010年和2016年3个时期的石漠化评价结果两两进行差值，将差值图像阈值设为：-1为石漠化减少、1为石漠化增加，结果如图2所示。

(a)2007—2010年 (b)2010—2016年

图2 石漠化变化检测

Fig.2ChangedetectionofKarstrockydesertification

由图2(a)可以看出，2007—2010年间，石漠化减少区域明显大于石漠化增加区域。在贵州省西部(六盘水市、兴义市)，中部(贵阳市、都匀市)和东北部(遵义市、铜仁市)等有少量石漠化增加现象，而贵州西北部(毕节市)石漠化明显减少。如图2(b)所示，2010—2016年间，石漠化趋势主要表现为减少，大部分减少区域位于贵州省西北部(毕节市、六盘水市)，东北部(遵义市、铜仁市)，中部、西南部有少许减少趋势；增加区域主要在贵州省西南部(兴义市、安顺市和六盘水市)，与2007—2010年相比，呈现持续增加趋势。

4.4 石漠化分级变化

通过比较分级区域的变化，可以检测出轻度石漠化、中度石漠化和重度石漠化相互间的转移变化情况，2007—2010年和2010—2016年期间石漠化转移情况如图3所示。

(a)2007—2010年 (b)2010—2016年

图3 石漠化转移情况

Fig.3Conversiondistributionofrockydesertification

由图3可以看出，在2007—2010年间，轻度石漠化总体增加趋势最明显，由中度和重度石漠化转移为轻度石漠化的比例分别为8.15%和0.25%，但仍有6.14%的轻度石漠化转化为中度石漠化；重度石漠化转为中度石漠化的比例仅有2.29%，表明这期间对石漠化的治理有一定成效，但还是存在石漠化程度加重的情况。在2010—2016年期间，石漠化现象得到明显改善，其中中度石漠化向轻度石漠化的转移比例为14.61%，转换为重度石漠化的比例仅为0.58%；轻度石漠化转化为中度石漠化的比例也只有2.59%。

对比2个阶段的石漠化转移情况，发现轻度转向中度石漠化的比例远远大于轻度转为重度石漠化的比例，重度石漠化转为中度石漠化的比例大于转为轻度石漠化的比例。该现象说明中度石漠化为重度石漠化与轻度石漠化相互转化的中间过程，从石漠化治理角度可以先由重度石漠化改善为中度石漠化，再由中度石漠化改善为轻度石漠化；中度石漠化与轻度石漠化相互之间的转化占比较大，说明这2种石漠化类型具有跳跃性；重度石漠化转移面积比较少，说明重度石漠化具有稳定性，难以改善。这与白晓永等[9]、赵丽苹[25]发现石漠化转移及演替的特征基本一致。

5 讨论

在熊康宁等[26]的研究中先将研究区划分为喀斯特区与非喀斯特区。本文考虑到贵州省碳酸盐岩几乎分布全境，石漠化区域也比较大，所以提取石漠化信息时没有先剔除非喀斯特区，而是将整个研究区作为喀斯特区域。由于选择的研究区范围较大，选择MODIS中分辨率成像光谱仪产品较为适合，并且选用MODIS地表反射率8 d合成产品，避免复杂的数据预处理过程。但由于MODIS影像的空间分辨率较低，不可避免的存在混合像元，加之研究区域范围较大，且没有实地考察来验证本文研究的提取精度。因此，模型计算结果只能宏观表示研究区石漠化概况与变化趋势。评价因子计算中，李霞等[27]曾基于SPOT影像，采用归一化土壤指数(normalized difference soil index，NDSI)和不透水面指数(normalized difference impervious surface index，NDISI)双指数法确定一定经验阈值区分裸土与建筑用地。但是，MODIS影像尚不具备区分建筑用地与裸土的空间分辨率，故本文只采用BSI单一指数估算了基岩裸露度。在今后的研究中可以尝试结合其他方法提高估算精度。

6 结论

基于MODIS地表反射率8 d合成产品，选取植被覆盖度、基岩裸露度、坡度和人口密度作为评价因子，采用层次分析法计算各因子权值，通过一致性检验后建立石漠化评价模型，确定分级标准划分石漠化程度，得到了不同时期的石漠化评价结果；对2007年、2010年和2016年3个时期的石漠化评价结果两两进行变化检测，获得石漠化的总体变化以及石漠化分级变化状况；完成了对贵州省石漠化状况的评价，并分析了石漠化变化趋势及时空演化特征。取得的主要结论如下：

1)研究区石漠化在近10 a间已得到了较大的改善，2010—2016年期间石漠化治理效果尤为突出。由2010年轻度石漠化与中度石漠化分别占全省总面积的55.67%和40.53%，到2016年轻度石漠化与中度石漠化分别占全省总面积的79.71%和17.35%，中度石漠化同比减少了23.18%，轻度石漠化增加了24.04%。重度石漠化有明显的改善，从2007年占全省总面积的6.99%降到了2016年的2.87%。

2)在石漠化等级转化过程中，中度石漠化与轻度石漠化相互之间的转化率较高，说明这2种石漠化类型具有跳跃性；中度石漠化为重度石漠化与轻度石漠化相互转化的中间过程，从石漠化治理角度可以先由重度石漠化改善为中度石漠化，再由中度石漠化改善为轻度石漠化。

3)对比2007—2010年和2010—2016年期间的石漠化转移矩阵，重度石漠化转移面积最小，表明重度石漠化较为稳定，改善难度较大。2007—2010年间由重度石漠化转化为中度石漠化和轻度石漠化所占比例分别为0.25%和2.29%，2010—2016年间由重度石漠化转化为中度石漠化和轻度石漠化所占比例分别为0.59%和0.45%。