非負(fù)約束全變差正則化方法在動(dòng)態(tài)光散射顆粒粒度分布反演中的應(yīng)用

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

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

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

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

1 動(dòng)態(tài)光散射技術(shù)與全變差正則化反演

對(duì)隨機(jī)的散射光信號(hào)進(jìn)行時(shí)間自相關(guān)運(yùn)算并作歸一化處理，可以得到歸一化的散射光強(qiáng)自相關(guān)函數(shù)g（2）（τ），它與歸一化的電場(chǎng)自相關(guān)函數(shù)g（1）（τ）的關(guān)系為

式中，τ是延遲時(shí)間，A是實(shí)驗(yàn)基線，β是相干因子．對(duì)于多分散顆粒系，歸一化的電場(chǎng)自相關(guān)函數(shù)

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

式（2）為第一類Fredholm積分方程，從帶有噪聲的g（1）求解顆粒粒徑屬于病態(tài)問(wèn)題（ill－posed problem）．原始數(shù)據(jù)的微小擾動(dòng)將導(dǎo)致解的嚴(yán)重偏離．求解時(shí)，將式（2）離散化為

采用全變差正則化方法求解式（3）．全變差正則化方法采用非線性正則化算子，實(shí)質(zhì)是基于全變差懲罰的最小二乘方法［16]，即指求解如下最優(yōu)化問(wèn)題

式中α為正則化參數(shù)；TV（x）為正則化算子，起著恢復(fù)解的穩(wěn)定性的作用．

由式（4）所得的結(jié)果中往往出現(xiàn)負(fù)解，因此對(duì)式（4）進(jìn)行非負(fù)約束，即非負(fù)TV正則化

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

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

式中β0為正的常數(shù)．

采用固定點(diǎn)迭代法求解非線性最優(yōu)化問(wèn)題（7），求解步驟如下：

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

（2）計(jì)算L（通過(guò)離散全變差方程計(jì)算得到）．

（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）計(jì)算Sk＋1＝－Hessian［Jα［x，g] ]－1grad［Jα［x，g] ]k．

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

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

2 數(shù)值模擬

利用非負(fù)約束TV正則化，對(duì)動(dòng)態(tài)光散射反演問(wèn)題進(jìn)行數(shù)值模擬．計(jì)算過(guò)程中采用Johnson’s SB分布函數(shù)作為模擬的真實(shí)粒徑分布［18]，該函數(shù)的表達(dá)式為

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

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

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

2．1 無(wú)噪聲時(shí)單峰和雙峰顆粒粒度分布反演

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

圖1 無(wú)噪聲單峰分布反演

圖2 無(wú)噪聲雙峰分布反演

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

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

為了驗(yàn)證算法的抗噪聲能力，分別在上述兩種顆粒的無(wú)噪聲模擬相關(guān)函數(shù)中加入噪聲水平為0．01的噪聲，其中噪聲相關(guān)函數(shù)為

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

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

圖3 噪聲水平為0．01時(shí)單峰分布反演結(jié)果

圖4 噪聲水平為0．01時(shí)雙峰分布反演結(jié)果

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

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

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

圖5 150nm聚苯乙烯顆粒反演結(jié)果

圖6 60nm和200nm混合顆粒反演結(jié)果

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

4 結(jié)束語(yǔ)

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

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

［2] 肖庭延，于慎根，王彥飛．反問(wèn)題的數(shù)值解法［M] ．北京：科學(xué)出版社，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積分方程正則化的數(shù)值分析［J] ．山東理工大學(xué)學(xué)報(bào)：自然科學(xué)版，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，F(xiàn)rontini 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，F(xiàn)rontiers in Applied Mathematics，2002，130－136．

［16] 王彥飛．反演問(wèn)題的計(jì)算方法及其應(yīng)用［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] 李紹新．動(dòng)態(tài)光散射測(cè)量粒徑分布的格雷碼編碼遺傳算法反演運(yùn)算［J] ．計(jì)算物理，2008，25（3）：323－329．