结合线性回归的离轴数字全息去载波相位恢复算法*

单明广 刘翔宇 庞成 钟志 于蕾 刘彬 刘磊†

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

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

为实现仅用一幅离轴数字全息图便能直接恢复相位，提出一种利用空间载波相移技术(spatial carrier phase shift，SCPS)和线性回归相结合的离轴数字全息去载波相位恢复算法.首先，利用SCPS 将一幅离轴数字全息图分为四幅含有载波相移的全息图，其中载波相移由沿行、列两个方向的正交载波所引入;然后，将四幅载波相移全息图作为输入，将所求物体相位和两个正交的载波作为未知量，结合最小二乘法和线性回归同时求出载波和相位信息.相较于已有的去载波技术，本算法无需背景全息图作为参考，便可准确地去除载波，实现高质量的相位重建.本文结合数值仿真和具体实验结果验证本算法的有效性和优越性.

1 引言

数字全息技术[1-6]因其非接触、分辨率高等优点，在光学测量领域得到了广泛的应用.按照恢复算法，可将数字全息分为同轴数字全息[2，3]和离轴数字全息[4-6]两类.其中，同轴数字全息可充分利用相机的空间带宽利用率，具有较高的分辨率，但是通常需要牺牲视场利用率[3，7]或时间利用率[8，9]采集两幅及以上的全息图，并且需要复杂的相移设备引入相移.而离轴数字全息技术虽然对相机的空间带宽利用率不高，但是仅需一幅全息图便可恢复物体信息[10-12]，提高了系统的测量效率，也降低了环境对多帧测量的影响.

傅里叶变换法(Fourier transform algorithm，FTA)[10]是离轴数字全息最常用的相位恢复算法.该方法简单快速，但恢复过程中使用的带通滤波会造成频谱泄露，同时带通滤波器截断频谱的操作会在恢复结果引入吉布斯效应，从而影响到恢复质量;同时，FTA 需精确确定实像频谱信息强度峰值的位置才能有效地去除载波，但由于全息图的数字化特点，实像频谱信息强度峰值很难恰好位于整数像素位置，从而造成载波残留等影响恢复质量的问题.为去除载波的影响，目前最常用也是最有效的方法是双曝光法(double exposure algorithm，DEA)[10，13]，即在实验之前采集一幅不含样品的全息图，以此记录系统的详细载波信息，通过简单的相减或相除便可准确地去除载波，但是DEA 要求系统具有极高的稳定性，保证系统参数在测量过程中保持不变;同时，某些反射式全息系统的无样品全息图很难采集，从而加大了实际使用的难度.为了实现亚像素级的载波估计，Fan 等[14]提出了质心法(spectrum centroid method，SCM)，通过计算峰值所在部分区域的质心估算载波，但是SCM 的估算精度严重依赖所选区域，其估算精度在实际应用中严重受限.Du 等[15]和Hincapié-Zuluaga 等[16]提出了补零法(zero padding，ZP)，通过对全息图补零来提高载波估算精度，但是估算精度与计算量成正比，计算量巨大;同时，还需反复尝试才能找到合适的补零倍数，使得该算法在实际应用中难以取得良好效果.近期，基于全息图翻转变换[17-19]的去载波算法得到了极大的发展，但是该类算法对物体的形状和在全息图中的分布位置有着严格的要求，并不具有普适性，制约了该类算法的适用范围.

为了避免带通滤波操作对相位恢复造成影响，空间载波相移技术(spatial carrier phase shift，SCPS)被应用到离轴全息的相位恢复领域.SCPS 是指将一幅离轴全息图在空域沿一个方向进行像素移位，获取两幅以上全息图，并将提取出的全息图按照相移全息图对待，使用同轴数字全息的恢复算法计算出含有载波和物体相位的相位信息，再去除其中载波信息，最终得到准确的物体相位信息.Guo 等[20]在SCPS 的基础上，提出了采用最小二乘的相位恢复算法，该算法首先选取局部区域并估算其所含有的载频信息，再通过最小二乘法计算物体相位，但是该算法的恢复质量与局部区域的位置选择和尺寸密切相关，而如何准确地选择合适的区域限制了该算法使用.Stykyu 和Patorski[21]指出该方法并没有解决背景和调制度造成的误差，提取的相位是不稳定的.Xu 等[22]针对Guo 等提出算法的缺点，提出了基于最小二乘的空间载波相移算法，提取了高精度的相位，并且提高了算法的稳定性.Huang等[23]提出在频谱中初步估算峰值信息以减少迭代次数.Liu 等[24]基于SCPS，利用利萨如图形和椭圆拟合实现了相位恢复，提升了算法的抗噪能力和恢复质量.由此可知，相较于相移法，SCPS 仅需一幅全息图便可完成相位恢复，并可避免吉布斯效应和频谱泄露等问题.但是上述基于SCPS 的相位恢复算法仍需要准确的载波信息才能恢复出物体相位信息.综上，FTA 和SCPS 在相位恢复时，均面临着如何准确、有效地去除载波的问题.

因此，基于SCPS，本文提出一种结合线性回归的离轴数字全息去载波相位恢复算法，实现从一幅离轴全息图同时直接恢复出载波与物体相位.本算法将物体相位、载波信息和全息图的直流量与调制量均作为未知量，结合线性回归等技术，实现载波和物体相位的同时直接恢复.本文将给出本算法的原理及计算过程，通过仿真比较几种相位恢复算法的性能，并使用实验数据验证本算法的有效性和优越性.

2 理论方法

一般情况下，离轴数字全息图可以用下式表示:

其中A(x，y)，B(x，y)分别表示全息图的直流量和调制量;φ(x，y)代表待测物体相位信息;kx，ky分别代表全息图中沿x方向、y方向的载波频率.

在实际全息图中，相较于载波，直流量A、调制量B和相位信息φ的变化非常缓慢.将I(x，y)在空域上进行移位，即下移、右移以及右下移各一个像素，从而得到四幅相移全息图:

其中由于直流量、调制量和相位变化缓慢，A2，A3，A4和A近似相等，B2，B3，B4和B近似相等，φ2，φ3，φ4和φ近似相等.因此，经过空间移位，从一幅离轴全息图提取出四幅含有载波相移的全息图，其中相移值分别为0，kx，ky，kx+ky.将(2)式整理为一般表达式，如下:

式中，下标m=1，2，3，4 代表图像索引，上标t代表理论值，δxm和δym代表沿x，y两个方向的相移量.设沿行、列方向的载波相移分别为Δym(y)，Δxm(x)，可以得到:

设坐标(x，y)的总载波相移为Δm(x，y)，为了进一步阐述载波相移和沿x，y两个方向的载波和相移之间的关系，Δm(x，y)可由下式表达:

在完成上述基本分析之后，接下来将详细地描述算法流程，该算法的每一次迭代都包含三个步骤:

步骤1利用行、列的载波相移计算相位φ

定义a(x，y)A(x，y)，b(x，y)B(x，y)cos[φ(x，y)]，c(x，y)-B(x，y)sin[φ(x，y)]，则(3)式可改写为

步骤1 中，假设kxm，kym和相移δxm，δym是已知的，因此Δm(x，y)是已知的，而a(x，y)，b(x，y)，c(x，y)是所要求的量，根据最小二乘法[25]，可得

其中M代表干涉图的数量;Δm代表第m张干涉图的载波相移;Im代表第m张全息图的实际强度.

则物体相位信息可通过下式计算:

步骤2用物体相位和行载波相移确定列载波相移，进而使用线性回归算法计算列载波kx

由(3)式和(4)式可知，列载波相移Δxm(x)、行载波相移Δym(y)分别代表着x，y方向的载波相移量，二者相互独立，且具有各自对应的线性函数.步骤2 将利用相位φ(x，y)和行载波相移Δym(y)计算列载波相移Δxm(x).

式中，Y代表总图像矩阵第x列的长度，即图像矩阵的行高.

则第x列载波相移Δxm(x)可以通过下式计算:

在离轴全息图中，载波频率通常都较大，上式通过反正切求取Δxm(x)必然出现包裹.又因为Δxm(x)是关于x的线性函数，因此，可以通过线性回归[26]得到载波频率kxm，即

最后，对所求得的kxm求取均值，得到列载波频率kx.

步骤3用物体相位和列载波相移确定行载波相移，进而使用线性回归算法计算行载波ky

X代表总图像矩阵第y行的长度，即图像矩阵的列宽.第y行载波相移Δym(y)可以通过下式决定:

与步骤2 同理，对求得的Δym(y)通过线性回归得到载波频率kym，最后对kym求取均值，得到列载波频率ky.

上述三个步骤构成了一次迭代过程.在迭代次数达到预设值或是恢复结果达到预设精度，则计算结束.本文中，将迭代过程中前后两次求得的kx和ky的误差值作为迭代是否停止的标准，计算过程如下:

其中阈值δ是预设值;i代表迭代次数.在本文仿真及实验中δ被设置为10—6.

为了更清楚地描述恢复过程，每次迭代的流程如图1 所示.

图1 算法的基本流程Fig.1.Flowchart of the proposed algorithm.

最后，对本算法进行总结，算法流程如下:

1) 使用SCPS 技术，从一幅离轴全息图中提取四幅载波相移全息图;

2) 将两个方向的载波相移作为已知量，通过最小二乘法求出物体相位;

3) 将物体相位、行载波相移作为已知量，通过最小二乘法和线性回归求出列载波相移和列载波;

4) 将物体相位、列载波相移作为已知量，通过最小二乘法和线性回归求出行载波相移和行载波;

5) 重复步骤(2)—(4)直至满足迭代停止条件.

当满足迭代条件计算结束时，便可同时得到载波信息和物体相位信息.由上述算法描述可知，不同于目前已有的最小二乘相位恢复算法，本算法利用线性回归技术建立了载波相移与载波之间的联系，实现了待测物体相位信息和载波信息的同时直接获取.

3 仿真实验

为了验证本文所提算法的有效性和优越性，首先利用仿真对本算法进行验证，同时与FTA，SCM，ZP 和DEA 进行比较，仿真过程中使用波长为632.8 nm 的光源.首先，仿真生成一个光程差为100 nm 的相位型半球，像素大小为256×256，像素尺寸为4.4 μm;其次，引入沿x方向载波kx=0.76，沿y方向载波ky=0.64，从而生成一幅离轴全息图，如图2(a)所示.分别使用FTA，SCM，10 倍ZP，20 倍ZP，DEA 和本算法对该全息图进行相位恢复，得到6 个恢复结果如图2(b)—(f)所示.由仿真结果可知，6 个恢复结果中，仅有本算法和DEA 能去除载波，得到准确的相位恢复结果.

图2 (a) 光程差100 nm 相位型半球生成的离轴全息图;(b)—(f) 使用FTA，SCM，10 倍ZP，20 倍ZP，DEA 和本算法的恢复结果Fig.2.(a) Off-axis hologram of a phase hemisphere with 100 nm optical path difference;(b)—(f) retrieved phase maps by FTA，SCM，ZP with 10 times zero-padding，ZP with 20 times zero-padding，DEA and the proposed algorithm.

经多次仿真实验可知，本算法对于载波数值的计算精度可达到10—7级;而ZP 法欲达到10—7的精度，需要将图像进行107倍的补零，数据量已经远远超过普通计算机的内存容量，致使该算法无法运行;同时，SCM 法的定位精度严重依赖于所选区域，很难在实际应用中取得良好的效果.因此，虽然本算法的迭代过程耗时较长，但是不需要提前采集背景全息图就能准确地去除载波，具有极高地实用价值.

为了验证本算法的抗噪能力，向生成全息图中加入均值为0、方差为0.005 的高斯白噪声，进行仿真验证.鉴于DEA 算法的优越性，仅使用其和本算法对含噪声的全息图进行恢复，得到恢复结果分别如图3(a)和图3(b)所示.为了便于分析比较，求出恢复结果相较于初始值的残差，如图3(c)和图3(d)所示，并将残差图的峰谷值(peak-valley，PV)和标准差(standard deviation，SD)标注在残差图的右上角.从图中可以看出，相较于DEA，本算法的PV 和SD 更小，具有更好的恢复质量.

图3 确定高斯白噪声影响下的(a)本算法和(b)DEA 的恢复结果，以及(c)本算法和(d)DEA 的残差图Fig.3.Retrieved phase maps from the hologram with Gaussian white noise by (a) the proposed algorithm and (b)DEA，and the corresponding residue maps by (c) the proposed algorithm and (d) DEA.

为了更深入地对比本算法和DEA 的抗噪能力，使用本算法和DEA 对含有0 均值和不同标准差高斯白噪声的全息图进行相位恢复，分别计算两种算法恢复结果与初始值间残差的PV 和SD，获得结果如图4 所示.从图4 中可以看出，随着噪声的放大，两种算法恢复结果的PV 和SD 都在增大，标志着恢复质量也在逐渐降低.但是，即使噪声方差增至0.025 的情况下，本算法仍能较高质量的恢复出物体相位信息，从而证明了本算法具有较强的抗噪干扰能力.

图4 不同噪声情况下本算法和DEA 算法对应残差值的PV 和SDFig.4.PV and SD of the residue maps by the proposed algorithm and DEA with different variance.

由此得知，本文提出的算法在不需要任何背景全息图的条件下，仅仅利用单幅离轴数字全息图即可在去除载波的同时，实现待测相位准确直接恢复.

4 实验结果及分析

为了进一步验证所提算法的有效性，利用FTA，10 倍ZP，SCM，DEA 和本算法分别对实际获取的一幅全息图进行恢复.在实验中，使用分光瞳离轴数字全息系统[27]测量一块反射式的硅片样品，所测得的全息图如图5(a)所示，相位恢复结果分别如图5(b)—(f)所示.从图5(b)—(f)中的恢复结果可知，如同仿真结果，FTA，10 倍ZP 和SCM 均无法准确实现相位恢复，而只有DEA 和本算法得到了较为准确的恢复结果.为了进一步验证本算法的恢复质量，对图5(e)和图5(f)实线处做一维剖面，其剖面数据如图6 所示.本算法高、低位置的标准差分别为0.0613 rad 和0.1218 rad，而DEA 高、低位置的标准差分别为0.1051 rad 和0.1303 rad.由此可见，本算法在不需要额外记录背景全息图的情况下，就能达到和DEA 相近的恢复质量.

图5 硅片实验结果 (a)全息图;利用(b) FTA，(c)10 倍ZP，(d) SCM，(e) DEA 和(f)本算法的恢复结果Fig.5.Experimental results for silicon wafer:(a) Hologram;retrieved phase maps by (b)FTA，(c)ZP with 10 times zero-padding，(d) SCM，(e) DEA and (f) the proposed algorithm.

图6 图5(e)和图5(f)中白线所标剖面数据Fig.6.1D phase profile along the white lines in Fig.5(e)and Fig.5(f).

上述实验已经验证了本算法在边缘突变物体上的有效性，接下来将对酒精蒸发过程中的一帧数据进行处理，以验证本算法对于边缘连续的物体仍具有有效性.图7(a)为酒精蒸发过程中采集到的一帧数据，图7(b)—(f)是分别使用FTA，10 倍ZP，SCM，DEA 和本算法的恢复结果.可见，相比于其他算法，本算法和DEA 仍能更准确的恢复出物体相位信息.为了更直观地展现本算法的恢复质量，将图7(e)和图7(f)中白线标记的剖面数据提取出来，如图8 所示.从图7(b)—(f)所示的恢复结果可以看出，本算法可以实现连续形貌高质量恢复.

图7 酒精实验结果 (a)全息图;利用(b) FTA，(c)10 倍ZP，(d) SCM，(e) DEA 和(f)本算法的恢复结果Fig.7.Experimental results for alcohol:(a) Hologram;retrieved phase maps by (b)FTA，(c)ZP with 10 times zero-padding，(d) SCM，(e) DEA and (f) the proposed algorithm.

图8 图7(e)和(f)中白线标注一维剖面图Fig.8.1D phase profiles along the white lines in Fig.7(e)and Fig.7(f).

5 结论

针对离轴数字全息相位恢复存在的问题，基于SCPS，本文基于SCPS 提出了一种结合线性回归的直接提取载波与物体相位的相位恢复方法.该算法利用SCPS，将一幅离轴全息图分解成4 幅含有载波相移的全息图，再利用线性回归建立起载波和相移之间的关系，最后通过迭代过程恢复出物体相位信息和沿x，y方向的载波信息.仿真与实验表明，提出的方法不需要背景全息图作参考，即可精确地提取载波信息和物体相位信息，且具有与最常用的DEA 相近的恢复效果.该方法不仅不需要先验信息，而且相较于以往的去载波方法具有更高的精度，对于离轴数字全息的实际应用具有重要的意义.