连续太赫兹波双物距叠层定量相衬成像*

王大勇 李兵 戎路† 赵洁 王云新 翟长超

1) (北京工业大学应用数理学院,北京 100124)

2) (北京市精密测控技术与仪器工程技术研究中心,北京 100124)

针对连续太赫兹波叠层成像重建算法收敛较为迟滞的问题,提出一种连续太赫兹波双物距叠层成像方法及相关重建算法,使用不同记录距离形成的差异化衍射图幅值作为重建算法记录平面的约束条件,增加了记录平面数据多样性和衍射信息冗余度.仿真结果表明,本方法可以加快算法收敛速率,有效减少迭代次数,提高连续太赫兹波定量相衬成像计算效率.随后构建了基于2.52 THz光泵连续太赫兹激光器的双物距叠层成像实验装置,应用双物距记录方法及改进算法重建获得了聚丙烯基字母图案样品的幅值和相位分布,结果表明改进方法可以减少算法迭代次数,提升计算效率,同时改善相位像重建结果保真度.

1 引 言

太赫兹(terahertz,THz)波是指位于电磁波谱中红外与微波之间的电磁波,波长范围为30 μm—3 mm,频率范围为 0.1—10.0 THz.相较于微波,太赫兹波具有更高的成像分辨率.更重要的是,太赫兹波所具备的穿透性、惧水性、低光子能量等独特传播特性,使其在医学成像[1−3]、反恐安检[4,5]和工业产品质量检测[6,7]和无损检测[8,9]等相关领域得到了广泛应用.连续太赫兹波相衬成像是近年来新兴的太赫兹成像研究方向,主要包括单点探测器逐点扫描采集方式和面阵探测器全场记录方式,前者数据采集速率较低,无法对大尺寸样品进行快速成像,同时重建像分辨率受太赫兹波束直径限制[10];后者包括太赫兹数字全息成像与叠层成像.其中,连续太赫兹波数字全息成像方法[11−13]需要利用物光波与参考光波干涉形成全息图,通过重建算法重建样品的复振幅信息.离轴数字全息成像光路较为复杂,同轴数字全息对穿过样品的直透光和衍射光能量比要求较高,同时样品尺寸受光斑直径及探测器靶面尺寸限制.连续太赫兹波叠层成像[14,15]是相干衍射成像方式的一种,通过移动照明探针或探针光束,探测器可记录样品一系列相邻位置的衍射图样信息,利用相位复原算法重建样品的复振幅信息.该方法不仅具有光路简单、结构紧凑、充分利用探测器空间带宽积的优势,相较于同轴数字全息,叠层成像系统可对大尺寸样品扫描成像,非常适合于大尺寸生物样品及非透明高分子材料(例如聚四氟乙烯、高密度聚乙烯等)的检测.

在重建过程中,共轭梯度法 (conjugate gradient)[16]、差分映射法 (difference map)[17]及扩展叠层衍射成像迭代引擎算法(extended ptychographic iterative engine,ePIE)[18]被用来同时重建样品复振幅分布与探针函数分布.其中,ePIE算法程序相较于其他算法具有流程简单、可高效使用计算机内存、算法运行过程中调用参数少,以及可使用基于图形显卡的并行计算方式实现高效计算的优点而被广泛应用在可见光、X射线、太赫兹等波段的叠层衍射成像中[14,19,20].但是ePIE算法在太赫兹波段重建过程中有时会出现无法收敛或需要几百次迭代,导致收敛迟滞的情况发生.叠层重建算法是相位复原迭代算法的一种,其收敛速度和重建精度不仅受子孔径扫描过程中的线性交叠率影响,也与算法在不同平面的约束条件有关.探针光束照明样品形成物平面约束条件,而探测器采集的衍射图样可以看作记录平面强度约束条件,强度约束条件的复杂程度影响相位复原算法的重建速度及重建像质量[21,22].通常为了提高计算效率与成像精度,多波长照明[23]、多角度照明[24]及多物距记录[25,26]的方式都可以增加衍射图样数据的多样性,使得相位复原算法更快的收敛至最优解,从而更快的获得样品的复振幅分布情况.将多角度照明和多波长照明方式引入到基于现阶段太赫兹源和探测器的成像系统中较为困难,而多物距记录衍射图样强度的方式只需要改变探测器的记录距离,不改变光路结构,较为适用于连续太赫兹波叠层成像.

本文提出一种连续太赫兹波双物距叠层相衬成像方法,保证系统结构简单紧凑的情况下,在采集数据过程中沿照明光束传播方向移动热释电面阵探测器,对同一照明探针光束获得两组不同物距的衍射图样作为记录平面差异化约束条件,增加记录平面实际数据的冗余度与多样性,提高连续太赫兹波叠层成像重建算法收敛速度及重建像质量,并通过仿真和实验验证.

2 理论分析及算法流程

将双物距记录方法引入ePIE算法中,改进算法称为双物距记录叠层迭代引擎算法(dual-plane ePIE,dp-ePIE),改进算法原理如图1 所示.设相邻探针扫描位置间隔为 Rj,j=1,2,···,J 表示扫描中的探针位置.经过小孔调制平行光照明光束在距离d处形成衍射图样,称为探针 Pj(r),物面坐标向量为 r=(x,y).探针光束透过样品衍射传播至记录平面I和II,分别记录样品的衍射场的强度分布,其中衍射传播过程利用角谱衍射理论(angular spectrum propagation,ASP)计算.连续太赫兹波叠层双物距成像方法与传统叠层成像理论分辨率计算方法相同,最小可分辨尺寸受限于探测器可探测最大衍射角的空间频率成分.所以双物距叠层成像理论分辨率约为

其中λ为太赫兹激光器中心波长,θ为最大衍射角.衍射极限分辨率不超过波长.重建算法具体过程如下.

图1 双物距叠层成像的原理图Fig.1.Schematic diagram of dual plane ptychographic imaging.

步骤1假设样品的复振幅透过率函数为On,j(r),其中,n=1,2,···,N 表示算法迭代次数,j表示探针光束位置序号.当 j+1 时,对 应 的On,1(r)为猜测物体的复振幅透过率函数,因此,探针经过猜测物体后的出射光场为

步骤2通过角谱衍射理论计算方式,出射场ψn,j(r)传播距离 z1至探测器平面,得到记录面I的复振幅值为即

其中 F {·}和F-1{·} 分别表示傅里叶变换和傅里叶逆变换,qx和qy为空间频率.

步骤3利用探测器在记录面I采集到的衍射图的均方根代替步骤2中的幅值,得到新的复振幅分布

步骤4采用角谱衍射理论计算方法将更新得到的记录面的复振幅逆向回传到物面,得到新的物面光场分布通过两个更新函数分别更新猜测的物体和探针,物体的更新函数为

这里,α 是权重系数,取值一般在 [0.9,1.0]之间,实验中取值为0.98.探针的更新函数

这里,β 是权重系数,取值一般在 [0.9,1.0]之间,实验中取值为0.98.

步骤5采用更新后的物函数和探针函数的乘积作为通过物体后的出射场,扫描第 j+1 位置,更新步骤1中的复振幅 ψj+1,n(r),重复步骤1至步骤5,直至扫描完成第J位置,即扫描完成整个物面,得到完整的更新物函数.

步骤6将步骤5得到的更新物体复振幅函数,即利用记录位置I处衍射强度图样重建得到的样品复振幅分布,作为利用记录位置II处衍射强度图样重建样品信息的初始样品函数.将步骤1得到的出射光场传播距离 z2至探测器位置,记录该位置样品衍射场信息.重复步骤3至步骤5,得到新的样品复振幅透过率,即完成双物距算法的一次迭代.

步骤7重复步骤1至步骤6,完成N次重建迭代直至算法收敛,得到样品最终的幅值分布与相位分布.

由于在仿真实验中真实样品的分布是已知的,所以重建算法的收敛度可以直接由方均根误差(root mean square error,RMSE)评价函数衡量:

其中 On(r) 是第n次迭代后重建的样品复振幅分布.

3 仿真分析

首先通过模拟实际系统参数评估连续太赫兹波双物距叠层成像方法的有效性.设置太赫兹激光器中心波长 λ 为 118.83 μm,对应频率为 2.52 THz,模拟面阵探测器像元尺寸为 100 μm×100 μm,小孔直径为3 mm.连续太赫兹波双物距叠层成像仿真实验如图2所示,样品放置在小孔后 d=4 mm处,设置照明探针光束扫描交叠率为75%,分别在记录面I,II记录6×6幅衍射图样作为叠层衍射成像重建数据.记录面I距离样品20 mm,记录面II距离样品23 mm.

图2 仿真实验原理图Fig.2.Schematic diagram of the terahertz dual-plane ptychographic simulation.

通过仿真实验对比传统单物距叠层算法ePIE与双物距叠层重建算法dp-ePIE的重建效果.图3(a)和图3(b)分别作为原始图像的幅值与相位分布,尺寸为 512 pixel×512 pixel.幅值信息对应强度范围为0—1,相位信息在0—π之间分布.

单物距叠层成像实验将只利用记录平面位置I处的衍射图样对输入图像进行恢复,而dp-ePIE同时利用记录位置I,II的衍射信息重建样品复振幅分布.为了量化衡量dp-ePIE算法收敛度改进情况,将 ePIE算法及 dp-ePIE算法分别迭代500次,获得振幅及相位重建图像,并利用E0(n)对重建像中心 288 pixel×288 pixel大小的中心区域进行收敛度评价,该区域大小由交叠率、小孔尺寸及扫描步长决定,评价结果如图4所示.蓝色虚线为ePIE算法收敛度评价曲线,红色实线为dp-ePIE算法收敛曲线.

图3 仿真实验中的样品 (a) 振幅分布,(b) 相位分布Fig.3.Sample for ptychographic simulation:(a) Input amplitude distribution;(b) input phase distribution.

对比两算法收敛度曲线可以看出,在500次算法迭代过程中,dp-ePIE算法在相同迭代次数时收敛度一直低于ePIE算法收敛值.在第100次迭代时,ePIE 算法收敛度值 E0(100)=0.002,此时认为算法接近收敛,而dp-ePIE算法达到该值是迭代次数仅为52次.dp-ePIE算法在迭代400次时下降至 1×10–3,此时认为算法完全收敛,然而ePIE迭代至500次时仍没有完全收敛.收敛度评价曲线说明双物距叠层重建算法收敛速度明显高于ePIE算法,可以减少算法迭代次数,提高计算效率.

为了量化评价dp-ePIE算法重建振幅分布与相位分布质量,可利用相关系数作为评价依据[15]:

其中 m (x,y) 表示输入图像强度或相位分布的中心区域,n (x,y) 表示对应区域的重建图像,μm(x,y),μn(x,y)分别表示第m,n张图像的平均值.相关系数取值区间为0—1,相关值越高代表输入图像与重建图像相关性越强,图像重建质量越高.

图4 收敛度评价曲线Fig.4.Convergence of the ePIE and dp-ePIE algorithms.

由收敛度评价曲线可知,dp-ePIE算法在50次左右时接近收敛.为了对比dp-ePIE算法重建效果,分别对ePIE算法和dp-ePIE进行4次、20次、及50次迭代,获得重建幅值分布与相位分布如图5所示.CC表示重建图像中心位置与输入图像相关系数.算法迭代开始前,设置初始样品为512×512 个像元的随机矩阵.对比图5(a1)和图5(a2),dp-ePIE算法完成第四次迭代是可以获得相关系数CC=0.9418的振幅重建结果,重建图像质量明显优于利用ePIE算法恢复的相关系数为 CC=0.8366 振幅分布结果;对比图5(b1)和图5(b2)可以看出,ePIE算法重建相位分布模糊不清,相关系数为 CC=0.6393,而 dp-ePIE 算法同样可以获得重建效果更好的相位重建分布,相关系数为 CC=0.7368.对比图5(c1)和图5(c2)可以看出ePIE算法在第20次迭代时振幅重建图像相关系数为CC=0.9737,重建图像仍含有相位分布中的混叠信息,选定区域边缘较为模糊;而dpePIE算法在第20次迭代时振幅重建图像相关系数为 CC=0.9909,与原始输入图像几乎一致,选定区域边缘及中心十分清晰.对比图5(d1)和图5(d2)可以看出ePIE算法重建相位分布的相关系数为CC=0.7886,而 dp-ePIE 算法重建相位分布的相关系数为 CC=0.9493.相同的迭代次数下,dp-ePIE算法表现出更好的相位分布重建像.对比图5(e1)和图5(e2)可以看出,ePIE算法获得的幅值分布相关系数为 CC=0.9919,dp-ePIE 算法重建像相关系数为CC=0.9978,两算法均可获得相同质量的幅值分布重建像.对比图5(f1)和图5(f2)可以看出利用ePIE算法重建的相位分布相关系数为CC=0.8557,dp-ePIE 算法重建像相关系数为CC=0.8800,dp-ePIE 算法重建像相关系数高于ePIE算法重建结果,ePIE算法重建结果左下方有明显条纹,而dp-ePIE算法重建结果幅面清晰锐利.对比仿真结果可知,获得相同相关系数重建结果,dp-ePIE算法迭代次数明显少于ePIE算法迭代次数,说明dp-ePIE算法能够利用更少的时间获得叠层衍射成像重建图像.

图5 叠层仿真重建结果 (a1)、(b1),(c1)、(d1)和 (e1)、(f1)分别表示 ePIE 算法迭代 4 次、20 次及 50 次迭代重建振幅分布与相位分布;(a2)、(b2),(c2)、(d2)和(e2)、(f2)分别表示dp-ePIE算法迭代4次、20次及50次迭代重建振幅分布与相位分布Fig.5.Comparison of simulation results with ePIE and dp-ePIE:(a1) and (b1),(c1) and (d1),(e1) and (f1) are the amplitude and phase reconstruction with ePIE algorithm;(a2) and (b2),(c2) and (d2),(e2) and (f2) are the amplitude and phase reconstruction with dp-ePIE algorithm.

为探究dp-ePIE算法收敛速度改进情况,分别利用两算法在相同仿真参数情况下迭代200次,并以图像相关度作为评价指标,可得到dp-ePIE算法及ePIE算法收敛速度对比,结果如图6所示.图6(a)中红色实线表示dp-ePIE算法在200次迭代过程中收敛情况,黑色虚线表示传统ePIE算法收敛情况.曲线斜率显示出dp-ePIE算法收敛速度较ePIE算法有明显提升,同时最终收敛至趋于1,明显高于ePIE算法收敛值.图6(b)表示两算法在200次迭代过程中恢复相位信息时的收敛情况.曲线斜率显示出dp-ePIE算法收敛速度相较于ePIE算法有大幅度提升,同时重建相位像与输入图像相关系数明显高于ePIE算法重建结果,说明连续太赫兹波双物距叠层成像方法不仅能够提升算法收敛度,而且能够改善重建像质量,其中dp-ePIE算法针对相位信息改善效果更加明显.

图6 重建结果的相关系数比较 (a) dp-ePIE 与 ePIE 算法幅值重建结果与迭代次数关系;(b) dp-ePIE与ePIE算法相位重建结果与迭代次数关系Fig.6.Comparison of correlation coefficients of reconstruction results:(a) Comparison of correlation coefficients of amplitude reconstruction results with dp-ePIE and ePIE;(b) comparison of correlation coefficients of phase reconstruction results with dp-ePIE and ePIE.

双物距记录叠层成像方法的一个重要参数为两次采集衍射强度图信息的间隔距离,即探测器放置间隔距离.为了确定最理想的探测器位置间隔距离,以记录面I为参考位置,记录间隔取值范围为1—15 mm,记录间隔 1 mm,分别设置记录 II,选取第200次迭代时相关系数作为评价指标.相关系数与间隔距离关系如图7所示.红色实线表示幅值重建像相关系数曲线,蓝色虚线表示相位重建像相关数曲线.幅值重建像相关系数在不同间隔距离下均接近1,振幅重建像相关系数在0.9附近,说明间隔距离的选取对重建像质量及算法收敛速度无明显影响.

4 实验与分析

4.1 实验装置

为了进一步实验验证该方法的有效性,搭建了连续太赫兹波双物距叠层相衬成像系统,成像装置如图8所示.实验利用FIR295型连续太赫兹波激光器作为辐射源,太赫兹波中心波长为118.83 μm,中心频率为 2.52 THz,最大输出功率为 500 mW.光束经由两个表面镀金的离轴抛物面镜进行准直扩束.样品固定在三维平移台上(MT3-Z8,精度为± 0.1 μm),并放置在直径为 3.3 mm 的小孔后的4 mm处.实验利用热释电面阵探测器(Pyrocam-III,Spiricon)采集衍射图样,探测器像素尺寸 Δ=100 μm,像素个数 N=124 pixel×124 pixel,斩波频率为48 Hz,并放置在样品后20 mm位置处并设定该位置为第一个探测器记录平面,23 mm设置为第二个探测器记录平面.沿每一记录平面的x,y方向分别扫描10次,即记录100幅衍射图样.实验过程中,扫描步长为 0.8 mm,交叠率为78%.视场范围为 10.025 mm×10.025 mm.

图8 连续太赫兹波叠层成像实验装置示意图Fig.8.Setup of continuous-wave terahertz dual-plane ptychography.

实验样品选择为聚丙烯“可回收利用”环保三角标志.这种材料是由丙烯聚合而成的一种热塑性树脂材料,对于太赫兹波的吸收较少,可以近似为纯相位样品.标志刻蚀在500 μm厚的聚丙烯基底上,图案模型如图9(a)所示.数字“5”表示聚丙烯材料,外围标志刻蚀宽度为 190 μm,内部数字刻蚀宽度范围为 190—220 μm,刻蚀深度 50 μm,精度为 0.1 μm,图9(b)为样品实物图.

图9 聚丙烯可回收三角标志图案样品 (a)三角标志模型;(b)样品实物图Fig.9.Sample of recyclable polypropylene triangle pattern:(a) Model of the sample;(b) picture of the sample.

多波长照明、多角度照明和双物距记录都是通过提升测量数据冗余度、增加采集数据多样性而提升系统重建像质量的方法.双物距较适合于基于现阶段太赫兹源和探测器搭设的叠层成像光路.一方面,FIR 295型光泵可调谐太赫兹激光器体积庞大,操作复杂,需要调节二氧化碳激光器泵浦源谱线才能对太赫兹激光器波长进行调节,难以做到波长快速切换,实现多波长照明方法较为困难.另一方面,考虑到太赫兹波传输损耗,要求实验系统结构简单紧凑.多角度照明方法如果采用分立元件,复杂的光路会增加太赫兹波的传输损耗,不可见光也会增加调节难度.连续太赫兹波采用光纤耦合现阶段也比较困难.与上述方法相比,双物距记录方法不需要调整太赫兹激光器输出参数,也无需额外的光学元件调制太赫兹照明光束方向,不会增加光路系统复杂度及太赫兹波的传输损耗.

4.2 重建结果及其分析

连续太赫兹波双物距叠层成像可以认为是远场衍射成像,理论分辨率受限于最大可探测衍射角,即由样品至探测器记录距离决定,所以叠层成像应将探测器尽量靠近样品放置.探测器在不同记录位置处衍射图样包含不同的衍射信息,不是几何缩放关系.通过自动聚焦算法精确计算样品到记录位置I距离为20.2 mm,记录位置II距离为23 mm.根据(1)式计算该系统理论分辨率约为182 μm.单物距叠层重建时应只利用记录位置I处的衍射图样.

分别利用ePIE、dp-ePIE算法对高斯预处理后的衍射图进行处理,实验过程中采集的数据量较大,双物距数据采集过程尽管牺牲了数据的采集时间,但增加了样品相同部位的衍射信息.两算法重建结果对比如图10所示.为了量化衡量dp-ePIE算法与ePIE算法重建图像质量,利用无参考结构清晰度评价方法 (no-reference structural sharpness,NRSS)比较两算法在相同迭代次数时重建的振幅分布与相位分布的清晰程度.无参考结构清晰度评价方法的理论在文献[27]中进行了详细推导,NRSS评价系数在0—1之间分布,NRSS值越大,说明重建图像越清晰,边缘更锐利.

图10 两种叠层重建算法分别迭代10次后可回收标志的重建结果 (a1),(b1)分别表示ePIE算法重建振幅分布及相位分布;(a2),(b2)表示dp-ePIE算法重建振幅分布及相位分布Fig.10.The reconstructed results after 10 iterations by two different reconstruction algorithms:(a1),(b1) Represent the amplitude and phase reconstructed based on ePIE algorithm:(a2),(b2) represent the amplitude and phase reconstructed based dp-ePIE algorithm.

对比图10(a1)与图10(a2)可以看出在10次迭代下,dp-ePIE算法重建振幅分布NRSS值为0.7879,ePIE算法重建振幅分布NRSS值为0.8211,dp-ePIE算法重建振幅分布轮廓较ePIE算法重建振幅分布更加清晰.对比图10(b1)与图10(b2)可以看出第10次迭代时ePIE算法无法重建样品基本轮廓信息,也无法对样品进行分辨,NRSS值为0.8233;而dp-ePIE在第10次迭代时可以获得样品清晰的相位重建结果,能够清晰分辨在聚丙烯材料上的三角形“回收标志”刻蚀边界及分类数字“5”,NRSS 值为 0.9831.对比重建相位分布可以看出dp-ePIE算法收敛速度提升明显,图像保真度大幅改善.NRSS评价系数反映的结论与仿真实验一致,即dp-ePIE算法能够在少量迭代次数时完成对待测样品的重建工作,可以快速获得相位型样品的高质量相位分布结果.

计算效率方面,运行ePIE和dp-ePIE算法所用计算机性能相同,内存 (RAM)大小 8 GB,处理器为英特尔酷睿 i5-6400 (四核,2.7 GHz).对实验收集的聚丙烯三角形“可回收”标志数据进行重建,dp-ePIE算法单次迭代时间为5.380 s,而ePIE算法单次迭代时间3.469 s.dp-ePIE算法单次迭代时间更长,原因在于dp-ePIE单次重建迭代过程需要在物平面与两个记录面之间进行衍射计算.但是ePIE算法迭代20次,即完成20次物平面与记录面之间的更新计算的总耗时为57.117 s,相位重建图像NRSS值为0.9519,而dp-ePIE算法迭代10次总耗时为 45.086 s,相位重建结果 NRSS值为0.9831.综上所述,dp-ePIE算法能够有效改善重建像质量,降低算法总耗时,提升连续太赫兹波叠层相衬成像计算效率.

实验结果证明了基于双物距的连续太赫兹波叠层成像方法利用差异化、多样性衍射信息为记录平面的约束条件,可以大幅度改善算法收敛速度,解决了ePIE算法在太赫兹波段收敛迟滞的问题,同时提升了相位重建像质量.两种算法的振幅重建结果中出现了网格状伪影,该伪影产生原因是扫描重建所致,与引入物平面约束条件及记录平面约束条件无关.

5 结 论

针对连续太赫兹波叠层成像重建过程中ePIE算法收敛迟滞问题,本文提出双物距叠层成像方法及dp-ePIE算法.在传统单物距基础上,对同一照明光束位置,沿光轴不同位置采集两次的衍射图样作为记录平面约束条件,利用差异化衍射图样重建样品复振幅透过率函数,提升记录平面衍射信息冗余度及复杂度.相较于多波长照明和多角度照明提升测量数据冗余度、增加采集数据多样性的方法,双物距记录方法不需要调整太赫兹激光器输出参数,也无需额外的光学元件调制太赫兹照明光束,因此双物距叠层成像方法不会增加光路系统复杂程度,能够保证实验元件紧凑放置以降低太赫兹波在传输过程中的功率损耗.计算机数值仿真和搭设的连续太赫兹波双物距叠层实验系统证明将双物距记录衍射图样作为记录面约束条件的方法可有效减少重建算法迭代次数,降低了连续太赫兹波叠层成像重建过程总耗时,加快算法收敛,有效解决了传统叠层成像重建算法收敛迟滞问题.