金 川，何 渝，唐 燕，刘俊伯，孙海峰，胡 松*

1中国科学院大学，北京 100049；2中国科学院光电技术研究所，四川 成都 610209

1 引言

近年来，数字全息显微作为一种定量相位成像技术，以其宽视场、非接触、免标记等优点，日益成为生物样品探测、微纳结构检测、流体场分析等领域有力的测量工具[1-5]。根据物波与参考波的夹角，数字全息显微系统一般可分为三种不同的干涉结构：同轴、离轴、微离轴[6-8]。作为同轴和离轴的中间结构，微离轴结构产生的数字全息图的直流像分别与实像和孪生像互相重叠，但实像和孪生像之间在空间频域上是互不重叠的。与同轴结构相比，微离轴数字全息因相移步数的减少，提供了更快的图像采集效率。与此同时，由于载频的降低，微离轴结构比离轴结构能更好地利用图像传感器的空间带宽。但是，传统的微离轴数字全息系统需要连续采集两幅全息图，不适用于样品的动态测量。为了提高微离轴数字全息系统的实时性，多种基于单次曝光的同步相移方法被提出。其中，最常用的同步相移方案通过视场复用技术实现，该技术使用不同的分光元件，例如：光栅、分束器、沃拉斯顿棱镜等，在单传感器的不同位置记录多幅具有相移差的全息图[9-11]。然而，该技术由于获取的全息图空间相对位置是未知的，存在空间失配问题，这会影响后续全息重建的精度。

为了解决该问题，Kiire 等人[12]通过移除全息图的相移使条纹可见度一致，然后通过强度差完成标定。Millerd 等人[13]通过计算同一目标的质心位置以实现标定。Hahn 等人[14]使用遗传算法优化基于强度信息的评价函数。以上这些标定方法都依赖于图像的强度信息，但实验过程中会产生的散斑噪声、寄生条纹或光束强度不一致等一系列成像误差，严重影响标定的准确性。基于相位的标定方法则可以避免上述问题，Li 等人[15]引入了圆形载波来获取标定参数，但是这种方法不能用于一些无法引入圆形载波的测量场景。Zheng 等人[16]报道了一种具有像素级精度的部分相位相关标定方法，标定过程中图像之间逐像素移动以计算空间位置。但是，以上这几种技术都没有考虑传感器平面倾斜导致的图像纵向位置误差。

本文提出了一种新颖的同步微离轴数字全息空间失配标定方法。整个标定过程被描述为非线性多变量问题，精准的相对位置误差在非线性优化过程中获得。标定的建模过程构建了纯相位波前，其相位是通过对提取系统相位畸变进行Zernike 拟合而来，以此排除全息图强度和相位噪声对标定精度的影响。与传统的标定方法[12-16]相比，我们提出的标定方法具有亚像素级的精度，并且易于实施。除此之外，相比于其他标定方法更多地考虑了由传感器平面倾斜引起的纵向位置误差。本文组织如下：第2 节描述了空间失配的原理和标定过程，第3 节展示了仿真效果，第4 节提供了实验结果。

2 理论推导

2.1 同步微离轴数字全息系统

基于偏振分光棱镜[17]的同步相移微离轴数字全息系统的光路设置如图1 所示。45°线偏振光束被分束器分成两束光。由显微镜、针孔和透镜组成的准直扩束器分别对两束光束进行扩束准直。其中一束光透射过1/4 波片，将线偏振光转换为圆偏振光作为参考光，中性密度衰减器用于调整光强比。另一束光作为物光透射通过样品，被10×显微物镜和管透镜组成的显微成像系统收集。物光和参考光由分光棱镜2 以微小角度重新组合，偏振分束棱镜用于将入射光束分成垂直S 偏振分量和水平P偏振分量。最后，两幅存在π/2 相移差的全息图由单个CCD (分辨率1544× 2064，像素大小3.45 μm)通过一次曝光同时捕获。

图1 同步微离轴数字全息显微系统光路图Fig.1 Optical setup of the parallel phase-shifting SODHM system

成像系统采集到的两幅相移量为 π/2的全息图可以表示为

全息图的直流像可以通过I1减去I2来移除：

通过在傅里叶频谱中执行滤波运算，可以得到物光复振幅信息:

其中：FT 和IFT 分别表示傅里叶变换和逆傅里叶变换，RD代表用于相位畸变校正的数字参考波，频谱窗W被用于执行滤波操作。

2.2 空间失配问题分析

从上述原理中可知，提取物光复振幅的关键是确保两幅全息图处于相对空间位置一致。然而，由于全息图空间失配问题的存在，上述步骤往往无法精确实现。根据2.1 小节分析，通过CCD 的单次曝光，两幅相移差为 π/2的全息图被同时采集，如图2(a)所示。在提取物光复振幅信息之前，需要对该全息图进行图像预分割，红色虚线的方框内得到两幅分离的全息图，其全息图相对的坐标系分别为(x1,y1)和(x2,y2)。因此，两幅全息图存在横向位置误差(Δx,Δy)=(x1-x2,y1-y2)。除此之外，很难完全保证光束垂直于CCD 平面入射，CCD 平面的轻微倾斜会导致纵向位置误差 Δz，如图2(b)所示。

图2 空间相对误差的生成原理。(a) 横向误差；(b) 纵向误差Fig.2 The generation principle of spatial relative position errors.(a) The longitudinal position error;(b) The transverse position error

2.3 空间误匹配标定方法

假设两幅全息图的光轴相互重叠，则全息图平面的位置关系如图3 所示。在笛卡尔坐标系中，位于参考平面H2的第二幅全息图与Z轴垂直，位于坐标系(x2,y2)。在坐标系(x1,y1)下，H1平面中的第一幅全息图相对于H2平面的第二幅全息图有Δx,Δy,Δz的平移量。全息图的空间相对位置误差导致两幅全息图中的光场信息不一致，如果已知位置误差参数，则可以利用基于角谱法的空间平面间的波场传播数值计算，将波场h1(x1,y1)从平面H1衍射到平面H2以消除差异：

图3 全息图平面的位置关系Fig.3 The position relationship of holograms

其中：F为整体的空间变换函数，u和v是x和y方向的空间频率。整个过程包括两个步骤：首先，执行从H1平面到M平面的平移过程，以消除横向位置误差。相对于像素点在空间域的离散特性造成的不连续平移，该过程转化为空间频域频谱的移相操作来实现，以完成高精度的亚像素级平移。第二步是纵向位置误差的补偿，利用数值计算将M平面衍射至H2平面，以消除纵向位置误差引起的物光场的相位信息变化，其中该衍射过程是通过角谱法实现的[18]。

在实际测量过程中，直接采集到的两幅全息图会存在散斑噪声、寄生条纹或光强不稳定等问题，若以图像的强度信息作为标定标准，则这些干扰会严重影响标定精度。为了避免强度噪声的影响，我们利用补偿后全息图和参考全息图中的相位畸变差异为标定标准来计算相对位置误差参数(Δx,Δy,Δz)。首先，调整系统反光镜的角度以增加载频，将微离轴数字全息系统转换为离轴全息结构，并且记录无样品的两幅全息图，可表示为

其中：Rn为参考光波，On为 物光波。在傅里叶频域进行滤波便能提取物光复振幅，可以表示为

其中：An(xn,yn)是幅度分布，相位分布 φn(xn,yn)由相位噪声φnn(xn,yn)和相位畸变φan(xn,yn)两部分组成。基于相位畸变的波前可以通过Zernike 多项式拟合的操作来构建，以避免相位噪声影响：

式中：Zi(xn,yn)是笛卡尔形式下的Zernike 多项式，是多项式的对应系数，相位畸变可以用1～9 低阶Zernike 多项式来描述。将波前 ψ1(x1,y1)代入等式(4)得到补偿后的波前：

综上所述，同步微离轴数字全息显微系统空间一致性标定共分为四个步骤，如图4 所示。

图4 标定算法的流程Fig.4 The flow of proposed calibration method

第一步：增加系统载频，将微离轴光路调整为离轴光路，获取无样本离轴数字全息图；

第二步：采用频谱滤波的方法提取物光波，构建基于Zernike 拟合的纯相位波前；

第三步：将标定问题转化为式(8)中的非线性多变量优化问题，其评价函数被设置为补偿后波前相位的均方根误差；

第四步：求解优化问题以得到相对位置误差。

2.4 粒子群算法原理

作为一种经典的启发式算法，粒子群算法(PSO)是根据鸟群的捕食行为提出的，其中鸟群被抽象为无质量和体积的粒子，每个粒子的运动方向和速度会受到个体和种群历史最优解的影响，以完成整个优化过程。相比于常用的基于梯度的算法[19-21]，该算法具有结构简单、收敛效率高和全局搜索能力强等特点。由于其优良的优化性能，粒子群算法被广泛应用于多个应用场景，如作业车间调度[22]、人工神经网络[23]、多峰函数优化[24]等。因此，线性惯性权重递减策略的粒子群优化算法(LDW-PSO)[25]被用于求解标定过程中的优化问题。

粒子群优化算法的整体流程如图5 所示，在标准粒子群算法中，会先初始化一组随机粒子，包括粒子的位置m和速度v两个变量。在标定模型中，粒子的位置m为相对位置误差，它决定了评价函数的大小。通常来说，评价函数的值称为适应度，适应度越小表面标定精度越高。粒子速度v控制着粒子每次迭代的步长和迭代方向，通过观察适应度最小的两个最优粒子进行迭代优化。这两个最优值包括多个迭代过程中粒子自身位置的个体最优值Pbest和整个种群位置的全局最优值Gbest。通过计算得到的这两个最优值，对粒子的两个速度和位置变量进行更新，其更新的表达式为

图5 LDW-PSO 算法的流程Fig.5 The whole flow of LDW-PSO algorithm

其中：t为当前迭代步数，c1和c2分别表示自我认知因素和社会认知因素。r1和r2两者都表示随机数，r1,r2∈(0,1)，w是惯性权重系数，它决定了粒子当前的适应度对下一个搜索方向和步长的影响。本算法采用惯性权线性递减策略(LDW)，该策略使惯性权重随迭代次数的增加从最大值wmax递减到最小值wmin，权重w表达式为

式中：tmax代 表最大的迭代次数。

一旦获得位置误差参数，我们再减小系统载频，将其恢复到微离轴的配置，并收集带有样品的全息图Ins(xn,yn)。全息图I1s(x1,y1) 从H1平面传播到H2平面，利用等式得到校正后的全息图I′1s(x2,y2)。最后，通过式(3)中传统微离轴的重建方法，可以计算出物光波的复振幅信息。

3 仿真分析

为了证明所提方法的有效性和高精度，我们进行了数值模拟。系统的相位畸变用ZEMAX 分类中的1～9 低阶Zernike 多项式来描述，如表1 所示的各阶Zernike 多项式以及其拟合系数设置，最后产生的相位畸变如图6(a)所示。模拟中所使用相位图的分辨率为800 pixels×600 pixels，像素大小为3.45 μm。图6(b)中两幅无样品的相移离轴全息图被生成用于标定，两幅全息图之间引入了相对位置误差，其具体参数位于表2 的第二行。与此同时，这些全息图都被添加了最大为10%的正态分布的随机强度和相位噪声。

表2 相对位置误差的数据比较Table 2 Data comparison of relative position error

图6 仿真输入。(a) 仿真的相位畸变；(b) 无样本的仿真π/2 相移量离轴全息图Fig.6 Simulation input.(a) Simulated phase aberration;(b) Simulated π/2 phase shifting off-axis holograms without samples

表1 仿真中ZEMAX 类型的Zernike 多项式及其系数Table 1 Zernike polynomials of ZEMAX type and their coefficients in the simulation

部分相位相关(partial phase correlation，PPC)标定方法[16]是一种像素级精度的空间位置配准方法，将图像逐像素点移动来寻找相位差最小点的位置，该方法已经被广泛应用于干涉图的空间一致性标定问题中。同时采用部分相位相关标定方法和本文标定方法在上述设置的条件下进行模拟仿真，分别获取图像间的相对位置误差。PPC 标定方法得到的误差参数如表2 第三行所示。由于PPC 方法在原理上的局限性，仅得到了具有像素精度的横向位置误差，并且忽略了纵向位置误差。在本文标定方法的仿真过程中，为保证优化速度和精度，粒子维度为3，粒子群规模被设置为30。自我认知因子和社会认知因子均设为2.05，线性惯性权重因子的最大值和最小值分别为0.9 和0.4。考虑图像的粗对准精度和粒子群算法的收敛效率，严格限制横向位置误差的位置和速度范围为(—50 pixels,50 pixels)，设置纵向位置误差的位置和速度范围为(—3 mm,3 mm)。LDW-PSO 算法经300 次迭代之后，本文方法的标定结果如表2 第四行所示。真实引入和标定的横向位置误差之差小于0.01 pixel，纵向位置误差之差为0.002 mm。毫无疑问，本文方法的标定数据与引入的位置误差几乎一致。

另外，基于LDW-PSO 算法的优化曲线如图7(a)所示，300 次迭代后相位分布的RMSE 从0.265 rad降低到1.798×10-3rad。采用文本方法标定后，两个波前的整体相位差分布如图7(b)所示，表明了两个波前的高度一致性。值得注意的是，相位差分布边缘出现波浪形不规则条纹，这是因为仿真中的频谱滤波操作会产生振铃伪影。即使在这种情况下，这种标定仍然具有极高的准确性。然而，采用PPC 方法得到的相位分布的RMSE 为0.0407 rad，图7(c)给出了两个波前的整体相位差分布。仿真结果表明，本文方法相比于PPC 方法具有更高的精度，甚至达到亚像素级。

图7 仿真结果。(a) 基于LDW-PSO 的相位分布RMSE 收敛曲线；(b) 使用本文方法标定后的整个相位差分布；(c) 使用PPC 方法标定后的整个相位差分布Fig.7 Simulation results.(a) RMSE convergence curve of phase distribution based on LDW-PSO;(b) The whole phase difference distribution after calibration using the method proposed in this paper;(c) The whole phase difference distribution after calibration using PPC method

4 实验结果

首先，对图1 所示的系统做空间一致性标定。为了提取系统的相位畸变，调整反射镜2 将系统转变为离轴结构。无样本离轴全息图被记录在图8(a)中，其被红色矩形粗略地分割为两幅全息图。由于系统的成像特性，两幅全息图沿图像中心暗区对称分布，因此需要对其中一幅全息图进行翻转操作。将图像左侧的全息图设置为参考全息图，并使用本文标定方法。在标定过程中，基于LDW-PSO 的目的函数的收敛曲线如图10(a)所示。经过300 次迭代，相位分布的RMSE 下降至 9.831×10-3rad，相位差分布如图10(b)所示。与仿真结果相比，由于相位噪声产生的拟合误差，因此波前差异较大，但不可否认，两个波前仍然具有良好的一致性。最终标定结果如表3 所示。

表3 相对位置误差的实验数据Table 3 Experimental data of relative position error

图8 实验系统采集的全息图 。(a) 无样品的离轴全息图；(b) 带样品的微离轴全息图Fig.8 Holograms collected by the experimental system.(a) The off-axis hologram without sample;(b) The slightly off-axis holograms with sample

图9 轮廓仪测量的样品轮廓Fig.9 Profile of the sample measured by the profilometer

图10 实验结果。(a) 基于LDW-PSO 的相位分布RMSE 收敛曲线；(b) 校准后的整个相位差分布；(c) 重建的样本相位分布；(d) 沿着图10(c)中标记的白线的轮廓Fig.10 Experimental results.(a) The RMSE convergence curve of phase distribution based on LDW-PSO;(b) The entire phase difference distribution after calibration;(c) The phase distribution of reconstructed sample;(d) The outline of the white line marked in Figure 10(c)

获得相对位置误差参数后，恢复微离轴光路结构，选择由熔融石英制成的微透镜阵列(折射率为1.458，平均高度为2.042 μm)作为样品，其轮廓经探针式轮廓仪测试，如图9 所示。样品的微离轴全息图如图8(b) 所示，全息图被相同的红色矩形分割成两幅，并利用所提标定方法进行校正。通过执行微离轴全息图的物波重建算法，微透镜的相位得到了很好的恢复，如图10(c)所示。图10(d)为微透镜白色虚线处的轮廓，测得的平均高度为2.025 μm，相应的轮廓结果偏差仅为17 nm。上述结果证明了所提出的同步微离轴数字全息的空间一致性标定方法的有效性。

5 结论

为了实现同步微离轴数字全息显微系统中的精准全息重建，本文提出一种基于相位差最小化的空间失配标定方法，该标定方法适用于同步微离轴数字全息显微系统。在标定过程中，利用系统的相位畸变构造基于Zernike 拟合的纯相位波前，以有效抑制了图像中影响标定精度的强度噪声和相位噪声。更重要的是，本文方法考虑了传感器平面倾斜引起的纵向位置误差。通过最小化构造波前相位分布的均方根误差来提取相对位置误差参数。仿真和实验结果证明了所提出的校准方法的亚像素级精度和可用性。在未来的工作中，可以进一步考虑将本文方法扩展到其他无需成像的干涉仪系统中。