基于经验模态分解的激光散斑噪声抑制方法

定量相位成像（QPI）［1-2］是一种新兴的无标记光学成像方法，可以实现反射样品三维形貌的可视化，测量透明和半透明样品的内部结构和折射率分布，其在透明生物样品上的应用是当前生物光学成像的研究热点之一。定量相位成像技术主要依赖全息［3］和干涉技术记录透明物体的相位信息，或是记录光场通过透明物体的衍射强度图，从记录的强度图中提取物体的相位。它可以与压缩感知［4］、深度学习［5-7］、自动对焦［8］、拓扑调制［9］或补偿技术［10］相结合，实现更高精度的测量结果。当用激光照射待测样品时，物体表面相对于激光波长量级而言相当粗糙，由于激光的高相干性，各个物点的散射光之间可能发生干涉，从而产生呈颗粒状无规则分布的散斑噪声。其存在会显著降低成像结果的分辨率、减少图像的信噪比，无法满足高精度相位重建的要求，因此有必要对激光散斑噪声进行抑制。

目前，常用的散斑噪声抑制的方法主要分成两类［11］，一类是添加或替换硬件降低光源的相干性，如使用LED等低相干性的光源［12］等。Farrokhi［13］等人采用一动一静双散射片系统来获得高速散斑场照明，以抑制定量相位成像中的相干散斑。李煊［14］等人提出一种基于旋转双散射片的全息散斑噪声抑制方法，其光学信号的非相关周期大，为获得大量非相关全息图求均值去噪提供了依据。总体来说，这类硬件方法对光程差进行精确控制的要求更高，增加了搭建光路的复杂度和难度。另外一类则是通过图像处理手段实现散斑噪声抑制。Uzan等人［15］对位于图像块中心的有噪声像素采用非局部均值（NLM）滤波，对散斑噪声起到了较好的抑制效果。但该方法在去噪的同时平滑了部分边缘处的信息，会导致相位出现误差。刘吉等人［16］提出一种基于正余弦分解结合自适应全变分的去噪方法，该方法提高了图像的峰值信噪比，减少了去噪后相位图的波动性。但该方法中自适应参数的选取是图像去噪和保留边缘信息的关键，需要通过缩小自适应参数的取值范围来达到最优去噪的重建效果。Montresor等人［17］应用残差网络（DNCNN）对数字全息干涉测量的相位图进行去噪，具有非高斯统计和非平稳特性，并表现出空间相关长度和出色的去噪性能。但该方法需要使用大量的含散斑噪声的全息图数据集对多等级神经网络进行训练，实验条件复杂。牛瑞等人［18］提出利用二维高斯窗口对包裹相位进行散斑噪声抑制，不仅保留了包裹相位的跳变边缘，还提高了相位重建的精度。但该方法中二维高斯窗口大小以及阈值会影响相位重建精度，需要研究最佳窗口大小及阈值。吴育民等人［19］提出一种基于Canny算子改进各向异性扩散（P-M）方程的全息散斑抑制方法，在去噪的同时有利于保留图像的细节信息。但该方法的再现图像仍然存在轻微的二次干涉条纹，且图像相较于原始图像而言存在一定的模糊。这类图像处理手段虽然在定量相位成像中达到了一定的散斑抑制效果，但对于相位重建精度还有待提高。

针对以上问题，本文提出了一种结合经验模态分解（EMD）的Canny算子改进P-M方程的定量相位成像散斑噪声抑制方法。首先对所记录的全息图进行EMD得到一系列频率由高到低的固有模态函数（IMF）分量，仅用高频的IMF分量重构图像以增强图像细节信息。随后引入Canny算子对细节信息突出的重构图像进行边缘检测，从而更好地控制扩散去噪过程，提高了P-M方程抑制散斑噪声的性能，并通过实验验证了该散斑噪声抑制方法的可行性。

2 基本原理

2.1　结合EMD的Canny算子改进P-M方程法在全息散斑噪声抑制中的应用

在全息图中的散斑噪声是一种乘性噪声［20-21］，其数学模型如式（1）所示：

H'(x,y)=H(x,y)N(x,y)

 ， （1）

其中：H'(x,y)为电荷耦合器件（CCD）所记录的图像，其包含着散斑噪声；H(x,y)为不含噪声的理想全息图；N(x,y)为散斑噪声。

结合EMD的Canny算子改进P-M方程法首先将图像H'(x,y)经EMD处理，得到一系列频率从高到低的IMF分量。通过这些不同频率的IMF进行选取并重构图像，从而得到效果不同的图像。比如需要对图像进行图像退化从而降低噪声的影响，则只需要选取低频的IMF分量进行重构新的图像；需要突出图像细节从而增强图像，则只需要选取高频的IMF分量进行重构图像。本文所提出的方法通过选取EMD分解出来的高频IMF分量重构出细节更加突出的图像，随后用Canny算子对重构图像进行边缘检测，从而引导P-M方程实现更好的扩散去噪性能。该散斑噪声抑制方法的主要处理过程如下：

（1）对图像H'(x,y)进行EMD分解，式（1）可以改写为：

H'(x,y)=(∑i=1nCi+r)N(x,y)

 . （2）

由于全息图的细节信息主要包含在高频分量C1中，因此只使用高频分量C1来重构新的全息图以突出细节信息，式（2）可以改写为：

H''(x,y)=C1N(x,y)

 . （3）

（2）对重构的图像H''(x,y)进行Canny边缘检测，得到边缘检测的结果C(x,y)，并记录上阈值k1。

（3）根据边缘检测结果引导各向异性扩散方程的扩散程度，进行去噪，得到散斑抑制后的全息图H'''(x,y)。

（4）随后对散斑抑制后的全息图H'''(x,y)进行数值重建，从再现像中提取出相位信息进行定量相位成像。由于此时的相位信息包裹在（-π，π）之中，需要使用相应的算法进行相位展开以获得真实的相位信息。最后根据相位误差的来源，采用相对应的相位补偿算法对误差进行补偿。

2.2　EMD

二维经验模态分解的目的是将一个二维信号f(x,y)分解成一系列频率从高到低的IMF分量，高频部分包含物体细节信息，低频部分包含物体轮廓信息，具体过程如下：

（1）找出f(x,y)的极大值和极小值，通过插值法构造出极大值曲面u(x,y)和极小值曲面v(x,y)，极大值曲面和极小值曲面的均值为：

e1(x,y)=u(x,y)+v(x,y)2

 . （4）

均值e1(x,y)与f(x,y)的差记为h1(x,y)，有：

h1(x,y)=f(x,y)−e1(x,y)

 . （5）

（2）将h1(x,y)作为新的输入重复上述过程k次，直到满足式（6）时方可停止循环。其中，SD的值一般在0.1~0.5之间，α一般为0.2。

SD=∑x=1X∑y=1Y∣∣h1(k−1)(x,y)−h1k(x,y)∣∣2h1k2(x,y)≤α

 . （6）

值得注意的是，SD的表达式实际上是两次相邻迭代计算结果之间的标准差计算公式，主要反映了两次迭代计算之间的偏差程度。当SD≤α时，说明两次相邻迭代计算结果h1(k−1)(x,y)与h1k(x,y)之间相差很小，达到了可接受的计算精度，如果再继续循环迭代，对于计算精度的提高作用不明显。此时计算得到的h1k(x,y)是一个IMF分量（C1），包含着图像f(x,y)的最高频信息。

（3）令r1(x,y)表示图像f(x,y)去除最高频信息C1之后剩余的部分，则有：

r1(x,y)=f(x,y)−C1

 . （7）

（4）令r1(x,y)作为新的待分析图像，并重复步骤（1）和步骤（2）的计算过程之后，得到第2个IMF分量C2。

（5）重复上述过程n次，当rn(x,y)和Cn小于预定误差或者rn(x,y)为单调函数时，则不能再从图像f(x,y)中提取出IMF分量。此时，图像f(x,y)如式（8）表示：

f(x,y)=∑i=1nCi+r

 . （8）

2.3　Canny算子改进P-M方程法

基于偏微分方程的图像去噪方法得到了广泛的应用，从线性均匀扩散最终发展到各向异性扩散。在各向异性扩散中，应用最为广泛的是Perona和Malik共同提出的Perona-Malik（P-M）方程［22］，许多方法都是基于P-M方程发展起来的［23-24］。

然而传统的P-M方程法存在以下问题［25］：

（1）由于梯度算子∇无法去除大的孤立噪声点，在大的噪声点处的梯度值会很大，同时扩散系数c与梯度值成反比，这会导致此处的扩散系数c的值变小，达不到扩散去噪的效果。

（2）梯度算子∇对噪声的敏感度高，抗噪性能不强，不能识别伪边缘。

为了克服上述问题，通过引入Canny边缘检测算子改进P-M方程从而克服梯度算子∇抗噪能力不强的问题［19］。其数学模型如式（9）所示：

∂I(x,y,t)∂t=div[(1−K) ⋅c(|∇I(x,y,t)|)⋅∇I(x,y,t)]

 ， （9）

其中：t为时间，I(x,y,0)为初始状态，∇为梯度算子，div为散度算子，K为Canny算子。值得注意的是，当该像素点是图像边缘时，K趋近于1，这会使得P-M方程的扩散速度减慢，有利于保留图像细节信息；否则，当K趋近于0时，会使得P-M方程的扩散速度加快，从而更好地去除噪声。c为扩散系数，用于控制P-M方程的扩散程度，常用的计算公式有两种，分别如式（10）和式（11）所示：

c(x)=exp(−xk2)

 ， （10）

c(x)=11+(xk)2

 ， （11）

在式（10）和式（11）中，k为控制扩散过程的一个系数。

使用Canny算子K改进P-M方程法的优势在于：一是Canny算子K相比于梯度算子∇对噪声的敏感度更低，抗噪能力更好。二是Canny边缘检测中的双阈值筛选一般涵盖图像70%的非边缘像素点，令扩散系数c中的参数k等于上阈值k1，显然这是一个合理的选择。该方法具体过程如下所示：

（1）对图像g(x,y)进行Canny边缘检测，得到边缘检测的结果C(x,y)，以及记录上阈值k1。

（2）按上、下、左、右4个方向求解图像g(x,y) 4个方向上的梯度，即∇N、∇S、∇E和∇W。

（3）令扩散系数c中的参数k等于上阈值k1，同时结合4个方向的梯度求解得到cN、cS、cE和cW。

（4）由于Canny边缘检测结果C(x,y)中只有0和1两种值，当C(x,y)=1时，令Canny算子K=0.01；当C(x,y)=0时，令Canny算子K=0.99。

（5）散度算子div可以按照式（12）计算：

div=14(cN×

∇N

+cS×

∇S

+cE×

∇E

+cW×

∇W)

 . （12）

（6）图像g(x,y)经过Canny算子改进P-M方程法处理后，可以表示为：

g'(x,y)=g(x,y)+(1−K)div

 . （13）

3 实验

3.1　实验装置及样本

本文搭建了一套基于Mach-Zehnder干涉方法的数字记录光路，如图1所示。将一个波长为632.8 nm的He-Ne激光器作为光源，其发出的光首先经过由一块焦距为15 mm的平凹透镜和一块焦距为250 mm的平凸透镜组成的扩束准直器，随后被分束镜BS1分成两束等光强的光。一束经过样品之后被反射镜M1反射到分束镜BS2，称之为物光；另外一束光则直接被反射镜M2反射到BS2上，称之为参考光。最后，物光和参考光在BS2处合成一束光，在CCD（BFLY-U3-23S6M-C，相机分辨率为1 920×1 200，像素大小为5.86 μm）的记录面上干涉形成全息图，并由CCD记录下来。