BFO-PSO算法下的弹性波数值模拟

赵平起， 何书梅， 倪天禄， 赵明， 张家良， 吴吉忠,2，魏朋朋， 李闻达， 白文磊*

1 中国石油大港油田分公司， 天津 300280 2 东北石油大学， 大庆 163318 3 中国科学院地质与地球物理研究所， 北京 100029 4 中国科学院地球科学研究院， 北京 100029 5 中国科学院油气资源研究院重点实验室， 北京 100029 6 中国科学院大学， 北京 100049

0 引言

地震波正演模拟是地震勘探反演和成像的基础和关键，提高弹性波数值模拟的精度和效率具有非常重要的理论和工程意义.目前广泛运用于地震波数值模拟的方法主要是有限差分法(Chu and Stoffa，2012；Di Bartolo et al.，2012；Yan et al.，2016；Yang et al.，2017；Liang et al.，2018；He et al.，2019；杜泽源等，2019；Miao et al.，2020)，其适用于GPU并行计算，计算速度快并且易于编程实现.

应用有限差分算法进行弹性波数值模拟时，使用差分算子来代替微分算子，因而必然会产生数值误差，如果忽略这一误差，将会严重影响数值模拟的精度，导致弹性波数值模拟过程中出现数值频散问题(Chu and Stoffa，2012；王之洋等，2015).为了解决差分算子频散问题，可以在数值模拟时选择主频较低的子波，较小的离散网格间距，或者优化差分算子.然而，随着子波主频的降低，高频成分缺失，由此也将导致地震反演和成像的分辨率大大降低；而减小离散网格的间距则会增加计算量，对于大模型来说，由此带来的计算量、存储量问题进一步制约勘探技术的发展(王之洋等，2015；He et al.，2019；刘立彬等，2020).优化差分算子，提高低阶数差分算子的精度，可以在不增加计算量，保证计算效率的同时，减小差分算子对微分算子的逼近误差，从而避免数值模拟过程中的数值计算误差，因此，优化有限差分算子在弹性波动方程的数值求解中具有非常重要的意义.

Chu和Stoffa(2012)指出，根据Taylor级数展开，可以使用二项式窗函数截断伪谱法的空间褶积序列推导得出有限差分算子，然而由于截断效应存在，频谱泄露不可避免，为了减弱频谱泄露，可以选择优化的窗函数.应用窗函数法优化有限差分算子，其主要目标是设计一种通带窄，阻带衰减大的窗函数，以尽可能的减弱截断所导致的频谱泄露(王之洋等，2015；Wang et al.，2017；Ren et al.，2018)，然而，不可忽略的问题是，窗函数的通带宽度与阻带衰减是相互矛盾的，很难同时满足通带窄、阻带衰减大的要求.因此，另一种更常用的优化有限差分算子的方法是最优化方法(Liu，2013；Zhang and Yao，2013；Yang et al.，2017；He et al.，2019；Miao and Zhang，2020)，通过构造包含有限差分系数的目标函数，将有限差分算子的优化问题转化为多参数优化问题，进而应用最优化算法进行优化求解.Zhang和Yao(2013)首次使用模拟退火算法优化有限差分算子.He等(2019)将Remez交换算法应用于常规网格和交错网格的有限差分算子优化.Miao和Zhang(2020)基于L1范数构造了包含有限差分系数的目标函数，并应用交替方向乘子算法(ADMM)获得了优化的有限差分算子.

应用最优化方法优化有限差分算子，其目标是在保证精度误差的同时，获得最大的谱覆盖范围(Liu，2013；Zhang and Yao，2013；He et al.，2019)，常规有限差分优化方法往往难以兼顾精度和谱覆盖范围两个方面，尤其是在高波数情况下，难以满足数值模拟精度的要求.因此，需要寻求一种更适用于高维多参数问题的快速收敛且易于跳出局部极值的高效优化算法.

启发式优化算法是目前广泛应用的一种优化算法(雷秀娟等，2012；Beheshti and Shamsuddin，2013；Pillay，2016；Qu et al.，2015；Li et al.，2016；Hussain et al.，2018)，不同于传统优化算法，启发式算法通过在解空间范围内随机搜索以获得最优结果，是一种基于概率计算的领域随机搜索算法，该类算法主要包括：模拟退火算法(SA)，粒子群算法(PSO)，遗传算法(GA)，等等.相比于其他算法，PSO算法在迭代更新时利用了当前最优解的信息，从而可以更快收敛到最优解(Kennedy and Eberhart，1995；Hazra and Sinha，2011；Jadoun et al.，2015)，然而，PSO算法的全局搜索能力对初始参数的要求较高，当初始参数随机性较大时，算法很容易陷入局部极值，尤其是对于一些参数维度较高的问题、多峰问题和病态问题，会出现过早收敛甚至难以收敛的问题.细菌觅食算法(BFO)是一种仿生类算法(Passino，2002；李珺等，2013；Daryabeigi and Dehkordi，2014)，其模拟大肠杆菌觅食行为，算法主体结构包括三层循环，分别对应趋化、复制、驱散三个步骤，具有较强的并行搜索能力，同时，驱散操作加强了搜索的随机性，从而有助于提高算法的全局搜索能力，然而，BFO算法包含三层循环嵌套，使得算法结构相对复杂，引入了较多参数，导致算法的收敛速度较慢.因此，结合PSO算法和BFO算法的思想，使其优势互补，提升性能，改进算法可以在快速收敛的同时搜索到全局最优解，以同时满足有限差分算子对精度和频谱覆盖范围的要求.

本文提出一种BFO-PSO算法下的有限差分优化方法，并采用优化的有限差分算子进行弹性波数值模拟.针对PSO算法全局搜索能力较弱的问题，引入BFO算法中的趋化、复制、驱散三个步骤，形成BFO-PSO混合优化算法；构造包含有限差分系数的目标函数，并运用BFO-PSO混合优化算法求取全局最优解，获得优化的有限差分算子，通过理论频散分析，比较优化性能；进而应用此优化的差分方法分别在层状介质模型和复杂模型上进行弹性波数值模拟，对比分析合成地震记录.

1 理论分析

1.1 BFO-PSO算法下的有限差分算子优化

根据采样定理，连续信号f(x)的一阶空间导数可以表示为(Chu and Stoffa，2012)：

亚当·斯密认为，人类的“爱”是“有层次的”，而自由市场通过平等自愿的交换，让人们能够“以自私为目的，达成利他的结果，最终使每个人都互惠互利”。从制度经济学的角度看，破除以地方政府为代表的政府机会主义，重塑政治生态，有耐于持之以恒、与时俱进地优化“规则”。正如布坎南所指出的那样：“要改变一种游戏或竞赛的结果，改变参加竞赛的人并不重要，而改变竞赛规则最为重要。”

(1)

其中，Δx为空间采样间隔，Δx/π为Nyquist波数，fn=f(nΔx).

将x=0代入公式(1)并截断，可推导得出一阶空间导数的有限差分算子：

(2)

对公式(2)左右两端同时应用Fourier变换，可推导出：

(3)

其中，kx为波数.

根据公式(3)，可得出有限差分算子的频散关系：

(4)

式中，考虑波数趋于0的情况，满足：

进一步，应用最优化算法优化有限差分系数，可以构造目标函数为：

(5)

因此，本文结合BFO算法的全局寻优能力，将BFO算法的趋化、复制、驱散三个步骤引入PSO算法中的粒子速度和位置更新策略，以形成BFO-PSO混合优化算法，进一步应用BFO-PSO算法对包含有限差分系数的目标函数(公式(5))进行求解，从而得到精度误差较小且频谱覆盖范围较大的有限差分算子.

PSO算法是一种基于概率计算的领域随机搜索算法，种群中的每一个粒子都代表目标函数的一组解，在该算法中，每个粒子都以随机的方式向当前最优粒子的方向移动，从而使得整个种群向最优解的方向运动，进而实现在目标函数解空间范围内以随机搜索的方式求解最优解.常规PSO算法的粒子速度和位置更新策略如下：

(6)

为了克服PSO算法对初始参数的高度依赖性以及难以跳出局部极值的缺陷，本文在PSO算法的粒子速度和位置更新方式中增加BFO算法的趋化、复制、驱散三种步骤，以提高算法在邻域搜索时的随机性，形成BFO-PSO混合优化算法，从而提高PSO算法的全局搜索能力，降低算法对初始参数的依赖性，同时进一步提高算法的收敛速度.

首先引入趋化步骤，在粒子速度和位置更新过程中，粒子种群向着全局最优解的方向运动，但是，对于一些参数维度较高的问题、多峰问题和病态问题，粒子种群很容易趋向局部最优，从而导致算法全局搜索能力较差.为了避免这一问题，在粒子位置更新时增加随机方向性，使粒子随机地向一个方向运动或者翻转，从而增加粒子种群在全局范围内的搜索随机性，使其具有跳出局部最优的能力，增加算法收敛到全局最优解的概率.引入趋化步骤的粒子位置更新公式为：

(7)

其次是复制步骤，将BFO算法中细菌优胜劣汰的繁殖过程引入到PSO算法粒子速度和位置更新策略中，对种群中的每个粒子计算适应度值，并根据目标函数要求按照从优到劣的顺序对其排序，并将排序在末尾的一半种群粒子淘汰，而将剩余的种群粒子进行完全复制，使得种群大小不变.通过这一步骤，可以保证粒子种群在每一次迭代更新时都向着部分较优解的方向运动，从而有效提高粒子种群向最优解方向移动的速度，提高算法的收敛速度.

最后是驱散步骤，为了进一步放置粒子种群过早的收敛到局部最优位置，避免种群粒子聚集在局部极值位置处，按照一定的概率将种群粒子随机移动到解空间中的新区域，以此提高算法在邻域内搜索的随机性，从而增加粒子种群的随机性，避免种群过早的收敛到局部极值，提高全局搜索能力.

根据上述分析，将BFO算法的趋化、复制、驱散三种步骤引入PSO算法，改进PSO算法的粒子速度和位置更新策略，形成BFO-PSO混合优化算法，算法流程图如图1所示.

图1 BFO-PSO算法流程图Fig.1 Flowchart of the BFO-PSO algorithm

1.2 理论频散分析

根据公式(3)所示有限差分算子的频散关系式，分析比较BFO-PSO算法，PSO算法以及Remez交换算法(He et al.，2019)下的优化有限差分算子的数值频散曲线，比较其精度误差.

分别应用BFO-PSO算法和PSO算法对包含有限差分系数的目标函数(公式(5))进行优化求解，以获得BFO-PSO算法下的优化有限差分算子.本文中，BFO-PSO算法和PSO算法的初始参数设定如下：粒子位置表示有限差分系数，并进行随机初始化；粒子速度表示有限差分系数的调整方向；种群大小为N=300，最大迭代次数为T=1500，权重系数为ω=0.7，粒子驱散概率为0.25.

如图2、3所示，分别为BFO-PSO算法和PSO算法优化得到的一阶空间导数有限差分算子的数值频散曲线，其中，图2应用BFO-PSO算法，图3应用PSO算法，图中黑色实线表示不同阶数下常规有限差分算子的数值频散曲线.对比分析，可以明显看出，采用BFO-PSO算法得到的优化有限差分算子在保证精度误差的同时，其频谱覆盖范围远大于常规有限差分算子.8阶BFO-PSO算法优化的有限差分算子的频散曲线甚至接近16阶常规有限差分算子的频散曲线，两者的频谱覆盖范围几乎一致，而12阶BFO-PSO算法优化的有限差分算子的频散曲线则远远优于20阶常规有限差分算子的频散曲线，在12阶时就可以满足高波数的精度要求.

图2 不同阶数的BFO-PSO算法下的优化有限差分算子的频散曲线Fig.2 Dispersion curves of the optimized FD operator based on the BFO-PSO algorithm for different values of N (different orders)

图3 不同阶数的PSO算法下的优化有限差分算子的频散曲线Fig.3 Dispersion curves of the optimized FD operator based on the PSO algorithm for different values of N (different orders)

进一步分析比较BFO-PSO混合优化算法的全局搜索能力和收敛速度，图4给出了BFO-PSO算法与PSO算法优化获得的有限差分算子的数值频散曲线对比.可以看出，在8阶、12阶、16阶、20阶的情况下，BFO-PSO算法和PSO算法优化得到的有限差分算子，其数值频散曲线几乎完全一致，即两种优化算法得到的有限差分算子是一致的.然而，当差分算子的阶数继续增加到24阶时，应用BFO-PSO算法得到的优化有限差分算子在保持精度误差不变的同时，具有更大的频谱覆盖范围.基于这一分析，可以说明，引入趋化、复制、驱散三种步骤后形成的BFO-PSO混合优化算法很好的结合了PSO算法的快速收敛性和BFO算法的全局搜索能力，对于多参数优化问题，BFO-PSO混合优化算法可以快速收敛到较好的全局最优解，就本文而言，应用BFO-PSO算法可以很快的获得精度误差小且适应高波数情况的优化有限差分算子.

图4 BFO-PSO算法与PSO算法优化的有限差分算子的频散曲线对比Fig.4 Comparison of numerical dispersion curves between optimized FD operator based on the BFO-PSO algorithm and the PSO algorithm

图5为分别应用BFO-PSO算法和Remez交换算法得到的优化有限差分算子的数值频散曲线对比，其中，彩色实线分别对应不同阶数的BFO-PSO算法优化的有限差分算子的频散曲线，彩色虚线分别对应不同阶数的基于Remez交换算法的有限差分算子的频散曲线.对比可得，在保证精度误差控制在一定范围内的同时，BFO-PSO算法优化后的有限差分算子具有更大的频谱覆盖范围，更适用于高波数情况下的数值模拟.表1列出了应用BFO-PSO算法得到的优化有限差分算子的系数.

图5 BFO-PSO算法与Remez交换算法(He et al.，2019)优化的有限差分算子的频散曲线对比Fig.5 Comparison of numerical dispersion curves between the optimized FD operator based on the BFO-PSO algorithm and the Remez algorithm (He et al., 2019)

表1 BFO-PSO算法下的优化有限差分算子系数Table 1 Optimized FD coefficients based on the BFO-PSO algorithm

应用有限差分法进行数值模拟时，差分格式的稳定性是一个非常重要的问题.这里假定空间步长均匀，即Δx=Δz，介质最大波速为v，时间步长为Δt，给出常规网格一阶有限差分算子的简单稳定性条件(Lines et al.，1999)：

(8)

2 数值模拟与分析

2.1 层状介质模型

首先在应用BFO-PSO算法优化的有限差分算子在各向同性双层介质模型上进行弹性波数值模拟测试，模型示意图如图6所示，模型上下两层的纵波速度、横波速度、密度分别在图中标出，模型宽和高均为4000 m，空间步长设置为Δx=Δz=8 m，时间采样间隔为Δt=0.0005 s.震源坐标为(2000 m,1000 m)，采用主频为25 Hz的Ricker子波.

图6 双层介质模型Fig.6 Double-layer velocity model

图7和图8为分别应用BFO-PSO算法优化得到的有限差分算子和常规有限差分算子进行弹性波数值模拟所得到的合成波场快照的X分量和Z分

量，其中，图7a和图8a应用8阶常规有限差分算子，图7b和图8b应用BFO-PSO算法优化得到的8阶有限差分算子，图7c和图8c应用16阶常规有限差分算子.对比分析发现，应用8阶常规有限差分算子获得的合成波场快照，无论是X分量还是Z分量，均存在明显的数值频散现象，如图中黑色箭头指示.而应用BFO-PSO算法优化得到的有限差分算子进行弹性波数值模拟，所获得合成波场快照中几乎没有数值频散，数值模拟的精度甚至达到了16阶常规差分算子的效果.因此，采用BFO-PSO算法可以有效的提高有限差分算子的精度和效率，优化获得的低阶有限差分算子在逼近微分算子时就具有较小的数值误差，进而可以在不增加计算量，保证计算效率的同时，达到常规高阶有限差分算子的精度，满足高精度数值模拟的要求，有效提高弹性波数值模拟的精度和效率.

图7 应用不同有限差分算子得到的合成波场快照(X分量)(a) 8阶常规有限差分算子； (b) 8阶BFO-PSO算法优化的有限差分算子； (c) 16阶常规有限差分算子.Fig.7 Synthetic wavefield snapshots (X component) using different FD operators(a) Eighth-order conventional FD operator; (b) Eighth-order FD operator optimized based on the BFO-PSO algorithm; (c) Sixteenth-order conventional FD operator.

图8 应用不同有限差分算子得到的合成波场快照(Z分量)(a) 8阶常规有限差分算子； (b) 8阶BFO-PSO算法优化的有限差分算子； (c) 16阶常规有限差分算子.Fig.8 Synthetic wavefield snapshots (Z component) using different FD operators(a) Eighth-order conventional FD operator; (b) Eighth-order FD operator optimized based on the BFO-PSO algorithm; (c) Sixteenth-order conventional FD operator.

2.2 Marmousi模型

之后，我们在较为复杂的Marmousi模型上进行弹性波数值模拟，以进一步对比分析BFO-PSO算法优化得到的有限差分算子在压制数值频散、提升弹性波数值模拟效率方面的性能.如图9所示，为Marmousi纵波速度模型，在本文数值模拟实验中，模型横波速度和纵波速度的比值为vP/vS=1.7，密度设置为ρ=2000 kg·m-3.模型空间离散网格大小设置767×376，空间步长设置为Δx=12 m,Δz=8 m，时间采样间隔为Δt=0.001 s，记录时长为t=5 s.震源坐标为(4600 m,1 m)，采用主频为25 Hz的Ricker子波.

图9 Marmousi模型Fig.9 Marmousi model

图10为分别应用BFO-PSO算法优化得到的有限差分算子和常规有限差分算子在Marmousi模型上进行弹性波数值模拟所得到的合成地震记录.其中，图10a、b、c分别为应用8阶常规有限差分算子、8阶BFO-PSO算法优化得到的有限差分算子、16阶常规有限差分算子得到的合成地震记录，图10d、e分别为区块1和2的放大记录.根据合成地震记录的对比分析，可以发现，应用BFO-PSO算法优化的有限差分算子得到的地震记录具有更好的精度和效果，几乎与16阶常规有限差分算子的精度一致，而8阶常规算子得到的地震记录上则可以清晰看到数值频散存在.因此，通过在层状介质模型和复杂的Marmousi模型上进行的弹性波数值模拟测试，可以充分说明，BFO-PSO算法优化得到的有限差分算子具有较高的精度，可以有效压制用差分算子替代微分算子时导致的数值频散，优化后的低阶差分算子就可以达到常规高阶算子的精度，从而可以有效提高对弹性波动方程求取数值解时的计算效率，降低对于计算资源的要求.

图10 在Marmousi模型应用不同有限差分算子进行弹性波数值模拟得到的合成地震记录(Z分量)(a) 8阶常规有限差分算子； (b) 8阶BFO-PSO算法优化的有限差分算子； (c) 16阶常规有限差分算子； (d) 区块1放大记录； (e) 区块2放大记录.图(d)、(e)中，左边为8阶常规有限差分算子，中间为8阶BFO-PSO算法优化的有限差分算子，右边为16阶常规有限差分算子.Fig.10 Synthetic shot records (Z component) obtained by numerical modelling of elastic waves using different FD operators on Marmousi model(a) Eighth-order conventional FD operator; (b) Eighth-order FD operator optimized based on the BFO-PSO algorithm; (c) Sixteenth-order conventional FD operator; (d) Zoomed view of block 1; (e) Zoomed view of block 2. In (d) and (e), the left is eighth-order conventional FD operator; middle is eighth-order FD operator based on the BFO-PSO algorithm; and right is sixteenth-order conventional FD operator.

3 结论

本文提出一种BFO-PSO算法下的有限差分算子优化方法，并应用此优化的有限差分算子分别在层状介质模型和复杂模型上进行弹性波数值模拟，对比分析合成地震记录.首先，针对PSO算法对初始参数依赖性大及算法易于陷入局部极值的缺点，引入BFO算法的趋化、复制、驱散三种步骤，以改进PSO算法的粒子速度和位置更新策略，形成BFO-PSO混合优化算法.其次，构造包含有限差分系数的目标函数，并应用BFO-PSO混合优化算法求取最优解，获得优化的有限差分算子，根据理论频散曲线分析，BFO-PSO算法可以在快速收敛的同时获得较好的全局最优解，对于低阶差分算子，在保证精度误差的同时，有效扩大了频谱覆盖范围，满足弹性波数值模拟对高波数的要求；对于高阶差分算子，BFO-PSO算法优化得到的差分系数精度更好，即BFO-PSO算法更适用于复杂的多参数优化问题.最后，应用BFO-PSO算法优化后的有限差分算子进行弹性波数值模拟，对比分析合成波场快照及合成地震记录，进一步验证了BFO-PSO混合优化算法获得的优化有限差分算子在压制数值频散方面的有效性，降低了弹性波动方程数值求解的误差，提高数值模拟的精度和效率.