CN114114403A - Fractional order Laplace operator-based anisotropic attenuation medium simulation method - Google Patents

Fractional order Laplace operator-based anisotropic attenuation medium simulation method Download PDF

Info

Publication number
CN114114403A
CN114114403A CN202111578115.5A CN202111578115A CN114114403A CN 114114403 A CN114114403 A CN 114114403A CN 202111578115 A CN202111578115 A CN 202111578115A CN 114114403 A CN114114403 A CN 114114403A
Authority
CN
China
Prior art keywords
wave
medium
vti
equation
attenuation
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.)
Granted
Application number
CN202111578115.5A
Other languages
Chinese (zh)
Other versions
CN114114403B (en
Inventor
王宁
石颖
王维红
赵海波
郭雪豹
王子墨
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northeast Petroleum University
Original Assignee
Northeast Petroleum University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Northeast Petroleum University filed Critical Northeast Petroleum University
Priority to CN202111578115.5A priority Critical patent/CN114114403B/en
Publication of CN114114403A publication Critical patent/CN114114403A/en
Application granted granted Critical
Publication of CN114114403B publication Critical patent/CN114114403B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

The application provides an anisotropic attenuation medium simulation method based on fractional order Laplace operators, which comprises the following steps: constructing a frequency dispersion relation of an elastic wave based on a preset frequency dispersion relation, wherein the elastic wave comprises a P wave and an S wave; obtaining the complex modulus of the elastic attenuation medium; obtaining an expression of the constitutive relation of the VTI attenuation medium in a frequency-wavenumber domain; expanding the complex modulus of the obtained elastic attenuation medium to VTI attenuation anisotropy, and converting the expression of the constitutive relation of the VTI attenuation medium in a frequency-wavenumber domain into an expression in a time-space domain; obtaining a DFL (design distribution model) viscoelasticity equation of the VTI medium based on an expression, a geometric equation and a motion balance equation of a constitutive relation of the VTI attenuation medium in a time-space domain; s150, determining underground medium parameters, and substituting the parameters into a VTI medium DFL viscoelasticity equation to calculate and obtain seismic wave field simulation values; the invention can make the earthquake simulation more accurate.

Description

Fractional order Laplace operator-based anisotropic attenuation medium simulation method
Technical Field
The application relates to the field of oil and gas detection, in particular to an anisotropic attenuation medium simulation method based on fractional order Laplace operators.
Background
In the development of oil and gas exploration, seismic exploration technology plays an extremely important role. At present, the development of seismic wave field numerical simulation aiming at complex media is a consensus in the industry at present, and the fractional order Decoupling (DFL) viscous wave equation developed in recent years has a plurality of unique advantages. However, numerical solution of the DFL viscous wave equation often faces a difficult problem. The Laplace operator, whose order varies with space, contained in this equation presents difficulties in the simulation of wavefields in inhomogeneous media. Currently, although anisotropic viscoelastic wave and acoustic wave equations in the form of DFLs already exist, the mixed-domain laplacian contained therein requires additional processing to be suitable for the simulation of inhomogeneous media.
Patent document (CN113341455A) discloses a method for numerically simulating viscous anisotropic seismic waves, which can accurately describe velocity and attenuation anisotropy at the same time, and facilitate numerical calculation equations, but only considers the attenuation and velocity anisotropy of longitudinal waves, and has a problem that the method is not suitable for actual underground conditions, so that the seismic simulation is inaccurate because seismic waves are transmitted in the underground as longitudinal waves and transverse waves at the same time.
Disclosure of Invention
In order to solve the technical problems, the anisotropic attenuation medium simulation method based on the fractional order Laplace operator is provided, anisotropy of longitudinal waves and anisotropy of transverse waves are considered, underground actual conditions can be better fitted, and seismic simulation is more accurate.
The embodiment of the invention adopts the technical scheme that:
the embodiment of the invention provides an anisotropic attenuation medium simulation method based on fractional order Laplace operators, which comprises the following steps:
s100, constructing a frequency dispersion relation of an elastic wave based on a preset frequency dispersion relation, wherein the elastic wave comprises a P wave and an S wave;
s110, obtaining the complex modulus of the elastic attenuation medium based on the constructed frequency dispersion relation of the elastic wave;
s120, obtaining an expression of the constitutive relation of the VTI attenuation medium in a frequency-wavenumber domain;
s130, expanding the complex modulus of the obtained elastic attenuation medium to VTI attenuation anisotropy, and converting the expression of the constitutive relation of the VTI attenuation medium in a frequency-wavenumber domain into an expression in a time-space domain;
s140, obtaining a DFL (doubly salient distortion) viscoelasticity equation of the VTI medium based on an expression, a geometric equation and a motion balance equation of a constitutive relation of the VTI attenuation medium in a time-space domain;
determining underground medium parameters, and substituting the underground medium parameters and the seismic wave parameters into the VTI medium DFL viscoelasticity equation to calculate and obtain seismic wave field simulation values;
the DFL viscoelasticity equation of the VTI medium is as follows:
Figure RE-GDA0003484890480000021
wherein t represents time, vxAnd vzThe vibration velocity components of the particles of the underground medium along the x direction and the z direction respectively, rho is the density of the underground medium, and sigma isxx、σzzAnd σxzThe horizontal, longitudinal and shear components of the stress tensor, respectively;
Figure RE-GDA0003484890480000031
wherein the content of the first and second substances,
Figure RE-GDA0003484890480000032
Q11and Q33Quality factors, Q, of the P-wave in the transverse and vertical directions, respectively55Is a quality factor of the S-wave,
Figure RE-GDA0003484890480000033
Figure RE-GDA0003484890480000034
vpis P wave edgeSpeed in the direction of the axis of symmetry; v. ofsIs the velocity of the S-wave, ε and δ being the velocity anisotropy parameters, δQIs a damping anisotropy parameter;
Figure RE-GDA0003484890480000035
velocity of the P-wave along the horizontal axis of the array, v2＝vpThe velocity of the P wave along the vertical array axis; omega0Is a reference angular frequency.
The invention has at least the following technical effects: the anisotropic attenuation medium simulation method based on the fractional order Laplace operator provided by the invention considers the anisotropies of longitudinal waves and transverse waves, can better fit the underground practical situation, and can accurately depict the attenuation mechanism of seismic waves in the underground propagation process.
Drawings
In order to more clearly illustrate the technical solutions in the embodiments of the present application, the drawings needed to be used in the description of the embodiments are briefly introduced below, and it is obvious that the drawings in the following description are only some embodiments of the present application, and it is obvious for those skilled in the art to obtain other drawings based on these drawings without creative efforts.
Fig. 1 is a flowchart of an anisotropic attenuation medium simulation method based on fractional order lagrange operators according to an embodiment of the present application;
fig. 2(a) and fig. 2(b) are schematic diagrams for verifying the simulation accuracy of the DFL viscoelasticity equation of the VTI medium provided by the embodiment of the present application by comparing with an analytic solution;
fig. 3 is a schematic diagram for verifying a decoupling characteristic of a DFL viscoelasticity equation of a VTI medium provided in an embodiment of the present application;
4 a-c are three layer models for testing the ability of the DFL viscoelastic equation of the VTI medium to process the medium with the Q value and the violent change of speed provided by the embodiment of the application;
FIG. 5 is a diagram of an observation system for estimating a Q value by 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 a spectral ratio method according to an embodiment of the present disclosure;
FIG. 7 is a graph of inverse comparison of quality factors provided by an embodiment of the present application;
fig. 8a to 8d are schematic diagrams of 2007BP TTI models provided in the embodiments of the present application;
FIG. 9a is a snapshot of a wave field in an anisotropic elastic medium according to an embodiment of the present application;
FIGS. 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 illustration of a seismic recording in an anisotropic elastic medium provided by an embodiment of the present application;
FIGS. 10b and 10c are schematic diagrams of 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 model of a Solton cell;
FIG. 11b is a schematic diagram of a resampled P-wave velocity model;
FIGS. 11c and 11d are schematic cross-sectional views of a different latitude and a different longitude of the Soluton's slot model, respectively;
12a to 12f are wavefield snapshots at different time points calculated by a DFL viscoelastic wave equation provided in an embodiment of the present application;
FIGS. 13a to c are schematic diagrams of three-dimensional shot-shared gather calculation using the viscoelastic wave equation provided in the embodiment of the present application and the conventional time-fractional order viscoelastic wave equation elastic wave equation;
FIG. 14 shows a viscoelastic wave equation and v provided in an example of the present applicationxComponent vyComponent sum vzAnd the schematic diagram is used for tracking and comparing a video envelope objective function, a time fractional order viscoelastic wave equation and a time-frequency phase error in the seismic record of the component.
Detailed Description
The technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the drawings in the embodiments of the present application, and it is obvious that the described embodiments are only a part of the embodiments of the present application, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present application.
As shown in fig. 1, the method for simulating an anisotropic attenuation medium based on a fractional order lawsonia operator provided by the embodiment of the present invention includes the following steps:
s100, constructing a frequency dispersion relation of the elastic wave based on a preset frequency dispersion relation, wherein the elastic wave comprises a P wave and an S wave.
In the embodiment of the present invention, the predetermined dispersion relation is an approximate dispersion relation of a constant Q model. The approximate dispersion relation of the constant Q model represents the variation relation of frequency and speed, and the formula is as follows:
ω2＝d1k+d2k2+d3k3+d4(iω)k+d5(iω)k2 (1)
(d1,d2,d3,d4,d5)＝(-γvω0,v2,γv3ω0 -1,πγv,πγ2v2ω0 -1) (2)
where k is the wave number, ω is the angular frequency, ω0For reference angular frequency, γ is space-variant fractional order, γ is arctanQ-1The Q is a quality factor of the elastic wave; v is a reference speed, v ═ v0cos(πγ/2)，v0At a reference frequency omega0The reference phase velocity of (c).
The elastic wave contains P wave and S wave components, and both of them satisfy the approximate dispersion relation, so the dispersion relation of the constructed constant Q model of the elastic wave is as follows:
Figure RE-GDA0003484890480000061
wherein v isαTo be shotVelocity of the wave, gammaα＝arctanQα -1/π，Qαα is a quality factor of the elastic wave, P or S.
And S110, obtaining the complex modulus of the elastic attenuation medium based on the constructed frequency dispersion relation of the elastic wave.
In an embodiment of the present invention, the complex modulus of the elastic damping medium may be:
Figure RE-GDA0003484890480000062
and S120, obtaining an expression of the constitutive relation of the VTI attenuation medium in a frequency-wavenumber domain.
In the embodiment of the present application, a simpler anisotropic medium is considered: a transversely isotropic media (VTI) with a perpendicular axis of symmetry. The stress-strain relationship is determined by an elastic stiffness matrix closely related to the elastic property of the medium, and for a VTI elastic medium, the expression of the elastic stiffness matrix in the Vogit form is
Figure RE-GDA0003484890480000063
The elastic stiffness coefficient is:
Figure RE-GDA0003484890480000064
where ρ represents the density of the subterranean medium; v. ofpRepresenting the speed of the P wave along the direction of the symmetry axis; v. ofsIs the velocity of the S-wave, and ε and δ are velocity anisotropy parameters.
For VTI attenuating media, its Q matrix can be expressed as
Figure RE-GDA0003484890480000071
Wherein Q11And Q33Corresponding to the quality factors of the P-wave in the transverse and vertical directions, respectivelyQ55Quality factor, Q, corresponding to S-wave12＝Q13,Q55＝Q66. Furthermore, a damping anisotropy parameter ε similar to the Thomsen coefficient is usedQAnd deltaQThus, Q is established11，Q33And Q13The relationship between them, specifically,
Figure RE-GDA0003484890480000072
Figure RE-GDA0003484890480000073
the anisotropic attenuation medium corresponds to a complex stiffness coefficient (or complex modulus, recorded as
Figure RE-GDA0003484890480000074
) It can be converted from the following relation through a rigidity matrix
Figure RE-GDA0003484890480000075
The constitutive relation of the VTI attenuation medium is expressed in a frequency-wavenumber domain as
Figure RE-GDA0003484890480000076
Figure RE-GDA0003484890480000077
Figure RE-GDA0003484890480000078
Wherein the content of the first and second substances,
Figure RE-GDA0003484890480000079
and
Figure RE-GDA00034848904800000710
the expressions in frequency for the horizontal, longitudinal and shear components of the stress tensor respectively,
Figure RE-GDA00034848904800000711
and
Figure RE-GDA00034848904800000712
the horizontal, longitudinal and shear components of the strain tensor are respectively expressions in frequency,
Figure RE-GDA00034848904800000713
and
Figure RE-GDA00034848904800000714
respectively, the complex modulus of the anisotropic attenuation medium.
S130, expanding the complex modulus of the obtained elastic attenuation medium to VTI attenuation anisotropy, and converting the expression of the constitutive relation of the VTI attenuation medium in a frequency-wavenumber 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
Figure RE-GDA0003484890480000081
And
Figure RE-GDA0003484890480000082
wherein the content of the first and second substances,
Figure RE-GDA0003484890480000083
Figure RE-GDA0003484890480000084
s132, since formula (14) is only related to v1Accordingly, the formula (4) can be used to derive the formula (14)
Figure RE-GDA0003484890480000085
In particular, v is1V in alternative (4)αTo obtain a compound of formula (14)
Figure RE-GDA0003484890480000086
The approximate expression of (c) is:
Figure RE-GDA0003484890480000087
s133, since formula (15) is only related to vsAccordingly, the formula (4) can be used to derive the formula (15)
Figure RE-GDA0003484890480000088
And
Figure RE-GDA0003484890480000089
in particular, v issV in alternative formula (1)αTo obtain a compound of formula (15)
Figure RE-GDA00034848904800000810
And
Figure RE-GDA00034848904800000811
the approximate expression of (c) is:
Figure RE-GDA00034848904800000812
Figure RE-GDA00034848904800000813
s134, substituting formula (16) for formula (14) and substituting formula (17) and formula (18) for formula (15), obtains the expression of formula (11) in the time-space domain:
Figure RE-GDA00034848904800000814
s135, similar operations to those of S131 to S134 are performed on expressions (12) and (13), respectively, that is, decomposed into P-wave and S-wave components, respectively, and then expressions of expressions (12) and (13) in the time-space domain are obtained based on expression (4):
Figure RE-GDA00034848904800000815
σxz＝L55exz (21)
wherein e isxx、ezzAnd exzThe horizontal, longitudinal and shear components of the strain tensor, respectively.
And S140, obtaining a DFL (doubly salient distortion) viscoelasticity equation of the VTI medium based on the expression, the geometric equation and the motion balance equation of the constitutive relation of the VTI attenuation medium in the time-space domain.
In the embodiment of the present invention, the geometric equation and the motion balance equation are:
Figure RE-GDA0003484890480000091
wherein i and j are x or z.
In the embodiment of the present invention, the DFL viscoelasticity equation of the VTI medium based on the new dispersion relation obtained by combining equations (19) to (22) is:
Figure RE-GDA0003484890480000092
wherein t represents time, vxAnd vzThe vibration velocity components, sigma, of the particles of the underground medium in the x-direction and the z-direction, respectivelyxx、σzzAnd σxzHorizontal, longitudinal and shear components of the stress tensor, respectively；
Figure RE-GDA0003484890480000093
Wherein the content of the first and second substances,
Figure RE-GDA0003484890480000101
vpthe speed of the P wave along the direction of the symmetry axis; v. ofsIs the velocity of the S-wave, epsilon and delta are velocity anisotropy parameters,
Figure RE-GDA0003484890480000102
velocity of the P-wave along the horizontal axis of the array, v2＝vpAnd is the velocity of the P wave along the vertical alignment axis.
And S150, determining underground medium parameters, and substituting the underground medium parameters into the VTI medium DFL viscoelasticity equation to calculate and obtain seismic wave field simulation values.
On the basis of obtaining the VTI medium DFL viscoelasticity equation, the seismic wave field simulation is realized by determining the parameters of the underground medium, including the underground medium density, P wave and S wave quality factors, P wave velocity and S wave velocity, attenuation anisotropy parameters, velocity anisotropy parameters and the like, and substituting the parameters into the VTI medium DFL viscoelasticity equation, so that the seismic wave field simulation is obtained.
(test example 1)
In this test example, the accuracy of the DFL viscoelastic wave equation shown in the above equation (23) is first verified by comparison with an analytical solution, which can be derived from the green's function in a homogeneous medium. The homogeneous medium model used has a size of 512 × 512 grid points and a spatial sampling interval of 10 m. The dominant frequency of a Rake wavelet seismic source arranged in the center of the uniform medium is 25Hz, and the reference speeds of the reference frequency and the seismic source which are consistent are c0p3000m/s and c0s2000m/s, a medium density of 2200 kg/m3. Record v at (3160m )xThe components, fig. 2a and 2b, correspond to the fully elastic medium and the damping medium (Q), respectivelyp＝32，Qs20). In the figure, the thick line, the dotted line and the thin black line correspond to the analytic solution, the numerical solution and the residual error. It can be found that the proposed new equation can be matched with the analytic solution well no matter in the attenuation medium or the complete elastic medium, and the accuracy of the DFL viscous elastic wave equation is proved.
(test example 2)
In the test example, the decoupling characteristic of the DFL viscoelastic wave equation is further verified. In fig. 3, the four wavefield snapshots correspond to the full elastic wave equation, the damped dominant equation, the dispersive dominant equation, and the proposed DFL visco-elastic wave equation, respectively. It can be observed that the attenuation principal equations describe well the attenuation of amplitude energy compared to elastic waves, the dispersion principal equations characterize the velocity dispersion caused by formation attenuation (phase lag can be observed), and the wavefield modeled by the DFL viscoelastic wave equation has both effects.
(test example 3)
In the test example, the ability of the DFL viscoelastic wave equation to handle non-uniform attenuating media was further verified. In this test example, three layer models were designed to examine their ability to handle strongly varying media in terms of Q and speed. As shown in fig. 4(a), the velocity in model 1 is a laminar medium and the Q value is uniform, the velocity in model 2 is uniform and the Q value is laminar, and both the velocity and Q in model 3 change in a laminar manner. The size of the group of models is 800 multiplied by 800 grid points, the space sampling interval is 10m, the horizontal interface is located at the depth 4400m, the adopted time step is 1ms, and the dominant frequency of a Rake wavelet source arranged in the center of a medium is 25 Hz. However, in the case of sharp Q-contrast (model 2 and model 3), the snapshot residual between the schemes is more significant than the reference snapshot, which is the result of the phase distortion additive effect (see fig. 4 c). The third column in fig. 4b has only slight residual errors and the two zone 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 attenuating media.
(test example 4)
In this test example, the VTI-DFL viscoelasticity equation presented for verification describes the quasi-variation of attenuation with directionIndeed, a set of anisotropic homogeneous model experiments was performed. Firstly, solving the proposed VTI-DFL viscoelasticity equation to obtain a corresponding snapshot. On one hand, the Q value of each direction can be estimated by using a spectral ratio method; on the other hand, by solving the christofel 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 and the estimated Q value. Fig. 5 shows a simulated geometry in which the 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 homogeneous anisotropic media models whose attenuation anisotropy parameters are shown in table 1. The model size is 400 × 400 grid points, with a spatial sampling interval of 10m and a time step of 1 ms. The remaining parameters are vp＝3.0km/s，vs1.5km/s, Thomsen velocity anisotropy parameters 0.16, 0.08, medium density 2000kg/m 3. To quantitatively express the intensity of the medium attenuation with direction change, an intensity coefficient Strength (Q) is definedp,s)＝1-min(Qp,s)/max(Qp,s)。
TABLE 1 set of elastic damping medium models for damping the variation of anisotropic parameters
Figure RE-GDA0003484890480000111
Figure RE-GDA0003484890480000121
The common shot gather recorded by the observation system used in the spectral ratio method can be shown in fig. 6, in which the left side is the inner concentric circle demodulator probe and the right side is the outer concentric circle demodulator probe in fig. 6. The inverse comparison of the obtained quality factors can be shown in fig. 7, wherein the circles represent the Q values estimated by the spectral ratio method, and the black solid line represents the theoretical Q values obtained by the christofel equation.
(test example 5)
In this test example, the proposed VTI-DFL stiction was verified with a complex modelThe simulation accuracy of the elastic equation was modeled as the 2007BP TTI model (without considering the tilt angle). FIGS. 8a to 8d show P-wave velocity and P-wave quality factor Q, respectivelypThomsen velocity anisotropy parameters δ and ε. S-wave velocity and quality factor QsRespectively through vs＝vp/1.73，Qs＝QpA/1.2 transformation to a damping anisotropy parameter of δQ＝2δ，ε Q3 epsilon. The size of the model is 500 multiplied by 270 grid points, the space sampling interval is 10m, the time sampling interval is 1ms, the simulation propagation time is 2.5s, the adopted excitation seismic source is a Rake wavelet with the main frequency of 25Hz, the position of the seismic source is (2500m and 30m), and the position of the wave detection point is positioned at each grid position with the depth of 50m along the horizontal direction. To evaluate the simulation accuracy, we calculated a reference solution from the existing VTI-TF viscoelastic wave equation (Zhu,2017) that can describe the direction-dependent attenuation, but requires significant computational and memory overhead. FIG. 9a shows a wave field snapshot in an anisotropic elastic medium, FIG. 9b shows a wave field snapshot obtained by simulation of a VTI-DFL viscoelastic equation, and FIG. 9c shows a wave field snapshot obtained by simulation of a VTI medium time fractional order viscoelastic wave equation, which does not need approximation when expanding to anisotropy, so that the change of formation attenuation along with direction can be well described, but because a time fractional order operator contained in the equation has the problems of calculation and storage in a numerical solving process, the equation is not suitable for large-scale numerical simulation; and comparing with the simulation result of the equation, the simulation precision of the DFL viscosity equation of the VTI medium 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. After analysis and comparison, the stratum attenuation effect obviously weakens the amplitude energy compared with a non-attenuation medium; the weak difference in fig. 9d and 10d furthermore indicates that the proposed VTI-DFL viscoelastic equation still has a 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 three-dimensional true model of Solton Hatch, Calif., southThe model simulates the propagation of seismic waves. Location map of solton sea as shown in fig. 11a, the initial P-wave velocity model was provided by (Ajala, 2019) and included 19 grids in depth from 0 to 9 km, 101 grids in longitude 115.7 ° W to 116.7 ° W, and 91 grids in latitude from 33.3 ° N34.2 ° N. Since the original data is not uniform, linear interpolation and resampling can obtain a new p-wave velocity model, which contains 505 × 455 × 48 grid points with a uniform grid spacing of 200m (see fig. 11 b). Fig. 11c and 11d show some sections at different latitudes and longitudes, respectively. Velocity of S wave is represented by vs＝vpThe velocity model of/1.73 (Li,1993), the Q model being derived from QP＝QS＝3.516(vp/1000)22(Li, 1993). A source of rake wavelets with a dominant frequency of 5hz is placed at (50000m,46000m, 1600m), with receivers at each grid point in the x-y plane of the earth's surface. The time is 15s and the time step is 5.0 ms. Fig. 12a to 12f show snapshots of components simulated with the proposed DFL viscoelastic wave equation at different times. It is important to note that some energy is observed to be trapped in the elliptically marked region, where seismic waves also travel at a slower rate than in the surrounding area. The reason for this phenomenon can be found in the model of the solton sea chest, and it can be seen from fig. 11b that there are some low-speed, highly attenuating deposits in the shallow part of this region. In order to evaluate the accuracy of the proposed DFL viscoelasticity equation of the solton heyday model, a time-fractional order viscoelasticity wave equation is used as a reference solution. Since computing the time fractional order visco-elastic wave equation requires a large amount of computer memory, in order for the simulation to be feasible, the model is resampled to 252 × 227 × 24 grid points, with grid spacing of 20m each. The seismic source dominant frequency is increased to 25Hz and the time step is reduced to 1.0 ms. FIG. 13 shows a 3D common shot gather, where the horizontal slice shows a snapshot at 1.1s and the two vertical sections show the two common shot gathers in the x and y directions. FIGS. 13a-c correspond to the DFL viscoelasticity wave equation, time fractional order viscoelasticity wave equation, and elasticity wave equation, respectively, as set forth in the present invention. The common shot gathers in FIGS. 13a and 13b appear very similar, their amplitude is significantly weaker than the elastic case due to energy attenuation(FIG. 13 c). For better comparison, the seismic recordings of (1500m,1500m,20m) are further shown in FIG. 14. The columns from left to right correspond to vx，vyAnd vzThe component (c). 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 for the three components (middle row in fig. 14) match well, with TFEM and TFPM within 2%, indicating that the proposed DFL viscoelastic wave equation generates seismic records that are nearly identical to the reference solution.
In summary, the method provided by the embodiment of the present application provides an acoustic wave equation based on a fractional laplacian, can accurately depict the characteristics of seismic wave velocity and attenuation changing with the azimuth at the same time, and does not include the fractional laplacian changing with the space, so that the propagation of seismic waves in a non-uniform medium can be naturally and accurately calculated, and the dimension reduction processing (from three-dimensional to two-dimensional, for example) of an anisotropic medium acoustic wave equation can be realized by simplifying the acoustic wave equation, so that the method can be widely applied to acoustic wave field simulation, absorption attenuation compensation reverse time migration imaging and full waveform inversion in two-dimensional (2D) and three-dimensional (3D) complex attenuation anisotropic media. And moreover, due to the fact that anisotropy of longitudinal waves and transverse waves is considered at the same time, the underground practical situation can be better attached, and the attenuation mechanism of 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 can be disposed in an electronic device to store at least one instruction or at least one program for implementing a method of the method embodiments, where the at least one instruction or the at least one program is loaded into and executed by a processor to implement the method provided by the above embodiments.
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 means for causing an electronic device to carry out the steps of the method according to various exemplary embodiments of the present application described above in the present description, 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 illustration, it should be understood by those skilled in the art that the above illustration is only for purposes of illustration and is not intended to limit the scope of the present application. It will also be appreciated by those skilled in the art that various modifications may be made to the embodiments without departing from the scope and spirit of the present application. The scope of the disclosure of the present application is defined by the appended claims.

Claims (9)

1. An anisotropic attenuation medium simulation method based on fractional Laplace operator is characterized by comprising the following steps:
s100, constructing a frequency dispersion relation of an elastic wave based on a preset frequency dispersion relation, wherein the elastic wave comprises a P wave and an S wave;
s110, obtaining the complex modulus of the elastic attenuation medium based on the constructed frequency dispersion relation of the elastic wave;
s120, obtaining an expression of the constitutive relation of the VTI attenuation medium in a frequency-wavenumber domain;
s130, expanding the complex modulus of the obtained elastic attenuation medium to VTI attenuation anisotropy, and converting the expression of the constitutive relation of the VTI attenuation medium in a frequency-wavenumber domain into an expression in a time-space domain;
s140, obtaining a DFL (doubly salient distortion) viscoelasticity equation of the VTI medium based on an expression, a geometric equation and a motion balance equation of a 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 viscoelasticity equation to calculate and obtain seismic wave field simulation values;
the DFL viscoelasticity equation of the VTI medium is as follows:
Figure FDA0003426036550000011
wherein t representsTime, vxAnd vzThe vibration velocity components of the particles of the underground medium along the x direction and the z direction respectively, rho is the density of the underground medium, and sigma isxx、σzzAnd σxzThe horizontal, longitudinal and shear components of the stress tensor, respectively;
Figure FDA0003426036550000021
wherein the content of the first and second substances,
Figure FDA0003426036550000022
Q11and Q33Quality factors, Q, of the P-wave in the transverse and vertical directions, respectively55Is a quality factor of the S-wave,
Figure FDA0003426036550000023
Figure FDA0003426036550000024
vpthe speed of the P wave along the direction of the symmetry axis; v. ofsIs the velocity of the S-wave, ε and δ being the velocity anisotropy parameters, δQIs a damping anisotropy parameter;
Figure FDA0003426036550000025
velocity of the P-wave along the horizontal axis of the array, v2＝vpThe velocity of the P wave along the vertical array axis; omega0Is a reference angular frequency.
2. The method of claim 1, wherein the predetermined dispersion relation is an approximate dispersion relation of a constant Q model;
the frequency dispersion relation of the constructed elastic wave is as follows:
Figure FDA0003426036550000031
where k is the wave number, ω is the angular frequency, vαIs the velocity of the elastic wave, gammaα＝arctanQα -1/π，Qαα is a quality factor of the elastic wave, P or S.
3. The method of claim 2, wherein the elastic damping medium has a complex modulus of
Figure FDA0003426036550000032
4. The method of claim 3, wherein the expression of the constitutive relation of the VTI attenuating medium in the frequency-wavenumber domain satisfies the following condition:
Figure FDA0003426036550000033
Figure FDA0003426036550000034
Figure FDA0003426036550000035
wherein the content of the first and second substances,
Figure FDA0003426036550000036
and
Figure FDA0003426036550000037
the expressions in frequency for the horizontal, longitudinal and shear components of the stress tensor respectively,
Figure FDA0003426036550000038
and
Figure FDA0003426036550000039
are respectively strain tensionThe expression in frequency of the horizontal component, the longitudinal component and the shear component of the quantity,
Figure FDA00034260365500000310
and
Figure FDA00034260365500000311
respectively, the complex modulus of the anisotropic attenuation medium.
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
Figure FDA00034260365500000312
And
Figure FDA00034260365500000313
wherein the content of the first and second substances,
Figure FDA00034260365500000314
Figure FDA00034260365500000315
s132, mixing v1V in alternative formula (1)αTo obtain a compound of formula (6)
Figure FDA00034260365500000316
The approximate expression of (c) is:
Figure FDA00034260365500000317
s133, mixing vsV in alternative formula (1)αTo give compounds of formula (7)
Figure FDA00034260365500000318
And
Figure FDA00034260365500000319
the approximate expression of (c) is:
Figure FDA00034260365500000320
Figure FDA0003426036550000041
s134, substituting formula (8) for formula (6) and substituting formula (9) and formula (10) for formula (7), obtains the expression of formula (3) in the time-space domain:
Figure FDA0003426036550000042
s135, similar operations to S131 to S134 are performed on equations (4) and (5), respectively, resulting in expressions of equations (4) and (5) in the time-space domain, respectively:
Figure FDA0003426036550000043
σxz＝L55exz (13)
wherein e isxx、ezzAnd exzThe horizontal, longitudinal and shear components of the strain tensor, respectively.
6. The method of claim 1, wherein the geometric equation and the motion balance equation are:
Figure FDA0003426036550000044
wherein 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, the at least one instruction or the at least one program being loaded and executed by a processor to implement the method of any of claims 1-6.
8. An electronic device comprising a processor and the non-transitory computer readable storage medium of claim 7.
9. A computer program product comprising a computer program, characterized in that the computer program is executed by a processor for implementing the method according to any of claims 1-6.
CN202111578115.5A 2021-12-22 2021-12-22 Anisotropic attenuation medium simulation method based on fractional order Laplace operator Active CN114114403B (en)

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 true CN114114403A (en) 2022-03-01
CN114114403B 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)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115392090A (en) * 2022-09-13 2022-11-25 中国矿业大学 Method for predicting seismic wave frequency dispersion and attenuation characteristics based on three-dimensional digital core
CN117471531A (en) * 2023-11-01 2024-01-30 中国矿业大学 VTI viscoelastic wave equation numerical simulation method based on power law frequency-varying Q effect

Citations (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20080298174A1 (en) * 2007-06-01 2008-12-04 Baker Hughes Incorporated Method for determining seismic anisotropy
US20090133476A1 (en) * 2007-05-01 2009-05-28 Paul Michaels Determination of permeability from damping
US20110222370A1 (en) * 2010-03-12 2011-09-15 Jon Downton Methods and systems for performing azimuthal simultaneous elastic inversion
CN102508296A (en) * 2011-11-14 2012-06-20 中国石油天然气股份有限公司 Method and device for analyzing dispersion and attenuation of unsaturated double-porosity medium earthquake waves
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
CN104360383A (en) * 2014-11-12 2015-02-18 中国石油大学(华东) Method and system for predicting seismic wave attenuation
CN104570084A (en) * 2013-10-29 2015-04-29 中国石油化工股份有限公司 Cross-scale seismic rock physical attenuation model and method for predicating attenuation and dispersion
CN105044771A (en) * 2015-08-05 2015-11-11 北京多分量地震技术研究院 3D TTI double-phase medium seismic wave field value simulation method based on finite difference method
CN105137486A (en) * 2015-09-01 2015-12-09 中国科学院地质与地球物理研究所 Elastic wave reverse-time migration imaging method and apparatus in anisotropic media
US20150362622A1 (en) * 2014-06-17 2015-12-17 Huseyin Denli Fast Viscoacoustic and Viscoelastic Full Wavefield Inversion
US20170371050A1 (en) * 2016-06-28 2017-12-28 Junzhe Sun Reverse time migration in anisotropic media with stable attenuation compensation
CN110058310A (en) * 2019-04-16 2019-07-26 中国石油大学(华东) The seimic wave velocity frequency dispersion and decaying prediction technique and device of earthquake rock
CN110542928A (en) * 2018-05-28 2019-12-06 中国石油化工股份有限公司 Seismic response simulation method based on VTI anisotropic propagation matrix
US20200348430A1 (en) * 2019-05-03 2020-11-05 Exxonmobil Upstream Research Company Inversion, Migration, and Imaging Related to Isotropic Wave-Mode-Independent Attenuation
CN113189642A (en) * 2021-04-28 2021-07-30 中国石油化工集团有限公司 Seismic source linear scanning signal design method based on forced vibration
CN113341455A (en) * 2021-06-24 2021-09-03 中国石油大学(北京) Viscous anisotropic medium seismic wave numerical simulation method, device and equipment

Patent Citations (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090133476A1 (en) * 2007-05-01 2009-05-28 Paul Michaels Determination of permeability from damping
US20080298174A1 (en) * 2007-06-01 2008-12-04 Baker Hughes Incorporated Method for determining seismic anisotropy
US20110222370A1 (en) * 2010-03-12 2011-09-15 Jon Downton Methods and systems for performing azimuthal simultaneous elastic inversion
CN102508296A (en) * 2011-11-14 2012-06-20 中国石油天然气股份有限公司 Method and device for analyzing dispersion and attenuation of unsaturated double-porosity medium earthquake waves
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
CN104570084A (en) * 2013-10-29 2015-04-29 中国石油化工股份有限公司 Cross-scale seismic rock physical attenuation model and method for predicating attenuation and dispersion
US20150362622A1 (en) * 2014-06-17 2015-12-17 Huseyin Denli Fast Viscoacoustic and Viscoelastic Full Wavefield Inversion
CN104360383A (en) * 2014-11-12 2015-02-18 中国石油大学(华东) Method and system for predicting seismic wave attenuation
CN105044771A (en) * 2015-08-05 2015-11-11 北京多分量地震技术研究院 3D TTI double-phase medium seismic wave field value simulation method based on finite difference method
CN105137486A (en) * 2015-09-01 2015-12-09 中国科学院地质与地球物理研究所 Elastic wave reverse-time migration imaging method and apparatus in anisotropic media
US20170371050A1 (en) * 2016-06-28 2017-12-28 Junzhe Sun 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
CN110058310A (en) * 2019-04-16 2019-07-26 中国石油大学(华东) The seimic wave velocity frequency dispersion and decaying prediction technique and device of earthquake rock
US20200348430A1 (en) * 2019-05-03 2020-11-05 Exxonmobil Upstream Research Company Inversion, Migration, and Imaging Related to Isotropic Wave-Mode-Independent Attenuation
CN113189642A (en) * 2021-04-28 2021-07-30 中国石油化工集团有限公司 Seismic source linear scanning signal design method based on forced vibration
CN113341455A (en) * 2021-06-24 2021-09-03 中国石油大学(北京) Viscous anisotropic medium seismic wave numerical simulation method, device and equipment

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
J. J. S. DE FIGUEIREDO: ""Shear wave anisotropy from aligned inclusions: ultrasonic frequency dependence of velocity and attenuation"", 《GEOPHYSICAL JOURNAL INTERNATIONAL》 *
赵海波: ""松辽盆地古龙页岩跨频带岩石物理测量"", 《大庆石油地质与开发》 *
高宇航: ""粘弹VTI介质数值模拟"", 《中国优秀硕士学位论文全文数据库》 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115392090A (en) * 2022-09-13 2022-11-25 中国矿业大学 Method for predicting seismic wave frequency dispersion and attenuation characteristics based on three-dimensional digital core
CN117471531A (en) * 2023-11-01 2024-01-30 中国矿业大学 VTI viscoelastic wave equation numerical simulation method based on power law frequency-varying Q effect
CN117471531B (en) * 2023-11-01 2024-04-19 中国矿业大学 VTI viscoelastic wave equation numerical simulation method based on power law frequency-varying Q effect

Also Published As

Publication number Publication date
CN114114403B (en) 2023-06-27

Similar Documents

Publication Publication Date Title
US8560242B2 (en) Pseudo-analytical method for the solution of wave equations
CN114114403B (en) Anisotropic attenuation medium simulation method based on fractional order Laplace operator
Withers et al. Memory-efficient simulation of frequency-dependent Q
Johnson Source mechanisms of induced earthquakes at The Geysers geothermal reservoir
MX2011003850A (en) Image domain signal to noise estimate.
CN109143351B (en) Pre-stack anisotropy characteristic parameter inversion method and computer readable storage medium
US9753166B2 (en) P-wave and S-wave separation of seismic data in the presence of statics and irregular geometry
CN111123354A (en) Method and equipment for predicting dense gas layer based on frequency-dependent reflection amplitude attenuation
CN110286410B (en) Fracture inversion method and device based on diffracted wave energy
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
NO20150821A1 (en) Efficient Wavefield Extrapolation in Anisotropic Media
Golubev et al. Continuum model of fractured media in direct and inverse seismic problems
US11199641B2 (en) Seismic modeling
Tang et al. Single-point dispersion measurement of surface waves combining translation, rotation and strain in weakly anisotropic media: theory
CN113552624B (en) Porosity prediction method and device
Chang et al. Elastic full waveform inversion for salt model building in Gulf of Mexico
AU2015290248A1 (en) System and method for rock property estimation of subsurface geologic volumes
CN109521470B (en) Method for analyzing influence of geological structure on seismic inversion crack density
De Matteis et al. BISTROP: Bayesian inversion of spectral‐level ratios and P‐wave polarities for focal mechanism determination
Sidler et al. Wave reflection at an anelastic transversely isotropic ocean bottom
CN114139335A (en) Viscous sound wave simulation method based on single relaxation time lattice Boltzmann model
MX2011003852A (en) Time reverse imaging attributes.
CN102096100A (en) Logging curve and seismographic record full-wave matching method and device
Zhu Seismic modeling, inversion, and imaging in attenuating media

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