CN102663686B - 基于Treelet变换和高斯尺度混合模型的图像去噪方法 - Google Patents

Abstract

本发明公开了一种基于Treelet变换和高斯尺度混合模型的图像去噪方法，主要解决现有去噪方法对被高斯白噪声腐蚀的自然图像去噪效果不佳的问题。其实现步骤如下，(1)输入一幅待去噪图像；(2)对图像块进行分类；(3)组成矩阵；(4)获得基矩阵；(5)投影；(6)估计无噪声系数；(7)计算去噪后图像块矩阵；(8)判断是否处理完图像内所有图像块，若处理完，则进行步骤(9)，否则转入步骤(3)；(9)归一化所有图像块，结果输出。本发明具有对含有高斯白噪声的自然图像去噪效果好的优点，能恢复出图像原有的特征，可用于图像分割、目标识别，变换检测等对图像的预处理。

Description

基于Treelet变换和高斯尺度混合模型的图像去噪方法

技术领域

本发明属于图像处理技术领域，更进一步涉及一种自然图像的Treelet变换和高斯尺度混合模型的去噪方法。该方法可用于开展森林资源调查、灾害评估、城市规划、医学影像、天文学影像等领域的数字图像预处理。

背景技术

由于受各种条件的限制，图像在获取、编码、传输过程中会受到的噪声影响，这给图像分割，目标识别等其他一些后续处理工作带来了不利，因此从噪声图像中尽量恢复不含噪声的图像很有必要。图像去噪解决了图像受到噪声干扰图像质量下降的问题，抑制了噪声影响，提高了图像质量，是图像处理和计算机视觉中非常重要的基本问题。

在过去的几十年里，人们提出了许多图像去噪方法，其中频域滤波是一种非常重要的去噪方法。频域滤波就是对原始图像进行多尺度变换，再用阈值，或者对系数分布进行建模的方法对变换后的系数进行处理。常见的频域滤波方法有小波域的阈值收缩方法，小波域的高斯尺度混合方法，还有一些基于后小波的图像去噪方法等。

Javier Portilla等人在“Image Denoising using Scale Mixtures of Gaussians inthe Wavelet Domain.IEEE Transactions on Image Processing，vol.12，no.11，Nov.2003，1338-1351.”中提出了一种小波域的高斯尺度混合模型的去噪方法。该方法首先对图像进行多分辨率的金字塔分解，再将相邻尺度内对应的系数和当前尺度内的系数的分布用高斯尺度混合模型(Gaussian Scale Mixture，GSM)来拟合，GSM由一个高斯向量和一个隐的正尺度缩放因子这两个独立随机变量的乘积组成，用Bayes最小二乘估计去噪后小波系数，最后进行反变换就得到最终的去噪结果。该方法的不足是，对图像进行金字塔分解需要足够多的方向和尺度，导致了计算复杂度的增加。而且，该方法对变化域系数进行处理，去噪的效结果中存在很多划痕，伪吉布斯效应十分明显，导致去噪效果不理想。

西安电子科技大学在其申请的专利“基于Treelet变换和非局部均值的图像去噪方法”(专利申请号201110001952.1，公开号CN102063708A)中公开了一种利用相似图像块灰度值Treelet变换的系数来计算计算图像块之间的相似度的非局部均值去噪方法。该方法在高噪声情况下能提高去噪效果，但是该方法存在的不足是，对图像中细节信息保持不好，去噪后图像过平滑，导致处理后的图像结构信息(边缘、纹理、点)会被模糊或滤除。

发明内容

本发明的目的在于克服现有技术中的不足，提出一种基于Treelet变换和高斯尺度混合模型的图像去噪方法。将图像中结构和纹理相似的图像块归为一类，再对每一类的图像块的灰度特征向量组成的矩阵进行Treelet变换。再用高斯尺度混合模型来拟合每一类图像块Treelet变换后的系数分布情况，用Bayes方法估计无噪声系数。

为实现上述目的，本发明的去噪方法主要包括以下步骤：

(1)输入一幅待去噪图像；

(2)对图像块进行分类

2a)在待去噪图像中，以任一像素点为中心，以固定长度为边长，确定一个正方形的图像块，用相同方法，对图像中所有像素都做同样操作；

2b)用k-means方法对步骤2a)中取得的所有图像块进行分类；

(3)将同一类的所有图像块的像素灰度向量组成矩阵Y；

(4)获得基矩阵B

4a)按照下式，计算灰度特征向量组成的矩阵Y的协方差矩阵和相关系数矩阵

${\widehat{\mathrm{\Σ}}}_{\mathrm{ij}}=E\left[(y_i-{\mathrm{Ey}}_i){(y_j-{\mathrm{Ey}}_j)}^T\right]$

${\hat{M}}_{\mathrm{ij}}={\widehat{\mathrm{\Σ}}}_{\mathrm{ij}}/\sqrt{{\widehat{\mathrm{\Σ}}}_{\mathrm{ii}}{\widehat{\mathrm{\Σ}}}_{\mathrm{jj}}}$

其中，表示协方差矩阵中的元素，i和j分别表示元素在协方差矩阵中的行和列的位置，E表示数学期望，y_i和y_j表示步骤(3)中任意两个列向量，T表示矩阵转置操作，表示相关系数矩阵中的元素，和分别为y_i和y_j的自协方差；

4b)初始化Treelet分解的层数l＝1，...，L；L为最大分解层数，其值为48；当l＝0时，初始化相似度矩阵为初始化基矩阵B₀为与协方差矩阵相同大小的单位矩阵，初始化和变量的下标集：δ＝{1，2，…，p}，其中p为矩阵Y的维数；

4c)按照下式，找出相似度矩阵中最相似的两个变量：

$(\mathrm{\α},\mathrm{\β})=\arg \max {\hat{M}}^{(l-1)}$

其中，α和β分别表示矩阵中最相似的两个变量的行和列位置，arg max表示在矩阵中寻找最大值，是l-1层相似度矩阵；

4d)对步骤4c)中得到的α和β两个变量由式下式计算雅克比旋转角度θ_l：

${\widehat{\mathrm{\Σ}}}_{\mathrm{\α\β}}^{\left(l\right)}={\widehat{\mathrm{\Σ}}}_{\mathrm{\β\α}}^{\left(l\right)}=0$

其中，表示第l层协方差矩阵，α和β分别表示步骤4c)计算得到的最相似两个变量的行和列位置；

4e)按下式计算雅克比旋转矩阵J：

其中，J表示雅克比旋转矩阵，cos和sin分别表示余弦和正弦函数，θ_l为雅克比旋转角度；

按照下式更新基矩阵B_l、和

B_l＝B_l-1J

${\widehat{\mathrm{\Σ}}}^{\left(l\right)}=J^T{\widehat{\mathrm{\Σ}}}^{(l-1)}J$

${\hat{M}}^{\left(l\right)}=J^T{\hat{M}}^{(l-1)}J$

其中，B_l第l层级的基矩阵，l表示分解层数，B_l-1第l层级的基矩阵，J为雅克比旋转矩阵，为第l层协方差矩阵，J^T为雅克比旋转矩阵J的转置，为l层相似度矩阵；

4f)重复步骤4b)、步骤4c)、步骤4d)、步骤4e)直至分解层数l＝L层，L为分解最大层数，其值为48，得到Treelet分解的基矩阵B：

B＝[φ_l，ψ₁，...ψ_l]^T

其中，B为基矩阵，φ_l表示矩阵B_l的第α列向量，ψ₁表示矩阵B_l的第β列向量，l表示分解层数，T表示矩阵转置；

(5)投影

5a)按照下式计算矩阵Y去均值的矩阵

$\stackrel{\‾}{Y}=Y-E\left[Y\right]$

其中，为Y去均值的矩阵，Y为步骤(3)得到的矩阵，E[Y]为Y的行均值矩阵；

5b)按照下式将步骤5a)得到的图像块去均值矩阵投影到步骤(4)得到的基矩阵B上：

$c=B\stackrel{\‾}{Y}$

其中，c是矩阵Y在基矩阵B上的投影系数，B是步骤(4)获得的基矩阵，是矩阵Y去均值的结果；

(6)估计无噪声系数

6a)根据下面公式计算含噪系数协方差矩阵C_c：

$C_c=\frac{1}{N}c*c$

其中，C_c向表示协方差矩阵，N表示投影系数c的维数，c表示投影系数；

按照下式估计无噪声系数的协方差矩阵C_u：

C_u＝(C_c-C_n)/E[z]

其中，C_u表示无噪声系数的协方差矩阵，C_c表示含噪系数的协方差矩阵，C_n是噪声的协方差矩阵，E表示数学期望，z表示高斯尺度混合模型中尺度随机因子；

6b)根据下式估计条件期望E[x|c，z]和p(c|z)：

E[x|c，z]＝zC_u(zC_u+C_n)^-1c

$p\left(c\right|z)=\frac{\exp ({-c}^T{({\mathrm{zC}}_u+C_n)}^{-1}c/2)}{\sqrt{{\left(2\mathrm{\π}\right)}^N|{\mathrm{zC}}_u+C_n|}}$

其中，E[x|c，z]是无噪声系数x关于含噪系数c和尺度因子z的条件期望，x表示无噪声系数，c是含噪系数，z是尺度随机因子，C_u和C_n分别是无噪声系数和噪声的协方差矩阵；p(c|z)表示含噪系数c关于z的条件概率，exp()表示指数，T表示矩阵转置，N表示信号的维数，||表示绝对值；

6c)根据下面公式估计p(z|c)：

$p\left(z\right|c)=\frac{p\left(c\right|z)p_z\left(z\right)}{{\∫}_0^{\∞}p\left(c\right|\mathrm{\α})p_z\left(\mathrm{\α}\right)\mathrm{d\α}}$

其中，p(z|c)表示表示尺度因子z关于含噪声系数c的条件概率，p(c|z)表示步骤6b)得到的含噪系数c关于z的条件概率，p_z(z)表示z的分布，表示积分操作，p(c|α)表示含噪系数c关于α的条件概率，p_z(α)表示α的概率密度；

6d)用下式估计无噪声系数：

$E\left[\hat{x}\right|c]={\∫}_0^{\∞}p\left(z\right|c)E[x|c,z]\mathrm{dz}$

其中，表示估计的无噪声系数关于含噪系数c的条件期望，表示估计得到的无噪声系数，c表示含噪系数，表示积分操作，p(z|c)是步骤6c)计算得到的z关于含噪系数c的条件概率，z表示尺度因子，c表示含噪声系数，E[x|c，z]是步骤6b)计算得到无噪声系数x关于含噪系数c和尺度因子z的条件期望，x表示无噪声系数；

(7)按照下式计算去噪后图像块矩阵

$\hat{Y}=B^T\hat{x}+E\left[Y\right]$

其中，是去噪后图像块矩阵，B^T是基矩阵的转置，表示估计的无噪声系数，E[Y]是Y的行均值矩阵；

(8)判断是否处理完图像内所有图像块，若处理完，则进行步骤(9)，否则转入步骤(3)；

(9)对去噪后的所有图像块灰度进行归一化处理，得到去噪结果。

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

第一，本发明将待去噪图像中的图像块进行分类，对相同类别的图像块组成的矩阵进行Treelet变化，由于Treelet变换对噪声具有鲁棒性，克服了现有技术受噪声影响较大的缺点，使得本发明的去噪结果对图像中的细节信息得到较为完整的保持。

第二，本发明在Treelet变化后，用高斯尺度混合模型对系数进行建模，能比较好的拟合Treelet变换后系数的分布情况。用贝叶斯最小均方准则对无噪声系数进行估计，克服了现有技术去噪结果中存在划痕的缺点，较好的保持了图像中的细节信息，平滑同质区域的噪声，提高了图像的去噪效果。

第三，本发明选择中心图像块时进行下采样操作，克服了现有技术对像素点处理存在的计算时间长的缺点，减少了运算时间，实现过程简单，且可以并行实现。

附图说明

图1为本发明的流程图；

图2为本发明与现有技术BLS_GSM的去噪效果比较图；

图3为本发明与输入的噪声标准差为20的图像比较图。

具体实施方式

下面结合附图1对本发明实现的步骤作进一步的详细描述。

步骤1，输入一幅待去噪图像。

步骤2，对图像块进行分类。

2a)在待去噪图像中，以任一像素点为中心，以固定长度为边长，确定一个正方形的图像块，用相同方法，对图像中所有像素都做同样操作，本发明实施例中选取7×7大小的图像块。

2b)从图像中选取的所有图像块随机选取q个图像块灰度矩阵作为聚类中心，本发明实施例中q取值为70。

2c)按下式分别计算图像块到所有聚类中心的距离d(y_p，m_k)：

$d(y_p,m_k)={\left|\right|y_p-m_k\left|\right|}_2^2$

其中，d(y_p，m_k)表示第p个图像块的灰度矩阵，m_k表示第k个聚类中心的灰度矩阵，表示二范数的平方。

根据图像块到聚类中心的距离，将它们分配给与其距离最近聚类中心所代表的类别。

2d)按照下式更新每个聚类中心：

$m_k=\frac{1}{\left|S_k\right|}\underset{y_p\∈S_k}{\mathrm{\Σ}}{y}_p$

其中，m_k表示第k个聚类中心的灰度矩阵，S_k表示第k个类别的图像块集合，|s_k|表示第k个类别中图像块的个数，y_p∈s_k表示属于第k类的图像块。

2e)判断聚类中心是否变化，若变化，则转入步骤2c)，否则转入步骤2f)。

2f)输出聚类结果。

步骤3，将属于同一类的图像块的灰度向量组成一个矩阵Y。

步骤4，对矩阵Y进行Treelet变换，得到基B。

4a)按照下式，计算灰度特征向量组成的矩阵Y的协方差矩阵和相关系数矩阵

${\widehat{\mathrm{\Σ}}}_{\mathrm{ij}}=E\left[(y_i-{\mathrm{Ey}}_i){(y_j-{\mathrm{Ey}}_j)}^T\right]$

${\hat{M}}_{\mathrm{ij}}={\widehat{\mathrm{\Σ}}}_{\mathrm{ij}}/\sqrt{{\widehat{\mathrm{\Σ}}}_{\mathrm{ii}}{\widehat{\mathrm{\Σ}}}_{\mathrm{jj}}}$

其中，表示协方差矩阵中的元素，E表示数学期望，y_i和y_j表示步骤(3)中任意两个列向量，i和j分别表示列向量y_i和y_j的在矩阵中的位置。T表示矩阵转置操作，表示相关系数矩阵中的元素，和分别为y_i和y_j的协方差。

4b)初始化Treelet分解的层数l＝1，...，L，本发明实施例中选取分解的最高层数L＝48。当l＝0时，初始化相似度矩阵为初始化基矩阵B₀为与协方差矩阵相同大小的单位矩阵；初始化和变量的下标集：δ＝{1，2，…，p}，其中p为矩阵Y的维数。

4c)按照下式，找出相似度矩阵中最相似的两个变量：

$(\mathrm{\α},\mathrm{\β})=\arg \max {\hat{M}}^{(l-1)}$

其中，α和β分别表示矩阵中最相似的两个变量的位置，arg max表示在矩阵中寻找最大值，是l-1层相似度矩阵。

4d)对步骤4c)中得到的两个变量进行主成分变换，得到去相关后的第一主成分和第二主成分。这里，主成分变换是通过雅克比旋转实现。旋转角度θ_l由且θ_l≤π/4计算得到，其中协方差矩阵雅克比旋转矩阵J由下式计算：

其中，θ_l为旋转角度。

按照下式更新基矩阵B_l和

B_l＝B_l-1J

${\hat{M}}^{\left(l\right)}=J^T{\hat{M}}^{(l-1)}J$

其中，B_l第l层级的基矩阵，B_l-1第l层级的基矩阵，J为雅克比旋转矩阵，J^T为雅克比旋转矩阵的转置，为l层相似度矩阵。

4e)旋转角度后使即α对应的变量为主成分分析后的第一主成分，β对应的变量为第二主成分。尺度向量φ_l和细节向量ψ_l分别为此分解层级下的基矩阵B_l的第α和β列。

4f)重复步骤(4b)至步骤(4d)直至l＝L层，得到第L层Treelet分解的基矩阵，记为B：

B＝[φ_l，ψ₁，...ψ_l]^T

其中，φ_l表示最高层的尺度，ψ_l表示第l层的细节向量，T表示矩阵转置操作。

步骤5，将步骤(3)得到的图像块矩阵Y投影到步骤(4)得到的基矩阵B上，得到了Treelet变换后的系数c。

5a)按照下式计算相似图像块矩阵Y去均值的矩阵

$\stackrel{\‾}{Y}=Y-E\left[Y\right]$

其中，为Y去均值的矩阵，Y为步骤(3)得到的矩阵，E[Y]为Y的行均值矩阵。

5b)将步骤5a)得到的图像块去均值矩阵投影到步骤(4)得到的基矩阵B上。

按照下式将矩阵投影到基矩阵B上：

$c=B\stackrel{\‾}{Y}$

c是图像块在基矩阵B上的投影系数，B是获得的基矩阵，是图像块矩阵去均值的结果。

步骤6，估计无噪声系数。

6a)根据下面公式计算含噪系数协方差矩阵C_c：

$C_c=\frac{1}{N}c*c$

其中，C_c向表示协方差矩阵，N表示投影系数c的维数，c表示投影系数；

按照下式估计无噪声系数的协方差矩阵C_u：

C_u＝(C_c-C_n)/E[z]

其中，C_u表示无噪声系数的协方差矩阵，C_c表示含噪系数的协方差矩阵，C_n是噪声的协方差矩阵，E表示数学期望，z表示高斯尺度混合模型中尺度随机因子；

6b)根据下式估计条件期望E[x|c，z]和p(c|z)：

E[x|c，z]＝zC_u(zC_u+C_n)^-1c

$p\left(c\right|z)=\frac{\exp ({-c}^T{(z{C}_u+C_n)}^{-1}c/2)}{\sqrt{{\left(2\mathrm{\π}\right)}^N|z{C}_u+C_n|}}$

其中，E[x|c，z]是无噪声系数x关于含噪系数c和尺度因子z的条件期望，x表示无噪声系数，c是含噪系数，z是尺度随机因子，C_u和C_n分别是无噪声系数和噪声的协方差矩阵；p(c|z)表示含噪系数c关于z的条件概率，exp()表示指数，T表示矩阵转置，N表示信号的维数，||表示绝对值；

6c)根据下面公式估计p(z|c)：

$p\left(z\right|c)=\frac{p\left(c\right|z)p_z\left(z\right)}{{\∫}_0^{\∞}p\left(c\right|\mathrm{\α})p_z\left(\mathrm{\α}\right)\mathrm{d\α}}$

其中，p(z|c)表示表示尺度因子z关于含噪声系数c的条件概率，p(c|z)表示步骤6b)得到的含噪系数c关于z的条件概率，p_z(z)表示z的分布，表示积分操作，p(c|α)表示含噪系数c关于α的条件概率，p_z(α)表示α的概率密度；

6d)用下式估计无噪声系数：

$E\left[\hat{x}\right|c]={\∫}_0^{\∞}p\left(z\right|c)E[x|c,z]\mathrm{dz}$

其中，表示估计的无噪声系数关于含噪系数c的条件期望，表示估计得到的无噪声系数，c表示含噪系数，表示积分操作，p(z|c)是步骤6c)计算得到的z关于含噪系数c的条件概率，z表示尺度因子，c表示含噪声系数，E[x|c，z]是步骤6b)计算得到无噪声系数x关于含噪系数c和尺度因子z的条件期望，x表示无噪声系数。

步骤7，按照下式计算去噪后图像块矩阵

$\hat{Y}=B^T\hat{x}+E\left[Y\right]$

其中，是去噪后图像块矩阵，B^T是基矩阵的转置，表示估计的无噪声系数，E[Y]是Y的行均值矩阵。

步骤8，判断是否处理完图像内所有类别图像块，若处理完，则进行步骤(9)，否则转入步骤(3)。

步骤9，对去噪后的所有图像块灰度进行归一化处理，得到去噪结果。

为了减少算法运算时间，在取当前要处理的图像块的时候进行采样处理，每四个像素选取一个图像块。

下面结合附图2、附图3对本发明的效果做进一步的描述。

附图2为本发明与现有技术BLS_GSM的去噪效果比较图。

附图2(a)，附图2(b)和附图2(c)为本发明的去噪效果图，附图2(d)，附图2(e)和附图2(f)为现有技术基于BLS_GSM的去噪方法的效果图，附图2(d)，附图2(e)和附图2(f)去噪效果图来源于文献“Image Denoising using Scale Mixtures ofGaussians in the Wavelet Domain.IEEE Transactions on Image Processing，vol.12，no.11，Nov.2003，1338-1351.”。

附图3是本发明与输入的噪声标准差为20的图像的比较图。

附图3(a)，附图3(b)和附图3(c)是本发明输入的待去噪的图像，图像加了标准差为20的随机高斯白噪声，附图3(a)为Lena图像，附图3(b)为Barbara图像，附图3(c)为Cameraman图像，所有图像都是灰度图像，灰度级为256级，其中Lena和Barbara为512×512像素大小，Cameraman为256×256像素大小。附图3(d)，附图3(e)和附图3(f)是本发明的对附图3(a)，附图3(b)和附图3(c)三幅图像去噪效果图。

将本发明的去噪效果图和现有技术去噪效果图进行比较，对本发明的效果作进一步描述。从图中可以看出，本发明的去噪效果，噪声得到有效的去除而且细节信息也得到较完整的保留，而现有方法的去噪效果中，存在一些划痕，失真明显。

下表分别列出了本发明与现有技术基于BLS_GSM方法，对附图3(a)，附图3(b)和附图3(c)三幅待去噪图像去噪效果图的峰值信噪比。表中第一列分别表示三幅待去噪的图像：Lena、Barbara和Cameraman；第二列分别表示待去噪图像中噪声标准差为15、20、25和35；第三列分别表示BLS_GSM方法和本发明去噪结果图的峰值信噪比。

从上表可以看出，对于Barbara和Cameraman两幅图像在不同噪声情况下，本发明去噪结果PSNR值均高于BLS_GSM方法，显示出了本发明对同质区域噪声抑制能力和保护点、线、边缘等高频信息的能力，去噪结果图像中细节信息得到较为完整的保持。从附图3(d)，附图3(e)和附图3(f)可以看出本发明的去噪结果细节信息保持完整，图像清晰，失真小，说明本发明与基于BLS_GSM方法相比，在对细节信息，像线和纹理进行恢复时表现出了更好的效果，对同质区域的平滑性能很好，而且对图像的边缘细节信息也保持的比较完整。

Claims (4)

1.一种基于Treelet变换和高斯尺度混合模型的图像去噪方法，包括如下步骤：

(1)输入一幅待去噪图像；

(2)对图像块进行分类

2a)在待去噪图像中，以任一像素点为中心，以固定长度为边长，确定一个正方形的图像块，用相同方法，对图像中所有像素都做同样操作；

2b)用k-means方法对步骤2a)中取得的所有图像块进行分类；

(3)将同一类的所有图像块的像素灰度向量组成矩阵Y；

(4)获得基矩阵B

4a)按照下式，计算灰度特征向量组成的矩阵Y的协方差矩阵和相关系数矩阵

${\widehat{\mathrm{\Σ}}}_{\mathrm{ij}}=E\left[(y_i-{\mathrm{Ey}}_i){(y_j-{\mathrm{Ey}}_j)}^T\right]$

${\widehat{\mathrm{\Σ}}}_{\mathrm{ij}}={\widehat{\mathrm{\Σ}}}_{\mathrm{ij}}/\sqrt{{\widehat{\mathrm{\Σ}}}_{\mathrm{ii}}{\widehat{\mathrm{\Σ}}}_{\mathrm{jj}}}$

其中，表示协方差矩阵中的元素，i和j分别表示元素在协方差矩阵中的行和列的位置，E表示数学期望，y_i和y_j表示步骤(3)中任意两个列向量，T表示矩阵转置操作，表示相关系数矩阵中的元素，和分别为y_i和y_j的自协方差；

4b)初始化Treelet分解的层数l＝0，1，...L；L为最大分解层数，其值为48；当l＝0时，初始化相似度矩阵为初始化基矩阵B₀为与协方差矩阵相同大小的单位矩阵，初始化和变量的下标集：δ＝{1，2，…，p}，其中p为矩阵Y的维数；

4c)按照下式，找出相似度矩阵中最相似的两个变量：

$(\mathrm{\α},\mathrm{\β})=\arg \max {\hat{M}}^{(l-1)}$

其中，α和β分别表示矩阵中最相似的两个变量的列位置，arg max表示在矩阵中寻找最大值，是l-1层相似度矩阵；

4d)对步骤4c)中得到的α和β两个变量由下式计算雅克比旋转角度θ_l：

${\widehat{\mathrm{\Σ}}}_{\mathrm{\α\β}}^{\left(l\right)}={\widehat{\mathrm{\Σ}}}_{\mathrm{\β\α}}^{\left(l\right)}=0$

其中，表示第l层协方差矩阵，α和β分别表示步骤4c)计算得到的最相似两个变量的列位置；

4e)按下式计算雅克比旋转矩阵J：

其中，J表示雅克比旋转矩阵，cos和sin分别表示余弦和正弦函数，θ_l为雅克比旋转角度；

按照下式更新基矩阵B_l、和

B_l＝B_l-1J

${\widehat{\mathrm{\Σ}}}^{\left(l\right)}=J^T{\widehat{\mathrm{\Σ}}}^{(l-1)}J$

${\hat{M}}^{\left(l\right)}=J^T{\hat{M}}^{(l-1)}J$

其中，B_l第l层级的基矩阵，l表示分解层数，B_l-1第l层级的基矩阵，J为雅克比旋转矩阵，为第l层协方差矩阵，J^T为雅克比旋转矩阵J的转置，为l层相似度矩阵；

4f)重复步骤4b)、步骤4c)、步骤4d)、步骤4e)直至分解层数l＝L层，L为最大分解层数，其值为48，得到Treelet分解的基矩阵B：

B＝[φ_l，ψ₁，...ψ_l]^T

其中，B为基矩阵，φ_l表示矩阵B_l的第α列向量，ψ₁表示矩阵B_l的第β列向量，l表示分解层数，T表示矩阵转置；

(5)投影

5a)按照下式计算矩阵Y去均值的矩阵

$\stackrel{\‾}{Y}=Y-E\left[Y\right]$

其中，为Y去均值的矩阵，Y为步骤(3)得到的矩阵，E[Y]为Y的行均值矩阵；

5b)按照下式将步骤5a)得到的图像块去均值矩阵投影到步骤(4)得到的基矩阵B上：

$c=B\stackrel{\‾}{Y}$

其中，c是矩阵Y在基矩阵B上投影得到的含噪系数，B是步骤(4)获得的基矩阵，是矩阵Y去均值的结果；

(6)估计无噪声系数

6a)根据下面公式计算含噪系数协方差矩阵C_c：

$C_c=\frac{1}{N}c*c$

其中，C_c向表示协方差矩阵，N表示投影系数c的维数，c表示含噪系数；

按照下式估计无噪声系数的协方差矩阵C_u：

$C_u=(C_c-C_n)/E\left[z\right]$

其中，C_u表示无噪声系数的协方差矩阵，C_c表示含噪系数的协方差矩阵，C_n是噪声的协方差矩阵，E表示数学期望，z表示高斯尺度混合模型中尺度随机因子；

6b)根据下式估计条件期望E[x|c，z]和p(c|z)：

E[x|c，z]＝zC_u(zC_u+C_n)^-1c

$p\left(c\right|z)=\frac{\exp ({-c}^T{({\mathrm{zC}}_u+C_n)}^{-1}c/2)}{\sqrt{{\left(2\mathrm{\π}\right)}^N|{\mathrm{zC}}_u+C_n|}}$

其中，E[x|c，z]是无噪声系数x关于含噪系数c和尺度随机因子z的条件期望，x表示无噪声系数，c是含噪系数，z是尺度随机因子，C_u和C_n分别是无噪声系数和噪声的协方差矩阵；p(c|z)表示含噪系数c关于z的条件概率，exp()表示指数，T表示矩阵转置，N表示信号的维数，||表示绝对值；

6c)根据下面公式估计p(z|c)：

$p\left(z\right|c)=\frac{p\left(c\right|z)p_z\left(z\right)}{{\∫}_0^{\∞}p\left(c\right|\mathrm{\α})p_z\left(\mathrm{\α}\right)\mathrm{d\α}}$

其中，p(z|c)表示尺度随机因子z关于含噪系数c的条件概率，p(c|z)表示步骤6b)得到的含噪系数c关于z的条件概率，p_z(z)表示z的分布，表示积分操作，p(c|α)表示含噪系数c关于α的条件概率，p_z(α)表示α的概率密度；

6d)用下式估计无噪声系数：

$E\left[\hat{x}\right|c]={\∫}_0^{\∞}p\left(z\right|c)E[x|c,z]\mathrm{dz}$

其中，表示估计的无噪声系数关于含噪系数c的条件期望，表示估计得到的无噪声系数，c表示含噪系数，表示积分操作，p(z|c)是步骤6c)计算得到的z关于含噪系数c的条件概率，z表示尺度随机因子，c表示含噪声系数，E[x|c，z]是步骤6b)计算得到无噪声系数x关于含噪系数c和尺度随机因子z的条件期望，x表示无噪声系数；

(7)按照下式计算去噪后图像块矩阵

$\hat{Y}=B^T\hat{x}+E\left[Y\right]$

其中，是去噪后图像块矩阵，B^T是基矩阵的转置，表示估计的无噪声系数，E[Y]是Y的行均值矩阵；

(8)判断是否处理完图像内所有图像块，若处理完，则进行步骤(9)，否则转入步骤(3)；

(9)对去噪后的所有图像块灰度进行归一化处理，得到去噪结果。

2.根据权利要求1所述的基于Treelet变换和高斯尺度混合模型的图像去噪方法，其特征在于：步骤2a)中所述的固定长度为7个像素。

3.根据权利要求1所述的基于Treelet变换和高斯尺度混合模型的图像去噪方法，其特征在于：步骤2b)中所述的图像块分类的类别数为70。

4.根据权利要求1所述的基于Treelet变换和高斯尺度混合模型的图像去噪方法，其特征在于：步骤2b)中所述的k-means方法具体步骤如下：

第1步，随机确定70个图像块灰度向量作为初始聚类中心；

第2步，计算每一图像块到每一聚类中心的距离，分别将它们分配给与其距离最相近聚类中心所代表的类别；

第3步，更新每个聚类中心；

第4步，判断聚类中心是否变化，若变化，则转入第2步，否则转入第5步；

第5步，输出聚类结果。