CN111999770A - Precise beam offset imaging method and system for converting TTI medium into PS wave - Google Patents

Precise beam offset imaging method and system for converting TTI medium into PS wave Download PDF

Info

Publication number
CN111999770A
CN111999770A CN202010915044.2A CN202010915044A CN111999770A CN 111999770 A CN111999770 A CN 111999770A CN 202010915044 A CN202010915044 A CN 202010915044A CN 111999770 A CN111999770 A CN 111999770A
Authority
CN
China
Prior art keywords
wave
anisotropic
medium
shot
point
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
CN202010915044.2A
Other languages
Chinese (zh)
Other versions
CN111999770B (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.)
Chinese Academy of Geological Sciences
Original Assignee
Chinese Academy of Geological Sciences
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 Chinese Academy of Geological Sciences filed Critical Chinese Academy of Geological Sciences
Priority to CN202010915044.2A priority Critical patent/CN111999770B/en
Publication of CN111999770A publication Critical patent/CN111999770A/en
Application granted granted Critical
Publication of CN111999770B publication Critical patent/CN111999770B/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
    • G01V1/30—Analysis
    • G01V1/301—Analysis for determining seismic cross-sections or geostructures
    • 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
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/362—Effecting static or dynamic corrections; Stacking
    • 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/51—Migration
    • G01V2210/512—Pre-stack
    • 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
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/70—Other details related to processing
    • G01V2210/74—Visualisation of seismic data

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 discloses a method and a system for converting a PS wave into a precise beam migration imaging by a TTI medium, which comprises the following steps of firstly, obtaining converted PS wave common shot gather seismic record data, anisotropic medium parameters and migration ray parameters of the TTI medium to be subjected to migration imaging; carrying out ray tracing of P wave and S wave in TTI medium respectively; then calculating the forward wave field of each shot point and the reverse wave field of each wave detection point; finally, imaging the forward wave field of each shot point and the reverse wave field of each wave detection point; and superposing the converted PS wave prestack depth migration profiles of the TTI medium of each shot point to obtain converted PS wave depth domain migration imaging profiles of the TTI medium. The method fully considers the influence of the anisotropy factor on the converted PS wave seismic wave field, can accurately shift and return the underground structure with the TTI medium, overcomes the influence of the anisotropy on the converted PS wave shift, and improves the precision of the acquired prestack depth shift imaging profile of the converted PS wave in the TTI medium.

Description

Precise beam offset imaging method and system for converting TTI medium into PS wave
Technical Field
The invention relates to the technical field of seismic exploration, in particular to a method and a system for accurate beam migration imaging of PS (packet switched) waves converted by TTI (transmission time interval) media.
Background
Compared with single longitudinal wave seismic exploration, the combined conversion of the PS wave seismic data can obtain more underground medium information, and better exploration effects are obtained on underground structure imaging, lithology estimation, fluid detection and reservoir monitoring. Due to the fact that a propagation path from a seismic source to a demodulator probe is asymmetric, imaging processing of PS wave conversion is very difficult, and accurate imaging results cannot be obtained through conventional superposition processing and a post-stack migration method. Although relatively good imaging effect can be obtained by converting the prestack time migration of the PS waves, in areas with complex geological structures and large speed change, imaging errors often exist in the prestack time migration, and prestack depth migration is an effective technical means for solving accurate imaging of the converted PS waves under the condition of the underground complex geological structures.
Most sedimentary rocks in the underground have seismic anisotropy, if a seismic migration method based on isotropic assumption is used for processing seismic data in an anisotropic medium, obvious imaging errors can be caused, and phenomena such as energy unfocusing and the like can occur in a migration section, so that subsequent seismic interpretation work can be seriously influenced. Therefore, with the increasing requirements for seismic exploration precision, it is more and more important to consider the influence of anisotropy on seismic migration imaging. Transversely Isotropic (TI) media are a common model of anisotropic media in practical applications, and most sedimentary rocks can be described as TI media, where TTI media with tilted axes of symmetry are a common form of TI media. The research on the migration technology in the TTI medium has important application value for eliminating the influence of anisotropy in the seismic wave propagation process and realizing high-precision imaging of a complex structure. The influence of the anisotropy on the S wave is more obvious than that of the P wave, and the anisotropy is more important to offset imaging of the converted PS wave and cannot be ignored. At present, the research on the converted PS wave prestack time migration in the anisotropic medium is mature, but the existing prestack time migration method cannot obtain an ideal converted PS wave migration imaging profile for the underground complex structure.
Disclosure of Invention
The invention aims to provide a method and a system for converting a PS wave into a precise beam offset imaging by a TTI medium, so as to overcome the influence of anisotropy on the offset of the converted PS wave and improve the accuracy of an acquired prestack depth offset imaging profile of the converted PS wave in the TTI medium.
In order to achieve the purpose, the invention provides the following scheme:
a TTI medium-converted PS wave precise beam offset imaging method, comprising the following steps:
acquiring converted PS wave common shot gather seismic record data, anisotropic medium parameters and migration ray parameters of a TTI medium to be subjected to migration imaging;
respectively carrying out ray tracing on the P wave and the S wave in the TTI medium by utilizing a two-dimensional anisotropic ray tracing equation system according to the anisotropic medium parameter and the offset ray parameter to obtain the anisotropic Gaussian beam complex value amplitude and complex value time of the P wave emitted by each shot point and the anisotropic Gaussian beam complex value amplitude and complex value time of the S wave emitted by each wave detection point;
constructing a forward wave field of each shot point according to the anisotropic Gaussian beam complex amplitude and the complex time of the P wave emitted by each shot point;
carrying out reverse continuation on the converted PS wave common shot point gather seismic record data of each shot according to the anisotropic Gaussian beam complex value amplitude and the complex value time of the S wave emitted by each wave detection point to obtain a reverse wave field of each wave detection point;
imaging the forward wave field of each shot point and the reverse wave field of each demodulator probe based on the cross-correlation migration imaging condition to obtain a converted PS wave prestack depth migration profile of the TTI medium of each shot point;
and superposing the converted PS wave prestack depth migration profiles of the TTI medium of each shot point to obtain converted PS wave depth domain migration imaging profiles of the TTI medium.
Optionally, the obtaining anisotropic medium parameters of the TTI medium to be subjected to offset imaging specifically includes:
obtaining Thomsen parameters (V) of an anisotropic medium model of a TTI medium to be imaged offsetP0,VS0B), (c)); wherein, VP0、VP0The vertical velocities of the P wave and the S wave are respectively, and the vertical velocities are a first dimensionless factor and a second dimensionless factor which represent the anisotropic strength of the VTI medium;
according to the relation between Thomsen parameters and VTI medium elasticity parameters, using a formula
Figure BDA0002664719550000021
Computing VTI medium density normalization elastic matrix
Figure BDA0002664719550000031
Elastic parameter element of
Figure BDA0002664719550000032
And
Figure BDA0002664719550000033
normalizing elastic parameter elements in an elastic matrix according to VTI medium density
Figure BDA0002664719550000034
Using Bond transformation formula
Figure BDA0002664719550000035
Computing TTI media density normalized elastic matrix
Figure BDA0002664719550000036
The anisotropic elastic parameter element a in11,a33,a55,a13,a15And a35；
Wherein,
Figure BDA0002664719550000037
and
Figure BDA0002664719550000038
are all elastic parameter elements in the VTI medium density normalization elastic matrix,
Figure BDA0002664719550000039
a11、a12、a13、a15、a22、a23、a25、a33、a35、a44、a46、a55、a66are all normalized elastic moments of TTI media densityAnisotropic elastic parameter element in array, theta0Is the symmetry axis tilt angle of the TTI medium.
Optionally, the performing ray tracing on the P-wave and the S-wave in the TTI medium respectively by using a two-dimensional anisotropic ray tracing equation set according to the anisotropic medium parameter and the offset ray parameter to obtain the complex value amplitude and the complex value time of the anisotropic gaussian beam of the P-wave emitted from each shot point and the complex value amplitude and the complex value time of the anisotropic gaussian beam of the S-wave emitted from each wave detection point specifically includes:
solving equation of equation function by using Newton iteration method
Figure BDA0002664719550000041
Obtaining an initial first ray parameter p10And a third ray parameter p30First component p, which is the initial slowness vector of the two-dimensional anisotropic ray tracing equation set1And a third component p3；
Wherein, κ1,κ2,κ3,κ4,κ5Respectively representing a first coefficient, a second coefficient, a third coefficient, a fourth coefficient and a fifth coefficient of the equation function,
Figure BDA0002664719550000042
the first component p of the initial slowness vector1And a third component p3Inputting the initial ray position coordinate into a two-dimensional anisotropic ray tracing equation set, solving the two-dimensional anisotropic ray tracing equation set by using a Runge Kutta method, and obtaining path information and travel time information of central rays of P waves and S waves emitted by each shot point and each wave detection point;
solving an anisotropic dynamic ray equation set according to the path information and the travel time information of the central rays of the P waves and the S waves emitted by each shot point and each wave detection point, and obtaining the dynamic parameters of the central rays of the P waves and the S waves emitted by each shot point and each wave detection point;
and calculating the anisotropic Gaussian beam complex value amplitude and complex value time of the P wave emitted by each shot point and the anisotropic Gaussian beam complex value amplitude and complex value time of the S wave emitted by each wave detection point according to the path information, travel time information and dynamic parameters of the central rays of the P wave and the S wave emitted by each shot point and each wave detection point.
Optionally, the constructing a forward wave field of each shot point according to the complex-valued amplitude and the complex-valued time of the anisotropic gaussian beam of the P wave emitted by each shot point specifically includes:
constructing a forward wave field of each shot seismic source according to the anisotropic Gaussian beam complex amplitude and the complex time of the P wave emitted by each shot point as follows:
Figure BDA0002664719550000043
wherein, G (x, x)sω) is shot point xsThe forward wave field of (a) is,
Figure BDA0002664719550000044
and
Figure BDA0002664719550000045
are respectively shot point xsComplex amplitude and complex time of the emergent P wave anisotropic Gaussian beam, omega is angular frequency,
Figure BDA0002664719550000051
is shot point xsThe emergent P-wave ray is shifted by a parameter vector,
Figure BDA0002664719550000052
shifting the horizontal component of the parameter vector for the P-wave ray;
Figure BDA0002664719550000053
is the vertical component of the P-wave ray migration parameter vector; x represents the position vector of any point in the subsurface.
Optionally, the performing reverse continuation on the converted PS wave common shot gather seismic record data of each shot according to the complex value amplitude and the complex value time of the anisotropic gaussian beam of the S wave emitted from each probe point to obtain a reverse wave field of each probe point specifically includes:
according to the anisotropic Gaussian beam complex amplitude and complex time of S wave emitted from each wave detection point, using formula
Figure BDA0002664719550000054
Carrying out reverse continuation on the converted PS wave common shot point gather seismic record data of each shot to obtain a reverse wave field of each demodulator probe;
wherein R (x, x)rω) denotes the detection point xrOf the backward wave field uPS(xr,xsω) is shot point xsAnd a point x of detectionrThe converted PS wave common shot gather seismic record data,
Figure BDA0002664719550000055
and
Figure BDA0002664719550000056
are respectively the detection point xrThe complex amplitude and complex time of the emergent S-wave anisotropic Gaussian beam,
Figure BDA0002664719550000057
is a detection point xrThe emergent S-wave rays are shifted by a parameter vector,
Figure BDA0002664719550000058
shifting the horizontal component of the parameter vector for the S-wave ray;
Figure BDA0002664719550000059
shifting the vertical component of the parameter vector for the S-wave ray; x represents the position vector of any point in the subsurface.
Optionally, the imaging the forward wave field and the reverse wave field of the seismic source of each shot based on the cross-correlation migration imaging condition to obtain the converted PS wave prestack depth migration profile of the TTI medium of each shot specifically includes:
using a formula based on cross-correlation offset imaging conditions
Figure BDA00026647195500000510
Imaging the forward wave field of each shot point and the reverse wave field of each demodulator probe to obtain a converted PS wave prestack depth migration profile of the TTI medium of each shot point;
wherein,
Figure BDA00026647195500000511
representing shot point xsThe converted PS wave offset profile of (1); sgn (x)r-xs) A sign function representing a correction for the phenomenon of polarity inversion of converted PS-wave seismic recordings,
Figure BDA0002664719550000061
uPS(xr,xsω) is shot point xsAnd a point x of detectionrThe converted PS wave common shot gather seismic record data,
Figure BDA0002664719550000062
and
Figure BDA0002664719550000063
are respectively the detection point xrThe complex amplitude and complex time of the emergent S-wave anisotropic Gaussian beam,
Figure BDA0002664719550000064
and
Figure BDA0002664719550000065
are respectively shot point xsComplex amplitude and complex time of the emitted P-wave anisotropic Gaussian beam, i represents a complex unit,
Figure BDA0002664719550000066
for the horizontal component of the P-wave ray offset parameter vector,
Figure BDA0002664719550000067
shifting the horizontal component of the parameter vector for the S-wave ray; x represents the position vector of any point in the subsurface, and ω is the angular frequency.
Optionally, the superimposing the converted PS wave prestack depth migration profile of the TTI medium of each shot point to obtain the converted PS wave depth domain migration imaging profile of the TTI medium specifically includes:
using formulas
Figure BDA0002664719550000068
Superposing the converted PS wave prestack depth migration profiles of the TTI medium of each shot point to obtain converted PS wave depth domain migration imaging profiles of the TTI medium;
where N represents the number of shots converting a PS-wave common shot gather seismic record, IPS(x) And (3) representing a final TTI medium conversion PS wave depth domain offset imaging profile, wherein x represents a position vector of any point in the underground.
A TTI media switched PS wave precision beam shifting imaging system, the imaging system comprising:
the parameter acquisition module is used for acquiring converted PS wave common shot gather seismic record data, anisotropic medium parameters and migration ray parameters of the TTI medium to be subjected to migration imaging;
the ray tracing module is used for respectively tracing the P wave and the S wave in the TTI medium by utilizing a two-dimensional anisotropic ray tracing equation set according to the anisotropic medium parameter and the offset ray parameter to obtain the anisotropic Gaussian beam complex value amplitude and the complex value time of the P wave emitted by each shot point and the anisotropic Gaussian beam complex value amplitude and the complex value time of the S wave emitted by each wave detection point;
the forward wave field construction module is used for constructing a forward wave field of each shot point according to the anisotropic Gaussian beam complex amplitude and complex time of the P wave emitted by each shot point;
the reverse wave field construction module is used for performing reverse continuation on the converted PS wave common shot point gather seismic record data of each shot according to the anisotropic Gaussian beam complex value amplitude and the complex value time of S waves emitted by each wave detection point to obtain a reverse wave field of each wave detection point;
the single shot point conversion PS wave prestack depth migration profile calculation module is used for imaging the forward wave field of each shot point and the reverse wave field of each demodulator probe based on the cross-correlation migration imaging conditions to obtain a conversion PS wave prestack depth migration profile of the TTI medium of each shot point;
and the conversion PS wave depth domain migration imaging section calculation module is used for superposing the conversion PS wave prestack depth migration sections of the TTI medium of each shot point to obtain the conversion PS wave depth domain migration imaging sections of the TTI medium.
Optionally, the parameter obtaining module specifically includes:
a Thomsen parameter acquisition submodule for acquiring Thomsen parameters (V) of an anisotropic media model of the TTI media to be imaged with offsetP0,VS0B), (c)); wherein, VP0、VP0The vertical velocities of the P wave and the S wave are respectively, and the vertical velocities are a first dimensionless factor and a second dimensionless factor which represent the anisotropic strength of the VTI medium;
an elasticity parameter calculation submodule of the VTI medium, which is used for utilizing a formula according to the relationship between Thomsen parameters and elasticity parameters of the VTI medium
Figure BDA0002664719550000071
Computing VTI medium density normalization elastic matrix
Figure BDA0002664719550000072
Elastic parameter element of
Figure BDA0002664719550000073
And
Figure BDA0002664719550000074
an elasticity parameter calculation submodule of the VTI medium for normalizing the elasticity parameter elements in the elasticity matrix according to the density of the VTI medium
Figure BDA0002664719550000075
Using Bond transformation formula
Figure BDA0002664719550000081
Computing TTI media density normalized elastic matrix
Figure BDA0002664719550000082
The anisotropic elastic parameter element a in11,a33,a55,a13,a15And a35；
Wherein,
Figure BDA0002664719550000083
and
Figure BDA0002664719550000084
are all elastic parameter elements in the VTI medium density normalization elastic matrix,
Figure BDA0002664719550000085
a11、a12、a13、a15、a22、a23、a25、a33、a35、a44、a46、a55、a66are all anisotropic elastic parameter elements in the TTI medium density normalized elastic matrix, theta0Is the symmetry axis tilt angle of the TTI medium.
Optionally, the ray tracing module specifically includes:
an initial parameter acquisition submodule of the two-dimensional anisotropic ray tracing equation set for solving the equation of the equation function by using the Newton iteration method
Figure BDA0002664719550000086
Obtaining an initial first ray parameter p10And a third ray parameter p30First component p, which is the initial slowness vector of the two-dimensional anisotropic ray tracing equation set1And a third component p3；
Wherein, κ1,κ2,κ3,κ4,κ5Respectively representing a first coefficient, a second coefficient, a third coefficient, a fourth coefficient and a fifth coefficient of the equation function,
Figure BDA0002664719550000087
a path information and travel time information acquisition submodule for converting the first component p of the initial slowness vector1And a third component p3Inputting the initial ray position coordinate into a two-dimensional anisotropic ray tracing equation set, solving the two-dimensional anisotropic ray tracing equation set by using a Runge Kutta method, and obtaining path information and travel time information of central rays of P waves and S waves emitted by each shot point and each wave detection point;
the dynamic parameter acquisition submodule is used for solving a dynamic ray equation set according to the path information and the travel time information of the central rays of the P waves and the S waves emitted by each shot point and each wave detection point, and acquiring the dynamic parameters of the central rays of the P waves and the S waves emitted by each shot point and each wave detection point;
and the complex value amplitude and complex value time calculation submodule is used for calculating the anisotropic Gaussian beam complex value amplitude and complex value time of the P wave emitted by each shot point and the anisotropic Gaussian beam complex value amplitude and complex value time of the S wave emitted by each wave detection point according to the path information, travel time information and dynamic parameters of the central rays of the P wave and the S wave emitted by each shot point and each wave detection point.
According to the specific embodiment provided by the invention, the invention discloses the following technical effects:
the invention discloses a TTI medium conversion PS wave accurate beam offset imaging method, which comprises the following steps: acquiring converted PS wave common shot gather seismic record data, anisotropic medium parameters and migration ray parameters of a TTI medium to be subjected to migration imaging; respectively carrying out ray tracing on the P wave and the S wave in the TTI medium by utilizing a two-dimensional anisotropic ray tracing equation system according to the anisotropic medium parameter and the offset ray parameter to obtain the anisotropic Gaussian beam complex value amplitude and complex value time of the P wave emitted by each shot point and the anisotropic Gaussian beam complex value amplitude and complex value time of the S wave emitted by each wave detection point; constructing a forward wave field of each shot point according to the anisotropic Gaussian beam complex amplitude and the complex time of the P wave emitted by each shot point; carrying out reverse continuation on the converted PS wave common shot point gather seismic record data of each shot according to the anisotropic Gaussian beam complex value amplitude and the complex value time of the S wave emitted by each wave detection point to obtain a reverse wave field of each wave detection point; imaging the forward wave field of each shot point and the reverse wave field of each demodulator probe based on the cross-correlation migration imaging condition to obtain a converted PS wave prestack depth migration profile of the TTI medium of each shot point; and superposing the converted PS wave prestack depth migration profiles of the TTI medium of each shot point to obtain converted PS wave depth domain migration imaging profiles of the TTI medium. The method fully considers the influence of the anisotropy factor on the converted PS wave seismic wave field, can accurately shift and return the underground structure with the TTI medium, overcomes the influence of the anisotropy on the converted PS wave shift, and improves the precision of the acquired prestack depth shift imaging profile of the converted PS wave in the TTI medium.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings needed to be used in the embodiments will be briefly described below, and it is obvious that the drawings in the following description are only some embodiments of the present invention, and it is obvious for those skilled in the art to obtain other drawings without inventive exercise.
FIG. 1 is a flowchart of a precise beam offset imaging method for converting PS waves by TTI media according to the present invention;
FIG. 2 is a diagram of a horizontal layered TTI medium model provided in the first embodiment of the present invention;
FIG. 3 is a single shot converted PS wave seismic record of a horizontal layered TTI media model provided in accordance with an embodiment of the present invention;
FIG. 4 is a cross-section of a single shot to PS wave prestack depth migration of a horizontal layered TTI medium model according to an embodiment of the present invention; wherein, fig. 4(a) is an offset profile of a horizontal layered TTI medium model obtained by an isotropic offset method, and fig. 4(b) is an offset profile of a horizontal layered TTI medium model obtained by a method of the present invention;
FIG. 5 is a diagram of a reverse-impact TTI medium model provided in the second embodiment of the present invention;
FIG. 6 is a cross-section of a multi-shot superposition transform PS-wave prestack depth migration of a thrust TTI medium model according to a second embodiment of the present invention; fig. 6(a) is a shifted cross section of the thrust TTI medium model obtained by the isotropic shift method, and fig. 6(b) is a shifted cross section of the thrust TTI medium model obtained by the method of the present invention.
Detailed Description
The invention aims to provide a method and a system for converting a PS wave into a precise beam offset imaging by a TTI medium, so as to overcome the influence of anisotropy on the offset of the converted PS wave and improve the accuracy of an acquired prestack depth offset imaging profile of the converted PS wave in the TTI medium.
In order to make the aforementioned objects, features and advantages of the present invention comprehensible, embodiments accompanied with figures are described in further detail below.
As shown in fig. 1, the present invention provides a TTI medium switched PS wave precise beam offset imaging method, which includes the following steps:
step 101, acquiring converted PS wave common shot gather seismic record data, anisotropic medium parameters and migration ray parameters of a TTI medium to be subjected to migration imaging.
The acquiring of the anisotropic parameter of the TTI medium to be subjected to offset imaging specifically includes:
reading Thomsen parameters (V) of an anisotropic medium model to be offset imagedP0,VS0And) obtaining a VTI medium density normalized elasticity matrix characterized by Thomsen parameters according to the relation between the Thomsen parameters and the VTI medium elasticity parameters
Figure BDA0002664719550000111
Elastic parameter element of
Figure BDA0002664719550000112
The method specifically comprises the following steps:
Figure BDA0002664719550000113
in formula (1), Thomsen parameter VP0、VP0The vertical velocities of the P wave and the S wave are respectively, and are dimensionless factors representing the anisotropic strength of the VTI medium;
obtaining the TTI medium density normalization elastic matrix through Bond transformation according to the elastic parameter of VTI medium density normalization
Figure BDA0002664719550000114
The anisotropic elastic parameter element a in11,a33,a55,a13,a15,a35The method specifically comprises the following steps:
Figure BDA0002664719550000115
Figure BDA0002664719550000116
Figure BDA0002664719550000117
Figure BDA0002664719550000118
Figure BDA0002664719550000119
Figure BDA00026647195500001110
in the formula (2), θ0Is the symmetry axis tilt angle of the TTI medium.
Wherein,
Figure BDA0002664719550000121
and
Figure BDA0002664719550000122
are all elastic parameter elements in the VTI medium density normalization elastic matrix,
Figure BDA0002664719550000123
a11、a12、a13、a15、a22、a23、a25、a33、a35、a44、a46、a55、a66are all anisotropic elastic parameter elements in the TTI medium density normalized elastic matrix, theta0Is the symmetry axis tilt angle of the TTI medium. In the operation of the invention only use
Figure BDA0002664719550000124
And
Figure BDA0002664719550000125
a11、a33、a55、a13、a15and a35And no other elements are needed for calculation.
And step 102, respectively performing ray tracing on the P wave and the S wave in the TTI medium by using a two-dimensional anisotropic ray tracing equation system according to the anisotropic medium parameter and the offset ray parameter, and obtaining the anisotropic Gaussian beam complex value amplitude and the complex value time of the P wave emitted by each shot point and the anisotropic Gaussian beam complex value amplitude and the complex value time of the S wave emitted by each wave detection point.
102, according to the anisotropic medium parameters and the offset ray parameters, respectively performing ray tracing on the P-wave and the S-wave in the TTI medium by using a two-dimensional anisotropic ray tracing equation set, and obtaining the anisotropic gaussian beam complex value amplitude and the complex value time of the P-wave emitted from each shot point and the anisotropic gaussian beam complex value amplitude and the complex value time of the S-wave emitted from each wave detection point, specifically comprising:
the two-dimensional anisotropic ray tracing equation set is:
Figure BDA0002664719550000131
Figure BDA0002664719550000132
Figure BDA0002664719550000133
Figure BDA0002664719550000134
in formula (3), τ is the travel time along the ray, xi(i is 1,3) is a coordinate in a rectangular coordinate system,
Figure BDA0002664719550000135
is a component of the slowness vector, akl＝cklRho (k, l ═ 1,3, or 5) is the elasticity parameter for density normalization of the TTI medium,
Figure BDA0002664719550000136
are the components of the Christoffel matrix eigenvector. The two-dimensional anisotropic ray tracing equation system (equation (3)) has the same form for P-waves and S-waves in a TTI medium, where m-1 and m-2 represent P-waves and S-waves, respectively;
the method for respectively realizing ray tracing of the P wave and the S wave in the TTI medium so as to obtain the complex amplitude and the complex time of the corresponding P wave and S wave anisotropic Gaussian beams specifically comprises the following steps:
(a) according to the equation G of the equation of the P wave and the equation of the S wave in the anisotropic mediumm(xi0,pi0) Obtaining the initial ray parameters p of corresponding wave mode as 110And p30Satisfies the following relation:
Figure BDA0002664719550000137
in the formula (4)
Figure BDA0002664719550000141
κ2＝(a15a33-a13a35)p10
Figure BDA0002664719550000142
Figure BDA0002664719550000143
Figure BDA0002664719550000144
(b) For the initial ray parameters p stated in (a)10And p30The equation of the fourth degree (equation (4)) is solved by using the Newton iteration method, and each initial ray parameter p is10Obtaining corresponding initial ray parameters p30Further, obtaining initial conditions corresponding to P-wave and S-wave anisotropic ray tracing, specifically:
the initial conditions of the two-dimensional anisotropic ray tracing equation set are specifically:
Figure BDA0002664719550000145
in the formula (6), (x)10,x30) Is the coordinates of the initial position of the ray-tracing, which is known for the P-wave and S-wave ray-tracing;
(c) according to the initial conditions corresponding to the P wave and S wave anisotropic ray tracing in the step (b), calculating a two-dimensional anisotropic ray tracing equation set by using a Runge Kutta method, solving the P wave and S wave ray tracing in the corresponding TTI medium, and obtaining the path and travel time information of the corresponding P wave and S wave central rays;
(d) on the paths of the central rays of the P wave and the S wave in the step (c), the dynamic parameters of the corresponding rays are obtained by utilizing an anisotropic dynamic ray equation system, and the method specifically comprises the following steps:
the system of kinetic ray equations in anisotropic media is:
Figure BDA0002664719550000146
in formula (7), P and Q are kinetic parameters of the ray, B1、B2、B3Is the characteristic value G of the corresponding wave-form Christoffel matrixmFor ray normal direction n and normal direction ray parameter pnA derivative of (a);
(e) acquiring complex amplitude A and complex time T of the corresponding P-wave and S-wave anisotropic Gaussian beams based on the travel time information of the P-wave and S-wave central rays in the step (c) and the dynamic parameters of the corresponding rays in the step (d), and specifically:
Figure BDA0002664719550000151
in the formula (8), the reaction mixture is,
Figure BDA0002664719550000152
is the initial amplitude, where V0Is the group velocity, L, of the mode corresponding to the initial position of the ray0Is the initial beam width, ω, of the corresponding mode anisotropic Gaussian beamrIs a reference frequency; v is the group velocity of the corresponding mode along the ray path and n is the vertical distance from the near-ray location to the central ray.
And 103, constructing a forward wave field of each shot point according to the anisotropic Gaussian beam complex amplitude and the complex time of the P wave emitted by each shot point.
103, constructing a forward wave field of each shot point according to the anisotropic gaussian beam complex amplitude and the complex time of the P wave emitted by each shot point, specifically comprising:
reading the seismic source position coordinates of each shot according to the converted PS wave common shot point gather seismic records, and determining the position of the seismic source;
constructing a forward wave field at the seismic source position by using a P wave anisotropic Gaussian beam emitted at the seismic source position, which specifically comprises the following steps:
Figure BDA0002664719550000153
in formula (9), G (x, x)sω) is the seismic source position xsThe forward wave field of (a) is determined,
Figure BDA0002664719550000154
and
Figure BDA0002664719550000155
respectively the complex amplitude and complex time of the P wave anisotropic Gaussian beam emitted by the seismic source, omega is angular frequency,
Figure BDA0002664719550000156
and the P wave ray parameter vector is emitted by the seismic source.
And 104, performing reverse continuation on the converted PS wave common shot gather seismic record data of each shot according to the anisotropic Gaussian beam complex value amplitude and the complex value time of the S wave emitted from each wave detection point to obtain a reverse wave field of each wave detection point.
Step 104, performing reverse continuation on the converted PS wave common shot gather seismic record data of each shot according to the complex value amplitude and the complex value time of the anisotropic gaussian beam of the S wave emitted from each probe point to obtain a reverse wave field of each probe point, specifically including:
the accurate wave field backward continuation formula of the wave detection point is as follows:
Figure BDA0002664719550000161
in the formula (10), R (x, x)rω) represents the backward continuation wave field, uPS(xr,xsω) is the converted PS-wave common shot gather seismic record spectrum in the TTI medium,
Figure BDA0002664719550000162
and
Figure BDA0002664719550000163
are respectively the detection point xrThe complex amplitude and complex time of the emergent S-wave anisotropic Gaussian beam,
Figure BDA0002664719550000164
and the S-wave ray parameter vector is emitted from the wave detection point.
And 105, imaging the forward wave field of each shot point and the reverse wave field of each demodulator probe based on the cross-correlation migration imaging conditions to obtain a converted PS wave prestack depth migration profile of the TTI medium of each shot point.
The imaging formula of the single shot converted PS wave accurate beam offset in the TTI medium is as follows:
Figure BDA0002664719550000165
in the formula (11), the reaction mixture is,
Figure BDA0002664719550000166
representing a single shot converted PS wave migration profile; using sign function sgn (x)r-xs) Correcting the polarity reversal phenomenon of the converted PS wave seismic record, wherein the sign function satisfies the following relation:
Figure BDA0002664719550000167
and 106, superposing the converted PS wave prestack depth migration profiles of the TTI medium of each shot point to obtain converted PS wave depth domain migration imaging profiles of the TTI medium.
Specifically, all the single-shot migration profiles are superposed to obtain a final converted PS wave depth domain migration imaging profile in the TTI medium, which specifically comprises the following steps:
Figure BDA0002664719550000168
in the formula (13), N represents the number of shots of the converted PS-wave common shot gather seismic record, IPS(x) Representing the final TTI media switched PS wave depth domain shifted imaging profile.
The invention also provides a TTI medium conversion PS wave accurate beam offset imaging system, which comprises:
and the parameter acquisition module is used for acquiring the converted PS wave common shot gather seismic record data, anisotropic medium parameters and migration ray parameters of the TTI medium to be subjected to migration imaging.
The parameter obtaining module specifically includes:
a Thomsen parameter acquisition submodule for acquiring Thomsen parameters (V) of an anisotropic media model of the TTI media to be imaged with offsetP0,VS0B), (c)); wherein, VP0、VP0The vertical velocities of the P-wave and the S-wave, respectively, and are a first dimensionless factor and a second dimensionless factor representing the anisotropic strength of the VTI medium.
An elasticity parameter calculation submodule of the VTI medium, which is used for utilizing a formula according to the relationship between Thomsen parameters and elasticity parameters of the VTI medium
Figure BDA0002664719550000171
Computing VTI medium density normalization elastic matrix
Figure BDA0002664719550000172
Elastic parameter element of
Figure BDA0002664719550000173
And
Figure BDA0002664719550000174
an elasticity parameter calculation submodule of the VTI medium for normalizing the elasticity parameter elements in the elasticity matrix according to the density of the VTI medium
Figure BDA0002664719550000175
Using Bond transformation formula
Figure BDA0002664719550000176
Computing TTI media density normalized elastic matrix
Figure BDA0002664719550000181
The anisotropic elastic parameter element a in11,a33,a55,a13,a15And a35。
Wherein,
Figure BDA0002664719550000182
and
Figure BDA0002664719550000183
are all elastic parameter elements in the VTI medium density normalization elastic matrix,
Figure BDA0002664719550000184
a11、a12、a13、a15、a22、a23、a25、a33、a35、a44、a46、a55、a66are all anisotropic elastic parameter elements in the TTI medium density normalized elastic matrix, theta0Is the symmetry axis tilt angle of the TTI medium. In the operation of the invention only use
Figure BDA0002664719550000185
And
Figure BDA0002664719550000186
a11、a33、a55、a13、a15and a35And no other elements are needed for calculation.
And the ray tracing module is used for respectively tracing the P wave and the S wave in the TTI medium by utilizing a two-dimensional anisotropic ray tracing equation set according to the anisotropic medium parameter and the offset ray parameter to obtain the anisotropic Gaussian beam complex value amplitude and the complex value time of the P wave emitted by each shot point and the anisotropic Gaussian beam complex value amplitude and the complex value time of the S wave emitted by each wave detection point.
The ray tracing module specifically comprises:
an initial parameter acquisition submodule of the two-dimensional anisotropic ray tracing equation set for solving the equation of the equation function by using the Newton iteration method
Figure BDA0002664719550000187
Obtaining an initial first ray parameter p10And a third ray parameter p30First component p, which is the initial slowness vector of the two-dimensional anisotropic ray tracing equation set1And a third component p3。
Wherein, κ1,κ2,κ3,κ4,κ5Respectively representing a first coefficient, a second coefficient, a third coefficient, a fourth coefficient and a fifth coefficient of the equation function,
Figure BDA0002664719550000188
a path information and travel time information acquisition submodule for converting the first component p of the initial slowness vector1And a third component p3And inputting the initial ray position coordinate into a two-dimensional anisotropic ray tracing equation set, solving the two-dimensional anisotropic ray tracing equation set by using a Runge Kutta method, and obtaining the path information and the travel time information of the central rays of the P wave and the S wave emitted by each shot point and each wave detection point.
And the dynamic parameter acquisition sub-module is used for solving an anisotropic dynamic ray equation set according to the path information and the travel time information of the central rays of the P wave and the S wave emitted by each shot point and each wave detection point, and acquiring the dynamic parameters of the central rays of the P wave and the S wave emitted by each shot point and each wave detection point.
And the complex value amplitude and complex value time calculation submodule is used for calculating the anisotropic Gaussian beam complex value amplitude and complex value time of the P wave emitted by each shot point and the anisotropic Gaussian beam complex value amplitude and complex value time of the S wave emitted by each wave detection point according to the path information, travel time information and dynamic parameters of the central rays of the P wave and the S wave emitted by each shot point and each wave detection point.
And the forward wave field construction module is used for constructing a forward wave field of each shot point according to the anisotropic Gaussian beam complex amplitude and the complex time of the P wave emitted by each shot point.
And the reverse wave field construction module is used for performing reverse continuation on the converted PS wave common shot gather seismic record data of each shot according to the anisotropic Gaussian beam complex value amplitude and the complex value time of the S wave emitted by each wave detection point to obtain a reverse wave field of each wave detection point.
And the single shot point conversion PS wave prestack depth migration profile calculation module is used for imaging the forward wave field of each shot point and the reverse wave field of each demodulator probe based on the cross-correlation migration imaging conditions to obtain the conversion PS wave prestack depth migration profile of the TTI medium of each shot point.
And the conversion PS wave depth domain migration imaging section calculation module is used for superposing the conversion PS wave prestack depth migration sections of the TTI medium of each shot point to obtain the conversion PS wave depth domain migration imaging sections of the TTI medium.
The first embodiment is as follows:
figure 2 is a model of the horizontal layered TTI medium according to the invention, with model anisotropy parameters as shown in figure 2, with a model grid of 301 x 301 and a longitudinal and transverse grid spacing of 10 m. A single explosion source (shot point) is arranged in the middle of the surface of the model, the source wavelet is a Ricker wavelet, the dominant frequency is 30Hz, the sampling time of seismic records is set to be 3s, and the sampling interval is 2 ms. And a middle blasting and two-side receiving observation system is adopted, and the track interval is 10 m. FIG. 3 is a single shot converted PS wave seismic record of the horizontal layered TTI media model shown in FIG. 2, from which FIG. 3 the phenomenon of polarity reversal of the converted PS wave can be clearly seen. FIG. 4 is a single shot transformed PS wave prestack depth migration profile of the horizontal layered TTI medium model shown in FIG. 2, where FIG. 4(a) is a migration profile obtained using an isotropic migration method and FIG. 4(b) is a migration profile obtained using the method of the present invention. As can be seen from fig. 4(a), the isotropic shift method cannot accurately return the reflection interface in the model due to neglecting the influence of anisotropy, and there is a significant imaging error in the profile. As can be seen from fig. 4(b), the method of the present invention can completely reposition the reflecting interface, resulting in an accurate offset imaging profile. The effectiveness of the method is verified by performing single shot conversion PS wave migration test on the horizontal layered TTI medium model.
The second embodiment is as follows:
fig. 5 shows the thrust TTI medium model provided by the present invention, the anisotropic parameters of the model are as shown in fig. 5, the model grid is 401 x 201, the longitudinal and transverse grid intervals are all 10m, and the model has a thrust rock slice structure with different symmetric axis inclination angles. 77 explosive seismic sources are arranged on the surface of the model, the shot spacing is 50m, the seismic source wavelet is a Ricker wavelet with the main frequency of 25Hz, 401 channels of each shot are received, and the channel spacing is 10 m. FIG. 6 is a multi-shot stack transform PS-wave prestack depth migration profile of the thrust TTI media model shown in FIG. 5, where FIG. 6(a) is a migration profile obtained using the isotropic migration method and FIG. 6(b) is a migration profile obtained using the method of the present invention. As can be seen from fig. 6, for the isotropic offset profile as shown in fig. 6(a), the horizontal interface image underlying the thrust slice structure is lifted upwards, the imaging position of the inclined interface in the thrust slice structure is not accurate, and there is significant divergent energy and strong noise interference near the reflecting interface as shown by the arrows in fig. 6 (a). As can be seen from FIG. 6(b), the method of the present invention obtains accurate focusing imaging for the thrust TTI medium model, and eliminates noise interference, and compared with the isotropic migration imaging effect shown in FIG. 6(a), the method of the present invention has a significant improvement, further verifying that the method of the present invention is an accurate and effective migration method suitable for converting the PS wave seismic data by the TTI medium.
The invention comprises the following steps: acquiring TTI medium conversion PS wave common shot gather seismic record data, TTI medium anisotropy parameters and migration ray parameters; realizing ray tracing of P waves and S waves in the TTI medium, and obtaining complex value amplitude and complex value time of corresponding P waves and S waves anisotropic Gaussian beams; constructing a forward wave field at each shot seismic source by using the obtained P wave anisotropy Gaussian beam complex amplitude and complex time; carrying out reverse continuation on the converted PS wave seismic record of each gun TTI medium to obtain an accurate beam reverse continuation wave field of a corresponding wave detection point of each gun; imaging the forward wave field and the reverse wave field of the wave detection point of each shot seismic source to obtain a prestack depth migration profile of the medium-to-PS wave of each shot TTI; and superposing all the single-shot migration profiles to obtain a final converted PS wave depth domain migration imaging profile in the TTI medium. By adopting the method, the influence of anisotropy on the migration of the converted PS wave can be effectively solved, and the high-precision prestack depth migration imaging profile of the converted PS wave in the TTI medium can be obtained.
According to the specific embodiment provided by the invention, the invention discloses the following technical effects:
1) the influence of anisotropic factors on the converted PS wave seismic wave field is fully considered, and the underground structure with the TTI medium can be accurately shifted and restored; 2) the method is not limited by the strength of anisotropy and is suitable for strong anisotropic media; 3) the invention adopts an accurate beam backward continuation wave field formula, does not perform approximate processing in the wave field continuation calculation process, and can obtain a high-precision conversion PS wave offset imaging profile; 4) the method is not limited by an observation system, and is suitable for irregularly acquired converted PS wave seismic data; 5) the invention can be widely applied to the field of converted PS wave seismic exploration and has more obvious imaging effect on a complex structure with an anisotropic medium.
The embodiments in the present description are described in a progressive manner, each embodiment focuses on differences from other embodiments, and the same and similar parts among the embodiments are referred to each other.
The principle and the implementation manner of the present invention are explained by applying specific examples, the above description of the embodiments is only used to help understanding the method of the present invention and the core idea thereof, the described embodiments are only a part of the embodiments of the present invention, not all embodiments, and all other embodiments obtained by a person of ordinary skill in the art based on the embodiments of the present invention without creative efforts belong to the protection scope of the present invention.

Claims (10)

1. A TTI medium conversion PS wave precise beam offset imaging method is characterized in that the imaging method comprises the following steps:
acquiring converted PS wave common shot gather seismic record data, anisotropic medium parameters and migration ray parameters of a TTI medium to be subjected to migration imaging;
respectively carrying out ray tracing on the P wave and the S wave in the TTI medium by utilizing a two-dimensional anisotropic ray tracing equation system according to the anisotropic medium parameter and the offset ray parameter to obtain the anisotropic Gaussian beam complex value amplitude and complex value time of the P wave emitted by each shot point and the anisotropic Gaussian beam complex value amplitude and complex value time of the S wave emitted by each wave detection point;
constructing a forward wave field of each shot point according to the anisotropic Gaussian beam complex amplitude and the complex time of the P wave emitted by each shot point;
carrying out reverse continuation on the converted PS wave common shot point gather seismic record data of each shot according to the anisotropic Gaussian beam complex value amplitude and the complex value time of the S wave emitted by each wave detection point to obtain a reverse wave field of each wave detection point;
imaging the forward wave field of each shot point and the reverse wave field of each demodulator probe based on the cross-correlation migration imaging condition to obtain a converted PS wave prestack depth migration profile of the TTI medium of each shot point;
and superposing the converted PS wave prestack depth migration profiles of the TTI medium of each shot point to obtain converted PS wave depth domain migration imaging profiles of the TTI medium.
2. The method for precise beam offset imaging of the TTI medium switching PS wave according to claim 1, wherein obtaining anisotropic medium parameters of the TTI medium to be offset imaged specifically comprises:
obtaining Thomsen parameters (V) of an anisotropic medium model of a TTI medium to be imaged offsetP0,VS0B), (c)); wherein, VP0、VS0The vertical velocities of the P wave and the S wave are respectively, and the vertical velocities are a first dimensionless factor and a second dimensionless factor which represent the anisotropic strength of the VTI medium;
according to the relation between Thomsen parameters and VTI medium elasticity parameters, using a formula
Figure FDA0002664719540000011
Computing VTI medium density normalization elastic matrix
Figure FDA0002664719540000021
Elastic parameter element of
Figure FDA0002664719540000022
And
Figure FDA0002664719540000023
normalizing elastic parameter elements in an elastic matrix according to VTI medium density
Figure FDA0002664719540000024
Using Bond transformation formula
Figure FDA0002664719540000025
Computing TTI media density normalized elastic matrix
Figure FDA0002664719540000026
The anisotropic elastic parameter element a in11,a33,a55,a13,a15And a35；
Wherein,
Figure FDA0002664719540000027
and
Figure FDA0002664719540000028
are all elastic parameter elements in the VTI medium density normalization elastic matrix,
Figure FDA0002664719540000029
a11、a12、a13、a15、a22、a23、a25、a33、a35、a44、a46、a55、a66all anisotropic elastic parameter elements in the density normalized elastic matrix of the TTI medium are provided, and theta is the inclination angle of the symmetry axis of the TTI medium.
3. The precise beam offset imaging method for converting the PS wave into the TTI medium according to claim 2, wherein the ray tracing of the P wave and the S wave in the TTI medium is performed respectively by using a two-dimensional anisotropic ray tracing equation system according to the anisotropic medium parameter and the offset ray parameter, so as to obtain the anisotropic gaussian beam complex value amplitude and the complex value time of the P wave emitted from each shot point and the anisotropic gaussian beam complex value amplitude and the complex value time of the S wave emitted from each wave detection point, specifically comprising:
solving equation of equation function by using Newton iteration method
Figure FDA0002664719540000031
Obtaining an initial first ray parameter p10And a third ray parameter p30First component p, which is the initial slowness vector of the two-dimensional anisotropic ray tracing equation set1And a third component p3；
Wherein, κ1,κ2,κ3,κ4,κ5Respectively representing a first coefficient, a second coefficient, a third coefficient, a fourth coefficient and a fifth coefficient of the equation function,
Figure FDA0002664719540000032
the first component p of the initial slowness vector1And a third component p3Inputting the initial ray position coordinate into a two-dimensional anisotropic ray tracing equation set, solving the two-dimensional anisotropic ray tracing equation set by using a Runge Kutta method, and obtaining path information and travel time information of central rays of P waves and S waves emitted by each shot point and each wave detection point;
solving an anisotropic dynamic ray equation set according to the path information and the travel time information of the central rays of the P waves and the S waves emitted by each shot point and each wave detection point, and obtaining the dynamic parameters of the central rays of the P waves and the S waves emitted by each shot point and each wave detection point;
and calculating the anisotropic Gaussian beam complex value amplitude and complex value time of the P wave emitted by each shot point and the anisotropic Gaussian beam complex value amplitude and complex value time of the S wave emitted by each wave detection point according to the path information, travel time information and dynamic parameters of the central rays of the P wave and the S wave emitted by each shot point and each wave detection point.
4. The TTI medium-converted PS wave precise beam migration imaging method according to claim 1, wherein the constructing a forward wave field of each shot according to the anisotropic Gaussian beam complex amplitude and complex time of the P wave emitted by each shot specifically comprises:
constructing a forward wave field of each shot seismic source according to the anisotropic Gaussian beam complex amplitude and the complex time of the P wave emitted by each shot point as follows:
Figure FDA0002664719540000033
wherein, G (x, x)sω) is shot point xsThe forward wave field of (a) is,
Figure FDA0002664719540000041
and
Figure FDA0002664719540000042
are respectively shot point xsComplex amplitude and complex time of the emergent P wave anisotropic Gaussian beam, omega is angular frequency,
Figure FDA0002664719540000043
is shot point xsThe emergent P-wave ray is shifted by a parameter vector,
Figure FDA0002664719540000044
shifting the horizontal component of the parameter vector for P-wave raysAn amount;
Figure FDA0002664719540000045
is the vertical component of the P-wave ray migration parameter vector; x represents the position vector of any point in the subsurface, and i is a complex unit.
5. The TTI medium-to-PS wave precise beam migration imaging method according to claim 1, wherein the method for performing reverse continuation on seismic record data of a gather of common shot points of the converted PS waves for each shot according to the complex value amplitude and the complex value time of the anisotropic Gaussian beam of the S wave emitted from each detection point to obtain a reverse wave field of each detection point specifically comprises the following steps:
according to the anisotropic Gaussian beam complex amplitude and complex time of S wave emitted from each wave detection point, using formula
Figure FDA0002664719540000046
Carrying out reverse continuation on the converted PS wave common shot point gather seismic record data of each shot to obtain a reverse wave field of each demodulator probe;
wherein R (x, x)rω) denotes the detection point xrOf the backward wave field uPS(xr,xsω) is shot point xsAnd a point x of detectionrThe converted PS wave common shot gather seismic record data,
Figure FDA0002664719540000047
and
Figure FDA0002664719540000048
are respectively the detection point xrThe complex amplitude and complex time of the emergent S-wave anisotropic Gaussian beam,
Figure FDA0002664719540000049
is a detection point xrThe emergent S-wave rays are shifted by a parameter vector,
Figure FDA00026647195400000410
shifting the horizontal component of the parameter vector for the S-wave ray;
Figure FDA00026647195400000411
shifting the vertical component of the parameter vector for the S-wave ray; x represents the position vector of any point in the subsurface, i is the complex unit, and ω is the angular frequency.
6. The TTI medium conversion PS wave accurate beam migration imaging method according to claim 1, wherein the imaging is performed on a forward wave field of each shot point and a reverse wave field of each demodulator probe based on a cross-correlation migration imaging condition to obtain a conversion PS wave prestack depth migration profile of the TTI medium of each shot point, and specifically comprises:
using a formula based on cross-correlation offset imaging conditions
Figure FDA00026647195400000412
Imaging the forward wave field of each shot point and the reverse wave field of each demodulator probe to obtain a converted PS wave prestack depth migration profile of the TTI medium of each shot point;
wherein,
Figure FDA0002664719540000051
representing shot point xsThe converted PS wave offset profile of (1); sgn (x)r-xs) A sign function representing a correction for the phenomenon of polarity inversion of converted PS-wave seismic recordings,
Figure FDA0002664719540000052
uPS(xr,xsω) is shot point xsAnd a point x of detectionrThe converted PS wave common shot gather seismic record data,
Figure FDA0002664719540000053
and
Figure FDA0002664719540000054
are respectively the wave detection pointxrThe complex amplitude and complex time of the emergent S-wave anisotropic Gaussian beam,
Figure FDA0002664719540000055
and
Figure FDA0002664719540000056
are respectively shot point xsComplex amplitude and complex time of the emitted P-wave anisotropic Gaussian beam, i represents a complex unit,
Figure FDA0002664719540000057
for the horizontal component of the P-wave ray offset parameter vector,
Figure FDA0002664719540000058
and x is a position vector of any point in the underground, and omega is an angular frequency.
7. The TTI medium conversion PS wave accurate beam offset imaging method according to claim 6, wherein the converting PS wave prestack depth offset profiles of TTI media of each shot point are superposed to obtain a converting PS wave depth domain offset imaging profile of TTI media, which specifically comprises:
using formulas
Figure FDA0002664719540000059
Superposing the converted PS wave prestack depth migration profiles of the TTI medium of each shot point to obtain converted PS wave depth domain migration imaging profiles of the TTI medium;
where N represents the number of shots converting a PS-wave common shot gather seismic record, IPS(x) And (3) representing a final TTI medium conversion PS wave depth domain offset imaging profile, wherein x represents a position vector of any point in the underground.
8. A TTI-media-switched PS-wave precise beam shift imaging system, the imaging system comprising:
the parameter acquisition module is used for acquiring converted PS wave common shot gather seismic record data, anisotropic medium parameters and migration ray parameters of the TTI medium to be subjected to migration imaging;
the ray tracing module is used for respectively tracing the P wave and the S wave in the TTI medium by utilizing a two-dimensional anisotropic ray tracing equation set according to the anisotropic medium parameter and the offset ray parameter to obtain the anisotropic Gaussian beam complex value amplitude and the complex value time of the P wave emitted by each shot point and the anisotropic Gaussian beam complex value amplitude and the complex value time of the S wave emitted by each wave detection point;
the forward wave field construction module is used for constructing a forward wave field of each shot point according to the anisotropic Gaussian beam complex amplitude and complex time of the P wave emitted by each shot point;
the reverse wave field construction module is used for performing reverse continuation on the converted PS wave common shot point gather seismic record data of each shot according to the anisotropic Gaussian beam complex value amplitude and the complex value time of S waves emitted by each wave detection point to obtain a reverse wave field of each wave detection point;
the single shot point conversion PS wave prestack depth migration profile calculation module is used for imaging the forward wave field of each shot point and the reverse wave field of each demodulator probe based on the cross-correlation migration imaging conditions to obtain a conversion PS wave prestack depth migration profile of the TTI medium of each shot point;
and the conversion PS wave depth domain migration imaging section calculation module is used for superposing the conversion PS wave prestack depth migration sections of the TTI medium of each shot point to obtain the conversion PS wave depth domain migration imaging sections of the TTI medium.
9. The system according to claim 8, wherein the parameter obtaining module specifically includes:
a Thomsen parameter acquisition submodule for acquiring Thomsen parameters (V) of an anisotropic media model of the TTI media to be imaged with offsetP0,VS0B), (c)); wherein, VP0、VP0The vertical velocities of the P wave and the S wave respectively, and the vertical velocities of the P wave and the S wave represent the VTI mediumA first dimensionless factor and a second dimensionless factor of the qualitative anisotropic strength;
an elasticity parameter calculation submodule of the VTI medium, which is used for utilizing a formula according to the relationship between Thomsen parameters and elasticity parameters of the VTI medium
Figure FDA0002664719540000061
Computing VTI medium density normalization elastic matrix
Figure FDA0002664719540000062
Elastic parameter element of
Figure FDA0002664719540000063
And
Figure FDA0002664719540000064
an elasticity parameter calculation submodule of the VTI medium for normalizing the elasticity parameter elements in the elasticity matrix according to the density of the VTI medium
Figure FDA0002664719540000071
Using Bond transformation formula
Figure FDA0002664719540000072
Computing TTI media density normalized elastic matrix
Figure FDA0002664719540000073
The anisotropic elastic parameter element a in11,a33,a55,a13,a15And a35；
Wherein,
Figure FDA0002664719540000074
and
Figure FDA0002664719540000075
are all elastic parameter elements in VTI medium density normalization elastic matrixThe content of the element is as follows,
Figure FDA0002664719540000076
a11、a12、a13、a15、a22、a23、a25、a33、a35、a44、a46、a55、a66all anisotropic elastic parameter elements in the density normalized elastic matrix of the TTI medium are provided, and theta is the inclination angle of the symmetry axis of the TTI medium.
10. The TTI media-switched PS wave precise beam offset imaging system of claim 8, wherein the ray tracing module specifically comprises:
an initial parameter acquisition submodule of the two-dimensional anisotropic ray tracing equation set for solving the equation of the equation function by using the Newton iteration method
Figure FDA0002664719540000077
Obtaining an initial first ray parameter p10And a third ray parameter p30First component p, which is the initial slowness vector of the two-dimensional anisotropic ray tracing equation set1And a third component p3；
Wherein, κ1,κ2,κ3,κ4,κ5Respectively representing a first coefficient, a second coefficient, a third coefficient, a fourth coefficient and a fifth coefficient of the equation function,
Figure FDA0002664719540000081
a path information and travel time information acquisition submodule for converting the first component p of the initial slowness vector1And a third component p3Inputting the initial ray position coordinate into a two-dimensional anisotropic ray tracing equation set, solving the two-dimensional anisotropic ray tracing equation set by using a Runge Kutta method, and obtaining path information and travel time information of central rays of P waves and S waves emitted by each shot point and each wave detection point;
the dynamic parameter acquisition submodule is used for solving an anisotropic dynamic ray equation set according to the path information and the travel time information of the central rays of the P waves and the S waves emitted by each shot point and each wave detection point, and acquiring the dynamic parameters of the central rays of the P waves and the S waves emitted by each shot point and each wave detection point;
and the complex value amplitude and complex value time calculation submodule is used for calculating the anisotropic Gaussian beam complex value amplitude and complex value time of the P wave emitted by each shot point and the anisotropic Gaussian beam complex value amplitude and complex value time of the S wave emitted by each wave detection point according to the path information, travel time information and dynamic parameters of the central rays of the P wave and the S wave emitted by each shot point and each wave detection point.
CN202010915044.2A 2020-09-03 2020-09-03 TTI medium conversion PS wave precise beam offset imaging method and system Active CN111999770B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010915044.2A CN111999770B (en) 2020-09-03 2020-09-03 TTI medium conversion PS wave precise beam offset imaging method and system

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010915044.2A CN111999770B (en) 2020-09-03 2020-09-03 TTI medium conversion PS wave precise beam offset imaging method and system

Publications (2)

Publication Number Publication Date
CN111999770A true CN111999770A (en) 2020-11-27
CN111999770B CN111999770B (en) 2024-01-16

Family

ID=73466179

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010915044.2A Active CN111999770B (en) 2020-09-03 2020-09-03 TTI medium conversion PS wave precise beam offset imaging method and system

Country Status (1)

Country Link
CN (1) CN111999770B (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115903032A (en) * 2022-10-25 2023-04-04 中国矿业大学(北京) Underground fluid storage space and migration channel detection method and device
US11953633B1 (en) * 2022-10-24 2024-04-09 China University Of Petroleum (East China) Method, device and computer device for decoupling anisotropic elastic wave

Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB8805489D0 (en) * 1987-03-31 1988-04-07 Amoco Corp Method for depth imaging multicomponent seismic data
US5293352A (en) * 1993-01-07 1994-03-08 Western Atlas International, Inc. Method for removing noise due to near surface scatterers
WO2010014379A2 (en) * 2008-07-30 2010-02-04 Chevron U.S.A. Inc. Method for propagating pseudo acoustic quasi-p waves in anisotropic media
CN102116870A (en) * 2011-02-12 2011-07-06 中国石油大学(华东) Elastic wave gaussian beam pre-stack depth migration technology
CN103995288A (en) * 2014-05-13 2014-08-20 中国石油天然气集团公司 Gauss beam prestack depth migration method and device
CN105549081A (en) * 2016-01-29 2016-05-04 中国石油大学(华东) Anisotropic medium common shot domain Gaussian beam migration imaging method
US20160306058A1 (en) * 2015-04-16 2016-10-20 The Board Of Trustees Of The Leland Stanford Junior Univerisity Reverse time migration based on geometric mean for imaging seismic sources
CN106291687A (en) * 2016-07-21 2017-01-04 中国地质科学院地质研究所 Anisotropy many ripples Gaussian beam pre-stack depth migration imaging method
CN106842300A (en) * 2016-12-21 2017-06-13 中国石油大学(华东) A kind of high efficiency multi-component seismic data true amplitude migration imaging method
CN108363101A (en) * 2018-02-02 2018-08-03 西安石油大学 A kind of inclined shaft crosshole seismic Gaussian beam pre-stack depth migration imaging method
CN109856679A (en) * 2019-03-26 2019-06-07 中国海洋石油集团有限公司 A kind of anisotropic medium elastic wave Gaussian beam offset imaging method and system
US20200124755A1 (en) * 2018-10-22 2020-04-23 Pgs Geophysical As Correction of source motion effects in seismic data recorded in a marine survey using a moving source

Patent Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB8805489D0 (en) * 1987-03-31 1988-04-07 Amoco Corp Method for depth imaging multicomponent seismic data
US5293352A (en) * 1993-01-07 1994-03-08 Western Atlas International, Inc. Method for removing noise due to near surface scatterers
WO2010014379A2 (en) * 2008-07-30 2010-02-04 Chevron U.S.A. Inc. Method for propagating pseudo acoustic quasi-p waves in anisotropic media
CN102116870A (en) * 2011-02-12 2011-07-06 中国石油大学(华东) Elastic wave gaussian beam pre-stack depth migration technology
CN103995288A (en) * 2014-05-13 2014-08-20 中国石油天然气集团公司 Gauss beam prestack depth migration method and device
US20160306058A1 (en) * 2015-04-16 2016-10-20 The Board Of Trustees Of The Leland Stanford Junior Univerisity Reverse time migration based on geometric mean for imaging seismic sources
CN105549081A (en) * 2016-01-29 2016-05-04 中国石油大学(华东) Anisotropic medium common shot domain Gaussian beam migration imaging method
CN106291687A (en) * 2016-07-21 2017-01-04 中国地质科学院地质研究所 Anisotropy many ripples Gaussian beam pre-stack depth migration imaging method
CN106842300A (en) * 2016-12-21 2017-06-13 中国石油大学(华东) A kind of high efficiency multi-component seismic data true amplitude migration imaging method
CN108363101A (en) * 2018-02-02 2018-08-03 西安石油大学 A kind of inclined shaft crosshole seismic Gaussian beam pre-stack depth migration imaging method
US20200124755A1 (en) * 2018-10-22 2020-04-23 Pgs Geophysical As Correction of source motion effects in seismic data recorded in a marine survey using a moving source
CN109856679A (en) * 2019-03-26 2019-06-07 中国海洋石油集团有限公司 A kind of anisotropic medium elastic wave Gaussian beam offset imaging method and system

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
俞岱;杨飞龙;孙渊;边瑞峰;王颖;: "井间地震高斯束叠前深度偏移成像方法", 物探化探计算技术, no. 02, pages 74 - 82 *
江玉乐，雷宛: "各向异性介质地震波传播与成像", 中国石油大学出版社, pages: 74 - 82 *
韩建光等: "二维各向异性多波高斯束叠前深度偏移", 《地球物理学进展》, 15 June 2017 (2017-06-15), pages 1130 - 1139 *
黄建平等: "双复杂条件下弹性波非倾斜叠加精确束偏移方法研究", 石油物探, vol. 54, no. 1, pages 56 - 63 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20240134075A1 (en) * 2022-10-23 2024-04-25 China University Of Petroleum (East China) Method, device and computer device for decoupling anisotropic elastic wave
US11953633B1 (en) * 2022-10-24 2024-04-09 China University Of Petroleum (East China) Method, device and computer device for decoupling anisotropic elastic wave
CN115903032A (en) * 2022-10-25 2023-04-04 中国矿业大学(北京) Underground fluid storage space and migration channel detection method and device
CN115903032B (en) * 2022-10-25 2023-12-15 中国矿业大学(北京) Underground fluid storage space and migration channel detection method and device

Also Published As

Publication number Publication date
CN111999770B (en) 2024-01-16

Similar Documents

Publication Publication Date Title
US10705238B2 (en) Method and apparatus for processing seismic data
US9632192B2 (en) Method of processing seismic data by providing surface offset common image gathers
CN102937721B (en) Limited frequency tomography method for utilizing preliminary wave travel time
CN109856679B (en) Method and system for imaging elastic wave Gaussian beam offset of anisotropic medium
CN106291687A (en) Anisotropy many ripples Gaussian beam pre-stack depth migration imaging method
CN101630016B (en) Method for improving imaging quality of vertical seismic profile
CN103995288B (en) Gauss beam prestack depth migration method and device
US20030018435A1 (en) System for estimating azimuthal variations in seismic data
CN101105537A (en) High accuracy depth domain prestack earthquake data inversion method
CN106970419A (en) A kind of non-homogeneous bent ripple 3D seismic data method for reconstructing based on linear Bregman algorithms
CN102103216A (en) Prestack migration method of two-dimensional Gaussian ray bundle
CN111999770A (en) Precise beam offset imaging method and system for converting TTI medium into PS wave
CN105629299A (en) Travel-time table and angle table acquisition method for angle domain prestack depth migration and imaging method
CN104570116A (en) Geological marker bed-based time difference analyzing and correcting method
CN110780341B (en) Anisotropic seismic imaging method
CN112034520A (en) Anisotropic medium dynamic focusing beam offset imaging method and system
CN106353798A (en) Multi-component joint Gaussian beam pre-stack reverse-time migration imaging method
CN111352153B (en) Microseism interference positioning method based on instantaneous phase cross-correlation weighting
CN105093327B (en) The vector mean filter method of multi-component earthquake data
CN115469362B (en) Energy flow density vector calculation method in seismic exploration
CN111999769B (en) Complex surface anisotropy multicomponent seismic data prestack depth migration method
CN112379431B (en) PS wave seismic data migration imaging method and system under complex surface condition
US20090052279A1 (en) Method and system for determining azimuth parameters for seismic data processing
CN109613614B (en) Method for selecting vertex of VSP (vertical seismic profiling) inclination filter
CN112379430B (en) Multi-component offset imaging method in angle domain

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