基于多源遥感数据的石漠化信息提取研究

陈 希，兰安军，熊康宁，廖建军

(贵州师范大学 喀斯特研究院/国家喀斯特石漠化防治工程技术研究中心，贵州 贵阳 550001)

0 引言

近年来，石漠化问题已经成为制约喀斯特地区发展的最关键问题．石漠化是指人类所实施的违反自然规律的社会生产经济活动所导致的喀斯特地区人地矛盾尖锐、植被破坏、水土流失、岩石逐渐裸露、土地生产力衰退丧失[1]，地表在视觉上呈现类似于荒漠化景观的演变过程[2]．其中，贵州省的喀斯特石漠化程度在全国范围内最为严重．目前，贵州省的石漠化土地面积达到4.2×104km2，占贵州省总面积的比例为22.9%，占贵州省喀斯特面积的37.1%[3]．因此，为了治理石漠化以及解决该现象所衍生出的诸多生态经济等难题，各界需要进一步加深对贵州省喀斯特石漠化的研究力度．

石漠化治理需要大量的生态环境数据、自然背景数据、社会经济数据等作为支撑，依托精准的遥感影像数据进行石漠化地区土地利用类型图、土壤侵蚀图等信息解译和石漠化等级判定是石漠化治理过程中必不可少的工作环节．但是贵州省喀斯特石漠化地区多云多雨的气候特征导致了在该地区难以获得有用的高质量的光学遥感数据．与光学遥感相比较，合成孔径雷达具备不受大雾或者云层影响的优点，穿透性和抗干扰性都较强，对纹理和地形的解译效果较好[4]．现阶段，国内对于石漠化信息提取的研究多为尝试利用光学遥感数据去开展，但很少涉及雷达遥感数据的利用，如2010年李丽等[5]利用中巴卫星光学遥感数据提出了一种基于植被覆盖度的石漠化信息提取方法，对云南隆林各族自治县完成了石漠化信息提取；2016年朱大运等[6]利用GF-1与Landsat-OLI光学遥感数据，对不同影像的植被指数分类能力进行了定量评价．

同时，国内对于利用雷达遥感数据去提取石漠化信息多基于单一时相的数据，解译效果并不理想．利用多时相雷达遥感数据的后向散射系数去构建时间序列的研究多集中于农作物信息提取、地表检测等，尝试其进行石漠化信息提取的研究较少，如2009年钱峻屏等[7]利用Radarsat雷达数据时间序列，通过密度异常检测的方法对南方多云雨地区的地面形变进行了检测；2015年钟礼山等[8]利用SAR时间序列对江苏省徐州市进行了耕地信息的提取研究.所以，目前利用雷达遥感数据的后向散射系数构建时间序列去提取石漠化信息和利用光学遥感数据与雷达遥感数据的融合进行石漠化强度定级的研究都尚属于探索阶段．基于上述，本研究尝试利用合成孔径雷达数据，通过构建后向散射系数值时间序列的方法进行研究区石漠化信息的提取，然后结合Landsat8遥感数据以及石漠化等级分类标准进行石漠化强度定级.

1 研究区概况与数据选择

1.1 研究区概况

本研究选取了“毕节撒拉溪喀斯特高原山地轻—中度石漠化综合治理示范区”和“关岭—贞丰花江喀斯特高原峡谷中—强度石漠化综合治理示范区”两个区域作为研究区，研究区分别位于贵州省的西北部和西南部地区，都具备了喀斯特地貌发育典型且喀斯特面积分布广泛的特点[9]，且各自代表了喀斯特高原山地轻-中度石漠化和喀斯特高原峡谷中-强度石漠化这两种典型的地貌，两个示范区的位置分布见图1.

图1 研究区位置图Fig. 1 Location of study areas

毕节撒拉溪喀斯特高原山地轻—中度石漠化综合治理示范区(简称“撒拉溪示范区”)位于毕节市西部的撒拉溪镇和野角乡境内，经纬跨度为105°1′00″E至105°09′30″E、27°39′11″N至27°18′00″N．研究区行政区划涉及2个乡镇，9个行政村，52个村民组，包括撒拉溪镇的朝营村、钟山村、冲锋村、永丰村、龙凤村、沙乐村、杨柳、撒拉溪和野角乡的茅坪村[10]．研究区属北亚热带湿润季风气候，植被类型以旱地、荒草地及稀疏灌木林地为主，土壤类型以沙壤和黄壤为主．区内土地利用类型以耕地和林地为主，水土流失整体状况以微度和中度为主．多年平均降雨量为984.40 mm，全年平均气温约为12 ℃．海拔高度为1500至2300 m，整体地势东高西缓，地表呈现斜坡丘陵自然延伸的形态，海拔差距较大，是典型的喀斯特高原山地轻度石漠化地区．

关岭—贞丰花江喀斯特高原峡谷中-强度石漠化综合治理示范区(简称“花江示范区”)位于贵州安顺市关岭县与黔西南自治州贞丰县交界处的北盘江峡谷花江段，经纬跨度为105°36′30″至105°46′30″E、25°39′13″至25°41′00″N．研究区行政区划涉及3个乡镇，8个行政村，50个村民组，包括北盘江镇的水淹坝、查尔岩、云洞湾和板贵乡的木工、坝山、三家寨、孔落箐以及花江镇的五里村[11]．研究区属亚热带季风湿润气候，植物类型多以次生灌木林为主，土壤类型以黄壤、黄色石灰土为主．区内土地利用类型以林地和耕地为主，水土流失整体状况以微度和轻度为主．多年平均降雨量为1100 mm，全年平均气温约为18.4 ℃，海拔高度为450至1450 m，整体地势陡峭，是典型的喀斯特高原峡谷强度石漠化地区．

1.2 雷达数据选择

本研究的SAR数据选择Sentinel-1所载的C波段SAR．Sentinel-1是GMES的重要组成部分，它由2颗卫星构成，重访周期为6 d，空间分辨率较高．它是搭载了C波段(5.405 GHz)的SAR传感器，它有4种成像模式：条带模式、干涉宽幅模式、超宽幅模式以及波模式．本研究中撒拉溪示范区的SAR数据获取时间分别为2016年1月12日、2月3日、3月15日、4月26日、5月23日、6月7日、7月13日、8月19日、9月30日、10月14日、11月13日和12月22日；花江示范区的SAR数据获取时间分别为2016年1月9日、2月14日、3月9日、4月14日、5月11日、6月15日、7月18日、8月21日、9月13日、10月17日、11月8日和12月11日.本研究的SAR数据均采用VV极化方式．

1.3 光学数据选择

本研究的光学数据选择Landsat8卫星所成的影像．Landsat8由美国国家航空航天局研发，它携带有两个主要荷载，分别是陆地成像仪(operational land imager)和热红外传感器(thermal infrared sensor)．Landsat8的陆地成像仪一共有9个波段，除全色的波段8空间分辨率是15 m外，其他波段的空间分辨率均为30 m．本研究中撒拉溪示范区和花江示范区分别选取2016年6月14日和10月4日的Landsat8陆地成像仪所成的影像作为光学遥感数据源．

2 提取方法

2.1 数据预处理

由于24景SAR数据的图像属性以及时相特征不一致，所以不能够直接把他们的图像信息有规律地联系起来．为了让各图像的原始信号幅度值能够准确量化地表达出来，使得不同时相的SAR数据具备可比性，必须要对影像进行辐射定标，将原始振幅信息转换为后向散射系数[12]．其公式为：

A·σo+Pn,

(1)

地形辐射校正的目的是通过消除遥感影像的几何变形，使多时相的SAR遥感数据能够互相匹配，从而提高数据的可对比性和可用性以及分析的准确性．本研究中的研究区生态环境复杂，地形破碎，加上SAR采用倾斜观测的成像方式，导致影像中地形畸变较大．本研究采用ASTER GDEM(先进星载热发射和反射辐射仪全球数字高程模型)进行地形辐射校正，其空间分辨率为30 m．本研究采用薛东剑等2017年提出的校正模型[13]，其公式为：

(2)

式中：β0为成像面散射单元对应的后向散射系数，α1为距离向坡度角，α2为方位向坡度角，δ为参考入射角，θ为当地入射角．

本研究采用精致Lee算子方法对SAR进行斑点噪声的滤波，这是一种基于边缘检测的算法．其原理是若检测范围内像元值较为一致，则可以认为该地区为匀质区域，此时利用低通滤波方法．反之，则采用高通滤波方法，这样可以最大程度上保留影像的纹理细节信息．

2.2 SAR后向散射系数时间序列的构建

SAR遥感数据后向散射系数时间序列构建是基于影像中同一个像元在不同时相数据中像元值的排列，为了能使各时相SAR数据的配准精度能够控制在0.5个像元误差之内[14]，需要在之前数据预处理的基础上进行几何精校配的步骤．几何精校配在ArcMap中进行，该操作以人工配准的方法，从两个示范区的12时相影像中选取质量较好的影像作为主影像，其余11时相影像作为辅影像，依次完成基于主影像的配准．设第m个点的配准误差为R，R的计算公式为[20]：

(3)

式中：x和y为第m个控制点的在主影像上的经度和纬度坐标，xi和yi是完成几何精校配之后第m个控制点在第i个时相在辅影像上的经度和纬度坐标．将配准后的24景影像按照示范区分区和时间顺序合成两幅影像，使得不同时相影像中的像元值一一对应．合成的影像中，X和Y轴分别代表经纬度，Z轴作为时间轴，那么每个像元则构成了后向散射系数的时间序列．SAR时间序列表达式为：

T={[(x1,y1),V1],[(x1,y1),V1],

…,[(xm×n,ym×n),Vm×n]},

(4)

Vi={((v1,t1),(v1,t1),…,(vn,tn)},

(5)

式中：(xi，yi)为第i个像元的行列坐标位置，m和n为SAR影像的行数和列数，Vi是基于第i个像元构建的时间序列，Vj为第i个像元处在时刻为tj的SAR影像后向散射系数值．

2.3 动态时间调整算法提取石漠化信息

DTW(dynamic time warping，动态时间调整)算法，能够对时间序列中的时间轴实现动态弯曲，对像元的后向散射系数在各个时间节点的变化具有很高的识别和匹配精度，它在解决尺度偏移的问题上具备一定的优势，而且在筛选和过滤畸变数值方面有较好的效果[15]．DTW算法提取石漠化信息的步骤为：首先，利用典型石漠化地物像元的后向散射系数时间序列作为石漠化地物标准时间序列；接着计算研究区内每个待分类像元的时间序列曲线与石漠化地物标准时间序列曲线之间的DTW距离，其算法公式为[16]：

D(1,1)=α11,

(6)

D(i,j)=αij+min{D(i-1,j-1),D(i,j-1),

D(i,j-1),D(i-1,j)}.

(7)

式中：i=2，3，…，m；j=2，3，…，n；D(i，j)为待分类像元的时间序列曲线与标准时间序列曲线距离的最小累加值，即DTW值，又称为最佳弯曲路径．然后，以目视解译的方法选取研究区石漠化地区边缘位置的若干像元，我们称之为混合像元[17]．接着计算混合像元在不同时间节点的后向散射系数值的平均值，这样就构建了一条新的基于12时相的混合像元时间序列曲线．最后计算混合像元时间序列和石漠化地物标准时间序列的DTW值，以该值作为研究区待分类地物类型判断的阈值[18]．若待分类像元的时间序列曲线与石漠化地物标准时间序列曲线的DTW值小于该阈值，则可以认为该像元属于石漠化地物类型；若待分类像元的时间序列曲线与石漠化地物标准时间序列曲线的DTW值大于该阈值，则可以认为该像不属于石漠化地物类型．

2.4 SAR与遥感数据融合的石漠化强度定级

将SAR与光学遥感数据在ENVI中进行融合操作，然后以3S(遥感技术、地理信息系统、全球定位系统)技术作为支撑，结合坡度、基岩裸露率、土壤信息、植被覆盖等作为主要判断因子，辅以平均降水量、土层厚度、水土流失状况等因子创建石漠化强度等级分类标准(表1)，最后利用ArcMap的空间分析功能对研究区的石漠化强度进行基于图斑识别的定级．

表1 石漠化强度等级分类标准Tab. 1 Classification standard of rock desertification intensity grade

3 结果与讨论

3.1 基于像元的后向散射系数分析

对于已构建的SAR时间序列，在影像中的石漠化地区选取训练样本，训练样本的选择应该遵循均匀地落在研究区内，且尽量能满足对整个研究区的代表性的选择原则．本研究在各示范区都选择了50个训练样本，均为6像元×6像元的正方形区域．下图为撒拉溪示范区部分训练样本在12时相(T1-T12)的后向散射系数值时间序列曲线图(见图2，篇幅有限不能全部列出)．

图2 部分训练样本后向散射系数值时间序列曲线图Fig. 2 Parts of training samples’backscatteringcoefficients time series trace

从数据上分析，撒拉溪示范区训练样本12个时相中的后向散射系数值大概集中在-24 db到-18 db之间，花江示范区训练样本12个时相中的数值集中在-27 db到-20 db之间，整体上数值大小范围为1月到6月数值较低，7月到9月数值较高，10月到12月数值再次呈现降低的趋势，其中，7月和8月数值最大，1月和12月达到最小值．究其原因，SAR后向散射系数值大小的主要影响因子主要有土壤水分、地表粗糙度以及植被覆盖等，而训练样本所在区域基岩裸露情况严重，植被覆盖较少，导致石漠化区域后向散射系数值一整年都比较平稳，变化趋势较小．花江示范区训练样本后向散射系数值较之撒拉溪示范区偏高一点的原因是该地区石漠化强度以中度和强度为主．而数值变化的原因在于数值高的月份太阳辐射较强，植被覆盖有所增强，土壤湿度也随之上升，后向散射系数值偏高；数值低的月份太阳辐射较弱，温度降低，植被稀疏，土壤水分冻结，导致土壤介电常数下降，后向散射系数值偏低．

3.2 SAR提取石漠化信息结果

在24时相的影像中目视选择质量最好的2016年7月13日的撒拉溪示范区雷达数据和2016年9月5日的花江示范区示范区雷达数据与选择的光学遥感数据进行影像融合，分别以51和42为DTW算法的阈值对撒拉溪示范区以及花江示范区进行石漠化信息提取，并对提取的结果进行二值化处理，结果见图3和图4.石漠化信息统计见表2.

表2 基于SAR的撒拉溪示范区和花江示范区石漠化信息提取结果Tab. 2 Information extraction of rocky desertification basedon SAR data in Salaxi and Huajiang demonstration areas

图3 基于SAR的撒拉溪示范区石漠化信息提取结果Fig. 3 Information extraction of rocky desertificationbased on SAR data in Salaxi demonstration area

图4 基于SAR的花江示范区石漠化信息提取结果Fig. 4 Information extraction of rocky desertification ofHuajiang demonstration area based on SAR data

3.3 数据融合的石漠化强度定级结果

基于SAR和光学遥感数据融合的研究区石漠化强度定级结果见图5和图6.

图5 基于SAR和光学遥感数据融合的撒拉溪示范区石漠化强度定级结果Fig. 5 Result of defining the levels of rock desertification intensitybased on SAR and optical data fusion in Salaxi demonstration areas

图6 基于SAR和光学遥感数据融合撒拉溪示范区石漠化强度定级结果Fig. 6 Result of defining the levels of rock desertification intensity ofHuajiang demonstration area based on SAR and optical data fusion

基于SAR和光学遥感数据融合的研究区各等级石漠化面积提取结果见表3.

表3 基于SAR和光学遥感数据融合的撒拉溪示范区各等级石漠化面积提取结果Tab. 3 Area extraction of different levels of rock desertification intensitybased on SAR and optical data fusion in Salaxi demonstration areas

3.4 DTW算法及提取结果精度验证

为了验证DTW算法的精度，本研究利用皮尔森相关系数算法和欧式距离算法对像元同样进行阈值分割处理，然后提取石漠化和非石漠化面积.为了验证提取结果的精度，本研究利用光学遥感数据对研究区石漠化强度进行等级确定，数据和定级方法与前文一致，分别选择2016年6月14日、10月4日的Landsat8数据，定级方法同样采用基于图斑识别的方法，定级结果见表4.

表4 基于光学遥感数据的撒拉溪示范区和花江示范区各等级石漠化面积提取结果Tab. 4 Area extraction of different levels of rock desertificationintensity based on optical data in Salaxi and Huajiang demonstration areas

对于DTW算法及提取结果精度的验证，本研究采用野外实地验证的方法．首先，在两个示范区分别选取200个精度验证训练样本．训练样本的选择遵循均匀分布且能尽量能代表整个研究区的原则．选取方式可以在地图上以目视的方式进行．然后利用GPS记录下整个野外实地验证的行程路径，用数码相机记录训练样本的照片并完成照片编号和该点石漠化等级的确定．最后，利用HOLUX ezTour软件，把行程路径和拍摄的照片文件转化为kml文件，再在谷歌地球中打开kml文件，作为DTW算法及提取结果精度验证的依据[19]．本研究选取石漠化界定正确率(该训练样本是否属于石漠化地物)和石漠化强度等级正确率两个指标作为精度验证的结果．结果见表5和表6.

从精度验证结果中分析，DTW算法在两个示范区的石漠化界定精度均要高于皮尔森相关系数算法和欧式距离算法；撒拉溪示范区基于SAR后向散射系数时间序列的提取方法的石漠化界定正确率为91%，比光学遥感数据提取方法的精度高7.69%；花江示范区基于SAR后向散射系数时间序列的提取方法的石漠化界定正确率为91%，比光学遥感数据提取方法的精度高9.41%；撒拉溪示范区SAR与光学遥感数据融合提取方法的石漠化界定正确率为88.5%，比光学遥感数据提取方法的精度高4.73%，石漠化强度等级界定正确率为83.5%，比光学遥感数据提取方法的精度高5.99%；花江示范区SAR与光学遥感数据融合提取方法的石漠化界定正确率为92.5%，比光学遥感数据提取方法的精度高8.82%，石漠化强度等级界定正确率为88%，比光学遥感数据提取方法的精度高5.39%.

表5 撒拉溪示范区和花江示范区不同算法的精度验证Tab. 5 Accuracy verification of different algorithms inSalaxi and Huajiang demonstration areas

表6 撒拉溪示范区和花江示范区石漠化信息提取结果精度验证Tab. 6 Accuracy verification of information extraction of rockydesertification in Salaxi and Huajiang demonstration areas

4 结论

本研究利用研究区多时相SAR数据集合，利用构建后向散射系数时间序列曲线的方式进行石漠化信息的提取，该方法使用DTW算法确定阈值从而对石漠化地物像元和非石漠化像元进行分割，利用SAR和光学遥感数据，结合石漠化强度等级分类标准表进行石漠化强度定级．实验结果表明：1)石漠化地物的后向散射系数值一年内趋于稳定，总体变化趋势较小，这也反证了利用SAR的后向散射系数时间序列提取石漠化信息的可行性与精准性；2)对比基于光学遥感数据的提取方法，基于SAR后向散射系数时间序列的提取方法在石漠化地物和非石漠化地物分割上精度更高；3)对比基于光学遥感数据的提取方法，基于多时相SAR与光学遥感数据融合的提取方法在石漠化强度定级上精度更高．

研究方法在一定程度上克服了使用单一时相SAR数据的缺陷，同时集雷达遥感和光学遥感之所长，兼顾雷达遥感成像中蕴含的丰富的纹理地形信息以及光学遥感成像中蕴含的丰富的光谱信息，使得石漠化信息提取精准度有所提高．