从发育过程来看, 泥石流主要发育于山坡或者山沟中, 发育中的泥石流表面植被覆盖差, 在图像中主要呈浅色调, 在SPOT5图像中一般表现为较亮的图斑, 与周围植被发育较好的区域有较明显的反差。因此在泥石流自动提取中, 去除植被是非常关键的一个步骤, 根据SPOT5图像和研究区的特征, 选择合适的植被指数进行去植被处理是非常重要的。

目前, 使用最多的植被指数有归一化差值植被指数(NDVI)和修正的土壤调节植被指数(modified soil adjustment vegetation index, MSAVI)。NDVI是Rouse等^[13]在1973年提出的, 是目前最常用的植被指数, 能够较有效地反映植被信息, 其针对SPOT5图像的计算公式为

NDVI= NIR-REDNIR+RED, (1)

式中: NIR为近红外波段的亮度值; RED为红波段的亮度值。

NDVI的计算结果如图1(a)所示, 图中白色和灰白色图斑代表植被发育的地区, 黑色和灰黑色图班代表没有植被或植被不发育的地区。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 SPOT5图像本底值指数计算结果Fig.1 Index calculation results of the background value of SPOT5 image

发育中的泥石流主要是裸露的土壤和植被稀少的地表, 在SPOT5图像中的光谱特征主要表现为较亮的图斑。为了将其与周围地物区分开来, 选用土壤亮度指数作为一个可选指数, 常用的土壤亮度指数有归一化差值土壤亮度指数(NDSI)和基于缨帽变换(tasseled cap transformation)的土壤亮度指数(brightness index, BI)等。本文经过反复研究发现, NDSI能够很好地反映泥石流在SPOT5图像中的光谱特征, 其计算公式为

NDSI= RED-GREENRED+GREEN, (2)

式中GREEN为绿波段的亮度值。

NDSI的计算结果如图1(b)所示, 图中白色和灰白色图斑代表裸露土壤和植被稀少土壤的地区, 其亮度值比较高。

由于多波段数据经常是高度相关的, 主成分变换通过坐标轴旋转使数据的方差达到最大, 生成互不相关的波段, 从而达到消除噪声和减少数据集维数的目的。SPOT5图像的4个波段从不同的光谱特征和空间特征反映了不同地物信息之间的区别与联系, 并且其主成分变换后的第一主成分包含了90%以上的地物信息, 因此选用主成分变换提取的第一主成分(PC1)作为一个遥感指数, 其能够较好地反映泥石流的遥感图像本底值。PC1的计算结果如图1(c)所示, 图中各类地物的边界更为清晰。

王一川等^[3]提出过一种基于DEM的沟谷中心线提取算法, 该算法从数学形态学的角度出发, 利用移动窗口(3× 3或者5× 5)在DEM上移动进行沟谷点的判断, 最终提取出沟谷中心线。图2为3× 3移动窗口示意图。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 3× 3移动窗口Fig.2 3× 3 moving window

其判断条件为:

1)若(g[i, j-1]-g[i, j])(g[i, j+1]-g[i, j])> 0, 且(g[i+1, j]-g[i, j])(g[i-1, j]-g[i, j])< 0,

当g[i, j-1] > g[i, j] 时, g[i, j]为山谷点;

当g[i, j-1] < g[i, j] 时, g[i, j]为山脊点。

2)若(g[i-1, j]-g[i, j])(g[i+1, j]-g[i, j])> 0, 且(g[i, j-1]-g[i, j])(g[i, j+1]-g[i, j])< 0,

当g[i-1, j] > g[i, j] 时, g[i, j]为山谷点;

当g[i-1, j] < g[i, j] 时, g[i, j]为山脊点;

当且仅当 (g[i-1, j]-g[i, j])(g[i+1, j]-g[i, j])=0 时, g[i, j]在此方向上是平地。

本文在此基础上, 针对DEM数据本身的误差和沟谷中心线提取结果的准确性做了2点改进:

1)由于DEM数据在采集和加工时难免带来的误差, 导致DEM数据中存在一些孤立点(和周围点的值相差较大的点), 因此在DEM的预处理过程中通过低通滤波来去除噪声点的干扰, 使DEM数据在一定范围内保持地形的连续性特征, 从而使提取的沟谷中心线尽可能少地出现短枝和离散点。

2)文献[3]中建议使用3× 3或者5× 5的小窗口进行移动匹配, 利用小窗口虽然计算速度较快, 但是本文通过试验发现, 窗口太小提取出来的沟谷中心线会出现较多的断裂现象。鉴于此, 本文将移动窗口大小设置为17× 17。另外, 由于本文是基于IDL语言实现泥石流自动提取, 而IDL语言对矩阵的执行效率相当高, 所以利用较大窗口来匹配DEM数据, 其计算速度并没有明显变化, 且沟谷中心线的提取效果也比较理想。

提取沟谷中心线的步骤包括: 首先对DEM数据做低通滤波去噪处理, 处理时直接利用ENVI软件系统二次开发所带的系统函数(CONV_DOIT), 函数参数METHOD设为3, 可以根据需要设置平滑核的大小; 然后利用一个17× 17的窗口在DEM数据中按照2个方向(即水平和垂直方向)进行移动匹配、判断符合条件的DEM数据值, 其判断条件是:

①水平方向, 若(g[i-1, j]-g[i, j]> 0且g[i+1, j]-g[i, j]> 0)且(g[i-2, j]-g[i, j]> 0且g[i+2, j]-g[i, j]> 0)且(g[i-3, j]-g[i, j]> 0且g[i+3, j]-g[i, j]> 0)且(g[i-4, j]-g[i, j]> 0且g[i+4, j]-g[i, j]> 0)且(g[i-5, j]-g[i, j]> 0且g[i+5, j]-g[i, j]> 0)且(g[i-6, j]-g[i, j]> 0且g[i+6, j]-g[i, j]> 0)且(g[i-7, j]-g[i, j]> 0且g[i+7, j]-g[i, j]> 0)且(g[i-8, j]-g[i, j]> 0且g[i+8, j]-g[i, j]> 0), 则g[i, j]为沟谷点。

②垂直方向, 若(g[i, j-1]-g[i, j]> 0且g[i, j+1]-g[i, j]> 0)且(g[i, j-2]-g[i, j]> 0且g[i, j+2]-g[i, j]> 0)且(g[i, j-3]-g[i, j]> 0且g[i, j+3]-g[i, j]> 0)且(g[i, j-4]-g[i, j]> 0且g[i, j+4]-g[i, j]> 0)且(g[i, j-5]-g[i, j]> 0且g[i, j+5]-g[i, j]> 0)且(g[i, j-6]-g[i, j]> 0且g[i, j+6]-g[i, j]> 0)且(g[i, j-7]-g[i, j]> 0且g[i, j+7]-g[i, j]> 0)且(g[i, j-8]-g[i, j]> 0且g[i, j+8]-g[i, j]> 0), 则g[i, j]为沟谷点。

需要说明的是, IDL中数组索引下标采用列行顺序。

尽管已经对DEM数据进行了上述预处理, 但是难免提取出的沟谷中心线仍然存在一些短枝和离散点(图3(a))。因此, 需要对提取结果中残存的短枝和离散点进一步剔除, 并对提取结果进行数学形态学闭运算, 以提高沟谷中心线的连续性。先膨胀后腐蚀的过程称为“ 闭运算” (先通过膨胀运算将不连续的沟谷中心线连接起来, 再通过腐蚀运算将因膨胀而变“ 肥大” 的沟谷中心线“ 瘦身” 而恢复原貌), 得到较为准确、连续的沟谷中心线(图3(b)); 最后对连续的沟谷中心线作形态学膨胀运算(膨胀核大小为5× 5)后, 得到沟谷范围(图3(c))。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 沟谷中心线与沟谷范围提取结果Fig.3 Extraction result for central lines of valley and valley ranges

将上述计算的NDVI, NDSI和PC1这3个遥感指数灰度图像合成, 并再次进行主成分变换提取新的第一主成分; 对所得到的灰度图像需要选择合适的阈值提取非植被信息和裸土信息。研究发现, 经上述计算后得到的第一主成分图像的直方图存在多个波峰和波谷, 为此, 本文对Tsai^[14]提出的基于双峰直方图的阈值选取算法进行了改进。改进后的算法针对新的第一主成分图像具有植被的灰度值较低、裸地的灰度值较高的特性, 寻找多峰直方图中的多个波谷, 并计算出对应于多个波谷的多个阈值。

由于计算后得到的新的第一主成分图像的灰度值为负数和小数, 为了计算方便, 首先要对图像的直方图做预处理。在实验区图像的直方图中, 灰度值在0和255处都有较多像元。根据波峰波谷的定义, 灰度值在0和255两个端点的直方图的波峰不存在, 而实际上这2个值恰恰对应于植被和裸地。为此, 在原始灰度直方图头尾各增加一个值为0的像元, 使直方图两端多了2个波峰或波谷。这样既不违反波峰的定义, 又可以把灰度值0和255处当作波峰处理。经预处理后, 就可以利用阈值自动选取算法对直方图进行迭代滤波和阈值计算, 实现非植被信息和裸地信息的自动提取。

本文基于ENVI+IDL软件平台开发的泥石流自动提取系统的技术流程如图4所示。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 泥石流自动提取流程Fig.4 Flow chart of automatic extraction of debris flow

包括如下4个重要部分:

1)泥石流信息遥感图像本底值计算(图4①)。通过对植被指数(NDVI)、土壤亮度指数(NDSI)和遥感图像经主成分变换得到的第一主成分(PC1)图像的集成计算和新一轮主成分变换, 获取新的第一主成分图像。该灰度图像中, 泥石流与周围地物的光谱值差异被增大, 便于利用阈值自动选择算法准确有效地提取泥石流候选区域。

2)阈值自动选取(图4②)。对获取的新的第一主成分图像, 利用改进后的基于多峰直方图的阈值自动选取算法进行泥石流阈值自动获取, 并根据获取的阈值提取出泥石流候选区域。由于新的灰度图像中存在部分不纯净的像元, 在提取的泥石流候选区域的二值图像中存在离散点和空洞点, 因此要利用形态学滤波的闭运算(先膨胀后腐蚀)对离散点和空洞点进行剔除。

3)沟谷范围提取(图4③)。首先, 基于DEM数据和改进的沟谷中心线算法提取沟谷中心线。由于地形总会存在一些洼地、鞍点和局部的一些小平原^[3], 且DEM数据在采集和处理过程中会带来各种误差, 导致提取的沟谷中心线出现短枝和离散点, 因此在提取沟谷中心线之前要对DEM数据进行平滑处理, 在沟谷中心线提取之后还要对沟谷中心线提取结果进行形态学滤波处理, 以消除短枝和离散点。然后, 对消除短枝和离散点后的沟谷中心线再作形态学闭运算, 以获得连续性较好的沟谷中心线(图3(b))。最后, 利用形态学滤波对已获取的连续沟谷中心线进行一定核大小(5× 5)的膨胀运算, 得到沟谷范围(图3(c))。经上述处理后的沟谷范围比较符合沟谷的连续性特征, 便于与泥石流候选区进行匹配。

4)泥石流信息筛选(图4④)。基于泥石流候选区的二值图像和沟谷范围二值图像做图像匹配处理, 得到疑似泥石流栅格图斑, 并将其矢量化; 然后以泥石流形成的条件和泥石流的空间特征为筛选条件对疑似泥石流矢量图斑进行面积、坡度和顺坡性筛选, 最终得到泥石流矢量图斑。

实验区位于福建省永泰县同安镇地区。采用2004年12月12日获取的SPOT5多波段图像数据(空间分辨率为10 m)。实验中假彩色合成图像采用B2(R), B1(G), B3(B)的波段组合(图5(a))。该区地势山高坡陡, 残坡积土体厚宽大, 土层结构松散, 又由于常年降雨丰富, 地表人为活动剧烈, 该地区成为泥石流多发区之一。

根据上述泥石流自动提取流程, 初步提取的疑似泥石流栅格结果如图5(b)所示, 其中有一部分误提的图斑和一些噪声点, 因此还要将其矢量化以后根据面积、坡度和顺坡性进行筛选, 最后提取的泥石流结果与原始的假彩色合成图像(图5(a))的叠加显示见图5(c)(图中蓝色图斑为提取的泥石流信息)。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 泥石流提取结果Fig.5 Extraction result of debris flow

通过对泥石流的自动提取结果和人工目视解译结果的对比分析与验证, 说明本文方法提取的结果具有较高的准确度。在实验结果的分析验证中, 为了更加直观准确地进行人工目视解译并保证解译结果的准确性, 将提取的泥石流结果与该地区三维透视立体影像图叠加(局部效果如图6所示, 红色圆圈内为提取的泥石流), 根据经验从泥石流的形成条件和植被发育程度、坡度坡向和顺坡性等泥石流解译标志进行人工目视解译和对比分析。结果表明, 本文方法提取的泥石流信息较为准确合理, 与基于野外经验的人工目视解译的结果较为吻合。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 泥石流提取结果与三维透视立体影像图叠加 (观察方位角为60° , 高度角为30° ; 红色圆圈内为提取的泥石流)Fig.6 Superposition of debris flow extraction result and 3D perspective image