基于模拟退火的随机取向含核椭球颗粒粒径反演

陈 都，邓小玖，李耀东，王 锋，张永明＊

（1.合肥工业大学电子科学与应用物理学院，合肥，230009；2.中国科学技术大学火灾科学国家重点实验室，合肥，230026）

0 引言

大气中悬浮颗粒物大多数为非球形，且并非仅由单种匀质成分组成，常见的含核颗粒即是其中一种，例如细菌、雾滴、城市水凝物气溶胶颗粒［1－3］，以及细水雾灭火抑制烟气和湿式洗涤塔除尘过程中生成的水烟凝并颗粒等。随着工业文明的进程，总体上城市大气环境中可吸入颗粒物逐渐增多，对人类生活环境和大气气候产生了影响，在一定程度上危害了人体健康及人类生活质量［4－7］。所以，控制并减少可吸入颗粒物的产生和排放是当下迫在眉睫的重大任务。

研究表明，细水雾对火灾烟气具有一定的抑制作用，在细水雾灭火抑制烟气过程中，粒径较小的雾滴与悬浮在空气中的火灾烟颗粒发生凝聚、合并，形成以碳烟颗粒作为凝结核、水为外壳的含核颗粒［8，9］；工业脱除燃烧颗粒物领域中，碳烟颗粒通过湿式洗涤塔，与塔腔内的过饱和水蒸汽凝聚、合并，形成以碳烟颗粒作为凝结核、水为外壳的含核颗粒。以上两种过程中形成的含核颗粒在重力和空气阻力的双重作用下外壳近似椭球形状。湿式洗涤塔烟颗粒脱除实验研究表明：碳烟颗粒粒径直接影响湿式洗涤塔的脱除效果［10］。对这类含核椭球颗粒粒径和形状准确的检测，可对细水雾灭火抑制烟气和湿式洗涤塔颗粒脱除过程进行更好的控制，以期实现更有效地灭火和对碳烟颗粒更好的脱除效果。

目前市场上大多数激光粒度仪都是一种基于Lorenz－Mie散射理论由颗粒的散射光强来反演颗粒粒径的探测仪器，而Lorenz－Mie散射理论是基于球形模型的理论［11］。本文将通过仿真模拟实验研究分析基于Lorenz－Mie散射理论的激光粒度仪对含核椭球颗粒的粒径探测的适用性。

粒径是颗粒的重要物理参数。在颗粒光散射过程中，散射光的角分布会随着颗粒粒径变化而变化，因此对散射光信号数据进行反演可获取颗粒粒径的信息［12］。本文以含核椭球颗粒为模型，模拟了含球形碳烟核的椭球状水滴颗粒光散射特性。利用DDA对随机取向含核椭球颗粒模型进行数值模拟计算得出Mueller矩阵，并基于Lorenz－Mie散射理论，采用全局搜索能力很强的模拟退火算法，实现了单分散系含核椭球颗粒的粒径反演。

1 理论与方法

1.1 Mueller矩阵

对于任意颗粒的光散射，入射光与散射光Stokes向量的变换关系表示为［13］：

式中I表示光强，Q，U，V表示光波的偏振状态，矩阵 ｛Fij（θ）（i，j＝1，2，3，4）｝为光散射 Mueller矩阵，包含6个独立变量，描述光散射从初态到待测态的变化关系，Mueller矩阵元随散射角的分布可以完整地描述散射光的强度与偏振特征［14］。

处理数据时，通常对Mueller矩阵进行归一化：矩阵元F11（θ）对应于00散射角的值归一化为F11（θ）／F11（0），其他矩阵元 Fij（θ）对应于 F11（θ）归一化为 Fij（θ）／F11（θ）。归一化的 Mueller矩阵与光强绝对值大小无关，能更好地反映光散射过程的物理特征［12，13］。

1.2 DDA方法

DDA方法可以计算任意形状颗粒模型的光散射和吸收特性［15，16］。DDA算法的基本思想是利用有限个离散的偶极子近似任意形状的颗粒，偶极子对局域电场的响应可以获得偶极矩，偶极子间通过各自产生的电场相互作用。

设将散射体离散为N个偶极子，各偶极子的坐标为rj（j＝1，2，···，N），极化率为αj，各偶极子的极化强度为Pj＝αjEj，其中Ej是rj处的电场，Ej是入射场Einc，j和其它N－1个偶极子所激发的电场的叠加：

应用DDA算法时必须满足以下两个条件［17，18］：

2）d必须足够小（N必须足够大），方能较好地描述目标颗粒的形貌特征。

1.3 模拟退火算法

模拟退火算法（Simulate Anneal Arithmetic，SAA）的思想最早由 Metropolis等提出，1983年Kirkpatrick S.，Gelatt C.D.和Vecchi M.P.将其成功用于组合优化。SAA是基于Monte Carlo迭代的一种全局概率型搜索算法，用来在一个大搜寻空间内寻找命题的最优解，其基本思想源于物理学中金属的退火过程。SAA算法流程如图1所示［19］：

SAA用一组冷却进度表参数来控制算法的进程，控制参数“降温”并趋于零时，最终求得组合优化问题的相对全局最优解。冷却进度表参数包括控制参数T的初始值T0、衰减函数、终止准则及Markov链长度Lk，优化问题的一个解i及其目标函数f（i）分别与固体的一个微观状态i及其能量E（i）相对应。从初始值X0和初始温度T0开始，对当前解重复进行“在解空间中随机产生新解 → Metropolis准则判断 → 接受或舍弃新解”的迭代，并按照一定的规则逐步降低T值，算法迭代终止时的当前解既是近似最优解。

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

2 计算结果与分析

本文采用匀质球体作为内核、水为外壳的轴对称的含核椭球颗粒模型，椭球纵横比（ε＝b／a，其中旋转轴的长度为a，短轴为b）设定为0.8［1］，入射光的真空波长为632.8nm，水的真空相对折射率为1.33，颗粒尺寸参数x0（x＝2πaeff／λ）设定为20，即其有效半径（有效半径定义为颗粒物的等体积 球的半径，即aeff＝ （3V／4π）1／3）aeff0为2.0143um。本文利用DDA方法模拟计算了两种情况下的含核椭球颗粒模型的光散射特性，并讨论分析了两种情况下的粒径反演结果和误差。

第一种情况是设定含核椭球颗粒模型的内核真空相对折射率n为1.57＋0.42i，讨论内核半径与外壳有效半径比q＝rcore／rshell＝rcore／aeff（变量q自0至1.0以0.1间隔递增）的变化对粒径反演结果和误差的影响。

第二种情况是设定含核椭球颗粒模型的内核半径与外壳有效半径比q为0.2，讨论内核真空相对折射率n（变量n自1.55＋0.42i至2.00＋0.42i，实部以0.05间隔递增）的变化对粒径反演结果和误差的影响。

应用DDA方法，利用运行于Linux平台上的ddscat7.1模拟计算软件，计算得到两种情况下的含核椭球颗模型的光散射Mueller矩阵，将矩阵元F11（θ）、F12（θ）、F33（θ）和 F34（θ）分别归一化为F11（θ）／F11（θ）、F12（θ）／F11（θ）、F33（θ）／F11（θ）和F34（θ）／F11（θ）。然后利用 Matlab仿真计算软件运行模拟退火反演程序对归一化矩阵元进行反演，获得含核椭球颗粒模型尺寸参数的反演值x。根据aeff＝λx／2π计算得出含核椭球颗粒模型的有效半径aeff，并计算得出aeff的相对误差。表1表示的是n为定值1.57＋0.42i，q的变化对粒径反演结果和误差的影响。表2表示的是q为定值0.2，n的变化对粒径反演结果和误差的影响。

表2 n对粒径反演结果和误差的影响Table 2 The influence of nto the result and error of the inversion of particle size

2.1 不同q的粒径反演结果与误差分析

图2和图3分别给出了有效半径aeff0为2.0143um的含核椭球颗粒模型在n为定值1.57＋0.42i时，反演半径aeff及其误差随q变化的关系曲线。由图2和图3分析知：当q值小于等于0.4时，反演半径与设定尺寸半径之差的绝对值较小，即aeff的值紧靠aeff0，对应的误差均小于0.5%；当q值较大时，反演半径与设定尺寸半径之差的绝对值比较大，aeff的值出现了较大的浮动，对应的误差比较大，均超过0.5%。由此可以说明，对于轴对称的含核椭球颗粒，q值较小时，基于Lorenz－Mie散射理论反演较为精确，反演结果较好符合仿真模拟实验的预期结果，即基于Lorenz－Mie散射理论激光粒度仪可以较为精确的探测含核颗粒的粒径；当q值较大时，基于Lorenz－Mie散射理论反演的结果出现较大误差，已不符合精度要求，此时利用基于Lorenz－Mie散射理论激光粒度仪探测含核颗粒的粒径出现较大误差。

2.2 不同n的粒径反演结果与误差分析

图4和图5分别给出了有效半径aeff0为2.0143um的含核椭球颗粒模型在内核半径与外壳有效半径比q为定值0.2时，反演半径aeff及其误差随n变化的关系曲线。从图4和图5可以看出，随着内核折射率n从1.65＋0.42i开始不断增加，粒径反演值逐渐偏离设定值aeff0，误差逐渐增大，由此看出，内核折射率n的改变在一定程度上影响颗粒的粒径反演；总体上，反演半径与设定尺寸半径之差的绝对值均较小，反演误差均小于0.5%，说明Lorenz－Mie散射理论适用于这种情况的粒径反演，反演结果较好符合仿真模拟实验的预期结果，即说明含核颗粒的内核折射率对基于Lorenz－Mie散射理论的激光粒度仪的颗粒探测结果在较小的误差范围内。

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

图3 aeff的相对误差Fig.3 The relative error of aeff

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

图5 aeff的相对误差Fig.5 The relative error of aeff

3 结论

本文采用随机取向、轴对称的含核椭球颗粒模型模拟了由球形碳烟颗粒作为内核、水为外壳的水凝物气溶胶颗粒的光散射特性。利用DDA方法计算其Mueller矩阵，采用全局搜索能力很强的模拟退火算法，同时使用4个归一化矩阵元数据进行粒径反演，实现了对单分散系含核椭球颗粒模型的粒径反演，分析了q及n的变化对粒径反演结果和误差的影响。验证了基于Lorenz－Mie散射理论的模拟退火反演算法对于含核椭球颗粒模型的粒径反演在一定程度上的有效性与实用性，并验证了基于Lorenz－Mie散射理论的激光粒度仪对含核椭球颗粒的粒径探测的适用性。从粒径反演结果分析表明：内核半径与外壳有效半径比q对含核椭球颗粒的粒径反演结果及误差影响较大，当q值较大时，基于Lorenz－Mie散射理论反演的结果出现较大误差，已不符合精度要求；内核折射率n的改变虽在一定程度上影响颗粒的粒径反演，但总体上，反演半径与设定尺寸半径之差的绝对值均较小，反演误差均小于0.5%，说明Lorenz－Mie散射理论适用于这种情况的粒径反演，反演结果较好符合仿真模拟实验的预期结果。

［1］孙贤明，等.随机取向双层椭球颗粒偏振散射特性研究［J］.物理学报，2011，60（11）：378－384.

［2］韩一平，孙贤明.冰水混合云对可见光的吸收和散射特性［J］.物理学报，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］汪宏七，等.大气气溶胶对用红外窗区通道遥感陆面温度的影响［J］.红外与毫米波学报，2001，20（4）：275－278.

［5］阴俊燕，等.影响激光大气传输因素分析［J］.红外与激光工程，2008，37（3）：399－404.

［6］黄朝军，等.大气气溶胶颗粒散射相函数的数值计算［J］.红外与激光工程，2012，41（3）：580－585.

［7］孙贤明，等.含有密集随机分布内核的椭球颗粒光散射特性研究［J］.物理学报，2009，58（9）：6222－6226.

［8］房玉东，等.细水雾作用下火灾烟气典型特性参数的变化规律［J］.燃烧科学与技术，2008，14（4）：372－377.

［9］房玉东.细水雾与火灾烟气相互作用的模拟研究［D］.合肥：中国科学技术大学，2006.

［10］颜金培，等.应用蒸汽相变机理脱除燃煤可吸入颗粒物实验研究［J］.中国电机工程学报，2007，27（35）：12－16.

［11］张艳春，等.基于米氏散射理论的粒度测试算法研究［J］.国外电子测量技术，2009，28（11）：24－25.

［12］李耀东，等.基于模拟退火的火灾烟颗粒粒径反演［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］张青，等.火灾烟颗粒分形模型和球形模型光散射的比较研究［J］.物理学报，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，Flatau PJ.Discrete dipole approximation for scattering calculations［J］.Journal of the Optical Society of America A，1994，11（4）：1491－1499.

［18］饶瑞中.随机取向立方颗粒光散射的数值分析［J］.物理学报，1998，47（11）：1790－1797.

［19］李明.详解MATLAB在最优化计算中的应用（上卷）［M］.北京：电子工业出版社，2011，397－399.