王保珺，申晉，李鑫強，王欽，劉偉，王雅靜，明虎

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

0 引言

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

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

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

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

對于振幅滿足高斯分布的散射場，光強ACFG（2）（τ）與電場ACFg（1）（τ）間滿足Siegert 關系，即

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

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

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

式（2）的離散形式為

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

式（4）的矩陣形式為

式中，g為歸一化電場ACF 的向量形式，其元素為g（1）（τj），A為電場ACF 數(shù)據(jù)對應的核矩陣，其元素為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)定性和準確性。

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

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

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

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

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

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

式中，λa、λb為產(chǎn)生交叉的兩個個體，β為交叉發(fā)生的概率。每個個體都有發(fā)生突變的概率，設第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 的離散點數(shù)設定為M=100。

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

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

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

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

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

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

圖1 為不同噪聲水平下采用L-curve 準則和MDP-GA 兩種正則參數(shù)選取方法對400 nm 窄分布顆粒的反演結(jié)果。圖2 為400 nm 寬分布顆粒的反演結(jié)果。表2 給出了對應圖1 和圖2 反演結(jié)果的性能指標和兩種算法的運行時間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 反演的性能指標Table 2 Performance index of PSD inversion

從圖1和表2可以看出，對于400 nm 窄分布顆粒，兩種參數(shù)選取方法得到的反演結(jié)果相同。對于400 nm 寬分布顆粒，由圖2 可知，在δ=10-5的噪聲水平下，L-curve 準則反演結(jié)果的峰值位置向粒徑增大方向偏移。隨著噪聲增加，在δ=10-4和δ=10-3噪聲水平下，該方法所得分布中小粒徑位置出現(xiàn)了虛假峰，且主峰分布變窄，峰值高于實際峰高。當噪聲增加到通常認為的極限水平δ=10-2時，兩種方法反演得到的PSD 都出現(xiàn)峰寬變窄且峰值位置向粒徑減小方向偏移的現(xiàn)象。對照表2 可以看出，在δ=10-5、δ=10-4和δ=10-3等較低或常規(guī)噪聲水平下，MDP-GA 方法均給出優(yōu)于L-curve 準則的峰值誤差和分布誤差指標。同時，在表2 中也能看到，相對于L-curve 準則方法，MDP-GA 方法的運算時間有明顯增加。考慮到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ù)字相關器組成。測量的散射角度為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）標準聚苯乙烯乳膠顆粒制備，所用分散劑為去離子蒸餾水，折射率和粘度系數(shù)分別為1.330 和0.891 mPa·s。兩種樣品顆粒的光強ACF、電場ACF、電場ACF 中的噪聲分量和PSD 反演結(jié)果如圖4、圖5。表3 為反演結(jié)果的性能指標。

圖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 標準聚苯乙烯乳膠顆粒反演的性能指標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）可以看出，目測良好的光強ACF 數(shù)據(jù)，經(jīng)Siegert 關系式中的開平方運算，會在求取的電場ACF 中明顯放大。圖4（c）和圖5（c）表明，放大后的噪聲隨延遲時間增加而增大，從中可以看出光強ACF 中的噪聲對反演結(jié)果的作用和根據(jù)噪聲調(diào)節(jié)正則化強度的重要性。

4.2 寬PSD 反演

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

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

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

表4 寬PSD 反演的性能指標Table 4 Performance index of wide PSD inversion

從表4 可以看出，對于納米級金剛石微粉和鈦酸鋇顆粒，MDP-GA 方法均給出了低于L-curve 準則的峰值重復性誤差，在圖6（d）給出的L-curve 準則方法反演納米級金剛石微粉所得PSD 中，小粒徑位置出現(xiàn)了虛假峰。結(jié)合模擬數(shù)據(jù)的反演結(jié)果可以看出，對于窄分布顆粒，兩種方法選取正則參數(shù)的反演結(jié)果沒有明顯差別，但對于寬分布顆粒，MDP-GA 方法所得反演結(jié)果的峰值誤差、分布誤差和重復性誤差均低于L-curve 準則方法。需要說明的是，圖6（a）和圖7（a）給出的光強ACF 中的縱軸截距較高，這是由于實測工業(yè)顆粒中存在混入大顆?；蝾w粒聚集等引起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é)論

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