基于求解光强传输方程的数字全息图解调方法

摘要

本发明公开了一种基于求解光强传输方程的数字全息图解调方法，首先重建数字全息图计算得到物平面上的复振幅分布，将该物平面上的复振幅分布用其振幅进行归一化，将归一化后的复振幅采用角谱数值传播法传播一小段距离，得到离焦光强图像，将离焦光强图像作为求解光强传输方程的输入数据，再利用基于快速傅里叶变换的求解法求解光强传输方程，得到物体的相位分布。本发明能够从离轴全息图中直接解调出物体的连续相位，不但有效地避免了传统数字全息中繁杂的去包裹过程，还可有效去除由于物参光路的离轴结构引入的倾斜像差与物光光路中显微物镜引入的二次相位畸变。

说明书

技术领域

本发明属于光学测量中的数字全息成像中的全息图处理技术，特别是一种基于求解 光强传输方程的数字全息图解调方法。

背景技术

相位恢复是光学测量与成像技术的一个重要课题，无论在生物医学还是工业检测领 域，相位成像技术都在发挥着重要的作用。特别是在生物医学显微领域，光学显微镜已 经成为生物学和医学研究中必不可少的工具，是细胞学和细胞生物学得以建立和发展的 重要基础。然而，由于大部分生物细胞的对光波的低吸收作用，采用普通明场光学显微 技术难以实现细胞的清晰成像。另一方面，生物细胞相位分布却对于细胞显微观察、结 构信息的提取以及动力学行为研究十分重要。所以相位测量、或者定量相位成像是解决 该项问题重要技术。

纵观光学测量近半个世纪的进展，最经典的相位测量方法应该非干涉测量法莫属。 数字全息作为干涉法中最具代表性的方法，由于其原理简单、成像速度快、记录再现灵 活、测量精度高，广泛应用于光学检测、材料科学、生物医学成像等领域。然而，数字 全息中重建得到的相位是包裹在(-π,π]之间的，对应于反正切函数的主值区间。为了获 得连续相位分布，必须对其进行相位解包裹。目前有许多相位解包裹的算法被提出([1]D. C.Ghiglia and M.D.Pritt,Two-dimensional phase unwrapping:theory,algorithms,and  software(Wiley New York:,1998).)，但当相位图中存在噪声以及相位“坏点”(phase  residue)时，如当遇到存在大面积相位间断区的复杂、非连续物体时，这些算法的表现 都不尽理想。由于测量过程中不可避免的各种干扰，有时不仅局部区域会出现相位瑕疵， 误差甚至还会沿去包裹方向传播扩散。除此以外以上方法还存在着诸多缺点，如算法复 杂，计算量大和去包裹耗时过长等，很难应用在高速实时测量的场合中。另一方面，由 于离轴倾斜与显微物镜所引入的二次相位畸变叠加在物体的相位上。为了准确获得物体 本身的相位信息，需要对这些像差进行物理或者数值补偿([2]T.Colomb,F.Montfort,J. Kühn,N.Aspert,E.Cuche,A.Marian,F.Charrière,S.Bourquin,P.Marquet,and C. Depeursinge,"Numerical parametric lens for shifting,magnification,and complete aberration  compensation in digital holographic microscopy,"J.Opt.Soc.Am.A 23,3177-3190 (2006).)，这通常难以操作且十分耗时。

发明内容

本发明的目的在于提供一种基于求解光强传输方程的数字全息图解调方法，能够从 离轴全息图中直接解调出物体的连续相位，不但有效地避免了传统数字全息中繁杂的去 包裹过程，还可有效去除由于物参光路的离轴结构引入的倾斜像差与物光光路中显微物 镜引入的二次相位畸变。

实现本发明目的的技术解决方案为：一种基于求解光强传输方程的数字全息图解调 方法，步骤如下：

第一步，重建数字全息图计算得到物平面上的复振幅分布U₀(x,y)；

第二步，将该物平面上的复振幅分布U₀(x,y)按式(7)进行归一化，得到归一化后 的物平面上的复振幅U₀′(x,y)，

$U_0^{\prime }(x,y)=\frac{U_0(x,y)}{\left|U_0(x,y)\right|}---\left(7\right)$

式中||代表取复数的绝对值操作，该步骤将待测物体改变为纯相位物，并假设物面 的光强分布为常数1；

第三步，将归一化后的复振幅U₀′(x,y)采用角谱数值传播法进行传播，得到离焦光强 图像I_Δz；

第四步，将离焦光强图像I_Δz作为求解光强传输方程的输入数据，通过式(11)计算 得到光强轴向微分

$\frac{\partial I(x,y)}{\partial z}\approx \frac{I_{\mathrm{\Delta z}}(x,y)-1}{\mathrm{\Delta z}}---\left(11\right)$

第五步，利用基于快速傅里叶变换的求解法求解光强传输方程，得到物体的相位分 布，即通过式(12)，采用快速傅里叶变换求解光强传输方程，得到最终的相位分布

$\phi (x,y)=F^{-1}\left\{\frac{k(u^2+v^2)F\{\partial I(x,y)/\partial z\}}{4{\pi }^2[{(u^2+v^2)}^2+\gamma ]}\right\}---\left(12\right)$

其中F代表傅里叶变换，γ是正则化参数，正则化参数γ有两种选择，如果需要去 除相位分布中的二次相位畸变，则其取值为γ＝(max(D_x,D_y))^-1，其中D_x与D_y是图像的 物理尺寸；如果不需要去除相位分布中的二次相位畸变，则正则化参数可取大于0的无 穷小量。

本发明与现有技术相比，其显著优点：不但有效地避免了传统数字全息中繁杂的去 包裹过程，还为去除由于物参光路的离轴结构引入的倾斜像差与物光光路中显微物镜引 入的离焦像差提供了一种简单而有效的手段。由于其仅仅采用四次傅里叶变换实现绝对 相位的获取与像差补偿，其相比传统相位解包裹方法简单高效的多，且不会出现由于去 包裹不当产生的相位瑕疵。具有广阔的应用前景。

下面结合附图对本发明作进一步详细描述。

附图说明

图1(a)是共聚焦显微镜测量得到的拓扑图。

图1(b)是如图1(a)所示的经过凸起中心处的一线剖面显示图。

图1(c)是采用单个分束镜的共路全息镜测量得到的包裹相位图。

图1(d)是路径相关相位解包裹法得到的连续相位图。

图1(e)是质量图引导法(quality-map guided)路径无关相位解包裹法得到的连续 相位图。

图1(f)是通过本发明方法得到的连续相位图。

图1(g)是表面形貌三维显示图。

图1(h)是对应于图1(f)中虚线位置的线剖面显示图。

图2(a)是通过频谱平移法得到的包裹相位分布图。

图2(b)是经过将图2(a)中相位解包裹得到的连续相位分布图。

图2(c)是没有经过倾斜校正的包裹相位分布图。

图2(d)是本发明方法得到的(数值反算)包裹相位分布图。

图2(e)是本发明方法结合Tikhonov正则化后得到的连续相位分布图。

图2(f)是图2(e)相位分布的三维显示图。

具体实施方式

本发明基于求解光强传输方程的数字全息图解调方法中，将数字全息图记录时的物 光与参考光的复振幅简化记作O与R，则全息图的强度分布表示为

I_H(x,y)＝|O|²+|R|²+RO^*+R^*O               (1)

式中*代表复数共轭，式(1)的前两项成为零级像，是物光与参考光强度的直接相 加。第三项RO^*是实像，或者-1级像，它包含了物光共轭像的信息，表现为物体本身 的一个镜像图像。第四项R^*O是虚像，或者+1级像，包含了物光本身的原始信息，与 原始物体直接成比例，是最利于直接观察的。实像与虚像二者统称为孪生像。对于离轴 记录结构的全息图，假设参考光波是一个理想平面波，且相对于物光具有一个夹角θ， 即

R(x,y)＝|R|exp(-jk sinθx)               (2) 则式(1)可以表示为

I_H(x,y)＝|O|²+|R|²+|R|Oexp(-jk sinθx)+|R|O^*exp(jk sinθx)      (3) 式中j为虚数单位，k为波数，θ为物参光的夹角。如果考虑全息图的傅里叶变换，那 么相位因子exp(-jk sinθx)的作用即是将虚像与实像的傅里叶频谱在频域进行了镜像平 移。当夹角θ足够大时，零级像的中心位于频域原点，而两项干涉项就会在频谱面上分 离，相对频域原点呈镜像分布。

下面详细描述本发明基于求解光强传输方程的数字全息图解调方法的具体实施步 骤：

第一步，重建数字全息图I_H(x,y)得到物平面上的复振幅分布。具体操作过程为： 首先将全息图I_H(x,y)进行傅里叶变换，然后采用一个频域带通滤波器去单独筛选出虚 像所对应的频谱(由前所述，对于离轴全息图，虚像、实像与零级项在频谱面上是分离 的，所以采用频域带通滤波器是很容易将其筛选出来的)，再逆傅里叶变换回空域获得 其全息图平面的复振幅分布U_z(x,y)。这个过程可以表述为式(4)

U_z(x,y)＝|R|Oexp(-jk sinθx)＝F^-1{F{I_H(x,y)}×H(u,v)}      (4) 式中傅里叶变换与反变换分别被记作F与F^-1，下标z代表了全息图平面与物平面之 间的距离；(u,v)代表相对于(x,y)的频域坐标；H(u,v)代表了频域带通滤波器的传输函 数，能够准确选筛出虚像所在的频谱区域。当然为了避免高频细节的丢失，其尺寸应该 尽可能的大，并且利用相对平滑的传输函数，如采用高斯(Gaussian)或汉宁(Hanning) 窗，可以有效避免重建图像中由于频谱截断而产生的吉布斯(Gibbs)效应。

然后将全息图平面的复振幅分布U_z(x,y)，采用式(5)将其“反传播”-z的距离，从 而就回到了其本身所在的物平面上，得到物平面上的复振幅分布U₀(x,y)

U₀(x,y)＝F^-1{F{U_z(x,y)}H_-z(u,v)},             (5) 式中傅里叶变换与反变换分别被记作F与F^-1，H_-z(u,v)是角谱传输函数，其形式为

$H_{-z}(u,v)=\exp [-j\frac{2\mathrm{\pi z}}{\lambda }\sqrt{1-{\left(\mathrm{\lambda u}\right)}^2-{\left(\mathrm{\lambda v}\right)}^2}]---\left(6\right)$

式中λ为光波波长，j为虚数单位，z代表了全息图平面与物平面之间的距离。

第二步，将第一步所得的物平面上的复振幅分布U₀(x,y)按式(7)进行归一化，得 到归一化后的物平面上的复振幅U₀′(x,y)，

$U_0^{\prime }(x,y)=\frac{U_0(x,y)}{\left|U_0(x,y)\right|}---\left(7\right)$

式中||代表取复数的绝对值操作，该步骤将待测物体改变为纯相位物，并假设物面的光 强分布为常数1。

第三步，将归一化后的物平面上的复振幅U₀′(x,y)采用式(8)与式(9)所表示的角 谱数值传播法传播Δz的距离(这里的Δz要取得很小(推荐为1μm)，以保证第四步中 有限差分逼近的准确度)，得到离焦光强图像I_Δz。

U_Δz(x,y)＝F^-1{F{U₀′(x,y))}H_Δz(u,v)},           (8)

I_Δz(x,y)＝|U₀(x,y)|²              (9)

式中傅里叶变换与反变换分别被记作F与F^-1，H_-z(u,v)是角谱传输函数，其形式为

$H_{\mathrm{\Delta z}}(u,v)=\exp \left[j\frac{2\mathrm{\pi \Delta z}}{\lambda }\sqrt{1-{\left(\mathrm{\lambda u}\right)}^2-{\left(\mathrm{\lambda v}\right)}^2}\right]---\left(10\right)$

第四步，将离焦光强图像I_Δz通过式(11)计算得到光强轴向微分

$\frac{\partial I(x,y)}{\partial z}\approx \frac{I_{\mathrm{\Delta z}}(x,y)-1}{\mathrm{\Delta z}}---\left(11\right)$

第五步，通过式(12)，采用快速傅里叶变换求解光强传输方程，得到最终的相位 分布

$\phi (x,y)=F^{-1}\left\{\frac{k(u^2+v^2)F\{\partial I(x,y)/\partial z\}}{4{\pi }^2[{(u^2+v^2)}^2+\gamma ]}\right\}---\left(12\right)$

其中F代表傅里叶变换，γ是正则化参数。正则化参数γ有两种选择，如果需要去除相 位分布中的二次相位畸变，则其取值为γ＝(max(D_x,D_y))^-1，其中D_x与D_y是图像的物理 尺寸。如果不需要去除相位分布中的二次相位畸变，则正则化参数可取一个很小的数值 (大于0的无穷小量)，如小于10^-10即可。

通过式(12)所得到的相位分布直接就是连续的，不存在相位包裹，从而也就避免 了复杂的相位解包裹操作以及解包裹操作本身可能导致的一系列相位瑕疵问题。此外由 于倾斜像差本身并不会导致光强轴向微分信号的改变，本发明方法可以自动地去除倾斜 像差，这是其另一大优势。当倾斜像差被准确去除后，残余的二次相位畸变则可通过设 置正则化参数γ有效消除。此外，由于本发明仅仅采用单幅离轴全息图，算法本身借助 于快速傅里叶变换实现，十分简单快捷，使其非常适于动态样品的测量。

实施例

下面用一实施例来说明本发明所提出的基于求解光强传输方程的数字全息图解调方 法的有效性。待测样品为一规则的反射凸起阵列。该阵列是在硅衬底上采用激光诱导非 侵入式结构加工技术在0.75个大气压下的纯氧环境中加工而成的。选取该样品的目的是 为了表现本发明算法在复杂相位与噪声情况下的鲁棒性。该微结构形貌的准确表征对于 其在集成光学、微机电系统、以及太阳能电池领域的应用至关重要。该微结构可以通过 共聚焦显微镜或者原子力显微镜测量得到，但是这些方法均需要大量的扫描时间。

图1(a)是采用共聚焦显微镜测量得到的拓扑图；图1(b)显示了经过凸起中心位置的 一剖面图。本实验采用了基于单个分束器的共路全息镜采集了离轴全息图。由于物参共 路结构本身可以自动补偿物镜引入的二次相位畸变，所以只需要采用频谱平移法对倾斜 误差进行补偿即可。图1(c)给出了去除倾斜之后得到的包裹相位，可以看出其中包括大 量的“坏点”(围绕其以任意路径进行回路积分结果不为0)，导致传统方法难以对该相位 图进行正确去包裹。图1(d)是采用一种简单的与路径相关的去包裹方法(逐行展开法) 得到的连续相位、由于误差累积与扩散，结果中出现了严重的“拉线”效应。采用更为复 杂的质量图引导法可以防止误差的传播与扩散，但是其仍然不能从根本上避免错误的出 现。如图1(e)所示，一些相位变化较快的区域(特别是在微小凸起接近顶部的区域，显 示于放大区域)出现了明显的去包裹错误。而采用发明方法解调的全息图给出了完整无 缺的连续相位，不存在任何瑕疵区域与误差传播现象，如图1(f)所示。通过相位-高度转 换后，可以得到样品的三维拓扑，显示于图1(g)。图1(h)给出了对应于图1(f)中虚线位 置的线剖面显示，从结果上来看与共聚焦显微镜符合得很好。由于本发明方法仅采用4 次快速傅里叶变换(两次用于数值传播，另外两次用于求解光强传输方程)，其运算复 杂度仅为4N²logN。该实施例中全息图的尺寸为1280×960，采用2.5GHz主频的笔记 本电脑在MATLAB环境下进行处理。方法处理总共所需时间为321.3ms，仅为传统的 质量图引导相位解包裹法的六分之一(1.971s)。

本发明采用基于迈克尔逊干涉仪结构的全息显微镜对固定在载玻片上的巨噬细胞进 行了测量。物参光相位曲率失配导致的二次相位畸变展宽了±1级的频谱，这让传统频 谱平移法难以对倾斜像差进行正确补偿。图2(a)给出了通过频谱平移法补偿后的包裹相 位分布；可以看出相位图中出现了许多偏心的同心圆环，这代表了倾斜误差没有完全去 除。图2(b)中给出了采用质量图引导法去包裹得到的连续相位分布。图中可以明显看到 一些去包裹误差。此外由于二次相位畸变的存在，很难准确得到物体本身的相位信息。 为了显示本发明算法对于像差的校正能力，我们直接将其应用于没有经过倾斜补偿(不 经过任何频谱平移)的滤波后的全息图中。此时包裹相位中将包含一个非常高频的倾斜 载波相位，如图2(c)所示。由于本发明方法可以直接获得连续相位，为了方便显示与比 较，将其反算回包裹相位，如图2(d)所示。相比于图2(a)，可以看出本发明的方法准确 地自动去除了高频倾斜像差，仅仅留下了由于二次相位畸变所致的位于图像中心的同心 圆环。如果在求解光强传输方程的过程中引入正则化(正则化因子为γ＝2×10^-3μm^-1)， 弯曲的背景将会有效地变平，如图2(e)与图2(f)所示。结果表明本发明方法不仅有效去 除了相位畸变，结果中也不存在任何由于去包裹不当所导致的相位瑕疵。