基于优化参数Adam算法的全波形反演

王倩倩, 宋鹏,2,3*, 华清峰, 刘保华, 李官保, 王绍文, 都国宁

1 中国海洋大学海洋地球科学学院, 青岛 266100 2 崂山实验室, 海洋矿产资源评价与探测技术功能实验室, 青岛 266100 3 中国海洋大学海底科学与探测技术教育部重点实验室, 青岛 266100 4 自然资源部第一海洋研究所, 青岛 266061

0 引言

Tarantola(1984)在20世纪80年代提出全波形反演(Full Waveform Inversion, FWI)理论,其可用于获得地下地层速度模型信息.基于最小二乘反演理论,FWI将模拟地震记录与实际地震记录之间的残差平方和作为反演的目标函数,使用伴随状态法计算目标函数的梯度,并对给定的初始速度模型进行迭代更新直至得到高分辨率的速度反演结果.由于FWI具有准确刻画地下复杂结构速度模型的潜力,该方法引起了业界的广泛关注并且成为地球物理领域的研究热点.为克服周期跳跃的影响,Bunks等(1995)提出了多尺度FWI策略,即将地震数据划分为由低到高的不同频率成分.除此之外,低波数信息对于避免局部极小值问题至关重要,因此需要从地震数据中提取低波数信息并研究多种方法来有效避免局部极小值(Wu et al., 2014; 张文生等, 2015; 陈生昌和陈国新, 2017; Lian et al., 2018; Wu and Alkhalifah, 2018; Yao et al., 2019a,b; Sun and Alkhalifah, 2019; He and Liu, 2020; Li and Alkhalifah, 2021; 李振春等, 2022).FWI已经被应用于黏声介质(Kamath et al., 2021; 邹鹏和程玖兵, 2020; Wang and Qu, 2021),弹性介质(Mora, 1987, 1988; Xu and McMechan, 2014; Ren and Liu, 2015, 2016; Qu et al., 2020),VTI介质(Zhang and Alkhalifah, 2017; Feng et al., 2021; 唐杰等, 2022),各向异性介质(张盼等, 2018; Oh and Alkhalifah, 2019),黏声各向异性介质(Qu et al., 2017)等.

经过几十年的发展,FWI在反演精度和计算效率方面都取得了巨大的进展,已有学者发展了反射波波形反演(Yao and Wu, 2017; Yao et al., 2019a,b; 胡勇等, 2021; 孙思宇等, 2021)等方法来提高反演精度.然而,常规FWI仍然存在着一个问题:目前FWI采用的大多数梯度求解算法,如共轭梯度算法(Polyak, 1969; Witte et al., 2018),L-BFGS算法(Nocedal and Wright, 2006)和截断牛顿算法(Nash, 2000)等,即使其梯度再经过预处理,在深部地层的收敛效率和精度依然不高(Alkhalifah et al., 2018).此外,这些梯度类算法在每次更新模型时还需要额外的波场计算来获得一个合适的迭代步长,进一步降低了FWI的计算效率.

以卷积神经网络(Hinton and Salakhutdinov, 2006)为代表的深度学习技术发展迅速,其为解决常规FWI问题提供了新的解决方案(李燕梅等, 2022; 王竟仪等, 2023).神经网络的优化过程本质上也是一个基于梯度优化算法参数的迭代更新过程,对于含有数百万个待优化参数的神经网络,其优化必须依靠高效的梯度优化算法.目前已有部分学者(Sun et al., 2020; Bernal-Romero and Iturrarán-Viveros, 2021)尝试将RMSProp(Hinton et al., 2012),Adam和AdaMax(Kingma and Ba, 2017)等具有椭球约束的一阶置信域(Adolphs et al., 2020)梯度优化算法引入FWI中,实验结果表明这些算法在反演精度和收敛效率方面均具有明显优势.

Adam算法为目前较为常用的深度学习梯度优化算法,其迭代步长和超参数{β1,β2}的取值对于梯度的优化效果至关重要(Ma and Yarats, 2019),然而深度学习框架所提供的默认参数并没有针对FWI进行优化.鉴于此,Sun等(2020)分析了迭代步长的取值对于梯度的优化影响,并给出了迭代步长的取值原则,Kuldeep和Shekar(2021)分析了{β1,β2}等超参数对于反演的影响,并给出了一个相对保守的取值范围来使损失曲线保持平稳.本文通过Marmousi模型数据和Overthrust模型数据的实验,系统分析了{β1,β2}等超参数的取值对于FWI的优化效果,并给出了优化的参数值.

1 基于Adam算法的时间域FWI

在密度恒定的二维各向同性介质中,时域声波方程为

(1)

FWI是以合成数据和观测数据的整体差异最小作为约束条件,基于梯度迭代方法来优化模型参数,最终重构能够表征地下介质构造的模型参数.最小二乘意义下的目标泛函可以表示为

(2)

其中,m为模型参数,在常密度声波方程中为地下介质的声波速度v(x);E(m)为FWI需要优化的目标函数,dobs和dsyn分别代表观测数据和模拟数据.

基于伴随状态法可得到FWI的梯度计算公式为

(3)

其中,下标k为第k次迭代,mk为第k次迭代的模型参数,Ns为震源总数,Nt为离散时间采样点数,ψ(x,t;xr)表示基于mk模型,在xr点处以合成地震记录与观测地震记录之差的残差地震记录为扰动的逆时波场,u(x,t;xs)表示基于mk模型,在xs点处通过波动方程数值正演模拟获得的正传波场.

对于常规的FWI,通过公式(3)计算得到的梯度通常还需要基于共轭梯度法(CG)、L-BFGS法(Liu and Nocedal, 1989; Moghaddam et al., 2013; Huang et al., 2017)、截断牛顿法(Yang et al., 2016; Matharu and Sacchi, 2019)等算法进行梯度优化,以达到更高精度的收敛效果.此外,一些梯度预处理方法,如基于波场能量的梯度预处理算法等(Zhang et al., 2011, 2012; Song et al., 2019)也常用来进一步提高梯度精度,进而提高深层部分的反演精度.

基于常规算法的FWI,其模型参数更新过程可以表示为

mk=mk-1-αkPgk,

(4)

其中,P表示梯度预处理算子,αk表示第k次迭代的迭代步长,其通常需要用线搜索策略(Gauthier et al., 1986; Virieux and Operto, 2009)来获得满足Wolfe条件的步长取值(Bernal-Romero and Iturrarán-Viveros, 2021).基于常规算法的FWI反演流程如图1所示.

图1 常规FWI计算流程Fig.1 The workflow of conventional FWI

传统FWI中线搜索和梯度预处理都需要额外的正演模拟计算,其降低了FWI的计算效率.此外,即使每次迭代中都使用梯度预处理算子进行预处理,其梯度值也并未得到足够优化,这会降低FWI的收敛效率与反演精度(Alkhalifah et al., 2018).

深度学习领域其神经网络的优化过程本质上也是一个基于梯度优化算法参数的迭代更新过程,对于含有数百万个待优化参数的神经网络,其优化必须依靠高效的梯度优化算法.Adam算法是目前较为常用的深度学习梯度优化算法,其结合了AdaGrad(Duchi et al., 2011)和RMSProp(Tieleman and Hinton, 2012)的优点,具有算法稳定、收敛效率高等优点,目前已有部分学者实现了基于Adam算法梯度优化的FWI(Sun et al., 2020; Bernal-Romero and Iturrarán-Viveros, 2021).

基于Adam算法的FWI,其迭代步长均由算法本身计算得到,因此在计算步长时不需要额外的波场延拓计算;此外,基于Adam算法的优化梯度均是由历史梯度的指数加权平均来获取,其不需要应用基于波场能量的梯度预处理等算法进行后续处理,这意味着非常可观的计算量节约.

在Adam算法中,其梯度的一阶矩估计的计算公式为

(5)

(6)

基于Adam算法的模型参数更新量的公式为

(7)

式中,Δmk为第k次迭代的模型更新量.ε为保持数值稳定的常数,其值一般为1×10-10.αk为第k次迭代的迭代步长.本文通过常规算法计算了第一次迭代步长值,并设定初始迭代步长值为α1,然后在后续的迭代中对初始步长施加一个衰减系数为0.99的指数衰减来计算其余迭代步长,如公式(8)所示.通过公式(8)计算所得的步长与Sun等(2020)所提出的步长有效范围一致:

αk=max{‖Δm0‖}/2*γk-1,

(8)

其中,‖Δm0‖为模型第一次迭代时基于常规算法计算所得的模型更新量的绝对值,γ=0.99为迭代步长的指数衰减系数.

基于Adam算法的FWI计算流程图如图2所示.与常规算法对比(如图1所示),Adam算法可以直接从梯度获得优化后的模型更新量,避免了迭代步长和梯度预处理算子P的求取,从而提高FWI的计算效率.

图2 基于Adam优化算法的FWI计算流程Fig.2 Workflow of FWI based on Adam algorithm

2 Adam算法参数测试与分析优化

Adam算法的迭代步长和超参数{β1,β2}的取值对于梯度的优化效果至关重要,当使用Adam算法的默认参数进行反演的时候,其反演精度是相近的(Bernal-Romero and Iturrarán-Viveros, 2021).由于这些默认参数是基于深度学习的任务进行设计的,因此其在FWI应用中还有改进的空间.Sun等(2020)分析了迭代步长的取值对于梯度的优化影响,并给出了迭代步长的取值原则,但目前尚未有对于超参数{β1,β2}的取值对于FWI影响的系统分析.本文基于Marmousi模型来系统测试超参数{β1,β2}的取值影响,并最终给出最优的参数取值,来最大限度的提升基于Adam算法的FWI的收敛效率和反演精度.

本文所采用的Marmousi模型的水平距离和模型深度分别为9.2 km和2.4 km,横向网格间距和纵向网格间距均为20 m,真实速度模型和初始速度模型分别如图3a、b所示.雷克子波的主频为10 Hz,其记录长度为4 s,采样间隔为1 ms.观测数据共有100炮,每炮共有461道进行接收,由时间二阶、空间十阶的有限差分进行模拟得到,边界条件使用混合吸收边界条件(Xie et al., 2020).炮点和检波点的深度为20 m,其中检波点的位置是固定的,炮点位置为从左向右进行移动.

图3 Marmousi模型(a) 真实模型; (b) 初始模型.Fig.3 Marmousi model(a) Real model; (b) Initial model.

考虑到梯度优化与数据频带无关,本文测试选取了0～5 Hz的频带范围进行测试.为系统测试Adam参数对于反演的影响,首先分别给出如表1和表2的β1和β2参数取值,其中β1由0.9均匀降至0.1共9个值,β2从0.999逐渐减小至0.1共10个值.测试分析时,以表1中β1的值为基准,共测试9组;每组中β1的取值固定,β2分别取自表2中的各个值,共10次反演测试,分别记为Adam1～Adam10,各组试验的目标函数曲线分别如图4a—i所示.

表1 各组β1的值Table 1 Values of β1 of different sets

表2 每组的β2值Table 2 Values of β2 in each set

图4 9组实验(a—i)的目标函数曲线Fig.4 The objective functions for 9 sets of experiments (a—i)

图4中各条曲线的平滑度一定程度上可反映反演的稳定性,最终目标函数值反映出反演的精度,因此相对平滑且最终值较小目标函数曲线可被认为对应于1组相对较优的参数{β1,β2}.鉴于以上分析,由于图4f—i(即β1≤0.4)中的目标函数曲线震荡较为剧烈且最终的值大于图4a—e的结果,因此更适合于FWI的Adam参数应该存在于1组到5组实验中.

在第1组到5组中选择最优的1～2个参数组合,重新绘制目标函数曲线如图5所示,其所对应的参数如表3所示.

表3 第1组至第5组的参数值Table 3 Suboptimal values of parameters from group 1 to group 5

图5中的蓝色曲线代表β1=0.9、β2=0.999,其为大多数深度学习训练过程中经常使用的参数值(Chollet, 2018),然而由图5可以看出,这组参数并不是FWI的最优参数,当β1=0.6、β2≥0.5时,其反演结果是所有测试中最优的,若想提高反演收敛的稳定性,可将β2设置为接近1的值,若想提高反演的收敛效率,可将β2设置为接近0.5的值.因此参数β1=0.6、β2≥0.5这个范围是适合于FWI的最优Adam算法参数范围.为了提高反演的收敛效率,在本文的模型实验中,将β2的值均设置为0.5.

3 模型实验

本文分别基于L-BFGS算法、默认参数Adam算法和优化参数的Adam算法,应用Marmousi模型(与第2节中使用的模型一致)和Overthrust模型进行FWI实验,反演频带范围依然为0～5 Hz.

由于基于Adam梯度优化算法的FWI可以直接获得优化后的模型,不需要计算迭代步长和基于波场能量的梯度预处理,其计算效率高于基于L-BFGS算法FWI的计算效率,按照波场延拓次数比较基于L-BFGS算法的FWI其计算成本约为基于Adam算法FWI的1.5倍.因此为对比相同计算效率下的反演收敛效率和精度,本文不比较二者迭代次数相同时的反演结果,而是比较相同波场延拓次数下的反演结果,因此基于L-BFGS算法的FWI进行了12次迭代,基于Adam算法的FWI进行了20次迭代.

通过计算,基于Adam算法FWI的初始迭代步长为45.73,且每次迭代的步长均在Sun等(2020)给出的参数范围之内.图6a、b、c分别为基于含梯度预处理的L-BFGS算法、基于默认参数Adam算法与基于最优参数Adam算法的FWI在经过6000次波场延拓计算后的反演结果.此时,基于L-BFGS算法的FWI进行了12次迭代,而基于Adam算法的FWI进行了20次迭代.由图6可知,基于Adam算法的FWI反演精度明显优于基于L-BFGS算法FWI的反演精度,尤其是深层地区(如黑色虚线圈所示).此外,与基于默认参数Adam算法的FWI相比,基于最优参数Adam算法的FWI在浅部区域具有更高的反演精度(如图6中白色虚线圈所示).

图6 不同算法的Marmousi模型反演结果(a) L-BFGS算法; (b) 默认参数Adam算法; (c) 最优参数Adam算法.Fig.6 Inversion results of Marmousi model based on(a)The L-BFGS algorithm result; (b) The Adam algorithm with default parameters; (c) The Adam algorithm with optimal parameters.

将图6中5.2 km处速度反演结果进行抽道并与真实速度曲线进行对比(如图7所示).在图7中,黑色曲线代表真实速度模型,蓝色曲线代表基于L-BFGS算法FWI的反演结果,红色曲线代表基于默认参数Adam算法FWI的反演结果,绿色曲线代表基于最优参数Adam算法FWI的反演结果.图7表明,与基于L-BFGS算法和基于默认参数Adam算法的FWI相比,基于最优参数Adam算法FWI的速度曲线更接近于真实速度,说明我们上述选择的最优参数Adam算法FWI与其他两种算法相比,具有更高的反演精度和收敛效率.

图7 5.2 km处速度纵向切片对比Fig.7 Velocity profiles at location 5.2 km

图8和图9分别显示了三种反演方法的目标函数对数迭代曲线与归一化误差迭代曲线.在图8和图9中,蓝色曲线代表基于L-BFGS算法的FWI反演结果,红色曲线代表基于默认参数Adam算法的FWI反演结果,绿色曲线代表基于最优参数Adam算法的FWI反演结果.图8和图9表明,当反演目标函数值相同时,基于最优参数Adam算法的FWI与其他两种算法相比,其收敛速度更快(波场模拟次数更少).因此,基于最优参数Adam算法的FWI在这三种算法中具有最快的收敛效率.

图8 目标函数对数迭代曲线(Simulations为波场模拟次数)Fig.8 The curves of logarithmic objective function (Simulations are the number of wave field simulations)

图9 归一化误差迭代曲线(Simulations为波场模拟次数)Fig.9 The curves of normalized data error (Simulations are the number of wave field simulations)

本文基于一张Telsa-V100卡分别对基于L-BFGS算法、基于默认参数Adam算法、基于最优参数Adam算法的FWI进行计算,基于L-BFGS算法12次迭代计算时间与基于Adam算法20次迭代计算时间如表4所示.由表4可知,基于Adam算法的FWI其迭代20次所需时间与基于L-BFGS算法的FWI迭代12次所需时间基本一致.由此可知,基于Adam算法的FWI其计算效率要高于基于L-BFGS算法的FWI,这是其不需要进行基于地震波能量预处理和计算梯度步长而节省的2次波场延拓计算时间所致.

表4 不同算法相同波场延拓次数下的FWI迭代时间Table 4 Computation time of different algorithms and same simulations

为了进一步验证基于最优参数Adam算法的收敛效率,本文还采用深度为4.0 km、宽度为12.5 km的Overthrust模型进行测试,其真实速度模型和初始速度模型分别如图10a、b所示.时间采样间隔为0.002 s,网格间距为25 m,我们通过100炮的合成地震数据,基于L-BFGS算法、默认参数Adam算法以及最优参数Adam算法进行迭代反演.通过计算,基于Adam算法FWI的初始迭代步长为56.28,且每次迭代的步长均在Sun等(2020)给出的参数范围之内.

图10 Overthrust速度模型(a) 真实模型; (b) 初始模型.Fig.10 The Overthrust velocity model(a) The true model; (b) The initial model.

本文同样比较经过6000次波场模拟计算后的反演结果,其中L-BFGS算法进行了12次迭代,Adam算法进行了20次迭代.图11a为基于梯度预处理L-BFGS算法FWI的反演结果,图11b、c分别为基于默认参数Adam算法和基于最优参数Adam算法FWI的反演结果.与基于L-BFGS算法相比,基于Adam算法的反演结果(图11b、c)具有更高的反演精度,尤其是在深部地层(如图11中的黑色虚线圈所示).

图11 Overthrust模型迭代反演结果(a) L-BFGS算法; (b) 默认参数Adam算法; (c) 最优参数Adam算法.Fig.11 Inversion results of Overthrust model(a) L-BFGS; (b) Adam with default parameters; (c) Adam with optimal parameters.

图11中4.8 km和10.1 km处反演结果纵向切片对比如图12所示.为了便于比较,将真实速度曲线用黑色表示.由图12可知,基于Adam算法的反演结果比基于L-BFGS算法的反演结果反演精度更高,此外,与基于默认参数Adam算法相比,基于最优参数Adam算法进一步提高了深部地层的反演精度.

图12 速度纵向切片对比(a) 4.8 km; (b) 10.1 km.Fig.12 Velocity profiles at locations

图13和图14为Overthrust模型的目标函数对数迭代曲线与归一化误差迭代曲线.从图中可以看出,与L-BFGS算法和默认参数Adam算法相比,优化参数Adam算法仍然具有最好的收敛效率和反演精度,由此可知,基于Marmousi模型所选择的最优参数也同样适用于其他复杂模型.

图13 目标函数对数迭代曲线(Simulations为波场模拟次数)Fig.13 The curves of logarithmic objective function (Simulations are the number of wave field simulations)

图14 归一化误差迭代曲线(Simulations为波场模拟次数)Fig.14 The curves of normalized model error (Simulations are the number of wave field simulations)

4 结论

本文在详细论述基于Adam梯度优化算法的FWI原理和实现步骤的基础上,通过大量的模型数据测试,系统分析了超参数{β1,β2}的取值对FWI的收敛性和反演精度的影响,并给出了更适合于FWI的优化参数.通过参数测试与分析优化,得到当β1=0.6、β2≥0.5时,其反演结果是所有测试中最优的,且若想提高反演收敛的稳定性,可将β2设置为接近1的值,若想提高反演的收敛效率,可将β2设置为接近0.5的值.基于Marmousi模型和Overthrust模型的数值实验表明,相比于基于L-BFGS算法和基于默认参数Adam算法的FWI,基于最优参数Adam算法的FWI其具有更高的收敛效率和反演精度.

当然,考虑到实际数据的反演还会受到模型低频信息缺失,初始模型精度不够,地下构造过于复杂(如盐丘模型和我国西部的逆冲断层模型)等众多因素的影响,在这种情况下基于优化参数的Adam算法在实际数据资料应用中不一定可以达到最优,因此将该算法引入实际数据的反演中,通过实际数据反演测试进一步提高该算法的实用性与精度将是我们下一步的工作重点.