基于經(jīng)驗(yàn)?zāi)B(tài)分解的彩虹法測(cè)量研究

羅道斌，秦毅盼，岳宗敏，騫來(lái)來(lái)，師 博，謝嬌嬌

（陜西科技大學(xué) 文理學(xué)院，陜西 西安 710021）

引言

折射率和直徑是能表征微小粒子特性的重要參數(shù)，高精度測(cè)量微小粒子的折射率和直徑有許多應(yīng)用，如噴霧液滴蒸發(fā)、燃燒等過(guò)程的測(cè)量[1-2]。彩虹測(cè)量技術(shù)由于其非接觸、高精度、快速、大信息量、可定量計(jì)算等優(yōu)點(diǎn)，在測(cè)量微小粒子的折射率和直徑中有廣泛的應(yīng)用[3-4]。隨著彩虹技術(shù)的發(fā)展，小粒子的形狀不再局限于球形，彩虹法同樣適用于扁圓、橢球粒子的測(cè)量[5-7]。此外，彩虹技術(shù)也可測(cè)量液柱或液滴折射率梯度和溫度梯度[8]，或測(cè)量噴霧場(chǎng)群粒子平均直徑和平均折射率[9-10]。在提高彩虹測(cè)量精度方面，有許多研究在測(cè)量方法及計(jì)算方法的優(yōu)化上。通過(guò)計(jì)算折射率或直徑的微小變化量與彩虹Ripple 結(jié)構(gòu)相位偏移之間的關(guān)系，可對(duì)折射率或直徑的相對(duì)變化進(jìn)行測(cè)量[11]；利用高階彩虹條紋對(duì)于橢球度的敏感度非常高的特征，通過(guò)高階彩虹可以更準(zhǔn)確反演扁圓液滴的橢球度[12]；通過(guò)研究高斯光束入射下的彩虹分布規(guī)律，采用合適的高斯光束腰大小和光束位置可提高信號(hào)強(qiáng)度，提高反演精度[13]；結(jié)合最小二乘優(yōu)化算法和彩虹技術(shù)反演算法，可對(duì)粒徑分布及折射率實(shí)現(xiàn)高精度測(cè)量[14]。

經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition,EMD)，是一種非線性、非平穩(wěn)信號(hào)分解方法[15]。這種方法無(wú)須預(yù)設(shè)任何基函數(shù)，本質(zhì)是通過(guò)數(shù)據(jù)的特征時(shí)間尺度或特征空間尺度，獲得信號(hào)的本征波動(dòng)模式，然后再分解數(shù)據(jù)，能將一個(gè)頻率不規(guī)則的信號(hào)波分解為多個(gè)單一頻率的波，和一個(gè)平均趨勢(shì)項(xiàng)波的形式。廣泛應(yīng)用于分析非線性、非平穩(wěn)信號(hào)序列[16-17]。本文結(jié)合經(jīng)驗(yàn)?zāi)B(tài)分解方法與液滴一階彩虹強(qiáng)度分布理論，研究了一種新的針對(duì)彩虹測(cè)量的反演算法，并用于液滴的折射率與直徑的同時(shí)測(cè)量。本方法簡(jiǎn)化了測(cè)量過(guò)程，具有去噪聲能力強(qiáng)、計(jì)算用時(shí)少，折射率與直徑的測(cè)量精度高的特征。

1 彩虹技術(shù)原理

一階彩虹強(qiáng)度分布主要包含Airy 結(jié)構(gòu)和Ripple 結(jié)構(gòu)，圖1為折射率m=1.332 4，直徑D=700 μm的一階彩虹強(qiáng)度分布。理論和實(shí)驗(yàn)研究表明，彩虹強(qiáng)度分布與入射光波長(zhǎng)、液滴折射率和液滴直徑有關(guān)[18]。

圖1 一階彩虹散射強(qiáng)度分布Fig.1 Distribution of first-order rainbow scattering intensity

由幾何光學(xué)理論，幾何光學(xué)彩虹角θrg與折射率m之間的關(guān)系為

根據(jù)彩虹的Airy 理論，彩虹角附近散射光虛擬波前的波動(dòng)方程為

式中z為無(wú)量綱參量，表示與彩虹角的歸一化角偏差，可表示為

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

幾何彩虹角位于z=0 處，由(2)、(3)式可求出幾何彩虹角處的光強(qiáng)為最大光強(qiáng)（一級(jí)Airy 峰值）的43.93%。

根據(jù)彩虹Airy 理論，第一個(gè)Airy 峰的角位置與折射率和直徑之間的關(guān)系表示為[19]

式中α為無(wú)量綱直徑參量α=πD/λ。根據(jù)圖2中的幾何關(guān)系，彩虹Ripple 結(jié)構(gòu)的頻率FRipple，可通過(guò)計(jì)算對(duì)應(yīng)于彩虹角時(shí)無(wú)限遠(yuǎn)處外反射光線與內(nèi)反射光線的光程差獲得[20]：

圖2 彩虹Ripple 結(jié)構(gòu)的形成原理Fig.2 The formation principle of rainbow Ripple structure

2 反演算法研究

2.1 彩虹信號(hào)的EMD 分析

EMD 是一種尺度分離算法，用I(θ)表示信號(hào)強(qiáng)度，首先確定I(θ)的所有極值點(diǎn)，再用三次樣條曲線構(gòu)造由極大值點(diǎn)組成的上包絡(luò)線u(θ)和由極小值點(diǎn)組成的下包絡(luò)線v(θ)，計(jì)算上下包絡(luò)線均值序列：

求原信號(hào)序列I(θ)與均值序列p(θ)的差值序列：

通過(guò)判斷q(θ)是否滿足：1）極值點(diǎn)數(shù)量與過(guò)零點(diǎn)數(shù)量相等或最多只差一個(gè)；2）在任意點(diǎn)q(θ)的上下包絡(luò)線均值為0。若滿足則q(θ)為分解出的第一個(gè)高頻本征模態(tài)函數(shù)(intrinsic mode function,IMF)，若不滿足，對(duì)q(θ)重復(fù)取上下包絡(luò)線均值序列，再計(jì)算差值序列，直至差值序列滿足停止條件，滿足停止條件的差值序列為分解出的第一個(gè)高頻模態(tài)IMF1。用I(θ)減去IMF1，對(duì)剩余分量繼續(xù)重復(fù)以上過(guò)程，直到分離出全部模態(tài)序列。信號(hào)最終可表示為

式中:IMFi(θ)表示第i個(gè)本征模態(tài)函數(shù)；R(θ)表示彩虹強(qiáng)度總體趨勢(shì)。

圖3為模擬m=1.332 4，D=700 μm 彩虹信號(hào)的EMD 分析結(jié)果，IMF1和IMF2代表彩虹的Ripple分量，IMF3代表彩虹的Airy 分量，最后一項(xiàng)R代表了彩虹信號(hào)的平均趨勢(shì)。從EMD 分析結(jié)果中可有效分離Airy 分量和Ripple 分量，并可得到第一個(gè)Airy 峰的角位置θ1。

圖3 彩虹信號(hào)經(jīng)驗(yàn)?zāi)B(tài)分解Fig.3 EMD of rainbow signal

實(shí)測(cè)彩虹信號(hào)含有噪聲，對(duì)同一液滴連續(xù)采集的2 幀彩虹信號(hào)做EMD 分析，求2 組EMD 結(jié)果對(duì)應(yīng)模態(tài)間的互相關(guān)系數(shù)，由于噪聲的隨機(jī)性，代表噪聲分量的對(duì)應(yīng)模態(tài)不相關(guān)，代表彩虹有效分量的對(duì)應(yīng)模態(tài)應(yīng)高度相關(guān)，可由此判斷彩虹信號(hào)的EMD 結(jié)果中的噪聲模態(tài)，對(duì)彩虹信號(hào)進(jìn)行去噪。互相關(guān)系數(shù)可表示為

式中：N為數(shù)據(jù)點(diǎn)個(gè)數(shù)；X和Y分別代表2幀采集數(shù)據(jù)的對(duì)應(yīng)模態(tài)序列。當(dāng)r＞0.8 時(shí)判定該模態(tài)高度相關(guān)，即該模態(tài)為有效信號(hào)；當(dāng)r＜0.3 時(shí)，判定為不相關(guān)，即該模態(tài)為噪聲分量。

為了驗(yàn)證上述去噪方法，理論模擬了m=1.332 4，D=700 μm 的液滴的一階彩虹信號(hào)，并對(duì)該信號(hào)加上2 組隨機(jī)噪聲，得到2 組加了不同隨機(jī)噪聲的彩虹信號(hào)，如圖4所示。

圖4 彩虹信號(hào)Fig.4 Rainbow signal

圖5為圖4中2 組加入隨機(jī)噪聲后的彩虹信號(hào)的EMD 結(jié)果。2 組EMD 結(jié)果對(duì)應(yīng)模態(tài)間的互相關(guān)系數(shù)計(jì)算結(jié)果如表1所示，由表1中數(shù)據(jù)分析，IMF1為噪聲分量，IMF2和IMF3為Ripple 結(jié)構(gòu)分量，IMF4為Airy 結(jié)構(gòu)分量，R反映光強(qiáng)趨勢(shì)。

圖5 彩虹信號(hào)的經(jīng)驗(yàn)?zāi)B(tài)分解結(jié)果Fig.5 EMD results of rainbow signal

表1 直徑和折射率的修正結(jié)果Table1 The inversion results of diameter and refractive index

結(jié)合表1分析結(jié)果，對(duì)圖5(a)加入隨機(jī)噪聲的彩虹信號(hào)減去IMF1，得到濾噪結(jié)果，濾噪結(jié)果與原信號(hào)相關(guān)系數(shù)為0.82，說(shuō)明了彩虹信號(hào)EMD 去噪技術(shù)的有效性。

2.2 反演算法的設(shè)計(jì)

折射率在1.3～1.5，直徑在100 μm～2 000 μm 范圍內(nèi)，根據(jù)彩虹分布在z=0 處強(qiáng)度與最大峰值強(qiáng)度間的關(guān)系，可通過(guò)在一級(jí)Airy 峰小角度側(cè)搜索強(qiáng)度最接近最大強(qiáng)度的43.93%的角位置，代入(3)式，得到折射率初值m0。利用FRipple對(duì)折射率和直徑的變化表現(xiàn)出的不同強(qiáng)度的變化規(guī)律將m0代入(6)式，得到直徑D1。再由θ1對(duì)折射率和直徑的變化表現(xiàn)出的不同強(qiáng)度的變化規(guī)律，將D1代入(5)式，得到折射率m1。通過(guò)反復(fù)將得到的折射率和直徑分別代入(6)式和(5)式不斷修正結(jié)果，再通過(guò)求解(11)式，得到多次修正結(jié)果的最優(yōu)值。

式中：I0θ為實(shí)測(cè)彩虹信號(hào)；Ikθ為每組修正結(jié)果模擬的理論彩虹信號(hào)。

3 實(shí)驗(yàn)結(jié)果與分析

3.1 測(cè)量裝置及測(cè)量方法

實(shí)驗(yàn)系統(tǒng)如圖6所示。由注射泵（LINZ-8A，雷恩，上海）定量控制壓出液體，懸掛在光纖末端形成球形液滴，由激光器（λ=632.8 nm，P=2 mW）發(fā)出的平行光束經(jīng)平面鏡反射后照射液滴，經(jīng)液滴散射形成彩虹條紋，用CCD（SG-14-01k80-00-R，1 024像素）記錄彩虹分布。顯微相機(jī)（LP-530HS，倍特嘉，蘇州）用來(lái)監(jiān)控液滴和輔助測(cè)量。本實(shí)驗(yàn)在光學(xué)暗室進(jìn)行，實(shí)驗(yàn)溫度22.4 ℃。

圖6 實(shí)驗(yàn)裝置圖Fig.6 The experimental diagram

通過(guò)移動(dòng)和旋轉(zhuǎn)平面鏡，角度精度為0.017°，使CCD 恰好能接收到一階彩虹信號(hào)，旋轉(zhuǎn)角度為φ，CCD 單個(gè)像素大小為14 μm，與液滴中心距離為30.4 cm，可確定CCD 測(cè)量角寬度Δθ=2.704°，由φ和Δθ可確定測(cè)量的散射角范圍。CCD 每次采集512 幀數(shù)據(jù)。實(shí)驗(yàn)液滴樣品為濃度0%～20%的乙醇溶液。

3.2 實(shí)驗(yàn)結(jié)果與分析

圖7為實(shí)驗(yàn)中線陣CCD 和面陣CCD 采集的純凈水的一階彩虹，圖7(a)為線陣CCD 連續(xù)采集的512 幀彩虹分布；圖7(b)為面陣CCD 采集的一階彩虹分布。從同次采集的512 幀數(shù)據(jù)中取連續(xù)的2 幀數(shù)據(jù)，對(duì)2 幀數(shù)據(jù)進(jìn)行EMD 分析，且對(duì)EMD 結(jié)果中對(duì)應(yīng)模態(tài)間互相關(guān)系數(shù)進(jìn)行計(jì)算，結(jié)果如表2所示，由表中數(shù)據(jù)分析，IMF1為噪聲分量，IMF2、IMF3、IMF4為有用信號(hào)與噪聲分量的混合信號(hào)，由于部分噪聲頻率接近Ripple 頻率，噪聲分量和信號(hào)分量會(huì)在同一個(gè)IMF上混疊，混疊信號(hào)只影響FFT 強(qiáng)度，不影響Ripple 頻率。對(duì)原信號(hào)減去高頻模態(tài)，得到Airy 分布，與原信號(hào)的對(duì)比如圖8所示，從圖中可以看出，濾波后的Airy 分量與原彩虹信號(hào)吻合。由Airy 分量得到θ1=138°，θrg=137.77°。對(duì)原信號(hào)去噪結(jié)果如圖9所示。對(duì)去噪信號(hào)進(jìn)行FFT，結(jié)果如圖10 所示，得到FRipple=38.81。根據(jù)反演算法，折射率和直徑的多次修正結(jié)果如表3所示。最后將修正結(jié)果代入(11)式得到最優(yōu)結(jié)果為m=1.332 3，D=1 881 μm。液滴直徑的參照值由顯微相機(jī)確定為1 883 μm；由文獻(xiàn)[21]可知，在該實(shí)驗(yàn)溫度及632.8 nm 波長(zhǎng)下，水的折射率為1.332 0，與實(shí)驗(yàn)反演結(jié)果吻合。

圖7 一階彩虹條紋圖Fig.7 First-order rainbow distribution

表2 2 幀數(shù)據(jù)EMD 結(jié)果中對(duì)應(yīng)模態(tài)間的互相關(guān)系數(shù)Table2 Correlation of two frames of measurement data

圖8 EMD 濾波后的Airy 分量與原信號(hào)的對(duì)比Fig.8 Comparison between detected signal and signal filtered by EMD

圖9 EMD 方法去噪后的彩虹信號(hào)圖Fig.9 Rainbow signal after wave filtering by EMD

表4是不同濃度不同直徑的乙醇液滴的測(cè)量結(jié)果與參照值的對(duì)比，折射率參照值由文獻(xiàn)[22]給出，直徑參照值為顯微相機(jī)測(cè)量值。

圖10 去噪后的彩虹信號(hào)頻譜Fig.10 Frequency spectrum of rainbow signal after denoised

表3 直徑和折射率的修正結(jié)果Table3 Inversion results of diameter and refractive index

表4 不同濃度乙醇液滴直徑和折射率的反演結(jié)果Table4 Inversion results of diameter and refractive index

4 結(jié)論

本文基于經(jīng)驗(yàn)?zāi)B(tài)分解法對(duì)彩虹信號(hào)分解和重構(gòu)，結(jié)合一級(jí)Airy 峰的角位置和Ripple 結(jié)構(gòu)與折射率、直徑間的關(guān)系，設(shè)計(jì)了彩虹-經(jīng)驗(yàn)?zāi)B(tài)分解算法對(duì)液滴的直徑和折射率測(cè)量。利用相關(guān)運(yùn)算，可快速分辨經(jīng)驗(yàn)?zāi)B(tài)分解后的各個(gè)模態(tài)屬性，能準(zhǔn)確分離彩虹強(qiáng)度分布中Airy、Ripple 結(jié)構(gòu)及噪聲。實(shí)驗(yàn)測(cè)量了水滴和不同濃度乙醇液滴的一階彩虹分布，并通過(guò)彩虹-經(jīng)驗(yàn)?zāi)B(tài)分解法反演獲得了液滴的直徑和折射率信息。研究結(jié)果表明，折射率在1.3～1.5，直徑在100 μm～2 000 μm，彩虹-經(jīng)驗(yàn)?zāi)B(tài)分解法的折射率的測(cè)量精度達(dá)到10-4，直徑的測(cè)量誤差不超過(guò)1%。本研究的方法具有設(shè)計(jì)簡(jiǎn)單、計(jì)算用時(shí)少、測(cè)量結(jié)果可靠、適用范圍更廣、抗噪聲強(qiáng)的特征。