CN102426701A - 基于双树复小波变换和pca的水下声纳图像的去噪方法 - Google Patents

Abstract

本发明的目的在于提供基于双树复小波变换和PCA的水下声纳图像的去噪方法，分为以下步骤：对一幅水下声纳图像应用双树复小波变换，将图像由空间域变换到复小波域，保持图像经三层双树复小波变换后获得的低频近似分量不变，对图像的高频分量进行处理，采用PCA方法估计高频子带中噪声的能量，从而确定阈值并采用硬阈值函数对复小波系数进行处理，对处理后的复小波系数进行双树复小波反变换，获得最终去噪后的图像。本发明克服了传统二维小波缺乏平移不变性和方向选择性的缺点，更好地捕捉图像的方向性信息，能够在去除噪声的同时，更好地保护图像的边缘、纹理等细节信息。

Description

基于双树复小波变换和PCA的水下声纳图像的去噪方法

技术领域

本发明涉及的是一种图像处理技术领域的去噪方法。

背景技术

目前图像去噪的方法的种类很多，常用的方法例如均值滤波、中值滤波、小波去噪等方法。均值滤波是线性的空间滤波器，它用掩膜确定的邻域像素的平均灰度值代替图像中每个像素点的值，这样的处理减少了图像灰度的“尖锐”变化，由于典型的随机噪声由灰度级的尖锐变换组成，因此经过均值滤波的图像噪声将会被减少。中值滤波是一种非线性滤波方式，它是将邻域内像素灰度的中值代替该像素的值，它的应用非常普遍，对很多种的随机噪声都有良好的去噪能力。小波去噪也是应用十分广泛的去噪方法，它将图像转换至小波域，采用阈值函数对小波系数进行处理，达到去噪的目的。

但是在图像的边缘，纹理等方向性信息保留方面，这些方法的处理效果并不是十分理想。

发明内容

本发明的目的在于提供能够在有效去除图像噪声的同时较好地保留图像的边缘、纹理等细节信息的基于双树复小波变换和PCA的水下声纳图像的去噪方法。

本发明的目的是这样实现的：

本发明基于双树复小波变换和PCA的水下声纳图像的去噪方法，其特征是：

(1)对水下声纳图像进行双树复小波变换：树A和树B分别代表复数小波的实部和虚部，它们采用不同的滤波器组，h₀(n)是树A的低通滤波器，h₁(n)是树A的高通滤波器，n是滤波器的长度，与之对应的尺度函数φ_h(t)和小波函数ψ_h(t)分别为：

g₀(n)是树B的低通滤波器，g₁(n)是树B的高通滤波器，n是滤波器的长度，与之对应的尺度函数φ_g(t)和小波函数ψ_g(t)分别为：

图像经双树复小波变换分解成3层，获得三层的子带系数，y₀代表经过3层分解后图像的低频近似分量的系数矩阵，

分别表示第一、二、三层的高频分量系数，他们各自包含6个系数矩阵，分别代表双树复小波变换分解的6个方向的高频细节信息；

每一层6个方向的变换系数实数部分为：

式中

其中ψ_1，i(x，y)和ψ_2，i(x，y)为：

x和y代表水平方向和垂直方向，φ_h和ψ_h分别是树A的尺度函数和小波函数，φ_g和ψ_g分别是树B的尺度函数和小波函数；

变换系数的虚数部分为：

式中

其中ψ_3，i(x，y)和ψ_4，i(x，y)为：

每一层分解后的6个方向的高频系数y_k为：

其中

(2)保持图像经三层双树复小波变换后获得的低频近似分量的系数矩阵y₀不变，对图像的高频分量y₁、y₂、y₃进行去噪处理：

包含6个系数矩阵，分别是每层分解后获得的6个方向的高频细节信息，对第一层的高频分量进行处理：

表示第一层高频分量y₁中任意一个方向的系数矩阵，将分解成大小为16λ16的子块，以子块为单位进行处理，用X表示任意子块，按以下步骤对任意子块进行处理：

计算X的协方差矩阵C：

矩阵X为16个16维的行向量X_j组成的矩阵，它的协方差矩阵为：

其中

求取协方差矩阵C的特征向量和特征值：

协方差矩阵特征值分解具有如下形式：

其中

是协方差矩阵的特征向量矩阵，μ₁，μ₂，......，μ₁₆是协方差矩阵C的特征向量，是协方差矩阵的特征值构成的对角阵，λ₁，λ₂，......，λ₁₆是协方差矩阵C的特征值，这些特征值按从大到小的顺序排列，也就是满足λ₁μλ₂μ......μλ₁₆，此时与之对应的特征向量μ₁，μ₂，......，μ₁₆便构成了特征空间的一组基；

根据特征向量对矩阵X进行重构，确定去噪阈值T：

通过计算特征值的方差贡献率η，

其中trace(C)是协方差矩阵C的迹，λ₁，λ₂，......，λ_m是协方差矩阵C的特征值，m的取值为1，2，...，16，选取前8个特征值，通过与之对应的特征向量对X进行重构，重构矩阵Y为：

其中

它是选取的前8个特征值所对应的特征向量组成的矩阵，μ₁，μ₂，......，μ₈是选取的8个特征向量，将Y矩阵中系数幅值的均值确定为去噪的阈值T；

采用硬阈值函数对子块X中的系数进行去噪处理：

硬阈值函数表达式如下：

其中X(i，j)为子块内的变换系数，|X(i，j)|是系数的模值，将子块中的每个系数都采用硬阈值函数进行处理，如果系数的幅值大于或者等于阈值T，那么系数保持不变；如果系数的幅值小于阈值T，则将该系数置为零，第二、三层的高频分量y₂和y₃也采用同样的方法进行处理；

(3)对处理后的系数进行双树复小波反变换，获得最终的去噪图像。

本发明的优势在于：本发明克服了传统二维小波缺乏平移不变性和方向选择性的缺点，更好地捕捉图像的方向性信息，能够在去除噪声的同时，更好地保护图像的边缘、纹理等细节信息。

附图说明

图1为本发明的流程图；

图2是一维双树复小波分解滤波器结构图；

图3是本发明与小波硬阈值去噪、小波软阈值去噪的效果对比图，图3(a)是含有噪声的水下声纳图像，噪声方差为

图3(b)是采用小波硬阈值去噪处理后的图像，图3(c)是采用小波软阈值去噪处理后的图像，图3(d)是采用本发明的去噪方法处理后的图像；

图4是三种去噪方法的峰值信噪比(PSNR)值对比图；

图5是三种去噪方法的均方值误差(MSE)值对比图。

具体实施方式

下面结合附图举例对本发明做更详细地描述：

结合图1～5，具体步骤是：

(1)对一幅水下声纳图像应用双树复小波变换，将图像由空间域变换到复小波域。

(2)保持图像经三层双树复小波变换后获得的低频近似分量不变。

(3)对图像的高频分量进行处理。采用PCA方法估计高频子带中噪声的能量，从而确定阈值并采用硬阈值函数对复小波系数进行处理。

(4)对处理后的复小波系数进行双树复小波反变换，获得最终去噪后的图像。

本发明的核心内容在于采用双树复小波变换将图像由空间域变换至复小波域，在复小波域进行去噪处理，并且采用PCA的方法确定阈值，通过硬阈值函数对复小波系数进行处理。

本发明包括的对图像进行双树复小波变换的主要内容为：对图像进行二维双树复小波变换，从而将图像的高低频信息分离。双树复小波变换变换层数为三层，变换后得到一个低频近似分量的系数矩阵，在每一层获得6个高频系数矩阵，分别表示6个方向的高频细节信息。

本发明包括的PCA去噪方法的主要内容为：保留图像的低频近似分量，采用PCA的方法处理图像的高频分量。将各个高频子带分解成大小为16λ16的若干个子块，以每一子块作为单位进行处理。采用PCA的方法确定每一子块系数的去噪阈值，采用硬阈值函数对每一子块的系数进行处理。对处理后的系数进行双树复小波反变换，获得最终的去噪图像。

本发明的优点在于采用双树复小波变换，克服了传统二维小波缺乏平移不变性和方向选择性的缺点，更好地捕捉图像的方向性信息，能够在去除噪声的同时，更好地保护图像的边缘、纹理等细节信息。

采用水下声纳图像，图像尺寸大小为256λ256。实验仿真中分别采用小波硬阈值去噪，小波软阈值去噪和本发明的去噪方法对含有相同噪声的图像进行去噪处理。

图1所示的是本发明的流程意图，参照图1本发明的具体实施步骤如下：

(1)选取一幅水下声纳图像，对该图像进行双树复小波变换对选取的声纳图像进行双树复小波变换，双树复小波变换是基于实数小波变换实现的复数小波变换，是一种可分离的复数小波变换。它通过两组并行的实数滤波器组实现。图2所示的是一维双树复小波变换的分解滤波器组的结构图。图中树A和树B分别代表复数小波的实部和虚部，它们采用不同的滤波器组。h₀(n)是树A的低通滤波器，h₁(n)是树A的高通滤波器，n是滤波器的长度，与之对应的尺度函数φ_h(t)和小波函数ψ_h(t)分别为：

g₀(n)是树B的低通滤波器，g₁(n)是树B的高通滤波器，n是滤波器的长度，与之对应的尺度函数φ_g(t)和小波函数ψ_g(t)分别为：

对图像进行双树复小波变换，是分别对图像进行水平和垂直方向的一维双树复小波变换。图像经双树复小波变换分解成3层，获得三层的子带系数，y₀代表经过3层分解后图像的低频近似分量的系数矩阵，分别表示第一、二、三层的高频分量系数，他们各自包含6个系数矩阵，分别代表双树复小波变换分解的6个方向的高频细节信息。

每一层6个方向的变换系数实数部分为：

式中

其中ψ_1，i(x，y)和ψ_2，i(x，y)为：

上式中的x和y代表的是水平方向和垂直方向。φ_h和ψ_h分别是树A的尺度函数和小波函数，φ_g和ψ_g分别是树B的尺度函数和小波函数。

变换系数的虚数部分为：

式中

其中ψ_3，i(x，y)和ψ_4，i(x，y)为：

上式中的x和y代表的是水平方向和垂直方向。φ_h和ψ_h分别是树A的尺度函数和小波函数，φ_g和ψ_g分别是树B的尺度函数和小波函数。因此每一层分解后的6个方向的高频系数y_k为：

其中

(2)保持图像经三层双树复小波变换后获得的低频近似分量的系数矩阵y₀不变。

(3)对图像的高频分量y₁、y₂、y₃进行去噪处理。

包含6个系数矩阵，分别是每层分解后获得的6个方向的高频细节信息。首先对第一层的高频分量进行处理。表示第一层高频分量y₁中任意一个方向的系数矩阵，，我们将

分解成大小为16λ16的若干个子块，以子块为单位进行处理。我们用X表示任意子块，按以下步骤对任意子块进行处理。

1)计算X的协方差矩阵C；

矩阵X可以看作是由16个16维的行向量X_j组成的矩阵，它的协方差矩阵为：

其中

本发明中

2)求取协方差矩阵C的特征向量和特征值；

协方差矩阵特征值分解具有如下形式：

其中

是协方差矩阵的特征向量矩阵，μ₁，μ₂，......，μ₁₆是协方差矩阵C的特征向量。

是协方差矩阵的特征值构成的对角阵，λ₁，λ₂，......，λ₁₆是协方差矩阵C的特征值。这些特征值按从大到小的顺序排列，也就是满足λ₁μλ₂μ......μλ₁₆，此时与之对应的特征向量μ₁，μ₂，......，μ₁₆便构成了特征空间的一组基。

3)根据特征向量对矩阵X进行重构，确定去噪阈值T。

通过计算特征值的方差贡献率η，

其中trace(C)是协方差矩阵C的迹，λ₁，λ₂，......，λ_m是协方差矩阵C的特征值，m的取值为1，2，...，16，随着n的增加，方差贡献率η也随之增大，当方差贡献率η足够大时(比如达到80％或90％)，前m个特征向量可以表征子块的主要特征。根据计算结果，前8个特征值的贡献率达到了90％，表明前8个特征向量可以表征子块的主要特征，因此选取前8个特征值，通过与之对应的特征向量对X进行重构，重构矩阵Y为：

其中

它是选取的前8个特征值所对应的特征向量组成的矩阵，μ₁，μ₂，......，μ₈是选取的8个特征向量。通过PCA方法重构的矩阵Y反应了由于噪声影响系数幅值的特征，将Y矩阵中系数幅值的均值确定为去噪的阈值T。

4)采用硬阈值函数对子块X中的系数进行去噪处理。

硬阈值函数表达式如下：

其中T是3)中确定的去噪阈值，X(i，j)为子块内的变换系数，|X(i，j)|是系数的模值。将子块中的每个系数都采用硬阈值函数进行处理，如果系数的幅值大于或者等于阈值T，那么系数保持不变；如果系数的幅值小于阈值T，则将该系数置为零。第二、三层的高频分量y₂和y₃也采用以上的步骤进行处理。(4)对处理后的系数进行双树复小波反变换，获得最终的去噪图像。并采用小波硬阈值、软阈值去噪的方法对同一幅图像进行处理，分别计算三种去噪方法的峰值信噪比(PSNR)和均方根误差(MSE)。

图3(a)是含有噪声的水下声纳图像，噪声方差为

图3(b)是采用小波硬阈值去噪处理后的图像，图3(c)是采用小波软阈值去噪处理后的图像，图3(d)是采用本发明的去噪方法处理后的图像。从图3中可以看出本发明的去噪效果要优于小波硬阈值和软阈值去噪，在有效地去除噪声的同时能够更好地保留图像的细节和边缘信息。

为了更加具体的比较三种去噪方法的去噪效果，采用PSNR和MSE作为评价指标，分别计算三种方法的PSNR和MSE，图4是三种去噪方法的PSNR对比，图5是三种去噪方法的MSE对比。PSNR是衡量处理后的图像与原始图像的差异程度的评价指标，PSNR值越高表明处理后的图像与原始图像的相似程度越大。从图4可以看出，采用本发明去噪后图像的PSNR的提升幅度要优于小波硬阈值和小波软阈值去噪，当信噪比较高的情况下，本发明的去噪方法的PSNR仍有较大的提高。MRSE是用来衡量图像之间的差异程度的另一个常用指标，MRSE的值越小，表明两幅图像越接近。从图5中可以看出，采用本发明的去噪方法处理后的图像的MSE最小，与原始图像更为接近，去噪效果更好。

Claims (1)

1.基于双树复小波变换和PCA的水下声纳图像的去噪方法，其特征是：

(1)对水下声纳图像进行双树复小波变换：树A和树B分别代表复数小波的实部和虚部，它们采用不同的滤波器组，h₀(n)是树A的低通滤波器，h₁(n)是树A的高通滤波器，n是滤波器的长度，与之对应的尺度函数φ_h(t)和小波函数ψ_h(t)分别为：

${\mathrm{\φ}}_h\left(t\right)=\sqrt{2}\underset{n}{\mathrm{\Σ}}{h}_0\left(n\right){\mathrm{\φ}}_h(2n-t),$

${\mathrm{\ψ}}_h\left(t\right)=\sqrt{2}\underset{n}{\mathrm{\Σ}}{h}_1\left(n\right){\mathrm{\φ}}_h(2n-t),$

g₀(n)是树B的低通滤波器，g₁(n)是树B的高通滤波器，n是滤波器的长度，与之对应的尺度函数φ_g(t)和小波函数ψ_g(t)分别为：

${\mathrm{\φ}}_g\left(t\right)=\sqrt{2}\underset{n}{\mathrm{\Σ}}{g}_0\left(n\right){\mathrm{\φ}}_g(2n-t),$

${\mathrm{\ψ}}_g\left(t\right)=\sqrt{2}\underset{n}{\mathrm{\Σ}}{g}_1\left(n\right){\mathrm{\φ}}_g(2n-t),$

图像经双树复小波变换分解成3层，获得三层的子带系数，y₀代表经过3层分解后图像的低频近似分量的系数矩阵，y_k(k＝1，2，3)分别表示第一、二、三层的高频分量系数，他们各自包含6个系数矩阵，分别代表双树复小波变换分解的6个方向的高频细节信息；

每一层6个方向的变换系数实数部分为：

ψ_i(x，y)＝ψ_1，i(x，y)-ψ_2，i(x，y)，

ψ_i+3(x，y)＝ψ_1，i(x，y)+ψ_2，i(x，y)，

式中i＝1，2，3，其中ψ_1，i(x，y)和ψ_2，i(x，y)为：

ψ_1，1(x，y)＝φ_h(x)ψ_h(y)，ψ_2，1(x，y)＝φ_g(x)ψ_g(y)

ψ_1，2(x，y)＝ψ_h(x)φ_h(y)，ψ_2，2(x，y)＝ψ_g(x)φ_g(y)

ψ_1，3(x，y)＝ψ_h(x)ψ_h(y)，ψ_2，3(x，y)＝ψ_g(x)ψ_g(y)

x和y代表水平方向和垂直方向，φ_h和ψ_h分别是树A的尺度函数和小波函数，φ_g和ψ_g分别是树B的尺度函数和小波函数；

变换系数的虚数部分为：

ψ_i(x，y)＝ψ_3，i(x，y)-ψ_4，i(x，y)，

ψ_i+3(x，y)＝ψ_3，i(x，y)+ψ_4，i(x，y)，

式中i＝1，2，3，其中ψ_3，i(x，y)和ψ_4，i(x，y)为：

ψ_3，1(x，y)＝φ_g(x)ψ_h(y)，ψ_4，1(x，y)＝φ_h(x)ψ_g(y)

ψ_3，2(x，y)＝ψ_g(x)φ_h(y)，ψ_4，2(x，y)＝ψ_h(x)φ_g(y)；

ψ_3，3(x，y)＝ψ_g(x)ψ_h(y)，ψ_4，3(x，y)＝ψ_h(x)ψ_g(y)

每一层分解后的6个方向的高频系数y_k为：

y_k(i)＝ψ_1，i(x，y)-ψ_2，i(x，y)+j[ψ_3，i(x，y)-ψ_4，i(x，y)]，

y_k(i+3)＝ψ_1，i(x，y)+ψ_2，i(x，y)+j[ψ_3，i(x，y)+ψ_4，i(x，y)]，

其中k＝1，2，3；

(2)保持图像经三层双树复小波变换后获得的低频近似分量的系数矩阵y₀不变，对图像的高频分量y₁、y₂、y₃进行去噪处理：

y_k(k＝1，2，3)包含6个系数矩阵，分别是每层分解后获得的6个方向的高频细节信息，对第一层的高频分量进行处理：y₁(l)(l＝1，2，...，6)表示第一层高频分量y₁中任意一个方向的系数矩阵，将y₁(l)(l＝1，2，...，6)分解成大小为16×16的子块，以子块为单位进行处理，用X表示任意子块，按以下步骤对任意子块进行处理：

计算X的协方差矩阵C：

矩阵X为16个16维的行向量X_j组成的矩阵，它的协方差矩阵为：

$C=\frac{1}{K}\underset{j=1}{\overset{K}{\mathrm{\Σ}}}(X_j-\stackrel{\‾}{X}){(X_j-\stackrel{\‾}{X})}^T,$ 其中 $\stackrel{\‾}{X}=\frac{1}{K}\underset{j=1}{\overset{K}{\mathrm{\Σ}}}{X}_j,K=16;$

求取协方差矩阵C的特征向量和特征值：

协方差矩阵特征值分解具有如下形式：

C＝U^TVU，

其中U＝[μ₁，μ₂，......，μ₁₆]是协方差矩阵的特征向量矩阵，μ₁，μ₂，......，μ₁₆是协方差矩阵C的特征向量，V＝diag[λ₁，λ₂，......，λ₁₆]是协方差矩阵的特征值构成的对角阵，λ₁，λ₂，......，λ₁₆是协方差矩阵C的特征值，这些特征值按从大到小的顺序排列，也就是满足λ₁≥λ₂≥......≥λ₁₆，此时与之对应的特征向量μ₁，μ₂，......，μ₁₆便构成了特征空间的一组基；

根据特征向量对矩阵X进行重构，确定去噪阈值T：

通过计算特征值的方差贡献率η，

$\mathrm{\η}=\frac{{\mathrm{\λ}}_1+{\mathrm{\λ}}_2+...+{\mathrm{\λ}}_m}{\mathrm{trace}\left(C\right)}$

其中trace(C)是协方差矩阵C的迹，λ₁，λ₂，......，λ_m是协方差矩阵C的特征值，m的取值为1，2，...，16，选取前8个特征值，通过与之对应的特征向量对X进行重构，重构矩阵Y为：

Y＝U_d ^TXU_d，

其中U_d＝[μ₁，μ₂，......，μ₈]，它是选取的前8个特征值所对应的特征向量组成的矩阵，μ₁，μ₂，......，μ₈是选取的8个特征向量，将Y矩阵中系数幅值的均值确定为去噪的阈值T；

采用硬阈值函数对子块X中的系数进行去噪处理：

硬阈值函数表达式如下：

$X(i,j)=\left\{\begin{array}{cc}X(i,j)& \left|X(i,j)\right|\&times;T\\ 0& \left|X(i,j)\right|<T\end{array}\right.,$

其中X(i，j)为子块内的变换系数，|X(i，j)|是系数的模值，将子块中的每个系数都采用硬阈值函数进行处理，如果系数的幅值大于或者等于阈值T，那么系数保持不变；如果系数的幅值小于阈值T，则将该系数置为零，第二、三层的高频分量y₂和y₃也采用同样的方法进行处理；

(3)对处理后的系数进行双树复小波反变换，获得最终的去噪图像。

Images