CN104616266B - 一种基于广义自回归异方差模型的噪声方差估计方法 - Google Patents

Abstract

一种基于广义自回归异方差模型的噪声方差估计方法，包括以下步骤：步骤1，读取被噪声污染的含噪图像；步骤2，对含噪图像做非子采样轮廓波变换；步骤3，对步骤2中的每个高频子带系数矩阵做去均值滤波处理；步骤4，将去均值滤波处理后的高频子带系数矩阵转换为一维序列数据；步骤5，对一维序列数据建立自回归模型，并求得其残差序列；步骤6，对残差序列建立统计学模型；步骤7，根据步骤5中求得的残差序列和步骤6的统计学模型，采用极大似然估计的方法计算统计学模型的参数；步骤8，求得含噪图像中的噪声的方差。该方法可以提高噪声方差估计的准确性，且适用于各种噪声水平的退化图像，为去噪、复原、特征提取等后续图像处理提供支持。

Description

一种基于广义自回归异方差模型的噪声方差估计方法

技术领域

本发明涉及图像处理技术领域，特别是涉及一种基于广义自回归异方差模型的噪声方差估计方法。

背景技术

由于成像环境和传输信道等各种因素的影响，获得的数字图像常常被噪声污染。噪声不仅降低了图像的质量和视觉效果，还会影响后续各种图像处理和分析过程。比如：图像复原中正则化参数的选取、稀疏表示中平衡因子的取值、图像压缩中最佳量化次数的确定等都依赖于噪声方差的先验知识，噪声方差估计得准确与否会显著地影响图像去噪、复原、表示、压缩、分割、特征提取、目标识别等的性能。

现有的图像噪声方差估计的方法主要有三大类：一类是基于滤波的方法，这类方法先通过滤波来提取图像的高频信息，再从高频信息中来估计噪声的方差，然而，该类方法容易将图像的边缘和纹理等高频细节信息视为噪声，导致噪声方差的过估计。第二类是基于分块的方法，该方法依赖“同质块的局部方差可以作为该块的噪声方差的估计量”这一假设，先按一定准则把噪声图像分成很多小块，然后从诸多小块中寻找同质子块，再基于同质子块来估计图像噪声的方差，该类方法性能的好坏取决于如何科学有效的确定图像的同质块，并保证有足够多的同质块参与到噪声方差的估计中。第三类是基于小波域统计建模的方法，该类方法中最著名的是Donoho方法，Donoho方法用统计量作为噪声方差的估计，其中median(·)表示取中值,HH为含噪图像小波变换域的第一层对角高频子带系数矩阵，该方法在噪声污染严重时估计的比较准确，当噪声水平较低时，会出现过估计。

因此，针对现有技术中的存在问题，亟需提供一种准确性高、且适用于各种噪声水平的退化图像，能够为去噪、复原、特征提取等后续图像处理提供支持的噪声方差估计方法显得尤为重要。

发明内容

本发明的目的在于针对现有含噪图像噪声方差估计技术的缺陷，提出一种 基于非子采样轮廓波变换和广义自回归异方差模型的噪声方差估计方法，该方法可以提高噪声方差估计的准确性，且适用于各种噪声水平的退化图像，为去噪、复原、特征提取等后续图像处理提供支持。

本发明的目的通过以下技术方案实现：

提供一种基于广义自回归异方差模型的噪声方差估计方法，包括以下步骤：

步骤1，读取被噪声污染的含噪图像；

步骤2，对含噪图像做非子采样轮廓波变换；

步骤3，对步骤2中的每个高频子带系数矩阵做去均值滤波处理；

步骤4，将去均值滤波处理后的二维高频子带系数矩阵转换为一维序列数据；

步骤5，对一维序列数据建立自回归模型，并求得其残差序列；

步骤6，对残差序列建立统计学模型；

步骤7，根据步骤5中求得的残差序列和步骤6的统计学模型，采用极大似然估计的方法计算统计学模型的参数；

步骤8，求得含噪图像中的噪声的方差。

其中，所述步骤1的具体步骤为：

读取被噪声污染的含噪图像g_x,y＝f_x,y+n_x,y，

其中：f_x,y为清晰图像，n_x,y为噪声图像；含噪图像g_x,y、清晰图像f_x,y、噪声图像n_x,y的大小均为N×M；n_x,y服从均值为0、方差为σ²的正态分布；下标x,y分别表示图像的行坐标与列坐标，1≤x≤N,1≤y≤M；

所述步骤2的具体步骤为：

对含噪图像g_x,y做L级非子采样轮廓波变换，即NSCT，其中，在第i级的方向数为2^{i -1}，得到大小为N×M的低频子带系数矩阵和各级各方向上的高频子带系数矩阵高频子带系数矩阵的个数为j为在第i级上的方向数，L为大于2的自然数；

其中，为f_x,y的NSCT的系数矩阵，为n_x,y的NSCT的系数矩阵。

其中，所述步骤3的具体步骤为：

采用以下公式对每个高频子带系数矩阵去均值滤波处理，

其中，为系数矩阵去均值滤波处理后得到的矩阵。

其中，所述步骤4的具体步骤为：

按照从左到右、从上到下的顺序依次回形读取去均值滤波后的二维系数矩阵 中的元素，得到一维序列数据序列数据既包含有噪声图像的信息，又包含有清晰图像的信息。

其中，所述步骤5的具体步骤为：

对一维序列数据建立r阶自回归模型，r为大于3的自然数，

为均值为0的残差序列，

采用Gram-Schmidt正交法确定自回归模型的阶数r，采用最小二乘估计的方法求得自回归模型的自回归系数进而得到均值为0的残差序列

所述残差序列包括两部分：其中，为清晰图像的残差信息，为噪声图像的残差信息。

其中，所述步骤6的具体步骤为：

用均值为0、方差为(σ^(i,j))²的正态分布模型对噪声图像的残差信息建模，得到模型：

用(1,1)阶广义自回归异方差模型GARCH(1,1)对清晰图像的残差信息建模，得到模型：

其中：α₀＞0,α₁≥0,β₁≥0，α₁+β₁＜1；

从而，残差序列服从的统计学模型为：

所述统计学模型中包含的参数有

其中，所述步骤7的具体步骤为：

采用极大似然估计的方法，求得估计统计模型的参数的目标函数：

其中，

基于步骤5中求得的残差序列并以极大化为目标，采用线性规划的方法，求得模型参数进而估计出噪声图像n_x,y在第i级第j方向数子带上的NSCT系数矩阵的方差(σ^(i,j))²。

其中，所述步骤8的具体步骤为：

分别对所有高频子带系数矩阵分别执行步骤3、4、5、6、7的操作，分别得到噪声n_x,y分解在每一个高频子带上的分量的方差(σ^(i,j))²；

噪声n_x,y的方差σ²的计算公式为：γ^(i,j)为加权系数，γ^(i,j)的计算公式为：

其中，所述L的最佳取值为3或4。

其中，所述r的最佳取值为4。

本发明的有益效果在于：

本发明的基于广义自回归异方差模型的噪声方差估计方法采用非子采样轮廓波变换和广义自回归异方差模型，设计出一套完善的噪声方差估计方法，该方法可以提高噪声方差估计的准确性，且适用于各种噪声水平的退化图像，为去噪、复原、特征提取等后续图像处理提供支持。

附图说明

利用附图对本发明做进一步说明，但附图中的内容不构成对本发明的任何限制。

图1为本发明的基于广义自回归异方差模型的噪声方差估算方法的流程示意图。

图2为步骤4中将二维高频子带系数矩阵转换为一维序列数据的方法图示。

图3为图像NSCT高频子带系数矩阵残差序列统计直方图与GARCH模型的概率密度函数的拟合情况图示。

图4(a)是无噪声污染的清晰图像的示例。

图4(b)是噪声标准差为σ₁＝10的含噪图像的示例。

图4(c)是噪声标准差为σ₂＝20的含噪图像的示例。

图4(d)是噪声标准差为σ₃＝30的含噪图像的示例。

图4(e)是噪声标准差为σ₄＝40的含噪图像的示例。

具体实施方式

为了使发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。

具体实施方式

本发明的一种基于广义自回归异方差模型的噪声方差估算方法，如图1至图3所示，包括如下步骤：

步骤1，读取被噪声污染的含噪图像g_x,y。

读取被噪声污染的含噪图像g_x,y＝f_x,y+n_x,y，其中f_x,y表示清晰图像，n_x,y表示噪声图像，含噪图像g_x,y、清晰图像f_x,y、噪声图像n_x,y的大小均为N×M，n_x,y服从均值为0、方差为σ²的正态分布，即n_x,y～N(0,σ²)下标x,y分别表示图像的行坐标与列坐标，1≤x≤N,1≤y≤M。

步骤2，对被噪声污染的含噪图像g_x,y做非子采样轮廓波变换，即NSCT。

对被噪声污染图像g_x,y做L级非子采样轮廓波变换(本实施示例中L最佳取值为3或4)，其中在第i(1≤i≤L)级的方向数为2^i-1，得到大小为N×M的低频子带系数矩阵和各级各方向上的高频子带系数矩阵 j为在第i级上的方向数)，高频子带系数矩阵的个数为 因为g_x,y包含f_x,y和n_x,y两个部分，可将子带系数矩阵表示为 其中为f_x,y的NSCT的系数矩阵，为n_x,y的NSCT的系数矩阵。

步骤3，对步骤2中的每个高频子带系数矩阵做去均值滤波处理。

对高频子带系数矩阵去均值滤波处理的过程可用下面的公式表示：

其中，为系数矩阵去均值滤波处理后得到的矩阵。

步骤4，将去均值滤波处理后的二维高频子带系数矩阵转换为一维序列数据。

按照从左到右、从上到下的顺序依次回形读取去均值滤波后的二维系数矩阵 中的元素，得到一维序列数据序列数据既包含有噪声图像的信息，又包含有清晰图像的信息。图2表示了将二维矩阵转为一维序列数据的读取顺序，每一个方格代表一个像素，箭头表示将二维系数矩阵转为一维序列数据的读取顺序。

步骤5，对一维数据序列建立r阶自回归模型，并求得其残差序列

对一维序列数据建立r阶自回归模型，r为大于3的自然数，r的最佳值为4；

采用Gram-Schmidt正交法确定自回归模型的阶数r(本实施示例中取r＝4)，模型中为自回归系数，为残差项，采样最小二乘估计的方法求得自回归模型的自回归系数序列减去自回归项即得到均值为0的残差序列残差序列包括两个部分： 其中为清晰图像的残差信息，为噪声图像的残差信息。

步骤6，对残差序列建立统计学模型。

用均值为0、方差为(σ^(i,j))²的正态分布模型对噪声图像的残差信息建模，得到模型：用(1,1)阶广义自回归异方差模型GARCH(1,1)对清晰图像的残差信息建模，得到模型：

其中；α₀＞0,α₁≥0,β₁≥0，α₁+β₁＜1，因为噪声图像与清晰图像相互独立，可得残差序列服从的统计模型为：模型中包含的参数有

步骤7，基于步骤5中求得的残差序列和步骤6中建立的统计模型，采用极大似然估计的方法计算统计学模型的参数并估计出噪声图像n_x,y在第i级第j方向数子带上的NSCT系数矩阵的方差(σ^(i,j))²。

采用极大似然估计的方法，求得估计统计模型的参数的目标函数：

其中，式中E(·)表示求 数学期望。基于步骤5中求得的残差序列并以极大化为目标，采用线性规划的方法，可以求得模型参数从而估计出噪声图像n_x,y在第i级第j方向数子带上的NSCT系数矩阵的方差(σ^(i,j))²。

步骤8，求得含噪图像中的噪声的方差。

分别对所有高频子带系数矩阵分别执行步骤3、4、5、6、7的操作，分别得到噪声n_x,y分解在每一个高频子带上的分量的方差(σ^(i,j))²，(高频子带系数矩阵的具体个数为：个)。

待估计的噪声n_x,y的方差σ²的计算公式即为：γ^(i,j)为加权系数，γ^(i,j)的计算公式为： 在步骤2中有说明，为子带系数矩阵。

非子采样轮廓波变换是一种真正意义上的多尺度的图像二维表示方法，与小波变换相比，它能更细致的刻画和区分图像的高频细节信息。此外，清晰图像(无噪声)的非子采样轮廓波变换的子带系数矩阵数据经过去均值滤波(步骤3)和去自回归项处理(步骤5)后的残差序列数据具有明显的“尖峰厚尾”统计特性，而广义自回归异方差模型(GARCH)是一种能很好地刻画“尖峰厚尾”统计特性的模型，如图3所示，清晰图像残差序列统计直方图(黑色柱状为残差序列统计直方图)与GARCH模型的概率密度函数曲线(虚线为GARCH模型的概率密度函数曲线)拟合得很好，说明GARCH模型可以很好的刻画清晰图像残差序列的统计特性。本发明用GARCH来刻画清晰图像的残差信息用高斯分布刻画噪声图像的残差信息(步骤6)，进而得到含噪图像残差序列的统计模型求解该模型参数即估计出噪声图像n_x,y在第i级第j方向数上的NSCT系数矩阵的方差(σ^(i,j))²，最后综合各子带的噪声方差，即可求得含噪图像中的噪声的方差。

具体实验与结果分析：

如图4(a)、图4(b)、图4(c)、图4(d)、图4(e)所示，其中，图4(a)是无噪声污染的清晰图像示例，通过仿真的方法在图4(a)上分别添加均值为0，标准差(标准差的平方即为方差)取一系列值的高斯噪声，加噪后的图像见图4(b)、图4(c)、图4(d)、图4(e)，其中图4(b)中噪声标准差为σ₁＝10，图4(c)中噪声标准差为σ₁＝20，图4(d)中噪声标准差为σ₁＝30，图4(e)中噪声标准差为σ₁＝40。

为说明本发明方法的有效性，我们将本发明方法与其他两种在图像噪声方差估计领域被广泛认可的方法做比较，这两个方法分别是：Donoho的方法和Ghazal的方法。本发明方法与Donoho的方法、Ghazal的方法分别对图4(b)、图4(c)、图4(d)、图4(e)中所含噪声的方差的估计结果见表1。

表1：三种方法实验统计结果

表1中的实验统计结果显示：无论在噪声水平较低(比如σ＝10)时，还是噪声水平较高(比如σ＝40)时，本发明方法估计结果的误差比Donoho方法、Ghazal方法估计结果的误差都要小(实验结果中，Donoho方法仅在σ＝30时与本发明方法的估计误差相当，但总体误差比本发明方法误差大)。本发明方法估计的相对误差都控制在百分之几，且噪声污染越严重，本发明方法的性能越好，可以达到1％左右的相对误差。

综上所述，非子采样轮廓波变换能细致的刻画和区分图像的高频细节信息，广义自回归异方差模型(GARCH)能很好地刻画图像非子采样轮廓波变换高频系数矩阵“尖峰厚尾”的统计规律，本发明充分利用上述特性，提出的方法能够有效的估计出噪声方差，且与当前国际上被广泛认可的方法相比，本发明估计的结果具有更小的误差、更高的精度，并适用于各种噪声水平的图像，可为去噪、复原、特征提取等后续图像处理提供有效的支持。

最后应当说明的是，以上实施例仅用于说明本发明的技术方案而非对本发明保护范围的限制，尽管参照较佳实施例对本发明作了详细说明，本领域的普通技 术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的实质和范围。

Claims (10)

1.一种基于广义自回归异方差模型的噪声方差估计方法，其特征在于，包括以下步骤：

步骤1，读取被噪声污染的含噪图像；

步骤2，对含噪图像做非子采样轮廓波变换；

步骤3，对步骤2中的每个高频子带系数矩阵做去均值滤波处理；

步骤4，将去均值滤波处理后的二维高频子带系数矩阵转换为一维序列数据；

步骤5，对一维序列数据建立广义自回归异方差模型，并求得其残差序列；

步骤6，对残差序列建立统计学模型；

步骤7，根据步骤5中求得的残差序列和步骤6的统计学模型，采用极大似然估计的方法计算统计学模型的参数；

步骤8，根据步骤7中统计学模型的参数求得含噪图像中的噪声的方差。

2.根据权利要求1所述的一种基于广义自回归异方差模型的噪声方差估计方法，其特征在于，

所述步骤1的具体步骤为：

读取被噪声污染的含噪图像g_x,y＝f_x,y+n_x,y，

其中：f_x,y为清晰图像，n_x,y为噪声图像；含噪图像g_x,y、清晰图像f_x,y、噪声图像n_x,y的大小均为N×M；n_x,y服从均值为0、方差为σ²的正态分布；下标x,y分别表示图像的行坐标与列坐标，1≤x≤N,1≤y≤M；

所述步骤2的具体步骤为：

对含噪图像g_x,y做L级非子采样轮廓波变换，即NSCT，其中，在第i级的方向数为2^i-1，得到大小为N×M的低频子带系数矩阵和各级各方向上的高频子带系数矩阵高频子带系数矩阵的个数为1≤i≤L,1≤j≤2^i-1，j为在第i级上的方向数，L为大于2的自然数；

其中，为f_x,y的NSCT的系数矩阵，为n_x,y的NSCT的系数矩阵。

3.根据权利要求2所述的一种基于广义自回归异方差模型的噪声方差估计方法，其特征在于，所述步骤3的具体步骤为：

采用以下公式对每个高频子带系数矩阵去均值滤波处理，

其中，为系数矩阵去均值滤波处理后得到的矩阵。

4.根据权利要求3所述的一种基于广义自回归异方差模型的噪声方差估计方法，其特征在于，所述步骤4的具体步骤为：

按照从左到右、从上到下的顺序依次回形读取去均值滤波后的二维系数矩阵中的元素，得到一维序列数据序列数据既包含有噪声图像的信息，又包含有清晰图像的信息。

5.根据权利要求4所述的一种基于广义自回归异方差模型的噪声方差估计方法，其特征在于，所述步骤5的具体步骤为：

对一维序列数据建立r阶广义自回归异方差模型，r为大于3的自然数，

为均值为0的残差序列，

采用Gram-Schmidt正交法确定广义自回归异方差模型的阶数r，采用最小二乘估计的方法求得广义自回归异方差模型的自回归系数进而得到均值为0的残差序列

所述残差序列包括两部分：其中，为清晰图像的残差信息，为噪声图像的残差信息。

6.根据权利要求5所述的一种基于广义自回归异方差模型的噪声方差估计方法，其特征在于，所述步骤6的具体步骤为：

用均值为0、方差为(σ^(i,j))²的正态分布模型对噪声图像的残差信息建模，得到模型：

用(1,1)阶广义自回归异方差模型GARCH(1,1)对清晰图像的残差信息建模，得到模型：

其中：α₀＞0,α₁≥0,β₁≥0，α₁+β₁＜1；

从而，残差序列服从的统计学模型为：

所述统计学模型中包含的参数有

7.根据权利要求6所述的一种基于广义自回归异方差模型的噪声方差估计方法，其特征在于，所述步骤7的具体步骤为：

采用极大似然估计的方法，求得估计统计学模型的参数的目标函数：

其中，

基于步骤5中求得的残差序列并以极大化为目标，采用线性规划的方法，求得统计学模型的参数进而估计出噪声图像n_x,y在第i级第j方向数子带上的NSCT系数矩阵的方差(σ^(i,j))²。

8.根据权利要求7所述的一种基于广义自回归异方差模型的噪声方差估计方法，其特征在于，所述步骤8的具体步骤为：

分别对所有高频子带系数矩阵分别执行步骤3、4、5、6、7的操作，分别得到噪声n_x,y分解在每一个高频子带上的分量的方差(σ^(i,j))²；

噪声n_x,y的方差σ²的计算公式为：γ^(i,j)为加权系数，γ^(i,j)的计算公式为：

9.根据权利要求2所述的一种基于广义自回归异方差模型的噪声方差估计方法，其特征在于，所述L的最佳取值为3或4。

10.根据权利要求8所述的一种基于广义自回归异方差模型的噪声方差估计方法，其特征在于，所述r的最佳取值为4。