基于Landsat8巫山石漠化区域监测识别方法探讨
2021-06-14陈艳英徐彦平杨琴游扬声
陈艳英 徐彦平 杨琴 游扬声
（1 重庆市气象科学研究所，重庆 401147；2 重庆市开州区突发事件预警信息发布中心，重庆 405400；3 重庆市气候中心，重庆 401147；4重庆大学土木学院，重庆 400045）
0 引言
石漠化也称石山化，指石灰岩山区植被强烈退化后岩石裸露的现象。岩溶石漠化是在热带、亚热带湿润、半湿润气候条件和岩溶极其发育的自然背景下，地表植被受人为活动干扰遭受破坏，造成土壤严重侵蚀，基岩大面积裸露，土地退化的表现形式。随着我国对生态文明建设的日益重视，石漠化也逐渐为公众所知晓。石漠化被喻为“地球的癌症”，反映了其治理的艰巨性。
20 世纪90年代以来，遥感技术在石漠化演变过程、石漠化的分布、石漠化提取等研究中得到广泛应用。李丽等以云南隆林各族自治县为研究区，结合土地利用现状，根据植被覆盖度提取了研究区石漠化分布状况；凌成星等基于中巴资源一号02B星数据，用改进的增强植被指数（EVI）提取石漠化信息，并使用水体指数（NDWI）修正水体信息，结合空间坡度图层进行决策树分析，对贵州省务川仡佬族苗族自治县进行石漠化提取，抽样检验表明石漠化提取精度达到88.2%；李伟等以植被类型、梯土化旱地因子、地形坡度、基岩裸露度、植被综合盖度5个因素为石漠化判识因子，建立石漠化识别条件，将研究区分为石漠化区域、潜在石漠化区域及非石漠化区域三类，并对岩基裸露度、地形坡度、植被综合盖度（结合旱地分布）、植被类型4类因子进行评分，从低到高将石漠化区域分为轻度石漠化区域、中度石漠化区域、重度石漠化区域、极重度石漠化区域4类，对广西百色市凌云县进行石漠化监测识别，以野外实地调查为佐证，验证了该方法具有较高准确性及可行性；苏欣利用地质数据、归一化植被指数、归一化岩石指数、植被覆盖度、土壤厚度等数据，构造石漠化提取指标，对重庆涪陵地区的10个乡镇街道进行了石漠化提取；黄岩等基于美国陆地卫星4～5号专题制图仪（thematic mapper） 所获取的多波段扫描影像提取植被覆盖度数据，并结合坡度数据、岩性背景资料及土壤厚度数据等，建立石漠化提取指标，对重庆南川进行了石漠化的提取，并进行了精度的验证；朱大运等以GF-1和Landsat-OLI为源数据，运用欧氏距离，以贵州关岭—贞丰花江一带为研究区，对多种植被指数在石漠化提取过程中的可分性和类型识别能力进行了定量的对比评价，结果表明，石漠化信息提取中推荐使用的最优植被指数均为归一化植被指数（NDVI），其次为长势监测指数（GRNDVI）。
在进行石漠化监测识别过程中，通常根据数据的获取情况，采用不同的数据进行石漠化提取，选取对石漠化识别比较敏感的数据进行研究，例如土地利用现状、植被覆盖度、归一化植被指数、增强型植被指数、地理国情地表覆盖、坡度、100 m 数字高程数据、岩石裸露度、岩土性质等数据。由此可见，在石漠化识别中，指标的选取并不固定。本文选取植被覆盖度、土地利用类型、岩石裸露度、岩土类型数据对巫山石漠化区域进行识别。
1 研究区概况
重庆岩溶地区有37个区县，辖888个乡镇、101个街道办事处。截至2006年末，岩溶地区总人口共3055万人，其中农业人口占76.23%，全年实现地区生产总值约2923.1亿元，粮食总产量为967.7万吨。重庆石漠化典型的区域主要分布于渝东北及渝东南地区，其余地区零星分布。
巫山县位于重庆市东北部（图1），幅员面积为2958 km，地处三峡库区腹心，介于30.767°—31.467°N，109.55°—110.183°E。长江横贯东西，大宁河、抱龙河等7 条支流呈南北纵向分布。深谷和中低山相间分布，地形起伏大，最高海拔达2698 m，最低海拔仅有63 m，属典型的喀斯特地貌。
全县山地占总面积约96%，丘陵平坝仅占4%。全县农业人口约占93.3%，人地关系紧张，土地过垦较严重。坡旱地分布广，大于25°的陡坡地占比大；森林覆盖率小，仅为12.97%；荒山草坡面积大，且土薄、坡陡，成片天然草场少，多呈块状分布。大量的良田沃土被交通、居住等建筑用地占用，农业生产与土地资源间存在较大的反差。
图1 巫山位置及高程图Fig. 1 Wushan location and elevation map
巫山属于贫困山区农业县，是长江三峡库区重点淹没县之一，喀斯特地貌遍布全县，地表破碎。境内水土流失严重，地质灾害及气象次生灾害频繁。石漠化是制约巫山发展的主要因素，碳酸盐岩面积约为991 km，占巫山县总面积的33.5%。2017年，巫山石漠化土地面积达到811.7 km。
2 资料与方法
2.1 研究资料的整理
1）Landsat8数据处理
Landsat8数据来源于地理空间数据云网站（http://www.gscloud.cn/）。
巫山县包含Landsat8数据中的第38行、126列和第39行、126列。由于第39行、126列在巫山只有红椿乡南部及庙宇镇东南角，范围较小。为了便于数据处理，选择了第38行、126列的数据后，通过ENVI软件对LandSat8数据进行大气校正、辐射定标处理，最终挑选了2013年10月24日、2014年10月11日、2015年10月14日、2017年5月28日、2017年11月4日共5组数据对巫山进行石漠化提取。
2）植被净初级生产力数据
整理国家卫星中心提供的植被净初级生产力（NPP）数据，计算巫山2000—2018年NPP及植被固碳量。
3）植被覆盖度
对植被指数数据逐月产品MYD13A3（来源于NASA官网）进行投影、格式转换、裁剪，合成年归一化植被指数（NDVI）数据，并利用NDVI计算出年植被覆盖度（VC），裁剪出巫山区域，形成2000—2018年巫山植被覆盖度数据序列。
4）地理信息数据
收集整理2016年土地利用数据和2004年岩土类型数据（来源于中科院北京数字空间科技有限公司）、100 数字高程数据（DEM）和重庆行政区划数据（来源于中国气象局信息中心）。
2.2 研究方法
2.2.1 植被覆盖度、岩石裸露度的计算
在ENVI软件中可实现基于Landsat数据的归一化植被指数NDVI和植被覆盖度VC的计算。
归一化植被指数NDVI采用可见光红光波段（波长在0.76～0.90 μm ）和近红外波段（波长在0.63～0.69 μm）计算，公式如下：
植被覆盖度计算公式如下：
VC
为植被覆盖度，N
_为格点植被指数，N
_为裸地或无植被覆盖区域植被指数，N
_为完全植被覆盖区域植被指数。此处N
_和N
_分别取NDVI累计百分率分别为5%和95%时的对应取值。当N
_＜N
_时，VC
=0；N
_＞N
_，VC
=1；N
_≤N
_≤N
_，采用公式（2）计算植被覆盖度。基于ENVI软件，利用公式（1）计算巫山不同日期植被指数NDVI。表1为NDVI累计百分率为5%和95%时对应的NDVI数值。
表1 5%和95%累计百分率对应NDVI值Table 1 NDVI values at cumulative percentages of 5%and 95%
根据表1及公式（2）计算巫山5次卫星过境时的植被覆盖度。
基于植被覆盖度计算岩土裸露度R
，R
=1-VC
，VC
为植被覆盖度。R
为0～1，值越大，说明岩石裸露越多。2.2.2 石漠化区域判别标准的建立
岩石裸露度、植被覆盖度、平均土层厚度、坡度、岩土类型、地质分布等数据常被用到，研究选用的指标略有区别，本文在参照已有研究的基础上，根据现有数据，选取了植被覆盖度、岩土裸露度、岩土类型、土地利用类型建立石漠化识别指标。
对植被覆盖度、岩石裸露度两个指标进行重分类并评分，表2给出了植被覆盖度、岩石裸露度的评分结果。表2中，植被覆盖度大（或者水体、硬化地等人工地表），评分较低；其余植被覆盖度较小的非水体和建筑地表，评分较高。岩石裸露度小（或者岩石裸露度接近1的水体、硬化地等人工地表），评分较低；其余岩石裸露度较大的非水体和建筑地表，评分较高。
导乐分娩时一种全新的分娩模式,有专业导乐人员为产妇提供专业化、人性化的服务,无创伤,给予产妇持续性的生理上的支持帮助及心理上的安慰鼓励,使其顺利完成分娩过程,有利于提高产科服务质量,降低剖宫产率,是国际上比较推崇的服务理念。自2017年1月我院开展导乐分娩以来,取得较好的收效,先总结如下。
表2 植被覆盖度、岩石裸露度评分标准Table 2 Scoring standards for vegetation coverage and rock exposure
表3 给出了土地利用类型的评分。评分值高，则现石漠化的概率高、石漠化程度重；评分值低，则出现石漠化的概率低、石漠化程度轻或无石漠化。表4中的地貌类型评分同理。
表3 土地利用类型评分表Table 3 Land use type scoring standard
表4 地貌类型评分表Table 4 Landform type scoring standard
对以上4个指标的评分求和，依据求和结果给出石漠化分级，表5给出巫山石漠化区域评分。
表5 巫山石漠化分级评分表Table 5 Grading of Wushan rocky desertification
在此基础上，将植被覆盖度、岩石裸露率和植被类型数据叠加运算，设置下列7个条件：①基岩裸露度＜30%；②有林地、灌木林地、疏林地、牧草地；③水田；建设用地；水域；④基岩裸露度≥30%；⑤植被综合盖度≥50%的有林地、灌木林地；⑥植被综合盖度≥70%的牧草地；⑦梯土化旱地。非石漠化和潜在石漠化区域分别叠加运算，参与叠加运算的所有数据投影、分辨率、行列数要求一致。
非石漠化区域判别：满足①和②；或满足③。
潜在石漠化区域判别：满足④和⑤；或满足④和⑥；或满足④和⑦。
最后，结合表5及非石漠化区域、潜在石漠化区域的判别条件，将巫山分为非石漠化区域、潜在石漠化区域、轻度石漠化区域、中度石漠化区域、重度石漠化区域、极重度石漠化区域6个等级。
2.2.3 生态植被监测分析
利用MODIS/AQUA逐月植被指数产品MYD13A3，采用最大值合成法合成季节植被指数，对季节植被指数取均值，得到格点年平均植被指数。
NDVI
为格点上的年植被指数；NDVI
为格点上的季节植被指数；n
为1、2、3、4，代表4个季节。格点年植被净初级生产力NPP
计算如下：
NPP
为格点上的年植被净初级生产力；NPP
为格点上的月植被净初级生产力；n
为1～12的整数，代表月份。计算区域内年平均植被覆盖度、植被净初级生产力，见公式（5）。
AVG
为植被覆盖度或植被净初级生产力的区域均值，Grid
为植被覆盖度或植被净初级生产力格点值，count
为计算区域内的格点点数。3 石漠化监测结果分析
利用选取的5次过境数据提取石漠化区域，给出了非石漠化区域、潜在石漠化区域、轻度石漠化、中度石漠化、重度石漠化、极重度石漠化6个等级的面积及每次提取的石漠化区域的总面积，见表6。
由表6可见，石漠化6个分区中非石漠化区域面积最大，其次是潜在石漠化区域，极重度和重度石漠化区域面积分居第3及第4位，轻度和中度石漠化区域分居第5及第6位。取5次监测结果均值作为标准值，计算每次监测结果的相对误差，见表7。
从表7可见，监测面积呈波动变化。5次监测中，非石漠化区域及潜在石漠化区域面积较稳定，分别在0.81%～1.1%、0.31%～17.18%波动；极重度石漠化区域变化最大，在17.09%～31.59%波动；轻度、中度、重度石漠化区域面积分别在9.94%～29.19%、0.06%～24.95%、1.66%～18.53% 波动，其中轻度石漠化区域面积波动范围仅次于极重度石漠化区域。
监测面积变化的影响因素主要来自遥感数据和地表信息的变化。卫星过境时的姿态使挂载仪器产生倾斜，仪器的倾斜导致传感器在同一个区域接收的地表反射信息产生变化（在地表起伏大的区域，传感器接收的地表反射信息变化更加明显），从而影响遥感监测指标；卫星过境时的光强对地表反射率也会产生影响，从而对大范围的监测指标产生影响。其次，地表反射信息影响遥感数据计算的石漠化分级指标，从而影响石漠化提取结果。同时，农作物的收割、休作及植被的不同生长季也会影响遥感监测指标及石漠化提取结果。
对5次石漠化计算指标进行加权平均并分级，得到石漠化分布结果，见图2。将该结果作为石漠化提取结果（见表8），并依此对石漠化区域进行分析。
加权后监测石漠化区域面积为875.9 km，非石漠化和潜在石漠化区域面积分别为1302.8 km、771.1 km。其中红椿乡西南部由于缺乏遥感数据，该乡南部石漠化分类结果有待进一步确认。
由图2可见，巫山石漠化区域分布与地形及土地利用相关，产生6个集中分布区：1）集中分布于大宁河流域及其支流小三峡区域；2）集中分布在七耀山山脉边缘至巫山县城东部呈东北—西南带状的区域；3）集中在省道105西部—官渡河—长江女峰下游区域，呈东北—西南带状分布；4）集中分布在巫山东南笃坪乡南部及雪花乡；5）集中分布在巫山东北部当阳乡、庙堂乡北部区域；6）河流右岸的狭长区域也有石漠化分布。非石漠化主要分布在县城、水域、林地等区域，潜在石漠化区域主要分布在草地、灌木地、农业用地区域。
图2 巫山石漠化监测结果Fig. 2 Wushan rocky desertification monitoring results
表8 巫山石漠化分级加权面积统计结果Table 8 Statistical results of Wushan rocky desertification weighted area
4 石漠化区域生态环境分析
将监测结果分为三类，即非漠化区域、潜在石漠化区域、石漠化区域（包含轻度石漠化、中度石漠化、重度石漠化、极重度石漠化4个分级），计算三个区域2001—2018年MODIS植被覆盖度、植被净初级生产力NPP、植被固碳量，见图3、图4。
表6 巫山石漠化分级面积统计Table 6 Statistics of Wushan rocky desertification classification
表7 巫山石漠化分级相对误差Table 7 Relative error of Wushan rocky desertification classification
从图3可见，巫山石漠化区域、潜在石漠化区域、非石漠化区域植被覆盖均呈增加趋势，非石漠化区域植被覆盖度最大，在60.5%～69.5%；潜在石漠化区域植被覆盖度略低于非石漠化区域植被覆盖度，在56.5%～66.0%；石漠化区域植被覆盖度最低，在54.5%～64.3%。植被覆盖度最小、最大值分别出现在2001和2017年。
从图4可见，石漠化区域、潜在石漠化区域、非石漠化区域的植被净初级生产力及固碳量差异明显，其中植被净初级生产力及固碳量两个指标在石漠化区域最小，在非石漠化区域最大；三个分区的植被净初级生产力及固碳量在2001—2018年呈现波动上升的趋势，且变化趋势基本一致。
5 总结与展望
5.1 总结
1）利用Landsat8遥感数据计算的植被覆盖度及岩石裸露度、土地利用数据、岩土类型数据建立了石漠化监测识别指标，提取巫山石漠化区域，结果分为6个等级。石漠化区域面积为875.9 km（包含轻度石漠化、中度石漠化、重度石漠化、极重度石漠化4个分级），非石漠化和潜在石漠化区域面积分别为1302.8和771.1 km。
图3 巫山非漠化区域、潜在石漠化区域、石漠化区域2001—2018年植被覆盖度Fig.3 The vegetation coverage of Wushan nondesertification area, potential rocky desertification area,and rocky desertification area from 2001 to 2018
图4 巫山2001—2018年植被净初级生 产力（NPP）及植被固碳量Fig.4 Vegetation net primary productivity and vegetation carbon sequestration in Wushan from 2001 to 2018
2）巫山非石漠化区域植被生态环境最好，其次是潜在石漠化区域，石漠化区域植被生态环境逊色于非石漠化区域及潜在石漠化区域生态环境。
5.2 不足与展望
在利用Landsat8数据、土地分类数据、地貌数据对巫山进行石漠化监测识别过程中，发现以下问题有待讨论：1）土地利用数据的精度及时效性对监测结果产生的影响不容忽视；2）在水域、建筑物及道路等硬化地、待开发用地等人为活动引起的植被破坏的区域，植被覆盖少，当土地利用数据得不到及时更新时，容易与石漠化区域混淆；3）在作物间歇期的农田区域，植被覆盖度偏低、岩土裸露度偏高，容易与石漠化区域混淆而造成误判。针对以上问题，后续可在同一年内收集中多组数据，对不同时期的数据融合，利用融合后的数据进行滚动识别、相互交叉验证，从而提升识别的精度。
