CN115407394A - 一种基于波动方程的地震记录合成方法 - Google Patents

Abstract

本发明涉及地震数值模拟技术，具体公开了一种基于波动方程的地震记录合成方法。该方法包括：步骤1，引入极值求取函数，求取波峰波谷位置；步骤2，对地震数据进行分层，求取对应分层的纵波速度、横波速度及密度平均值；步骤3，根据矢量化反射率法的公式，求取慢度‑频率域反射系数；步骤4，对求取的慢度‑频率域反射系数进行重采样，获得角度‑频率域反射系数；步骤5，提取频率域子波；步骤6，得到角度‑时间域合成地震记录。本发明以实际地震资料为基础数据，引入了改进之后的自动寻峰算法，实现地下介质的层数的合理划分，在保证波场特征的同时减少层数，结合利用单道地震记录提取的地震子波，可以稳定地实现地下介质波场模拟。

Description

一种基于波动方程的地震记录合成方法

技术领域

本发明涉及地震数值模拟技术领域，具体包含一种基于波动方程的地震记录合成方法。

背景技术

地震数值模拟是根据对地下模型的构建来模拟波在地下介质中传播的方法和手段。通过地震数值模拟，可以构建地震波在地下实际介质中传播的机制，同时有助于了解地下的真实形态。目前来说，在地震勘探、地震学以及多种学科中，数值模拟都起到了很重要的作用。数值模拟的方法有很多种，基于波动方程的方法是最为实用且高效的方法之一，而矢量化反射率法以其精确的波场模拟特质及卓越的计算效率成为波动方程类方法的首选。

矢量化反射率法是一种在一维假设下求解弹性波方程从而实现层状介质波场模拟的方法，它能够提供波场的完全解，准确模拟多次波、转换波，并考虑了地震波在传播过程中的能量损耗等影响。对于各向同性介质，矢量化反射率法不涉及到纵横波解耦的问题，纵横波都能完全分离开。它是实现层状垂向非均匀介质全波场模拟的最有效的方法之一。对于精细的波场模拟以及叠前地震反演，矢量化反射率法的应用都起到极为关键的作用。

在利用矢量化反射率法进行波场模拟时，针对于已有的地震资料，需要人为估算分层的段数及地下介质的厚度。但是在实际工作中发现，针对于数据量庞大的地震资料进行人为的估算分层具有强烈的主观影响，而且对于整个工区的数据分层工作量极大，几乎不可能实现。此外，如若整个工区利用井上信息的同一个分层，不能良好地表征地下介质的横向变化。而且根据井信息计算的层数远远大于地震层数，极大地增加了计算负荷，增加了数倍的计算时间。因此，对地震数据进行合理的分层是提升矢量化反射率法的运算速度的关键，也是准确快速利用矢量化反射率法进行数值模拟亟待解决的问题。

现有的矢量化反射率法分层方法主要有两种：第一种是利用声波测井资料进行划分，但是由于测井分辨率高于地震，层数划分极多且分层较薄，致使运算速度大幅度降低且不能良好地表征地下介质的横向变化。第二种是根据先验认识人为划分地层，其缺点是主观因素极强且面对于实际数据工作量太大，不易实现。上述的分层方法毋庸置疑地降低了矢量化反射率法在正演时的效率，阻碍了该方法在地球物理勘探中的应用。

发明内容

为了降低分层的主观性并提升利用矢量化反射率法正演的运算效率，本发明提出了一种基于波动方程的地震记录合成方法。本方法以实际地震资料为基础数据，结合以精确的峰值求取算法，实现地下介质的层数划分，在保证波场特征的同时减少层数，实现地下介质的准确、快速模拟。

为解决上述技术问题，本发明采用的技术方案是：一种基于波动方程的地震记录合成方法，包括：

步骤1，抽取单道观测地震记录，引入极值求取函数，求取波峰波谷位置；

步骤2，根据求取的波峰波谷位置对地震数据进行分层，代入测井数据中的纵波速度、横波速度及密度，求取对应分层的纵波速度、横波速度及密度平均值；

步骤3，将求取的对应分层的纵波速度、横波速度及密度平均值代入矢量化反射率法的公式，求取慢度-频率域反射系数；

步骤4，对求取的慢度-频率域反射系数进行重采样，获得角度-频率域反射系数；

步骤5，从单道观测地震记录中提取频率域子波；

步骤6，将获得的角度-频率域反射系数和提取的频率域子波相乘，得到角度-频率域合成地震记录，进而转换到角度-时间域，得到角度-时间域合成地震记录。

优选的，所述步骤1包括：抽取单道观测地震记录，代入峰值求取函数计算波峰的位置，而后计算波形最大值，并与原始波形对应相减，求取波谷的位置；

P＝f(G(t)) (2)，

G'(t)＝G(t)-MAX(G(t)) (3)，

L＝f(G'(t)) (4)，

式(2)、(3)、(4)中，f表示极值求取函数，G(t)是单道观测地震记录，P是波峰的位置，L是波谷的位置，MAX代表求取最大值函数。

优选的，所述步骤2中，求取对应分层的纵波速度、横波速度及密度平均值的方式为：

式(5)中，

表示第n层的参数平均值，N表示第n个极值点的位置，M表示第n+1个极值点的位置，r_i表示这一段范围内任意参数值，所述参数包括纵波速度、横波速度以及密度。

优选的，所述步骤3中，求取慢度-频率域反射系数的方式为：

式(6)中，R_pp(p,ω)表示慢度-频率域反射系数，p是水平慢度，ω是角频率，υ₀为六元向量，υ₀(4)和υ₀(1)分别表示υ₀的第4个和第1个元素；

υ₀的计算方式为：

υ₀＝Q₀Q₁···Q_N-1υ_N, (7)，

υ_n＝Q_nυ_n+1, (8)，

υ_N＝[1 0 0 0 0 0]^T. (10)，

式(6)、(7)、(8)、(10)、(11)中，υ₀表示介质顶层的总反射响应，n代表任一层介质，N表示介质的总层数，υ_n+1表示地下第n+1层的响应，υ_N表示介质最底层的初始响应，Q_n表示传递矩阵，

和

表示第n层介质对地震波振幅的影响，

上标的+表示该层界面的下界面，

上标-表示该层界面的上界面，

和

各自都包含16个独立分量；矩阵E_n表示地震波经过第n层介质时其相位的变化。

优选的，所述步骤4包括：

引入p_n＝sinθ/α_n，其中p_n代表第n层的水平慢度，θ代表入射角，α_n代表第n层的纵波速度；

对求取的慢度-频率域反射系数进行重采样，获得角度-频率域的反射系数：

R_pp(θ,ω)＝resample(R_pp(p,ω)) (12)，

式(12)中，resample表示重采样过程，R_pp(θ,ω)表示角度-频率域反射系数。

优选的，所述步骤5中，从单道观测地震记录中提取频率域子波的方式为：

S(ω)＝waveletextract(G(t)) (13)，

式(13)中，S(ω)表示频率域子波，G(t)表示单道观测地震记录，waveletextract表示子波提取函数；

所述子波提取函数的方式为：

Z＝log(fft(G(t))) (14)，

wl＝ifft(Z)*W (15)，

S(ω)＝exp(fft(wl))； (16)，

式(14)、(15)、(16)中，log表示取对数，fft表示傅里叶变换，ifft表示傅里叶反变换，exp表示取指数，W表示时窗范围。

优选的，所述步骤6包括：将获得的角度-频率域反射系数和子波在频率域相乘，得到角度-频率域合成地震记录，通过傅里叶反变换至角度-时间域，得到角度-时间域合成地震记录：

式(17)中，d(θ,t)为角度-时间域合成地震记录，S(ω)表示频率域子波，R_pp(θ,ω)为角度-频率域反射系数,ω是角频率。

本发明方法以实际地震资料为基础数据，引入了改进之后的自动寻峰算法，实现地下介质的层数的合理划分，在保证波场特征的同时减少层数，结合利用单道地震记录提取的地震子波，可以稳定地实现地下介质波场模拟，从而提高正演的计算效率和准确度。

附图说明

图1为本发明提供的一种基于波动方程的地震记录合成方法一实施例的流程图。

图2为某井实际测得的纵波速度、横波速度及密度示意图。

图3为单道地震记录求取的波峰波谷位置结果示意图。

图4为某地层根据波峰波谷位置划分的参数分层结果示意图。

图5是从单道地震记录中提取的子波(时间域)示意图。

图6a是传统方法合成地震记录示意图。

图6b是本发明方法合成地震记录示意图。

图6c是传统方法合成地震记录与本发明方法合成地震记录的残差示意图。

具体实施方式

为便于对本发明的方法及达到的效果有进一步的了解，现结合附图实例详细说明如下，然而附图仅提供参考与说明用，并非用来对本发明加以限制。

本发明构建了一种基于波动方程的地震记录合成方法。由于地震记录单道波形的波峰、波谷代表了地震波在传播过程中经历的地层界面，所以其波形的极值点可以作为划分地层的依据，在水平层状假设情况下，两个界面之间的地层拥有相同的地震波传播速度和密度。通过引入波峰寻取算法并加以改进，可以得到准确的波形极值搜索结果。根据极值点的位置求取相对应的地层厚度，并对纵波速度、横波速度和密度数据进行分层，代入矢量化反射率法求取反射系数，并在频率域与求取的地震子波相乘，实现地下介质的正演模拟。

式(1)中，R_pp(p,ω)表示慢度-频率域反射系数，S(ω)是频率域子波，d(p,t)表示慢度-时间域合成地震记录，e是自然常数，ω是角频率，i是虚数单位。由式(1)可知，合成地震记录的精确程度和反射系数及子波具有密切的关系。

进一步的，根据斯奈尔定律，慢度P等于sinθ/α,其中θ是角度，α是地震波速度，可以将慢度-时间域合成地震记录转换成角度-时间域合成地震记录，只有将合成地震记录转换到角度域才能与实际地震记录进行对比、运算。

如图1所示，本发明提供的一种基于波动方程的地震记录合成方法一实施例的流程包括：

步骤1，抽取单道观测地震记录，引入极值求取函数，求取波峰波谷位置；

步骤2，根据求取的波峰波谷位置对地震数据进行分层，代入测井数据中的纵波速度、横波速度及密度，求取对应分层的纵波速度、横波速度及密度平均值；

步骤3，将求取的对应分层的纵波速度、横波速度及密度平均值代入矢量化反射率法的公式，求取慢度-频率域反射系数；

步骤4，对求取的慢度-频率域反射系数进行重采样，获得角度-频率域反射系数；

步骤5，从单道观测地震记录中提取频率域子波；

步骤6，将获得的角度-频率域反射系数和提取的频率域子波相乘，得到角度-频率域合成地震记录，进而转换到角度-时间域，得到角度-时间域合成地震记录。

步骤1，抽取单道观测地震记录，引入极值求取函数，求取波峰波谷位置，具体包括：抽取单道观测地震记录，代入峰值求取函数计算波峰的位置，而后计算波形最大值，并与原始波形对应相减，求取波谷的位置；

P＝f(G(t)) (2)，

G'(t)＝G(t)-MAX(G(t)) (3)，

L＝f(G'(t)) (4)，

式(2)、(3)、(4)中，f表示极值求取函数，G(t)是单道观测地震记录(即实际地震记录)，P是波峰的位置，L是波谷的位置，MAX代表求取最大值函数。

步骤2，根据所求波峰波谷位置及其对应深度，对地震数据进行分层，代入测井数据中的速度和密度值，计算每层的参数值(纵波速度、横波速度及密度)，具体计算过程可表示如下：

式(5)中，

表示第n层的参数平均值，N表示第n个极值点的位置，M表示第n+1个极值点的位置，r_i表示这一段范围内任意参数值，具体包括纵波速度、横波速度以及密度；

测井数据是由测井实测而来，在地球物理之中当作准确的数据，其测得的时间与深度的关系是验证合成地震记录准确性的标准。

步骤3，根据矢量化反射率法的公式，求取慢度-频率域反射系数，具体包括：

式(6)中，R_pp(p,ω)表示慢度-频率域反射系数，p是水平慢度，ω是角频率，该反射系数是由一个六元向量υ₀的两个元素的比值组成，υ₀(4)和υ₀(1)分别表示六元向量υ₀的第4个和第1个元素，该六元向量的计算方法如下：

υ₀＝Q₀Q₁···Q_N-1υ_N (7)，

υ_n＝Q_nυ_n+1 (8)，

Q_n＝E_nF_n (9)，

υ_N＝[1 0 0 0 0 0]^T (10)，

在以上各式中，所有的下标皆代表其所在的层数，n代表任一层介质，N表示介质的总层数；Q_n表示第n层的传递矩阵，六元向量υ_n表示地下第n层的响应，υ_n+1表示地下第n+1层的响应，而υ₀表示介质顶层的总反射响应，υ_N表示介质最底层的初始响应，其上标T代表该向量的转置。矩阵E_n表达了地震波经过第n层介质时其相位的变化，F_n表述了第n层介质对地震波振幅的影响；

为了方便于后续计算，重新将Q_n写成如下形式

和

表示第n层介质对地震波振幅的影响，其上标的+表示该层界面的下界面，而-表示该层界面的上界面，各自都包含16个独立分量；矩阵E_n表达了地震波经过第n层介质时其相位的变化。

步骤4，在求取完慢度-频率域反射系数之后，通过引入水平慢度和入射角的关系p_n＝sinθ/α_n，此处p_n代表第n层的水平慢度，θ代表入射角，α_n代表第n层的纵波速度，对R_pp(p,ω)进行重采样，获得角度-频率域的反射系数。

R_pp(θ,ω)＝resample(R_pp(p,ω)) (12)，

式(12)中，resample表示重采样过程，R_pp(θ,ω)表示角度-频率域反射系数。

步骤5，从单道观测地震记录中提取频率域子波：

S(ω)＝waveletextract(G(t)) (13)，

式(13)中，S(ω)表示频率域子波，G(t)表示单道观测地震记录，waveletextract表示子波提取函数，其具体过程为：

Z＝log(fft(G(t))) (14)，

wl＝ifft(Z)*W (15)，

S(ω)＝exp(fft(wl))； (16)，

式(14)、(15)、(16)中，log表示取对数，fft表示傅里叶变换，ifft表示傅里叶反变换，exp表示取指数，而W表示时窗范围。

步骤6，将获得的角度-频率域反射系数和子波在频率域相乘，得到角度-频率域合成地震记录，通过傅里叶反变换至角度-时间域，得到角度-时间域合成地震记录：

式(17)中，d(θ,t)为角度-时间域合成地震记录，S(ω)表示频率域子波，R_pp(θ,ω)是由式(12)求得的角度-频率域的反射系数，e是自然常数，ω是角频率，i是虚数单位，π是圆周率，至此便完成了整个正演过程。

选定工区中某块地震数据和相应区域内的井资料，利用本发明方法进行地震资料分层并完成正演。图2展示了测井曲线中的纵波速度、横波速度及密度。图3表示对单道地震记录进行波峰波谷求取的结果，空心点代表波峰，实心点代表波谷。可以发现本方法精准地计算了每个波峰波谷的位置，可以很明确地区分出地层界面的位置。根据波峰波谷对应的位置，代入测井的相应数据，得到速度和密度的分层结果，如图4所示，此过程可以根据地质资料的先验认识调整分层的精细程度，本图仅以实际地震记录分层的结果划分，未做调整。而后，根据图3所示的单道地震记录可以提取地震子波，其结果如图5所示。为验证本方法的准确性，将传统的测井分层的合成记录与本方法进行对比，得到图6a、图6b、图6c所示的结果，通过对比其残差可以发现，本方法合成单个道集(30°以内)具有几乎和传统方法一样的精度。但是，通过记录传统方法和本方法的运算时间(如表1)，本方法将速度提升了约45倍。

表1传统方法和本方法的运算时间对比

名称	传统方法	本方法
			运算时间(S)	4.5716	0.1002

从以上的图表中可以看出，本发明方法不仅具有极其高的准确性，更在速度上有了极大的提升。这证明了本发明方法的合理性和正确性。

以上所述实施例仅是为充分说明本发明而所举的较佳的实施例，本发明的保护范围不限于此。本技术领域的技术人员在本发明基础上所作的等同替代或变换，均在本发明的保护范围之内。

Claims (7)

1.一种基于波动方程的地震记录合成方法，其特征在于，包括：

步骤1，抽取单道观测地震记录，引入极值求取函数，求取波峰波谷位置；

步骤2，根据求取的波峰波谷位置对地震数据进行分层，代入测井数据中的纵波速度、横波速度及密度，求取对应分层的纵波速度、横波速度及密度平均值；

步骤3，将求取的对应分层的纵波速度、横波速度及密度平均值代入矢量化反射率法的公式，求取慢度-频率域反射系数；

步骤4，对求取的慢度-频率域反射系数进行重采样，获得角度-频率域反射系数；

步骤5，从单道观测地震记录中提取频率域子波；

步骤6，将获得的角度-频率域反射系数和提取的频率域子波相乘，得到角度-频率域合成地震记录，进而转换到角度-时间域，得到角度-时间域合成地震记录。

2.如权利要求1所述一种基于波动方程的地震记录合成方法，其特征在于，

所述步骤1包括：抽取单道观测地震记录，代入峰值求取函数计算波峰的位置，而后计算波形最大值，并与原始波形对应相减，求取波谷的位置；

P＝f(G(t)) (2)，

G'(t)＝G(t)-MAX(G(t)) (3)，

L＝f(G'(t)) (4)，

式(2)、(3)、(4)中，f表示极值求取函数，G(t)是单道观测地震记录，P是波峰的位置，L是波谷的位置，MAX代表求取最大值函数。

3.如权利要求1所述一种基于波动方程的地震记录合成方法，其特征在于，

所述步骤2中，求取对应分层的纵波速度、横波速度及密度平均值的方式为：

式(5)中，

表示第n层的参数平均值，N表示第n个极值点的位置，M表示第n+1个极值点的位置，r_i表示这一段范围内任意参数值，所述参数包括纵波速度、横波速度以及密度。

4.如权利要求1所述一种基于波动方程的地震记录合成方法，其特征在于，

所述步骤3中，求取慢度-频率域反射系数的方式为：

式(6)中，R_pp(p,ω)表示慢度-频率域反射系数，p是水平慢度，ω是角频率，υ₀为六元向量，υ₀(4)和υ₀(1)分别表示υ₀的第4个和第1个元素；

υ₀的计算方式为：

υ₀＝Q₀Q₁···Q_N-1υ_N, (7)，

υ_n＝Q_nυ_n+1, (8)，

υ_N＝[1 0 0 0 0 0]^T. (10)，

式(6)、(7)、(8)、(10)、(11)中，υ₀表示介质顶层的总反射响应，n代表任一层介质，N表示介质的总层数，υ_n+1表示地下第n+1层的响应，υ_N表示介质最底层的初始响应，Q_n表示传递矩阵，

和

表示第n层介质对地震波振幅的影响，

上标的+表示该层界面的下界面，

上标-表示该层界面的上界面，

和

各自都包含16个独立分量；矩阵E_n表示地震波经过第n层介质时其相位的变化。

5.如权利要求1所述一种基于波动方程的地震记录合成方法，其特征在于，

所述步骤4包括：

引入p_n＝sinθ/α_n，其中p_n代表第n层的水平慢度，θ代表入射角，α_n代表第n层的纵波速度；

对求取的慢度-频率域反射系数进行重采样，获得角度-频率域的反射系数：

R_pp(θ,ω)＝resample(R_pp(p,ω)) (12)，

式(12)中，resample表示重采样过程，R_pp(θ,ω)表示角度-频率域反射系数。

6.如权利要求1所述一种基于波动方程的地震记录合成方法，其特征在于，

所述步骤5中，从单道观测地震记录中提取频率域子波的方式为：

S(ω)＝waveletextract(G(t)) (13)，

式(13)中，S(ω)表示频率域子波，G(t)表示单道观测地震记录，waveletextract表示子波提取函数；

所述子波提取函数的方式为：

Z＝log(fft(G(t))) (14)，

wl＝ifft(Z)*W (15)，

S(ω)＝exp(fft(wl))； (16)，

式(14)、(15)、(16)中，log表示取对数，fft表示傅里叶变换，ifft表示傅里叶反变换，exp表示取指数，W表示时窗范围。

7.如权利要求1所述一种基于波动方程的地震记录合成方法，其特征在于，

所述步骤6包括：将获得的角度-频率域反射系数和子波在频率域相乘，得到角度-频率域合成地震记录，通过傅里叶反变换至角度-时间域，得到角度-时间域合成地震记录：

式(17)中，d(θ,t)为角度-时间域合成地震记录，S(ω)表示频率域子波，R_pp(θ,ω)为角度-频率域反射系数,ω是角频率。

Images