粒径分布及光学性质的准确测量对于环境监测、清洁能源利用以及生物医疗等方面具有重要意义。近年来，参与性介质颗粒系粒径测量越来越朝着自动化和微小颗粒粒径测量的方向发展，研究内容大多数集中于测量的在线实现和参与性介质颗粒系粒径分布的反演两部分[5-6]。光全散射法测量颗粒系的粒径分布，在对数据进行处理时会碰到求解第一类Fredholm积分方程的问题。鉴于目前该方程尚不能全面理论求解，因此一般会用到反演的方法解决这一问题。光学常数包含实部n和虚部k两部分，分别表征参与性介质颗粒系的折射行为和吸收行为。即使是相同种类的参与性介质颗粒系，其光学常数也可能因为温度不同、粒子的几何尺寸不同等因素而不同。虽然参与性介质颗粒系的光学常数的求解不易，但是由于其在解决颗粒系辐射传输机理上具有重要的作用，国内外****从未间断对其数值快速获得的研究[7-8]。
对于非线性领域系统物性的测量，需要结合相应的反问题算法。群体智能优化算法可以在一定的条件下，通过随机的优化方式找到反问题的解。目前广泛用于各个领域的反问题智能优化算法包括微粒群算法、蚁群算法等[9-12]。2010年Tan和Zhu提出了一种新型智能优化算法——烟花算法[13]。烟花算法主要受到烟花爆炸现象的启发，其基本原理为将烟花以及烟花爆炸产生的火花看作反问题的一个可能解，然后通过烟花的爆炸行为产生火花，从而实现在问题的解空间内搜索反问题的最优解。由于烟花算法的高效性和鲁棒性，提出后迅速在各个领域引起了广泛关注，如产品设计以及多目标优化问题方面[14]。
本文提出了一种改进的烟花算法，并将其应用于颗粒系粒径分布以及光学常数的反演中。主要原理为通过测量脉冲激光辐照下参与性介质不同角度的时域光学响应，结合反问题算法，可以同时获得颗粒系的粒径分布以及光学常数。
1 正问题模型 本文研究的是平行光入射一维吸收、散射、非发射的颗粒系，如图 1所示。
 |
| R-反射率；T-透射率；z-样品厚度方向；θ-透射角 图 1 实验原理图 Fig. 1 Experimental schematic |
| 图选项 |
瞬态激光在一维参与性介质中的辐射传输过程可以通过求解式(1)得到相应的各个角度辐射信号[15]。
 | (1) |
式中：c0为光速；I为在z位置θ方向上的辐射强度；βλ和σsλ分别为颗粒系的消光系数和散射系数，下角标λ为入射激光的波长；散射相函数Φλ(θ′, θ)为从θ′方向入射的激光从θ方向射出去的概率。对于不同角度的透射及反射信号的测量，可以采用图 1(b)的实验装置。
参与性介质颗粒系的辐射特性常用吸收截面、散射截面和衰减截面来表示[16]：
 | (2) |
式中：G为入射方向散射体的几何投影面积；Qext、Qsca和Qabs分别为粒子的衰减因子、散射因子和吸收因子。参与性介质颗粒系的光谱吸收系数κλ和光谱散射系数σsλ能够通过式(3)和式(4)计算出来：
 | (3) |
 | (4) |
式中：D为粒子的直径；Caλ(D)和Csλ(D)分别为直径为D的粒子的吸收截面和散射截面；N(D)为直径为D的粒子数密度。对于球形粒子来说，吸收截面Caλ和散射截面Csλ可以通过Mie理论求得[15]。
实际的颗粒系粒径测量中发现，大多颗粒系的粒径分布情况符合某种函数分布的规律，较为常用的参数分布R-R、S-N、L-N的分布函数分别为
 | (5) |
 | (6) |
 | (7) |
式中：D为颗粒特征尺寸；σ为特征参数。
2 反问题模型 2.1 烟花算法 烟花算法的基本原理是通过模拟烟花的2种爆炸行为：①火花较多，爆炸范围较小；②火花较少，爆炸半径较大，来实现在问题的解空间内搜索反问题的最优解。在烟花算法中，适应度较小的烟花(距离最优解更接近)会采用第1种爆炸行为，以更好地搜索邻近区域。适应度较大的烟花(距离最优解较远)会采用第2种爆炸行为，以更好地搜索较远距离的空间。
对于一个k维问题来说，火花数量以及爆炸半径定义为[13]
 | (8) |
 | (9) |
式中：ni和Ai分别为第i个烟花产生的火花数和爆炸半径；Ntot和Amax分别为产生火花的总数和最大爆炸半径；m为烟花的数量；fmax和fmin分别为最大和最小的适应度函数值；f(Xi)为第i个烟花Xi的适应度函数；ξ为一个非零极小值, 用来防止除数为零。
通过式(8)和式(9)可以看出, 如果某一个烟花距离最优解很近的时候，爆炸半径会很小，然而产生的火花数量却很大，这是没有必要的。因此，需要更进一步地限制火花的产生数量和爆炸的半径范围，其公式分别为
 | (10) |
 | (11) |
 | (12) |
式中：nmin和nmax分别为允许的最小和最大火花数量；Ainit和Afinal分别为初始和最终的爆炸半径；t和tmax分别为迭代次数和最大迭代次数。
如何通过烟花爆炸生成火花以及如何选择下一代的烟花对于烟花算法来说是极为重要的步骤。对于一个k维问题，火花是随机由烟花在解空间内产生的，火花的位置可以由式(13)计算:
 | (13) |
式中：xij为由第i个烟花产生的新的火花的第j维位置；xij为第i个烟花产生的火花的第j维位置；r1为一个在[-1, 1]之间均匀分布的随机数。为了更好地适应多维问题，减少计算量。在一次爆炸过程中，并不是所有的方向都受影响。k维问题中，k0个方向随即被选择按照式(13)进化。如果新产生的火花超出了搜索区域，则按照式(14)进行更新：
 | (14) |
式中：xminj和xmaxj分别为第j维的最小和最大值；r2为一个均匀分布在(0, 1)之间的随机数。
为了保持烟花的多样性，防止算法陷入局部收敛，一部分烟花被选择来通过另外一种方式生成火花。
 | (15) |
式中：g为一个高斯分布随机数。
在生成了所有的火花之后，需要选择下一代的烟花位置。在烟花算法中，新一代的烟花位置是从上一代的烟花以及火花中选择的。上一代中的适应度函数最低的烟花或者火花肯定会被选择成为下一代的烟花位置。其他的烟花位置则从上一代的烟花和火花中根据它们距离最优位置的不同随机选取。被选中的概率以及距离计算公式分别为
 | (16) |
 | (17) |
式中：K为上一代所有烟花以及火花的位置集合；d(Xi, Xj)为颗粒i与其他颗粒j之间的距离。
2.2 改进的烟花算法 通过2.1节的介绍发现，烟花的爆炸半径是和其适应度值保持正相关的关系。在搜索的后期，这个性质对于找到更精确的解是很有帮助的。然而在初期，如果某代最优烟花距离最优解也较远时，这将会极大地降低算法的收敛速度(如图 2所示)。除此之外，下一代烟花的选择方式也会降低其进化速率。因此，本文提出了一种改进的烟花算法，每一个烟花被添加一个速度，用来提高其搜索效率(如图 2所示)。下一代烟花的选择方式为
 |
| 图 2 烟花算法的搜索过程示意图 Fig. 2 Schematic of searching process for firework algorithm |
| 图选项 |
 | (18) |
 | (19) |
式中：Vi(t)为第t代第i个烟花的搜索速度；Xib(t)为第i个烟花产生的适应度值最低的火花；r3为一个均匀分布在[0, 1]之间的随机数；Pg为每一代中最好的位置；w为一个权重系数，决定了上一代的速度对下一代的影响程度。在本文中，w=0.9-0.5t/tmax，改进烟花算法与标准的烟花算法的计算结果见文献[17]，本文不再赘述。
3 结果与讨论 本文提出一种利用瞬态激光入射参与性介质颗粒系，然后通过测量特定方向的时间分辨的透射率同时反演参与性介质颗粒系的粒径分布和光学常数的方法。此方法需要在2块厚度不同的一维吸收、散射、非发射的半透明参与性介质样品左侧，分别入射2个不同波段的连续激光。
介质的光学常数通常是与波长相关的。由于在本文中，涉及2种不同波长的入射激光，因此需要反演的光学常数为2组。也就是说，需要同时反演6个参数(2组光学常数以及1组粒径分布参数)。综合考虑到所需测量的信息个数以及实验实现过程的复杂程度，使用双波长双厚度法对颗粒系的光学常数及粒径分布进行重建。并使用上述条件下2个不同角度的透射信号作为已知信息。在本文以下的算例中，透射信号的测量角度θ设置为18.5°以及56°。探测器接收角为7.5°。样品厚度分别为0.02和0.04 m。激光的脉冲宽度为0.2 m。实验的基本操作流程如下：首先测量波长为λ1的激光作用条件下，厚度为L1的颗粒系的2种角度下的透射信号；然后更换厚度为L2的样品测量相同条件下的透射信号；最后更换波长为λ2的激光器，重复上述操作。
本文对于粒径反演以及光学常数的反演是通过使2个不同角度的透射信号的测量值与模拟值匹配来完成的。也就是通过改进的烟花算法得到目标函数的最小值，目标函数如下：
 | (20) |
式中：下角标a、b和c分别代表不同厚度、不同入射波长和不同透射角度的辐射信号；t0为测量时间；Test和Tmea分别为计算和实际测量得到的透射信号。反演流程如图 3所示。首先使用实际的物性通过正问题模型计算出不同角度的时域透射信号，然后与假设的物性计算得出的结果进行对比，如果式(20)的值小于某一预先设定的ε=10-6的话，则假设的物性参数即为所需测量的物性，反演结束。否则，使用本文提出的改进的烟花算法重新假设新的物性参数。重复以上流程直到目标函数值Fobj < ε为止。
 |
| 图 3 物性反演流程图 Fig. 3 Flowchart of parameter retrieval |
| 图选项 |
烟花算法的参数设置如表 1所示，ng为高斯烟花的个数，Vmax为烟花的最大搜索速度。其粒径分布及光学常数设置如表 2所示。其中假设所需测量的颗粒系为微藻悬浮液，因此光学常数的值设置为微藻的实际测量值[18]，且光学常数与粒径分布的搜索范围均设置为实际微藻颗粒的取值范围，并假设微藻颗粒近似为球形。粒子数密度设置为定值1012 m-3。对光学常数实部n值、虚部k值、颗粒特征尺寸D和特征参数σ进行同时重建，规定其搜索范围为：n∈(1.3, 1.5)、k∈(10-3, 10-2)、D∈(0, 10)、σ∈(0, 10)。测量时间t0设置为0.2 s。
表 1 改进的烟花算法参数设置 Table 1 Parameter setting of improved firework algorithm
| 参数 | Amax | Ainit | Afinal | Ntot | ng | nmax | nmin | tmax | Vmax | w |
| 数值 | 5.0 | 0.1 | 0.001 | 100 | 5 | 20 | 5 | 3 000 | 10.0 | 0.9-0.5t/tmax |
表选项
表 2 粒径分布以及复折射率真值 Table 2 True value of particle size distribution and complex refractive index
| 参数分布 | (D, σ) | (n, k) | |
| 波长为λ1 | 波长为λ2 | ||
| R-R | (3.4, 1.8) | (1.359, 4.9×10-3) | (1.363, 2.5×10-3) |
| S-N | (4.1, 2.3) | (1.352, 4.9×10-3) | (1.358, 2.5×10-3) |
| L-N | (3.6, 4.2) | (1.352, 5.0×10-3) | (1.353, 2.5×10-3) |
表选项
反演结果如表 3所示，可以看出，对于参与性介质颗粒系光学常数的反演效果较好。最大误差小于1%。尽管对于粒径分布参数的反演偏差相对于光学常数来说较大(最大2.22%)。然而，实际的粒径分布曲线吻合较好，如图 4所示。
表 3 光学常数和粒径分布反演结果 Table 3 Retrieval results of optical constant and particle size distribution
| 参数分布 | (n, k) | (εrel(n), εrel(k))/% | (D, σ) | (εrel(), εrel(σ))/% | |||
| 波长为λ1 | 波长为λ2 | 波长为λ1 | 波长为λ2 | ||||
| R-R | (1.358 9, 4.85 × 10-3) | (1.364 5, 2.52 × 10-3) | (0.01, 1.04) | (0.11, 0.67) | (3.45, 1.84) | (1.47, 2.22) | |
| S-N | (1.352 5, 4.94 × 10-3) | (1.356 0, 2.52 × 10-3) | (0.03, 0.80) | (0.11, 0.95) | (4.17, 2.33) | (1.71, 1.30) | |
| L-N | (1.351 6, 4.99 × 10-3) | (1.352 7, 2.51 ×10-3) | (0.03, 0.13) | (0.02, 0.50) | (3.57, 4.14) | (0.83, 1.43) | |
表选项
 |
| 图 4 粒径分布反演结果 Fig. 4 Retrieval results of particle size distribution |
| 图选项 |
4 结论 本文采用一种改进的烟花算法同时反演了粒径分布与光学常数，得到：
1) 使用一种波长的激光入射待测颗粒系，无法准确地通过测量透射信号来同时得到其粒径分布与光学常数。
2) 通过对3种常见的粒径分布进行测试发现，本文提出的方法可以准确得到颗粒系粒径分布与光学常数。
参考文献
| [1] | 王亚辉, 王强, 张伯川, 等. 高超声速飞行器红外窗口热辐射特性试验[J].红外与激光工程, 2015, 44(6): 1716–1720. WANG Y H, WANG Q, ZHANG B C, et al. Experiment of the thermo-radiation characteristic of infrared window of hypersonic vehicles[J].Infrared and Laser Engineering, 2015, 44(6): 1716–1720.(in Chinese) |
| [2] | FRANKEL J I. Cumulative variable formulation for transient conductive and radiative transport in participating media[J].Journal of Thermophysics & Heat Transfer, 2015, 9(2): 210–218. |
| [3] | QIAO Y, HONG Q, QIN C, et al. Multi-start iterative reconstruction of the radiative parameter distributions in participating media based on the transient radiative transfer equation[J].Optics Communications, 2015, 351: 75–84.DOI:10.1016/j.optcom.2015.04.048 |
| [4] | COQUARD R, ROCHAIS D, BAILLIS D. Conductive and radiative heat transfer in ceramic and metal foams at fire temperatures[J].Fire Technology, 2012, 48(3): 699–732.DOI:10.1007/s10694-010-0167-8 |
| [5] | 汪雪, 苏明旭, 蔡小舒. 超声衰减谱法颗粒粒径测量中遗传算法参数优化[J].上海理工大学学报, 2016, 38(2): 148–153. WANG X, SU M X, CAI X S. Parameters optimization of genetic inversion algorithm for particle sizing by ultrasonic attenustion spectroscopy[J].Journal of University of Shanghai for Science and Technology, 2016, 38(2): 148–153.(in Chinese) |
| [6] | DREYER A, KIRCHGEORG T, WEINBERG I, et al. Particle-size distribution of airborne poly-and perfluorinated alkyl substances[J].Chemosphere, 2015, 129: 142–149.DOI:10.1016/j.chemosphere.2014.06.069 |
| [7] | REN Y, QI H, CHEN Q, et al. Simultaneous retrieval of the complex refractive index and particle size distribution[J].Optics Express, 2015, 23(15): 19328–19337.DOI:10.1364/OE.23.019328 |
| [8] | LUK'YANCHUK B S, VOSHCHINNIKOV N V, PANIAGUA-DOMíNGUEZ R, et al. Optimum forward light scattering by spherical and spheroidal dielectric nanoparticles with high refractive index[J].ACS Photonics, 2015, 2(7): 993–999.DOI:10.1021/acsphotonics.5b00261 |
| [9] | ZHANG B, XU C L, WANG S M. An inverse method for flue gas shielded metal surface temperature measurement based on infrared radiation[J].Measurement Science and Technology, 2016, 27(7): 074002.DOI:10.1088/0957-0233/27/7/074002 |
| [10] | DEEPA O, SENTHILKUMAR A. Swarm intelligence from natural to artificial systems:Ant colony optimization[J].Applied Optics, 2016, 8(1): 9–17. |
| [11] | EATON J, YANG S, MAVROVOUNIOTIS M. Ant colony optimization with immigrants schemes for the dynamic railway junction rescheduling problem with multiple delays[J].Soft Computing, 2016, 20(8): 1–16. |
| [12] | QI H, REN Y T, CHEN Q, et al. Fast method of retrieving the asymmetry factor and scattering albedo from the maximum time-resolved reflectance of participating media[J].Applied Optics, 2015, 54(16): 5234–5242.DOI:10.1364/AO.54.005234 |
| [13] | TAN Y, ZHU Y.Fireworks algorithm for optimization[C]//International Conference on Advances in Swarm Intelligence.Heidelberg:Springer, 2010, 6145:355-364. |
| [14] | IMRAN A M, KOWSALYA M. A new power system reconfiguration scheme for power loss minimization and voltage profile enhancement using fireworks algorithm[J].International Journal of Electrical Power & Energy Systems, 2014, 62: 312–322. |
| [15] | MODEST M F. Radiative heat transfer[M].New York: Academic Press, 2013: 265-270. |
| [16] | HODKINSON J, GREENLEAVES I. Computations of light-scattering and extinction by spheres according to diffraction and geometrical optics, and some comparisons with the Mie theory[J].Journal of the Optical Society of America, 1963, 53(5): 577–588.DOI:10.1364/JOSA.53.000577 |
| [17] | REN Y T, QI H, HE M J, et al. Application of an improved firework algorithm for simultaneous estimation of temperature-dependent thermal and optical properties of molten salt[J].International Communications in Heat & Mass Transfer, 2016, 77: 33–42. |
| [18] | PILON L, BERBEROG.LU H, KANDILIAN R. Radiation transfer in photobiological carbon dioxide fixation and fuel production by microalgae[J].Journal of Quantitative Spectroscopy & Radiative Transfer, 2011, 112(17): 2639–2660. |
