DOI: 10.11817/j.issn.1672-7207.2016.11.040

基于RQPSO的颗粒粒径分布反演算法

张彪¹，李舒²，许传龙¹，王式民¹

(1. 东南大学 能源热转换及其过程测控教育部重点实验室，江苏 南京，210096；

2. 南京市计量监督检测院科技发展部，江苏 南京，210037)

摘要：针对标准量子微粒群算法(BQPSO)，提出改进量子微粒群算法(IQPSO)和含正则化项的改进量子微粒群算法(RQPSO)，并将它们引入到粒径分布的反演中，利用光全散射法在独立模式下，通过测量可见光波段内不同波长下的光谱消光值反演几种粒径分布，其中正问题利用反常衍射近似(ADA)计算得到估计值，测量值则通过Mie理论计算得到。研究结果表明：与BQPSO相比，IQPSO在计算效率和稳定性上得到很大提升；在粒径分布的反演中，RQPSO提高IQPSO的维数极限，并具有更高的反演精度、稳定性和抗噪性，为粒径分布的反演提供一种新的方法。

关键词：粒径分布；量子微利群算法；正则化；光全散射法；独立模式

中图分类号：TP212         文献标志码：A         文章编号：1672-7207(2016)11-3922-07

Retrieval of particle size distribution based on RQPSO algorithm

ZHANG Biao¹, LI Shu², XU Chuanlong¹, WANG Shimin¹

(1. Key Laboratory of Energy Thermal Conversion and Control of Ministry of Education,

Southeast University, Nanjing 210096, China;

2. Technology Development Department, Nanjing Institute of Measurement and Testing Technology,

Nanjing 210037, China)

Abstract: An improved quantum behavior particle swarm optimization (IQPSO) and regularized quantum behavior particle swarm optimization (RQPSO) were developed based on basic quantum behavior particle swarm optimization (BQPSO). Furthermore, these three algorithms were introduced in retrieval of particle size distribution (PSD). Several types of PSDs were retrieved by measuring the spectral extinction values in the visible spectrum, which used total light scattering method under independent mode. In the direct problem, the anomalous diffraction approximation was used to calculate the estimation values, and Mie theory was used for measurement values. The results show that the efficiency and stability of the IQPSO algorithm was proved to be more greatly improved than the BQPSO algorithm. For the retrieval of these PSDs, the RQPSO algorithm has better performances on dimension limit, accuracy, stability and noise immunity than the IQPSO algorithm. Thus, this algorithm provides a new method for retrieving of PSDs.

Key words: particle size distribution; quantum behavior particle swarm optimization; regularization; total light scattering; independent mode

颗粒粒径及其尺寸分布是颗粒技术最重要的参数和技术指标之一，在许多情况下，颗粒粒径不仅直接影响到产品的性能和质量，而且与能源的高效利用、环境污染的防治、生产工艺的优化、人身健康的保障等等都密切相关^[1-3]，因此，准确、快速、稳定的颗粒粒径测量技术对国防、气象、煤矿、动力、化工等方面有重要意义。光散射法近30年来随着激光和计算机等技术的迅速发展而得到了普遍应用，是一种重要的燃烧诊断手段。它通过测量颗粒的散射光强来反演颗粒粒径分布，能够实现非接触式在线测量，具有适用性广、测量范围宽、响应快、易于实现自动化等特点^[4-6]。该方法的最大难点在于如何通过衰减光谱的反演获得准确的粒径分布，在理论上归结于第一类Fredholm积分方程的求解。这是一个典型的不适定问题，直接求解具有很大的困难，因此，粒径分布反演算法的研究受到了众多研究者的重视^[7-9]。光全散射粒径测量方法分为独立模式算法和非独立模式算法，独立模式算法事先无须假定粒径分布，通过求解离散线性方程组得到粒径分布；非独立模式算法需要假设被测颗粒系的颗粒尺寸分布满足某个已知的分布函数，通过一定的优化方法确定分布函数中的待定参数。事实上，在绝大多数实际应用中，往往不知道被测颗粒系的粒径分布规律，或者颗粒系的尺寸分布无法简单地用某个分布来描述，这就使非独立模式求得的结果变得不可靠^[10]。目前已经发展了多种独立模式下的粒径分布反演算法，每种算法都有各自的应用背景和局限性。传统算法如：Phillips-Twomey算法、Chahine迭代算法、共轭梯度法等，这些算法存在目标函数及导数计算复杂、抗噪声能力差及多峰粒径重建困难等问题。近年来，遗传算法、模拟退火法等智能算法也已用于解决粒径测量问题，与传统方法相比，智能算法具有很好的全局搜索和抗噪声能力，但进化速度慢，反演结果不稳定等缺点突出^[11-12]。本文作者利用光全散射法在独立模式下反演了几种粒径分布，其中正问题利用反常衍射近似(ADA)计算得到估计值，测量值则通过Mie理论计算得到；反问题方法采用量子微粒群算法，并在标准量子微粒群算法的基础上提出了改进提出改进量子微粒群算法，在目标函数中利用Markov随机场建立正则化项，提升了算法的反演能力，并进行了数值反演验证。

1  正问题模型

1.1  光全散射法测量原理

光全散射法以光的散射理论为基础，当一束光强为I₀，波长为的平行单色光照射到厚度为L的悬浮待测颗粒系时，由于颗粒对入射光的吸收和散射作用，穿过颗粒系透射光的光强将减弱。根据Lambert-Beer定律，如果假设颗粒系为服从一定粒径分布范围的多分散球形粒子系，并且颗粒间满足不相关单散射的条件(忽略多次散射效应)，则多分散球形颗粒系在波长为时的消光值可以表示为^[13]

    (1)

式中：I₀为入射激光强度；I为透射光强；为波长为的消光值，它可以通过人工实验测量得到；N为待测颗粒系的颗粒总数；为颗粒系的体积频度分布；D_max和D_min分别为颗粒粒径分布的上、下限；为消光系数，它是波长、介质复折射率m以及颗粒粒径D的函数。

将粒径分布区间划分为M个子区间后，积分项可以变成离散项求和：

    (2)

式中：c_j为数值积分系数；M为整个待测粒径范围 [D_min, D_max]内划分的子区间个数；D_j为各子区间的等效粒径。

由于消光系数可以通过Mie理论或者一些近似方法计算得到，则可以通过测量多个波长下的消光值来反演颗粒系的体积频率分布。

1.2  Mie理论计算公式

Mie理论是一种球形颗粒散射特性的严格计算方法，Mie散射利用光的电磁波性质，应用Maxwell方程对散射颗粒形成的边界条件进行求解，以得到光散射的物理量。可用于计算在单色平行光照射下任意尺寸和任意成分的球形颗粒的散射场^[1]。利用Mie理论推导得到的颗粒消光系数的计算公式为

       (3)

式中：，为无因次粒径参数；a_n和b_n为Mie系数，是与Bessel函数和Hankel函数有关的函数，

        (4)

        (5)

           (6)

           (7)

式中：为颗粒相对于周围介质的复折射率；和分别为半奇阶的第1类Bessel函数和第2类Hankel函数。

1.3  ADA计算公式

由于Mie理论公式计算繁琐，特别是计算中涉及到无穷级数求和问题，占用大量内存资源和计算时间，不适合作为在线检测的正问题求解模型。反常衍射近似(ADA)假设颗粒的消光主要是吸收以及透射光与衍射光之间的干涉引起的，它的计算速度要远快于Mie理论的计算速度。最初的ADA只能适用于满足＞＞ 1和＜＜    1条件的粒子，ZHAO等^[14]引入边界效应，提出了修正ADA算法，扩展了的适用范围。

对于球形粒子，标准的反常衍射近似衰减因子和边界效应的表达式如下：

    (8)

             (9)

式中：，。

当粒径为0.1~10 μm时，修正消光系数为

 (10)

2  反问题模型

2.1  标准量子微粒群算法(BQPSO)

量子微粒群算法是由SUN等^[15]提出。量子微粒群算法将PSO系统看成是一个量子系统，每个粒子具有量子行为，量子的状态由波函数决定，为粒子的位置的概率密度。在第t次迭代中，粒子i在N维搜索空间内以粒子的局部吸引因子为中心在领域内搜索。

用蒙特卡洛法模拟后，位置可由下式表示：

            (11)

         (12)

            (13)

            (14)

式中：R，R₁和R₂为[0,1]区间内服从均匀分布的随机数；C₁和C₂为加速系数；p_ij为粒子i的个体历史最优位置在第j维上的坐标；为群体历史最优位置在第j维上的坐标；为吸引扩散系数，在＜1.781时可以保证量子微粒群的全局收敛。

2.2  改进量子微粒群算法(IQPSO)

对于一般的优化过程来说，初始粒子位置相对分散有利于群体的寻优，在标准量子微粒群算法中，初始化中是利用随机数产生位置的，这样有可能在初始化的过程中，产生位置相近的粒子，本文将典型的logistic映射应用于产生混沌的信号，这样初始化得到的粒子位置会相对分散。

  (15)

式中：为混沌系数，当时，logistic映射将处于一个完全的混沌状态。

在标准量子微粒群算法中，局部吸引因子是1个非常重要的参数，1个不利的局部吸引因子容易让算法陷入局部收敛，因此，本文提出了一种服从正态分布的局部吸引因子，它以群体历史最优位置为平均值，以为标准偏差的正态分布随机数。可以通过下式来定义：

            (16)

如同其他群体智能优化算法一样，标准量子微粒群算法在进化后期容易丧失种群的多样性，这样也不利于优化，本文引进了基因变异机制，它可以帮助粒子逃出局部最优位置而增加种群的多样性。位置更新可由下式来确定：

   (17)

              (18)

式中：sgn( )为符号函数；为在区间内服从均匀分布的随机数；为服从标准正态分布的随机数；R₃为[0,1]区间内服从均匀分布的随机数；为变异因子，为1个小于0.5的正数。

2.3  含正则化项的量子微粒群算法(RQPSO)

当反演参数过多时，标准量子微粒群算法和改进量子微粒群算法容易达到维数极限，导致反演结果不稳定而失效，这时通常需要对目标函数进行光滑化处理。本文在改进量子微粒群算法的基础上，利用Markov随机场理论认为每个待测粒径范围[D_min, D_max]的子区间上的体积频率密度仅仅与它邻近子区间的体积频率密度有关，与其他子区间内的体积频率密度无关，整个待测粒径范围内的体积频率密度为一个Markov随机场，因而建立起目标函数的正则化项，其中目标函数的表达式如下：

            (19)

  (20)

式中：M为测量值；E为估计值；A为光顺矩阵；为向量的2-范数。

3  结果与讨论

3.1  正问题验证

为了验证ADA的可靠性，本文通过Mie理论的计算结果来进行检验。在可见光范围内，光全散射法的最佳粒径测量范围为0.1~10.0 μm^[1]，本文粒子的复折射率参照典型飞灰粒子，它的范围分别为和^[6]。因此，分别利用Mie理论和ADA计算复折射率为的粒子在波长为0.45 μm下消光因子和复折射率为的粒子在波长为0.73 μm下消光因子，计算结果如图1所示。

从图1可以看出：ADA的计算结果与Mie理论的计算结果十分吻合，但ADA的计算效率是Mie理论的1 758倍。因此，ADA适合作为粒径反演中正问题的求解方法，而Mie理论的计算结果可以作为粒径反演中的测量值。

3.2  函数反演

为了比较标准量子微粒群算法(BQPSO)和改进量子微粒群算法(IQPSO)的性能，文中通过选用Sphere函数、Rastrigin函数和Schaffer函数3种标准测试函数对这2个算法进行研究，其具体参数如表1所示。

图1  ADA与Mie理论的计算结果对比

Fig. 1  Comparison of calculation results by ADA and Mie theory

在这2种算法中粒子的总数都为N_p=50，最大迭代次数都为N_t=50 000，吸引扩散系数随着迭代从0.9线性下降到0.4，在改进量子微粒群算法中混沌系数为3.9，变异因子为0.05，本文中若未作特别说明，算法的参数设置保持不变。为了对比这3种算法的计算效率和鲁棒性，将每种算法都独立执行100次，计算结果如表2所示。从表2可以看出：随着测试函数维数的增加，这2种算法的迭代次数均增加，标准偏差也在增大；相比于标准量子微粒群算法而言，改进量子微粒群算法在计算效率和稳定性上得到了很大的提升，说明了改进量子微粒群算法在性能上具有很大的优势。在后续的算例中，本文都采用改进量子微粒群算法作为反问题的求解方法。

表1  3个测试函数的具体参数

Table 1  Parameters of three test functions

表2  在固定精度下3种算法执行100次的迭代次数

Table 2  Iterations of three algorithms in 100 independent runs under fixed accuracy

3.3  单峰分布反演

为了考察量子微粒群算法在粒径反演中的效果，本文利用改进量子微粒群算法(IQPSO)和含正则化项的改进量子微粒群算法(RQPSO)反演1个服从单峰R-R分布的颗粒系，表达式为

  (21)

已知颗粒系的粒径范围为0.01~1.00 μm，颗粒的复折射率为。首先将粒径范围均匀的分成10个子区间，通过在可见光波段内选择10个波段来测量粒子系的消光值，反演出10个子区间内的体积频率密度，反演结果如图2(a)所示，其中算法的搜索范围为。为了衡量反演结果的优劣，定义误差为

        (22)

式中：f_ret和f_ori分别为颗粒体积频率密度的反演值和真值。

它们的反演误差分别为和，从图2(a)可以看出：2种算法反演出的分布都基本复合分布规律；从误差水平上看，RQPSO算法比IQPSO算法的误差要小，说明加上正则化后，反演结果更接近真实分布。

为了考察2种算法的维数极限，将粒径分布范围的子区间划分成20份，反演的误差水平分别为和，反演结果如图2(b)所示。从图2(b)可以看出：IQPSO算法的反演结果已经严重偏离真实分布，而RQPSO算法的反演结果能很好地符合真实分布。

为了考察RQPSO算法能力的提升程度，将粒径分布范围的子区间分别划分成30份和50份，它们的反演误差水平分别为和，反演结果分别如图2(c)~(d)所示。从图2可以看出：反演的分布都基本符合真实分布。说明RQPSO的反演能力得到了很大提升，不仅提高了反演的维数极限，而且提高了反演的精度。

3.4  双峰分布反演

为了进一步的考察含正则化项的改进量子微粒群算法的性能，本文利用其反演了1个服从双峰R-R分布的颗粒系，表达式为

图2  单峰分布的反演结果

Fig. 2  Retrieval results of monomodal distribution

 (23)

颗粒系的粒径范围和复折射率与上述颗粒系相同，分别将粒径范围划分成20份和30份，它们的反演误差水平分别为和，反演结果如图3所示。从图3可以看出：反演的分布都基本符合真实分布，进一步证明了RQPSO算法的可靠性。

3.5  误差分析

为了考察测量误差对反演结果的影响，本文以单峰分布为例，将粒径范围划分成30份，分别加入0，2%，5%和10%的测量误差后，反演误差水平分别为，，和，反演结果如图4所示。从图4可以看出：随着测量误差的增大，反演的误差水平有所增加，但反演的分布基本符合真实分布，说明了含正则化的改进量子微粒群算法具有很好的抗噪性能。

图3  双峰分布的反演结果

Fig. 3  Retrieval results of bimodal distribution

图4  含测量误差时的反演结果

Fig. 4  Retrieval results with measurement error

4  结论

1) 针对标准量子微粒群算法，提出了改进量子微粒群算法和含正则化项的改进量子微粒群算法，通过函数反演证明了改进量子微粒群算法相对于标准量子微粒群算法具有更高的效率和稳定性。

2) 将改进的量子微粒群算法引入到粒径分布的反演当中，利用光全散射法在独立模式下反演了服从单峰分布和双峰分布的颗粒系粒径分布，得出含正则化项的改进量子微粒群算法提高了改进量子微粒群算法的维数极限，并具有更高的精度和可靠性。

3) 含正则化项的改进量子微粒群算法具有很好的抗噪性能，为独立模式下的粒径分布反演提供了一个新的方法。

参考文献：

[1] 唐红. 光全散射法颗粒粒径分布反演算法的研究[D]. 哈尔滨: 哈尔滨工业大学电气工程及自动化学院, 2008: 1-8.

TANG Hong. Study of inversion algorithm of particle size distribution using total light scattering method[D]. Harbin: Harbin Institute of Technology. School of Electrical Engineering and Automation, 2008: 1-8.

[2] QI Hong, RUAN Liming, WANG Shenggang, et al. Application of multi-phase particle swarm optimization technique to retrieve the particle size distribution[J]. Chinese Optics Letters, 2008, 6(5): 346-349.

[3] TANG Hong. Retrieval of spherical particle size distribution with an improved Tikhonov iteration method[J]. Thermal Science, 2012, 16(5): 1400-1404.

[4] 菅立川, 齐宏, 阮立明, 等. 微粒群算法反演粒径分布的研究[J]. 节能技术, 2007, 141(25): 7-13.

JIAN Lichuan, QI Hong, RUAN Liming, et al. Inverse analysis of particles size distribution by particle swarm optimizer algorithm[J]. Energy Conservation Technology, 2007, 141(25): 7-13.

[5] HE Zhenzong, QI Hong, WANG Yuqing, et al. Inverse estimation of the spheroidal particle size distribution using Ant Colony Optimization algorithms in multispectral extinction technique[J]. Optics Communications, 2014, 328(19): 8-22.

[6] QI Hong, NIU Chunyang, GONG Shuai, et al. Application of the hybrid particle swarm optimization algorithms for simultaneous estimation of multi-parameters in a transient conduction- radiation problem[J]. International Journal of Heat and Mass Transfer, 2015, 83(4): 428-440.

[7] 孙晓刚, 唐红, 原桂彬. 非独立模式算法下粒径分布反演及分类的研究[J]. 光谱学与光谱分析, 2008, 28(5): 1111-1114.

SUN Xiaogang, TANG Hong, YUAN Guibin. Study of inversion and classification of particle size distribution under dependent model algorithm[J]. Spectroscopy and Spectral Analysis, 2008, 28(5): 1111-1114.

[8] QI Hong, HE Zhenzong, GONG Shuai, et al. Inversion of particle size distribution by spectral extinction technique using the attractive and repulsive particle swarm optimization algorithm[J]. Thermal Science, 2015, 19(6): 2151-2160.

[9] HE Zhenzong, QI Hong, YAO Yuchen, et al. Inverse estimation of the particle size distribution using the fruit fly optimization algorithm[J]. Applied Thermal Engineering, 2015, 88(9): 306-315.

[10] 徐峰, 蔡小舒, 苏明旭, 等. 独立模式算法求解颗粒粒径分布的研究[J]. 中国激光, 2004, 31(2): 223-228.

XU Feng, CAI Xiaoshu, SU Mingxu, et al. Study of independent model algorithm for determination of particle size distribution[J]. Chinese Journal of Lasers, 2004, 31(2): 223-228.

[11] 王丽, 孙晓刚. 基于模式搜索的光谱消光粒径分布反演算法的研究[J]. 光谱学与光谱分析, 2013, 33(3): 618-622.

WANG Li, SUN Xiaogang. Research on pattern search method for inversion of particle size distribution in spectral extinction technique[J]. Spectroscopy and Spectral Analysis, 2013, 33(3): 618-622.

[12] HE Zhenzong, QI Hong, YAO Yuchen, et al. An effective inversion algorithm for retrieving bimodal aerosol particle size distribution from spectral extinction data[J]. Journal of Quantitative Spectrascopy and Radiative Transfer, 2014, 149(18): 117-127.

[13] QI Hong, ZHANG Biao, REN Yatao, et al. Retrieval of spherical particle size distribution using ant colony optimization algorithm[J]. Chinese Optics Letters, 2013, 11(11): 112901.

[14] ZHAO Jianqi, HU Yinqiao. Bridging technique for calculating the extinction efficiency of arbitrary shaped particles[J]. Applied Optics, 2003, 42(24): 4937-4945.

[15] SUN Jun, FENG Bin, XU Wenbo. Particle swarm optimization with particles having quantum behavior[C]// Pedro Isasi. CEC2004, Proceedings of the 2004 Congress on Evolutionary Computation. Piscataway N J, US: IEEE, 2004: 325-331.

(编辑  陈爱华)

收稿日期：2016-01-02；修回日期：2016-03-20

基金项目(Foundation item)：国家自然科学基金资助项目(51506030, 51376049)；江苏省自然科学基金资助项目(BK20150622)；国家质量监督检验检疫总局科技计划项目(2012QK176) (Projects(51506030, 51376049) supported by the National Natural Science Foundation of China; Project(BK20150622) supported by the Natural Science Foundation of Jiangsu Province; Project(2012QK176) supported by the State Administration of Quality Supervision, Inspection and Quarantine Science and Technology)

通信作者：许传龙，博士，教授，博士生导师，从事多相流测试方面的研究；E-mail: chuanlongxu@seu.edu.cn