光谱消光法基于Lambert-Beer(LB)定律, 当一束光强为I₀的平行单色光, 穿过离散颗粒物两相体系会发生散射和吸收并导致透射光的衰减, 将入射光强I₀与透射光强I之比的对数形式称为消光值, 可以描述颗粒物对某一波长光吸收/散射的强弱与颗粒物浓度及光程的关系, 即

ln(I0/I)=πcnLR2kext=cnLCext(1)

式(1)中, c_n为颗粒数目浓度, R为被测颗粒半径, L为光程, k_ext为消光系数, C_ext为消光截面。

LB定律适用于单一颗粒系的单分散情况, 在不同波长条件下运用式(1), 可获得消光谱信息并确定颗粒系粒径分布。 不过, 对于混合颗粒体系, 该定律适用性遇到困难, 不同类型颗粒光散射特性不同, 其消光系数/截面也不一致, 在模型中很难表达并使得理论消光值无法计算。 此时, 引入蒙特卡罗方法, 构建新的消光谱预测模型。

根据消光法原理, 蒙特卡罗方法核心思想是将物理上的入射光束按照“ 光子” 概念抽象并做离散化处理, 通过光子与体系中颗粒的相互作用将光子历程分为透射、 散射与吸收, 分类统计各去向的光子数, 进而得到其消光值。 根据其特点, 易于引入到颗粒及混合颗粒多相体系的模拟计算。

如图1所示, 两相体系中混合了两种球形颗粒: 黑色颗粒Ⅰ 、 蓝色颗粒Ⅱ 。 几何参数L为样品池厚度(光程), d_T、 d_R分别为发射器与接收器的直径, 2H为样品池上下边界距离, S为发射器或者接收器与样品池的距离, l为随机散射自由程。 通过统计获取接收器光子数, 即可模拟在一定颗粒粒径、 浓度和光波长条件下混合颗粒系中光波行为并计算消光值。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 光子在混合颗粒体系中传播过程示意图Fig.1 Schematic of photon propagation through a mixed particle system

图2给出了蒙特卡罗模型计算流程, 根据输入初始参数颗粒半径R、 波长λ 、 体积浓度c_v, 进行消光值计算。 假设光束沿着图1中x轴方向传播, 平行光入射。 以发射器中心点为原点, 光子的初始出射坐标(x₀, y₀)

x0=0y0=(ε1-0.5)dT(2)

式(2)中, ε ₁为[0, 1]区间内服从均匀分布的随机数。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 蒙特卡罗模型算法流程Fig.2 Flow chart of Monte Carlo model

光子刚进入样品池时坐标为(S, y₀), 若光子和颗粒发生碰撞, 则光子首次发生散射的坐标(x₁, y₁)可以表示为

x1=S+ly1=y0(3)

l=-lnεlτ(4)

式中, ε _l是[0, 1]范围内均匀分布随机数, τ 为浊度, 其在混合颗粒系中, 仍可定义为^[12]

τ=cn×Cext(5)

其中, 对于混合颗粒体系, 需要判断颗粒类型

ε2> ϕ颗粒Ⅰε2≤ϕ颗粒Ⅱ(6)

式(6)中, ε ₂是在[0, 1]区间内服从均匀分布的随机数, 混合比ϕ 指颗粒Ⅱ 在整个混合颗粒体系中所占的数目比, 例如: ϕ =0时全为颗粒Ⅰ , ϕ =1时则全为颗粒Ⅱ 。 当前颗粒的消光截面C_ext可由Mie散射理论计算^[1]。 颗粒系数目浓度c_n为

cn=cv(1-ϕ)43πR13+ϕ43πR23(7)

通常进入两相体系的光子可能被颗粒吸收、 被接收器直接接收(透射)、 未被接收器接收(逃逸)或在样品池里再次散射, 定义散射截面和消光截面比值为反照率a

a=Csca/Cext(8)

其中: 消光截面C_ext和散射截面C_sca由Mie散射理论计算得出。 根据设定条件判断光子的下一事件

(ε3> a)(n≥1)吸收(x> S+L)(n< 1)透射(S< x< S+L)(|y|< H)再次散射其他逃逸(9)

式(9)中, ε ₃是[0, 1]区间服从均匀分布的随机数, n为散射次数。

对于多次散射, 光子与颗粒碰撞后的空间散射角分布可以由Henyey-Greenstein相函数^[13]确定, 散射角θ ₀的抽样表示为

θ0=arcos12g1+g2-1-g21-g+2gεθ02(10)

式(10)中, εθ0是另一个[0, 1]范围内随机数, g是不对称因子。 发生第一次散射时, 光子散射方向θ ₁=θ ₀。 重复式(4)—式(10), 根据散射自由程与散射角, 第n次散射时的光子散射方向为θ _n=θ _n-1+θ ₀(n> 1), 光子散射后的坐标为

xn+1=xn-lcosθnyn+1=yn+lsinθn(11)

根据式(9)条件判断并统计其最终去向, 进行下一个光子历程, 直到完成所有光子计算, 统计各个物理过程的光子数, 得到散射介质的消光特性。 消光谱可表示为

ln(I0/I)=ln(Nset/N)(12)

式(12)中, N为接收透射光子总数, N_set为设定光子数。

按照前述模型编制计算程序, 模型尺寸(参见图1)、 颗粒性质、 可变参数均列于表1中, 介质为水(折射率m=1.33), 颗粒Ⅰ 选取聚苯乙烯(折射率m₁), 颗粒Ⅱ 选取高密度玻璃(折射率m₂), 忽略颗粒的吸收特性。 根据前期初步数值结果分析, 设定光子数为10⁵进行计算。

表1

Table 1

表1(Table 1)

表1 计算参数 Table 1 Parameters in the calculation 模型尺寸 S/mm H/mm L/mm dT/mm dR/mm 45.0 15.0 10.0 1.0 1.0 颗粒性质 m1 R1/μm m2 R2/μm 1.591 {0.2, 0.8, 2.0} 1.910 {0.15, 0.2, 0.3, 0.8, 2.0} 可变参数 λ /μm ϕ cv [0.4, 0.8] [0, 1] 2× 10-5

表1 计算参数 Table 1 Parameters in the calculation

使用蒙特卡罗模型计算程序模拟R=0.2 μm的颗粒Ⅰ 的光子去向, 将经过一次与多次散射后被接收器接收到的光子分为单散射与复散射, 定义从前向边界出射的其他逃逸光子为前向逃逸, 反之为后向。

图3可以看出, 在给定粒径下, 随着波长增大透射光子数逐渐增多。 按光散射理论, 波长增大, 亚微米区颗粒的无因次参数α (α =2π R/λ )减小, 消光系数减小, 再结合式(4)计算出射光子随机自由程, 光子准直透射的概率增大。 在当前光学厚度下, 散射光较微弱, 且以单散射为主, 其随光波长增大而减小。 由于接收器位置设定离样品池较远, 接收器尺寸小, 最终进入接收器散射光有限(计算消光时限定只考虑透射光子), 逃逸发生的光子以前向为主, 数目随光波长增大逐渐减小。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 光子事件统计Fig.3 Statistics of photon events

为了验证蒙特卡罗方法预测消光谱的准确性, 首先分别考虑颗粒Ⅰ 和颗粒Ⅱ 的单一颗粒系情况, 采用蒙特卡罗方法仿真计算, 研究颗粒系的消光谱变化, 并和Lambert-Beer模型进行对比。

图4给出两种颗粒各自的消光谱, 颗粒半径分别为0.2、 0.8、 2.0 μm。 可以看出, 随着颗粒粒径的增大, 消光谱有明显变化, 表明其对于粒径是非常敏感, 有利于粒径的反演; 当粒径超出亚微米区后(如2.0 μm), 消光系数随着无因次参数α 的增大而呈现振荡现象; 而相同粒径下, 由于颗粒Ⅰ 和颗粒Ⅱ 自身折射率不同, 导致消光谱趋势不同。 通过结果对比看出蒙特卡罗方法预测与Lambert-Beer模型在数值上吻合, 相对误差均在± 2%以内, 表明蒙特卡罗方法适用于颗粒系的消光谱预测。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 蒙特卡罗方法与Lambert-Beer模型的消光谱对比 (a): 聚苯乙烯; (b): 高密度玻璃Fig.4 Comparison of extinction spectra of Monte Carlo and Lambert-Beer models (a): Polystyrene; (b): High-density glass

单一颗粒系体积浓度c_v与消光值存在线性关系[见式(1)], 其并不影响消光随光波长的变化趋势即归一化的消光谱, 对于混合颗粒系, 同样通过数值模拟进行了验证。 后续计算中, 仅以c_v=2× 10^-5进行计算, 不考虑c_v对消光谱的影响。

将模型拓展到由颗粒Ⅰ 与颗粒Ⅱ 组成的混合颗粒系中, R₁=0.2 μm不变, R₂不同时, 改变两种颗粒的数目混合比ϕ , 进行消光谱预测。

图5(a)、 (b)、 (c)给出了不同混合比条件下颗粒系的消光谱。 随着入射光波长的增大, 消光谱呈递减趋势, 颗粒Ⅱ 的粒径对消光谱曲线有明显影响。 当R₂为0.15 μm, 图5(a)中消光谱随波长的增大逐渐趋缓, 而当R₂增至0.3 μm, 由于消光系数差异, 图5(c)中消光谱随波长增大趋势更近于线性。 对于给定波长和粒径, 从图5(d)可以看出: 颗粒系的消光值随混合比增大而递增, 而波长减小时, 消光值由线性趋势向非线性趋势转变, 且两种颗粒消光特性差别越大, 非线性趋势越明显, 因为两种颗粒类型不同, 当入射光波长改变时, 颗粒的消光特性相应变化。 上述结果说明, 在一定浓度下混合颗粒系的消光由颗粒类型、 混合比、 颗粒粒径和光波长等共同决定, 如混合比对消光谱的影响与不同类型颗粒的散射特性差异有关, 其受制于粒径和波长取值, 二者同时影响颗粒的消光和散射效应, 据此可根据消光谱同步反演出颗粒粒径与混合比。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 混合颗粒系消光预测 (a): R1=0.2 μm, R2=0.15 μm; (b): R1=0.2 μm, R2=0.2 μm; (c): R1=0.2 μm, R2=0.3 μm; (d): λ =0.4 μm与λ =0.8 μmFig.5 Predicted light extinction of mixed particles (a): R1=0.2 μm, R2=0.15 μm; (b): R1=0.2 μm, R2=0.2 μm; (c): R1=0.2 μm, R2=0.3 μm; (d): λ =0.4 μm and λ =0.8 μm

最优化算法是消光谱反演求解颗粒粒径时最为常见方法^[14], 混合颗粒系可借助蒙特卡罗模型获得理论消光谱, 不过, 由于两种甚至多种颗粒类型的存在, 增加了混合比等待定参数, 反演难度会增加。 采用消光谱误差构建目标函数

F=∑j=1M[ln(I/I0)sim(R1, R2, ϕ, λj)-ln(I/I0)set, j]2(13)

式(13)中: ln(I/I₀)_sim是理论消光谱, ln(I/I₀)_set为设定值的消光谱, j为计算消光谱的波长数, j=1, 2, …, M。 在设定反演参数的上下限范围内寻优搜索最佳粒径与混合比, 使目标函数达到最小, 从而实现参数的反演。

为更好地验证反演效果, 采用三种全局最优算法对所构造的目标函数进行寻优, 分别是优化遗传算法(improved genetic algorithm, IGA)^[15]、 粒子群算法(particle swarm optimization, PSO)^[16]、 改进差分进化算法(improved differential evolution, IDE)^[17]。 其中, IGA具有较好全局收敛性, 经优化设定最大迭代次数为500、 种群尺度为50、 交叉概率为0.8、 变异概率为0.045, 以提高算法稳定性。 PSO采用自适应惯性权重, 对适应度值小于平均值的粒子进行小波变异, 增强群体多样性, 并使粒子在解空间的其他区域进行搜索, 防止粒子群算法陷入局部最优, 设定种群尺度为500, 迭代次数在50~60之间。 IDE引入了自适应控制参数因子对缩放因子和交叉因子进行优化, 选定交叉因子为0.85, 缩放因子上限和下限分别为0.5和0.3, 种群尺度为50, 并设定最大迭代次数为50, 以提高优化效率。

如图6所示, 进行颗粒和混合比的单参数、 双参数和三参数同步反演。 分别为: ①已知混合颗粒系粒径, 对混合比进行反演; ②已知颗粒系混合比, 对两种颗粒粒径进行反演; ③两种颗粒类型粒径相同或不同时, 对粒径与混合比进行同步反演。 分别选取R₁=R₂=0.2 μm、 ϕ =0.8和R₁=0.2 μm、 R₂=0.3 μm、 ϕ =0.4两种情形, 设置粒径参数上下限为R∈ [0.05 μm, 1 μm]、 混合比ϕ ∈ [0, 1], 计算波长数M=9。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 反演算例Fig.6 Inversion examples

图7(a)为已知颗粒粒径R₁和R₂, 求解其混合比, 可以看出: 无论颗粒粒径是否相同, 三种优化算法反演混合比与设定值相对误差均小于1.5%。 图7(b)给定颗粒系的混合比ϕ =0.4, 反演两种颗粒的粒径, 为双参数反演, 粒径相对误差在3%以内, 均较为精确。 为验证算法稳定性, 以IGA为例, 进行三次重复反演, 其粒径R₁、 R₂的标准差为8.2× 10^-4与1.1× 10^-3。 图7(c)、 (d)为三参数同步反演, R₁=R₂时, ϕ 的反演误差仍小于2%, 但R₁反演误差可至约8%(PSO算法); 当R₁≠ R₂时, 三参数的反演误差增大, R₁、 R₂的误差接近9%。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 混合颗粒系参数反演结果 (a): 混合比反演结果; (b): 粒径反演结果; (c): 粒径和混合比同步反演结果(R1=R2); (d): 粒径和混合比同步反演结果(R1≠ R2)Fig.7 Inversion results of parameters in mixed particle systems (a): Inversion of the mixing ratio; (b): Inversion of particle size; (c): Simultaneous inversion of particle size and mixing ratio (R1=R2); (d): Simultaneous inversion of particle size and mixing ratio (R1≠ R2)

此外, 三种算法的反演时间不尽相同, 双参数反演时, IGA和IDE的一次计算时间约需40 min, 而PSO反演耗时为其数倍。 总体上, IGA算法效果较好且更加稳定, 就本文算例来看, 所有算例下反演误差均在10%以内。