魏淑玉^1,2, 陈欣^1,2, 文永富^1,2, 程灏波^1,2
1. 北京理工大学 深圳研究院, 广东 深圳 518057;
2. 北京理工大学 光电学院, 北京 100081
2020-09-24 收稿, 2020-11-13 录用
基金项目: 深圳市科技创新项目（JCYJ20170817114558026，JCYJ20200109143424977）
^*通讯作者: 程灏波

摘要: 数字全息技术由于其高灵敏度、高准确度、分辨能力强，且再现、存储及传输方便灵活等特点，在微纳结构与生物细胞的测量领域中得到越来越广泛的应用。然而，记录数字全息图受孔径有限等因素的影响，其重构图像会产生一些类似于波纹的边缘误差，导致成像质量降低和CCD视场不能充分利用。本文首先基于数字全息理论对误差产生原因进行了理论分析，并对边缘误差产生影响的数值参数进行了详细讨论，在此基础上，提出了一种周期延拓迭代方法对全息图进行预处理，并使用图像均方差值对迭代结果进行评价。通过计算机仿真与实验，结果表明，该方法能够显著降低边缘误差，提高再现图像质量。
关键词: 数字全息边缘误差周期延拓
Research on Edge Error Suppression Method of Reconstructed Image of Digital Holography
WEI Shuyu^1,2, CHEN Xin^1,2, WEN Yongfu^1,2, CHENG Haobo^1,2
1. Shenzhen Research Institute, Beijing Institute of Technology, Shenzhen 518057, Guangdong, P. R. China;
2. School of Optics and Photonics, Beijing Institute of Technology, Beijing 100081, P. R. China
^*Corresponding author: CHENG Haobo
Abstract: Due to its high sensitivity, high accuracy, strong resolving power, convenient and flexible reproduction, storage and transmission, digital holography technology has been widely used in the measurement of micro-nano structures and biological cells. However, the recording of digital holograms is affected by factors such as limited apertures, and the reconstructed image will produce some edge errors that resemble ripples, resulting in the degradation of image quality and insufficient utilization of CCD field of view. Based on the digital holographic theory, this paper firstly analyzes the causes of errors, and discusses in detail the numerical parameters that affect edge errors. On this basis, a period continuation iterative method is proposed to preprocess the hologram and use the mean square difference of images to evaluate the iteration results. Through computer simulation and experiments, the results show that this method can significantly reduce edge errors and improve the quality of the reproduced image.
Key words: digital holographyedge errorperiod extension
近年来，随着光电成像器件与计算机能力的飞速发展，数字全息技术凭借其高灵敏度、高准确度、分辨能力强，且再现、存储及传输方便灵活等优点^[1-3]，在微元结构与生物细胞的测量领域展现出广阔的应用前景。但是，数字全息图的数值孔径十分有限^[4]，使得其重构图像边缘出现了类似于波纹的误差，导致重构图像质量降低和CCD视场不能充分利用。在光学系统中，往往采用切趾法来降低孔径衍射，即选用切趾窗函数与全息图数据相乘，目前常用的切趾窗函数有Gaussian窗^[5]、汉宁窗^[6]、Tukey窗^[7]、Blackman窗等，利用边缘处透过率从1突变为0改换成由1平滑变为0，实质是降低旁瓣以达到削弱孔径衍射的目的。切趾处理后，虽然能在一定程度上降低孔径衍射影响, 使靠近重构图像中央区域的质量提高，却是以牺牲最边缘信息为代价，而在一定场合，尤其是微生物或者微结构的动态在线测量当中，充分利用CCD获取的数据获得有效信息是十分重要的。因此，既能改善孔径衍射所带来的边缘误差，又能不损失原始数据信息是相当重要的一个研究方向。
为解决上述问题，通常可直接采用补零操作，即将原全息图或者全息图频谱放入更大的零矩阵中，利用新的全息图或者全息图频谱再进行数值计算^[8-10]。由于全息图周围延拓填充的内容均为零，原始边缘处数据仍是在边界处突变成零，所以补零操作对改善孔径衍射所带来的边缘误差效果有限。2001年，美国****Massig等^[11]为降低由离散傅里叶变换产生的边缘误差提出了一种延拓方法，即对全息图的边缘部分数据人为地向外延拓，延拓后在原全息图外缘部分产生新的数据信息，该部分数据信息尽可能与原全息图信息相似，这样在理想状态下只会在新的全息图重构图像边缘产生误差，从而抑制了有效数据边缘误差的产生^[12]。2002年，Dubois等^[13]根据将最高空间频率最小化原则，用一系列复杂数据逐行逐列扩展数字全息图的方法，保留了原始全息图振幅分布，因此能削弱孔径衍射所带来的边缘误差，比较好地保留了边界部分数据，但是，该方法因添加一些不存在信息数据可能会干扰重建。本文在分析此类边缘误差产生基础上，对重构边缘误差进行理论分析，讨论影响重构边缘精度的参数，并分析这些参数与边缘误差之间的变化规律并提出一种周期延拓迭代方法对全息图进行预处理，首先使用补零操作将原始全息图放入更大的零矩阵中，之后对新全息图光强信息反复周期迭代，最后得到的全息图延拓区域处数据与原始数据相似，从而解决了上述问题，并通过仿真与实验证明此方法在抑制边缘误差方面有非常明显的效果。
1 基本理论1.1 边缘误差产生原因理论分析数字全息记录与重构过程如图 1所示。
图 1

图 1 数字全息记录与重构过程

根据菲涅尔衍射理论，图 1中O点物光波在记录平面P点光波的复振幅为O_P=(A₀/r)exp(jkr)，参考光在P点的复振幅分布为R_p=A_Rexp(jkρsinθ)，该点物光波与参考光波干涉之后叠加的干涉条纹强度分布见式(1)。

(1)

用原始参考光R_I照射全息图，假设R_I的复振幅为R_I=A_Rexp(jkρ′sinθ)，根据角谱衍射理论^{[14, 15]}，记传递函数为，使用频谱滤波的方式去除重构物光波共轭项^{[16, 17]}，记频谱滤波截取边界，可以得到P点在像平面重构的离散物光强度分布A_obj，见式(2)。

(2)

式中，Δx、Δy为CCD的记录感光单元间隔，αΔx、βΔy为感光单元大小，α和β为感光单元尺寸填充因子，M、N为感光单元总数；sinc(MΔx*f_x, NΔy*f_y)为CCD截取图像的边界函数，sinc(αΔx*f_x, βΔy*f_y)与comb(Δx*f_x, Δy*f_y)是CCD内部感光单元的特征函数。忽略CCD内部感光单元特征对整体的重构精度影响，记：
则物光波强度可表示为：

(3)

由式(3)可见，在忽略CCD本身内部因素对重构图像影响的前提下，物光波强度的重构精度可能会受到传递函数、CCD截取图像的边界、频谱滤波截取边界的影响。而CCD截取图像的边界与频谱滤波截取边界在重构过程中实际上为sinc函数的傅里叶变换与逆傅里叶变换，因此可能会产生孔径衍射效应，从而造成边缘误差的产生；传递函数在频域中与sinc函数直接相乘，可视为影响sinc函数的一变量，因此可能会对衍射引起振荡波纹的大小产生影响，此外，重构距离d作为传递函数中的一次变量也可能会造成一定影响。综上，本文认为理想情况下重构误差的产生主要跟重构距离d、CCD截取图像矩形边界、频域中截取矩形边界有关。
1.2 重构误差仿真分析从前文可知，在忽略CCD本身内部因素对重构图像影响的前提下，可能影响重构边缘精度的有3项，分别是重构距离d、CCD截取图像边界、频谱滤波截取边界, 只有找到这些参数与边缘误差之间的变化规律，才能选择合适的值对边缘误差进行调控。因此，本文使用控制变量的方法，分别以此3项为变量进行仿真讨论，利用每幅重构图像的中心区域与边缘区域的残差均方值作为边缘误差来评价重构图像边缘质量，残差均方值越低，边缘误差越小。
(1) 其他变量固定，只改变重构距离d对图像进行仿真重构(图 2)。取每个重构图像边界一定数据归一化之后的均方差值作为误差，边缘误差与其重构距离的关系见图 3，可见重构距离的不同会对重构图像边缘造成一定的影响，重构距离越大边缘误差越明显。因此，在设计光学系统时，物体到CCD距离设计尽可能小才能在重构过程中最大程度降低边缘误差。
图 2

图 2 不同距离重构的图像(a) 10 mm; (b) 20 mm; (c) 40 mm

图 3

图 3 重构距离d对重构边缘误差影响

(2) 其他变量固定不变，改变CCD截取图像矩形边界，即改变截取图像矩形边长的像素数，对图像进行仿真重构，如图 4所示。随着矩形边长的逐步增加，边缘误差变化起伏比较大，没有明显的线性关系，因此，这部分影响很难从系统上降低。
图 4

图 4 CCD截取全息图大小对重构图像边缘误差的影响

(3) 其他变量固定不变，改变频域中截取矩形边界。设定截取的形状为正方形，即只改变矩形边长大小，对图像进行仿真重构。截取频谱的区域大小跟重构图像清晰程度与再现细节有很大的联系，如图 5所示，可见选择频谱时的矩形框大小会对重构图像边界造成一定的影响，频谱越大越好，实际情况下，为重构图像细节更加清晰，也是尽可能选择较大的频谱。
图 5

图 5 频域中截取矩形边界大小对边缘误差影响

2 全息光强图周期延拓算法2.1 全息光强图周期延拓算法原理该算法的基本原理(图 6)是直接对全息图的光强进行傅里叶变换，在频域中对全息图高频频谱截断，在吉布斯效应的作用下，旁瓣频谱能量泄漏到补零区域，通过反复的迭代使周围填充区域的能量增强，原始全息图的条纹在此部分得到延拓，因为内容相近，这样重构出来的波纹误差在新图像边缘出现，便能有效抑制边缘误差对有效信息数据的干扰。
图 6

图 6 全息光强图周期延拓算法原理

假设对干涉图A_P进行延拓，步骤如下：
(1) 先对干涉全息图周围进行补零，记此时图像为I₀。可记后一项为影响条纹强制以及相位部位。

(4)

(2) 对I₀傅里叶变换，理想情况下干涉图的频谱主要分为一个零级频谱，两个一级频谱，且零级频谱的强度大于一级频谱。进行频移之后记为(5)。

(5)

(3) 将小于0.01A(f_x, f_y)的频谱均设为零，即对高频频谱截断，再进行逆傅里叶变换得到图像I，在吉布斯效应的作用下，旁瓣频谱能量泄漏到原始边界之外，即补零区域。
(4) 将原图图像A_P重新填充现在的图像，得到I₁。
(5) 再重复以上步骤。随着迭代次数增加，补零区域的能量逐步增加，得到延拓的条纹就越清晰。
2.2 全息光强图延拓算法仿真分析考虑到物象与共轭像的分开，以及消除共轭像问题，本文采用离轴全息方法进行仿真实验。设定仿真环境，激光为635 nm波长，重建距离为25 mm，CCD光敏面间隔为4.65 μm。物光与参考光夹角为2°，生成离轴干涉全息图，使用对全息图光强延拓的方法对全息图进行预处理后，再用角谱法进行重现。仿真得到干涉全息图后对全息图周围进行补零，再用对全息图光强延拓的方法对补完零的全息图进行延拓迭代，得到新的全息图，如图 7所示。为方便对比，分别对(a)、(b)、(c)全息图进行重构，得到3幅重构图像(图 8)，从图 8可以看出，全息图光强延拓后全息图重构出来的图像(c)边缘波纹误差降低到很小，边缘信息比较清晰。
图 7

图 7 仿真得到的全息图(a) 原始全息图，大小为512×512；(b) 补零后得到的全息图，大小为712×712；(c) 通过对全息图光强迭代100次得到的全息图；(d)、(e) 分别为(b)、(c)的局部放大图

图 8

图 8 仿真重构图像结果(a) 直接重构得到的图像; (b) 补零全息图重构出来的图像; (c) 通过全息图光强延拓后全息图重构的图像; (d)、(e)、(f) 分别为(a)、(b)、(c)的局部放大图

为了得到更加直观的数据，以及分析延拓像素数对边缘误差的影响，分别在不同延拓像素数情况下，利用SSIM公式对仿真数据进行处理。SSIM公式是一种衡量两幅图像结构相似度的公式^[18]，利用两幅图的平均值、方差、协方差分别作为亮度、对比度、结构相似度的估计，当两张图一样时其值为1，即最大值为1。表达式见式(6)，其中μ代表平均值，σ代表方差。分别取边缘部分数据用SSIM公式求相似值，结果如图 9所示，可见当延拓像素数超过100时，通过对全息图光强延拓方法重构出来图像SSIM系数最接近于1，对边缘误差抑制效果最好。

(6)

图 9

图 9 不同延拓像素数下3种方法的SSIM值对比

根据前面所分析的内容，随着迭代次数的增加，延拓区域的能量逐步增加，延拓的条纹越加清晰，但是，迭代的能量到达一定值后，再增加迭代次数多余的计算步骤会对计算机造成计算负担而降低计算效率。相应的，延拓像素数的增加也会大大增加计算时间。鉴于此，本文对最佳迭代次数及最佳延拓像素数进行了初步探索。
以延拓像素数为变量，分别观测以迭代次数为变量造成的误差变化，处理数据时仍采用取边缘部分数据用SSIM公式求相似值，考虑到SSIM性质，本文认为SSIM大于0.900得到的图像质量为佳。选取两组数据，先设迭代次数为100，SSIM值大于0.900为终止条件，找到相应的延拓像素数，之后，使用求得的延拓像素数，以及SSIM值大于0.903且相邻两个值小于0.0005为迭代终止条件。另外一组，先设迭代次数为100，SSIM值0.906为终止条件，之后，求最佳迭代次数时设置SSIM值大于0.910且相邻两个值小于0.0005时迭代终止。
通过图 10可知，延拓像素数及迭代次数与SSIM值的关系并不是简单的线性关系，不同的延拓像素数区间，相应的SSIM值及最佳迭代次数不是相同的，但是随着延拓像素数的增加，在相应最佳迭代次数下SSIM值呈增长趋势。(a)、(b)组中，迭代像素数为320时SSIM值大于0.900，为了得到更好的效果，取延拓像素数略大于320的400，即一边延拓约为原图像尺寸的2/5，可以看到迭代次数在30~35之间便满足迭代终止要求；(c)、(d)组中，延拓像素数为620，即一边延拓约为原图像尺寸的3/5，从(d)图中可以看到，迭代次数在30~35以内便满足要求。因此，一般情况下，一边延拓像素最多为原全息图最大尺寸的3/5，并且迭代次数在30~50之间几乎都能满足要求。
图 10

图 10 迭代次数与延拓像素数对延拓后重构边缘误差影响(a) SSIM值大于0.900为终止条件，延拓像素数与SSIM之间的变化关系；(b) 延拓400像素数，迭代次数与SSIM之间的变化关系；(c) SSIM值大于0.906为终止条件，延拓像素数与SSIM之间的变化关系；(d) 延拓620像素数，迭代次数与SSIM之间的变化关系

3 实验分析进行实验分析验证方法的可行性。采用基于迈克逊干涉显微原理反射式数字全息显微实验平台，光路原理图如图 11所示，波长为635 nm的红色激光经过扩束准直后，在分光棱镜的作用下分成两束，一束光到达反射镜M3后被反射回来后，作为参考光到达CCD感光面，另一束照射到物体上被反射回来，到达CCD感光面后，作为物光与参考光发生干涉。
图 11

图 11 反射式数字全息显微实验平台光路原理图

实验样品是Umasch公司生产的TGXYAZ02高精度刻蚀箭头的晶元表面和‘MIN’字母字样的反射式样品。经过调节数字全息显微平台后，得到离轴全息图(图 12)，大小为960×1280，将全息图利用角谱法直接重构后得到图像为图 13(a)、(d)，将原全息图补零到1560×1880后重构得到的图像为图 13(b)、(e)，且光强延拓迭代50次后重构得到的图像为图 13(c)、(f)，很明显(c)、(f)的边缘效果最好，孔径衍射效应得到了很好的抑制。
图 12

图 12 离轴全息图(a)刻蚀箭头晶元；(b)‘MIN’字样的反射式样品

图 13

图 13 实验结果(a)、(d) 直接重构得到的图像; (b)、(e) 补零全息图重构出来的图像; (c)、(f)通过全息图光强延拓后全息图重构出来的图像

考虑到迭代次数在数值较小的30~50便能满足，相比于延拓像素数变化造成的对计算的负担，其影响比较小，因此，对延拓不同像素数造成的边缘误差变化进一步实验。以刻蚀箭头的晶元为样本，取迭代次数固定为50，分别延拓50~2250像素数来观察重构图像边缘误差的变化。取边缘部分数据并分别计算其均方差值评价边缘误差，均方差数值越小，则此组图像边缘质量越好。如图 14所示，总体上，边缘误差随延拓像素数增加而呈降低的趋势，重构图像中类似图像边缘波纹也在逐渐消失，单边延拓500之后的曲线总体变化比较缓慢，大至趋于平稳，也验证了前文结论，一般情况下，一边的延拓像素最多为原全息图最大尺寸的3/5，迭代次数30~50之间几乎都能满足要求。
图 14

图 14 延拓不同像素数造成的边缘误差变化

4 结论本文针对数字全息术中对全息图进行重构时，受孔径衍射效应影响，重构图像边缘产生一些类似于波纹的边缘误差而造成的成像质量降低和CCD视场不能充分利用问题，对边缘误差产生的原因进行了理论分析，并仿真验证了边缘误差的产生与重构距离、CCD截取边框等因素有关；提出了一种周期延拓迭代对全息图进行预处理的方法，通过仿真及实验证明此方法在抑制边缘误差方面有非常明显的效果。

参考文献

[1]	邓丽军. 数字全息显微术及若干应用技术研究[D]. 天津: 南开大学, 2014.
[2]	Cuche E, Marquet P, Depeursinge C. Aperture apodization using cubic spline interpolation: application in digital holographic microscopy[J]. Optics Communications, 2000, 182(1-3): 59-69. DOI:10.1016/S0030-4018(00)00747-1
[3]	Zygmund A, Bary N K. A treatise on trigonometric series[J]. Mathematics of Computation, 1966, 20(93): 191-195. DOI:10.2307/2004319
[4]	Hadley G. Diffraction by apodized apertures[J]. IEEE Journal of Quantum Electronics, 1974, 10(8): 603-608. DOI:10.1109/JQE.1974.1068217
[5]	Tukey J W. An introduction to the calculations of numerical spectrum analysis[J]. Spectra Analysis of Time, 1967, 25-46.
[6]	Hamilton O H, Achieser N I, Hyman C J. Theory of approximation[J]. The American Mathematical Monthly, 1957, 64(9): 686. DOI:10.2307/2309990
[7]	张延曹, 赵建林, 张伟, 等. Tukey窗函数用于数字全息图的切趾研究[J]. 光子学报, 2007, 36(12): 2256-2260.
[8]	Zhang F, Yamaguchi I, Yaroslavsky L P. Algorithm for reconstruction of digital holograms with adjustable magnification[J]. Optics Letters, 2004, 29(14): 1668. DOI:10.1364/OL.29.001668
[9]	Li J C, Tankam P, Peng Z J, et al. Digital holographic reconstruction of large objects using a convolution approach and adjustable magnification[J]. Optics Letters, 2009, 34(5): 572-574. DOI:10.1364/OL.34.000572
[10]	Shin J G, Kim J W, Eom T J, et al. Frequency domain zero padding for accurate autofocusing based on digital holography[J]. Journal of the Korean Physical Society, 2018, 72(1): 57-65. DOI:10.3938/jkps.72.57
[11]	Massig J H, Heppner J. Fringe-pattern analysis with high accuracy by use of the fourier-transform method: theory and experimental tests[J]. Applied Optics, 2001, 40(13): 2081-2088. DOI:10.1364/AO.40.002081
[12]	王雨雷, 吕志伟, 何伟明, 等. 采用一种新的外推延拓法重构激光波前[C]. 第十七届全国激光学术会议论文集, 2005.
[13]	Dubois F, Monnom O, Yourassowsky C, et al. Border processing in digital holography by extension of the digital hologram and reduction of the higher spatial frequencies[J]. Applied Optics, 2002, 41(14): 2621-2626. DOI:10.1364/AO.41.002621
[14]	Guo C S, Zhang L, Rong Z Y, et al. Effect of the fill factor of CCD pixels on digital holograms: comment on the papers[J]. Optical Engineering, 2003, 42(9): 2768-2772. DOI:10.1117/1.1599841
[15]	Kreis T M. Frequency analysis of digital holography with reconstruction by convolution[J]. Optical Engineering, 2002, 41(8): 573-607.
[16]	叶春玲. 数字全息图记录及再现算法的研究[D]. 哈尔滨: 哈尔滨工业大学, 2006.
[17]	刘诚, 刘志刚, 程笑天, 等. 数字滤波法再现电子全息图[J]. 光学学报, 2003, 23(2): 150-154. DOI:10.3321/j.issn:0253-2239.2003.02.005
[18]	Alain Horé, Ziou D. Image quality metrics: PSNR vs. SSIM[C].20th International Conference on Pattern Recognition, ICPR 2010. IEEE Computer Society, 2010.