基于Landsat-7 ETM SLC-OFF影像的岩溶区石漠化研究

马友平

(湖北民族大学 林学园艺学院，湖北 恩施 445000)

石漠化是一种石质荒漠化，有地球癌症之称[1]，是土壤受到侵蚀后的直接后果与岩溶地区生态系统逆向演替的顶级阶段，一旦发生，恢复治理相当困难，导致生态环境的崩溃，且恢复速率极慢，已成为我国三大生态自然灾害之一[2-4].在我国主要分布在以云贵高原为中心的西南八省、市、自治区的451个县(市)中，这样大范围内有关石漠化基础数据的获取一直是石漠化研究的核心问题之一，其中遥感影像已成为石漠化研究的主要数据源，遥感影像主要依据的是地物的波谱特性，易受成像时间、地点、地物类别等的影响较大，现还没有一个大范围的统一标准，为此作者拟在该方面进行一些探讨，丰富以遥感影像为主要数据源的石漠化研究.

图1 恩施市位置图Fig.1 Enshi city position map

在石漠化的研究中,多以特定的行政区划范围为单位研究石漠化的发生机理、评定方法和绘制石漠化的等级分布图.如周忠发等[5]以贵州省清镇市为例，应用TM遥感影像制作了清镇市石漠化等级分布图;周涛[6]以贵州省关岭县为例，以ETM+、TM遥感影像为数据源，采用CA_MarkoV模型研究了2007年喀斯特石漠化分布情况，并对模型模拟情况进行评价;胡宝清等[7]以广西都安瑶族自治县为例，应用TM影像，建立数据库，分析影响因子，认为石漠化产生的主导原因是人类不合理的经营活动；朱林富等[8]应用TM影像和ArcGIS软件平台对重庆巫山县喀斯特石漠化与土地覆被关系进行了研究；这些都是以县域为单位.然而还有以省域或更大区域为对象的.如白晓永等[9]用三个时期TM影像对贵州省的石漠化演变过程进行了分析和评价；陈起伟等[10]以TM、ASTER和ALOS 三期卫星影像为遥感动态监测数据源对贵州省的石漠化进行了遥感监测；靖娟利等[11]利用MODIS遥感影像对广西自治区的喀斯特石漠化演变特征进行了研究；童立强等[12]应用ETM、TM、CBERS-02B、ALOS和ASTER卫星遥感数据对中国南方地区石漠化进行了遥感调查与演变分析.上述研究其基础数据的获取主要来自遥感影像，以遥感影像的光谱特性为基础，以不同波段之间的线性与非线性组合来提取能反映石漠化状况的指标.为此本研究拟在借鉴上述成果时以湖北省恩施市为例，讨论石漠化范围的确定和等级的划分标准；对于指标提取时影像的“缺失”区通过插值法获得，有研究证明遥感影像各点的灰度值(DN)，并不相互独立，它们之间存在一定的相关性[13-14]，为“缺失”区石漠化信息的地统计学估计提供了可能.

1 研究地区概况及数据来源

1.1 研究区概况

恩施市位于湖北省的西南部,如图1所示，是恩施土家族苗族自治州的首府，东经109°4′48″～109°58′42″、北纬29°50′33″～30°39′30″，面积397 2 km2,有世界硒都之美称，小气候复杂多变，石灰岩分布广泛，是南方岩溶发育区的一部分；虽经过石漠化的专项治理，石漠化现状有所改善，但问题依然存在，治理任务依然十分艰巨.

1.2 数据来源

美国Landsat系列数据现已免费开放，其中1999-04-15发射了Landsat-7，有8个波段，在波段设置上与Landsat-5非常相似，不同的是波段数上增加了一个全色波段，空间分辨率为15 m，在第6波段热红外波段上空间分辨率从120 m提高到了60 m.但在2003-05-31发生了故障，导致了一景影像中约有25%的数据丢失，因此这个时段以后的Landsat-7遥感数据称之为Landsat-7 ETM SLC-OFF数据，该研究主要采用有条带缺失的Landsat-7 ETM SLC-OFF影像，时间2015-10-22，轨道号/行号为126/39，空间分辨率30 m，云量0，而通过实际的图片观察恩施市上空有少量的云层；全市的DEM数据采用SRTM1(Shuttle Radar Topography Mission)，对应的分辨率精度为30 m，来源地理空间数据云(网址：http://www.gscloud.cn/).同时也收集一些有关恩施市的社会、经济、文化资料，主要为《恩施市统计年鉴》、《恩施州统计年鉴》.收集到恩施市的行政边界和区划图后，利用GPS采集地面控制点坐标进行几何校正，校正后坐标系均为WGS_1984_UTM_Zone_49N，与遥感影像、DEM为同一坐标系.

2 研究方法

对于大区域石漠化的研究，遥感影像作为基础数据源是许多专家学者的选择之一，该研究选择了Landsat-7 ETM SLC-OFF数据源，探讨在有云、阴影、条带缺失区的情形下与石漠化研究相关的数据修复问题，以弥补在应用过程中因云量、时段等导致的数据源缺乏的问题，下面同时结合其他数据源阐明遥感影像在石漠化应用中的研究方法.

2.1 Landsat-7 ETM SLC-OFF数据处理

按行政边界进行裁剪，分别将第1波段B1(Bands 1)、第4波段B4(Bands 4)、第7波段B7(Bands 7)放入红绿蓝(RGB)通道进行假彩色合成，提取影像数据的缺失区域，灰度级为1bit，正常区为1，缺失区数据为0，为机载扫描行校正器SLC(Scan Lines Corrector) 突然发生故障，导致获取的数据丢失部分；云覆盖的位置彩用梯度法[15-16]确定；阴影采用图像RGB中的G通道，HIS变换后的I通道进行阈值标定后，进行逻辑与运算所得[17]；在获取了图像的掉线、阴影、云的位置后进行彭胀一个象元即30 m的处理，如图2所示.图2(a)为第一、四、七波段合成的原图，图2(b)为去掉不正常区域后的一、四、七波段合成图形，经统计 “缺失”区面积1 503.829 0 km2，占整个国土面积的37.860 8%.

(a) B147原图 (b) B147“缺失”图图2 原图与“缺失”图Fig.2 The original image and the “missing” map

2.2 基于遥感影像的石漠化指数与基岩裸露率

石漠化信息的提取精度与方法，直接关系到研究区域石漠化的成果，经过查阅文献[18-20]、地面实际调查和比较，采用增强型植被指数法建立研究区域的石漠化指数模型，公式如下：

(1)

式(1)中：Di为第i像元的石漠化指数值，DNiB5、DNiB4分别为第i个像元在第5波段、第4波段中的灰度值，GB4、GB5分别为遥感影像第4、第5波段所有像元灰度值的平均值.

由式(1)提取了石漠化指数后，还要转化为基岩裸露率，基岩裸露率是石漠化等级划分的一个非常重的条件指标，计算公式为[20-21]：

(2)

式(2)中:Dgi为研究区域中第i个像元的基岩裸露率，Dmax、Dmin分别研究区域中石漠化指数的最大值与最小值，该研究中不包含“缺失”区的数据值.

2.3 石漠化等级的划分

石漠化主要分布在我国广大的南方地区，面积大，地域差异明显，由于研究方法和数据的差异，现在还没有一个统一的划分标准，作者在参考前辈的研究成果和在实际调查的基础上[12-13，22]，石漠化等级的划分标准为二项指标，即基岩裸露率和坡度；等级划分为四级，分别为无石漠化(I)、轻度石漠化(II)、中度石漠化(III)和重度石漠化(IV)，详见表1.

表1 岩溶区石漠化分级标准Tab.1 Classification standard of karst rock desertification

图3 图2(b)石漠化指数分布图Fig.3 Distribution map of karst rocky desertification index figure 2(b)

2.4 克里格法(Kriging)

克里格法(Kriging)又称为空间局部插值法，是在有限区域内对区域化变量的取值进行线性无偏最优估计的一种方法，是地统计学的重要研究内容之一，以变异函数来研究空间的结构性与随机性[14，23-24]，不同于经典统计学考虑了变量间的自相关性；而遥感影像为栅格数据结构，像元在图像上呈规则网状排列，特征值(石漠化指数)的提取只考虑了地物的灰度值，而灰度值存在一定的空间相关性，显然应用统计学中的克里格内插法可对“缺失”区域的数据进行预估，提高精度[25-26].

石漠化指数作为一项石漠化状况评价的基础指标，现还没有统一的评定方法和严格的定义，在大尺度研究时，多以遥感影像为数据源，通过光谱特征来计算，对于“缺失”区预测数据，无法同地形因子如高程等，与地面实测值相较，评价拟合结果，仅能将参与拟合的数据按随机的方法抽选一部分作为检验点，该研究中抽取5%的数据，而剩余95%的数据参与克里格法预测建模.建模精度拟采用三个指标评价：①决定系数(r)，用以反映观测值与估计值的关系紧密程度，值越大估测精度越高；②平均绝对误差(MAE)，反映估测值与检验值的偏离大小，其值越接近0效果越好；③均方根误差(RMSE)，介于0与1之间，取值越接近0，估测越精确.MAE和RMSE的计算公式如下:

(3)

(4)

式中：O实测值，P预测值，n为检验数据的个数.

3 结果与分析

3.1 石漠化指数的统计特征

石漠化指数的提取以ERDAS软件中的建模能功为工具，以图2(b)为遥感数据源，由式(1)可得图3.

图3中白色颜表示数据“缺失”区，数据值为NoData，将其转换为点状的矢量数据结构，属性值为石漠化指数，对其进行是否服从正态布检验.经试验和检验原始数据做对数变换后服从正态分布，满足克里格分析中简单克里格、普通克里格等要求数据服从正态分布的要求[27]，否则可能会使变异函数产生比例效应.

3.2 石漠化指数“缺失”的克里格估计

在使用遥感影像作为数据源时，由于地形起伏、成像时间、大气状况和成像传感器自身等因素，导致所获影像部分不是地物光谱特征的真实反映，有的甚至为数据缺失，如2003年后的Landsat-7影像，对这些区域的研究，统一称为数据“缺失”区，由于灰度值的自相关性[13-14]，用灰度值推演出的石漠化指数也应存在自相关性，因而对“缺失”区数据的估计采用地统计中的克里格法预测.

从图1中可以发现恩施市的行政边界极不规则，东西、南北、西北-东南、西南-东北方向的间距在60～100 km之间，石漠化指数的采样间距除去“缺失”区外为30 m，可认为是规则的采样点，可将采样间距或倍数设为滞后距离[28]，因此滞后距离设为3 km，滞后组数为12，变程为R=36 km，说明像素之间的距离小于该值时，它们之间空间相关，反之则认为无关，即相互独立.

由于原始数据经对数变换后服从正态分布，分别采用普通克里格法(Ordinary Kriging)、简单克立里格法(Simple Kriging)对“缺失”区进行数据估算，变程均为36 km，滞后距离3 km，滞后组数12，数据按随机原则从总体中抽取的训练数据，得半变异函数模型为指数模型.为了对其拟合精度检验证，利用5%的没有参与建模的数据，按照式(3)、(4)计算得图4、图5，可以发现无论决定系数(r)、平均绝对误差(MAE)、均方根误差(RMSE)均以普通克里格法最优，块金值C0=0.024 1，拱高C=0.008 1，变程R=36 km，半变异指数函数模型，如图6所示，由此可得石漠化指数的Kriging插值，如图7所示.

图4 普通克里格法预测石漠化指数精度图 图5 简单克里格法预测石漠化指数精度图

图6 石漠化指数半变异函数图 图7 石漠化指数的Kriging插值图

图8 石漠化等级分布图Fig.8 Distribution map of karst rocky desertification(K.R.D)

由上面的分析和图6可以发现基台值C0+C=0.032 2，C0/(C0+C)=74.844 7%，表示随机因素所引起的空间异质性占总变异的比例，该数值大，说明石漠化指数的空间相关性弱，反之则相关性强，由结构性引起的空间变异度较大.经研究表明该数值小于25%，空间相关性强，25%～75%空间相关性中等，大于75%中间相关性弱[13]，显然74.844 7%小于75%，属于中等空间相关.结构性因素岩性、气候和随机性因素土地耕作、林木采伐等对石漠化指数变异的影响差异不大.

3.3 石漠化等级的空间分布

图7为石漠化指数分布图，依据式(2)将其转变基岩裸露率分布图；坡度的获取从SRTM1数据中提取，然后据石漠化等级划分标准(表1)，可得恩施市2015年的石漠化等级分布(图8)，统计得轻度石漠化面积850.531 5 km2,中度石漠化面积27.267 3 km2，重度石漠化面积0.223 2 km2，石漠化的总面积为878.022 0 km2.

2015-10-22的Landsat-7遥感影像由于传感器自身的原因存在条带缺失，如果仅对条带缺失区域的数据采用克里格插值法进行完善，不考虑云、影阴的影响，可分为：情形①，不对影阴区进行克里格插值，用原影像的光谱信息分析则轻度石漠化为108.895 5 km2，中度石漠化为0.082 8 km2，而处理后的结果为轻度石漠化为130.309 2 km2，中度为6.585 3 km2；情形②，对于云区，不进行克里格处理，结果为轻度石漠化为4.471 2 km2,中度石漠化为0.234 0 km2，处理后轻度石漠化为2.686 5 km2，中度石漠化为0.042 3 km2,无重度,可以发现对于阴影区、云区二者结果差异较大.

4 结语

该研究选取了传感器存在缺陷的Landsat-7数据源，条带缺失明显，存在大面积的数据缺失，云层受大气状况的影响，所选影像全景云量很低近似为零；阴影由于地形起伏、太阳高度角、成像时间等的影响，作者在研究时对云、阴影的位置和范围进行了确定，作为数据“缺失”区处理.

地统计学分析时，得半变异函数为指数模型，变程为36 km，说明在36 km的范围内石漠化指数是存在空间相关性的，可以认为所有“缺失”区的图斑均在该范围内，应用克里格空间插值是存在意义的；对其预测精度的评估，现多采用在建模时预留一部分检证数据和参与建模的数据与预测值的误差大小来评价，而实际调查只能得到石漠化的状况和面积大小，而无法获得石漠化指数值.

对于条带缺失区的数据修复是必做的，云区、阴影区作者做了一个对比试验，即使用原有的光谱辐射值和克里格空间插值两种方法对比，结果波动性较大.对石漠化指数的空间结构分析表明，由于受气候、岩性、坡度和人们经营措施等的影响，研究区域内以石漠化指数为指标的空间结构呈中等空间自相关，恩施市的石漠化以轻度石漠化为主，面积850.531 5 km2，占国土总面积的21.413 2%，主要分布在研究区域的西北和东南地区，若采取有效的措施石漠化将得到有效的控制，转为无石漠化地区，反之则相反.