CN102289800A - 基于Treelet的Contourlet域图像去噪方法 - Google Patents

Abstract

本发明公开了一种基于Treelet的Contourlet域图像去噪方法，主要解决现有去噪方法去噪效果不佳的问题。其实现步骤是：(1)由含噪图像NI求出差变量矩阵；(2)对NI进行非局部均值预滤波，得到滤波后图像LI；(3)对LI进行平移，得到平移后图像PI，对PI进行Contourlet分解；(4)对分解出的各个高频子带分别进行去噪，并对去噪后的子带进行Contourlet逆变换，得到去噪图像FI；(5)将FI逆平移，得到去噪图像；(6)将步骤(3)-(5)重复八次去噪，对去噪的八幅图像取平均后输出。本发明能有效去除含有高斯白噪声的自然图像中的噪声，可用于变化检测，目标识别时对图像的预处理。

Description

基于Treelet的Contourlet域图像去噪方法

技术领域

本发明属于图像处理技术领域，涉及被高斯白噪声腐蚀的自然图像的去噪，可用于开展土地利用与覆盖的变化检测、环境变化评估、城市规划、医学影像等领域的数字图像预处理。

背景技术

图像去噪的主要目的是解决实际图像由于噪声干扰导致的图像质量下降问题。通过去噪可以提高图像质量，增大信噪比，更好地体现图像所携带的信息，因此图像去噪技术在很多领域中都占据着重要位置。

根据图像的特点和噪声的统计特征，多年来已经提出了很多图像去噪方法，现有的去噪方法主要分为空间域滤波和频率域滤波。图像的空间域滤波方法直接对图像的灰度做运算，图像的频率域滤波方法则是在变换域内对图像的变换系数进行运算，然后反变换到图像空间域。小波变换是频率域滤波方法中最具代表性的，但是小波变换使用固定的正方形基函数去逼近原始图像，只能有效地表示点奇异信息，不能有效表示图像的二维奇异信息，如线、轮廓等重要信息，从而制约了小波降噪方法的性能。而Contourlet变换不仅继承了小波变换的多分辨率时频分析特征，具备了小波变换所没有的多方向性和各向异性，并采用可调的长方形基函数去逼近原始图像，解决了小波变换不能有效表示二维奇异性的缺点，能准确地将图像中的边缘捕获到不同尺度、不同方向的子带中，实现了对图像信号的稀疏分离，因此选择合适的阈值进行去噪就能获得比小波变换更好的效果。

现有的Contourlet域图像去噪方法主要有以下几种：

(1)Wiener滤波：Zhou等在“Contourlet-based Image Denoising Algorithm usingAdaptive Windows.ICIEA，2009：3654-3657.”中提出了一种基于自适应窗口的Contourlet域Wiener滤波去噪方法。该方法首先在高频子带内选择各向异性窗口，然后对高频系数进行该窗口的Wiener滤波。该方法的缺点是噪声大量残留且没有考虑到Contourlet缺乏平移不变性，导致去噪结果有严重的伪Gibbs现象。

(2)系数建模：D D.-Y.Po等在“Directional multiscale modeling of images usingthe Contourlet transform.IEEE Transactions on Image Processing，2006，15(6)：1610-1620.”中提出了一种对Contourlet系数进行HMT建模的去噪方法。该方法考虑了高频系数及其父系数，邻域系数的相关性，对系数进行HMT建模，根据联合概率分布进行去噪。该方法的缺点是复杂度高且没有考虑到Contourlet缺乏平移不变性，导致去噪结果有严重的伪Gibbs现象。

(3)Cycle Spinning：梁栋等在“一种基于Contourlet递归Cycle Spinning的图像去噪方法.电子学报，2005，33(11)：2044-2046.”中提出了一种基于Cycle Spinning的Contourlet域去噪方法。该方法采用递归循环平移(Cycle Spinning)来克服采用Contourlet去噪时产生的伪Gibbs现象，采用了窗口硬阈值对自然图像进行去噪。该方法的缺点是破坏了图像的细节。

发明内容

本发明的目的在于克服上述去噪方法的缺点，提出一种基于Treelet的Contourlet域图像去噪方法，以降低计算复杂度，去除伪Gibbs现象，保持图像细节，减少噪声残留，提高去噪后图像的清楚度。

实现本发明的技术方案是通过引入一种新的自适应多尺度分析和表示方法-Treelet，对高维Contourlet系数进行逐层降维分解，反映出系数潜在的结构以及系数之间的相关性，对系数进行了准确的分析，从而估计出Contourlet高频子带的自适应阈值，其具体步骤包括如下：

(1)由大小为A×B的含噪图像矩阵NI求出自协方差矩阵S_s和自相关系数矩阵M_s；

(2)设Treelet分解的最高层数为L，每次分解所在层为lev，lev取值为1，.....，L，在lev＝1层，即Treelet的最底层，令Treelet的初始输入S⁽⁰⁾＝S_s，M⁽⁰⁾＝M_s，X⁽⁰⁾＝NI，L＝B-1，B是NI的列数，A是NI的行数；

(3)将上述输入S⁽⁰⁾，M⁽⁰⁾，X⁽⁰⁾和L代入Treelet，对X⁽⁰⁾进行Treelet变换分解到最高层，得到差变量矩阵{dif_i}，dif_i为Treelet在第i层分解得到的差变量，i取值为1，...，L；

(4)将图像矩阵NI投影在步骤(3)所得的差变量矩阵{dif_i}上，估计出非局部均值滤波的平滑因子h，然后对含噪图像NI进行搜索窗为7×7，相似窗为3×3，平滑因子为h的非局部均值滤波，得到滤波后图像LI；

(5)设步骤(6)-(9)中所用的循环平移的初始值a＝0，b＝0，cycle＝1，a为平移的行数，b为平移的列数，cycle为平移的次数；

(6)对滤波后图像LI进行循环平移，下移a行，右移b列，得到平移后的图像PI，对PI进行Contourlet变换，分解为3层，第一层分解16个高频子带，第二层分解8个高频子带，第三层分解8个高频子带，分别对各个高频子带进行步骤(7)-(8)的处理；

(7)由大小为P×N的高频子带Y求出自协方差矩阵

和自相关系数矩阵

得到Treelet变换的初始输入：

X⁽⁰⁾＝Y，L＝N-1，其中N为Y的列数，P是Y的行数；

(8)将步骤(7)得到的Treelet的初始输入S⁽⁰⁾，M⁽⁰⁾，X⁽⁰⁾及L代入步骤(3)中，得到差变量矩阵{dif_i}，根据{dif_i}估计阈值Th并计算权重g，然后对高频子带Y依次进行软阈值去噪和独立点噪声去除；

(9)将所有经过步骤(7)-(8)处理过的高频子带进行Contourlet逆变换，得到平移后的去噪图像FI，将FI左移b列，上移a行，得到去噪图像DI(cycle)，令cycle＝cycle+1，若cycle＜＝8，则a＝a+3，b＝b+3，返回步骤(6)，否则进入步骤(10)；

(10)将循环八次得到的去噪图像DI(1)，...，DI(8)进行平均，得到最终的去噪结果SI。

本发明与现有技术相比具有以下优点：

(a)本发明利用Treelet变换估计非局部均值的平滑因子，进行非局部均值预滤波，滤除了图像中的大噪声，保持了图像细节，减少了噪声残留；

(b)本发明利用Treelet分解的差变量提取含噪数据的噪声信息，得到该数据的主要噪声信息表示，提出了去噪的最优阈值，对含噪数据进行软阈值去噪，降低了计算复杂度，提高了去噪图像的清楚度；

(c)本发明利用Cycle Spinning在含噪图像的Contourlet分解过程中对图像进行了循环平移，改善了Contourlet不具有平移不变的缺点，有效的减少了去噪图像的伪Gibbs现象。

附图说明

图1是本发明的实现流程图；

图2是本发明使用的未加噪图像和加噪图像；

图3是采用背景技术中提到的三种方法和本发明方法去噪后的图像。

具体实施方式

参照图1，本发明的实现步骤如下：

步骤1，由大小为A×B的含噪图像NI求出自协方差矩阵S_s和自相关系数矩阵M_s。

对于大小为A×B的含噪图像NI，用第i列的像素灰度构成A×1维的列向量NI_i，A是NI的行数，i取值为1，....，B，B是NI的列数，为了判断列向量之间的相关性，求NI的自协方差矩阵S_s：

$S_s=\frac{1}{B}\underset{i=1}{\overset{B}{\mathrm{\Σ}}}({\mathrm{NI}}_i-\stackrel{\‾}{\mathrm{NI}}){({\mathrm{NI}}_i-\stackrel{\‾}{\mathrm{NI}})}^T$

其中是所有列向量的平均，B是NI的列数，上标^T为矩阵转置；

将自协方差矩阵进行归一化，得到NI的自相关系数矩阵M_s，其中M_s在点(si，sj)处的系数值计算如下：

$M_s(\mathrm{si},\mathrm{sj})=\frac{S_s(\mathrm{si},\mathrm{sj})}{\sqrt{S_s(\mathrm{si},\mathrm{si})S_s(\mathrm{sj},\mathrm{sj})}}$

其中M_s(si，sj)为M_s在点(si，sj)处的系数值，S_s(si，sj)为S_s在点(si，sj)处的系数值，S_s(si，si)为S_s在点(si，si)处的系数值，S_s(sj，sj)为S_s在点(sj，sj)处的系数值，si取值为1，...，B，sj取值为1，...，B，为开根号操作。

步骤2，设Treelet分解的最高层数为L，每次分解所在层为lev，lev取值为1，.....，L，在lev＝1层，即Treelet的最底层，令Treelet的初始输入S⁽⁰⁾＝S_s，M⁽⁰⁾＝M_s，X⁽⁰⁾＝NI，L＝B-1，B是NI的列数。

步骤3，将上述输入S⁽⁰⁾，M⁽⁰⁾，X⁽⁰⁾和L代入Treelet，对X^(lev-1)进行Treelet分解计算出第lev层的差变量dif_lev，得到差变量矩阵{dif_lev}，lev取值为1，...，L。

(3a)令lev＝1，...，L，由相关系数矩阵M^(lev-1)找出最相似的两个变量η和μ：

(η，μ)＝argmaxM^(lev-1)(ii，jj)，ii＜jj

其中M^(lev-1)(ii，jj)为M^(lev-1)在点(ii，jj)处的系数值，ii取值为1，...，L+1，jj取值为1，...，L+1，η和μ分别代表Treelet在第lev-1层分解时的第一和第二主成分，取为令M^(lev-1)(ii，jj)最大的ii和jj；

(3b)求出雅克比旋转矩阵J：

$J=\left[\begin{array}{ccccccc}1& L& 0& L& 0& L& 0\\ M& O& M& & M& & M\\ 0& L& \mathrm{cc}& L& -\mathrm{ss}& L& 0\\ M& & M& O& M& & M\\ 0& L& \mathrm{ss}& L& \mathrm{cc}& L& 0\\ M& & M& & M& O& M\\ 0& L& 0& L& 0& L& 1\end{array}\right]$

矩阵中cc＝cos(θ_lev)，ss＝sin(θ_lev)，其中旋转角度 ${\mathrm{\θ}}_{\mathrm{lev}}=\frac{1}{2}\arctan \frac{2*S^{(\mathrm{lev}-1)}(\mathrm{\η},\mathrm{\μ})}{S^{(\mathrm{lev}-1)}(\mathrm{\η},\mathrm{\η})-S^{(\mathrm{lev}-1)}(\mathrm{\μ},\mathrm{\μ})},$ |θ_lev|≤π/4，更新变量：X^(lev)＝J^TX^(lev-1)，更新协方差矩阵和相关系数矩阵：S^(lev)＝J^TS^(lev-1)J，M^(lev)＝J^TS^(lev-1)J(此处上标^T为矩阵转置)；

(3c)第lev层的差变量

此处代表X^(lev)第μ列系数值。

步骤4，将图像矩阵NI投影在步骤(3)所得的差变量矩阵{dif_lev}上，估计出非局部均值滤波的平滑因子h，然后对含噪图像NI进行搜索窗为7×7，相似窗为3×3，平滑因子为h的非局部均值滤波，得到滤波后图像LI。

平滑因子h计算如下：

$h={(\frac{1}{3}*\mathrm{mean}\left(\right|I*{\mathrm{Difs}}^T\left|\right))}^2$

其中NI为图像矩阵，Difs＝{dif_lev}，是由NI得到的差变量矩阵，上标^T为矩阵转置，||为取绝对值操作，mean()为取均值操作。

步骤5，设步骤6-9中所用的循环平移的初始值a＝0，b＝0，cycle＝1，a为平移的行数，b为平移的列数，cycle为平移的次数。

步骤6，对滤波后图像LI进行循环平移，下移a行，右移b列，得到平移后的图像PI，对PI进行Contourlet变换，其中变换参数拉普拉斯滤波器组选择“9-7”，方向滤波器组选择“pkva”，PI总共分解为3层，第一层分解16个高频子带，第二层分解8个高频子带，第三层分解8个高频子带，分别对各个高频子带进行步骤7-8的处理。

步骤7，由大小为P×N的高频子带Y求出自协方差矩阵

和自相关系数矩阵

得到Treelet变换的初始输入：

X⁽⁰⁾＝Y，L＝N-1，其中N为Y的列数，P是Y的行数。

对大小为P*N的高频子带，其第r列的Contourlet系数值构成了P*1维的列向量Y_r(r＝1，....，N)，则系数矩阵为Y＝{Y₁，...，Y_N}，为了判断列向量之间的相关性，求Y的自协方差矩阵：

$\widehat{\mathrm{\Σ}}=\frac{1}{N}\underset{r=1}{\overset{N}{\mathrm{\Σ}}}(Y_r-\stackrel{\‾}{Y}){(Y_r-\stackrel{\‾}{Y})}^T$

其中 $\stackrel{\‾}{Y}=\frac{1}{N}\underset{r=1}{\overset{N}{\mathrm{\Σ}}}{Y}_r;$

将自协方差矩阵进行归一化，得到X的自相关系数矩阵在点(ti，tj)处的系数值计算如下：

$\hat{M}(\mathrm{ti},\mathrm{tj})=\frac{\widehat{\mathrm{\Σ}}(\mathrm{ti},\mathrm{tj})}{\sqrt{\widehat{\mathrm{\Σ}}(\mathrm{ti},\mathrm{ti})\widehat{\mathrm{\Σ}}(\mathrm{tj},\mathrm{tj})}}$

其中

为在点(ti，tj)处的系数值，

为在点(ti，tj)处的系数值，

为

在点(ti，ti)处的系数值，

为

在点(tj，tj)处的系数值，ti取值为1，...，N，tj取值为1，...，N。

步骤8，将步骤(7)得到的Treelet的初始输入S⁽⁰⁾，M⁽⁰⁾，X⁽⁰⁾及L代入步骤(3)中，得到差变量矩阵{dif_lev}，根据{dif_lev}估计阈值Th并计算权重g，然后对高频子带Y依次进行软阈值去噪和独立点噪声去除。

8a)将高频子带Y投影在差变量矩阵上，估计阈值Th：

Th＝mean(|Y*Detail^T|)

其中差变量矩阵Detail＝{dif_lev}，上标^T为矩阵转置，||为取绝对值操作，mean()为取均值操作；

8b)计算高频子带Y的权重g：

$g=\frac{\sqrt{2\ln (P*N)}}{\sqrt{\ln (k+1)}}*(1-2*e\_r)*s$

这里P为Y的行数，N为Y的列数，k为子带所在的分解层，在1，2，3中取值；

为各子带与所在层之间的系数能量比，其中q表示子带为第k层的第q个子带，Q为第k层的子带个数，

为第k层第q个子带内的所有系数幅值的平方和；s＝argmin(GCV)为调整因子，按照广义交叉验证准则在集合{0.1，0.2，...，1}中选取一个具体元素，GCV为广义交叉验证的公式；

s的选取过程如下：

设ll为s的选取过程中循环的次数，G(ll)为第ll次循环中得到的GCV公式的值，令初始值s＝0.1，ll＝1，计算G(ll)：

$G\left(\mathrm{ll}\right)=\mathrm{GCV}\left(s\right)=\frac{\frac{1}{P*N}{\left|\right|Y-Y_s\left|\right|}^2}{{\left(\frac{{\mathrm{LY}}_s}{P*N}\right)}^2}$

其中P为高频子带Y的行数，N为高频子带Y的列数，Y_s＝sign(Y)max(0，|Y|-Th*g)，为Y经软阈值去噪后的系数矩阵，LY_s为Y_s中幅值为0的系数个数，|| ||²为取均方值操作；

令s＝s+0.1，ll＝ll+1，重新计算G(ll)，重复此过程直到s＝1，ll＝10，可得到G＝{G(1)，G(2)，....G(10)}，再由s＝argmin(G)选出相应的s；

8c)根据估计的阈值Th和权重g，对高频子带Y进行软阈值去噪，得到阈值去噪后的子带Y_s：

Y_s＝sign(Y)max(0，|Y|-Th*g)

式中|Y|为Y的绝对值矩阵；

8d)对Y_s进行独立点噪声去除，得到最终的去噪子带Y_z：

$Y_z(\mathrm{li},\mathrm{lj})=\left\{\begin{array}{c}0,\mathrm{neigh}=0\\ Y_s(\mathrm{li},\mathrm{lj}),\mathrm{else}\end{array}\right.$

其中Y_s(li，lj)为Y_s在点(li，lj)处的Contourlet系数值，Λ为点(li，lj)的8邻域， $\mathrm{neigh}=\underset{(\mathrm{li},\mathrm{lj})\∈\mathrm{\Λ}}{\mathrm{\Σ}}\left|Y_s(\mathrm{li},\mathrm{lj})\right|$ 为Y_s在邻域Λ内的系数绝对值之和。

步骤9，将所有经过步骤(7)-(8)处理过的高频子带进行Contourlet逆变换，得到平移后的去噪图像FI，将FI左移b列，上移a行，得到去噪图像DI(cycle)，令cycle＝cycle+1，若cycle＜＝8，则a＝a+3，b＝b+3，返回步骤(6)，否则进入步骤(10)。

步骤10，将循环八次得到的去噪图像DI(1)，...，DI(8)进行平均，得到最终去噪结果SI。

$\mathrm{SI}=\frac{1}{8}*\underset{\mathrm{cycle}=1}{\overset{8}{\mathrm{\Σ}}}\mathrm{DI}\left(\mathrm{cycle}\right)$

本发明的效果可通过以下仿真实验结果与分析进一步说明：

1.实验数据

实验所使用的输入图像如图2所示，其中图2(a)为未加噪的Lena图像，图2(b)为未加噪的Barbara图像，图2(c)为未加噪的Boats图像，图2(d)为加入噪声水平为20的高斯随机白噪声的lena图像，四幅图像均为灰度图像，大小均为512×512，灰度级为256。

2.实验评价指标

图像去噪效果的评价分为主观和客观两个方面。在主观上评价一幅图像去噪效果的优劣主要是通过人眼的视觉特性来衡量，图像质量非常好，感觉很清晰则去噪效果好，反之效果较差。在客观上评价一幅图像的去噪效果本发明采用峰值信噪比PSNR来衡量。令真实图像为V，去噪后的图像为V’，V(t)为V中第t个像素的灰度值，V’(t)为V’中第t个像素的灰度值，Ω为图像中的像素集合，V_max＝max{V(t)，t∈Ω}。则峰值信噪比按如下公式计算：

$\mathrm{PSNR}=10\log \left[\frac{V_{\max }^2}{\frac{1}{\left|\mathrm{\Ω}\right|}\underset{t\∈\mathrm{\Ω}}{\mathrm{\Σ}}{[V^{\′}\left(t\right)-V\left(t\right)]}^2}\right]$

评价去噪效果的方法是在一幅清晰的图像上加入随机的高斯白噪声，然后在加噪的图像上进行去噪的实验。

3.对比实验

本发明使用的对比实验方法如下所述：

对比方法1：ZF Zhou等在文献“Contourlet-based Image Denoising Algorithm usingAdaptive Windows.ICIEA，2009：3654-3657.”中提出的Contourlet域自适应窗口Wiener滤波的图像去噪方法，简称C-W。

对比方法2：D D.-Y.Po，M N Do等在文献“Directional multiscale modeling ofimages using the Contourlet transform.IEEE Transactions on Image Processing，2006，15(6)：1610-1620.”中提出的图像Contourlet域多尺度HMT建模去噪方法，简称C-HMT。

对比方法3：梁栋，沈敏等在文献“一种基于Contourlet递归Cycle Spinning的图像去噪方法.电子学报，2005，33(11)：2044-2046.”中提出的基于Cycle Spinning的Contourlet硬阈值去噪方法，简称C-CS。

4.实验结果与分析

为了验证本发明方法的有效性，采用上述三种对比实验方法和本发明方法分别对图2(d)所示的图像进行去噪，去噪结果如图3所示，其中图3(a)是用现有C-W方法得到的去噪结果，图3(b)是用现有C-HMT方法得到的去噪结果，图3(c)是用现有C-CS方法得到的去噪结果，图3(d)是用本发明方法得到的去噪结果，从图3可以看出，本发明的主观去噪效果明显优于现有的三种去噪方法。

表1列出了噪声水平为10～50时四种方法的实验评价指标，表中sigma表示噪声水平。

表1不同的去噪方法对不同图像在不同噪声水平下的峰值信噪比PSNR

表1中的粗体部分为4组数据最大值。

从表1的PSNR可以直观的看出，本发明方法在各种噪声水平下的PSNR都优于其他三种方法，充分验证了本发明方法对自然图像去噪的有效性。

Claims (4)

1.一种基于Treelet的Contourlet域图像去噪方法，包括如下步骤：

(1)由大小为A×B的含噪图像矩阵NI求出自协方差矩阵S_s和自相关系数矩阵M_s；

(2)设Treelet分解的最高层数为L，每次分解所在层为lev，lev取值为1，.....，L，在lev＝1层，即Treelet的最底层，令Treelet的初始输入S⁽⁰⁾＝S_s，M⁽⁰⁾＝M_s，X⁽⁰⁾＝NI，L＝B-1，B是NI的列数，A是NI的行数；

(3)将上述输入S⁽⁰⁾，M⁽⁰⁾，X⁽⁰⁾和L代入Treelet，对X⁽⁰⁾进行Treelet变换分解到最高层，得到差变量矩阵{dif_lev}，dif_lev为Treelet在第i层分解得到的差变量，lev取值为1，...，L；

(4)将图像矩阵NI投影在步骤(3)所得的差变量矩阵{dif_i}上，估计出非局部均值滤波的平滑因子h，然后对含噪图像NI进行搜索窗为7×7，相似窗为3×3，平滑因子为h的非局部均值滤波，得到滤波后图像LI；

(5)设步骤(6)-(9)中所用的循环平移的初始值a＝0，b＝0，cycle＝1，a为平移的行数，b为平移的列数，cycle为平移的次数；

(6)对滤波后图像LI进行循环平移，下移a行，右移b列，得到平移后的图像PI，对PI进行Contourlet变换，其中变换的参数拉普拉斯滤波器组选择“9-7”，方向滤波器组选择“pkva”，PI总共分解为3层，第一层分解16个高频子带，第二层分解8个高频子带，第三层分解8个高频子带，分别对各个高频子带进行步骤(7)-(8)的处理；

(7)由大小为P×N的高频子带Y求出自协方差矩阵

和自相关系数矩阵

得到Treelet变换的初始输入：

X⁽⁰⁾＝Y，L＝N-1，其中N为Y的列数，P是Y的行数；

(8)将步骤(7)得到的Treelet的初始输入S⁽⁰⁾，M⁽⁰⁾，X⁽⁰⁾及L代入步骤(3)中，得到差变量矩阵{dif_i}，根据{dif_i}估计阈值Th并计算权重g，然后对高频子带Y依次进行软阈值去噪和独立点噪声去除；

(9)将所有经过步骤(7)-(8)处理过的高频子带进行Contourlet逆变换，得到平移后的去噪图像FI，将FI左移b列，上移a行，得到去噪图像DI(cycle)，令cycle＝cycle+1，若cycle＜＝8，则a＝a+3，b＝b+3，返回步骤(6)，否则进入步骤(10)；

(10)将循环八次得到的去噪图像DI(1)，...，DI(8)进行平均，得到最终的去噪结果SI。

2.根据权利要求1所述的图像去噪方法，其中步骤(4)所述的非局部均值滤波的平滑因子h，计算如下：

$h={(\frac{1}{3}*\mathrm{mean}\left(\right|\mathrm{NI}*{\mathrm{Difs}}^T\left|\right))}^2$

其中NI为图像矩阵，Difs＝{dif_lev}，是由NI得到的差变量矩阵，上标^T为矩阵转置，||为取绝对值操作，mean()为取均值操作。

3.根据权利要求1所述的图像去噪方法，其中步骤(8)所述的阈值Th和权重g计算如下：

8a)将高频子带Y投影在差变量矩阵上，估计阈值Th：

Th＝mean(|Y*Detail^T|)

其中Detail＝{dif_lev}，是由子带Y得到的差变量矩阵；

8b)计算高频子带Y的权重g：

$g=\frac{\sqrt{2\ln (P*N)}}{\sqrt{\ln (k+1)}}*(1-2*e\_r)*s$

这里P为Y的行数，N为Y的列数，k为子带所在的分解层，在1，2，3中取值；

为子带Y与所在层之间的系数能量比，其中q表示子带Y为第k层的第q个子带，Q为第k层的子带个数，为第k层第q个子带内的所有系数幅值的平方和；s＝argmin(GCV)为调整因子，按照广义交叉验证准则在集合{0.1，0.2，...，1}中选取一个具体元素，GCV为广义交叉验证的公式；

8c)根据估计的阈值Th和权重g，对高频子带Y进行软阈值去噪，得到阈值去噪后的子带Y_s：

Y_s＝sign(Y)max(0，|Y|-Th*g)

式中|Y|为Y的绝对值矩阵；

8d)对Y_s进行独立点噪声去除，得到最终的去噪子带Y_z；

Y_z在点(li，lj)处的系数值计算如下：

$Y_z(\mathrm{li},\mathrm{lj})=\left\{\begin{array}{c}0,\mathrm{neigh}=0\\ Y_s(\mathrm{li},\mathrm{lj}),\mathrm{else}\end{array}\right.$

其中Y_s(li，lj)为子带Y_s在点(li，lj)处的系数值，Λ为点(li，lj)的8邻域， $\mathrm{neigh}=\underset{(\mathrm{li},\mathrm{lj})\∈\mathrm{\Λ}}{\mathrm{\Σ}}\left|Y_s(\mathrm{li},\mathrm{lj})\right|$ 为Y_s在邻域Λ内的系数绝对值之和。

4.根据权利要求3所述的权重g的计算方法，其中步骤(8b)所涉及的调整因子s，是按照广义交叉验证准则在集合{0.1，0.2，...，1}中选取的一个具体元素，选取过程如下：

8b1)设ll为s的选取过程中循环的次数，G(ll)为第ll次循环中得到的GCV公式的值，令初始值s＝0.1，ll＝1，计算G(ll)：

$G\left(\mathrm{ll}\right)=\mathrm{GCV}\left(s\right)=\frac{\frac{1}{P*N}{\left|\right|Y-Y_s\left|\right|}^2}{{\left(\frac{{\mathrm{LY}}_s}{P*N}\right)}^2}$

其中P为高频子带Y的行数，N为高频子带Y的列数，Y_s＝sign(Y)max(0，|Y|-Th*g)，为Y经软阈值去噪后的子带，LY_s为Y_s中幅值为0的系数个数，|| ||²为取均方值操作；

8b2)令s＝s+0.1，ll＝ll+1，返回(8b.1)重新计算G(ll)，重复此过程直到s＝1，ll＝10，可得到G＝{G(1)，G(2)，....G(10)}，再由s＝argmin(G)选出相应的s。

Images