混合二阶全变分的抗核辐射图像降噪方法

黄小莉，陈春梅,2，刘桂华,2

(1.西南科技大学 信息工程学院，四川 绵阳 621000；2.特殊环境机器人技术四川省重点实验室，四川 绵阳 621000)

0 引 言

核环境中，图像采集器采集到的图像含有大量核噪声亮斑，加大了后续对核辐射场景感知的难度。目前，对于核辐射干扰图像的降噪研究相对较少，对核环境下的图像进行降噪处理具有极高的学术研究价值和较强的创新性。

传统的图像降噪方法，如中值滤波法[1]、均值滤波法[2]等在进行降噪时会对图像的细节信息造成损失。文献[3]首次提出全变分(total variation，TV)算法，很好地解决了传统降噪方法造成的图像边缘模糊的问题。文献[4]将分裂布雷格曼方法应用于全变分模型中，解决了全变分模型的正则化问题。文献[5-6]提出了基于邻近算子的全变分方法(proximity approximation operator total variation，PATV)，通过线性变换计算凸函数组成的近似算子提升降噪性能。

为了缓解全变分降噪引起的阶梯伪影问题，很多文献进行了大量研究。文献[7]提出了广义全变分(generalized total variation，TGV)，在一定程度上能抑制阶梯伪影的产生。文献[8]提出了一种自适应TGV模型的降噪算法，通过自适应参数使得边缘和区域平滑。文献[9]提出了一种重叠组稀疏全变分(overlapping group sparse total variation，OGSTV)方法，通过突出图像平滑区域与边缘区域的差异性，有效抑制了阶梯伪影。文献[10]将高阶非凸正则化加入到OGSTV上，使图像纹理更为平滑。文献[11-12]提出的高阶方程法在抑制阶梯伪影的同时很好地保留了图像边缘信息。

为了解决高阶全变分造成的图像过度平滑问题，很多学者提出了混合技术。文献[13]将一阶全变分同二阶全变分相结合，既保持了图像的边缘又抑制了阶梯效应。文献[14]提出利用交替方向乘子法(alternating direction method of multiplier, ADMM)作为快速求解器有效地解决L1-TV正则化问题。文献[15]提出一种针对脉冲噪声去噪的带L1保真项的混合变分模型，增强模型的去噪性能。

受到混合算法的启发，本文提出一种混合去噪算法，将重叠组稀疏正则化与非凸二阶全变分方法相结合，采用ADMM和增广拉格朗日乘子法进行交替循环求解，将降噪后的图像进行相减，在初始图像上对得到的差值进行叠加操作与迭代处理，优化最终降噪效果。与对比算法的结果相比，本算法保留了更多的边缘细节信息，能更大程度地还原图像。

1 混合二阶全变分

图像降噪可以看成一个线性逆问题的求解过程，即从待处理图像中推导出原图像，常见的去噪模型为

I=Kg+γ

(1)

(1)式中：g∈RMN×1为未受噪声污染的原始图像；γ表示噪声；K为点扩散函数；I∈RMN×1为含有噪声的待处理图像。

求解图像降噪的关键在于对(1)式进行正则化处理，即将原始图像g中的先验知识融入线性问题的求解过程中，由此获得具有正则化解，同时抑制噪声。

1.1 二维重叠组稀疏正则化

传统的组稀疏正则化器[16]方程表示为

(2)

图像g的二维M×M点群表示为

(3)

(4)

1.2 重叠组稀疏正则化与二阶全变分混合

将非凸二阶全变分与重叠组稀疏正则化混合构成本文提出方法的前半部分，可以简单表示为

(5)

通过假设图像周期性边界条件，已知∇g与∇2g为一阶和二阶差分算子，引入辅助变量将(5)式转换为约束优化问题，表示为

(6)

1.3 ADMM和增广拉格朗日乘子法优化求解

ADMM是一种求解优化问题的计算方法[19]，将受约束的全局问题分解为多个局部子问题。在本算法中，ADMM通过引入辅助变量将全局问题分裂为多个子问题，交替迭代优化多个子问题的解从而得到全局问题的解。使用增广拉格朗日函数和ADMM对(6)式进行求解，增广拉格朗日函数表示为

(7)

(8)

(8)式中，求解gm是一个最小二乘问题。

在图像的周期边界条件下，由于矩阵均为块循环，因此，使用二维离散傅里叶变换对(8)式进行求解，通过将时域乘法问题转换为频域问题，极大地减少计算复杂度，表示为

(9)

(9)式中，F-1为傅里叶逆变换。T的计算表达式为

T=

(10)

1.3.1 变量ζ1的优化

变量ζ1的优化是一个重叠组稀疏性问题。可以表示为

(11)

使用优化-最小化(MM)[20]方法解决重叠组稀疏的问题，在每次迭代中使用一个二次函数φ(ζ1)。在MM算法的每次迭代中，二阶问题被最小化，(11)式的收敛式为

(12)

(12)式中，重叠组稀疏全变分正则化φ(ζ1)的定义为

(13)

(13)式中，M为重叠组大小。

1.3.2 变量ζ2的优化

变量ζ2的优化是一个非凸二阶去噪问题，可以表示为

(14)

采用分离变量法将(14)式改写为

(15)

采用迭代重加权算法(IRL1)[21]对(15)式进行求解，每一次迭代，(15)式都可以表示为

(16)

(16)式中，每次迭代更新的权重t表达式为

(17)

(17)式中，δ是一个避免被除数为0的小数字。

利用收缩阈值函数对(16)式进行计算，可得

(18)

(18)式中，shrink(·)表示一维收缩运算符。

1.3.3 变量ζ3的优化

变量ζ3是一个确保像素值保持在[0,255]的投影，表示为

(19)

采用投影的方法来求解(19)式，使用更新规则，将迭代次数大大减小，最终的拉格朗日乘子迭代为

(20)

1.4 差值迭代

将上述混合非凸二阶全变分算法降噪过程设为H，降噪后的图像与原图像进行迭代降噪，多次操作后得到图像的灰度值接近于原始图像的灰度值，提高了降噪的效率。具体流程如下。

1)将图像g通过上述提出的混合非凸二阶全变分算法进行降噪，得到处理后的图像g1为

g1=H(g,1)

(21)

2)将处理后的图像g1再次进行混合非凸二阶全变分降噪得到g2为

g2=H(g1,1)

(22)

将g1与g2相减得到gΔ，表示为

gΔ=g1-g2=H(g,1)-H(g1,1)

(23)

设定β为迭代次数的倒数，修改β的值来提高算法的降噪鲁棒性。令β∈(0,1)，将差值gΔ进行β倍的缩小并与原图像相加，得到图像g3，g3不仅包含有初始图像的框架，还包含有前两次降噪后的结果，得到的g3为

g3=g+β(g1-g2)=g+βgΔ

(24)

3)对g3进行连续迭代降噪，设迭代次数为Ter，得到最终图像gn为

gn=H(g3,Ter)

(25)

2 模拟核噪声

2.1 核噪声的特性

在核辐射场景下，将图像采集镜头用镜盖盖住，此时无可见光输入，核辐射粒子直接作用于传感器，生成噪声亮斑图像如图1所示。从图1a可以看出，核噪声为多个块状的亮斑，每个亮斑的像素值相近。增强核辐射强度，得到图1b所示图像，对比图1a、图1b能看出，当核辐射强度增强时，放射性粒子产生的脉冲电流的数量与强度也随之增加，同时图像上形成的噪声亮斑的数量和面积也会随之增加。由图1c和图1d得出，灰度值的幅值服从高斯分布且灰度级在50～150出现的概率最大。

图1 无可见光下不同辐射强度采集的图像及对应灰度直方图

2.2 核噪声的模拟

根据核噪声的随机性，采用随机行走方法来模拟核噪声斑块的形状。将图像中的像素点比作核辐射粒子，随机选择一个像素点作为随机行走的起点，以起点为中心向周围8领域以概率P进行扩散。通过调整行走概率P和迭代次数来调整模拟的噪声斑块面积，噪声斑块的数量则由斑块的分布率决定。将随机行走得到的矩阵与高斯矩阵相乘得到模拟核噪声斑块图。

综上分析，模拟核噪声斑块的具体流程如下。

1)生成一个只含有0和1的二值随机矩阵，矩阵中1的个数为泊松函数得到的泊松随机数，将生成的泊松随机数作为模拟生成的核噪声个数；

2)以矩阵中的1为起点向周围8领域以概率P进行随机行走，通过调整概率P和迭代次数来调整模拟的噪声面积，接下来设置斑块的分布率-斑块噪声的像素个数与图像中所有像素个数的比值，从而模拟得到不同的斑块数量与不同辐射强度中的噪声斑块；

3)将随机行走得到的矩阵与高斯矩阵相乘得到最终模拟的核噪声斑块；

4)将模拟得到的核噪声斑块加入到不同场景的自然图像中，得到不同场景下的核辐射噪声图像。

图2所示为不同参数模拟的不同强度核噪声斑块，设定图2a中的模拟弱核斑块行走概率P为0.1，迭代次数为5，斑块分布率为0.3；图2b中的模拟中度核斑块行走概率P为0.1，迭代次数为10，斑块分布率为0.5；图3c中的模拟辐射较强时的核斑块行走概率P为0.1，迭代次数为25，斑块分布率为0.8。由核噪声特性可知，随着核辐射强度的增加，产生的核辐射斑块噪声的面积变大，斑块的密度也随之增加。

图2 不同参数模拟的不同强度核噪声斑块

将模拟的核斑块加入到数字图像处理常用的数据集set12中，图3为其中一张图像的效果。本文通过模拟得到不同场景下的核辐射噪声图片，为后续实验提供条件。

图3 加入不同强度模拟噪声后的图像

3 实验参数的选取

3.1 差值迭代次数的确定

在实验过程中需要设置差值迭代次数Ter，迭代次数过大会消耗过多的计算时间，迭代次数过小则无法达到最佳降噪效果。用不同迭代次数对获得的模拟数据集进行实验，以峰值信噪比(peak signal-to-noise ratio, PSNR)和结构相似性(structural similarity, SSIM)为指标来获得最佳迭代次数。图4为不同迭代次数处理后的PSNR和SSIM折线图。通过图4看出，差值迭代次数Ter为15时，PSNR和SSIM数值达到最高，继续增加迭代次数，PSNR和SSIM无明显变化。

图4 不同迭代次数处理后的PSNR和SSIM折线图

3.2 重叠组M大小的确定

在实验过程中需要设置重叠组M的大小。固定其他参数，选取模拟辐射强度为弱的数据集进行实验，以PSNR和SSIM为指标来获得最佳重叠组M的大小。图5为不同重叠组M处理后的PSNR和SSIM折线图。从图5可知，当M=4时，PSNR和SSIM数值最大。由此，设置重叠组M的大小为4。

图5 不同重叠组M处理后的PSNR和SSIM折线图

3.3 其他参数的调整

首先，对变量λ的范围进行调整。变量λ控制数据保真项，在图像去噪上具有重要作用。由于在实验中使用不同模拟噪声强度的数据集，因此，根据经验设置λ∈[1×10-4，4]。

然后，对混合算法的内部迭代次数进行调整。设置差值迭代次数分别为1、5、10、15、20，对不同情况下的PSNR进行对比，以此得到最佳的内部迭代次数。选取模拟辐射强度为弱的数据集“house”和“woman”进行实验，结果如表1所示。当迭代次数为10的时候，PSNR数值达到最大值。由于迭代次数过高会导致计算时间增加，综合考虑，选取内部迭代最佳次数为10。

表1 不同迭代次数的降噪结果

4 实验结果

为了验证算法的有效性和可靠性，在两类不同数据集上进行实验，数据集分为模拟得到的核噪声数据集及真实核环境下的数据集。本算法与传统中值滤波、邻近算子全变分方法(PATV)[5]、广义全变分方法(TGV)[8]、重叠组稀疏的全变分方法(OGSTV)[9]、L1全变分方法[14]等多个算法进行比较，采用峰值信噪比、结构相似性以及图像粗糙度作为客观评判算法性能的指标，同时从主观角度出发，通过图像细节部分对图像最终处理结果进行评定。本实验在处理器为AMD Ryzen 72700X Eight-Core Processor at 3.70 GHz，内存为16.0 GByte，系统为Windows 10专业版，MATLAB为R2015b的环境下进行实验。

4.1 对模拟得到的图像进行实验

对获得的多张不同辐射强度的图片进行处理，表2为不同去噪算法的结果对比。由表2可见，在不同的模拟核辐射强度中，本文算法得到的PSNR和SSIM值最高。将不同算法得到的最终结果的红框处进行放大，如图6—图7所示。由图6—图7能看出，本算法在降噪的同时能更好地保留图像的边缘信息。

表2 不同去噪算法的结果对比

图6 “peppers”处理后的局部放大图

图7 “woman”处理后的局部放大图

4.2 对真实核环境下的核图像进行实验

为了进一步验证算法的有效性，对真实核环境下采集到的噪声图像进行降噪处理，同时将处理结果与不同算法进行对比。

图8为不同算法对核环境下噪声图片进行处理后的结果。图8a为广州番禺事故及河南杞县事故处理时在真实核环境下采集到的图像，整个图像被噪声包围且能见度极低，严重影响了核环境的识别及后续操作。将图8b—图8g红框处部分进行放大，可得到“铁柱”处理后的局部放大图如图9所示，得到“小车”处理后的局部放大图如图10所示。从图9可知，本算法结果的细节处更加清晰，引起的图像缺陷也更少，裂缝衔接处纹理更加流畅。从图10可知，本算法比其他算法的处理结果损失更少，结果中含有的噪声含量较少，同时物体边缘更加清晰。

图8 不同算法对核环境下噪声图片进行处理的结果

图9 “铁柱”处理后的局部放大图

图10 “小车”处理后的局部放大图

5 结束语

本文通过将非凸二阶全变分与重叠组稀疏结合起来，在降噪的同时抑制阶梯效应，从而极大地保留了图像的细节信息。在此混合算法的基础上，将处理后得到的图像进行差值迭代，使得图像的灰度更加偏向原图灰度，进一步提升了图像的质量。实验结果表明，本混合算法在视觉主观效果上和客观数据评价指标上优于对比算法，证实了算法的优越性和可行性。