基于密度聚类与灰度变换的NSST 域声呐图像去噪

2022-09-02刘光宇曾志勇曹禹赵恩铭邢传玺

湖南大学学报（自然科学版） 2022年8期

刘光宇，曾志勇†，曹禹，赵恩铭，邢传玺

（1.大理大学工程学院，云南大理 671003；2.哈尔滨工程大学物理与光电工程学院，黑龙江哈尔滨 150001；3.云南民族大学电气信息工程学院，云南昆明 650031）

随着经济社会的发展和相关技术的成熟，人们探索海洋的欲望越来越强烈.声呐探测技术作为海洋探索的基本手段，在海底地形测绘、水下目标物体探测以及生物种群监测方面都发挥着重要作用［1］.由于受海风、洋流、水温、杂质以及成像设备等影响，使用声呐探测技术得到的图像中通常包含各种类型的噪声，主要表现为颗粒状的斑点，尤其在浅海区域，该现象更为严重［2］.斑点噪声的存在导致图像质量严重下降，对后续的图像处理产生不利影响［3］.因此，如何有效去除图像中的噪声成了人们研究的热点.与光学图像一样，声呐图像中常见的去噪方法主要包括空间域去噪法（如均值滤波、中值滤波和维纳滤波等）和变换域去噪法（如小波变换、脊波变换和轮廓波变换等）［4］.其中变换域去噪法因其复杂程度低，且具有良好的边缘保持能力而更受相关研究者们的关注［5］.

近年来，剪切波作为变换域中一种新的多尺度几何分析工具在各种图像去噪方法中脱颖而出，以其结构简单，对多维数据的各向异性特征近乎最优的稀疏表示而深受人们的青睐［6］.由于剪切波变换后图像存在一定程度的伪吉布斯现象［7］，因此，Easley 等人在其基础上进行改进，提出了非下采样剪切波变换［8］.2020年，Routray等把非下采样剪切波和双边滤波进行结合，构建了一种新的图像去噪模型［9］；同年，Morteza 等提出了基于TLS（t-locations scale）和非下采样剪切波的噪声滤波方法［10］.这些方法在保持剪切波原有优点的同时，提高了图像的峰值信噪比和结构相似度，降低了均方误差，克服了图像中出现的伪吉布斯现象，使图像的质量得到进一步提升.

由于声呐图像的噪声特性与光学图像、红外图像等存在差异，直接运用光学图像和红外图像等的去噪方法将难以使声呐图像达到理想的去噪效果.斑点噪声在声呐图像中的分布呈现密度不一的特点，主要集中在混响区域和目标区域，在阴影区域密度较小［11］，与常见的噪声类型相比，其复杂性更强.聚类算法在图像处理中已有成功应用.其中，DBSCAN（Density-based Spatial Clustering of Applications with Noise）算法可根据数据的分布特性自主计算邻域值，利用数据内部的聚集性和连通性可实现更为准确和高效的数据分类［12］.

本文把非下采样剪切波变换与DBSCAN 密度聚类相结合，利用灰度变换对声呐图像进行处理，构建一种新的图像去噪模型.该模型在去除斑点噪声和缓解伪吉布斯现象的同时，进一步保持了图像的边缘细节特征，对提高图像峰值信噪比、结构相似度和降低均方误差等拥有明显的效果.

1 基本理论

1.1 噪声模型

根据噪声与信号的关系，图像的噪声类型可分为加性噪声和乘性噪声［13］.声呐图像中常见的噪声类型为乘性斑点噪声，其模型可描述为：

式中：f（x，y）为含有噪声的图像；n（x，y）为原始图像；q（x，y）为噪声信号.乘性斑点噪声通常服从瑞利分布，与光学图像中常见的加性高斯噪声相比，具有起伏剧烈、均匀度较低以及存在较大不确定性等特点，处理起来更为困难［14］.加性高斯噪声和乘性斑点噪声在图像中的分布如图1所示.

图1 噪声分布模型图Fig.1 Model diagram of noise distribution

1.2 剪切波变换

剪切波变换（Shearlet Transform，ST）是在合成小波变换基础上发展的数学算法，可通过对基函数实行伸缩、平移和剪切等操作生成所需的剪切波函数，是一种接近最优的多维函数稀疏表示方法［15］.

当维数n=2时，剪切波的仿射系统为：

式中：ψ∈L2（R2）；A、B均为二维可逆矩阵，且|detB|=1；j和l分别为尺度参数和几何参数；k为平移参数.如果AAB（ψ）满足紧支撑框架，即

则AAB（ψ）系统中的元素ψ就称为合成小波，其中膨胀矩阵Aj代表尺度变换，矩阵Bl代表几何变换.通过上述仿射变换，系统可以构造各个尺度和各个方向上的紧支撑框架［16］.令式（2）中的A为各向异性的膨胀矩阵，B为剪切矩阵，通常情况下：

当a=4，s=1时，

此时的合成小波被称为剪切波［17］.剪切波是一个带有尺度、方向和位置的函数集合.在频域中，不同尺度的剪切波在其频域支撑下，表现为近似22j×2j大小的梯形形状.剪切波基函数的支撑区域可以随着尺度的长宽比发生变化，能够最优地表示图像的内在几何结构，其频域剖分和频域支撑示意图如图2 所示.图2（a）中的ξ1和ξ2为剪切波频域剖分的尺度参数.

图2 剪切波频域剖分及其频域支撑示意图Fig.2 Shearlet frequency domain subdivision diagram and schematic diagram of frequency domain support

1.3 非下采样剪切波变换

传统剪切波变换在多尺度分解过程中存在下采样操作，容易导致细节信息丢失，重构后的图像极易产生伪吉布斯现象.非下采样剪切波变换（Nonsubsampled Shearlet Transform，NSST）的出现解决了上述问题，在保持传统剪切波良好的局部性、方向敏感性、最佳稀疏性等特性的同时，有效缓解了频谱混叠现象［18］.

非下采样剪切波变换的离散化过程分为多尺度分解和方向局部化两部分［19］：

1）多尺度分解：利用非下采样金字塔滤波器（Non-subsampled Pyramid Filter，NSPF）将图像分解为k+1 个与原图像大小相同的子带图像，其中包括1个低频子带图像和k个高频子带图像.

2）方向局部化：标准的剪切波变换中使用的剪切滤波器（Shear Filter，SF）是在伪极化坐标系中通过窗函数的平移操作实现的，此过程中包含下采样操作，不具备平移不变性.而非下采样剪切波变换将标准的剪切滤波器从伪极化网络系统映射回到笛卡尔坐标系统，摒弃下采样操作，实现了平移不变性.

非下采样剪切波的具体实现过程：首先，采用非下采样金字塔滤波器将图像分解为低频子带图像和高频子带图像.其次，对高频子带图像构建Meyer窗，进行方向局部化，获得不同的方向子带.最后，对每一个方向子带进行傅立叶逆变换，以获得非下采样剪切波系数.NSST 二级分解示意图如图3所示.

图3 NSST二级分解示意图Fig.3 Schematic diagram of NSST level two decomposition

1.4 密度聚类

DBSCAN 密度聚类作为聚类中常见的算法，因其具有识别任意形状的簇、不需要事先设置簇的数量以及能够识别噪声点等优点，受到人们的青睐［20］.它的基本思想是定义一个簇为所有紧密连接的数据点的最大集合，在以某个数据点为圆心，以ε为半径的区域内，DBSCAN 通过区域查询，判断该区域包含的数据点个数是否大于设定的密度阈值δ，如果大于该数值，则创建一个新的簇，然后把相关联的簇进行合并.如果某个对象不在任何簇中，则该对象被视为噪声点，从而达到分类的目的［21］.基于参数ε和δ，给出数据集D=｛x1，x2，…，xm｝，其中xi∈D，xj∈D.DBSCAN密度聚类算法有以下基本概念［22］：

1）核心对象：若以xj为圆心，ε为半径的区域内至少包含δ个样本，则称xj为核心对象；

2）密度直达：若xj位于以xi为圆心，ε为半径的区域内，且xi为核心对象，则称xj由xi密度直达；

3）密度可达：若存在样本序列K1，K2，…，Kn，其中K1=xi，Kn=xj且Ki+1由Ki密度直达，则称xj由xi密度可达；

4）密度相连：若存在xh均可使xi和xj由xh密度可达，则称xi和xj密度相连.

当δ=3 时，DBSCAN 密度聚类示意图如图4 所示.其中，虚线圈内的区域是以ε为半径的区域，x1为核心对象，x2由x1密度直达，x3由x1密度可达，x3与x4密度相连.

图4 DBSCAN密度聚类示意图Fig.4 Schematic diagram of DBSCAN density clustering

1.5 灰度变换

灰度变换是指将原图像中的灰度值映射到另一灰度值，其中比例线性变换、分段线性变换是常见的映射函数［23］.由于分段线性变换在一定的灰度值区间内具有不同的变换函数，比单一的线性变换函数更为灵活和更具针对性，可根据需求拉伸图像每一部分的灰度值，达到增强图像质量，提高对比度的效果.分段线性变换函数示意图如图5所示.

图5 分段线性变换函数示意图Fig.5 Schematic diagram of piecewise linear transformation function

2 NSST域声呐图像去噪

本文结合非下采样剪切波变换、DBSCAN 密度聚类和灰度变换对声呐图像进行去噪，具体去噪步骤如下：

1）利用非下采样金字塔滤波器对含噪图像进行4 层分解，分解为1 层低频子带图像和4 层高频子带图像，同时通过剪切滤波器把4 层高频子带图像分别分为23、23、24和24个方向.含噪图像及其高低频子带图像如图6所示.

图6 含噪图像及其高低频子带图像Fig.6 Noisy image and its high and low frequency subband images

从图6 中可以看出，含噪图像经过分解后，原始信息被分解到不同层次.其中，低频子带图像包含了原始信息的主要能量，高频子带图像包含了大量细节信息和噪声信号.随着分解层数的增加，细节信息和噪声信号越来越集中.因此，将更高层次的子带图像分成多个方向，以精准捕获每个方向上的细节特征.

2）根据斑点噪声在声呐图像中的分布特性，利用DBSCAN 密度聚类对每一层的高频子带图像进行聚类处理，δ的取值通常为N+1，其中N为数据的维度.通过调节ε值，将信号划分为细节信号和噪声信号，对噪声信号进行零值处理，保留细节信号.

3）低频子带图像中通常不包含噪声干扰，因此，只需对其对比度进行调整.利用灰度变换中的分段函数对低频子带图像进行灰度拉伸，进一步提升图像的整体视觉效果.

4）利用非下采样剪切波逆变换对处理后的高频子带图像和低频子带图像进行重构，得到恢复后的图像.

整体去噪流程图如图7所示.

图7 整体去噪流程图Fig.7 Overall denoising flow chart

3 实验结果与分析

实验1选取大小为256×256的沉船图像［（图8（a）］作为实验对象，对其添加方差为0.5 的斑点噪声，运用本文所提方法对其进行处理.得到重构后的图像与变换域中的小波变换（Wavelet Transform，WT）、脊波变换（Ridgelet Transform，RT）、轮廓波变换（Contourlet Transform，CT）、剪切波变换（Shearlet Transform，ST）和非下采样剪切波变换（Non-subsampled Shearlet Transform，NSST）去噪法的去噪结果进行对比.对比结果如图8所示.

图8 沉船图像Fig.8 Images of shipwreck

从图8（c）（d）（e）可以看出，通过小波变换、脊波变换和轮廓波变换去噪后，图像中的斑点噪声得到一定程度的去除，但同时图像整体出现了模糊，细节信息丢失严重.在图8（f）中，通过剪切波变换去噪后，图像中的噪声明显减少，但其边缘保持能力不足，影响了图像的质量.在图8（g）中，通过非下采样剪切波变换处理后，图像中的伪吉布斯现象得到明显改善，边缘保持能力相比剪切波去噪法有了较大的提升.在图8（h）中，采用DBSCAN 密度聚类、灰度变换与非下采样剪切波变换相结合的方法进行处理后，图像中的噪声得到有效滤除，边缘轮廓和细节特征得到准确捕捉，失真程度明显小于小波变换（WT）、脊波变换（RT）、轮廓波变换（CT）、剪切波变换（ST）和非下采样剪切波变换（NSST）去噪法，从直观视觉的角度证明了该方法的优越性.

为进一步客观具体地评估所提出的方法，对图8（b）中的含噪图像（Noisy Image，NI）重复实验20次，计算去噪后图像的质量评价指标［24］，即均方误差（Mean Squared Error，MSE）、信噪比（Signal to Noise Ratio，SNR）、峰值信噪比（Peak Signal to Noise Ratio，PSNR）和结构相识度（Structural Similarity，SSIM）的平均值，并通过数据对图像质量进行判断.MSE 值、SNR值、PSNR值和SSIM值的计算公式分别如下：

式中：QMSE表示MSE 的值；m×n表示图像的大小；I（i，j）为原始图像；R（i，j）为重构后的图像.MSE值越小表明图像质量越好.

式中：QSNR表示SNR的值.

式中：QPSNR表示PSNR 的值；Lmax表示原始图像I（i，j）的最大灰度值.与SNR 值一样，PSNR 值越大表明图像失真程度越小，图像质量越好.

式中：QSSIM表示SSIM 的值；ux、uy分别为图像x和y的均值；σx、σy分别为图像x和y的标准差.SSIM 是衡量两幅图像相似度的指标，其值越大表明两幅图像越相似.

经过多次实验后，沉船图像的均方误差（MSE）、信噪比（SNR）、峰值信噪比（PSNR）和结构相识度（SSIM）的平均值如表1所示.

表1 沉船图像质量评价指标Tab.1 Quality evaluation index of shipwreck image

从表1 中的数据可以明显看出，采用小波变换（WT）、脊波变换（RT）和轮廓波变换（CT）去噪后，图像的SNR 值和PSNR 值有一定程度的提升，但SSIM值效果不理想，图像失真较为严重.采用剪切波变换（ST）去噪后，SNR 值、PSNR 值和SSIM 值分别提升了232.87%、59.85%和60.73%，MSE 值下降了86.46%.采用非下采样剪切波变换（NSST）去噪后，SNR 值、PSNR 值和SSIM 值分别提升了253.05%、65.04%和65.39%，MSE 值下降了88.40%，去噪效果好于剪切波变换.运用密度聚类、灰度变换和非下采样剪切波变换相结合的方法处理后，图像的SNR 值从3.702 8 dB 提升至13.855 7 dB，提升了274.20%；PSNR 值从14.406 8 dB 提升至24.559 7 dB，提升了70.47%；SSIM 值从0.328 2 提升至0.565 2，提升了72.21%；MSE 值从0.036 2 下降至0.003 5，下降了90.33%.本文方法的SNR 值、PSNR 值和SSIM 值提升的幅度以及MSE值下降的幅度均大于小波变换（WT）、脊波变换（RT）、轮廓波变换（CT）、剪切波变换（ST）和非下采样剪切波变换（NSST）去噪法，从客观数据表明，本文方法在上述6 种变换域去噪法中去噪效果最好，失真程度最小，边缘保持能力最强.

实验2为进一步验证本文方法在不同噪声密度条件下的客观性和优越性，利用公共图像数据集进行实验，并与现有文献方法进行对比.将图9 中的公共图像数据集添加方差分别为0.3、0.4、0.5 和0.6的斑点噪声，运用文献［25］的BM3D 算法、文献［26］的PGPD 算法、文献［27］的EPLL-GMM 算法、文献［28］的ACVA 算法和本文所提方法进行去噪处理，得到去噪后图像的MSE、SNR、PSNR 和SSIM 的平均值如表2所示.

图9 公共图像数据集Fig.9 Public image dataset

表2 公共图像数据集质量评价指标Tab.2 Quality evaluation index of public image dataset

从表2 中的数据可知，本文所提方法在噪声方差分别为0.3、0.4、0.5和0.6的情况下，MSE 值分别下降了89.76%、91.21%、91.40% 和91.64%；SNR 值、PSNR 值和SSIM 值分别提升了164.34%、197.74%、223.60%、247.41%，86.42%、97.90%、103.25%、110.20%和293.32%、345.45%、387.73%、422.05%.各项指标的变化幅度随噪声方差的增大而逐渐提升，并且其值均优于文献［25-28］的去噪方法，表明本文方法在去除声呐图像斑点噪声中拥有明显优势.通过折线图进一步展示各项评价指标的变化趋势，趋势图如图10所示.

图10 图像质量评价指标趋势图Fig.10 Trend charts of image quality evaluation indexes

图10 中，随着噪声方差的增大，BM3D、PGPD、EPLL-GMM、ACVA 和本文所提方法的MSE 值呈上升趋势，SNR 值、PSNR 值和SSIM 值呈下降趋势.其中，BM3D 的MSE 值在5 种去噪方法中保持最大，SNR值和PSNR 值保持最小，表明其去噪效果不太理想；ACVA 的SNR 值、PSNR 值和SSIM 值均大于PGPD 和EPLL-GMM 算法，表明其去噪效果好于PGPD 和EPLL-GMM 算法；本文去噪法的各项指标变化平稳，始终优于其他4 种去噪方法，并且其优势随噪声方差的增大而更加明显，证明了本文方法更适用于具有高密度噪声的声呐图像去噪.不同方法的去噪效果如图11所示.

图11 不同方法的去噪效果Fig.11 Denoising effect of different methods

在图11 中，含噪图像分别通过BM3D、PGPD、EPLL-GMM、ACVA 和本文方法进行去噪处理后，图像质量得到了不同程度的提升.其中BM3D 算法去除了图像中的部分斑点噪声，但效果还不够理想；PGPD 和EPLL-GMM 算法噪声去除明显，但细节信息丢失较为严重，影响了图像的质量；ACVA 算法的边缘保持能力得到较大的提升，但与本文方法相比还存在一定差距.

除了MSE、SNR、PSNR 和SSIM 等常用指标，时间复杂度也是评价图像去噪算法的重要标准.运用搭载Intel Core i5-6200 处理器和4 GB 运行内存的计算机在MATLAB R2019b 平台上进行多次试验，得到不同算法的平均运行时间如表3所示.