摘要:散射介质对光的散射是当前限制光学成像深度或距离的一个严重的问题. 本文首先数值模拟比较了光透过随机散射介质成像研究中常用的基于光学记忆效应(memory effect, ME)和自相关(autocorrelation, AC)方法的HIO&ER算法和乒乓(Ping-Pang, PP)算法的优缺点. 通过对HIO&ER算法和PP算法的恢复效果和迭代次数进行比较, 发现PP算法在保持较高恢复效果的前提下拥有更快的运行速度. 实验中, 利用连续He-Ne激光器和旋转毛玻璃产生赝热光源, 通过物镜对随机散射介质后数毫米距离内的不同形状物体进行了单帧成像, 并采用PP算法成功地恢复出微米量级物体的实际图像. 这一研究结果将进一步促进ME和AC方法在深层生物组织医学成像研究上的应用. 最后, 实验研究了不同的物镜和散射介质的间距对成像恢复的放大率、分辨率和图像强度的影响特性, 并进行了详细研究.
关键词: 记忆效应/
强度自相关/
随机散射/
相位恢复算法

English Abstract

全文HTML

--> --> -->

随着纳米技术、光子学技术和生物医学技术的发展, 活体光学成像在生物光子及现代医疗科学领域的应用越来越广泛^[1-3]. 然而, 由于生物组织对光具有吸收和散射作用, 激发光和信号光都难以穿透足够深度的生物组织, 并且探测光在生物组织内的散射会引起其相位信息的破坏, 因而难以直接通过相机来对深层生物组织进行高分辨成像^[4-6]. 散射介质除了引入以上问题, 研究也发现光经过随机散射介质可以增加实际成像系统的数值孔径、增大成像系统的视场角、接收来自物体表面的倏逝波并将其散射为行波在远场传播, 因而散射介质也被一些研究小组用来进行光学超衍射成像^[7-9]. 多年来, 国内外的一些研究小组利用波前调制等反馈控制调节法^[10,11]及相位共轭时间反演法^[12-14], 来实现光经过散射介质的深度成像. 此外, 在早期鬼成像的研究中, 就已经通过光学互相关的原理来减小散射介质对成像的影响^[15-17]. 但这些方法的光路系统相对复杂, 容易受环境影响, 难以适应生物组织深度成像的实际应用.
最近的研究惊喜地发现, 利用赝热光源基于“光学记忆效应”(memory effect, ME)^[18-22]和自相关(autocorrelation, AC)方法^[23-25], 可以从探测器探测到的杂乱无章的散斑图中直接恢复出物体的二维图像信息. 其中, 非常重要的一点是用到了Gerchberg-Saxton (GS)算法^[26-29]来进行物体的相位恢复. 本文对GS算法中的HIO&ER算法^[30-32]和PP (Ping-Pang, PP)算法^[33]的恢复效果和迭代次数进行了比较, 结果表明PP算法在保持较高恢复效果的前提下拥有更快的运行速度. 实验上, 采用PP算法, 并通过前置物镜将赝热光源聚焦到特定的成像目标上, 通过后置物镜将经过散射介质的散射光收集到相机上, 其物体和散射介质间距可达毫米量级. 同时, 通过平移台移动成像物体即可对其不同部位进行扫描成像, 本文成功恢复出了不同形状的微米量级物体的实际图像.

2

-->

2.1.随机散射介质特性

光波经过随机散射介质会发生散射效应, 即一部分光在散射介质中随机传播, 导致其传播的相位面发生畸变, 不能清晰成像, 最终形成散斑图, 如图1所示. 散斑指光束经过随机散射介质, 并在介质内发生多重散射后, 产生的随机无序的颗粒状图样, 本质上是一种干涉现象. 因为观察点的光场是散射介质颗粒上各点发出的相干子波的叠加, 且光波波长小于散射介质颗粒尺寸, 所以到达观察点的各个相干子波的相位是随机分布的, 随机相位的相干叠加就产生了散斑的随机强度图样.
图 1 散斑产生示意图
Figure1. Schematic of speckle generation.

但实质上, 携带物体信息的光波经过随机散射介质作用后, 物体信息并没有完全丢失, 它们只是在散射过程中进行了随机组合. 因此, 本文假设散射介质是一个光学透镜, 利用它的光学传播特性, 从散射光中重组并恢复出物体的原有信息. 例如波前整形方法便可有效地补偿扭曲的波阵面, 但该方法需要空间光调制器对扭曲波形面的相位进行逐一调制, 其操作过程耗时复杂, 不适于实时成像.
2 -->

2.2.记忆效应

当入射光照射到散射介质上时, 其空间相位信息将会被打乱因而具有随机性, 但其中也会包含相应的相关性. 对于入射光束, 其相对散射介质的入射角度$\theta = \left( {{\theta _x},{\theta _y}} \right)$可控, 物体经过散射介质后所成的像与入射角度有关,

$I\left( \theta \right) = \int_{ - \infty }^\infty {O\left( r \right)} S\left( {r - \theta d} \right){{\rm{d}}^2}r = \left[ {O*S} \right]\left( \theta \right),$

式中$I\left( \theta \right)$表示相机探测到的强度图, $O\left( r \right)$表示物体的实际像, $S\left( r \right)$表示散斑强度, $r$表示空间坐标矢量, $\theta $表示入射光束的入射角, $d$表示物体与散射介质的间距, $*$表示卷积操作, 即物体透过散射介质的像为物体实际像和散斑强度的卷积.
为了将物体实际像信息从随机散斑中提取出来, 对探测到的光强信息做AC操作,

$\begin{split}& \left\langle {I \otimes I} \right\rangle \left( {\text{δ} \theta } \right) = \left\langle {O*S} \right\rangle \otimes \left\langle {O*S} \right\rangle \\ & = \left\langle {O \otimes O} \right\rangle *\left\langle {S \otimes S} \right\rangle = \left[ {O \otimes O} \right]*\left\langle {S \otimes S} \right\rangle ,\end{split}$

式中$ \otimes $表示互相关操作, $\left\langle \cdot \right\rangle $表示取平均操作.
假定入射光束宽度为$w$, 则

$\begin{split}\left\langle {I \otimes I} \right\rangle \left( \theta \right) = & {\left( {\frac{{k\left| \theta \right|L}}{{\sinh \left( {k\left| \theta \right|L} \right)}}} \right)^2} \\ & \times \left[ {\left[ {O \otimes O} \right]*{{\left( {\frac{{2{{\rm{J}}_1}\left( {k\left| \theta \right|w} \right)}}{{k\left| \theta \right|w}}} \right)}^2}} \right]\left( \theta \right),\end{split}$

式中, ${{\rm{J}}_1}$是一阶贝塞尔函数; $L$是散射介质厚度, 远大于平均自由程; $k = 2{\text{π}}/\lambda $是波数. (3)式中第一部分可以用来表示图案相关性程度$C\left( \theta \right) = $${\left( {\dfrac{{k\left| \theta \right|L}}{{\sinh \left( {k\left| \theta \right|L} \right)}}} \right)^2}$, 对于ME, 其小角度条件需满足$\left| \theta \right| \leqslant \lambda /2{\text{π}}L$; 第二部分中${\left( {\dfrac{{2{{\rm J}_1}\left( {k\left| \theta \right|w} \right)}}{{k\left| \theta \right|w}}} \right)^2}$表示平均散斑大小, 通过增大$w$可无限接近于衍射极限.
假定散斑在空间分布是随机的, 则散斑的AC结果是一个${\text{δ}}$函数, 即$\left\langle {S \otimes S} \right\rangle = {\text{δ}}\left( r \right)$, 由此可将散斑信息从成像AC中消除, 从而得到物体实际像AC的近似值$\left\langle {I \otimes I} \right\rangle \left( \theta \right) \approx \left[ {O \otimes O} \right]$. 根据卷积定理, 可得

$F\left\{ {\left\langle {I \otimes I} \right\rangle } \right\} \!=\! F\left\{ {O \otimes O} \right\} \!=\! F\{ O\} F{\{ O\} ^*} \!=\! {\left| {F\{ O\} } \right|^2}.$

2 -->

2.3.相位恢复算法

根据相机采集到的光强图像信息, 首先对其进行滤波等预处理操作, 得到相机光强图像$I\left( {x,y} \right)$. 根据卷积定理, 对其做AC处理可得

$A\left( {x,y} \right) \!=\! I\left( {x,y} \right) \otimes I\left( {x,y} \right) \!=\! {\rm{F}}{{\rm{T}}^{ - 1}}\left\{ {{{\left| {{\rm{FT}}\left\{ {I\left( {x,y} \right)} \right\}} \right|}^2}} \right\}.$

根据维纳-辛钦定理, 物体的能量谱等于物体AC的傅里叶变换振幅大小, 因此可以通过矩形窗口$W\left( {x,y} \right)$截取AC的中心部分, 再对其进行二维傅里叶变换, 即可得到物体的能量谱

$S\left( {{k_x},{k_y}} \right) \!=\! \left| {{\rm{FT}}\left\{ {W\left( {x,y} \right)A\left( {x,y} \right)} \right\}} \right|\! \approx \!{\left| {{\rm{FT}}\left\{ {O\left( {x,y} \right)} \right\}} \right|^2}.$

由于前面所述的图像信息均是强度信息, 因此在计算推导傅里叶变换过程中丢失了图像的相位信息. 由此, 将采用GS算法对丢失的相位信息进行恢复, 在迭代过程中需要将得到的物体能量谱开根作为替代模量不断替换傅里叶变换的模量. 基本的GS算法分为如下五步:
1)对随机相位值${g_k}\left( {x,y} \right)$做傅里叶变换

${G_k}\left( {{k_x},{k_y}} \right) = {\rm{FT}}\left\{ {{g_k}\left( {x,y} \right)} \right\};$

2)求得傅里叶变换域的角度

${\theta _k}\left( {{k_x},{k_y}} \right) = \arg \left\{ {{G_k}\left( {{k_x},{k_y}} \right)} \right\};$

3)用测量得到的物体能量谱开根替换傅里叶变换的模量

$G_k'\left( {{k_x},{k_y}} \right) = \sqrt {S\left( {{k_x},{k_y}} \right)} {{\rm{e}}^{{\rm{i}}{\theta _k}\left( {{k_x},{k_y}} \right)}};$

4)对(9)式做傅里叶逆变换

$g_k'\left( {x,y} \right) = {\rm{F}}{{\rm{T}}^{ - 1}}\left\{ {G_k'\left( {{k_x},{k_y}} \right)} \right\};$

5)根据物理约束条件, 生成新的迭代相位值

$g_k'\left( {x,y} \right) \to {g_{k + 1}}\left( {x,y} \right).$

这里, ${g_k}\left( {x,y} \right)$表示物体的估计相位值, ${G_k}\left( {{k_x},{k_y}} \right)$表示对${g_k}\left( {x,y} \right)$进行傅里叶变换得到的模量, ${\theta _k}\left( {{k_x},{k_y}} \right)$表示${G_k}\left( {{k_x},{k_y}} \right)$所对应的角度信息, $G_k'\left( {{k_x},{k_y}} \right)$表示新得到的模量信息, $g_k'\left( {x,y} \right)$表示对$G_k'\left( {{k_x},{k_y}} \right)$进行逆傅里叶变换得到的相位, ${g_{k + 1}}\left( {x,y} \right)$表示下一次迭代的输入相位值.
根据(11)式的不同迭代公式, 可将其分为Error Reduction (ER)和Hybrid Input-Output (HIO)两种算法:
1) ER算法

${g_{k + 1}}\left( {x,y} \right) = \left\{ {\begin{aligned}& {g_k'\left( {x,y} \right),}\quad {\left( {x,y} \right) \notin \varGamma ,}\\& {0,}\quad\quad\quad\quad\, {\left( {x,y} \right) \in \varGamma ;}\end{aligned}} \right.$

2) HIO算法

${g_{k + 1}}\left( {x,y} \right) = \left\{ {\begin{aligned}&{g_k'\left( {x,y} \right),} \quad\quad\quad\quad\quad\quad\; {\left( {x,y} \right) \notin \varGamma ,}\\&{{g_k}\left( {x,y} \right) - \beta g_k'\left( {x,y} \right),}\quad {\left( {x,y} \right) \in \varGamma.}\end{aligned}} \right.$

这里$\varGamma $表示不满足非负的实数的集合.
2 -->

2.4.Ping-Pang算法

由上可见, ER和HIO方法各有优缺点, ER方法是目前唯一在数学上被证明的解决相位问题的方法, 但该方法收敛速度慢、易受噪声影响; HIO方法则是目前应用最广泛的一种方法, 其算法简单、运行效率高. 鉴于上述两种算法的优缺点, 将HIO和ER两种算法结合, 衍生出更具潜力的PP算法. 如图2所示, 首先运用HIO算法对随机预测的初始相位进行迭代, 将HIO算法迭代的结果${g_{{\rm{middle}}}}\left( {x,y} \right)$作为ER算法的输入, 利用ER算法进行迭代, 最后根据物理约束条件来确定是否需要继续迭代. 在PP算法中, 每次迭代均先后进行HIO算法和ER算法, 利用两种算法的优点, 使得总迭代次数较小, 可以极大地缩短恢复算法运行时间.
图 2 相位恢复算法框图
Figure2. Schematic of phase retrieval algorithm.

图3是按照图2所示的相位恢复算法进行数值模拟的过程, 该方法截取$20 \times 20$像素大小的数字图像和$3600 \times 3600$像素大小的散射介质点扩散函数进行卷积, 模拟物像经过散射介质后的散斑图. 图3(a)为标准数字“5”图像; 图3(b)为模拟的散射介质的点扩散函数; 图3(c)为用$1000 \times 1000$像素大小窗口截取点扩散函数后的卷积合成的带有数字“5”的散斑图; 图3(d)为散射介质的点扩散函数的三维AC结果, 可以看出其是一个${\text{δ}}$函数, 说明在AC操作下散射介质对物体图像信息没有影响; 图3(e)为标准数字“5”的AC结果; 图3(f)为散斑图的AC结果, 与标准数字的AC相同; 图3(g)为优化的能量谱开根; 图3(h)为最终重建出来的数字“5”.
图 3 成像过程的数值模拟　(a) 物体; (b) 点扩散函数; (c) 散斑图; (d) 点扩散函数AC; (e) 物体AC; (f) 散斑AC; (g) 能量谱开根; (h) 重建结果
Figure3. Simulations of imaging process: (a) Object; (b) point diffusion function; (c) speckle pattern; (d) AC of point diffusion function; (e) AC of object; (f) AC of speckle pattern; (g) square root of power spectrum; (h) result of reconstruction.

图3的成像恢复数值模拟使用的是基于GS算法的HIO算法和ER算法的顺序叠加, 其最终的迭代参数设定为物理约束次数$N = 30$, HIO约束回归系数$\beta = 1: - 0.04:0$, 即先对得到的能量谱开根图像进行$30 \times \left\{ {\left( {1/0.04} \right) + 1} \right\} = 780$次的HIO迭代, 再对HIO的迭代结果进行30次的ER迭代, 最终结果如图3(h)所示.
为了进一步探究迭代次数以及不同迭代算法对恢复结果的影响, 图4还利用matlab对其不同情形进行仿真. 图4(a)—(c)是HIO & ER算法恢复结果, 图4(d)—(f)是PP算法恢复结果, 参数设定如表1所列.

Algorithm	Physical contraint N	Physical contraint $\beta$	Interation times
HIO&ER	30	$1: - 0.02:0$	1560
	30	$1: - 0.04:0$	810
	30	$1: - 0.05:0$	660
PP	30	$3: - 0.02:1$	202
	30	$3: - 0.05:1$	82
	30	$3: - 0.1:1$	42

表1不同情形下算法迭代次数
Table1.Interation times of algorithm in different conditions.

图 4 不同迭代次数下的恢复效果　(a)—(c) HIO&ER算法的恢复结果, 其中, (a) $\beta = 1: - 0.02:0$, (b) $\beta = 1: - 0.04:0$, (c) $\beta = 1: - 0.05:0$; (d)—(f) PP算法的恢复结果, 其中, (d) $\beta = 3: - 0.02:1$, (e) $\beta = 3: - 0.05:1$, (f) $\beta = 3: - 0.1:1$
Figure4. Retrieval results in different interation times: (a)?(c) Retrieval results of HIO&ER algorithm when (a) $\beta = 1: - 0.02:0$, (b) $\beta = 1: - 0.04:0$, (c) $\beta = 1: - 0.05:0$; (d)?(f) Retrieval results of PP algorithm when (d) $\beta = 3: - 0.02:1$, (e) $\beta = 3: - 0.05:1$, (f) $\beta = 3: - 0.1:1$.

通过对比图4(a)—(c)和表1可以看出, 随着约束回归系数$\beta $步长的增加, 总迭代次数减少, 表示其运行速度加快, 但图像恢复效果也逐渐下降. 图4(d)是采用PP算法进行恢复的结果, 可以看出其恢复效果接近图4(a)的效果, 但总迭代次数仅为202次, 相较于图4(a)有近8倍的速度提升, 进一步证明PP算法的优势. 同时, 随着约束回归系数$\beta $步长的增加, 迭代次数进一步减少, 但恢复效果更差, 图像出现模糊, 如图4(e)和图4(f)所示, 因而PP算法也难以无限制地缩短运行时间.

本文通过散射介质分别对标准分辨率板上的不同数字进行成像, 并通过PP算法对其进行图像恢复, 验证其对透过散射介质的成像能力.
2

-->

4.1.不同物体的恢复效果

实验光路图如图5(a)所示, 实验中采用He-Ne激光器(632.8 nm, $\phi $0.48 mm, ThorLabs)作为光源, 经空间扩束系统(透镜1焦距为25 mm, 透镜2焦距为200 mm)进行8倍的扩束后入射到旋转的散射片(600砂, $\phi $2口径, ThorLabs)上, 该散射片利用电机以20 Hz的转速驱动会将扩束之后连续光的空间相干性打乱, 产生实验所需的赝热光, 效果如图5(b)所示.
图 5 通过散射介质成像的光学装置　(a)实验光路图; (b)赝热光的产生; (c)散斑的产生
Figure5. Optical setups used for imaging through the scattering media: (a) Optical path in experiment; (b) generation of pseudothermal light source; (c) generation of speckle pattern.

产生的赝热光经由一个前置显微物镜(× 20, 0.4NA, 工作距1.2 mm, Olympus)聚焦到标准分辨率板(1951 USAF负测试靶, $\phi $1口径, ThorLabs)的特定数字上, 其中, 分辨率板可以利用一手动平移台移动位置, 实现前置显微镜对分辨率板上不同数字图像的聚焦提取. 透过分辨率板上不同数字的光, 又经过一片距离分辨率板8 mm远的散射片(220砂, $\phi $1口径, ThorLabs), 使得物体像信息模糊紊乱形成散斑图, 经过散射片后的实际效果如图5(c)所示. 这些散斑图最后通过后置显微物镜2 (× 10, 0.25NA, 工作距10.6 mm, Olympus)收集到相机sCMOS(Quantalux?黑白相机, 1920×1080像素, USB 3.0接口, ThorLabs)上采集成像, 获得的散斑数据输送到PC机上进行算法恢复.
图6是根据图5所述光学装置进行实验的结果, 其中图6(a)—(e)分别为(a)原始物像, (b)相机接收到的散斑图, (c)散斑图的AC结果, (d)光谱能量的根值, (e)通过PP相位恢复算法恢复出来的物体幅度信息. 在恢复过程中, 首先截取相机采集到的散斑图中心$1050 \times 1050$像素大小的区域, 对该区域图像做预处理(滤波、归一化处理等, 第二列图), 然后对该区域进行AC运算操作(第三列图)并取得近似的物体能量谱信息(第四列图), 最后通过PP算法进行相位恢复(第五列图), 为了得到较好的恢复效果, 在此设定的物理约束条件为$N = 30$, $\beta = 3: - 0.01:1$, 总迭代次数为402次. 图6(f)—(t)是对数字“3”,“5”,“6”的恢复过程, 这些数据也验证了该方法对不同形状物体成像具有普适性.
图 6 不同数字的实验结果　(a)—(e)数字“1”的恢复过程, 其中, (a)物体, (b) sCOMS成像, (c)散斑AC, (d)能量谱开根, (e)重建结果; (f)—(t)数字“3”, “5”, “6”的恢复过程
Figure6. Experimental results for different numbers: (a)?(e) Retrieval process of number “1”, namely, (a) object, (b) sCOMS image, (c) autocorrelaction of speckle pattern, (d) square root of power spectrum, (e) result of reconstruction; (f)?(t) retrieval processes of number “3”, “5” and “6”.

2 -->

4.2.物镜和散射介质间距的影响

本实验系统采用前置物镜将赝热光源聚焦到特定的成像目标上, 并通过后置物镜将经过散射介质的散射光收集到相机上. 为了探究物镜和散射介质间距对物体成像质量的影响, 实验调节后置显微镜焦平面与散射介质的间距, 利用PP算法对sCMOS采集到的散斑图进行恢复.
如图7所示, 其中图7(a)—(f)表示物镜焦平面与散射介质的间距分别为$700$, $900$, $1100$, $1300$, $1500$和$1700\;{\text{μ}\rm{m}}$时对应的sCMOS采集到的散斑图的AC结果, 图7(g)—(l)表示相对应的最终恢复效果. 通过对比可以发现, 随着物镜焦平面和散射介质间距的增加, 所得到的AC图和恢复结果图都有相应的放大, 这表明物镜对成像具有放大作用, 且放大程度与物镜焦平面和散射介质间距有关. 此外, 随着间距的增加, sCOMS上收集到的物体采样点信息增加, 分辨率得到进一步提高. 但由于物镜的远离, 其所收集到的经过散射介质后的散射光强变小, 最后恢复出来的图像强度有所减弱. 因此, 实际情况中, 应根据需要来选取合适的散射片与物镜的距离, 以平衡成像分辨率和图像强度, 实现最好的效果.
图 7 不同物镜和散射介质间距对成像效果的影响　(a)—(f)不同间距下的散斑AC结果; (g)—(l)不同间距下的恢复结果
Figure7. Effects of different distance between objective and diffuser: (a)?(f) AC results of speckle pattern in different distance; (g)?(l) retrieval results in different distance.

本文研究了一种基于赝热光照明的散斑成像方法, 该方法利用ME和AC原理消除了散射介质对物体成像的影响, 利用相位恢复算法通过单幅散斑图即可实现物体成像的恢复. 同时, 针对传统基于GS算法的HIO和ER算法迭代次数多的缺陷, 设计出一种快速高效的PP算法, 该算法在保持较高分辨率的同时可以有效缩短运行时间, 有利于实现生物组织的实时成像.
采用$20 \times 20$像素大小的数字“5”图像和$3600 \times 3600$像素大小的散射介质点扩散函数进行卷积模拟散斑图, 通过数值模拟验证了基于相位恢复算法成像的可行性. 对HIO&ER算法和PP算法的成像效果进行了模拟对比, 结果表明PP算法相对于HIO&ER算法有近8倍的速度提升, 证明PP算法具有更高的效率. 通过实验对不同形状的物体进行了散斑成像, 利用PP算法获得微米量级物体的图像恢复结果. 同时, 通过改变物镜焦平面与散射介质的间距, 使其从700 ${\text{μ}\rm{m}}$增加到1700 $ {\text{μ}\rm{m}}$, 发现最终物体散斑成像的放大率和分辨率有了相应的提高, 而图像强度有所下降, 这对今后实际生物组织医学成像的发展具有重要指导意义.