CN109946742B - Pure qP wave seismic data simulation method in TTI medium - Google Patents
Pure qP wave seismic data simulation method in TTI medium Download PDFInfo
- Publication number
- CN109946742B CN109946742B CN201910250705.1A CN201910250705A CN109946742B CN 109946742 B CN109946742 B CN 109946742B CN 201910250705 A CN201910250705 A CN 201910250705A CN 109946742 B CN109946742 B CN 109946742B
- Authority
- CN
- China
- Prior art keywords
- wave
- representing
- field
- grid
- simulation
- 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
Images
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
The embodiment of the specification discloses a pure qP wave seismic data simulation method in a TTI medium. According to the scheme provided by the embodiment of the specification, the qP waves coupled together in the TTI medium are decoupled from the qSV waves, compared with the forward modeling of the anisotropic full elastic wave fluctuation equation, the decoupling pure qP wave simulation method naturally solves the problem of difficult wave field separation, and the practical application proves that the pure qP wave simulation is well applied to longitudinal wave exploration; the simulation process has the advantages that pseudo transverse wave interference does not exist at all and the anisotropy parameter is stable; when the anisotropic parameters are changed violently, the wave field can still be stably transmitted in the simulation process, and the stability is better.
Description
Technical Field
The invention relates to a pure qP wave seismic data simulation method in a TTI medium, and belongs to the field of exploration geophysics.
Background
The earth medium develops widely anisotropic properties, among which are The Transverse Isotropic (TTI) media with tilted symmetry axis.
In data processing for such media, acoustically approximated coupled (qP waves propagate with residual spurious qSV waves coupled) pseudoacoustic wave equations can describe more accurately the kinematic characteristics of qP waves, but they suffer from spurious shear wave interference and can produce instability at the anisotropy parameter < n >. Meanwhile, the conventional pure qP wave method mostly adopts a pseudo-spectrum method to carry out numerical simulation, and the pseudo-spectrum method has the defects of low calculation efficiency and high requirement on stability conditions.
Based on this, a more stable and efficient method for simulating seismic data is needed.
Disclosure of Invention
The invention aims to provide an accurate and stable pure qP wave seismic data simulation method in a TTI medium.
In order to solve the technical problems, the invention adopts the following technical scheme:
a pure qP wave seismic data simulation method in a TTI medium comprises the following steps:
obtaining an initial seismic geological anisotropy parameter field, the anisotropy parameters including, phi, and vpWherein the sum is the Thomsen parameter of the medium,. phi.pRepresenting the longitudinal wave velocity of the medium;
dividing the parameter field into a plurality of grids;
for each grid, obtaining a simulation equation coefficient λ associated with the anisotropic parameter1,λ2,λ3；
Adopting a preset wave field propagation operator, arranging shot points and detection point positions based on the simulation equation coefficient related to the grid position and the anisotropic parameter field, and then carrying out numerical simulation to obtain a stress field value received by the detector, wherein the wave field propagation operator comprises:
c1＝(1+λ3+λ1cos4φ+2cos2φ+λ2cos2φ)，
c2＝(1+λ3+λ1sin4φ+2sin2φ+λ2sin2φ)，
c3＝(2+2λ3+λ2+2+1.5λ1sin22φ)，
c4＝(-2sin2φ-λ2sin2φ-2λ1sin2φcos2φ)，
c5＝(-2sin2φ-λ2sin2φ-2λ1sin2φsin2φ)；
wherein u ispRepresenting the value of the stress field, uqShowing the value of auxiliary stress field, phi is the included angle between the anisotropic symmetry axis and the vertical direction, t shows the propagation time, i and j respectively show the positions of the vertical and horizontal grid points, and k shows the discrete time value△ x denotes the lateral grid spacing, △ z denotes the longitudinal grid spacing, △ t denotes the temporal sampling interval, a0And annRepresenting difference coefficients, FFT representing fast Fourier transform, FFT-1Representing the inverse of the fast Fourier transform, kxAnd kzRespectively, represents the wavenumber values in the vertical and horizontal directions, N represents the difference order, and x and z represent the abscissa and ordinate, respectively.
According to the scheme provided by the embodiment of the specification, the qP waves coupled together in the TTI medium are decoupled from the qSV waves, compared with the forward modeling of the anisotropic full elastic wave fluctuation equation, the decoupling pure qP wave simulation method naturally solves the problem of difficult wave field separation, and the practical application proves that the pure qP wave simulation is well applied to longitudinal wave exploration; the simulation process has the advantages that pseudo transverse wave interference does not exist at all and the anisotropy parameter is stable; when the anisotropic parameters are changed violently, the wave field can still be stably transmitted in the simulation process, and the stability is better.
Drawings
Fig. 1 is a schematic flowchart of a pure qP wave seismic data simulation method in a TTI medium according to an embodiment of the present disclosure;
FIG. 2 is a schematic diagram comparing simulation results of conventional coupled pseudoacoustic wave equations provided herein;
FIG. 3 is another schematic diagram comparing simulation results with conventional coupled pseudoacoustic wave equations provided herein;
FIG. 4 is a diagram illustrating the relationship between the anisotropy parameter and the position of a model provided in an embodiment of the present disclosure;
FIG. 5 is a diagram illustrating a comparison between a wave field snapshot obtained by a conventional method and a method provided by an embodiment of the present disclosure;
FIG. 6 is a schematic diagram of a comparison of the method provided by the embodiments of the present disclosure with seismic shot recordings from a conventional method.
Detailed Description
In order to make the objects, technical solutions and advantages of the present application more apparent, the technical solutions of the present application will be described in detail and completely with reference to the following specific embodiments of the present application and the accompanying drawings. It should be apparent that the described embodiments are only some of the embodiments of the present application, and not all of the embodiments. All other embodiments obtained by a person skilled in the art based on the embodiments in the present specification without any inventive step are within the scope of the present application.
The wavefield propagation operators used in the simulation process of the present application will first be described.
In the seismic field, the qP wave phase velocity formula equation is:
to facilitate solving the wave equation, the present invention approximates equation (1) with the approximate qP-wave phase velocity equation (2). Wherein v ispRepresenting the qP wave propagation phase velocity, theta representing the phase angle, and the Thomsen parameter value of the medium, vp0Representing the medium longitudinal wave velocity.
The formula (2) is:
wherein the coefficient value λ in the formula (2)1(x,z),λ2(x,z),λ3(x, z) is a position-dependent function, λ1(x,z),λ2(x,z),λ3The value of (x, z) is expressed by the sum of anisotropy parameters, which can be obtained by approximating equation (1), specifically, by the best-squares approximation. For example, the specific principle of maximizing equation (2) using the principle of the best-squares approximation to approximate equation (1) is as follows:
if f (x) is the interval [ a, b ]]The continuous function of (a) to (b),
is [ a, b ]]A linear independent function system of above, and
in [ a, b ]]All are continuous, doDefinite generalized polynomial
Coefficient a of0,a1,…,anTo make
Function thus obtained
Is referred to as f (x) in [ a, b ]]The best square approximation of (a). In the present embodiment, f (x) corresponds to equation (1),
i.e. corresponds to equation (2).
The solution to the optimal square approximation requires the use of an integral formula, here a complex Simpson formula, with the following specific principles:
will be interval [ a, b]n is divided equally: division point xk＝a+kh，
In the interval [ xk,xk+1]Simpson's formula for k-0, 1, … n-1
Wherein the content of the first and second substances,
based on determined lambda1(x,z),λ2(x,z),λ3(x, z), mathematically transforming equation (2) to obtain the following decoupled qP wave equation.
The qP wave fluctuation equation is:
wherein u isp(x, z, t) is the stress value, vp0Is the particle longitudinal wave velocity, λ1,λ2,λ3The approximation coefficient is obtained from the front, phi is the included angle between the anisotropy symmetry axis and the vertical direction, t represents the propagation time, and x and z respectively represent the horizontal and vertical coordinates.
Equation (3) can only be calculated by pseudo-spectrometry, resulting in an excessive amount of calculation. Equation (3) is rewritten here as equation (4), and equation (4) consists of a modified acoustic wave equation and poisson's equation.
Wherein u isq(x, z, t) are auxiliary stress field values and have no practical physical significance.
At this time, it can be seen that the first part in equation (4) is a modified acoustic wave equation, and the second part is a poisson equation, which are separated. The numerical simulation can be performed on equation (4) using a finite difference-pseudo-spectral mixture method. That is, the modified acoustic wave equation part is solved by applying the temporal second-order spatial high-order difference format, and the poisson equation part is calculated by using the pseudo-spectral method, then the qP wave field propagation operator of the above formula can be expressed as:
wherein, c1＝(1+λ3+λ1cos4φ+2cos2φ+λ2cos2φ)，
c2＝(1+λ3+λ1sin4φ+2sin2φ+λ2sin2φ)，
c3＝(2+2λ3+λ2+2+1.5λ1sin22φ)，
c4＝(-2sin2φ-λ2sin2φ-2λ1sin2φcos2φ)，
c5＝(-2sin2φ-λ2sin2φ-2λ1sin2φsin2φ)，
Wherein u ispRepresenting the value of the stress field, uqRepresenting auxiliary stress field values, i, j respectively representing longitudinal and transverse grid point positions, k representing discrete time values, △ x representing transverse grid spacing, △ z representing longitudinal grid spacing, △ t representing time sampling intervals, a0And annRepresenting difference coefficients, FFT representing fast Fourier transform, FFT-1Representing the inverse of the fast Fourier transform, kxAnd kzRespectively representing the values of waves in the longitudinal and transverse directions, λ1(x,z),λ2(x,z),λ3(x, z) are values expressed by anisotropy parameters and the magnitude of which is determined when equation (2) approximates equation (1), phi is the angle between the anisotropy symmetry axis and the vertical direction, m and N are intermediate variables in calculation, there is no practical physical meaning, N represents the difference order, and when N is 5, 10-order difference accuracy is obtained.
The foregoing section explains and explains the wave field propagation operator used in the embodiment of the present specification, and a specific use manner is shown in fig. 1, where fig. 1 is a schematic flow chart of a pure qP wave seismic data simulation method in a TTI medium provided in the embodiment of the present specification, including:
s101, acquiring an initial seismic geological anisotropy parameter field, wherein anisotropy parameters comprise phi, phi and vpWherein the sum is the Thomsen parameter of the medium,. phi.pRepresenting the longitudinal wave velocity of the medium;
in particular, a geological model map may be drawn using geological geophysical models, field geological outcrops, drilling, geophysical prospecting, well logging data, and an anisotropic parameter field associated with the initial seismic data may be obtained by model filling.
S103, dividing the parameter field into a plurality of grids.
The grid division is generally uniform, that is, the grid distance can be randomly specified. In practical application, selection is performed in order to realize the principles that the smaller calculation amount is small, the memory requirement is small, the grid distance dx and the time sampling interval dt meet the stability condition, and the numerical dispersion is small, so that the forward modeling result is efficient, stable and accurate. The stability conditions to be met here are:
where v denotes the velocity field, △ t is the time sample interval, △ d is the grid spacing.
S105, acquiring a simulation equation coefficient lambda related to the anisotropic parameter for each grid1,λ2,λ3。
As previously described, the simulation equation coefficients may be determined based on the anisotropy parameter field that has been determined, using the principle of best squares approximation for equation (1) based on equation (2).
And S107, determining a shot point and a detection point position by adopting a preset wave field propagation operator based on the simulation equation coefficient related to the grid position and the anisotropic parameter field, and then performing numerical simulation to obtain a stress field value received by the detector.
The wavefield propagation operator has been described in detail above. In practical application, the presetting can be performed in the form of a functional module or an algorithm module. The shot and wave detector locations may be located at any position on the grid, typically on the first layer of the grid (i.e., the earth's surface).
As described above, in the numerical simulation process, the modified acoustic wave equation part is solved by applying the temporal second-order spatial high-order difference format, and the poisson equation part is calculated by using the pseudo-spectral method, so that numerical simulation based on the decoupled pure qP wave is realized, and stress field data about wave field propagation is obtained.
Fig. 2 is a diagram comparing simulation results with conventional coupled pseudoacoustic wave equations provided in the present specification. The part a is a qP wave field snapshot at 0.4s time (0.24, ═ 0.12,. phi. 30 °) obtained by a traditional coupled sound wave equation, and the part b is a qP wave field snapshot at 0.4s time (0.24,. phi. 0.12,. phi. 30 °) provided by the present application. As can be seen from FIG. 2, the qP wave simulation result of the invention has no pseudo-shear wave interference compared with the simulation result of the conventional coupled sound wave equation. Namely, the invention completely decouples the qP wave and the qSV wave, and the traditional coupling simulation result of the sound wave simulation equation still couples the residual qSV wave. In the figure, the abscissa represents the length x, and the ordinate represents the depth z.
Fig. 3 is another comparison diagram of simulation results of the conventional coupled pseudoacoustic wave equation provided in the present specification. The part a is a qP wave field snapshot at 0.4s time (0.12, ═ 0.24,. phi. 30 °) obtained by a traditional coupled sound wave equation, and the part b is a qP wave field snapshot at 0.4s time (0.12,. phi.24,. phi. 30 °) provided by the present application. As can be seen from fig. 3, the method provided by the present application is still stable at the anisotropy parameter < while the traditional coupled pseudoacoustic wave equation is not stable. In the figure, the abscissa represents the length x, and the ordinate represents the depth z.
According to the scheme provided by the embodiment of the specification, the qP waves and the qSV waves which are coupled together are decoupled, compared with the forward modeling of the anisotropic full elastic wave fluctuation equation, the decoupling pure qP wave simulation method naturally solves the problem of difficult wave field separation, and the practical application proves that the pure qP wave simulation is well applied to longitudinal wave exploration; the simulation process has the advantages that pseudo transverse wave interference does not exist at all and the anisotropy parameter is stable; when the anisotropic parameters are changed violently, the wave field can still be stably transmitted in the simulation process, and the stability is better.
In order to verify the applicability of the qP wave simulation method to the complex model, the scheme provided by the invention is used in a BP model. The parameters of this model (part of the BP model) are shown in fig. 4, from which it can be seen that there is a strong anisotropic dip change in the model. Fig. 4 is a schematic diagram of a relationship between an anisotropy parameter and a position of a model provided in an embodiment of the present disclosure. In the figure, the abscissa represents the length x, and the ordinate represents the depth z.
The conventional coupled pseudoacoustic wave equation is chosen for reference herein. Fig. 5 is a schematic diagram of a comparison between the wave field snapshot obtained by the method provided in the embodiment of the present disclosure and the wave field snapshot obtained by the conventional method, where the abscissa is length x and the ordinate is depth z. FIG. 6 is a schematic diagram of a comparison between simulated shot records obtained by the method of the present disclosure and conventional methods, where the abscissa is length x and the ordinate is time t. It can be seen that the solution proposed by this patent is stable in the regions of strong variation of the anisotropic tilt angle. The forward result obtained by the traditional coupled sound wave simulation equation is unstable in the area with the violent change of the inclination angle. The effectiveness of the seismic simulation method based on the optimal square approximation provided by the patent is proved through BP model test.
Correspondingly, the embodiment of the present application further provides a computer device, which includes a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the pure qP wave seismic data simulation method in the TTI medium.
The embodiments in the present specification are described in a progressive manner, and the same and similar parts among the embodiments are referred to each other, and each embodiment focuses on the differences from the other embodiments. Especially, as for the device, apparatus and medium type embodiments, since they are basically similar to the method embodiments, the description is simple, and the related points may refer to part of the description of the method embodiments, which is not repeated here.
Claims (4)
1. A forward modeling method for pure qP wave seismic data in a TTI medium comprises the following steps:
obtaining an initial seismic geological anisotropy parameter field, the anisotropy parameters including, phi, and vpWherein the sum is the Thomsen parameter of the medium,. phi.pRepresenting the longitudinal wave velocity of the medium;
dividing the parameter field into a plurality of grids;
for each grid, obtaining a simulation equation coefficient λ associated with the anisotropic parameter1,λ2,λ3；
Adopting a preset wave field propagation operator, determining a shot point and a detection point position based on the simulation equation coefficient related to the grid position and the anisotropic parameter field, and then carrying out numerical simulation to obtain a stress field value received by the detector, wherein the wave field propagation operator comprises:
c1＝(1+λ3+λ1cos4φ+2cos2φ+λ2cos2φ)，
c2＝(1+λ3+λ1sin4φ+2sin2φ+λ2sin2φ)，
c3＝(2+2λ3+λ2+2+1.5λ1sin22φ)，
c4＝(-2sin2φ-λ2sin2φ-2λ1sin2φcos2φ)，
c5＝(-2sin2φ-λ2sin2φ-2λ1sin2φsin2φ)；
wherein u ispRepresenting the value of the stress field, uqRepresenting auxiliary stress field values, t representing propagation time, i representing transverse grid point positions, j representing longitudinal grid point positions, k representing discrete time values, Δ x representing transverse grid spacing, Δ z representing longitudinal grid spacing, Δ t representing time sampling interval, a0And annRepresenting difference coefficients, FFT representing fast Fourier transform, FFT-1Express fastInverse of fast Fourier transform, kxRepresenting the value of the wave in the transverse direction, kzRepresents the wave number value in the vertical direction, N represents the difference order, and x and z represent the abscissa and ordinate, respectively.
2. The method of claim 1, dividing the parameter field into a plurality of grids, comprising:
simulating stability conditions based on numerical values
And determining the space between grids, and carrying out grid division on the speed field, wherein v represents the speed field, delta t is a time sampling interval, and delta d is the grid space.
3. The method of claim 1, obtaining, for each grid, a simulation equation coefficient λ associated with the anisotropic parameter1,λ2,λ3The method comprises the following steps:
using the best square approximation method to make the equation
Equation of approximate qP wave phase velocity
Determining lambda1(x,z),λ2(x,z),λ3(x, z) wherein vp0Is the particle longitudinal velocity.
4. A computer device comprising a memory, a processor and a computer program stored on the memory and executable on the processor, wherein the processor implements the method of any of claims 1 to 3 when executing the program.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910250705.1A CN109946742B (en) | 2019-03-29 | 2019-03-29 | Pure qP wave seismic data simulation method in TTI medium |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910250705.1A CN109946742B (en) | 2019-03-29 | 2019-03-29 | Pure qP wave seismic data simulation method in TTI medium |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109946742A CN109946742A (en) | 2019-06-28 |
| CN109946742B true CN109946742B (en) | 2020-09-11 |
Family
ID=67012313
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910250705.1A Active CN109946742B (en) | 2019-03-29 | 2019-03-29 | Pure qP wave seismic data simulation method in TTI medium |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109946742B (en) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111025386B (en) * | 2019-12-13 | 2020-11-17 | 成都理工大学 | Vertical and horizontal wave separation method without separation false image |
| CN112764105B (en) * | 2020-10-16 | 2022-07-12 | 中国石油大学(华东) | HTI medium quasi-longitudinal wave forward simulation method and device, storage medium and processor |
| CN115685337A (en) * | 2022-10-24 | 2023-02-03 | 中国石油大学(华东) | Anisotropic elastic wave decoupling method and device and computer equipment |
| CN116559941B (en) * | 2023-04-07 | 2024-03-12 | 中国地质调查局油气资源调查中心 | Norris-KG model-based earthquake longitudinal wave simulation analysis method |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104216011B (en) * | 2013-06-05 | 2017-08-04 | 上海青凤致远地球物理地质勘探科技有限公司 | A kind of stable qP ripple reverse-time migration methods of TTI media |
| CN105044771B (en) * | 2015-08-05 | 2017-10-27 | 北京多分量地震技术研究院 | Three-dimensional TTI two-phase medias seismic wave field method for numerical simulation based on finite difference calculus |
| US11237283B2 (en) * | 2017-07-13 | 2022-02-01 | Exxonmobil Upstream Research Company | Visco-pseudo-elastic TTI FWI/RTM formulation and implementation |
-
2019
- 2019-03-29 CN CN201910250705.1A patent/CN109946742B/en active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN109946742A (en) | 2019-06-28 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109946742B (en) | Pure qP wave seismic data simulation method in TTI medium | |
| EP3028071B1 (en) | Method and device for the generation and application of anisotropic elastic parameters in horizontal transverse isotropic (hti) media | |
| RU2523734C2 (en) | Seismic data collection system and method | |
| EP2335093B1 (en) | Estimation of soil properties using waveforms of seismic surface waves | |
| US9081115B2 (en) | Convergence rate of full wavefield inversion using spectral shaping | |
| US8296069B2 (en) | Pseudo-analytical method for the solution of wave equations | |
| EP3063562B1 (en) | Methods of subsurface exploration, computer program product and computer-readable storage medium | |
| CN107272058B (en) | Imaging method, imaging apparatus, and computer storage medium | |
| EP2497043B1 (en) | Seismic imaging systems and methods employing a 3d reverse time migration with tilted transverse isotropy | |
| CN112327358B (en) | Forward modeling method for acoustic seismic data in viscous medium | |
| US11327195B2 (en) | Correction of source motion effects in seismic data recorded in a marine survey using a moving source | |
| CN103149585A (en) | Elastic migration seismic wave field construction method and elastic migration seismic wave field construction device | |
| MX2011003850A (en) | Image domain signal to noise estimate. | |
| CN113341455B (en) | Viscous anisotropic medium seismic wave numerical simulation method, device and equipment | |
| CN110542928A (en) | Seismic response simulation method based on VTI anisotropic propagation matrix | |
| CN111665556B (en) | Stratum acoustic wave propagation velocity model construction method | |
| CN112285778B (en) | Reverse time migration imaging method for pure qP waves in sticky sound TTI medium | |
| US11199641B2 (en) | Seismic modeling | |
| CN109239776A (en) | A kind of seimic wave propagation the Forward Modeling and device | |
| CN111505714A (en) | Elastic wave direct envelope inversion method based on rock physical constraint | |
| CN113866823A (en) | Forward imaging method in visco-acoustic anisotropic medium | |
| Wapenaar et al. | The wavelet transform as a tool for geophysical data integration | |
| MX2011003852A (en) | Time reverse imaging attributes. | |
| CN115704913A (en) | Shale reservoir stratum crustal stress prediction method and device | |
| CN117272717A (en) | Source numerical simulation method suitable for forward modeling of elastic wave field |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |