基于模擬退火的隨機(jī)取向含核橢球顆粒粒徑反演

陳 都，鄧小玖，李耀東，王 鋒，張永明＊

（1.合肥工業(yè)大學(xué)電子科學(xué)與應(yīng)用物理學(xué)院，合肥，230009；2.中國科學(xué)技術(shù)大學(xué)火災(zāi)科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，合肥，230026）

0 引言

大氣中懸浮顆粒物大多數(shù)為非球形，且并非僅由單種勻質(zhì)成分組成，常見的含核顆粒即是其中一種，例如細(xì)菌、霧滴、城市水凝物氣溶膠顆粒［1－3］，以及細(xì)水霧滅火抑制煙氣和濕式洗滌塔除塵過程中生成的水煙凝并顆粒等。隨著工業(yè)文明的進(jìn)程，總體上城市大氣環(huán)境中可吸入顆粒物逐漸增多，對(duì)人類生活環(huán)境和大氣氣候產(chǎn)生了影響，在一定程度上危害了人體健康及人類生活質(zhì)量［4－7］。所以，控制并減少可吸入顆粒物的產(chǎn)生和排放是當(dāng)下迫在眉睫的重大任務(wù)。

研究表明，細(xì)水霧對(duì)火災(zāi)煙氣具有一定的抑制作用，在細(xì)水霧滅火抑制煙氣過程中，粒徑較小的霧滴與懸浮在空氣中的火災(zāi)煙顆粒發(fā)生凝聚、合并，形成以碳煙顆粒作為凝結(jié)核、水為外殼的含核顆粒［8，9］；工業(yè)脫除燃燒顆粒物領(lǐng)域中，碳煙顆粒通過濕式洗滌塔，與塔腔內(nèi)的過飽和水蒸汽凝聚、合并，形成以碳煙顆粒作為凝結(jié)核、水為外殼的含核顆粒。以上兩種過程中形成的含核顆粒在重力和空氣阻力的雙重作用下外殼近似橢球形狀。濕式洗滌塔煙顆粒脫除實(shí)驗(yàn)研究表明：碳煙顆粒粒徑直接影響濕式洗滌塔的脫除效果［10］。對(duì)這類含核橢球顆粒粒徑和形狀準(zhǔn)確的檢測(cè)，可對(duì)細(xì)水霧滅火抑制煙氣和濕式洗滌塔顆粒脫除過程進(jìn)行更好的控制，以期實(shí)現(xiàn)更有效地滅火和對(duì)碳煙顆粒更好的脫除效果。

目前市場(chǎng)上大多數(shù)激光粒度儀都是一種基于Lorenz－Mie散射理論由顆粒的散射光強(qiáng)來反演顆粒粒徑的探測(cè)儀器，而Lorenz－Mie散射理論是基于球形模型的理論［11］。本文將通過仿真模擬實(shí)驗(yàn)研究分析基于Lorenz－Mie散射理論的激光粒度儀對(duì)含核橢球顆粒的粒徑探測(cè)的適用性。

粒徑是顆粒的重要物理參數(shù)。在顆粒光散射過程中，散射光的角分布會(huì)隨著顆粒粒徑變化而變化，因此對(duì)散射光信號(hào)數(shù)據(jù)進(jìn)行反演可獲取顆粒粒徑的信息［12］。本文以含核橢球顆粒為模型，模擬了含球形碳煙核的橢球狀水滴顆粒光散射特性。利用DDA對(duì)隨機(jī)取向含核橢球顆粒模型進(jìn)行數(shù)值模擬計(jì)算得出Mueller矩陣，并基于Lorenz－Mie散射理論，采用全局搜索能力很強(qiáng)的模擬退火算法，實(shí)現(xiàn)了單分散系含核橢球顆粒的粒徑反演。

1 理論與方法

1.1 Mueller矩陣

對(duì)于任意顆粒的光散射，入射光與散射光Stokes向量的變換關(guān)系表示為［13］：

式中I表示光強(qiáng)，Q，U，V表示光波的偏振狀態(tài)，矩陣 ｛Fij（θ）（i，j＝1，2，3，4）｝為光散射 Mueller矩陣，包含6個(gè)獨(dú)立變量，描述光散射從初態(tài)到待測(cè)態(tài)的變化關(guān)系，Mueller矩陣元隨散射角的分布可以完整地描述散射光的強(qiáng)度與偏振特征［14］。

處理數(shù)據(jù)時(shí)，通常對(duì)Mueller矩陣進(jìn)行歸一化：矩陣元F11（θ）對(duì)應(yīng)于00散射角的值歸一化為F11（θ）／F11（0），其他矩陣元 Fij（θ）對(duì)應(yīng)于 F11（θ）歸一化為 Fij（θ）／F11（θ）。歸一化的 Mueller矩陣與光強(qiáng)絕對(duì)值大小無關(guān)，能更好地反映光散射過程的物理特征［12，13］。

1.2 DDA方法

DDA方法可以計(jì)算任意形狀顆粒模型的光散射和吸收特性［15，16］。DDA算法的基本思想是利用有限個(gè)離散的偶極子近似任意形狀的顆粒，偶極子對(duì)局域電場(chǎng)的響應(yīng)可以獲得偶極矩，偶極子間通過各自產(chǎn)生的電場(chǎng)相互作用。

設(shè)將散射體離散為N個(gè)偶極子，各偶極子的坐標(biāo)為rj（j＝1，2，···，N），極化率為αj，各偶極子的極化強(qiáng)度為Pj＝αjEj，其中Ej是rj處的電場(chǎng)，Ej是入射場(chǎng)Einc，j和其它N－1個(gè)偶極子所激發(fā)的電場(chǎng)的疊加：

應(yīng)用DDA算法時(shí)必須滿足以下兩個(gè)條件［17，18］：

2）d必須足夠?。∟必須足夠大），方能較好地描述目標(biāo)顆粒的形貌特征。

1.3 模擬退火算法

模擬退火算法（Simulate Anneal Arithmetic，SAA）的思想最早由 Metropolis等提出，1983年Kirkpatrick S.，Gelatt C.D.和Vecchi M.P.將其成功用于組合優(yōu)化。SAA是基于Monte Carlo迭代的一種全局概率型搜索算法，用來在一個(gè)大搜尋空間內(nèi)尋找命題的最優(yōu)解，其基本思想源于物理學(xué)中金屬的退火過程。SAA算法流程如圖1所示［19］：

SAA用一組冷卻進(jìn)度表參數(shù)來控制算法的進(jìn)程，控制參數(shù)“降溫”并趨于零時(shí)，最終求得組合優(yōu)化問題的相對(duì)全局最優(yōu)解。冷卻進(jìn)度表參數(shù)包括控制參數(shù)T的初始值T0、衰減函數(shù)、終止準(zhǔn)則及Markov鏈長度Lk，優(yōu)化問題的一個(gè)解i及其目標(biāo)函數(shù)f（i）分別與固體的一個(gè)微觀狀態(tài)i及其能量E（i）相對(duì)應(yīng)。從初始值X0和初始溫度T0開始，對(duì)當(dāng)前解重復(fù)進(jìn)行“在解空間中隨機(jī)產(chǎn)生新解 → Metropolis準(zhǔn)則判斷 → 接受或舍棄新解”的迭代，并按照一定的規(guī)則逐步降低T值，算法迭代終止時(shí)的當(dāng)前解既是近似最優(yōu)解。

圖1 模擬退火算法流程圖Fig.1 The flow chart of simulated annealing algorithm

2 計(jì)算結(jié)果與分析

本文采用勻質(zhì)球體作為內(nèi)核、水為外殼的軸對(duì)稱的含核橢球顆粒模型，橢球縱橫比（ε＝b／a，其中旋轉(zhuǎn)軸的長度為a，短軸為b）設(shè)定為0.8［1］，入射光的真空波長為632.8nm，水的真空相對(duì)折射率為1.33，顆粒尺寸參數(shù)x0（x＝2πaeff／λ）設(shè)定為20，即其有效半徑（有效半徑定義為顆粒物的等體積 球的半徑，即aeff＝ （3V／4π）1／3）aeff0為2.0143um。本文利用DDA方法模擬計(jì)算了兩種情況下的含核橢球顆粒模型的光散射特性，并討論分析了兩種情況下的粒徑反演結(jié)果和誤差。

第一種情況是設(shè)定含核橢球顆粒模型的內(nèi)核真空相對(duì)折射率n為1.57＋0.42i，討論內(nèi)核半徑與外殼有效半徑比q＝rcore／rshell＝rcore／aeff（變量q自0至1.0以0.1間隔遞增）的變化對(duì)粒徑反演結(jié)果和誤差的影響。

第二種情況是設(shè)定含核橢球顆粒模型的內(nèi)核半徑與外殼有效半徑比q為0.2，討論內(nèi)核真空相對(duì)折射率n（變量n自1.55＋0.42i至2.00＋0.42i，實(shí)部以0.05間隔遞增）的變化對(duì)粒徑反演結(jié)果和誤差的影響。

應(yīng)用DDA方法，利用運(yùn)行于Linux平臺(tái)上的ddscat7.1模擬計(jì)算軟件，計(jì)算得到兩種情況下的含核橢球顆模型的光散射Mueller矩陣，將矩陣元F11（θ）、F12（θ）、F33（θ）和 F34（θ）分別歸一化為F11（θ）／F11（θ）、F12（θ）／F11（θ）、F33（θ）／F11（θ）和F34（θ）／F11（θ）。然后利用 Matlab仿真計(jì)算軟件運(yùn)行模擬退火反演程序?qū)w一化矩陣元進(jìn)行反演，獲得含核橢球顆粒模型尺寸參數(shù)的反演值x。根據(jù)aeff＝λx／2π計(jì)算得出含核橢球顆粒模型的有效半徑aeff，并計(jì)算得出aeff的相對(duì)誤差。表1表示的是n為定值1.57＋0.42i，q的變化對(duì)粒徑反演結(jié)果和誤差的影響。表2表示的是q為定值0.2，n的變化對(duì)粒徑反演結(jié)果和誤差的影響。

表2 n對(duì)粒徑反演結(jié)果和誤差的影響Table 2 The influence of nto the result and error of the inversion of particle size

2.1 不同q的粒徑反演結(jié)果與誤差分析

圖2和圖3分別給出了有效半徑aeff0為2.0143um的含核橢球顆粒模型在n為定值1.57＋0.42i時(shí)，反演半徑aeff及其誤差隨q變化的關(guān)系曲線。由圖2和圖3分析知：當(dāng)q值小于等于0.4時(shí)，反演半徑與設(shè)定尺寸半徑之差的絕對(duì)值較小，即aeff的值緊靠aeff0，對(duì)應(yīng)的誤差均小于0.5%；當(dāng)q值較大時(shí)，反演半徑與設(shè)定尺寸半徑之差的絕對(duì)值比較大，aeff的值出現(xiàn)了較大的浮動(dòng)，對(duì)應(yīng)的誤差比較大，均超過0.5%。由此可以說明，對(duì)于軸對(duì)稱的含核橢球顆粒，q值較小時(shí)，基于Lorenz－Mie散射理論反演較為精確，反演結(jié)果較好符合仿真模擬實(shí)驗(yàn)的預(yù)期結(jié)果，即基于Lorenz－Mie散射理論激光粒度儀可以較為精確的探測(cè)含核顆粒的粒徑；當(dāng)q值較大時(shí)，基于Lorenz－Mie散射理論反演的結(jié)果出現(xiàn)較大誤差，已不符合精度要求，此時(shí)利用基于Lorenz－Mie散射理論激光粒度儀探測(cè)含核顆粒的粒徑出現(xiàn)較大誤差。

2.2 不同n的粒徑反演結(jié)果與誤差分析

圖4和圖5分別給出了有效半徑aeff0為2.0143um的含核橢球顆粒模型在內(nèi)核半徑與外殼有效半徑比q為定值0.2時(shí)，反演半徑aeff及其誤差隨n變化的關(guān)系曲線。從圖4和圖5可以看出，隨著內(nèi)核折射率n從1.65＋0.42i開始不斷增加，粒徑反演值逐漸偏離設(shè)定值aeff0，誤差逐漸增大，由此看出，內(nèi)核折射率n的改變?cè)谝欢ǔ潭壬嫌绊戭w粒的粒徑反演；總體上，反演半徑與設(shè)定尺寸半徑之差的絕對(duì)值均較小，反演誤差均小于0.5%，說明Lorenz－Mie散射理論適用于這種情況的粒徑反演，反演結(jié)果較好符合仿真模擬實(shí)驗(yàn)的預(yù)期結(jié)果，即說明含核顆粒的內(nèi)核折射率對(duì)基于Lorenz－Mie散射理論的激光粒度儀的顆粒探測(cè)結(jié)果在較小的誤差范圍內(nèi)。

圖2 不同q的半徑反演值Fig.2 Radius of inversion values of differentq

圖3 aeff的相對(duì)誤差Fig.3 The relative error of aeff

圖4 不同n的半徑反演值Fig.4 Radius of inversion values of differentn

圖5 aeff的相對(duì)誤差Fig.5 The relative error of aeff

3 結(jié)論

本文采用隨機(jī)取向、軸對(duì)稱的含核橢球顆粒模型模擬了由球形碳煙顆粒作為內(nèi)核、水為外殼的水凝物氣溶膠顆粒的光散射特性。利用DDA方法計(jì)算其Mueller矩陣，采用全局搜索能力很強(qiáng)的模擬退火算法，同時(shí)使用4個(gè)歸一化矩陣元數(shù)據(jù)進(jìn)行粒徑反演，實(shí)現(xiàn)了對(duì)單分散系含核橢球顆粒模型的粒徑反演，分析了q及n的變化對(duì)粒徑反演結(jié)果和誤差的影響。驗(yàn)證了基于Lorenz－Mie散射理論的模擬退火反演算法對(duì)于含核橢球顆粒模型的粒徑反演在一定程度上的有效性與實(shí)用性，并驗(yàn)證了基于Lorenz－Mie散射理論的激光粒度儀對(duì)含核橢球顆粒的粒徑探測(cè)的適用性。從粒徑反演結(jié)果分析表明：內(nèi)核半徑與外殼有效半徑比q對(duì)含核橢球顆粒的粒徑反演結(jié)果及誤差影響較大，當(dāng)q值較大時(shí)，基于Lorenz－Mie散射理論反演的結(jié)果出現(xiàn)較大誤差，已不符合精度要求；內(nèi)核折射率n的改變雖在一定程度上影響顆粒的粒徑反演，但總體上，反演半徑與設(shè)定尺寸半徑之差的絕對(duì)值均較小，反演誤差均小于0.5%，說明Lorenz－Mie散射理論適用于這種情況的粒徑反演，反演結(jié)果較好符合仿真模擬實(shí)驗(yàn)的預(yù)期結(jié)果。

［1］孫賢明，等.隨機(jī)取向雙層橢球顆粒偏振散射特性研究［J］.物理學(xué)報(bào)，2011，60（11）：378－384.

［2］韓一平，孫賢明.冰水混合云對(duì)可見光的吸收和散射特性［J］.物理學(xué)報(bào)，2006，55（2）：682－687.

［3］Quirantes A.Light scattering properties of spheroidal coated particles in random orientation［J］.Journal of Quantitative Spectroscopy and Radiative Transfer，1999，63（2）：263－275.

［4］汪宏七，等.大氣氣溶膠對(duì)用紅外窗區(qū)通道遙感陸面溫度的影響［J］.紅外與毫米波學(xué)報(bào)，2001，20（4）：275－278.

［5］陰俊燕，等.影響激光大氣傳輸因素分析［J］.紅外與激光工程，2008，37（3）：399－404.

［6］黃朝軍，等.大氣氣溶膠顆粒散射相函數(shù)的數(shù)值計(jì)算［J］.紅外與激光工程，2012，41（3）：580－585.

［7］孫賢明，等.含有密集隨機(jī)分布內(nèi)核的橢球顆粒光散射特性研究［J］.物理學(xué)報(bào)，2009，58（9）：6222－6226.

［8］房玉東，等.細(xì)水霧作用下火災(zāi)煙氣典型特性參數(shù)的變化規(guī)律［J］.燃燒科學(xué)與技術(shù)，2008，14（4）：372－377.

［9］房玉東.細(xì)水霧與火災(zāi)煙氣相互作用的模擬研究［D］.合肥：中國科學(xué)技術(shù)大學(xué)，2006.

［10］顏金培，等.應(yīng)用蒸汽相變機(jī)理脫除燃煤可吸入顆粒物實(shí)驗(yàn)研究［J］.中國電機(jī)工程學(xué)報(bào)，2007，27（35）：12－16.

［11］張艷春，等.基于米氏散射理論的粒度測(cè)試算法研究［J］.國外電子測(cè)量技術(shù)，2009，28（11）：24－25.

［12］李耀東，等.基于模擬退火的火災(zāi)煙顆粒粒徑反演［J］.中國激光，2011，38（1）：209－214.

［13］Mishchenko MI，et al.Scattering，absorption，and emission of light by small particles［M］.UK：Cambridge University Press，2002，115－139.

［14］張青，等.火災(zāi)煙顆粒分形模型和球形模型光散射的比較研究［J］.物理學(xué)報(bào)，2010，59（10）：7442－7446.

［15］Purcell EM，Pennypacker CR.Scattering and absorption of light by nonspherical dielectric grains［J］.Astrophysical Journal，1973，186：705－714.

［16］Draine BT.The discrete－dipole approximation and its application to interstellar graphite grains［J］.Astrophysical Journal，1988，333：848－872.

［17］Draine BT，F(xiàn)latau PJ.Discrete dipole approximation for scattering calculations［J］.Journal of the Optical Society of America A，1994，11（4）：1491－1499.

［18］饒瑞中.隨機(jī)取向立方顆粒光散射的數(shù)值分析［J］.物理學(xué)報(bào)，1998，47（11）：1790－1797.

［19］李明.詳解MATLAB在最優(yōu)化計(jì)算中的應(yīng)用（上卷）［M］.北京：電子工業(yè)出版社，2011，397－399.