利用掩膜主动轮廓模型提取水库水面面积

苗立新,唐守正,李 霞,纪中奎,文 强,熊志明

（1.中国林业科学研究院资源信息研究所,北京100091;2.二十一世纪空间技术应用股份有限公司,北京 100096）

密云水库作为北京惟一的地表水饮用水源,承担着北京城区60%的自来水供应任务,对保障北京稳定与发展起着至关重要的作用。近年来,密云水库入库水量连年递减、水面逐年萎缩,对城市的用水安全造成影响,成为人们关注的焦点。水面面积是衡量水库储水量的一个重要指标[1]。遥感技术具备有宏观、实时及成本低等特点,从卫星遥感影像中快速、准确地提取水体信息已成为水资源调查、水资源宏观监测及湿地保护的重要手段。北京一号小卫星多光谱数据具有时间分辨率高的优点,可以满足对密云水库进行密集连续动态监测的需求。但是由于监测周期比较短,水库水面变化不明显,同时水面提取的精度受系统误差和人为误差的影响,水面的实际变化趋势不能很好地反映出来。为了最大限度地减少系统误差和人为误差,提高水库边界提取的精度,有必要寻找一种稳定的水体提取方法。

针对几种常见分辨率遥感影像的光谱特点,学者们提出了单波段阈值法[1-2]、谱间关系法[3-4]、色彩空间变换法等多种水体提取方法。例如,刘建波等[1]利用密度分割法从TM影像中提取水体的分布范围。陆家驹等[2]分别利用阀值法、色度判别法、比率测算法从TM影像中识别水体。杨存建等[3]和钟春棋等[4]分别以TM影像中TM2+TM3和TM4+TM5的差值及二者的比值辅以适合的阈值来提取水体。王刚等[5]提出两种水体提取方法:一种是将谱间关系法与IHS彩色空间构建模型相结合的方法,另一种是基于卫星图像数据的LBV变换与归一化植被指数NDVI的提取方法。都金康等[6]及邓劲松等[7]利用决策树方法分别从SPOT 4和SPOT 5影像中提取水体信息;赵书河等人[8]利用迭代混合分析法从中巴资源一号卫星数据中提取水体。韩栋等人[9]针对北京一号小卫星多光谱数据的特点,引入特征波段PRWI=(B3+B1)/(B3-B1),通过设置不同的PRWI阈值提取水体信息。在以上各种水体提取方法中,一般都存在一个必不可少的步骤,即阈值的选取。由于阈值的选取一般受人的主观影响比较大,不仅需要反复试验来确定,而且当阈值选择不合适时,提取的水体边界准确度会受到很大影响。同以上方法不同,本文提出的利用掩膜主动轮廓模型提取水体边界的方法,在水体提取过程中不需要人为选取任何阈值,而是充分利用影像自身内在的信息(梯度信息以及对象间的异质性信息),利用水平集演化的方法提取水库水面边界,为水库水面面积提取方法的探索提供了一种新的思路。

1 掩膜主动轮廓模型

1.1 主动轮廓模型

主动轮廓模型(Active Contour Model)最早由Kass[10]等人提出,用来解决图像边缘检测问题。主动轮廓模型的基本原理是在满足给定图像约束的情况下,为了检测出该图像中的对象,而让一条曲线沿着一定的规则进行变形和移动。这种为了恢复对象的形状而在数字图像内进行变形的曲线称为主动轮廓(Active Contour)。在Kass等人提出主动轮廓模型之后,国内外学者从多个方面对该模型进行了改进,形成多种形式的主动轮廓模型[10-19]。根据主动轮廓模型中停止项是否依赖于图像的梯度可以将主动轮廓模型分为基于边缘的主动轮廓模型和基于区域的主动轮廓模型。Ron Kimme[16]将边缘与区域信息结合起来进行目标提取,提出了一种结合了鲁棒对齐项(Robust Alignment Term)、加权面积(Weighted Area)、测地主动轮廓(Geodesic Active Contour)和最小方差(Minimal Variance)的主动轮廓模型。

1.2 掩膜主动轮廓模型

本文提出的掩膜主动轮廓模型,即一种带有掩膜的主动轮廓模型。它主要依据Ron Kimme[16]所提出的主动轮廓模型,并根据实际需要进行了4个方面的改进:

(1)Ron Kimme模型只能从灰度图像中提取目标,为了使模型能够处理彩色或多光谱图像,依据Chen等人的方法对最小方差项进行扩展[17-18],并对所有涉及图像梯度的项,采用Di Zenzo梯度[20]代替传统的灰度梯度,从而将模型的适用范围从灰度图像扩展到彩色或多光谱图像。

(2)根据Chunming Li[21]的不需要重新初始化的水平集思想,增加了一个惩罚项,使水平集函数接近符号距离函数。因此,可以允许采用比较大的时间步长。

(3)Ron Kimme模型及其它各种主动轮廓一般是在迭代开始前确定一个初始的时间步长,这个时间步长在各次迭代过程中保持不变。掩膜主动轮廓模型以一个逐渐收敛的可变时间步长代替固定时间步长,在迭代开始时可以将时间步长设置得比较大,随着迭代的进行,步长逐渐缩小,从而在加快收敛速度的同时,保证了边界提取的精度。

(4)在模型中加入了掩膜的概念。由于模型中最小方差项要求待分割的图像是近似双峰的,为了满足这个条件,引入了掩膜的概念,从而得到掩膜主动轮廓模型。掩膜主动轮廓模型需要对待分割图像和水平集函数同时做掩膜运算,即Im=I◦Mask,φm=φ◦Mask。初始掩膜取整个图像范围,在每次迭代过程中,取水平集内部区域(φ小于0的部分)为前景对象(水库水面),对它进行数学形态学膨胀,从而得到当前迭代中使用的掩膜,代入水平集函数的演化方程进行下一次迭代。随着迭代的进行,经过掩膜后的待分割图像越来越接近双峰图像,从而满足最小方差项对图像的要求。同时,引入掩膜后,由于每次迭代都有一部分明显不属于水库的像元不参加计算,从而计算的速度也大大加快。

改进后的掩膜主动轮廓模型可写成如下形式:

式中:φm——经过掩膜后的水平集函数;Im——掩膜后的待分割的图像,可以为灰度、彩色或多光谱图像;α,β,γ,μ,ν——均为常系数;div ——散度 ;Δ——拉普拉斯操作符;gm(x,y)——边缘指标函数(或称边缘停止函数),这里取

式中:n——迭代次数;tn——第n次迭代的时间步长。

为简洁起见,记

则掩膜主动轮廓模型的迭代方程可以简写为如下形式:

式中:GAC——测地主动轮廓(Geodesic Active Contour)项,它对应的泛函∫g(Cs)ds可以理解为一个加权了的弧长。其中,Cs为参数化曲线,gm(x,y)为边缘指标函数。当gm为常数1时该泛函表示曲线的弧长。作为对数据敏感的规则项,GAC用来控制主动轮廓线的光滑程度;AR——鲁棒对齐项,它的作用是使图像梯度方向与曲线法线方向尽可能一致,其中sign为符号函数,对于参数大于0、小于 0、等于0的情况,其返回值分别为1,-1,0;MV——最小方差项,它的作用是使区域内的方差尽可能小,其中cm1和cm2分别为闭合曲线(零水平集)内部和外部的图像像元亮度中值(向量);WA——加权面积项,其对应的能量泛函EW(C)=∫∫Ω Cf(x,y)dxdy是曲线围成的封闭图形内部区域的加权面积,其中f(x,y)可以是任何标量函数,当简单地取 f(x,y)=1时,加权面积项实际上变成了气球力(Balloon Force)[16,21];P——惩罚项,用来描述一个函数φ接近符号距离函数的程度。根据掩膜主动轮廓模型的迭代公式(1),采用IDL(Interactive Data Language)语言[22]编程实现了掩膜主动轮廓模型,用以支持后续的试验。

2 案例分析

2.1 试验数据

虽然当前可供选择的遥感数据源已经相当丰富,其中不乏分辨率在米级以及亚米级的数据,如RapidEye数据、Spot5 数据 、QuickBird 数据 、WorldView数据等。但高分辨率数据一般幅宽较窄,重复覆盖较大区域的周期比较长,并且多数卫星不属于国内自主运营,数据的可获取性较差,因此,难以满足以月度或更短时间周期对密云水库水面面积进行连续动态监测的要求。北京一号小卫星及运营系统,是国家“十五”科技攻关计划和高技术研究发展计划(863计划)联合支持的研究成果,卫星上装有4 m全色和32 m多光谱双传感器,其32 m多光谱数据包括绿波段(520～600 nm),红波段(630～690 nm),近红外波段(760～900 nm)3个波段,幅宽达600 km,重访周期3～5 d,能够满足对密云水库进行月度监测的需求。4 m全色数据成像幅宽24 km,虽然无法在每个月完全覆盖密云水库一次,但可选择完全覆盖密云水库的4 m全色数据作为同期多光谱数据提取结果的验证数据,验证多光谱数据上水库面积的提取精度。研究使用的遥感数据情况见表1。

表1 研究使用的遥感数据

2.2 模型参数设置

在模型迭代过程中,鲁棒对齐项实际上对应着图像中边缘或梯度所起到的作用;最小方差项对应着对象内部同质性大小;测地主动轮廓项对应着曲线的光滑程度;加权面积项或气球力起到的主要是加快曲线收缩或膨胀的作用。在这几个参数中,最重要的是控制边缘的鲁棒对齐项和控制对象内部同质性的最小方差项在所有项中所占的权重,在确定模型参数时,可根据具体图像的不同特征,给两者赋予适当的权重系数。

通过取不同组权重系数进行反复试验,发现模型中的鲁棒对齐项系数、最小方差项系数、测地主动轮廓系数、气球力系数和初始时间步长的参数值分别设置为5 000,100,100,100,0.01,对试验中的绝大多数图像可以取得比较好的提取结果,而且,除非各项系数的值发生非常大的变化,否则,系数的改变对结果的影响非常小,表现出非常强的稳定性。因此,取这组参数作为模型的输入参数,对包括上述6期多光谱影像在内的一系列多光谱影像进行水库水面提取。以2007年6月26日数据为例,全色影像上人机交互提取结果和多光谱影像上模型法提取结果对比见图1。

2.3 试验结果及分析

基于上述6期北京一号小卫星32 m多光谱影像数据,利用掩膜主动轮廓模型提取水库水面边界并统计水面面积。同时,采用人机交互解译方法对同期的北京一号小卫星4 m全色影像提取水库水面边界并统计水面面积,作为精度验证的参考数据。将掩膜主动轮廓模型法提取的面积值(简称模型法面积)与4 m全色影像提取水库水面面积(简称全色面积)对比,结果见表2。

图1 2007年6月26日北京一号全色影像人机交互提取结果(a)与多光谱影像模型法提取结果(b)对比

表2 模型法面积与全色面积对比

从表2中可以看出,掩膜主动轮廓模型从6期影像上提取结果的平均精度可达97%以上。说明采用模型法提取的水库水面面积具有较高的精度。通过计算精度的标准差可知,模型法精度的标准差为1.14%,说明模型法提取的水面面积比较稳定,波动较小,比较适合连续监测。另外,模型法面积与全色面积的差值全部为负值,说明模型中仍然存在着一定的系统误差,使模型的结果比参考结果稍有偏小,经过目视检查认为,该误差主要来源于水库边界处的混合像元,当需要进一步提高水面面积监测精度时,可以考虑采用混合像元分解方法对面积进行进一步修正。

将模型法面积与全色面积(验证值)生成散点图,并进行线性拟合,模型法面积与验证值非常接近,其线性拟合方程为y=1.166x-15.13,R2为0.896 1。

2.4 模型应用

经过验证,采用掩膜主动轮廓模型从多光谱数据上提取密云水库水面面积能够达到比较高的精度,满足业务应用的精度要求,同时,各期提取结果的稳定性较好。因此,采用该方法对2006年2月到2008年10月质量良好的北京一号多光谱影像进行处理,提取水库水面边界并统计面积,实现对密云水库的月度动态变化遥感监测,从而定量地了解其变化和发展趋势,为密云水库的科学管理、高效利用和有效保护提供理论依据和基础数据。2006年2月到2008年10月各月份密云水库水面面积变化曲线如图2所示。从图中可以看出密云水库水面面积从2006年到2008年对应月份总的变化趋势是逐年减少的。

图2 2006年 2月至 2008年10月密云水库各月份水面面积变化图

3 结论

采用北京一号小卫星多光谱数据等重访周期较短的遥感数据对密云水库水面面积进行连续密集监测,要求相邻时相的监测结果具有较高的稳定性和可比性。尤其是水体边界变动比较小,要求水体提取方法的误差相对于不同时相间水体变化必须足够小,才能真实地反映密云水库水体面积的变化情况。从试验结果可以看出,掩膜主动轮廓模型法能比较准确地提取出水库水面的边界,满足连续监测对精度的要求;各期数据提取结果与参考数据相比,精度变化波动小,具有很高的相关性,因此能够保证连续监测的稳定性。

掩膜主动轮廓模型提取水体,只需要在程序运行前设置好参数,提取过程中不需要人工干预,消减了人为误差的干扰。采用本方法进行连续动态监测,结果具有较高的稳定性和可比性,能够准确地反映水库水面面积变化的趋势,为水库水面面积自动获取和连续监测方法研究探索了一条新路。另外,由于遥感图像上不可避免地存在混合像元,会对提取面积产生一定的影响,如何采用混合像元分解方法进一步提高水面面积监测精度是后续研究需要解决的问题。

[1]刘建波,戴昌达.TM图像在大型水库库情监测管理中的应用[J].环境遥感,1996,11(1):53-58.

[2]陆家驹,李士鸿.TM资料水体识别技术的改进[J].环境遥感,1992,7(1):17-23.

[3]杨存建,徐美.遥感信息机理的水体提取方法的探讨[J].地理研究,1998,17(增刊):86-89.

[4]钟春棋,曾从盛,柳铮铮.基于谱间特征与比值型指数的水体影像识别分析[J].地球信息科学,2008,10(5):663-669.

[5]王刚,李小曼,田杰.几种TM影像的水体自动提取方法比较[J].测绘科学,2008,33(3):141-142.

[6]都金康,黄永胜,冯学智,等.卫星影像的水体提取方法及分类研究[J].遥感学报,2001,5(3):214-219.

[7]邓劲松,王珂,邓艳华,等.SPOT-5卫星影像中水体信息自动提取的一种有效方法[J].上海交通大学学报:农业科学版,2005,23(2):198-201.

[8]赵书河,冯学智,都金康.中巴资源一号卫星水体信息提取方法研究[J].南京大学学报:自然科学版,2003,39(1):106-112.

[9]韩栋,杨晓梅,纪凯.小卫星遥感影像自动提取水体方法研究[J].测绘科学,2008,33(1):51-54.

[10]Kass M,Witkin A,Terzopoulos D.Snake:Active contour models[J].International Journal of Computer Vision,1988,1(4):321-331.

[11]Cohen L D.On active contour models and balloons[J].CVGIP:Image Understanding,1991,53(2):211-218.

[12]Xu C,Prince J L.Snakes,shapes,and gradient vector flow[J].IEEE Trans.Imag.Proc.,1998,7(3):359-369.

[13]Caselles V,Catte F,Coll T,et al.A Geometric Model for Active Contours in Image Processing[J].Numerische Mathematik,1993,66(1):1-31.

[14]Caselles R V,Kimmel G.Sapiro,Geodesic active contours[J].International Journal of Computer Vision,1997,22(1):61-79.

[15]Chan T F,Vese L A.Active contours without edges,IEEE Trans[J].Image Processing,2001,10(2):266-277.

[16]Ron Kimme.Fast Edge Integration.Geometric Level Set Methods in Imaging[M].New York:Springer,2003:59-77.

[17]Chan T F,Sandberg B Y,Vese L A.Active contours without edges for vector-valued Images[J].Journal of Visual Communication and Image Representations,2000,11(2):130-141.

[18]Rousson M,Deriche R.A Variational Framework for Active and Adaptative Segmentation of Vector Valued Images[M].Geometric Level Set Methods in Imaging,Vision,and Graphics.New York:Springer,2003:195-205.

[19]Cohen L D.On active contour models and balloons[J].CVGIP:Image Understanding,1991,53(2):211-218.

[20]Silvano di Zenzo.A note on the gradient of a multi-image[J].Computer Vision,Graphics and Image Processing,1986,33(1):116-125.

[21]Li Chunming,Xu Chenyang,Gui Changfeng,et al.Level Set Evolution Without Reinitialization:A New Variational Formulation[C]//Proceedings of the 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition.Washington DC,USA:IEEE Computer Society,2005:430-436.

[22]IDL-Data Visualization Solutions.ITT Visual Information Solutions.[2009-12-25].http://www.ittvis.com/idl.