CN105629696A - 一种基于迭代去噪收缩阈值算法的数字全息重构方法 - Google Patents

Abstract

本发明公开了一种基于迭代去噪收缩阈值算法的数字全息重构方法：1、读入原始图像，其振幅为O₀(ξ,η)，利用此原始图像的物光信息获得衍射光复振幅O(x,y)；2、由菲涅尔衍射得到的衍射光复振幅O(x,y)与加入了0和π/2相移量后的参照光R(x,y)相干涉到达全息干板上的光波复振幅分别为U₁(x,y)和U₂(x,y)，经过数值处理形成两幅全息图I₍₁₎和I₍₂₎；3、对两幅全息图I₍₁₎和I₍₂₎进行叠加生成一幅相移全息图，构造传感矩阵A，使得其中Φ为已知的测量矩阵，R表示稀疏基矩阵，y表示得到的观测数据；4、利用迭代去噪收缩阈值算法(IDNST)求解获得重构原始物体的图像。本发明的IDNST算法在(TwIST)算法的基础上引入阈值和正则化参数的双收缩，信噪比得以提高，并且加快了收敛速度，提高了重构精度，使得再现质量达到更优水平。

Description

一种基于迭代去噪收缩阈值算法的数字全息重构方法

技术领域

本发明涉及数字全息技术领域，特别涉及一种基于迭代去噪收缩阈值算法的数字全息重构方法。

背景技术

压缩感知技术应用于数字全息技术中，和传统的光学全息术相比，具有以下优点：1)是由于用电子成像器件代替全息干板来记录全息图，所需的曝光时间可以很短，而没有化学银盐干版的湿处理过程，可连续记录运动物体的各个瞬间过程，有利于实现实时在线全息记录；2)是数字计算得到的再现图样可通过计算机进行定量分析。数字再现全息图得到的是物场的复振幅分布，并可通过计算得到强度和相位分布，易于实现两个或多个全息图的相加减、增减背景图像、叠加图像[1]等操作，而这些在光学全息中很难做到，由此也显现了数字全息技术更具实用性；3)将压缩感知理论与数字全息图再现相结合，消除了传统算法对数字全息重构过程中出现的0极像、共轭像干扰[2，3]的问题，对实现计算全息实时性具有重要的意义。

在后续工作中针对稀疏表示中还可以考虑冗余字典上的稀疏分解理论，进一步减少测量值提高重构效果。目前国内外已经有一系列的工作针对基于压缩感知的数字全息成像问题展开。文献[4]研究了2种不同变换域对全息图稀疏化的影响，根据压缩感知理论直接获取图像的压缩表示，并且利用傅里叶稀疏域下使用正交匹配追踪(OMP)重构算法能获得1.5％的压缩率。文献[5]将压缩感知理论与数字全息图再现相结合提出基于全变差的重构算法，并应用于数字全息图压缩感知全息图重建，消除了传统算法对数字全息重构过程中会出现0级像、共轭像干扰的问题。文献[6]将压缩感知应用在相移全息图中，通过利用可变密度采样结构将更多信号积聚在中心位置，很少的信号分布在全息图的边缘，利用全变分求解最优值，从而获得再现像。文献[7]验证了压缩感知技术从低维数据中产生高维图像的重要性，并且阐明了利用压缩感知2D全息数据中重构出一幅3D层析图。文献[8]详细阐明了基于全变分的两步迭代收缩阈值的重构算法，尤其说明了对于病态的问题，相比IST算法，此法对于目标函数求解最优化值有很快的收敛速度。文献[9]介绍了基于压缩感知的次曝光同轴全息结构，能够捕捉微生物的行为动态，相比在线数字全息术，此方法只需要一次曝光特性。文献[10]主要将压缩感知技术应用在获得多个视图全息投影，并且通过对生成的全息图二次抽样重构3D场景。

现有技术中，数字全息图像重构过程中，收敛速度慢，重构精度低。

参考文献

[1]T.-C.PoonandP.P.Banerjee,ContemporaryOpticalImageProcessingwithMATLAB(Elsevier,2001).

[2]YairRivensonandAdrianStern.”ConditionsforpracticingCompressiveFresnelholography,”Opt.Lett.36,3365-3367(2010).

[3]V.Mico,J.Garc1a,Z.Zalevsky,andB.Javidi,“Phase-shiftingGaborholography,”Opt.Express17,1492–1494(2009).

[4]李科,李军.基于压缩传感的全息图压缩研究[J].华南师范大学学报,2012,44(4):61-65.

[5]简献忠,周海,乔静远,李莹,王佳.基于全变差重构算法的数字全息研究[J].激光技术,

2014,38(2):236-239.

[6]YairRivenson,AdrianSternaandBahram,Javidi.”CompressiveFresnelHolography,”J.DisplayTechnol,vol.6,NO.10,pp.506-509,Oct.2010.

[7]DavidJ.Brady,KerkilChoi,DanielL.Marks,RyoichiHorisakiandSehoonLim.”CompressiveHolography,”Opt.Express17(15),13040-13049.

[8]J.M.Bioucas-DiasandM.A.T.Figueiredo,“AnewTwIST:two-stepiterativeshrinkage/thresholdingalgorithmsforimagerestoration,”IEEETrans.ImageProcess.16,2992–3004(2007).

[9]Y.Rivenson,A.SternandB.Javidi,“Improveddepthresolutionbysingle-exposurein-linecompressiveholography,”Appl.Opt.52,A223–A229(2013).

[10]YairRivenson,AdrianSternandJosephRosen.”Compressivemultipleviewprojectionincoherentholography,”Opt.Express19(7),6109-6118(2011).

发明内容

本发明的目的在于提供一种基于迭代去噪收缩阈值算法的数字全息重构方法；在两步迭代收缩阈值算法(TwIST)的基础上，引入阈值和正则化参数的双收缩，提出了迭代去噪收缩阈值(IDNST)算法；利用前两次迭代的值和不断更新的参数t以及不断收缩的正则化参数来获得新的迭代值，提高了算法的重构精度,缩短了迭代时间，并且加强了抗噪能力。本发明提出的重构算法IDNS,加快了收敛速度，提高了图像的重构精度，使得再现质量达到更优水平。

为了实现上述目的，本发明采用如下技术方案：

一种基于迭代去噪收缩阈值算法的数字全息重构方法，包括以下步骤：

步骤一：读入原始图像，其振幅为O₀(ξ,η)，利用此原始图像的物光信息进行菲涅尔衍射得到距离孔径平面为z的观察平面上任意一点P(x,y)的衍射光复振幅O(x,y)；

步骤二：由菲涅尔衍射得到的衍射光复振幅O(x,y)与加入了0和π/2相移量后的参照光R(x,y)相干涉到达全息干板上的光波复振幅分别为U₁(x,y)和U₂(x,y)，经过数值处理形成两幅全息图I₍₁₎和I₍₂₎；

步骤三：对两幅全息图I₍₁₎和I₍₂₎进行叠加得到相移全息图构造传感矩阵A，使得y＝ΦI＝ΦRO'＝AO'，其中Φ为已知的测量矩阵，R表示稀疏基矩阵，y表示得到的观测数据；

步骤四：利用迭代去噪收缩阈值算法重构出原始图像

进一步的，步骤一中：

$O(x,y)=\frac{1}{i\mathrm{\λ}z}\exp \left(ikz\right)\exp \[\frac{i\mathrm{\π}}{2z}(x^2+y^2)\]F\{O_0(\mathrm{\ξ},\mathrm{\η})\exp \[\frac{ik}{2z}({\mathrm{\ξ}}^2+{\mathrm{\η}}^2)\]\}$

其中，λ是光的波长，k是波数，z为物体所在平面与全息面的距离。

进一步的，步骤二具体包括以下步骤：

步骤2.1：衍射光O(x,y)与参考光R(x,y)相干涉在全息记录平面xy上所叠加的总光场分布U(x,y)为：

U(x,y)＝O(x,y)+R(x,y)

步骤2.2：利用全息记录平面xy上的总光场U(x,y)计算形成的全息图光强分布为：

I＝U(x,y)·U^*(x,y)＝(O+R)(O^*+R^*)

＝|O|²+|R|²+O·R^*+O^*·R

式中I代表全息图I₍₁₎和I₍₂₎；上式中U(x,y)为U₁(x,y)时，计算获得全息图I₍₁₎；U(x,y)为U₂(x,y)时，计算获得全息图I₍₂₎；

步骤2.1中，在参考光R(x,y)的相移量为0时，获得U₁(x,y)；在参考光R(x,y)的相移量为π/2时，获得U₂(x,y)。

进一步的，步骤三中：传感矩阵A：

步骤3.1：将步骤二中得到的两幅全息图I₍₁₎和I₍₂₎叠加得

$\hat{I}(x,y)=\frac{1}{2r}(I_{\left(1\right)}+{\mathrm{iI}}_{\left(2\right)})$

步骤3.2：根据步骤二，CCD上的光强分布I(x,y)具体为

$\begin{array}{c}\hat{I}(x,y)=\mathrm{\α}|U(x,y){|}^2+I_0(x,y)\\ =(\mathrm{\α}+I_0)2\mathrm{\α}\mathrm{Re}(o*r)+\mathrm{\α}(o*r){(o*r)}^{*}\end{array}$

式中，α表示CCD量子系数，I₀(x,y)表示CCD背景噪音。式中第一项(α+I₀)表示全息图的直流分量，可以忽略；第二项Re表示实部，2αRe(o*r)表示物体干涉形成的全息图，其中o和r分别为衍射光和参考光的振幅大小；第三项α(o*r)(o*r)^*表示物体之间对干涉的影响；

步骤3.3：α(o*r)(o*r)^*＜＜2αRe(o*r)，因而可以忽略。因而全息图可以简化为：

$\hat{I}(x,y)=\mathrm{Re}(o*r)=\mathrm{Re}\{F^{-1}\[F\left(o\right)F\left(h\right)\]\}$

式中F、F^-1分别表示二维傅里叶变换和傅里叶逆变换。

步骤3.4：假设CCD的分辨力为N_x×N_y，像素大小为Δx×Δy，被记录的目标空间o(x₀,y₀)被划分为步长为Δx×Δy的N_x×N_y个采样区间，对上式进行离散化，得

$\begin{array}{c}\hat{I}(k\mathrm{\Δ}x,l\mathrm{\Δ}y)=\underset{n}{\overset{N_y}{\Σ}}\underset{m}{\overset{N_x}{\Σ}}o(m\mathrm{\Δ}x,n\mathrm{\Δ}y)*(k\mathrm{\Δ}x-m\mathrm{\Δ}x,l\mathrm{\Δ}y-n\mathrm{\Δ}y)\\ =F^{-1}\{F\[o(m\mathrm{\Δ}x,n\mathrm{\Δ}y)\]F\[r(k\mathrm{\Δ}x-m\mathrm{\Δ}x,l\mathrm{\Δ}y-n\mathrm{\Δ}y)\]\}.\end{array}$

步骤3.5：为了应用压缩感知理论，将目标空间二维矩阵与全息图化为一维向量，定义 $O_{(n-1)\×N_x+m}^{\′}=o(m,n),r_{(n-1)\×N_x+m}^{\′}=r(m\mathrm{\Δ}x,n\mathrm{\Δ}y),I_{(l-1)\×N_x+k}^{\′}=\hat{I}(k,l),$ 则上式简化为

I′＝T_2DWQBO′＝RO′

其中，和O'分别为式中和的简写；式中B＝F_2D是大小为(N_x×N_y)×(N_x×N_y)二维分块对角阵，表示二维离散傅里叶变换；BO'对应式中F[o(mΔx,nΔy,qΔz)]；Q为(N_x×N_y)×(N_x×N_y)的对角矩阵，其[(n-1)×N_x+m]为F[r(kΔx-mΔx,lΔy-nΔy)]矩阵第n行m列元素；F[r(kΔx-mΔx,lΔy-nΔy)]表示菲涅尔变换核r(mΔx,nΔy)的傅里叶变换；W＝I₁,为(N_x×N_y)×(N_x×N_y)单位矩阵；T_2D大小为(N_x×N_y)×(N_x×N_y)，表示二维傅里叶逆变换；R＝T_2DWQB为压缩感知理论的稀疏矩阵，大小为(N_x×N_y)×(N_x×N_y)，则A＝ΦR＝ΦT_2DWQB。

进一步的，步骤四具体包括以下步骤：

步骤4.1：初始化；初始图像O'₀＝0，初始两步迭代次数TwIST_iters＝0，迭代次数t₁＝1，最大迭代次数Mt＝1000，算法终止条件值tolA＝10^-3，正则化收缩因子u＝0.9，计数器s＝1，以及各种临时参数O'₁＝O'₀,O'₂＝O'₀；

步骤4.2：定义目标函数为

$X\left(I\right)=\frac{1}{2}\left|\right|y-\mathrm{\Φ}I|{|}_2^2+{\mathrm{\τ\Θ}}_{TV}\left(I^{\′}\right)$

根据步骤三得

$X\left(I^{\′}\right)=\frac{1}{2}\left|\right|y-{\mathrm{AO}}^{\′}|{|}_2^2+{\mathrm{\τ\Θ}}_{TV}\left(I^{\′}\right)$

其中A是传感矩阵，y表示得到的观测数据，I'代表是步骤三获得的叠加生成的相移全息图；τ代表正则化参数，Θ_TV(I')为总变分函数，其公式如下：

${\mathrm{\Θ}}_{TV}\left(f\right)=\underset{i}{\Σ}\sqrt{{\left({\mathrm{\Δ}}_i^h{I}^{\′}\right)}^2+{\left({\mathrm{\Δ}}_i^v{I}^{\′}\right)}^2}$

和代表对于图像像素值I'的一阶水平和垂直差分操作

${\mathrm{\Δ}}_i^h{I}^{\′}={I^{\′}}_{i,j,k}-{I^{\′}}_{i+1,j,k}$

${\mathrm{\Δ}}_i^v{I}^{\′}={I^{\′}}_{i,j,k}-{I^{\′}}_{i,j+1,k}$

步骤4.3：在最优化问题中将式 $X\left(O^{\′}\right)=\frac{1}{2}\left|\right|y-{\mathrm{AO}}^{\′}|{|}_2^2+{\mathrm{\τ\Θ}}_{TV}\left(I^{\′}\right)$ 转化为最小化式

ψ_λ(O')＝argmin{X}

并计算O'₀对应的目标函数X(O'₀)，并且赋值pref＝X(O'₀)；

步骤4.4：对O'₀进行去噪处理，O'₁＝Γ_λ(O'₀),其中，

Γ_λ(O')＝Ψ_λ(O'+Φ^T(y-ΦO'))

Ψ_λ为去噪软收缩函数，在此引入迭代因子t，则

$t_{k+1}=\frac{1+\sqrt{1+4t_k^2}}{2}$

从而

${O^{\′}}_{k+1}={O^{\′}}_k+\left(\frac{t_k-1}{t_{k+1}}\right)({O^{\′}}_k-{O^{\′}}_{k-1})$

并计算O'_k+1对应的目标函数X(O'_k+1)；

步骤4.5：判断算法是否已经进行过两步迭代；

如果TwIST_iters≠0不成立，则计算O'的目标函数值X(O')；之后再判断条件

X(O'_k)>X(O'_k+1)是否成立；

若X(O'_k)>X(O'_k+1)成立，则s＝s×2；并判断条件s>1000是否成立，若成立，则终止程序，最终输出估计值O'_k，若不成立，则赋值TwIST_iters＝0并转至步骤4.4；

若X(O'_k)>X(O'_k+1)不成立，则赋值TwIST_iters＝TwIST_iters+1，并执行步骤4.6；

如果TwIST_iters≠0成立，则两步估计O'_k+1＝(1-α)O'_k-1+(α-β)O'_k+βΓλ(O'_k)，同时计算O'_k+1的目标函数值X(O')，之后再判断X(O')>pref是否成立；若成立则赋值TwIST_iters＝0，并转至步骤4.4；若不成立则赋值TwIST_iters＝TwIST_iters+1，O'＝O'_k+1,并执行步骤4.6；

步骤4.6：判断算法终止条件；判断条件迭代次数s>1000是否成立，若成立则终止算法，输出最终估计值O'；若不成立则赋值O'_k＝O'_k-1，O'_k+1＝O'_k，并将迭代次数加1，并计算算法终止条件并赋值X(O'_k)＝X(O'_k-1)；

步骤4.7：判断C(O'_k,O'_k-1)≤tolA和t>Mt是否成立：

若都不成立，则收缩正则化参数，τ＝τ×u，转至步骤4.4；

只要有一条件成立，则终止算法；最终重构出原始图像

相对于现有技术，本发明具有以下有益效果：

为了解决全息图像数据在传输过程中占用大量内存并在一定程度上增加设计成本的问题，本发明方法基于压缩感知对全息图压缩采样再通过通信信道传输后重构恢复出原始图像。

因此首先由物体光场分布的复振幅O₀(ξ,η)利用菲涅尔衍射，得到距离孔径平面为z的观察平面衍射光场中任意一点P(x,y)的复振幅O(x,y)；由得到的物光波衍射光O(x,y)与加入了0和π/2相移量后的参照光R(x,y)相干涉到达全息干板上的光波复振幅分别为U₁(x,y)和U₂(x,y)，经过数值处理形成两幅全息图I₍₁₎和I₍₂₎；对两幅全息图I₍₁₎和I₍₂₎进行叠加生成一幅相移全息图构造传感矩阵A，使得其中Φ为已知的测量矩阵，R表示稀疏基矩阵，y表示得到的观测数据；引入总变分函数Θ_TV(x)去掉噪声影响，同时利用迭代去噪收缩阈值(IDNST)重构算法求解获得最后的复原图像本发明在实际应用中大量减少占用内存，实用性强。仿真结果表明，本发明所提出的基于迭代去噪收缩阈值压缩全息方法，能够高概率的恢复出原始图像。

本发明用电子成像器件代替全息干板来记录全息图，所需的曝光时间可以很短，同时将压缩感知理论与数字全息图再现相结合，消除了传统算法对数字全息重构过程中出现的0极像、共轭像干扰的问题，对实现计算全息实时性具有重要的意义；另一方面运用我们所提出的基于迭代去噪收缩阈值(IDNST)重构算法成功地由二维的编码像复原出了原始图像的数据，通过引入阈值和正则化参数的双收缩减小了迭代运行时长，提高了算法的重构精度和抗噪能力，根据图8所示，与TwIST算法相比重构误差明显减小，重构效果达到更优。

附图说明

图1(a)为原始图像，图1(b)和图1(c)为利用原有的TwIST算法生成的全息图1和全息图2；

图2(a)为原始图像，图2(b)和图2(c)为利用本发明IDNST算法生成的全息图1和全息图2；

图3为本发明方法的流程图；

图4为IDNST重构算法流程图；

图5(a)、(b)和(c)为利用IDNST算法重建光的幅值与相位图和复原的图像；

图6(a)、(b)和(c)为利用原有的TwIST算法重建光的幅值与相位图和复原的图像；

图7为采用本发明IDNST算法与原有TwIST算法的信噪比对比图；

图8为采用本发明IDNST算法与原有TwIST算法的均方误差对比图。

图9为数字全息记录与再现的直观立体图。

具体实施方式

一种基于迭代去噪收缩阈值算法的数字全息重构方法，包括以下步骤：

步骤一：由已知物体所在平面，物体自身的物光波由该平面衍射出，设该平面为x₀y₀平面，物体光场分布的振幅O₀(ξ,η)利用菲涅尔衍射，根据惠更斯-菲涅尔子波干涉原理得到距离孔径平面为z的观察平面上任意一点P(x,y)的衍射光复振幅O(x,y)，并将其二次相位因子展开，得到菲涅尔衍射的傅里叶变换表达式；读入原始图像，其振幅为O₀(ξ,η)，利用此原始图像的物光信息进行菲涅尔衍射得到在观测平面上的复振幅O(x,y)；

步骤二：由得到的物光波衍射光复振幅O(x,y)与加入了0和π/2相移量后的参照光R(x,y)相干涉到达全息干板上的光波复振幅分别为U₁(x,y)和U₂(x,y)，经过数值处理形成两幅全息图I₍₁₎和I₍₂₎；

步骤三：对两幅全息图I₍₁₎和I₍₂₎进行叠加生成一幅相移全息图构造传感矩阵A，使得其中Φ为已知的测量矩阵，R表示稀疏基矩阵，y表示得到的观测数据；

步骤四：引入总变分函数Θ_TV(x)去掉随机噪声对重构图像的影响，同时利用基于迭代去噪收缩阈值(IDNST)重构算法求解获得重构出原始图像

步骤一中物体物光波的振幅O₀(ξ,η)利用菲涅尔衍射，得到距离孔径平面为z的观察平面上任意一点P(x,y)的衍射光复振幅O(x,y)，具体包括以下步骤：

步骤1.1：衍射光场在观察平面中任意一点的复振幅表示成在孔径平面S上的面积积分形式：

$O(x,y)=\frac{1}{i\mathrm{\λ}}\underset{S}{\∫\∫}{O}_0(\mathrm{\ξ},\mathrm{\η})K\left(\mathrm{\θ}\right)\frac{\exp \left(ikr\right)}{r}dS$

其中，O₀(ξ,η)是自身物体物光波的振幅，O(x,y)是观察平面衍射光的复振幅，K(θ)是倾斜因子，λ是光的波长，k是波数，z是物体所在平面与全息面的距离，r是物体所在平面上的点到全息面上的点的距离。

步骤1.2：在基尔霍夫边界条件下，孔径平面上孔径S以外的区域内光场为0，因此可以将上式积分扩展到无穷。在傍轴条件下K(θ)近似为1，再利用二项式近似

$r=\sqrt{{(x-\mathrm{\ξ})}^2+{(y-\mathrm{\η})}^2+z^2}\≈z\[1+\frac{1}{2}{\left(\frac{x-\mathrm{\ξ}}{z}\right)}^2+\frac{1}{2}{\left(\frac{y-\mathrm{\η}}{z}\right)}^2\]$

从而可以得到菲涅尔衍射公式

$O(x,y)=\frac{1}{i\mathrm{\λ}z}\exp \left(ikz\right)\underset{\mathrm{\∞}}{\∫\∫}{O}_0(\mathrm{\ξ},\mathrm{\η})\exp (\frac{ik}{2z}\[{(x-\mathrm{\ξ})}^2+{(y-\mathrm{\η})}^2\])d\mathrm{\ξ}d\mathrm{\η}$

步骤1.3：根据傅立叶变换的定义，可以将上式表达为

$O(x,y)=\frac{1}{i\mathrm{\λ}z}\exp \left(ikz\right)\exp \[\frac{i\mathrm{\π}}{2z}(x^2+y^2)\]F\{O_0(\mathrm{\ξ},\mathrm{\η})\exp \[\frac{ik}{2z}({\mathrm{\ξ}}^2+{\mathrm{\η}}^2)\]\}$

其中称 $h(x,y,z)=\frac{1}{i\mathrm{\λ}z}\exp \left(ikz\right)\exp \[\frac{i\mathrm{\π}}{2z}(x^2+y^2)\]$ 是参与记录过程的系统传递函数。

步骤二中分别在参考光R(x,y)中加入了0和π/2相移量后与衍射光复振幅O(x,y)干涉形成两幅同轴全息图I₍₁₎和I₍₂₎的具体步骤如下：

步骤2.1：在参考光R(x,y)的相移量分别为0和π/2时，则衍射光O(x,y)与参考光R(x,y)相干涉在全息记录平面xy上所叠加的总光场分布U(x,y)(相移量为0时通过下式计算获得U₁(x,y)，相移量为π/2时通过下式计算获得U₂(x,y))为：

U(x,y)＝O(x,y)+R(x,y)

步骤2.2：利用全息记录平面xy上的总光场U(x,y)计算形成的全息图光强分布为：

I(x,y)＝U(x,y)·U^*(x,y)＝(O+R)(O^*+R^*)

＝|O|²+|R|²+O·R^*+O^*·R

式中I(x,y)代表全息图I₍₁₎和I₍₂₎(上式中U(x,y)为U₁(x,y)时，计算获得全息图I₍₁₎；U(x,y)为U₂(x,y)时，计算获得全息图I₍₂₎)。上式中，O和R为O(x,y)和R(x,y)的简写。第三项和第四项包含了参考光波与物光波相互干涉后所形成干涉条纹的信息，干涉条纹不仅包含了物光波的振幅信息同时也将相位信息记录在了其中。

步骤2.3：在整个全息图物光波的波前记录过程中，按照一定的数量关系使得入射光光强能够透过全息记录平面，则称该数量关系为记录平面的透过率函数：

T_H(x,y)＝T₀+β·I(x,y)

其中T₀和β都是常数，所以在分析中，我们将无关紧要的T₀和β省略。

步骤三中构造传感矩阵A，具体步骤如下：

步骤3.1：将步骤二中得到的两幅全息图I₍₁₎和I₍₂₎叠加得

$\hat{I}(x,y)=\frac{1}{2r}(I_{\left(1\right)}+{\mathrm{iI}}_{\left(2\right)})$

步骤3.2：根据步骤二，CCD上的光强分布具体为

$\begin{array}{c}\hat{I}(x,y)=\mathrm{\α}|U(x,y){|}^2+I_0(x,y)\\ =(\mathrm{\α}+I_0)2\mathrm{\α}\mathrm{Re}(o*r)+\mathrm{\α}(o*r){(o*r)}^{*}\end{array}$

式中，α表示CCD量子系数，I₀(x,y)表示CCD背景噪音。式中第一项(α+I₀)表示全息图的直流分量，可以忽略；第二项Re表示实部，2αRe(o*r)表示物体干涉形成的全息图，其中o和r分别为衍射光和参考光的振幅大小；第三项α(o*r)(o*r)^*表示物体之间对干涉的影响；

步骤3.3：α(o*r)(o*r)^*＜＜2αRe(o*r)，因而可以忽略。因而全息图可以简化为：

$\hat{I}(x,y)=\mathrm{Re}(o*r)=\mathrm{Re}\{F^{-1}\[F\left(o\right)F\left(h\right)\]\}$

式中F、F^-1分别表示二维傅里叶变换和傅里叶逆变换。

步骤3.4：假设CCD的分辨力为N_x×N_y，像素大小为Δx×Δy，被记录的目标空间o(x₀,y₀)被划分为步长为Δx×Δy的N_x×N_y个采样区间，对上式进行离散化，得

$\begin{array}{c}\hat{I}(k\mathrm{\Δ}x,l\mathrm{\Δ}y)=\underset{n}{\overset{N_y}{\Σ}}\underset{m}{\overset{N_x}{\Σ}}o(m\mathrm{\Δ}x,n\mathrm{\Δ}y)*(k\mathrm{\Δ}x-m\mathrm{\Δ}x,l\mathrm{\Δ}y-n\mathrm{\Δ}y)\\ =F^{-1}\{F\[o(m\mathrm{\Δ}x,n\mathrm{\Δ}y)\]F\[r(k\mathrm{\Δ}x-m\mathrm{\Δ}x,l\mathrm{\Δ}y-n\mathrm{\Δ}y)\]\}.\end{array}$

步骤3.5：为了应用压缩感知理论，将目标空间二维矩阵与全息图化为一维向量，定义 $O_{(n-1)\×N_x+m}^{\′}=o(m,n),r_{(n-1)\×N_x+m}^{\′}=r(m\mathrm{\Δ}x,n\mathrm{\Δ}y),I_{(l-1)\×N_x+k}^{\′}=\hat{I}(k,l),$ 则上式简化为

I′＝T_2DWQBO′＝RO′

其中，和O'分别为式中和的简写。式中B＝F_2D是大小为(N_x×N_y)×(N_x×N_y)二维分块对角阵，表示二维离散傅里叶变换；BO'对应式中F[o(mΔx,nΔy,qΔz)]；Q为(N_x×N_y)×(N_x×N_y)的对角矩阵，其[(n-1)×N_x+m]为F[r(kΔx-mΔx,lΔy-nΔy)]矩阵第n行m列元素；F[r(kΔx-mΔx,lΔy-nΔy)]表示菲涅尔变换核r(mΔx,nΔy)的傅里叶变换；W＝I₁,为(N_x×N_y)×(N_x×N_y)单位矩阵；T_2D大小为(N_x×N_y)×(N_x×N_y)，表示二维傅里叶逆变换；R＝T_2DWQB为压缩感知理论的稀疏矩阵，大小为(N_x×N_y)×(N_x×N_y)，则A＝ΦR＝ΦT_2DWQB。

步骤四中利用IDNST算法重构出两幅全息图的具体步骤如下：

步骤4.1：初始化；初始图像O'₀＝0，初始两步迭代次数TwIST_iters＝0，迭代次数t₁＝1，最大迭代次数Mt＝1000，算法终止条件值tolA＝10^-3，正则化收缩因子u＝0.9，计数器s＝1，以及各种临时参数O'₁＝O'₀,O'₂＝O'₀；

步骤4.2：定义目标函数为

$X\left(I\right)=\frac{1}{2}\left|\right|y-\mathrm{\Φ}I|{|}_2^2+{\mathrm{\τ\Θ}}_{TV}\left(I^{\′}\right)$

根据步骤三得

$X\left(I^{\′}\right)=\frac{1}{2}\left|\right|y-{\mathrm{AO}}^{\′}|{|}_2^2+{\mathrm{\τ\Θ}}_{TV}\left(I^{\′}\right)$

其中A是传感矩阵，y表示得到的观测数据，I'代表是步骤二获得的待重构的数据全息图；τ代表正则化参数，Θ_TV(I')为总变分函数，其公式如下：

${\mathrm{\Θ}}_{TV}\left(f\right)=\underset{i}{\Σ}\sqrt{{\left({\mathrm{\Δ}}_i^h{I}^{\′}\right)}^2+{\left({\mathrm{\Δ}}_i^v{I}^{\′}\right)}^2}$

和代表对于图像像素值I'的一阶水平和垂直差分操作

${\mathrm{\Δ}}_i^h{I}^{\′}={I^{\′}}_{i,j,k}-{I^{\′}}_{i+1,j,k}$

${\mathrm{\Δ}}_i^v{I}^{\′}={I^{\′}}_{i,j,k}-{I^{\′}}_{i,j+1,k}$

步骤4.3：在最优化问题中将式转化为最小化式

ψ_λ(O')＝argmin{X}

并计算O'₀对应的目标函数X(O'₀)，并且赋值pref＝X(O'₀)；

步骤4.4：对O'₀进行去噪处理，O'₁＝Γ_λ(O'₀),其中，

Γ_λ(O')＝Ψ_λ(O'+Φ^T(y-ΦO'))

Ψ_λ为去噪软收缩函数，在此引入迭代因子t，则

$t_{k+1}=\frac{1+\sqrt{1+4t_k^2}}{2}$

从而

${O^{\′}}_{k+1}={O^{\′}}_k+\left(\frac{t_k-1}{t_{k+1}}\right)({O^{\′}}_k-{O^{\′}}_{k-1})$

并计算O'_k+1对应的目标函数X(O'_k+1)；

步骤4.5：判断算法是否已经进行过两步迭代；

如果TwIST_iters≠0不成立，则计算O'的目标函数值X(O')；之后再判断条件X(O'_k)>X(O'_k+1)是否成立；

若X(O'_k)>X(O'_k+1)成立，则s＝s×2；并判断条件s>1000是否成立，若成立，则终止程序，最终输出估计值O'_k，若不成立，则赋值TwIST_iters＝0并转至步骤4.4；

若X(O'_k)>X(O'_k+1)不成立，则赋值TwIST_iters＝TwIST_iters+1，并执行步骤4.6；

如果TwIST_iters≠0成立，则两步估计O'_k+1＝(1-α)O'_k-1+(α-β)O'_k+βΓ_λ(O'_k)，同时计算O'_k+1的目标函数值X(O')，之后再判断X(O')>pref是否成立；若成立则赋值TwIST_iters＝0，并转至步骤4.4；若不成立则赋值TwIST_iters＝TwIST_iters+1，O'＝O'_k+1,并执行步骤4.6；

步骤4.6：判断算法终止条件；判断条件迭代次数s>1000是否成立，若成立则终止算法，输出最终估计值O'；若不成立则赋值O'_k＝O'_k-1，O'_k+1＝O'_k，并将迭代次数加1，并计算算法终止条件并赋值X(O'_k)＝X(O'_k-1)；

步骤4.7：判断C(O'_k,O'_k-1)≤tolA和t>Mt是否成立：

若都不成立，则收缩正则化参数，τ＝τ×u，转至步骤4.4；

只要有一条件成立，则终止算法；最终重构出原始图像

请参阅图1至图4所示，本发明利用PC平台仿真代替干涉仪实验装置获取全息图。硬件为acerAL1916W3GHz、内存3G，软件为Windows7旗舰版，MATLABR2014a版本。仿真中光源采用He-Ne激光器，其波长为632.8nm；原始图像像素为256pixel×256pixel。在IDNST重构算法中，参数设定：在编程过程中，根据经验得 ${\mathrm{\λ}}_1=1e-4,{\mathrm{\λ}}_m=1.k={\mathrm{\λ}}_1/{\mathrm{\λ}}_m,\mathrm{\ρ}\≡1-\sqrt{k}/1+\sqrt{k},$ 从而得α＝ρ²+1和β＝2α/λ₁+λ_m。将图5和图6的(c)图进行对照可以看出在相同的初始化条件下本文提出的IDNST算法恢复效果更佳。根据图7和图8，可以看出IDNST算法迭代时间简短，信噪比更大并且均方误差更小。

Claims (5)

1.一种基于迭代去噪收缩阈值算法的数字全息重构方法，其特征在于，包括以下步骤：

步骤一：读入原始图像，其振幅为O₀(ξ,η)，利用此原始图像的物光信息进行菲涅尔衍射得到距离孔径平面为z的观察平面上任意一点P(x,y)的衍射光复振幅O(x,y)；

步骤二：由菲涅尔衍射得到的衍射光复振幅O(x,y)与加入了0和π/2相移量后的参照光R(x,y)相干涉到达全息干板上的光波复振幅分别为U₁(x,y)和U₂(x,y)，经过数值处理形成两幅全息图I₍₁₎和I₍₂₎；

步骤三：对两幅全息图I₍₁₎和I₍₂₎进行叠加生成一幅相移全息图构造传感矩阵A，使得其中Φ为已知的测量矩阵，R表示稀疏基矩阵，y表示得到的观测数据；

步骤四：利用迭代去噪收缩阈值算法重构出原始图像

2.根据权利要求1所述的一种基于迭代去噪收缩阈值算法的数字全息重构方法，其特征在于，步骤一中：

$O(x,y)=\frac{1}{i\mathrm{\λ}z}\exp \left(ikz\right)\exp \[\frac{i\mathrm{\π}}{2z}(x^2+y^2)\]F\{O_0(\mathrm{\ξ},\mathrm{\η})\exp \[\frac{ik}{2z}({\mathrm{\ξ}}^2+{\mathrm{\η}}^2)\]\}$

其中，λ是光的波长，k是波数。

3.根据权利要求1所述的一种基于迭代去噪收缩阈值算法的数字全息重构方法，其特征在于，步骤二具体包括以下步骤：

步骤2.1：衍射光O(x,y)与参考光R(x,y)相干涉在全息记录平面xy上所叠加的总光场分布U(x,y)为：

U(x,y)＝O(x,y)+R(x,y)

步骤2.2：利用全息记录平面xy上的总光场U(x,y)计算形成的全息图光强分布为：

I＝U(x,y)·U^*(x,y)＝(O+R)(O^*+R^*)

＝|O|²+|R|²+O·R^*+O^*·R

式中I代表全息图I₍₁₎和I₍₂₎；上式中U(x,y)为U₁(x,y)时，计算获得全息图I₍₁₎；U(x,y)为U₂(x,y)时，计算获得全息图I₍₂₎；

步骤2.1中，在参考光R(x,y)的相移量为0时，获得U₁(x,y)；在参考光R(x,y)的相移量为π/2时，获得U₂(x,y)。

4.根据权利要求1所述的一种基于迭代去噪收缩阈值算法的数字全息重构方法，其特征在于，步骤三中：构造传感矩阵A的步骤为：

步骤3.1：将步骤二中得到的两幅全息图I₍₁₎和I₍₂₎叠加得

$\hat{I}(x,y)=\frac{1}{2r}(I_{\left(1\right)}+{\mathrm{iI}}_{\left(2\right)})$

步骤3.2：根据步骤二，CCD上的光强分布具体为

$\begin{array}{c}\hat{I}(x,y)=\mathrm{\α}|U(x,y){|}^2+I_0(x,y)\\ =(\mathrm{\α}+I_0)+2\mathrm{\α}\mathrm{Re}(o*r)+\mathrm{\α}(o*r){(o*r)}^{*}\end{array}$

式中，α表示CCD量子系数，I₀(x,y)表示CCD背景噪音；式中第一项(α+I₀)表示全息图的直流分量；第二项Re表示实部，2αRe(o*r)表示物体干涉形成的全息图，其中o和r分别为衍射光和参照光的振幅大小；第三项α(o*r)(o*r)^*表示物体之间对干涉的影响；

步骤3.3：α(o*r)(o*r)^*＜＜2αRe(o*r)，将全息图简化为：

$\hat{I}(x,y)=\mathrm{Re}(o*r)=\mathrm{Re}\{F^{-1}\[F\left(o\right)F\left(h\right)\]\}$

式中F、F^-1分别表示二维傅里叶变换和傅里叶逆变换；

步骤3.4：CCD的分辨力为N_x×N_y，像素大小为Δx×Δy，被记录的目标空间o(x₀,y₀)被划分为步长为Δx×Δy的N_x×N_y个采样区间，对上式进行离散化，得

$\begin{array}{c}\hat{I}\left(k\mathrm{\Δ}x,l\mathrm{\Δ}y\right)=\underset{n}{\overset{N_y}{\Σ}}\underset{m}{\overset{N_x}{\Σ}}o\left(m\mathrm{\Δ}x,n\mathrm{\Δ}y\right)*r\left(k\mathrm{\Δ}x-m\mathrm{\Δ}x,l\mathrm{\Δ}y-n\mathrm{\Δ}y\right)\\ =F^{-1}\left\{F\[o\left(m\mathrm{\Δ}x,n\mathrm{\Δ}y\right)\]F\[r\left(k\mathrm{\Δ}x-m\mathrm{\Δ}x,l\mathrm{\Δ}y-n\mathrm{\Δ}y\right)\]\right\}\end{array}$

步骤3.5：将目标空间二维矩阵与全息图化为一维向量，定义 $r_{(n-1)\×N_x+m}^{\′}=r(m\mathrm{\Δ}x,n\mathrm{\Δ}y),$ $I_{(l-1)\×N_x+k}^{\′}=\hat{I}(k,l)$ 则上式简化为

I′＝T_2DWQBO′＝RO′

其中，和O'分别为式中和的简写；B＝F_2D是大小为(N_x×N_y)×(N_x×N_y)二维分块对角阵，表示二维离散傅里叶变换；BO'对应式中F[o(mΔx,nΔy,qΔz)]；Q为(N_x×N_y)×(N_x×N_y)的对角矩阵，其[(n-1)×N_x+m]为F[r(kΔx-mΔx,lΔy-nΔy)]矩阵第n行m列元素；F[r(kΔx-mΔx,lΔy-nΔy)]表示菲涅尔变换核r(mΔx,nΔy)的傅里叶变换；W＝I₁,为(N_x×N_y)×(N_x×N_y)单位矩阵；T_2D大小为(N_x×N_y)×(N_x×N_y)，表示二维傅里叶逆变换；R＝T_2DWQB为压缩感知理论的稀疏矩阵，大小为(N_x×N_y)×(N_x×N_y)，则传感矩阵A＝ΦR＝ΦT_2DWQB。

5.根据权利要求1所述的一种基于迭代去噪收缩阈值算法的数字全息重构方法，其特征在于，步骤四具体包括以下步骤：

步骤4.1：初始化；初始图像O'₀＝0，初始两步迭代次数TwIST_iters＝0，迭代次数t₁＝1，最大迭代次数Mt＝1000，算法终止条件值tolA＝10^-3，正则化收缩因子u＝0.9，计数器s＝1，以及各种临时参数O'₁＝O'₀,O'₂＝O'₀；

步骤4.2：定义目标函数为

$X\left(I\right)=\frac{1}{2}\left|\right|y-\mathrm{\Φ}I|{|}_2^2+{\mathrm{\τ\Θ}}_{TV}\left(I^{\′}\right)$

根据步骤三得

$X\left(I^{\′}\right)=\frac{1}{2}\left|\right|y-{\mathrm{AO}}^{\′}|{|}_2^2+{\mathrm{\τ\Θ}}_{TV}\left(I^{\′}\right)$

其中A是传感矩阵，y表示得到的观测数据，I'代表是步骤三获得的叠加生成的相移全息图；τ代表正则化参数，Θ_TV(I')为总变分函数，其公式如下：

${\mathrm{\Θ}}_{TV}\left(f\right)=\underset{i}{\Σ}\sqrt{{\left({\mathrm{\Δ}}_i^h{I}^{\′}\right)}^2+{\left({\mathrm{\Δ}}_i^v{I}^{\′}\right)}^2}$

和代表对于图像像素值I'的一阶水平和垂直差分操作

${\mathrm{\Δ}}_i^h{I}^{\′}={I^{\′}}_{i,j,k}-{I^{\′}}_{i+1,j,k}$

${\mathrm{\Δ}}_i^v{I}^{\′}={I^{\′}}_{i,j,k}-{I^{\′}}_{i,j+1,k}$

步骤4.3：在最优化问题中将式转化为最小化式

ψ_λ(O')＝argmin{X}

并计算O'₀对应的目标函数X(O'₀)，并且赋值pref＝X(O'₀)；

步骤4.4：对O'₀进行去噪处理，O'₁＝Γ_λ(O'₀),其中，

Γ_λ(O')＝Ψ_λ(O'+Φ^T(y-ΦO'))

Ψ_λ为去噪软收缩函数，在此引入迭代因子t，则

$t_{k+1}=\frac{1+\sqrt{1+4t_k^2}}{2}$

从而

${O^{\′}}_{k+1}={O^{\′}}_k+\left(\frac{t_k-1}{t_{k+1}}\right)({O^{\′}}_k-{O^{\′}}_{k-1})$

并计算O'_k+1对应的目标函数X(O'_k+1)；

步骤4.5：判断算法是否已经进行过两步迭代；

如果TwIST_iters≠0不成立，则计算O'的目标函数值X(O')；之后再判断条件X(O'_k)>X(O'_k+1)是否成立；

若X(O'_k)>X(O'_k+1)成立，则s＝s×2；并判断条件s>1000是否成立，若成立，则终止程序，最终输出估计值O'_k，若不成立，则赋值TwIST_iters＝0并转至步骤4.4；

若X(O'_k)>X(O'_k+1)不成立，则赋值TwIST_iters＝TwIST_iters+1，并执行步骤4.6；

如果TwIST_iters≠0成立，则两步估计O'_k+1＝(1-α)O'_k-1+(α-β)O'_k+βΓ_λ(O'_k)，同时计算O'_k+1的目标函数值X(O')，之后再判断X(O')>pref是否成立；若成立则赋值TwIST_iters＝0，并转至步骤4.4；若不成立则赋值TwIST_iters＝TwIST_iters+1，O'＝O'_k+1,并执行步骤4.6；

步骤4.6：判断算法终止条件；判断条件迭代次数s>1000是否成立，若成立则终止算法，输出最终估计值O'；若不成立则赋值O'_k＝O'_k-1，O'_k+1＝O'_k，并将迭代次数加1，并计算算法终止条件 $C({O^{\′}}_k,{O^{\′}}_{k-1})=\frac{|X\left({O^{\′}}_k\right)-X\left({O^{\′}}_{k-1}\right)|}{X\left({O^{\′}}_{k-1}\right)},$ 并赋值X(O'_k)＝X(O'_k-1)；

步骤4.7：判断C(O'_k,O'_k-1)≤tolA和t>Mt是否成立：

若都不成立，则收缩正则化参数，τ＝τ×u，转至步骤4.4；

只要有一条件成立，则终止算法；最终重构出原始图像