全波形反演中波场存储量的优化分析

全波形反演中波场存储量的优化分析

0波场重建策略地球物理勘探的主要目的是获得地球上的精确图像。近年来,随着计算硬件性能的提升以及GPU(GraphicProcessingUnit)等高性能计算设备在地震数据成像中的成功应用,具有高分辨能力的全波形反演(FullWaveformInversion,FWI)方法受到了地球物理学家们青睐时间域的FWI方法最先由Tarantola实现,随后,Pratt将FWI发展到频率域除了初始模型的问题会干扰FWI的应用效果,巨大的计算量和存储量同样也成为影响FWI实用化推广的棘手问题。在FWI迭代计算过程中,反演梯度的计算与逆时偏移(ReverseTimeMigration,RTM)方法流程类似,因此,采用互相关成像方法时,震源波场的存储量决定着RTM和FWI方法对计算硬件存储性能的要求。同理,针对RTM方法的存储策略同样可借鉴应用至FWI方法中。Clapp提出在速度外围构建随机速度层作为边界条件在前人研究工作的基础上,本文提出了一种基于波场重建的时域FWI方法。借助拉格朗日乘子法,首先推导了FWI的梯度表达式,并引入解析步长方法和混合共轭梯度优化方法求解该反问题。针对FWI成像梯度计算所需的巨大存储量问题,我们引入了一种能大幅降低存储成本的波场重建策略,为时域FWI的实用化推广提供一套切实可行的理论方案。最后通过匀速模型和Marmousi2模型的波场重建和FWI测试验证了本文方法的有效性。1方法原理1.1震源波场模拟时空域二阶声波常密度波动方程为:式中,x是空间坐标;x式中,x全波形反演(FWI)通过最小化模拟数据和野外观测数据之间的误差,推导地下的速度参数,其目标函数通常定义为模拟数据和观测数据残差的L2范数式中,ns和ng分别为炮点和检波点数目;p该目标函数受到模拟数据满足的波动方程的限制,我们采用拉格朗日乘子法法求解该约束性优化问题,建立拉格朗日乘子符为:式中,λ(x,t;x该伴随方程与声波方程(1)具有相同的表达形式,因此只需要将震源项替换为模拟数据和观测数据的残差,即可利用相同的数值求解方法模拟震源波场和伴随波场。因此,传统的声波全波形反演的流程可以简述为:1)沿时间正向求解方程(1),得到正向传播的震源波场p2)沿时间反方向求解方程(6),得到反向传播的检波点波场λ(x,t;x3)从计算机内存中提取保存的波场,应用公式(5),构建FWI的梯度。4)采用合适的反演算法,计算优化的步长,更新速度参数v。5)重复1~4步骤,直到目标函数值小于设定的值或满足最大迭代次数。1.2快速搜索步长方法上一节,我们已经推导了FWI的梯度,给出了传统FWI的实现流程。然而,为了得到合理的反演结果,最优化的搜索步长也尤其重要,常规的步长求取方法有很多,如:固定步长将第k+1次迭代的目标函数进行二阶泰勒展开,可得:式中,α式中,<,>表示两个矢量之间的内积;我们近似了Hessian矩阵H式中,f为式(1)所示的正演算子,ε为试探步长,它满足条件1.3嘴唇交叉试验常规最速下降方法在迭代寻优的过程,收敛速度非常缓慢,会带来锯齿折叠效应。为了在迭代过程中,求解一个合理的下降方向d其中,β在得到共轭梯度方向和优化步长后,速度参数则能够迭代地更新,更新公式可以表示为:2逆时波场重建在传统FWI的实现过程中,需要对背景震源波场和伴随波场进行零延迟互相关来构建梯度。但二者在时间方向上互逆,因此可预先保存某一方向所有时刻的波场值,当另一个波场传播时,提取相应时刻保存的波场进行互相关运算,但由此产生巨大的存储需求和I/O(Input/Output)耗时通常是难以接受的。因此,本文提出一种适用于FWI的逆时重建震源波场的方法,旨在降低了FWI方法的存储需求和I/O耗时。本文采用高阶有限差分法求解波动方程(1):波动方程差分公式如图1所示,假设真实波场范围为图1中的A区域,如果采用空间12阶差分计算格式求解波动方程,仅需保存A区域外侧6层波场数据,如图1中B区域,即可满足差分条件,C区域为吸收边界区域。为了实现震源波场的逆时重建,震源波场沿正向传播时,我们需要记录每时刻的B区域的波场信息以及最后时刻的A区域的波场信息。利用这些保存的信息,我们可以实现逆时波场重建:首先通过波动方程的逆向运算实现震源波场的逆时传播,同时,将对应时刻的B区域内的波场使用正传过程中保存的波场更新,以满足差分运算条件。在震源波场逆时重建过程中,结合公式5即可进行梯度的构建,直到所有时刻和所有炮完成叠加。其余反演步骤与传统方法一致,综上所述,本文提出的基于波场重建的时域FWI流程如图2所示。文中首先采用简单的均匀各向同性的介质进行震源波场逆时重建方法测试,该模型的水平和垂直方向网格点均为301,沿水平和垂直方向的网格尺寸均为10m。该模型的纵波传播速度为2800m/s,为了吸收来自由于边界截断产生的反射波,我们应用了完美匹配层吸收边界(PML)图3为均匀各向同性介质不同时刻的波场快照,其中图3a和图3b分别显示了正向传播和本文方法逆时重建的震源波场,从左自右的时刻分别为0.15,0.30,0.45,0.60s。我们不难发现,模型区域内的震源波场已经被很好地重建出来,由于PML边界内的波场不参与梯度的计算,因此并没有重建。我们进一步对波场快照抽取单道进行对比,正传波场和重建波场的单道数据对比如图4所示,抽取位置在t=0.45s,x=2km处。观察可知,正传的单道图4a和逆时重建的单道图4b几乎没有差别。两者的误差单道图4c所示,可知单道误差10为了证明文章中所述方法能够节约存储,我们对比了传统全波场存储的方法和文章所提出方法。假设,该模型的传播总时间间隔为5000,波场数据在GPU端的数据类型为4字节浮点型,全波场存储的计算方法为:nx·hz·nt·4byte,本文的方法的存储为z·(nx+nz)·nt·6·4byte,其中6表示为有限差分方法的空间差分阶数的1/2,nx和nz为模型沿水平和垂直方向的空间网格点数。存储消耗对比如表1所示,这验证了本文所述方法能够大幅降低存储。特别地,当模型为3D时,文中所述的方法将大幅降低存储消耗,加速反演的计算速度。3模型测试为了验证文中所述波场重建方法和全波形反演方法的有效性,分别对Marmousi2模型进行的波场重建测试和FWI方法测试。3.1震源波场重建首先,进行了基于Marmousi2模型的波场重建测试,速度模型如图5所示。模型沿水平和垂直方向的网格点数分别为500和200,沿水平和垂直方向的网格间距均为10m,为了压制了来自边界的虚假反射,我们在模型四周加载了厚度为50层PML吸收边界。震源位置位于水平方向2.5km,垂直方向0.65km处,震源函数为主频20Hz的雷克子波,时间采样间隔为1ms。分别截取正传至t=0.20,0.40,0.60,0.80s时刻的震源波场,如图6a所示,当地震波传播至最大时刻5s时,进行波场逆时重建,重建所得的t=0.20,0.40,0.60,0.80s时刻的波场如图6b所示。根据图6对比可知,本文所述的逆时重建方法能够很好地恢复出在模型区域内的震源波场。进一步对波场快照抽取单道数据进行对比,所选数据为t=0.40s,x=3km处的波场数据,对比结果如图7所示。正传波场和逆时重建波场的的单道数据分别如图7a和7b所示,二者并无肉眼可见差别,对两者做相减运算,分析所得误差如图7c所示,误差为103.2迭代次数对fwi方程力学性能的影响为了验证本文提出的FWI方法对复杂模型的有效性,对Marmousi2模型进行了FWI测试,精确速度模型如图5所示。模型的基本信息与波场重建测试相同。震源和检波器均置于地面,震源的范围在0~5km,炮间距为250m,共20炮,采用主频为10Hz的雷克子波,共500个检波器,检波点间距为10m。时间采样间隔为1ms,总时间长度为5s。由于实际地震勘探无法获得精确的速度模型,因此在进行FWI测试时,首先将图5所示的精确速度进行平滑,如图8所示,以此作为初始速度模型。不同迭代次数的反演结果如图9所示。随着迭代次数的增加,反演结果逐渐接近真实速度,对比不同迭代次数的反演结果可进一步发现,初期的反演结果对中浅部的速度构造更新效果较好,当迭代次数超过50次后,深层的速度构造也得以更新,趋近真实速度。分别提取了水平方向1,2,3,4km处的反演结果的单道数据,如图10所示,黑色实线表示初始速度,蓝色实线表示真实速度,绿色和红色实线分别表示迭代20次和100次的反演结果。曲线对比的结果与前文所述结论一致,执行20次迭代计算后,文中所述的方法可以很好地反演出中浅层的速度,随着迭代次数的增加,深层构造也得到了有效更新。图11显示了目标函数随迭代次数的下降情况,随着迭代计算,理论模拟数据与观测数据误差逐渐减少,进一步验证了本文方法的有效性。4重建策略与模型试算波场数据的存储策略对全波形反演方法的应