CN111913217A - Fracture type stratum seismic scattering wave field characteristic simulation method based on imaginary number domain Rytov approximation - Google Patents

Fracture type stratum seismic scattering wave field characteristic simulation method based on imaginary number domain Rytov approximation Download PDF

Info

Publication number
CN111913217A
CN111913217A CN202010794592.4A CN202010794592A CN111913217A CN 111913217 A CN111913217 A CN 111913217A CN 202010794592 A CN202010794592 A CN 202010794592A CN 111913217 A CN111913217 A CN 111913217A
Authority
CN
China
Prior art keywords
wave field
field
underground
grid
medium
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.)
Pending
Application number
CN202010794592.4A
Other languages
Chinese (zh)
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.)
Beijing Research Center of CNOOC China Ltd
CNOOC China Ltd
Original Assignee
Beijing Research Center of CNOOC China Ltd
CNOOC China Ltd
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 Beijing Research Center of CNOOC China Ltd, CNOOC China Ltd filed Critical Beijing Research Center of CNOOC China Ltd
Priority to CN202010794592.4A priority Critical patent/CN111913217A/en
Publication of CN111913217A publication Critical patent/CN111913217A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/282Application of seismic models, synthetic seismograms
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30Analysis
    • G01V1/301Analysis for determining seismic cross-sections or geostructures
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/60Analysis
    • G01V2210/64Geostructures, e.g. in 3D data cubes
    • G01V2210/642Faults
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/60Analysis
    • G01V2210/66Subsurface modeling

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)

Abstract

The invention relates to a fracture type stratum scattered wave field characteristic simulation method with an imaginary number domain Rytov approximation, which obtains an analytic expression of a background wave field and a scattered wave field through an imaginary number domain Rytov approximation scattered wave field equation and constructs a Green function integral operator on the basis; then, deducing a total wave field and a backscattering field phase screen propagation operator of forward propagation of the three-dimensional seismic waves by using the sheet approximation and small angle approximation conditions; and finally, dividing the underground crack medium by adopting a variable grid, defining the space sampling interval of the underground crack medium, determining an encryption boundary, performing forward simulation on the region containing the crack in the underground crack medium by adopting an encryption grid Green function integral, and performing forward simulation on the region without the crack by adopting a coarse grid Green function integral. The invention can simulate the seismic scattering wave field characteristics caused by the underground fracture, guide the earthquake to identify the fracture type reservoir, better solve the problems of simulation efficiency and stability and has higher advantages of improving the signal-to-noise ratio and precision of the simulation result.

Description

Fracture type stratum seismic scattering wave field characteristic simulation method based on imaginary number domain Rytov approximation
Technical Field
The invention relates to a forward modeling method for underground fracture seismic scattering wave field characteristics, in particular to a fracture type stratum seismic scattering wave field characteristic modeling method based on imaginary number domain Rytov approximation, and belongs to the field of exploration geophysics.
Background
With the continuous development of geophysical exploration technology, exploration targets are changed from conventional oil and gas reservoirs to unconventional oil and gas reservoirs, and the requirements on seismic wave forward modeling are higher and higher. The underground homogeneous stratum is influenced by the movement of structures and the like, high-angle seams in directional arrangement are easy to develop to form a fracture type stratum, and for the unconventional exploration field, the fractures in the stratum have the functions of oil-gas migration and storage, so that the fracture type stratum becomes one of important exploration targets. Therefore, the forward simulation of the high-angle fracture medium with the oriented arrangement characteristics has important significance in the new-period oil and gas reservoir exploitation. The integral method is a main method for seismic wave equation forward simulation, and the differential method is to use a differential method to obtain a differential form solution of the wave equation by derivation, and high-order approximation is adopted in the derivation process. Although all information of wave field propagation can be guaranteed, for a small-scale heterogeneous body such as a crack, the problems of stability, dispersion and the like exist, and the stability and the accuracy of calculation are influenced.
The Green function integral forward modeling method can effectively model and depict the wavelength characteristics of seismic waves in the small-scale heterogeneous medium, and is widely concerned. Aki, the wake wave is the backscattering wave caused by the heterogeneity of the medium, Martin and Flatte establish a phase screen method suitable for the propagation of seismic waves in any medium, Wu and R.S utilize different methods to obtain the scattering wave field of the Green function and the background scattering field of elastic waves, and the forward modeling theory of the integral equation is gradually mature. In Qinxiefei et al and the bud, deep research is carried out on forward modeling research of a scattering wave field for seismic exploration application. However, on the premise of ensuring the accuracy, the improvement of the efficiency of the Green function integration method needs to be deeply researched.
Moczo firstly proposes a variable grid idea, and changes corresponding space and time step length for a region with large local change so as to improve efficiency on the premise of meeting precision. The Tessmer introduces the idea of changing the grid step length in the finite difference forward modeling of the wave equation, so that the seismic record of the depicted underground medium is more accurate. In order to meet the precision requirement of simulating hole seam media, Sunlienjie proposes a multi-stage grid-changing idea, so that the step length of the grid-changing is increased from ten times to more than one hundred times. However, currently variable meshes are mainly applied to differential simulation methods.
Disclosure of Invention
Aiming at the problems, the invention aims to provide a fracture type stratum scattering wave field characteristic simulation method based on imaginary number domain Rytov approximation, which can simulate the seismic scattering wave field characteristics caused by underground fractures and guide the earthquake to identify a fracture type reservoir, can better solve the problems of simulation efficiency and stability, and has higher advantages for improving the signal-to-noise ratio and the precision of a simulation result.
In order to achieve the purpose, the invention adopts the following technical scheme: a fracture type stratum scattered wave field characteristic simulation method of imaginary number domain Rytov approximation comprises the following steps:
step 1: establishing a ground observation system to determine the position of a ground seismic source excitation point and a receiving position, and designing a proper seismic source signal to be used as an initial wave field of a forward simulation system;
step 2: performing two-dimensional Fourier transform on the seismic source signal in the real number domain, converting the seismic source signal in the real number domain into the seismic source signal in the imaginary number domain, and taking the seismic source signal in the imaginary number domain as an incident wave field u of the underground fracture mediumi=u0(x, y, z; ω) in which u0Is a background wave field, x, y and z are coordinate axes of a Cartesian coordinate system, and omega is an angular frequency;
and step 3: and establishing a scattered field extrapolation equation based on the Rytov approximation:
obtaining the following result according to a Rytov transformation formula:
Figure BDA0002625065830000021
secondly, obtaining an imaginary number domain Rytov approximate scattered field equation U based on the Rytov approximationsComprises the following steps:
Figure BDA0002625065830000022
in the formula, i is an imaginary number unit; s0Uniform background slowness; k is a radical ofx、kyAnd kzWave numbers in x, y and z directions, respectively, and
Figure BDA0002625065830000023
wherein k is0Is the background wave number; delta z is the thickness of a thin plate which divides the underground fracture medium into equal thickness along the vertical direction; delta h is the horizontal slowness disturbance quantity between the thin plates;
and thirdly, based on an imaginary number domain Rytov approximate scattered field equation, the forward scattered field and the back scattered field are respectively characterized as follows:
Figure BDA0002625065830000024
Figure BDA0002625065830000025
in the formula, FT-1Is a two-dimensional inverse fourier transform operator.
And 4, step 4: firstly, respectively calculating the backscattering field of the upper-layer underground fracture medium by using virtual number domain Rytov approximate scattering wave propagation equations (3) and (4)
Figure BDA0002625065830000026
And forward scattered field
Figure BDA0002625065830000027
And back scattering field to upper underground crack medium
Figure BDA0002625065830000028
Storing; then based on Rytov approximation, the background wave field u of the upper-layer underground fracture medium is determined0(z1) And forward scattered field
Figure BDA0002625065830000029
Substitution into
Figure BDA00026250658300000210
Taking the calculation result as the incident wave field of the next layer of underground fracture medium, substituting the new incident wave field into the virtual number domain Rytov approximate scattering wave propagation equations (3) and (4), and respectively calculating the forward scattering field of the next layer of underground fracture medium
Figure BDA00026250658300000211
And a back scattered field
Figure BDA00026250658300000212
And back scattering field to the next layer of underground crack medium
Figure BDA00026250658300000213
Continuing to store;
and 5: calculating a virtual number domain coarse grid Rytov approximate scattering wave field according to the step 3 and the step 4, when the simulated wave field is transmitted to an encryption boundary of the underground crack medium, carrying out interpolation encryption processing on the coarse grid wave field in the horizontal direction by using an interpolation principle until the encryption multiple of the grid is met, segmenting the coarse grid in the vertical direction, and carrying out extrapolation calculation on the wave field of the underground crack medium in the encryption grid layer by layer based on a wave field extrapolation theory;
step 6: when the wave field of the underground crack medium in the encryption grid is extrapolated to the encryption boundary of the underground crack medium, transforming the wave field of the encryption grid in the horizontal direction into the coarse grid based on a conversion algorithm, and continuously extrapolating the wave field of the underground crack medium in the coarse grid under the encryption boundary layer by layer based on a wave field extrapolation theory in the vertical direction;
and 7: when the wave field is transmitted to the bottom of the underground fracture medium, the wave field extrapolation of the underground fracture medium is calculated; at this time, from the bottom of the underground fracture medium, the backward scattering field calculated by using the bottom of the underground fracture medium is used as the incident wave field of the underground fracture medium of the upper layer, and the backward propagation backward scattering field is calculated by using the formula (4)
Figure BDA0002625065830000031
And 8: backward propagating back scattered field
Figure BDA0002625065830000032
And (4) carrying out superposition processing on the backward scattering field of the n-1 layer underground fracture medium stored in the step (4), taking the superposition result as an incident wave field of the previous layer underground fracture medium, and calculating a backward propagation forward scattering field by using the formula (3)
Figure BDA0002625065830000033
And step 9: repeating the step 4-8, backward propagating the wave field, and finishing the backward extrapolation calculation of the wave field of the underground fracture medium when the wave field reaches the receiving position of the underground fracture medium; and finally, converting the imaginary number domain seismic scattering record into a real number domain seismic scattering record by utilizing two-dimensional Fourier inverse transformation, characterizing the wave field characteristics of a time-space domain, storing the wave field result after inverse transformation, and finally obtaining the seismic scattering record of the received underground fracture medium.
Preferably, in the step 5, the algorithm for obtaining the encryption grid by interpolation of the wave field through the coarse grid is as follows:
Figure BDA0002625065830000034
the algorithm for transforming the wavefield from the encryption grid to the coarse grid is:
Figure BDA0002625065830000041
in the formula, k is the multiple of grid encryption in the horizontal direction; m is the order of coarse grid points; j is the order of the encrypted mesh points; n is the number of coarse grid points; u shapen mThe wave field value of the m point in the coarse grid; u shapekn kmIs an encrypted grid wave field; u shapekn km+jA wavefield in the encryption grid for km + j points; u shapekn k(m+1)Wavefields in the encryption grid for k (m +1) points;
Figure BDA0002625065830000042
is the wavefield at point l in the encryption grid.
Due to the adoption of the technical scheme, the invention has the following advantages: 1. the method is based on Rytov approximation of weak change of the wave field, the imaginary number field wave field is calculated more stably, and the accumulative error when the Born approximate calculates the wave field is avoided. 2. On the premise of ensuring the precision and the stability, the invention simulates the seismic scattering response characteristics of the fracture stratum by utilizing a multilevel variable grid Green function integration method on the basis of the constructed imaginary number domain Rytov approximate scattered field equation, and the precision can reach millimeter scale. 3. For the same-scale crack model, the simulation efficiency of the multi-stage variable grid Green function integration method is tens of times higher than that of the multi-stage variable grid finite difference method, and the memory consumption of a computer is much lower. 4. The method is not influenced by direct waves, multiple waves and the like, and crack seismic scattering response can be clearly distinguished.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly described below, and the drawings in the following description are only some embodiments of the present invention.
FIG. 1 is a schematic representation of a mesh subdivision of a fracture-bearing formation;
FIG. 2 is a schematic flow diagram of the present invention;
FIG. 3 is a schematic diagram of a fracture-bearing formation velocity model for forward simulation;
FIGS. 4(a) - (f) are surface receiver logs for different fracture dip angles;
fig. 5(a) - (d) are surface reception recordings of fractures at different locations.
Detailed Description
The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, 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 invention.
The basic idea of the invention is that firstly, an analytic expression of a background wave field and a scattering wave field is obtained by approximating a scattering field equation through a virtual number field Rytov, and a Green function integral operator is constructed on the basis; then, deducing a total wave field and a backscattering field phase screen propagation operator of forward propagation of the three-dimensional seismic waves by using the sheet approximation and small angle approximation conditions; finally, the variable grids are adopted to subdivide the underground crack medium (as shown in figure 1), the encryption boundary is determined, the area containing the cracks in the underground crack medium is subjected to forward modeling by adopting the encryption grid Green function integral, and the area without the cracks is subjected to forward modeling by adopting the coarse grid Green function integral.
As shown in fig. 2, the method for simulating the fracture-type formation scattered wave field characteristics with an imaginary number domain Rytov approximation provided by the invention comprises the following specific steps:
step 1: and establishing a ground observation system to determine the position of the ground seismic source excitation point and the receiving position, and designing a proper seismic source signal to be used as an initial wave field of the forward simulation system.
Step 2: performing two-dimensional Fourier transform on the seismic source signal in the real number domain (namely representing the frequency and wave number of the seismic source signal in the real number domain), converting the seismic source signal in the real number domain into the seismic source signal in the imaginary number domain, and taking the seismic source signal in the imaginary number domain as an incident wave field u of the underground fracture mediumi=u0(x, y, z; ω) in which u0Is a background wave field, x, y and z are coordinate axes of a Cartesian coordinate system, and omega is an angular frequency; according to the wave field superposition principle, the total wave field u of the underground fracture medium can be regarded as the incident wave field uiAnd the scattered wave field UsSuperposition of (2): u-ui+Us
And step 3: scattered wave field UsCan be equivalent to a forward scattered field
Figure BDA0002625065830000051
And the back scattering field
Figure BDA0002625065830000052
Linear superposition of (a):
Figure BDA0002625065830000053
the forward scattering field and the back scattering field can be quantitatively characterized on the basis of a scattering wave propagation equation by taking the background velocity of the underground fracture medium as the propagation velocity. In the wave field extrapolation calculation, the Born approximation is mostly adopted in the traditional method, but because the approximation is based on the weak scattering approximation to establish quantitative characterization, the calculation process necessarily generates accumulative errors. Aiming at the problems, the Rytov approximation is introduced to establish a scattered field extrapolation equation. The underground cracks can cause phase change, the Rytov approximation is based on weak change approximation of the scattered field, and in order to improve the stability of calculation, an extrapolation equation of the scattered field is established by the following method:
obtaining the following result according to a Rytov transformation formula:
Figure BDA0002625065830000054
secondly, based on the Rytov approximation, an imaginary number domain Rytov approximation scattering field equation U can be obtainedsComprises the following steps:
Figure BDA0002625065830000055
in the formula, i is an imaginary number unit; s0To be uniformA background slowness; k is a radical ofx、kyAnd kzWave numbers in x, y and z directions, respectively, and
Figure BDA0002625065830000061
wherein k is0Is the background wave number; delta z is the thickness of a thin plate which divides the underground fracture medium into equal thickness along the vertical direction; delta h is the horizontal slowness disturbance quantity between the thin plates;
and thirdly, based on an imaginary number domain Rytov approximate scattered field equation, the forward scattered field and the backward scattered field can be respectively characterized as follows:
Figure BDA0002625065830000062
Figure BDA0002625065830000063
in the formula, FT-1Is a two-dimensional inverse fourier transform operator.
And 4, step 4: firstly, respectively calculating the backscattering field of the upper-layer underground fracture medium by using virtual number domain Rytov approximate scattering wave propagation equations (3) and (4)
Figure BDA0002625065830000064
And forward scattered field
Figure BDA0002625065830000065
And back scattering field to upper underground crack medium
Figure BDA0002625065830000066
Storing; then based on Rytov approximation, the background wave field u of the upper-layer underground fracture medium is determined0(z1) And forward scattered field
Figure BDA0002625065830000067
Substitution into
Figure BDA0002625065830000068
Taking the calculated result as the next underground layerSubstituting the new incident wave field into the virtual number field Rytov approximate scattering wave propagation equations (3) and (4) to respectively calculate the forward scattering field of the next layer of underground fracture medium
Figure BDA0002625065830000069
And a back scattered field
Figure BDA00026250658300000610
And back scattering field to the next layer of underground crack medium
Figure BDA00026250658300000611
And continuing to store.
And 5: calculating the Rytov approximate scattering wave field of the virtual number domain coarse grid according to the step 3 and the step 4, when the simulated wave field is transmitted to the encryption boundary of the underground crack medium, carrying out interpolation encryption processing on the coarse grid wave field in the horizontal direction by using an interpolation principle until the encryption multiple of the grid is met, segmenting the coarse grid in the vertical direction, and carrying out extrapolation calculation on the wave field of the underground crack medium in the encryption grid layer by layer based on a wave field extrapolation theory:
the algorithm for obtaining the encryption grid by coarse grid interpolation of the wave field is as follows:
Figure BDA00026250658300000612
the algorithm for transforming the wavefield from the encryption grid to the coarse grid is:
Figure BDA0002625065830000071
in the formula, k is the multiple of grid encryption in the horizontal direction; m is the order of coarse grid points; j is the order of the encrypted mesh points; n is the number of coarse grid points; u shapen mThe wave field value of the m point in the coarse grid; u shapekn kmIs an encrypted grid wave field; u shapekn km+jA wavefield in the encryption grid for km + j points; u shapekn k(m+1)Wavefields in the encryption grid for k (m +1) points;
Figure BDA0002625065830000072
is the wavefield at point l in the encryption grid.
Step 6: and when the wave field of the underground crack medium in the encryption grid is extrapolated to the encryption boundary of the underground crack medium, transforming the wave field of the encryption grid in the horizontal direction into the coarse grid based on a conversion algorithm, and continuously extrapolating the wave field of the underground crack medium in the coarse grid under the encryption boundary layer by layer based on a wave field extrapolation theory in the vertical direction.
And 7: when the wave field is transmitted to the bottom of the underground fracture medium, the wave field extrapolation of the underground fracture medium is calculated; at this time, from the bottom of the underground fracture medium, the backward scattering field calculated by using the bottom of the underground fracture medium is used as the incident wave field of the underground fracture medium of the upper layer, and the backward propagation backward scattering field is calculated by using the formula (4)
Figure BDA0002625065830000073
And 8: backward propagating back scattered field
Figure BDA0002625065830000074
And (4) carrying out superposition processing on the backward scattering field of the n-1 layer underground fracture medium stored in the step (4), taking the superposition result as an incident wave field of the previous layer underground fracture medium, and calculating a backward propagation forward scattering field by using the formula (3)
Figure BDA0002625065830000075
And step 9: and repeating the step 4-8, backward propagating the wave field, and finishing the backward extrapolation calculation of the wave field of the underground fracture medium when the wave field reaches the receiving position of the underground fracture medium. And finally, converting the imaginary number domain seismic scattering record into a real number domain seismic scattering record by utilizing two-dimensional Fourier inverse transformation, characterizing the wave field characteristics of a time-space domain, storing the wave field result after inverse transformation, and finally obtaining the seismic scattering record of the received underground fracture medium.
FIG. 3 shows a velocity model of a fracture-bearing formation for forward modeling, which is used by the present invention to calculate synthetic seismic records with fracture angles varying from 0-180 (the angle is the angle between the horizontal line and the fracture clockwise), and other parameters are kept unchanged. Partial results are shown in fig. 4(a) - (f), which are synthetic seismic records with fracture angles of 0 °, 30 °, 60 °, 90 °, 120 ° and 150 ° in the order of fig. 4(a) - (f). When the trend of the crack is horizontal, the earthquake in-phase axis generated by the existence of the crack is most continuous, and when the crack is vertical, the continuity is the worst and the energy is relatively concentrated.
5(a) - (d) show the surface reception records of the fractures at different positions of the stratum boundary, wherein the fractures have three groups and seven pieces in each group, and it can be seen from the figure that the surface reception records of the fractures at different positions of the stratum boundary in the model are different. The seismic response of the bed boundaries is increased when the fracture is above the bed reflection boundary and is diminished when the fracture is in the middle and below the bed boundary.
The above embodiments are only for illustrative purposes and are not limited to the above embodiments, and any modifications, equivalent substitutions, improvements and the like within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims (2)

1. A fracture type stratum scattered wave field characteristic simulation method of imaginary number domain Rytov approximation is characterized by comprising the following steps:
step 1: establishing a ground observation system to determine the position of a ground seismic source excitation point and a receiving position, and designing a proper seismic source signal to be used as an initial wave field of a forward simulation system;
step 2: performing two-dimensional Fourier transform on the seismic source signal in the real number domain, converting the seismic source signal in the real number domain into the seismic source signal in the imaginary number domain, and taking the seismic source signal in the imaginary number domain as an incident wave field u of the underground fracture mediumi=u0(x, y, z; ω) in which u0For the background wavefield, x, y and z are the coordinates of a Cartesian coordinate systemAxis, ω is angular frequency;
and step 3: and establishing a scattered field extrapolation equation based on the Rytov approximation:
obtaining the following result according to a Rytov transformation formula:
Figure FDA0002625065820000011
secondly, obtaining an imaginary number domain Rytov approximate scattered field equation U based on the Rytov approximationsComprises the following steps:
Figure FDA0002625065820000012
in the formula, i is an imaginary number unit; s0Uniform background slowness; k is a radical ofx、kyAnd kzWave numbers in x, y and z directions, respectively, and
Figure FDA0002625065820000013
wherein k is0Is the background wave number; delta z is the thickness of a thin plate which divides the underground fracture medium into equal thickness along the vertical direction; delta h is the horizontal slowness disturbance quantity between the thin plates;
and thirdly, based on an imaginary number domain Rytov approximate scattered field equation, the forward scattered field and the back scattered field are respectively characterized as follows:
Figure FDA0002625065820000014
Figure FDA0002625065820000015
in the formula, FT-1A two-dimensional inverse Fourier transform operator;
and 4, step 4: firstly, respectively calculating the backscattering field of the upper-layer underground fracture medium by using virtual number domain Rytov approximate scattering wave propagation equations (3) and (4)
Figure FDA0002625065820000016
And forward scattered field
Figure FDA0002625065820000017
And back scattering field to upper underground crack medium
Figure FDA0002625065820000018
Storing; then based on Rytov approximation, the background wave field u of the upper-layer underground fracture medium is determined0(z1) And forward scattered field
Figure FDA0002625065820000019
Substitution into
Figure FDA00026250658200000110
Taking the calculation result as the incident wave field of the next layer of underground fracture medium, substituting the new incident wave field into the virtual number domain Rytov approximate scattering wave propagation equations (3) and (4), and respectively calculating the forward scattering field of the next layer of underground fracture medium
Figure FDA00026250658200000111
And a back scattered field
Figure FDA00026250658200000112
And back scattering field to the next layer of underground crack medium
Figure FDA00026250658200000113
Continuing to store;
and 5: calculating a virtual number domain coarse grid Rytov approximate scattering wave field according to the step 3 and the step 4, when the simulated wave field is transmitted to an encryption boundary of the underground crack medium, carrying out interpolation encryption processing on the coarse grid wave field in the horizontal direction by using an interpolation principle until the encryption multiple of the grid is met, segmenting the coarse grid in the vertical direction, and carrying out extrapolation calculation on the wave field of the underground crack medium in the encryption grid layer by layer based on a wave field extrapolation theory;
step 6: when the wave field of the underground crack medium in the encryption grid is extrapolated to the encryption boundary of the underground crack medium, transforming the wave field of the encryption grid in the horizontal direction into the coarse grid based on a conversion algorithm, and continuously extrapolating the wave field of the underground crack medium in the coarse grid under the encryption boundary layer by layer based on a wave field extrapolation theory in the vertical direction;
and 7: when the wave field is transmitted to the bottom of the underground fracture medium, the wave field extrapolation of the underground fracture medium is calculated; at this time, from the bottom of the underground fracture medium, the backward scattering field calculated by using the bottom of the underground fracture medium is used as the incident wave field of the underground fracture medium of the upper layer, and the backward propagation backward scattering field is calculated by using the formula (4)
Figure FDA0002625065820000021
And 8: backward propagating back scattered field
Figure FDA0002625065820000022
And (4) carrying out superposition processing on the backward scattering field of the n-1 layer underground fracture medium stored in the step (4), taking the superposition result as an incident wave field of the previous layer underground fracture medium, and calculating a backward propagation forward scattering field by using the formula (3)
Figure FDA0002625065820000023
And step 9: repeating the step 4-8, backward propagating the wave field, and finishing the backward extrapolation calculation of the wave field of the underground fracture medium when the wave field reaches the receiving position of the underground fracture medium; and finally, converting the imaginary number domain seismic scattering record into a real number domain seismic scattering record by utilizing two-dimensional Fourier inverse transformation, characterizing the wave field characteristics of a time-space domain, storing the wave field result after inverse transformation, and finally obtaining the seismic scattering record of the received underground fracture medium.
2. A fracture type formation scattering wave field characteristic simulation method according to claim 1, wherein in the step 5, the algorithm for obtaining the dense grid by coarse grid interpolation of the wave field is as follows:
Figure FDA0002625065820000024
the algorithm for transforming the wavefield from the encryption grid to the coarse grid is:
Figure FDA0002625065820000031
in the formula, k is the multiple of grid encryption in the horizontal direction; m is the order of coarse grid points; j is the order of the encrypted mesh points; n is the number of coarse grid points; u shapen mThe wave field value of the m point in the coarse grid; u shapekn kmIs an encrypted grid wave field; u shapekn km+jA wavefield in the encryption grid for km + j points; u shapekn k(m+1)Wavefields in the encryption grid for k (m +1) points;
Figure FDA0002625065820000032
is the wavefield at point l in the encryption grid.
CN202010794592.4A 2020-08-10 2020-08-10 Fracture type stratum seismic scattering wave field characteristic simulation method based on imaginary number domain Rytov approximation Pending CN111913217A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010794592.4A CN111913217A (en) 2020-08-10 2020-08-10 Fracture type stratum seismic scattering wave field characteristic simulation method based on imaginary number domain Rytov approximation

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010794592.4A CN111913217A (en) 2020-08-10 2020-08-10 Fracture type stratum seismic scattering wave field characteristic simulation method based on imaginary number domain Rytov approximation

Publications (1)

Publication Number Publication Date
CN111913217A true CN111913217A (en) 2020-11-10

Family

ID=73283455

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010794592.4A Pending CN111913217A (en) 2020-08-10 2020-08-10 Fracture type stratum seismic scattering wave field characteristic simulation method based on imaginary number domain Rytov approximation

Country Status (1)

Country Link
CN (1) CN111913217A (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115576014A (en) * 2022-10-26 2023-01-06 江苏科技大学 Intelligent identification method for fractured reservoir based on acoustic wave remote detection imaging

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108508482A (en) * 2018-05-25 2018-09-07 中国海洋石油集团有限公司 A kind of subterranean fracture seismic scattering response characteristic analogy method

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108508482A (en) * 2018-05-25 2018-09-07 中国海洋石油集团有限公司 A kind of subterranean fracture seismic scattering response characteristic analogy method

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
ZHANG KAI 等: "Wave equation tomographic velocity inversion method based on the Born/Rytov approximation", 《APPLIED GEOPHYSICS》 *
侯凯: "非均匀介质地震波正演模拟方法研究", 《中国优秀硕士学位论文全文数据库 基础科学辑》 *
潘成磊 等: "正演模拟倾斜裂缝对应的散射波场", 《中国石油和化工标准与质量》 *
陈生昌 等: "基于Rytov近似的叠前深度偏移方法", 《石油地球物理勘探》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115576014A (en) * 2022-10-26 2023-01-06 江苏科技大学 Intelligent identification method for fractured reservoir based on acoustic wave remote detection imaging
CN115576014B (en) * 2022-10-26 2023-07-11 江苏科技大学 Crack type reservoir intelligent identification method based on acoustic wave far detection imaging

Similar Documents

Publication Publication Date Title
CN111239802B (en) Deep learning speed modeling method based on seismic reflection waveform and velocity spectrum
AU2012220584B2 (en) Sensitivity kernel-based migration velocity analysis in 3D anisotropic media
Vinje et al. Traveltime and amplitude estimation using wavefront construction
CN108508482A (en) A kind of subterranean fracture seismic scattering response characteristic analogy method
CN112883564B (en) Water body temperature prediction method and prediction system based on random forest
CN107765308B (en) Reconstruct low-frequency data frequency domain full waveform inversion method based on convolution thought Yu accurate focus
CN107894618B (en) A kind of full waveform inversion gradient preprocess method based on model smoothing algorithm
RU2451951C2 (en) Method of searching for hydrocarbon deposits confined to fissured-cavernous collectors
WO2016032353A1 (en) Method of searching for hydrocarbon deposits confined to fractured-cavernous reservoirs
CN113031068A (en) Reflection coefficient accurate base tracking prestack seismic inversion method
CN111948708B (en) Seismic wave field forward modeling method for dipping in undulating surface of boundary
CN115600373A (en) Viscous anisotropic medium qP wave simulation method, system, equipment and application
Fan et al. A discontinuous-grid finite-difference scheme for frequency-domain 2D scalar wave modeling
CN111913217A (en) Fracture type stratum seismic scattering wave field characteristic simulation method based on imaginary number domain Rytov approximation
Furumura et al. Parallel PSM/FDM hybrid simulation of ground motions from the 1999 Chi-Chi, Taiwan, earthquake
Yu et al. A minimum traveltime ray tracing global algorithm on a triangular net for propagating plane waves
CN104123449A (en) Subregion local variable-density non-equidistant dual mesh division method for complex mountainous region
Zhao et al. Combined inversion of first‐arrival travel times and reflection travel times
Wu et al. Mesh-free radial-basis-function-generated finite differences and their application in reverse time migration
Deng et al. Hydrocarbon accumulation conditions and key exploration and development technologies for PL 19–3 oilfield
CN109459790A (en) For coal measure strata seismic velocity field method for building up and system
CN107807392A (en) A kind of piecemeal space-time of adaptive anti-frequency dispersion is double to become reverse-time migration method
US20220026593A1 (en) Implicit property modeling
CN116819616B (en) Method for determining thickness of ultrathin high-quality shale reservoir
CN104062680A (en) Method for calculating wave resistance and backstepping gradient of objective function

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
RJ01 Rejection of invention patent application after publication

Application publication date: 20201110

RJ01 Rejection of invention patent application after publication