CN112505765B - Method for scanning travel time of seismic waves by using Lax Friedrichs - Google Patents

Method for scanning travel time of seismic waves by using Lax Friedrichs Download PDF

Info

Publication number
CN112505765B
CN112505765B CN202011294721.XA CN202011294721A CN112505765B CN 112505765 B CN112505765 B CN 112505765B CN 202011294721 A CN202011294721 A CN 202011294721A CN 112505765 B CN112505765 B CN 112505765B
Authority
CN
China
Prior art keywords
travel time
equation
lax
friedrichs
factor
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
Application number
CN202011294721.XA
Other languages
Chinese (zh)
Other versions
CN112505765A (en
Inventor
黄光南
李红星
张华�
李泽林
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
East China Institute of Technology
Original Assignee
East China Institute of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by East China Institute of Technology filed Critical East China Institute of Technology
Priority to CN202011294721.XA priority Critical patent/CN112505765B/en
Publication of CN112505765A publication Critical patent/CN112505765A/en
Application granted granted Critical
Publication of CN112505765B publication Critical patent/CN112505765B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/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
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/282—Application of seismic models, synthetic seismograms

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

The invention relates to a method for scanning travel time of seismic waves by using Lax Friedrichs, which comprises the following steps: by introducing a slowness vector, establishing a program function equation of the seismic waves in the anisotropic TI medium; decomposing the solution of the equation into products of two factors, and bringing the products of the two factors into the equation to be converted into corresponding Hamiltonian forms; solving group velocity and phase velocity according to the slowness vector; and solving a equation by using Lax-Friedrichs scanning format values, and determining the travel time of the seismic waves. The method, the device and the storage medium for scanning the travel time of the seismic waves by the Lax Friedrichs adopt a factorization and rapid scanning mode, ensure the accuracy and the rapidity of travel time calculation, and facilitate effective geological investigation.

Description

Method for scanning travel time of seismic waves by using Lax Friedrichs
Technical Field
The invention relates to the technical field of geological exploration, in particular to a method for scanning travel time of seismic waves by using Lax Friedrichs.
Background
The travel time of a seismic wave refers to the time that it takes for the seismic wave to travel from the source to the observation point, and is also known as travel time. Usually, the travel time of seismic waves is recorded through a ground sensor, the real position of the underground geologic body is truly inverted, and then the physical properties of underground rock are obtained through inversion. Therefore, the calculation of travel time of the seismic waves can truly reflect the geological conditions, and is a key technology in the technical field of geological exploration.
Seismic waves have very different kinematic and dynamic characteristics when propagating in isotropic and anisotropic media. Only longitudinal and transverse waves are present in the isotropic medium, which involves three modes in the anisotropic medium, including one quasi-longitudinal wave (qP) and two quasi-transverse waves (qSV and qSH). Each mode has a respective mode and polarization direction. In the prior art, a finite difference-based program solver is the most stable and efficient travel time calculation method, however, a method for using an anisotropic medium to perform the finite difference program solver is lacked, and meanwhile, the source singularities cause larger errors at source points, the errors can be spread to the whole area, and the accuracy is lower.
In fact, more and more commercial seismic processing software (such as Landmark, CGG, omega, etc.) takes into account the anisotropy. In addition, the equation for controlling weak anisotropy is much simpler than the equation for controlling strong anisotropy, and thus, weak anisotropy assumptions are widely used in various techniques of seismic data processing. This assumption makes the techniques of seismic data processing easier to implement. However, in the prior art, the method for performing the travel time calculation by using weak anisotropy often lacks accuracy and rapidity, and how to propose a method for performing the travel time calculation rapidly and accurately is a problem to be solved.
Disclosure of Invention
In view of the foregoing, it is necessary to provide a method for scanning travel time of seismic waves by using the Lax Friedrichs, so as to solve the problem of how to provide a method for rapidly and accurately calculating travel time.
The invention provides a method for scanning travel time of seismic waves by using Lax Friedrichs, which comprises the following steps:
by introducing a slowness vector, establishing a program function equation of the seismic waves in the anisotropic TI medium;
decomposing the solution of the equation into products of two factors, and bringing the products of the two factors into the equation to be converted into corresponding Hamiltonian forms;
solving group velocity and phase velocity according to the slowness vector;
and solving the equation of the equation by utilizing Lax-Friedrichs scanning format values, and determining the travel time of the seismic waves.
Further, decomposing the solution of the equation of the program into a first factor and a second factor, the first factor and the second factor being expressed as:
Figure BDA0002784893810000021
wherein x, y, z are coordinates in the computational domain, T (x, y, z) is the first factor, s (x, y, z) is the second factor, T 0 (x, y, z) is a predetermined capture source singularity, u (x, y, z) is an unknown factor smoothed around the source point, s 0 (x, y, z) is a phase slowness model, and α (x, y, z) is the phase slowness model s 0 (x, y, z) a corresponding location factor;
bringing the first and second factors into the equation and based on the unknown factors u (x, y, z) and T 0 (x, y, z) performing formula conversion on the derivatives of the coordinates x, y, z, and determining the corresponding hamiltonian form as:
Figure BDA0002784893810000022
f(x,y,z)＝s m (P,Q,R),(m＝1,2,3)
wherein ,(ux ,u y ,u z ) Is the derivative of the unknown factor u (x, y, z) with respect to the coordinates x, y, z, H (x, y, z, u) x ,u y ,u z ) Is Hamiltonian, f (x, y, z) is corresponding Hamiltonian derivative function, s m (P, Q, R) is the phase slowness model after it is brought into P, Q, R, P, Q and R are the derivatives of the first factor T (x, y, z) with respect to coordinates x, y, z.
Further, the solving for the group velocity and the phase velocity based on the slowness vector comprises:
setting a phase slowness direction and a ray direction;
and determining the group velocity and the phase velocity according to the relation between the slowness vector and the group velocity vector.
Further, in the corresponding hamiltonian form, the determination of the derivative of the travel time with respect to the coordinate axes x, y, z comprises:
determining an analytical expression of the derivative of the travel time with respect to the coordinate axes x, y, z, expressed as:
Figure BDA0002784893810000031
wherein ,v0 As a function of the phase velocity at the source point, (p) 0 ,q 0 ,r 0 )＝(T 0x ,T 0y ,T 0z ) Is the derivative of the travel time with respect to x, y, z;
sequentially extracting derivatives p of the travel time with respect to x, y, z from square roots of the analytical expression 0 ,q 0 ,r 0 And determining the derivative p of said travel time with respect to x, y, z by algebraic operation 0 ,q 0 ,r 0 Expressed as the corresponding absolute value of:
Figure BDA0002784893810000032
Figure BDA0002784893810000033
Figure BDA0002784893810000034
wherein ,
Figure BDA0002784893810000035
is the inclination angle and azimuth angle of the phase slowness direction.
Further, the solving the equation by using the Lax-Friedrichs scan format value includes:
establishing a first-order Lax-Friedrichs format of the corresponding Hamiltonian form;
replacing travel time values on the grid points according to the third-order windward bias WENO approximate value corresponding to each grid point;
writing the static Hamiltonian-Jacobian equation in the WENO-based third-order Lax-Friedrichs format into a compact form, and updating the travel time value on the grid point according to the compact form.
Further, the updating travel time values at the grid points according to the compact form includes:
initializing and assigning grid points in the first-order Lax-Friedrichs format;
solving the compact form by using Gauss-Seidel of 8 alternate direction sweep frequency iterations, and updating the grid;
setting a plurality of corresponding maximum and minimum value conditions on six planes of an inner boundary and an outer boundary, and determining a numerical solution of an out-of-domain point by adopting a linear extrapolation method;
and judging whether the L1 norm difference of the solution between two continuous iterations meets a preset precision condition, and if so, terminating the iteration.
Compared with the prior art, the invention has the beneficial effects that: the first arrival travel time of the three modes, namely qP-, qSV-, qSH-waves, can be calculated and the upwind singularities of the finite difference program solver are processed by a factorization method to decompose the solution of the general program equation into the product of two factors, one of which is the travel time of the homogeneous medium, this factor captures the upwind singularities and the other factor (the basis function) is a necessary modification near the source point. Further, a Lax-Friedrichs finite difference format is adopted, iteration convergence can be guaranteed no matter whether the Hamiltonian is convex or concave, a rapid sweep frequency format based on a third-order weighted intrinsic non-oscillation format (WENO) is adopted to carry out numerical calculation on the basis function, the iteration times are irrelevant to the size of the grid, and as the size of the grid approaches 0, the numerical solution converges on the expected weak solution. In conclusion, the method, the device and the storage medium for scanning the travel time of the seismic waves by the Lax Friedrichs adopt the mode of factorization and rapid scanning, ensure the accuracy and the rapidity of travel time calculation, and facilitate effective geological investigation.
Drawings
FIG. 1 is a flow chart of a method for Lax Friedrichs to scan travel time of seismic waves;
FIG. 2 is a schematic flow chart of factorization provided by the present invention;
FIG. 3 is a schematic diagram of a flow chart for solving group velocity and phase velocity provided by the present invention;
FIG. 4 is a schematic diagram of a Lax-Friedrichs scanning flow provided by the invention;
FIG. 5 is a flowchart of updating a grid according to the present invention;
FIG. 6 is a schematic diagram of the travel time of a qP wave provided by the present invention;
FIG. 7 is a schematic diagram of the travel time of qSV wave provided by the present invention;
FIG. 8 is a schematic diagram of the travel time of qSH wave provided by the present invention;
FIG. 9 is a graph showing the travel time of a qP wave provided by the present invention;
FIG. 10 is a timing diagram of qSV provided by the present invention;
FIG. 11 is a timing diagram of qSH provided by the present invention;
FIG. 12 is a schematic diagram of anisotropic parameters provided by the present invention;
FIG. 13 is a schematic diagram showing the travel time of qP waves under Lax-Friedrichs scanning provided by the invention;
FIG. 14 is a schematic diagram showing the travel time of qSV waves under Lax-Friedrichs scanning provided by the invention;
FIG. 15 is a schematic diagram showing the travel time of qSH waves under Lax-Friedrichs scanning provided by the invention;
FIG. 16 is a three-dimensional travel time contour plot of a qP wave provided by the present invention;
FIG. 17 is a three-dimensional travel time contour plot of qSV waves provided by the present invention;
FIG. 18 is a three-dimensional travel time contour plot of qSH waves provided by the present invention;
FIG. 19 is a schematic diagram of a device for Lax Friedrichs to scan travel time of seismic waves.
Detailed Description
Preferred embodiments of the present invention will now be described in detail with reference to the accompanying drawings, which form a part hereof, and together with the description serve to explain the principles of the invention, and are not intended to limit the scope of the invention.
Example 1
The embodiment of the invention provides a method for scanning travel time of seismic waves by using Lax Friedrichs, and referring to FIG. 1, FIG. 1 is a flow chart of the method for scanning travel time of seismic waves by using Lax Friedrichs, wherein the method for scanning travel time of seismic waves by using Lax Friedrichs comprises steps S1 to S4, and the method comprises the following steps:
in step S1, establishing a program function equation of the seismic waves in the anisotropic TI medium by introducing a slowness vector;
in step S2, decomposing the solution of the equation into products of two factors, and bringing the products of the two factors into the equation to be converted into a corresponding Hamiltonian formula;
in step S3, solving group velocity and phase velocity according to the slowness vector;
in step S4, the equation is solved by using Lax-Friedrichs scan format values.
Preferably, the step S1 specifically includes the following steps:
according to the kristolochia equation, the slowness vector is expressed as:
Figure BDA0002784893810000061
wherein ,
Figure BDA0002784893810000062
v is the phase velocity of the seismic waves, which include qP waves, qSV waves and qSH waves, as the normal vector of the wavefront;
from the slowness vector, a program function equation is established expressed as:
Figure BDA0002784893810000063
wherein, T is the travel time,
Figure BDA0002784893810000064
is the reciprocal of travel time and is equal to the slowness vector, v m (m=1, 2, 3) is the phase velocities of qP wave, qSV wave, and qSH wave, expressed by the following formulas, respectively:
Figure BDA0002784893810000071
Figure BDA0002784893810000072
wherein ,v1,2 The phase velocity of qP wave and qSV wave, v 3 A is the phase velocity of qSH wave 44 and a66 For the elastic modulus of the anisotropic TI medium,
Figure BDA0002784893810000073
is formed by the phase slowness direction and the symmetry axis of the anisotropic TI medium, expressed as:
Figure BDA0002784893810000074
wherein ,
Figure BDA0002784893810000075
tilt angle and azimuth angle of symmetry axis for anisotropic TI medium, θ and +.>
Figure BDA0002784893810000076
The tangential angle and the azimuth angle of the phase slowness direction are expressed as follows:
Figure BDA0002784893810000077
Figure BDA0002784893810000078
wherein ,P＝Tx ，R＝T z and Q＝Ty ，T x 、T z and Ty Is the partial derivative of the travel time T.
Wherein M and N are represented by the following formula:
M＝0.5(Q 1 +Q 2 )
N＝Q 1 Q 2 -Q 3 (7)
wherein ,Q1 、Q 2 、Q 3 Expressed by the following formula:
Figure BDA0002784893810000079
wherein ,(a11 ，a 22 ，a 33 ，a 44 ，a 55 ，a 66 ) Is the elastic modulus of the anisotropic TI medium.
Preferably, as seen in fig. 2, fig. 2 is a schematic flow chart of factorization provided in the present invention, and step S2 includes steps S21 to S22, wherein:
in step S21, decomposing the solution of the equation of the program into a first factor and a second factor;
in step S22, the first and second factors are brought into a program function equation and are based on the unknown factors u (x, y, z) and T 0 And (x, y, z) carrying out formula conversion on the derivatives of the coordinates x, y and z to determine the corresponding Hamiltonian form.
In a specific embodiment of the present invention, the solution of formula (2) is decomposed into the product of two factors, expressed as:
Figure BDA0002784893810000081
wherein x, y, z are coordinates in the computational domain, T (x, y, z) is the first factor, s (x, y, z) is the second factor, T 0 (x, y, z) is a predetermined capture source singularity, u (x, y, z) is an unknown factor smoothed around the source point, s 0 (x, y, z) is a phase slowness model, and α (x, y, z) is the phase slowness model s 0 (x, y, z), the first factor capturing source singularities by parsing or numerical calculation, the second factor being the necessary correction to the first factor, being smooth near the source point. Considering the decomposition of equation (2), in a homogeneous medium, the equation satisfies the following relationship:
Figure BDA0002784893810000082
wherein ,T0 Is predetermined to capture source singularities, u is an unknown factor that is smoothed around the source. s is a phase slowness model, s 0 Is the source point (x 0 ,y 0 ,z 0 ) Is a phase slowness of (a). They are all related to the phase slowness direction
Figure BDA0002784893810000083
Is a function of (2).
Substituting the above formula (9) into the above formula (2), the obtained result is expressed as the following formula:
Figure BDA0002784893810000084
with p=t x ，R＝T z and Q＝Ty The calculation was carried out and the result was expressed as follows:
P＝T 0x u+T 0 p (12-a)
Q＝T 0y u+T 0 q (12-b)
R＝T 0z u+T 0 r (12-c)
wherein, (p, q, r) = (u) x ,u y ,u z ) Is the derivative of u with respect to x, y, z, T 0x ，T 0y and T0z Is T 0 Regarding the derivatives of x, y, z, further, equation (11) can be written in hamiltonian, expressed as:
Figure BDA0002784893810000091
f(x,y,z)＝s m (P,Q,R),(m＝1,2,3) (13-b)
wherein ,(ux ,u y ,u z ) Is the derivative of the unknown factor u (x, y, z) with respect to the coordinates x, y, z, H (x, y, z, u) x ,u y ,u z ) Is Hamiltonian, f (x, y, z) is corresponding Hamiltonian derivative function, s m (P, Q, R) is a phase slowness model after being brought into P, Q, R, P, Q and R is the derivative of the first factor T (x, y, z) with respect to coordinates x, y, z, T 0 (x, y, z) is a predetermined capture source singularities, u (x, y, z) is an unknown factor smoothed around the source points.
Preferably, as seen in conjunction with fig. 3, fig. 3 is a schematic flow chart for solving group velocity and phase velocity provided in the present invention, and step S3 includes steps S31 to S32, wherein:
in step S31, a phase slowness direction and a radial direction are set;
in step S32, the group velocity and the phase velocity are determined from the relationship between the slowness vector and the group velocity vector.
Thus, the group velocity and the phase velocity are effectively solved by the phase slowness direction and the ray direction so as to solve the equation of the program.
In a specific embodiment of the present invention, step S3 specifically includes the following steps:
t according to formula (13) 0 Is still unknown and can be calculated as
Figure BDA0002784893810000092
wherein ,
Figure BDA0002784893810000093
to calculate the coordinates in the domain +.>
Figure BDA0002784893810000094
For the coordinates of the source position +.>
Figure BDA0002784893810000095
To be along the radial direction
Figure BDA0002784893810000096
Group velocity of (c).
Wherein the group velocity
Figure BDA0002784893810000101
Expressed by the following formula:
Figure BDA0002784893810000102
when the reciprocal of the phase velocity of the qP wave, qSV wave, and qSH wave are expressed as follows:
Figure BDA0002784893810000103
Figure BDA0002784893810000104
wherein the derivatives of M and N are represented by the following formula:
Figure BDA0002784893810000105
Figure BDA0002784893810000106
phase slowness direction
Figure BDA0002784893810000107
And (c) angle->
Figure BDA0002784893810000108
Are implicitly defined, they depend on the solution t=t 0 U. In the following description we will give a numerical program to calculate the group velocity and phase velocity along a given ray direction.
In the formula (13), T 0 Derivatives (T) about x, y, z 0x ,T 0y ,T 0z ) And is also an unknown. Unlike isotropic media, high-order finite difference methods cannot be used for computation because high-frequency noise may be present when the wavefront is not a standard sphere/sphere. Thus, we present an analytical expression for calculating these three variables, expressed in a uniform anisotropic medium as:
Figure BDA0002784893810000109
/>
wherein ,v0 As a function of the phase velocity at the source point, (p) 0 ,q 0 ,r 0 )＝(T 0x ,T 0y ,T 0z ) Is T 0 Derivatives with respect to x, y, z. Extracting p sequentially from square root of equation (18) 0 ,q 0 ,r 0 And performing algebraic operations, wherein the result obtained after the operations is expressed as follows:
Figure BDA00027848938100001010
Figure BDA0002784893810000111
Figure BDA0002784893810000112
for a specific ray direction, θ and θ may be determined using a solution method for group velocity demonstrated later
Figure BDA0002784893810000113
According to theta and->
Figure BDA0002784893810000114
We can determine the value ranges of these three variables (p 0 ,q 0 ,r 0 ) Is a symbol of (c).
Wherein, in calculating T 0 When we need to calculate the group velocity along a given ray direction. However, the group velocity is a function of the slowness vector, not the direction of the rays (19). If the slowness vector of the line direction can be calculated, the group velocity can be easily calculated. Earlier work studied how to calculate slowness vectors for a particular ray direction.
If the phase slowness direction is given
Figure BDA0002784893810000115
And phase velocity v m The slowness vector may be expressed as:
Figure BDA0002784893810000116
slowness vector
Figure BDA0002784893810000117
And group velocity vector>
Figure BDA0002784893810000118
The following relationship should be satisfied:
Figure BDA0002784893810000119
phase slowness direction
Figure BDA00027848938100001110
And ray direction->
Figure BDA00027848938100001111
Expressed by the following formula:
Figure BDA00027848938100001112
dividing both sides of (28) by
Figure BDA00027848938100001113
The following formula is obtained:
Figure BDA00027848938100001114
equation (25) will be used to calculate the phase velocity and group velocity along a given ray direction. Let the ray direction be
Figure BDA00027848938100001115
The phase slowness direction is +.>
Figure BDA00027848938100001116
Formula (25) is rewritten as:
Figure BDA0002784893810000121
the equation may be solved using pseudo-position methods or equal-division methods, etc. Once in the radial direction
Figure BDA0002784893810000122
Calculate the phase slowness direction +.>
Figure BDA0002784893810000123
Group velocity U can be calculated from equations (19) and (3), respectively m Sum phase velocity v m . The method is applicable to three modes in any anisotropic medium.
Preferably, as seen in fig. 4, fig. 4 is a schematic flow chart of the Lax-Friedrichs scanning provided by the present invention, and step S4 includes steps S41 to S42, where:
in step S41, a corresponding first-order Lax-Friedrichs format of Hamiltonian form is established;
in step S42, according to the third-order windward bias WENO (third-order weighted intrinsic non-oscillation format) approximation value corresponding to each grid point, replacing the travel time value on the grid point, and establishing a static hamilton-jacobian equation based on the third-order Lax-Friedrichs format of the WENO;
in step S43, the static hamiltonian-jacobian equation in the WENO-based third-order Lax-Friedrichs format is written in a compact form, and travel time values at grid points are updated according to the compact form.
Therefore, by adopting the Lax-Friedrichs finite difference format, iteration convergence can be ensured no matter whether the Hamiltonian is convex or concave. And (3) carrying out numerical calculation on the basis function by using a fast sweep frequency format based on a third-order weighted intrinsic non-oscillation format (WENO), so as to achieve the purpose of fast scanning, wherein the iteration times are irrelevant to the size of the grid, and the numerical solution converges to the expected weak solution as the size of the grid approaches 0.
In a specific embodiment of the present invention, formula (13) is represented as follows:
Figure BDA0002784893810000124
wherein ,αx ,α y and αz Is the viscosity coefficient.
In the formula (27), u i-1,j,k ,u i+1,j,k ,u i,j-1,k ,u i,j+1,k ,u i,j,k-1 and ui,j,k+1 Replaced by
Figure BDA0002784893810000131
wherein ,
Figure BDA0002784893810000132
and />
Figure BDA0002784893810000133
Is u x Is similar to the third-order windward bias WENO; />
Figure BDA0002784893810000134
and />
Figure BDA0002784893810000135
Is u y Is similar to the third-order windward bias WENO; />
Figure BDA0002784893810000136
and />
Figure BDA0002784893810000137
Is u z Is approximated by a third order windward bias, WENO, wherein:
Figure BDA0002784893810000138
when (when)
Figure BDA0002784893810000139
and
Figure BDA00027848938100001310
When (when)
Figure BDA00027848938100001311
Similarly, we can define
Figure BDA00027848938100001312
and />
Figure BDA00027848938100001313
Is a third order WENO approximation. Epsilon is a small positive number to avoid dividing by 0.
Then, the static hamiltonian-jacobian equation based on the WENO third-order Lax-Friedrichs format can be written in a compact form:
Figure BDA00027848938100001314
Figure BDA0002784893810000141
wherein ,
Figure BDA0002784893810000142
for the numerical solution of the grid points (i, j, k) to be updated, < >>
Figure BDA0002784893810000143
Is the old numerical solution of the same point.
Preferably, as seen in conjunction with fig. 5, fig. 5 is a schematic flow chart of updating a grid provided by the present invention, and step S43 includes steps S431 to S434, wherein:
in step S431, initializing and assigning values to grid points in the first-order Lax-Friedrichs format;
in step S432, the compact form is solved iteratively with Gauss-Seidel (gaussian-Seidel iteration) of 8 alternating direction sweeps, and grid updating is performed;
in step S433, a plurality of maximum and minimum conditions are set on six planes of the inner boundary and the outer boundary, and a numerical solution of the off-domain point is determined by using a linear extrapolation method;
in step S434, it is determined whether the L1-norm difference of the solution between two successive iterations satisfies a preset accuracy condition, and if so, the iteration is terminated.
Therefore, the Gauss-Seidel iterative solution of 8 times of alternate direction sweep frequency is utilized to achieve the purpose of rapid scanning solution, and meanwhile, the numerical solution of the out-of-domain points is accurately calculated to ensure the high efficiency of grid updating.
In an embodiment of the invention, the hamilton-jacobian equation (13) is solved using a Lax-Friedrichs format (19) based on third-order WENO. The procedure can be summarized as follows (Zhang et al, 2006; luo et al, 2012):
algorithm diagram of three-order Lax-Friedrichs scanning format of anisotropic factor-path function equation
Initialization, assigning values on grid points at a distance of less than or equal to 2h from the source point, because the third order WENO approximation requires at least three grid points. These values will be fixed during the iteration and all other grid points are assigned larger positive values.
Gauss-Seidel iteration with 8 alternating direction sweeps solved the discretization equation (33):
Figure BDA0002784893810000151
at each grid point (i, j, k), an update is made according to the method described above.
Termination if the L1-norm difference of the solution between two successive iterations is less than the specified accuracy requirement, the iteration is terminated.
In performing the third-order Lax-Friedrichs scanning scheme, linear extrapolation is employed near and at boundary points. On the six planes of the inner boundary, the following conditions are set:
Figure BDA0002784893810000152
then, on the six planes of the outer boundary, the following conditions are applied:
Figure BDA0002784893810000153
Figure BDA0002784893810000161
and calculating a numerical solution of the point outside the calculation domain by using an extrapolation method, a maximum value method and a minimum value method. In the gaussian-saidel iteration of the fast sweep method, the solution equations must be solved (33) efficiently and the unknown factors calculated.
In a specific embodiment of the present invention, four test cases are given for illustration:
the first example is a test of a sinusoidal anisotropy model. The calculated domain is [ -0.5,0.5] × [0.0,0.5] km rectangular domain, the point source is located at x=0.0 km, y=0.0 km. The anisotropic parameters (a 11 (x), a13 (x), a33 (x), a44 (x), a66 (x), θ0 (x)) of the model are listed in table 1, table 1 being as follows:
TABLE 1
Anisotropic parameters Generating a function of an anisotropy parameter
a 11 5.2×{1.5+0.5sin(0.4π(z+0.01))sin(3π(x+0.05))}
a 13 0.98×{1.5+0.5sin(0.4π(z+0.01))sin(3π(x+0.05))}
a 33 4.0×{1.5+0.5sin(0.4π(z+0.01))sin(3π(x+0.05))}
a 44 1.0×{1.5+0.5sin(0.4π(z+0.01))sin(3π(x+0.05))}
a 66 1.0×{1.5+0.5sin(0.4π(z+0.01))sin(3π(x+0.05))}
θ 0 0 0
The reference solution is computed on a very fine grid 1601 x 801 using the SPM model. The model was sampled on a 101 x 51 grid for the Lax-Friedrichs scanning method. Table 2 lists the L1 norm errors and convergence order for four different grid sizes. The convergence order is greater than the first order, verifying the effectiveness of the Lax-Friedrichs scan method to solve for source singularities, table 2 is shown below:
TABLE 2
Figure BDA0002784893810000171
As can be seen from fig. 6, 7 and 8, fig. 6 is a schematic diagram of the travel time of the qP wave provided by the present invention, fig. 7 is a schematic diagram of the travel time of the qSV wave provided by the present invention, and fig. 8 is a schematic diagram of the travel time of the qSH wave provided by the present invention. In fig. 6, 7, 8 we compare travel times of the three modes generated by the SPM and LAX methods.
The second example tested the Hess VTI model with a rectangular area of [ -16.0,16.0] × [0.0,12.0] km. The point source is located at x=0.0 km, and y=0.0 km. There are five elastic models (a 11 (x), a13 (x), a33 (x), a44 (x), a66 (x)), and θ0 (x) =00 (fig. 4). Steep faults and high velocity salts are present in the model. The model is resampled using a grid 401 x 151. As can be seen from fig. 9 to 11, fig. 9 shows the intention of the travel time of the qP wave provided by the present invention, fig. 10 shows the intention of the travel time of the qSV provided by the present invention, and fig. 11 shows the intention of the travel time of the qSH provided by the present invention, and it can be seen from the figure that the Lax-Friedrichs sweep method can well calculate the travel time table of the qP, qSV and qSH waves.
A third example is a test performed on a anticline formation with a high inclination on both sides of the formation. This model is referred to as part of the BP TTI model. Referring to fig. 12, fig. 12 is a schematic diagram of parameters of anisotropy provided by the present invention, and it can be seen from the figure that (a 11 (x), a13 (x), a33 (x), a44 (x), a66 (x)) and an inclination angle θ0 (x), wherein the calculated area is a [ -5.0,5.0] × [0.0,10.0] km rectangular area, and the source point is located at x=0.0 km, and y=0.0 km. The model uses 251 x 251 grid sampling. Referring to fig. 13 to 15, fig. 13 is a schematic diagram of the travel time of the qP wave under the scan of the Lax-Friedrichs provided by the present invention, fig. 14 is a schematic diagram of the travel time of the qSV wave under the scan of the Lax-Friedrichs provided by the present invention, and fig. 15 is a schematic diagram of the travel time of the qSH wave under the scan of the Lax-Friedrichs provided by the present invention, and it can be seen that the process of the Lax-Friedrichs scan successfully generates the travel time of three wave modes.
A fourth example is to test the proposed method in a three-layer anisotropic medium. There is an inclined interface in the model. The anisotropic parameters (a 11 (x), a13 (x), a33 (x), a44 (x), a66 (x)), tilt angle θ0 (x) and azimuth angle Φ0 (x) are listed in table 3, table 3 being as follows:
TABLE 3 Table 3
Figure BDA0002784893810000181
As can be seen from fig. 16 to 18, fig. 16 is a three-dimensional travel time contour map of the qP wave provided by the present invention, fig. 17 is a three-dimensional travel time contour map of the qSV wave provided by the present invention, and fig. 18 is a three-dimensional travel time contour map of the qSH wave provided by the present invention. From the figure, the three-dimensional expansion of the two-dimensional model is used to calculate the travel time and ray path of the elastic wave, the calculation domain is [ -2.5,2.5] × [ -1.0,1.0] × [0.0,5.0] km rectangular domain, the source point is located at x=0.0 km, z=0.0 km, and the grid of the Lax-Friedrichs scanning method is 101×41×101. Three-dimensional travel time contours of qP, qSV and qSH waves were plotted on an anisotropic model (a 11 (x) model) and a three-dimensional travel time cube (fig. 16-18), respectively.
Example 2
The embodiment of the invention provides a device for scanning travel time of seismic waves by using Lax Friedrichs, and referring to FIG. 19, FIG. 19 is a schematic structural diagram of the device for scanning travel time of seismic waves by using Lax Friedrichs, where the device 1900 includes:
a modeling unit 1901 for creating a program function equation of the seismic wave in the anisotropic TI medium by introducing a slowness vector;
a conversion unit 1902, configured to decompose a solution of the equation into a product of two factors, and bring the product of the two factors into the equation, and convert the product into a corresponding hamiltonian form;
a solving unit 1903 for solving the group velocity and the phase velocity according to the slowness vector;
and a scanning unit 1904, configured to determine travel time of the seismic wave by using the Lax-Friedrichs scanning format values and solving a equation.
Example 3
The embodiment of the invention provides a computer readable storage medium, wherein a computer program is stored on the computer readable storage medium, and when the computer program is executed by a processor, the method for scanning the travel time of the seismic waves by the Lax Friedrichs is realized.
The invention discloses a method for scanning seismic wave travel time by Lax Friedrichs, which can calculate the first arrival travel time of three wave types, namely qP-, qSV-, qSH-waves, and adopts a factor decomposition method to treat the upwind singularity of a finite difference program function solver, and decompose the solution of a general program function equation into the product of two factors, wherein one factor is the travel time of a uniform medium, the factor captures upwind singularities, and the other factor (basic function) is one necessary modification near a source point. Further, a Lax-Friedrichs finite difference format is adopted, iteration convergence can be guaranteed no matter whether the Hamiltonian is convex or concave, a rapid sweep frequency format based on a third-order weighted intrinsic non-oscillation format (WENO) is adopted to carry out numerical calculation on the basis function, the iteration times are irrelevant to the size of the grid, and as the size of the grid approaches 0, the numerical solution converges on the expected weak solution. Therefore, the invention provides Hamiltonian form of the factor anisotropic equation of qP-, qSV-and qSH-waves in the TTI medium by utilizing an effective factor decomposition method, designs a rapid scanning method combining Lax-Friedrichs finite difference format, and ensures the high efficiency and accuracy of calculating travel time.
According to the technical scheme, a factor solving and quick scanning mode is adopted, and the very accurate first arrival time can be calculated through a Hamiltonian form of a three-order weighted basic non-oscillation approximation and a Lax-Friedrichs format discretization anisotropic process function equation and then a quick scanning method for iteratively solving a discrete system. The anisotropic equation does not have weak anisotropy assumption on phase and group velocity, so the method is suitable for any anisotropic oblique transverse isotropic (TTI) medium, has strong flexibility, high generalization degree and correctness and effectiveness, ensures accuracy and rapidity of travel time calculation, and is convenient for effective geological investigation.
The present invention is not limited to the above-mentioned embodiments, and any changes or substitutions that can be easily understood by those skilled in the art within the technical scope of the present invention are intended to be included in the scope of the present invention.

Claims (5)

1. A method for LaxFriedrichs to scan travel time of seismic waves, comprising:
by introducing a slowness vector, establishing a program function equation of the seismic waves in the anisotropic TI medium;
decomposing the solution of the equation into products of two factors, and bringing the products of the two factors into the equation to be converted into corresponding Hamiltonian forms;
solving group velocity and phase velocity according to the slowness vector;
solving the equation of the equation by utilizing Lax-Friedrichs scanning format values, and determining the travel time of the seismic waves;
the decomposing the solution of the equation into products of two factors, and bringing the products of the two factors into the equation, the converting into the corresponding Hamiltonian form includes:
decomposing the solution of the equation of the program into a first factor and a second factor, the first factor and the second factor being expressed as:
Figure QLYQS_1
wherein ,
Figure QLYQS_3
to calculate the coordinates in the domain +.>
Figure QLYQS_5
For said first factor,/o>
Figure QLYQS_8
As a result of the said second factor,
Figure QLYQS_4
is to predetermine the travel time field coordinates of the capture source singularities,/->
Figure QLYQS_6
For an unknown factor that is smoothed near the source point,
Figure QLYQS_7
for the phase slowness model, +.>
Figure QLYQS_9
For the phase slowness model +.>
Figure QLYQS_2
A corresponding location factor;
bringing the first factor and the second factor into the equation of the program function and based on the unknown factor
Figure QLYQS_10
And +.about.of the predetermined capture source singularities>
Figure QLYQS_11
Performing formula conversion on the derivatives of the coordinates x, y and z, and determining the corresponding Hamiltonian form as follows:
Figure QLYQS_12
Figure QLYQS_13
, />
Figure QLYQS_14
wherein ,
Figure QLYQS_16
the unknown factors +.>
Figure QLYQS_20
Derivatives with respect to coordinates x, y, z, < >>
Figure QLYQS_23
Is a travel time field of predetermined capture source singularities, +.>
Figure QLYQS_17
Is Hamiltonian, japan>
Figure QLYQS_19
For the corresponding Hamiltonian derivative function, +.>
Figure QLYQS_22
For bringing in->
Figure QLYQS_24
Posterior phase slowness model, +.>
Figure QLYQS_15
、/>
Figure QLYQS_18
And R is said first factor +.>
Figure QLYQS_21
Derivative with respect to coordinates x, y, z。
2. The method of claim 1, wherein said solving for group velocity and phase velocity based on said slowness vector comprises:
setting a phase slowness direction and a ray direction;
and determining the group velocity and the phase velocity according to the relation between the slowness vector and the group velocity vector.
3. A method of LaxFriedrichs scanning seismic travel time according to claim 2, wherein in the corresponding hamiltonian form the determination of the derivatives of the travel time with respect to coordinate axes x, y, z comprises:
determining an analytical expression of the derivative of the travel time with respect to the coordinate axes x, y, z, expressed as:
Figure QLYQS_25
, />
Figure QLYQS_26
.
wherein ,
Figure QLYQS_27
as a function of the phase velocity at the source point +.>
Figure QLYQS_28
，/>
Figure QLYQS_29
Derivatives of the travel time with respect to x, y, z, respectively;
sequentially extracting derivatives of the travel time with respect to x, y, z from square roots of the analytical expression
Figure QLYQS_30
And determining the travel time by algebraic operationDerivative of x, y, z->
Figure QLYQS_31
Expressed as the corresponding absolute value of:
Figure QLYQS_32
, />
Figure QLYQS_33
Figure QLYQS_34
, />
Figure QLYQS_35
Figure QLYQS_36
Figure QLYQS_37
wherein ,
Figure QLYQS_38
is the tilt angle and azimuth angle of the phase slowness direction.
4. A method of scanning travel time of seismic waves according to claim 3, wherein said solving said equation of the equation using the Lax-Friedrichs scan format values comprises:
establishing a first-order Lax-Friedrichs format of the corresponding Hamiltonian form;
according to the approximate value of the third-order windward bias WENO corresponding to each grid point, replacing the travel time value on the grid point, and establishing a static Hamiltonian-Jacobian equation based on the third-order Lax-Friedrichs format of the WENO;
writing the static Hamiltonian-Jacobian equation in the WENO-based third-order Lax-Friedrichs format into a compact form, and updating the travel time value on the grid point according to the compact form.
5. The method of claim 4, wherein updating travel time values at the grid points based on the compact form comprises:
initializing and assigning grid points in the first-order Lax-Friedrichs format;
solving the compact form by using Gauss-Seidel of 8 alternate direction sweep frequency iterations, and updating the grid;
setting a plurality of corresponding maximum and minimum value conditions on six planes of an inner boundary and an outer boundary, and determining a numerical solution of an out-of-domain point by adopting a linear extrapolation method;
and judging whether the L1 norm difference of the solution between two continuous iterations meets a preset precision condition, and if so, terminating the iteration.
CN202011294721.XA 2020-11-18 2020-11-18 Method for scanning travel time of seismic waves by using Lax Friedrichs Active CN112505765B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011294721.XA CN112505765B (en) 2020-11-18 2020-11-18 Method for scanning travel time of seismic waves by using Lax Friedrichs

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011294721.XA CN112505765B (en) 2020-11-18 2020-11-18 Method for scanning travel time of seismic waves by using Lax Friedrichs

Publications (2)

Publication Number Publication Date
CN112505765A CN112505765A (en) 2021-03-16
CN112505765B true CN112505765B (en) 2023-05-09

Family

ID=74956908

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011294721.XA Active CN112505765B (en) 2020-11-18 2020-11-18 Method for scanning travel time of seismic waves by using Lax Friedrichs

Country Status (1)

Country Link
CN (1) CN112505765B (en)

Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5062086A (en) * 1990-08-27 1991-10-29 Conoco Inc. Calculation of raypaths and wavepaths from traveltime tables for the tomographic estimation of transmission velocities
US5394325A (en) * 1993-04-07 1995-02-28 Exxon Production Research Company Robust, efficient three-dimensional finite-difference traveltime calculations
AU2013206119A1 (en) * 2007-10-30 2013-06-20 University Of Utah Research Foundation Fast iterative method for processing Hamilton-Jacobi equations
CN105388518A (en) * 2014-09-04 2016-03-09 中国石油化工股份有限公司 Centroid frequency and spectral ratio integrated borehole seismic quality factor inversion method
CN106154325A (en) * 2016-06-20 2016-11-23 吉林大学 Relief surface based on ray theory combination source wavefield orientation method
CN107966729A (en) * 2016-10-19 2018-04-27 中国石油化工股份有限公司 A kind of three-dimensional TTI media ray-tracing procedure and system
CN108363101A (en) * 2018-02-02 2018-08-03 西安石油大学 A kind of inclined shaft crosshole seismic Gaussian beam pre-stack depth migration imaging method
CN109784724A (en) * 2019-01-15 2019-05-21 常州大学 A kind of subsea production system equipment fault diagnosis method based on Bayesian network
CN110568497A (en) * 2019-09-26 2019-12-13 核工业北京地质研究院 Accurate solving method for seismic first-motion wave travel time under complex medium condition
CN110674893A (en) * 2019-10-30 2020-01-10 江苏方天电力技术有限公司 Self-adaptive correction method for diagnosis experience in rotary machine fault diagnosis knowledge base

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6748330B2 (en) * 2002-04-10 2004-06-08 Schlumberger Technology Corporation Method and apparatus for anisotropic vector plane wave decomposition for 3D vertical seismic profile data
US20160202375A1 (en) * 2013-09-20 2016-07-14 Westerngeco Llc Eikonal Solver for Quasi P-Waves in Anisotropic Media

Patent Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5062086A (en) * 1990-08-27 1991-10-29 Conoco Inc. Calculation of raypaths and wavepaths from traveltime tables for the tomographic estimation of transmission velocities
US5394325A (en) * 1993-04-07 1995-02-28 Exxon Production Research Company Robust, efficient three-dimensional finite-difference traveltime calculations
AU2013206119A1 (en) * 2007-10-30 2013-06-20 University Of Utah Research Foundation Fast iterative method for processing Hamilton-Jacobi equations
CN105388518A (en) * 2014-09-04 2016-03-09 中国石油化工股份有限公司 Centroid frequency and spectral ratio integrated borehole seismic quality factor inversion method
CN106154325A (en) * 2016-06-20 2016-11-23 吉林大学 Relief surface based on ray theory combination source wavefield orientation method
CN107966729A (en) * 2016-10-19 2018-04-27 中国石油化工股份有限公司 A kind of three-dimensional TTI media ray-tracing procedure and system
CN108363101A (en) * 2018-02-02 2018-08-03 西安石油大学 A kind of inclined shaft crosshole seismic Gaussian beam pre-stack depth migration imaging method
CN109784724A (en) * 2019-01-15 2019-05-21 常州大学 A kind of subsea production system equipment fault diagnosis method based on Bayesian network
CN110568497A (en) * 2019-09-26 2019-12-13 核工业北京地质研究院 Accurate solving method for seismic first-motion wave travel time under complex medium condition
CN110674893A (en) * 2019-10-30 2020-01-10 江苏方天电力技术有限公司 Self-adaptive correction method for diagnosis experience in rotary machine fault diagnosis knowledge base

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Lax Friedrichs;何雷宇等;《地球物理学进展》;第32卷(第03期);第1140-1148页 *
Lax–Friedrichs sweeping scheme for static Hamilton–Jacobi equations;Kao等;《Journal of Computational physics》;第196卷(第1期);第367-391页 *
一种改进的利用程函方程计算旅行时的方法;李文杰 等;《石油地球物理勘探》(第05期);第589-594+487页 *
基于因式分解形式程函方程的地震走时层析成像方法研究;胡秋萍;《中国优秀硕士学位论文全文数据库基础科学辑》(第01期);第A011-799页 *

Also Published As

Publication number Publication date
CN112505765A (en) 2021-03-16

Similar Documents

Publication Publication Date Title
Rawlinson et al. Seismic ray tracing and wavefront tracking in laterally heterogeneous media
Leung et al. An adjoint state method for three-dimensional transmission traveltime tomography using first-arrivals
CN106199742B (en) A kind of dimension of Frequency-domain AEM 2.5 band landform inversion method
AU2008239658B2 (en) Inverse-vector method for smoothing dips and azimuths
Waheed et al. A fast sweeping algorithm for accurate solution of the tilted transversely isotropic eikonal equation using factorization
Cameron et al. Seismic velocity estimation from time migration
JP7142968B2 (en) FULL WAVEFORM INVERSION METHOD, APPARATUS AND ELECTRONICS
CN110058307B (en) Full waveform inversion method based on fast quasi-Newton method
CN111239819B (en) Direct envelope inversion method with polarity based on seismic channel attribute analysis
Qian et al. A level set-based Eulerian approach for anisotropic wave propagation
Le Bouteiller et al. A discontinuous Galerkin fast-sweeping eikonal solver for fast and accurate traveltime computation in 3D tilted anisotropic media
CN111290019B (en) L-BFGS initial matrix solving method applied to least square reverse time migration
Huang et al. Local algorithm for computing complex travel time based on the complex eikonal equation
CN111580163B (en) Full waveform inversion method and system based on non-monotonic search technology
Fang et al. A hybrid approach to solve the high-frequency Helmholtz equation with source singularity in smooth heterogeneous media
WO2021116800A1 (en) System and method for using a neural network to formulate an optimization problem
Hou et al. Approximate Gauss-Newton iteration for full-waveform inversion
CN112505765B (en) Method for scanning travel time of seismic waves by using Lax Friedrichs
Zhou et al. An iterative factored topography-dependent eikonal solver for anisotropic media
Wang et al. A fast sweeping method for P-wave traveltimes in attenuating VTI media with an irregular surface
Leung et al. A level-set adjoint-state method for transmission traveltime tomography in irregular domains
Martinelli An application of the level set method to underwater acoustic propagation
Zhou et al. A topography-dependent eikonal solver for accurate and efficient computation of traveltimes and their derivatives in 3D heterogeneous media
Padhi et al. Accurate quasi-p traveltimes in 3d transversely isotropic media using a high-order fast-sweeping-based eikonal solver
Zhang et al. Optimal transport with a new preprocessing for deep-learning full waveform inversion

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