基于改进差分进化算法的超声衰减谱反演计算

蒋 瑜， 贾 楠， 苏明旭

（上海理工大学 能源与动力工程学院，上海 200093）

近年来随着与颗粒紧密关联的行业如能源、化工、食品、环境、医药等的飞速发展，颗粒两相流测量技术获得了广泛的关注[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为6 阶方阵；An，Bn，Cn分别表示压缩波、热波以及剪切波；e为散射系数向量。

求解方程组（1）可得到系数An，声波总衰减由系数An决定，最终的衰减系数

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

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

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

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

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

对不同的优化算法适用性进行研究，设计算例，对半径为5 μm 的单峰R-R 分布玻璃微珠颗粒和设定参数双峰分布系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 原理的反演算法完全可行。V（D）为体积频度分布。

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

2 差分进化算法

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

2.1 标准差分进化算法

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

式中：N为种群大小；分别为xj取值的下限和上限，rand（0,1）为在（0,1）区间均匀分布的随机数。

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

式中：F为缩放因子，本文取F=0.35。

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

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

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

式中，F为目标函数值。

反复执行式（9）～（12），直到满足终止条件为止。

2.2 算法改进

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

式中：τ1，τ2为调节因子；FL，FU分别为缩放因子的下限和上限。

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

将改进后的算法称为IDE（improved differential evolution）算法，算法流程如图3 所示。

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

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

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、正态和对数正态分布的颗粒系进行反演计算。R，K（或σ）分别为特征尺寸以及分布特征参数。可以看出，标准DE 算法和改进算法反演得到的结果均与原始分布相同，改进算法IDE 计算一次需要时间10 s 左右，比DE 算法多2 s 左右。

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

3.1.2 双峰分布

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

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

对双峰分布的反演结果如图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 求解精度得到明显提高。

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

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

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

3.2 抗噪性分析

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

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

图 5 在不同噪声幅度下K，值的变化Fig.5 Changes of the values of and K with different noise magnitudes

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

4 结 论

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