CN114428292B - Method for constructing near-surface velocity model and storage medium - Google Patents
Method for constructing near-surface velocity model and storage medium Download PDFInfo
- Publication number
- CN114428292B CN114428292B CN202011004206.3A CN202011004206A CN114428292B CN 114428292 B CN114428292 B CN 114428292B CN 202011004206 A CN202011004206 A CN 202011004206A CN 114428292 B CN114428292 B CN 114428292B
- Authority
- CN
- China
- Prior art keywords
- ray
- point
- angle
- determining
- ray path
- 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 58
- 238000003860 storage Methods 0.000 title claims abstract description 8
- 238000004364 calculation method Methods 0.000 claims abstract description 25
- 238000001514 detection method Methods 0.000 claims abstract description 20
- 238000004590 computer program Methods 0.000 claims description 4
- 230000003094 perturbing effect Effects 0.000 claims description 3
- 238000010276 construction Methods 0.000 claims description 2
- 238000004587 chromatography analysis Methods 0.000 abstract description 8
- 230000014509 gene expression Effects 0.000 description 20
- 238000007796 conventional method Methods 0.000 description 3
- 238000003384 imaging method Methods 0.000 description 3
- 238000012937 correction Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 230000003068 static effect Effects 0.000 description 2
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000007918 pathogenicity Effects 0.000 description 1
- 230000001575 pathological effect Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 238000003325 tomography Methods 0.000 description 1
- 238000012876 topography Methods 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. analysis, for interpretation, for correction
- 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. analysis, for interpretation, for correction
- G01V1/282—Application of seismic models, synthetic seismograms
-
- 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
-
- 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/66—Subsurface modeling
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Abstract
The application provides a method for constructing a near-surface velocity model and a storage medium. The method comprises the following steps: s100: determining the azimuth angle and incidence inclination angle of the rays; s200: determining the position of each selected point on the ray path; s300: if each point selected on the ray path meets the preset angle condition, the ray does not meet the preset calculation stopping criterion, S400 is executed, otherwise S500 is executed; s400: if the distance from each point selected on the ray path to the detection point is smaller than or equal to the preset distance threshold, executing S600, otherwise, executing S500; s500: selecting a new set of incidence dip angle and azimuth angle, and returning to execute S200; s600: determining plane wave solution, sphere wave solution, gaussian beam, travel time chromatography kernel function and near-surface velocity updating quantity; s700: obtaining an updated speed model; s800: and determining a difference value between the measured travel time and the calculated travel time, wherein the difference value does not meet a preset time condition, and executing S100, otherwise, taking the updated speed model as a final near-surface speed model.
Description
Technical Field
The invention relates to the technical field of geophysical exploration, in particular to a method for constructing a near-surface velocity model and a storage medium.
Background
With the continuous deep development of seismic exploration, the surface topography of the complex mountain with severe change of near-surface speed seriously affects the quality of static correction and imaging. To improve the static correction effect and the imaging accuracy of complex subsurface structures, the accuracy of the near-surface velocity model is critical to solving the above problems. The method for modeling the surface velocity is that first arrival travel time information recorded by a shot gather is picked up firstly, an initial gradient velocity model is established, ray paths are calculated by ray tracing and first arrival time based on the model, a chromatographic equation is established according to the relation between residual errors between the picked first arrival and the calculated first arrival and the ray paths, the velocity update quantity is solved, and a final near-surface velocity model is obtained through multiple iterations, so that a foundation is laid for subsequent seismic processing and imaging.
The method of near-surface velocity modeling includes the following: the first is ray chromatography, which has high calculation efficiency and low dependence on an initial model, but ray theory is limited by problems of focal dispersion, shadow areas and the like, has poor effect in complex structures, and the ray chromatography equation contains a large amount of 0 and is a pathological equation and is difficult to process; the second is fat ray chromatography, the sparsity of a chromatography inversion matrix is reduced through the expansion of the path of a correlation line, and the stability of a speed inversion result is improved, but the method lacks strict theoretical evidence; the third is full waveform inversion chromatography, based on a full waveform inversion theoretical frame, the early-arrival waveform inversion can theoretically complete near-surface high-precision speed modeling by utilizing the kinematics and dynamics information of the early-arrival characteristic waves, but has high calculation complexity and high dependence on an initial model, and is difficult to apply in actual production; and the fourth is Gaussian beam chromatography, and strict theory proves that the theoretical basis of the Gaussian beam method is a high-frequency approximation of the wave equation concentrated in a ray area by combining a dynamic ray tracing and paraxial approximation method, the calculation speed is high, the pathogenicity caused by too sparse of the chromatographic equation is reduced by constructing a plurality of chromatographic equations near the central ray, and the accuracy of the chromatographic inversion is improved.
In solving the Gaussian beam tomography kernel function, ray tracing is a key step, in the two-dimensional case, the ray tracing can quickly find a correct path through a dichotomy or other methods, and in the three-dimensional ray tracing, other quick ray tracing algorithms such as dichotomy are not applicable any more due to the existence of two parameters of azimuth angle and inclination angle. The conventional method for three-dimensional Gaussian beam chromatography center ray tracing is to conduct global scanning on the initial inclination angle and azimuth angle of a ray to obtain a ray path from a shot point to a detector point. However, on one hand, global scanning of the initial tilt angle and azimuth angle consumes a lot of time, and most of the calculation amount is invalid, so that the efficiency is low and the method is not reasonable; on the other hand, when the ray exceeds the velocity model boundary or reaches the set time, the ray stops the calculation, which also brings about a lot of ineffective calculation amount. Therefore, conventional three-dimensional ray tracing methods are not effectively applied in industry.
Disclosure of Invention
The invention mainly aims to provide a construction method and a storage medium of a near-surface velocity model, which are used for solving the problem of quickly constructing a near-surface Gaussian beam chromatography kernel function and the near-surface velocity model.
In a first aspect, embodiments of the present application provide a method for constructing a near-surface velocity model, including the following steps: s100: determining the surface initial speed of the ray according to a given initial speed model, determining an included angle between the ray and a specified reference line and a distance between the shot point and the detector in the horizontal direction according to the position of the shot point and the position of the detector, taking the included angle between the ray and the specified reference line as the azimuth angle of the ray, and determining the incidence inclination angle of the ray according to the given initial speed model and the distance between the shot point and the detector in the horizontal direction; s200: selecting a plurality of points on a ray path from a shot point to a detector point according to a preset time interval, and determining the position of each selected point on the ray path by utilizing a kinematic ray tracing method according to the initial speed, incidence inclination angle and azimuth angle of the surface of the ray; s300: judging whether each point selected on the ray path meets a preset angle condition according to the positions of the shot points and the positions of the detection points and the positions of each point selected on the ray path, if so, judging that the ray does not meet a preset stop calculation criterion, executing S400, otherwise, executing S500; s400: for each point selected on the ray path, judging whether the distance from each point selected on the ray path to the detection point is smaller than or equal to a preset distance threshold value according to the position of the detection point and the position of each point selected on the ray path, if so, executing S600, otherwise, executing S500; s500: respectively perturbing the incidence inclination angle and the azimuth angle of the ray by utilizing a preset perturbation angle, selecting a new group of incidence inclination angle and azimuth angle from the perturbed incidence inclination angle and azimuth angle, and returning to execute S200; s600: determining the ray velocity of each point according to the position of each point selected on a ray path by using an initial velocity model, determining plane wave solutions and spherical wave solutions corresponding to the moment of reaching each point selected on the ray path according to the ray velocity of each point selected on the ray path by using a dynamic ray tracing method, determining Gaussian beams according to the plane wave solutions and spherical wave solutions at each moment, determining a travel time chromatographic kernel according to the Gaussian beams, and determining near-surface velocity updating quantity according to the travel time chromatographic kernel; s700: updating a given initial velocity model according to the near-surface velocity updating quantity to obtain an updated velocity model; s800: determining a difference value between the travel time of the measured rays from the shot point to the detector point and the travel time determined according to a preset time interval, taking the updated velocity model as a new initial velocity model when the difference value does not meet a preset time condition, and returning to execute S100, and taking the velocity model updated this time as a final near-surface velocity model when the difference value meets the preset time condition.
In one embodiment, the initial velocity model is:
v(z)＝v 0 +gz
wherein z represents depth direction, g represents velocity gradient, v 0 Representing an initial surface velocity;
the azimuth angle of the ray is determined using:
wherein, 
representing the azimuth angle of the ray, (x) s ，y s ) Representing the plane coordinates of the shot point, (x) r ，y r ) Representing plane coordinates of the detector points;
the incidence angle of the ray is determined using:
wherein θ 0 Represents the incidence inclination of the ray, offset represents the horizontal distance between the shot point and the detector, g represents the velocity gradient, v 0 Representing the initial surface velocity.
In one embodiment, determining the location of each selected point on the ray path from the initial velocity, incidence angle, and azimuth angle of the surface of the ray using kinematic ray tracing includes: the position and ray velocity of each selected point on the ray is determined using:
wherein v (x 0 ,y 0 ,z 0 ) Representing the ray path coordinates (x) 0 ,y 0 ,z 0 ) Velocity at a point of (x), coordinates (x) 0 ,y 0 ,z 0 ) I.e., the point preceding the point on the ray path with coordinates (x, y, z),
represents the azimuth angle, theta, of the ray 0 Representing the incidence inclination of the rays, τ representing a preset time interval, p x 、p y And p z Slowness in the x, y and z directions, respectively.
In one embodiment, the preset angle condition includes: the included angle between the straight line direction from the shot point to each point selected on the ray path and the straight line direction from each point selected on the ray path to the detection point is smaller than or equal to a preset angle threshold.
In one embodiment, determining whether each selected point on the ray path satisfies a preset angle condition includes: and judging whether cosine values of included angles between the straight line direction from the shot point to each point selected on the ray path and the straight line direction from each point selected on the ray path to the detection point are larger than or equal to cosine values of a preset angle threshold value.
In one embodiment, using a dynamic ray tracing method, determining a plane wave solution and a spherical wave solution corresponding to a moment of arrival at each selected point on a ray path according to a ray velocity of each selected point on the ray path includes: the plane and sphere solutions corresponding to the moment of arrival at each selected point on the ray path are determined using:
wherein m and n respectively represent two mutually perpendicular directions in a plane perpendicular to the ray direction, V (a) represents the velocity at the selected point a on the ray path, V (a) represents the second partial derivative of the selected point a on the ray path with respect to the central ray coordinate system, and P and Q respectively represent the plane wave solution and the spherical wave solution.
In one embodiment, the Gaussian beam is determined from the plane and sphere solutions at each time instant using the following equation:
wherein u (Q, a, ω) represents a gaussian beam, a represents a selected point on the ray path, T (a) represents a propagation time of the wave along the central ray, Q represents a vector of the point a in the central ray coordinate system, i represents a complex unit, v (a) represents a velocity at the selected point a on the ray path, P and Q are solutions of dynamic ray tracing representing a plane wave solution and a spherical wave solution, respectively, ω represents a frequency;
determining a travel time chromatographic kernel from the gaussian beam using:
wherein omega 1 And omega 2 Respectively represent the lower limit and the upper limit of the frequency, P (omega) is a weight coefficient, and satisfies the relation
u (q, a, ω) representsA Gaussian beam;
determining a near-surface velocity update amount from the travel-time chromatographic kernel comprises: determining a slowness update amount by using the following method, wherein the reciprocal of the slowness update amount is taken as a near-surface speed update amount:
K T Δs＝Δt
wherein K is T Representing the travel time tomograph kernel, deltas representing the slowness update and deltat representing the difference between the travel time of the measured ray from the shot to the detector and the travel time determined from the preset time interval.
In one embodiment, the preset time condition includes: and for the difference value between the travel time of the measured rays from the shot point to the detection point and the travel time determined according to the preset time interval, the difference value calculated in the iteration and the difference value calculated last time is smaller than or equal to a preset time threshold value.
In one embodiment, S100 further comprises: and determining the elevation difference between the shot point and the demodulation point according to the position of the shot point and the position of the demodulation point, and adjusting the incidence dip angle of the ray according to the elevation difference between the shot point and the demodulation point.
In a second aspect, embodiments of the present application provide a storage medium storing a computer program which, when executed by a processor, implements the steps of the method of constructing a near-surface velocity model as described above.
The modeling method of the Gaussian beam near-surface velocity model based on the rapid ray tracing can be realized, the method utilizes the analytic solution of the ray path under the gradient velocity model to calculate the incidence angle at the shot points, carries out fine adjustment on the incidence angle according to the elevation difference between the shot points, calculates the azimuth relation of the shot points at the same time, finally obtains the initialization parameters of the correct ray incidence angle and azimuth angle, effectively reduces the scanning range of the two parameters, and finishes the calculation of the incorrect ray path according to the calculation stopping criterion, thereby rapidly and accurately obtaining the central ray path and travel time information from the seismic source to the detection points, and improving the calculation efficiency of the near-surface Gaussian beam chromatographic kernel function and the modeling speed of the near-surface velocity model.
Drawings
The accompanying drawings, which are included to provide a further understanding of the invention and are incorporated in and constitute a part of this specification, illustrate embodiments of the invention and together with the description serve to explain the invention and do not constitute a undue limitation on the invention, wherein:
FIG. 1 is a flow chart of a method of constructing a near-surface velocity model according to an embodiment of the present application;
FIG. 2 is a diagram of a given initial velocity model according to an embodiment of the present application;
FIG. 3A is a cross-sectional view of a Gaussian beam at Inline100 with a finite frequency of 1-50Hz according to a conventional method;
FIG. 3B is a cross-sectional view of a three-dimensional Gaussian beam at line100 obtained by a method of constructing a near-surface velocity model according to an embodiment of the application;
FIG. 4 is a graph comparing conventional ray tracing with improved method calculation efficiency.
Detailed Description
It should be noted that, in the case of no conflict, the embodiments and features in the embodiments may be combined with each other. The invention will be described in detail below with reference to the drawings in connection with embodiments.
Example 1
In the three-dimensional space, for a given shot point and a given detector point, the correct path is found, so that the inclination angle and the azimuth angle need to be searched in an omnibearing manner, but the process is a time-consuming and labor-consuming process, so that the initial incidence angle and the initial azimuth angle are reasonably set, and the calculation time of ray tracing can be effectively reduced.
FIG. 1 is a flow chart of a method of constructing a near-surface velocity model according to an embodiment of the present application. As shown in fig. 1, the present embodiment provides a method for constructing a near-surface velocity model, which includes the following steps:
s100: according to the given initial speed model, determining the surface initial speed of the ray, determining the included angle between the ray and the appointed datum line and the distance between the shot point and the datum line in the horizontal direction according to the position of the shot point and the position of the datum point, taking the included angle between the ray and the appointed datum line as the azimuth angle of the ray, and determining the incidence inclination angle of the ray according to the given initial speed model and the distance between the shot point and the datum line in the horizontal direction.
Wherein, the given initial velocity model may be a gradient velocity model, as shown in expression (1):
v(z)＝v 0 +gz (1)
wherein z represents depth direction, g represents velocity gradient, v 0 Representing the initial surface velocity.
For a given initial velocity model, for any point excitation, the ray path received by another point, we have analytical expression (2):
wherein x and z represent parameters in the analytical expression, θ, respectively 0 Representing the incidence angle of the ray.
Using the trigonometric relationship, the incidence inclination of a ray can be determined by expression (3):
wherein θ 0 Represents the incidence inclination of the ray, offset represents the horizontal distance between the shot point and the detector, g represents the velocity gradient, v 0 Representing the initial surface velocity.
Since the heights of the shot point and the demodulation point may be inconsistent, the height difference between the shot point and the demodulation point may be determined according to the position of the shot point and the position of the demodulation point, and the incidence inclination angle of the ray may be adjusted according to the height difference between the shot point and the demodulation point. If the shot point is higher than the detector, increasing the incidence dip angle, and if the shot point is lower than the detector, reducing the incidence dip angle, and specifically adjusting the incidence dip angle of the rays by using expressions (4) and (5):
θ＝θ 0 +Δθ (5)
wherein z is s Representing the shot point elevation, z r Representing the elevation of the detector.
The specified reference line may be a line pointing in the direction of north. Determining the azimuth of the ray by means of the shot-to-detector azimuth relationship using expression (6):
wherein, 
representing the azimuth angle of the ray, (x) s ，y s ) Representing the plane coordinates of the shot point, (x) r ，y r ) Representing the plane coordinates of the detector points. />
S200: and selecting a plurality of points on a ray path from a shot point to a detector point according to a preset time interval, and determining the position of each selected point on the ray path by utilizing a kinematic ray tracing method according to the initial speed, incidence inclination angle and azimuth angle of the surface of the ray.
The method comprises the steps of starting from a shot point, selecting points on a ray path reached by rays every preset time interval, and selecting a plurality of points on the ray path according to the method. Determining the position of each selected point on the ray path by using a kinematic ray tracing method according to the initial speed, incidence inclination angle and azimuth angle of the surface of the ray, wherein the method comprises the following steps: determining the position and ray velocity of each selected point on the ray path using expression (7):
wherein v (x 0 ,y 0 ,z 0 ) Representing the ray path coordinates (x) 0 ,y 0 ,z 0 ) At a point of (x) and the ray is at a coordinate (x 0 ,y 0 ,z 0 ) After a predetermined time interval, the point reaches a point with coordinates (x, y, z) of (x 0 ,y 0 ,z 0 ) I.e., the point preceding the point with coordinates (x, y, z),
represents the azimuth angle, theta, of the ray 0 Representing the incidence inclination of the rays, τ representing a preset time interval, p x 、p y And p z Slowness in the x, y and z directions, respectively. The kinematic ray tracing equation is calculated through the fourth-order Dragon-Geku tower, and the position coordinates of each point selected on the ray path from the shot point to the detector point according to the preset time interval can be obtained in sequence.
S300: for each point selected on the ray path, judging whether each point selected on the ray path meets the preset angle condition according to the positions of the shot point and the detector point and the positions of each point selected on the ray path, if so, judging that the ray does not meet the preset stop calculation criterion, executing S400, otherwise, executing S500.
The preset angle condition may include: the included angle between the straight line direction from the shot point to each point selected on the ray path and the straight line direction from each point selected on the ray path to the detection point is smaller than or equal to a preset angle threshold.
The analytical expression represented by the expression (2) is a circle, the center of the circle is higher than the connecting line between the detector point and the shot point, and according to the given gradient speed model (1), the path from the shot point to the detector point is a semicircle in the most extreme case. Therefore, if the included angle between any point in the shot-to-ray and the point-to-detector is smaller than or equal to 90 degrees, in an acceptable range, if the included angle is larger than 90 degrees, the calculation of the current ray needs to be stopped, and the incident inclination angle and/or the azimuth angle of the ray are changed and then recalculated. At this time, the preset angle threshold may be set to 90 °. However, the preset angle threshold may be appropriately increased for the case of a complex near-surface in the actual data, and may be set to 95 ° or 100 ° or the like, for example.
Specifically, determining whether each point selected on the ray path satisfies a preset angle condition may include: and judging whether cosine values of included angles between the straight line direction from the shot point to each point selected on the ray path and the straight line direction from each point selected on the ray path to the detection point are larger than or equal to cosine values of a preset angle threshold value.
When the preset angle threshold is 90 °, the calculation can be performed using the following formula:
wherein, 
representing the vector from the point on the ray to the shot point,/->
Representing the vector from the detector point onto the ray.
When the preset angle threshold is greater than 90 °, for example, the cosine value of the preset angle threshold is-0.2, the following formula may be used for calculation:
s400: for each point selected on the ray path, judging whether the distance from each point selected on the ray path to the detection point is smaller than or equal to a preset distance threshold value according to the position of the detection point and the position of each point selected on the ray path, executing S600 if the distance from each point selected on the ray path to the detection point is smaller than or equal to the preset distance threshold value, otherwise executing S500.
Specifically, the distance between the selected point on the ray path and the detector point can be calculated according to the position coordinates of each selected point on the ray path and the position coordinates of the detector point.
S500: and respectively perturbing the incidence inclination angle and the azimuth angle of the ray by utilizing a preset perturbation angle, selecting a new group of incidence inclination angle and azimuth angle from the perturbed incidence inclination angle and azimuth angle, and returning to the execution S200.
For S300, when there is a point on the ray path that does not satisfy the preset angle condition, that is, satisfies the stopping calculation criterion, or for S400, although each point selected on the ray path satisfies the preset angle condition, there is a point on the ray path that has a distance to the detector point greater than the preset distance threshold, the incident inclination angle and the azimuth angle of the ray are respectively disturbed by using the preset disturbance angle, and the disturbance range may be that
And (2) obtaining the incidence inclination angle range and the azimuth angle range of the rays, selecting a group of new incidence inclination angles and azimuth angles from the incidence inclination angle range and the azimuth angle range, and returning to the S200 until the incidence inclination angles and the azimuth angles of the rays which simultaneously meet the preset angle condition and the preset distance condition are determined.
S600: determining the ray velocity of each point according to the position of each point selected on the ray path by using an initial velocity model, determining a plane wave solution and a spherical wave solution corresponding to the moment of reaching each point selected on the ray path according to the ray velocity of each point selected on the ray path by using a dynamic ray tracing method, determining Gaussian beams according to the plane wave solution and the spherical wave solution at each moment, determining a travel time chromatographic kernel according to the Gaussian beams, and determining a near-surface velocity update quantity according to the travel time chromatographic kernel.
For a ray path determined by an incidence inclination angle and an azimuth angle of rays simultaneously meeting a preset angle condition and a preset distance condition, determining a plane wave solution and a spherical wave solution corresponding to a moment of reaching each point selected on the ray path according to a ray speed of each point selected on the ray path by utilizing a dynamic ray tracing method, wherein the method comprises the following steps:
determining a plane wave solution and a spherical wave solution corresponding to a moment of arrival at each selected point on the ray path using expression (10):
wherein m and n respectively represent two mutually perpendicular directions in a plane perpendicular to the ray direction, V (a) represents the speed at the selected point a on the ray path, V (a) represents the second partial derivative of the selected point a on the ray path relative to the central ray coordinate system, and P and Q are solutions of dynamic ray tracing respectively represent a plane wave solution and a spherical wave solution, and the values of P and Q corresponding to each selected point on the ray path can be calculated through a fourth-order longgrid tower. Wherein,
the expression (12) can be used to translate from rectangular to central ray coordinate system:
wherein t represents the tangential direction of the ray, and m and n represent two mutually perpendicular directions in a plane perpendicular to the direction of the ray, respectively.
Determining Gaussian beams according to plane wave solutions and spherical wave solutions at various moments by using an expression (13):
where u (Q, a, ω) represents a gaussian beam, a represents a selected point on the ray path, T (a) represents a propagation time of the wave along the central ray, Q represents a vector of the point a in the central ray coordinate system, i represents a complex unit, v (a) represents a velocity at the selected point a on the ray path, P and Q represent a plane solution and a spherical solution, respectively, and ω represents a frequency.
Gaussian beams are high frequency progressive solutions where the elastohydrodynamic equations are concentrated near the rays, which have a uniform form in the frequency domain as shown in expression (13), both in two and in three dimensions.
Determining a travel-time chromatographic kernel from the gaussian beam using expression (14):
wherein omega 1 And omega 2 Respectively representing the lower limit and the upper limit of the frequency, P (omega) is a weight coefficient, and P (omega) satisfies the relation
u (q, a, ω) represents a gaussian beam. Expression (13) is a single-frequency Gaussian beam, and expression (14) is a band-limited Gaussian beam travel-time chromatographic kernel function.
Determining a near-surface velocity update amount from the travel-time chromatographic kernel comprises: determining a slowness update amount by using the expression (15), and taking the reciprocal of the slowness update amount as a near-surface velocity update amount:
K T Δs＝Δt (15)
wherein K is T Representing the travel time tomograph kernel, deltas representing the slowness update and deltat representing the difference between the travel time of the measured ray from the shot to the detector and the travel time determined from the preset time interval.
S700: and updating the given initial velocity model according to the near-surface velocity updating quantity to obtain an updated velocity model.
Specifically, the near-surface velocity update amount may be added to the right side of expression (1), resulting in expression (16):
and simultaneously storing and calculating to obtain the correct incidence inclination angle and azimuth angle as the initial angle of the next iteration.
S800: determining a difference value between the travel time of the measured rays from the shot point to the detector point and the travel time determined according to a preset time interval, taking the updated velocity model as a new initial velocity model when the difference value does not meet a preset time condition, and returning to execute S100, and taking the velocity model updated this time as a final near-surface velocity model when the difference value meets the preset time condition. The length of each preset time interval on the ray path is accumulated, and the travel time is obtained.
Wherein, the preset time condition includes: and for the difference value between the travel time of the measured rays from the shot point to the detection point and the travel time determined according to the preset time interval, the difference value calculated in the iteration and the difference value calculated last time is smaller than or equal to a preset time threshold value.
According to the invention, by determining tomographic inversion parameters (an initial gradient model, an azimuth angle and a disturbance maximum angle), calculating an initial inclination angle of rays by utilizing an analytic solution of rays of an initial gradient velocity model, adjusting the initial inclination angle according to an elevation relation between a shot point and a detector, calculating the initial azimuth angle of the rays by utilizing a coordinate relation of the shot point and the detector, if the rays meet a calculation stopping criterion, carrying out disturbance on the azimuth angle or the inclination angle, and if the rays do not meet the calculation stopping criterion, and the distance from any point on a ray path to the detector meets a preset distance condition, obtaining a correct ray path, and storing the current incident inclination angle and azimuth angle of the rays as the initial angle of the next iteration.
According to the method, the initial azimuth angle and the initial inclination angle are calculated, the calculation range is reduced, meanwhile, calculation is stopped in time aiming at an incorrect ray path, the correct azimuth angle and the correct inclination angle are found through local optimization, the calculation efficiency is effectively improved, and the modeling method of the Gaussian beam near-surface velocity model based on rapid ray tracing can be realized.
Example two
The initial velocity model in which the initial surface velocity of the ray is 1500m/s, the initial inclination angle and azimuth angle of the ray from the shot point to the detector point are calculated using expressions (4), (5) and (6), and then the Gaussian beam tomographic kernel function is calculated, as shown in FIG. 2, in which the ray has no velocity change in the horizontal direction and the velocity gradient in the depth direction is 1. Fig. 3A is a cross-sectional view of a gaussian beam obtained by a conventional method at an line100, and fig. 3B is a cross-sectional view of a gaussian beam obtained by fast ray tracing according to the present application at the line 100. Compared with the traditional method, the method has the advantages that the result obtained by calculation is almost the same as that obtained by the traditional method. As shown in fig. 4, the method of the present application calculates faster.
Example III
The present embodiment provides a storage medium storing a computer program which, when executed by a processor, implements the steps of the method of constructing a near-surface velocity model as described above.
It is noted that the terms used herein are used merely to describe particular embodiments and are not intended to limit exemplary embodiments in accordance with the present application and when the terms "comprises" and/or "comprising" are used in this specification they specify the presence of stated features, steps, operations, devices, components, and/or combinations thereof.
It should be noted that the terms "first," "second," and the like in the description and claims of the present application and in the drawings are used for distinguishing between similar objects and not for describing a particular sequential or chronological order. It is to be understood that the terms so used are interchangeable under appropriate circumstances such that the embodiments of the application described herein are, for example, capable of operation in sequences other than those illustrated or otherwise described herein.
It should be understood that the exemplary embodiments in this specification may be embodied in many different forms and should not be construed as limited to only the embodiments set forth herein. These embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the concept of these exemplary embodiments to those skilled in the art, and should not be construed as limiting the invention.
Claims (10)
1. The construction method of the near-surface velocity model is characterized by comprising the following steps of:
s100: determining the surface initial speed of the ray according to a given initial speed model, determining an included angle between the ray and a specified reference line and a distance between the shot point and the detector in the horizontal direction according to the position of the shot point and the position of the detector, taking the included angle between the ray and the specified reference line as the azimuth angle of the ray, and determining the incidence inclination angle of the ray according to the given initial speed model and the distance between the shot point and the detector in the horizontal direction;
s200: selecting a plurality of points on a ray path from a shot point to a detector point according to a preset time interval, and determining the position of each selected point on the ray path by utilizing a kinematic ray tracing method according to the initial speed, incidence inclination angle and azimuth angle of the surface of the ray;
s300: judging whether each point selected on the ray path meets a preset angle condition according to the positions of the shot points and the positions of the detection points and the positions of each point selected on the ray path, if so, judging that the ray does not meet a preset stop calculation criterion, executing S400, otherwise, executing S500;
s400: for each point selected on the ray path, judging whether the distance from each point selected on the ray path to the detection point is smaller than or equal to a preset distance threshold value according to the position of the detection point and the position of each point selected on the ray path, if so, executing S600, otherwise, executing S500;
s500: respectively perturbing the incidence inclination angle and the azimuth angle of the ray by utilizing a preset perturbation angle, selecting a new group of incidence inclination angle and azimuth angle from the perturbed incidence inclination angle and azimuth angle, and returning to execute S200;
s600: determining the ray velocity of each point according to the position of each point selected on a ray path by using an initial velocity model, determining plane wave solutions and spherical wave solutions corresponding to the moment of reaching each point selected on the ray path according to the ray velocity of each point selected on the ray path by using a dynamic ray tracing method, determining Gaussian beams according to the plane wave solutions and spherical wave solutions at each moment, determining a travel time chromatographic kernel according to the Gaussian beams, and determining near-surface velocity updating quantity according to the travel time chromatographic kernel;
s700: updating a given initial velocity model according to the near-surface velocity updating quantity to obtain an updated velocity model;
s800: determining a difference value between the travel time of the measured rays from the shot point to the detector point and the travel time determined according to a preset time interval, taking the updated velocity model as a new initial velocity model when the difference value does not meet a preset time condition, and returning to execute S100, and taking the velocity model updated this time as a final near-surface velocity model when the difference value meets the preset time condition.
2. The method of claim 1, wherein the initial velocity model is:
v(z)＝v 0 +gz
wherein z represents depth direction, g represents velocity gradient, v 0 Representing an initial surface velocity;
the azimuth angle of the ray is determined using:
wherein, 
representing the azimuth angle of the ray, (x) s ，y s ) Representing the plane coordinates of the shot point, (x) r ，y r ) Representing plane coordinates of the detector points;
the incidence angle of the ray is determined using:
wherein θ 0 Represents the incidence inclination of the ray, offset represents the horizontal distance between the shot point and the detector, g represents the velocity gradient, v 0 Representing the initial surface velocity.
3. The method of claim 1, wherein determining the location of each selected point on the ray path using a kinematic ray tracing method based on the initial surface velocity, incidence angle, and azimuth angle of the ray comprises:
the position and ray velocity of each selected point on the ray is determined using:
wherein v (x 0 ,y 0 ,z 0 ) Representing the ray path coordinates (x) 0 ,y 0 ,z 0 ) Velocity at a point of (x), coordinates (x) 0 ,y 0 ,z 0 ) I.e., the point preceding the point on the ray path with coordinates (x, y, z),
represents the azimuth angle of the ray, θ represents the incidence inclination angle of the ray, τ represents a preset time interval, p x 、p y And p z Slowness in the x, y and z directions, respectively.
4. The method for constructing a near-surface velocity model according to claim 1, wherein the preset angle condition includes: the included angle between the straight line direction from the shot point to each point selected on the ray path and the straight line direction from each point selected on the ray path to the detection point is smaller than or equal to a preset angle threshold.
5. The method for constructing a near-surface velocity model according to claim 1, wherein determining whether each point selected on the ray path satisfies a preset angle condition comprises:
and judging whether cosine values of included angles between the straight line direction from the shot point to each point selected on the ray path and the straight line direction from each point selected on the ray path to the detection point are larger than or equal to cosine values of a preset angle threshold value.
6. The method for constructing a near-surface velocity model according to claim 1, wherein determining a plane wave solution and a spherical wave solution corresponding to a moment of arrival at each selected point on the ray path from the ray velocity of each selected point on the ray path by using a dynamic ray tracing method, comprises:
the plane and sphere solutions corresponding to the moment of arrival at each selected point on the ray path are determined using:
wherein v (a) represents the velocity at the selected point a on the ray path, P and Q are solutions for dynamic ray tracing, respectively represent a plane wave solution and a spherical wave solution,
v (a) represents the second partial derivative of the selected point a on the ray path with respect to the central ray coordinate system, m and n representing two mutually perpendicular directions in a plane perpendicular to the ray direction, respectively.
7. The method of constructing a near-surface velocity model according to claim 1, wherein the gaussian beam is determined from the plane and spherical solutions at each time by using the following formula:
wherein u (Q, a, ω) represents a gaussian beam, a represents a selected point on the ray path, T (a) represents a propagation time of the wave along the central ray, Q represents a vector of the point a in the central ray coordinate system, i represents a complex unit, v (a) represents a velocity at the selected point a on the ray path, P and Q represent a plane solution and a spherical solution, respectively, ω represents a frequency;
determining a travel time chromatographic kernel from the gaussian beam using:
wherein omega 1 And omega 2 Respectively represent the lower limit and the upper limit of the frequency, P (omega) is a weight coefficient, and satisfies the relation
u (q, a, ω) represents a gaussian beam;
determining a near-surface velocity update amount from the travel-time chromatographic kernel comprises: determining a slowness update amount by using the following method, wherein the reciprocal of the slowness update amount is taken as a near-surface speed update amount:
K T Δs＝Δt
wherein K is T Representing the travel time tomograph kernel, deltas representing the slowness update and deltat representing the difference between the travel time of the measured ray from the shot to the detector and the travel time determined from the preset time interval.
8. The method for constructing a near-surface velocity model according to claim 1, wherein the preset time condition includes: and for the difference value between the travel time of the measured rays from the shot point to the detection point and the travel time determined according to the preset time interval, the difference value calculated in the iteration and the difference value calculated last time is smaller than or equal to a preset time threshold value.
9. The method of constructing a near-surface velocity model according to claim 1, wherein S100 further comprises: and determining the elevation difference between the shot point and the demodulation point according to the position of the shot point and the position of the demodulation point, and adjusting the incidence dip angle of the ray according to the elevation difference between the shot point and the demodulation point.
10. A storage medium storing a computer program, wherein the computer program, when executed by a processor, implements the steps of the method of constructing a near-surface velocity model according to any one of claims 1 to 9.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011004206.3A CN114428292B (en) | 2020-09-22 | 2020-09-22 | Method for constructing near-surface velocity model and storage medium |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011004206.3A CN114428292B (en) | 2020-09-22 | 2020-09-22 | Method for constructing near-surface velocity model and storage medium |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114428292A CN114428292A (en) | 2022-05-03 |
| CN114428292B true CN114428292B (en) | 2023-06-02 |
Family
ID=81308855
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011004206.3A Active CN114428292B (en) | 2020-09-22 | 2020-09-22 | Method for constructing near-surface velocity model and storage medium |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114428292B (en) |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109581495A (en) * | 2018-10-24 | 2019-04-05 | 中国石油天然气集团有限公司 | Mountainous region surface velocity model construction method and system |
Family Cites Families (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9013956B2 (en) * | 2009-10-27 | 2015-04-21 | Chevron U.S.A Inc. | Method and system for seismic imaging and earth modeling using beam tomography |
| US10162070B2 (en) * | 2012-04-05 | 2018-12-25 | Westerngeco L.L.C. | Converting a first acquired data subset to a second acquired data subset |
| CN102841376A (en) * | 2012-09-06 | 2012-12-26 | 中国石油大学(华东) | Retrieval method for chromatography speed based on undulating surface |
| CN106842313B (en) * | 2015-12-04 | 2021-04-16 | 中国石油化工股份有限公司 | Anisotropic parameter inversion method based on azimuth pre-stack seismic data |
| CN105549081B (en) * | 2016-01-29 | 2018-09-11 | 中国石油大学(华东) | Anisotropic medium is total to big gun domain Gaussian beam offset imaging method |
| KR101656862B1 (en) * | 2016-03-15 | 2016-09-13 | 한국지질자원연구원 | Apparatus and method for performing stochastic modeling of earthquake fault rupture |
| CN106597540B (en) * | 2016-12-30 | 2017-10-31 | 中国科学院地质与地球物理研究所 | Gaussian beam offset imaging method and device |
| CN106772583B (en) * | 2017-01-10 | 2018-09-04 | 中国科学院地质与地球物理研究所 | A kind of earthquake diffracted wave separation method and device |
| CN108303736B (en) * | 2017-12-07 | 2020-11-17 | 东华理工大学 | Ray tracing forward method for shortest path of anisotropic TI medium |
| CN111638551A (en) * | 2019-03-01 | 2020-09-08 | 中国石油天然气集团有限公司 | Seismic first-motion wave travel time chromatography method and device |
-
2020
- 2020-09-22 CN CN202011004206.3A patent/CN114428292B/en active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109581495A (en) * | 2018-10-24 | 2019-04-05 | 中国石油天然气集团有限公司 | Mountainous region surface velocity model construction method and system |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114428292A (en) | 2022-05-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106054134B (en) | A kind of method for rapidly positioning based on TDOA | |
| JP6896180B2 (en) | Wind flow detection system and method for finding the velocity field of wind flow | |
| CN106912105B (en) | Three-dimensional positioning method based on PSO _ BP neural network | |
| KR101914550B1 (en) | Method for tracking target position of radar | |
| CN105136054A (en) | Fine structure deformation monitoring method and system based on ground three-dimensional laser scanning | |
| CN110132281B (en) | Underwater high-speed target high-precision autonomous acoustic navigation method based on inquiry response mode | |
| CN108614268A (en) | The acoustics tracking of low altitude high speed airbound target | |
| WO2019242045A1 (en) | Method for calculating virtual source two-dimensional wavefront construction seismic wave travel time | |
| CN108896053B (en) | Planet landing optical navigation optimal road sign selection method | |
| CN109459787B (en) | coal mine underground structure imaging method and system based on seismic channel wave full-waveform inversion | |
| CN109814110A (en) | The method of structuring the formation of deep-sea Long baselines positioning formation topological structure | |
| CN108303736B (en) | Ray tracing forward method for shortest path of anisotropic TI medium | |
| CN103307968A (en) | Method for detecting posture of robot carrying platform | |
| WO2020192182A1 (en) | Indoor positioning method and system, and electronic device | |
| Spencer | Closed-form analytical solutions of the time difference of arrival source location problem for minimal element monitoring arrays | |
| CN114428292B (en) | Method for constructing near-surface velocity model and storage medium | |
| CN105573963A (en) | Reconstruction method for horizontal nonuniform structure of ionized layer | |
| JP2020063958A (en) | Position estimation device and method | |
| CN112596113A (en) | Method for identifying field source position based on intersection points of characteristic values of different gradients of gravity | |
| CN116660980A (en) | Microseism positioning method based on improved path function equation | |
| US20150316671A1 (en) | Migration Velocity Analysis Method for VSP Data | |
| CN109581494B (en) | Pre-stack migration method and device | |
| CN115826038A (en) | VTI medium adjoint state method travel time multi-parameter tomography method and system | |
| CN113608262B (en) | Seismic data processing method and device for calculating rotation component by using translation component | |
| CN112257241B (en) | Triangular net Fresnel time difference tomography inversion method |
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 | ||
| TR01 | Transfer of patent right | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20231106 Address after: No. 219 Shanggao Road, Dongshan Street, Jiangning District, Nanjing City, Jiangsu Province, 211103 Patentee after: Sinopec Petroleum Geophysical Exploration Technology Research Institute Co.,Ltd. Address before: 100728 No. 22 North Main Street, Chaoyang District, Beijing, Chaoyangmen Patentee before: CHINA PETROLEUM & CHEMICAL Corp. Patentee before: SINOPEC GEOPHYSICAL Research Institute |