引　言

在我们的生活中，粒子无处不在，如沙土、粉尘、雨滴、血细胞等。在生物医学领域，红细胞的计数是评估患者健康状况的关键指标，细胞的浓度对于精确控制所需化学物质的用量至关重要^[1]。工业生产中，通过检测燃烧后产生的颗粒尺寸，可以判断燃料是否完全燃烧^[2]。社会生活中，基于光散射的粒度分析可以改善许多行业产品的质量控制，包括制药、食品、化妆品等^[3]。因此，对粒子的研究具有重大的科学意义及广阔的应用前景。

已提出的光学粒子测量方法主要有激光多普勒测量法、激光衍射测量法、数字全息法等^[4]。其中，数字全息技术具有记录大景深、三维信息获取、可实时处理和分析等优点，在光学粒子测量中具有重要的应用潜力。数字全息法通过对采集到的全息图在不同传播距离进行数字重建，获得粒子的三维位置。然而，当待测粒子位于散射介质中，由于光线在介质中的散射效应，会导致获得的全息图信噪比（signal-to-noise ratio，SNR）严重降低。在此场景下，数字全息反向传播（back propagation，BP）数值重建将遇到困难，给粒子三维信息的获取带来挑战^[5]。

近年来，压缩感知（compressed sensing，CS）^[6]理论的提出为解决数字全息中的低信噪比问题带来了新的机遇。该理论突破了奈奎斯特采样定理的限制，通过将全息技术与压缩感知相结合。利用CS优化算法重构原始信号，能够精确重建三维粒子场信息^[7]。基于全变分（total variation，TV）^[8]正则化的CS方法为降低全息图中的散射和孪生像噪声提供了有效的解决方案。CS方法已经证明了其在非均匀介质成像中抑制散射噪声的能力，并已成功应用于稀疏性约束下二维全息图重建三维分层物体^[9]。2012年，美国麻省理工大学的Liu等^[10]将压缩全息引入边缘稀疏约束，以亚像素精度对物体进行定量定位测量，实现了1/45像素的定位精度。2017年，美国西北大学的Wang等^[11]将压缩全息术和编码曝光技术结合起来，通过数字微镜实现了像素级的时间调制，根据单幅振动头发全息图来重建其四维运动场景。2020年，上海大学的伍小燕等^[12]将变密度减采样的方法与压缩感知重建相结合，实现了较低采样率条件下高质量层析重建。压缩全息术的出现为散射介质中粒子信息的恢复提供了一种可行的方法，提高了粒子检测和分析的准确性。

本研究利用压缩全息和自动聚焦的方法，恢复了清晰的散射介质中粒子全息图，并定位了粒子的三维位置。该方法的轴向精度达到3 μm，有效地解决了数字全息图重建中的逆问题，降低了杂散光噪声，减弱了孪生像干扰，并抑制了散射介质中物体信息的散射。

1 原　理

同轴数字全息是一种常用的三维粒子场成像方法，具有结构简单，抗干扰能力强等优点^[13]。图1为同轴全息成像原理图，其中相干光束穿过粒子样本，与未穿过样本的参考光束发生干涉，形成全息图，并在成像面被电荷耦合器件（charge coupled device，CCD）所捕捉^[14]。

参考光和物光形成干涉，在CCD平面的干涉条纹光强I可表示为

$\begin{split} I\left(x,y\right)=&{\left|O\left(x,y\right)+R\left(x,y\right)\right|}^{2}\\ =&{\left|O\left(x,y\right)\right|}^{2}+{\left|R\left(x,y\right)\right|}^{2}+\\ & {R}^{*}\left(x,y\right)O\left(x,y\right)+R\left(x,y\right){O}^{*}\left(x,y\right) \end{split}$

(1)

式中：O（x, y）表示物光波；$ R\left(x，y\right) $表示参考光；R*（x, y）表示参考光波的共轭；O*（x, y）表示物光波的共轭。激光穿过粒子场后，经过粒子的调制，其强度信息包含在O（x, y）中，通过数值方法求解强度信息，得到粒子的信息。

采用角谱数值重建方法，在重建距离z处，物光波复振幅分布为

$ \begin{array}{c}{u}_{ {\textit{z}}}\left(x,y\right)={F}^{-1}\left\{F\left[{u}_{0}\left(x,y\right)\right]H\left(\xi ,\gamma \right)\right\} \end{array} $

(2)

式中：F表示傅里叶变换；$ {F}^{-1} $表示傅里叶逆变换；$ {u}_{0}\left(x,y\right) $代表z = 0处的物光波复振幅；$ H\left(\xi ,\gamma \right) $为传递函数，其表达式为

$ \begin{array}{c}H\left(\xi ,\gamma \right)=\mathrm{exp}\left[\mathrm{j}\dfrac{2{\text{π}}{\textit{z}}}{\lambda }\sqrt{1-{\left(\lambda \xi \right)}^{2}-{\left(\lambda \gamma \right)}^{2}}\right] \end{array} $

(3)

式中：$ \lambda $为光波波长；频域$ \xi $与空域$ \gamma $的坐标对应为

$ \begin{array}{c}\xi =\dfrac{x}{\lambda {\textit{z}}}\text{，}\gamma =\dfrac{y}{\lambda {\textit{z}}} \end{array} $

(4)

其中，式（1）归一化后仍然具有适用性，可得

$ \begin{array}{c}H\left(x,y\right)=2{\mathrm{Re}}\left\{O\left(x,y\right)\right\}+e\left(x,y\right) \end{array} $

(5)

式中，e（x, y）表示系统引入的误差。图像大部分区域内的梯度分布通常趋向于零，这表明在梯度域内物体可以被认为是稀疏的，而孪生像产生的扩散模式在梯度域非稀疏。压缩全息通过稀疏性约束集成了数字全息和压缩传感，主要求解方法有两步：迭代收缩阈值算法（two-step iterative shrinkage/thresholding algorithm, TwIST）和阈值迭代算法（iterative shrinkage/thresholding algorithms，IST）。TwIST算法相对于IST算法具有更快的收敛速度^[15]，因此本研究利用TwIST算法来重建图像抑制噪声。重建平面中物体的估计复场为

$ \begin{array}{c}\widehat{f}={\mathrm{argmin}}\dfrac{1}{\gamma }\parallel g-2{\mathrm{Re}}\left\{Hf\right\}{\parallel }_{2}+k{{\varGamma} }_{\mathrm{T}\mathrm{V}}\left(f\right) \end{array} $

(6)

式中：g表示测量图像之间的误差；k>0为常数；Γ_TV（f）为全变分（total variation，TV）正则化算子^[16]，其计算公式为

$ \begin{array}{c}\varGamma_{T V}(f)=\iint|\nabla f| \mathrm{d} x \mathrm{~d} y=\iint (f_x+f_y) {\mathrm{d}} x {\mathrm{d}} y\end{array} $

(7)

式中：f_x，f_y分别表示图像在x，y方向上的复振幅值差的绝对值。TV值的含义为图像梯度，该值越小则图像越平滑，噪声越低。因此，通过限制图像的全变分可以达到降低噪声的效果。

为了从不同深度的重建中自动检测粒子的三维位置，使用了Intensity（I），Laplace（L），Brenner（B），Roberts（R）^[17]等自动对焦评价函数来确定粒子的深度。它们各自的方程为

$ \begin{array}{c}I={\displaystyle\sum }_{x}{\displaystyle\sum} _{y} \left\{\left[f\right(x,y)-\mu {]}^{2}\right\} \end{array} $

(8)

$ \begin{array}{c}L={\displaystyle\sum }_{x}{\displaystyle\sum} _{y}{\left[ \begin{array}{c}q\left(x + 1,y\right) + q\left(x - 1,y\right) + \\q\left(x,y + 1\right) + q(x,y - 1) - 4q(x,y{)}^{2}\end{array} \right]}^{2} \end{array} $

(9)

$ \begin{array}{c}B={\displaystyle\sum} _{x}{\displaystyle\sum} _{y}\left\{{\left[q\left(x+2,y\right)-q\left(x,y\right)\right]}^{2}\right\} \end{array} $

(10)

$ \begin{array}{c}R={\displaystyle\sum} _{x}{\displaystyle\sum} _{y} \left\{\begin{array}{c}{\left[q\left(x+1,y+1\right)-q\left(x,y\right)\right]}^{2}+\\ {\left[q\left(x+1,y\right)-q\left(x,y+1\right)\right]}^{2}\end{array}\right\} \end{array} $

(11)

式中：$ \mu $表示图片的平均强度；q（x, y）表示图像中坐标（x, y）处的像素强度值。

2 模拟结果和分析

当粒子处于散射环境时，使用同轴全息采集的粒子衍射图像信噪比很低，这给传统的数值重建带来了挑战。通过模拟生成一幅包含5个不同深度粒子的数字全息图，并加入阈值为0.4的高斯噪声。噪声将模拟粒子全息图的衍射环淹没，这模拟了捕捉散射介质中的粒子的成像结果。图2（a）和图2（b）显示了模拟粒子在空间上的分布。颗粒深度分别为462 μm，526 μm，582 μm（2个），614 μm和734 μm。

图3展示了模拟粒子在不同深度中的重建结果，选取其中3个聚焦平面，分别采用角谱法反向传播（BP法）和CS法对模拟粒子进行重构。在模拟粒子深度前后5 μm处对粒子进行重建，间隔为1 μm。图3（a）～（c）为用CS法重建的模拟粒子全息图。CS中总变分正则化条件的约束有效地减弱了模拟粒子重构中固有的零阶项和孪生图像噪声。当一个平面上的粒子聚焦时，其他平面上的粒子就会散焦。压缩全息法还可以有效地抑制不同层间粒子的干涉。使用MATLAB中的imfindcircle函数，识别出边缘锐利的粒子，并在图中用红色圆圈表示。应用上文说明的4种自动对焦方法对角谱法重建后的图像中的粒子进行定位，自动对焦结果如图3（d）～（f）所示。对图3（b）中的粒子，4种方法并不能找到极小值点，其自动聚焦曲线单调递增，说明在背景噪声较强的情况下，角谱法的重构不能有效地恢复粒子信息。图3（g）～（i）为CS重建全息图的自动对焦曲线。可以看出，4种方法能够有效地识别出共同聚焦平面。红色虚线表示设置的粒子深度位置。4种自动聚焦方法的极值点与设定的粒子深度位置误差均小于3 μm，5种粒子的平均深度误差为2.2 μm，证明了CS方法在噪声环境下处理粒子全息图能够获得较高的轴向精度。

3 实验结果和分析

为进一步验证本研究方法的有效性，进行了水和明胶介质中的粒子识别实验，采用如图4（a）所示的同轴数字全息系统进行数据采集。实验系统如图4（b）所示，使用的激光器为富和科技的632 nm红色光源，型号为FU650ADX-PXG1670，输出功率为5 mW，相干长度为20 cm。全息图使用奥林巴斯倒置显微镜，用10×物镜放大，用分辨率为1280×1024，像元尺寸为4.8 μm的Wattentech HTSUA34GC相机记录。样品为聚苯乙烯微粒，直径为30 μm，分别分散于清水和明胶当中，并放置于载玻片上方便观察。含有聚苯乙烯微粒的明胶实验样本如图4（c）所示。

相机记录的粒子全息图如图5所示。图5（a）为水中粒子全息图；图5（b）为明胶中的粒子全息图。水中的颗粒表现出清晰的衍射环，而明胶中的颗粒由于介质的散射而在衍射环内表现出散射效应，这种散射效应会给后续的重建带来严重挑战。

图6（a）为角谱法重建的水中粒子，图6（b）为CS法重构的水中粒子。水中粒子的实验数据与模拟结果接近，重建质量高，有效地减弱了背景噪声和孪生像。图6（c）描绘了用角谱法重建的明胶中的颗粒。很明显，使用角谱法从明胶中恢复粒子是具有挑战性的。图6（d）为CS重建的明胶中的颗粒，可以看出明胶中的颗粒结构清晰可见，明显减少了明胶介质造成的散射。红色圆圈表示通过MATLAB中的imfindcircle函数识别的该平面内的聚焦粒子。

表1给出了CS和BP 重构的水及明胶中颗粒的峰值信噪比（PSNR）、结构相似指数（SSIM） 、均方误差（MSE） 以及信噪比（SNR）。信噪比的值V_SNR =（I_avg/I_o），其中I_avg 为粒子图像的平均强度，I_o 为粒子图像的背景强度。MSE 是用于测量两幅图像的像素值之间差异的平均平方误差。MSE的值越小，两幅图像之间的差异就越小，图像质量就越好。PSNR 用来衡量两幅图像之间的信号强度与噪声强度之比，该值越高，图像质量越好。SSIM 用于评估两幅图像之间的结构相似性，主要考虑3个因素：亮度、对比度和结构。与PSNR 相比，SSIM 能更好地反映人眼对图像质量的感知，该值越大，两幅图像之间的相似性越好。4种评价指标都通过将背景图像作为参考来计算。表1中4种评价指标比较的结果显示，CS 算法处理的颗粒，无论是在水中还是在明胶中，都得到显著改善。这表明CS 与TV 正则化相结合可显著提高重建图像的质量。

为了验证对散射介质中粒子大小的测量结果，使用Hough变换识别明胶中不同直径的颗粒数量，通过设置Hough变换的直径，统计出不同直径粒子的分布，结果如图7所示。实验中添加了20个直径为30 μm的粒子，CS法识别的直径为30 μm的粒子占全部粒子数的75%，而角谱法识别的30 μm直径的粒子只占30%。说明光线在介质中散射以及背景噪声会导致角谱法难以准确识别明胶散射介质中的粒子大小，而通过CS方法可以很好地解决这个问题，从而验证了所提出的方法的有效性。

4 结　论

本文介绍了一种基于压缩全息的检测方法，它能够用于观察和分析散射介质中微粒的三维分布，判断粒子直径。将CS与TV正则化相结合，能够有效抑制杂散光噪声，消除孪生图像干扰，使重建图像具有较高的轴向精度。实验结果表明，本文提出的方法在处理明胶中粒子全息图时，显著提高了图像质量，在PSNR，SSIM，MSE和SNR等评价指标上均有提升。与传统的角谱重建方法相比，CS处理后的粒子图像质量明显得到提高。然而，由于明胶介质的独特性质，在某些深度位置的粒子识别仍然存在挑战，表明该方法需要进一步改进。未来的工作将集中在优化算法以提高散射介质中粒子识别的准确性，并探索该技术在检测其他生物材料方面的潜在应用。