CN114839675B - Method for establishing three-dimensional speed model - Google Patents
Method for establishing three-dimensional speed model Download PDFInfo
- Publication number
- CN114839675B CN114839675B CN202110132180.9A CN202110132180A CN114839675B CN 114839675 B CN114839675 B CN 114839675B CN 202110132180 A CN202110132180 A CN 202110132180A CN 114839675 B CN114839675 B CN 114839675B
- Authority
- CN
- China
- Prior art keywords
- point
- ray
- shot
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 43
- 238000001514 detection method Methods 0.000 claims abstract description 64
- 238000012937 correction Methods 0.000 claims abstract description 10
- 238000007781 pre-processing Methods 0.000 claims abstract description 4
- 238000004364 calculation method Methods 0.000 claims description 13
- 239000013598 vector Substances 0.000 claims description 12
- 230000008569 process Effects 0.000 claims description 7
- 230000001902 propagating effect Effects 0.000 claims description 6
- 230000005855 radiation Effects 0.000 claims description 6
- 230000009466 transformation Effects 0.000 claims description 4
- 238000011160 research Methods 0.000 claims description 3
- 238000012545 processing Methods 0.000 abstract description 2
- 238000003384 imaging method Methods 0.000 description 7
- 238000013508 migration Methods 0.000 description 6
- 230000005012 migration Effects 0.000 description 6
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000003325 tomography Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000002474 experimental method 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
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 the common shot point gather and the common detector point gather to form an observation data space; establishing a discrete speed model, forming a model space, and giving an initial discrete speed model; in the current speed model, respectively calculating the initial ray directions of the shot point and the wave detection point which are transmitted to the underground by using the picked slope data, and obtaining the ray paths, the travel time and the corresponding ray parameters of the shot point and the wave detection point respectively by solving a three-dimensional function equation; establishing an inversion objective function and an inversion equation, solving and obtaining correction quantity of the speed model, and updating the speed model; judging whether the current speed model meets the requirements, if so, outputting the current speed model, wherein the invention 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 pre-stack depth migration imaging of seismic survey data is greatly dependent on the accuracy of the subsurface velocity model used. Seismic travel time tomographic inversion is an effective technique for building subsurface low-frequency velocity models using seismic travel time data. Conventional reflection wave travel-time tomographic inversion needs to continuously track each phase axis and pick up travel-time data of each phase axis, and particularly needs to establish corresponding relations between the phase axes and the underground stratum interface in advance, which is often difficult to realize for complicated geological structure areas or low signal-to-noise ratio seismic data. In addition, in the reflection wave travel time tomographic inversion, the seismic wave travel time of the current velocity model and a kernel function matrix in an inversion equation need to be calculated through ray tracing. A ray path from the shot point to the detector point is determined, and the exit angle and the incident angle of the ray are known in addition to the velocity model, the shot point and the detector point position. The conventional reflection wave travel time tomographic inversion determines the exit angle and the incident angle through a plurality of experiments, or adopts a mode of firstly solving a wavefront surface and then determining rays by the wavefront surface, so that the calculated amount of ray tracing is large, and the travel time tomographic inversion has low efficiency.
Summary of the invention
The invention aims at solving at least one of the technical problems in the prior art, and provides a method for establishing a three-dimensional speed model by using seismic wave travel time and slope data.
The method for establishing the three-dimensional speed model provided by the embodiment of the invention comprises the following steps:
step one, preprocessing original three-dimensional seismic data;
step two, picking up the earthquake travel time and slope data of the shot point and the wave detection point to form an observation data space, which comprises the following steps: tracking the local related in-phase axis in the common shot point gather and the common detection point gather respectively by utilizing a slope and travel time picking algorithm of the local related in-phase axis based on curvelet transformation, picking corresponding travel time and slope data, and obtaining the following observation data space after operating the three-dimensional seismic data along a main line (x direction) and a cross line (y direction) respectively:
wherein ,TSR_obs Reflecting the travel time for the observed seismic waves from the shot point S to the detection point R; p (P) SX_obs and PSY_obs Seismic slopes along x and y directions at the shot point for observation, respectively; p (P) RX_obs ,、P RY_obs The observed seismic slopes in the x and y directions at the geophones are shown, respectively, with N being the number of pairs of the geophones, i.e., the number of pairs of the geophones, from which data was picked up.
Step three, a discrete speed model is established, a model space is formed, and an initial discrete speed model is given, and the method specifically comprises the following steps: according to the elevation data of the observation system data, determining an undulating observation surface, and then according to the range of a research area, determining the sizes of the three-dimensional speed model in the x, y and z directions; determining the size of a discrete grid according to the complexity of the underground structure and the resolution requirement of a speed model to be built; discretizing the discrete speed model, setting the speed value of each grid node, wherein the speed value of any point in the discrete speed model is represented by a cubic spline function of the speed value of the grid node, the speed value of each grid node is a model parameter, namely an unknown quantity, and the model space is:
wherein ,Vj The speed value of the jth grid node is given, M is the number of the grid nodes, and the speed of each grid point is assigned to form an initial speed model;
determining the emergent directions of rays at the shot point and at the detector point according to the picked shot point slope data, namely obtaining an initial value of the ray parameters, and then obtaining the ray path, 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:
calculating an exit direction angle of a ray at the shot point, and a ray parameter vertical component at the shot point, using the picked shot point slope data, with the following formula:
wherein ,the included angle between the shot point ray direction and the x-axis passing plane; θ S The included angle between the projection of the ray direction on the x-axis passing plane and the z-axis; v (V) 0S Is the velocity at the shot point; p (P) SX_obs and PSY_obs Seismic slopes along x and y directions at the picked shots, respectively; p (P) SZ0 Is the vertical slowness component at the shot point; sin (sin) -1 Representing an arcsine function; cos represents the cosine function, the initial ray parameters P at the shot point S0 Is that
P S0 ＝[P SX_obs ,P SY_obs ,P SZ0 ]
Using the pick-up dip slope data, the ray's exit direction angle at the dip, and the ray parameter vertical component at the dip, are calculated with the following formula:
wherein ,the included angle between the ray direction of the wave detection point and the x-axis passing plane is set; θ R The included angle between the projection of the ray direction on the x-axis passing plane and the z-axis; v (V) 0R Is the velocity at the detector point; p (P) RX_obs and PRY_obs Seismic slopes along x and y directions at the picked-up pickup points, respectively; p (P) RZ0 The vertical slowness component at the detection point; sin (sin) -1 Representing an arcsine function; cos represents a cosine function. At this time, the initial ray parameter P at the detector point R0 Is that
P R0 ＝[P RX_obs ,P RY_obs ,P RZ0 ]
Starting from the shot point S and the detector point R respectively, using the shot point and the detector point initial ray parameters obtained in the steps S3 and S4, solving the following three-dimensional equation by adopting a Runge-Kutta method,
wherein x= (X, Y, Z) is a spatial position vector; p (x) = [ P ] X (x),P Y (x),P Z (x)]Is a vector composed of ray parameters or slowness components in three coordinate axis directions; v (x) is the speed at x; t is the travel time of the seismic wave from the shot point or the wave detection point; and represents a gradient operator.
Solving the differential equation to continuously obtain the path of the ray propagating from the shot point S to the underground and the single-pass travel time t S And a ray path and a single travel time t propagating from the detector point R into the subsurface R When t S And t R The sum is equal to the observed double pass time T of the offset pair SR_obs When the calculation is stopped, the position of the ray end point of the shot point S is marked as x S The position of the ray end point of the detector R is marked as x R Storing the rays of the shots and detectors in respective grid cellsRay length and ray parameters of (a);
and fifthly, calculating reflection point errors and errors of observation reflection travel time and calculation reflection travel time respectively determined from the shot point and the detection point, judging whether the current speed model meets the requirement, and determining whether the inversion process is terminated.
Further, the fifth step of calculating the reflection point error and the error between the observation reflection travel time and the calculation reflection travel time determined from the shot point and the detection point respectively, judging whether the current speed model meets the requirement, and determining whether the inversion process is terminated, specifically comprising the following steps:
s1, if x S And x R When the single-pass travel time sum of rays along the shot point and the detector point coincides with the observed double-pass reflection travel time, the discrete speed model is correct, the reflection point error determined from the shot point and the detector point is represented by the distance between the shot point and the detector point of the shot-detection pair, and the reflection point error and the reflection travel time error are calculated by the following formulas:
wherein ,(XSi ，Y Si ，Z Si) and (XRi ，Y Ri ，Z Ri ) Coordinates of ray ends of the shot point and the geophone from the ith offset pair, respectively; t is t Si and tRi The travel time from the shot point and the detection point of the ith offset pair to the tail end of each ray; t (T) SRi obs Reflected wave travel time for the i-th offset pair observed;
s2, giving an error limit, when the err is calculated D and errT Outputting the current speed model when the error time is smaller than the given error time limit, and stopping operation; otherwise, the following steps are executed;
s3, establishing a specific objective function and an inversion equation, solving to obtain correction quantity of the current speed model, and updating the speed model;
s4, establishing the following objective function:
order of (A)
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 offset pairs; t is t Si and tRi The travel time from the shot point and the detection point of the ith offset pair to the tail end of each ray; t (T) SRi_obs Double-pass time-shifting of the reflected wave for the ith offset pair of observation; t (T) SRi The difference between the sum of the two single-trip times of the ith offset pair and the observed double-trip time; x is x Si ＝(X Si ，Y Si ，Z Si) and xRi ＝(X Ri ，Y Ri ，Z Ri ) Coordinate vectors of ray ends of the shot point and the detector point of the ith offset pair respectively; d (D) SRi For ray end point x Si And x Ri A distance therebetween; v is a model parameter vector consisting of grid node velocities; lambda is a weight coefficient;
s5, solving a minimum solution of the nonlinear objective function by adopting a Gaussian-Newton method, and performing Taylor series expansion on the objective function to obtain a first-order term to obtain a linear inversion equation as follows:
G·Δv＝Δd
wherein ,
wherein i=1, 2,..n, j=1, 2,..m;
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 The travel time from the shot point and the detection point of the ith offset pair to the tail end of each ray; t (T) SRi_obs Reflected wave travel time for the i-th offset pair observed; (X) Si ，Y Si ，Z Si) and (XRi ，Y Ri ，Z Ri ) Coordinates of ray ends of a shot point and a detector point of the ith offset pair respectively; v (V) i Is the speed value of the ith discrete grid; v (V) Si and VRi The speeds of the units where the shot point of the ith offset pair and the tail end of the wave detection point are respectively located; p (P) SXi ，P SYi and PSZi The parameters of the ray end of the shot point of the ith offset pair along the x, y and z directions are respectively; p (P) RXi ，P RYi and PRZi The radiation parameters of the radiation end of the detection point of the ith offset pair along the x, y and z directions are respectively; l (L) Sij and LRij The lengths of the shot point rays and the wave detection point rays of the ith offset pair in the j unit are respectively; l (L) SXij ，L SYij and LSZij The projection lengths of shot rays of the ith offset pair in the x, y and z directions of the ray segments of the jth unit are respectively; l (L) RXij ，L RYij and LRZij The wave-point rays of the ith offset pair are respectively atProjection lengths of the ray segments of the jth cell in x, y and z directions; n is the total number of gun-receiver 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 speed model update quantity Deltav, and modifying the current discrete speed model by using the model update quantity to obtain a modified model, wherein the calculation formula is as follows:
v (k+1) ＝v (k) +Δv
wherein ,v(k) The method is a current discrete speed model, namely a discrete speed model of the last (kth) iteration, and Deltav is a discrete speed model correction amount obtained by inversion of the iteration; v (k+1) The modified discrete velocity model obtained in the iteration is used as a new initial model, and the step four is carried out.
The beneficial effects of the invention are as follows: by adopting the method for establishing the three-dimensional speed model by using the travel time and the slope data of the seismic reflection waves, the slope data of the phase axis of the local correlation reflection waves is used for calculating the parameters of the propagation direction from the shot point and the wave detection point to the underground and the initial ray, so that the efficiency of calculating the ray path and the travel time is improved; the used local correlation reflected wave phase axis is easy to automatically track, the travel time and slope data are easy to automatically pick up, and the method has strong adaptability to complex underground structures and low signal-to-noise ratio seismic data. The embodiment of the invention is easy for computer automation realization, and provides an effective three-dimensional speed modeling method for seismic imaging of a region with a complex structure.
Drawings
The invention is further described below with reference to the drawings and examples;
FIG. 1 is a flow chart of a method for establishing a three-dimensional velocity model according to an embodiment of the present invention;
FIG. 2 is a slice of the three-dimensional velocity model created using travel time and slope data for the embodiment of FIG. 1;
FIG. 3 is a slice of a three-dimensional velocity model constructed by conventional tomographic inversion methods using travel time data;
FIG. 4 is a graph of a common imaging point gather after shifting using a three-dimensional velocity model constructed in accordance with an embodiment of the present invention;
FIG. 5 is a plot of a trace of co-imaging points after displacement using a three-dimensional velocity model established by conventional travel-time tomography.
Detailed Description
The present invention will be described in further detail with reference to the accompanying drawings, in order to make the objects and technical solutions of the present invention more apparent.
In the seismic data, the common shot gather and the common detector gather are easy to track the local correlation phase axis, and the change rate of the travel time of the local correlation phase axis along with the offset distance is called a 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 and the detector are known, the exit angle of the reflected ray at the shot and the angle of incidence at the detector can be calculated using the slope parameters. In addition, the local coherence phase axis is easy to track, the travel time and the slope of the local coherence phase axis are easy to pick up, and the local coherence phase axis has strong adaptability to low signal-to-noise ratio seismic data. If a smooth speed model is adopted, rays are constrained through seismic slopes, and only locally related in-phase axis travel time and slope data are needed, so that the corresponding relation between each in-phase axis and stratum interface is not needed to be established in advance, and the automation degree and the practicability of the reflection wave travel time tomographic inversion are improved. Therefore, if the seismic slope data is introduced into the reflection wave travel-time tomographic inversion, the ray emergence angle and the incident 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 pre-stack depth migration velocity analysis and waveform inversion.
According to the embodiment of the invention, the earthquake slope data are introduced into the three-dimensional earthquake reflection travel time tomographic inversion, the ray emergence angle and the incident angle are determined by using the earthquake slope, rays which respectively propagate from the shot point and the detection point to the underground are encountered at the reflection point according to the condition that the velocity model is accurate and the travel time sum of the shot point and the detection point is equal to the actually measured two-way reflection travel time, a method for establishing a three-dimensional velocity model by using the earthquake reflection travel time and the slope data is provided, and a technology for establishing the underground three-dimensional velocity model is provided for the migration imaging of earthquake exploration or low signal-to-noise ratio earthquake data in a region with a complex structure.
Referring to fig. 1, a method for establishing a three-dimensional velocity model according to an embodiment of the present invention includes the following steps:
step one, preprocessing original three-dimensional seismic data;
step two, picking up the earthquake travel time and slope data of the shot point and the wave detection point to form an observation data space, which comprises the following steps: tracking the local related in-phase axis in the common shot point gather and the common detection point gather respectively by utilizing a slope and travel time picking algorithm of the local related in-phase axis based on curvelet transformation, picking corresponding travel time and slope data, and obtaining the following observation data space after operating the three-dimensional seismic data along a main line (x direction) and a cross line (y direction) respectively:
wherein ,TSR_obs Reflecting the travel time for the observed seismic waves from the shot point S to the detection point R; p (P) SX_obs and PSY_obs Seismic slopes along x and y directions at the shot point for observation, respectively; p (P) RX_obs ,、P RY_obs The observed seismic slopes in the x and y directions at the geophones are shown, respectively, with N being the number of pairs of the geophones, i.e., the number of pairs of the geophones, from which data was picked up.
Step three, a discrete speed model is established, a model space is formed, and an initial discrete speed model is given, and the method specifically comprises the following steps: according to the elevation data of the observation system data, determining an undulating observation surface, and then according to the range of a research area, determining the sizes of the three-dimensional speed model in the x, y and z directions; determining the size of a discrete grid according to the complexity of the underground structure and the resolution requirement of a speed model to be built; discretizing the discrete speed model, setting the speed value of each grid node, wherein the speed value of any point in the discrete speed model is represented by a cubic spline function of the speed value of the grid node, the speed value of each grid node is a model parameter, namely an unknown quantity, and the model space is:
wherein ,Vj The speed value of the jth grid node is given, M is the number of the grid nodes, and the speed of each grid point is assigned to form an initial speed model;
determining the emergent directions of rays at the shot point and at the detector point according to the picked shot point slope data, namely obtaining an initial value of the ray parameters, and then obtaining the ray path, 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:
calculating an exit direction angle of a ray at the shot point, and a ray parameter vertical component at the shot point, using the picked shot point slope data, with the following formula:
wherein ,the included angle between the shot point ray direction and the x-axis passing plane; θ S The included angle between the projection of the ray direction on the x-axis passing plane and the z-axis; v (V) 0S Is the velocity at the shot point; p (P) SX_obs and PSY_obs Seismic slopes along x and y directions at the picked shots, respectively; p (P) SZ0 Is the vertical slowness component at the shot point; sin (sin) -1 Representing an arcsine function; cos represents the cosine function, the initial ray parameters P at the shot point S0 Is that
P S0 ＝[P SX_obs ,P SY_obs ,P SZ0 ]
Using the pick-up dip slope data, the ray's exit direction angle at the dip, and the ray parameter vertical component at the dip, are calculated with the following formula:
wherein ,the included angle between the ray direction of the wave detection point and the x-axis passing plane is set; θ R The included angle between the projection of the ray direction on the x-axis passing plane and the z-axis; v (V) 0R Is the velocity at the detector point; p (P) RX_obs and PRY_obs Seismic slopes along x and y directions at the picked-up pickup points, respectively; p (P) RZ0 The vertical slowness component at the detection point; sin (sin) -1 Representing an arcsine function; cos represents a cosine function. At this time, the initial ray parameter P at the detector point R0 Is that
P R0 ＝[P RX_obs ,P RY_obs ,P RZ0 ]
Starting from the shot point S and the detector point R respectively, using the shot point and the detector point initial ray parameters obtained in the steps S3 and S4, solving the following three-dimensional equation by adopting a Runge-Kutta method,
wherein x= (X, Y, Z) is a spatial position vector; p (x)＝[P X (x),P Y (x),P Z (x)]Is a vector composed of ray parameters or slowness components in three coordinate axis directions; v (x) is the speed at x; t is the travel time of the seismic wave from the shot point or the wave detection point;representing the gradient operator.
Solving the differential equation to continuously obtain the path of the ray propagating from the shot point S to the underground and the single-pass travel time t S And a ray path and a single travel time t propagating from the detector point R into the subsurface R When t S And t R The sum is equal to the observed double pass time T of the offset pair SR_obs When the calculation is stopped, the position of the ray end point of the shot point S is marked as x S The position of the ray end point of the detector R is marked as x R Storing the ray length and the ray parameters of the rays of the shot point and the detector point in the corresponding grid cells;
and fifthly, calculating reflection point errors and errors of observation reflection travel time and calculation reflection travel time respectively determined from the shot point and the detection point, judging whether the current speed model meets the requirement, and determining whether the inversion process is terminated.
Further, the fifth step of calculating the reflection point error and the error between the observation reflection travel time and the calculation reflection travel time determined from the shot point and the detection point respectively, judging whether the current speed model meets the requirement, and determining whether the inversion process is terminated, specifically comprising the following steps:
s1, if x S And x R When the single-pass travel time sum of rays along the shot point and the detector point coincides with the observed double-pass reflection travel time, the discrete speed model is correct, the reflection point error determined from the shot point and the detector point is represented by the distance between the shot point and the detector point of the shot-detection pair, and the reflection point error and the reflection travel time error are calculated by the following formulas:
wherein ,(XSi ，Y Si ，Z Si) and (XRi ，Y Ri ，Z Ri ) Coordinates of ray ends of the shot point and the geophone from the ith offset pair, respectively; t is t Si and tRi The travel time from the shot point and the detection point of the ith offset pair to the tail end of each ray; t (T) SRi_obs Reflected wave travel time for the i-th offset pair observed;
s2, giving an error limit, when the err is calculated D and errT Outputting the current speed model when the error time is smaller than the given error time limit, and stopping operation; otherwise, the following steps are executed;
s3, establishing a specific objective function and an inversion equation, solving to obtain correction quantity of the current speed model, and updating the speed model;
s4, establishing the following objective function:
and is also provided with
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 offset pairs; t is t Si and tRi The travel time from the shot point and the detection point of the ith offset pair to the tail end of each ray; t (T) SRi_obs Double-pass time-shifting of the reflected wave for the ith offset pair of observation; t (T) SRi The difference between the sum of the two single-trip times of the ith offset pair and the observed double-trip time; x is x Si ＝(X Si ，Y Si ，Z Si) and xRi ＝(X Ri ，Y Ri ，Z Ri ) Respectively the ith cannonCoordinate vectors of ray ends of the shot point and the detector point to be detected; d (D) SRi For ray end point x Si And x Ri A distance therebetween; v is a model parameter vector consisting of grid node velocities; lambda is a weight coefficient;
s5, solving a minimum solution of the nonlinear objective function by adopting a Gaussian-Newton method, and performing Taylor series expansion on the objective function to obtain a first-order term to obtain a linear inversion equation as follows:
G·Δv＝Δd
wherein ,
wherein i=1, 2,..n, j=1, 2,..m;
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 The travel time from the shot point and the detection point of the ith offset pair to the tail end of each ray; t (T) SRi_obs Reflected wave travel time for the i-th offset pair observed; (X) Si ，Y Si ，Z Si) and (XRi ，Y Ri ，Z Ri ) Cannon of ith offset pair respectivelyCoordinates of the ray ends of the point and the detector; v (V) i Is the speed value of the ith discrete grid; v (V) Si and VRi The speeds of the units where the shot point of the ith offset pair and the tail end of the wave detection point are respectively located; p (P) SXi ，P SYi and PSZi The parameters of the ray end of the shot point of the ith offset pair along the x, y and z directions are respectively; p (P) RXi ，P RYi and PRZi The radiation parameters of the radiation end of the detection point of the ith offset pair along the x, y and z directions are respectively; l (L) Sij and LRij The lengths of the shot point rays and the wave detection point rays of the ith offset pair in the j unit are respectively; l (L) SXij ，L sYij and LSZij The projection lengths of shot rays of the ith offset pair in the x, y and z directions of the ray segments of the jth unit are respectively; l (L) RXij ，L RYij and LRZij The projection lengths of the ray segments of the ray of the wave-detecting point of the ith offset pair in the x, y and z directions of the jth unit are respectively; n is the total number of gun-receiver 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 speed model update quantity Deltav, and modifying the current discrete speed model by using the model update quantity to obtain a modified model, wherein the calculation formula is as follows:
v (k+1) ＝v (k) +Δv
wherein ,v(k) The method is a current discrete speed model, namely a discrete speed model of the last (kth) iteration, and Deltav is a discrete speed model correction amount obtained by inversion of the iteration; v (k+1) The modified discrete velocity model obtained in the iteration is used as a new initial model, and the step four is carried out.
As a specific example of the embodiment of the invention, the method of the embodiment of the invention is tested by three-dimensional oil-gas seismic exploration data in a region of southwest of China and compared with the result of reflection tomography by using reflection travel time data only.
Step 100, the length of the work area is 16km, and the width is 8km. And uniformly arranging 25 gun lines and 59 wave detection lines, wherein the distance between the gun lines is 300m, the distance between the wave detection lines is 150m, the distance between the gun lines is 100m and the distance between the wave detection points is 50m. When each gun is excited, 8 detection lines are received, and each detection line is received by 240 detection points. The original seismic data is preprocessed by noise suppression, amplitude compensation, deconvolution, static correction and the like.
Step 110, picking up travel time and slope data in a common shot point domain and a common detector point domain along the directions of a detection line and a shot line respectively by using a local correlation phase axis slope and travel time picking-up method based on curvelet transformation, obtaining reflected travel time and slope data of 236152 shot-detection pairs, and forming a data space.
Step 120, a speed model with the size of 16km×8km×4km is built, and the three-dimensional speed model is discretized by using a grid with the size of 250m×250m, so as to form a model space. The model surface builds up relief features with corresponding elevation data. The initial speed of the model is set to be linearly increased along with the depth, the highest speed 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 130, obtaining the directions of the seismic rays which respectively propagate from the shot point and the wave point to the underground by using the initial speed and the slope data of the shot point and the wave point which are picked up, and giving initial conditions of ray emergence.
And 140, solving a three-dimensional equation by adopting a range-Kutta method according to the initial velocity model and the ray initial conditions to obtain the ray end positions and the travel time of each gun-receiver pair and the corresponding ray paths and ray parameters of each grid unit.
And 150, calculating a reflection point error by using the distances between the shot point and the ray end positions of the detection points of each shot-detection pair, and calculating a reflection travel time error by using the travel time of the shot point and the ray end of the detection point of each shot-detection pair and the actual reflection wave travel time.
And 160, comparing the calculated reflection point error and the reflection time error with error limits set in advance respectively, and judging whether the current speed model meets the requirement. If the calculated reflection point error and the reflection travel time error are smaller than the set error limit, outputting a current speed model, and stopping calculation; otherwise, the following steps are performed.
Step 170, establishing inversion equations using the observed reflection travel time and slope data, the ray end positions of each of the gun-receiver pairs calculated in step S5, the travel time, and the corresponding ray paths and ray parameters of each of the grid cells. And solving an inversion equation to obtain the correction quantity of the current speed model.
And 180, correcting the current model by using the obtained model correction amount to obtain an updated model. The updated model is taken as a new initial model and the process goes to step S4.
FIG. 2 is a vertical slice of the subsurface velocity model along the middle line after 30 iterations of the invention, and FIG. 3 is a vertical slice of the subsurface velocity model along the middle line built using only reflected travel time data. Comparing these two models, the velocity model established by using the travel time and slope data in the present invention is greatly different from the conventional velocity model established by using only the travel time data in the lateral direction and the vertical direction.
In order to further illustrate the beneficial effects of the invention, the pre-stack depth migration imaging is performed on seismic data by utilizing the velocity model established by using travel time and slope data and the conventional velocity model established by using travel time data only. FIG. 4 shows the migration of the velocity model established by the present invention to obtain a common-image-point gather along the intermediate line, and FIG. 5 shows the migration of the velocity model established by conventional travel-time tomographic inversion to obtain a common-image-point gather along the intermediate line. The phenomenon of pull-down and upwarp of the same phase shaft in fig. 5 is obvious, the same phase shaft basically presents a non-horizontal form, and the same phase shaft in fig. 4 is relatively horizontal and smooth. In the common imaging point trace set, the accuracy of the velocity model used at the time of offset is often judged by the level of the in-phase axis. The more horizontal the same phase axis is, 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, indicating a good application effect of the present invention.
The foregoing description of the preferred embodiments of the invention is not intended to be limiting, but rather is to cover all modifications, equivalents, and alternatives falling within the spirit and principles of the invention.
Claims (2)
1. A method for creating a three-dimensional velocity model using seismic data, comprising the steps of:
step one, preprocessing original three-dimensional seismic data;
step two, picking up the earthquake travel time and slope data of the shot point and the wave detection point to form an observation data space, which comprises the following steps: tracking the local related in-phase axis in the common shot point gather and the common detection point gather respectively by utilizing a slope and travel time picking algorithm of the local related in-phase axis based on curvelet transformation, picking corresponding travel time and slope data, and operating the three-dimensional seismic data along a main line and a cross line, namely an x direction and a y direction respectively to obtain the following observation data space:
wherein ,TSR_obs Reflecting the travel time for the observed seismic waves from the shot point S to the detection point R; p (P) SX_obs and PSY_obs Seismic slopes along x and y directions at shot taken respectively; p (P) RX_obs 、P RY_obs The seismic slopes of the picked-up positions along the x and y directions are respectively shown, and N is the number of shot detection point pairs for picking up data, namely the number of shot detection pairs;
step three, a discrete speed model is established, a model space is formed, and an initial discrete speed model is given, and the method specifically comprises the following steps: according to the elevation data of the observation system data, determining an undulating observation surface, and then according to the range of a research area, determining the sizes of the three-dimensional speed model in the x, y and z directions; determining the size of a discrete grid according to the complexity of the underground structure and the resolution requirement of a speed model to be built; discretizing the discrete speed model, setting the speed value of each grid node, wherein the speed value of any point in the discrete speed model is represented by a cubic spline function of the speed value of the grid node, the speed value of each grid node is a model parameter, namely an unknown quantity, and the model space is:
wherein ,Vj The speed value of the jth grid node is given, M is the number of the grid nodes, and the speed of each grid point is assigned to form an initial speed model;
determining the emergent directions of rays at the shot point and at the detector point according to the picked shot point slope data, namely obtaining an initial value of the ray parameters, and then obtaining the ray path, 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:
calculating an exit direction angle of a ray at the shot point, and a ray parameter vertical component at the shot point, using the picked shot point slope data, with the following formula:
wherein ,the included angle between the shot point ray direction and the x-axis passing plane; θ S The included angle between the projection of the ray direction on the x-axis passing plane and the z-axis; v (V) 0S Is the velocity at the shot point; p (P) SX_obs and PSY_obs Seismic slopes along x and y directions at the picked shots, respectively; p (P) SZ0 Is the vertical slowness component at the shot point; sin (sin) -1 Representing an arcsine function; cos represents the cosine function, the initial ray parameters P at the shot point S0 Is that
P S0 ＝[P SX_obs ,P SY_obs ,P SZ0 ]
Using the pick-up dip slope data, the ray's exit direction angle at the dip, and the ray parameter vertical component at the dip, are calculated with the following formula:
wherein ,the included angle between the ray direction of the wave detection point and the x-axis passing plane is set; θ R The included angle between the projection of the ray direction on the x-axis passing plane and the z-axis; v (V) 0R Is the velocity at the detector point; p (P) RX_obs and PRY_obs Seismic slopes along x and y directions at the picked-up pickup points, respectively; p (P) RZ0 The vertical slowness component at the detection point; sin (sin) -1 Representing an arcsine function; cos represents a cosine function; at this time, the initial ray parameter P at the detector point R0 Is that
P R0 ＝[P RX_obs ,P RY_obs ,P RZ0 ]
Starting from the shot point S and the detector point R respectively, using the initial ray parameters of the shot point and the detector point obtained in the third step and the fourth step, adopting a Runge-Kutta method to solve the following three-dimensional equation,
wherein x= (X, Y, Z) is a spatial position vector; p (x) = [ P ] X (x),P Y (x),P Z (x)]Is a vector composed of ray parameters or slowness components in three coordinate axis directions; v (x) is the speed at x; t is the travel time of the seismic wave from the shot point or the wave detection point;representing a gradient operator;
solving the differential equation to continuously obtain the path of the ray propagating from the shot point S to the underground and the single-pass travel time t S And a ray path and a single travel time t propagating from the detector point R into the subsurface R When t S And t R The sum is equal to the observed double pass time T of the offset pair SR_obs When the calculation is stopped, the position of the ray end point of the shot point S is marked as x S The position of the ray end point of the detector R is marked as x R Storing the ray length and the ray parameters of the rays of the shot point and the detector point in the corresponding grid cells;
and fifthly, calculating reflection point errors and errors of observation reflection travel time and calculation reflection travel time respectively determined from the shot point and the detection point, judging whether the current speed model meets the requirement, and determining whether the inversion process is terminated.
2. The method for building three-dimensional velocity model according to claim 1, wherein the fifth step is to calculate the reflection point errors and the errors between the observation of the reflection travel time and the calculation of the reflection travel time determined from the shot and the detector, respectively, to determine whether the current velocity model meets the requirement, and to determine whether the inversion process is terminated, and the method specifically comprises the following steps:
s1, if x S And x R When the single-pass travel time sum of rays along the shot point and the detector point coincides with the observed double-pass reflection travel time, the discrete speed model is correct, the reflection point error determined from the shot point and the detector point is represented by the distance between the shot point and the detector point of the shot-detection pair, and the reflection point error and the reflection travel time error are calculated by the following formulas:
wherein ,(XSi ,Y Si ,Z Si) and (XRi ,Y Ri ,Z Ri ) Coordinates of ray ends of the shot point and the geophone from the ith offset pair, respectively; t is t Si and tRi The travel time from the shot point and the detection point of the ith offset pair to the tail end of each ray; t (T) SRi_obs Reflected wave travel time for the i-th offset pair observed;
s2, giving an error limit, when the err is calculated D and errT Outputting the current speed model when the error time is smaller than the given error time limit, and stopping operation; otherwise, the following steps are executed;
s3, establishing a specific objective function and an inversion equation, solving to obtain correction quantity of the current speed model, and updating the speed model;
s4, establishing the following objective function:
and is also provided with
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 offset pairs; t is t Si and tRi The travel time from the shot point and the detection point of the ith offset pair to the tail end of each ray; t (T) SRi_obs Double-pass time-shifting of the reflected wave for the ith offset pair of observation; t (T) SRi The difference between the sum of the two single-trip times of the ith offset pair and the observed double-trip time; x is x Si ＝(X Si ,Y Si ,Z Si) and xRi ＝(X Ri ,Y Ri ,Z Ri ) Coordinate vectors of ray ends of the shot point and the detector point of the ith offset pair respectively; d (D) SRi For ray end point x Si And x Ri A distance therebetween; v is the model space; lambda is a weight coefficient;
s5, solving a minimum solution of the nonlinear objective function by adopting a Gaussian-Newton method, and performing Taylor series expansion on the objective function to obtain a first-order term to obtain a linear inversion equation as follows:
G·Δv＝Δd
wherein ,
where i=1, 2, …, N, j=1, 2, …, M;
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 Ti ) 2 +(X 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 The travel time from the shot point and the detection point of the ith offset pair to the tail end of each ray; t (T) SRi_obs Reflected wave travel time for the i-th offset pair observed; (X) Si ,Y Si ,Z Si) and (XRi ,Y Ri ,Z Ri ) Shot and detection for the ith offset pair respectivelyCoordinates of the ray end of the point; v (V) i Is the speed value of the ith discrete grid; v (V) Si and VRi The speeds of the units where the shot point of the ith offset pair and the tail end of the wave detection point are respectively located; p (P) SXi ，P SYi and PSZi The parameters of the ray end of the shot point of the ith offset pair along the x, y and z directions are respectively; p (P) RXi ，P RYi and PRZi The radiation parameters of the radiation end of the detection point of the ith offset pair along the x, y and z directions are respectively; l (L) Sij and LRij The lengths of the shot point rays and the wave detection point rays of the ith offset pair in the j unit are respectively; l (L) SXij ，L SYij and LSZij The projection lengths of shot rays of the ith offset pair in the x, y and z directions of the ray segments of the jth unit are respectively; l (L) RXij ，L RYij and LRZij The projection lengths of the ray segments of the ray of the wave-detecting point of the ith offset pair in the x, y and z directions of the jth unit are respectively; n is the total number of offset 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 speed model update quantity Deltav, and modifying the current discrete speed model by using the model update quantity to obtain a modified model, wherein the calculation formula is as follows:
v (k+1) ＝v (k) +Δv
wherein ,v(k) The method is a current discrete speed model, namely a kth iteration discrete speed model, and Deltav is a discrete speed model correction amount obtained by the iteration inversion; v (k+1) The modified discrete velocity model obtained in the iteration is used as a new initial model, and the step four is carried out.
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 CN114839675A (en) | 2022-08-02 |
| CN114839675B true 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 (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| 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 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CA2838713A1 (en) * | 2013-01-11 | 2014-07-11 | 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 |
-
2021
- 2021-01-31 CN CN202110132180.9A patent/CN114839675B/en active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| 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 (1)
| Title |
|---|
| 地震立体层析成像的实现方法及效果分析;金昌昆等;《CT理论与应用研究》;20141130;第23卷(第06期);第939-950页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114839675A (en) | 2022-08-02 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| RU2693495C1 (en) | Complete wave field inversion with quality factor compensation | |
| US8406081B2 (en) | Seismic imaging systems and methods employing tomographic migration-velocity analysis using common angle image gathers | |
| CN105277978B (en) | A kind of method and device for determining near-surface velocity model | |
| US9869783B2 (en) | Structure tensor constrained tomographic velocity analysis | |
| CN109444956B (en) | Three-dimensional undulation observation surface seismic slope tomography method | |
| CA2683618C (en) | Inverse-vector method for smoothing dips and azimuths | |
| MX2011006036A (en) | Using waveform inversion to determine properties of a subsurface medium. | |
| CN105319589B (en) | A kind of fully automatic stereo chromatography conversion method using local lineups slope | |
| US8731838B2 (en) | Fresnel zone fat ray tomography | |
| CN111856577B (en) | Method for reducing calculation amount of reverse-time migration earth surface offset gather | |
| CN105137479B (en) | A kind of computational methods and device of bin degree of covering | |
| CN111665556B (en) | Stratum acoustic wave propagation velocity model construction method | |
| CN113466933B (en) | Depth weighting-based seismic slope tomography method | |
| CN114839675B (en) | Method for establishing three-dimensional speed model | |
| Popov et al. | Reverse time migration with Gaussian beams and velocity analysis applications | |
| 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 | |
| CN115061197A (en) | Two-dimensional sea surface ghost wave water body imaging measurement method, system, terminal and flow measurement equipment | |
| CN109655888B (en) | Quantitative selection method and system for smooth floating reference surface in seismic data processing | |
| 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 | |
| CN110873894A (en) | Shot record obtaining method and system based on Gaussian beam anti-migration | |
| CN111665544A (en) | Acoustic parameter acquisition method for underground goaf detection | |
| Tavakoli F et al. | Adjoint stereotomography |
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 |