（浙江理工大学 理学院，浙江 杭州 310018）

0 引言

周期性噪声和准周期性噪声可能源于图像采集与传输过程中的机电干扰，也可能来自电子干扰，例如现代数码相机的电子电路和电荷耦合器件等。当多成像系统中检测器的灵敏度不匹配时，经常会产生一种特殊的条纹状拟周期性噪声，这种噪声被称为摩尔纹（moiré）噪声。摩尔纹噪声经常出现在磁共振成像、显微镜图像及遥感图像等实时图像的应用过程中［1］。由于摩尔纹噪声在空间上是相关的，并散布在整个图像中，因此使用空间域滤波算法去除摩尔纹噪声可能会在实体周围产生伪影，有时还会使恢复的图像变得模糊。但在傅里叶变换图像频谱中，可以很容易观测到由摩尔纹噪声引起的高振幅频谱分量。因此，通过频域操作恢复这些频谱分量能够抑制噪声带来的影响，从而恢复原始图像信息。

为去除摩尔纹噪声，可采用的一种方法是使用频域滤波器。如Moallem 等［2］提出一种新颖的自适应高斯陷波滤波器，可有效降低周期与准周期噪声；Varghese［3］使用交换中值滤波器（Switching Median Filter）去除高密度的周期噪声。在此基础上，Charkraborty 等［4-5］利用双阈值法确认波峰位置，并使用自适应高斯恢复滤波器（Adaptive Gaussian Restoration Filter）去除周期噪声，并利用自动陷波抑制滤波器对Gabor 变换后图像的特定周期信号进行滤波，从而去除对应周期区间内的准周期噪声。在医学领域，Ionita 等［6］采用自适应阈值，能更好地确定受周期噪声影响的傅里叶频域波峰检测阈值，然后使用窗口高斯陷波滤波器去除出现在显微镜图像中的周期噪声。此外，Ionita 等［7］通过分析一阶Haar 小波变换频域，确定了与噪声对应的频率区域，然后通过校正显微镜图像的傅里叶幅度谱分量去除噪声。另一种方法是使用最小化模型去除摩尔纹噪声。如Chang等［8］提出一种针对遥感影像的图像恢复算法，该算法首先对条纹信号进行分类，然后利用全变差（Total variation）模型去除噪声信号，以保持影像光谱与空间的一致性。由于该模型需要目标图像的特征梯度小于条纹梯度，因此对于一般的自然图像，使用该算法可能会破坏原始图像边缘信息。

上述研究虽然能较为有效地去除摩尔纹噪声，但恢复的傅里叶频域信号依赖于其领域信号。因此，若图像在傅里叶低频区域内受噪声影响的系数范围较大时，上述研究无法准确地恢复原始图像低频信息。为更好地恢复原始图像低频信息，本文提出一种基于紧小波框架的图像恢复算法。算法本质是利用紧框架下图像小波系数的稀疏性特征，根据压缩感知原理建立优化模型，通过求解该模型最优解得到原始图像稀疏恢复结果。该算法有助于在机电干扰较为严重的情况下进行图像采集。

1 图像恢复原理

1.1 稀疏信号恢复

压缩感知理论表明，利用目标（如信号、图像等）的稀疏性，可使用远远少于香农采样定理（Shannon-Nyquist Sampling Theorem）所需的样本数重构原始信号。以带高斯噪声的信号恢复问题为例，若n维列向量x为需要恢复的原始稀疏信号，则实际信号采集过程可表示为：

其中，y∈ℝm为采集到的信号，A∈ℝm×n为测量矩阵，w∈ℝm为独立同分布的高斯噪声项。因此，稀疏信号恢复问题可转换为一个线性方程组的求解问题。在实际问题中，矩阵A的条件数通常较大，因此无法直接采用最小二乘法。一种可行的方法是借助原始信号x的稀疏性使得问题具有唯一解。

目前普遍采用l0范数和l1范数衡量目标稀疏度。对于n维向量x=(x1,x2,…,xn)，l0范数为x中的非零元素个数，用‖x‖0表示。虽然‖x‖0不满足范数的定义，但为了统一起见，仍称其为范数。l1范数为x中所有元素绝对值之和，用‖x‖1表示。如果问题（1）中的噪声项w满足‖w‖2＜ε（ε为一个近似于零的正阈值），且信号x∈ℝn×1是稀疏的，即x的l0范数远小于信号长度，则问题（1）可转化为带约束的非凸优化模型（模型1）：

求解l0范数的过程是NP（Non-deterministic Polynomi⁃al-time）难问题。由于l1范数可以很好地逼近l0范数，因此可考虑用l1项替换模型1 中的l0项，从而得到一个凸的优化模型（模型2）：

Donoho 等［8-9］证明了当矩阵A满足一定的零空间性质时，模型1 与模型2 的解是等价的；Candès 等［10-11］证明了当矩阵A满足一定的约束等距性质（Restricted Isometry Property，RIP）时，模型1 与模型2 的解是稳定且等价的。此外，Candès 等［12］还证明了若矩阵A∈ℝm×n为随机高斯矩阵，且A的行列满足m=O()kln(n/k)，矩阵A满足k阶约束等距性质的概率不低于。在此基础上，Rul⁃son 等［13-14］证明了当随机高斯矩测量次数k满足一定条件时，通过模型2 能够精确地恢复原始信号。由于l1项是凸的，因此可利用凸优化算法快速求解模型2。为便于求解，大部分算法先将模型2 转化为无约束形式（模型3）：

其中，λ＞0 是一个正则化参数。模型2 到模型3 的转换在数学上可表述为：对于任意的ε＞0，存在一个λ使得模型2 与模型3 的最优解相同。

1.2 基于框架的图像恢复原理

数字图像的退化主要包括像素丢失、图像变模糊及像素点含有噪声等。一张灰度图像在数学上通常用向量x=(x1,x2,…,xn) 表示，其中xi∈[0,255] 为对应像素点的值。进入21 世纪之后，基于框架（特别是紧小波框架）的图像恢复方法得到了迅速发展。对一个给定的矩阵D∈ℝn×m(n≤m)，如果D满足DDT=I，其中I是单位矩阵，DT是D的转置，则D是一个紧框架。因此，对任何向量x∈ℝn，有：

此外，构成向量u=DTx的元素是用来表示x的规范系数。如果紧框架D的元素是通过有限多个小波函数的扩展与转换生成的，则D被称为一个紧小波框架。

利用压缩感知理论的一个重要前提是测量目标是稀疏的。自然图像虽无法满足稀疏条件，但自然图像的小波系数通常是稀疏或者可压的。小波分析采用细分函数构造小波函数，并通过小波函数的伸缩与平移构成函数空间的一组正交基或框架，空间中任一函数均可由小波函数生成的基或框架展开。在数字图像处理过程中，大多使用高通和低通滤波器作为图像分解与重构的基本工具，小波系数可呈现出图像在不同频域中的特征。现实中的大规模数据往往具有某种意义上的低维结构或特征。对于一般的自然图像，小波系数具有稀疏性特征，例如200 万像素的数字图像仅需10 万个小波系数便能获取几乎所有可见信息，剩余的190 万小波系数仅贡献少量、可忽略的细节信息。在数学上，小波系数的稀疏性可通过‖DTx‖p0≤p≤1 进行刻画。关于多尺度分析、细分函数、小波函数、小波滤波器等严格的数学理论和应用可参考文献［15］、［16］等。在稀疏性假设下，由文献［17］可知，图像恢复问题可建模为模型4：

同理，模型4 可转化为无约束形式（模型5）：

模型5 也被称为ALASSO 模型［18］。Lin 等［18-19］证明了如果矩阵A满足一定的D-RIP 条件，则模型5 的最优解是稳定的。

一般情况下，自然图像在空间域上具有一定连续性，因此自然图像所产生的小波系数一般较大，而高斯噪声经过小波变换后的系数具有较强的随机性，且噪声对应的大部分小波系数都很小。若设噪声的小波系数方差为σ，则根据切比雪夫不等式，绝大部分噪声系数都位于[-3σ,3σ]区间内。因此，信号经过小波变换后，使用阈值法可以较好地去除噪声信号。

1.3 摩尔纹噪声

在数码影像中，如果主体中有密纹的纹理，通常会出现类似水波的条纹，该条纹即为摩尔纹。图1 表示含有摩尔纹噪声的小丑图像去噪结果。

Fig.1 Denoising result of image“Clown”图1 小丑图像去噪结果

其中，图1（a）为含有摩尔纹噪声的小丑图像，其傅里叶频谱如图1（b）所示。图1（b）中正中心的高亮部分表示直流通分量，图像中心附近的白色十字状光谱特征即为摩尔纹噪声在傅立叶频谱中的尖峰状分量。因此，通过校正这些由噪声改变的频谱分量，可以有效抑制噪声的影响。图1（c）表示使用软阈值模型的处理结果，图1（d）表示提取的摩尔纹噪声。具体而言，对于一幅尺寸为M×N的光滑自然图像f，其傅立叶频谱为：

其中，H(k,l) 为傅立叶频谱H的第(k,l)个系数。傅立叶频谱中波峰位置可通过人工标注或使用阈值法加以确认。对于传统阈值法，判定频谱中第(k,l)点为波峰的条件为：该点系数Hkl与其周围I×J个点系数的中值MEDI×J(Hkl)满足不等式：

其中，CL×F为大小为L×F的低频置信区域，T2为给定阈值。在确认波峰所在位置之后，选用不同类型滤波器校正由噪声改变的频谱分量过程可表示为：

图2 表示高斯陷波滤波器的频域处理过程。

Fig.2 Gaussian notch filtering illustration图2 高斯陷波滤波器频域处理过程

对于周期性噪声，可选用高斯陷波滤波器［20］进行滤波，即滤波器GI×J满足：

2 基于紧小波框架的最小化模型

对于包含摩尔纹噪声的图像f，通过阈值法可确认其由噪声引起频谱中的波峰位置，然后根据这些波峰位置可建立一个掩膜矩阵（Mask matrix）。在数学上，掩膜矩阵建立过程可表示为：

其中，Pi,j为掩膜矩阵P的第(i,j)个元素；s为给定阈值，用于规范矩阵大小。掩膜矩阵的建立确定了需要恢复的频域系数范围。由于在紧小波框架下，这部分系数是稀疏的，根据模型5，可得到优化模型（模型6）：

其中，ℱ 表示傅立叶卷积过程，P为利用公式（13）求得的掩膜矩阵，y=Pℱf，λ是正则化参数，D是一个紧小波框架，x为测量到的图像，x为原始图像。为便于表述，本文称Θ=1 时的模型6 为软阈值模型，Θ=0 时的模型6 为硬阈值模型。

求解模型6 的难点在于其中的正则项‖Dx‖Θ是不可分离的。DFISTA（Decomposition-basedFast Iterative Shrink⁃age Thresholding Algorithm，DFISTA）算法［21］提供了一个可行的分解思路：通过引入辅助变量，从而以简单和显式的方式应用迭代收缩阈值算法求解模型。具体而言，对于最小化问题，引入辅助变量τ，得到损失函数（lose function），如模型7 所示：

由于l1范数是不光滑的，在τ为零处不可微，因此需要使用近端梯度法（Proximal Gradient Method，PG）。PG 算法思想是使用近端算子作为近似梯度，从而进行梯度下降。近端梯度下降法及其加速算法的收敛性分析可参见文献［21］。

则函数∇J的利普希茨常数满足不等式：

其中，ε是一个近似于零的正阈值。利用DFISTA 算法求解模型7 的过程如下：

由以上过程可知，影响恢复结果的主要因素包括：①检测波峰的阈值T2；②规范矩阵大小的阈值s；③紧小波框架函数选取；④小波阈值函数选取与阈值T1。

3 实验结果与分析

3.1 实验说明

为了验证本文算法的可行性与有效性，利用Matlab R2019a 软件进行仿真实验。实验对象选取Barbara、Boat、Cameramen 等4 幅灰度图像以及2 幅使用智能手机拍摄的彩色自然图像。对于纯净的灰度图像，首先在对应傅里叶频域中靠近中心的低频信号区域添加干扰噪声，然后对使用加速算法与未使用加速算法求解软阈值模型6 的PSNR（Peak Signal to Noise Ratio）值变化情况，以及分别使用中值滤波器［22］、窗口高斯陷波滤波器、由掩模矩阵直接生成的滤波器与模型7 恢复结果的PNSR 值进行对比。由于在傅里叶频域中，靠近中心的低频信号区域决定了大部分图像细节信息，因此添加的摩尔纹噪声主要集中在低频区域。对于彩色自然图像，将其红色、绿色与蓝色通道图像分别作为灰度图进行处理，最后再将去噪结果重新合成彩色图像［23］。

图3 表示对灰度图像去除摩尔纹噪声的实验流程。首先使用快速离散傅里叶逆变换函数得到与输入图像尺寸相同的摩尔纹噪声图像，然后向图像添加摩尔纹噪声并使用式（8）得到标记的傅里叶频域波峰位置，最后分别利用频域滤波器与模型6 对含有摩尔纹噪声的图像进行处理，并对比实验结果的峰值信噪比。其中，频域滤波器包括高斯陷波滤波器、中值滤波器以及直接使用掩膜矩阵生成的滤波器。

Fig.3 Experimental procedures图3 实验流程

为便于比较实验结果，设式（8）中，L=F=8，I=J=11，T2=5；式（10）中，α=1，β=0.1。根据波峰检测结果，设s为大于1 且小于6 的整数。模型7 中，设当Θ=1 时，μ=0.7；当Θ=0 时，μ=1。使用分段样条紧小波框架，设小波函数分解层次为4 层。

为能客观地比较算法结果，采用以下定义的PSNR 值对图像恢复质量进行评价：

其中，x表示恢复的图像，表示原始图像。

3.2 实验结果

为验证加速近端梯度算法的有效性与稳定性，在使用含有噪声的Cameramen 图像作为实验对象时，对比使用加速算法与未使用加速算法的收敛速度。即在利用DFISTA算法求解模型7 的过程中，设为加速算法，tk+1=tk为未加速算法。使用加速算法求解软阈值模型6 的迭代结果如图4 所示。

Fig.4 Iterative results of noised image“Cameramen”图4 含噪声图像Cameramen 迭代结果

由图4 可以看出，当迭代次数达到50 次之后，人眼已难以观察出实验结果与原始图像的差别。

图5 给出了随着迭代次数增加，含噪声图像Camera⁃men 的PSNR 值变化情况。其中，红色圈状轨迹表示未使用加速算法的软阈值模型6 的处理结果，蓝色的“+”号轨迹表示使用加速算法的软阈值模型6 的处理结果。从图中可以看出，图像Cameramen 在使用加速算法求解软阈值模型6 时能更快地收敛到合适的结果。

图6 给出了4 种去噪算法的对比结果。

Fig.5 Variation of the PSNR for the different algorithms on nosed image“Cameramen”图5 含噪声图像Cameramen PSNR 值变化情况

Fig.6 Comparison of the result for fouralgorithms图6 4 种去噪算法结果对比

由图6 可以看出，采用不同算法处理Barbara 与Cam⁃eramen 图像的结果各不相同。其中，使用模型6 处理图像，结果的边缘震荡效果更弱，更符合人的视觉审美。由于高斯陷波滤波器以及掩膜矩阵滤波器对傅里叶频域系数的处理只是进行相应缩放或截断，因此结果的频域系数相较于原始频域系数可能会有较大误差，具体表现为实验结果中实物边缘产生的震荡效应。实验图片处理结果的PSNR 值如表1 所示。

表1 给出了所有实验对象处理结果的PSNR 值。由表1 可以看出，使用模型6 处理图像的结果相比频域滤波器有了较大提升，验证了利用稀疏性逼近的傅里叶频域系数相较于高斯陷波滤波器的滤波结果更接近原始频域系数。由于l1范数的凸性质，在大部分情况下，使用软阈值模型6的结果比使用硬阈值模型6 的结果更逼近全局最优解。

图7 表示2 张含有摩尔纹噪声的彩色自然图像及其处理结果，图片来源于使用智能手机摄像头分别拍摄IPS面板与TN 面板计算机液晶显示器得到的彩色图像。由图7 可以看出，频域滤波器无法在保护原始图像信息的前提下去除摩尔纹噪声。通过观察软阈值模型6 的处理结果可以发现，虽然图像部分细节有所丢失，但对摩尔纹噪声的总体去除效果较为显著。

Table 1 Comparison of PSNR of image restoration results表1 图像恢复结果PSNR 值对比

Fig.7 Denoising results of color images图7 彩色图像去噪结果

4 结语

本文通过分析摩尔纹噪声的傅里叶频域特性，针对传统频域滤波器无法准确恢复低频信号的缺陷，提出一种基于稀疏表示的摩尔纹噪声去除方法。首先通过阈值法确定频域中受噪声影响的区域，然后根据这些区域位置建立掩膜矩阵。由于紧小波框架下自然图像的小波系数具有稀疏性，因此根据压缩感知原理可建立优化模型稀疏逼近原始图像的傅里叶频域系数。通过使用DFISTA 算法的优化模型，求解过程能够收敛到稳定的结果。数值实验结果表明，当噪声影响区域涉及低频区域时，本文提出的图像恢复算法相比中值滤波方法与高斯陷波滤波方法不仅具有更高的峰值信噪比，而且能保留更清晰的细节信息。以含有摩尔纹噪声的彩色自然图像作为实验对象的图像恢复数值实验也证明了该算法的有效性与可行性。