潮滩通常指平均大潮高潮线和平均大潮低潮线之间的地带, 也称潮间带^[1], 本文研究的潮滩是复杂、敏感、动态不稳定、可进入条件差、冲淤多变的淤泥质潮滩。潮滩水边线是指潮汐波动起伏下海陆的瞬时交界线^[2], 用于反演潮滩地形和进行潮滩冲淤定量计算, 是提取海岸线的关键步骤之一。

常见的水边线提取技术包括阈值法^[3-4]、边缘检测算子法^[5]、数学形态学法^[6]、面向对象法^[7]等。遥感图像上, 淤泥质潮滩由于特殊的物质组成、平坦的地形和悬浮泥沙的影响, 表现出水沙特征区别度小、边界模糊的特点, 水边线提取难度极大。相关研究表明^[8-10], 小波变换能够有效检测图像边缘, 有较好的边缘定位性质^[11], 可以在保证图像边缘细节提取精度情况下较好地抑制噪声。故本文提出一种基于à trous小波分解与自适应阈值自动提取相结合的方法, 提取淤泥质潮滩水边线。这种算法既可以保证多尺度分析, 又避免了人工确定阈值的不确定性, 在水陆边界模糊条件下, 可以自动进行去噪、融合和边缘检测, 快速准确提取水边线, 为其他类似的水中泥沙含量高、水沙特征区别不明显、水陆边界模糊条件下的水边线提取问题的解决提供借鉴。

1 研究区与数据 1.1 研究区概况

本文的研究区是腰沙(图 1)。腰沙位于江苏省南通市东南侧海域, 长江口北岸^[12], 东安港口和遥望港口之间^[13]。腰沙形似半岛, 西侧与岸相连, 东北侧为三沙水道, 南侧为小庙洪水道^[13]。腰沙滩面平坦, 表层物质以粉砂质黏土及细颗粒泥沙为主, 是典型的粉砂淤泥质潮滩, 一旦被浪掀动悬扬, 就不容易沉降, 长期处于悬浮状态^[13]。由于粉砂中含有水且含水量不均匀, 水体、含水淤泥及含水粉沙表现出相近的色泽, 以致水沙边界不清晰。

1.2 数据来源

本文使用的影像的成像时间是2014年5月8日, 数据来源于高分一号卫星, 成像地点位于苏北辐射沙洲南部南通市通州区东部岸外滩涂(121°29′~ 122°30′ E, 31°53′~32°09′ N)。此高分影像全色分辨率2 m, 多光谱分辨率8 m, 包括4个波段, 覆盖范围从0.45 μm到0.9 μm, 坐标系为WGS84, 地图投影为经纬度投影Geographic Lat/Lon。水陆边界十分模糊, 具有一定的典型性(图 1)。

2 研究方法与实现过程 2.1 研究方法

本文采用基于à trous小波分解与自适应阈值自动提取相结合的方法, 提取淤泥质潮滩水边线。

2.1.1 À trous小波分解

本研究设置离散化尺度参数为2, 分解层次为4。基于à trous小波的分解原理如下:假设原图像数据为C₀(k), 经过尺度函数φ(x)的一次滤波后所得数据为{C₁(k)}, 由{C₀(k)-C₁(k)}得到两尺度图像间的细节信息(小波面)^[2], 滤波后得到不同层次的图像信息:

$ \sum {\phi {C_j}(k) = \frac{1}{{{2^j}}}} \left\langle {f(x), \varphi (\frac{{x - k}}{{{2^j}}})} \right\rangle , $

(1)

该过程也可用卷积运算代替, 其公式为

$ {C_j}(k) = \sum\limits_1^l {h(l){C_{j - 1}}(k + {2^{j - 1}})l} , $

(2)

h(l)代表低通滤波算子, f(x)代表尺度函数, 二者的关系为:

$ \frac{1}{2}\phi \left( {\frac{x}{2}} \right) = \sum\limits_l^4 {h\left( l \right)\phi \left( {x - l} \right)} , $

(3)

得到小波面, 如下:

$ {W_j}\left( k \right) = \left\{ {{C_{j - 1}}\left( k \right)} \right\} - \left\{ {{C_j}\left( k \right)} \right\}, $

(4)

则小波函数的定义如下:

$ \frac{1}{2}\varphi \left( {\frac{x}{2}} \right) = \phi \left( x \right) - \frac{1}{2}\phi \left( {\frac{x}{2}} \right), $

(5)

À trous小波算法的重构公式为:

$ {C_0}\left( k \right) = {C_p} + \sum\limits_{j = 1}^p {{w_j}\left( k \right)} , $

(6)

À trous算法分解和重构的滤波器会随着尺度的变化做膨胀和压缩, 不需要抽取偶数样本, 也不需要插零, 变换后也不会出现信号或图像压缩得现象, 经过对图像分解和聚焦图像上的细节, 提取图像上的低频信息, 再融合具有纹理信息的高频部分, 可达到水边线提取的目的。

2.1.2 自适应阈值的自动提取

自动阈值确定在Matlab中编程实现, 通过基于小波分解得到每一层的的模值和相角值, 进行模极大值检测。将局部区域内的最大的模值对应的小波系数留存, 其余的点对应的小波系数赋予0值剔除, 本质上是一种实现小波消噪的阈值选取方案, 最终得到边缘图像p(x, y)。

2.1.3 技术路线

本文采用的基于à trous小波分解与自适应阈值自动提取相结合的方法在水陆模糊条件下提取水边线, 并和经典算子提取结果对比进行精度验证, 技术路线见图 2。

2.2 边缘检测的实现过程

1) 图像去噪

选择分解层次j=4, 在Matlab中对原图像进行à trous小波变换和多分辨率分析。得到低频和高频子图像, 突出提取图像的低频信息, 保留图像的纹理等高频信息, 达到图像去噪的目的。

2) 模值和梯度计算

根据à trous小波变换对图像进行多层分解信息, 求小波变换各层的模值$W_s^1f\left( {x, y} \right)$和$W_s^2f\left( {x, y} \right)$, 计算每一尺度下小波变换的模值${M_s}f\left( {x, y} \right)$和梯度方向${A_s}f\left( {x, y} \right)$, 幅角s在–π/2到π/2,

$ \begin{array}{l} {M_s}f\left( {x, y} \right) = \sqrt {{{\left| {\mathop \sum \limits_j W_s^jf\left( {x, y} \right)} \right|}^2}} , 1 \le j \le 4, \\ - {\rm{ \mathsf{ π} }}/2 \le s \le {\rm{ \mathsf{ π} }}/2, \end{array} $

(7)

$ {A_s}f\left( {x, y} \right) = {\rm{arctg}}\left[ {\frac{{W_s^yf(x, y)}}{{W_s^xf(x, y)}}} \right], $

(8)

3) 模极大值检测

检验梯度方向内局部区域的模极大值, 幅角范围为–π/2到π/2之间, 进行局部模极大值检测。沿着梯度的4个方向, 比较相邻点之间的模值, 将局部最大值留存用于下一流程, 局部内其余的点赋予0值, 作为非边缘点剔除, 得到可能的边缘图像p(x, y)。

4) 自适应阈值自动选择

在可能的边缘图像p(x, y)中, 将不间断的点构成为线。同时将每一条线中横坐标和纵坐标作为最大值和最小值, 围成矩形小分块。在矩形小分块中, 统计每个点的模极大值, 计算出均值, 这个值就是自适应阈值。如果小分块内的点的模值大于等于阈值, 将此点输出, 其余点赋予0值剔除。

5) 水边线检测

将不同层次上的边缘图像融合并输出, 实现对水边线的检测。

3 研究结果分析 3.1 本文小波变换多尺度检测结果分析

利用à trous自动阈值结合算法, 在Matlab中编程, 实现遥感图像水边线边缘的提取, 结果如图 3所示。一二层检测尺度小, 提取的水边线细节突出, 清晰连续, 定位准确, 偏移较小(图 3a, 图 3b)。三四层检测尺度大, 水边线轮廓光滑连续, 但水边线有向陆一侧移动, 部分信息丢失(图 3c, 图 3d)。

由于受不同尺度平滑效果的影响, 水陆过渡较为模糊的区域, 也就是图中中间部分, 曲折较多的部分, 水边线检测结果有较大差异, 作进一步具体分析(图 4)。如图 4所示, 水边线在一层检测时比较清晰, 但呈点状分布, 这是因为本文使用的是样条小波变换, 是基于离散点的变换, 在尺度较小时, 会显示出离散点的特征。二层检测时水边线能被清晰地检测出来, 效果佳。但在三层检测和四层检测时, 水边线被进一步平滑, 水边线附近的边缘信息和阴影信息融合, 水边线变粗, 细节体现不明显。

3.2 小波变换结果与其他算子提取结果定性对比分析

水边线提取的核心是要保证精度。在对高分遥感图像进行预处理之后, 又分别采用Sobel算子, Log算子, Canny算子对图像进行边缘检测(图 5), 通过定性对比各种经典算子处理效果, 对本文探究的方法精度进行验证。

由于该图像中部水边线沟壑和曲折丰富, 具有典型的代表性, 所以本文对图像中间部分进行重点比较, 对影像中部进行的检测处理结果如图 6所示。

由图 6看出, Sobel算子提取的水边线断断续续呈点状分布, 这是因为Sobel是一个离散的基于一阶导数的局域窗口梯度算子, 且Sobel算子没有严格的模拟人的视觉生理特征。传统的Log算子采用中心对称函数计算, 具有各向同性的特点, 对不同方向的边缘检测效果不佳^[14]。Canny算子边缘检测原理类似于Log算子, 属于先平滑后求导的方法, 且Canny算子容易把噪点误判为边界。在水陆边界模糊, 边界噪声突出条件下, Canny算子水边线检测效果十分不理想(图 6)。虽然Sobel算子检测结果与本文采用的à trous自动阈值结合法的检测效果相近, 但Sobel算子得到的水边线边缘较粗, 这是因为Sobel算子在去除噪声的同时, 将边缘的基本轮廓保留下来, 造成轮廓线和真正水边线的重叠, 同时无法剔除噪声点。故本文算法检测出的水边线结果更加理想。

3.3 小波变换结果与其他算子提取结果定量对比分析

对水边线提取结果图像进行边缘清晰度(image edge definition, IED)检测, IED越大, 图像越清晰, IED计算公式如下, 其中P(i, j)表示图像坐标为(i, j)的像素, N, M对应图像矩阵的行数和列数^[15]:

$ {\rm{IED}} = 10\left| {\frac{{\sum\limits_{\begin{array}{*{20}{c}} {2 \le i \le N - 1}\\ {2 \le j \le M - 1} \end{array}} {{p^2}\left( {i,j} \right)} }}{{\left( {N - 2} \right)\left( {M - 2} \right)}}} \right|, $

(9)

IED计算结果表明, Sobel算子为37.254 3, Log算子为28.775 2, Canny算子为27.265 7, à trous自动阈值结合法为38.997 4。通过定量计算, 表明à trous自动阈值结合法提取的水边线边缘清晰度数值比Sobel算子、Log算子、Canny算子提取结果的边界清晰度数值高, 证明采用à trous自动阈值结合法提取效果比较理想。

4 结论

1) 水陆模糊条件下, 采用基于à trous小波函数和自适应阈值自动提取相结合的方法对高分影像进行多层次分析, 既不会出现信号压缩, 又有利于获取图像的细节特征, 并且避免了人为确定阈值带来的不确定性误差, 可以快速准确提取水边线, 具有一定的方法优势。

2) 实验结果表明对于本文方法提取的水边线结果无论在视觉效果定性对比上, 还是在图像边缘清晰度定量指标计算上, 均优于Sobel算子、Canny算子和Log算子方法提取的结果, 表明此方法可以为其他类似的水中泥沙含量高、水沙特征区别不明显、水陆边界模糊条件下的水边线提取困难的问题的解决提供方法学上的借鉴。