荧光光谱和荧光探针在检测和测量领域有着广泛的应用^[1-3]。为了使光学显微镜的分辨率达到电子显微镜的水平是一个长远的目标^[4]，提高空间分辨率的一个常用方法就是超分辨荧光显微成像中的单分子定位方法^[5]，其横向分辨率通常能够达到2~25 nm^[6]。此方法已广泛应用于二维活体组织细胞的研究，如在细胞膜中跟踪独立蛋白^{[7- 8]}。在活体细胞动力学过程的研究，如细胞内的转换、传输等过程，都是三维 (Three Dimension, 3D)的情形，所以大量科研工作者对三维的单分子定位进行了深入研究，得到一系列成果：如柱透镜引入像散 (Cylindrical Lens)^[9-10]、通过改变相位形成双螺旋点扩展函数 (Double Helix Point Spread Function, DH-PSF)^[11-13]、双焦面探测 (dDual Focus Imaging) 和荧光干涉等^[14-15]。但上述方法通常都存在一定的限制，例如要求成像分子密度不能太高。后来，科研人员以提高单帧图像里荧光分子数的识别率为目标，即荧光分子密度，又发展了一些三维高密度荧光分子定位算法，如3D-DAOSTORM（Three Dimension Dominion Astrophysical Observatory Stochastic Optical Reconstruction Microscopy）。从横向投影密度看，3D-DAOSTORM可以识别1 μm⁻²以下的荧光分子图像，其横向定位精度达到20 nm，轴向定位精度达到50 nm，当图像里的分子密度超过1 μm⁻²时，其定位精度会显著降低^[16-17]。BARSICA^[18]采用DH-PSF方法，对荧光分子图像用匹配追踪(Matching Pursuit, MP)法进行粗略定位，然后对粗略定位的结果再用凸最优化 (Convex Optimization, CO)进行精确的分子定位。此方法虽然采用了MP方法进行提速，但是由于CO算法计算速度太慢，导致算法总体上计算速度过慢，对于重构成百上千幅荧光分子图像而言，耗时太久。2014年，Lusheng Gu等人提出采用双层面成像进行三维的单分子定位。其重构算法采用的是基于MATLAB工具箱的CVX算法进行重构定位的，此方法对于重构一张大的超分辨图像来说，速度太慢^[19-20]。2016年，Christy F. Landes课题组采用DH-PSF方法进行荧光成像，并采用并行计算来缩短算法的重构时间，但该算法对电脑配置要求比较高^[21]。

针对上述问题，本文通过像散成像，即在成像光路中引入柱透镜形成像散，然后对图像采集过程建立三维压缩感知模型，并采用L1范数最小化的同伦算法（L1-Homotopy, L1-H）来重构三维超分辨图像，以提高算法的运行速度和荧光分子识别密度，减少图像采集时间和算法的重构时间，提高时空分辨率。

基于压缩感知的超分辨荧光显微成像方法，根据光学系统采用的点扩展函数来建立测量矩阵，根据成像过程建立压缩感知模型，然后采用合适的算法对模型进行重构。

在成像光路中引入一个弱的柱透镜可以在轴向产生像散，通过改变荧光分子的轴向位置，可以在探测平面形成不同形状的椭圆光斑，从而可以实现轴向定位，以及提高轴向分辨能力^[9]。包含轴向像散的单个辐射源的强度分布为^{[10, 20]}：

其中，椭圆率$\varepsilon = \sqrt {{\sigma _y}/{\sigma _x}} $，一个广义的半峰全宽(Full Width at Half Maximum, FWHM)值为$\sigma _r^2 = \sqrt {\sigma _x^2\sigma _y^2} $。式（1）中的（${\mu _x},{\mu _y}$）代表荧光分子的横向位置坐标。${\sigma _x}$和${\sigma _y}$分别表示在x和y方向强度分布的FWHM值。除了在x和y两个焦面的中间位置外（z = 0），${\sigma _x}$和${\sigma _y}$都不相等。${\sigma _x}$，${\sigma _y}$与轴向位置z有以下关系^{[20, 22]}：

其中，${\sigma _0}$表示点光源衍射极限的FWHM值，${\sigma _0} = 0.51\dfrac{\lambda }{{NA}}$，$NA$为光学系统的数值孔径，$\lambda $为辐射点源的波长。所以，对于同样波长的荧光分子，${\sigma _0}$值唯一。像散值γ由光学系统的柱透镜、管镜和物镜的焦距来确定。而${{\textit{z}}_r }$由成像系统决定^[23]：

其中，${w_{x{\textit{z}}}}$和${w_{y{\textit{z}}}}$等于XZ与YZ两个对应焦平面的FWHM值。

利用压缩感知方法分析大量重叠的荧光分子图像，在二维单分子定位中取得了一定的成功^[24-25]。在本文中，需采用三维的点扩展函数才能得到三维超分辨图像。这表示分子的轴向坐标z不同，其辐射的光子在探测器上所成的像也不同，即二维点扩展函数的形式不同。可从三维荧光分子分布与三维点扩展函数的卷积中得到探测器图像I。为便于采用三维压缩感知方法进行处理，做如下设定：探测器图像I记作矢量b、3D-PSF记作矩阵A，荧光分子三维分布记作矢量s。

如果已知光学系统参数或者采集的荧光分子图像非常稀疏，那么可以通过计算或拟合得到点扩展函数，进而可以通过点扩展函数来构造测量矩阵A。在得到高密度荧光分子图像形成的列矢量b的情况下，可以建立压缩感知模型，而荧光分子的位置分布s，作为包含了荧光位置信息的矢量，可以被重构。考虑在获得图像的过程中存在噪声，无约束优化问题可表示如下^[26]：

其中，稀疏性与求解误差之间的平衡由$\lambda $的值进行控制。本文采用L₁范数的同伦算法对其进行求解。

压缩感知算法用于三维单分子定位时，采用的是二维网格方法^[24]，如果要实现三维的单分子定位，需要运用三维网格和三维的点扩展函数，如图1所示。在STORM中，分子位置与图像间具有线性平移不变性，使用如图1所示的三维网格，4 × 4 × 10个网格点正好是单位像素对应的三维网格。80 nm × 80 nm为单个像素所具有的有效像元尺寸。将每个像素横向细分成4 × 4的网格后，横向就可以达到期望的20 nm的精度，轴向分成10层轴向就可以达到期望40 nm的定位精度。

图  1  荧光图像对应样本进行三维网格细分图。 （a）模拟生成荧光分子成像；（b）荧光图像对应的荧光分子三维空间细分；（c）一个像素和整个图像对应网格细分数

Figure 1.  3D mesh subdivision process of fluorescence image corresponding sample. (a) Simulated fluorescence molecular image; (b) the three-dimensional subdivision of sample corresponding to fluorescent image; (c) the number of the subdivision corresponding to a pixel and the whole image

可使用压缩感知的方法从一幅单一的高密度图像中恢复出粒子的三维空间位置。粒子的位置在最终的超分辨图像中不是用坐标点来表示的，而是用一些离散的三维网格来表示的。相应坐标轴方向所期望的精度必须大于网格的轴向距离与横向距离。利用2.1节中的方法生成模拟图像（图1(a)），假设图像大小为46 pixel × 46 pixel，分成184 × 184 × 10个网格点（图1（b））。使用的三维网格如图1(c)所示。由三维网格点[s]_mp依次连接构成一维矢量[s_j], 图1(b)给出了三维网格点的连接顺序。在矢量s中，坐标网格中分子发射的光子数是用元素值的大小表示的，当网格点中没有分子时元素值大小等于零。利用系统的三维点扩展函数计算矩阵[A_ij]，当且仅当只有一个分子在位置s_j产生的图像排列成列向量时，可得到矩阵A 的第j^th列。许多荧光粒子产生的图像排列而成的列向量可构成矢量b。当然，也可通过已知三维点扩展函数恢复分子的三维空间位置坐标，即使荧光分子形成的光斑图像出现重叠实现。如果直接对整幅图像进行处理，需要对整个样品进行三维空间划分，这样会导致由三维点扩展函数形成的观测矩阵比较大，从而导致重构时间长，矩阵存储量较大。因此，考虑到点扩展函数的展宽比较窄，对原始图像进行分块处理不会影响定位精度，具体步骤如下：

（1）探测图像I大小为46 pixel × 46 pixel。超分辨图像为S，184 × 184 × 10个网格点。定义R为对横向像素进行细分的比例，R = 4，轴向的分层数L为10，见图1(c)。

（2）考虑重构时间和点扩展函数所占像素的情况，将原图像分成大小为16 pixel×16 pixel的小图像块I₀（图2(a)），处理后的小块图像，相邻块之间共享6个像素，这是由系统的点扩展函数的半高全宽决定的，如图2(b)所示。

图  2  对整个图像进行分块处理。（a）从一帧图像里取一小块；（b）重合边缘；（c）对取出的小块图像进行边缘扩展；（d）对重构结果保留有效区域

Figure 2.  The original image subdivided into smaller patches. (a) A small patch taken from one frame of image； (b) superposition edges; (c) edge expansion of the small block of image; (d) the retaining valid areas of reconstruction results

（3）对于每一小块I₀，横向细分(16 × 4) × (16 × 4)，分块边缘外的样本分子决定了边缘的像素强度值，所以在细分后的网格点周围加上6个网格点，网格点总数为(16 × 4 + 12) × (16 × 4 + 12) × 10，如图2(c)所示。

（4）对小块进行处理，使得分块图像I₀的每列与超分辨图像S的每列依次连接同一个列向量，I₀为256 × 1，s为57760 × 1；s为仅第j个元素有分子发光时得到的图像，即为点扩展函数矩阵A的第j列，A大小为256 × 57 760。利用已知的I₀与A得到s，然后把s转换成超分辨图像。为了避免相邻块边缘对每个小块的影响，使其外围18个网格点为零。这18个网格点是由12个相邻块重叠的网格和在第（3）步中考虑边缘外影响的6个网格点组成的。如图2(d)所示。

（5）最后组合而成的一幅三维超分辨图像是由每一小块重构出的三维矩阵而得到的。

通过分析一系列高密度荧光分子图像可知道三维压缩感知重构方法的性能。在本节中，一个样本体积为2.76 μm × 2.76 μm × 0.4 μm，分子随机分布在其中。虽然本次模拟生物样品时取的样本体积比较小，但不会影响算法的性能分析，因为文中对图像进行分块处理。在处理整幅图像时，对每次处理的一个小块都做了分块操作，所以精度并不会随着体积横向投影的增加而衰减。又由于超出轴向范围的都当作背景处理了，所以轴向大小会对精度有一定的影响。相机的有效像元大小为80 nm，设v为实际探测器的尺寸， M为探测光路的放大率，则其有效像元尺寸为v/M。通过一个100倍物镜和实际像元尺寸为8 μm EMCCD的组合可以实现80 nm × 80 nm的有效像元尺寸。将点扩展函数的半高全宽τ₀设置为200 nm，相当于2.5个像素的大小。由衍射极限公式知，其数值孔径大小为1.4，波长λ为560 nm。设置焦深z_r为500 nm，设置像散参数γ为200 nm。根据大量荧光分子图像的分析结果，设置每个荧光分子所发射的光子数服从对数正态分布， 峰值为5000， 标准差为1700。探测器所能接收的光子数由荧光分子的曝光时间、辐射效率以及探测效率共同决定。实际上，信噪比越好代表采集的光子数越多，这就意味着拟合结果以及定位精度更好。将坐标原点设置为样品中心。在模拟生成荧光分子图像时，让32个分子(横向投影密度为4 μm⁻²) 随机分布在三维空间内，如图3（彩图见期刊电子版）所示。给图像添加一个高斯噪声（均值为零，方差为2个光子）和50个光子数的均匀背景噪声。为了使模拟结果更接近实际情况，在整个图像上添加泊松化的噪声。添加噪声后，图像的峰值信噪比 (PSNR) 为24.2746 dB，由公式（1）产生的图像如图3中底面所示，图像的尺寸应略大于样品的成像区域，以避免图像边缘对重构算法的影响。

图  3  模拟生成的荧光分子图像（绿色星号为三维空间随机分布的荧光分子；三维空间的底面图像为荧光分子模拟生成的图像）

Figure 3.  Simulated fluorescent molecule image (green asterisks are randomly distributed fluorescent molecules in 3D space; the bottom image in 3D space is the image generated by the simulated fluorescent molecules)

为了验证基于压缩感知的L1-H算法的有效性，通过比较CVX方法，OMP算法与L1-H三种算法的性能，模拟生成不同荧光分子的密度图像，其密度从0.5到16变化，每张图像大小为46 pixel × 46 pixel，有效像元大小设置为80 nm。为了定量分析几种算法性能，采用模拟生成的不同方法对生成的荧光分子图像进行重构，分别从恢复率，横向与轴向定位精度，算法运行时间三方面对3种算法进行分析（图4）。恢复率和定位精度的定义与计算方法可参考文献[27]。从图4可以看出：随着分子密度的增加，恢复率会降低，但密度小于16 时，3种算法的恢复率都大于70%，并且比较接近（图4(a)）；随着荧光分子密度的增加，横向定位精度与轴向定位精度会变差，CVX方法和L1-H两种方法的定位性能相近，OMP算法相对其它两种算法，定位精度差很多（图4(b)，4(c)）；算法的运行时间是分析重构一张46 pixel × 46 pixel的荧光分子图像所用的时间，从图4(d)可以看出，采用OMP算法和L1-H算法的运行速度比CVX方法快了两个量级。综合考虑，L1-H算法最适合用于三维荧光分子的定位。

图  4  CVX，OMP和L1-H 三种算法对模拟生成的荧光分子图像的重构结果。（a）恢复率；（b）横向XY用标准差（Stdev）显示定位精度；（c）Z向用标准差显示定位精度；（d）算法运行时间

Figure 4.  Reconstructed results of simulated fluorescence molecule image by three kinds of algorithms. (a) Recovery rate; (b) localization accuracy of horizontal XY shown with the standard deviation; (c) localization accuracy of axial Z shown with the standard deviation; (d) algorithm running time

本文对有结构的生物样品和实际生物样品进行了超分辨图像重构，以验证算法的实用性。通过上文分析知，L1-H与CS-CVX的重构密度和定位精度相近，都比OMP算法好，而且L1-H的执行时间比CS-CVX算法快近两个量级，所以接下来仅采用L1-H算法进行分析。

根据荧光显微成像原理，采用模拟软件使八条射线束上随机分布荧光分子，并且使八条射线逆时针分布在不同的三维空间内，如图5（彩图见期刊电子版）所示。轴向范围为−200~200 nm (图5(a))，分子处在不同的轴向位置时，其颜色不同，图5(a)中的颜色条表示分子的轴向位置。为准确对比清晰度，将−250 nm设置为所有的背景值，黑色为图像的背景颜色。为实现荧光分子的稀疏分布，对三维空间内的荧光分子进行随机激发。图5(b)即为经过光学系统后所成的像，其中白色比例尺为800 nm，随机分布的图像里分子数为54。为了让射线束连续，随后模拟随机生成200帧不同分子分布的图像。然后，进行图像叠加，这200帧图像的分子叠加结果如图5(a)所示。图5(c)为直接进行叠加的宽场图像。对每张图像进行重构，图5(d)是把所有的结果进行叠加而得到的三维超分辨图像。其轴向范围由图5(a)，5(d)中的颜色条所表示，单位为nm。光子数由图5(b)，5(c)的颜色标尺表示。从8条射线形状可以看出横向定位基本正确，从颜色可以看出轴向定位结果非常接近真实结果，由此验证了算法的实用性。

图  5  模拟生物实验的三维重构结果，图中的比例尺为800 nm。 （a）模拟生成随机分布的荧光分子累加结果，横向范围为6.9 μm × 6.9 μm，轴向范围为−200~200 nm；（b）随机生成一帧图像；（c）随机生成200幅图像进行累加的结果；（d）用算法分别对200帧图像进行重构定位出分子整合成的一张三维超分辨图像。

Figure 5.  3D reconstructed results of simulated bioexperiment (scale bar is 800 nm). (a) Accumulation results of randomly distributed fluorescent molecules with a transverse range of 6.9 μm×6.9 μm and an axial range from −200 nm to 200 nm. (b) A randomly generated frame of image. (c) Accumulation results of randomly generated 200 frames of images. (d) A 3D super-resolution image obtained by reconstructing 200 frames of images.

通过上述的模拟实验可知,超分辨成像可以应用于模拟的生物样品中。接下来，选取实际的生物样品对此算法进行验证。选取非洲绿猴肾细胞进行了荧光显微成像实验，然后对采集到的非洲绿猴肾细胞荧光图像进行重构处理，如图6（彩图见期刊电子版）所示。宽场图像是由2 000帧图像直接累加而成，在图6(a)右边的颜色标尺表示归一化的光子数；对每帧荧光图像用L1-H算法进行定位，然后对所有帧定位之后的图像进行累加，最后获得的三维超分辨图像如图6(b)所示。图右边的颜色标尺表示分子所处在的轴向深度。从图6(a)，6(b)分别取出图6(a)白线位置对应的像素值，然后对其进行归一化，如图7所示。从图可以看出横向分辨率有了很大的提高。从而验证此算法具有一定的实用性。

图  6  三维超分辨图像的重构结果。（a）由2 000帧荧光图像累加得到的宽场图像；（b）本文算法重构的三维超分辨图像

Figure 6.  Reconstructed results of 3D super-resolution image (scale bar is 2 μm). (a) Wide-field image obtained by accumulating 2 000 frames of fluorescent images. (b) 3D super-resolution image reconstructed by the proposed algorithm.

图  7  图6（a）中白线所在位置的宽场图像与对应位置在图6（b）中的重构结果

Figure 7.  Comparison of reconstructed results and wide-field results along the white line in Fig.6

本文提出将L1范数最小化的同伦算法应用于三维分子定位中，并在传统的荧光显微镜中引入柱透镜，形成三维的点扩展函数，然后建立压缩感知模型。此模型可以应用于各种三维点扩展函数中，其随压缩感知模型中观测矩阵的变化而变化。根据模拟和实验结果可得出以下结论：基于L1-H的三维压缩感知方法比采用MATLAB 工具箱CVX方法的计算速度快了两个数量级，并且其定位精度优于OMP算法。L1-H算法可对高密度的STORM数据实现三维的超分辨成像，同时还可以实现纳米精度的快速定位。