非負約束全變差正則化方法在動態光散射顆粒粒度分布反演中的應用

王田田，申 晉，劉 偉，孫賢明

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

20世紀60年代以來，動態光散射（Dynamic Light Scattering，DLS）技術已成為亞微米、納米級顆粒粒徑測量的主要手段［1]．動態光散射技術的測量原理建立在顆粒的布朗運動基礎上，在分散介質中，顆粒的運動使得散射光強隨時間不斷地起伏漲落，這種漲落與顆粒的粒徑有關，粒徑越小，漲落越快．測量顆粒散射光信號的相關函數，再由相關函數反演出顆粒的粒徑分布，因此，動態光散射也稱為光子相關光譜（Photo Correlation Spectroscopy，PCS）技術．由相關函數反演顆粒粒度分布需求解第一類Fredholm積分方程，其高度病態性決定了原始數據微小的擾動都可能導致所求解與真實解的巨大偏離［2]．為解決這個問題，人們提出和應用了多種顆粒粒度反演算法，包括累積分析法、拉普拉斯變換法、指數采樣法、CONTIN算法、非負最小二乘法、Tikhonov正則化方法及截斷奇異值分解法［3－9]，但每種方法都有其適用范圍和局限性．一些改進的算法和顆粒粒度反演新方法不斷地被提出，包括改進的累計分析法［10]、非負最小二乘截斷奇異值分解法［11]、改進的正則化方法［12]及遺傳算法［13]、神經網絡算法［14]等智能優化算法．其中正則化方法因其較簡單的理論框架而得到廣泛的應用．

本文采用全變差（Total Variation，TV）正則化方法將離散的積分方程化為非線性最優化問題，并用固定點迭代法（Lagged Diffusivity Fixed Point Method）求解．實際測量中，粒度分布的誤差水平是不可知的，而反演過程中所用的系數矩陣A＇A是正定矩陣，固定點迭代法全局收斂［15]，所以算法的終止不取決于分布的誤差水平．非負約束條件的引入有效地消除了分布中無意義的負值．

1 動態光散射技術與全變差正則化反演

對隨機的散射光信號進行時間自相關運算并作歸一化處理，可以得到歸一化的散射光強自相關函數g（2）（τ），它與歸一化的電場自相關函數g（1）（τ）的關系為

式中，τ是延遲時間，A是實驗基線，β是相干因子．對于多分散顆粒系，歸一化的電場自相關函數

式中∫∞0G（Γ）dΓ＝1，是歸一化的衰減線寬分布函數Γ為衰減線寬DT為平移擴散系數；q為散射矢量；kB是波爾茲曼常數；T代表絕對溫度；η為溶液粘度系數；d為顆粒直徑；n為溶液的折射率；λ0為激光在真空中的波長；θ為散射角．

式（2）為第一類Fredholm積分方程，從帶有噪聲的g（1）求解顆粒粒徑屬于病態問題（ill－posed problem）．原始數據的微小擾動將導致解的嚴重偏離．求解時，將式（2）離散化為

采用全變差正則化方法求解式（3）．全變差正則化方法采用非線性正則化算子，實質是基于全變差懲罰的最小二乘方法［16]，即指求解如下最優化問題

式中α為正則化參數；TV（x）為正則化算子，起著恢復解的穩定性的作用．

由式（4）所得的結果中往往出現負解，因此對式（4）進行非負約束，即非負TV正則化

TV（x）為函數x定義在區間I＝［0，1] 上的全變差函數

要求得式（5）的極小值，首先要求Jα［x，g] 的梯度，由于TV（x）關于x是不可微的，TV（x）的梯度很難獲得，因此用一個光滑的泛函逼近TV（x），最優化問題化為

式中β0為正的常數．

采用固定點迭代法求解非線性最優化問題（7），求解步驟如下：

（1）設定迭代次數k＝50，假定初始粒度分布z＝0．

（2）計算L（通過離散全變差方程計算得到）．

（3）求解Jα［x，g] 的梯度，grad［Jα［x，g] ] ＝A＇（Ax－g）＋αLx．

（4）求解Jα［x，g] 的Hessian矩陣，Hessian［Jα［x，g] ] ＝A＇A＋αL．

（5）計算Sk＋1＝－Hessian［Jα［x，g] ]－1grad［Jα［x，g] ]k．

（6）更新分布xk＋1＝xk＋Sk＋1，直到滿足設定的迭代次數．

本文采用L曲線準則［17]選擇正則化參數α，即以log－log尺度描述‖xk‖與‖Axk－g‖的曲線，該曲線呈明顯的L狀，曲線拐角處的值即取為正則化參數．

2 數值模擬

利用非負約束TV正則化，對動態光散射反演問題進行數值模擬．計算過程中采用Johnson’s SB分布函數作為模擬的真實粒徑分布［18]，該函數的表達式為

式中，t＝（d－dmin）／（dmax－dmin）為歸一化尺寸，dmax和dmin是粒子的最大和最小的粒徑；u和σ是分布參數，通過改變這兩個參數，能模擬出不同的分布形式．

模擬實驗條件為入射光在真空中的波長λ0＝632．8nm，溶液的折射率n＝1．331，絕對溫度T＝298K，波爾茲曼常數kB＝1．3087×10－23J·K－1，溶液粘度系數η＝0．89×10－3Pa·s，散射角θ＝90°．

反演結果中，用分布誤差和峰值誤差來衡量反演得到的顆粒粒度分布（Particle Size Distribution，PSD）與真實粒度分布的匹配程度．分布誤差越小，峰值誤差越小，反演得到的分布越接近真實的粒度分布．本文給出的所有粒度分布均為五次反演結果的平均值．

2．1 無噪聲時單峰和雙峰顆粒粒度分布反演

單峰分布情況下，模擬初始粒徑分布參數為u＝0．8，σ＝1．5，αmax＝600nm，αmin＝100nm；雙峰分布顆粒為兩個等份額的單峰Johnson’s SB函數之和，參數為u1＝－3．2，σ1＝2．8，u2＝3．4，σ2＝2．1，αmax＝900nm，αmin＝100nm．單峰分布、雙峰分布顆粒的粒徑反演范圍初始值分別?。?，700] ，［0，1 000] ．單峰分布和雙峰分布的無噪聲相關函數反演結果分別如圖1、圖2所示．

圖1 無噪聲單峰分布反演

圖2 無噪聲雙峰分布反演

由圖1和圖2可以看出，非負TV正則化方法能很好地反演出單峰分布和雙峰分布，其中，單峰分布的分布誤差為0．013 3，峰值誤差為0；雙峰分布的分布誤差為0．061 5，第一個峰的峰值誤差為4．65%，第二個峰的峰值誤差為1．39%．在無噪聲情況下，非負TV正則化反演顆粒粒度分布是可行的．

2．2 加入噪聲時單峰和雙峰顆粒粒度分布反演

為了驗證算法的抗噪聲能力，分別在上述兩種顆粒的無噪聲模擬相關函數中加入噪聲水平為0．01的噪聲，其中噪聲相關函數為

g＿noise＝g＋noise＿level＊rand（N，1）；式中noise＿level為噪聲水平．

粒徑反演范圍初值同上．由非負TV正則化反演得到的單峰分布和雙峰分布如圖3和圖4所示．

圖3 噪聲水平為0．01時單峰分布反演結果

圖4 噪聲水平為0．01時雙峰分布反演結果

由圖3和圖4看出，加入噪聲水平為0．01的噪聲后，反演出的分布較大的偏離真實分布．單峰分布的分布誤差為0．0344 9，峰值誤差為6．56%；雙峰分布的分布誤差為0．089 9，第一個峰的峰值誤差為9．31%，第二個峰的峰值誤差為2．8%．雖然加入噪聲后，分布誤差和峰值誤差加大，但均在可接受范圍內．所以非負TV正則化具有一定的抗噪聲能力．

3 實驗結果

實驗溫度為25℃．光電倍增管在90接收散射光信號．實驗樣品為粒徑150nm（3150A）、60nm（3060A）和200nm（3200A）的聚苯乙烯乳膠顆粒．分散介質（水）的折射率為1．331．采用非負TV正則化方法反演測得的兩種顆粒系的相關函數，結果如圖5、圖6所示．

圖5 150nm聚苯乙烯顆粒反演結果

圖6 60nm和200nm混合顆粒反演結果

由非負TV正則化反演的150nm聚苯乙烯乳膠顆粒系的峰值為159．3nm，誤差為6．20%；反演的60nm和200nm混合顆粒的第一個峰值為61．3 nm，第二個峰值為221nm，第一個峰的峰值誤差為2%，第二個峰的峰值誤差為10．7%．非負TV正則化可以反演出實際的單分散和雙分散顆粒系．

4 結束語

本文采用全變差正則化方法對模擬和實驗數據進行了分析．結果表明，在無噪聲情況下，TV正則化反演的單峰分布和雙峰分布的分布誤差小于0．0615，峰值誤差小于4．65%．在相關函數中加入噪聲水平為0．01的噪聲后，分布誤差小于0．089 9，峰值誤差小于9．31%，所以，全變差正則化能較好地反演出單峰分布和雙峰分布．通過反演150nm標準顆粒和由60nm與200nm組成的混合顆粒驗證了這一結論．

［1] Pecora R．Dynamic light scattering：application of photon correlation spectroscopy［M] ．New York：Plenum Press，1985．

［2] 肖庭延，于慎根，王彥飛．反問題的數值解法［M] ．北京：科學出版社，2003：8－17．

［3] Koppel D E．Analysis of macromolecular polydispersity in intensity correlation spectroscopy：the method of cumulants［J] ．The Journal of Chemical Physics，1972，57：4 814－4 820．

［4] McWhirter J G and Pike E R．On the numerical inversion of the Laplace transform and similar Fredholm integral equations of the first kind［J] ．Phys．A Math．Theor，1978，11：1 729－1 745．

［5] Ostrowsky N，Sornette D，Pike E R．Exponential sampling method for light－scattering polydispersity analysis［J] ．Mod Opt，1981，28：1 059－1 070．

［6] Provencher S W．A general purpose constrained regularization program for inverting noisy linear algebraic and integral equations［J] ．Computer Physics Communication，1982，27：229－242．

［7] Iqbal M．On photon correlation measurements of colloidal size distributions using bayesian strategies［J] ．J Comput Appl Math，2000，126：77－89．

［8] 李功勝，于杰．一種新的求解第一類Fredholm積分方程正則化的數值分析［J] ．山東理工大學學報：自然科學版，2003，17（2）：5－8．

［9] Ubera J V，Aguilar J F，Gale D M．Reconstruction of particlesize distribution from light－scattering patterns using three inversion methods［J] ．Appl Opt，2007，40：124－132．

［10] Frisken B J．Revisiting the method of cumulants for analysis of dynamic light－scattering data［J] ．Appl．Opt，2001，40：4 087－4 091．

［11] Zhu X J，Shen J．Nonnegative least－square truncated singular value decomposition to particle size distribution inversion from dynamic light scattering data［J] ．Appl Opt，2010，49（34）：6 591－6 596．

［12] Arias M L，Frontini G L．Particle size distribution retrieval from elastic light scattering measurement by a modified regularization method［J] ．Part．Part．Syst．Charact，2006 23：374－380．

［13] Ye M，Wang S．Inversion of particle－size distribution from angular light－scattering data with genetic algorithms［J] ．Appl Opt，1999，38（12）：2 677－2 685．

［14] Gugliotta L M，Stegmayer G S．A neural network model for estimating the particle size distribution of dilute latex from multiangle dynamic light scattering measurements［J] ．Part Part Syst Charact，2009 26：41－52．

［15] Vogel C R．Computational methods for inverse problems［M] ．Siam，Frontiers in Applied Mathematics，2002，130－136．

［16] 王彥飛．反演問題的計算方法及其應用［M] ．北京：高等教育出版社，2007，38－49．

［17] Agarwal V．Total variation regularization and L－curve method for the selection of regularization parameter［D] ．2003，4－7．

［18] 李紹新．動態光散射測量粒徑分布的格雷碼編碼遺傳算法反演運算［J] ．計算物理，2008，25（3）：323－329．