1 引言

折射率和直径是能表征微小粒子特性的重要参数，高精度测量微小粒子的折射率和直径有许多应用，如喷雾液滴蒸发、燃烧等过程的测量^[1-2]。彩虹测量技术由于其非接触、高精度、快速、大信息量、可定量计算等优点，在测量微小粒子的折射率和直径中有广泛的应用^[3-4]。随着彩虹技术的发展，小粒子的形状不再局限于球形，彩虹法同样适用于扁圆、椭球粒子的测量^[5-7]。此外，彩虹技术也可测量液柱或液滴折射率梯度和温度梯度^[8]，或测量喷雾场群粒子平均直径和平均折射率^[9-10]。在提高彩虹测量精度方面，有许多研究在测量方法及计算方法的优化上。通过计算折射率或直径的微小变化量与彩虹Ripple结构相位偏移之间的关系，可对折射率或直径的相对变化进行测量^[11]；利用高阶彩虹条纹对于椭球度的敏感度非常高的特征，通过高阶彩虹可以更准确反演扁圆液滴的椭球度^[12]；通过研究高斯光束入射下的彩虹分布规律，采用合适的高斯光束腰大小和光束位置可提高信号强度，提高反演精度^[13]；结合最小二乘优化算法和彩虹技术反演算法，可对粒径分布及折射率实现高精度测量^[14]。

经验模态分解(empirical mode decomposition, EMD)，是一种非线性、非平稳信号分解方法^[15]。这种方法无须预设任何基函数，本质是通过数据的特征时间尺度或特征空间尺度，获得信号的本征波动模式，然后再分解数据，能将一个频率不规则的信号波分解为多个单一频率的波，和一个平均趋势项波的形式。广泛应用于分析非线性、非平稳信号序列^[16-17]。本文结合经验模态分解方法与液滴一阶彩虹强度分布理论，研究了一种新的针对彩虹测量的反演算法，并用于液滴的折射率与直径的同时测量。本方法简化了测量过程，具有去噪声能力强、计算用时少，折射率与直径的测量精度高的特征。

1 彩虹技术原理

一阶彩虹强度分布主要包含Airy结构和Ripple结构，图1为折射率m=1.332 4，直径D=700 μm的一阶彩虹强度分布。理论和实验研究表明，彩虹强度分布与入射光波长、液滴折射率和液滴直径有关^[18]。

图 1. Distribution of first-order rainbow scattering intensity

Fig. 1. Distribution of first-order rainbow scattering intensity

下载图片 查看所有图片

由几何光学理论，几何光学彩虹角θ_rg与折射率m之间的关系为

1${\theta _{{\rm{rg}}}} = 4{\rm{arccos}}\sqrt {\frac{{4 - {m^2}}}{{3{m^2}}}} - 2{\rm{arccos}}\sqrt {\frac{{4 - {m^2}}}{3}} $

根据彩虹的Airy理论，彩虹角附近散射光虚拟波前的波动方程为

2$\Omega ({\textit{z}}) = \int_0^\infty {\cos \frac{1}{2}} {\text{π}} ({\textit{z}}t - {t^3}){\rm{d}}t$

式中z为无量纲参量，表示与彩虹角的归一化角偏差，可表示为

3${\textit{z}} = - (\theta - {\theta _{{\rm{rg}}}}){\Bigg(\frac{{16{D^2}}}{{h{\lambda ^2}}}\Bigg)^{1/3}}$

式中h为彩虹三次波前曲率，表示为

4$h = \frac{9}{4}\frac{1}{{{m^2} - 1}} \cdot {\Bigg(\frac{{4 - {m^2}}}{{{m^2} - 1}}\Bigg)^{1/2}}$

几何彩虹角位于z=0处，由(2)、(3)式可求出几何彩虹角处的光强为最大光强（一级Airy峰值）的43.93%。

根据彩虹Airy理论，第一个Airy峰的角位置与折射率和直径之间的关系表示为^[19]

5${\theta _1} = {\theta _{{\rm{rg}}}} + 1.087\;376{\Bigg(\frac{{h{{\text{π}} ^2}}}{{12}}\Bigg)^{1/3}}{\alpha ^{ - 2/3}}$

式中α为无量纲直径参量α=πD/λ。根据图2中的几何关系，彩虹Ripple结构的频率F_Ripple，可通过计算对应于彩虹角时无限远处外反射光线与内反射光线的光程差获得^[20]：

6${F_{{\rm{Ripple}}}} = \frac{D}{{2\lambda }}\Bigg(\cos {i_{{\rm{rg}}}} + \cos \frac{{{\theta _{{\rm{rg}}}}}}{2}\Bigg)\frac{{\text{π}} }{{180}}$

图 2. The formation principle of rainbow Ripple structure

Fig. 2. The formation principle of rainbow Ripple structure

下载图片 查看所有图片

2 反演算法研究

2.1 彩虹信号的EMD分析

EMD是一种尺度分离算法，用I(θ)表示信号强度，首先确定I(θ)的所有极值点，再用三次样条曲线构造由极大值点组成的上包络线u(θ)和由极小值点组成的下包络线v(θ)，计算上下包络线均值序列：

7$p(\theta ) = \frac{{u(\theta ) + v(\theta )}}{2}$

求原信号序列I(θ)与均值序列p(θ)的差值序列：

8$q(\theta ) = I(\theta ) - p(\theta )$

通过判断q(θ)是否满足：1） 极值点数量与过零点数量相等或最多只差一个；2） 在任意点q(θ)的上下包络线均值为0。若满足则q(θ)为分解出的第一个高频本征模态函数(intrinsic mode function, IMF)，若不满足，对q(θ)重复取上下包络线均值序列，再计算差值序列，直至差值序列满足停止条件，满足停止条件的差值序列为分解出的第一个高频模态IMF₁。用I(θ)减去IMF₁，对剩余分量继续重复以上过程，直到分离出全部模态序列。信号最终可表示为

9$I(\theta ) = \sum\limits_{i = 1}^n {IM{F_i}} (\theta ) + R(\theta )$

式中:IMF_i(θ)表示第i个本征模态函数；R(θ)表示彩虹强度总体趋势。

图3为模拟m=1.332 4，D=700 μm彩虹信号的EMD分析结果，IMF₁和IMF₂代表彩虹的Ripple分量，IMF₃代表彩虹的Airy分量，最后一项R代表了彩虹信号的平均趋势。从EMD分析结果中可有效分离Airy分量和Ripple分量，并可得到第一个Airy峰的角位置θ₁。

图 3. EMD of rainbow signal

Fig. 3. EMD of rainbow signal

下载图片 查看所有图片

实测彩虹信号含有噪声，对同一液滴连续采集的2帧彩虹信号做EMD分析，求2组EMD结果对应模态间的互相关系数，由于噪声的随机性，代表噪声分量的对应模态不相关，代表彩虹有效分量的对应模态应高度相关，可由此判断彩虹信号的EMD结果中的噪声模态，对彩虹信号进行去噪。互相关系数可表示为

10$r = \left| {\frac{{\displaystyle\sum {XY} - \frac{{\displaystyle\sum X \sum Y }}{N}}}{{(\displaystyle\sum {{X^2}} - \frac{{{{(\displaystyle\sum X )}^2}}}{N})(\displaystyle\sum {{Y^2}} - \frac{{{{(\displaystyle\sum Y )}^2}}}{N})}}} \right|$

式中：N为数据点个数；X和Y分别代表2帧采集数据的对应模态序列。当r>0.8时判定该模态高度相关，即该模态为有效信号；当r<0.3时，判定为不相关，即该模态为噪声分量。

为了验证上述去噪方法，理论模拟了m=1.332 4，D=700 μm的液滴的一阶彩虹信号，并对该信号加上2组随机噪声，得到2组加了不同随机噪声的彩虹信号，如图4所示。

图 4. Rainbow signal

Fig. 4. Rainbow signal

下载图片 查看所有图片

图5为图4中2组加入随机噪声后的彩虹信号的EMD结果。2组EMD结果对应模态间的互相关系数计算结果如表1所示，由表1中数据分析，IMF₁为噪声分量，IMF₂和IMF₃为Ripple结构分量，IMF₄为Airy结构分量，R反映光强趋势。

图 5. EMD results of rainbow signal

Fig. 5. EMD results of rainbow signal

下载图片 查看所有图片

结合表1分析结果，对图5(a)加入随机噪声的彩虹信号减去IMF₁，得到滤噪结果，滤噪结果与原信号相关系数为0.82，说明了彩虹信号EMD去噪技术的有效性。

2.2 反演算法的设计

折射率在1.3~1.5，直径在100 μm~2 000 μm范围内，根据彩虹分布在z=0处强度与最大峰值强度间的关系，可通过在一级Airy峰小角度侧搜索强度最接近最大强度的43.93%的角位置，代入(3)式，得到折射率初值m₀。利用F_Ripple对折射率和直径的变化表现出的不同强度的变化规律将m₀代入(6)式，得到直径D₁。再由θ₁对折射率和直径的变化表现出的不同强度的变化规律，将D₁代入(5)式，得到折射率m₁。通过反复将得到的折射率和直径分别代入(6)式和(5)式不断修正结果，再通过求解(11)式，得到多次修正结果的最优值。

11$f\left( {m,D} \right) = \min \left\{ {\sum\nolimits_{\theta = 1}^N {{{({I_{0\theta }} - {I_{k\theta }})}^2}} } \right\}$

式中：I_0θ为实测彩虹信号；I_kθ为每组修正结果模拟的理论彩虹信号。

3 实验结果与分析

3.1 测量装置及测量方法

实验系统如图6所示。由注射泵（LINZ-8A，雷恩，上海）定量控制压出液体，悬挂在光纤末端形成球形液滴，由激光器（λ=632.8 nm，P=2 mW）发出的平行光束经平面镜反射后照射液滴，经液滴散射形成彩虹条纹，用CCD（SG-14-01k80-00-R，1 024像素）记录彩虹分布。显微相机（LP-530HS，倍特嘉，苏州）用来监控液滴和辅助测量。本实验在光学暗室进行，实验温度22.4 ℃。

图 6. The experimental diagram

Fig. 6. The experimental diagram

下载图片 查看所有图片

通过移动和旋转平面镜，角度精度为0.017°，使CCD恰好能接收到一阶彩虹信号，旋转角度为φ，CCD单个像素大小为14 μm，与液滴中心距离为30.4 cm，可确定CCD测量角宽度Δθ=2.704°，由φ和Δθ可确定测量的散射角范围。CCD每次采集512帧数据。实验液滴样品为浓度0%~20%的乙醇溶液。

3.2 实验结果与分析

图7为实验中线阵CCD和面阵CCD采集的纯净水的一阶彩虹，图7(a)为线阵CCD连续采集的512帧彩虹分布；图7(b)为面阵CCD采集的一阶彩虹分布。从同次采集的512帧数据中取连续的2帧数据，对2帧数据进行EMD分析，且对EMD结果中对应模态间互相关系数进行计算，结果如表2所示，由表中数据分析，IMF₁为噪声分量，IMF₂、IMF₃、IMF₄为有用信号与噪声分量的混合信号，由于部分噪声频率接近Ripple频率，噪声分量和信号分量会在同一个IMF上混叠，混叠信号只影响FFT强度，不影响Ripple频率。对原信号减去高频模态，得到Airy分布，与原信号的对比如图8所示，从图中可以看出，滤波后的Airy分量与原彩虹信号吻合。由Airy分量得到θ₁=138°，θ_rg=137.77°。对原信号去噪结果如图9所示。对去噪信号进行FFT，结果如图10所示，得到F_Ripple=38.81。根据反演算法，折射率和直径的多次修正结果如表3所示。最后将修正结果代入(11)式得到最优结果为m=1.332 3，D=1 881 μm。液滴直径的参照值由显微相机确定为1 883 μm；由文献[21]可知，在该实验温度及632.8 nm波长下，水的折射率为1.332 0，与实验反演结果吻合。

图 7. First-order rainbow distribution

Fig. 7. First-order rainbow distribution

下载图片 查看所有图片

图 8. Comparison between detected signal and signal filtered by EMD

Fig. 8. Comparison between detected signal and signal filtered by EMD

下载图片 查看所有图片

图 9. Rainbow signal after wave filtering by EMD

Fig. 9. Rainbow signal after wave filtering by EMD

下载图片 查看所有图片

表4是不同浓度不同直径的乙醇液滴的测量结果与参照值的对比，折射率参照值由文献[22]给出，直径参照值为显微相机测量值。

图 10. Frequency spectrum of rainbow signal after denoised

Fig. 10. Frequency spectrum of rainbow signal after denoised

下载图片 查看所有图片

4 结论

本文基于经验模态分解法对彩虹信号分解和重构，结合一级Airy峰的角位置和Ripple结构与折射率、直径间的关系，设计了彩虹−经验模态分解算法对液滴的直径和折射率测量。利用相关运算，可快速分辨经验模态分解后的各个模态属性，能准确分离彩虹强度分布中Airy、Ripple结构及噪声。实验测量了水滴和不同浓度乙醇液滴的一阶彩虹分布，并通过彩虹−经验模态分解法反演获得了液滴的直径和折射率信息。研究结果表明，折射率在1.3~1.5，直径在100 μm~2 000 μm，彩虹−经验模态分解法的折射率的测量精度达到10⁻⁴，直径的测量误差不超过1%。本研究的方法具有设计简单、计算用时少、测量结果可靠、适用范围更广、抗噪声强的特征。