河流是生态系统的重要组成部分,也是人类生存所需的重要物质资源^[1]。河流信息的准确获取对水土资源管理、灾害预测以及环境评估极具价值^[2]。针对不同数据源影像国内外学者提出了许多水体信息提取方法,其中包括光谱特征阈值法^[3,4]、水体指数法^[5,6]、监督/非监督分类法^[7,8,9]等。然而上述方法多针对面状水体的提取,对于细小水体,其提取效果并不理想。细小水体主要指受遥感影像空间分辨率影响,水体总体宽度较窄、深度较浅,在遥感影像中呈现细长线状空间特征的地物^[10]。该类水体其光谱信息与周围环境较为相似,导致水体信息提取时不可避免存在河流间断现象的出现。

目前,针对细小河流提取已有许多研究,Geling等^[11]对Duda算子(Duda road operator,DRO)进行改进,提出用于线边缘检测的乘性DRO算子(multiply DRO,MDRO),该算子具有良好的暗、亮线边缘检测能力; 为避免MDRO算子在阶梯状边缘产生响应的缺陷,姜浩等^[12]对其算子进行改进,提出了线状特征增强算子(linear feature enhancement,LFE),较好地突出了遥感影像中线状地物信息,提高了算法对蜿蜒曲折自然河流的较好适应性; 杨树文等^[13]利用TM影像的多波段谱间关系提出一种山区细小水体提取方法,有效地去除了阴影噪声干扰; 周艺等^[14]为克服城镇建筑用地与细小水体混淆的缺陷,提出了伪归一化水体指数(false normalized difference water index,FNDWI),实现了城镇周边细小河流的提取; 李艳华等^[15]针对山区线状细小河流,使用面向对象的方法实现了精确化提取。此类结合细小河流光谱特性及其线状特征的河流识别算法,较准确地提取了遥感影像中的细小河流。但由于河流走向复杂,干扰较多,该类算法无法避免河流提取不连续的现象。

本文针对TM遥感影像中细小河流提取不连续问题,提出一种利用水体提取模型^[16](water extracting model,WEM)结合结构相似(structural similarity,SSIM)区域搜索的河流提取方法。该方法充分运用了河流像元间的高度相关性,构建启发式河流提取模型,发展面向细小河流的遥感水体提取方法,较好地解决了传统方法以地物特征差异为主的不足,实现复杂环境下的细小河流完整提取。

传统细小河流提取算法由于受混合像元效应影响,当水体与背景像元值难以区分时,所提取的河流往往会发生不连续现象,不能保证其提取完整性。本文基于传统算法改进,提出基于SSIM区域搜索的细小河流提取算法,具体流程为: 首先,利用TM影像中水体与其他地物之间亮度值差异,使用WEM区分水体与其他无关地物; 在此基础上,结合TM各波段地物间的光谱特征差异以及剩余噪声的面积大小,进行噪声去除; 最终,在去噪后的不连续细小河流结果上,对待连接区域的河流断点进行评估选取,并基于SSIM区域搜索方法对细小河流进行连接。完整的河流提取流程如图1所示。

TM影像中,水体光谱特征在第2—5波段上与其他地物具有明显区别^[17]。第2波段为绿光波段,对水体具有一定的穿透力,可反映水下特征,如图2(a)所示; 第3波段为红光波段,对水中悬浮泥沙以及裸露地表反映敏感,可为地物分类提供丰富信息,如图2(b)所示; 第4波段为近红外波段,该波段处于水体强吸收区,对突出水体信息独具优势,如图2(c)所示; 第5波段为中红外波段,该波段中的水体反射率继续下降,与无关背景信息间差距增大,可大大抑制背景噪声,如图2(d)所示。

对这4个波段组合运算,可抑制水体外的其他无关信息,突出显示水体区域。本文采用WEM进行水体信息提取,即

式中TM₂,TM₃ ,TM₄和TM₅分别对应TM影像第2—5波段的像素值。经该模型所提取的影像,水体部分较为明亮,非河流地物较为灰暗。

植被、河流以及高山峡谷等受光照、成像过程和预滤波等影响,在水体提取过程中会产生大量噪声。这些噪声主要体现在2个方面: 山体阴影和线型噪声。

文献[18]认为,在TM影像绿光波段中噪声与水体的灰度值差异最大。可将绿光波段灰度值小于T_shadow1的像素点标记为噪声。由于TM影像近红外波段中水体具有较强吸收能力^[19],因此,近红外波段灰度值大于T_shadow2的位置也可标记为噪声。即在水体信息图W上,将TM₂和TM₄满足式(2)的位置点(i,j)的W值设为0,即

对于线型噪声,其连通域面积较小,通过设定面积阈值便可删除这些噪声。满足式(3)设定的区域被认为是噪声,即

式中: S_i为第i个连通区域的面积; N为连通区域的总数目; T_α为面积阈值。

TM影像中定义细小河流宽度为1~3个像素。在沟谷较深、地物光谱与河流相近的地方,常常存在水体像元值与周围地物像元值非常接近的情况,使得经水体信息提取、噪声去除后的局部河流出现不连续现象,造成河流结果呈现残缺不完整状态。因此,为连接河流不连续区域,需首先确定河流待连接端点位置,以及待连接端点之间的连接关系。设p(r,c)为河流图中第r行第c列的端点,当列号和行号已知或可以忽略时,将p(r,c)简写为p。因此,对于端点p,要求点集合E,即

有|E|=2成立。

由所有的端点p组成端点集合P,即

对于端点集P中所有点,需要评估其是否为河流待连接点,同时需要找到其对应待连接点。二者组成一个待连接点对,即在所有河流待连接点之间建立一个一一对应的映射。本文用来评估待连接点p和q之间的连接属性的函数为

式中: D(·)为距离函数; S(·)为相似性函数。

距离函数采用欧式距离,p(r,c)和q(k,l)2点的欧式距离为

相似性函数采用SSIM^[20]度量,即

式中: x和y为2个信号序列; μ_x,μ_y和 σx2, σy2分别为x和y的均值和方差; σ_xy为x和y的协方差; c₁和c₂分别为用来维持函数稳定的常数,该算法中取c₁=0.01,c₂=0.001。

SSIM常用于图像质量评价,其取值范围为[0,1],值越大表明图像失真越小,图像相似度越高。为提高间断河流的连接准确性,本文将SSIM应用于待连接像元点的评估,当2个河流像元点靠得越近,其在亮度、对比度和结构度上差异越小,即SSIM越大。算法中针对不同连通域的2个待连接河流端点p和q,选择可组成符合人类视觉习惯的假彩色影像的TM影像的3、4、5波段,分别以p和q点为中心,根据其周围3像元×3像元区域的影像块,组成27维向量,代入式(8)进行计算得到其平均SSIM度量值。

连接属性值L越小,表明2个断点之间的欧氏距离越小,SSIM越大,属于连接点对的可能性越大。根据式(6),计算P中除自身之外的所有点与p的相似度L。p的待连接对象q'应满足

对于连接点对p和q',其连接结果应尽可能地与原图中河流位置相一致。借鉴启发式搜索算法的思想,构建河流连接点对之间搜索算法。估价函数F(n)一般定义为

式中: G(n)为代价函数,即从当前节点n到初始节点p(或目标节点q')所需付出的代价; H(n)为启发函数,即从节点n到目标节点q'所需付出的预估代价; λ为其所占权重。

TM影像中,河面光谱信息较为一致,当进行河流连接时,每一步搜索先计算3—5波段中当前点8领域与起点、终点的SSIM,再选择距离最近的点与当前点计算SSIM,每次选择相似度最大值位置进行连接。即SSIM取决于当前节点n_i的位置,因此G(n)可采用SSIM来构造,即

启发函数H(n)可直接使用距离函数。

搜索过程如图3(a)所示。图中p为起点,q'为终点,根据式(11)计算p的8邻域的估价函数值,如图3(b)。根据计算,此时p点下方的估价函数值最小,将此点选定为最佳候选点。移动到最佳候选点,并对新的点重新计算G(n)和H(n),再次计算估价函数值,直至q'为最佳候选点时停止搜索。

采用Landsat TM/ETM+影像进行算法验证,为保持影像的原始光谱特征,所有影像均未做任何预处理操作,文中假彩色影像均采用B5(R),B4(G),B3(B)组合方案合成。为验证算法的精确度、稳定性和鲁棒性,分别针对不同区域设计了不同试验、试验区共有4个: 云南省德宏傣族景颇族自治州怒江试验区; 云南省丽江市玉龙纳西族自治县金沙江试验区; 四川省凉山彝族自治州木里藏族自治县理塘河试验区; 云南省怒江傈僳族自治州独龙江试验区。

2.1.1 阴影噪声去除

算法中所涉及到的阴影噪声去除,与绿光波段TM₂和近红外波段TM₄的阈值设定紧密相关,不同阈值的设定,影响着算法最终提取结果的呈现。由于TM₂与TM₄可互补地去除不同区域阴影,因此利用式(2),可将TM₂中小于T_shadow1或者TM₄中大于T_shadow2的像素点当做阴影噪声去除。为进一步探索参数对细小河流提取结果的影响,以云南省德宏傣族景颇族自治州怒江局部区域为试验区(图4),该区多高山陡坡,地域环境复杂。在式(2)的阴影去除步骤中进行如下对比试验: ①只考虑TM₂分量的阈值T_shadow1; ②只考虑TM₄分量的阈值T_shadow2; ③同时考虑TM₂和TM₄2个分量的阈值。

试验①只考虑T_shadow1,其噪声去除效果明显较差,阈值设置过小,噪声几乎无法去除; 阈值设置过大,河流信息丢失严重。试验②只考虑T_shadow2,其阈值设定规则与试验①完全相反,即阈值越小河流信息丢失越严重,并且噪声去除效果与试验①呈现明显互补趋势。试验③设定T_shadow1=0.1,T_shadow2=0.2时,细小河流信息几乎被完整保留,但剩余噪声面积较大,需设置较大阈值T_α才能去除剩余噪声,但这将导致不连续河流片段大面积丢失。类似地,当T_shadow1=0.11,T_shadow2=0.2时,河流严重不连续现象,以至于影响后续细小河流启发式搜索方向的确定。当T_shadow1=0.1,T_shadow2=0.17和T_shadow1=0.11,T_shadow2=0.17时,河流信息虽有丢失,但细小河流整体轮廓存在,且剩余噪声面积较小,便于后续去噪时较完整地保留细小河流信息。其余阈值去噪所呈现的结果中,河流信息丢失严重,因此不予考虑进一步去噪处理。

2.1.2 小面积噪声去除

小面积噪声去除过程中,阈值T_α的大小关乎着噪声去除效果和不连续河流保留问题。当阈值T_α设置过大,会在去除噪声的同时,将不连续河流信息误去除,造成后续搜索困难; 当T_α设置过小,则会留下太多的噪声,一定程度上影响待连接河流断点的判断。因此,继续以阴影噪声去除后的怒江局部区域为试验区展开对比分析(图5)。

本试验中,当阴影去除阈值最佳时(T_shadow1=0.1,T_shadow2=0.17),阈值设定小面积噪声阈值T_α=3为理想情况,此时去噪效果和河流保留问题可得到较好平衡。此外,T_α=2和T_α=4,T_α=5分别代表了阈值过小和阈值过大的情况,噪声去除效果均不理想。

2.1.3 启发函数权重

无关噪声去除后(T_shadow1=0.1,T_shadow2=0.17, T_α=3),仍需结合SSIM区域搜索算法进行不连续河流连接,此时,相关启发式函数的权重λ发挥着决定性作用。因此,在小面积噪声去除的基础上仍以怒江局部区域为试验区,分析启发式权重的取值(图6)。

本文试验的启发式函数权重最优取值区间为[1.2,2.5],该实验中最佳权重为λ=2,此时,河流连接正确率较高。当λ取值过小时,河流结果受混合像元效应影响出现过搜索现象; λ=1时河流以相似性最大方向作为主搜索方向,导致细小河流部分结果区域掺杂着些许无关块状信息; 若λ取值过大,河流结果易呈现出偏直线化不自然结果; λ=5和λ=10时整个搜索过程中距离函数作为主导,导致河流搜索方向不受相似性函数引导,出现“畸形”河流结果。

2.2.1 本文方法提取结果与分析

混合像元效应是导致细小河流难以准确提取的重要原因,选取混合像元效应明显且地形复杂的云南省丽江市玉龙纳西族自治县金沙江地区的4个局部区域作为研究对象,利用本文方法提取细小河流(T_shadow1=0.1,T_shadow2=0.09,T_α=0,λ=3),以验证本方法的提取准确性,主要步骤结果见图7。

在金沙江4个局部区域的TM5(R),4(G),3(B)假彩色合成影像中,图7(a)—(b)河流流经区域峡谷较多,部分河流信息受其影响,不易与其他地物较好区分; 图7(c)—(d)地形相对统一,受其影响的河流区域较少。经WEM处理后的影像中,细小河流与峡谷相较于其他地物亮度明显提高,河流信息基本完全保留,但山体阴影仍然存在。噪声去除后的细小河流,受山体阴影及混合像元效应影响,许多河流弯道几乎未被提出,局部信息丢失严重,其中,图7(i)中共发生6处断裂,图7(j)中有9处区域发生断裂,图7(k)中共丢失13段河流信息,图7(l)中有7处断裂。经SSIM区域搜索方法所提取的细小河流结果中,河流完整且连续。与图7(i)—(l)对比,4个区域中即使河流局部被山体阴影遮挡,发生多处断裂,通过SSIM区域搜索,仍能准确地将其进行连接,使其与原河流轨迹一致。

2.2.2 细小河流局部对比

为了验证本文方法的普适性以及相较于其他算法的优势,以四川省凉山彝族自治州木里藏族自治县理塘河地区(河流深切,岭谷相对高差大)和云南省怒江傈僳族自治州独龙江地区(多高山陡坡,地域环境复杂)为试验区开展对比分析,分别采用区域生长算法、MDRO算法、LFE算法和本文算法针对地貌差异明显的不同地区细小河流进行提取,其局部原始影像和4种算法的提取结果如图8所示。

从图8中可以看出,区域生长下的细小河流,虽已提前设定6个种子点,但提取结果仍存在欠生长现象,甚至部分区域由于混合像元效应出现明显过生长现象。MDRO算法所提取的细小河流结果,河流不连续现象明显,虽保留了河流绝大部分基本结构信息,但信息丢失比较严重。基于LFE算法的河流提取结果由MDRO改善而得,所提取细小河流较为准确,信息丢失相对没有那么严重,但提取结果不稳定,无法确保其完整连续性。本文方法提取的细小河流取结果中,准确地保留了河流的每一个细节信息,河流结果与原河流完全吻合。

2.2.3 细小河流全局对比

完整的河流提取大多是在较大尺寸(地理范围)的影像上进行,在此过程中河流会产生多处不连续,导致待连接的河流端点较多,容易造成错误连接。为此,选取相对完整的独龙江支流区域为试验对象,来进一步验证本文方法的稳定性。河流全局原始影像和4种算法提取结果对比如图9所示。

提前设定13个种子点的区域生长提取结果中,细小河流信息相对完整,但图9(b)中右上角区域由于河流较细且与周围地物间灰度级差异较小,导致出现局部无关信息过生长现象。MDRO算法提取结果中河流主体信息保留相对完整,但残留信息相对散碎,河流断点较多,且图9(c)中右下方河流信息基本未被保留。LFE算法提取结果中,整体河流轮廓清晰,信息保留相对完整,但仍不可避免河流不连续现象。本文方法较好地克服了多断点易错误连接的障碍,精准地处理了河流弯道。

为分析4个试验区域的不连续细小河流断点连接情况,对相关阴影噪声去除影像进行统计,得到其间断河流应连接点对的参考数据,并对比假彩色影像和本文方法提取结果,进而获得断点连接情况,如表1所示。从表1中可见,4个试验区域的河流间段情况各不相同,但连接情况几乎一致,连接正确率均在85%以上,实现了绝大部分间断河流准确连接。本次试验中断点最多的全局河流中,共出现6处错误连接,其中5处为连接偏移,仅有1处为待连接点对寻找错误。此外,试验中连接最佳的河流,未出现错误连接点对,实现了完美连接。

本文方法从细小河流像元间相关关系进行考虑,充分利用了河流像元间的SSIM,大大降低了无关信息的干扰,使得河流即使局部区域被阴影遮挡,发生间断,仍能被准确连接,较好地解决了传统方法的不足。通过我国滇西和川西4个试验区的结果表明,相较于传统细小河流提取方法,本文方法较好地实现了间断河流的准确连接,在一定程度上,保证了河流信息提取的完整性。

本研究的不足之处在于,启发式搜索依赖于间断河流起始点与目标点之间距离的确定,当河流弯度较大或聚集复杂时,河流待连接点对往往会出现评估异常,导致河流连接出现“偏移”现象。因此,如何进一步自动纠正河流走向,成为今后的研究目标。