基于改进海洋捕食者优化算法和瑞雷波频散曲线的近地表地层参数反演

于涵，刘财，王典，赵鹏飞，鹿琪，李鹏

吉林大学地球探测科学与技术学院，长春 130026

0 引言

面波勘探是一种无损地对近地表进行探测的方法，在确定小尺度地质体的结构特征、对近地表介质进行描述、探测城市浅层地下空间等领域应用潜力巨大(Park et al., 1999; 周熙襄和王振国, 2004; Hayashi and Suzuki, 2004; 夏江海等, 2015；程光华等，2019; Li et al., 2021).随着近地表地震勘探的目标逐渐转向“精细化”，利用反演方法识别近地表地下介质的意义更加突出，因此，如何有效获得近地表场地特征，提取近地表地层结构信息已经成为研究热点.瑞雷波频散曲线隐含着地下地层的信息，利用瑞雷波频散曲线可以反演得到地下不同深度介质的属性(刘昊楠等，2020; 胡书凡等，2021).地球物理反演是典型的高度非线性、多参数、多极值的反问题和最优化控制问题(赵鹏飞和刘财，2021).从瑞雷波频散曲线反演地层参数也是一个典型的多目标优化求解问题.通过构建与瑞雷波频散曲线相应的反演问题，从而得到浅层的速度图像，有效获得近地表场地特征，可以更好地了解近地表结构(夏江海，2015).

目前，瑞雷波频散曲线的反演方法可分为局部线性化方法(Local Search Method)和非线性全局优化算法(Global Optimization Algorithms)两大类.局部线性化方法是基于梯度的计算方法，反演结果在很大程度上取决于初始模型的选择以及目标函数对模型参数的偏导数计算的准确性(Ganji et al., 1998; Xia et al., 1999).蒙特卡洛算法和广义模式识别算法是非线性全局优化算法的典型代表(Socco and Boiero, 2008；宋先海，2008；张志厚等，2022).近年来，元启发式优化方法广泛应用于频散曲线的提取过程中，主要分为两大类别：以遗传算法(GA)为代表的进化算法和以粒子群算法(PSO)为代表的群智能优化算法.多位学者在此基础上将更多的元启发式算法应用到瑞雷波频散曲线的反演中来，如模拟退火(SA)、萤火虫算法(FA)及其改进的优化算法(如算法步长改进、多种优化算法之间相互结合).通过分析近地表介质参数的反演结果，可以发现以上算法各有特点.遗传算法计算耗时较长，粒子群算法容易存在种群多样性趋同，多参数存在串扰现象，易陷入局部极小值(Beaty et al., 2002；Dal Moro et al., 2007; Song et al., 2015；于东凯等, 2018；蔡伟等, 2018).

海洋捕食者算法(Marine Predators Algorithm, MPA)是由Faramarzi等(2020)提出，是一种新的群体智能优化算法.它源于对海洋捕食者的觅食行为的简化和模拟，在反演过程中结合Lévy飞行及布朗运动，具有参数少、设置简单、易于实现、计算准确等优点，该反演方法对多参数反演具有极大的应用潜力，在解决实际工程问题方面表现出了优越的性能.

本文提出一种改进的海洋捕食者优化算法(ACMPA算法).针对初始化随机不均匀的现象，引入基于Cubic映射的混沌序列初始化，加强了搜索的随机性，产生的混沌映射函数值对海洋捕食者种群位置均匀初始化，从而有助于增强种群多样性和算法的全局搜索能力；引入自适应函数因子进一步提高迭代前期全局搜索能力及其在搜索空间的分布质量；引入精英等级制度，避免算法陷入局部极值点，以取得最佳处理效果，加强海洋捕食者算法的寻优精度和稳定性.模型及实际数据反演测试结果表明，与其他现有的优化算法相比，改进的海洋捕食者优化算法可以有效地提取夹低速地层构造的近地表介质信息，对于浅地表地层信息拾取较为精细，可以有效地反演近地表地层结构.

1 方法原理

1.1 海洋捕食者优化算法的基本原理及其改进

海洋捕食者算法是模拟自然界中物竞天择的生物学特性发展而来的.在一个种群中，根据海洋中捕食者的觅食行为，在算法优化过程中引入捕食者和猎物概念.海洋捕食者算法的核心主要是将布朗运动与Lévy运动相互结合，并引入海洋旋涡阶段影响，拥有捕食者存储记忆特性，均有效地避免了收敛“早熟”而陷入局部最优解现象.

(1)种群初始化：在模型搜索空间范围中随机生成模型X0：

X0=Xmin+rand(Xmax-Xmin)，

(1)

式中，Xmax、Xmin为搜索空间范围；rand为服从均匀分布[0,1]的随机数.

该算法分别引入了捕食者和猎物矩阵，通过公式(1)生成初始种群定义为猎物矩阵Prey，而最优解被命名为捕食者矩阵Elite，分别表示为：

(2)

式中，n为种群数量，d为优化问题的变量个数，Xi,j(j=1,2,…,d)为第i个猎物第j维空间位置，X1,j*(j=1,2,…,d)称为最优捕食者向量，精英矩阵由此向量重复n次形成.

(2)优化过程：调整种群位置更新猎物位置

海洋捕食者优化算法采用分类寻优策略，将游走行为主要分为三个阶段来更新捕食者矩阵位置：

第一阶段(前1/3迭代过程)：假定捕食者比猎物移动速度快，探索阶段进行空间搜索并捕获猎物，整体运动行为为布朗运动；

第二阶段：假定捕食者与猎物以相同速度进行移动，在搜索的过程中，猎物种群用于开发阶段，对搜索空间进行Lévy运动，跳出该捕食圈层；捕食者种群用于探索阶段，对探索空间进行布朗运动，对搜索空间进行随机游走；

第三阶段(后1/3迭代过程)：假定捕食者比猎物移动速度慢，捕食者矩阵进行Lévy运动.

随迭代次数增加移动步长变换公式如式(3)，即S(Iter)矩阵的第i行向量为：

(3)

随迭代次数增加，猎物矩阵Prey(Iter)的第i行的变更公式如下：

(4)

式中：P=0.5为常数量，用于优化产生随机数对移动步长的影响.R表示在[0,1]区间内抽取的一组d维随机向量，其他变量及其参数的含义与式(4)相同.CF是自适应步长参数，用于控制捕食者运动步长：

(5)

(3)涡流及鱼类聚集效应

为避免算法陷入局部最优解，引入涡流及鱼类聚集效应(FADs)的环境影响，即捕食者在80%搜索空间内游荡于猎物附近，剩余20%游动于不同区域从而提高全局搜索，增强算法的收敛能力，并将其表示为：

(6)

式中，FADs=0.2，是影响涡流优化过程的概率，U为包含0和1的n维向量，当随机解小于0.2时，U变为0，否则变为1.R为[0,1]内随机数.Preyrand 1和Preyrand 2为猎物矩阵中的两个随机索引向量.

1.2 改进的海洋捕食者优化算法(ACMPA)

1.2.1 混沌初始化

MPA算法在初始过程阶段采用随机分布，会使初始种群个体分布不均，多样性较差，对算法的求解效率有很大的影响，甚至导致寻优失败.为增强种群分布的均匀性和遍历性，本文将混沌序列进行初始化.不同混沌映射的搜索能力在收敛精度上存在差异，相比Logistic混沌映射，Cubic混沌映射具较好的混沌遍历性，有寻优速度快、精度高的特点(王光义和袁方，2013；Feng et al., 2017).根据Cubic混沌映射，在保证较优个体的同时，令个体尽可能均匀分布，可有效平衡局部寻优和全局寻优.

Cubic混沌映射(图1)将函数值映射到(0,1)的区间，表达式如下：

图1 Cubic混沌映射

(7)

式中，i为种群个数，Ci∈(0,1)为第i个混沌映射函数值；取C0=0.3和μ=2.595的混沌参数，Cubic混沌映射在(0.1)之间具有较好的遍历性(张孟健等，2021).

通过Cubic混沌映射函数生成的初始种群，将公式(1)改进如下：

X0=Xmin+Ci(Xmax-Xmin).

(8)

1.2.2 自适应步长因子

在优化过中为提高算法在迭代前期及中期的全局搜索能力，引入自适应步长因子A，并将其与CF(式(5))进行比较(图2)，自适应步长因子A的变化较为敏感，可以自适应地随迭代次数增加变换移动速度，公式描述如下：

图2 步长变化比较图

(9)

利用自适应步长因子对于猎物矩阵的布朗运动进行改进，将部分公式(4)改进如下：

(10)

1.2.3 精英竞赛制度

为改善种群后期多样性降低问题，避免出现早熟现象，引入一种精英等级制度.比较前三个最优解均值与当前最优解的适应度函数值，择优选取最优解，引导其余种群个体进行变换，帮助算法跳出局部最优，克服早熟的不足，表示为：

(11)

1.3 改进的海洋捕食者算法的实现流程

改进的海洋捕食者算法可以归纳如下四个步骤，流程图如图3所示：

图3 计算流程框架

(1)初始化种群：设定猎物总数，初始混沌序列化每个猎物的初始位置.

(2)计算猎物的适应值：根据目标函数结合精英等级制度，计算并比较前三个最优解均值与最优解之间的适应函数值，不断对捕食者进行替换获得更优解，最后将当前迭代最优解与迭代前进行比较，更新猎物群体及捕食者个体位置.

(3)进入迭代过程，算法在迭代初期时，加大全局搜索范围，结合自适应步长因子改进的布朗运动，增强猎物的位置扰动范围；在迭代过程中，运用精英竞赛选择机制，避免陷入局部最优解；在迭代后期，结合布朗运动与Lévy运动，提升算法的精度和可靠性.

(4)终止判断：判断目标函数是否满足终止条件，输出最终结果.

2 反演试算

2.1 目标函数的建立

本文利用快速矢量传递算法计算瑞雷波频散曲线，包含非线性隐函数关系的瑞雷波频散曲线方程表示为F(f,VR,VS,VP,ρ,h)=0，式中：f为频率,VR为瑞雷波相速度，VS,VP,ρ,h分别代表地层模型中的横波速度、纵波速度、密度及厚度参数，各参数对于瑞雷波频散曲线的变化均有不同程度的影响(Renalier et al., 2010；Thota et al., 2021).利用频散曲线方程的介质参数反演是一个求解多目标函数的优化问题，采用改进的海洋捕食者优化算法进行处理，找到一组最优的地下地层模型.最优化问题的目标函数为：

(12)

2.2 理论模型设计及反演

实际近地表地下结构与所覆盖沉积层密切相关，近地表勘探工程的主要目标为土层，地表盖层厚度在30 m范围以内的平均剪切波速度是面波勘探的主目标.研究岩土层的介质参数，在工程应用中具有重要意义.在实际工程勘探应用中，通常遇到检测近地表异常地质结构(如含有淤泥质软土层、层间断裂、空洞等含低速夹层地质目标体)、含矿层、高速的盐丘体等含高速夹层的地质目标体等实际需求(王栩等，2021; 徐力群等，2022).根据近地表主要介质物性参数关系及实际面波勘探应用总结，瑞雷波频散曲线与地层横波速度及厚度密切相关，将地层密度及泊松比作为已知参数进行设置(Feng et al., 2005；夏江海等2015；Shen et al., 2016; 林志平等，2015；赵浩等，2018；秦彤威等，2021).本文利用表1给定的地层模型参数设计含有不同夹高(低)速层的地质模型.

考虑到在实际面波勘探过程中，面波基阶模态频散能量较强，应用更为广泛，故本文采用基阶模态瑞雷波频散曲线反演水平层状介质的物性参数.在建立模型过程中，各模型包含多种地质情况，以表1为参数依据，建立四层水平地质模型：图4a1的速度递增模型A；图4a2的第三层为低速夹层或高速表层速度异常模型B；图4a3的第二层为低速夹层模型C；图4a4的含高速夹层模型D.这四种模型对应的瑞雷波频散曲线如图4b所示.运用ACMPA优化算法用不同模型的频散曲线反演地层参数以验证优化方法的效果.

图4 不同地质模型示意图及其基阶模式瑞雷波频散曲线

表1 近地表地层参数

在模型数据B的测试中，为了对比ACMPA算法的效果，本文使用标准MPA方法、粒子群算法(PSO)、遗传算法(GA)、蝴蝶优化算法(BOA)、秃鹫优化算法(BES)、基于螺旋上升机制的鲸鱼优化算法(WOA)、基于四种等级寻优政策的灰狼优化算法(GWO)、基于链状群体模型的樽海鞘优化算法(SSA)、基于振荡反应的黏菌算法(SMA)等多种优化算法(Mirjalili et al., 2014, 2017; Mirjalili, 2016; Arora and Singh, 2019; Li et al., 2020)进行反演，种群数量设为40，分别比较迭代次数在50、200次的反演效果.从收敛曲线图5可以看出，ACMPA优化算法在反演过程中可以随着迭代次数的增加不断有效寻优，其精度更高，迭代收敛效果更加稳定，体现该算法在全局寻优和局部寻优的平衡性，迭代次数较少亦可以进行有效寻优.运用其特有的精英矩阵与猎物矩阵交互的过程，可以有效跳出局部最优解，加入自适应步长改进的布朗运动，可增加其在迭代前期过程中的全局寻优效果.

图5 模型B不含噪声数据ACMPA与多种优化算法收敛曲线对比

本文使用迭代次数200，种群个数40的ACMPA优化算法，在后面所有的理论模型测试中，使用相同的反演参数对理论模型的无噪、含噪数据以及实测数据进行反演测试，了解改进的ACMPA算法的稳定性与抗噪性.为了定量说明改进后ACMPA反演算法的效果，在迭代200次时，对MPA优化算法及ACMPA优化算法反演模型B的地层参数的相对误差及标准差进行统计，示于表2，从模型B的反演数值误差可以看出，MPA优化算法对于厚度参数反演误差值差异较大，随机度较高，ACMPA优化算法具有较强的鲁棒性和稳定性，明显优于MPA优化算法.

表2 地质模型B真值、反演对比结果

图6分别展示4种不同地层模型参数设置、基于ACMPA优化算法的瑞雷波频散曲线反演拟合效果及地层模型拟合效果图.表3为不同地质模型通过ACMPA优化算法反演得到介质参数及相对误差统计表.在近地表地层参数反演过程中，ACMPA优化算法对于模型速度的变化更敏感，可以有效反演至参数个位数值，在速度及地层厚差异较小的情况下，亦可以有效地进行反演，反演求解更高，反演误差值较小.在不同模型反演测试过程中，可以发现，当含有高速异常层并且速度与上下地层速度差异较大时，以模型4为例，相对于其他地质模型而言，厚度反演误差相对较大，对反演效果影响较为明显.

表3 地质模型真值、反演结果及相对误差统计表

图6 ACMPA优化算法对四种模型的反演结果

由于采集数据中包含的部分噪声会影响频散曲线的拾取，导致反演结果的多解性，反演方法的抗噪性是反演方法的一种重要的能力.为了测试本方法的抗噪性能，在频散曲线数据中加入噪声，使频散曲线产生扰动，然后测试ACMPA反演方法反演效果.低速异常地层是工程勘探过程中重要的监测目标体，研究含低速夹层模型更具有实用价值，故在模型B的基础上进行改进，分别设置第二层、第三层为低速层的地质模型，并基于此模型进行正演，在得到的频散曲线数据中加入10%的噪声.从图7可以看到，ACMPA优化算法对于含有10%噪声的频散数据也可以有效地提取出地层的参数.ACMPA优化算法对于含噪频散数据的反演仍具有较好的稳定性，可以较好地拟合理论频散曲线，各模型之间误差值较小.夹低速层的较强的参数反演能力使该方法具有较好的应用前景.

图7 含噪声情况下的低速层模型反演效果

2.3 实际数据处理

为了进一步验证ACMPA优化算法反演效果，利用实际数据进行反演测试.实际数据的测量地点位于杭州东火车站附近(Mi et al., 2022).对采集的背景噪声地震数据进行了分析.沿两条交叉道路部署东北线和东南线，密集线性阵列间距约为 5 m，有47个地震道，长度超过 200 m，探测深度达100 m，时间采样间隔为0.016 s，得到由NE-142至NE-189台站(以NE-142作为虚拟源)生成的地震记录(图8a左)以及SE-136至SE-183台站(以SE-136作为虚拟源)生成的地震记录(图8a右).图8b—d分别为对来自NE方向及SE方向的多通道表面波地震记录运用频率分解法、相移法、高分辨率Radon方法进行频散能量成像图，其中从高分辨率Radon变换频散能量成像图中有效地提取基阶模态瑞雷波频散曲线，频散能量主要分布于5～15 Hz，运用ACMPA算法反演近地表地层介质信息.

图8 NE方向(左)及SE方向(右)实际地震记录与频散能量谱

运用2号井测井数据将初始勘探空间划分为10层，分别从NE方向及SE方向的地震数据中，运用频散能量谱成像的方法进行处理，从中提取出的瑞雷波频散曲线(白线)，如图9b所示，NE方向及SE方向提取出的瑞雷波频散曲线整体趋势较为一致，并运用NE方向及SE方向的实测频散曲线与ACMPA优化算法反演后的频散曲线进行验证.可以看出，在没有足够先验信息的情况下，ACMPA优化算法进行反演后获得的横波速度模型和测井资料基本吻合，在浅层地表部分可较为精细地体现出不同层厚横波速度的变化，根据该方法反演出的横波速度与实测测井数据的拟合效果较好(图9a)，反演计算得到的相速度频散曲线与拾取频散曲线亦有较好的拟合效果(图9b).由此可见，运用ACMPA优化算法由瑞雷波频散曲线反演地层参数具有较好的应用效果.

图9 应用ACMPA算法由实际数据进行反演的结果

3 结论与展望

本文提出了一种改进混沌初始化、自适应步长和精英竞赛机制优化的海洋捕食者算法，并将其应用于瑞雷波频散曲线反演地层参数中.通过对四层地层模型合成数据、含噪声频散曲线和实际面波数据进行反演测试，验证了ACMPA算法在处理由瑞雷波频散曲线反演地下介质参数的问题时，具有较好的种群多样性、收敛精度较高.本文得到以下结论与展望.

(1)改进的ACMPA优化算法相较于传统的ACMPA优化算法，修改种群初始化，引入自适应步长和精英竞赛机制优化，从而在求解反问题时表现出更好的稳定性和更高的求解精度.

(2)改进的ACMPA优化算法不依赖于初始模型，具有较强的全局搜索能力，从局部与全局寻优两个层面提高反演算法的收敛速度，在利用瑞雷波频散曲线反演地表物性参数的问题中，可获得较高精度的反演结果.

(3)在实际数据的测试中，改进的ACMPA优化算法的反演结果与测井数据结果较好吻合，且对于夹低速层的探测能力具有较好的应用效果，稳定性更强，反演精度较高，能够更好地提取出近地表地层参数信息.

(4)在今后研究频散曲线反演过程中，应考虑多参数反演可能性，例如地层泊松比、近地表土壤密度等参数一同反演.