CN114114403B - Anisotropic attenuation medium simulation method based on fractional order Laplace operator - Google Patents
Anisotropic attenuation medium simulation method based on fractional order Laplace operator Download PDFInfo
- Publication number
- CN114114403B CN114114403B CN202111578115.5A CN202111578115A CN114114403B CN 114114403 B CN114114403 B CN 114114403B CN 202111578115 A CN202111578115 A CN 202111578115A CN 114114403 B CN114114403 B CN 114114403B
- Authority
- CN
- China
- Prior art keywords
- medium
- wave
- elastic
- vti
- equation
- 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
- 238000004088 simulation Methods 0.000 title claims abstract description 35
- 238000000034 method Methods 0.000 title claims abstract description 30
- 238000013016 damping Methods 0.000 claims abstract description 26
- 230000014509 gene expression Effects 0.000 claims abstract description 25
- 239000006185 dispersion Substances 0.000 claims abstract description 23
- 239000002245 particle Substances 0.000 claims description 3
- 150000001875 compounds Chemical class 0.000 claims 2
- 238000012360 testing method Methods 0.000 description 15
- 238000010586 diagram Methods 0.000 description 12
- 239000000243 solution Substances 0.000 description 9
- 238000005070 sampling Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 4
- 239000011159 matrix material Substances 0.000 description 4
- 230000003595 spectral effect Effects 0.000 description 4
- 230000015572 biosynthetic process Effects 0.000 description 3
- 230000008859 change Effects 0.000 description 3
- 230000006870 function Effects 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 239000012088 reference solution Substances 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 238000012952 Resampling Methods 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000004590 computer program Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000004807 localization Effects 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000001615 p wave Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000004613 tight binding model Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/40—Transforming data representation
- G01V2210/48—Other transforms
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/50—Corrections or adjustments related to wave propagation
- G01V2210/58—Media-related
- G01V2210/584—Attenuation
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/50—Corrections or adjustments related to wave propagation
- G01V2210/58—Media-related
- G01V2210/586—Anisotropic media
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
The application provides an anisotropic attenuation medium simulation method based on fractional order Laplace operator, which comprises the following steps: constructing a dispersion relation of elastic waves based on a preset dispersion relation, wherein the elastic waves comprise P waves and S waves; obtaining the complex modulus of the elastic damping medium; acquiring an expression of constitutive relation of the VTI attenuation medium in a frequency-wave number domain; expanding complex modulus of the obtained elastic damping medium to VTI damping anisotropy, and converting an expression of constitutive relation of the VTI damping medium in a frequency-wave number domain into an expression in a time-space domain; based on an expression, a geometric equation and a motion balance equation of the constitutive relation of the VTI attenuation medium in a time-space domain, obtaining a DFL viscous-elastic equation of the VTI medium; s150, determining parameters of the underground medium, and substituting the parameters into a viscous-elastic equation of a DFL (distributed-energy-density-distribution-coefficient) medium to calculate and obtain a seismic wave field simulation value; the invention can enable the earthquake simulation to be more accurate.
Description
Technical Field
The application relates to the field of oil gas detection, in particular to an anisotropic attenuation medium simulation method based on fractional order Laplace operator.
Background
In oil and gas exploration and development, the seismic exploration technology plays an extremely important role. At present, the development of seismic wave field numerical simulations for complex media has become a consensus in the industry today, and fractional Decoupling (DFL) viscous wave equations developed in recent years have a number of unique advantages. However, the numerical solution of the DFL viscous wave equation often faces challenges. The Laplacian, which contains orders that vary spatially, presents difficulties in the wavefield simulation of non-uniform media. Currently, although anisotropic viscoelastic waves and acoustic wave equations exist in the form of DFLs, the hybrid domain laplace operator contained therein requires additional processing to be suitable for modeling of non-uniform media.
Patent document (CN 113341455 a) discloses a numerical simulation method of seismic waves of viscous multidirectional medium, and the technical scheme can accurately describe speed and attenuation anisotropy at the same time and is convenient for numerical calculation equation, but only takes attenuation of longitudinal waves and anisotropy of speed into consideration, so that the problem of non-conforming to actual underground conditions exists, and the seismic simulation is inaccurate, because the seismic waves are transmitted in the underground in longitudinal waves and transverse waves at the same time.
Disclosure of Invention
Aiming at the technical problems, the application provides the anisotropic attenuation medium simulation method based on the fractional order Laplace operator, and meanwhile, the anisotropy of longitudinal waves and transverse waves is considered, so that the method can better fit the underground actual situation, and the earthquake simulation is more accurate.
The technical scheme adopted by the embodiment of the invention is as follows:
the embodiment of the invention provides an anisotropic attenuation medium simulation method based on fractional order Laplace operator, which comprises the following steps:
s100, constructing a dispersion relation of elastic waves based on a preset dispersion relation, wherein the elastic waves comprise P waves and S waves;
s110, obtaining complex modulus of the elastic attenuation medium based on the dispersion relation of the constructed elastic waves;
s120, acquiring an expression of constitutive relation of the VTI attenuation medium in a frequency-wave number domain;
s130, expanding complex modulus of the obtained elastic damping medium to VTI damping anisotropy, and converting an expression of constitutive relation of the VTI damping medium in a frequency-wave number domain into an expression in a time-space domain;
s140, obtaining a DFL viscous elastic equation of the VTI medium based on an expression, a geometric equation and a motion balance equation of the constitutive relation of the VTI attenuation medium in a time-space domain;
determining underground medium parameters, substituting the underground medium parameters and seismic wave parameters into a VTI medium DFL viscous-elastic equation to calculate and obtain a seismic wave field simulation value;
wherein, the VTI medium DFL viscous elastic equation is:
wherein t represents time, v x And v z Vibration velocity components of particles of an underground medium along the x direction and the z direction respectively, wherein rho is the density and sigma of the underground medium xx 、σ zz Sum sigma xz A horizontal component, a longitudinal component and a shear component of the stress tensor, respectively;
The invention has at least the following technical effects: the anisotropic attenuation medium simulation method based on fractional order Laplace operator provided by the invention considers the anisotropism of longitudinal waves and transverse waves, can better fit the underground actual condition, and can more accurately describe the attenuation mechanism of seismic waves in the underground propagation process.
Drawings
In order to more clearly illustrate the technical solutions of the embodiments of the present application, the drawings that are needed in the description of the embodiments will be briefly introduced below, and it is obvious that the drawings in the following description are only some embodiments of the present application, and that other drawings may be obtained according to these drawings without inventive effort for a person skilled in the art.
FIG. 1 is a flow chart of an anisotropic attenuation medium simulation method based on fractional order Laplace operator provided in an embodiment of the present application;
FIGS. 2 (a) and 2 (b) are schematic diagrams for verifying the simulation accuracy of the VTI medium DFL viscoelasticity equation provided by an embodiment of the present application by comparison with an analytical solution;
FIG. 3 is a schematic diagram for verifying the decoupling characteristics of the DFL viscoelastic equation of the VTI medium according to an embodiment of the present application;
FIGS. 4 a-c are three layered models for examining the ability of the VTI media DFL viscoelasticity equation provided by an embodiment of the present application to handle Q-value and velocity strongly varying media;
FIG. 5 is a diagram of an observation system for estimating Q value by using a spectral ratio method according to an embodiment of the present application;
FIG. 6 is a schematic diagram of a common shot gather recorded by an observation system used in the spectral ratio method according to the embodiment of the present application;
FIG. 7 is a reciprocal comparison chart of quality factors provided in embodiments of the present application;
fig. 8a to 8d are schematic diagrams of 2007BP TTI models according to embodiments of the present application;
FIG. 9a is a snapshot of wavefield in an anisotropic elastic media provided in an embodiment of the present application;
fig. 9b and 9c are schematic diagrams of simulation results of an anisotropic DFL viscoelastic wave equation and an anisotropic time fractional order viscoelastic wave equation, respectively;
FIG. 9d is the difference between FIG. 9b and FIG. 9 c;
FIG. 10a is a schematic diagram of seismic recording in an anisotropic elastic media provided in an embodiment of the present application;
fig. 10b and 10c are schematic diagrams showing simulation results of an anisotropic DFL viscoelastic wave equation and an anisotropic time fractional order viscoelastic wave equation, respectively;
FIG. 10d is the difference between FIG. 10b and FIG. 10 c;
FIG. 11a is a position diagram of a Sorton trough model;
FIG. 11b is a schematic diagram of a resampled P-wave velocity model;
FIGS. 11c and 11d are schematic cross-sectional views of a Sortton groove model with different latitudes and different longitudes, respectively;
FIGS. 12 a-12 f are wave field snap graphs of DFL viscoelastic wave equation calculations for different moments provided in embodiments of the present application;
FIGS. 13a-c are schematic diagrams of calculating a three-dimensional common shot gather using the viscoelastic wave equation provided by the embodiments of the present application, the existing time fractional order viscoelastic wave equation elastic wave equation;
FIG. 14 shows a viscoelastic wave equation and v provided in an embodiment of the present application x Component, v y Component v z And a schematic diagram for tracking and comparing the video envelope objective function, the time fractional order viscoelastic wave equation and the time-frequency phase error in the seismic record of the component.
Detailed Description
The following description of the embodiments of the present application will be made clearly and fully with reference to the accompanying drawings, in which it is evident that the embodiments described are only some, but not all, of the embodiments of the present application. All other embodiments, which can be made by those skilled in the art based on the embodiments herein without making any inventive effort, are intended to be within the scope of the present application.
As shown in fig. 1, the anisotropic attenuation medium simulation method based on fractional order Laplace operator provided by the embodiment of the invention comprises the following steps:
s100, constructing a dispersion relation of elastic waves based on a preset dispersion relation, wherein the elastic waves comprise P waves and S waves.
In the embodiment of the present invention, the preset dispersion relation is an approximate dispersion relation of a constant Q model. The approximate dispersion relation of the constant Q model represents the change relation of frequency and speed, and the formula is as follows:
ω 2 ＝d 1 k+d 2 k 2 +d 3 k 3 +d 4 (iω)k+d 5 (iω)k 2 (1)
(d 1 ,d 2 ,d 3 ,d 4 ,d 5 )＝(-γvω 0 ,v 2 ,γv 3 ω 0 -1 ,πγv,πγ 2 v 2 ω 0 -1 ) (2)
where k is the wave number, ω is the angular frequency, ω 0 For reference angular frequency, γ is the space-variant fractional order, γ=arctanq -1 Pi, Q is the quality factor of the elastic wave; v is the reference speed, v=v 0 cos(πγ/2)，v 0 For at reference frequency omega 0 Reference phase velocity at.
The elastic wave contains both P wave and S wave components, and both satisfy the approximate dispersion relation, so the dispersion relation of the constructed constant Q model of the elastic wave is:
wherein v is α Is the velocity of elastic wave, gamma α ＝arctanQ α -1 /π，Q α α=p or S, which is the quality factor of the elastic wave.
S110, obtaining the complex modulus of the elastic attenuation medium based on the dispersion relation of the constructed elastic waves.
In an embodiment of the present invention, the complex modulus of the elastic damping medium may be:
s120, acquiring an expression of constitutive relation of the VTI attenuation medium in a frequency-wave number domain.
In the embodiments of the present application, consider a simpler anisotropic medium: a transversely isotropic medium (VTI) having a vertical symmetry axis. The stress-strain relationship is determined by an elastic stiffness matrix closely related to the elastic property of the medium, and for the VTI elastic medium, the expression of the elastic stiffness matrix in the form of Vogit is
The elastic rigidity coefficient is:
wherein ρ represents the density of the subsurface medium; v p Representing the speed of the P wave along the direction of the symmetry axis; v s Is the velocity of the S wave, epsilon and delta being the velocity anisotropy parameters.
For a VTI attenuation medium, its Q matrix may be expressed as
Wherein Q is 11 And Q 33 Corresponding to the quality factors of the P wave in the transverse direction and the vertical direction respectively, and Q 55 Quality factor corresponding to S wave, Q 12 ＝Q 13 ,Q 55 ＝Q 66 . Furthermore, the decaying anisotropy parameter ε, which is similar to the Thomsen coefficient, is used Q And delta Q Thereby establishing Q 11 ，Q 33 And Q 13 The relationship between the two, specifically,
unlike the anisotropic elastic medium whose stiffness coefficient is a real number, the anisotropic attenuation medium corresponds to a complex stiffness coefficient (or complex modulus, noted as
) It can be represented by the stiffness matrix as followsRelationship conversion
The constitutive relation of the VTI attenuation medium is expressed as the following in the frequency-wave number domain
Wherein,,
and->
Expression in frequency of horizontal, longitudinal and shear components of stress tensor, respectively,/->
And->
Expression in frequency of horizontal, longitudinal and shear components, respectively, of strain tensor,/->
And->
The complex moduli of the anisotropic damping medium respectively.
S130, expanding complex modulus of the obtained elastic damping medium to VTI damping anisotropy, and converting an expression of constitutive relation of the VTI damping medium in a frequency-wave number domain into an expression in a time-space domain.
The method specifically comprises the following steps:
s131, decomposing the formula (3) into P wave and S wave components to obtain corresponding components
And->
Wherein,,
s132, because the formula (14) is equal to v only 1 In relation, therefore, from equation (4), the equation (14) can be derived
In particular, will v 1 V in alternative (4) α Obtaining +.about.in formula (14)>
Is expressed as:
s133 because the formula (15) is equal to v only s In relation, therefore, the formula (15) can be deduced from the formula (4)
And->
In particular, will v s V in alternative (1) α Obtaining +.about.in formula (15)>
And->
Is expressed as:
s134, substituting formula (16) into formula (14) and substituting formulas (17) and (18) into formula (15), to obtain an expression of formula (11) in the time-space domain:
s135, operations similar to S131 to S134 are performed on the formulas (12) and (13), respectively, that is, decomposed into P-wave and S-wave components, respectively, and then expressions of the formulas (12) and (13) in the time-space domain are obtained, respectively, based on the formula (4):
σ xz ＝L 55 e xz (21)
wherein e xx 、e zz And e xz A horizontal component, a longitudinal component and a shear component of the strain tensor, respectively.
S140, obtaining a DFL viscous-elastic equation of the VTI medium based on an expression, a geometric equation and a motion balance equation of the constitutive relation of the VTI attenuation medium in a time-space domain.
In the embodiment of the invention, the geometric equation and the motion balance equation are as follows:
where i and j are x or z.
In the embodiment of the invention, the VTI medium DFL viscoelastic equation based on the new dispersion relation obtained by combining the formulas (19) to (22) is:
wherein t represents time, v x And v z Vibration velocity components, sigma, of particles of an underground medium along x direction and z direction respectively xx 、σ zz Sum sigma xz A horizontal component, a longitudinal component and a shear component of the stress tensor, respectively;
S150, determining underground medium parameters, and substituting the underground medium parameters into the VTI medium DFL viscous elasticity equation to calculate and obtain a seismic wave field simulation value.
On the basis of obtaining a VTI medium DFL viscous-elastic equation, the underground medium parameters including underground medium density, P wave and S wave quality factors, P wave speed and S wave speed, attenuation anisotropy parameters, speed anisotropy parameters and the like are determined, and the parameters are substituted into the VTI medium DFL viscous-elastic equation, so that the simulation of a seismic wave field is realized, and a seismic wave field simulation value is obtained.
Test example 1
In this test example, the accuracy of the DFL viscoelastic wave equation shown in the above equation (23) was first verified by comparison with an analytical solution in which the analysis was performed in a uniform mediumThe solution may be derived from a green function. The uniform media model size used was 512 x 512 grid points with a spatial sampling interval of 10m. The main frequency of the Rake wavelet source arranged in the center of the uniform medium is 25Hz, and the reference frequencies are respectively c with the same reference speed as the source 0p =3000 m/s and c 0s =2000 m/s, medium density 2200 kg/m 3 . V at record (3160 m ) x The components, part 2a and part 2b correspond to the fully elastic medium and the damping medium (Q p ＝32，Q s =20). In the figure, the thick line, the dotted line, and the thin black line correspond to the analytical solution, the numerical solution, and the two residuals, respectively. It can be found that the proposed new equation fits very well with the analytical solution, both for the damping medium and the completely elastic medium, proving the accuracy of the DFL viscoelastic wave equation.
Test example 2
In this test example, the decoupling characteristics of the DFL viscoelastic wave equation were further verified. In fig. 3, the four wavefield snapshots correspond to the full elastic wave equation, the damping master equation, the dispersion master equation, and the proposed DFL viscoelastic wave equation, respectively. It can be observed that the damping master equation describes well the damping of amplitude energy compared to elastic waves, the dispersion master equation describes the velocity dispersion (phase lag observable) caused by the formation damping, while the DFL viscoelastic wave equation simulates the wave field with both effects.
Test example 3
In this test example, the ability of the DFL viscoelastic wave equation to handle non-uniform attenuation media was further verified. In this test example, three layered models were designed to examine their ability to handle medium of dramatic changes in Q value and velocity. As shown in fig. 4 (a), the velocity in the model 1 is a laminar medium and the Q value is uniform, the velocity in the model 2 is uniform and the Q value is laminar, and both the velocity and Q in the model 3 are laminar. The model size of the group is 800×800 grid points, the space sampling interval is 10m, the horizontal interface is positioned at the depth 4400m, the adopted time step is 1ms, and the main frequency of the Rake wavelet source arranged at the center of the medium is 25Hz. However, in case of significant Q contrast (model 2 and model 3), the snapshot residual between schemes is more pronounced than the reference snapshot, which is the result of the phase distortion superposition effect (see fig. 4 c). The third column in fig. 4b has only subtle residuals and the two lines in fig. 4c match well, which demonstrates the good ability of the DFL viscoelastic wave equation provided by the present invention to handle non-uniform attenuation media.
Test example 4
In this test case, a set of anisotropic homogeneous model tests were performed to verify the accuracy of the proposed VTI-DFL viscoelasticity equation describing the decay as a function of direction. Firstly, solving a proposed VTI-DFL viscous elastic equation to obtain a corresponding snapshot. On one hand, the Q value of each direction can be estimated by using a spectrum ratio method; on the other hand, by solving the Christoffel equation of the plane wave equation, the theoretical Q value can be obtained. The accuracy of the VTI-DFL viscoelasticity equation provided by the invention is evaluated by comparing the compatibility of the theoretical Q value with the estimated Q value. Fig. 5 is a simulation geometry in which a star is a rake wavelet with a dominant frequency of 25Hz and two concentric circles of triangles are receivers, each angularly spaced 600m apart. The test used a set of uniform anisotropic media models whose decay anisotropy parameters are shown in table 1. The model size is 400 x 400 grid points, the spatial sampling interval is 10m, and the time step is 1ms. The remaining parameters are v p ＝3.0km/s，v s The Thomsen velocity anisotropy parameter was epsilon=0.16, delta=0.08, medium density 2000kg/m3, =1.5 km/s. To quantify the intensity of the medium decay with direction, the intensity coefficient Strength (Q p,s )＝1-min(Q p,s )/max(Q p,s )。
TABLE 1 elastic damping medium model for damping anisotropic parameter variation
Common shot gathers recorded using the observation system used for spectral ratio method can be shown in fig. 6, wherein in fig. 6, the left side is the inner concentric circle detector point and the right side is the outer concentric circle detector point. The reciprocal comparison of the obtained quality factors can be shown in fig. 7, wherein circles represent the Q values estimated by the spectral ratio method, and black solid lines represent theoretical Q values obtained by Christoffel equations.
Test example 5
In this test example, the simulation accuracy of the proposed VTI-DFL viscoelasticity equation was verified with a complex model, which was 2007BP TTI model (without considering tilt angle). FIGS. 8a to 8d show the P-wave velocity and the P-wave quality factor Q, respectively p Thomsen velocity anisotropy parameters δ and ε. S-wave velocity and quality factor Q s Respectively through v s ＝v p /1.73，Q s ＝Q p Converted to 1.2, the damping anisotropy parameter is delta Q ＝2δ，ε Q =3ε. The model has the size of 500×270 grid points, a spatial sampling interval of 10m, a time sampling interval of 1ms and a simulation propagation time of 2.5s, the adopted excitation source is a Rake wavelet with a main frequency of 25Hz, the source position is (2500 m,30 m), and the detector point position is positioned at each grid position with a depth of 50m along the horizontal direction. To evaluate simulation accuracy, we calculated a reference solution from the existing VTI-TF viscoelastic wave equation (Zhu, 2017) that is capable of describing the direction dependent attenuation, but requires significant computational and memory overhead. FIG. 9a shows a snapshot of the wavefield in an anisotropic elastic medium, FIG. 9b a snapshot of the wavefield simulated by the VTI-DFL viscoelastic equation, FIG. 9c a snapshot of the wavefield simulated by the VTI medium time fractional order viscoelastic wave equation, which does not need to be approximated when expanding to anisotropy, so the variation of formation attenuation with direction can be well described, but because of the computational and memory problems of the time fractional operators contained therein during the numerical solution, it is not suitable for large scale numerical simulation; and comparing with the simulation result of the equation, the simulation precision of the proposed VTI medium DFL viscous equation can be verified. Fig. 9d is the result of the subtraction of fig. 9b and 9 c. Fig. 10 a-10 d are common shot gather seismic records corresponding to fig. 9 a-9 d, respectively. As can be seen from the comparison of the analysis,the formation attenuation effect significantly attenuates amplitude energy as compared to the non-attenuating medium; in addition, the weak difference between fig. 9d and 10d shows that the proposed VTI-DFL viscoelastic equation still has high numerical simulation accuracy for complex models.
Test example 6
In this test example, to demonstrate the versatility and feasibility of the proposed DFL viscoelastic wave equation, the propagation of seismic waves was simulated in a three-dimensional real model of the solton sea chest, south california. The localization map of the solton sea chest is shown in fig. 11a, where the initial P-wave velocity model is provided by (Ajala, 2019) which contains 19 grids within a depth of 0 to 9 km, 101 grids within a longitude 115.7 ° W to 116.7 ° W, and 91 grids within a latitude of 33.3 ° N34.2 °. Since the raw data is not uniform, linear interpolation and resampling can result in a new p-wave velocity model, comprising 505×455×48 grid points with a grid pitch uniformity of 200m (see fig. 11 b). Fig. 11c and 11d show some sections of different latitudes and longitudes, respectively. S wave velocity is v s ＝v p A velocity model of 1.73 (Li, 1993), Q model is obtained from Q P ＝Q S ＝3.516(v p /1000) 22 (Li, 1993). At (50000 m,46000m, 1600 m) a Rake wavelet source with a dominant frequency of 5hz is placed and a receiver is located at each grid point in the x-y plane of the earth's surface. The time is 15s and the time step is 5.0ms. Fig. 12a to 12f show snapshots of components simulated at different times using the proposed DFL viscoelastic wave equation. It is particularly noted that some energy is observed to be trapped in the elliptical marked area where the seismic waves travel slower than in the surrounding area. The reason for this phenomenon can be found in the Sorton sea-trough model, and it can be seen from FIG. 11b that there are some low-velocity, high-attenuation deposits in the shallow portion of the zone. In order to evaluate the accuracy of the proposed solton sea-groove model DFL viscoelasticity equation, a time fractional order viscoelasticity wave equation was used as a reference solution. Since computing the time-fractional order viscoelastic wave equation requires a lot of computer memory, the model is resampled to 252 x 227 x 24 grid points, with a grid spacing of 20m for simulation to be feasible. Seismic sourceThe dominant frequency is raised to 25Hz and the time step is reduced to 1.0ms. Fig. 13 shows a 3D co-shot gather, wherein the horizontal slice shows a snapshot at 1.1s and the two vertical sections show two co-shot gathers in x and y directions. Fig. 13a-c correspond to the DFL viscoelastic wave equation, the time fractional order viscoelastic wave equation, and the elastic wave equation, respectively, proposed by the present invention. The common shot gathers in fig. 13a and 13b look very similar, their amplitude being significantly weaker than the elastic case due to energy decay (fig. 13 c). For better contrast, the seismic recordings (1500 m ) are further shown in FIG. 14. Columns from left to right correspond to v x ，v y And v z Is a component of (a). The upper and lower rows represent the time-frequency envelope objective function (TFEM) and the time-frequency phase error (TFPM, kristekova et al, 2006, 2009), respectively. The seismic records of the three components (middle row of fig. 14) match well, with TFEM and TFPM within 2%, indicating that the proposed DFL viscoelastic wave equation generates a seismic record that is nearly identical to the reference solution.
In summary, the method provided by the embodiment of the application provides the acoustic wave equation based on the fractional Laplace operator, can accurately describe the characteristics of the speed and the attenuation of the seismic wave along with the change of the azimuth at the same time, does not contain the fractional Laplace operator along with the change of the space, so that the propagation of the seismic wave in a non-uniform medium can be naturally and accurately calculated, and the dimension reduction processing (such as three-dimensional to two-dimensional) of the acoustic wave equation of the anisotropic medium can be realized through the simplified processing of the acoustic wave equation, thereby being widely applied to the simulation of the acoustic wave field, the absorption attenuation compensation reverse time migration imaging and the full waveform inversion of the acoustic wave field in the two-dimensional (2D) and three-dimensional (3D) complex attenuation anisotropic medium. In addition, due to the fact that anisotropy of longitudinal waves and transverse waves is considered, the underground actual situation can be well attached, and an attenuation mechanism of the seismic waves in the underground propagation process can be accurately described.
Embodiments of the present application also provide a non-transitory computer readable storage medium that may be disposed in an electronic device to store at least one instruction or at least one program for implementing one of the methods embodiments, the at least one instruction or the at least one program being loaded and executed by the processor to implement the methods provided by the embodiments described above.
Embodiments of the present application also provide an electronic device comprising a processor and the aforementioned non-transitory computer-readable storage medium.
Embodiments of the present application also provide a computer program product comprising program code for causing an electronic device to carry out the steps of the method according to various exemplary embodiments of the present application as described in the present specification, when said program product is run on the electronic device.
Although some specific embodiments of the present application have been described in detail by way of example, it should be understood by those skilled in the art that the above examples are for illustration only and are not intended to limit the scope of the present application. Those skilled in the art will also appreciate that various modifications might be made to the embodiments without departing from the scope and spirit of the present application. The scope of the present disclosure is defined by the appended claims.
Claims (8)
1. An anisotropic attenuation medium simulation method based on fractional order Laplace operator is characterized by comprising the following steps:
s100, constructing a dispersion relation of elastic waves based on a preset dispersion relation, wherein the elastic waves comprise P waves and S waves;
s110, obtaining complex modulus of the elastic attenuation medium based on the dispersion relation of the constructed elastic waves;
s120, acquiring an expression of constitutive relation of the VTI attenuation medium in a frequency-wave number domain;
s130, expanding complex modulus of the obtained elastic damping medium to VTI damping anisotropy, and converting an expression of constitutive relation of the VTI damping medium in a frequency-wave number domain into an expression in a time-space domain;
s140, obtaining a DFL viscous elastic equation of the VTI medium based on an expression, a geometric equation and a motion balance equation of the constitutive relation of the VTI attenuation medium in a time-space domain;
s150, determining underground medium parameters, and substituting the underground medium parameters into the VTI medium DFL viscous elastic equation to calculate and obtain a seismic wave field simulation value;
wherein, the VTI medium DFL viscous elastic equation is:
wherein t represents time, v x And v z Vibration velocity components of particles of an underground medium along the x direction and the z direction respectively, wherein rho is the density and sigma of the underground medium xx 、σ zz Sum sigma xz A horizontal component, a longitudinal component and a shear component of the stress tensor, respectively;
2. The method of claim 1, wherein the predetermined dispersion relation is an approximate dispersion relation of a constant Q model;
the dispersion relation of the constructed elastic wave is as follows:
wherein k is the wave number, ω is the angular frequency, v α Is the velocity of elastic wave, gamma α ＝arctanQ α -1 /π，Q α α=p or S, which is the quality factor of the elastic wave.
3. The method of claim 2, wherein the elastic damping medium has a complex modulus of
4. A method according to claim 3, wherein the constitutive relation of the VTI attenuation medium in the frequency-wavenumber domain is expressed as:
wherein,,
and->
Expression in frequency of horizontal, longitudinal and shear components of stress tensor, respectively,/->
And->
The expressions in frequency for the horizontal component, the longitudinal component and the shear component of the strain tensor respectively,
and->
The complex moduli of the anisotropic damping medium respectively.
5. The method of claim 4, wherein S130 further comprises:
s131, decomposing the formula (3) into P wave and S wave components to obtain corresponding components
And->
Wherein,,
s132, v 1 V in alternative (1) α Obtaining the compound of formula (6)
Is expressed as:
s133, v s V in alternative (1) α Obtaining the compound of formula (7)
And->
Is expressed as:
s134, substituting formula (8) into formula (6) and substituting formulas (9) and (10) into formula (7), to obtain an expression of formula (3) in the time-space domain:
s135, performing operations similar to S131 to S134 on the formulas (4) and (5), respectively, to obtain expressions of the formulas (4) and (5) in the time-space domain, respectively:
σ xz ＝L 55 e xz (13)
wherein e xx 、e zz And e xz A horizontal component, a longitudinal component and a shear component of the strain tensor, respectively.
6. The method of claim 5, wherein the geometric equation and the motion balance equation are:
where i and j are x or z.
7. A non-transitory computer readable storage medium having stored therein at least one instruction or at least one program loaded and executed by a processor to implement the method of any one of claims 1-6.
8. An electronic device comprising a processor and the non-transitory computer-readable storage medium of claim 7.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111578115.5A CN114114403B (en) | 2021-12-22 | 2021-12-22 | Anisotropic attenuation medium simulation method based on fractional order Laplace operator |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111578115.5A CN114114403B (en) | 2021-12-22 | 2021-12-22 | Anisotropic attenuation medium simulation method based on fractional order Laplace operator |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114114403A CN114114403A (en) | 2022-03-01 |
| CN114114403B true CN114114403B (en) | 2023-06-27 |
Family
ID=80361940
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111578115.5A Active CN114114403B (en) | 2021-12-22 | 2021-12-22 | Anisotropic attenuation medium simulation method based on fractional order Laplace operator |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114114403B (en) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115392090B (en) * | 2022-09-13 | 2023-05-23 | 中国矿业大学 | Prediction method of seismic wave frequency dispersion and attenuation characteristics based on three-dimensional digital rock core |
| CN117471531B (en) * | 2023-11-01 | 2024-04-19 | 中国矿业大学 | VTI viscoelastic wave equation numerical simulation method based on power law frequency-varying Q effect |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102508296A (en) * | 2011-11-14 | 2012-06-20 | 中国石油天然气股份有限公司 | Unsaturated double-pore medium seismic wave frequency dispersion attenuation analysis method and device |
| CN102590859A (en) * | 2011-12-31 | 2012-07-18 | 中国石油集团西北地质研究所 | Anisotropic reverse time migration method for quasi-P wave equation in transverse isotropy with a vertical axis of symmetry (VTI) medium |
| CN105137486A (en) * | 2015-09-01 | 2015-12-09 | 中国科学院地质与地球物理研究所 | Elastic wave reverse-time migration imaging method and apparatus in anisotropic media |
Family Cites Families (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7930926B2 (en) * | 2007-05-01 | 2011-04-26 | Boise State University | Determination of permeability from damping |
| US7924652B2 (en) * | 2007-06-01 | 2011-04-12 | Baker Hughes Incorporated | Method for determining seismic anisotropy |
| WO2011112932A1 (en) * | 2010-03-12 | 2011-09-15 | Cggveritas Services (Us) Inc. | Methods and systems for performing azimuthal simultaneous elatic inversion |
| CN104570084B (en) * | 2013-10-29 | 2018-01-05 | 中国石油化工股份有限公司 | Across yardstick earthquake rock physicses attenuation model and the method for prediction decay and frequency dispersion |
| WO2015199800A1 (en) * | 2014-06-17 | 2015-12-30 | Exxonmobil Upstream Research Company | Fast viscoacoustic and viscoelastic full-wavefield inversion |
| CN104360383B (en) * | 2014-11-12 | 2017-03-22 | 中国石油大学(华东) | Method and system for predicting seismic wave attenuation |
| 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 |
| WO2018004789A1 (en) * | 2016-06-28 | 2018-01-04 | Exxonbil Upstream Research Company | Reverse time migration in anisotropic media with stable attenuation compensation |
| CN110542928A (en) * | 2018-05-28 | 2019-12-06 | 中国石油化工股份有限公司 | Seismic response simulation method based on VTI anisotropic propagation matrix |
| CN110058310B (en) * | 2019-04-16 | 2021-07-09 | 中国石油大学(华东) | Seismic wave velocity frequency dispersion and attenuation prediction method and device for seismic rock |
| US11360224B2 (en) * | 2019-05-03 | 2022-06-14 | Exxonmobil Upstream Research Company | Inversion, migration, and imaging related to isotropic wave-mode- independent attenuation |
| CN113189642B (en) * | 2021-04-28 | 2024-04-05 | 中国石油化工集团有限公司 | Seismic source linear scanning signal design method based on forced vibration |
| CN113341455B (en) * | 2021-06-24 | 2024-02-09 | 中国石油大学(北京) | Viscous anisotropic medium seismic wave numerical simulation method, device and equipment |
-
2021
- 2021-12-22 CN CN202111578115.5A patent/CN114114403B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102508296A (en) * | 2011-11-14 | 2012-06-20 | 中国石油天然气股份有限公司 | Unsaturated double-pore medium seismic wave frequency dispersion attenuation analysis method and device |
| CN102590859A (en) * | 2011-12-31 | 2012-07-18 | 中国石油集团西北地质研究所 | Anisotropic reverse time migration method for quasi-P wave equation in transverse isotropy with a vertical axis of symmetry (VTI) medium |
| CN105137486A (en) * | 2015-09-01 | 2015-12-09 | 中国科学院地质与地球物理研究所 | Elastic wave reverse-time migration imaging method and apparatus in anisotropic media |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114114403A (en) | 2022-03-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114114403B (en) | Anisotropic attenuation medium simulation method based on fractional order Laplace operator | |
| Song et al. | Modeling of pseudoacoustic P-waves in orthorhombic media with a low-rank approximation | |
| US8560242B2 (en) | Pseudo-analytical method for the solution of wave equations | |
| EP2542920B1 (en) | Method for local attribute matching in seismic processing | |
| US20120314538A1 (en) | System and method for seismic data inversion | |
| AU2012268720A1 (en) | System and method for data inversion with phase extrapolation | |
| Withers et al. | Memory-efficient simulation of frequency-dependent Q | |
| US9753166B2 (en) | P-wave and S-wave separation of seismic data in the presence of statics and irregular geometry | |
| CN113341455B (en) | Viscous anisotropic medium seismic wave numerical simulation method, device and equipment | |
| MX2011003850A (en) | Image domain signal to noise estimate. | |
| CN110286410B (en) | Fracture inversion method and device based on diffracted wave energy | |
| Johnson | Source mechanisms of induced earthquakes at The Geysers geothermal reservoir | |
| CN105891884A (en) | Micro-earthquake focus mechanism inversion method and micro-earthquake focus mechanism inversion device | |
| NO20150821A1 (en) | Efficient Wavefield Extrapolation in Anisotropic Media | |
| CN114139335A (en) | Viscous sound wave simulation method based on single relaxation time lattice Boltzmann model | |
| Raknes et al. | A numerical study of 3D elastic time-lapse full-waveform inversion using multicomponent seismic data | |
| CN109541682A (en) | Isotropic elasticity parameter protects width inversion method and device | |
| CN109946742A (en) | The pure rolling land qP shakes digital simulation method in a kind of TTI medium | |
| CN105425298A (en) | Method and device for eliminating numerical frequency dispersion in finite difference forward process | |
| US20160018541A1 (en) | System and method for rock property estimation of subsurface geologic volumes | |
| Barnier et al. | Full-waveform inversion by model extension: Practical applications | |
| US11199641B2 (en) | Seismic modeling | |
| Chang et al. | Elastic full waveform inversion for salt model building in Gulf of Mexico | |
| CN109738944B (en) | Wide-angle reflection-based seismic acquisition parameter determination method and device | |
| De Matteis et al. | BISTROP: Bayesian inversion of spectral‐level ratios and P‐wave polarities for focal mechanism determination |
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 |