数字差分-积分快速相位解包裹算法研究*

王子硕 刘磊2)† 刘晨博 刘珂 钟志2) 单明广2)‡

1) (哈尔滨工程大学信息与通信工程学院,哈尔滨 150001)

2) (哈尔滨工程大学,先进船舶通信与信息技术工信部重点实验室,哈尔滨 150001)

数字全息技术是目前应用最广泛的定量相位成像技术之一,但是当测量相位较大的物体时,需要解包裹算法才能计算出正确的相位信息.目前,已有的解包裹算法均面临计算量巨大、耗时的问题.为了解决上述问题,本文基于傅里叶变换相位恢复算法,利用复振幅相位信息的完整性,提出针对薄相位和连续大相位的基于数字差分-积分的快速直接解包裹算法.该算法首先通过基本的傅里叶变换相位恢复算法操作后,得到含有物体完整相位的复振幅信息;随后,从中提取两幅子复振幅信息,并将二者相除,再相位提取出其中信息,便可得到一个物体真实相位差分信息;最后,沿差分方向对提取的相位差分信息进行积分,便可得到解包裹后的相位信息.同时通过仿真与具体实验对该算法进行了验证.结果表明,本文算法可以实现快速准确的解包裹相位直接恢复.

1 引言

数字全息技术是目前应用最广泛的定量相位成像技术之一,它具有无需标记、大视场、非接触和无损等优点.目前为止,国内外众多课题组针对数字全息技术开展了大量的研究工作,并取得了丰富的成果[1-11].其中,相位重建算法[12,13]是数字全息领域的一个极重要的研究方向.

数字全息技术按照重建算法的区别可分为同轴数字全息[14-16]和离轴数字全息[9-11]两类全息技术.同轴数字全息一般需要采集多幅同轴全息图消除共轭像,因此耗时较长,而且对实验装置要求很高,实验过程中实现难度较大.而离轴数字全息技术仅需一张离轴数字全息图便可完成相位恢复,使其成为最具效率的数字全息技术.学者曾提出一系列基于傅里叶变换的相位恢复算法[12,13]对离轴数字全息进行快速且准确相位恢复.基于傅里叶变换的相位恢复算法通常都包含二维傅里叶变换、带通滤波、二维逆傅里叶变换和求取相位等步骤.其中带通滤波会损失很多高频信息,降低了成像的分辨率.因此,也有学者研究同-离轴混合数字全息恢复方法,以兼顾离轴数字全息和同轴数字全息的优势[17,18].但是,无论采用哪种数字全息技术进行测量,由于相位的周期性,仅能提取包裹的相位信息.因此,那些算法恢复出的相位存在着包裹问题,需要解包裹算法对包裹相位进行解包裹才能得到准确的相位信息.目前,解包裹算法主要分为路径跟踪法[19-21]、最小范数法[22-24]和基于深度学习方法[25-27]三大类.其中,枝切法是最经典的路径跟踪解包裹算法,虽然已有若干改进的基于枝切法的相位解包裹算法被提出,但是这些解包裹算法易受局部噪声影响,会出现拖影的解包裹失败的现象.而最小二乘算法利用全局最小化来进行相位解包裹,该方法可以得到整幅解包裹结果,但是会平滑相位信息,降低重建结果的精度.最小二乘算法可分为无加权值[22]和有加权值[23]两种,其中权值代表噪声对像素的干扰程度,权值越大,干扰程度越低;反之亦然.为了提升解包裹的精度,Wei 等[24]将枝切法和最小二乘法相结合,提出了留数点校准最小二乘解包裹算法,该方法利用留数点作为最小二乘算法的加权值,对含噪声相位图表现出卓越的抗噪性,但需多次迭代运算才可求取最优结果,运算时间较长.需要注意的是,不论是路径跟踪法还是最小范数法,重建过程都很复杂,计算量很大,无法满足快速相位重建的实际需求.此外,基于深度学习的相位解包裹算法发展迅速,但是这类算法不仅需要事先采集处理大量数据作为数据集,还需要搭建网络进行长时间的训练才能生效,增加了算法实现的难度,也很难对不同全息系统采集到的数据实现高质量的解包裹.

为了解决数值解包裹算法计算量大、耗时长的问题,本文提出了一种针对薄相位和连续大相位的基于数字差分-积分的快速直接解包裹算法.该算法利用基于傅里叶变换相位恢复算法在逆傅里叶变换所得二维复数矩阵中含有完整物体相位信息的特点,通过数字位移提取和复数相除等操作,得到一幅含有物体相位差分信息的复振幅矩阵,求取其相位信息便可得到无包裹的相位差分信息,最后,沿差分方向对提取的相位差分信息进行积分,便可得到解包裹后的相位信息.该算法避免了局部或全局优化过程,因此在保证解包裹质量的同时,大幅降低了计算量,提高了重建速度.

2 理论方法

假设采集到一幅离轴数字全息图,表达如下式:

其中,A(x,y)表示全息图所含直流量,B(x,y) 表示全息图调制量,φ(x,y)表示物体相位信息,k1和k2表示沿x和y方向的载波.

对全息图做二维傅里叶变换、带通滤波和逆傅里叶变换,表达式如下:

其中,FT[·]表示二维傅里叶变换操作,IFT[·]表示二维逆傅里叶变换操作,BPF 表示带通滤波器.已有的傅里叶变换的相位恢复算法是直接提取复数矩阵的相位信息,得到包裹的相位信息,计算过程如下:

其 中,exp[i(k1x+k2y)] 是对无样品全息图进行(2)式操作得到的系统所含载波信息和背景相位;ANGLE(·)表示取相位操作.

根据傅里叶变换的基本特性可知,在(2)式求取的复振幅数据所含相位信息应为完整的物体相位信息,但是由于相位的周期性和相位提取算法的局限性,(3)式仅能恢复出(-π,π]的包裹相位φ′(x,y),必须使用解包裹算法才能完成相位重建.但是,目前已有的解包裹算法大多计算量大、耗时长.

为了简化解包裹过程,提高相位重建速度,在经(2)式求出复振幅信息exp{i[φ(x,y)+k1x+k2y]}

后,并不直接求取其相位信息,而是从中提取两幅具有一个像素位移量的尺寸略小的子复振幅信息: 从第1 列像素开始提取得到复振幅信息

从第2 列像素开始提取尺寸与II′相同的复振幅信息

并将I′(x,y) 与II′′(x,y) 相除,可得含待测相位沿y方向差分相位信息 ∆yφ(x,y),计算如下:

对(4)式所得复振幅信息求取其相位信息:

其中exp(-ik2)是对一幅无样品全息进行相同操作得到的一幅参考复振幅矩阵.由于相位的周期性和相位提取算法的局限性,所求∆yφ(x,y)也处于(-π,π],因此只要所求相位φ(x,y)相位中相邻相位间差值处于(-π,π],∆yφ(x,y)便能表达真实的相位变化,因此通过将(5)式计算所得∆yφ(x,y)沿y方向积分,便可得真实的相位信息.

根据上述的分析,本文算法的计算流程(图1)总结如下:

图1 快速解包裹算法流程图Fig.1.Algorithm flow diagram for fast phase unwrapping.

步骤1二维傅里叶变换.对全息图进行二维傅里叶变换,得其频谱图.

步骤2带通滤波、逆傅里叶变换.对步骤1中得到的傅里叶频谱进行带通滤波、逆傅里叶变化,得到相位中仅含物体相位信息的复振幅矩阵.

步骤3提取生成两幅新的具有数字位移的子复振幅信息.

步骤4复振幅相除及相位提取.将步骤2 和步骤3 中得到的复振幅矩阵相除,得到一幅新的复振幅信息,并求出其相位信息.

步骤5积分操作.对步骤4 中得到的相位信息沿步骤3 中像素平移方向反向积分,即可得到解包裹的相位信息.

3 仿真结果及分析

为了验证所提方法的可行性,基于CPU i7-87 50H@2.2 GHz,32 GB 内存的硬件平台对该方法进行计算机模拟验证.仿真生成尺寸为500 pixel×501 pixel 的Peak 的相位信息,如图2(a)所示.使用波长为532 nm 的光源,生成离轴全息图如图2(b)所示.使用传统的复振幅相除的傅里叶变换相位恢复算法,得到恢复相位如图2(c)所示.可见,使用传统算法恢复出的相位包裹,必须使用解包裹算法才能得到准确的解包裹信息.分别使用枝切法和最小二乘法对图2(c)进行解包裹,得到解包裹相位为图2(d)和图2(e)所示.然后,使用本文算法得到恢复相位如图2(f)所示.可见,本文算法恢复出来的相位已经是解包裹相位.为了显示统一,仅在文中显示尺寸为500 pixel×500 pixel 的区域.

图2 (a) Peak 值相位型物体及其(b)离轴全息图;(c)傅里叶变换法的恢复结果;使用(d)枝切法和 (e)最小二乘法的解包裹结果;(f)本文所提算法的恢复结果Fig.2.(a) Phase object with peak values and (b) its off-axis hologram;(c) the retrieved result of the Fourier transform;the unwrapped results by (d) branch-cut and (e) least square;(f) the retrieved result of this proposed algorithm.

将3 个恢复出的相位转换成对应的光程差,并提取黑色虚线标注位置的数据展示在图3 中,三者数据完全重合.经计算3 种恢复算法恢复结果标准差均为0.5371 nm.

图3 图2(d)–(f)中黑色虚线所标剖面数据Fig.3.One-dimensional (1D) phase profile along the black dashed lines in Figs.2 (d)–(f).

此外,为了更全面评价所提算法在相位恢复方面的速度优势,使用枝切法、最小二乘法和该算法对多种尺寸(256 pixel×256 pixel,512 pixel×512 pixel 和1024 pixel×1024 pixel)的全息图进行恢复,所能达到的恢复速度如表1 所列.在处理尺寸为256 pixel×256 pixel 的全息图时,相较于枝切法和最小二乘法,本文算法分别提升了7 倍和2 倍,并且随着全息图尺寸的增大,速度提升越大.可见本文算法在保证恢复质量的同时,大幅缩短计算时间.

表1 三种算法比较Table 1. Comparison of three algorithms.

上述数值仿真实验验证了本文所提算法的有效性,不仅可以直接得到解包裹的真实相位信息,而且恢复质量和精度也与传统算法一致.此外,由于仅需一幅数字位移的复振幅矩阵和一维积分等操作,整个算法所需时间远低于目前已有的单波长解包裹算法.

4 实验结果及分析

使用课题组最近提出的离轴数字全息系统[28]采集全息图,其中使用波长为632.8 nm 的 He-Ne激光器作为光源,并使用 1280 pixel×1024 pixel(像素尺寸为 4.8 μm×4.8 μm)的 CCD 相机记录全息图.为了验证该技术的有效性,首先测量高度约为 580.22 nm(由BRUKER 原子力显微镜测量所得)的相位台阶物体,该相位台阶由折射率为1.5168 的BK7 玻璃制成.在本次实验中,采集的离轴全息图如图4(a)所示,使用基于傅里叶变换的相位恢复算法进行相位重建,恢复出的相位图如图4(b)所示.由于引起光程差小于使用波长,所以并未引起包裹.使用本文所提的算法,得到恢复结果如图4(c)所示.

图4 (a) 全息图;(b)傅里叶变换法恢复结果;(c)本文算法恢复结果Fig.4.(a) A hologram;(b) the retrieved result by Fourier transform;(c) the retrieved result by the proposed algorithm.

将两个恢复结果中白色虚线标注的相位信息提取出来,计算出对应的高度信息,一维高度剖面如图5 所示,二者恢复结果完全重合.高度剖面的高度差为583.24 nm,与原子力显微镜测量结果基本吻合.此外,一维高度剖面的标准差(standard deviation,SD)分别为20.997 nm 和20.997 nm.本实验验证了该技术可以实现精确测量.

图5 图4(b)和图4(c)中白色虚线所标剖面数据Fig.5.1D phase profile along the black dashed lines in Fig.4(b) and (c).

为了进一步验证本文算法对连续大相位物体的相位重建能力,使用本文算法对酒精蒸发过程中某一时刻的状态进行处理.图6(a)显示了记录的酒精蒸发过程中的全息图之一,图6(b),(c)显示了经过枝切法和最小二乘解包裹后的两个相位信息,图6(d)显示了本文算法直接恢复的相位信息.为了更好地验证本文算法的有效性,将图6(b)–(d)中白色虚线标注的数据提取出来,展示在图7中,二者高度一致.

图7 图6(b)–(d)中白色虚线所标剖面数据Fig.7.1D phase profile along the black dashed lines in Fig.6 (b)-(d).

通过上述两个实验验证了所提算法能够对具体实验中所采集到富含噪声的全息图进行有效正确的相位解包裹,充分验证了本文算法的有效性和抗噪能力.

5 结论

为了解决单波长数字全息解包裹算法计算量大、速度慢等问题,本文利用复振幅中相位信息完整性的特点,通过数字位移结合复振幅相除,提出数字差分-积分的快速解包裹算法.本文从数值仿真和具体实验两方面验证了所提算法,不仅恢复质量和精度与已有算法高度一致,而且恢复速度也得到了大幅提升.并且,本文算法简单便于实现,为数字全息的快速、高质量重建提供有利的条件.