基于改進(jìn)Morozov 偏差原理的動態(tài)光散射粒度反演

王?，B，申晉，李鑫強(qiáng)，王欽，劉偉，王雅靜，明虎

（山東理工大學(xué) 電氣與電子工程學(xué)院，淄博 255049）

0 引言

動態(tài)光散射（Dynamic Light Scattering，DLS）技術(shù)是測量顆粒粒度分布（Particle Size Distribution，PSD）的有效方法［1-2］，該技術(shù)以其快速、準(zhǔn)確和非接觸等優(yōu)點(diǎn)［3］在材料［4-5］、化工［6］、食品［7-8］、醫(yī)藥［9-11］等領(lǐng)域得到廣泛應(yīng)用。DLS 技術(shù)通過對散射光強(qiáng)信號進(jìn)行自相關(guān)運(yùn)算獲得光強(qiáng)自相關(guān)函數(shù)（Autocorrelation Function，ACF），通過反演光強(qiáng)ACF 獲得待測顆粒的PSD。反演ACF 需要求解第一類Fredholm 積分方程，該方程是一典型的病態(tài)方程，測量數(shù)據(jù)中的微小誤差或噪聲均會導(dǎo)致所求PSD 的顯著變化。為提高PSD 反演的準(zhǔn)確性，研究人員已提出了多種反演方法，包括奇異值分解法［12］、CONTIN 法［13，14］、非負(fù)約束最小二乘法［15，16］和Tikhonov 正則化方法［17］等。這些方法各有其特點(diǎn)，其中的Tikhonov 正則化方法由于不受粒度分布的限制，具有良好的適應(yīng)性，在DLS 測量中得到了廣泛的應(yīng)用。該方法通過引入穩(wěn)定泛函改善病態(tài)性，由正則化參數(shù)控制穩(wěn)定泛函的修正程度，因此對反演結(jié)果有重要影響［18］。如果選擇的正則化參數(shù)過小，反演的PSD中會出現(xiàn)振蕩和虛假峰，而正則化參數(shù)過大，則會出現(xiàn)過于平滑的反演結(jié)果［19］。選擇正則化參數(shù)的主要方法包括L-curve 準(zhǔn)則［20-21］、Morozov 偏差原理（Morozov's Discrepancy Principle，MDP）［22-23］以及廣義交叉驗證準(zhǔn)則（Generalized Cross Validation，GCV）［24，25］等。

在DLS 測量技術(shù)中，進(jìn)行正則化反演通常采用L-curve 準(zhǔn)則選取正則參數(shù)。2012 年，劉曉艷［26］采用Lcurve 準(zhǔn)則進(jìn)行正則化反演顆粒粒度分布，研究了DLS 技術(shù)對散射角的依賴關(guān)系。2016 年，林承軍等［27］將多參數(shù)正則化用于前向散射測量中的粒度反演，降低了正則化算法帶來的振蕩和負(fù)值。2022 年，LIU Z 等［28］和韓錦壯等［29］分別采用L-curve 準(zhǔn)則選取正則參數(shù)進(jìn)行了流動顆粒DLS 粒度反演，文獻(xiàn)［28］對L-curve 準(zhǔn)則和GCV 兩種正則參數(shù)選取方法進(jìn)行了比較，認(rèn)為在噪聲水平較低時，GCV 方法選取正則參數(shù)的反演結(jié)果優(yōu)于L-curve 準(zhǔn)則的結(jié)果，而當(dāng)噪聲水平較高時，L-curve 準(zhǔn)則的反演效果則優(yōu)于GCV 方法。文獻(xiàn)［29］通過條件預(yù)優(yōu)處理降低了正則化對流速增加的敏感性，從而改善了流動顆粒DLS 反演的性能指標(biāo)。

L-curve 準(zhǔn)則是通過解的向量模間接引入穩(wěn)定性分析，但其在寬粒度分布條件下通常不能得到準(zhǔn)確的PSD 反演結(jié)果。MDP 方法可依據(jù)電場ACF 的噪聲水平選取正則參數(shù)以適應(yīng)不同的粒度分布，但原始數(shù)據(jù)噪聲水平通常是未知的，這一條件限制使其難以在DLS 的實際測量中得以應(yīng)用。與窄粒度分布顆粒的DLS反演相比，寬分布顆粒反演難以獲取與之相適應(yīng)的正則參數(shù)。為解決寬粒度分布反演時正則參數(shù)不易準(zhǔn)確選取的問題，本文提出基于改進(jìn)Morozov 偏差原理的MDP-GA 正則參數(shù)求取方法，即采用小波包分解（Wavelet Packet Decomposition，WPD）結(jié)合MDP 的方式建立迭代判據(jù)，利用遺傳算法（Genetic Algorithm，GA）迭代求取正則參數(shù)，較好地解決了寬分布顆粒DLS 反演中正則參數(shù)不易準(zhǔn)確選取的難題。

1 動態(tài)光散射與正則化反演

對于振幅滿足高斯分布的散射場，光強(qiáng)ACFG（2）（τ）與電場ACFg（1）（τ）間滿足Siegert 關(guān)系，即

式中，τ為延遲時間，B為測量基線，β為相干因子，且β≤1。對于多分散顆粒體系，g（1）（τ）表示為

式中，G（Γ）為衰減線寬分布函數(shù)，Γ為衰減線寬，可由Stokes-Einstein 公式求得。

式中，d為顆粒粒徑，KB、T、η、λ0、n和θ分別為Boltzmann 常數(shù)、絕對溫度、介質(zhì)粘度系數(shù)、入射光波長、懸浮介質(zhì)的折射率和散射角度。

式（2）的離散形式為

式中，τj為離散的延遲時間，N和i分別為離散粒度的總點(diǎn)數(shù)和第i個離散點(diǎn)，M和j分別為光子相關(guān)器的總通道數(shù)和第j通道，f（di）為的粒度分布，滿足通過衰減線寬—平移擴(kuò)散系數(shù)—顆粒粒度關(guān)系，可求得顆粒粒度分布f（d）。

式（4）的矩陣形式為

式中，g為歸一化電場ACF 的向量形式，其元素為g（1）（τj），A為電場ACF 數(shù)據(jù)對應(yīng)的核矩陣，其元素為exp（-Γiτj），f是由離散的PSD 組成的向量，其元素為f（di）。

式（5）是病態(tài)方程，通常采用Tikhonov 正則化方法將其求解轉(zhuǎn)化為未知函數(shù)的約束優(yōu)化問題，即

式中，M、α、‖·‖和‖Lf‖分別為穩(wěn)定泛函、正則參數(shù)、歐幾里得范數(shù)和懲罰因子。正則矩陣L選用二階差分矩陣，α控制解的穩(wěn)定性和準(zhǔn)確性。

2 MDP-GA 正則參數(shù)選取方法

本文提出的基于改進(jìn)Morozov 偏差原理的MDP-GA 方法，是一種以GA 全局尋優(yōu)求取正則參數(shù)的方法。該方法在正則參數(shù)經(jīng)驗范圍生成初始種群，利用WPD 求出電場ACF 的噪聲分量，并將該分量的噪聲水平代入MDP 判據(jù)，建立適應(yīng)值函數(shù)。通過保留每代種群中適應(yīng)值最優(yōu)的個體，來保證迭代的進(jìn)化方向，進(jìn)而對種群做交叉和突變運(yùn)算，生成下一代種群。

作為待求正則參數(shù)的隨機(jī)初始值，種群中的初始個體目標(biāo)變量可表示為

式中，λmax和λmin分別為目標(biāo)變量允許的最大和最小值，在DLS 正則參數(shù)選取中對應(yīng)值分別取20 和0，Rand為［0，1］區(qū)間內(nèi)的隨機(jī)數(shù)。

為評價種群中個體的優(yōu)劣，根據(jù)MDP 判據(jù)‖Afα-g‖2-δg2=0 建立適應(yīng)值函數(shù)f，表示為

式中，δg為g的噪聲水平，由WPD 給出。在本代種群做交叉和突變運(yùn)算前，需要保留本代的最優(yōu)個體λbest，即本代種群中適應(yīng)值最大的個體。每兩個個體都有發(fā)生交叉的概率，發(fā)生交叉的兩個個體由式（9）求得。

式中，λa、λb為產(chǎn)生交叉的兩個個體，β為交叉發(fā)生的概率。每個個體都有發(fā)生突變的概率，設(shè)第t代種群為（λ1…λi…λm），若僅λi發(fā)生突變，則t+1 代種群變?yōu)椋é?…λi'…λm），λi'表示為

3 模擬數(shù)據(jù)反演

采用半對數(shù)正態(tài)分布模型模擬單峰PSD，模擬實驗條件為：分散介質(zhì)折射率n=1.331 6，入射光在真空中的波長λ=632.8 nm，絕對溫度T=298.15 K，介質(zhì)粘度系數(shù)η=0.89 mPa·s，Boltzmann 常數(shù)KB=1.380 7×10-23J/K，相干因子β=0.7，基線B=1，PSD 的離散點(diǎn)數(shù)設(shè)定為M=100。

式中，Di、Dg和σ分別為離散的顆粒粒徑、標(biāo)稱粒徑和標(biāo)準(zhǔn)偏差，相應(yīng)的G（2）（τ）數(shù)據(jù)通過式（1）和式（5）獲得，表1 為兩種不同分布寬度的400 nm 顆粒的PSD 模擬參數(shù)。

表1 PSD 的分布參數(shù)和屬性Table 1 Parameters and properties of the simulated PSD

為接近實測情況，在模擬的光強(qiáng)ACF 中加入了高斯噪聲，含噪光強(qiáng)ACF 可表示為

式中，n（τ）和δ分別為高斯隨機(jī)噪聲和噪聲水平。

粒度反演時，首先由式（1）求得g（1）（τ），通過反演g（1）（τ）求取顆粒粒度。為評估反演結(jié)果，引入峰值位置相對誤差（EP）和分布誤差（VE）兩個性能指標(biāo)。

式中，P表示峰值位置處的粒徑，下標(biāo)true 和meas 分別表示真實值和反演值。

圖1 為不同噪聲水平下采用L-curve 準(zhǔn)則和MDP-GA 兩種正則參數(shù)選取方法對400 nm 窄分布顆粒的反演結(jié)果。圖2 為400 nm 寬分布顆粒的反演結(jié)果。表2 給出了對應(yīng)圖1 和圖2 反演結(jié)果的性能指標(biāo)和兩種算法的運(yùn)行時間T。

圖1 窄分布400 nm 的反演結(jié)果Fig.1 Inversion results of 400 nm narrow PSD

圖2 寬分布400 nm 的反演結(jié)果Fig.2 Inversion results of 400 nm broad PSD

表2 PSD 反演的性能指標(biāo)Table 2 Performance index of PSD inversion

從圖1和表2可以看出，對于400 nm 窄分布顆粒，兩種參數(shù)選取方法得到的反演結(jié)果相同。對于400 nm 寬分布顆粒，由圖2 可知，在δ=10-5的噪聲水平下，L-curve 準(zhǔn)則反演結(jié)果的峰值位置向粒徑增大方向偏移。隨著噪聲增加，在δ=10-4和δ=10-3噪聲水平下，該方法所得分布中小粒徑位置出現(xiàn)了虛假峰，且主峰分布變窄，峰值高于實際峰高。當(dāng)噪聲增加到通常認(rèn)為的極限水平δ=10-2時，兩種方法反演得到的PSD 都出現(xiàn)峰寬變窄且峰值位置向粒徑減小方向偏移的現(xiàn)象。對照表2 可以看出，在δ=10-5、δ=10-4和δ=10-3等較低或常規(guī)噪聲水平下，MDP-GA 方法均給出優(yōu)于L-curve 準(zhǔn)則的峰值誤差和分布誤差指標(biāo)。同時，在表2 中也能看到，相對于L-curve 準(zhǔn)則方法，MDP-GA 方法的運(yùn)算時間有明顯增加。考慮到DLS 測量時間通常為數(shù)分鐘，因此，PSD 反演時間增加不足1.5 s，對測量不會產(chǎn)生影響。

4 實測數(shù)據(jù)反演

實驗數(shù)據(jù)取自自行研制的DLS 測量實驗平臺（如圖3），測量裝置由波長為532 nm 的固體激光器（MGL-III-532 nm-15 mW）、光電探測器（CH326）和512 通道的數(shù)字相關(guān)器組成。測量的散射角度為90°，采用單模光纖探針接收散射光，樣品池溫度控制在298.15 K，實驗樣品顆粒按分布的寬窄分為兩類。其中光纖探針由數(shù)值孔徑0.12、纖芯直徑3.5 μm 的單模保偏光纖和一個0.25 節(jié)距自聚焦透鏡耦合組成，透鏡發(fā)散角為2 mrad，光束直徑為0.4 mm。

圖3 實驗裝置示意Fig.3 Schematic of the experimental apparatus

4.1 窄分布粒度反演

兩個窄分布樣品分別為由（31±3） nm（Duke Scientific，3 030 A）和（203±5） nm（Duke Scientific，3 200 A）標(biāo)準(zhǔn)聚苯乙烯乳膠顆粒制備，所用分散劑為去離子蒸餾水，折射率和粘度系數(shù)分別為1.330 和0.891 mPa·s。兩種樣品顆粒的光強(qiáng)ACF、電場ACF、電場ACF 中的噪聲分量和PSD 反演結(jié)果如圖4、圖5。表3 為反演結(jié)果的性能指標(biāo)。

圖4 31 nm 聚苯乙烯乳膠顆粒的G（2）（τ）、g（1）（τ）、噪聲分量以及反演結(jié)果Fig.4 G（2）（τ）， g（1）（τ）， noise component and inversion results of 31 nm standard polystyrene latex particles

圖5 203 nm 聚苯乙烯乳膠顆粒的G（2）（τ）、g（1）（τ）、噪聲分量以及反演結(jié)果Fig.5 G（2）（τ）， g（1）（τ）， noise component and inversion results of 203 nm standard polystyrene latex particles

表3 標(biāo)準(zhǔn)聚苯乙烯乳膠顆粒反演的性能指標(biāo)Table 3 Performance index for inversion of standard polystyrene latex particles

從圖4（d）和圖5（d）可以看出，對于31 nm 和203 nm 兩種窄分布顆粒，采用兩種正則參數(shù)選取方法所得的反演結(jié)果幾乎無差異，表3 給出的兩種方法所得的平均粒徑相同。對比圖4（a）、5（a）與圖4（b）、5（b）可以看出，目測良好的光強(qiáng)ACF 數(shù)據(jù)，經(jīng)Siegert 關(guān)系式中的開平方運(yùn)算，會在求取的電場ACF 中明顯放大。圖4（c）和圖5（c）表明，放大后的噪聲隨延遲時間增加而增大，從中可以看出光強(qiáng)ACF 中的噪聲對反演結(jié)果的作用和根據(jù)噪聲調(diào)節(jié)正則化強(qiáng)度的重要性。

4.2 寬PSD 反演

兩個寬分布樣品分別為由納米級金剛石微粉（Nanodiamond）和鈦酸鋇（BaTiO3）制備。為評估反演結(jié)果，采用重復(fù)性誤差（δr）作為被測顆粒的性能評價指標(biāo)。

式中，xi表示第i次測量得到的峰值粒徑，xˉ代表平均峰值粒徑。由兩種樣品得到的光強(qiáng)ACF、電場ACF、電場ACF 中的噪聲分量和PSD 反演結(jié)果如圖6、7。表4 為進(jìn)行6 次測量的峰值重復(fù)性誤差。

圖6 納米金剛石微粉的G（2）（τ）、g（1）（τ）、噪聲分量以及反演結(jié)果Fig.6 G（2）（τ）， g（1）（τ）， noise component and inversion results of nanodiamond particle system

表4 寬PSD 反演的性能指標(biāo)Table 4 Performance index of wide PSD inversion

從表4 可以看出，對于納米級金剛石微粉和鈦酸鋇顆粒，MDP-GA 方法均給出了低于L-curve 準(zhǔn)則的峰值重復(fù)性誤差，在圖6（d）給出的L-curve 準(zhǔn)則方法反演納米級金剛石微粉所得PSD 中，小粒徑位置出現(xiàn)了虛假峰。結(jié)合模擬數(shù)據(jù)的反演結(jié)果可以看出，對于窄分布顆粒，兩種方法選取正則參數(shù)的反演結(jié)果沒有明顯差別，但對于寬分布顆粒，MDP-GA 方法所得反演結(jié)果的峰值誤差、分布誤差和重復(fù)性誤差均低于L-curve 準(zhǔn)則方法。需要說明的是，圖6（a）和圖7（a）給出的光強(qiáng)ACF 中的縱軸截距較高，這是由于實測工業(yè)顆粒中存在混入大顆粒或顆粒聚集等引起ACF 基線抬升所致，在圖6（b）和圖7（b）給出的電場ACF 的長延遲時段可清晰看到。

圖7 鈦酸鋇顆粒的G（2）（τ）、g（1）（τ）、噪聲分量以及反演結(jié)果Fig.7 G（2）（τ）， g（1）（τ）， noise component and inversion results of BaTiO3 particle system

5 結(jié)論

本文提出了基于改進(jìn)Morozov 偏差原理的MDP-GA 正則參數(shù)選取方法，該方法可以顯著改善寬粒度分布顆粒體系的DLS 正則化反演效果。首先，通過小波包分解求出電場ACF 的噪聲分量，得到其噪聲水平。然后，通過MDP 建立適應(yīng)值函數(shù)，在正則參數(shù)經(jīng)驗范圍內(nèi)生成初始種群。最后，將適應(yīng)值函數(shù)與初始種群帶入遺傳算法，通過全局尋找最優(yōu)適應(yīng)值對應(yīng)的參數(shù)值作為正則參數(shù)。模擬與實測數(shù)據(jù)的反演結(jié)果表明，對于窄分布顆粒體系，所提MDP-GA 方法與L-curve 準(zhǔn)則方法的反演結(jié)果無顯著差異，對于寬分布顆粒體系，MDP-GA 方法反演結(jié)果的性能指標(biāo)優(yōu)于L-curve 準(zhǔn)則，避免了正則參數(shù)選取不當(dāng)導(dǎo)致的PSD 中出現(xiàn)虛假峰和峰值偏移情況，得到了更為準(zhǔn)確的寬分布顆粒體系的反演結(jié)果。鑒于復(fù)雜粒度分布的準(zhǔn)確測量需要采用多角度DLS 方法，利用MDP-GA 方法在多角度DLS 反演中進(jìn)行正則參數(shù)選取，將是本研究后續(xù)的主要工作。