水边线是水体与陆地的分界线^[1], 在遥感影像上表现为传感器成像时刻的瞬时水边线^[2]。瞬时水边线提取是基于遥感技术提取海岸线的前提^[3], 其结果对于海岸带科学管理具有重要意义。

近些年来, 米级和亚米级星载高分遥感技术的发展, 为瞬时水边线精细提取提供了新的技术手段^[4]。尽管目视解译结果的精度较高, 但受限于目视解译的效率, 科研人员一直致力于利用计算机进行瞬时水边线的自动化或半自动化提取。利用计算机进行瞬时水边线的自动化或半自动化提取, 其实质是依据水陆两侧地物光谱或纹理等差异进行图像分割或分类, 并由此产生了多种分割和分类方法, 如经典的阈值分割法、边缘检测法, 以及近些年发展起来的深度学习^[5]等算法。经典的阈值分割法有密度分割法^[6]、大津法^[7]和泛洪算法^[8]等, 在此基础上, 自适应阈值化方法^[9], 以及大津法与影像分割技术相结合方法^[10]也被应用于瞬时水边线或海岸线提取。阈值分割法中, 阈值的选取直接影响分割的效果^[11]。当两种地物的像素值相近, 对比度不强时, 使用阈值分割方法难以做到有效分割。此外, 由于阈值分割法是对影像进行逐像素分类, 高分遥感数据中的噪声会导致阈值分割出的瞬时水边线不连续。而常用的边缘检测法是利用边缘检测算子, 如Laplacian^[12]、Prewitt^[13]、Roberts^[14]和Canny算子^[15-16]进行瞬时水边线提取。利用边缘检测算子在进行复杂海岸带瞬时水边线提取时, 同样容易受高分遥感影像中礁石和白浪等噪声影响, 导致瞬时水边线提取结果产生断裂, 因此后续常采用数学形态学对边缘检测算子提取结果进行优化, 以消除孔洞和毛刺, 使提取结果连续, 但同时也带来瞬时水边线提取结果进一步出现偏差, 以及提取效率显著降低等问题。此外, 最近发展起来的深度学习算法也被应用于遥感影像瞬时水边线提取, 如通过升级残差网络为主干网络, 经过全局卷积网络、密集连接网络和解码器网络配合, 产生海陆分割结果^[17], 但深度学习过于依赖样本的数量, 且受计算机硬件条件限制, 耗费的时间成本较高。

主动轮廓模型最初由Kass提出, 可以实现指定地物初始轮廓自动演化运动到目标边界, 称之为Snake模型^[18]。但Snake模型容易忽略地物边缘的微小特征, 该模型各项参数的调整过程比较复杂^[19-20], 且初始化结果对地物边缘收敛的准确性影响较大^[21]。Malladi等^[22]首次将水平集方法引入主动轮廓模型中, 提出了水平集活动轮廓模型。由Chan和Vese提出的CV模型^[23]也通过了水平集方法对模型进行改进。该类模型能够表示复杂拓扑结构的轮廓和处理拓扑结构变化, 但其在演化的过程中需要不断地重新进行初始化水平集函数。李纯明等人为了避免上述弊端, 提出距离正则化的水平集演化模型(distance regularized level set evolution, DRLSE)^[24], DRLSE模型对水平集公式重新进行定义, 在零水平集区域使水平集函数规则化, 解决了模型中初始轮廓重复进行初始化的问题, 使得影像分割效率大大提高^[25]。但是, 使用DRLSE模型提取出的地物轮廓线, 包含噪声轮廓线(如沿海养殖场和货轮等轮廓), 需要后期人工辅助处理, 才能得到瞬时水边线, 导致基于DRLSE模型的瞬时水边线提取效率较低。

综上所述, 本文在兼顾精度和效率的基础上, 提出一种顾及轮廓信息和DRLSE模型的瞬时水边线快速提取方法, 可有效地实现瞬时水边线分离提取, 并保证瞬时水边线提取结果的连续性。

1 研究区概况和数据预处理 1.1 研究区概况

本文的研究区是福建泉州, 位于福建东南部, 北部与莆田、福州接壤, 南部紧邻漳州、厦门, 西部紧靠龙岩和三明, 东与台湾岛隔海相望, 地理范围为117°25′~119°05′E, 24°30′~25°56′N, 研究区位置示意图如图 1所示。

泉州境内溪流密布, 沿岸的入海河流众多。河流所携泥沙常年堆积, 给泉州沿海带来了丰富的淤泥质海岸线。淤泥质岸线土质肥沃, 为沿海滩涂养殖、开辟盐场等生产活动提供了极为便利的条件。除淤泥质岸线外, 经实地调研, 基岩岸线、砂质和人工岸线也多有分布, 自然岸线和人工岸线共同构成了泉州丰富的岸线资源。

1.2 数据获取和预处理

本文选取了研究区2016年3月6日高分二号(GF-2)多光谱数据, 空间分辨率为4 m, 图幅大小为45 km×45 km。首先, 根据文件自带的参数进行辐射定标, 将DN(Digital Number)值转化为辐射亮度; 其次, 基于ENVI FLAASH大气校正模块进行大气校正, 以减弱由于大气散射和吸收等造成的衰减; 最后, 利用ASTER GDEM数据(第二版)基于C-correction模型^[26]对大气校正后的数据进行进一步地形校正, 以弱化地形影响。

2 研究方法

本文提出的顾及轮廓信息和DRLSE模型的瞬时水边线快速提取方法, 主要包括: 首先对预处理后的遥感数据, 基于DRLSE模型, 在目标水域标定矩形初始轮廓, 据此向地物轮廓边缘处演化, 得到瞬时水边线粗提取轮廓(包括影像的外边缘轮廓), 之后利用初始轮廓信息(周长和中心位置), 消除噪声点轮廓和部分离岸岛礁对粗提取结果的影响, 获得瞬时水边线的精确结果。

2.1 基于DRLSE模型水边线信息粗提取

基于DRLSE模型的水边线信息粗提取主要包括以下内容:

(1) 增强对比度

DRLSE模型将图像梯度信息作为边界停止的条件, 对于弱边界的地物目标边缘信息, 容易产生边缘泄露问题。因此使用DRLSE模型对地物边缘轮廓提取前, 强化感兴趣地物与周围地物的对比度十分重要。本文利用归一化差异水体指数^[27](normalized difference water index, NDWI)增强水体和其他地物的对比度:

$ I_{\mathrm{NDWI}}=\frac{\left(B_{\mathrm{Green}}-B_{\mathrm{NIR}}\right)}{\left(B_{\mathrm{Green}}+B_{\mathrm{NIR}}\right)} , $

(1)

式中, I_NDWI代表归一化差异水体指数, B_Green和B_NIR分别对应于GF-2数据的绿光波段和近红外波段。

(2) 选取初始矩形轮廓

在DRLSE模型中, 初始轮廓是一个封闭矩形, 该封闭矩形在DRLSE模型中作为输入条件, DRLSE模型的演化过程可以认为是由初始矩形轮廓演化到地物边缘轮廓的一个过程。本文将初始矩形轮廓选定在水体区域任意范围内, 以便后续利用初始矩形轮廓的周长和中心位置对提取出的地物轮廓进行筛选, 从而直接得到目标地物(瞬时水边线)边缘信息。

(3) 提取水边线信息

由于DRLSE模型^[24]的本质是以图像的梯度信息作为边缘检测的标准, 因此, DRLSE模型也是边缘检测算子的一种方法, 其能量泛函公式如下所示:

$ \begin{aligned} E(\phi)= & \mu \int_{\Omega} p(|\nabla \phi|) \mathrm{d} x+\lambda \int_{\Omega} g \delta(\phi)(|\nabla \phi|) \mathrm{d} x \\ & +\alpha \int_{\Omega} g H(-\phi) \mathrm{d} x, \end{aligned} $

(2)

式中, g(x, y)表示边缘停止函数, 可用如下公式表示:

$ g(x, y) = \frac{1}{{1 + |\nabla {G_\delta } \times I(x, y){|^2}}} , $

(3)

式中, G_δ是参数δ的高斯核; I(x, y)为图像函数。

式(2)中, 能量泛函第1项为惩罚能量项, 也称之为距离正则化项。该项是DRLSE模型区别于传统水平集的主要内容。其中μ、λ > 0和α∈R分别为各能量项的常系数, P(|▽Φ|)为双阱势函数, 如下所示:

$ p(|\nabla \phi|)= \begin{cases}\frac{1}{(2 \pi)^2}(1-\cos (2 \pi|\nabla \phi|)), & \text { 如果 }|\nabla \phi| \leqslant 1 \\ \frac{1}{2}(|\nabla \phi|-1)^2, & \text { 如果 }|\nabla \phi| \geqslant 1 .\end{cases} $

(4)

能量泛函第2项为长度能量项, 其中δ(ּ)为Dirac函数, 其计算在水平集演化过程中的轮廓线积分。长度能量项可以光滑水平集演化过程中的主动轮廓, 减少及避免了零水平集尖点、角点和奇异点的产生。当长度能量项达到最小时, 运动演化的曲线会稳定停留在地物边界附近。

能量泛函第3项为面积能量项, 也称为区域能量项。其中H(ּ)为Heaviside函数。该项计算主动轮廓内的加权面积。当g = 1时, 该项为演化的零水平集内部面积, 该项的主要目的是零水平集距离目标边界较远时, 可以加快演化速度^[28]。

由于不存在理想化的Dirac函数和Heaviside函数, 因此通常使用如下公式进行表示^[29-30]:

$ \delta(x)=\left\{\begin{array}{ll} \frac{1}{(2 \pi)^2}(1-\cos (2 \pi x)), & |x| \leqslant \varepsilon \\ 0, & |x|>\varepsilon \end{array}, \right. $

(5)

$ H(x)= \begin{cases}\frac{1}{2}\left(1+\frac{x}{\varepsilon}+\frac{1}{\pi}\left(\frac{\pi x}{\varepsilon}\right)\right) & |x| \leqslant \varepsilon \\ 1, & |x|>\varepsilon \\ 0, & |x|>-\varepsilon\end{cases}, $

(6)

式中, ε为光滑参数, 其值一般为1.5。

使用DRLSE模型提取各个地物的轮廓, 由于引入了惩罚能量项, 有效避免了主动轮廓在演化过程中的不规则运动, 同时该模型不需进行重复初始化, 大大提高了图像分割效率^[31]。为方便水平集方法在计算机中具体实现, 使用的数值求解方法为有限差分法。在有限差分法实现数值求解过程中, 时间步长Δt应满足CFL^[32](Courant Friedrichs Lewy)收敛条件, 以保证函数的收敛性。Toure等人在使用DRLSE模型进行海岸线提取时, 将参数设置为μ=0.2, λ=5, α=3, Δt=1^[33]。考虑到两者研究背景一致性, 本文使用相同参数进行实验。观察DRLSE模型主动轮廓的运动演化过程, 当DRLSE主动轮廓运动到指定地物边界, 且在地物边缘附近无明显移动, 则停止迭代。由此可得到水体和地物之间的边缘轮廓线。相较于传统水平集方法, DRLSE模型可以随着目标的演化自动进行合并或者分裂, 并且其演化速度也得到很大提高。此外惩罚能量项的增加, 解决了DRLSE模型重新初始化的问题。

2.2 水边线信息精确提取

基于DRLSE模型提取出的轮廓中, 不但包含目标水边线, 还受噪声点轮廓、货轮轮廓、内陆湖泊轮廓以及沿海养殖场轮廓等影响。在DRLSE模型中, 水体的初始矩形轮廓是DRLSE模型的输入条件。为此, 本文对DRLSE模型的初始矩形轮廓的作用进行拓展, 一方面利用最小初始矩形轮廓周长可以删除噪声轮廓与周长数值较小的非水体轮廓, 加快解算速度; 另一方面根据该初始矩形轮廓的几何中心与地物轮廓面域之间的关系, 可以在诸多轮廓中选取出目标水体轮廓。

(1) 基于初始轮廓的几何中心的筛选

若矩形轮廓的几何中心X_j位于地物轮廓C_i构成的面域外, 表明该轮廓并非水体轮廓, 剔除C_i, 反之, 则保留此轮廓。

(2) 基于初始轮廓的周长轮廓的筛选

首先需要剔除噪声点轮廓和小周长轮廓。取矩形轮廓的周长最小值C_min=min{C}, 将其作为轮廓剔除的标准。若第i个边缘轮廓的周长小于矩形轮廓最小值, 则剔除该轮廓; 反之, 则保留该轮廓。一般来说, 使用轮廓的几何中心这一属性即可完成筛选过程, 但是使用选定的影像数据试验时, 结果显示该算法的运行效率低、时间长。因此, 先使用初始轮廓的轮廓周长筛选, 再基于初始轮廓的几何中心的筛选, 可有效降低数据冗余度, 提高算法运行速度。

经过以上的判断筛选, 将得到的地物轮廓与影像的外框矩形轮廓相重叠的地物轮廓线剔除, 进而得到最终的目标瞬时水边线。

3 实验结果与分析 3.1 瞬时水边线提取

为了验证顾及轮廓信息和DRLSE模型的遥感影像瞬时水边线快速提取方法效果, 采用研究区预处理后的GF-2数据进行实验分析。实验硬件环境: CPU为Inter(R)Core(TM)i5-7300HQ, 主频为2.50 GHz, 内存为16 GB; 软件环境为ENVI 5.3, 集成开发工具Visual Studio(2019), 以及在此基础上开发的DRLSE模块。

(1) 基于DRLSE模型水边线粗提取

利用归一化差异水体指数(NDWI)增强水体和其他地物的对比度, 得到归一化差异水体指数效果图如图 2(a)。DRLSE模型初始轮廓需要手动选取, 为了进行多片不连通海域瞬时水边线的提取, 本文开发的DRLSE模型可以同时进行多个不连通海域的初始矩形轮廓选取。手动选取DRLSE模型的初始矩形轮廓如图 2(b)所示, 该轮廓在影像中使用黑色矩形表示。选取方式为使用矩形在水体内部任意处进行标定即可。该遥感影像中有两块彼此隔离的水体面域, 使用初始矩形轮廓分别标定, 记录初始矩形轮廓的中心点坐标, 可以确定水体面域位置。同时这两块初始矩形轮廓属性信息可用作后续其他轮廓筛选的条件。

DRLSE模型的运行结果如图 2(c)所示, 得到水边线粗提取结果如图 3(a)所示, 可以看出, 图中不但包含目标瞬时水边线信息, 还包括噪声点如轮船等、海陆分割线轮廓、岛屿轮廓、内陆湖泊轮廓以及沿海养殖场轮廓。

(2) 顾及轮廓信息的瞬时水边线精确提取

DRLSE模型提取过程中的轮廓部分属性如表 1所示, 该轮廓属性包括了轮廓几何中心像元点坐标和轮廓周长, 其中轮廓周长为轮廓像元点个数。初始矩形轮廓1和轮廓2分别对应图 2(b)中上、下两个轮廓, 沿海养殖场轮廓和岛屿轮廓是算法运行之后产生的部分水体轮廓。初始矩形轮廓的周长始终小于水体真实轮廓周长, 因此小于该周长的噪声点轮廓和部分岛屿轮廓可有效剔除。养殖场等零星地物轮廓可使用初始矩形轮廓的几何中心像元位置进行判断, 若初始矩形轮廓的几何中心位于该轮廓包围的面域内, 可确定该轮廓为水体轮廓, 若初始矩形轮廓的几何中心不在轮廓包围的面域内, 则该轮廓为非水体轮廓进行剔除。通过上述的筛选, 得到最终的水边线信息如图 3(b)。

3.2 结果评价

为了验证本文方法的有效性, 选取泛洪算法、Canny算子和CV模型三种方法和本方法对比方法进行比较分析, 如图 4所示。研究区中涉及四类典型海岸带, 图 5(a)、图 5(f)、图 5(k)和图 5(p)框选区域分别对应人工岸线、砂质岸线、基岩岸线和淤泥质岸线, 其后分别对应4种方法的瞬时水边线提取结果。将目视解译提取出的水边线作为基准线, 对水边线提取结果进行定性分析。

对于边界清晰、形状规则的人工海岸区域[图 5(b)—图 5(e)], 4种方法均能够以较高的精度提取人工岸线区域的水边线, 误差在1~2个像元。对于亮度较强的砂质海岸[图 5(g)—图 5(j)], 除Canny算子外, 泛洪算法和CV模型与本文方法所提取的水边线与基准线较为接近, 误差在2~3个像元, 提取出的水边线结果符合一般误差要求。对于边界清晰度较低、岸线曲折的基岩海岸[图 5(l)—图 5(o)], 泛洪算法和Canny算子不能较为准确地提取出水边线, 而CV模型和本文方法仍能保证误差在两个像元以内浮动。对于边界不清晰, 且与水体颜色相近的淤泥质海岸[图 5(q)—图 5(t)], 其他3种对比方法提取出的水边线均较差, 而本文方法仍然能较准确地提取淤泥质瞬时水边线。

此外, 使用Canny算子获取到的研究区地物边缘轮廓数量大(总共2 049个轮廓), 后期需要人工从诸多边缘信息中目视判读出海岸带瞬时水边线位置, 之后使用矢量编辑器对其进行删除操作, 人工作业量大, 且提取出的地物轮廓边缘非连续如图 4(b), 水边线提取效率低; 而泛洪算法和CV模型只能处理一片连通水域。采用本文提出的瞬时水边线快速提取方法, 一方面解决了边缘提取算法中边缘非连续性的问题, 另一方面集成了地物轮廓边缘信息的提取以及轮廓筛选的功能, 可实现对指定瞬时水边线边缘轮廓的快速提取, 还可实现多水域的提取。

最后, 在现有运行平台下, 记录各方法模型对应研究区一景GF-2瞬时水边线提取时间, 并将提取的水边线按300 m间距进行抽样, 计算这些点到基准线的欧式距离, 采用均方根误差对提取结果位置精度进行定量评价, 如表 2所示。

从4种海岸水边线提取结果来看, 在人工海岸上, 本文和其他3种方法提取结果的均方根误差均小于10 m, 说明这几种方法均对人工海岸水边线提取效果较好。针对砂质海岸, Canny算子提取结果的均方根误差为133 m, 效果较差, 而泛洪算法、CV模型和本文方法分别为13 m、5 m和4 m。对于基岩海岸, 泛洪算法和Canny算子的均方根误差较大, 这是因为水和岩石混合、边界模糊, 而CV模型和本文方法的均方根误差较小, 均为5 m。对于淤泥质海岸, 其他三种方法提取结果均方根误差均较大, 但本文方法的均方根误差相对较小, 均方根误差为31 m, 总体来说, 4种方法对淤泥质海岸瞬时水边线提取结果均较差。从运行时间上来看, CV模型运行的时间最长, 原因是该模型需要不断地进行初始化。而泛洪模型、Canny算子和本文方法运行时间较短。由提取精度、运行时间以及后期处理的工作量上可以看出, 本文方法提取水边线能在保证精度的基础上有更高的效率。

4 结论和讨论

针对经典算法用于高分影像瞬时水边线提取结果不连续、效率不高和无法同时提取多片水域等问题, 提出了顾及轮廓信息和DRLSE模型的遥感影像瞬时水边线快速提取方法, 并以福建泉州一景高分2号4 m多光谱数据为例进行了验证。结果表明, 将初始轮廓信息(周长和中心位置)引入DRLSE中, 可有效消除噪声点轮廓等对瞬时水边线提取结果的影响, 效率更高, 且保证了瞬时水边线提取结果的连续性。目视解译的结果表明, 采用本方法提取结果的精度要优于泛洪算法、Canny算子, 而CV模型虽然在精度上有较好的表现, 但是耗时较长, 在本文选取的研究区实验中运行时间为28 801.4 s, 而本文方法只耗费了7.31 s, 在效率上远远高于CV模型。

尽管本方法对于给定的研究区中人工岸线、砂质岸线和基岩岸线瞬时水边线提取结果较好, 但在淤泥质岸线区域的提取精度仍有待提高。在其他地区的研究发现, 当影像空间分辨率进一步提高时(如优于1 m时), 影像中的白浪或泡沫等对瞬时水边线提取结果影响较大, 导致提取的瞬时水边线偏差较大, 因此对于更高空间分辨率的遥感数据, 应探讨其他白浪去除方法, 以获取较好的提取结果。此外, 本文给出的是大陆海岸线瞬时水边线提取方法, 对于岛屿海岸线瞬时水边线提取, 可以对本方法和程序进行修改, 以适应岛屿海岸线瞬时水边线提取工作。