CN114839675A - Method for establishing three-dimensional velocity model - Google Patents
Method for establishing three-dimensional velocity model Download PDFInfo
- Publication number
- CN114839675A CN114839675A CN202110132180.9A CN202110132180A CN114839675A CN 114839675 A CN114839675 A CN 114839675A CN 202110132180 A CN202110132180 A CN 202110132180A CN 114839675 A CN114839675 A CN 114839675A
- Authority
- CN
- China
- Prior art keywords
- shot
- point
- ray
- model
- travel time
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 34
- 238000001514 detection method Methods 0.000 claims abstract description 90
- 238000012937 correction Methods 0.000 claims abstract description 10
- 238000007781 pre-processing Methods 0.000 claims abstract description 4
- 239000000523 sample Substances 0.000 claims description 29
- 239000013598 vector Substances 0.000 claims description 13
- 238000004364 calculation method Methods 0.000 claims description 8
- 230000008569 process Effects 0.000 claims description 7
- 238000012360 testing method Methods 0.000 claims description 5
- 239000006185 dispersion Substances 0.000 claims description 4
- 230000009466 transformation Effects 0.000 claims description 4
- 230000001902 propagating effect Effects 0.000 claims description 3
- 238000012545 processing Methods 0.000 abstract description 2
- 230000000644 propagated effect Effects 0.000 abstract 1
- 238000003384 imaging method Methods 0.000 description 10
- 238000003325 tomography Methods 0.000 description 8
- 238000013508 migration Methods 0.000 description 5
- 230000005012 migration Effects 0.000 description 5
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000001427 coherent effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 230000001629 suppression Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
-
- 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/303—Analysis for determining velocity profiles or travel times
- G01V1/305—Travel times
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/622—Velocity, density or impedance
- G01V2210/6222—Velocity; travel time
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 belongs to the technical field of seismic data processing, and provides a method for establishing a three-dimensional velocity model, which comprises the following steps: preprocessing the acquired original three-dimensional seismic data; picking up local related reflected wave travel time and slope data in a common shot point gather and a common detection point gather to form an observation data space; establishing a discrete velocity model, forming a model space, and giving an initial discrete velocity model; in the current speed model, the picked slope data is used for respectively calculating the initial ray directions propagated from the shot point and the wave detection point to the underground, and the ray paths, the travel time and the corresponding ray parameters of the shot point and the wave detection point are obtained by solving a three-dimensional function equation; establishing an inversion target function and an inversion equation, solving to obtain a correction quantity of the velocity model, and updating the velocity model; and judging whether the current speed model meets the requirements, if so, outputting the current speed model, and the method has high efficiency and strong adaptability to complex geological structures and low signal-to-noise ratio data.
Description
Technical Field
The invention belongs to the technical field of seismic data processing, and particularly relates to a method for establishing a three-dimensional velocity model.
Background
The effect of prestack depth migration imaging of seismic survey data is highly dependent on the accuracy of the subsurface velocity model used. Seismic travel time analytic inversion is an effective technology for establishing an underground low-frequency velocity model by utilizing seismic wave travel time data. Conventional reflection wave travel time tomography inversion needs to continuously track each in-phase axis and pick travel time data of each in-phase axis, particularly needs to establish a corresponding relation between the in-phase axes and an underground stratum interface in advance, and is often difficult to realize for complex geological structure areas or low signal-to-noise ratio seismic data. In addition, in the reflected wave travel time tomography inversion, the seismic travel time of the current velocity model and a kernel function matrix in an inversion equation need to be calculated through ray tracing. Determining a ray path from the shot to the demodulator probe requires knowledge of the ray exit angle and incidence angle, in addition to the velocity model, the shot and demodulator probe positions. The conventional reflection time-lapse tomography inversion determines an emergent angle and an incident angle through a plurality of tests, or adopts a mode of firstly solving a wave front and then determining rays by the wave front, so that the calculation amount of ray tracing is large, and the time-lapse tomography inversion efficiency is low.
Disclosure of the invention
The invention aims to solve at least one of the technical problems in the prior art, and provides a method for establishing a three-dimensional velocity model by using seismic wave travel time and slope data, which is easy to realize by computer automation and provides an effective three-dimensional velocity modeling method for seismic imaging of areas with complex structures.
The method for establishing the three-dimensional speed model provided by the embodiment of the invention comprises the following steps:
firstly, preprocessing original three-dimensional seismic data;
picking up earthquake travel time and slope data of a shot point and a demodulator probe to form an observation data space, which specifically comprises the following steps: tracking local related homophase axes by using a slope and travel time picking algorithm of local related homophase axes based on curvelet transformation, respectively tracking the local related homophase axes in a common shot point gather and a common survey point gather, picking corresponding travel time and slope data, and operating three-dimensional seismic data along a main survey line (x direction) and a cross survey line (y direction) respectively to obtain the following observation data space:
wherein ,TSR_obs Reflecting travel time for the observed seismic wave from the shot point S to the wave detection point R; p SX_obs and PSY_obs The observed seismic slopes in the x and y directions at the shot point, respectively; p RX_obs ,、P RY_obs Respectively, the observed seismic slopes in the x and y directions at the waypoints, and N is the number of shot point waypoint pairs, i.e., shot pair number, from which data is picked up.
Step three, establishing a discrete velocity model, forming a model space, and giving an initial discrete velocity model, which specifically comprises the following steps: determining an undulating observation surface according to elevation data of observation system data, and determining the sizes of the three-dimensional speed model in the x direction, the y direction and the z direction according to the range of a study area; determining the size of a discrete grid according to the complexity of an underground structure and the resolution requirement of a speed model to be established; discretizing the discrete velocity model, setting velocity values of all grid nodes, wherein the velocity value of any point in the discrete velocity model is represented by a cubic spline function of the velocity values of the grid nodes, the velocity value of each grid node is a model parameter, namely an unknown quantity, and the model space is as follows:
wherein ,Vj Assigning a value to the speed of each grid point to form an initial speed model, wherein M is the number of the grid nodes and the speed of the jth grid node is the speed value of the jth grid node;
determining the emergent directions of the rays at the shot point and the wave detection point according to the picked shot point slope data, namely obtaining initial values of ray parameters, and then obtaining ray paths, travel time and corresponding ray parameter data of the current speed model by solving a three-dimensional equation, wherein the method specifically comprises the following steps of:
using the picked shot slope data, calculating the emergent direction angle of the ray at the shot and the ray parameter vertical component at the shot by the following formula:
wherein ,
the included angle between the shot point ray direction and the plane passing through the x axis is formed; theta S The included angle between the projection of the ray direction on the plane passing through the x axis and the z axis is shown; v 0S The velocity at the shot point; p SX_obs and PSY_obs Seismic slopes in the x and y directions at the picked shot point, respectively; p SZ0 Is the vertical slowness component at the shot point; sin for medical use -1 Representing an arcsine function; cos represents the cosine function, the initial ray parameter P at the shot point S0 Is composed of
P S0 ＝[P SX_obs ,P SY_obs ,P SZ0 ]
Using the picked slope data of the detection point, the exit direction angle of the ray at the detection point and the ray parameter vertical component at the detection point are calculated by the following formula:
wherein ,
the included angle between the ray direction of the wave detection point and the plane passing through the x axis is shown; theta R The included angle between the projection of the ray direction on the plane passing through the x axis and the z axis is shown; v 0R Is the velocity at the pickup point; p RX_obs and PRY_obs Seismic slopes in x and y directions at the picked-up wave detection points, respectively; p RZ0 Is the vertical slowness component at the detection point; sin for medical use -1 Representing an arcsine function; cos denotes a cosine function. At this time, the initial ray parameter P at the detection point R0 Is composed of
P R0 ＝[P RX_obs ,P RY_obs ,P RZ0 ]
Starting from the shot point S and the demodulator probe R respectively, solving the following three-dimensional equation by using the initial ray parameters of the shot point and the demodulator probe obtained in the steps S3 and S4,
wherein X ═ (X, Y, Z) is a spatial position vector; p (x) ═ P X (x),P Y (x),P Z (x)]The vector is composed of ray parameters or slowness components in three coordinate axis directions; v (x) is the velocity at x; t is slave gunSeismic travel time of points or geophone points; ^ represents the gradient operator.
Solving the differential equation to continuously obtain the path of the ray from the shot point S to the underground and the single-pass travel time t S And the path of the ray propagating from the point of detection R into the subsurface and the single travel time t R When t is S And t R The sum of which is equal to the observed two-way travel time T of the shot-check pair SR_obs Stopping calculation, and recording the ray end point position of the shot point S as x S And the ray end point position of the wave detection point R is recorded as x R Storing the ray length and the ray parameter of the ray of the shot point and the wave detection point in each corresponding grid unit;
and step five, calculating errors of the reflection points determined from the shot points and the wave detection points respectively and errors of observed reflection travel time and calculated reflection travel time, judging whether the current velocity model meets the requirements or not, and determining whether the inversion process is terminated or not.
Further, the fifth step of calculating the errors of the reflection points and the errors of the observed reflection travel time and the calculated reflection travel time respectively determined from the shot point and the demodulator probe, judging whether the current velocity model meets the requirements, and determining whether the inversion process is terminated, specifically comprises the following steps:
s1, if x S And x R When the single-pass travel time sum of the rays respectively along the shot point and the demodulator probe is consistent with the observed double-pass reflection travel time, the used discrete velocity model is correct, the distance between the tail end positions of the rays of the shot point and the demodulator probe of the shot-geophone pair is used for representing the error of the reflection point determined from the shot point and the demodulator probe, and the error of the reflection point and the error of the reflection travel time are respectively calculated by the following formulas:
wherein ,(XSi ，Y Si ，Z Si) and (XRi ，Y Ri ，Z Ri ) Coordinates of ray ends of a shot point and a demodulator probe of the ith shot-check pair are respectively; t is t Si and tRi Respectively the travel time from the shot point and the wave detection point of the ith shot-detection pair to the tail end of each ray; t is SRi obs The reflected wave travel time of the ith shot-checking pair is observed;
s2, setting error limit, calculating err D and errT When the error is smaller than the given error limit, outputting the current speed model and stopping running; otherwise, executing the following steps;
s3, establishing a specific objective function and an inversion equation, solving to obtain a correction quantity of the current speed model, and updating the speed model;
s4, establishing the following objective function:
eyes of a user
D SRi ＝[(X Si -X Ri ) 2 +(Y Si -Y Ri ) 2 +(Z Si -Z Ri ) 2 ] 1/2
T SRi ＝t Si +t Ri -T SRi_obs
Wherein N is the number of shot detection pairs; t is t Si and tRi Respectively the travel time from the shot point and the wave detection point of the ith shot-detection pair to the tail end of each ray; t is SRi_obs The reflected wave of the ith shot-checking pair is subjected to double-pass travel; t is SRi The difference between the sum of two single-pass travel times of the ith shot-checking pair and the observed two-pass travel time; x is the number of Si ＝(X Si ，Y Si ，Z Si) and xRi ＝(X Ri ，Y Ri ，Z Ri ) Coordinate vectors of ray tail ends of shot points and wave detection points of the ith shot and survey pair respectively; d SRi Is the ray end point x Si And x Ri The distance between them; v is a model parameter vector consisting of mesh node velocities; λ is a weight coefficient;
s5, solving the minimum value solution of the nonlinear objective function by adopting a Gaussian-Newton method, performing Taylor series expansion on the objective function to obtain a first-order term, and obtaining a linear inversion equation as follows:
G·Δv＝Δd
wherein ,
B ij ＝λ(L Sij +L Rij )，i＝1，2，...，N，j＝1，2，...，M；
Δd＝[D SR1 D SR2 … D SRN T SR1 T SR2 … T SRN ] T ；
Δv＝[V 1 V 2 … V M ] T ；
D sRi ＝[(X Si -X Ri ) 2 +(Y Si -Y Ri ) 2 +(Z Si -Z Ri ) 2 ] 1/2 ，i＝1，2，...，N；
T SRi ＝t Si +t Ri -T SRi_obs ，i＝1，2，...，N；
t Si and tRi Respectively the travel time from the shot point and the wave detection point of the ith shot-detection pair to the tail end of each ray; t is SRi_obs Travel time for the reflected wave of the ith shot-detection pair; (X) Si ，Y Si ，Z Si) and (XRi ，Y Ri ，Z Ri ) Coordinates of ray tail ends of a shot point and a wave detection point of the ith shot and survey pair respectively; v i Is the velocity value of the ith discrete grid; v Si and VRi The speeds of the unit where the shot point and the wave detection point ray tail end of the ith shot and survey pair are located respectively; p SXi ，P SYi and PSZi The shot points of the ith shot detection pair respectivelyRay parameters of the ray end in x, y and z directions; p RXi ，P RYi and PRZi Ray parameters of the ray tail end of the wave detection point of the ith shot-examination pair along the x direction, the y direction and the z direction respectively; l is Sij and LRij The lengths of ray segments of shot point rays and wave detection point rays of the ith shot and survey pair in the jth unit are respectively set; l is SXij ，L SYij and LSZij Projection lengths of ray segments of shot point rays of the ith shot and examine pair in the jth unit in the x direction, the y direction and the z direction are respectively set; l is RXij ，L RYij and LRZij Projection lengths of ray segments of wave detection point rays of the ith shot-examination pair in the jth unit in the directions of x, y and z are respectively set; n is the total number of the shot-test pairs; m is the total number of grid nodes of the discrete model;
s6, solving the linear inversion equation by using an LSQR algorithm to obtain a discrete velocity model updating quantity delta v, modifying the current discrete velocity model by using the model updating quantity to obtain a modified model, wherein the calculation formula is as follows:
v (k+1) ＝v (k) +Δv
wherein ,v(k) The current discrete velocity model is the discrete velocity model of the last iteration (the kth time), and the delta v is the correction quantity of the discrete velocity model obtained by the current iteration inversion; v. of (k+1) And (4) taking the corrected discrete velocity model obtained by the iteration as a new initial model, and turning to the step four.
The invention has the beneficial effects that: by adopting the method for establishing the three-dimensional velocity model by using the travel time and the slope data of the seismic reflection waves, the slope data of the local related reflection wave in-phase axis is used for calculating the propagation direction from the shot point and the wave detection point to the underground and the initial ray parameters, so that the efficiency of calculating the ray path and the travel time is improved; the in-phase axis of the used local related reflected wave is easy to automatically track, the travel time and slope data of the in-phase axis are easy to automatically pick up, and the in-phase axis local related reflected wave has strong adaptability to complex underground structures and low signal-to-noise ratio seismic data. The embodiment of the invention is easy to realize by computer automation, and provides an effective three-dimensional velocity modeling method for seismic imaging of areas with complex structures.
Drawings
The invention is further described below with reference to the accompanying drawings and examples;
FIG. 1 is a flow chart of a method for building a three-dimensional velocity model according to an embodiment of the present invention;
FIG. 2 is a slice view of a three-dimensional velocity model created using travel time and slope data for the embodiment of FIG. 1;
FIG. 3 is a slice diagram of a three-dimensional velocity model built by a conventional tomographic inversion method using travel-time data;
FIG. 4 is a diagram of a common imaging trace set after the three-dimensional velocity model constructed according to the embodiment of the present invention is shifted;
FIG. 5 is a common imaging point trace set plot after migration of a three-dimensional velocity model built using conventional time-lapse tomographic inversion.
Detailed Description
In order to make the purpose and technical solution of the present invention more apparent, the present invention is further described in detail below with reference to the accompanying drawings.
In seismic data, a common shot gather and a common geophone gather easily track a local correlation in-phase axis, and the change rate of travel time of the local correlation in-phase axis along with offset is called seismic slope. According to ray theory, the seismic slope is a ray parameter, i.e., the horizontal component of the slowness vector. If the velocity values at the shot point and the demodulator point are known, the slope parameters can be used to calculate the exit angle of the reflected ray at the shot point and the incident angle at the demodulator point. In addition, the local coherent in-phase axis is easy to track, the travel time and the slope are easy to pick up, and the method has strong adaptability to the seismic data with low signal-to-noise ratio. If a smooth velocity model is adopted, rays are constrained through the seismic slope, only locally related event travel time and slope data are needed, and the corresponding relation between each event and a stratum interface does not need to be established in advance, so that the automation degree and the practicability of reflection wave travel time tomography inversion are improved. Therefore, if the seismic slope data is introduced into the reflection wave time-lapse tomography inversion, the ray emergence angle and the ray incidence angle are determined by the seismic slope, the ray tracing efficiency can be improved, and a practical method for establishing an underground low-frequency velocity model is formed, so that a low-frequency or initial velocity model is provided for prestack depth migration velocity analysis and waveform inversion.
The embodiment of the invention introduces seismic slope data into three-dimensional seismic reflection travel time analytic inversion, determines a ray emergence angle and an incident angle by using the seismic slope, provides a method for establishing a three-dimensional velocity model by using seismic reflection travel time and slope data according to the condition that rays respectively transmitted from a shot point and a wave detection point to the underground meet at a reflection point when a velocity model is accurate, and the sum of the travel time of the two rays is necessarily equal to the actually measured two-way reflection travel time, and provides a technology for establishing the underground three-dimensional velocity model for seismic exploration in a complex construction area or offset imaging of low signal-to-noise ratio seismic data.
Referring to fig. 1, a method for building a three-dimensional velocity model according to an embodiment of the present invention includes the following steps:
firstly, preprocessing original three-dimensional seismic data;
picking up earthquake travel time and slope data of a shot point and a demodulator probe to form an observation data space, which specifically comprises the following steps: tracking local related homophase axes by using a slope and travel time picking algorithm of local related homophase axes based on curvelet transformation, respectively tracking the local related homophase axes in a common shot point gather and a common survey point gather, picking corresponding travel time and slope data, and operating three-dimensional seismic data along a main survey line (x direction) and a cross survey line (y direction) respectively to obtain the following observation data space:
wherein ,TSR_obs Reflecting travel time for the observed seismic wave from the shot point S to the wave detection point R; p SX_obs and PSY_obs Observed seismic slopes in the x and y directions at the shot point, respectively; p RX_obs ,、P RY_obs Respectively, the observed seismic slopes in the x and y directions at the waypoints, and N is the number of shot point waypoint pairs, i.e., shot pair number, from which data is picked up.
Step three, establishing a discrete velocity model, forming a model space, and giving an initial discrete velocity model, which specifically comprises the following steps: determining an undulating observation surface according to elevation data of observation system data, and determining the sizes of the three-dimensional speed model in the x direction, the y direction and the z direction according to the range of a study area; determining the size of a discrete grid according to the complexity of an underground structure and the resolution requirement of a speed model to be established; discretizing the discrete velocity model, setting velocity values of all grid nodes, wherein the velocity value of any point in the discrete velocity model is represented by a cubic spline function of the velocity values of the grid nodes, the velocity value of each grid node is a model parameter, namely an unknown quantity, and the model space is as follows:
wherein ,Vj Assigning a value to the speed of each grid point to form an initial speed model, wherein M is the number of the grid nodes and the speed of the jth grid node is the speed value of the jth grid node;
determining the emergent directions of the rays at the shot point and the wave detection point according to the picked shot point slope data, namely obtaining initial values of ray parameters, and then obtaining ray paths, travel time and corresponding ray parameter data of the current speed model by solving a three-dimensional equation, wherein the method specifically comprises the following steps of:
using the picked shot slope data, calculating the emergent direction angle of the ray at the shot and the ray parameter vertical component at the shot by using the following formula:
wherein ,
the included angle between the shot point ray direction and the plane passing through the x axis is formed; theta S The included angle between the projection of the ray direction on the plane passing through the x axis and the z axis is shown; v 0S The velocity at the shot point; p SX_obs and PSY_obs Seismic slopes in the x and y directions at the picked shot point, respectively; p SZ0 Is the vertical slowness component at the shot point; sin for medical use -1 Representing an arcsine function; cos represents the cosine function, the initial ray parameter P at the shot point S0 Is composed of
P S0 ＝[P SX_obs ,P SY_obs ,P SZ0 ]
Using the picked slope data of the detection point, the exit direction angle of the ray at the detection point and the ray parameter vertical component at the detection point are calculated by the following formula:
wherein ,
the included angle between the ray direction of the wave detection point and the plane passing through the x axis is shown; theta R The included angle between the projection of the ray direction on the plane passing through the x axis and the z axis is shown; v 0R Is the velocity at the pickup point; p RX_obs and PRY_obs Seismic slopes in x and y directions at the picked-up wave detection points, respectively; p RZ0 Is the vertical slowness component at the detection point; sin for medical use -1 Representing an arcsine function; cos denotes a cosine function. At this time, the initial ray parameter P at the detection point R0 Is composed of
P R0 ＝[P RX_obs ,P RY_obs ,P RZ0 ]
Starting from the shot point S and the demodulator probe R respectively, solving the following three-dimensional equation by using the initial ray parameters of the shot point and the demodulator probe obtained in the steps S3 and S4,
wherein X ═ (X, Y, Z) is a spatial position vector; p (x) ═ P X (x),P Y (x),P Z (x)]The vector is composed of ray parameters or slowness components in three coordinate axis directions; v (x) is the velocity at x; t is the seismic travel time from the shot point or the demodulator probe;
a gradient operator is represented.
Solving the differential equation to continuously obtain the path of the ray from the shot point S to the underground and the single-pass travel time t S And the path of the ray propagating from the point of detection R into the subsurface and the single travel time t R When t is S And t R The sum of which is equal to the observed two-way travel time T of the shot-check pair SR_obs Stopping calculation, and recording the ray end point position of the shot point S as x S And the ray end point position of the wave detection point R is recorded as x R Storing the ray length and the ray parameter of the ray of the shot point and the wave detection point in each corresponding grid unit;
and step five, calculating errors of the reflection points determined from the shot points and the wave detection points respectively and errors of observed reflection travel time and calculated reflection travel time, judging whether the current velocity model meets the requirements or not, and determining whether the inversion process is terminated or not.
Further, the fifth step of calculating the errors of the reflection points and the errors of the observed reflection travel time and the calculated reflection travel time respectively determined from the shot point and the demodulator probe, judging whether the current velocity model meets the requirements, and determining whether the inversion process is terminated, specifically comprises the following steps:
s1, e.g.Fruit x S And x R When the single-pass travel time sum of the rays respectively along the shot point and the demodulator probe is consistent with the observed double-pass reflection travel time, the used discrete velocity model is correct, the distance between the tail end positions of the rays of the shot point and the demodulator probe of the shot-geophone pair is used for representing the error of the reflection point determined from the shot point and the demodulator probe, and the error of the reflection point and the error of the reflection travel time are respectively calculated by the following formulas:
wherein ,(XSi ，Y Si ，Z Si) and (XRi ，Y Ri ，Z Ri ) Coordinates of ray ends of a shot point and a demodulator probe of the ith shot-check pair are respectively; t is t Si and tRi Respectively the travel time from the shot point and the wave detection point of the ith shot-detection pair to the tail end of each ray; t is SRi_obs The reflected wave travel time of the ith shot-checking pair is observed;
s2, setting error limit, calculating err D and errT When the error is smaller than the given error limit, outputting the current speed model and stopping running; otherwise, executing the following steps;
s3, establishing a specific objective function and an inversion equation, solving to obtain the correction quantity of the current speed model, and updating the speed model;
s4, establishing the following objective function:
and is
D SRi ＝[(X Si -X Ri ) 2 +(Y Si -Y Ri ) 2 +(Z Si -Z Ri ) 2 ] 1/2
T SRi ＝t Si +t Ri -T SRi_obs
Wherein N is the number of shot detection pairs; t is t Si and tRi Respectively the travel time from the shot point and the wave detection point of the ith shot-detection pair to the tail end of each ray; t is SRi_obs The reflected wave of the ith shot-checking pair is subjected to double-pass travel; t is SRi The difference between the sum of two single-pass travel times of the ith shot-checking pair and the observed two-pass travel time; x is the number of Si ＝(X Si ，Y Si ，Z Si) and xRi ＝(X Ri ，Y Ri ，Z Ri ) Coordinate vectors of ray tail ends of shot points and wave detection points of the ith shot and survey pair respectively; d SRi Is the ray end point x Si And x Ri The distance between them; v is a model parameter vector consisting of mesh node velocities; λ is a weight coefficient;
s5, solving the minimum value solution of the nonlinear objective function by adopting a Gaussian-Newton method, performing Taylor series expansion on the objective function to obtain a first-order term, and obtaining a linear inversion equation as follows:
G·Δv＝Δd
wherein ,
B ij ＝λ(L Sij +L Rij )，i＝1，2，...，N，j＝1，2，...，M；
Δd＝[D SR1 D SR2 … D SRN T SR1 T SR2 … T SRN ] T ；
Δv＝[V 1 V 2 … V M ] T ；
D SRi ＝[(X Si -X Ri ) 2 +(Y Si -Y Ri ) 2 +(Z Si -Z Ri ) 2 ] 1/2 ，i＝1，2，...，N；
T SRi ＝t Si +t Ri -T SRi_obs ，i＝1，2，...，N；
t Si and tRi Respectively the travel time from the shot point and the wave detection point of the ith shot-detection pair to the tail end of each ray; t is SRi_obs The reflected wave travel time of the ith shot-checking pair is observed; (X) Si ，Y Si ，Z Si) and (XRi ，Y Ri ，Z Ri ) Coordinates of ray tail ends of a shot point and a wave detection point of the ith shot and survey pair respectively; v i Is the velocity value of the ith discrete grid; v Si and VRi The speeds of the unit where the shot point and the wave detection point ray tail end of the ith shot and survey pair are located respectively; p SXi ，P SYi and PSZi Ray parameters of the shot point ray tail end of the ith shot and survey pair along the x, y and z directions respectively; p RXi ，P RYi and PRZi Ray parameters of the ray tail end of the wave detection point of the ith shot-examination pair along the x direction, the y direction and the z direction respectively; l is Sij and LRij The lengths of ray segments of shot point rays and wave detection point rays of the ith shot and survey pair in the jth unit are respectively set; l is SXij ，L sYij and LSZij Projection lengths of ray segments of shot point rays of the ith shot and examine pair in the jth unit in the x direction, the y direction and the z direction are respectively set; l is RXij ，L RYij and LRZij Projection lengths of ray segments of wave detection point rays of the ith shot-examination pair in the jth unit in the directions of x, y and z are respectively set; n is the total number of the shot-test pairs; m is the total number of grid nodes of the discrete model;
s6, solving the linear inversion equation by using an LSQR algorithm to obtain a discrete velocity model updating quantity delta v, modifying the current discrete velocity model by using the model updating quantity to obtain a modified model, wherein the calculation formula is as follows:
v (k+1) ＝v (k) +Δv
wherein ,v(k) Is the current discrete velocity model, i.e. the discrete velocity model of the last (k) iteration, and Δ v is the distance obtained by the inversion of the current iterationA dispersion velocity model correction amount; v. of (k+1) And (4) taking the corrected discrete velocity model obtained by the iteration as a new initial model, and turning to the step four.
As a specific example of the embodiment of the invention, the method of the embodiment of the invention is tested by using three-dimensional oil and gas seismic exploration data in a certain area in the southwest of China and compared with the result of reflection tomography only using reflection travel time data.
And step 100, the length of the work area is 16km, and the width of the work area is 8 km. And 25 gun lines and 59 demodulation point lines are uniformly distributed, wherein the distance between the gun lines is 300m, the distance between the demodulation lines is 150m, the distance between the guns is 100m, and the distance between the demodulation points is 50 m. Each shot is fired with 8 lines of detectors and 240 detector points per line of detectors. The original seismic data is preprocessed by noise suppression, amplitude compensation, deconvolution, static correction and the like.
And step 110, acquiring time and slope data in a common shot point domain and a common geophone point domain by using a local correlation same-phase axis slope and travel time picking method based on curvelet transformation, and acquiring reflection travel time and slope data of 236152 shot-geophone pairs to form a data space.
Step 120, a velocity model with a size of 16km × 8km × 4km is established, and the three-dimensional velocity model is discretized by a grid with a size of 250m × 250m × 250m to form a model space. The model surface uses the corresponding elevation data to construct the relief. The initial speed of the model is set to linearly increase along with the depth, the speed at the highest position of the surface is 1.0km/s, and the speed at the bottom boundary of the model is 4.5km/s, so that an initial speed model is formed.
And step 130, acquiring the directions of the seismic rays respectively transmitted from the shot point and the demodulator probe to the underground by using the initial speed and the picked slope data of the shot point and the demodulator probe, and giving initial conditions for ray emission.
And step 140, solving a three-dimensional equation by adopting a Runge-Kutta method according to the initial velocity model and the initial ray conditions to obtain the ray tail end position and the travel time of each shot-inspection pair and the corresponding ray path and ray parameters of each grid unit.
And 150, calculating a reflection point error by using the distance between the shot point of each shot-detection pair and the tail end position of the ray of the wave detection point, and calculating a reflection travel time error by using the travel time of the shot point of each shot-detection pair, the tail end of the ray of the wave detection point and the actually measured travel time of the reflected wave.
And step 160, comparing the calculated reflection point error and the calculated reflection travel time error with preset error limits respectively, and judging whether the current speed model meets the requirements or not. If the calculated reflection point error and the calculated reflection travel time error are both smaller than the set error limit, outputting a current speed model and terminating the calculation; otherwise, the following steps are performed.
Step 170, establishing an inversion equation by using the observed reflection travel time and slope data, the ray end positions and the travel times of each shot-test pair calculated in step S5, and the corresponding ray paths and ray parameters of each grid cell. And solving an inversion equation to obtain the correction quantity of the current speed model.
And step 180, correcting the current model by using the obtained model correction quantity to obtain an updated model. The updated model is used as a new initial model, and the process proceeds to step S4.
FIG. 2 is a vertical slice along the mean line of the subsurface velocity model output after 30 iterations of the present invention, and FIG. 3 is a vertical slice along the mean line of the subsurface velocity model constructed using only reflection-time data. Comparing the two models, the velocity model established by the invention using the travel time and slope data is greatly different from the conventional velocity model established by using only travel time data in the transverse direction and the vertical direction.
In order to further illustrate the beneficial effects of the invention, the prestack depth migration imaging is carried out on the seismic data by respectively utilizing the velocity model established by using the travel time and the slope data and the conventional velocity model established by using the travel time data only. Fig. 4 shows that the common imaging point gather along the middle survey line is obtained by shifting the velocity model established by the invention, and fig. 5 shows that the common imaging point gather along the middle survey line is obtained by shifting the velocity model established by conventional travel time tomography inversion. The pulling-down and lifting phenomena of the in-phase axis in fig. 5 are more obvious, the in-phase axis basically presents a non-horizontal form, while the in-phase axis in fig. 4 is relatively horizontal and smooth. In co-imaging trace sets, the accuracy of the velocity model used in migration is often determined by the level of the in-phase axis. The more the in-phase axis goes to the horizontal, the more accurate the velocity model used to account for the offset. Comparing fig. 4 and fig. 5, it can be seen that the velocity model established by the present invention is more accurate, which shows the good application effect of the present invention.
The above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not to be construed as limiting the invention, but rather the intention is to cover all modifications, equivalents, and alternatives falling within the spirit and scope of the invention as defined by the appended claims.
Claims (2)
1. A method of building a three-dimensional velocity model using seismic data, comprising the steps of:
firstly, preprocessing original three-dimensional seismic data;
picking up earthquake travel time and slope data of a shot point and a demodulator probe to form an observation data space, which specifically comprises the following steps: tracking local related homophase axes by using a slope and travel time picking algorithm of local related homophase axes based on curvelet transformation, respectively tracking the local related homophase axes in a common shot point gather and a common survey point gather, picking corresponding travel time and slope data, and operating three-dimensional seismic data along a main survey line (x direction) and a cross survey line (y direction) respectively to obtain the following observation data space:
wherein ,TSR_obs Reflecting travel time for the observed seismic wave from the shot point S to the wave detection point R; p SX_obs and PSY_obs Observed seismic slopes in the x and y directions at the shot point, respectively; p RX_obs ,、P RY_obs Respectively, the observed seismic slopes in the x and y directions at the waypoints, and N is the number of shot point waypoint pairs, i.e., shot pair number, from which data is picked up.
Step three, establishing a discrete velocity model, forming a model space, and giving an initial discrete velocity model, which specifically comprises the following steps: determining an undulating observation surface according to elevation data of observation system data, and determining the sizes of the three-dimensional speed model in the x direction, the y direction and the z direction according to the range of a study area; determining the size of a discrete grid according to the complexity of an underground structure and the resolution requirement of a speed model to be established; discretizing the discrete velocity model, setting velocity values of all grid nodes, wherein the velocity value of any point in the discrete velocity model is represented by a cubic spline function of the velocity values of the grid nodes, the velocity value of each grid node is a model parameter, namely an unknown quantity, and the model space is as follows:
wherein ,Vj Assigning a value to the speed of each grid point to form an initial speed model, wherein M is the number of the grid nodes and the speed of the jth grid node is the speed value of the jth grid node;
determining the emergent directions of the rays at the shot point and the wave detection point according to the picked shot point slope data, namely obtaining initial values of ray parameters, and then obtaining ray paths, travel time and corresponding ray parameter data of the current speed model by solving a three-dimensional equation, wherein the method specifically comprises the following steps of:
using the picked shot slope data, calculating the emergent direction angle of the ray at the shot and the ray parameter vertical component at the shot by using the following formula:
wherein ,
the included angle between the shot point ray direction and the plane passing through the x axis is formed; theta S The included angle between the projection of the ray direction on the plane passing through the x axis and the z axis is shown; v 0S The velocity at the shot point; p SX_obs and PSY_obs Seismic slopes in the x and y directions at the picked shot point, respectively; p SZ0 Is the vertical slowness component at the shot point; sin for medical use -1 Representing an arcsine function; cos represents the cosine function, the initial ray parameter P at the shot point S0 Is composed of
P S0 ＝[P SX_obs ,P SY_obs ,P SZ0 ]
Using the picked slope data of the detection point, the exit direction angle of the ray at the detection point and the ray parameter vertical component at the detection point are calculated by the following formula:
wherein ,
the included angle between the ray direction of the wave detection point and the plane passing through the x axis is shown; theta R The included angle between the projection of the ray direction on the plane passing through the x axis and the z axis is shown; v 0R Is the velocity at the pickup point; p RX_obs and PRY_obs Seismic slopes in x and y directions at the picked-up wave detection points, respectively; p RZ0 Is the vertical slowness component at the detection point; sin for medical use -1 Indicates inversion or orthogenesisA chord function; cos represents the cosine function. At this time, the initial ray parameter P at the detection point R0 Is composed of
P R0 ＝[P RX_obs ,P RY_obs ,P RZ0 ]
Starting from the shot point S and the demodulator probe R respectively, solving the following three-dimensional equation by using the initial ray parameters of the shot point and the demodulator probe obtained in the steps S3 and S4,
wherein X ═ (X, Y, Z) is a spatial position vector; p (x) ═ P X (x),P Y (x),P Z (x)]The vector is composed of ray parameters or slowness components in three coordinate axis directions; v (x) is the velocity at x; t is the seismic travel time from the shot point or the demodulator probe; ^ represents the gradient operator.
Solving the differential equation to continuously obtain the path of the ray from the shot point S to the underground and the single-pass travel time t S And the path of the ray propagating from the point of detection R into the subsurface and the single travel time t R When t is S And t R The sum of which is equal to the observed two-way travel time T of the shot-check pair SR_obs Stopping calculation, and recording the ray end point position of the shot point S as x S And the ray end point position of the wave detection point R is recorded as x R Storing the ray length and the ray parameter of the ray of the shot point and the wave detection point in each corresponding grid unit;
and step five, calculating errors of the reflection points determined from the shot points and the wave detection points respectively and errors of observed reflection travel time and calculated reflection travel time, judging whether the current velocity model meets the requirements or not, and determining whether the inversion process is terminated or not.
2. The method of claim 1, wherein the step five comprises calculating errors of reflection points and errors of observed reflection travel times and calculated reflection travel times determined from the shot point and the geophone point, respectively, determining whether a current velocity model meets requirements, and determining whether an inversion process is terminated, and specifically comprises the steps of:
s1, if x S And x R When the single-pass travel time sum of the rays respectively along the shot point and the demodulator probe is consistent with the observed double-pass reflection travel time, the used discrete velocity model is correct, the distance between the tail end positions of the rays of the shot point and the demodulator probe of the shot-geophone pair is used for representing the error of the reflection point determined from the shot point and the demodulator probe, and the error of the reflection point and the error of the reflection travel time are respectively calculated by the following formulas:
wherein ,(XSi ,Y Si ,Z Si) and (XRi ,Y Ri ,Z Ri ) Coordinates of ray ends of a shot point and a demodulator probe of the ith shot-check pair are respectively; t is t Si and tRi Respectively the travel time from the shot point and the wave detection point of the ith shot-detection pair to the tail end of each ray; t is SRi_obs The reflected wave travel time of the ith shot-checking pair is observed;
s2, setting error limit, calculating err D and errT When the error is smaller than the given error limit, outputting the current speed model and stopping running; otherwise, executing the following steps;
s3, establishing a specific objective function and an inversion equation, solving to obtain a correction quantity of the current speed model, and updating the speed model;
s4, establishing the following objective function:
and is
D SRi ＝[(X Si -X Ri ) 2 +(Y Si -Y Ri ) 2 +(Z Si -Z Ri ) 2 ] 1/2
T SRi ＝t Si +t Ri -T SRi_obs
Wherein N is the number of shot detection pairs; t is t Si and tRi Respectively the travel time from the shot point and the wave detection point of the ith shot-detection pair to the tail end of each ray; t is SRi_obs The reflected wave of the ith shot-checking pair is subjected to double-pass travel; t is SRi The difference between the sum of two single-pass travel times of the ith shot-checking pair and the observed two-pass travel time; x is the number of Si ＝(X Si ,Y Si ,Z Si) and xRi ＝(X Ri ,Y Ri ,Z Ri ) Coordinate vectors of ray tail ends of shot points and wave detection points of the ith shot and survey pair respectively; d SRi Is the ray end point x Si And x Ri The distance between them; v is a model parameter vector consisting of mesh node velocities; λ is a weight coefficient;
s5, solving the minimum value solution of the nonlinear objective function by adopting a Gaussian-Newton method, performing Taylor series expansion on the objective function to obtain a first-order term, and obtaining a linear inversion equation as follows:
G·Δv＝Δd
wherein ,
B ij ＝λ(LS ij +LR ij )，i＝1,2,…,N,j＝1,2,…,M；
Δd＝[D SR1 D SR2 … D SRN T SR1 T SR2 … T SRN ] T ；
Δv＝[V 1 V 2 … V M ] T ；
D SRi ＝[(X Si -X Ri ) 2 +(Y Si -Y Ri ) 2 +(Z Si -Z Ri ) 2 ] 1/2 ,i＝1,2,…,N；
T SRt ＝t St +t Rt -T SRt_obs ,i＝1,2,…,N；
t Si and tRi Respectively the travel time from the shot point and the wave detection point of the ith shot-detection pair to the tail end of each ray; t is SRi_obs The reflected wave travel time of the ith shot-checking pair is observed; (X) Si ,Y Si ,Z Si) and (XRi ,Y Ri ,Z Ri ) Coordinates of ray tail ends of a shot point and a wave detection point of the ith shot and survey pair respectively; v i Is the velocity value of the ith discrete grid; v Si and VRi The speeds of the unit where the tail ends of the shot point and the wave detection point of the ith shot and wave detection pair are located respectively; p SXi ，P SYi and PSZi Ray parameters of the shot point ray tail end of the ith shot and survey pair along the x, y and z directions respectively; p RXi ，P RYi and PRZi Ray parameters of the ray tail end of the wave detection point of the ith shot-examination pair along the x direction, the y direction and the z direction respectively; l is Sij and LRij The lengths of ray segments of shot point rays and wave detection point rays of the ith shot and survey pair in the jth unit are respectively set; l is SXij ，L SYij and LSZij Projection lengths of ray segments of shot point rays of the ith shot and examine pair in the jth unit in the x direction, the y direction and the z direction are respectively set; l is RXij ，L RYij and LRZij Projection lengths of ray segments of wave detection point rays of the ith shot-examination pair in the jth unit in the directions of x, y and z are respectively set; n is the total number of the shot-test pairs; m is the total number of grid nodes of the discrete model;
s6, solving the linear inversion equation by using an LSQR algorithm to obtain a discrete velocity model updating quantity delta v, modifying the current discrete velocity model by using the model updating quantity to obtain a modified model, wherein the calculation formula is as follows:
v (k+1) ＝v (k) +Δv
wherein ,v(k) Is the current departureThe dispersion velocity model is a dispersion velocity model of the last iteration (the kth time), and delta v is the correction quantity of the dispersion velocity model obtained by the current iteration; v. of (k+1) And (4) taking the corrected discrete velocity model obtained by the iteration as a new initial model, and turning to the step four.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110132180.9A CN114839675B (en) | 2021-01-31 | 2021-01-31 | Method for establishing three-dimensional speed model |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110132180.9A CN114839675B (en) | 2021-01-31 | 2021-01-31 | Method for establishing three-dimensional speed model |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114839675A true CN114839675A (en) | 2022-08-02 |
| CN114839675B CN114839675B (en) | 2023-09-05 |
Family
ID=82560816
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110132180.9A Active CN114839675B (en) | 2021-01-31 | 2021-01-31 | Method for establishing three-dimensional speed model |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114839675B (en) |
Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20140200814A1 (en) * | 2013-01-11 | 2014-07-17 | Cgg Services Sa | Dip tomography for estimating depth velocity models by inverting pre-stack dip information present in migrated/un-migrated pre-/post-stack seismic data |
| CN105093279A (en) * | 2014-05-16 | 2015-11-25 | 中国石油化工股份有限公司 | Three-dimensional seismic primary wave Fresnel volume chromatography inversion method specific for Piedmont zones |
| CN105319589A (en) * | 2014-07-25 | 2016-02-10 | 中国石油化工股份有限公司 | Full-automatic three-dimensional tomography inversion method using a local event slope |
| EP3067718A1 (en) * | 2015-03-12 | 2016-09-14 | CGG Services SA | Boundary layer tomography method and device |
| CN107505651A (en) * | 2017-06-26 | 2017-12-22 | 中国海洋大学 | Seismic first break and back wave joint slope chromatography imaging method |
| WO2018021630A1 (en) * | 2016-07-29 | 2018-02-01 | 가톨릭대학교 산학협력단 | Method for manufacturing customized artificial tooth and customized artificial tooth |
| CN109387868A (en) * | 2018-09-28 | 2019-02-26 | 中国海洋石油集团有限公司 | A kind of three-dimensional chromatography imaging method based on seismic wave lineups slope information |
| CN109444956A (en) * | 2019-01-09 | 2019-03-08 | 中国海洋大学 | Three-dimensional fluctuating inspection surface earthquake slope chromatography imaging method |
| WO2019071504A1 (en) * | 2017-10-12 | 2019-04-18 | 南方科技大学 | Two-point ray tracing based seismic travel time tomography inversion method |
-
2021
- 2021-01-31 CN CN202110132180.9A patent/CN114839675B/en active Active
Patent Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20140200814A1 (en) * | 2013-01-11 | 2014-07-17 | Cgg Services Sa | Dip tomography for estimating depth velocity models by inverting pre-stack dip information present in migrated/un-migrated pre-/post-stack seismic data |
| CN105093279A (en) * | 2014-05-16 | 2015-11-25 | 中国石油化工股份有限公司 | Three-dimensional seismic primary wave Fresnel volume chromatography inversion method specific for Piedmont zones |
| CN105319589A (en) * | 2014-07-25 | 2016-02-10 | 中国石油化工股份有限公司 | Full-automatic three-dimensional tomography inversion method using a local event slope |
| EP3067718A1 (en) * | 2015-03-12 | 2016-09-14 | CGG Services SA | Boundary layer tomography method and device |
| WO2018021630A1 (en) * | 2016-07-29 | 2018-02-01 | 가톨릭대학교 산학협력단 | Method for manufacturing customized artificial tooth and customized artificial tooth |
| CN107505651A (en) * | 2017-06-26 | 2017-12-22 | 中国海洋大学 | Seismic first break and back wave joint slope chromatography imaging method |
| WO2019071504A1 (en) * | 2017-10-12 | 2019-04-18 | 南方科技大学 | Two-point ray tracing based seismic travel time tomography inversion method |
| CN109387868A (en) * | 2018-09-28 | 2019-02-26 | 中国海洋石油集团有限公司 | A kind of three-dimensional chromatography imaging method based on seismic wave lineups slope information |
| CN109444956A (en) * | 2019-01-09 | 2019-03-08 | 中国海洋大学 | Three-dimensional fluctuating inspection surface earthquake slope chromatography imaging method |
Non-Patent Citations (2)
| Title |
|---|
| 金昌昆等: "地震立体层析成像的实现方法及效果分析", 《CT理论与应用研究》 * |
| 金昌昆等: "地震立体层析成像的实现方法及效果分析", 《CT理论与应用研究》, vol. 23, no. 06, 30 November 2014 (2014-11-30), pages 939 - 950 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114839675B (en) | 2023-09-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US8760965B2 (en) | Time lapse marine seismic surveying employing interpolated multicomponent streamer pressure data | |
| US7523003B2 (en) | Time lapse marine seismic surveying | |
| US8792299B2 (en) | Method and device for processing seismic data | |
| US8406081B2 (en) | Seismic imaging systems and methods employing tomographic migration-velocity analysis using common angle image gathers | |
| US9869783B2 (en) | Structure tensor constrained tomographic velocity analysis | |
| MX2011006036A (en) | Using waveform inversion to determine properties of a subsurface medium. | |
| CN109444956B (en) | Three-dimensional undulation observation surface seismic slope tomography method | |
| EP1879052A2 (en) | Time lapse marine seismic surveying employing interpolated multicomponent streamer pressure data | |
| CN110187382B (en) | Traveling time inversion method for wave equation of reverse wave and reflected wave | |
| EA032186B1 (en) | Seismic adaptive focusing | |
| CN111665556B (en) | Stratum acoustic wave propagation velocity model construction method | |
| US9658354B2 (en) | Seismic imaging systems and methods employing correlation-based stacking | |
| CN113466933B (en) | Depth weighting-based seismic slope tomography method | |
| Zheng et al. | Double-beam stacking to infer seismic properties of fractured reservoirs | |
| CN114839675B (en) | Method for establishing three-dimensional speed model | |
| CN115598704A (en) | Method and device for generating amplitude-preserving angle gather based on least square reverse time migration and readable storage medium | |
| CN111665546B (en) | Acoustic parameter acquisition method for combustible ice detection | |
| AU2015201786A1 (en) | Methods and systems to separate wavefields using pressure wavefield data | |
| CN115061197A (en) | Two-dimensional sea surface ghost wave water body imaging measurement method, system, terminal and flow measurement equipment | |
| CN111665550A (en) | Underground medium density information inversion method | |
| CN111665549A (en) | Inversion method of stratum acoustic wave attenuation factor | |
| CN111665551B (en) | Acoustic parameter acquisition method for bridge substrate detection | |
| CN111665544A (en) | Acoustic parameter acquisition method for underground goaf detection | |
| CN111665548A (en) | Acoustic parameter acquisition method for seafloor detection | |
| CN111665543A (en) | Acoustic parameter acquisition method for subway danger detection |
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 |