CN110764146A - 一种基于声波算子的空间互相关弹性波反射波形反演方法 - Google Patents

Abstract

本发明提出了一种基于声波算子的空间互相关弹性波反射波形反演方法，采用空间相关函数实现对多波多分量反射地震数据运动学信息的有效利用，降低在数据匹配过程中对完整地震波形的高度依赖，缓解反演的非线性问题，从而准确、高效地反演纵、横波速度的背景模型。在实现过程中，本发明利用声波算子分别构建反射PP和PS的波路径，进行纵、横波速度的分步反演，这种分布反演策略可以有效压制纵、横波速度之间的耦合效应，提高两种参数的反演精度，同时显著降低算法对计算量和内存的需求。本发明得到的纵、横波速度反演结果可以为高分辨率、高精度的地震成像提供了准确的背景模型，因此对解决实际多分量地震数据的处理需求具有重要意义。

Description

一种基于声波算子的空间互相关弹性波反射波形反演方法

技术领域

本发明涉及一种弹性波多参数背景建模方法，属于地震波成像技术领域。

背景技术

弹性波全波形反演方法(EFWI)可以利用观测的多分量地震数据获取地下弹性参数模型，如纵、横波速度、密度、弹性阻抗、各向异性参数等，是探索地球内部构造、进行油气资源勘探开发等重要手段。然而，受实际观测孔径有限的影响，EFWI通常只能从观测的初至波中获取浅层速度模型的有效更新，而有关中深层模型信息则体现在反射波中，却很难利用。为此，有学者提出了弹性波反射波形反演方法(ERWI)。该方法通过提取较低波数的反射波路径，可以实现多种参数模型的中深层更新，是目前研究的热点。

然而，基于波形残差目标函数的ERWI仍然具有高度的反演非线性，导致在匹配反射波波形时出现周波跳跃现象，随着偏移距的增大，这种现象更为明显，从而产生错误的反演结果。另外，反射波形的精确模拟依赖于真振幅成像技术。事实上，这种保真成像本身就是目前成像问题中的一个巨大挑战。相较于利用反射数据的波形信息，反射数据的运动学信息更为稳健。其中，空间相关函数可以有效地衡量反射波的空间聚焦程度，可以很好地描述反射波的运动学特征，因此可以有效降低反演的非线性，同时规避对真振幅偏移的依赖。另外，弹性波场通常在时间序列上具有较长的延续，采用空间相关函数可以提高运动学信息提取的准确性，因而在复杂介质参数反演中具有更好的适应性。

这种空间相关目标函数目前仅在单分量声波RWI问题中取得了较好的实际应用案例，尚未应用于多分量多参数的弹性波参数反演问题中。值得注意的是，相对于单参数声波RWI反演问题，ERWI还需要考虑更为严重的参数耦合效应，同时负担指数增长的内存需求、计算代价，因此EFWI中常用的Hessian预条件、波场分解等技术并不适用于ERWI问题。可见，发展基于反射波运动学信息的弹性波反射波形反演方法，同时设计更为合理的多参数反演策略，对实现弹性参数背景建模至关重要。

发明内容

针对现有ERWI方法存在的周波跳跃问题和对真振幅偏移技术的依赖等问题，本发明引入空间相关函数，以充分利用地震反射数据的运动学信息，缓解周波跳跃问题和对真振幅偏移的过度依赖，获取纵横波速度的背景模型。同时，利用采用分步反演策略，基于声波算子分别构建反射PP和PS波路径，实现纵横波速度的有效解耦，提高参数反演精度，降低算法对计算量和内存的需求。

本发明的具体技术方案如下：

一种基于声波算子的空间互相关弹性波反射波形反演方法，该方法具体包括以下步骤：

(1)读取的纵、横波速度初始模型α₀和β₀以及地震观测的PP波和PS波数据d_p,obs,d_s,obs；

(2)利用声波算子基于在纵、横波速度的当前模型α_k和β_k下计算炮端正传P波场u_p,b；同时，反传多分量地震观测数据d_p,obs,d_s,obs，计算检端反传P波场w_p,b和S波场w_s,b；其中本申请涉及到的所有正、反传波场均是时间t的函数；

(3)利用步骤2中得到的炮、检端波场，获取PP波的成像结果I_piI_pp＝∫u_p,b·ψ_p,bdt，此处u_p,b和ψ_i,b均是时间t的函数；

同时，基于炮端正传P波场u_p,b和PP波的成像结果I_pp，分别获取炮端反偏移P波场u_p,r和反偏移S波场u_s,r，并在检波点处抽取合成PP波和PS波数据d_p,cal,d_s,cal；

(4)基于地震观测数据和合成数据以及空间相关函数，求取空间相关道集，其计算公式为C_p,h(t,λ)＝∫d_p,obs(t,h+λ)·d_p,cal(t,h)dh(此处均以PP波为例，PS波也同样计算)，其中，变量t代表时间序列，变量h代表观测道号，变量λ代表观测道的偏移量，h取值范围一般设置为[-50,+50]；

根据上述公式得到PP波的二维空间相关剖面C_p,h(t,λ)；对应t中每一个时刻，记录C_p,h(t,λ)取最大值时的空间偏移量λ，生成所有时刻的Δ_p(t)，并求取t中每一时刻上，检端反传的震源函数，

获得单炮情况下PP波的伴随震源g_p；

同理，计算PS波的相关道集C_s,h(t,λ)，最大相关性对应的偏移量Δ_s(t)，以及PS波的伴随震源g_s；

然后，根据空间相关目标泛函计算目标函数值，其中，s代表观测数据的炮信息；根据目标函数值，判断当前模型下反射波的空间聚焦性是否满足收敛准则，若不满足，进入步骤(5)，迭代更新纵波速度和横波速度的当前模型，直至目标函数满足收敛准则，进入步骤(9)；

(5)计算地下炮端正传波场u_p,b及其反偏移波场u_p,r，u_s,r，并根据步骤(4)中获得的PP 波和PS波的伴随震源g_p和g_s分别计算检端反传波场ψ_p,r，ψ_s,b及其反偏移波场ψ_p,r；

(6)求取炮、检端波场的零延迟互相关，并沿时间积分，即构建单炮情况下反射PP波和 PS波的波路径，

Path_pp＝∫(u_p,b·ψ_p,r+u_p,r·ψ_p,b)dt

Path_ps＝∫(u_s,r·ψ_s,b)dt

(7)将多炮反射PP、PS波的波路径进行累加，计算纵波速度和横波速度的梯度方向g_α和 g_β：

(8)根据纵波和横波速度的梯度利用给定的步长，对纵波和横波速度的当前模型α_k和β_k进行更新，

α_k+1＝α_k+γ_αg_α

β_k+1＝β_k+γ_βg_β

最后，返回步骤(2)进行下一次迭代，k＝k+1；

(9)输出纵、横波速度的当前模型α_k和β_k，作为最终的反演结果。

优选地，步骤(2)中，计算炮端正传P波场u_p,b公式如下：

式中，

代表二维空间偏导数，代表时间偏导数，f为地震子波；

计算检端反传P波场w_p,b和S波场w_s,b公式如下：

其中，dp代表P波的数据残差，ds代表S波的数据残差，符号*代表共轭算子。由于本申请采用空间相关值作为衡量标准，因而对地震数据的振幅和相位差具有较高的容忍度，因此使用声波算子实现P波和S波的快速计算，节省反射波形反演中对内存和计算量的过度需求。

优选地，步骤(3)中，获取炮端反偏移P波场u_p,r和反偏移S波场u_s,r的公式如下：

优选地，步骤(5)中，

a)对于纵波速度的当前模型更新，计算炮端正传P波场u_p,b的公式为：

计算炮端反偏移波P波场u_p,r的公式为：

计算检端反传P波场ψ_p,b的公式为：

检端反偏移P波场ψ_p,r的公式为：

b)对于横波速度的当前模型更新，需要计算炮端反偏移S波场u_s,r的公式为：

和计算检端反传S波场ψ_s,b的公式为：

由于PS像通常具有更多的偏移假象(耦合脚印和非物理波现象)，同时需要解决极性反转问题，因此，在计算反偏移的S波场时，本申请仍然采用PP波的像进行替代。

本发明首先利用观测的多分量地震数据和纵横波速度的初始模型获取地下界面信息；然后，基于成像结果利用反偏移技术预测反射数据，并根据空间相关函数估计运动学差异，求取反射PP波和PS波的伴随震源；最后，利用构建的反射PP和PS波路径合成梯度，更新参数模型。为了尽可能地压制纵横波速度之间的耦合效应，本发明需要优先反演纵波速度模型，当纵波速度较为准确后，再进行横波速度背景反演。另外，由于空间相关目标函数对地震数据的振幅、相位精度具有较高的容忍度，因此，本发明在计算弹性波场时均采用声波算子，已实现多分量地震数据的快速计算，以显著降低内存和计算量，同时便于分步反演策略的有序开展。

本发明相比现有技术具有如下有益效果：

本发明提的基于声波算子的空间互相关弹性波反射波形反演方法，采用空间相关函数突出对多分量反射数据运动学信息的有效利用，有效降低反演的非线性问题，从而准确、稳健地反演纵横波速度的背景模型。

在实现过程中，本发明采用声波算子分别构建反射PP和PS的波路径，实现纵、横波速度的分步反演，这种分布反演策略可以有效压制纵横波速度之间的耦合效应，提高参数的反演精度，同时显著降低算法对计算量和内存的需求。

本发明得到的反演结果可以为高分辨率地震成像和高精度的全波形反演提供了必要的参数建模，因此对实现参数高精度建模的实际应用具有重要意义。

附图说明

图1是本发明流程图(也作摘要附图)；

图2是Sigbee2A实验中真实纵(a)、横波(b)速度模型；

图3是Sigbee2A实验中初始纵(a)、横波(b)速度模型；

图4是地震观测反射数据水平(a)和垂直(b)分量；

图5是初始模型下得到的PP波(a)和PS波(b)成像结果；

图6是真实模型下得到的PP波(a)和PS波(b)成像结果；

图7是初始模型得到的像生成的反射PP波(a)和PS波(b)数据；

图8是初始模型下得到的反射PP波(a)和PS波(b)空间相关道集；

图9是单炮情况下的反射PP波炮(a)、检(b)端反射波路径；

图10是反演得到的纵(a)、横波(b)速度模型；

图11是反演得到的PP波(a)和PS波(b)空间相关道集；

图12是反演结果作为偏移速度得到的PP波(a)和PS波(b)成像结果；

图13是反演结果作为偏移速度得到的PP波(a)和PS波(b)空间相关道集；

图14是反演结果作为EFWI初始模型得到的纵(a)、横波(b)速度反演结果；

图15是水平位置2km处纵(a)、横波(b)速度模型的垂向抽线结果；

图16是初始模型作为EFWI初始模型得到的纵(a)、横波(b)速度反演结果。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步说明，

实施例：

本发明基于声波算子的空间互相关弹性波反射波形反演方法，下面以Sigbee2A模型实验为例，其真实模型如图2所示，该方法具体包括以下步骤：

(1)读取的纵、横波速度初始模型α₀和β₀，如图3所示，以及多分量地震观测，包括反射PP波和PS波d_p,obs,d_s,obs，如图4所示；

(2)根据声波算子以及纵、横波速度的当前模型α_k和β_k，计算炮端正传P波场u_p,b；根据观测的反射PP波数据d_p,obs，计算检端反传P波场w_p,b。本申请涉及到的所有正、反传波场均是时间t的函数；

(3)根据炮端正传P波场u_p,b和检端反传P波场w_p,b，以及成像公式I_pp＝∫u_p,b·w_{p, b}dt，获取PP波成像结果，如图5a所示，为了同时反应纵横波速度模型的准确性，图5也展示了反射PS波的成像结果，但在算法中仍然采用PP波成像结果来合成反射PS波；

根据真实模型得到的PP波成像结果，如图6所示。对比可以看出，在PP波的成像结果中，由于初始模型距离真实模型较远，主要层位定位不准确，同时底部散射点无法收敛。

根据炮端正传P波场u_p,b和PP波的成像结果I_pp，计算炮端反偏移P波场u_p,r和反偏移S波场 u_s,r，并在检波点处抽取反偏移得到的合成反射PP波和PS波数据d_p,cal,d_s,cal，如图7所示；由于反射PS波数据在零偏位置处存在极性转换，这是声波算子无法正确模拟的，因此在实际匹配之前需要对单炮PS数据在近偏移距进行置零处理)；

(4)根据地震观测P波数据和合成P波数据以及空间相关函数，其计算公式为 C_p,h(t,λ)＝∫d_p,obs(t,h+λ)·d_p,cal(t,h)dh，计算反射PP波空间相关道集，如图8所示，其中，变量t代表时间序列，变量h代表观测道号，变量λ代表空间偏移量，λ取值一般为[-50,+50] 内的所有整数；

根据上述公式得到PP波的二维空间相关剖面C_p,h(t,λ)；对应t中每一个时刻，记录C_p,h(t,λ)取最大值时的空间偏移量λ，生成所有时刻的Δ_p(t)，并求取t中每一时刻上，检端反传的震源函数，

获得单炮情况下PP波的伴随震源g_p；

根据地震观测S波数据和合成S波数据，计算反射PS波的相关道集C_s,h(t,λ)，最大相关性对应的偏移量Δ_s(t)，以及PS波的伴随震源g_s；

根据空间相关目标泛函，其公式如下

计算纵横波速度当前模型下的目标函数值，其中，shot代表观测数据的炮号；

根据目标函数值，判断当前模型下反射波的空间聚焦性是否满足收敛准则，如：是否小于初始目标函数的1％，若不满足，进入步骤(5)；若满足，进入步骤(9)；

(5)根据纵横波速度当前模型，计算炮端正传波场u_p,b及其反偏移波场u_p,r，u_s,r，并根据反射PP波和PS波的伴随震源g_p和g_s计算检端反传波场ψ_p,r，ψ_s,b及其反偏移波场ψ_p,r；

(6)根据炮端正传波场u_p,b及其反偏移波场u_p,r，检端反传波场ψ_p,r，ψ_s,b及其反偏移波场ψ_p,r，计算单炮情况下反射PP波和PS波的波路径，如图9所示，其公式如下

Path_pp＝∫(u_p,b·ψ_p,r+u_p,r·ψ_p,b)dt

Path_ps＝∫(u_s,r·ψ_s,b)dt

(7)根据单炮情况下反射PP波和PS波的波路径，计算计算纵波速度和横波速度的梯度g_α和g_β

(8)根据纵波和横波速度的梯度利用给定的步长，对纵波和横波速度的当前模型α_k和β_k进行更新，

α_k+1＝α_k+γ_αg_α

β_k+1＝β_k+γ_βg_β

最后，返回步骤(2)进行下一次迭代，k＝k+1；

(9)输出纵、横波速度的当前模型α_k和β_k，作为最终的反演结果，如图10所示。

步骤(2)中，计算炮端正传P波场u_p,b公式如下：

式中，

代表二维空间偏导数，

代表时间偏导数，f为地震子波；

计算检端反传P波场w_p,b和S波场w_s,b公式如下：

其中，dp代表P波的数据残差，ds代表S波的数据残差，符号*代表共轭算子。需要注意的是，由于本申请采用空间相关值作为衡量标准，因而对地震数据的振幅和相位差具有较高的容忍度，因此使用声波算子实现P波和S波的快速计算，节省反射波形反演中对内存和计算量的过度需求。

步骤(3)中，计算炮端反偏移P波场u_p,r和反偏移S波场u_s,r的公式如下：

步骤(5)中，

a)计算炮端正传P波场u_p,b的公式为：

计算炮端反偏移波P波场u_p,r的公式为：

计算检端反传P波场ψ_p,b的公式为：

检端反偏移P波场ψ_p,r的公式为：

b)计算炮端反偏移S波场u_s,r的公式为：

和计算检端反传S波场ψ_s,b的公式为：

从空间相关道集(如图11所示)可以看出，经过SC-ERWI方法反演后得到的空间聚焦性得到改善，偏离中心位置的程度明显降低。同时，基于反演得到的纵横波速度模型，分别求取PP波和PS波的像(如图12所示)，并提取共成像点道集(如图13所示)。可见，经过多次迭代后，界面位置更加准确，底部散射点收敛性更好，且与真实模型得到的像非常接近，这说明了利用SC-ERWI方法反演得到的纵横波速度模型已经可以较好地描述参数的背景，使得预测的反射波可以很好的匹配观测数据的运动学特征。

同时，可以将本发明反演得到的纵横波速度模型作为高精度EFWI的初始模型，得到的反演结果如图14所示，其垂直抽线如图15所示。从图中可以看出，本发明所研究的方法可以准确地恢复纵横波速度的宏观背景，再经过EFWI反演得到的纵横波速度模型则非常准确，可以很好地拟合真实模型，而如果以图2模型作为EFWI的初始模型，EFWI则会陷入局部极值，得到错误的反演结果(如图16所示)。可见，本发明方法，可以正确利用反射数据的运动学信息，恢复纵横波速度的背景模型，可以满足EFWI对初始模型精度的需求。

最后，需要指出的是，由于利用声波算子构建反射波路径来替代常规ERWI中使用的标准弹性波算子，因此，在内存上，本发明研究方法只需要存储压力分量，而ERWI至少需要存储水平和位移压力分量，因此在二维情况下内存上可以节省50％(三维情况下为66.7％)。而在计算量方面，求解两次声波方程和弹性波方程并实现波场分解所需的计算量比约为4：12，因此SC-ERWI方法在计算量上可以节省67％(如表1所示)。

表1 SC-ERWI和ERWI方法资源消耗对比

Claims (4)

1.一种基于声波算子的空间互相关弹性波反射波形反演方法，该方法具体包括以下步骤：

(1)读取的纵、横波速度初始模型α₀和β₀以及多分量地震观测，包括反射PP波数据和PS波数据d_p,obs,d_s,obs；

(2)根据声波算子以及纵、横波速度的当前模型α_k和β_k，计算炮端正传P波场u_p,b；根据观测的反射PP波数据d_p,obs，计算检端反传P波场w_p,b，其中正、反传波场均是时间t的函数；

(3)根据炮端正传P波场u_p,b和检端反传P波场w_p,b，以及成像公式I_pp＝∫u_p,b·w_p,bdt，获取PP波的成像结果；

根据炮端正传P波场u_p,b和PP波的成像结果I_pp，计算炮端反偏移P波场u_p,r和反偏移S波场u_s,r，并在检波点处抽取反偏移得到的合成反射PP波和PS波数据d_p,cal,d_s,cal；

根据地震观测PP波数据和合成PP波数据以及空间相关函数，其计算公式为C_p,h(t,λ)＝∫d_p,obs(t,h+λ)·d_p,cal(t,h)dh，计算反射PP波的空间相关道集，其中，变量t代表时间序列，变量h代表观测道号，变量λ代表空间偏移量，λ取值一般为[-50,+50]内的所有整数；

根据上述公式得到PP波的二维空间相关剖面C_p,h(t,λ)；对应t中每一个时刻，记录C_p,h(t,λ)取最大值时的空间偏移量λ，生成所有时刻的Δ_p(t)，并求取t中每一时刻上，检端反传的震源函数，

获得单炮情况下PP波的伴随震源g_p；

(4)根据地震观测PS波数据和合成PS波数据，计算反射PS波的相关道集C_s,h(t,λ)，最大相关性对应的偏移量Δ_s(t)，以及PS波的伴随震源g_s；

根据空间相关目标泛函，其公式如下计算纵横波速度当前模型下的目标函数值，其中，shot代表观测数据的炮号；

根据目标函数值，判断当前模型下反射波的空间聚焦性是否满足收敛准则，如：是否小于初始目标函数的1％，不满足，则进入步骤(5)；若满足条件，进入步骤(9)；

(5)根据纵横波速度当前模型，计算炮端正传波场u_p,b及其反偏移波场u_p,r，u_s,r，并根据反射PP波和PS波的伴随震源g_p和g_s计算检端反传波场ψ_p,r，ψ_s,b及其反偏移波场ψ_p,r；

(6)根据炮端正传波场u_p,b及其反偏移波场u_p,r，检端反传波场ψ_p,r，ψ_s,b及其反偏移波场ψ_p,r，计算单炮情况下反射PP波和PS波的波路径，其公式如下

Path_pp＝∫(u_p,b·ψ_p,r+u_p,r·ψ_p,b)dt

Path_ps＝∫(u_s,r·ψ_s,b)dt

(7)根据单炮情况下反射PP波和PS波的波路径，计算计算纵波速度和横波速度的梯度g_α和g_β：

(8)根据纵波和横波速度的梯度利用给定的步长，对纵波和横波速度的当前模型α_k和β_k进行更新，

α_k+1＝α_k+γ_αg_α

β_k+1＝β_k+γ_βg_β

最后，返回步骤(2)进行下一次迭代，k＝k+1；

(9)输出纵、横波速度的当前模型α_k和β_k，作为最终的反演结果。

2.根据权利要求1所述基于声波算子的空间互相关弹性波反射波形反演方法，步骤(2)中，计算炮端正传P波场u_p,b公式如下：

式中，

代表二维空间偏导数，

代表时间偏导数，f为地震子波；

计算检端反传P波场w_p,b和S波场w_s,b公式如下：

其中，dp代表反射PP波数据，ds代表反射PS波数据，符号*代表共轭算子。

3.根据权利要求1或2所述基于声波算子的空间互相关弹性波反射波形反演方法，步骤(3)中，计算炮端反偏移P波场u_p,r和反偏移S波场u_s,r的公式如下：

4.根据权利要求3所述基于声波算子的空间互相关弹性波反射波形反演方法，步骤(5)中，

a)计算炮端正传P波场u_p,b的公式为：

计算炮端反偏移波P波场u_p,r的公式为：

计算检端反传P波场ψ_p,b的公式为：

检端反偏移P波场ψ_p,r的公式为：

b)计算炮端反偏移S波场u_s,r的公式为：

和计算检端反传S波场ψ_s,b的公式为：

Images