A method, device, medium and product for determining travel time of seismic waves in three-dimensional medium

Abstract

The present invention discloses a method, device, medium and product for determining seismic wave travel time in a three-dimensional medium, and relates to the technical field of geophysical seismic data processing. The method comprises: dividing the three-dimensional medium to be measured into a plurality of cubic discrete units; using the finite difference method, according to the eikonal equation, determining the initial value of the seismic wave travel time; according to the initial value of the seismic wave travel time, based on the discontinuous finite element method, respectively determining the polynomial of each cubic discrete unit; updating the coefficients of multiple groups of polynomials, determining multiple groups of final coefficients; determining the final value of the seismic wave travel time according to the multiple groups of final coefficients. The present invention determines the final value of the seismic wave travel time by combining the finite difference method and the discontinuous finite element method, taking into account the accuracy of the seismic wave travel time results and the efficiency of calculation, while improving the accuracy of the seismic wave travel time results, improving the efficiency of calculation, and further providing an important data basis for detecting the structure and lithology of underground media.

Description

Technical Field

[0001] The present invention relates to the technical field of geophysical seismic data processing, and in particular to a method, device, medium and product for determining travel time of seismic waves in a three-dimensional medium. Background Art

[0002] The seismic exploration method in geophysics is to stimulate seismic waves through artificial sources, so that they propagate in the underground medium and transmit the earth structure information to the ground. The processing and analysis of this information is crucial for detecting the structure and lithology of the underground medium. Among them, seismic wave travel time is one of the key information, which is widely used in data processing processes such as pre-stack migration, seismic imaging, travel time tomography inversion, and seismic positioning. The travel time of seismic waves and the propagation speed of seismic waves in the medium are combined to obtain the propagation distance. Therefore, the calculation of seismic wave travel time can be applied to seismic imaging and migration, mainly to solve the time difference caused by offset distance and inclination to complete imaging or homing. Seismic wave travel time is also used for inversion tomography. Inversion tomography is mainly used for matching and comparison of direct problems and inverse problems, mapping and velocity update of travel kernel functions, and seismic positioning. Therefore, when seismic wave travel time is used for inversion tomography, it can help realize underground medium structural imaging.

[0003] At present, the methods for calculating seismic wave travel time mainly include two-point ray tracing method, shortest path method and numerical solution based on eikonal equation. Although the two-point ray tracing method is simple, it is inefficient and not suitable for complex media. The shortest path method has high accuracy, but it has low computational efficiency in complex models and requires a large amount of computing memory. The numerical method based on the eikonal equation does not need to calculate the ray path and has the characteristics of high efficiency, stability and easy implementation. Among them, the eikonal equation finite difference method proposed by Vidale (1988), the fast march method of Sethian et al. (1999), and the fast scanning method algorithm of Zhao et al. (2004) are typical representatives. However, in a three-dimensional complex model, when analyzing complex structures and velocity changes, the finite difference method alone is used to calculate the seismic wave travel time, which cannot well characterize the structure, resulting in low accuracy of the calculated seismic wave travel time; the discontinuous finite element method alone is used to calculate the seismic wave travel time, which will lead to low computational efficiency, and ultimately lead to the inability to balance the accuracy of the seismic wave travel time results and the efficiency of the calculation. Summary of the invention

[0004] The purpose of the present invention is to provide a method, device, medium and product for determining seismic wave travel time in a three-dimensional medium, so as to solve the problem in the prior art that, when analyzing complex structures and velocity changes under a three-dimensional complex model, the finite difference method is used alone to calculate the seismic wave travel time, which cannot well characterize the structure, resulting in low accuracy of the calculated seismic wave travel time; the discontinuous finite element method is used alone to calculate the seismic wave travel time, which leads to low calculation efficiency, and ultimately leads to the problem that it is impossible to take into account both the accuracy of the seismic wave travel time results and the efficiency of the calculation.

[0005] To achieve the above object, the present invention provides the following solutions:

[0006] A method for determining travel time of seismic waves in a three-dimensional medium, comprising:

[0007] Divide the three-dimensional medium to be measured into a plurality of discrete cubic units;

[0008] Using the finite difference method and the eikonal equation, the initial value of the seismic wave travel time is determined;

[0009] According to the initial value of the seismic wave travel time, based on the discontinuous finite element method, respectively determine the polynomial of each of the cube discrete units;

[0010] updating the coefficients of the multiple sets of polynomials to determine multiple sets of final coefficients;

[0011] The final value of the seismic wave travel time is determined according to the multiple groups of final coefficients.

[0012] Optionally, the eikonal equation is:

[0013]

[0014] Wherein, V(x,y,z) is the seismic wave propagation velocity of the three-dimensional medium to be measured at the target position (x,y,z) of the seismic wave propagation, and t(x,y,z) is the seismic wave travel time from the earthquake source to the target position (x,y,z) of the seismic wave propagation in the three-dimensional medium to be measured.

[0015] Optionally, the polynomial is:

[0016]

[0017] Wherein, t is the average travel time of the cube discrete unit, ξ is the travel time gradient of the cube discrete unit in the X direction of the spatial rectangular coordinate system, η is the travel time gradient of the cube discrete unit in the Y direction, γ is the travel time gradient of the cube discrete unit in the Z direction, wherein the spatial rectangular coordinate system is established with the center of the three-dimensional medium to be measured as the origin and the direction from the origin to the earthquake source as the X axis, t is the travel time of the seismic wave, h is the side length of the cube discrete unit, (x, y, z) is the target position of the seismic wave propagation, (xi ,y j , z k ) is the center coordinate of any cubic discrete unit, i is the i-th cubic discrete unit in the X direction, j is the j-th cubic discrete unit in the Y direction, and k is the k-th cubic discrete unit in the Z direction.

[0018] Optionally, the expressions of the coefficients are:

[0019]

[0020] in, is the travel time from the earthquake source to the first node of any cubic discrete unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the second node of any discrete cubic unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the third node of any discrete cubic unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the fourth node of any discrete cubic unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the fifth node of any cube discrete unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the sixth node of any discrete cubic unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the seventh node of any discrete cubic unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the eighth node of any cube discrete unit in the three-dimensional medium to be measured, and the first node, the second node, the third node, the fourth node, the fifth node, the sixth node, the seventh node and the eighth node are the eight vertices of the cube discrete unit.

[0021] Optionally, updating the coefficients of the multiple sets of polynomials to determine the multiple sets of final coefficients specifically includes:

[0022] Determining a set of travel time equations based on the polynomial and the weak form of the eikonal equation;

[0023] Based on the travel time equations, multiple sets of coefficients of the polynomials are updated to determine multiple sets of final coefficients.

[0024] Optionally, based on the travel time equations, updating multiple sets of coefficients of the polynomials to determine multiple sets of final coefficients specifically includes:

[0025] Based on the travel time equation group, the coefficients of the multiple sets of polynomials are updated, and when the coefficients no longer change, the current coefficients of the multiple sets of polynomials are determined to be final coefficients.

[0026] Optionally, the travel time equations are:

[0027]

[0028] in: The coordinates from the earthquake source to the center of the three-dimensional medium to be measured are (x i ,y j , z k ) is the average travel time of the cube discrete unit, ξ i,j,k The center coordinates are (x i ,y j , z k ) in the X direction, η i,j,k The center coordinates are (x i ,y j , z k ) is the travel time gradient of the cubic discrete unit in the Y direction, γ i,j,k The center coordinates are (x i ,y j , z k ) is the travel time gradient of the cube discrete unit in the Z direction; h is the side length of the cube discrete unit, R 1i,j,k is the first simplified formula, R 2i,j,k is the second simplified formula, R 3i,j,k is the third simplified formula, R 4i,j,k is the fourth simplified formula, A 1i,j,k is the fifth simplified formula, A 2i,j,k is the sixth simplified formula, A 3i,j,k is the seventh simplified formula, A 4i,j,k is the eighth simplified formula, A 5i,j,k is the ninth simplified formula, A 6i,j,k is the tenth simplified formula, A 7i,j,k is the eleventh simplified formula;

[0029]

[0030] The coordinates from the earthquake source to the center of the three-dimensional medium to be measured are (x i-1 ,y j , z k ) is the average travel time of the discrete cube unit, ξ i-1,j,k The center coordinates are (x i-1 ,y j , z k ) in the X direction, η i-1,j,k The center coordinates are (x i-1 ,y j , z k ) is the travel time gradient of the cubic discrete unit in the Y direction, γ i-1,j,k The center coordinates are (x i-1,y j , z k )’s travel time gradient in the Z direction of the cube discrete unit; The coordinates from the earthquake source to the center of the three-dimensional medium to be measured are (x i ,y j+1 , z k ) is the average travel time of the discrete cube unit, ξ i,j+1,k The center coordinates are (x i ,y j+1 , z k ) in the X direction, η i,j+1,k The center coordinates are (x i ,y j+1 , z k ) is the travel time gradient of the cubic discrete unit in the Y direction, γ i,j+1,k The center coordinates are (x i ,y j+1 , z k )’s travel time gradient in the Z direction of the cube discrete unit; The coordinates from the earthquake source to the center of the three-dimensional medium to be measured are (x i ,y j-1 , z k ) is the average travel time of the discrete cube unit, ξ i,j-1,k The center coordinates are (x i ,y j-1 , z k ) in the X direction, η i,j-1,k The center coordinates are (x i ,y j-1 , z k ) is the travel time gradient of the cubic discrete unit in the Y direction, γ i,j-1,k The center coordinates are (x i ,y j-1 , z k )’s travel time gradient in the Z direction of the cube discrete unit; The coordinates from the earthquake source to the center of the three-dimensional medium to be measured are (x i ,y j , z k+1 ) is the average travel time of the discrete cube unit, ξ i,j,k+1 The center coordinates are (x i ,y j , z k+1 ) in the X direction, η i,j,k+1 The center coordinates are (x i ,y j , z k+1 ) is the travel time gradient of the cubic discrete unit in the Y direction, γi,j,k+1 The center coordinates are (x i ,y j , z k+1 )’s travel time gradient in the Z direction of the cube discrete unit; The coordinates from the earthquake source to the center of the three-dimensional medium to be measured are (x i ,y j , z k-1 ) is the average travel time of the discrete cube unit, ξ i,j,k-1 The center coordinates are (x i ,y j , z k-1 ) in the X direction, η i,j,k-1 The center coordinates are (x i ,y j , z k-1 ) is the travel time gradient of the cubic discrete unit in the Y direction, γ i,j,k-1 The center coordinates are (x i ,y j , z k-1 )’s cubic discrete unit in the Z direction; α r ,α l ,α f ,α b ,α n ,α s are the flux coefficients of the six faces of the cubic discrete unit;

[0031]

[0032] A computer device includes: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the method for determining the travel time of seismic waves in three-dimensional media.

[0033] A computer-readable storage medium stores a computer program, which, when executed by a processor, implements the method for determining the travel time of seismic waves in a three-dimensional medium.

[0034] A computer program product comprises a computer program, which, when executed by a processor, implements the method for determining the travel time of seismic waves in a three-dimensional medium.

[0035] According to the specific embodiments provided by the present invention, the present invention discloses the following technical effects:

[0036] The method, device, medium and product for determining the travel time of seismic waves in three-dimensional media provided by the present invention are as follows: the three-dimensional medium to be measured is divided into a plurality of cubic discrete units; the initial value of the travel time of seismic waves is determined according to the eikonal equation using the finite difference method; the polynomial of each cubic discrete unit is determined respectively according to the initial value of the travel time of seismic waves based on the discontinuous finite element method; the coefficients of multiple groups of polynomials are updated to determine multiple groups of final coefficients; the final value of the travel time of seismic waves is determined according to the multiple groups of final coefficients. The present invention determines the final value of the travel time of seismic waves by combining the finite difference method and the discontinuous finite element method, taking into account the accuracy of the results of the travel time of seismic waves and the efficiency of calculation, improving the accuracy of the results of the travel time of seismic waves while improving the efficiency of calculation, and further providing an important data basis for detecting the structure and lithology of underground media. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0038] Figure 1 A flow chart of the method for determining the travel time of seismic waves in three-dimensional media provided by the present invention;

[0039] Figure 2 A schematic diagram of a discretized grid of a three-dimensional medium of the present invention;

[0040] Figure 3 It is a schematic diagram of a cubic discrete unit in the present invention;

[0041] Figure 4 It is a travel time contour map in the three-dimensional medium model of the present invention;

[0042] Figure 5 It is a travel time contour map on the plane of z=200m in the three-dimensional medium model of the present invention;

[0043] Figure 6 The error diagram between the seismic wave travel time calculated by the finite difference method in the present invention and the analytical solution;

[0044] Figure 7 This is the error diagram between the seismic wave travel time calculated by the discontinuous finite element method in the present invention and the analytical solution. DETAILED DESCRIPTION

[0045] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0046] The purpose of the present invention is to provide a method for determining the travel time of seismic waves in a three-dimensional medium, which aims to determine the final value of the travel time of seismic waves by combining a finite difference method and a discontinuous finite element method, taking into account the accuracy of the results of the travel time of seismic waves and the efficiency of calculation, while improving the accuracy of the results of the travel time of seismic waves, improving the efficiency of calculation, and further providing an important data basis for detecting the structure and lithology of underground media. The present invention has good applicability when solving problems under complex geological conditions, and provides more accurate underground structure information for the field of geophysics.

[0047] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0048] Example 1

[0049] like Figure 1 As shown, the present invention provides a method for determining the travel time of seismic waves in a three-dimensional medium, comprising:

[0050] Step 101: divide the three-dimensional medium to be measured into a plurality of discrete cubic units.

[0051] The three-dimensional medium model is divided using cubic units with a side length of h to form a discretized grid of the three-dimensional medium. Specifically, a spatial rectangular coordinate system can be established with the center of the three-dimensional medium to be measured as the origin and the direction from the origin to the source as the X-axis; the three-dimensional medium to be measured can also be approximated as a cuboid or a cube, with a vertex of the three-dimensional medium to be measured as the origin and the three sides corresponding to the vertex as the X-axis, Y-axis and Z-axis respectively to establish a spatial rectangular coordinate system, such as Figure 2 As shown in . If the number of discrete grids of the three-dimensional medium in the three directions of the spatial rectangular coordinate system X, Y, and Z is Nx, Ny, and Nz respectively, then the number of corresponding cubic discrete units is (Nx-1)×(Ny-1)×(Nz-1). The cubic discrete unit is as follows Figure 3 As shown, its center coordinates are (x i ,y i 、z i ), the node coordinate on the side close to the origin in the X direction is x i-1 / 2 , the node coordinate on the positive side of the X-axis is x i+1 / 2 ; The node coordinate on the Y direction close to the origin is y j-1 / 2，The node coordinate on the positive side of the Y axis is y j+1 / 2 ; The node coordinate on the side close to the origin in the Z direction is z k-1 / 2 , the node coordinate on the Z direction close to the positive direction of the Z axis is z k+1 / 2 .

[0052] Step 102, using the finite difference method and the eikonal equation, determine the initial value of the seismic wave travel time.

[0053] Specifically, the Eikonal equation is solved based on the fast scanning method of finite differences to obtain the initial value of the seismic wave travel time.

[0054] In a three-dimensional medium, the Eikonal equation is used to describe the travel time required for a seismic wave to travel from the source to the position (x, y, z). The Eikonal equation is:

[0055]

[0056] Wherein, V(x,y,z) is the seismic wave propagation velocity of the three-dimensional medium to be measured at the target position (x,y,z) of the seismic wave propagation, and t(x,y,z) is the seismic wave travel time from the earthquake source to the target position (x,y,z) of the seismic wave propagation in the three-dimensional medium to be measured.

[0057] The fast scanning method based on the finite difference method is used to solve the eikonal equation, and the grid points (x i ,y i 、z i ) is the initial value of the seismic wave travel time t i,j,k , the initial value of the seismic wave travel time is of first-order accuracy.

[0058] Step 103, according to the initial value of the seismic wave travel time, based on the discontinuous finite element method, respectively determine the polynomial of each cubic discrete unit.

[0059] Among them, the discontinuous finite element method uses the basis functions in the polynomial on the cubic element 1, The coefficients in the polynomial to approximate the travel time of seismic waves.

[0060] Specifically, the polynomial is:

[0061]

[0062] in, is the average travel time of the cube discrete unit, ξ is the travel time gradient of the cube discrete unit in the X direction of the spatial rectangular coordinate system, η is the travel time gradient of the cube discrete unit in the Y direction, and γ is the travel time gradient of the cube discrete unit in the Z direction, wherein the spatial rectangular coordinate system is established with the center of the three-dimensional medium to be measured as the origin and the direction from the origin to the source as the X axis, t is the travel time of the seismic wave, h is the side length of the cube discrete unit, (x, y, z) is the target position of the seismic wave propagation, (x i ,y j , z k ) is the center coordinate of any cubic discrete unit, i is the i-th cubic discrete unit in the X direction, j is the j-th cubic discrete unit in the Y direction, and k is the k-th cubic discrete unit in the Z direction.

[0063] The essence of this polynomial is to use the background travel time field t and the disturbance travel time field (the sum of the disturbances in different coordinate system directions) to represent the real travel time within the discrete unit of the cube.

[0064] When the seismic wave propagates to a certain position (x, y, z), the center coordinates of the discrete cube unit at that location are (x i ,y i 、z i ), then the basis function in the above polynomial is becomes a known number, so the solution of the seismic wave travel time is converted into the coefficients of the polynomial Solution for ξ,η,γ.

[0065] The coefficients of the polynomial corresponding to the polynomial in the cube discrete unit are initialized using the initial values ​​of the seismic wave travel time.

[0066] The expressions of each coefficient are:

[0067]

[0068] in, is the travel time from the earthquake source to the first node of any cubic discrete unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the second node of any discrete cubic unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the third node of any discrete cubic unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the fourth node of any discrete cubic unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the fifth node of any cube discrete unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the sixth node of any discrete cubic unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the seventh node of any discrete cubic unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the eighth node of any cube discrete unit in the three-dimensional medium to be measured, and the first node, the second node, the third node, the fourth node, the fifth node, the sixth node, the seventh node and the eighth node are the eight vertices of the cube discrete unit.

[0069] Specifically, the first node is the node of the cube discrete unit in the negative direction of the X axis (the node in the negative direction of the Y axis, the node in the negative direction of the Z axis), the second node is the node of the cube discrete unit in the negative direction of the X axis (the node in the positive direction of the Y axis, the node in the negative direction of the Z axis), the third node is the node of the cube discrete unit in the negative direction of the X axis (the node in the negative direction of the Y axis, the node in the positive direction of the Z axis), and the fourth node is the node of the cube discrete unit in the negative direction of the X axis (the node in the positive direction of the Y axis, the node in the positive direction of the Z axis). ), the fifth node is the node of the cube discrete unit in the positive direction of the X-axis (the node in the negative direction of the Y-axis, the node in the negative direction of the Z-axis), the sixth node is the node of the cube discrete unit in the positive direction of the X-axis (the node in the positive direction of the Y-axis, the node in the negative direction of the Z-axis), the seventh node is the node of the cube discrete unit in the positive direction of the X-axis (the node in the negative direction of the Y-axis, the node in the positive direction of the Z-axis), and the eighth node is the node of the cube discrete unit in the positive direction of the X-axis (the node in the positive direction of the Y-axis, the node in the positive direction of the Z-axis).

[0070] According to the above formula, the basis function coefficients of all discrete units are obtained, and the coefficients of the polynomials of each single discrete unit are expressed as ξ i,j,k ,η i,j,k ,γ i,j,k .

[0071] Step 104, updating multiple sets of polynomial coefficients to determine multiple sets of final coefficients.

[0072] Step 104 specifically includes: determining a travel time equation system based on the polynomial and the weak form of the eikonal equation; updating multiple sets of coefficients of the polynomial based on the travel time equation system to determine multiple sets of final coefficients. The gradient of the solution of the polynomial is fixed to the gradient calculated by the finite difference method to reduce the degree of freedom of the equation and avoid solving multiple parameters at the same time.

[0073] Based on the travel time equations, the coefficients of multiple sets of polynomials are updated to determine multiple sets of final coefficients, specifically including: based on the travel time equations, the coefficients of multiple sets of polynomials are updated, and when the coefficients no longer change, the current coefficients of the multiple sets of polynomials are determined to be the final coefficients.

[0074] Specifically, according to the discontinuous finite element standard weighted residual weak form representation of the eikonal equation, the travel time solution equation group can be constructed:

[0075] The travel time equations are:

[0076]

[0077] in: The coordinates from the earthquake source to the center of the three-dimensional medium to be measured are (x i ,y j , z k ) is the average travel time of the discrete cube unit, ξ i,j,k The center coordinates are (x i ,y j , z k ) in the X direction, η i,j,k The center coordinates are (x i ,y j , z k ) is the travel time gradient of the cubic discrete unit in the Y direction, γ i,j,k The center coordinates are (x i ,y j , z k ) is the travel time gradient of the cube discrete unit in the Z direction; h is the side length of the cube discrete unit, R 1i,j,k is the first simplified formula, R 2i,j,k is the second simplified formula, R 3i,j,k is the third simplified formula, R 4i,j,k is the fourth simplified formula, A 1i,j,k is the fifth simplified formula, A 2i,j,k is the sixth simplified formula, A 3i,j,k is the seventh simplified formula, A 4i,j,k is the eighth simplified formula, A 5i,j,k is the ninth simplified formula, A 6i,j,k is the tenth simplified formula, A 7i,j,k This is the eleventh simplified formula.

[0078] It should be noted that in order to simplify the formula, R 1i,j,k Defined as the first simplified formula, R 2i,j,k Defined as the second simplified formula, R 3i,j,k Defined as the third simplified formula, R 4i,j,k Defined as the fourth simplified formula, A 1i,j,k Defined as the fifth simplified formula, A 2i,j,k Defined as the sixth simplified formula, A 3i,j,k Defined as the seventh simplified formula, A 4i,j,k Defined as the eighth simplified formula, A 5i,j,k Defined as the ninth simplified formula, A 6i,j,kDefined as the tenth simplified formula, A 7i,j,k Defined as the eleventh simplified formula.

[0079]

[0080] The coordinates from the earthquake source to the center of the three-dimensional medium to be measured are (x i-1 ,y j , z k ) is the average travel time of the discrete cube unit, ξ i-1,j,k The center coordinates are (x i-1 ,y j , z k ) in the X direction, η i-1,j,k The center coordinates are (x i-1 ,y j , z k ) is the travel time gradient of the cubic discrete unit in the Y direction, γ i-1,j,k The center coordinates are (x i-1 ,y j , z k )’s travel time gradient in the Z direction of the cube discrete unit; The coordinates from the earthquake source to the center of the three-dimensional medium to be measured are (x i ,y j+1 , z k ) is the average travel time of the discrete cube unit, ξ i,j+1,k The center coordinates are (x i ,y j+1 , z k ) in the X direction, η i,j+1,k The center coordinates are (x i ,y j+1 , z k ) is the travel time gradient of the cubic discrete unit in the Y direction, γ i,j+1,k The center coordinates are (x i ,y j+1 , z k )’s travel time gradient in the Z direction of the cube discrete unit; The coordinates from the earthquake source to the center of the three-dimensional medium to be measured are (x i ,y j-1 , z k ) is the average travel time of the discrete cube unit, ξ i,j-1,k The center coordinates are (x i ,y j-1 , z k ) in the X direction, η i,j-1,k The center coordinates are (x i ,yj-1 , z k ) is the travel time gradient of the cubic discrete unit in the Y direction, γ i,j-1,k The center coordinates are (x i ,y j-1 , z k )’s travel time gradient in the Z direction of the cube discrete unit; The coordinates from the earthquake source to the center of the three-dimensional medium to be measured are (x i ,y j , z k+1 ) is the average travel time of the discrete cube unit, ξ i,j,k+1 The center coordinates are (x i ,y j , z k+1 ) in the X direction, η i,j,k+1 The center coordinates are (x i ,y j , z k+1 ) is the travel time gradient of the cubic discrete unit in the Y direction, γ i,j,k+1 The center coordinates are (x i ,y j , z k+1 )’s travel time gradient in the Z direction of the cube discrete unit; The coordinates from the earthquake source to the center of the three-dimensional medium to be measured are (x i ,y j , z k-1 ) is the average travel time of the discrete cube unit, ξ i,j,k-1 The center coordinates are (x i ,y j , z k-1 ) in the X direction, η i,j,k-1 The center coordinates are (x i ,y j , z k-1 ) is the travel time gradient of the cubic discrete unit in the Y direction, γ i,j,k-1 The center coordinates are (x i ,y j , z k-1 )’s cubic discrete unit in the Z direction; α r ,α l ,α f ,α b ,α n ,α s are the flux term coefficients of the six faces of the cubic discrete unit.

[0081]

[0082] Specifically, the update of the basis function coefficients of the cubic discrete unit requires the information of the adjacent cubic discrete units, and the constant α r ,α l ,α f ,α b ,α n ,α s The flux coefficients corresponding to the six faces of the cube discrete unit are used to determine whether the cube discrete unit exchanges information with the adjacent cube discrete unit. If the travel time of the seismic wave on a cube discrete unit increases monotonically from left to right, there is a t x >0, then with a r = 0. At this time, the travel time calculation of the cubic discrete unit needs to be completed using the information of the cubic discrete unit adjacent to the cubic discrete unit along the negative direction of the X axis, which is consistent with the propagation direction of the seismic wave.

[0083] Furthermore, when solving the travel time equations, the present invention uses the first-order characteristics of the gradient term to simplify the calculation and avoid solving multiple parameters at the same time. The solution of the travel time equations is smooth, and the initial value and gradient of the seismic wave travel time are both first-order accurate. Approximately (x i ,y j ,z k ) at t x ,t y ,t z (gradient), ξ is no longer updated during the solution process ijk ,η ijk γ ijk , but fix them to:

[0084] ξ ijk =t x,ijk h,η ijk =t y,ijk h,γ ijk =t z,ijk h.

[0085] At this point, we only need to use the first formula in the travel time equations:

[0086]

[0087] renew Specifically, the formula can be transformed into:

[0088]

[0089] in, The coefficients to be updated for the current cube discrete unit.

[0090] It should be noted that when updating the coefficients, the three-dimensional medium model needs to be scanned, and the cubic discrete units need to be traversed to update the coefficients. The traversal directions include X, Y and Z positive and negative directions, which are traversed alternately in sequence.

[0091] Step 105, determining the final value of the seismic wave travel time according to the multiple sets of final coefficients.

[0092] After traversing the cube discrete unit update coefficients, get the final And with i,j,k ,η i,j,k ,γ i,j,k Combinatorial, polynomial-based The seismic wave travel time at all positions in the three-dimensional medium can be obtained, thereby determining the final value of the seismic wave travel time, that is, the seismic wave travel time with second-order accuracy.

[0093] In summary, the present invention tests the accuracy of the solution of seismic wave travel time in a three-dimensional medium model. The three-dimensional medium model can be a three-dimensional uniform seismic medium model with a model size of 400m×400m and a velocity of V=1000m / s. The location of the earthquake source (x s ,y s ,z s ) is (200m, 200m, 200m), the calculation formula of the analytical solution of seismic wave travel time is:

[0094]

[0095] The calculation formula of travel time gradient is:

[0096]

[0097] Among them, x s is the coordinate of the earthquake source in the X direction, y s is the coordinate of the earthquake source in the Y direction, z s is the coordinate of the earthquake source in the Z direction.

[0098] First, the finite difference method is used to solve the seismic wave travel time in the three-dimensional uniform seismic medium model. The grid spacing Δh is set to 1m, 2m and 4m respectively. The calculated seismic wave travel time is shown in Table 1. Table 1 is the seismic wave travel time table obtained by the finite difference method, which includes the error and convergence order of the analytical solution. For L 1 Norm and L ∞ The travel time accuracy calculated by the finite difference method meets the first-order convergence.

[0099] Table 1 Seismic wave travel time obtained by finite difference method

[0100]

[0101] Table 2 is the travel time gradient table obtained by the finite difference method, which includes the error and convergence order table of the analytical solution gradient. 1 Norm and L ∞ norm, the travel time gradient calculated by the finite difference method also meets the first-order convergence accuracy. This result supports the use of the "superconvergence" phenomenon in the subsequent discontinuous finite element solver to simplify the solution steps.

[0102] Table 2 Seismic wave travel time gradient obtained by finite difference method

[0103]

[0104] The finite difference method is used to solve the eikonal equation to obtain the first-order accurate seismic wave travel time. The first-order accurate travel time is used to initialize the coefficients of the discontinuous finite element. Then, the local solution formula of the discontinuous finite element is used to update the coefficients of each cube discrete unit to obtain the final seismic wave travel time result.

[0105] Table 3 is the travel time table obtained by the discontinuous finite element method, which contains the error of the analytical solution. It can be seen from Table 3 that the seismic wave travel time finally solved meets the second-order convergence accuracy. In order to make a more intuitive comparison, the travel time contours calculated by the finite difference method and the discontinuous finite element method are compared, as shown in Figure 4-Figure 5 shown. Figure 4 is the travel time contour map in the three-dimensional medium model, Figure 5 is the travel time contour map on the two-dimensional plane at z = 200m. Figure 5 The contour details are magnified in the upper right corner of the figure, and the travel time calculated by the discontinuous finite element method is closer to the analytical solution.

[0106] Table 3 Seismic wave travel time obtained by discontinuous finite element method

[0107]

[0108] Figure 6 Error plot between seismic wave travel times calculated by the finite difference method and the analytical solution; Figure 7 This is the error diagram between the seismic wave travel time calculated by the discontinuous finite element method and the analytical solution. Figure 6 and Figure 7 It can be seen that the error calculated by the discontinuous finite element method is much smaller than the calculation result of the finite difference method, which proves that the seismic wave travel time calculated by the present invention has higher accuracy.

[0109] Example 2

[0110] A computer device includes: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the method for determining travel time of seismic waves in three-dimensional media in Example 1.

[0111] Example 3

[0112] A computer-readable storage medium stores a computer program, which, when executed by a processor, implements the method for determining the travel time of seismic waves in a three-dimensional medium in Example 1.

[0113] Example 4

[0114] A computer program product includes a computer program, which, when executed by a processor, implements the method for determining the travel time of seismic waves in a three-dimensional medium in embodiment 1.

[0115] The technical features of the above embodiments may be combined arbitrarily. To make the description concise, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0116] The principles and implementation methods of the present invention are described in this article using specific examples. The description of the above embodiments is only used to help understand the method and core idea of ​​the present invention. At the same time, for those skilled in the art, according to the idea of ​​the present invention, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting the present invention.

Claims

1. A method for determining travel time of seismic waves in a three-dimensional medium, characterized in that: include: Divide the three-dimensional medium to be measured into a plurality of discrete cubic units; Using the finite difference method and the eikonal equation, the initial value of the seismic wave travel time is determined; According to the initial value of the seismic wave travel time, based on the discontinuous finite element method, the polynomial of each of the cubic discrete units is determined respectively, and the polynomial is: in, is the average travel time of the cube discrete unit, ξ is the travel time gradient of the cube discrete unit in the X direction of the spatial rectangular coordinate system, η is the travel time gradient of the cube discrete unit in the Y direction, γ is the travel time gradient of the cube discrete unit in the Z direction, wherein the spatial rectangular coordinate system is established with the center of the three-dimensional medium to be measured as the origin and the direction from the origin to the earthquake source as the X axis, t is the seismic wave travel time, h is the side length of the cube discrete unit, (x, y, z) is the target position of seismic wave propagation, (x i ,y j , z k ) is the center coordinate of any cubic discrete unit, i is the i-th cubic discrete unit in the X direction, j is the j-th cubic discrete unit in the Y direction, and k is the k-th cubic discrete unit in the Z direction; Updating the coefficients of the multiple sets of polynomials to determine multiple sets of final coefficients specifically includes: Determining a set of travel time equations based on the polynomial and the weak form of the eikonal equation; Based on the travel time equations, updating multiple sets of coefficients of the polynomials to determine multiple sets of final coefficients; The final value of the seismic wave travel time is determined according to the multiple groups of final coefficients.

2. The method for determining travel time of seismic waves in three-dimensional medium according to claim 1, characterized in that: The eikonal equation is: Wherein, V(x,y,z) is the seismic wave propagation velocity of the three-dimensional medium to be measured at the target position (x,y,z) of the seismic wave propagation, and t(x,y,z) is the seismic wave travel time from the earthquake source to the target position (x,y,z) of the seismic wave propagation in the three-dimensional medium to be measured.

3. The method for determining travel time of seismic waves in three-dimensional medium according to claim 1, characterized in that: The expressions of the coefficients of the polynomials are: in, is the travel time from the earthquake source to the first node of any cubic discrete unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the second node of any discrete cubic unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the third node of any discrete cubic unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the fourth node of any discrete cubic unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the fifth node of any cube discrete unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the sixth node of any discrete cubic unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the seventh node of any discrete cubic unit in the three-dimensional medium to be measured, is the travel time from the earthquake source to the eighth node of any cube discrete unit in the three-dimensional medium to be measured, and the first node, the second node, the third node, the fourth node, the fifth node, the sixth node, the seventh node and the eighth node are the eight vertices of the cube discrete unit.

4. The method for determining travel time of seismic waves in three-dimensional medium according to claim 1, characterized in that: Based on the travel time equations, multiple sets of coefficients of the polynomials are updated to determine multiple sets of final coefficients, specifically including: Based on the travel time equation group, the coefficients of the multiple sets of polynomials are updated, and when the coefficients no longer change, the current coefficients of the multiple sets of polynomials are determined to be final coefficients.

5. The method for determining travel time of seismic waves in three-dimensional medium according to claim 1, characterized in that: The travel time equations are: in: The coordinates from the earthquake source to the center of the three-dimensional medium to be measured are (x i ,y j , z k ) is the average travel time of the discrete cube unit, ξ i,j,k The center coordinates are (x i ,y j , z k ) in the X direction, η i,j,k The center coordinates are (x i ,y j , z k ) in the Y direction, γ i,j,k The center coordinates are (x i ,y j , z k ) is the travel time gradient of the cube discrete unit in the Z direction; h is the side length of the cube discrete unit, R 1i,j,k is the first simplified formula, R 2i,j,k is the second simplified formula, R 3i,j,k is the third simplified formula, R 4i,j,k is the fourth simplified formula, A 1i,j,k is the fifth simplified formula, A 2i,j,k is the sixth simplified formula, A 3i,j,k is the seventh simplified formula, A 4i,j,k is the eighth simplified formula, A 5i,j,k is the ninth simplified formula, A 6i,j,k is the tenth simplified formula, A 7i,j,k is the eleventh simplified formula; The coordinates from the earthquake source to the center of the three-dimensional medium to be measured are (x i-1 ,y j , z k ) is the average travel time of the discrete cube unit, ξ i-1,j,k The center coordinates are (x i-1 ,y j , z k ) in the X direction, η i-1,j,k The center coordinates are (x i-1 ,y j , z k ) is the travel time gradient of the cubic discrete unit in the Y direction, γ i-1,j,k The center coordinates are (x i-1 ,y j , z k )’s travel time gradient in the Z direction of the cube discrete unit; The coordinates from the earthquake source to the center of the three-dimensional medium to be measured are (x i ,y j+1 , z k ) is the average travel time of the cube discrete unit, ξ i,j+1,k The center coordinates are (x i ,y j+1 , z k ) in the X direction, η i,j+1,k The center coordinates are (x i ,y j+1 , z k ) is the travel time gradient of the cubic discrete unit in the Y direction, γ i,j+1,k The center coordinates are (x i ,y j+1 , z k )’s travel time gradient in the Z direction of the cube discrete unit; The coordinates from the earthquake source to the center of the three-dimensional medium to be measured are (x i ,y j-1 , z k ) is the average travel time of the cube discrete unit, ξ i,j-1,k The center coordinates are (x i ,y j-1 , z k ) in the X direction, η i,j-1,k The center coordinates are (x i ,y j-1 , z k ) is the travel time gradient of the cubic discrete unit in the Y direction, γ i,j-1,k The center coordinates are (x i ,y j-1 , z k )’s travel time gradient in the Z direction of the cube discrete unit; The coordinates from the earthquake source to the center of the three-dimensional medium to be measured are (x i ,y j , z k+1 ) is the average travel time of the discrete cube unit, ξ i,j,k+1 The center coordinates are (x i ,y j , z k+1 ) in the X direction, η i,j,k+1 The center coordinates are (x i ,y j , z k+1 ) in the Y direction, γ i,j,k+1 The center coordinates are (x i ,y j , z k+1 )’s travel time gradient in the Z direction of the cube discrete unit; The coordinates from the earthquake source to the center of the three-dimensional medium to be measured are (x i ,y j , z k-1 ) is the average travel time of the discrete cube unit, ξ i,j,k-1 The center coordinates are (x i ,y j , z k-1 ) in the X direction, η i,j,k-1 The center coordinates are (x i ,y j , z k-1 ) in the Y direction, γ i,j,k-1 The center coordinates are (x i ,y j , z k-1 )’s cubic discrete unit in the Z direction; α r ,α l ,α f ,α b ,α n ,α s are the flux coefficients of the six faces of the cubic discrete unit; 6. A computer device comprising: A memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the method for determining the travel time of seismic waves in three-dimensional media as described in any one of claims 1-5.

7. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the method for determining the travel time of seismic waves in three-dimensional media as described in any one of claims 1 to 5 is implemented.

8. A computer program product, comprising a computer program, characterized in that When the computer program is executed by a processor, the method for determining the travel time of seismic waves in three-dimensional media as described in any one of claims 1 to 5 is implemented.

Images