1 引 言

水体分布信息的准确获取对于水资源调查、流域综合治理、水利规划、洪水监测与灾害评估等领域具有重要的意义；水体的形成与消失、扩张与萎缩对区域乃至全球的生态环境演化和气候变化过程产生重大影响；准确及时地获取水体的动态变化信息更是揭示区域、全球生态环境演变规律的重要技术手段。在陆地表层系统中，水体主要表现为冰川、湖泊、河流等形态。由于水体的反射光谱特征比较明显，而且在空间分布上表现为很强的连续性和均质性，水陆界线清晰，为利用遥感手段快速、准确地获取水体分布及其变化信息提供了可能。国内外学者很早就对此展开了研究，提出了许多有效的水体信息提取模型。如：比值法^[1]、光谱知识法^[2]、光谱分类法^{[3, 4]}、谱间关系法^{[5, 6]}、监督分类法^[7]、特征变化法^[8]、特征分割法^{[9, 10]}、决策树模型^[11]、指数模型法^{[12, 13, 14, 15]}等。

另外，为了提高提取的精度，一些学者对现有的方法进行了改进，如文献^[16]通过比较5种不同水体提取指数模型的原理分析，结合具体实例进行水体遥感提取，确定了在不同时期和不同用途时可采用最佳的水体提取模型；为了实现水体提取的自动化和精确化，文献^[17]将指数模型同分层分类相结合，提出了一种分步迭代的自动化湖泊水体提取方法，一定程度上实现了对湖泊水体信息精准、自动化的信息提取。

但是这些算法基本都是针对光谱性质相对均一的面状水体(湖泊)，对于线性水体(河流)的研究目前相对较少，整体精度不高，自动化水平不强。尤其是河流的支流的提取效果并不理想^[6]。主要原因在于，对线状的河道水体来说，一方面由于水体空间分布形态呈细小狭长状，受到流域环境的影响大，光谱混杂，另一方面由于河流丰水期和枯水期泥沙含量不同，河流上、中、下游的光谱信息差异显著。所以，相对于湖泊等面域水体来说，河道细小线性水体高精度自动化提取比较困难。

鉴于此，本文在ASTER DEM数据生成水系的辅助下，采用空间“全域－局部”分步迭代方法，提出了一种基于水体指数的自适应阈值分割和空间迭代优化相结合算法。并利用DEM数据提取的河网粗略位置信息，融入线状水体空间分布特征，降低信息提取过程中大的环境背景信息的干扰，建立了一种遥感河道线状水体信息提取的技术方案，实现了河流提取的精确化、自动化，提高了计算效率。

河流作为一种独特的地理单元，是在陆地表面凹槽处形成的比较狭窄而修长的水域，上游相对狭窄、下游相对开阔，支流相对狭窄、干流相对开阔。在遥感影像上，河流一方面具有水体光谱特征的共性，河流与地表背景之间光谱反差明显，边界清晰，局部区域内的河流水体光谱相似性强，易聚集成斑块状；另一方面，不同的流域环境下的河流水体(如浑浊度、悬浮物、冰冻、植被)具有较大的光谱差异，且很多支流河道宽度不足一个像元，在遥感影像上表现为混合像元，相对于背景光谱信息，河流光谱信息显得较弱。本文综合考虑了河流的地理空间分布特征和在遥感影像上光谱表征，根据“全域－局部”分步迭代思想，设计了DEM的辅助下图谱迭代反馈机制算法，通过局部水体光谱指数特征阈值的自适应的调整，实现河道水体的自动提取。其总体算法流程见图 1。

图 1 DEM辅助下的河流水体自动提取总体流程 Fig. 1 The flow chart of river water extraction assisted by DEM data

图选项

(1) 影像预处理。从Landsat遥感数据的头文件中获取辐射校正系数、成像时太阳角度等信息，对多光谱影像作大气校正，计算地表反射率，恢复遥感影像的物理属性。获取与遥感数据覆盖范围相同的数字高程数据(DEM)，相互配准，统一坐标体系。

(2) NDWI计算与DEM提取水系。根据归一化水体指数计算公式，计算遥感影像的归一化水体指数；利用DEM生成水系分布图，确定河流大致可能的分布范围。

(3) 空间分析与初始阈值分割。根据水系分布图提供的河流位置信息，根据周边的高程信息以及距离信息，向两边缓冲扩张，形成河流提取的目标区，减小整体背景光谱信息干扰。通过对目标靶区内的NDWI特征分析，确定初始分割阈值，完成NDVI初始分割。

(4) 局部自适应迭代。通过对局部窗口内部的河道水体单元缓冲操作，局部自适应阈值的判断，迭代反馈运算，提高河道水体提取的精度。根据上一步迭代中新增水体单元的数量同当前的数量比较，判断迭代计算是否终止。

(5) 提取结果后处理。根据上述的算法流程提取形成栅格专题层，存在一些小斑块和部分河流水体的断线，需要对其进行去除小斑块、对象化等后处理，获得最终的河流专题图。

在遥感信息提取过程中，最为常见的方法就是构建地物波谱指数模型，通过波谱特征的数学运算达到增强有效的地物诊断特征并抑制其他信息。水体指数通过波段差异比值来增强水体与背景之间的差异。目前关于水体指数的算法有很多，分别适用于不同类型的水体提取。本文采用文献^[12]定义的归一化化水体指数(normal difference water index,NDWI)

点击浏览公式

式中，ρ_green、ρ_NIR分别为绿波段和近红外波段的光谱反射率。

基于栅格DEM提取水系的方法有很多，在此不再赘述。文中采用最具代表性的基于地表径流漫流模型的坡面流模拟方法中的D8算法^[18]。其基本原理是依据水总是沿着斜坡最陡方向流动的原理，确定DEM中的高程数据点水流方向，根据水流方向计算每个高程数据点的上游集水区，通过阈值确定属于水系的高程数据点，追踪形成整个水系。该方法具体实现可以分4步：首先，对DEM数据填坑预处理；其次，确定DEM中每一个网格单元的水流方向，根据每一个网格单元的水流方向可以得到整个流域的排水网络；第三，利用网格单元的水流方向，可以计算出每一个网格单元的上游集水面积；最后，通过积水面积内网格单元的阈值设置，生成流域水系。

根据DEM提取的河网水系分布图，可以整体上确定河流分布的大概位置(即，可能存在河道的区域)。依据水系分布的空间位置信息，通过空间距离，对水系两边作缓冲区计算，文中初始缓冲区的大小设置为3.6 km，包括了最宽的河道；并通过的高程落差限制，对初始缓冲区修改，将高程落差大于60 m的像元去除，得到河流信息提取的计算目标区。通过信息计算目标区的确定，可以消除整个大背景光谱信息的干扰(如山体阴影、暗目标、冰川)，有利于提高信息提取的精度。

由于水体像元与周围背景像元组成的直方图呈双峰分布，而非水体像元则呈单峰分布^[19]。 根据NDWI的物理属性，通过对目标区域的NDWI直方图分布特征分析，选择初始阈值T₀。图 2描述了水体单元同背景单元的理想高斯双峰分布，(图 2中u_land表示背景区域直方图高斯分布的均值，u_water表示水体区域直方图高斯分布的均值；σ_land表示背景区域高斯分布的方差，σ_land表示水体区域高斯分布的方差。)通过判断双峰直方图的谷底(本文试验中取-0.12)，对NDWI图像进行初始分割，得到初步的河流水体信息。

图 2 初始水体分割最佳阈值选择示意图 Fig. 2 The sketch map of optimal threshold selection for water body initial segmentation

图选项

从遥感信息提取的角度而言，水体指数阈值分割法是实现水体快速制图的有效方法。但是，最佳阈值的选择往往比较困难。文中通过分步迭代，实现水体信息提取的逐步优化，达到最佳分割阈值的选择。局部自适应迭代水体算法流程见图 3。

图 3 局部自适应迭代水体提取流程图 Fig. 3 The flow chart of local adaptive iteration

图选项

根据全局分割的初步水体提取结果，以3倍水体大小作局部缓冲区扩展。对影像实行局部窗口运算，初始窗口的设置以水系为中心，窗口大小与局部缓冲区宽度一致。通过移动窗口，分别计算缓冲区的NDWI直方图；分析直方图特征确定分割阈值(T)；通过局部直方图阈值分割得到局部缓冲区内的水体像元(N)；通过同初始水体信息(N₀)比较，判断是否终止迭代。如果本次迭代提取的结果和上次提取的结果一样，说明水体单元基本稳定，迭代终止；否则将T的值赋给T₀，N的值赋给N₀(T₀=T,N₀=N)；进入下一次迭代计算。

为了验证本文提出算法的有效性，试验选取伊犁河的支流为研究区，如图 5所示。该区地形复杂，河流位于山谷洼地、河面狭窄，河道水体提取相对困难。试验选用的卫星遥感影像是2000年8月5日获取的该区域一景ETM+多光谱遥感数据，轨道号Path 146，Row 30；试验中用到的DEM是来自NASA的ASTER GDEM数据。

图 5 河流水体初始阈值分割提取结果 Fig. 5 The result of river water initial segmentation

图选项

根据文献^[18]的坡面径流模拟法提取DEM水系算法(即D8算法)，生成的区域水系分布图如图 4所示。图 4(a)为DEM影像；根据D8算法原理，每个单元格的水流有8种可能的流向，即东、东南、南、西南、西、西北、和北、东北，并分别用1、2、4、8、16、32、64和128这8个有效方位特征码表示，计算坡度最大的方向定义为中心栅格单元的水流方向，见图 4(b)；根据水流方向栅格矩阵，计算出每个栅格单元的上游集水面积(产流)，见图 4(c)；通过产流面积阈值分割，将上游集水区面积大于50的划为水系单元，连接水系单元生成流域水系网络，得到区域水系分布矢量图，见图 4(d)。

图 4 DEM 生成流域水系 Fig. 4 Using DEM to generate drainage network

图选项

依据DEM生成的流域水系，提供了河道存在的大概位置信息。通过空间分析，自动建立水系两侧缓冲区即河道水体提取的目标区域，见图 5(b)。通过水体指数(NDWI)计算(见图 5(c)和全局阈值分割，得到流域水体的初步分布信息(见图 5(d))。

在NDWI初始阈值分割的基础上(图 6(a))，根据局部自适应迭代算法，对初始的水体提取结果进行局部迭代计算，以提高水体信息识别的精度。通过局部自适应迭代计算，并对一些小斑水体(非河流)进行去除，得到局部迭代优化后的河流水体信息(图 6(b))。

从图 6(b)中可以看出，利用DEM辅助下的多光谱遥感影像河道线性水体自适应迭代提取算法，能够快速、准确地提取河道水体信息。但是由于桥梁、沿岸植被、河床地形等因素的影响，部分地区细小河流水体存在着断线现象。这是多光谱遥感影像上记录的客观现象。为此笔者需要通过数字图像后处理技术，实现对河流断线现象的自动连接。文章中采用数学形态学方法对河流提取的结果进行闭运算计算(close)^[20]

点击浏览公式

式中，X为提取结果，即二值影像；B为一个3×3的结构模板。

图 6 局部自适应迭代优化与后处理 Fig. 6 The result of extracted river water after local adaptive iteration optimization and post-processing

图选项

首先通过一个3×3的结构模板与提取河流结果二值影像进行“或”操作，如果都为0，则输出的结果为0，否则输出结果为1，实现河道水体膨胀处理，通过此步骤可以连接河流之间的线状断流区域(见图 6(c))。然后对膨胀后的河流区域腐蚀细化处理，理论上的腐蚀操作是通过结构模板的“与”操作，如果都是1，则输出的结果为1，否则输出结果为0。但是，实际应用中河道细化算法需要考虑的情况要比膨胀算法复杂得多。一方面要实现河道的细化，另一方面要从8邻域保证河道的连通性，为此文中通过比较采用了经典的Hilditch算法^[21]。最后自动跟踪矢量化，得到河流的水系分布图(见图 6(d))。

通过伊犁河上游山区细小河道的提取试验，结果表明：本文的算法在提高精度的同时，避免了遥感信息计算的盲目性，增强了遥感信息提取的针对性，大大地提高了计算效率(见表 1算法效率与精度分析)。从算法的执行效率看，本文的算法耗时差不多是全局-局部迭代算法的一半；特别是在河流断线方面，人工后处理编辑操作工作量很大，文本的方法提取的结果正常情况下几乎没有什么断线，只需要打开查看、检查错误，大大降低了人工后处理工作量。另外，为了验证本文提出算法对河道细小线性水体的识别精度，通过人工目视解译的方法在河道沿岸随机选取了80个验证点，其中有50个验证点分布在水体区域，30个验证点分布在非水体。通过目视解译验证，50个水体验证点，全局阈值分割法识别为水体的像素有36，全局-局部迭代识别为水体的有44个，本文提出的算法识别为水体的有47个；30个非水体验证点，判断为非水体的分别是：24、26、29个。从算法提取的精度看，全局阈值分割算法的总精度为75%，全局-局部迭代水体提取算法的精度为87.5%，而文中的方法河流水体提取总体精度达到了95%。尤其是非水体的辨伪率，由于本文引入空间位置信息、通过空间分析，消除了背景环境因素的影响，提高的尤为明显。

本文针对线性河流水体的光谱一致性差，易受到流域环境、地表覆被的影响，导致水体的自动化提取困难等问题。在“图-谱”耦合和“全局-局部”迭代算法思路的基础上，进一步提取了DEM空间知识辅助下的河道水体自适应迭代提取。该方法充分利用了DEM生成水系的“图”空间特征和多光谱遥感影像上水体的“谱”特征，通过NDWI的水体分割阈值选择和局部自适应迭代机制，实现了河流细小水体的自动化、精确化提取。试验结果表明，该方法能够精确地识别出河流水体，可以为流域水资源遥感调查、水利规划、洪涝灾害监测等提供技术支撑和决策支持。通过伊犁河支流的河流水体提取试验，研究得到以下结论：

(1) “图”信息和“谱”信息是遥感影像固有的特征，也是遥感影像区别于其他影像的重要特征，综合利用图谱信息，实现优势互补，可以提高遥感影像自动解译的精度。

(2) 波谱指数计算是水体快速提取的有效方式。通过融入空间位置信息，建立信息提取目标区，可有效地突出河道水体信息，减弱背景光谱的干扰，改善遥感信息计算的针对性和目的性；有效地避免了整体影像输入的盲算，提高了计算效率。

(3) 由于河流水体光谱的非均一性，受流经区域环境影响很大，河流的上中下游光谱特征差异明显；通过局部自适应迭代计算，可有效地挖掘局部细小水体，提高分类精度。