基于波动光学的显微光场成像点扩散函数

2019-09-10顾梦涛宋祥磊张彪唐志永许传龙

北京航空航天大学学报 2019年8期

顾梦涛，宋祥磊，张彪，唐志永，许传龙，*

（1.东南大学能源与环境学院能源热转换及其过程测控教育部重点实验室，南京210096；2.中国科学院上海高等研究院低碳转化科学与工程重点实验室，上海201210）

航空航天发动机由于换热量大，且内部结构复杂，传统换热器已经不能满足其对空间及换热性能的双重需求，而微通道换热器因换热性能强、结构紧凑、质量轻、体积小等优点成为研究热点。目前，微通道换热器仍处于发展阶段，研究其流动特性可指导微通道的构型优化，且流动特性的研究有助于换热特性的研究。因此，流动特性是微通道换热器重要的研究内容之一［1］。由于微通道流动特征尺度微小，流体分子间作用力、静电力等表面力效应增强，同时受到微通道构型、壁面粗糙度和浸润性等因素影响，微尺度流动特性非常复杂，目前还没有合适的理论模型来对微通道内复杂的流动特性进行解释［2-3］。传统的接触式测量会侵入微通道影响流动，无法完成精确测量，而显微粒子图像测速（Micro-PIV）技术可实现非接触、全场、瞬态、定量的微尺度速度场测量，是流动可视化的实验方法之一，特别是三维Micro-PIV技术可以揭示微通道内的三维流动，对研究微尺度流动有重要意义［4］。现有的三维Micro-PIV技术或设备复杂或局限于周期性流场的测量，无法利用单相机完成非定常和非周期性流场的测量［5-6］。

近年来，光场技术进一步发展，利用微透镜阵列，单相机一次曝光能同时完成光线的位置和方向信息（即光场）的记录。因此，在传统显微成像系统的基础上加装微透镜阵列构成显微光场成像系统［7］，可以实现单相机一次曝光记录物空间的光场信息，进一步利用图像处理技术可以得到物空间的三维坐标信息［8］，因此光场显微成像结合粒子图像测速技术（光场Micro-PIV），通过单个光场相机有望实现微尺度瞬时三维速度场的测量。微尺度流场中，示踪粒子的三维坐标信息的准确重建是光场Micro-PIV技术实现瞬时三维速度场测量的关键。示踪粒子的成像过程可表示为示踪粒子的光强分布与成像系统的点扩散函数（PSF）的卷积［9］，因此可通过反卷积进行示踪粒子三维重建。反卷积算法主要有线性反卷积法［10］、盲反卷积法［11］和Lucy-Richardson（L-R）算法［12］。线性反卷积法对噪声敏感度高且无法复原高频信息，重建精度低；盲反卷积法虽然精度高，但需在迭代重建的同时更新PSF，计算量大，耗时长；L-R算法能够按照泊松噪声统计标准对给定PSF的退化图像进行反卷积迭代推演，当给定的PSF足够精确时，可以准确而快速地完成三维重建。但准确获得PSF是L-R算法完成示踪粒子三维重建的关键。

现有的PSF模型可分为两大类：基于几何光学的PSF模型和基于波动光学的PSF模型。由于显微光场成像系统衍射无法忽略，同时系统中微透镜阵列的引入，导致现有的基于几何光学和波动光学的PSF模型均不适用于光场Micro-PIV中PSF的确定。因此，需要建立适合显微光场成像系统的PSF模型。目前，PSF可通过理论计算和实验测量方法确定。对于显微光场成像系统，实验测量方法要求粒子直径小于衍射限制，同时，当系统参数变化时，需要重新进行实验测量，导致PSF确定过程极为耗时。而理论计算方法没有上述限制，且系统参数变化时，只需改变PSF模型中相应的参数，即可计算出新的PSF。

针对光场Micro-PIV中PSF的确定，本文基于波动光学建立了显微光场成像PSF模型，进行了数值模拟，获得了理论的PSF图像。构建了光场Micro-PIV实验装置，实验确定了显微光场成像系统的PSF图像，进一步通过L-R算法实现了微尺度三维流动中示踪粒子的三维重建，对PSF理论计算模型的准确性进行了验证。

1 光场M icro-PIV

光场的概念最早由Gershun［13］提出，用来描述光在三维空间中的辐射传输特性。当光线在传播过程中保持能量不变时，可表示为四维光场。对于四维光场的表征，主要采用双平面法，表达形式为L（u，v，s，t），其中（u，v）和（s，t）分别表示光线与2个相互平行的平面的坐标交点。基于光场理论，在传统显微镜与相机传感器之间加装微透镜阵列，可组成显微光场成像系统，从而记录待测对象的光场信息。图1为显微光场成像原理示意图。其中，微透镜阵列位于筒镜后焦面处，相机传感器位于微透镜阵列后焦面处。测量体内任意一点S发出光线，经过显微镜内物镜与筒镜的作用，会聚于微透镜阵列C处，每个微透镜把接收的光线按其传播方向分散到相机传感器的不同像素点上，因此可通过单光场相机实现四维光场信息采集。

基于显微光场成像的Micro-PIV系统原理如图2所示。光场Micro-PIV 系统主要由激光器（激光波长λ0＝532 nm）、显微镜、微注射泵、光场相机系统和计算机组成。通过光场Micro-PIV系统获得时间间隔为Δt的2张示踪粒子的光场图像后，通过L-R算法反卷积重建出示踪粒子的三维坐标信息，再利用互相关算法获得瞬时三维速度场。L-R算法的迭代公式可表示为［11］

图1 显微光场成像原理示意图Fig.1 Schematic diagram ofmicroscopic light field imaging

图2 光场M icro-PIV系统原理图Fig.2 Schematic diagram of light field Micro-PIV system

从式（1）可见，显微光场成像系统的PSF是光场Micro-PIV通过L-R算法完成三维重建的前提。准确获得显微光场成像系统的PSF图像之后，反卷积问题将转换为线性复原问题。

2 显微光场成像PSF模型

点光源是最基本的光源单位，对于复杂的光源，可以将其视为众多点光源的叠加。现以点光源为例，分析其在显微光场成像系统中的成像过程。以物镜的焦点为原点，系统光轴为z轴，建立如图1所示坐标系。其中，z＝0表示物镜焦平面，z＞0表示离焦且靠近物镜，z＜0表示离焦且远离物镜。测量体内任一点光源发出的球面波在系统中的传播过程可分为以下3个部分：

1）球面波经过显微镜内物镜和筒镜到达微透镜阵列。显微镜可视为衍射受限光学系统，点光源S（x，y，z）发出的发散球面波投射到物镜上，经显微镜变换为筒镜上的会聚球面波，根据基于波动光学的标量衍射理论［14］，到达微透镜阵列的球面波的波阵面U（x，y，z）可表示为

式中：x、y、z为点光源的三维坐标；r和a分别为物面径向和轴向光学坐标；fo为物镜焦距；λ为荧光波长；M为物镜的放大倍率；α为物镜孔径角的半数，sinα＝NA，NA为物镜的数值孔径；ρ为物镜孔径的规格化径向坐标；P（ρ）为归一化光瞳函数；J0（·）为零阶贝塞尔函数；s、t为微透镜阵列面坐标。

2）球面波透过微透镜阵列。微透镜阵列可简化为疏状函数，本文采用的微透镜孔径为正方形，因此微透镜阵列的透过率函数T（x，y，z）可表示为

式中：rect（·）和comb（·）分别为矩形函数和梳状函数；fμ和D分别为微透镜的焦距和孔径。

球面波穿过微透镜阵列后，其波阵面U′（x，y，z）为

3）球面波传播到相机传感器面。球面波传播到相机传感器面可视为菲涅尔衍射，到达相机传感器面的波阵面U″（x，y，z）可表示为

U″（x，y，z）＝U′（x，y，z）*

式中：k为波矢，k＝2π／λ；u、v为传感器面坐标。

通过傅里叶变换，式（7）可表示为

式中：ξ和η分别为传感器面上坐标u和v的空间频率，ξ＝u／（λfμ），η＝v／（λfμ）；F（·）和F－1（·）分别为傅里叶变换和傅里叶逆变换。

因此，测量体内点光源S（x，y，z）发出的球面波到达相机传感器面的光强分布h（x，y，z）为

光学系统的PSF定义为输入物为一点光源时其输出像的光强分布，因此，h（x，y，z）即为光场Micro-PIV系统的PSF。

3 实验系统

为了验证本文建立的PSF模型，基于图2所示的光场Micro-PIV系统原理，搭建了光场Micro-PIV实验系统，如图3所示。其中，光场相机系统由微透镜阵列、1∶1中继镜和CCD相机组成，各部件通过笼板和笼杆加以固定。显微镜可控制载物台与物镜距离，精度为1μm，因此可确定测量体所处深度范围。显微镜内的二向色镜用于过滤被反射的激光，以提高光场图像的信噪比。显微光场成像系统参数如表1所示，改变系统参数主要通过更换显微镜物镜和相应的微透镜阵列来实现。实验所用微流控芯片为矩形层流芯片，其长直微通道截面宽度为150μm，深度为50μm。实验时，含有示踪粒子的溶液经微注射泵注入到微流控芯片中，示踪粒子被激光激发出峰值波长为584 nm的荧光，粒子激发的荧光经显微光场成像系统成像于相机传感器面，通过计算机可实时记录含有示踪粒子的光场图像。

图3 光场M icro-PIV实验系统Fig.3 Light field Micro-PIV experimental system

表1 显微光场成像系统参数Table 1 M icroscopic ligh t field im aging system param eters

4 结果分析

4.1 仿真结果

参照表1所示的显微光场成像系统参数，利用式（2）～式（9），计算不同深度的轴上点（x＝0，y＝0）的PSF，结果如图4所示。图4为深度z＝0～80μm时光轴上点的模拟PSF图像。可以发现，除焦面附近（z＝0～20μm）外，轴上点的PSF已经不再是圆斑，而是由离散的点斑组成。这是因为系统中加入了微透镜阵列，轴上的点发出的光线会被微透镜按照其传播方向折射到微透镜覆盖的不同像素点上，从而形成离散的点斑。当深度不同时，轴上点光源发出的球面波到达微透镜阵列时会有不同的波阵面，且会经过不同数量的微阵列，因此不同深度的轴上点会产生不同的PSF图像，由此可以通过PSF反演出深度信息。

图4 光轴上不同深度点的模拟PSF图像Fig.4 Simulated PSF image of points at different depth on optical axis

图5为点光源处于深度z＝50μm时不同水平位置处PSF仿真结果。图5（b）给出了位于同一深度（z＝50μm）不同水平位置的点的PSF图像，各点的具体位置如图5（a）所示，图5（c）为同一深度不同微透镜对应位置点的PSF图像。图5（a）中一个色块对应一个微透镜，边长为D／M，其中S0为轴上点，S1～S8为离轴点。通过比较图5（b）中S0～S8的PSF图像可以发现，当测量体内的点处于同一深度不同位置时，其PSF图像不一致，即PSF图像不满足平移不变性。通过比较图5（c）中S1和S′1的PSF图像可以发现，当点光源处于不同微透镜的对应位置时，如图5（a）中的S1和S′1，点光源形成的PSF图像一致。这是因为当点光源产生值为D／M的位移时，其产生的球面波到达微透镜阵列时会形成值为D的位移，而D恰为微透镜孔径的大小，因此同一深度的PSF图像满足周期分布，周期为D／M，即

图5 深度z＝50μm 时不同水平位置处PSF仿真结果Fig.5 PSF simulation results at different horizontal positions when depth z＝50μm

式中：m、n为正整数。

从式（10）可见，只需计算同一深度中心微透镜区域各点的PSF图像，再经过平移，即可获得同一深度任意点的PSF图像。通过PSF所处微透镜区域及其图像，可以反演出水平位置信息。

结合图4和图5，可见本文所建立的显微光场成像系统PSF模型可以模拟出点光源处于不同水平位置不同深度时的PSF图像，PSF反映了测量体内任意点的三维坐标信息。

4.2 模型验证

基于搭建的光场Micro-PIV实验系统，实验获取了粒径为2的示踪粒子的光场图像。深度处于z＝0～80μm 的轴上点的PSF图像如图6所示。

图6 光轴上不同深度点的实际PSF图像Fig.6 Actual PSF image of points at different depth on optical axis

从图6可以发现，当粒子深度逐渐增加时，PSF图像先呈圆斑（z＝0～20μm），然后变为离散的点斑，且PSF图像占据的微透镜阵列数随着离焦距离增大而增大，这与模拟的PSF图像表现一致。

通过比较图4和图6中对应的PSF图像，发现其形状与大小基本一致。

为了进一步评价PSF模型的准确性，本文从图像相似度与三维重建两方面进行了验证。首先通过结构相似性（SSIM）算法［15］进行模型验证，SSIM算法计算公式为

式中：SSIM（T1，T2）为图像相似度函数值，T1和T2分别为比较的2张图像矩阵；μ1和μ2分别为T1和T2的平均值；σ12为T1、T2的协方差；和分别为T1和T2的方差；c1＝（k1L）2，c2＝（k2L）2为用来维持稳定的常数，L为像素值的动态范围，k1＝0.01，k2＝0.03。相似度函数值越接近1，表明比较的2张图片相似度越高。

对图4和图6中对应的PSF图像进行图像相似度计算，结果如表2所示。可以看出，在相同条件下，模型得到的模拟PSF图像与光场Micro-PIV系统的实际PSF图像计算得到的图像相似度函数值均大于0.94，可以认为2组图像具有较好的一致性，证明了本文建立的基于波动光学的显微光场成像PSF模型的准确性。

进一步以计算获得的PSF为已知量，对实验获得的单个粒子和不同粒子浓度下的流场的光场图片，利用L-R算法进行反卷积重建，结果如图7～图9所示。图7为单个示踪粒子的实验光场图片和重建结果，示踪粒子的三维坐标为（61.5μm，61μm，50μm）。从图7（b）～（d）可见，重建粒子大致为橄榄球形，在水平和竖直方向均有一定的拉伸，水平拉伸为1μm，竖直拉伸为4.3μm，这是由于在拉伸区域内，任意点光源发出的光线经系统成像后形成的图像相似，系统无法准确识别，使得重建后出现拉伸。虽然重建存在拉伸，但重建粒子的中心点的三维坐标为（61μm，60.5μm，50.5μm），与实际三维坐标仅有一个像素的误差，在误差允许范围之内，即通过L-R算法和本文建立的PSF模型可以完成单个示踪粒子的三维坐标信息的重建。图8和图9分别为低粒子浓度（0.025 g／L）和高粒子浓度（0.25 g／L）的重建结果。低粒子浓度下，示踪粒子形成的PSF图像相互影响较小，而高粒子浓度下，示踪粒子形成的PSF图像相互影响较大，PSF图像相互叠加。但在高、低示踪粒子浓度下，L-R算法结合PSF模型，可实现示踪粒子的重建，获取示踪粒子的三维坐标信息，也进一步验证了基于波动光学的显微光场成像PSF模型的准确性。