基于改进Canny算子的水域边界自动提取方法研究

郑云柯，王俊霖，张世涛

（昆明理工大学 国土资源工程学院，云南 昆明 650093）

0 引言

水是人类生存必不可少的资源，是大自然赋予人类最宝贵的财富［1］。随着全球气候变化和人类活动的不断增加，地球上水体的分布和特征发生了显著变化，给人类社会和自然生态带来了严重挑战，水资源的可持续性问题引起人们越来越多的关注。因此，及时、准确、全面地监测和分析水体信息对于理解水循环过程、评估水资源状况、保障水资源的可持续开发利用至关重要。卫星遥感技术因其能够高效、准确地获取全球范围内的水体信息而备受关注。利用卫星遥感技术，可以获取大量水体信息，从多源、多时相、多尺度的遥感影像中提取和反演水体信息，包括水体范围、深度、透明度、色度、温度等参数，并结合地面观测数据和数值模型进行验证与分析，可以实现对大范围、复杂区域内水体变化情况的动态监测和定量评估，为制定合理、有效的水资源开发利用和管理措施提供科学依据［2］。如今，快速、准确地从卫星遥感影像上提取水体信息已成为水资源调查、水资源宏观监测及湿地保护的重要手段。

1 相关研究

目前，利用遥感技术监测湖泊主要集中在湖泊的水质［3］、水位［4］、面积［5］等方面。常见的水体信息提取方法有单波段阈值法［6］、多波段谱间关系法［7］和水体指数法等。其中，水体指数法是一种简单、快速、高效的水体信息识别方法，如Mcfeeters［8］在1996 年提出归一化差异水体指数（Normalized Difference Water Index，NDWI），利用绿光波段和近红外波段增强了影像中的水体特征；徐涵秋［9］在NDWI 的基础上，提出改进的归一化差异水体指数（Modified Normalized Difference Water Index，MNDWI），采用中红外波段代替近红外波段，取得了比NDWI 更好的提取效果；闫霈等［10］提出增强型水体指数（Enhanced Water Index，EWI），相较于NDWI 和MNDWI，能更有效地区分河道与背景噪音。针对不同研究区的水体特性，以上水体指数法都取得了较好效果。

湖泊作为具有一定面积的相对封闭水域，为了获取湖泊的矢量边界，常用的方法包括手工数字化方法、区域分割法［11］与边缘检测法［12］。当前国内外学者对水边界提取的研究主要集中在海岸线与湖岸线，如Alesheik 等［13］利用阈值分割法对乌尔米耶湖岸线进行提取；李秀梅等［14］利用Canny 算子对渤海湾海岸带进行提取以监测其时空变化；Karantzalos 等［15］基于遥感影像，分别用Laplacian 算子和 Canny 算子提取海岸线；申家双等［16］在分析了边缘检测算法用于影像水边界提取的优缺点后，提出将Canny 算子与GAC 模型结合提取影像水边界的方法；Abolhassani 等［17］利用边缘检测算法提取了美国东海岸线；魏东岚等［18］在MATLAB 平台上利用边缘检测算法、阈值分割法及小波变换法提取海岸线，结果表明边缘检测法中的Canny 算子能准确提取出围填区域的海岸线。

已有研究主要使用Canny 边缘检测算子提取海岸线，且Canny 算子的阈值多是通过大量实验获取，这是Canny算子自身的不足：高低阈值需人为设定，容易导致检测中出现大量虚假边缘［19］。因此，本文在已有研究基础上，以湖泊为研究对象，将Canny 边缘检测算子与遥感水体指数法相结合以提取湖岸边界。针对Canny 算子的不足，利用类间方差最大阈值分割法来计算其阈值，使Canny 算子具有较理想的阈值，从而实现自动化提取湖泊水边界。

2 研究区概况

洱海位于云南省大理市郊区，地跨大理市和洱源县。洱海是云贵高原九大湖泊之一，也是云南省第二大淡水湖，海拔1 980 m，经纬度为100°05'-100°18'E，25°36'-25°58'N。洱海北起洱源县，南至大理市下关镇，湖泊呈东西窄、南北长的条带状，因湖泊形状像耳朵而得名。洱海水域面积250 km2，平均水深10.5 m，平均水位1 974 m，湖水主要靠河流补给。

3 数据来源及预处理

本文采用的遥感数据是从USGS 网站上下载的大理市洱海的Landsat-8 OIL 影像，轨道号/行号为131/42，选取影像的云量少于3%，遥感影像获取时间为2019 年11 月29号，处于湖泊的枯水期。Landsat-8 OIL数据共有9个波段，除全色波段空间分辨率为15 m外，其余波段都为30 m。另外，为验证水体边界提取精度，从中国资源卫星应用中心获取了一幅时期相近的大理市洱海的2 m 全色/8 m 多光谱GF-6号影像，遥感影像获取时间为2019年12月7号。

对Landsat-8 OIL 影像进行辐射定标、大气校正、图像裁剪等一系列预处理，对GF-6 号影像进行2 m 全色波段和8 m 多光谱融合，形成空间分辨率为2 m 的融合影像。

4 研究方法

4.1 水体信息增强

国内外学者根据水体在蓝绿波段吸收率较低、近红外波段吸收率较高这一光谱特征，开展了大量基于遥感影像的表面水体自动提取算法研究。本文将利用不同的影像波段，分别采用归一化差分水体指数法、改进的归一化差分水体指数法和增强型水体指数法来增强洱海表面水体信息，如表1所示。

Table 1 Model of the water body index表1 水体指数模型

4.2 水体边缘检测算法

传统的边缘检测算法有Roberts 算子、Prewitt 算子、Sobel 算子等，这些算法简单，虽容易实现，但处理噪声的能力较差，裂纹边缘识别不完整，还容易出现伪边缘现象。能否选择合适的边缘检测算子，将直接影响到结果的精度。笔者对前人的研究结果进行比较，发现使用Canny 算子处理边缘和噪声的效果最佳［20］。因此，本文选用Canny算子对水体增强的遥感影像进行边缘检测，主要分为以下4个步骤：

（1）利用高斯函数对图像进行平滑处理。高斯函数表达式为：

平滑图像后的 结果I(x，y) 表达式为：I(x，y)=G(x，y)*f(x，y)，其中G(x，y)为高斯函数，f(x，y)为原图像。

（2）对平滑后的图像进行梯度幅值和方向的计算。基于2×2 模板，通过X、Y方向像素的一阶导数来确定梯度幅值。设Fx(x，y)为X方向的偏导数，Fy(x，y)为Y方向的偏导数，则梯度幅值M(x，y)和方向θ的表达式为：

（3）对梯度幅值的非极大值进行抑制处理，判断梯度幅值在其八邻域内是否为最大值，若为最大值则为边缘，否则置零。

（4）对梯度选取两个门限阈值，即高门限阈值TH和低门限阈值TL，两者通常的关系为［21］：

取出经过非极大值抑制后图像中的最大梯度幅值，标记大于高门限阈值的点（边缘点），将小于低门限阈值的点置0，即可提取出完整边缘。

4.3 Canny算子阈值设定

由Canny 算法步骤可知，门限阈值的选取是进行图像边缘提取的关键，而Canny 算子需要人为预先设定高低阈值，因此需要很多次反复实验才能找到合适的阈值。为解决阈值需要人为预先设定的问题，本文利用类间方差最大阈值分割算法来确定Canny算子的门限阈值。

类间方差最大阈值分割算法别称大津（Otsu）法［22］，是由日本学者大津在1979 年提出的一种自适应阈值确定法。其数学描述为对于给定的图像L(x，y)，设图像的灰度级为M，灰度值h像素点的像素数为n，则其出现的概率为［23］：

采用阈值t为界限把图像区分为两个区域，假设背景像素区域为A1，前景像素区域为A2，则两个区域的总概率为：

设图像整体的灰度均值为η，则前景像素区域A1与背景像素区域A2的灰度均值为：

即当σ2为最大值时，对应的t为选择的最佳阈值。

根据上述Otsu 算法原理确定一个分割阈值t，将选取的分割阈值t作为Canny 算子的高阈值TH，再用此高阈值乘以一个比例因子0.5 作为其低阈值TL［24-25］。得到阈值的表达式如下：

由此Canny 算子的高低阈值选择只与分割阈值t有关，不需要再人为多次反复实验设定，这就为Canny 算子中阈值的确定问题提供了一个较好的解决方法，也增强了Canny算子的分割能力和自适应性。

总体技术路线如图1所示。

Fig.1 Technical flow图1 技术流程

5 结果与分析

5.1 水体边界自动提取

5.1.1 水体灰度差增强

如图2 所示，NDWI、MNDWI 和EWI 水体指数法都能有效拉伸水体和非水体之间的灰度差异。从目视解译的效果来看，3 种方法拉伸灰度差异的效果各不相同。相比于MNDWI 法和EWI 法，NDWI 法对洱海南边的出水口即西洱河的提取有所欠缺。

Fig.2 Water body enhancement effect图2 水体增强效果

5.1.2 边缘矢量化

先对3 幅水体增强影像进行归一化处理，由大津法的数学原理可知其特性［26-27］：当目标地物与背景地物的面积比例悬殊时，灰度直方图可能会无明显“双峰”或呈“多峰”形态，此时使用大津法的分割效果不佳。为得到最优阈值，需满足灰度频率直方图呈“双峰”的情况［28］。通过图2 可看出目标地物水体与背景地物非水体的面积比例相差不大，研究发现3 种水体指数法的灰度频率直方图均呈“双峰”分布，如图3 所示。其中，横坐标Data Value 代表归一化后的灰度像素值。

Fig.3 Frequency histogram of water index method图3 水体指数法频率直方图

然后，使用Ostu 算法得到水体增强影像的分割阈值（保留两位有效数字），最后利用关系式（7）得到Canny 算子的高阈值TH 和低阈值TL，结果如表2所示。

Table 2 Thresholds of various water index methods表2 各类水体指数法阈值

将提取好的栅格影像结果转换为矢量数据，再与原始遥感影像叠加显示。如图4 所示，整体上，水体边界提取的边缘连续性及去噪效果都达到最佳。3 种水体指数法均可提取出洱海的水体边界，对于洱海的北部、西部、东部的非目标地物类，也均提取了其边缘。就北部的目视效果而言，MNDWI+Canny 算子提取的非目标地物略少一点。相比其它两种水体指数法，NDWI+Canny 算子还提取出了位于洱海西南部的苍山积雪边缘及南部的部分建筑边缘。因此，选取合适的边缘检测阈值，可减少目标地物的非必要干扰信息，为后续的处理节省时间。

Fig.4 Vector boundary overlay image of water body图4 水体矢量边界叠加影像

5.1.3 边缘细部提取效果

选取洱海湖岸的湿地、河滩及洱海出水口3 种典型地区进行分析，湿地是水体、植被与裸地等按不同比例混合组成的一种土地形式。洱海湖岸周围存在着不少湿地区，比较著名的有海舌湿地公园、洱海月湿地公园等。以海舌湿地为例，分析其湿地边界提取效果。如图5 所示，运用NDWI+Canny 算子提取出的海舌湿地边缘出现明显的不连续现象，且产生了细微的伪边缘现象。其将湿地的近水侧岸线和近陆侧岸线均提取了出来，整体上不能很好地提取出完整岸线。运用MNDWI+Canny 算子提取出的海舌湿地边缘连续、完好，可以看出其将湖岸湿地划分给了水域。运用EWI+Canny 算子提取出的海舌边缘同样出现了不连续现象及少量断点，位置不够准确，可以看出其将湖岸湿地划分给了陆地。

Fig.5 Sea Tongue wetland boundary图5 海舌湿地边界

河滩是由于泥沙沉积而形成的天然滩涂土地，选取洱海湖岸的某处河滩进行对比分析。如图6 所示，3 种方法所提取的边缘都是连续的。NDWI+Canny 算子在大河滩（右下方）处和小河滩（右上方）处均有向外膨胀的现象。MNDWI+Canny 算子在提取河滩边界时有较好效果，可以看出其将河滩也划分给了水域。EWI+Canny 算子则在大河滩处多提取了一处伪边缘，在小河滩处出现了向外膨胀的现象。

Fig.6 River beach boundary图6 河滩边界

位于洱海西南方向的湖岸出水口称为西洱河，洱海与西洱河以大理市的兴盛大桥为界，该地附近还有著名的洱海月湿地公园，选取该区域进行分析。如图7 所示，3 种方法都能有效、连续地提出洱海出水口边界，但NDWI+Canny算子和EWI+Canny 算子多提取了大桥边缘及西洱河（桥梁左侧）的部分边缘，且两种方法多提取的部分几乎是一致的。MNDWI+Canny 算子则恰好提取了出水口边界。宏观洱海月湿地的提取效果（位于出水口的左上方），3 种方法并无明显差异，均能有效、连续地提取出边缘；微观洱海月湿地的提取效果，NDWI+Canny 算子和EWI+Canny 算子都出现了向外膨胀的现象，而MNDWI+Canny 算子提取的边缘效果较好。

5.1.4 边界整体提取结果

由于影像分辨率的问题，混合像元的存在使Canny 算子无法完全克服噪声，导致提取的边界存在微小的不连续情况。叠加原始影像，将存在的微小边缘断裂进行连接，并剔除非目标地物的边缘信息，输出3 种水体指数法的洱海水域矢量图，如图8所示。

Fig.8 Automatic extraction results of water boundary vector data图8 水体边界矢量数据自动提取结果

5.2 水域边界精度评价

5.2.1 目视判别

将提取的洱海边界与时期相近的GF-6 号2m 融合影像进行叠加，如图9 所示。目视检查发现NDWI、MNDWI和EWI 都能有效提取出边缘，MNDWI 更接近精确的边缘位置，NDWI和EWI都略微有向外膨胀的现象。

Fig.9 Gf-6 2m fusion image of superimposed water boundary vector data图9 叠加水体边界矢量数据的GF-6号2m融合图像

5.2.2 空间统计分析

基于空间缓冲区统计法对提取精度进行量化分析，从统计学角度提出3 点准则：①向水缓冲区B 的最大、最小及平均值都要小于向陆缓冲区A；②B 内标准差要小于A；③两个缓冲区的标准差Δσ=σA-σB要尽可能大［29］。

因水体对近红外波段的吸收率较强，故这里选取Landsat-8 OIL 近红外波段影像作为统计的基础影像。以水边界为基准分别在水陆两侧建立缓冲区A 和B，缓冲范围为2 个像素值，即60 m。再统计缓冲区域的像素灰度特征（DN）值，如表3所示。

Table 3 Characteristic statistics of pixel DN value in buffer表3 缓冲区内像素DN值特征统计

基于表3 的特征值分析，洱海水体边界提取位置较准确的是MNDWI+Canny 算子。

6 结语

本文选用2019 年的OIL 影像作为研究数据，以洱海为研究对象，分别运用归一化差分水体指数、改进的归一化差分水体指数和增强型水体指数对遥感影像进行水体信息灰度增强。分析3 种水体指数法的灰度频率直方图，使用大津法得到水体灰度影像的分割阈值。将分割阈值作为Canny 算子的高门限阈值，结合Canny 算子高低门限的阈值关系式，得到低门限阈值以自动提取洱海水体边缘信息，进而叠加原始影像进行典型地区分析。将提取的水边界叠加在GF-6 号2 m 分辨率的融合影像上进行目视判别，再经空间缓冲区分析得到：在NDWI、MNDWI 与EWI 3种水体指数法中，MNDWI 与Canny 算子结合的提取效果最好。

另外，本文中存在些许边缘不连续的情况，这是由于水体指数法大多采用具有中红外波段的遥感影像，混合像元的存在对提取结果仍具有一定干扰性。Canny 算子无法完全克服噪声影响，今后可将混合像元分解或将蚁群算法等技术应用于研究中。