冯光胜
(中铁第四勘察设计院集团有限公司,湖北武汉 430063)
崩塌是最重要的地质灾害之一,也是工程地质不良地质体之一,是工程勘察重点调查对象。崩塌多发生在陡峻的山坡上,具有突发性,常随着暴雨、地震及人为破坏等发生,对工程建设,人类生产、居住等造成重大危害和影响。随着遥感技术的快速发展,基于遥感及其集成技术对崩塌、滑坡及泥石流等地质灾害的调查取得了一系列成果[1-8]。崩塌的遥感识别、调查虽然已有诸多成功实例,但是通过大量文献对比研究发现,目前的研究方法主要集中在图像增强结合人工经验目视解译,而通过算法从影像上自动提取崩塌的研究鲜有报道,且未实现阈值自动选取等核心算法。
“工程地质遥感解译信息系统”软件针对崩塌信息遥感识别、解译中存在的问题,以华南山区为研究区,在文献[9]等研究的基础上,通过改进去植被算法、阈值自动选取算法及数据处理流程等,实现了崩塌信息的遥感自动提取,且取得了较好的实验效果。
崩塌发育于比较陡峭的斜坡上,发展中的崩塌面、崩塌体植被覆盖差,多呈浅色调,在SPOT-5影像中一般为白至浅灰色调。而在植被发育的华南山区,崩塌面、崩塌体及堆积体与周围环境差异较大。因此,软件在崩塌的自动提取算法中,针对SPOT影像及研究区的特征,选择合适的去植被计算算法是非常关键的一步。
目前,可用于植被计算的指数有很多,常用的有比值植被指数(RVI)、归一化植被指数(NDVI)、修正的土壤调节植被指数(MSAVI)等。通过文献及实验比较研究,在研究中选择了 NDVI和MSAVI进行去植被实验。
NDVI是Rouse(1973)提出的目前最为常用的植被指数,能够比较有效地提取植被信息,其计算公式为
式中,NDVI为归一化植被指数,ρNIR为近红外的反射率,ρRED为红光的反射率(下同)。
NDVI植被指数在大多数情况下可用来进行植被信息的提取,但是NDVI对土壤背景的变化比较敏感[11],为了消除土壤背景的影响,Huete[12]提出 SAVI,QI J.[13]等在此基础上提出了 MSAVI。MSAVI 能够比较有效地消除土壤背景的干扰,其计算公式为
大量的实验研究发现,由于依据NDVI和MSAVI计算的阈值图像均不能有效地增大植被和非植被间的差异,不利于用阈值自动选取算法实现自动提取。因此,在反复试验的基础上,对 MSAVI进行了修改,(2ρNIR+1)取值100时,可将植被信息调整为较大的正值,非植被信息调整为较小的负值,常量化后的计算公式为
式(3)中,常量100是个经验值,该值范围在1~255,如果接近1或255时容易受图像大小影响,造成提取结果存在误差,故软件中选择100作为常量。
利用 CMSAVI(Constant modified soil-adjusted vegetation index)进行植被信息的提取,既有效地减少土壤背景的影响,又增大了植被与非植被间的差异。反之,在崩塌信息的提取中,要消除植被的干扰,突出裸土信息。因此,基于 CMSAVI指数法去植被更为有效。
根据CMSAVI计算的阈值图像需要指定合适的阈值来提取非植被信息,目前阈值的提取方法多以人工统计方式获取,自动化程度不够,极大地影响了效率和准确性。经过研究发现,计算的阈值图像存在多个波峰谷,为此,对 D.M.Tsai[14]提出的基于多峰直方图的阈值选取算法做了改进。改进算法针对CMSAVI计算的灰度图像植被的灰度值较高、非植被地物的灰度值较低的特征,寻找多峰直方图中的多个波谷,并计算出对应于多个波谷的多个阈值。由于非植被地物的灰度值较低,所以选取最小的阈值作为提取非植被信息的阈值。
在算法设计中,由于CMSAVI计算的灰度图的阈值存在负值和小数,为了计算的方便首先要对图像的直方图做预处理。具体处理如下:在实验区图像直方图中,灰度值0和255处都有较多像素,但根据波峰波谷的定义,灰度值0和255作为直方图的两个端点不存在波峰,可是灰度0和255处的像素对应于原始图像中的植被和裸地等非常重要的目标,理应作为波峰处理。为解决这个问题,软件在原始灰度直方图头尾各加一个0值,这样既可以把灰度值0和255处当作波峰处理,又不违反波峰的定义。计算过程如图1所示,图1(a)是经CMSAVI计算的灰度图像的原始直方图,存在多个波峰,不利于自动获取阈值;图1(b)是对灰度图像进行32次迭代滤波处理后灰度图像的直方图,处理后的直方图有4个波峰,其3个波谷的阈值分别为31,131和212。在预处理之后,就可以利用阈值自动选取算法进行阈值计算,实现非植被信息的自动提取。
图1 灰度图像的直方图预处理
基于ENVI+IDL二次开发的“工程地质遥感解译信息系统”软件中崩塌自动提取的流程如图2所示,该过程包括7个计算步骤。
发展中的崩塌面、崩塌体植被覆盖差,土体或岩体充分暴露,使得CMSAVI值降低,且明显低于周围植被发育区域,故可利用改进的多峰直方图法进行阈值的自动选取,从而有效地将非植被信息从影像中分离出来。
在利用CMSAVI指数法去植被的计算过程中,由于地形阴影的影响,提取的非植被信息中,仍存在部分阴影中的植被信息被误提出来,因此需利用近红外波段对阴影比较敏感的特征[9]将阴影中植被信息去除,阴影反射率的计算亦采用多峰直方图阈值自动选取算法实现。
在山地丘陵地区,崩塌发生在55°~70°斜坡上的几率较大,故通过坡度筛选计算,可将低坡度上的建筑物、耕地等非崩塌信息图斑去除。
在提取的目标图斑中存在大量的空洞和不连续,其原因在于裸地并不是非常纯净,其中夹杂少许未破坏的植被,因此需要利用形态学膨胀滤波算法进行填充,以还原崩塌面较为连续的实际情况,但是膨胀滤波后存在图斑扩大的问题,需要进一步利用腐蚀滤波进行处理,从而获得比较理想的目标信息。
为了实现后续的筛选计算,需要将前面计算的栅格图斑转为矢量图斑,并对图斑边缘进行追踪形成新的图斑,以消除栅格数据转换后的锯齿现象。
根据崩塌在SPOT5影像中可识别、提取的尺度分析[9,15],在 10 m 分辨率的影像上,面积小于 600 m2的小崩塌体无法准确提取出来,而面积大于50 000 m2崩塌较少,或是由伐林残留区等交错粘连形成的非崩塌信息,可以删除。
由于崩塌发生过程中是受重力控制的,虽然地形有一定的影响,但是崩塌总体上具有沿最大重力梯度下降的特性,利用顺坡性算法可有效去除坡地、道路等人类活动信息。崩塌自动提取的流程如图2所示。
软件实验区选为向莆铁路的永泰县地区,采用SPOT-5影像,多光谱分辨率为10 m,全色波段为2.5 m,该影像获取于2009年12月12日。该实验区是福建省崩塌等灾害的多发区之一,残坡积土体厚度大,土层结构松散,地表人为活动强烈,暴雨后不同强度的降雨入渗是诱发该区出现不同规模崩滑体的主要因素[16]。
图2 崩塌自动提取的流程
根据崩塌自动提取流程,首先进行去植被计算。图3是MSAVI和CMSAVI去植被结果的比较。图3(a)是实验地区SPOT合成假彩色影像截图(400像素×400像素);图3(b)是基于CMSAVI指数计算的阈值图像,该阈值图像与NDVI指数计算的结果比较相似;图3(c)是基于CMSAVI指数计算的阈值图像;图3(d)是利用NDVI和阈值自动选取算法计算的去植被结果;图3(e)是对图3(b)进行去植被自动计算的结果,图3(f)是对图3(c)进行去植被自动计算的结果。
通过实验比较分析,基于NDVI、MSAVI植被指数法提取的裸地图斑数较多,而基于CMSAVI去植被,能更为准确有效提取非植被信息,误提图斑极大减少,更利于用阈值自动选取算法进行处理以自动进行崩塌信息的提取。
图3 去植被算法效果比较
图4 崩塌提取结果比较验证
进一步的实验是在利用CMSAVI去植被获取非植被信息后,依次对提取的结果进行阴影区植被信息去除、坡度筛选、形态学滤波、栅矢转换及新矢量图斑生成、面积筛选及去规则图斑筛选、顺坡性计算等处理,最终提取出研究区的崩塌信息,如图4所示。图4(a)是基于NDVI去植被算法提取的崩塌结果与影像叠加显示,共提取崩塌信息19处;图4(b)是基于MSAVI去植被法提取的崩塌结果与影像叠加显示,共提取崩塌信息11处;图4(c)是基于CMSAVI去植被法提取的崩塌结果与影像叠加显示,共提取崩塌信息8处。
提取结果通过和2.5 m的全色波段影像的人工解译及现场验证资料进行对比,研究区人工解译的崩塌图斑有16处,现场验证为11处,其中面积大于600 m2的有10处。软件自动提取崩塌8处,与现场验证资料相吻合,提取精度约73%,如图5所示。经过对比分析,该方法误提图斑少,提取的精度较高,达70%以上。
针对崩塌体进行大量、反复的实验研究及现场验证,结果表明利用改进的CMSAVI算法去植被后提取崩塌效果更为理想。在崩塌提取过程中,对阈值自动选取算法进行了反复测试,基于该算法成功实现了崩塌体信息的自动提取,为铁路工程地质选线提供快速、准确的地质资料,有效指导现场勘察工作。
本文研究的不足之处在于:①CMSAVI算法改进过程中是基于实验经验的,对其机理的研究还有待深入;②形态学滤波的滤波算子大小是根据实验设定的,目前还无法根据图像自行确定。
图5 CMSAVI提取结果与现场验证资料对比
[1]杨则东,鹿献章.安徽境内长江岸带崩塌遥感调查[J].国土资源遥感,1998(3):22-25
[2]杨则东,徐小磊,谷丰.巢湖湖岸崩塌及淤积现状遥感分析[J].国土资源遥感,1999(4):1-7
[3]高克昌,赵纯勇.基于TM影像的万州主城区崩塌地质灾害研究[J].遥感技术与应用,2003,18(2):91-94
[4]吴泉源,姜春玲,邹敏,等.GIS技术支持下的招远市金矿区崩塌遥感调查[J].山地学报,2006,24(6):672-677
[5]陈文平,韩小明,范英霞.SPOT5影像数据在伊犁谷地地质灾害遥感调查中的应用[J].新疆地质,2008,26(4):396-398
[6]张景华,张建龙.Quickird卫星影像在汉源县小堡乡崩塌群识别中的应用[J].沉积与特提斯地质,2009,29(2):104-106
[7]邵虹波,唐川,李为乐.滇西某铁路拟选线路地质灾害遥感解译特征研究[J].防灾减灾工程学报,2009,29(3):342-346
[8]许冲,戴福初,陈剑,等.汶川Ms8.0地震重灾区次生地质灾害遥感精细解译[J].遥感学报,2009,13(4):754-762
[9]姚鑫,张永双,王献礼,等.基于地貌特征的浅层崩滑体遥感自动识别[J].地质通报,2008,27(11):1870-1874
[10]Rouse,J.W.,R.H.Haas,J.A.Schell,et al.Monitoring Vegetation Systems in the Great Plains with ERTS[J].Third ERTS Symposium,1973,NASA SP-351 I:309-317
[11]陈述彭,赵英时.遥感地学分析[M].北京:测绘出版社,1990
[12]Huete A R.A soil-adjusted vegetation index(SAV I)[J].Remote Sensing of Environment,1988,25(3):295-309
[13]QI J,CHEHBOUNI A,HUETE A R,et al.A modified soil adjusted vegetation index[J].Remote Sensing of the Environment,1994,48:119-126
[14]D.M.Tsai.A fast thresholding selection procedure for multimodal and unimodal histograms[J].Pattern Recognition Letters,1995,16:653-666
[15]仇大海,蒋炜,牛海波,等.遥感影像分辨率分析技术在滑坡研究中的应用[J].地质灾害与环境保护,2010,21(1):105-108
[16]袁东,池永翔,程刚.闽北地区不同植被类型下滑坡体土层入渗性能研究[J].长江科学院院报,2010,27(5):8-12