CN105572734B - 一种以逆时偏移算法为引擎的波动方程初至走时层析方法 - Google Patents
一种以逆时偏移算法为引擎的波动方程初至走时层析方法 Download PDFInfo
- Publication number
- CN105572734B CN105572734B CN201410548681.5A CN201410548681A CN105572734B CN 105572734 B CN105572734 B CN 105572734B CN 201410548681 A CN201410548681 A CN 201410548681A CN 105572734 B CN105572734 B CN 105572734B
- Authority
- CN
- China
- Prior art keywords
- wave
- equation
- chromatography
- preliminary
- travel time
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种以逆时偏移算法为引擎的波动方程初至走时层析方法。包括：输入地震子波，观测系统，深度域初始速度场，观测地震记录的初至走时，以及走时残差范数的最小值；利用声波波动方程正演方法，计算初始速度场中正演的模拟地震记录，求解波动方程；计算模拟地震记录的初至波到达时；计算初至波走时残差，判断走时残差的范数是否小于预先设定的精度要求；计算波动方程走时敏感度矩阵；求解层析线性方程组；更新速度模型；输出反演得到的速度模型。优点是：采用声波波动方程具有更高的反演分辨率；更适用于描述波场的前向散射现象；将波动方程层析用于初至波走时层析中，得到了高精度的背景速度反演方法，能够反演大尺度的宏观背景速度场。
Description
技术领域
本发明涉及地震层析成像技术领域的层析速度反演技术，尤其涉及一种基于双程声波方程的初至走时层析速度反演方法。
背景技术
根据利用的数据信息不同，地震层析成像技术主要可分为走时层析，振幅层析及波形层析。其中，走时层析方法最为稳健。按正过程不同，走时层析主要可分为射线层析，胖射线层析，菲涅尔体层析和波动方程层析。
射线层析(Bishop等，1985；Tomographic determination of velocity anddepth in laterally varying media，Geophysics，50，903-923)基于高频射线理论，对初始模型精度的要求较低，但反演精度不高。同时，由于层析矩阵十分稀疏，矩阵的零空间比较大，因而求解时收敛较慢，且反演结果受射线照明的影响严重(Hu和Marcinkovich，2012；A sensitivity controllable target-oriented tomography algorithm，82nd SEGAnnual Meeting，Expanded Abstracts，1-5)。
考虑到地震波的传播不仅受到射线路径上的速度的影响，也受到射线周围速度结构的影响。胖射线层析(Vasco等，1995；Beyond ray tomography：Wavepaths and Fresnelvolumes，Geophysics，60(6)，1790-1804)将射线加宽，其物理背景是波传播是有一定宽度的，但没有严格的数学物理推导。此外，射线宽度的选择也有技巧，如固定宽度(Michelena和Harris，1991；Michelena R.J.，and J.M.Harris，1991，Tomographic traveltimeinversion using natural pixels)或随着射线传播距离增大而逐渐增大射线宽度(Xu等，2006；Enhanced tomography resolution by a fat ray technique.76th SEG AnnualMeeting，Expanded Abstracts，3354-3358)。但这些做法不能定量描述波动现象。
相比而言，菲涅尔带(体)层析(Yomogida，1992；Fresnel zone inversion forlateral heterogeneities in the earth，Pure and Applied Geophysics，138(3)，391-406)在改善射线层析矩阵稀疏性的同时，还考虑了地震波在第一菲涅尔带(体)内的有限频效应。但单频菲涅尔带只有在常速介质中存在解析表达式(Cerveny和Soares，1992；Fresnel volume ray tracing，Geophysics，57，902-915)。在变速介质中，可以用常速介质中主频对应的菲涅尔带来近似有限频带地震波的菲涅尔带宽度(刘玉柱等，2009；初至波菲涅耳体地震层析成像，2009，52(9)，2310-2320)；或计算炮点和检波点的走时场，利用单频菲涅尔带的公式计算菲涅尔带的范围(Watanabe等，1999；Seismic traveltimetomography using Fresnel volume approach：69thAnnual International Meeting，SEG，ExpandedAbstracts，1402-1405)。但这些做法只能近似描述变速介质中的带限地震波的菲涅尔带(体)范围。
波动方程走时层析没有引入高频近似假设，在理论上比射线类的层析具有更高的反演分辨率。在一定条件下(波动方程线性化近似，如一阶Born/Rytov近似等)能够描述波场(或相位)扰动与模型扰动的定量关系。Woodward(1992；Wave-equation tomography：Geophysics，57，15-26)引入了一阶Born及Rytov近似下“波路径”的概念(即波形与相位敏感度核函数)。Luo和Schuster(1991；Wave-equation traveltime inversion.Geophysics，56，645-653)利用互相关函数测量观测数据与模拟数据的时差，并给出了互相关准则下走时扰动对模型扰动的Fréchet导数(基于一阶Born近似)。Marquering等(1999；Three--dimensional sensitivity kernels for finite--frequency traveltimes：the banana-doughnut paradox.Geophys.J.Int.，137，805-815)则给出了扰动场与走时扰动的关系，进而可以导出一阶Born近似下走时敏感度核函数的表达。
考虑到一阶Rytov近似更适用于描述波场的前向散射现象，且适用条件较Born近似更为广泛。因此有必要发展一套在一阶Rytov近似下以波动方程为传播算子的层析速度反演技术。
发明内容
本发明的目的是针对上述现有技术存在的不足，提出了一种以逆时偏移算法为引擎的波动方程初至走时层析方法。该方法是一种高精度的背景速度反演方法，能够反演大尺度的宏观背景速度场。此外，由于正问题采用声波波动方程，没有引入高频近似，在理论上比射线类的层析具有更高的反演分辨率，因而可以定量的描述有限频带地震波的传播效应，对小尺度的速度异常有更好的反演能力。
本发明从波动方程线性化出发，并考虑地震波场的前向散射近似(一阶Rytov近似)，导出了基于波动方程的走时敏感度核函数在时空域的显式计算方法，其具体算法与著名的逆时偏移方法相似。而走时敏感度核函数就是层析线性方程组中的稀疏矩阵，通过求解层析线性方程组，进行速度更新，可以实现波动方程初至走时层析方法。
本发明的核心在于基于波动方程的走时敏感度核函数的构造，技术方案包括以下步骤：
步骤1：输入地震子波f(t)，观测系统，深度域初始速度场u(x)，观测地震记录的初至走时Tobs(xr，xs)，走时残差L2范数的最小值ε。
步骤2：利用声波波动方程正演方法，计算初始速度场中正演的模拟地震记录ucal(xr，t|xs)，求解如下波动方程：
其中，x为地下空间任意一点的坐标，xr和xs分别代表观测系统中检波器和震源的坐标，t代表时间，Δ代表拉普拉斯算子。
步骤3：计算模拟地震记录的初至波到达时Tcal(xr，xs)。
步骤4：计算初至波走时残差ΔT(xr，xs)，按照如下公式：
ΔT(xr，xs)＝Tobs(xr，xs)-Tcal(xr，xs)。 (8)
其中，Tobs(xr，xs)代表输入的观测数据的初至波到达时。
同时，判断走时残差的L2范数是否小于预先设定的精度要求(即)。若满足要求，跳转至步骤8。否则，执行下一步。
步骤5：计算波动方程走时敏感度矩阵K(x|xr，xs)，按照如下公式：
其中，u0(x|xs；t)为震源在xs处的地震子波产生的地震波场，x为地下空间任意一点的坐标。上标代表时间导数。p(x|xr；t)为检波器处逆时传播的地震波场，表示为如下公式
p(x|xr；t)＝g0(x，-t|xr，0)*u0(xr|xs；t)。 (10)
其中，g0(x，-t|xr，0)为声波波动方程的格林函数，代表在震源xr处，从时间t＝0开始的逆时传播过程，u0(xr|xs；t)为正演波场u0(x|xs；t)在检波器位置xr处的地震记录，即
显然，波动方程走时敏感度矩阵K(x|xr，xs)的计算与业内著名的波动方程逆时偏移有相类似的计算流程，即：震源处正演的模拟波场与检波器处逆时传播的地震记录对时间的积分。不同之处在于，公式(9)中需要的是正演的模拟波场的时间导数，并且需要能量校正因子因此，公式(9)在已有的逆时偏移算法基础之上，可以很容易实现。
步骤6：求解层析线性方程组：
Kδv＝ΔT (12)
其中，δvn×1为速度更新量，矩阵Km×n为公式(9)中计算的波动方程走时敏感度矩阵，ΔTm×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)，观测地震记录的初至走时Tobs(xr，xs)，走时残差L2范数的最小值ε；
步骤2：利用(声波)波动方程正演方法，计算初始速度场中正演的模拟地震记录ucal(xr，t|xs)，求解如下波动方程：
其中，x为地下空间任意一点的坐标，xr和xs分别代表观测系统中检波器和震源的坐标，t代表时间，Δ代表拉普拉斯算子。
步骤3：自动计算模拟地震记录的初至波到达时Tcal(xr，xs)；
步骤4：计算初至波走时残差ΔT(xr，xs)，按照如下公式：
ΔT(xr，xs)＝Tobs(xr，xs)-Tcal(xr，xs)。 (14)
其中，Tobs(xr，xs)代表输入的观测数据的初至波到达时。
同时，判断走时残差的L2范数是否小于预先设定的精度要求(即)。若满足要求，跳转至步骤8。否则，执行下一步。
步骤5：计算波动方程走时敏感度矩阵K(x|xr，xs)，按照如下公式：
其中，u0(x|xs；t)为震源在xs处的地震子波产生的地震波场，x为地下空间任意一点的坐标。上标代表时间导数。p(x|xr；t)为检波器处逆时传播的地震波场，表示为如下公式
p(x|xr；t)＝g0(x，-t|xr，0)*u0(xr|xs；t)。 (16)
其中，g0(x，-t|xr，0)为声波波动方程的格林函数，代表在震源xr处，从时间t＝0开始的逆时传播过程，u0(xr|xs；t)为正演波场u0(x|xs；t)在检波器位置xr处的地震记录，即
步骤6：求解层析线性方程组：
Kδv＝ΔT (18)
其中，δvn×1为速度更新量，矩阵Km×n为公式(15)中计算的波动方程走时敏感度矩阵，ΔTm×1为公式(14)中计算的初至波走时残差。m为观测系统中检波器的数目，n为速度模型网格离散化后的网格数目。
步骤7：更新速度模型，进行下一次迭代反演。跳转步骤2。
步骤8：输出反演得到的速度模型。
结合附图2-4，更加证明了本发明的技术效果。
Claims (1)
1.一种以逆时偏移算法为引擎的波动方程初至走时层析方法，其特征是基于波动方程的走时敏感度核函数的构造，包括以下步骤：
步骤1：输入地震子波f(t)，观测系统，深度域初始速度场v(x)，观测地震记录的初至走时Tobs(xr,xs)，走时残差L2范数的最小值ε；
步骤2：利用声波波动方程正演方法，计算初始速度场中正演的模拟地震记录ucal(xr,t|xs)，求解如下波动方程：
其中，x为地下空间任意一点的坐标，xr和xs分别代表观测系统中检波器和震源的坐标，t代表时间，Δ代表拉普拉斯算子；
步骤3：计算模拟地震记录的初至波到达时Tcal(xr,xs)；
步骤4：计算初至波走时残差ΔT(xr,xs)，按照如下公式：
ΔT(xr,xs)＝Tobs(xr,xs)-Tcal(xr,xs) (2)
其中，Tobs(xr,xs)代表输入的观测数据的初至波到达时；
同时，判断走时残差的L2范数是否小于预先设定的精度要求，即若满足要求，跳转至步骤8，否则，执行下一步；
步骤5：计算波动方程走时敏感度矩阵K(x|xr,xs)，按照如下公式：
其中，u0(x|xs；t)为震源在xs处的地震子波产生的地震波场，x为地下空间任意一点的坐标，上标·代表时间导数，p(x|xr；t)为检波器处逆时传播的地震波场，表示为如下公式
p(x|xr；t)＝g0(x,-t|xr,0)*u0(xr|xs；t) (4)
其中，g0(x,-t|xr,0)为声波波动方程的格林函数，代表在震源xr处，从时间t＝0开始的逆时传播过程，u0(xr|xs；t)为正演波场u0(x|xs；t)在检波器位置xr处的地震记录，即
步骤6：求解层析线性方程组:
Kδv＝ΔT (6)
其中，δvn×1为速度更新量，矩阵Km×n为公式(3)中计算的波动方程走时敏感度矩阵，ΔTm×1为公式(2)中计算的初至波走时残差，m为观测系统中检波器的数目，n为速度模型网格离散化后的网格数目；
步骤7：更新速度模型，进行下一次迭代反演，跳转步骤2；
步骤8：输出反演得到的速度模型。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410548681.5A CN105572734B (zh) | 2014-10-16 | 2014-10-16 | 一种以逆时偏移算法为引擎的波动方程初至走时层析方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410548681.5A CN105572734B (zh) | 2014-10-16 | 2014-10-16 | 一种以逆时偏移算法为引擎的波动方程初至走时层析方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105572734A CN105572734A (zh) | 2016-05-11 |
| CN105572734B true CN105572734B (zh) | 2018-10-09 |
Family
ID=55883075
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201410548681.5A Active CN105572734B (zh) | 2014-10-16 | 2014-10-16 | 一种以逆时偏移算法为引擎的波动方程初至走时层析方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105572734B (zh) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107807387B (zh) * | 2017-10-31 | 2019-08-27 | 中国科学技术大学 | 基于神经网络的地震初至波走时获取方法 |
| CN113281808B (zh) * | 2021-04-22 | 2023-10-20 | 南方海洋科学与工程广东省实验室(湛江) | 一种抗频散地震波正演方法、系统、装置及介质 |
| CN113504566B (zh) * | 2021-06-01 | 2024-04-30 | 南方海洋科学与工程广东省实验室(湛江) | 基于波动方程的地震反演方法、系统、装置及介质 |
| CN113791447B (zh) * | 2021-10-12 | 2023-06-20 | 同济大学 | 一种反射结构导引的反射波层析反演方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6643590B2 (en) * | 2002-01-04 | 2003-11-04 | Westerngeco, L.L.C. | Method for computing finite-frequency seismic migration traveltimes from monochromatic wavefields |
| CN101980052A (zh) * | 2010-09-28 | 2011-02-23 | 中国科学院地质与地球物理研究所 | 叠前逆时偏移成像的方法及装置 |
| CN102937721A (zh) * | 2012-11-07 | 2013-02-20 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 利用初至波走时的有限频层析成像方法 |
-
2014
- 2014-10-16 CN CN201410548681.5A patent/CN105572734B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6643590B2 (en) * | 2002-01-04 | 2003-11-04 | Westerngeco, L.L.C. | Method for computing finite-frequency seismic migration traveltimes from monochromatic wavefields |
| CN101980052A (zh) * | 2010-09-28 | 2011-02-23 | 中国科学院地质与地球物理研究所 | 叠前逆时偏移成像的方法及装置 |
| CN102937721A (zh) * | 2012-11-07 | 2013-02-20 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | 利用初至波走时的有限频层析成像方法 |
Non-Patent Citations (3)
| Title |
|---|
| Sensitivity kernels for time-distance inversion based on the Rytov approximation;J.M.Jensen 等;《Astronomy&Astrophysics》;20031231;第257-265页 * |
| Wave equation tomographic velocity inversion method based on the Born/Rytov approximation;Zhang Kai 等;《APPLIED GEOPHYSICS》;20131231;第10卷(第3期);第314-322页 * |
| 初至波菲涅尔体地震层析成像;刘玉柱 等;《地球物理学报》;20090930;第52卷(第9期);第2310-2320页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105572734A (zh) | 2016-05-11 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105572734B (zh) | 一种以逆时偏移算法为引擎的波动方程初至走时层析方法 | |
| Lin et al. | L1-norm regularization and wavelet transform: An improved plane-wave destruction method | |
| Pleshkevich et al. | Sixth-order accurate pseudo-spectral method for solving one-way wave equation | |
| Chen et al. | Fractional Laplacian viscoacoustic wave equation low-rank temporal extrapolation | |
| Feng et al. | Joint acoustic full-waveform inversion of crosshole seismic and ground-penetrating radar data in the frequency domain | |
| Yuan et al. | Analysis of attenuation and dispersion of Rayleigh waves in viscoelastic media by finite-difference modeling | |
| Jun et al. | Weighted pseudo-Hessian for frequency-domain elastic full waveform inversion | |
| CN107315192B (zh) | 基于二维各向同性介质的弹性波波场数值的模拟方法 | |
| Yu et al. | Convergence analyses of different modeling schemes for generalized Lippmann–Schwinger integral equation in piecewise heterogeneous media | |
| Liu et al. | First-arrival phase-traveltime tomography | |
| Bakir et al. | Modeling seismic attributes of Pn waves using the spectral-element method | |
| Zhao et al. | Double plane wave reverse time migration in the time domain | |
| Son et al. | An efficient 3D traveltime calculation using coarse-grid mesh for shallow-depth source | |
| Kamath et al. | Elastic full-waveform inversion of transmission data in 2D VTI media | |
| Naveed et al. | Elastic reverse time migration with non-staggered grid finite difference modeling approach | |
| Plessix et al. | 3D full-waveform inversion with a frequency-domain iterative solver | |
| Xia et al. | A time domain full waveform inversion method based on well-constrained regularization | |
| Mengxuan et al. | Time-domain full waveform inversion using the gradient preconditioning based on seismic wave energy: Application to the South China Sea | |
| Yang et al. | Time-windowed frequency domain full waveform inversion using phase-encoded simultaneous sources | |
| Yang et al. | A second-order adjoint truncated newton approach to time-domain multiparameter full waveform inversion in viscoacoustic medium | |
| Sun et al. | 2D frequency-domain elastic waveform inversion using matrix decomposition method | |
| Zhang et al. | An approach of full waveform inversion on GPU clusters and its application on land datasets | |
| Qu* et al. | Pre-stack elastic wave reverse time migration of irregular surface based on layered mapping method | |
| Jin et al. | Stack-based full wave-field velocity tomography | |
| Zhang* et al. | Full waveform inversion on the CPU/GPU heterogeneous platform and its application on land datasets |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |