基于机器学习的椭球颗粒群消光系数预测
1 引言
颗粒群的粒径分布表征在大气、化工及医药等领域具有重要意义[1-4],目前已开发了多种粒径表征技术,其中,光谱消光技术具有测量范围大、非接触、可在线测量等优点,应用前景广阔。但与其他直接测量技术不同,光谱消光技术需根据不同波长的消光系数反演颗粒群的粒径参数,计算过程复杂,且需针对参数反演的严重病态问题设计算法以减小噪声干扰。
目前粒径参数反演算法可分为独立及非独立两种模式,独立模式无需颗粒群的先验知识,直接离散化Fredholm积分方程将参数反演简化为线性方程组的迭代求解,效率较高,理论上可反演任意颗粒群的粒径分布,求解方法包括最小二乘法,Tikhonov正则化算法[5-6]及Chahine迭代算法[7]等。但粒径分辨率有限,且通常需加入正则项保证粒径分布的平滑变化。非独立模式假设颗粒群粒径分布满足Rosin-Rammler(R-R)分布、正态分布或对数正态分布等规律[8-9],进一步基于遗传算法或粒子群优化等进化算法优化粒径参数以最小化测量与模拟消光系数偏差,具有很好的抗干扰性及全局优化能力。当实际粒径分布与假设相符时,可准确反演粒径参数,但需多次求解散射问题,效率较低,实时性差,因此提高优化效率是加强非独立模式实用性的核心问题。
进化算法的总优化时间在串行编程时近似为单次消光系数求解时间与总求解次数的乘积,其中,总求解次数与进化算法的收敛速率相关,本研究主要关注的是单次求解效率的提高问题。Sun等[10]用多项式拟合消光系数与粒径参数之间的函数关系,取得了较好的效果,但拟合过程与粒径分布的具体形式相关,表达形式不具有普适性。机器学习技术具有强大的非线性映射能力,近年来得到了广泛关注。Nascimento等[11-13]将机器学习技术应用于粒径反演领域,取得了一定的效果。这些技术一般将神经网络的输入及输出分别设置为光谱消光系数及粒径分布参数,通过预先建立的训练数据集更新神经网络参数后直接根据消光系数获得粒径分布参数,效率较高。但神经网络的参数不具备物理意义,难以调试,导致该方法抗干扰性能差。针对该问题,本文以粒径分布参数为输入,消光系数为输出,利用机器学习效率高的优势,缩短单次计算时间。由于正向建模与测量噪声无关,因此,后续可进一步结合进化算法实现粒径分布参数的准确高效反演。
2 基本原理
均匀球形颗粒的消光系数可由Mie散射理论严格求解获得,但实际测量中颗粒可能为椭球形状,导致Mie散射理论不再适用。可通过T矩阵法[14]、有限差分法[15]、有限元法[16]等多种数值计算方法准确求解不规则颗粒的电磁散射问题,但计算复杂、效率较低,难以满足实时反演粒径参数的需求。反常衍射近似理论[17-18]简化了消光系数的求解过程、计算效率高且原理简单,能适用于弱散射椭球颗粒。因此分别用Mie散射及反常衍射近似理论建立均匀球形及椭球形颗粒群消光系数的训练与测试数据集。
2.1 反常衍射近似
当椭球颗粒群浓度(单位体积内颗粒的数目)较低且不存在多重散射时,根据朗伯比尔定律可将颗粒群的消光系数τ(λ)由第一类Fredholm积分方程表示为
式中,N0为颗粒数目,λ为照明波长。
式中,下标ad、edge分别为椭球体及边缘,θ为光束与旋转轴之间夹角,w = 2ka(m-1)/g, g= (cos2θ+ε2sin2θ)1/2,k=2p/λ,c0为与折射率相关的常数,当颗粒折射率较小时,c0约为0.996163。S为颗粒在沿光束传播方向的投影面积,s为S边缘的弧长,P为颗粒在光束传播方向投影的边界,R为颗粒边缘处的曲率半径。
为了简化计算,一般假设颗粒群的粒径与形状分布相互独立[19],将f(ε, b)改写为
式中,fε(ε)及fb(b)分别为颗粒群的形状及粒径分布。
文献[ 20]将矿物气溶胶颗粒群的形状分布fε(ε)表示为
式中,n为形状分布调节参数,C为归一化因子,以确保fε(ε)的积分为1。
实验测试了两种常见的单峰粒径分布函数,包括
式中,u、σ为尺寸分布的两个特征参数,u为粒径分布均值,σ可以表示粒径的离散程度。
将(2)式~(10)式代入(1)式,可得到任意波长时颗粒群的消光系数。使用进化算法反演颗粒群粒径与形状分布参数向量F=[u, σ, n]时,需要根据一定规则不断改变F,减小测量与理论计算消光系数的误差。由(1)式可知,虽然
2.2 机器学习
由于消光系数与向量F之间为复杂的非线性函数关系,基于多项式拟合易出现欠拟合或过拟合现象,基于神经网络的机器学习方法具有优良的泛化能力,因此用神经网络技术拟合两者之间的函数关系。目前已开发了多种神经网络类型及训练算法,尽管深度学习近年来应用广泛,但其网络规模庞大、运算复杂,更适合于大数据。因此,实验仍使用传统多层感知器(MLP)结构,如
基于独立模式求解时,为了将(1)式简化为线性方程组形式,一般需假设形状分布参数n已知,如文献[
21]中假设n=0,可表示不同纵横比椭球均匀分布,导致反演存在一定误差。进化算法则将n作为优化变量,输入向量包含形状及尺寸分布的三个参数,输出层为消光系数。由于神经网络训练及测试具有一定随机性,难以确定隐藏层的最优结构,因此,使用两层神经元数目相同的隐藏层。理论上可将波长作为输入参数,但为了提高精度,将数据集容量、网络规模及训练时间均以数量级程度增加。由于实际应用中波长数目一般有限,因此针对每个波长单独训练一个网络,并存储于内存中,后续直接调用,即通过增加空间复杂度降低时间复杂度。现代计算机的内存一般较大,而本神经网络的规模较小,同时存储多个网络数据并不消耗太多的额外资源。由于网络针对特定颗粒群粒径与形状分布及折射率进行训练,当采用其他粒径分布形式及颗粒材料时,需更新相应训练及测试数据集并重新训练测试网络。
模拟均在Matlab平台上实现,首先,根据颗粒群选择粒径及形状分布并初始化相应参数空间;然后,根据反常衍射近似理论生成训练及测试数据集,为避免在(1)式中频繁计算消光系数,可预先计算各离散点处的消光系数并存储在内存中,数值积分时直接使用以提高计算效率;最后,训练神经网络及优化网络参数,其中,隐藏层及输出层神经元的激活函数分别为tansig及purelin函数,并用测试数据集验证神经网络的预测精度。
3 结果与分析
分别测试了两种粒径分布时机器学习的预测精度,粒径及形状参数空间如
根据反常衍射近似理论针对无吸收颗粒生成了训练及测试数据集,由于每个波长的消光系数均单独使用一个神经网络进行训练预测,因此将颗粒折射率设置为常数1.33,实际应用中折射率可根据真实参数进行设置,不影响训练及预测过程。光谱仪的测量范围为300~950 nm,分辨率为50 nm,共测量14个波长的消光系数。(1)式中的积分参数(bmin, bmax)=[0.5 μm, 20 μm]、(εmin, εmax)=[0.2, 5]。由于消光系数的绝对值受N0及b的影响,因此实验主要关注相对变化趋势,并对消光系数进行归一化处理,即直接将消光系数除以训练集中该波长消光系数的最大值。采用误差反向传播算法优化神经网络参数,训练最大次数为5000次,优化目标为e-10。计算机的CPU为Intel(R) Core(TM) i7-7700 @3.60 GHz,内存为8 G,当隐藏层神经元数目为20时,训练时间约持续2.5 h。
表 1. 粒径与形状参数空间
Table 1. Particle size and shape parameter space
|
图 3. R-R分布下的颗粒群参数。(a)消光系数,(b)神经元数目对预测误差的影响
Fig. 3. Particle group parameters under R-R distribution. (a) Extinction coefficient; (b) influence of the neuron number on the prediction error
为了验证基于机器学习对消光系数预测效率的提升,计算了不同神经元数目时的预测时间。由于Matlab执行for或while循环时引入了额外时间消耗、效率较低,无法反映实际预测时间,因此采用向量运算原理,同时输入了100000组粒径参数,计算了平均预测时间。基于进化算法反演粒径参数时,同样可以采取相似的向量计算原理提高效率。
图 4. 神经元数目对机器学习预测时间的影响
Fig. 4. Influences of the neuron number on prediction time for machine learning
采用相同的训练及测试过程计算了正态分布颗粒群的消光系数,粒径及形状参数向量(u, σ, n)为(1.534 μm, 0.139 μm, 3.508),由于预测时间只与神经元的数目有关,因此未重复计算预测时间,结果如
图 5. 正态分布下的颗粒群。(a)消光系数;(b)神经元数目对预测误差的影响
Fig. 5. Particle group under normal distribution. (a) Extinction coefficient; (b) influence of the neuron number on the prediction error
4 结论
针对基于进化算法反演椭球颗粒群粒径及形状参数效率较低的问题,提出了一种基于机器学习的计算方法,根据反常衍射近似理论计算了不同粒径及形状参数时颗粒群光谱消光系数,并基于多层感知器神经网络建立了粒径及形状参数与消光系数之间的联系。数值模拟结果表明,当隐藏层神经元数目为20时,测试集的平均预测误差均低于0.05%,预测时间约为0.6 μs。若实际测量中无折射率先验信息,可在神经网络输入向量中加入对应元素,使用相同训练及测试过程预测光谱消光系数。该方法可以显著提高消光系数的计算效率,满足进化算法实时高精度反演参数的要求,有望在颗粒特性在线实时表征中得到广泛应用。
[7] Ferri F, Bassini A, Paganini E. Modified version of the Chahine algorithm to invert spectral extinction data for particle sizing[J]. Applied Optics, 1995, 34(25): 5829-5839.
[8] Li M Z, Wilkinson D. Particle size distribution determination from spectral extinction using evolutionary programming[J]. Chemical Engineering Science, 2001, 56(10): 3045-3052.
[9] He Z Z, Qi H, Yao Y C, et al. Inverse estimation of the particle size distribution using the fruit fly optimization algorithm[J]. Applied Thermal Engineering, 2015, 88: 306-314.
[10] Sun X G, Tang H, Dai J M. Retrieval of particle size distribution in the dependent model using the moment method[J]. Optics Express, 2007, 15(18): 11507-11516.
[17] Fournier G R. Evans B T N. Approximation to extinction efficiency for randomly oriented spheroids[J]. Applied Optics, 1991, 30(15): 2042-2048.
[18] Sun X G, Tang H, Yuan G B. Anomalous diffraction approximation method for retrieval of spherical and spheroidal particle size distributions in total light scattering[J]. Journal of Quantitative Spectroscopy and Radiative Transfer, 2008, 109(1): 89-106.
张晓浩, 陈功叶, 李浩淼, 彭浩辰, 曹兆楼. 基于机器学习的椭球颗粒群消光系数预测[J]. 激光与光电子学进展, 2021, 58(6): 0629001. Zhang Xiaohao, Chen Gongye, Li Haomiao, Peng Haochen, Cao Zhaolou. Prediction of Extinction Coefficient of Ellipsoid Particle Group Based on Machine Learning[J]. Laser & Optoelectronics Progress, 2021, 58(6): 0629001.