上海理工大学学报  2020, Vol. 42 Issue (4): 332-338   PDF    
基于改进差分进化算法的超声衰减谱反演计算
蒋瑜, 贾楠, 苏明旭     
上海理工大学 能源与动力工程学院,上海 200093
摘要: 针对传统差分进化算法存在早熟收敛和求解精度低的缺点,研究了一种自适应控制参数的差分进化算法。通过引入自适应控制变量因子、自适应缩放因子和交叉因子使种群不断地向更新成功的个体学习,促进了后续种群的进化。对于颗粒粒径分布服从高斯分布、R-R(Rosin-Rammler)分布以及对数正态分布的3种典型颗粒系进行数值模拟,研究算例发现,改进差分进化算法反演得出分布参数值 $\bar R$ K的误差小于5%,体积中位径相比于设定分布的误差小于5%,因此,改进差分进化算法具有较强的稳定性与抗噪性。
关键词: 衰减谱     差分进化算法     反演     粒度分布    
Inverse calculation of ultrasonic attenuation spectrum based on an improved differential evolution algorithm
JIANG Yu, JIA Nan, SU Mingxu     
School of Energy and Power Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China
Abstract: The traditional differential evolution algorithm frequently encounters the problem of premature convergence and low accuracy. By introducing adaptive control variable factors, a differential evolution (DE) algorithm was developed to ensure that the population could be updated continuously and successfully. The individual learning, which is conducive to the evolution of subsequent populations, and the overall inversion yield more accurate results. As a validation, the particle systems obeying three typical distribution functions of Gaussian distribution, R-R distribution and lognormal distribution were numerically simulated. The resultant distribution parameter values $\bar R $ and K yield errors less than 5%, and the deviations of median volume diameter are within ±5% compared with the given distribution. The DE algorithm also reveals obvious stability and noise resistance.
Key words: attenuation spectrum     differential evolution algorithm     inversion     particle size distribution    

近年来随着与颗粒紧密关联的行业如能源、化工、食品、环境、医药等的飞速发展,颗粒两相流测量技术获得了广泛的关注[1-5]。超声颗粒测量技术具有穿透力强、非浸入式、易实现在线检测、可测粒径范围宽等特点,其中,超声衰减谱法[6-7]是一种具有广阔发展前景的在线颗粒测量方法,超声波在颗粒两相体系中传播并与颗粒之间发生物理作用,得到超声衰减谱中包含的大量颗粒粒径信息,结合BLBL(Bouguer-Lambert-Beer-Law)、ECAH(Epstein-Carhart-Allegra-Hawley)模型[8-9],进一步反演计算,可以快速分析颗粒粒径分布。

颗粒粒径测量的反演问题得到了广泛关注。栾志超[10]利用拟牛顿算法模拟反演粒径分布。苏明旭等[11]曾通过非线性最小二乘优化LM算法对目标函数进行反演计算。孔明等[12]用自动基线拟合的累计量算法对光子相关光谱函数进行反演。汪雪等[13]采用改进遗传算法(genetic algorithm)反演颗粒分布参数。与上述方法相比,差分进化算法(differential evolution, DE)[14]具有很好的全局搜索能力和鲁棒性。与标准遗传算法相比,差分进化算法所需控制参数较少,原理简单易行,已在许多优化问题中获得了良好效果[15-17]

针对超声法颗粒测量,本文研究了一种非独立模式条件下的超声衰减谱反演差分进化改进算法,针对标准差分进化算法容易陷入局部最优以及早熟收敛问题,提出了自适应控制参数因子并进行优化,通过对高斯分布、R-R(Rosin-Rammler)分布以及对数正态分布函数的反演过程模拟,将改进差分进化算法和其他算法比较,分析了改进差分进化算法的特点,尤其是不同设定峰值数情况下噪声水平影响及抑制。

1 超声衰减谱法及反演优化问题

当超声波通过含颗粒物介质时,由于颗粒与声波的作用而产生衰减,声波在两相之间的分界面处,会在固体颗粒的内部和外部产生压缩波,同时在内部产生热波、剪切波。在球坐标下,运用边界条件得到在颗粒内部和周围介质的3种波动的势函数,其线性方程组可简写为[18]

$ {{{M}}}{[{A_{{n}}},{C_{{n}}},A_{{n}}^{'},C_{{n}}^{'},{B_{{n}}},B_{{n}}^{'}]^{\rm{T}}} = {{e}} $ (1)

式中:M为6阶方阵;AnBnCn分别表示压缩波、热波以及剪切波;e为散射系数向量。

求解方程组(1)可得到系数An,声波总衰减由系数An决定,最终的衰减系数

${\alpha _{\rm{s}}} = - \frac{{3\varphi }}{{2k_{{c}}^2{R^3}}}\sum\limits_{{{n}} = 0}^\infty {(2n + 1){A_n}} $ (2)

式中:kc为连续介质中的波数;φ为颗粒相体积浓度;R为颗粒半径。

$k = \frac{\omega }{{{c_{\rm{s}}}(\omega )}} + {\rm{j}}{\alpha _{\rm{s}}}(\omega )$ (3)
${\left( {\frac{k}{{{k_{\rm{c}}}}}} \right)^2} = 1 + \frac{{3\varphi }}{{{\rm{j}}k_{\rm{c}}^2{R^3}}}\sum\limits_{{{n = }}0}^\infty {(2n + 1){A_n}} $ (4)

式中: $ {\rm{j = }}\sqrt {{\rm{ - 1}}} $ cs为混合介质的声速;k为复波数;ω为角频率。

${{\alpha}} = {{S}} \times {{W}}$ (5)

式中:α为声衰减分布列向量;S为衰减系数矩阵;W为被测颗粒的粒径分布。

为了便于采用最优化理论求解,往往根据超声频率、粒径分布以及模型参数结合ECAH模型计算理论声衰减谱,同时将其同实验声衰减谱比较,定义误差函数

${E_{\rm{SSD}}} = \sum\limits_{{{j}} = 1}^{{N}} {({\alpha _{\rm{m}}} - {\alpha _{\rm{s}}}} {)^2}$ (6)

式中:αs为理论声衰减预测谱;αm为实验声衰减谱。

通过反演算法最小化误差函数,即可求解得最优的颗粒粒径分布。对于最优化求解,通常颗粒系被描述为粒径频度分布服从一定的函数形式,例如,高斯分布、对数正态分布或R-R(Rosin-Rammler)分布。其中,R-R分布函数为

$f(R) = \frac{k}{R}{(R/\bar R)^{k - 1}}\exp \left( { - {{\left( {\frac{R}{{\bar R}}} \right)}^k}} \right)$ (7)

式中: $\overline {{R}} $ 为特征尺寸;K为分布特征参数

对不同的优化算法适用性进行研究,设计算例,对半径为5 μm的单峰R-R分布玻璃微珠颗粒和设定参数双峰分布系( $\overline {{R_1}} $ =5,K1=10, $\overline {{R_2}} $ =15,K2=25,W1=0.5)进行数值反演。W1为第一个峰所占的比例。按式(6)形式构造目标函数,分别采用DFP(Davidon-Fletcher-Powell),LM(Levenberg-Marquard)、BFGS(Broyden-Fletcher-Goldfarb-Shanno)和基本差分进化算法(DE)在0.1~50 μm粒径区间寻优。从图1中可以看出,无论是给定单峰还是双峰分布,3种局部最优化算法DFP,LM,BFGS的结果均明显偏离真值,而DE算法的结果总体上符合对真值范围的预期,双峰反演结果存有一定偏差,但明显优于DFP,LM,BFGS算法,由此可见,基于DE原理的反演算法完全可行。VD)为体积频度分布。


图 1 4种算法反演计算的颗粒粒径分布 Fig. 1 Inversion results by four algorithms for the given particle size distribution
2 差分进化算法

差分进化算法是一种基于群体智能理论的优化算法,通过群体内个体间的合作与竞争来对目标进行优化。近年来差分进化算法在最优化问题求解方面得到了广泛的应用。

2.1 标准差分进化算法

基本的差分进化算法主要由4个步骤组成。

a. 初始化种群。初始种群 ${\rm{\{ }}{x_i}{\rm{(0)}}|x_{j,i}^{\rm L}$ ${x_{j,i}}{\rm{(0)}}$ $x_{j,i}^{\rm U},i = 1,2, \cdots ,N,j = 1,2, \cdots ,D{\rm{\} }}$ 随机产生

${x_{j,i}}{\rm{(0) = }}x_{j,i}^{\rm L} + {\rm{rand}}(0,1)(x_{j,i}^{\rm U} - x_{j,i}^{\rm L})$ (8)

式中:N为种群大小; $ {x}_{j,i}^{{\rm{L}}} $ $ {x}_{j,i}^{{\rm{U}}} $ 分别为xj取值的下限和上限,rand(0,1)为在(0,1)区间均匀分布的随机数。

b. 变异操作。通过如下差分策略,产生变异中间体

${v_i}(g + 1) = {x_{{r_1}}}(g) + F({x_{{r_2}}}(g) - {x_{{r_3}}}(g))$ (9)
$i \ne r_1 \ne r_2 \ne r_3$

式中:F为缩放因子,本文取F=0.35。

c. 交叉操作。对第g代种群{xig)}及其变异的中间体{vig+1)}进行个体间的交叉操作

${u_{j,i}}(g + 1)\! =\! \left\{ \begin{aligned} &{v_{j,i}}(g + 1),\\ &\;\;\;\;{\text{若}}{\rm{rand}(0,1)} \leqslant CR{\text{或者}} j = {j_{{\rm{rand}}}}\\ &{x_{j,i}},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\text{其他}} \end{aligned} \right.$ (10)

式中:CR为交叉因子,本文取CR=0.85;jrand为随机取出解向量的位置。

d. 选择操作。差分进化算法采用贪婪算法来选择进入下一代种群的个体。

${x_i}(g + 1)\! = \!\left\{ \begin{array}{l} \!\!\!{u_i}(g + 1),\;{\text{若}}f({u_i}(g + 1)) \leqslant f({x_i}{\rm{(g)}})\\ \!\!\!{x_i}{\rm{(g)}},\;\;\;\;\;\;\;\,{\text{其他}} \end{array} \right.$ (11)

式中,F为目标函数值。

反复执行式(9)~(12),直到满足终止条件为止。

2.2 算法改进

DE算法主要有3个控制参数,即种群大小N、缩放因子F和交叉因子CR。缩放因子F控制偏差变量的放大作用,改变搜索的方向,若种群过早收敛,则F应增加,反之亦然。交叉因子CR改变种群多样性,较大的CR会加速收敛。标准DE算法中的FCR为固定值,在算法迭代后期会因其种群个体聚集,造成种群间差异性变小,使算法易陷入局部最优解、出现早熟收敛等问题。因此,在算法中引入自适应控制参数因子加以改进,改进后的缩放因子和交叉因子为

$\begin{aligned}\!\!\!\!\!{F_{i,\;G + 1}} \!=\! \left\{ \begin{array}{l} \!\!\!{F_{\rm{L}}} +{\rm{rand}(0,1)}{F_{\rm{U}}},\;{\text{若}}\;{\rm{rand}(0,1)} \leqslant {\tau _1}\\ \!\!\!{F_{i,\;G}},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\qquad{\text{其他}} \end{array} \right.\end{aligned}$ (12)
$C{R_{i,\;G + 1}} = \left\{ \begin{array}{l} {\rm{rand}(0,1)}\;\;\;,\;\;\;{\text{若}}\;{\rm{rand}(0,1)} < {\tau _1}\\ C{R_{i,\;G}}\;\;\;,\;\;\;{\text{其他}} \end{array} \right.$ (13)

式中:τ1τ2为调节因子;FLFU分别为缩放因子的下限和上限。

图2给出了不同的FLFU对应的目标函数值,为了提高优化效率,给定FL(0.1,0.4),FU(0.5,0.9)的范围,在多种粒径分布参数下进行反演寻优,选取目标函数值最小处对应的FLFU图2(b)FLFU对应的目标函数值等高线图,由图2(b)可确定FLFU的可选范围,对应目标函数值最小,将FLFU分别选取为0.3,0.8。


图 2 适应值随FUFL的变化 Fig. 2 Variation of the best fitness value along with FU and FL

将改进后的算法称为IDE(improved differential evolution)算法,算法流程如图3所示。


图 3 IDE算法流程 Fig. 3 Flowchart of the improved differential evolution algorithm

a. 确定DE算法参数,随机产生初始种群并计算个体的适应度(fitness)。

b. 判断是否最佳适应度基本不变。若是,则输出最佳个体为最优解;若否,结合自适应控制参数产生缩放因子F、交叉因子CR,进行变异和交叉操作,得到中间种群,在原种群和中间种群选择个体,得到新种群。

c. 若满足停止条件,算法结束;反之,进化代数t=t+1,转步骤b。

3 结果与分析 3.1 数值模拟

为了验证改进算法的性能,分别采用标准DE算法和改进算法IDE(improved differential evolution)对满足单峰、双峰分布的颗粒群进行数值模拟。

3.1.1 单峰分布

按照表1设定参数,对服从R-R、正态和对数正态分布的颗粒系进行反演计算。 $\overline {{R}} $ K(或σ)分别为特征尺寸以及分布特征参数。可以看出,标准DE算法和改进算法反演得到的结果均与原始分布相同,改进算法IDE计算一次需要时间10 s左右,比DE算法多2 s左右。


表 1 3种分布函数下的单峰参数反演结果 Table 1 Inversion results of three kinds of distribution functions by DE and IDE
3.1.2 双峰分布

双峰分布反演结果如表2所示。 $\overline {{R_1}} $ K1 $\overline {{R_2}} $ K2W1分别为第一个峰的特征尺寸、分布特征参数,第二个峰的特征尺寸、分布特征参数,及第一个峰所占的比例。可以看出,当设定分布变化时,IDE均表现出一定的鲁棒性,而DE算法虽然多数峰值位置偏离不大,但是,其分布宽度均存在较大的反演误差,例如, $\overline {{R_1}} $ =5,K1=10, $\overline {{R_2}} $ =15,K2=25,W1=0.5时,DE算法反演出K1=6.91,K2=19.87,结果偏差较大。


表 2 R-R分布函数下的5个参数反演值 Table 2 Inversion results of R-R distribution functions by DE and IDE

表2数据还可以看出,在上述算例将频谱范围从之前1~10 MHz增至1~15 MHz时,反演结果和设定值吻合。增加频宽的范围对于反演更有利,但频宽增加对于实验技术提出了更高要求。

根据反演结果求得的颗粒特征尺寸如表3所示,Dv10Dv50Dv90分别表示颗粒体积累积概率达到10%,50%,90%时对应的颗粒粒径,可以看出,改进差分进化算法反演得到的特征粒径更接近设定分布值,其相比于设定分布的误差均小于5%,而原始反演算法得到的误差甚至可以达到30%。


表 3 R-R分布函数下反演得到的颗粒特征尺寸 Table 3 Inversion results of R-R distribution functions

对双峰分布的反演结果如图4所示,可以看出,IDE反演结果均接近原始分布,而DE算法反演分布曲线相比原始分布均出现一定偏离。如图4(a)所示,其第一个峰的位置向右产生了偏移,同时两个峰的峰值均小于原始设定值。同步获得目标函数反演残差,原始DE算法谱误差ESSD=0.011,IDE算法反演得到的ESSD=1.09×10−14图4(b)中第二个峰的结果与设定基本吻合,但第一个峰的反演粒径偏小且峰宽偏窄,对应原始DE算法反演得到的ESSD=0.0027,IDE算法反演得到的ESSD=1.04×10−12,可见IDE求解精度得到明显提高。


图 4 双峰颗粒粒径分布反演结果 Fig. 4 Inversion results of bimodal particle size distribution
3.2 抗噪性分析

在超声衰减谱的实际测量过程中会不可避免地引入测量误差(如实验过程噪声或谱分析误差),为了模拟实际测量中的随机噪声和其他干扰,在超声衰减谱中人为加入3%,5%,10%的随机噪声,进行反演计算。为了讨论算法的随机性影响和抑制,分析了25组连续反演结果,并将其均值作为最终结果与设定颗粒粒径分布和无噪声情形进行对比。

图5(a)所示,设定R-R分布 $\overline {{R}} $ =5,K=7,加入不同幅度随机噪声时,IDE反演R-R分布特征参数K出现了一定波动,在加入3%的噪声时,K值呈现出微小震荡,平均值为6.91,相对误差为1.3%,95%置信水平下的置信区间为(6.58~7.23)。而当加入5%,10%的噪声时,会出现少量严重偏离真值的情况,但多数结果都分布在真值附近。如图5(b)所示,特征尺寸 $\overline {{R}} $ 则基本不会随噪声幅度而产生较大变化,3%,5%,10%噪声下的平均值分别为5.02,5.05,5.21,相对误差均在5%以内,置信区间分别为(4.98~5.05),(4.98~5.12),(5.1~5.34),即反演结果基本稳定在噪声水平以内,没有放大。


图 5 在不同噪声幅度下K $\overline {{R}} $ 值的变化 Fig. 5 Changes of the values of $\overline {{R}} $ and K with different noise magnitudes

不同噪声幅度下单峰R-R分布反演结果如图6(a)所示,可以看出,采用改进差分进化算法的反演曲线光滑平顺,未加入噪声时反演分布和设定曲线重合;超声谱加入3%,5%,10%的噪声后,反演曲线稍有偏移,但其 $\overline {{R}} $ 值偏差始终在5%以内,且基本保持了原有分布宽度。R-R双峰分布颗粒系的超声谱加入1%,2%噪声后的反演结果如图6(b),其反演出两个峰的分布宽度相对偏窄,但偏差均在5%及以内;当加入3%噪声后,偏差变大,第一个峰分布过宽,第二个峰也向左偏移且分布宽度变窄,其偏差在15%及以内。考虑到实验室条件下超声谱的测量随机误差为1%~2%,上述结果表明,IDE对于单峰和双峰分布颗粒系的反演均是有效的,对反演结果求平均值(本文采用了25次),可以较好地抑制此类噪声的影响,相对误差小于5%。


图 6 R-R分布函数加入噪声的反演计算结果 Fig. 6 Inversion results of R-R distribution functions with noise
4 结 论

超声衰减谱法粒径反演时,通过对差分进化算法优化,自适应控制缩放因子F、交叉因子CR,由本文数值模拟结果可知,改进差分进化算法反演得出的分布参数值 $\overline {{R}} $ K误差小于5%,而标准差分进化算法反演得到的 $\overline {{R}} $ K值误差达到20%及以上。改进差分进化算法得到的Dv10Dv50Dv90与设定分布的相对误差均小于5%,而原始差分进化算法得到的误差甚至可以达到30%。反演结果显示,改进差分进化算法具有很强的抗干扰能力,反演峰值位置偏差很小,粒度相对误差不超过5%。

参考文献
[1]
蔡小舒, 苏明旭, 沈建琪. 颗粒粒度测量技术及应用[M]. 北京: 化学工业出版社, 2010.
[2]
刘东红, 朱潘炜. 超声检测技术及在食品安全检测中的研究进展[J]. 食品工业科技, 2009, 30(9): 373-376.
[3]
ABDA F, AZBAID A, ENSMINGER D, et al. Ultrasonic device for real-time sewage velocity and suspended particles concentration measurements[J]. Water Science & Technology, 2009, 60(1): 117-125.
[4]
严玉鹏, 唐亚东, 万彪, 等. 颗粒尺寸对纳米氧化物环境行为的影响[J]. 环境科学, 2018, 39(6): 2982-2990.
[5]
付加, 祁贵生, 刘有智, 等. 超重力湿法脱除气体中细颗粒物研究[J]. 化学工程, 2015, 43(4): 6-10. DOI:10.3969/j.issn.1005-9954.2015.04.002
[6]
MOUGIN P, WILKINSON D, ROBERTS K J, et al. Characterization of particle size and its distribution during the crystallization of organic fine chemical products as measured in situ using ultrasonic attenuation spectroscopy[J]. The Journal of the Acoustical Society of America, 2001, 109(1): 274-282. DOI:10.1121/1.1331113
[7]
POVEY M J W. Ultrasound particle sizing: a review[J]. Particuology, 2013, 11(2): 135-147. DOI:10.1016/j.partic.2012.05.010
[8]
EPSTEIN P S, CARHART R R. The absorption of sound in suspensions and emulsions. I. Water fog in air[J]. The Journal of the Acoustical Society of America, 1953, 25(3): 553-565. DOI:10.1121/1.1907107
[9]
林春丹, 梁永燊, 张万松, 等. 超声衰减谱法测量含蜡原油中蜡晶粒度[J]. 声学技术, 2013, 32(4): 294-298.
[10]
栾志超. 激光粒度仪的无模式数据处理算法研究[D]. 天津: 天津大学, 2005.
[11]
苏明旭, 王夕华, 黄有贵, 等. LM优化算法在消光法测粒技术中的应用[J]. 仪器仪表学报, 2005, 26(S1): 63-65.
[12]
孔明, 赵军, 李春燕. 自动基线拟合累积量算法纳米颗粒反演[J]. 光电工程, 2009, 36(9): 52-55. DOI:10.3969/j.issn.1003-501X.2009.09.010
[13]
汪雪, 苏明旭, 蔡小舒. 超声衰减谱法颗粒粒径测量中遗传算法参数优化[J]. 上海理工大学学报, 2016, 38(2): 148-153.
[14]
STORN R, PRICE K. Differential evolution—a simple and efficient heuristic for global optimization over continuous spaces[J]. Journal of Global Optimization, 1997, 11(4): 341-359. DOI:10.1023/A:1008202821328
[15]
崔利刚, 邓洁, 王林, 等. 基于改进联合采购及配送模型的RFID投资决策研究[J]. 中国管理科学, 2018, 26(5): 86-97.
[16]
周华, 刘昱, 郭谨玮, 等. 基于差分进化算法和相关向量机的车辆油耗预测[J]. 汽车技术, 2018(12): 19-22.
[17]
简兆圣, 艾剑良. 差分进化算法在气动力参数辨识中的应用[J]. 复旦学报(自然科学版), 2017, 56(5): 545-550.
[18]
GREENWOOD M S, PANETTA P D, BAMBERGER J A, et al. System and technique for ultrasonic characterization of settling suspensions: USA, US7140239B2[P]. 2006-11-28.