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

摘要

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

说明书

技术领域

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

背景技术

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

现有的图像噪声方差估计的方法主要有三大类：一类是基于滤波的方法，这 类方法先通过滤波来提取图像的高频信息，再从高频信息中来估计噪声的方差， 然而，该类方法容易将图像的边缘和纹理等高频细节信息视为噪声，导致噪声方 差的过估计。第二类是基于分块的方法，该方法依赖“同质块的局部方差可以作 为该块的噪声方差的估计量”这一假设，先按一定准则把噪声图像分成很多小块， 然后从诸多小块中寻找同质子块，再基于同质子块来估计图像噪声的方差，该类 方法性能的好坏取决于如何科学有效的确定图像的同质块，并保证有足够多的同 质块参与到噪声方差的估计中。第三类是基于小波域统计建模的方法，该类方法 中最著名的是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(1≤i≤L) 级的方向数为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的具体步骤为：

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

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

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

按照从左到右、从上到下的顺序依次回形读取去均值滤波后的二维系数矩阵 中的元素，得到一维序列数据$\left\{s_t^{(i,j)}\right\}=\{s_1^{(i,j)},s_2^{(i,j)},s_3^{(i,j)}......s_t^{(i,j)}......s_{M\times N}^{(i,j)}\},$序列 数据既包含有噪声图像的信息，又包含有清晰图像的信息。

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

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

$s_t^{(i,j)}=b_0^{(i,j)}+b_1^{(i,j)}{s}_{t-1}^{(i,j)}+b_2^{(i,j)}{s}_{t-2}^{(i,j)}+......b_r^{(i,j)}{s}_{t-r}^{(i,j)}+e_t^{(i,j)},$为均值为0的残差序列，

采用Gram-Schmidt正交法确定自回归模型的阶数r，采用最小二乘估计的方 法求得自回归模型的自回归系数进而得到均值为0的残 差序列$e_t^{(i,j)}=s_t^{(i,j)}-b_0^{(i,j)}-b_1^{(i,j)}{s}_{t-1}^{(i,j)}-b_2^{(i,j)}{s}_{t-2}^{(i,j)}-b_3^{(i,j)}{s}_{t-3}^{(i,j)}-b_4^{(i,j)}{s}_{t-4}^{(i,j)};$

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

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

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

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

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

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

所述统计学模型中包含的参数有 $\{{\alpha }_0^{(i,j)},{\alpha }_1^{(i,j)},{\beta }_1^{(i,j)},{\sigma }^{(i,j)}\}.$

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

采用极大似然估计的方法，求得估计统计模型的参数 $\Gamma =\{{\alpha }_0^{(i,j)},{\alpha }_1^{(i,j)},{\beta }_1^{(i,j)},{\sigma }^{(i,j)}\}$的目标函数：

$\log \left(f\left(e_t^{(i,j)}\right)\right)=-0.5\underset{t=2}{\overset{M\times N}{\Sigma }}\log ({\left({\sigma }^{(i,j)}\right)}^2+h_{t|t-1}^{(i,j)})-\underset{t=2}{\overset{M\times N}{\Sigma }}({\left(e_t^{(i,j)}\right)}^2/2({\left({\sigma }^{(i,j)}\right)}^2+h_{t|t-1}^{(i,j)}))$

其中，$h_{t|t-1}^{(i,j)}={\alpha }_0^{(i,j)}+{\alpha }_1^{(i,j)}E\left({}_{f}e_{t-1}{}_{(i,j)}\right|\{e_{t-1}^{(i,j)},h_{t-1|t-2}^{(i,j)}\})+{\beta }_1^{(i,j)}{h}_{t-1|t-2}^{(i,j)},$

基于步骤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的低频 子带系数矩阵和各级各方向上的高频子带系数矩阵(1≤i≤L,1≤j≤2^i-1，j为在第i级上的方向数)，高频子带系数矩阵的个数为 因为g_x,y包含f_x,y和n_x,y两个部分，可将子带系数矩阵表示为 其中为f_x,y的NSCT的系数矩阵，为n_x,y的NSCT 的系数矩阵。

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

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

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

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

按照从左到右、从上到下的顺序依次回形读取去均值滤波后的二维系数矩阵 中的元素，得到一维序列数据$\left\{s_t^{(i,j)}\right\}=\{s_1^{(i,j)},s_2^{(i,j)},s_3^{(i,j)}......s_t^{(i,j)}......s_{M\times N}^{(i,j)}\},$序列 数据既包含有噪声图像的信息，又包含有清晰图像的信息。图2表示了将 二维矩阵转为一维序列数据的读取顺序，每一个方格代表一个像素，箭头表 示将二维系数矩阵转为一维序列数据的读取顺序。

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

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

$s_t^{(i,j)}=b_0^{(i,j)}+b_1^{(i,j)}{s}_{t-1}^{(i,j)}+b_2^{(i,j)}{s}_{t-2}^{(i,j)}+......b_r^{(i,j)}{s}_{t-r}^{(i,j)}+e_t^{(i,j)},$采用Gram-Schmidt正交法确定 自回归模型的阶数r(本实施示例中取r＝4)，模型中为自 回归系数，为残差项，采样最小二乘估计的方法求得自回归模型的自回归系 数序列减去自回归项即得到均值为0的残差序列 $e_t^{(i,j)}=s_t^{(i,j)}-b_0^{(i,j)}-b_1^{(i,j)}{s}_{t-1}^{(i,j)}-b_2^{(i,j)}{s}_{t-2}^{(i,j)}-b_3^{(i,j)}{s}_{t-3}^{(i,j)}-b_4^{(i,j)}{s}_{t-4}^{(i,j)};$残差序列包括两个部分： 其中为清晰图像的残差信息，为噪声图像的残差信息。

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

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

其中；α₀＞0,α₁≥0,β₁≥0，α₁+β₁＜1，因为噪声图像与清晰图像相互独立，可 得残差序列$e_t^{(i,j)}={}_{f}e_t{}_{(i,j)}+{}_{n}e_t{}_{(i,j)}$服从的统计模型为：模型中 包含的参数有

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

采用极大似然估计的方法，求得估计统计模型的参数 $\Gamma =\{{\alpha }_0^{(i,j)},{\alpha }_1^{(i,j)},{\beta }_1^{(i,j)},{\sigma }^{(i,j)}\}$的目标函数：

$\log \left(f\left(e_t^{(i,j)}\right)\right)=-0.5\underset{t=2}{\overset{M\times N}{\Sigma }}\log ({\left({\sigma }^{(i,j)}\right)}^2+h_{t|t-1}^{(i,j)})-\underset{t=2}{\overset{M\times N}{\Sigma }}({\left(e_t^{(i,j)}\right)}^2/2({\left({\sigma }^{(i,j)}\right)}^2+h_{t|t-1}^{(i,j)}))$

其中，$h_{t|t-1}^{(i,j)}={\alpha }_0^{(i,j)}+{\alpha }_1^{(i,j)}E\left({}_{f}e_{t-1}{}_{(i,j)}\right|\{e_{t-1}^{(i,j)},h_{t-1|t-2}^{(i,j)}\})+{\beta }_1^{(i,j)}{h}_{t-1|t-2}^{(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)能很好地刻画图像非子采样轮廓波变换高频 系数矩阵“尖峰厚尾”的统计规律，本发明充分利用上述特性，提出的方法能够 有效的估计出噪声方差，且与当前国际上被广泛认可的方法相比，本发明估计的 结果具有更小的误差、更高的精度，并适用于各种噪声水平的图像，可为去噪、 复原、特征提取等后续图像处理提供有效的支持。

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