基于改进萤火虫算法的波阻抗反演方法

滕明，王典,2，于涵，张敬权，徐雨歆

1.吉林大学 地球探测科学与技术学院，长春 130026；2.吉林大学 地质资源立体探测虚拟仿真实验教学中心，长春 130026

0 引言

萤火虫算法(firefly algorithm,简称FA)是由Yang[1]于2008年提出的一种仿生随机优化算法。该算法设计思路源于对萤火虫觅食、求偶等行为模拟，通过萤火虫个体之间相互吸引、移动从而到达寻优目的。与其他群智能算法相比，萤火虫算法的优点表现为：算法结构简单、所需调整参数少以及流程清晰易懂，已在数值优化、工程技术和聚类分析等[2--6]方面得到成功应用。和多数仿生随机算法类似，FA算法也存在一些不足，例如收敛“早熟”现象、后期收敛速度缓慢和易陷入局部极值。针对FA算法存在的不足，许多学者进行了研究与改进。主要分为以下两大类：一类是改进FA算法参数设置方式。如Lukasik et al.[7]对吸收系数和移动步长进行修改，改进后FA算法的求解精度得以提高，但求解速度较慢；欧阳喆等[8]采用动态方式调整步长；Gandomi et al.[9]采用混沌序列修改FA算法参数设置，提升FA算法的求解精度和收敛速度；徐丽华等[10]针对FA算法“早熟”收敛现象，提出了一种变尺度混沌光强吸收系数调整策略的萤火虫优化算法；刘畅等[11]利用指数分布等改进吸引项。另一类是将FA与粒子群算法(PSO)、遗传算法(GA)和差分演化算法(DE)等优化算法相结合，提高FA算法的性能。曹秀爽[12]提出当FA算法处于全局搜索时，将模拟退火搜索机制融入FA算法，当FA算法处于局部搜索时，把回火策略引入FA算法，以此解决萤火虫易陷入局部极值和后期收敛速度慢等问题； Xia et al.[13]基于萤火虫算法和粒子群算法，提出一种混合优化算法。

笔者提出基于混沌搜索和动态步长修改的萤火虫算法(CAFA)。在种群初值选取时，利用混沌序列代替随机方式对萤火虫种群进行初始化。混沌是确定性非线性系统内存在的特殊随机现象，具有以下特性：有界性、随机性和遍历性。有界性是指混沌变量运动轨迹始终局限于一个确定区域内，此区域被称为混沌吸引域；随机性是指某一时刻的混沌变量会受到前一状态影响而确定出现；遍历性是指混沌变量运动轨迹不会停留于某一状态而是不重复地遍历空间区域中所有可能存在状态。混沌初始化提高了种群多样性和搜索空间分布质量，避免算法陷入局部最优。另外，采用一种新型的双曲线递减步长代替传统的固定步长，实现了提高算法后期收敛速度、寻优精度以及算法稳定性的目的。

波阻抗参数是联系地震、测井及地质信息的纽带。通过波阻抗反演可以把界面型的地震剖面转换成岩层型的波阻抗剖面，可直接将地震资料与测井对比进行储层岩性解释和物性分析，因此波阻抗反演技术在储层参数估算[14]、油藏预测[15]等研究工作中发挥关键的作用。地震资料波阻抗反演是一个非线性反演问题，传统的地震波阻抗反演常用线性化算法，如最速下降法、牛顿法和共轭梯度法等，此类方法计算效率高、稳定性好，但过于依赖初始模型且算法收敛速度慢、易陷入局部极值等缺点。近些年，遗传算法、粒子群算法和模拟退火算法等智能优化算法已经在波阻抗反演中得到广泛应用。 聂茹等[16]将免疫算法和遗传算法相结合应用于波阻抗反演；许永忠等[17]将粒子群算法应用到火山岩的波阻抗反演，将反演结果与概率神经网络方法相对比，证明了粒子群算法在波阻抗反演的有效性。萤火虫算法在波阻抗反演方向研究较少，笔者将CAFA算法应用于波阻抗反演，通过理论模型试算分析CAFA算法的可行性。

1 传统萤火虫算法(FA)

1.1 生物学机制

自然界中大多数萤火虫具有发光行为，发出的闪光信号可以被一定范围内其他萤火虫所感知，萤火虫借助此独特的发光行为进行觅食、求偶等群体行为。如图1所示，将萤火虫个体用搜索空间范围内红点表示，萤火虫的荧光亮度取决于搜索范围内红点所在位置的目标函数值。吸引度与亮度有关，越亮的萤火虫拥有越强的吸引力，吸引视线范围内亮度弱的萤火虫向其所在位置方向移动，若发光亮度相同，则萤火虫各自随机移动。把萤火虫个体之间的不断飞行和吸引的行为模拟成算法迭代和优化过程，舍弃萤火虫发光的某些生物学上意义，对闪光特性做出理想化处理[17]。

图1 萤火虫算法仿生模拟图Fig.1 Bionic simulation diagram of firefly algorithm

①假设萤火虫没有雌雄之分。萤火虫相互吸引与萤火虫的性别无关。

②萤火虫之间的吸引度与荧光亮度和距离有关。亮度高的萤火虫会吸引周围亮度弱的萤火虫，随着距离增大吸引度逐渐减小，亮度相同或是亮度比邻域其他萤火虫高的个体在自身周围随机运动。

③萤火虫的亮度与具体研究问题的目标函数有关。

1.2 FA算法描述

萤火虫算法的核心思想是荧光亮度弱的萤火虫会被亮度强的萤火虫吸引，并根据位置更新公式进行位置更新，经过多次迭代最终收敛到最优值附近。下面定义相对亮度、吸引力、位置更新等核心公式。

相对荧光亮度[17]:

(1)

式中：I0表示自身(r=0处) 的荧光亮度，与目标函数值有关，荧光亮度越高目标值越优；γ为光强吸收系数，衡量荧光亮度在传播路径中的减弱程度；rij为萤火虫i和j之间的欧氏距离，对于其他优化问题，距离公式依据具体问题进行定义。

吸引度[17]：

(2)

式中：β0表示光源处(r=0处)的吸引度；γ为光强吸收系数；rij为萤火虫i到萤火虫j的欧氏距离[17]：

(3)

位置更新公式[17]：

Xi(t+1)=Xi(t)+β(rij){Xj(t)-Xi(t)}

+α(rand-1/2)

(4)

Xi(t+1)=Xi(t)+α(rand-1/2)

(5)

式中：Xi,Xj是t时刻萤火虫i和j的位置；α为移动步长；β为吸引度系数；rand为0-1的随机数。把t时刻萤火虫Xi,Xj代入公式(1)并比较二者荧光亮度的大小，如果I(Xi)

在公式(4)中，Xi(t)表示萤火虫当前位置对更新位置的影响，起到平衡全局和局部搜索的能力；β(rij){Xj(t)-Xi(t))}反映群体信息交流的影响，是萤火虫之间的信息共享，取决于吸引力的大小；α(rand-1/2)是带有特定系数的随机扰动项，避免萤火虫陷入局部极值点。通过查阅文献以及测试大量复杂函数发现FA算法主要存在两点不足之处：①传统FA算法的位置更新公式采用固定步长向较优个体移动，如果被吸引萤火虫与较优个体之间距离小于固定步长因子，被吸引萤火虫在移动时会跃过较优个体，并随着迭代进行反复在较优个体附近震荡；如果步长较小，则需要多次的迭代计算，导致收敛速度降低。②FA采用随机方式对萤火虫个体位置进行初始化，这样会造成萤火虫位置分布不均匀，导致个体间差异较大，降低种群的多样性，算法易陷入局部最优点。

2 改进的萤火虫算法(CAFA)

2.1 动态步长

在FA 算法搜索过程中，如果固定步长较大会使搜索精度降低且收敛后期出现震荡现象。如果固定步长较小导致搜索速度降低，但求解精度提高。CAFA算法利用动态步长代替固定步长，动态步长是一种随迭代次数增加，步长不断变化的函数。在算法迭代初期，利用较大步长促使萤火虫个体加速向较优值附近移动，随着迭代次数增加步长会逐渐减小直至零值，避免了算法后期被吸引萤火虫个体与较优个体之间距离小于固定步长值而出现振荡现象。以四峰函数为例，当坐标为(0, 0)、 (0， -4)时，函数取得极大值fmax=2。图2经过15次迭代，对30只萤火虫进行了450次估算得到收敛效果图，图2a自适应步长萤火虫成功搜索到全局最优值fmax=2，然而图2b固定步长FA算法收敛精度为1.957，增加迭代次数固定步长FA算法收敛精度未能提高。

图2 动态步长(a)及固定步长(b) Fig.2 Dynamic step(a) and fixed step(b)

在理想情况下，群体中所有个体随着飞行次数的增加逐渐向最优值附近移动，最终收敛到一点。

因此，对于任意的两个萤火虫i和j可以得到[18]：

(6)

(7)

式(6)表示搜索空间内萤火虫个体收敛到一点(最优值)，式(7)表示收敛结束后个体(解)不再发生移动(变化)。结合上述公式可以得到[18]：

(8)

从式(8)可以得到，随着算法进行多次迭代，移动步长α趋于零。在传统FA中，设置步长α是静态的，它不能真正模拟算法搜索过程。通常对FA算法来说，α值较大有利于探索新的搜索空对全局最优收敛没有帮助。如果步长α值较小，则结果相反，因此移动步长α对算法的探索和收敛有着较大影响。本文使用双曲线方程设计动态步长调整方案。双曲线递减动态步长：

(9)

式中：t为当前迭代次数；Tmax为最大迭代次数；k为控制动态步长减小速率的参数。双曲线递减步长的动态递减性质如图3所示。从图3可以看出，步长在初期阶段较大，然后随着迭代次数的增加逐渐减小，这有助于算法平衡全局探索和局部搜索的能力。表1列举其他学者设计的各种动态步长[19--22]和笔者提出的双曲线递减步长。以上述四峰函数作为目标函数，测试各种动态步长的收敛效果。四峰函数全局最优值为fmax=2。这里定义差值参数DP，它是全局最优值与自适应步长萤火虫算法求解的差值，DP越接近于0值，表示算法求解精度越高。在表1中，步长α=0.4是FA算法的固定步长，其余的步长均为动态递减步长且满足收敛公式(8)。测试结果表明动态步长在求解精度和收敛速度方面明显优于固定步长。本文提出的双曲线递减步长萤火虫算法迭代次数最少、收敛速度最快，在收敛性能方面明显优于上述文献改进的动态步长。

图3 双曲线动态步长收敛效果图Fig.3 Convergence effect diagram of hyperbolic dynamic step

表1 测试自适应步长收敛效果

2.2 混沌序列优化FA

传统FA算法采用随机方式对萤火虫种群个体位置加以初始化，导致萤火虫个体位置分布不均匀、个体之间差异性较大及种群的多样性降低，仅有少数优秀个体在移动和迭代过程中被保留下来。为了解决这个问题，笔者将Logistic混沌序列引入传统FA算法初始化阶段，并结合具体问题相应地修改从混沌空间到搜索空间的映射方式。混沌序列具有独特性质即有界性、随机性、规律性和遍历性，混沌变量类似于随机变量，但能够按照指定规律在搜索范围内不重复地遍历所有状态以确保生成搜索能力较强、分布较广且均匀的初始种群，增强初始种群多样性及其在搜索空间分布质量，达到避免算法陷入局部极值点，提高算法求解精度。

采用Logistic映射函数生成混沌序列为

(10)

混沌序列初始化萤火虫种群的基本步骤：

①随机对D维空间内萤火虫进行初始化得到一个混沌变量；

②按照式(10)递推M-1次，最终得到的M个混沌变量；

③将产生混沌变量按式(11)映射到目标函数的搜索空间，从而得到M只萤火虫初始位置，公式如下：

xi,d=min{|LB|,|UB|}·yi,d

(11)

式中：LB和UB分别表示搜索空间第d维的下限和上限，yi,d是第i只萤火虫相对应的第d维混沌变量，xi,d是第i只萤火虫在搜索空间第d维的坐标值。图4混沌序列初始化萤火虫算法的流程图；图5直观得出采用混沌序列优化的萤火虫算法，其种群分布的均匀性要明显优于随机方式初始化的萤火虫种群。引入霍普金斯统计量[23](Hopkins statistic)定量衡量初始种群分布均匀性，如果霍普金斯统计量H接近于0.5表示数据分布均匀，如果H接近于1表示数据有聚类趋势，定义霍普金斯统计量H：

(12)

分别计算两种初始化萤火虫的方式的霍普金斯统计量H，采用混沌序列初始化的H1=0.485 6，而随机方式初始化的H2=0.889 4，证明在种群分布均匀性方面上混沌序列初始化优于传统随机方式。

图4 混沌序列初始化萤火虫算法的流程图Fig.4 Flowchart of chaotic sequence initialization firefly algorithm

图5 随机方式与混沌搜索初始化萤火虫位置对比图Fig.5 Comparison chart of firefly position initialized by random method and chaotic search

图6 测试函数的三维图Fig.6 Three dimensional diagram of test function

2.3 改进萤火虫算法和传统萤火虫算法比较

3 CAFA反演波阻抗

3.1 波阻抗反演

以褶积模型为基础进行波阻抗反演，地震记录的数学模型为：

S(t)=W(t)×R(t)

(13)

式中:S为地震记录序列；W为地震子波序列是地震记录的基本单位；R为地震反射系数序列，第i层的反射系数如下：

(14)

式中：ρi、vi是第i层的密度和速度，zi是第i层的波阻抗。由公式(13)和(14)可知，地震记录与波阻抗为非线性关系，因此求取波阻抗是一个求解非线性组合优化问题。目前模拟退火算法[24]、遗传算法[25]、粒子算法[26]和差分进化算法[27]等智能优化算法已经在反演波阻抗方面得到了成功应用，为萤火虫算法在波阻抗反演研究提供一些理论支持。笔者以改进萤火虫算法(CAFA)研究为基础，将CAFA算法和地震记录波阻抗反演相结合，验证CAFA算法在地球物理反演方面的可行性。使用最小二乘法建立萤火虫算法反演波阻抗的目标函数，建立的目标函数如下：

(15)

式中：M为空间采样点；W为Ricker子波；R(t)为所求的反射序列；S(t)为实际地震记录。CAFA经过多次迭代后得到的反射系数和Ricker子波褶积与实际地震记录的残差平方和达到全局极小值，再利用CAFA算法得到反射系数来递推各层的波阻抗。在反射系数递推波阻抗过程中误差会不断积累，若反射系数存在较大误差会导致递推的波阻抗曲线越来越偏离实际的波阻抗曲线。

3.2 理论模型试算

为检验CAFA算法反演波阻抗的可行性，设计一个五层水平地层理论模型(图 8)，萤火虫种群规模n=50,光强吸收系数β=1,初始步长α=0.2，自适应步长选取双曲线递减步长。图8 理论模型各层介质速度为1 600、2 000、2 300、1 900、2 600 m/s，各层介质密度为1 500、1 900、2 000、1 800、2 400 kg/m3。图9水平地层介质波阻抗模型，该模型各层介质波阻抗为2.40、4.60、3.06、1.82、3.42，单位为106kg·m-2·s-1。图10是水平地层模型正演得到单道合成地震记录，其中子波选取Ricker子波，子波主频为50 Hz，采样率为2 ms。从图11可知，CAFA反演波阻抗与模型波阻抗吻合程度更好，能够反演理论模型的各个层位波阻抗变化情况。图12 中CAFA算法和FA算法分别反演合成地震记录与模型合成地震记录都具有较好吻合，且CAFA算法的求解精度明显高于FA算法(CAFA算法求解精度为0.000 284、FA算法求解精度为0.005 320)。分别在原始合成地震记录中加入3%、5%和10%随机噪声验证CAFA算法的抗噪能力。图13、图14分别是波阻抗和合成地震记录反演的结果，含3%、5%的随机噪声时，反演结果误差较小能够真实的反映有效地层信息。含10%的随机噪声时，误差较大，但反演的总体变化趋势仍然正确，这表明CAFA具有较强的抗噪能力。

图8 五层水平地层Fig.8 Five horizontal strata

图9 模型的波阻抗序列Fig.9 Wave impedance sequence of model

图10模型的合成地震记录Fig.10 Synthetic seismic record of model

图11 CAFA和FA反演波阻抗Fig.11 Inversion wave impedance of CAFA and FA

图12 CAFA和FA反演合成地震记录Fig.12 Inversion synthetic seismic records of CAFA and FA

图13 反演波阻抗效果(含3%，5%，10%随机噪声)Fig.13 Inversion of impedance effect(including 3%, 5%, and 10% random noise)

图14 反演合成地震记录(含3%，5%，10%随机噪声)Fig.14 Inverted synthetic seismic records(including 3%, 5%, and 10% random noise)

4 结论

(1)双曲线递减步长相比于传统固定步长，不但提高算法求解精度(7.684 9e-07)，而且大幅度地提升收敛速度(1.3 s)，避免了算法收敛后期反复振荡现象。

(2)采用混沌序列替换随机方式初始化萤火虫种群能够提高初始种群分布质量，有效地解决算法易陷入局部极值问题。

(3)利用CAFA算法反演水平地层的波阻抗模型，得到误差较小的反演波阻抗和反演合成地震记录。在理论模型中分别加入不同强度的随机噪声，验证本文提出的CAFA算法具有较好的抗噪能力。