0 引言

瑞雷波勘探(李庆春等, 2006)具有无损、高效、经济等特点，越来越受到浅地表地球物理和地质工程界的重视.与该技术密切相关的瑞雷波是一种由纵波和横波干涉形成，并沿着自由表面传播的波.具有衰减速度慢、信噪比高、抗干扰能力强、以及在层状介质中的频散特性等特点，已被广泛应用于获取浅地表、岩石圈或地幔的横波速度结构信息(胡家富等, 1999; 黄忠贤等, 2014; 潘佳铁等, 2014)，获取近地表地层的品质因子(Xia, 2014)，路基压实度(刘江平等, 2009)等.诸多学者也进行了相关研究，如瑞雷波反演(石耀霖和金文, 1995; Xia et al., 1999; 陈祥和孙进忠, 2006)，瑞雷波数值模拟研究(Xu et al., 2007; 周竹生等, 2007; 熊章强等, 2008)，利用群速度估计浅地表横波速度(鲁来玉等, 2014; 余大新等, 2014)，岩洞的探测和面波的绕射(赵成刚等, 2008)，地形对瑞雷波传播的影响(周红和陈晓非, 2007)等.

瑞雷波勘探可分为三个过程：面波数据采集(Zhang et al., 2004)、拾取频散曲线(Luo et al., 2008)和频散曲线反演(张碧星等, 2000; Lu et al., 2007).在正确地拾取瑞雷波频散曲线之后，频散曲线的反演是获取可靠的近地表横波速度剖面的关键步骤(Zhang and Chan, 2003).多数局部线性化优化方法已被用来反演瑞雷波频散曲线，并且在各种商业软件中占主导地位，如最速下降法，L-M算法的最小二乘法等.然而，与大多数其他地球物理问题一样，瑞雷波频散曲线反演是高度非线性、多参数、多极值的.局部线性化的方法易收敛到局部极小值，是否成功在很大程度上取决于初始模型的选择以及偏导数计算的准确性.因此，可以克服这种限制的全局优化方法已被广泛地应用于瑞雷波频散曲线的反演中，如蒙特卡洛法、遗传算法(石耀霖和金文, 1995)、模拟退火法(崔建文, 2004)、粒子群算法(Song et al., 2012)等，但这类算法易出现收敛速度慢、早熟收敛以及搜索精度低等问题.

蝙蝠算法(Bat Algorithm，BA)(Yang, 2010)和萤火虫算法(Firefly Algorithm，FA)(Yang, 2010)均为新颖的仿生群智能优化算法.自然界中蝙蝠利用一种声呐探测猎物、避免障碍物，BA由此生物学特性发展而来.该算法具收敛速度快、精度高等优点，但存在易陷入局部极小值等不足.FA是模拟自然界中萤火虫发光的生物学特性发展而来，该算法具有稳定，全局搜索能力强等优点，但存在后期搜索精度低等不足.鉴于此，本文结合二者的优点，设计了一种新的优化策略，称其为FBA，即在前期采用稳定、全局搜索能力强的FA对种群进行优化，后期采用收敛速度快、搜索精度高的BA对前期FA生成的种群继续进行优化.为进一步增强BA算法的局部搜索能力，本文借鉴粒子群算法(PSO)对BA做了改进，在BA中引入了惯性权重，可形成带有惯性权重的BA(WBA)，继而形成带有惯性权重的FBA(WFBA).

为了评价WFBA对瑞雷波数据反演的有效性与稳定性，本文首先对三个典型地质模型的含噪声与不含噪声的理论频散曲线进行反演，并分析了算法的抗噪能力.然后单独将WBA、FA、PSO以及不含惯性权重的FBA对瑞雷波数据进行反演，并与WFBA进行了对比分析，说明了WFBA相对于WBA、FA、PSO以及FBA具有更稳定、收敛速度更快，反演精度更高的优越性.最后反演了来自美国怀俄明州某地区的实测资料(Xia, 2014)，以检验WFBA对瑞雷波频散曲线反演的适用性.理论模型试算与实测资料分析的结果表明，WFBA具有稳定、快速、精度高、适用性强的特点，可有效地对瑞雷波频散曲线进行定量解释.

1 萤火虫算法和蝙蝠算法的基本原理 1.1 萤火虫算法基本原理

FA是以模拟萤火虫的群体行为而构造出的群智能优化算法.用模型参数模拟萤火虫个体，将搜索和优化过程模拟成萤火虫个体之间吸引或移动的过程，将所需优化问题的目标函数度量成萤火虫个体所处位置的优劣，将萤火虫个体的优胜劣汰过程类比为搜索和优化过程中用好的可行解取代较差可行解的迭代过程.FA主要包含两个要素：荧光亮度与吸引度.其中，荧光亮度体现了萤火虫所处位置的优劣并决定其移动方向，吸引度决定了萤火虫移动的距离，通过荧光亮度和吸引度的不断更新实现目标优化.萤火虫算法的优化原理描述如下(Yang, 2010)：

定义1：萤火虫的相对荧光亮度为

(1)

定义2：萤火虫的吸引度为

(2)

定义3：萤火虫i被吸引向萤火虫j移动的位置更新公式为

(3)

式中，I₀为萤火虫的最大荧光亮度，与目标函数值f(x)相关，目标函数值越优，I₀越大；γ_F为光强吸收系数，理论上，γ_F∈[0, ∞]，但在实际中一般设置为模型搜索空间长度S_k的平均值，模型搜索空间长度S_k表达式为S_k=up_k-low_k, k=1, …, d，其中up_k为x第k个分量的上界，low_k为x第k个分量的下界，d为x的变量个数；r_ij为萤火虫i与j之间的空间距离；β₀为最大吸引度，一般设β₀=1.x_i、x_j分别为萤火虫i和j所处的空间位置，x_{i, k}为x_i的第k个分量，x_{j, k}为x_j的第k个分量.rand为[0, 1]上均匀分布的随机因子.α_F为步长因子，α_F=cS_k，c∈[0, 1]，在本文中，c随迭代次数按指数规律递减：

(4)

其中iter为当前迭代次数，Maxiter为最大迭代次数，c₀∈[0, 1]，本文中c₀=0.2.

萤火虫优化算法流程如下：

(1) 种群初始化：在模型搜索空间中随机生成一组萤火虫初始位置x_i(i=1, 2, …, n)，x_i=(x_i1, …, x_id)其中n为萤火虫数目，d为种群的变量个数.参数初始化：

(2) 根据初始化萤火虫的位置，计算萤火虫的目标函数值作为各自最大亮度I₀.

(3) 由式(1)(2)计算群体中萤火虫的相对亮度I和吸引度β，并对群体的相对亮度I进行排序，找出当前最优值和最优解.

(4) 对所有萤火虫个体进行遍历，如果I_j＞I_i(i, j=1, …, n)，则根据式(3)更新萤火虫的空间位置，使萤火虫i向萤火虫j移动.

(5) 当达到最大搜索次数则转(6)；否则，搜索次数增加1，转(3)，进行下一次搜索.

(6) 输出全局最优值和最优解.

1.2 蝙蝠算法基本原理 1.2.1 基本的蝙蝠算法

蝙蝠算法的主要步骤可以描述如下(Yang, 2010)：

(1) 种群初始化：在模型搜索空间中随机生成一组蝙蝠初始位置x_i(i=1, 2, …, n)，x_i=(x_i1, …, x_id)和初始速度v_i(i=1, 2, …, n)，v_i=(v_i1, …, v_id)，其中n为蝙蝠的种群数量.参数初始化：目标函数f(x)，最大脉冲音量A₀，最大脉冲速率R₀，搜索脉冲频率范围[f_min, f_max], 音量衰减系数α_B，搜索频率的增强系数γ_B.

(2) 根据初始化蝙蝠的位置x_i和v_i计算适应度值，并找到当前最优解x^*；

(3) 蝙蝠的搜索脉冲频率、位置和速度更新.种群在进化过程中，每一代个体的搜索频率、速度和位置按如下公式进行变化：

(5)

(6)

(7)

其中，φ∈[0, 1]是一个随机数；f_i是蝙蝠i的搜索脉冲频率f_i∈[f_min, f_max]，v_i^t、v_i^t-1分别表示蝙蝠i在迭代次数为t和t-1时的速度，x_i^t、x_i^t-1分别表示蝙蝠i在迭代次数为t和t-1时的位置.

(4) 生成随机数Rand，Rand∈[0, 1]，如果Rand＞R_i^t(R_i^t为蝙蝠i在迭代次数为t时的脉冲速率)，则对当前最优解进行随机扰动，产生一个新解.在本文中，其扰动方式为x_i^t=x^*+a×(Rand-0.5)，其中，a与萤火虫算法中的步长因子α含义相同.

(5) 生成随机数Rand，如果Rand＜A_i^t(A_i^t为蝙蝠i在迭代次数为t时的脉冲音量)，且f(x_i)＜f(x^*)时，则接受步骤(4)产生的新解x_i^t，然后按照如下公式对R_i^t+1和A_i^t+1进行更新：

(8)

(9)

其中，0＜α_B＜1，γ_B＞0.当t→∞时，A_i^t→0, R_i^t→R₀，本文中α_B=0.9，γ_B=0.9，A₀=1，R₀=0.5.

(6) 对所有蝙蝠的适应度值f(x)进行排列，找出当前最优解x^*和最优值f(x^*).

(7) 当达到最大搜索次数则转(8)；否则，搜索次数增加1，转(3)，进行下一次搜索.

(8) 输出全局最优值和最优解.

1.2.2 改进的蝙蝠算法

在BA中，如果用随机的参数来替换变化的频率f_i，并设置A_i=0，R_i=1，BA的速度位置的更新方程则变为与基本的PSO类似的形式.在基本的PSO中加入惯性权重能增强PSO的局部搜索能力(Elegbede, 2005)，故本文对BA做了相应的改进，在速度位置更新方程中引入了惯性权重ω，可形成带有惯性权重的BA(WBA)，改进后的速度位置更新方程如式(10)(11)所示：

(10)

(11)

与在PSO中的作用类似，惯性权重ω能控制前面的速度对当前速度的影响，较大的ω可以加强BA的全局搜索能力，而较小的ω能加强局部搜索能力，基本的BA是ω=1的情况，因此在后期相对缺少局部搜索能力.本文提出的反演策略主要利用了BA的局部搜索能力，故采用了较小的ω，并使其线性递减，这样可进一步增强BA的局部搜索能力，提高求解精度.线性递减的惯性权重ω如式(12)所示，

(12)

其中，iter为当前的迭代次数；Maxiter为最大迭代次数；ω_start和ω_end分别为初始的和终止的惯性权重.在本文中ω_start为0.6，ω_end为0.4.

2 理论模型试算

考虑到大多数情况下，基阶波是能量较强、较易观测到，且在实际中也是较广泛的，故本文采用基阶波频散曲线来反演横波速度结构.Xia等(1999)的研究表明，与瑞雷波频散曲线特征变化关系最密切的是地层的横波速度，其次是地层厚度、纵波速度和地层密度.因此，为减少反演时的计算量及其他参数的影响，本文只反演横波速度和厚度，而泊松比、密度这些参数假定是已知的.另外，在实际应用中，通常是没有足够的先验信息来估计浅地表的横波速度和厚度的.为了更符合实际情况，本文使用了较宽的模型参数搜索范围.在后面所有的理论模型测试中，搜索区域的下限和上限设定为与真实值相差50%或更多的值.我们使用相同的WFBA反演参数对理论模型的无噪声、含噪声数据以及实测数据进行10次独立反演测试，用标准差来评价算法的稳定性.反演参数设置为：种群数量为5×d，总迭代次数为100，其中FA和WBA的迭代次数分别为40和60.

瑞雷波频散曲线反演的过程是一个求解目标函数最小值的优化问题.反演的目标函数为实测的与理论的相速度值的均方差(RMS)，即

(13)

式中，V_R^obs为实测瑞雷波相速度；V_R^cal为理论计算的瑞雷波相速度，M为频点数；频散曲线的正演模拟采用的是快速矢量传递算法(凡友华等, 2002).

为了检验和评价WFBA的有效性和稳定性，使用了三个理论地质模型.这些模型可用于模拟浅层工程勘察中经常遇到的理论地质模型.如表 1所示，模型A是一个横波速度递增型的四层地质模型；表 1所示的模型B是一个含有低速软弱夹层的四层地质模型，比如常见的路面结构；表 1所示的模型C是一个含有高速硬夹层的四层地质模型.对于模型A和C，我们在5~100 Hz的频率范围内计算基阶模式频散曲线.但对于模型B，我们只在5~70 Hz范围内计算基阶模式频散曲线，因为实际上在高于70 Hz的频率范围内，模型B的基阶波是观察不到的.

2.1 无噪声理论数据反演分析

首先采用WFBA对无噪声的理论数据进行反演，模型A、模型B和模型C的反演结果如图 1与表 2所示(表中数据均保留两位小数).可以发现在没有足够的横波速度和厚度的先验信息的情况下，WFBA反演模型的模型响应(图 a1，b1，c1中的实线)几乎和真实模型的频散曲线完全拟合(图a1，b1，c1中的实心点)，且各个模型的真实值均被WFBA精确地反演和重建，如图 (a1，b1，c1)所示.模型A、模型B和模型C这3个模型的模型参数平均相对误差分别为0.17%，0.13%，0.16%，平均标准差分别为1.85，0.79，0.49.由此可见，WFBA可以有效并精确地对理论瑞雷波频散曲线进行反演.

2.2 含噪声理论数据反演分析

在实际应用中，拾取的瑞雷波相速度不可避免地存在噪声(Zhang and Chan, 2003)，数据噪声会导致反演算法的不稳定，还会导致反演算法找到的解不是真正的解，为了检验和评价含噪声的瑞雷波数据对WFBA性能的影响，在前面设计的3个理论模型中分别加入10%的白噪声作为随机扰动，然后分别进行反演.加噪声后的频散曲线如图 2(a1, b1, c1)中的实线所示，对加噪声的数据进行反演后，获得的反演曲线如图 2(a1，b1, c1)中的虚线所示.从图中可以看出WFBA在对含有噪声的数据反演时仍具有较强的稳定性，反演模型的频散曲线均能较好地拟合理论值的频散曲线.模型A、模型B和模型C这3个模型的模型参数平均相对误差分别为5.27%，3.97%，4.47%，平均标准差分别为2.40，2.26，1.49，如表 2所示.由此可见，WFBA反演瑞雷波频散曲线具有良好的稳定性.

为了检验WFBA的收敛性，需要分析在算法迭代过程中，最小目标函数值是否随着迭代次数逐渐减小.以模型C为例，将模型C的不含噪声和含噪声的理论数据反演迭代过程中最小目标函数的变化进行了分析，如图 3所示.由图可见，最小目标函数值在前20次迭代中快速收敛，然后逐渐收敛到零值附近(实心圆点)或某一常数附近(实线)，这表明WFBA已经完成了对最优解的搜索.由此可见在无噪声和有噪声的情况下，WFBA反演瑞雷波频散曲线均具有很好的收敛性.

2.3 含多模式理论数据反演分析

考虑到在某些含低速软弱夹层等特殊地质结构中，高阶波较基阶波在某些频段(如高频段)有时可能会更加发育，此时利用多个模式来进行反演是有必要的，且在反演时需注意模式的正确识别，将理论曲线与实测曲线相对应的模式进行拟合.本节以含有低速软弱夹层的模型B为例，在原含有噪声的基阶波数据中加入高阶模式数据，如图 4a中实线所示.然后利用WFBA对此含有多个模式的数据进行反演，其反演结果如图 4和表 3所示.在图 4a中，反演模型对应的频散曲线与各模式频散曲线均拟合较好；且各个模型参数的真实值均被较好地反演和重建，如图 4b所示.相对于单纯地反演基阶波数据，反演多个模式的数据能降低瑞雷波频散曲线反演的多解性.本文将仅含基阶模式数据和含多模式数据的反演结果进行了对比，发现后者相对于前者，模型参数反演平均误差从3.97%降低至2.65%，尤其是第三层厚度，其反演误差从14.65%降低至5.79%；平均标准差从2.26减小至0.05.由此可验证，利用多个模式数据来进行反演的确能降低反演的多解性，使得反演结果更准确可靠；也说明了WFBA不仅可用于基阶模式的频散曲线反演，也适用于多个模式的频散曲线反演.

2.4 WFBA与WBA、FA、PSO和FBA的对比分析

为了验证本文提出的WFBA的搜索性能的确有所提升，本文将WFBA与单独使用WBA和FA时的搜索性能进行了比较，并与典型的群智能优化算法PSO以及不含惯性权重的FBA进行了对比.以模型C为例，本文采用了WFBA、WBA、FA、PSO和FBA对模型C的不含噪声的理论数据进行10次独立反演测试.为了更方便比较WFBA和其他算法的性能，本文使用了相同的搜索空间、群体大小和迭代次数，WBA、FA的相关参数设置与WFBA相同，FBA的参数设置相对于WFBA除不含惯性权重外，其余参数设置均相同，PSO中的自身经验加速系数c₁与社会经验加速系数c₂分别设置为2.8、1.3(Song et al., 2012).反演结果如图 5所示.由图 5b明显可见，在100次迭代内，WFBA比FA、PSO和FBA具有更快的收敛速度，尤其是在迭代40次使用了WBA之后更为明显，FA在迭代大约60次后无法进一步收敛，PSO在100次迭代内保持收敛趋势，但收敛速度缓慢.WBA与WFBA的收敛速度相当，但从总体上看，WFBA是优于WBA的.为了定量地说明问题，对WBA、FA、PSO、FBA和WFBA的反演的模型参数的平均相对误差、平均标准差以及反演成功概率做了统计，如表 4所示，由表可明显看出，无论从模型参数平均相对误差、平均标准差还是成功概率上看，WFBA均明显优于WBA、FA、PSO和FBA.

3 实测资料分析

为了进一步检验WFBA的适用性，利用WFBA对在美国怀俄明某地获得的地震数据(Xia, 2014)进行了分析.如图 6所示，图 6a为采集得到的地震记录，采用了48个8 Hz垂直分量检波器，道间距为0.9 m，最小偏移距为0.9 m，震源采用锤击震源.图 6b为由炮集记录提取的频散能量图，从基阶波能量上拾取实测数据(图中黑色实心圆点).

然后，采用WFBA对其进行反演，与反演理论数据类似，只反演横波速度和厚度，固定密度和泊松比.根据其测井数据，可将地层大致划分为5个层位，然后根据测井资料的横波速度范围粗略估计WFBA的模型搜索范围，并根据测井资料所反映的地下岩性来估计各层的泊松比和密度.模型搜索空间和泊松比以及密度的设置如表 5所示.

利用WFBA反演拾取的基阶波相速度曲线，其最小目标函数随迭代次数的变化如图 6b所示，由图可见，算法已在100次迭代内收敛.反演获得的该地区的横波速度模型(图 7c中实线)和测井资料(图 7c中带圆点的折线)基本吻合，尤其是浅地表部分，而且与Xia(2014)利用局部线性化优化方法反演的结果(图 7c中带方块的折线)相比，其在5~13 m的深度范围内与测井资料更为接近.图 7a中的实线表示根据反演结果计算的相速度，其与测量值(图 7a中的实心点)也十分吻合.由此可见，WFBA反演瑞雷波频散曲线具有良好的实用性.

4 结论

本文提出了一种新的基于萤火虫算法和带惯性权重的蝙蝠算法的优化策略(WFBA)，并将其应用到了瑞雷波频散曲线反演中，反演中设置了较宽的模型搜索范围来模拟实际应用中没有足够先验信息的情况.本文首先采用无噪声、含噪声的理论数据对WFBA进行了有效性与稳定性实验，然后将WFBA与WBA、FA、PSO以及不含惯性权重的FBA进行了对比，分析了WFBA的反演效果，最后用来自美国怀俄明地区的实测数据检验了WFBA的实用性.理论数据和实测资料的反演实验结果表明：

(1) WFBA兼备FA稳定，全局收敛性强和WBA收敛速度快，求解精度高的优点，是一个有效、稳定、适用性强的全局优化算法，具有良好的发展潜力，很适用于瑞雷波频散曲线反演或其他类似的地球物理反演问题.

(2) 在WFBA中，带惯性权重的BA比基本的BA在收敛速度和求解精度上均有所提高，以模型C为例，求解精度从0.96%提高到了0.16%.

(3) WFBA相对于单独使用WBA、FA以及典型的群智能优化算法PSO具有更好的稳定性、更快的收敛速度和更高的求解精度，从而可加快瑞雷波频散曲线反演的反演速度，并使反演结果更稳定、准确、可靠.