CN112505765A - Method for scanning seismic wave travel time by Lax Friedrichs - Google Patents

Abstract

The invention relates to a Lax Friedrichs scanning seismic wave travel time method. The method comprises: establishing a seismic wave function equation in anisotropic TI medium by introducing a slowness vector; decomposing the solution of the function equation into two factors product, and bring the product of the two factors into the functional equation, and convert it into the corresponding Hamiltonian form; according to the slowness vector, solve the group velocity and phase velocity; use the Lax-Friedrichs scan format value to solve the functional equation , to determine the travel time of seismic waves. The method, device and storage medium for the Lax Friedrichs scanning seismic wave travel time provided by the present invention adopt the method of factoring and fast scanning, which ensures the accuracy and fastness of travel time calculation, so as to carry out effective geological survey.

Description

Lax Friedrichs' method for scanning the travel time of seismic waves

technical field

The invention relates to the technical field of geological exploration, in particular to a method for scanning the travel time of seismic waves by Lax Friedrichs.

Background technique

Earthquake travel time refers to the time it takes for seismic waves to travel from the epicenter to the observation point, which is called the travel time of seismic waves, also known as travel time. Usually, we record the travel time of the seismic wave through the ground sensor, and invert the real position of the underground geological body, and then invert the physical properties of the underground rock. Therefore, the calculation of the travel time of seismic waves can truly reflect the geological situation, which is a key technology in the field of geological exploration technology.

Seismic waves have very different kinematics and dynamics when propagating in isotropic and anisotropic media. There are only longitudinal and shear waves in isotropic media, which involves three wave modes in anisotropic media, including one quasi-longitudinal wave (qP) and two quasi-shear waves (qSV and qSH). Each wave mode has its own wave mode and polarization direction. In the prior art, the finite-difference-based equation solver is one of the most stable and efficient travel time calculation methods. However, there is a lack of methods for finite-difference equation solvers using anisotropic media. At the same time, the singularity of the source This property causes a large error at the source point, and the error spreads over the entire area, resulting in lower accuracy.

In fact, more and more commercial seismic processing software (such as Landmark, CGG, Omega, etc.) takes anisotropy into account. In addition, the equations governing weak anisotropy are much simpler than those governing strong anisotropy, so the weak anisotropy assumption is widely used in various techniques of seismic data processing. This assumption makes various techniques of seismic data processing easier to implement. However, in the prior art, the method for calculating travel time by using weak anisotropy pairs often lacks accuracy and rapidity. How to propose a method for calculating travel time quickly and accurately is an urgent problem to be solved.

SUMMARY OF THE INVENTION

In view of this, it is necessary to provide a Lax Friedrichs method for scanning the travel time of seismic waves to solve the problem of how to propose a fast and accurate method for travel time calculation.

The present invention provides a Lax Friedrichs method for scanning the travel time of seismic waves, comprising:

By introducing the slowness vector, the function equation of seismic wave in anisotropic TI medium is established;

Decompose the solution of the functional equation into a product of two factors, and bring the product of the two factors into the functional equation, and convert it into a corresponding Hamiltonian form;

According to the slowness vector, solve the group velocity and the phase velocity;

The travel time of the seismic wave is determined by solving the functional equation using the Lax-Friedrichs scan format numerical values.

Further, the solution of the functional equation is decomposed into a first factor and a second factor, and the first factor and the second factor are expressed as:

Among them, x, y, z are the coordinates in the computational domain, T(x, y, z) is the first factor, s(x, y, z) is the second factor, T ₀ (x, y) ,z) is pre-determined to capture the source singularity, u(x,y,z) is an unknown factor smoothed around the source point, s ₀ (x,y,z) is the phase slowness model, α(x,y , z) is the position factor corresponding to the phase slowness model s ₀ (x, y, z);

The first factor and the second factor are brought into the functional equation, and the coordinates _x , y, The derivative of z is converted into a formula, and the corresponding Hamiltonian form is determined as:

f(x,y,z)=s _m (P,Q,R),(m=1,2,3)

where (u _x , 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, u _x , u _y , u _z ) is the Hamiltonian operator, f(x, y, z) is the corresponding Hamiltonian derivative function, s _m (P, Q, R) is the phase slowness model brought into P, Q, R, P, Q and R are derivatives of the first factor T(x, y, z) with respect to the coordinates x, y, z.

Further, calculating the group velocity and the phase velocity according to the slowness vector includes:

Set the phase slowness direction and ray direction;

The group velocity and the phase velocity are determined from the relationship 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, and z includes:

An analytical expression for determining the derivative of the travel time with respect to the coordinate axes x, y, and z is expressed as:

Wherein, v ₀ is the phase velocity function at the source point, (p ₀ , q ₀ , r ₀ )=(T _0x , T _0y , T _0z ), which is the derivative of the travel time with respect to x, y, z;

The derivatives p ₀ , q ₀ , r ₀ of the travel time with respect to x, y, z are sequentially extracted from the square root of the analytical expression, and through algebraic operations, the derivatives of the travel time with respect to x, y, z are determined The corresponding absolute values of p ₀ , q ₀ , r ₀ are expressed as:

为相位慢度方向的的倾角和方位角。in,

are the inclination and azimuth of the phase slowness direction.

Further, described utilizing Lax-Friedrichs scanning format numerical value, solving described functional equation comprises:

establishing a first-order Lax-Friedrichs scheme of the corresponding Hamiltonian form;

Replace the travel time value on the grid point according to the approximation value of the third-order upwind bias WENO corresponding to each grid point;

The static Hamiltonian-Jacobi equation based on the WENO third-order Lax-Friedrichs format is written in a compact form, and the travel time value on the grid point is updated according to the compact form.

Further, the updating the travel time value on the grid point according to the compact form includes:

initializing and assigning grid points in the first-order Lax-Friedrichs format;

Solve the compact form with 8 Gauss-Seidel iterations swept in alternating directions, and perform mesh updates;

On the six planes of the inner boundary and the outer boundary, set the corresponding maximum and minimum value conditions, and use the linear extrapolation method to determine the numerical solution of the point outside the domain;

It is judged whether the L1 norm difference of the solution between two consecutive iterations satisfies the preset accuracy condition, and if so, the iteration is terminated.

Compared with the prior art, the beneficial effects of the present invention include: the first arrival travel time of three wave modes, namely qP-, qSV-, and qSH-waves, can be calculated, and the factorization method is used to process the finite difference distance function solver. Headwind singularity, which decomposes the solution of the general equation of function into the product of two factors, one of which is the traveltime of the homogeneous medium, this factor captures the headwind singularity, and the other (the basis function) is a factor near the source point Modifications are necessary. The Lax-Friedrichs finite-difference scheme is further adopted, which can guarantee the iterative convergence regardless of whether the Hamiltonian function is convex or concave. The basis function is numerically calculated by the fast frequency sweep scheme based on the third-order weighted intrinsically non-oscillating scheme (WENO), and iteratively The degree is independent of the grid size, and as the grid size approaches 0, the numerical solution converges to the desired weak solution. To sum up, the method, device and storage medium for the LaxFriedrichs scanning seismic wave travel time provided by the present invention adopts the method of factoring and fast scanning, which ensures the accuracy and speed of travel time calculation, so as to carry out effective geological survey.

Description of drawings

Fig. 1 is the schematic flow chart of the method for the Lax Friedrichs scanning seismic wave travel time provided by the present invention;

Fig. 2 is the schematic flow sheet of factorization provided by the present invention;

3 is a schematic flow chart for solving group velocity and phase velocity provided by the present invention;

4 is a schematic flow chart of Lax-Friedrichs scanning provided by the present invention;

5 is a schematic flowchart of updating grid provided by the present invention;

6 is a schematic diagram of the travel time of the qP wave provided by the present invention;

Fig. 7 is the travel time schematic diagram of qSV wave provided by the present invention;

8 is a schematic diagram of the travel time of the qSH wave provided by the present invention;

9 is a schematic diagram of a travel time table of a qP wave provided by the present invention;

Fig. 10 is the travel time table schematic diagram of qSV provided by the present invention;

Fig. 11 is the travel time representation of qSH provided by the present invention;

12 is a schematic diagram of anisotropy parameters provided by the present invention;

13 is a schematic diagram of the travel time of the qP wave under the Lax-Friedrichs scan provided by the present invention;

14 is a schematic diagram of the travel time of the qSV wave under the Lax-Friedrichs scan provided by the present invention;

15 is a schematic diagram of the travel time of the qSH wave under the Lax-Friedrichs scan provided by the present invention;

16 is a three-dimensional travel time contour diagram of a qP wave provided by the present invention;

17 is a three-dimensional travel time contour map of a qSV wave provided by the present invention;

18 is a three-dimensional travel time contour map of a qSH wave provided by the present invention;

FIG. 19 is a schematic structural diagram of the Lax Friedrichs device for scanning the travel time of seismic waves provided by the present invention.

Detailed ways

The preferred embodiments of the present invention are specifically described below with reference to the accompanying drawings, wherein the accompanying drawings constitute a part of the present application, and together with the embodiments of the present invention, are used to explain the principles of the present invention, but are not used to limit the scope of the present invention.

Example 1

An embodiment of the present invention provides a method for scanning the travel time of seismic waves by Lax Friedrichs. Referring to FIG. 1 , FIG. 1 is a schematic flowchart of the method for scanning travel time of seismic waves by Lax Friedrichs provided by the present invention. The method comprises steps S1 to S4, wherein:

In step S1, by introducing the slowness vector, the function equation of the seismic wave in the anisotropic TI medium is established;

In step S2, the solution of the function equation is decomposed into the product of two factors, and the product of the two factors is brought into the function equation, and converted into the corresponding Hamiltonian formula;

In step S3, according to the slowness vector, the group velocity and the phase velocity are solved;

In step S4, the Lax-Friedrichs scanning format numerical value is used to solve the functional equation.

Preferably, step S1 specifically includes the following steps:

According to the Christoffel equation, the slowness vector is expressed as:

为波前法向向量，v为地震波的相速度，地震波包括qP波、qSV波和qSH波；in,

is the wavefront normal vector, v is the phase velocity of the seismic wave, and the seismic wave includes qP wave, qSV wave and qSH wave;

According to the slowness vector, the established program function equation is expressed as the following formula:

为旅行时间的倒数且与慢度向量相等，v_m(m＝1,2,3)为qP波、qSV波和qSH波的相速度，分别表示为下式：where T is the travel time,

is the inverse of the travel time and is equal to the slowness vector, v _m (m=1,2,3) is the phase velocity of the qP wave, qSV wave and qSH wave, respectively expressed as the following equations:

是由相慢度方向和各向异性TI介质的对称轴形成的，表示为下式：where v _1,2 are the phase velocities of the qP wave and qSV wave, v ₃ is the phase velocity of the qSH wave, a ₄₄ and a ₆₆ are the elastic moduli of the anisotropic TI medium,

is formed by the phase slowness direction and the symmetry axis of the anisotropic TI medium and is expressed as:

为各向异性TI介质的对称轴的倾斜角和方位角，θ和

分别是相位慢度方向的切斜角和方位角，表示为下式：in,

are the inclination and azimuth angles of the symmetry axis of the anisotropic TI medium, θ and

are the shear angle and the azimuth angle in the phase slowness direction, respectively, and are expressed as:

where P ₌ Tx, R ₌ Tz and Q _{= Ty} , Tx, _Tz and _Ty are the partial derivatives of the travel time T.

where M and N are expressed as:

M=0.5(Q ₁ +Q ₂ )

N=Q ₁ Q ₂ -Q ₃ (7)

Among them, Q ₁ , Q ₂ , Q ₃ are represented by the following formulas:

where (a ₁₁ , a ₂₂ , a ₃₃ , a ₄₄ , a ₅₅ , a ₆₆ ) are the elastic moduli of anisotropic TI media.

Preferably, with reference to FIG. 2, FIG. 2 is a schematic flowchart of the factorization provided by the present invention, and step S2 includes steps S21 to S22, wherein:

In step S21, the solution of the function equation is decomposed into a first factor and a second factor;

In step S22, the first factor and the second factor are brought into the functional equation, and the derivatives of the unknown factors u(x, y, z) and T ₀ (x, y, z) with respect to the coordinates x, y, z are used Convert the formula 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, which is expressed as the following formula:

Among them, x, y, z are the coordinates in the computational domain, T(x, y, z) is the first factor, s(x, y, z) is the second factor, T ₀ (x, y) ,z) is predetermined to capture the source singularity, u(x,y,z) is an unknown factor smoothed around the source point, s ₀ (x,y,z) is the phase slowness model, α(x,y ,z) is the position factor corresponding to the phase slowness model s ₀ (x, y, z), the first factor captures the singularity of the source through analytical or numerical calculation, and the second factor is the first factor The necessary correction of , is smooth near the source point. Considering the decomposition of Equation (2), in a homogeneous medium, the function equation satisfies the following relationship:

的函数。where T ₀ is pre-determined to capture the source singularity, and u is an unknown factor smoothed around the source point. s is the phase slowness model, and s ₀ is the phase slowness of the source point (x ₀ , y ₀ , z ₀ ). They are all about the phase slowness direction

The function.

Substituting the above formula (9) into the above formula (2), the result obtained is expressed as the following formula:

With P=T _x , R=T _z and Q=T _y brought into the calculation, the result obtained is expressed as the following formula:

P=T _0x u+T ₀ p (12-a)

Q=T _0y u+T ₀ q (12-b)

R=T _0z u+T ₀ r (12-c)

where (p,q,r)=(u _x ,u _y ,u _z ) is the derivative of u with respect to x, y, and z, and T _0x , T _0y and T _0z are the derivatives of T ₀ with respect to x, y, and z , and further, equation (11) can be written in Hamiltonian form, expressed as the following equation:

f(x,y,z)=s _m (P,Q,R),(m=1,2,3) (13-b)

where (u _x , 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, u _x , u _y , u _z ) is the Hamiltonian operator, f(x, y, z) is the corresponding Hamiltonian derivative function, s _m (P, Q, R) is the phase slowness model brought into P, Q, R, P, Q and R are the derivatives of the first factor T(x, y, z) with respect to the coordinates x, y, z, T ₀ (x, y, z) is predetermined to capture the source singularity, u(x , y, z) are unknown factors smoothed around the source point.

Preferably, referring to Fig. 3, Fig. 3 is a schematic flowchart of solving the group velocity and phase velocity provided by the present invention, and step S3 includes steps S31 to S32, wherein:

In step S31, the phase slowness direction and the ray direction are set;

In step S32, the group velocity and the phase velocity are determined according to the relationship between the slowness vector and the group velocity vector.

Thus, through the phase slowness direction and the ray direction, the group velocity and the phase velocity are effectively solved in order to solve the equation of function.

In a specific embodiment of the present invention, step S3 specifically includes the following steps:

According to Eq. (13), T ₀ is still unknown and can be calculated as

为计算域中的坐标，

为源位置的坐标，

为沿射线方向

的群速度。in,

are the coordinates in the computational domain,

are the coordinates of the source location,

along the ray direction

the group velocity.

表达为下式：where the group velocity

Expressed as:

When the inverses of the phase velocities of the qP wave, qSV wave and qSH wave are respectively expressed as the following equations:

where the derivatives of M and N are expressed as:

和夹角

是隐式定义的，它们依赖于解T＝T₀·u。在后续描述中，我们将给出一个数值程序来计算沿给定射线方向的群速度和相速度。Phase slowness direction

and angle

are implicitly defined and they depend on the solution T=T ₀ ·u. In the following description, we will give a numerical procedure to calculate the group and phase velocities along a given ray direction.

It should be noted that, in the formula (13), the derivatives of T ₀ with respect to x, y, and z (T _0x , T _0y , T _0z ) are also unknowns. Unlike isotropic media, high-order finite difference methods cannot be used for calculations due to the possibility of high-frequency noise when the wavefront is not a perfect circle/sphere. Therefore, we give analytical expressions for computing these three variables, in a homogeneous anisotropic medium, as:

Among them, v ₀ is the phase velocity function at the source point, (p ₀ , q ₀ , r ₀ )=(T _0x , T _0y , T _0z ) is the derivative of T ₀ with respect to x, y, and z. Extract p ₀ , q ₀ , r ₀ from the square root of equation (18) in turn and do some algebraic operations. The result after the operation is expressed as follows:

根据θ和

的取值范围，我们可以确定这三个变量(p₀,q₀,r₀)的符号。It should be noted that, for a specific ray direction, θ and θ can be determined using the group velocity solution method demonstrated in the following

According to θ and

The value range of , we can determine the sign of these three variables (p ₀ , q ₀ , r ₀ ).

Among them, when calculating _T0 , 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 ray direction (19). The group velocity can be easily calculated if the slowness vector in the ray direction can be calculated. Previous work investigated how to calculate slowness vectors for specific ray directions.

和相位速度v_m，则慢度矢量可以表示为下式：If the phase slowness direction is given

and the phase velocity v _m , then the slowness vector can be expressed as:

和群速度向量

应满足以下关系:slowness vector

and the group velocity vector

The following relationship should be satisfied:

和射线方向

分别表示为下式：Phase slowness direction

and ray direction

They are respectively expressed as the following formulas:

得到下式：Divide both sides of equation (28) by

get the following formula:

相位慢度方向为

将式(25)改写为：Equation (25) will be used to calculate the phase and group velocities along a given ray direction. Let the ray direction be

The phase slowness direction is

Rewrite equation (25) as:

计算相位慢度方向

则可分别由式(19)和式(3)计算群速度U_m和相速度v_m。该方法适用于任意各向异性介质中的三种波型。This equation can be solved by the pseudo-position method or the bisection method, etc. Once along the ray direction

Calculate the phase slowness direction

Then, the group velocity U _m and the phase velocity v _m can be calculated from equations (19) and (3), respectively. The method is applicable to three wave modes in any anisotropic medium.

Preferably, with reference to FIG. 4 , FIG. 4 is a schematic flowchart of Lax-Friedrichs scanning provided by the present invention, and step S4 includes steps S41 to S42, wherein:

In step S41, the first-order Lax-Friedrichs format of the corresponding Hamiltonian form is established;

In step S42, according to the approximation value of the third-order upwind bias WENO (third-order weighted intrinsic non-oscillation scheme) corresponding to each grid point, the travel time value on the grid point is replaced, and a third-order Lax-Friedrichs based on WENO is established The static Hamiltonian-Jacobi equation of the format;

In step S43, the static Hamiltonian-Jacobi equation based on the third-order Lax-Friedrichs format of WENO is written in a compact form, and the travel time value on the grid point is updated according to the compact form.

Thus, using the Lax-Friedrichs finite difference scheme, iterative convergence can be guaranteed regardless of whether the Hamiltonian function is convex or concave. The basis function is numerically calculated by the fast frequency sweep scheme based on the third-order weighted intrinsically non-oscillating scheme (WENO) to achieve the purpose of fast sweep. The number of iterations has nothing to do with the grid size, and as the grid size approaches 0, The numerical solution converges to the expected weak solution.

In a specific embodiment of the present invention, formula (13) is expressed as the following formula:

where α _x , α _y and α _z are the viscosity coefficients.

In 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 u _i,j,k+1 are replaced by

和

是u_x的三阶迎风偏置WENO近似；

和

是u_y的三阶迎风偏置WENO近似；

和

是u_z的三阶迎风偏置WENO近似，其中：in,

and

is the third-order upwind bias WENO approximation of u _x ;

and

is the third-order upwind bias WENO approximation of u _y ;

and

is the third-order upwind bias WENO approximation of u _z , where:

when

as well as

when

和

的三阶WENO近似。ε是一个很小的正数，以避免除以0。Similarly, we can define

and

The third-order WENO approximation of . ε is a small positive number to avoid division by 0.

Then, the static Hamiltonian-Jacobi equation based on WENO's third-order Lax-Friedrichs format can be written in compact form:

为待更新的网格点(i,j,k)的数值解，

为同一点的旧数值解。in,

is the numerical solution of the grid point (i, j, k) to be updated,

for the old numerical solution for the same point.

Preferably, with reference to FIG. 5, FIG. 5 is a schematic flowchart of updating grid provided by the present invention, and step S43 includes steps S431 to S434, wherein:

In step S431, initialize and assign the grid points in the first-order Lax-Friedrichs format;

In step S432, the compact form is solved iteratively with 8 Gauss-Seidel (Gauss-Seidel iteration) frequency sweeps in alternating directions, and the grid is updated;

In step S433, on the six planes of the inner boundary and the outer boundary, a plurality of corresponding maximum and minimum value conditions are set, and a linear extrapolation method is used to determine the numerical solution of the point outside the domain;

In step S434, it is judged whether the L1-norm difference of the solution between two consecutive iterations satisfies the preset accuracy condition, and if so, the iteration is terminated.

Therefore, the Gauss-Seidel iterative solution with 8 frequency sweeps in alternating directions is used to achieve the purpose of fast scanning and solution, and at the same time, the numerical solution of the out-of-domain point is accurately calculated to ensure the efficiency of grid update.

In a specific embodiment of the present invention, the Hamilton-Jacobi equation (13) is solved using the Lax-Friedrichs scheme (19) based on the third-order WENO. The implementation procedure can be summarized as follows (Zhang et al., 2006; Luo et al., 2012):

Algorithm Schematic: Third-Order Lax-Friedrichs Scan Format for Anisotropic Factorial Function Equations

Initialization: Assign values on grid points with a distance less than or equal to 2h from the source point, because at least three grid points are required for the third-order WENO approximation. These values will be fixed during the iteration and all other grid points will be assigned large positive values.

Gauss-Seidel Iteration: Solve the discretized equation (33) with 8 Gauss-Seidel iterations swept in alternating directions:

At each grid point (i,j,k), it is updated according to the method described above.

Terminate: Terminate an iteration if the L1-norm difference of the solution between two consecutive iterations is less than the specified accuracy requirement.

When performing a third-order Lax-Friedrichs scanning scheme, linear extrapolation is employed near and at boundary points. On the six planes of the inner boundary, set the following conditions:

Then, on the six planes of the outer boundary, the following conditions are imposed:

Numerical solutions for points outside the computational domain are calculated by extrapolation, maximum and minimum methods. During the Gauss-Seidel iteration of the fast frequency sweep method, the (33) decomposition equation must be solved efficiently and the unknown factors calculated.

In a specific embodiment of the present invention, four test cases are given for description:

The first instance is a test of a sinusoidal anisotropy model. The computational domain is [-0.5,0.5]×[0.0,0.5]km rectangular domain, and the point source is located at x=0.0km, y=0.0km. The anisotropic parameters of the model (a11(x), a13(x), a33(x), a44(x), a66(x), θ0(x)) are listed in Table 1, which is shown below:

Table 1

Anisotropy parameter function to generate anisotropic parameters a₁₁ 5.2×{1.5+0.5sin(0.4π(z+0.01))sin(3π(x+0.05))} a₁₃ 0.98×{1.5+0.5sin(0.4π(z+0.01))sin(3π(x+0.05))} a₃₃ 4.0×{1.5+0.5sin(0.4π(z+0.01))sin(3π(x+0.05))} a₄₄ 1.0×{1.5+0.5sin(0.4π(z+0.01))sin(3π(x+0.05))} a₆₆ 1.0×{1.5+0.5sin(0.4π(z+0.01))sin(3π(x+0.05))} θ₀ 0⁰

The reference solution is computed on a very fine grid 1601×801 using the SPM model. The model is sampled on a 101 × 51 grid for the Lax-Friedrichs scanning method. Table 2 lists the L1 norm error and convergence order for four different grid sizes. The convergence order is greater than the first order, which verifies the effectiveness of the Lax-Friedrichs scanning method to solve the source singularity. Table 2 is as follows:

Table 2

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 Figures 6, 7, and 8, we compare the travel times of the three modes produced by the SPM and LAX methods.

The second instance tests 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, y=0.0 km. There are five elastic models (a11(x), a13(x), a33(x), a44(x), a66(x)), with θ0(x) = 00 (Fig. 4). Steep faults and high velocity salt bodies are present in the model. The model is resampled with a grid of 401×151. 9 to 11, FIG. 9 is a schematic diagram of the travel time table of the qP wave provided by the present invention, FIG. 10 is a schematic diagram of the travel time table of the qSV provided by the present invention, and FIG. 11 is a schematic diagram of the travel time table of the qSH provided by the present invention. As can be seen from the figure, the Lax-Friedrichs frequency sweep method can calculate the travel time table of qP, qSV and qSH waves well.

The third example is a test performed on an anticline formation with large dips on both sides of the formation. This model is called part of the BP TTI model. Referring to Fig. 12, Fig. 12 is a schematic diagram of the parameters of anisotropy provided by the present invention. It can be seen from the figure that (a11(x), a13(x), a33(x), a44(x), a66(x)) and the inclination angle θ0(x), where the calculation area is [-5.0, 5.0]×[0.0, 10.0]km rectangular area, and the source point is located at x=0.0km, y=0.0km. The model is sampled with a 251×251 grid. 13 to 15, FIG. 13 is a schematic diagram of the travel time of the qP wave under the Lax-Friedrichs scan provided by the present invention, FIG. 14 is a schematic diagram of the travel time of the qSV wave under the Lax-Friedrichs scan provided by the present invention, and FIG. 15 is According to the travel time diagram of the qSH wave under the Lax-Friedrichs scanning provided by the present invention, it can be seen that the Lax-Friedrichs scanning method successfully generates the travel times of three wave modes.

The fourth example is to test the proposed method in a three-layer anisotropic medium. A sloped interface exists in the model. Anisotropy parameters (a11(x), a13(x), a33(x), a44(x), a66(x)), inclination angle θ0(x) and azimuth angle φ0(x) are listed in Table 3, Table 3 is shown below:

table 3

16 to 18, FIG. 16 is a three-dimensional travel time contour diagram of a qP wave provided by the present invention, FIG. 17 is a three-dimensional travel time contour diagram of a qSV wave provided by the present invention, and FIG. 18 is the present invention. Three-dimensional travel time contour plots of qSH waves are provided. It can be seen from the figure that the 3D extension of the adopted 2D model is used to calculate the travel time and ray path of elastic waves, and the computational 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.0km, x=0.0km, z=0.0km, and the grid of Lax-Friedrichs scanning method is 101×41×101. Three-dimensional travel time contour maps of qP, qSV, and qSH waves were plotted on the anisotropy model (a11(x) model) and the three-dimensional travel time cube (Figures 16-18), respectively.

Example 2

An embodiment of the present invention provides a Lax Friedrichs device for scanning the travel time of seismic waves. Referring to FIG. 19 , FIG. 19 is a schematic structural diagram of the device for scanning the travel time of seismic waves by Lax Friedrichs provided by the present invention. The above device 1900 includes:

The modeling unit 1901 is used to establish the function equation of seismic waves in anisotropic TI medium by introducing the slowness vector;

a conversion unit 1902, configured to decompose the solution of the functional equation into the product of two factors, and bring the product of the two factors into the functional equation, and convert it into a corresponding Hamiltonian form;

The solving unit 1903 is used to solve the group velocity and the phase velocity according to the slowness vector;

The scanning unit 1904 is configured to use the Lax-Friedrichs scanning format value to solve the functional equation and determine the travel time of the seismic wave.

Example 3

An embodiment of the present invention provides a computer-readable storage medium on which a computer program is stored. When the computer program is executed by a processor, the above-mentioned Lax Friedrichs scanning seismic wave travel time method is implemented.

The invention discloses a Lax Friedrichs scanning seismic wave travel time method, which can calculate the first arrival travel time of three wave types, namely qP-, qSV- and qSH-waves, and uses a factor decomposition method to process a finite difference equation solver The headwind singularity of A necessary modification. The Lax-Friedrichs finite-difference scheme is further adopted, which can guarantee the iterative convergence regardless of whether the Hamiltonian function is convex or concave. The basis function is numerically calculated by the fast frequency sweep scheme based on the third-order weighted intrinsically non-oscillating scheme (WENO), and iteratively The degree is independent of the grid size, and as the grid size approaches 0, the numerical solution converges to the desired weak solution. Therefore, the present invention uses the effective factorization method to give the Hamiltonian form of the factorial anisotropy equations of qP-, qSV- and qSH-waves in TTI medium, and designs a combination of the Lax-Friedrichs finite difference scheme. The fast scanning method ensures efficient and accurate calculation of travel time.

The technical scheme of the present invention adopts the method of factoring and fast scanning, discretizes the Hamiltonian form of the anisotropic functional equation through the third-order weighted basic non-oscillation approximation and Lax-Friedrichs format, and then uses the fast scanning method to iteratively solve the discrete system that can calculate very precise first arrival times. The anisotropic functional equation does not have weak anisotropy assumptions in phase and group velocity, so this method is suitable for inclined transversely isotropic (TTI) media with arbitrary anisotropy, with strong flexibility and generalization. High, both correctness and effectiveness, ensure the accuracy and rapidity of travel time calculation, so as to carry out effective geological survey.

The above description is only a preferred embodiment of the present invention, but the protection scope of the present invention is not limited to this. Substitutions should be covered within the protection scope of the present invention.

Claims (6)

1. a method for Lax Friedrichs scanning seismic wave travel time, is characterized in that, comprises: By introducing the slowness vector, the function equation of seismic wave in anisotropic TI medium is established; Decompose the solution of the functional equation into a product of two factors, and bring the product of the two factors into the functional equation, and convert it into a corresponding Hamiltonian form; According to the slowness vector, solve the group velocity and the phase velocity; The travel time of the seismic wave is determined by solving the functional equation using the Lax-Friedrichs scan format numerical values. 2. The method for scanning the travel time of seismic waves by Lax Friedrichs according to claim 1, wherein the solution of the functional equation is decomposed into the product of two factors, and the product of the two factors is banded Into the functional equation, the conversion to the corresponding Hamiltonian form includes: The solution of the functional equation is decomposed into a first factor and a second factor, and the first factor and the second factor are expressed as:

Among them, x, y, z are the coordinates in the computational domain, T(x, y, z) is the first factor, s(x, y, z) is the second factor, T ₀ (x, y) ,z) is predetermined to capture the source singularity, u(x,y,z) is an unknown factor smoothed around the source point, s ₀ (x,y,z) is the phase slowness model, α(x,y , z) is the position factor corresponding to the phase slowness model s ₀ (x, y, z); Bring the first factor and the second factor into the functional equation, and capture the source singularity T ₀ (x, y) according to the unknown factor u(x, y, z) and the predetermined , z) formula transformation is carried out on the derivatives of the coordinates x, y, and z, and it is determined that the corresponding Hamiltonian form is:

f(x,y,z)=s _m (P,Q,R),(m=1,2,3) Among them, (u _x , u _y , u _z ) are the derivatives of the unknown factor u(x, y, z) with respect to the coordinates x, y, z respectively, T ₀ is predetermined to capture the source singularity, H(x ,y,z,u,u _x ,u _y ,u _z ) is the Hamiltonian operator, f(x,y,z) is the corresponding Hamiltonian derivative function, s _m (P,Q,R) is the bring-in Phase slowness model after P, Q, R, where P, Q and R are derivatives of the first factor T(x, y, z) with respect to coordinates x, y, z. 3. The method for Lax Friedrichs scanning seismic wave travel time according to claim 2, wherein, according to the slowness vector, solving the group velocity and the phase velocity comprises: Set the phase slowness direction and ray direction; The group velocity and the phase velocity are determined from the relationship between the slowness vector and the group velocity vector. 4. The method for Lax Friedrichs scanning seismic wave travel time according to claim 3, wherein, in the corresponding Hamiltonian form, the determination of the derivative of the travel time with respect to the coordinate axes x, y, and z comprises: The analytical expression for determining the derivative of the travel time with respect to the coordinate axes x, y, and z is expressed as:

Wherein, v ₀ is the phase velocity function at the source point, (p ₀ , q ₀ , r ₀ )=(T _0x , T _0y , T _0z ), which are the derivatives of the travel time with respect to x, y, and z, respectively; The derivatives p ₀ , q ₀ , r ₀ of the travel time with respect to x, y, z are sequentially extracted from the square root of the analytical expression, and through algebraic operations, the derivatives of the travel time with respect to x, y, z are determined The corresponding absolute values of p ₀ , q ₀ , r ₀ are expressed as:

in,

are the tilt and azimuth angles in the direction of phase slowness. 5. the method for Lax Friedrichs scanning seismic wave travel time according to claim 4, is characterized in that, described utilizing Lax-Friedrichs scanning format numerical value, solving described function equation comprises: establishing a first-order Lax-Friedrichs scheme of the corresponding Hamiltonian form; According to the approximation value of the third-order upwind bias WENO corresponding to each grid point, the travel time value on the grid point is replaced, and the static Hamiltonian-Jacobi equation of the third-order Lax-Friedrichs format based on WENO is established; The static Hamiltonian-Jacobi equation based on the WENO third-order Lax-Friedrichs format is written in a compact form, and the travel time values on the grid points are updated according to the compact form. 6. The method for scanning travel time of seismic waves by Lax Friedrichs according to claim 5, wherein the updating the travel time value on the grid point according to the compact form comprises: initializing and assigning grid points in the first-order Lax-Friedrichs format; Solve the compact form with 8 Gauss-Seidel iterations swept in alternating directions, and perform mesh updates; On the six planes of the inner boundary and the outer boundary, set the corresponding maximum and minimum value conditions, and use the linear extrapolation method to determine the numerical solution of the point outside the domain; It is judged whether the L1 norm difference of the solution between two consecutive iterations satisfies the preset accuracy condition, and if so, the iteration is terminated.

Images