1 引 言

从遥感影像中提取水体信息，是土地覆被制图、水域分布数据更新、数字流域构建等技术的关键支撑，并为水资源调查、河流健康评价、灾害评估等应用提供快捷而准确的信息^[1]。较植被、城镇等土地覆被类型而言，水体的光谱反射率总体偏低，特征明显^[2]。常见的水体提取方法分为4类：单波段阈值法^[3]、多波段谱间关系法^[4]、水体指数法^{[5, 6]}、监督/非监督分类法^{[7, 8]}。然而，单纯依赖光谱信息，无法有效识别宽度接近影像分辨率的细小河流^[1]。原因在于，细小河流深度总体较浅，其像元可能成为包含水体、悬浮泥沙、河床、河岸等非水体地物的混合像元，易与其他地物类型像元混淆^[9]。

河流具有独特的空间特征，能够用以提高其识别精度^[10]。针对TM/ETM+影像，文献^{[1, 11]}利用谱间关系法初步识别河流，然后利用数学形态学方法进行结果修整。文献^[12]认为，数学形态学方法适于连接长度小于等于5个像元的间断，否则被连接的河流将发生变形。文献^[9]将数字高程模型引入河流识别过程。文献^[13]针对SPOT影像，利用DRO(Duda’s road operator)算子结合启发式路径搜索算法提取河流。文献^[14]针对高分辨率影像，利用基于空间像素模板的Adaboost方法，有效利用邻近像素的类别来辅助区分河流与其他地物。

本文提出一种利用改进的归一化差异水体指数(modified normalized difference water index，MNDWI)^[6]结合多种数字图像处理技术的细小河流提取方法。以TM遥感影像作为数据源，通过在我国北方4省中研究区的试验结果认为，该方法能够提取到90%以上的河流，具有较高的普适性。

本文定义的细小河流其宽度小于等于3个像元。因为对于更宽的河流，其中心像元一般为水体纯像元，具有面状水体的光谱特征，可以通过阈值法进行提取。

本文采用区域划分的方法以获取完整的水系分布(如图 1)：对于面状水体，利用全局阈值进行提取，并利用局部阈值确定其边缘，具体方法在2.1节介绍；对剩余区域，利用图 2中描述的方法提取细小河流，具体方法在2.2至2.4节介绍。这两种方法相互补充，通过合并得到最终结果，如图 1。

图 1 提取完整水系分布的流程 Fig. 1 Flow chart of water extraction

图选项

图 2 细小河流提取流程 Fig. 2 Flow chart of narrow river extraction

图选项

本文采用MNDWI表达水体信息，MNDWI的计算公式如下^[6]

点击浏览公式

式中，ρ_green与ρ_swir分别对应TM影像的第2和第5波段。经过试验认为：取0.3这个较高的阈值能有效识别面状水体，且避免绝大多数阴影、城镇的干扰，只是在水陆交界处会漏分部分像元。因此，在已提取像元邻域中，补充满足MNDWI≥0的像元作为水体边缘(图 1)。之后，对面状水体以外的区域(图 1)，利用接下来介绍的方法提取细小河流(图 2)。

线状特征增强能够突出影像中的线状地物特征，同时抑制其他形状地物特征。河流形状较为复杂，可能笔直(人工河道)也可能蜿蜒曲折(河流源头)，为提高普适性，本文对DRO算子进行了修改：①只考虑中心像元值大于左右两侧像元值的像元，避免了DRO算子在阶梯状边缘产生响应的缺陷；②只考虑3个像元，相比考虑9个像元DRO算子，对于蜿蜒曲折的自然河流效果更好。

该算子包含4个方向算子(图 3)，计算方法为：对每个方向算子利用式(2)求得方向性线状增强结果LFE_i，并且取LFE_i的最大值为线状特征增强结果

点击浏览公式

式中，i代表 4个方向算子的序号；而a、b和c分别对应各方向算子中相应位置的像元值(图 3)。这里将本步骤的滤波结果简称为LFE(linear feature enhancement)。

图 3 线状特征增强算子 Fig. 3 Linear feature enhancement operator

图选项

Canny 算法是数字图像边缘检测的经典方法^[15]。该方法用高阈值筛选强边缘，然后用低阈值追踪与强边缘相连的弱边缘。本文借鉴该方法：利用双阈值进行河流线段追踪。具体方法为：首先，把所有满足高阈值(公式(3))的像元标记为河流；再将剩余像元中与河流相连接，且满足低阈值(公式(4))的像元标记为河流。

点击浏览公式

点击浏览公式

在提取细小河流的过程中，容易混淆的噪声主要有3种类别，其成因和特点如表 1所示。

(1) 山体阴影：由于现有的ASTER DEM数据在水体区域质量不高^[16]，不利于采用地形校正方法去除阴影^[17]。文献^[18]认为，阴影和水体在TM绿光波段中灰度值差异最大。因此，以绿光波段灰度值小于T_shadow的像元作为阴影。而如果影像地形平坦，则设置0忽略此功能

点击浏览公式

(2) 道路：道路的MNDWI值接近于裸地、城镇，而高于植被。以植被为背景的道路作为线状信息也会被增强。而在短波红外波段，道路的反射率高于植被，而水体的反射率低于植被。因此，如果该像元满足公式(6)，则作为道路而排除

点击浏览公式

(3) 其他类型：产生原因较多，如呈线状分布的空地，或者机场、煤场等人工设施；之前消噪步骤残留的像元。这些噪声的共同特点是长度短。因此，判断其像元数量(公式(7))，删除较小斑块

点击浏览公式

本文根据不同的土地覆被类型，挑选我国北方地区的四景TM影像进行试验，其基本情况见表 2。通过人工解译得到参考结果，利用以下4个指标来评价提取效果：制图精度(producer’s accuracy)、用户精度(user’s accuracy)、Kappa系数(kappa coefficient)、完整度(completeness)^{[19, 20]}，结果见表 3。其中，前3个指标基于逐像元的对比，而完整度基于河流的长度，公式如下

点击浏览公式

式中，L_ref 、L_t分别为参考、正确提取的细小河流长度。

试验1展示了本文方法主要步骤的中间结果。图 4(a)为TM 543波段组合影像，其中较大城镇为漠河县城。其中4条河流分别是额木尔河(图 4(e)中标注为1)、大林河(3)、老槽河(2)和古莲河(4)。图 4(c)是阈值法提取结果，其中额木尔河和大林河的提取结果较完整，而老槽河、古莲河几乎完全无法提取到，制图精度极低。而对比图 4(d)和图 4(b)可以看出，经过线状特征增强后，河流和道路在影像中十分突出。图 4(e)为利用双阈值线段追踪方法的初步提取结果，4条河流被完整提取，且相互连接无间断，但其中混合有河滩和道路等噪声，最明显的是4段道路，图中均标以星号。而图 4(f)显示了去噪后的结果，不足之处是老槽河有一处间断，此处的河宽经Google Earth测算小于15 m。

图 4 试验1中原始影像和主要步骤结果(白色区域表示提取结果) Fig. 4 Experiment 1,comparison of original image and results of the main steps (white pixels indicate the extraction results)

图选项

试验2的影像中山体阴影较多，阈值法的用户精度在4个试验中最低，因为结果混有大片的阴影和城镇，如图 5(d)所示。而利用本文方法，用户精度显著提升，可以提取影像中两条河流的主要部分，然而其中也存在3处明显间断，其中河宽最宽处为20 m。

试验3的影像中有黄河和4条人工河道。其中，对于黄河等面状水体，两种方法提取效果差别不大，如图 5(e)和图 5(k)。但阈值法识别人工河道的效果不佳。而利用本文方法时，4条人工河大部分被提取，右侧一条水渠发生断裂，经测算其河床宽度25 m左右。

试验4的影像位于青海省果洛州优云乡，黄河与优尔曲的交汇处。由于靠近黄河源头，此处的河流蜿蜒曲折，分支交错，河流宽度变化大。黄河宽度大于40 m，而优尔曲的宽度在3 m至15 m不等。结果显示，利用阈值法可以提取黄河干流，但优尔曲只能提取到一些微小片段，如图 5(f)。而利用本文方法可以提取黄河和优尔曲的干流，如图 5(l)。

图 5 试验2—4，对比原始影像和提取结果(白色像元) Fig. 5 Experiment 2 to 4, comparison of original images and extraction results (white pixels)

图选项

混合像元效应是导致难以从TM影像中提取细小河流的主要原因。本文提出一种细小河流提取方法：①首先利用阈值分割，识别面状水体；②对水体指数进行线状特征增强；③利用双阈值线段追踪方法，提取影像中的细小河流；④去除山体阴影、道路和其他这3类噪声。通过我国北方4省试验区的结果显示，针对宽度小于等于3个像元的细小河流，该方法对比阈值法有明显改善，提取结果的制图精度高于82%，用户精度高于93%，Kappa系数高于0.993，完整度高于90%。

本研究的不足之处在于：①当河流宽度小于像元宽度时，其水体指数值接近于背景而难以提取，因此，关于算法能够提取的河流宽度，存在一个下限，本文根据试验认为，对于宽度小于等于25 m的河流，提取效果将会不稳定，但这一结论缺乏定量化的描述，本文只是根据试验结果作出简单估计；②文中的双阈值线段追踪方法依赖于给定的阈值，在拐弯、淤塞或者阴影处，往往因为个别像元难以识别而发生间断，如何利用启发式方法来自动连接间断，成为今后的研究内容。