0 引　言

在声呐信号采集和传输的过程中，数字信号经常被各种噪声所污染，且海洋环境复杂多变，混合的声音信号使声呐图像不同于光学图像。针对声呐图像分辨率低、背景噪声高、边缘对比度差及缺少有用信息等特性，结合光学图像的研究成果，提出了一些去噪增强理论和算法，包括线性与非线性函数^[1]、多尺度retinex^[2]、形态学方法^[3]等，由于声呐图像本质与光学图像有很大的差别，因此对声呐图像处理的应用受到限制，效果也并不理想。从20世纪末开始发展起来的多分辨率分析理论，例如傅里叶变换^[4]、小波变换^[5]、Ridgelet变换^[6]和Curvelet变换^[7]，为处理信号奇异性提供强大的分析工具。近年来，小波在声呐图像降噪和识别中取得了一定的成功，由于小波的局限性，并不能有效地解决线奇异性和有限的方向选择性。为了克服这些缺点，Donoho等^[8]提出一种新的图像表示方法Ridgelet，可以有效处理线或超平面奇点并以足够精细的比例以局部方式分析曲线。Da C等^[9]以Contourlet 变换为基础进行改进提出了具有多尺度、多方向和平移不变性的非下采样Contourlet变换（NSCT），使用非下采样拉普拉斯金字塔（LP）和方向滤波器组（DFB）构造分解结构。该过程可以避免下采样过程，从而使变换具有平移不变性，并进一步提高Contourlet变换在图像去噪领域的性能。基于传统贝叶斯原理的NSCT系数方差估计算法存在一个问题，即如何选择合适的原始图像先验概率系数模型和噪声系数^[10-11]。为了解决此问题，可以使用变贝叶斯方法。V-Bayes的基本思想^[12]是用更容易近似的分布近似真实的后验概率，通过最小化近似分布和实际后验概率分布之间的Kullback-Leibler Divergence（KL散度）^[13]差异来预测估计。

本文结合NSCT二元模型和V-Bayes估计，提出了一种图像去噪算法，对声呐图像进行了仿真分析，并与现有的经典方法进行了比较。结果表明，该算法可以有效去除图像中的噪声，获得最高的PSNR值和噪声抑制能力。

1 相关介绍 1.1 非下采样轮廓波NSCT

Contourtlet变换可以以最佳方式稀疏地表示图像，并且能有效捕获和刻画图像中弯曲和定向的几何结构，其结构中LP和DFB都进行了下采样操作，因此Contourtlet系数的冗余大大降低，以至于Contourlet 变换缺少平移不变性，使增强后的图像在不连续的点边缘产生伪“Gibbs”现象。为了解决上述问题，NSCT舍弃了Contourtlet变换中的下采样操作，将非下采样金字塔（NSP）和非下采样方向滤波器组（NSDFB）组合，在采样的比例上与原始图像比例相同。

NSP可确保NSCT的多尺度特性，NSP分解阶段生成2个子带：低频子带和高频子带，利用NSDFB对高频子带进行方向分析，从各个方向展示图像的轮廓细节和边缘特征。NSP对图像进行n级尺度分解，可得到n+1个子带图像。图1为2层NSCT变换的原理图。

1.2 双变量模型

小波变换系数存在一定程度的相关性，其中，双变量模型充分利用了图像小波系数的尺度相关性，使用非高斯密度函数对“父”和“子”系数的分布进行建模，并使用贝叶斯统计理论获得解析表达式最大后验估计。

NSCT分解相邻尺度之间的各级系数有很强的相关性，即NSCT变换域中的系数具有与小波系数相似的分布特性。因此，小波域中的双变量模型的处理思想可用于将其扩展到NSCT变换域中。

为了估计NSCT系数复方差 $ \sigma _{{n}}^2 $ ，可以使用贝叶斯估计方法，将图像观察信号表示为：

$ y = w + n, $

(1)

其中 $ w=[{w}_{1},{w}_{2}] $ ， $ {w}_{1} $ 为原始图像的NSCT系数，而 $ {w}_{2} $ 是在先前图中与 $ {w}_{1} $ 的父系数具有相同空间位置的NSCT系数； $ y=[{y}_{1},{y}_{2}] $ 和 $ n=[{n}_{1},{n}_{2}] $ 分别表示由观察图像的NSCT分解系数和噪声组成的向量。

在NSCT域，根据贝叶斯最大后验估计（MAP），从 $ y $ 获得 $ w $ 的估计 $ \widehat{w} $ ：

$ \hat w(y) = {\rm{argmax}}[{\rm{log}}({p_n}(y{\text{ - }}w)) + {\rm{log}}({p_n}(w))]。$

(2)

上面的方程式中，首先必须知道噪声的NSCT分解系数的概率分布。在此假设噪声密度函数 $ {p}_{n} $ 服从高斯分布，系数方差为 $ {\sigma }_{n}^{2} $ ：

$ {p_n}(n) = \frac{1}{{2{\text{π}} \sigma _n^2}}{\rm{exp}}\left( - \frac{{n_1^2 + n_2^2}}{{2\sigma _n^2}}\right) 。$

(3)

Sendur等提出双变量模型对小波上下级系数进行描述，NSCT内系数概率分布用联合概率密度函数 $ {p}_{w} $ 表示：

$ {p_w}(w) = \frac{3}{{2{\text{π}} {\sigma ^2}}}{\rm{exp}}( - \frac{{\sqrt 3 }}{\sigma }\sqrt {w_1^2 + w_2^2} ), $

(4)

式中： $ \sigma $ 表示每个NSCT系数的模型 $ {w}_{1} $ 和 $ {w}_{2} $ 的方差。将式(3)、式(4)代入式(2)中得NSCT系数 $ {w}_{1} $ 的MAP估计为：

$ {\hat w_1} = \frac{{\left(\sqrt {y_1^2 + y_2^2} - \dfrac{{\sqrt 3 \sigma _n^2}}{\sigma }\right)}}{{\sqrt {y_1^2 + y_2^2} }} + {y_1} 。$

(5)

式(5)是双变量模型中当前系数 $ {w}_{1} $ 的MAP估计，称为双变量萎缩函数。可以看出，如果要获得 $ {\widehat{w}}_{1} $ ，则必须首先估计噪声的NSCT系数方差 $ {\sigma }_{n}^{2} $ 。在本文中，使用变分贝叶斯估计(V-Bayes)方法来估计系数方差。

1.3 变分贝叶斯阈值估计

MAP估计方法不需要显示求出后的后验密度函数，只要得出最大概率密度的点值，属于点估计范畴，对模型参数没有给出一个很好的估计准则，而变分贝叶斯方法是对整个后验密度的求解，因此可以利用更多的信息.为了解决传统贝叶斯估计问题，变分贝叶斯方法提供了有效的解决方案。

对给定的原始图像w，噪声图像y，根据贝叶斯定理：

$ p\left( {\nabla \left. w \right|\nabla y} \right) = \frac{{p\left( {\nabla \left. y \right|\nabla w} \right)p\left( {\nabla w} \right)}}{{p\left( {\nabla y} \right)}}。$

(6)

其中： $\nabla y $ 为噪声图像的梯度， $\nabla w $ 为原始图像梯度，P(y|w)为似然密度函数；P(w)为原始图像的先验密度函数；P(y|w)为后验密度函数。

用易于处理的近似分布q( $\nabla w $ )近似真实的后验概率分布p( $\nabla w $ | $\nabla y $ )。KL散度可以测量2个分布之间的距离，因此可以通过最小化近似分布和实际后验概率分布之间的KL散度来实现估计：

$ \begin{split} &KL\left\{q(\nabla w,{\sigma }^{2})\Vert p(\nabla w|\nabla y)\right\}=\\ &\qquad \int q(\nabla w,{\sigma }^{2})\mathrm{ln}\dfrac{q(\nabla w,{\sigma }^{2})}{p(\nabla w|\nabla y)}{\rm{d}}\nabla w+\mathrm{ln}p(\nabla y)\geqslant 0，\end{split} $

(7)

当且仅当q（ $\nabla $ ω,σ2）= p（ $\nabla $ ω| $\nabla y $ ）时，以上等式成立。

$p(\nabla y) $ 在整个估计过程中始终是常数，因此可以定义成本函数 $ {C}_{KL} $ 以获得近似分布的最佳值：

$ \begin{split} {C}_{KL}=&KL\{q(\nabla w,{\sigma }^{2})\Vert p(\nabla w|\nabla y)\}-\mathrm{ln}p(\nabla y)=\\ & \int q(\nabla w,{\sigma }^{2})\mathrm{ln}\dfrac{q(\nabla w,{\sigma }^{2})}{p(\nabla w|\nabla y)}{\rm{d}}\nabla w+\\ & \int q(-{\sigma }^{2})\mathrm{ln}\dfrac{q(-{\sigma }^{2})}{p(-{\sigma }^{2})}{\rm{d}}(-{\sigma }^{2})。\end{split} $

(8)

贝叶斯变化估计通过迭代方法（即贝叶斯期望最大化定理）使成本函数最小，从而估计噪声函数。使用变分贝叶斯方法进行噪声参数估计，特点在于可以求出后验密度分布的逼近分布，存在闭合形式的解，可以更好地估计后续的高频子带噪声参数。

2 模糊对比度增强

NSCT分解后，其低频子带具有 $ L $ 个灰度级的 $ M\times N $ 图像，可以看作为一个模糊集，其每一个元素对应的模糊特征集合，记为模糊矩阵：

$ X = \mathop \bigcup \limits_{m = 1}^M \mathop \bigcup \limits_{n = 1}^N \frac{{{u_{mn}}}}{{{x_{mn}}}} ,$

(9)

模糊集的隶属度函数将低频带图像映射到模糊特征平面，从而使图像特征更容易识别和提取。采用一种改进的正余弦结合模糊隶属函数：

$ {u_{mn}} = \frac{1}{2}\left[{\rm{sin}}\left(\frac{{\text{π} ({x_{mn}} - {x_{{\rm{min}}}})}}{{2({x_{{\rm{max}}}} - {x_{{\rm{min}}}})}}\right) - {\rm{cos}}\left(\frac{{\text{π} ({x_{mn}} - {x_{{\rm{min}}}})}}{{2({x_{{\rm{max}}}} - {x_{{\rm{min}}}})}}\right) + 1\right] 。$

(10)

改进的隶属度函数曲线为“S”型，符合模糊理论对隶属度函数凹凸性的要求，且具有效果较好的去噪性。

运用模糊增强算子正弦函数的回归调用来修正隶属度：

$ \begin{split} u_{mn}^{'} =& {T^r}\left( {{u_{mn}}} \right) = T\left( {{T^{r - 1}}\left( {{u_{mn}}} \right)} \right) = \\ &\frac{{{\rm{sin}}(\text{π} ({u_{mn}}{{ - }}0.5)) + 1}}{2}。\end{split} $

(11)

式中，r= $1,2,3\cdots $ ， $ {T}^{r}\left({u}_{mn}\right) $ 为变换函数 $ T\left({u}_{mn}\right) $ 的第r次迭代，此变换以越渡点0.5为界限，根据非线性增强规则进行迭代变换。

将图像从模糊域映射到空间域，其逆变换为：

$ {x_{mn}} = {x_{{\rm{min}}}} + ({x_{{\rm{max}}}} - {x_{{\rm{min}}}})\frac{{{\rm{arcsin}}{{[1 - {{\left( {2{u_{mn}} - 1} \right)}^2}]}^{\sqrt r }}}}{{\text{π}}} 。$

(12)

2.1 具体步骤

NSCT具有多方向选择性和平移不变性的优点，而双变量模型可以充分利用图像小波系数的尺度间相关性。

步骤1　对噪声图像进行NSCT，可以设置每个尺度上的NSP分解层数和NSDFB分解方向子带数。

步骤2　估算噪声的NSCT系数方差 $ {\sigma }_{n}^{2} $ 。

步骤3　对于每个高频方向子带，使用局部自适应方法来估计每个系数的模型边缘方差，即 $ {\sigma }^{2}\left(x\right)= $ $ \mathrm{m}\mathrm{a}\mathrm{x} \left({\displaystyle\sum }_{v}^{2}\left(x\right)-{\sigma }^{2}\left(x\right),0\right) $ ,其中 $ {\displaystyle\sum }_{v}^{2}\left(x\right)=\dfrac{1}{MN} {\displaystyle\sum }_{m=1}^{M}{\displaystyle\sum }_{n=1}^{N} $ $ {y}_{m,n}^{2}\left(x\right) $ 。式中M，N为子带图像的长度和宽度，m和n代表系数所在位置。

步骤4　将 $ {\sigma }_{n}^{2} $ 和 $ {\sigma }^{2} $ 代入式(5)，以获得二元阈值 $ {\widehat{w}}_{1} $ 。

步骤5　建立模糊域特征平面矩阵，将低频分量带入到模糊域中，再利用改进的正余弦结合隶属度函数并对其特征增强，最后模糊逆变换为空间域矩阵。

步骤6　重建去噪的非下采样Contourlet系数矩阵，并执行逆NSCT变换以获得去噪够的图像。

3 实验对比分析

为了演示本文算法的去噪和增强性能，采用Matlab R2019b，计算机内存16 GB，处理器为intel(R)Core(TM)i7，主频为2.2 GHz作为实验环境。选择以灰度等级为512×512像素的声呐图像为测试样本，NSCT分解层数为3层。为了便于比较和分析，该算法与直方图均衡(HE)、离散小波(DWT)、双树复数小波(DTCWT)、非下采样轮廓波（NSCT）进行比较。表1比较了不同方法增强后图像计算信息熵以及峰值信噪比（PSNR）和平均绝对误差（MAE）三项性能指标。

水下恶劣环境下，声呐图像会产生伪Gibbs和目标边缘模糊的现象，图3为对前视声呐图像进行实验结果。

可以看出，NSCT变换能实现图像的最优稀疏分布逼近，其系数的双变量模型根据信息能量的不同，充分表示了尺度间以及尺度内的尺度相关性，实现变换系数自适应收缩，改进的隶属域函数不仅解决了传统模糊增强中参数因子计算复杂问题，同时采用正余弦结合“S”型函数，低灰度像素值有效保留。实验结果表明，该算法在去除噪声的同时尽可能的保留边缘细节，做到了整体和局部的增强，提高了恢复图像PSNR值，算法的主观视觉效果和客观评价标准都好于对比算法。

4 结　语

本文通过对含噪声呐图像进行NSCT分解，并在低频分量和高频分量上分别进行模糊集增强和变分贝叶斯阈值去噪增强。通过在模糊集上构造新的隶属度函数，并改进高频子带阈值自适应取值，使图像的纹理细节突出，图像的边缘更加清晰，对于舰船声呐图像的后续应用有重要意义。