High-precision finite difference seismic travel time calculation method under complex terrain condition

Abstract

The invention discloses a rapid iteration seismic travel time calculation method under a complex terrain condition, and relates to the field of seismic travel time calculation, and the method comprises the steps: firstly reading a geological model and related parameter information; then attribute declaration is carried out on the whole calculation grid points, and corresponding parameters are set; initializing a seismic source, and setting seismic source point attributes in calculation grid points according to provided seismic source point positions; then selecting all adjacent points around the seismic source point to form an initial activity list, and carrying out attribute setting; updating the travel time of the activity list; the travel time of the grid nodes is obtained through calculation according to an iterative calculation formula until elements in the activity list are empty, travel time values of all the grid nodes are obtained, and finally the calculation result is output. According to the rapid iterative seismic wave travel time calculation method under the complex terrain condition, high precision is guaranteed, meanwhile, the algorithm has high calculation efficiency, and the calculation efficiency is high. And the calculation efficiency of subsequent seismic exploration data processing is improved.

Description

Technical Field

[0001] The invention relates to the field of seismic wave travel time calculation, and in particular to a high-precision finite difference seismic wave travel time calculation method under complex terrain conditions. Background Art

[0002] Seismic wave travel time is an important parameter that describes the kinematic characteristics of seismic waves and is widely used in technologies such as migration and imaging. Existing travel time calculation methods are basically based on horizontal surfaces, while today's oil and gas exploration work is mostly carried out in complex surface areas. The travel time calculation under complex surface conditions requires the algorithm to not only have good adaptability to complex surfaces but also have high computational efficiency when dealing with irregular boundaries. Therefore, it is of great significance to study efficient and high-precision seismic wave travel time calculation methods that can adapt to complex models.

[0003] "Progress in Geophysical Exploration" Vol. 30, No. 5, 2007 published "Calculation of Seismic Wave Traveltimes on Undulating Surfaces by the Rapid Wavefront Marching Method" by Sun Zhangqing et al., which introduced the rapid wavefront marching method (FMM) as an accurate and stable method for calculating seismic wave traveltimes. Based on the FMM, a method for calculating seismic wave traveltimes under complex surface conditions was proposed, which achieved good results. The computational complexity of the algorithm is O(Nlog2N), where N represents the total number of calculation grid points.

[0004] Chinese invention patent application 201110397445.4 discloses a method for calculating the travel time of seismic waves in VTI media. The algorithm inputs the underground geological velocity model and observation system parameters, and then calculates the corresponding intermediate results. It mainly uses polynomial solution to obtain the travel time of seismic waves, and has achieved good results through verification of calculation results using different models.

[0005] Chinese invention patent application 201810077621.8 discloses a hybrid two-dimensional seismic wave travel time calculation method. The algorithm reads in the velocity model and related calculation parameters, tracks the ray information of a certain distance from the shot point in different directions, calculates the seismic wave travel time within the ray range using the wavefront construction method, and then calculates the travel time of the remaining grid nodes using the rapid advancement method. The good calculation effect of this method has been demonstrated through numerical simulation.

[0006] From the above examples, it can be seen that the above seismic wave travel time calculation method has good calculation effect for conventional horizontal models, but this type of method relies on sorted data structure, and the calculation efficiency is closely related to the complexity of calculation. When facing complex undulating surface problems, the actuarial efficiency is not high, or due to the limitations of the algorithm itself, it has no good adaptability to complex surface problems. Summary of the invention

[0007] The purpose of the present invention is to provide a high-precision finite-difference seismic wave travel time calculation method under complex terrain conditions. Based on the upwind difference scheme and narrowband technology, an iterative method is used to solve the eikonal equation to obtain the seismic wave travel time. The algorithm does not use the sorting technology in the narrowband to ensure the causal relationship of the calculation. All grid nodes in the active list can be updated at one time during the calculation, so that while ensuring high accuracy, the algorithm also has high computational efficiency.

[0008] To achieve the above object, the present invention provides the following technical solution: a method for calculating the travel time of seismic waves under complex terrain conditions, comprising the following steps:

[0009] Step 1: Read in the velocity model and related parameter information, wherein the parameters include the grid size, grid spacing, source location and complex surface information of the velocity model;

[0010] Step 2: Calculate grid point attribute declaration and parameter definition;

[0011] Step 2.1: Declare the attributes of all grid points, and divide them into fixed nodes, active nodes, and far-away points; fixed nodes refer to grid points whose travel time has been calculated, and their attributes are represented by z0; active nodes refer to grid points that are about to be updated, and their attributes are represented by z1; far-away points refer to grid points that have not yet been calculated, and their attributes are represented by z2;

[0012] Step 2.2: According to the information of the undulating surface, the grid point attribute of all points above the undulating surface in the calculation is z0; the grid point attribute in the active list is z1; and the attributes of the remaining grid points are set to z2;

[0013] Step 2.3: Parameter definition in the algorithm: X represents the calculation area, and the travel time when it does not participate in the update calculation is represented by U(X i,j ) means, then U(X i,j )=1000;t i,j represents a calculation point in the calculation area X, and its travel time during update calculation is g(t i,j ) indicates; L indicates the activity list;

[0014] Step 3: Source initialization;

[0015] Step 3.1: According to the information of the undulating surface, set the grid point attributes of all points above the undulating surface to z0; set the attribute of the earthquake source point to z0, and its travel time U=0; set the attributes of the remaining grid points to z2, and their travel time U=1000;

[0016] Step 3.2: Construct the initial activity list, select all neighboring points around the earthquake source point, change their attributes to z1, and construct the initial activity list based on all neighboring points;

[0017] Step 4: In this method, the travel time of each grid node is obtained by solving the Eikonal equation. The travel time of each grid node refers to the points in the initial activity list excluding the source point. In the two-dimensional case, the seismic wave propagation wavefront satisfies the Eikonal equation:

[0018]

[0019] Where t(x,z) is the travel time function, s(x,z) is the model medium slowness function;

[0020] This method is based on the upwind difference formula to discretize the travel time gradient term of the eikonal equation to obtain the seismic wave travel time, and the following formula is used for calculation:

[0021]

[0022] Where I, J represent the longitudinal and transverse lengths of the calculation model respectively, i=2,…,I-1, j=2,…,J-1, h represents the grid spacing, It represents the calculated travel time value of a point in the calculation area. The calculation formula is as follows:

[0023]

[0024] Combining can get The solution is as follows:

[0025]

[0026] in,

[0027] Step 5: Activity list update calculation;

[0028] Step 5.1: For each point L in the active list i,j Do the following:

[0029] 1) Get the travel time U(L) of each point in the activity list i,j ), and recorded as p;

[0030] 2) Get the calculated value g(L) of all points in the active list i,j ), and record it as q, and replace q with U(L i,j )value;

[0031] Step 5.2: If p and q satisfy: |pq|<e, e is a very small positive number, then for each point L in the active list i,j All neighboring points X nb Do the following:

[0032] 1) If X nb For a point that is not in the active list, get the travel time U(X) of all neighboring pointsnb ), and recorded as p;

[0033] 2) Get the calculated value g(X) of all neighboring points nb ), and recorded as q;

[0034] 3) If p and q satisfy: |pq|>0, replace q with the U(X nb ) value, and X nb Move it into the active list and change the grid point attribute to z1;

[0035] 4) After the above calculations are completed, L i,j Remove from the active list and change the grid point attribute to z0;

[0036] Step 6: Determine whether the activity list is empty, if not, return to step 5 to continue calculation.

[0037] Compared with the prior art, the present invention has the following beneficial effects:

[0038] The present invention provides an accurate and efficient method for calculating seismic wave travel time under complex terrain conditions. Based on the upwind difference scheme and narrowband technology, an iterative method is used to solve the eikonal equation to obtain the seismic wave travel time. The algorithm does not use the sorting technology in the narrowband to ensure the causal relationship of the calculation. All grid nodes in the active list can be updated at one time during the calculation, so that while ensuring high accuracy, the algorithm also has high calculation efficiency.

[0039] For the travel time calculation method under complex near-surface conditions, high requirements are placed on the stability of the method. The present invention has unconditional stability and strong adaptability. In addition, the timeliness requirement is relatively high in actual data processing. The present method has high computational efficiency and can improve the computational efficiency of subsequent seismic exploration data processing. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 This is a flow chart of a method for calculating the travel time of seismic waves through rapid iteration under complex terrain conditions according to the present invention;

[0041] Figure 2 It is a schematic diagram of expanding the activity list of the present invention;

[0042] Figure 3 is a relative error diagram of a uniform model in an embodiment of the present invention;

[0043] Figure 4 It is a contour map of seismic wave travel time under complex terrain conditions in an embodiment of the present invention. DETAILED DESCRIPTION

[0044] 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.

[0045] See also Figure 1-4 ,

[0046] The present invention provides a method for calculating the travel time of seismic waves under complex terrain conditions. Figure 1 As shown, the following steps are included:

[0047] Step 1: Read in the velocity model and related parameter information, wherein the parameters include the grid size, grid spacing, source location and complex surface information of the velocity model;

[0048] In this embodiment, the relative error is calculated using a uniform velocity model, the undulating surface model is a custom complex model, the model size is 12011201, the grid spacing is 5m, and the embodiment data described below is a uniform velocity model, the model size is 20002000, the grid spacing is 1m, the velocity is 2000m / s, and the source point is (1000m, 2m);

[0049] Step 2: Calculate grid point attribute declaration and parameter definition;

[0050] Step 2.1: Declare the attributes of all grid points, and divide them into fixed nodes, active nodes, and far-away points; fixed nodes refer to grid points whose travel time has been calculated, and their attributes are represented by z0; active nodes refer to grid points that are about to be updated, and their attributes are represented by z1; far-away points refer to grid points that have not yet been calculated, and their attributes are represented by z2;

[0051] Step 2.2: According to the information of the undulating surface, the grid point attribute of all points above the undulating surface in the calculation is z0; the grid point attribute in the active list is z1; and the attributes of the remaining grid points are set to z2;

[0052] Step 2.3: Parameter definition in the algorithm: X represents the calculation area, and the travel time when it does not participate in the update calculation is represented by U(X i,j ) means, then U(X i,j )=1000;t i,j represents a calculation point in the calculation area X, and its travel time during update calculation is g(t i,j ) indicates; L indicates the activity list;

[0053] Step 3: Source initialization;

[0054] Step 3.1: According to the information of the undulating surface, set the grid point attributes of all points above the undulating surface to z0; set the attribute of the earthquake source point to z0, and its travel time U=0; set the attributes of the remaining grid points to z2, and their travel time U=1000;

[0055] In this embodiment, the uniform velocity model is used for trial calculation, and the attribute of the source point (2000m, 2000m) is set to z0, and its travel time U=0, and the attributes of the remaining grid points are set to z2, and their travel time U=1000;

[0056] Step 3.2: Construct the initial activity list, select all neighboring points around the earthquake source point, change their attributes to z1, and construct the initial activity list based on all neighboring points;

[0057] In this embodiment, the attributes of all neighboring points around the earthquake source point (1000m, 2m), namely (1000m, 1m), (1000m, 3m), (999m, 2m), and (1001m, 2m), are modified to z1, and an initial activity list is formed based on the four neighboring points;

[0058] Step 4: In this method, the travel time of each grid node is obtained by solving the Eikonal equation. The travel time of each grid node refers to the points in the initial activity list excluding the source point. In the two-dimensional case, the seismic wave propagation wavefront satisfies the Eikonal equation:

[0059]

[0060] Where t(x,z) is the travel time function, s(x,z) is the model medium slowness function;

[0061] This method is based on the upwind difference formula to discretize the travel time gradient term of the eikonal equation to obtain the seismic wave travel time, and the following formula is used for calculation:

[0062]

[0063] Where I, J represent the longitudinal and transverse lengths of the calculation model respectively, i=2,…,I-1, j=2,…,J-1, h represents the grid spacing, It represents the calculated travel time value of a point in the calculation area. The calculation formula is as follows:

[0064]

[0065] Can get The solution is as follows:

[0066]

[0067] in,

[0068] Step 5: Activity list update calculation;

[0069] Step 5.1: For each point L in the active list i,j Do the following:

[0070] 1) Get the travel time U(L) of each point in the activity list i,j ), and recorded as p;

[0071] 2) Get the calculated value g(L) of all points in the active list i,j ), and record it as q, and replace q with U(L i,j )value;

[0072] Step 5.2: If p and q satisfy: |pq|<e, e is a very small positive number, then for each point L in the active list i,j All neighboring points X nb Do the following:

[0073] 1) If X nb For a point that is not in the active list, get the travel time U(X) of all neighboring points nb ), and recorded as p;

[0074] 2) Get the calculated value g(X) of all neighboring points nb ), and recorded as q;

[0075] 3) If p and q satisfy: |pq|>0, replace q with the U(X nb ) value, and X nb Move it into the active list and change the grid point attribute to z1;

[0076] 4) After the above calculations are completed, L i,j Remove from the active list and change the grid point attribute to z0;

[0077] Step 6: Determine whether the activity list is empty, if not, return to step 5 to continue calculation.

[0078] When used specifically, the present invention provides a high-precision finite-difference seismic wave travel time calculation method under complex terrain conditions. For a model under complex near-surface conditions, the analysis and calculation results show that the distribution of travel time contour lines conforms to the propagation law of seismic waves inside the medium. When faced with complex geological structures such as faults in the velocity model, the travel time line diagram clearly reflects the propagation law that seismic waves diffuse outward quickly in the high-speed layer and diffuse outward slowly in the low-speed layer. In addition, the travel time lines are still distributed continuously at the boundary between the high-speed layer and the low-speed layer, which also proves that the algorithm of the present invention has good adaptability to complex surface models.

[0079] Although the present invention has been described in detail with reference to the aforementioned embodiments, it is still possible for those skilled in the art to modify the technical solutions described in the aforementioned embodiments, or to make equivalent substitutions for some of the technical features therein. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A high-precision finite-difference seismic wave travel time calculation method under complex terrain conditions, characterized in that: The following steps are involved: Step 1: Read in the velocity model and related parameter information; Step 2: Calculate grid point attribute declaration and parameter definition; Step 3: Source initialization; Step 4: In this method, the travel time of each grid node is obtained by solving the Eikonal equation. In the two-dimensional case, the seismic wave propagation wavefront satisfies the Eikonal equation: |▽t(x,z)|=s(x,z) Where t(x,z) is the travel time function, s(x,z) is the model medium slowness function; Based on the upwind difference formula, the travel time gradient term of the discretized eikonal equation is used to obtain the seismic wave travel time, which is calculated using the following formula: Where I, J represent the longitudinal and transverse lengths of the calculation model respectively, i=2,…,I-1, j=2,…,J-1, h represents the grid spacing, It represents the calculated travel time value of a point in the calculation area. The calculation formula is as follows: Step 5: Activity list update calculation; The list update calculation steps are as follows: Step 5.1: For each point L in the active list i,j Do the following: 1) Get the travel time U(L) of each point in the activity list i,j ), and recorded as p; 2) Get the calculated value g(L) of all points in the active list i,j ), and record it as q, and replace q with U(L i,j )value; Step 5.2: If p and q satisfy: |pq|<e, then for each point L in the active list i,j All neighboring points X nb Do the following: 1) If X nb For a point that is not in the active list, get the travel time U(X) of all neighboring points nb ), and recorded as p; 2) Get the calculated value g(X) of all neighboring points nb ), and recorded as q; 3) If p and q satisfy: |pq|>0, replace q with the U(X nb ) value, and X nb Move it into the active list and change the grid point attribute to z1; 4) After the above calculations are completed, L i,j Remove from the active list and change the grid point attribute to z0; Step 6: Determine whether the activity list is empty, if not, return to step 5 to continue calculation.

2. The high-precision finite-difference seismic wave travel time calculation method under complex terrain conditions according to claim 1, characterized in that: In step 1, the parameters include the grid size of the velocity model, the grid spacing, the earthquake source location and the complex surface information.

3. The high-precision finite-difference seismic wave travel time calculation method under complex terrain conditions according to claim 1, characterized in that: The calculation of grid point attribute declaration and parameter definition in step 2 includes the following steps: Step 2.1: Declare the attributes of all grid points, and divide them into fixed nodes, active nodes, and far-away points; fixed nodes refer to grid points whose travel time has been calculated, and their attributes are represented by z0; active nodes refer to grid points that are about to be updated, and their attributes are represented by z1; far-away points refer to grid points that have not yet been calculated, and their attributes are represented by z2; Step 2.2: According to the information of the undulating surface, the grid point attribute of all points above the undulating surface in the calculation is z0; the grid point attribute in the active list is z1; and the attributes of the remaining grid points are set to z2; Step 2.3: Parameter definition in the algorithm: X represents the calculation area, and the travel time when it does not participate in the update calculation is represented by U(X i,j ) means, then U(X i,j )=1000;t i,j represents a calculation point in the calculation area X, and its travel time during update calculation is g(t i,j ) indicates; L indicates the activity list.

4. The high-precision finite-difference seismic wave travel time calculation method under complex terrain conditions according to claim 1, characterized in that: The steps of the source initialization in step 3 are as follows: Step 3.1: According to the information of the undulating surface, set the grid point attributes of all points above the undulating surface to z0; set the attribute of the earthquake source point to z0, and its travel time U=0; set the attributes of the remaining grid points to z2, and their travel time U=1000; Step 3.2: Construct the initial activity list. Select all neighboring points around the earthquake source point, change their attributes to z1, and construct the initial activity list based on all neighboring points.

5. The high-precision finite-difference seismic wave travel time calculation method under complex terrain conditions according to claim 1, characterized in that: In step 4, according to the formula as well as Can get The solution is as follows: in, 6. The high-precision finite-difference seismic wave travel time calculation method under complex terrain conditions according to claim 1, characterized in that: Each grid node in step 4 refers to a point in the initial activity list excluding the earthquake source point.

7. The high-precision finite-difference seismic wave travel time calculation method under complex terrain conditions according to claim 1, characterized in that: In the step 5.2, |pq|<e, and e is a very small positive number.

Images