一种以逆时偏移算法为引擎的波动方程初至走时层析方法

摘要

本发明公开了一种以逆时偏移算法为引擎的波动方程初至走时层析方法。包括：输入地震子波，观测系统，深度域初始速度场，观测地震记录的初至走时，以及走时残差范数的最小值；利用声波波动方程正演方法，计算初始速度场中正演的模拟地震记录，求解波动方程；计算模拟地震记录的初至波到达时；计算初至波走时残差，判断走时残差的范数是否小于预先设定的精度要求；计算波动方程走时敏感度矩阵；求解层析线性方程组；更新速度模型；输出反演得到的速度模型。优点是：采用声波波动方程具有更高的反演分辨率；更适用于描述波场的前向散射现象；将波动方程层析用于初至波走时层析中，得到了高精度的背景速度反演方法，能够反演大尺度的宏观背景速度场。

说明书

技术领域

本发明涉及地震层析成像技术领域的层析速度反演技术，尤其涉及一种基于 双程声波方程的初至走时层析速度反演方法。

背景技术

根据利用的数据信息不同，地震层析成像技术主要可分为走时层析，振幅层 析及波形层析。其中，走时层析方法最为稳健。按正过程不同，走时层析主要 可分为射线层析，胖射线层析，菲涅尔体层析和波动方程层析。

射线层析(Bishop等，1985；Tomographicdeterminationofvelocityand depthinlaterallyvaryingmedia,Geophysics,50,903-923)基于高频射 线理论，对初始模型精度的要求较低，但反演精度不高。同时，由于层析矩阵 十分稀疏，矩阵的零空间比较大，因而求解时收敛较慢，且反演结果受射线照 明的影响严重(Hu和Marcinkovich，2012；Asensitivitycontrollable target-orientedtomographyalgorithm,82ndSEGAnnualMeeting,Expanded Abstracts,1-5)。

考虑到地震波的传播不仅受到射线路径上的速度的影响，也受到射线周围速 度结构的影响。胖射线层析(Vasco等，1995；Beyondraytomography:Wavepaths andFresnelvolumes,Geophysics,60(6),1790-1804)将射线加宽，其物理 背景是波传播是有一定宽度的，但没有严格的数学物理推导。此外，射线宽度 的选择也有技巧，如固定宽度(Michelena和Harris，1991；MichelenaR.J., andJ.M.Harris,1991,Tomographictraveltimeinversionusingnatural pixels)或随着射线传播距离增大而逐渐增大射线宽度(Xu等，2006；Enhanced tomographyresolutionbyafatraytechnique.76thSEGAnnualMeeting, ExpandedAbstracts,3354-3358)。但这些做法不能定量描述波动现象。

相比而言，菲涅尔带(体)层析(Yomogida，1992；Fresnelzoneinversion forlateralheterogeneitiesintheearth,PureandAppliedGeophysics, 138(3),391-406)在改善射线层析矩阵稀疏性的同时，还考虑了地震波在第一 菲涅尔带(体)内的有限频效应。但单频菲涅尔带只有在常速介质中存在解析 表达式(Cerveny和Soares，1992；Fresnelvolumeraytracing,Geophysics, 57,902-915)。在变速介质中，可以用常速介质中主频对应的菲涅尔带来近似 有限频带地震波的菲涅尔带宽度(刘玉柱等，2009；初至波菲涅耳体地震层析 成像，2009,52(9),2310-2320)；或计算炮点和检波点的走时场，利用单频菲 涅尔带的公式计算菲涅尔带的范围(Watanabe等，1999；Seismictraveltime tomographyusingFresnelvolumeapproach:69thAnnualInternational Meeting,SEG,ExpandedAbstracts,1402–1405)。但这些做法只能近似描述 变速介质中的带限地震波的菲涅尔带(体)范围。

波动方程走时层析没有引入高频近似假设，在理论上比射线类的层析具有更 高的反演分辨率。在一定条件下(波动方程线性化近似，如一阶Born/Rytov近 似等)能够描述波场(或相位)扰动与模型扰动的定量关系。Woodward(1992； Wave-equationtomography:Geophysics,57,15-26)引入了一阶Born及Rytov 近似下“波路径”的概念(即波形与相位敏感度核函数)。Luo和Schuster(1991； Wave-equationtraveltimeinversion.Geophysics,56,645–653)利用互 相关函数测量观测数据与模拟数据的时差，并给出了互相关准则下走时扰动对 模型扰动的Fréchet导数(基于一阶Born近似)。Marquering等(1999； Three-dimensionalsensitivitykernelsforfinite-frequencytraveltimes: thebanana-doughnutparadox.Geophys.J.Int.,137,805-815)则给出了 扰动场与走时扰动的关系，进而可以导出一阶Born近似下走时敏感度核函数的 表达。

考虑到一阶Rytov近似更适用于描述波场的前向散射现象，且适用条件较 Born近似更为广泛。因此有必要发展一套在一阶Rytov近似下以波动方程为传 播算子的层析速度反演技术。

发明内容

本发明的目的是针对上述现有技术存在的不足，提出了一种以逆时偏移算 法为引擎的波动方程初至走时层析方法。该方法是一种高精度的背景速度反演 方法，能够反演大尺度的宏观背景速度场。此外，由于正问题采用声波波动方 程，没有引入高频近似，在理论上比射线类的层析具有更高的反演分辨率，因 而可以定量的描述有限频带地震波的传播效应，对小尺度的速度异常有更好的 反演能力。

本发明从波动方程线性化出发，并考虑地震波场的前向散射近似(一阶 Rytov近似)，导出了基于波动方程的走时敏感度核函数在时空域的显式计算方 法，其具体算法与著名的逆时偏移方法相似。而走时敏感度核函数就是层析线 性方程组中的稀疏矩阵，通过求解层析线性方程组，进行速度更新，可以实现 波动方程初至走时层析方法。

本发明的核心在于基于波动方程的走时敏感度核函数的构造，技术方案包 括以下步骤：

步骤1：输入地震子波，观测系统，深度域初始速度场，观测地震记录的初 至走时，以及走时残差范数的最小值；

步骤2：利用声波波动方程正演方法，计算初始速度场中正演的模拟地震记 录，求解如下波动方程；

步骤3：计算模拟地震记录的初至波到达时；

步骤4：计算初至波走时残差，判断走时残差的范数是否小于预先设定的精 度要求，若满足要求，跳转至步骤8；否则，执行下一步；

步骤5：计算波动方程走时敏感度矩阵；

步骤6：求解层析线性方程组；

步骤7：更新速度模型，进行下一次迭代反演，跳转回步骤2；

步骤8：输出反演得到的速度模型。

上述技术方案进一步优化为：

步骤1：输入地震子波f(t)，观测系统，深度域初始速度场v(x)，观测地震记 录的初至走时T^obs(x_r,x_s)，走时残差L2范数的最小值ε。

步骤2：利用声波波动方程正演方法，计算初始速度场中正演的模拟地震记 录u^cal(x_r,t|x_s)，求解如下波动方程：

$\frac{1}{v^2\left(x\right)}\frac{{\partial }^2{u}^{\mathrm{cal}}(x,t|x_s)}{\partial t^2}-\Delta u^{\mathrm{cal}}(x,t|x_s)=f\left(t\right)\delta (x-x_s)---\left(7\right)$

其中，x_r和x_s分别代表观测系统中检波器和震源的坐标，t代表时间，Δ代 表拉普拉斯算子。

步骤3：计算模拟地震记录的初至波到达时T^cal(x_r,x_s)。

步骤4：计算初至波走时残差ΔT(x_r,x_s)，按照如下公式：

ΔT(x_r,x_s)＝T^obs(x_r,x_s)-T^cal(x_r,x_s)。(8)

其中，T^obs(x_r,x_s)代表输入的观测数据的初至波到达时。

同时，判断走时残差的L2范数是否小于预先设定的精度要求(即)。

若满足要求，跳转至步骤8。否则，执行下一步。

步骤5：计算波动方程走时敏感度矩阵K(x|x_r,x_s)，按照如下公式：

$K\left(x\right|x_r,x_s)=\frac{-2}{v^3\left(x\right)}\frac{{\int \stackrel{·}{u}}_0\left(x\right|x_s;t)p\left(x\right|x_r;t)\mathrm{dt}}{\int u_0^2(x_r,t|x_s)\mathrm{dt}}.---\left(9\right)$

其中，u₀(x|x_s；t)为震源在x_s处的地震子波产生的地震波场，x为地下空间 任意一点的坐标。上标·代表时间导数。p(x|x_r；t)为检波器处逆时传播的地震波 场，表示为如下公式

p(x|x_r；t)＝g₀(x,-t|x_r,0)*u₀(x_r|x_s；t)。(10)

其中，g₀(x,-t|x_r,0)为声波波动方程的格林函数，代表在震源x_r处，从时间 t＝0开始的逆时传播过程，u₀(x_r|x_s；t)为正演波场u₀(x|x_s；t)在检波器位置x_r处 的地震记录，即

$u_0\left(x_r\right|x_s;t)=u_0\left(x\right|x_s;t){|}_{x=x_r}.---\left(11\right)$

显然，波动方程走时敏感度矩阵K(x|x_r,x_s)的计算与业内著名的波动方程逆 时偏移有相类似的计算流程，即：震源处正演的模拟波场与检波器处逆时传播 的地震记录对时间的积分。不同之处在于，公式(9)中需要的是正演的模拟波场 的时间导数，并且需要能量校正因子因此，公式(9)在已 有的逆时偏移算法基础之上，可以很容易实现。

步骤6：求解层析线性方程组:

Kδv＝ΔT(12)

其中，δv_n×1为速度更新量，矩阵K_m×n为公式(9)中计算的波动方程走时敏感 度矩阵，ΔT_m×1为公式(8)中计算的初至波走时残差。m为观测系统中检波器的数 目，n为速度模型网格离散化后的网格数目。

步骤7：更新速度模型，进行下一次迭代反演。跳转步骤 2。

步骤8：输出反演得到的速度模型，算法结束。

与现有技术相比，本发明的有益效果是：

传统的射线层析基于高频射线理论，反演精度不高。同时，由于层析矩阵 十分稀疏，矩阵的零空间比较大，因而求解时收敛较慢，且反演结果受射线照 明的影响严重。本发明将层析正问题用波动方程描述，没有引入高频近似假设， 在理论上比射线类的层析具有更高的反演分辨率。

同时，考虑到一阶Rytov近似更适用于描述波场的前向散射现象，且适用 条件较Born近似更为广泛。本发明给出了基于一阶Rytov近似的波动方程正问 题的表达，即利用业内著名的逆时偏移快速算法构造波动方程走时敏感度核函 数。因此，该方法比基于Born近似的传统的有限频层析方法的适用范围更广。

上述方法构成了基于一阶Rytov近似的波动方程走时层析的理论基础，本 发明将该理论具体应用于初至波走时层析方法中。

综上，与现有技术相比，本发明的一种以逆时偏移算法为引擎的波动方程 初至走时层析方法的所具有的优点是：

(1)波动方程层析正问题采用声波波动方程，没有引入高频近似假设，在 理论上比射线类的层析具有更高的反演分辨率；

(2)对波动方程进行一阶Rytov近似，更适用于描述波场的前向散射现象， 因而比基于Born近似的传统的有限频层析方法的适用范围更广；

(3)给出了基于逆时偏移算法的走时敏感度核函数的显式计算方法；

(4)将波动方程层析理论具体用于初至波走时层析中，得到了一种高精度 的背景速度反演方法，能够反演大尺度的宏观背景速度场。

附图说明

图1为本发明的波动方程初至走时层析方法的流程图。

图2为真实的速度模型。

图3为经过波动方程初至波走时反演得到的速度模型。

图4初至波走时对比。背景为观测的炮集。“圆点”，“加号”和“星号”分 别是真实模型、初始模型和反演得到的速度模型中的初至走时曲线。

具体实施方式

下面结合附图1对本发明的方案做进一步阐述。

步骤1：输入地震子波f(t)，观测系统，深度域初始速度场v(x)，观测地震 记录的初至走时T^obs(x_r,x_s)，走时残差L2范数的最小值ε；

步骤2：利用(声波)波动方程正演方法，计算初始速度场中正演的模拟地 震记录u^cal(x_r,t|x_s)，求解如下波动方程：

$\frac{1}{v^2\left(x\right)}\frac{{\partial }^2{u}^{\mathrm{cal}}(x,t|x_s)}{\partial t^2}-\Delta u^{\mathrm{cal}}(x,t|x_s)=f\left(t\right)\delta (x-x_s)---\left(13\right)$

其中，x_r和x_s分别代表观测系统中检波器和震源的坐标，t代表时间，Δ代 表拉普拉斯算子。

步骤3：自动计算模拟地震记录的初至波到达时T^cal(x_r,x_s)；

步骤4：计算初至波走时残差ΔT(x_r,x_s)，按照如下公式：

ΔT(x_r,x_s)＝T^obs(x_r,x_s)-T^cal(x_r,x_s)。(14) 其中，T^obs(x_r,x_s)代表输入的观测数据的初至波到达时。

同时，判断走时残差的L2范数是否小于预先设定的精度要求(即)。

若满足要求，跳转至步骤8。否则，执行下一步。

步骤5：计算波动方程走时敏感度矩阵K(x|x_r,x_s)，按照如下公式：

$K\left(x\right|x_r,x_s)=\frac{-2}{v^3\left(x\right)}\frac{{\int \stackrel{·}{u}}_0\left(x\right|x_s;t)p\left(x\right|x_r;t)\mathrm{dt}}{\int u_0^2(x_r,t|x_s)\mathrm{dt}}.---\left(15\right)$

其中，u₀(x|x_s；t)为震源在x_s处的地震子波产生的地震波场，x为地下空间 任意一点的坐标。上标·代表时间导数。p(x|x_r；t)为检波器处逆时传播的地震波 场，表示为如下公式

p(x|x_r；t)＝g₀(x,-t|x_r,0)*u₀(x_r|x_s；t)。(16)

其中，g₀(x,-t|x_r,0)为声波波动方程的格林函数，代表在震源x_r处，从时间 t＝0开始的逆时传播过程，u₀(x_r|x_s；t)为正演波场u₀(x|x_s；t)在检波器位置x_r处 的地震记录，即

$u_0\left(x_r\right|x_s;t)=u_0\left(x\right|x_s;t){|}_{x=x_r}.---\left(17\right)$

步骤6：求解层析线性方程组:

Kδv＝ΔT(18)

其中，δv_n×1为速度更新量，矩阵K_m×n为公式(9)中计算的波动方程走时敏感 度矩阵，ΔT_m×1为公式(14)中计算的初至波走时残差。m为观测系统中检波器的 数目，n为速度模型网格离散化后的网格数目。

步骤7：更新速度模型，进行下一次迭代反演。跳转步骤 2。

步骤8：输出反演得到的速度模型。

结合附图2-4，更加证明了本发明的技术效果。