CN103150706A - 一种改进的小波ica去噪方法 - Google Patents

Abstract

本发明公开了一种改进的小波ICA去噪方法，其特征在于，包括如下步骤：（1）对含噪声图像进行小波降维；（2）对小波降维后的图像进行小波ICA调整。本发明提出了一种改进的小波ICA滤波器，从图像源中将噪声与信号分离，该方法使用小波降维，通过ICA规范化降维后的信号，从而发现独立噪声特征，使用Morlet小波来解决非正交问题，并通过实验与PCA方法进行了对比，验证了该方法的效率。结果显示该方法降噪效果较PCA方法有大幅提高。

Description

一种改进的小波ICA去噪方法

技术领域

本发明涉及图像去噪技术领域，具体地讲，涉及一种改进的小波ICA去噪方法。 

背景技术

目前比较流行的低通滤波器方法主要通过平滑方法达到去噪的目的，而平滑方法会造成目标边缘模糊现象。在这些方法中，基于傅立叶变换的去噪方法主要工作在频域，基于小波变换的方法在频率和空间域上都有应用，但这两种方法的数据都不具备适应性。而在图像去噪过程中，如果滤波方法的数据具有适应性，会产生比较准确的结果，这主要是因为图像的去噪算法依赖于含噪图像的类型，因此算法的数据适应性在图像去噪中起到关键作用，而ICA方法恰恰具备数据适应性的本质特点。在武器系统中，光学设备获取的图像通常被白高斯噪声污染，这是不同的图像捕捉方法造成的，该噪声包含了可见光波段的所有频率成分。 

如何用ICA方法来消除比较复杂的噪声是目前研究界重点研究的问题，通常是对图像噪声的数量进行定量测量，先假设噪声图像的在特定范围内的变化已知，或通过不同的计算方法和优化方法得到。 

去噪算法通常对噪声的数量以及图像采样窗口的大小依赖性较强。如果采样窗口大，计算复杂度增加，算法的实施比较困难；反之，如果采样窗口较小，计算复杂度会降低，算法比较容易实现。同时，采样窗口的大小也会影响处理过程中引入的噪声数量。虽然大部份的算法也比较去噪前后的图像均方差BP以及信噪比SNR，但与视觉质量相比，这些定量计算的直观性要差一些。目前比较流行的 ICA方法主要包括以下几种： 

(1)适应性PCA方法 

主成份分析是线性数据适应变换技术，大家比较熟悉的是Hotelling变换。假定数据的维数为m，变换子空间的维数为n，则需要确定正交向量w_i，(i＝1,2,3…,n)。首先在上一步发现的方向上进行正交，然后找到数据最大差异的方向，再进行优化获得主成份，其数学表达式为： 

w_i＝arg maxw_iE{(w_i ^Tx)²}    (1) 

最大差异查找的限制条件为： 

w_i ^Tw_j＝0,j<I    (2) 

||w_i||=1    (3) 

通过计算相关矩阵C_x=E{x x^t}的特征向量来查找主成份方向矢量，该矩阵降低了相应特征值的序数，在这里特征值为特征向量方向上的差异。特征向量 

λ_i表示为 

λ_i＝E{(w_i ^Tx)²}    (4) 

通过得到的主成份向量w_i构造矩阵W，则向量w_i ^T为矩阵W的第i行，此时PCA变换矩阵为W。该变换可以通过现成的应用软件如计算无关向量、主成份降维等完成，从而作为ICA方法的预处理。可适应PCA技术主要用来产生2D的基集，其内的向量沿着边而不与其相交。选择2D局部可适应性基集，2D基函数通过主成分方法在采样窗口中产生，最大的特征向量在局部图像的边缘方向产生，通过局部适应主成份分析技术，经系数阈值过滤后，再进行逆变换重构降维后的图像。在图像的不同位置，局部基主成分需要重新计算。 

处理时首先估计图像中存在的噪声，噪声图像表示为： 

Y＝X＋n    （5) 

其中，Y是噪声图像观察数据，X为去噪后图像或者来源图像，n为图像存在的噪声，通常为零均值和方差为σ的白高斯噪声， 

图像被分割成重叠的块，每个块都包含训练域，去噪后的域和重叠域如图1 所示。Y的主要成份是局部适应基函数，Y是由局部噪声图像块组成的观测图像块，图像块的大小是M像素，即I=1,2,……M，矩阵Y中列数为M×M，矩阵Y的整体维数为M²×N，N为块的数量，主成份系数方差通过最大相似估计方法来计算。 

基降维比小波降维更高效，在可见光图像中，算法拥有固定的转换模式，若去噪幅度过大，即使用改善的PSNR算法，去噪后的图像包含的噪声也非常明显。为了使噪声在不同的区域均匀分布，通常使用重叠域进行处理。表1显示四种方法在噪声方差σ=50、25以及15时对四幅测试图像的效果。 

表1不同噪声方差的PSNR(dB)的比较 

(2)稀疏编码收缩方法 

稀疏编码收缩是一种数据适应性去噪技术，收缩参数通过数据本身进行计算，系数编码是一种线性变换，变换中每个成分量都起不到决定作用，在观测数据的最大相似估计基础上进行数据表示。假设变换成分与观测数据的非高斯方差相等，则从整体上通过不同的参数来计算最大相似估计。对该算法来说，收缩策略最为关键，系数分布在0处存在较重的峰值，通常称为超高斯性。收缩是一种独立成分去噪的阈值方法，通常采用两种成分阈值，即软阈值与硬阈值。在每次的计算中，只有少量的成分起关键作用，而其他绝对值较小的成分作为纯噪声，可以通过软阈值将不需要的分量置0，也可以将重要的图像数据模型参数经过阈值计算后替代原分量。 

令随机变化的观测数据为x＝(x₁，x₂，x₃,…x_n)^T，即数据或噪声图像窗口。s=(s₁,s₂s₃,…..s_n)^T为线性变换成分向量，若W为n×n矩阵，则s＝Wx。 

假设系数的分布不服从高斯分布，而噪声服从高斯分布，令非高斯随机变量 为s，高斯噪声表示为具有0均值与方差σ的f，此时噪声图像数据表示为 

Y＝s+f    (6) 

其中，s为去噪后的图像，若源用概率密度函数p来表示，f=-log(P)，则最大相似估计s’为： 

s’=argmin_u0.05σ^-2(y-u)²+f(u)    (7) 

其中f为一个完全可微的凸函数，问题的解为： 

g^-1(u)=u+σ²f’(u)    (8) 

该方法的关键是非线性估计，s的概率密度需要通过非线性估计进行建模，为使用非线性估计进行图像降噪。需要采用参数化方法，充分考虑可见光图像中遇到的最大密度，比较适合超高斯密度的公式为： 

p(s)=Cexp(-as²/2-b|s|)    (9) 

这里a、b为正数，C为比例常数或尺度常数。对于上述密度公式，其非线性为 

g(u)=1/(1+σ²a(sign(u)max(0，/u|-bσ²)))    (10) 

稀疏编码收缩方法的降噪过程为：首先使用噪声数据作为x的训练集，通过稀疏编码方法计算正交矩阵W，使s＝Wx中的成份s_i拥有尽可能的稀疏分布。对每个稀疏成分，用公式9估计其密度模型p_i(s_i)，对于x的每个观测数据x’，稀疏成分都表示为y＝Wx’，应用10给出的收缩非线性公式g_i()指示s’＝g_i(y_i)包含的成分。将s＝Wx进行逆运算得到去噪后的图像，其中W的行向量称为ICA滤波器。 

改进稀疏编码收缩方法通过补偿操作避免在滤波过程中丢失细节，采用滑动非方块窗口方法大大减少了计算复杂度，使用Hyvarinen提出的固定点快速ICA方法进行ICA变换。该方法可以概括为：首先将同类已经去噪的图像进行测试，通过与ICA变换矩阵正交来估计系数编码变换W，使用稀疏编码变换W而不是噪声图像进行变换，将每个成分都应用非线性估计g，同时都乘上补偿系数β，最后执行逆变换得到去噪后图像。 

改进的稀疏编码收缩算法对于随机图像的执行效率如表2所示，表2列出了不同尺度的窗口的PSNR均值。 

表2不同窗口尺寸与噪声级别的平均PSNR 

噪声层次	噪声图像	8×8	4×8	4×4	6×6
						0.1	32.27	37.10	37.17	37.00	37.03
0.2	26.25	33.16	33.22	33.06	33.17
						0.3	22.73	31.13	31.13	30.84	31.08
0.4	20.21	29.63	29.28	29.28	29.63
						0.5	18.26	28.62	28.57	28.10	28.51

(3)ICA混合模型： 

PCA属于二阶表示以及数据处理，但也可以用于高阶来观测图像的更多信息，因此ICA的图像表示要比PCA方法优秀。在ICA变换中，将数据表示为统计上独立的随机变量，图像数据看做包含白高斯噪声的非高斯随机变量。从而适合通过线性ICA进行高效求解。随机变量（噪声图像）通过X进行建模，表示为： 

X＝As+η    (11) 

s_i为统计上独立的随机向量，A为混合矩阵，η为独立向量s_i中的噪声，无噪数据通过重构的无噪图像表示：X’=As’。 

随机向量X与独立向量s_i的维度给出了ICA的三个不同条件：(1)在结束之前，随机向量X的维度高于独立向量s_i的维度；(2)结束情况下二者相同；(3)未确定的情况下随机向量X的维度低于独立向量s_i的维度。 

该方法的目的是从随机向量X中找到优化混合矩阵A以及独立源s。在观测向量与源向量维度相同的情况下，可以认为通过ICA方法完成了图像去噪，该过程通过让数据变白，找到正交非标准投影，去噪过程分为两步： 

a、研究噪声图像块，可以进行计算选取，然后存储在观察矩阵X的行中，再进行线性变换，得到g＝VX，这样g的成分具有单位方差并且线性无关，通过设置V＝C_x ^-1/2得到I＝E{gg^-1}，其中C_x＝E{XX¹}为数据的相关矩阵。 

b、计算正交矩阵B，该矩阵给出了代价函数的极大极小值，因此X可以通过X＝Bg计算。 

在基于ICA混合模型的图像去噪中，通过非高斯数据结构的参数化非线性函 给出每个类的数据密度。图像模型基于不同的数据密度函数构建，每个类的数据由高斯混合函数中的多变量高斯源使用高斯内核产生，用双曲线正切非线性快速ICA进行ICA估计，ICA变换通过正交化得到优于PCA的超高斯性。 

(4)独立子空间分析（Independent Subspace Analysis,ISA）: 

这是ICA的非线性方法，在子空间分析中，成分分成多个组，找出其独立性。该方法使得各组之间的成分相互独立，但组内成分相互依赖，组内的成分服从球对称稀疏密度： 

$P\left(s\right)={\mathrm{\&Pi;}}_k\exp (-f\left({\mathrm{\&Sigma;}}_{i\&times;G}{s}_i^2\right))---\left(12\right)$

其中，k=[1,2,…,G]为假设的子空间的数目，子空间相互独立，它们各自的概率密度乘积可以通过式(12)来表示。学习算法采用基于梯度上升的最大相似估计方法。图2显示组内成分具有相同的方向、空间频率以及位置，通常用该模型对局部噪声进行消除。 

(5)地形独立分析(Topographic Independent Analysis,ICA) 

在实际数据样本中，研究人员发现成分组的独立性达不到ICA的要求，相关成分间存在依赖，可以将这些成分放进1D或2D的网格中进行建模。假设网格中距离较近的成分具有能量相关性，距离较远的成分相互独立，则概率密度函数如下： 

$P\left(s\right)={\mathrm{\&Pi;}}_i\exp (-f\left({\mathrm{\&Sigma;}}_{j\&times;N_i}{s}_j^2\right))---\left(13\right)$

这里i与j两者的取值区间为[1,n]。 

该模型或多或少地与ISA模型相似，不同之处仅在于该方法具有重叠成分组，在该方法中，距离较近的特征具有相同的方向、空间频率以及位置。地形ICA的基在图3中给出，显示了靠近的特征具有相同的空间频率与位置。 

(6)多分辨率傅立叶变换(MFT-ICA) 

MFT-ICA组合MFT与ICA两种方法，使用二者优点如MFT的计算效率以及ICA的数据适应性。图像去噪中，方向信息是一个重要的成分，该方法使用方向基分析图像的统计结构。MFT基于多方向选择滤波器，而ICA比较适合将图像分解成一系 列方向基集，因此将两者联合的优势比较明显，MFT-ICA方法中，噪声图像v’、无噪图像v以及噪声η的关系为： 

v’=v+η    (14) 

ICA的分解信号存在如下关系： 

P(AnB)＝P(A).P(B)    (15) 

这里A与B为信号的独立成分，用ICA成份估计研究无关方向，将成分进行系数分布。首先，在拥有有限数量成分的MFT中，对每个子波段执行ICA；然后在成分基的傅立叶频谱中获取2D高斯滤波器。对大部分成分来说，无论在方向上还是频率上都是局部的，从而使逆变换得到无噪图像。 

发明内容

本发明要解决的技术问题是提供一种改进的小波ICA去噪方法，提高了含噪声图像的降噪效果。 

本发明采用如下技术方案实现发明目的： 

一种改进的小波ICA去噪方法，其特征在于，包括如下步骤： 

（1）对含噪声图像进行小波降维； 

（2）对小波降维后的图像进行小波ICA调整。 

作为对本技术方案的进一步限定，所述步骤（1）包括如下步骤： 

（1.1）将含噪图像标准化后进行输入； 

（1.2）确定分解的层次n，n为自然数； 

（1.3）使用Morlet小波对输入的含噪图像进行n层分解； 

（1.4）对Morlet小波分解后的系数进行量化。 

作为对本技术方案的进一步限定，所述步骤（2）包括如下步骤： 

（2.1）将小波系数进行线性变换，得到具有单位方差且线性无关的矩阵； 

（2.2）设置相关矩阵V； 

（2.3）计算正交矩阵B，得出代价函数的极大极小值； 

（2.4）一维小波重构，通过低频与高频的小波系数对图像进行重构。 

与现有技术相比，本发明的优点和积极效果是：本发明提出了一种改进的小波ICA滤波器，从图像源中将噪声与信号分离，该方法使用小波降维，通过ICA规范化降维后的信号，从而发现独立噪声特征，使用Morlet小波来解决非正交问题，并通过实验与PCA方法进行了对比，验证了该方法的效率。结果显示该方法降噪效果较PCA方法有大幅提高。 

附图说明

图1为PCA方法去噪后的域和重叠域的区域分割示意图。 

图2为稀疏矩阵编码方法基于独立子空间的分析图像。 

图3为地形ICA的图像基。 

图4为使用MFT-ICA对lena的去噪结果。 

图5为使用MFT-ICA对Jaguar去噪结果。 

图6中（a）为小波变换双通道滤波器示意图；（b）为小波变换3层滤波器示意图。 

图7不同频段内的独立分量的不同混合比展示图。 

图8为两种MICA算法与快速ICA算法的噪声标准偏差比较示意图。 

图9为对Lena图像去噪的效果比较图。 

图10为对包含不同比例噪声的Lena图像去噪的PSNR。 

图11为对包含不同比例噪声的Lena图像去噪的MAE。 

图12为对包含混合噪声Lady图像的恢复效果图。 

图13为本发明的处理流程图。 

具体实施方式

下面结合附图和优选实施例对本发明作更进一步的详细描述。 

一、改进的小波ICA滤波器 

(1)ICA模型 

ICA的本质是发现一个混合矩阵的逆，但由于没有提供混合矩阵的先验信息， 因而不能直接求解逆矩阵。根据中心极限理论，独立随机变量的总体分布服从高斯分布，从统计学观点来看，两个随机变量的总体分布比其中一个变量的总体分布更接近高斯分布，因此可以用非高斯性来显示收敛，著名的快速ICA算法通过负熵值来测量非高斯性，从而分离这些独立分量。负熵的定义为： 

J(x)＝H(x_gauss)-H(x)    (16) 

这里x_gauss为随机变量x中具有相同协方差矩阵的高斯随机变量集合，H(.)表示Shannon熵函数，定义为： 

H(x)＝-∫p(x)logp(x)dx    (17) 

根据信息理论的相关成果，一个高斯随机变量在所有的等方差随机变量中具有最高的熵值，因此熵总是非负的，仅仅在x服从高斯分布的情况下才取0值。通过概率密度分布估计计算熵值，近似表达式为： 

$J\left(x\right)\&ap;\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}{k}_i{\left[E\right\{G_i\left(x\right)\}-E\{G_i\left(v\right)\left\}\right]}^2---\left(18\right)$

这里E为理想期望值，k_i为正数常量，υ为0均值、单位方差的高斯随机变量，x为随机变量，也具有0均值与单位方差。对于一个给定的随机系列x，在使用公式(18)之前必须进行标准化。函数G_i为非二次函数，确保能够较好地近似。A.

和E.Oja提出了两个比较优秀的G_i函数得到较好的近似，分别为： 

$G_1\left(u\right)=\frac{1}{a_1}\log \cosh \left(a_{1}u\right),1\&le;a_1\&le;2---\left(19\right)$

G₂(u)＝-exp(-u²/2)    (20) 

比较流行的FastICA算法用这些G_i函数来实现熵最小化。 

(2)改进的小波降维方法 

信号与图像处理中的小波为研究人员提供了一种非常灵活的工具，可以用于不同的领域如语音与图像处理，在图像去噪领域，2D小波变换仅仅用来进行预处理，预处理减少了特征向量的维度并且可以去噪，然而，计算复杂度相对较高， 因此，本文提出了1D小波变换作为滤波器来降低特征向量的维度，可以大大降低计算复杂度。 

小波通过双通道滤波器进行构造，如图6所示，在进行1D信号降维的小波中，信号同时输入一个低通滤波器L与一个高通滤波器H，处理结果分别为低频成分A[n]与高频成分D[n]。 

小波滤波器通过回旋将信号分解为高频与低频成分 

$D\left[n\right]=\underset{k=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}s\left[k\right]\&CenterDot;H[n-k]\&DoubleLeftRightArrow;D=<s,H>---\left(21\right)$

$A\left[n\right]=\underset{k=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}s\left[k\right]\&CenterDot;L[n-k]\&DoubleLeftRightArrow;A=<s,L>---\left(22\right)$

为了构造多通道滤波器，我们可以采用流水通道滤波器，图6(b)显示了三级均衡流水结构，这对于多分辨率分析来说是一个重要概念。 

垂直投影技术的目的是为了减少图像的维度，变为1D能量轮廓信号，图像去噪的效率将会有极大提高，但是1D能量轮廓信号仍然非常巨大，为了克服该问题，我们借助小波变换来将1D能量信号作为特征向量来降低维度。 

(3)小波ICA滤波器 

ICA模型要求输入量必须是矩阵，不能输入向量，在小波ICA滤波器中，通常用小波降维构造输入矩阵，该方法的思想为： 

小波降维可以是不变的，但有的信号分量在降维后都保持相同的相位。为了简化算法，本文操作中不引入时间延迟。 

独立成分变量的混合比随着带通滤波器不同位置进行变化。图7假设s₁、s₂与s₃的分布为观测信号包含的三个独立分量。窗口T₁与T₂表示与两个不同的时不变线性变换函数相应的频率分布，很明显，s₁、s₂与s₃的混合比在窗口T₁与T₂内是不同的。 

s₁与s₂为信号x的两个独立分量，x_T1与x_T2为执行线性变换T₁与T₂后得到的结果，令s_1a与s_1b分别表示来自x_a与x_b中s₁的分布，s_2a与s_2b表示来自x_a与x_b中s₂的分布。在ICA的迭代过程中，信号中的分量需要处理成包含分量尽可能多的 数据，s_1a与s_1b也整合同样的分离成分。当线性变换采用小波变换时，该方法可以称为ICS滤波器，因为小波变换实质上可以看做一个带通滤波器组。 

但是，该方法的执行不如ICA的效果好，因为只有一组图像数据，并且图像混合成分中的独立分量互不相同。基于上述分析，为了显示小波ICA滤波器检测独立成分的可行形，必须证明两点： 

首先，必须证明s_1a或s_1b与s_2a或s_2b是独立的，令： 

$s_{1a}={\&Integral;}_{-\&infin;}^{\&infin;}{h}_a\left(\mathrm{\&tau;}\right)s_1(t-\mathrm{\&tau;})\mathrm{dt}$

$s_{1b}={\&Integral;}_{-\&infin;}^{\&infin;}{h}_b\left(\mathrm{\&tau;}\right)s_1(t-\mathrm{\&tau;})\mathrm{dt}$

$s_{2a}={\&Integral;}_{-\&infin;}^{\&infin;}{h}_a\left(\mathrm{\&tau;}\right)s_2(t-\mathrm{\&tau;})\mathrm{dt}$

$s_{2b}={\&Integral;}_{-\&infin;}^{\&infin;}{h}_b\left(\mathrm{\&tau;}\right)s_2(t-\mathrm{\&tau;})\mathrm{dt}$

h_a(t)与h_b(t)为域T₁与T₂相应的滤波函数，因此，s₁(t)与s₁(t)为两个独立分量，其中： 

$R_{s_{1a}{s}_{2a}}(t_1,t_2)={\&Integral;}_u{\&Integral;}_v{h}_a\left(u\right)h_a\left(v\right)E\left[s_1(t_1-u)s_2(t_2-v)\right]\mathrm{dudv}---\left(23\right)$

E[s₁(t₁)s₂(t₂)]＝0    (24) 

则： 

$R_{s_{1a}{s}_{2a}}(t_1,t_2)=0---\left(25\right)$

因为s₁(t)与s₂(t)为独立的，s_1a与s_2a在任何高阶序数下都不相关，同样，s_1b与s_2b也是如此。 

其次，要证明s_1a与s_1b或s_2a与s_2b是相关的，s_1a与s_1b间的关系函数为： 

$R_{s_{1a}{s}_{1b}}(t_1,t_2)={\&Integral;}_u{\&Integral;}_v{h}_a\left(u\right)h_b\left(v\right)E\left[s_1(t_1-u)s_1(t_2-v)\right]\mathrm{dudv}={\&Integral;}_u{\&Integral;}_v{h}_a\left(u\right)h_b\left(v\right)R_s(t_1-u,t_2-v)]\mathrm{dudv}---\left(26\right)$

R_s1(·)表示独立成分s₁的自相关函数，通常大于0，R_slas1b(t₁,t₂)不为0，当且仅当滤波器h_a(t)与h_b(t)是正交的，因此，可以证明R_s1as1b(t₁,t₂)>0，当且仅当线性滤波器为非正交。这就意味着s_1a与s_1b是时间相关的。 

对于非正交形，本文使用Morlet小波作为基小波，为了保留s_1a与s_1b，或s_2a与s_2b的相关性，线性滤波器h_a(t)与h_b(t)是相关的，来自Morlet小波线性滤波器可以是相关的，因为Morlet小波是非正交的，实际上，在h_a(t)与h_b(t)是线 性相关时，其他的基小波也可以使用，Morlet小波的定义为： 

为了分离无用的信号，小波变换定义为： 

$W(a,b)=\frac{1}{\sqrt{a}}\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}x\left(k\right)\mathrm{\&psi;}*\left(\frac{k-b}{a}\right)---\left(28\right)$

这里，*表示联合操作，a为尺度参数或扩张参数。 

本文使用改进的小波ICA消除附加噪声，被附加噪声污染的图像很难进行去噪，因为附加噪声不像标准白高斯噪声一样独立，比较合理的方法是在补偿情况下改变图像的数据统计，需要修改基于独立源的简单线性与非线性混合标准ICA模型。采用附加ICA（MICA）模型对图像进行去噪，同时使用二阶以及更高阶统计来解决该问题，被附加噪声污染的独立源混合模型为： 

Y_i=v_ix_i   fori=1,2,………..n    (29) 

x=As    (30) 

其中，Y＝[Y₁,Y₂,…Y_n]^T为向量表示的噪声图像块，s=[s₁,s₂,…s_n]^T为n个独立源向量，附加噪声成份为v＝[v₁,v₂,…v_n]^T，噪声消除数据由向量x给出，A为混合矩阵。混合矩阵的逆称为分离矩阵，也就是该问题的解。去噪过程分为两个步骤：首先是变白，该步骤对降维非常重要，同时提高了收敛速度。ICA的目标是发现一个非线性数据变换矩阵B，这样，输出就尽可能独立，矩阵B表示为： 

B＝UW    (31) 

这里W为变白的矩阵，U分离矩阵（归一矩阵），通过二阶、三阶、四阶图像统计结构来求解该问题。ICA问题将协方差、三阶与四阶累积量作为比较函数，通过线性与非线性最小化进行求解。这里基于三阶的方法称为Third order MICA(TMICA)，基于四阶的方法称为FMICA。FMICA比TMICA更加复杂，但由于TMICA的估计函数计算更准确，因而TMICA给出的结果比FMICA更准确。图8显示采用标准快速ICA方法对TMICA与FMICA的比较结果。 

二、实验及分析 

本文通过一系列512×512的测试图像进行试验来评估改进的小波ICA算法 的效率，所有的实验都采用3×3的过滤器窗口，使用峰值信噪比(peak signal-to-noise ratio,PSNR)与平均绝对误差(mean absolute error,MAE)对图像恢复能力与细节保留能力进行评估，PSNR定义为： 

$\mathrm{PSNR}=10\lg \left(\frac{{\mathrm{\&Sigma;}}_k{255}^2}{{\mathrm{\&Sigma;}}_k{(t\left(k\right)-y\left(k\right))}^2}\right)\mathrm{dB}---\left(32\right)$

这里255为图像的峰值灰度等级，t(k)表示期望输出值，y(k)为物理输出值，平均绝对误差(MAE)定义为： 

$\mathrm{MAE}=\frac{{\mathrm{\&Sigma;}}_k|t\left(k\right)-y\left(k\right)|}{H\&times;W}---\left(33\right)$

其中H与W为图像高度与宽度， 

噪声比为r的脉冲噪声模型描述为： 

其中g(k)与n(k)分别表示原始图像像素与噪声像素，在8位灰度图像中，脉冲噪声均与分布在[0,255]的范围内。另外，评估中也考虑了脉冲噪声与高斯噪声的混合噪声。 

首先将PCA与ICA进行，这里对PCA进行简单的回顾，在模式识别领域，广泛使用PCA对特征成分进行分离，假设特征向量x有两种成分x₁与x₂组成，通过有N个采样点进行采样，表示为： 

X₁₁,X₁₂,…,X_1N

X₂₁,X₂₂,…,X_2N

使用PCA方法查找新的坐标轴D₁与D₂，使观测数据在D₁的投影具有最大标准偏差，然后，D₁的方向将会携带原始数据集的最多信息，可以用来逼近两个原始成分，D₁方向称为主要方向，该方向的投影称为主成分，很明显，PCA是一种线性变换，旋转一个坐标轴到主成分方向。 

y₁＝a₁₁x₁+a₁₂x₂    (35) 

y₂＝a₂₁x₁+a₂₂x₂

为了简化，将其写为： 

Y＝AX    (36) 

则 

YYT＝AXXTAT    (37) 

代入期望值公式： 

E(YYT)＝A·E(XXT)·AT    (38) 

或 

C_y=Ac_xAT    (39) 

C_x与C_y为x与y的协方差矩阵，一个合适的矩阵A可以使C_y变成对角线矩阵。换句话说，所有的主元素都是无关的因为C_y的对角元素正好是每个主成分Y₁,Y₂,…,Y_m的方差，这就是为什么通常使用PCA进行特征分离，不同于PCA，ICA使这些成分独立，而不是仅仅在二阶累积量上分离，等同于在更高序数累积量上对他们进行分离，因此，ICA可以认为是PCA的一个泛化。 

现在使用PCA、FastICA以及小波ICA分别处理图像，都使用小波降维作为处理方法，对ICA方法，全部的处理过程构成小波ICA滤波器。通过对比试验验证了提出小波ICA滤波器的效率，使用被20%的脉冲噪声污染的图像作为实验图像。表3显示了三种方法对包含20%脉冲噪声图像去噪时PSNR的结果比较。可以看出，本文方法的效果远好于其他方法，可以得到较高的PSNR值。本文方法恢复的图像效果较好，包含较低的噪声，保留了更多的细节。并给出了最小的MAE值。为了证实本文算法的鲁棒性，图10与图11显示了PSNR与MAE的比较，相应地，对于包含2%-50%脉冲噪声的Lena图像。图10与图11清楚地显示了本文方法远胜过其他滤波器。 

表3处理20%脉冲噪声图像的PSNR比较 

虽然本文方法最初设计为消除脉冲噪声，但对于高斯与脉冲混合噪声的去噪 效果也非常好，本文也比较了本方法与其他方法对混合噪声的消除效果，表4显示本文方法的效果非常好，比其他方法提供更理想的结果。 

表4去除混合噪声的PSNR对比 

	Baboon	Bridge	Boats	Lake	Lena	Goldhill
							PCA	21.62	24.41	26.81	24.52	25.37	26.12
FastICA	20.15	23.78	23.56	23.92	24.12	24.05
							本文方法	19.85	22.24	22.12	22.49	23.53	23.79

图12显示了对于包含混合噪声（高斯噪声(σ=20)脉冲噪声(20%)）Lady图像的去噪效果，如图12所示，本文方法可以得到最好的可视效果，噪声斑点最少，MAE值也最低。 

当然，上述说明并非对本发明的限制，本发明也不仅限于上述举例，本技术领域的普通技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也属于本发明的保护范围。 

Claims (3)

1.一种改进的小波ICA去噪方法，其特征在于，包括如下步骤：

（1）对含噪声图像进行小波降维；

（2）对小波降维后的图像进行小波ICA调整。

2.根据权利要求1所述改进的小波ICA去噪方法，其特征在于，所述步骤（1）包括如下步骤：

（1.1）将含噪图像标准化后进行输入；

（1.2）确定分解的层次n，n取自然数；

（1.3）使用Morlet小波对输入的含噪图像进行n层分解；

（1.4）对Morlet小波分解后的系数进行量化。

3.根据权利要求1所述改进的小波ICA去噪方法，其特征在于，所述步骤（2）包括如下步骤：

（2.1）将小波系数进行线性变换，得到具有单位方差且线性无关的矩阵；

（2.2）设置相关矩阵V；

（2.3）计算正交矩阵B，得出代价函数的极大极小值；

（2.4）一维小波重构，通过低频与高频的小波系数对图像进行重构。

Images