一种基于压缩感知的太赫兹图像去噪方法

摘要

本发明涉及一种基于压缩感知的太赫兹图像去噪方法，可包括下列步骤：S1、原始太赫兹图像的信号稀疏表示；S2、原始太赫兹图像的信号感知测量，构造测量矩阵；S3、通过采用局部快速傅里叶变换求解最小化全变差TV的重建算法恢复太赫兹图像。本发明可以在最大限度的去除太赫兹图像噪声的同时，保护图像有用的边缘和纹理细节信息。

权利要求

1.一种基于压缩感知的太赫兹图像去噪方法，其特征在于，包括下列步骤：S1、原始太赫兹图像的信号稀疏表示，具体过程为：令 表示一幅像素为N(N＝n×n)的二维灰度太赫兹数字图像，称为原始太赫兹图像，u表示重建后的图像，原始太赫兹图像的信号稀疏表示为 其中，Ψ＝[ψ1,ψ2,......ψN]∈RN×N是RN的小波正交基，x∈RN是内积 的向量，||x||0＝K，K是向量x中的非零个数，K＜＜N，称 是在Ψ域的K稀疏；

S2、原始太赫兹图像的信号感知测量，构造测量矩阵，具体过程为：采用局部傅里叶矩阵作为测量矩阵Φ，采样测量原始太赫兹图像信号 Φ是M×N矩阵，M满足K

S3、通过采用局部快速傅里叶变换求解最小化全变差TV的重建算法恢复太赫兹图像，所述步骤S3的具体过程为：S31、将TV的定义为：TV(u)＝∑|ui+1,j-ui,j|+|ui,j+1-ui,j|，其中(i,j)为图像像素的位置坐标，获取 的局部傅里叶频率测量向量fp(b＝fp)： 其中Fp∈Cp×N表示一个局p部离散傅里叶变换，作为测量矩阵Φ，p为Φ矩阵的行数(M＝p)，ω∈C表示随机噪声；

S32、通过求解 从fp重建图像 其

中， 涵盖了所有像素，Di为2×N元矩阵，Diu表示u在像素i处水平和垂直的局部有限差分，表示图像u的TV的离散量，||ψTu||1表示小波基底Ψ下u的l1范数，τ,μ>0是均衡正则参数，T表示矩阵的转秩；

S33、重做(1)式到 其中， 对于任

意的i，均有wi＝Diu，wi∈R2，z＝ψTu，z∈RN，S34、引入线性约束变矢量v，约束u＝v，惩罚违例||u-v||2，使(2)式具有变量分离结构，为了解决线性约束问题，采用增广拉格朗日函数：其中， φ1(s ,t ,v)＝τ|s|-v(s-t)+(β/2)·|s-t|2 ,β>0是均衡正

则参数；

S35、引入定义： 进行迭代，算法的迭代基础是：采用交替方向求解方法，交替求解LA、最小算子wi、zi和u，具体求解如下：

1)对于固定的u和λ，最小算子zi的计算：

其中，s1(·,τ/β)是一维收敛操作数，其定义为：

2)对于固定的u和λ，最小算子wi的计算：

其中s2(·,1/β)是二维收敛操作数，其定义为：

3)对于固定的(w,z)和λ，LA关于u的最小化，通过求解等价方程(4)实现：其中， D(1),D(2)∈RN×N分别表示水平和垂直的全局有限差分算子，I是原始图像灰度矩阵，用局部快速傅里叶变换和反变换，求解该方程(4)。

说明书

一种基于压缩感知的太赫兹图像去噪方法

技术领域

[0001] 本发明涉及图像处理领域，具体地涉及一种基于压缩感知的太赫兹图像去噪方法。

背景技术

[0002] 太赫兹波是频率范围覆盖0.3THz～10THz(1THz＝1012Hz)的电磁波，该频段介于红外与微波之间，在20世纪80年代中期之前，由于缺乏太赫兹波段高效的发射源和灵敏的探测器，该频段的电磁辐射并未得到深入研究，太赫兹波段曾被称为电磁波谱的“太赫兹空隙”。然而，在过去的20多年中，随着材料科学和激光科学的发展，太赫兹波相关的产生与探测技术都出现较快的发展，太赫兹技术逐渐被成功应用于制药、生物探测、工业无损检测以及国防安全等领域。

[0003] 然而，当前由于太赫兹时域光谱系统的硬件限制，太赫兹成像的质量普遍比较差，诸如，图像存在着大量的噪声，图像的空间分辨率较低。为了提高太赫兹成像质量，有两种常用的方法。一种是从太赫兹成像系统的硬件入手，提高激光发射器的功率，研发相应的高精度光学传感器件，并提高光学系统的加工工艺，其优点是能够从物理层提高太赫兹图像的质量，缺点是受限于激光发射器的功率和光学系统的加工工艺，系统制造困难且造价昂贵；另一种方法是采用数字图像处理技术，对太赫兹图像进行后期的数字图像处理，改善和提高成像质量，满足相应领域对太赫兹图像的需求。一些典型的数字图像去噪的方法，例如，均值滤波、中值滤波、非局部均值滤波(NLM)方法在太赫兹图像去噪过程中会模糊图像的边缘和细节信息，去噪效果差，算法复杂，实时性不高。因此，研发有效的快速太赫兹图像去噪方法是太赫兹图像理论研究以及应用推广的关键因素。

发明内容

[0004] 本发明旨在提供一种基于压缩感知的太赫兹图像去噪方法，以解决现有图像去噪方法无法获得满意的太赫兹图像效果的问题。为此，本发明采用的具体技术方案如下：

[0005] 一种基于压缩感知的太赫兹图像去噪方法，可包括下列步骤：

[0006] S1、原始太赫兹图像的信号稀疏表示，具体过程为：令 表示一幅像素为N(N＝n×n)的二维灰度太赫兹数字图像，称为原始太赫兹图像，u表示重建后的图像，原始太赫兹图像的信号稀疏表示为 其中，Ψ＝[ψ1,ψ2,......ψN]∈RN×N是RN的小波正交基，x∈RN是内积 的向量，||x||0＝K，K是向量x中的非零个数，K＜＜N，称 是在Ψ域的K稀疏；

[0007] S2、原始太赫兹图像的信号感知测量，构造测量矩阵，具体过程为：采用局部傅里叶矩阵作为测量矩阵Φ，采样测量原始太赫兹图像信号 Φ是M×N矩阵，M满足K

[0008] S3、通过采用局部快速傅里叶变换求解最小化全变差TV的重建算法恢复太赫兹图像。

[0009] 进一步的，所述步骤S3的具体过程为：

[0010] S31、将TV的定义为：TV(u)＝∑|ui+1,j-ui,j|+|ui,j+1-ui,j|，其中(i,j)为图像像素p×N的位置坐标，获取u的局部傅里叶频率测量向量fp(b＝fp)： 其中Fp∈C 表示一个局部离散傅里叶变换，作为测量矩阵Φ，p为Φ矩阵的行数(M＝p)，ω∈Cp表示随机噪声；

[0011] S32、通过求解 从fp重建图像其中， 涵盖了所有像素，Di为2×N元矩阵，Diu表示u在像素i处水平和垂直的局部有限差分， 表示图像u的TV的离散量，||ψTu||1表示小波基底Ψ下u的l1范数，τ,μ>0是均衡正则参数，T表示矩阵的转秩；

[0012] S33、重做(1)式到 其中，对于任意的i，均有wi＝Diu，wi∈R2，z＝ψTu，z∈RN，

[0013] S34、引入线性约束变矢量v，约束u＝v，惩罚违例||u-v||2，使(2)式具有变量分离结构，为了解决线性约束问题，采用增广拉格朗日函数：

[0014]

[0015] 其中， φ1(s,t,v)＝τ|s|-v(s-t)+(β/2)·|s-t|2,β>0是均衡
正则参数；

[0016] S35、引入定义： 进行迭代，算法的迭代基础是：

[0017]

[0018] 采用交替方向求解方法，交替求解LA、最小算子wi、zi和u，具体求解如下：

[0019] 1)对于固定的u和λ，最小算子zi的计算：

[0020]

[0021] 其中，s1(·,τ/β)是一维收敛操作数，其定义为：

[0022]

[0023] 2)对于固定的u和λ，最小算子wi的计算：

[0024]

[0025] 其中s2(·,1/β)是二维收敛操作数，其定义为：

[0026]

[0027] 3)对于固定的(w,z)和λ，LA关于u的最小化，通过求解等价方程(4)实现：

[0028]

[0029] 其中， D(1),D(2)∈RN×N分别表示水平和垂直的全局有限差分算子，I是原始图像灰度矩阵，用局部快速傅里叶变换和反变换，求解该方程(4)。

[0030] 本发明采用上述技术方案，具有的有益效果是：本发明可以在最大限度的去除太赫兹图像噪声的同时，保护图像有用的边缘和纹理细节信息。

附图说明

[0031] 图1是本发明实施例的流程图；

[0032] 图2是用不同去噪方法的得到的纯度为99％的氧化铝陶瓷片的太赫兹图像，其中(a)原始的光学图像，(b)实验采集的太赫兹图像，(c)Donoho阈值去噪后的太赫兹图像，(d)广义小波阈值去噪后的太赫兹图像，(e)四阶偏微分方程去噪后的太赫兹图像(f)自适应阈值去噪后的太赫兹图像，(g)Bayes阈值去噪后的太赫兹图像，(h)本发明的方法去噪后的太赫兹图像；

[0033] 图3是用不同去噪方法的得到的纯度为96％的氧化铝陶瓷片(内部含有裂纹缺陷)的太赫兹图像，其中(a)原始的光学图像，(b)实验采集的太赫兹图像，(c)Donoho阈值去噪后的太赫兹图像，(d)广义小波阈值去噪后的太赫兹图像，(e)四阶偏微分方程去噪后的太赫兹图像(f)自适应阈值去噪后的太赫兹图像，(g)Bayes阈值去噪后的太赫兹图像，(h)本发明的方法去噪后的太赫兹图像；

[0034] 图4是用不同去噪方法的得到的纯度为96％的氧化铝陶瓷片(预置多个不同规格尺寸套孔缺陷)的太赫兹图像，其中(a)原始的光学图像，(b)实验采集的太赫兹图像，(c)Donoho阈值去噪后的太赫兹图像，(d)广义小波阈值去噪后的太赫兹图像，(e)四阶偏微分方程去噪后的太赫兹图像(f)自适应阈值去噪后的太赫兹图像，(g)Bayes阈值去噪后的太赫兹图像，(h)本发明的方法去噪后的太赫兹图像。

具体实施方式

[0035] 为进一步说明各实施例，本发明提供有附图。这些附图为本发明揭露内容的一部分，其主要用以说明实施例，并可配合说明书的相关描述来解释实施例的运作原理。配合参考这些内容，本领域普通技术人员应能理解其他可能的实施方式以及本发明的优点。

[0036] 现结合附图和具体实施方式对本发明进一步说明。如图1所示，一种基于压缩感知的太赫兹图像去噪方法，可包括下列步骤：

[0037] S1、原始太赫兹图像的信号稀疏表示，具体过程为：令 表示一幅像素为N(N＝n×n)的二维灰度太赫兹数字图像，称为原始太赫兹图像，u表示重建后的图像，原始太赫兹图像的信号稀疏表示为 其中，Ψ＝[ψ1,ψ2,......ψN]∈RN×N是RN的小波正交基，x∈RN是内积 的向量，||x||0＝K，K是向量x中的非零个数，K＜＜N，称是在Ψ域的K稀疏；

[0038] S2、原始太赫兹图像的信号感知测量，构造测量矩阵，具体过程为：采用局部傅里叶矩阵作为测量矩阵Φ，采样测量原始太赫兹图像信号 Φ是M×N矩阵，M满足K

[0039] S3、通过采用局部快速傅里叶变换求解最小化全变差TV的重建算法恢复太赫兹图像；所述步骤S3的具体过程为：

[0040] S31、将TV的定义为：TV(u)＝∑|ui+1,j-ui,j|+|ui,j+1-ui,j|，其中(i,j)为图像像素的位置坐标，获取 的局部傅里叶频率测量向量fp(b＝fp)： 其中Fp∈Cp×N表示一个局部离散傅里叶变换，作为测量矩阵Φ，p为Φ矩阵的行数(M＝p)，ω∈Cp表示随机噪声；

[0041] S32、通过求解 从fp重建图像其中， 涵盖了所有像素，Di为2×N元矩阵，Diu表示u在像素i处水平和垂直的局部有限差分， 表示图像u的TV的离散量，||ψTu||1表示小波基底Ψ下u的l1范数，τ,μ>0是均衡正则参数，T表示矩阵的转秩；

[0042] S33、重做(1)式到 其中，对于任意的i，均有wi＝Diu，wi∈R2，z＝ψTu，z∈RN，

[0043] S34、引入线性约束变矢量v，约束u＝v，惩罚违例||u-v||2，使(2)式具有变量分离结构，为了解决线性约束问题，采用增广拉格朗日函数：

[0044]

[0045] 其中， φ1(s,t,v)＝τ|s|-v(s-t)+(β/2)·|s-t|2,β>0是均衡
正则参数；

[0046] S35、引入定义： 进行迭代，算法的迭代基础是：

[0047]

[0048] 采用交替方向求解方法(ADM)，交替求解LA、最小算子wi、zi和u，具体求解如下：

[0049] 1)对于固定的u和λ，最小算子zi的计算：

[0050]

[0051] 其中，s1(·,τ/β)是一维收敛操作数，其定义为：

[0052]

[0053] 2)对于固定的u和λ，最小算子wi的计算：

[0054]

[0055] 其中s2(·,1/β)是二维收敛操作数，其定义为：

[0056]

[0057] 3)对于固定的(w,z)和λ，LA关于u的最小化，通过求解等价方程(4)实现：

[0058]

[0059] 其中， D(1),D(2)∈RN×N分别表示水平和垂直的全局有限差分算子，I是原始图像灰度矩阵，用局部快速傅里叶变换和反变换，求解该方程(4)。

[0060] 整个算法的计算机求解伪代码流程如下：

[0061] 设置离散空间步长γ＝1.618

[0062] 1)初始化参数u,λ1,λ2,k；

[0063] u＝u0，λ1＝(λ1)0，λ2＝(λ2)0，k＝0

[0064] 当“不收敛时”，进行下一步；

[0065] 2)通过下面的式子计算z和w：

[0066]

[0067]

[0068] 3)通过求解公式(4)，计算uk+1，这里y的值取决于(w,z)＝(wk+1,zk+1)和λ＝((λ1)k,(λ2)k)；

[0069] 4)通过下式更新λ1和λ2

[0070]

[0071]

[0072] 5)判断是否收敛，收敛的条件：||uk+1-uk||2≤ε(1+||uk||2)，ε>0为足够小，如果不收敛，把k←k+1，重复执行第2步，如果收敛，执行下一步；

[0073] 6)满足收敛条件，程序结束。

[0074] 效果验证

[0075] 为了便于实验比对，我们采用了太赫兹成像实验获取的三种不同规格的氧化铝陶瓷基复合材料具有代表性的太赫兹含噪图像，分别是99％纯度氧化铝标准片，图像大小为128×128；内部有裂纹缺陷的96％氧化铝瑕疵片，图像大小为128×128；预置不同规格套孔缺陷的96％氧化铝瑕疵片，图像大小128×128。为了便于客观比较去噪效果，每幅图像加高斯白噪声，噪声水平参数σ＝5。与此同时，我们选择目前具有代表性的Donoho阈值，广义小波阈值，四阶偏微分方程，贝叶斯(Bayes)阈值和自适应阈值这五种其他的不同图像去噪方法，进行去噪效果比对。实验采用MATLAB R2014a数学工具编程实现，其运行的计算机环境配置为：中央处理器：i5-2400，主频3.1GHz；内存：4096MB；操作系统：微软Windows 10，64位。对以上所提出的六种不同的去噪算法，每一种都进行一百次的实验计算。

[0076] 图像信号局部傅里叶快速变换采样的依据是其傅里叶域中中心放射线的数量，放射线的数量RLs与测量向量的数量M成正比。采样率r定义为：r＝M/N。具体实验参数设置见表1。

[0077] 表1压缩感知去噪方法参数设置

[0078]

[0079] 图2～4给出了实验的结果。从实验结果参数，我们可以看出Donoho阈值，广义小波阈值和四阶偏微分方程的方法由于丢失了一部分有用的图像信息，因此去噪效果比较差；自适应阈值和贝叶斯阈值的方法去噪效果稍强；本发明的方法去噪效果要明显优于其他五种方法，是这六种方法中最优的。

[0080] 为了定量化评价这六种图像去噪方法的效果，我们采用均方误差(MSE)，图像信噪比(PSNR)，均方根(RMS)，以及归一化均方误差(NMSE)这四种典型的评价方法，对不同的图像去噪方法的效果进行比对。此外，我们也比较了这六种不同图像去噪方法的计算复杂度(仿真计算程序运行时间)。以下给出了评价方法的数学公式：

[0081]

[0082]

[0083] 这里MAX是像点的最大值, PSNR的单位是dB。

[0084] 表2.PSNR,MSE,NMSE,RMS以及运行时间的执行效果对比

[0085]

[0086] 从表2可以明确的看到，本发明的去噪方法的MSE、NMSE和RMS值是所有方法中最小的，去噪后的图像与原始的图像差异最小。此外，本发明的方法的PSNR值是这几种方法中最大的(是其他方法的约3-5倍)，这也就意味着本发明的方法可以在去噪过程中最大程度的保留图像的边缘和纹理信息。因此，从实验结果，可以得到，虽然本发明的去噪方法所耗的程序仿真时间要稍长于其他方法，但是去噪效果要显著地优于其他方法。

[0087] 尽管结合优选实施方案具体展示和介绍了本发明，但所属领域的技术人员应该明白，在不脱离所附权利要求书所限定的本发明的精神和范围内，在形式上和细节上可以对本发明做出各种变化，均为本发明的保护范围。