First arrival travel time tomography inversion method based on double-difference constraint

Abstract

The invention discloses a first arrival travel time tomography inversion method based on double-difference constraint, and the method comprises the following steps: S1, building a travel time-based first arrival ray tomography inversion equation based on a high-frequency approximation theory and a perturbation theory; s2, introducing double-difference constraint when solving the travel-time-based first-arrival ray tomography inversion equation. According to the method, a new double-difference constraint condition is added to calculation of a tomography inversion equation set, the constraint enables the change of the velocity field within a certain range to be more accurate by minimizing the travel time residual error of the adjacent detection points, the travel time accumulative error of the adjacent detection points is eliminated to a certain extent, and the inversion precision of the near-surface velocity field is improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of seismic data processing for oil and gas exploration, and particularly relates to a first-arrival travel-time tomography inversion method based on double-difference constraint. Background Art

[0002] At the present stage, in the field of oil and gas exploration, the main methods for obtaining a more accurate underground velocity model are travel-time tomography inversion and waveform inversion, etc. The former mainly includes ray-based travel-time tomography inversion and travel-time tomography inversion based on the wave equation. In practical applications, ray-based tomography inversion has shown extremely strong practical value due to its high computational efficiency and strong robustness. And due to the advantages of strong energy and easy picking of travel time of the first arrival wave, the ray tomography inversion method based on the first arrival has been widely used in the field of constructing the near-surface velocity model. However, in practical applications, it is often affected by defects such as low ray tracing accuracy and strong ill-conditioning of the tomography equation set, resulting in low accuracy of near-surface velocity inversion and affecting the final inversion velocity accuracy.

[0003] The conventional first-arrival ray tomography inversion method updates the velocity using absolute travel time, aiming to minimize the error between the observed travel time and the calculated travel time. However, it will be affected by factors such as travel time measurement error and cumulative travel time calculation error. For example, absolute travel time measurement is affected by factors such as the first arrival jump point selection criterion and noise interference, resulting in the absolute travel time update criterion reducing the accuracy and reliability of the inversion result. On the contrary, replacing "absolute time measurement" with "relative time measurement", the measured seismic travel time is also called double-difference travel time, which can weaken the error brought by the absolute travel time update criterion to a certain extent and improve the inversion accuracy of first-arrival tomography.

[0004] Therefore, aiming at the problem of absolute travel time in first-arrival tomography inversion, the present invention designs a first-arrival tomography inversion method with double-difference travel time constraint. Summary of the Invention

[0005] The problem to be solved by the present invention is to provide a first-arrival travel-time tomography inversion method based on double-difference constraint. This method introduces the double-difference idea into the field of seismic exploration for oil and gas, and improves the inversion accuracy of first-arrival ray tomography by constraining the travel time of adjacent underground rays.

[0006] To solve the above technical problems, the technical solution adopted by the present invention is: a first-arrival travel-time tomography inversion method based on double-difference constraint, including the following steps,

[0007] S1: Based on the high-frequency approximation theory and perturbation theory, establish a first-arrival ray tomography inversion equation based on travel time;

[0008] S2: Introduce double-difference constraint when solving the first-arrival ray tomography inversion equation based on travel time.

[0009] Further, S1 includes the following steps:

[0010] S11: Based on the high-frequency approximation theory, seismic waves propagate along the ray path, and the required travel time is:

[0011] t = ∫ L(s) s(x,z)dl (Formula 1)

[0012] In the formula, s(x,z) is the slowness field, which is the reciprocal of the velocity field v(x,z), x and z represent spatial coordinates, and L(s) is the ray path under the slowness field s(x,z);

[0013] S12: When there is a perturbation relationship as shown in Formula 2 in the slowness field, a small perturbation δL occurs in the ray path, and at this time the travel time perturbation is δt, as shown in Formulas 3 and 4:

[0014] s = s0 + δs + O(δs 2 ) (Formula 2)

[0015] In the formula, δs is the first-order perturbation quantity, and O(δs 2 ) is the high-order perturbation small quantity.

[0016]

[0017] t = t0 + δt (Formula 4)

[0018] Formula 3 is expanded to obtain the following Formula 5:

[0019]

[0020] S13: Combine Formula 4 and Formula 5, and omit the high-order small quantity O(δs 2 ), to obtain the following Formula 6:

[0021]

[0022] S14: Under discrete conditions, combine the tomography equations corresponding to all rays during the propagation process to obtain the first-arrival ray tomography inversion equation based on travel time.

[0023] Further, the first-arrival ray tomography inversion equation based on travel time is as follows:

[0024] LΔs = Δt (Formula 7)

[0025] In the formula, L represents the ray path matrix, which is a large-scale sparse matrix, the number of rows is the number of rays, the number of columns is the total number of discrete grids of the underground velocity (slowness) field, and each element in the matrix is the ray path length of the corresponding ray at the corresponding discrete grid.

[0026] Furthermore, the slowness update Δs of the first arrival ray tomography inversion equation based on travel time is obtained by the LSQR algorithm.

[0027] Furthermore, the S2 includes the following steps:

[0028] S21: When minimizing the absolute travel time residual at the geophone point in the first arrival travel time tomography, the following conditions need to be satisfied:

[0029]

[0030] In the formula, t obs,i represents the observed travel time of the i-th channel, and t cal,i represents the calculated travel time of the i-th channel;

[0031] S22: For adjacent i and i+1 channels corresponding to the same source, the formula is as follows:

[0032]

[0033] In the formula, and are the observed travel time difference and the calculated travel time difference of adjacent i and i+1 channels corresponding to the same source, respectively;

[0034] S23: The double difference between and

[0035]

[0036] needs to satisfy the following conditions:

[0037]

[0038] S24: Perform double difference constrained tomography inversion.

[0039] Furthermore, the double difference constrained tomography inversion is to solve a constrained problem that satisfies the following conditions:

[0040]

[0041] In the formula, L represents the ray path matrix of all rays, and L i+1 (L i ) represents the distribution of the ray path length corresponding to the geophone points of the i+1(i)-th channel in the underground discrete grid, and s cal represents the updated slowness.

[0042] Furthermore, in the S2, the double difference travel time is defined by the channel spacing of adjacent 5 channels.

[0043] Furthermore, the present invention provides a device for running the above data processing method.

[0044] Furthermore, the present invention provides a device, including a memory, a processor, and an algorithm stored in the memory and executable on the processor. When the processor executes the computer program, the above data processing method is implemented.

[0045] Furthermore, the present invention provides a computer-readable storage medium storing a computer algorithm, which implements the above data processing when executed by a processor.

[0046] The advantages and positive effects of the present invention are as follows:

[0047] The present invention adds a new double-difference constraint condition to the calculation of the tomographic inversion equations. By minimizing the travel-time residuals of adjacent geophones, this constraint makes the variation of the velocity field within a certain range more accurate, eliminates the travel-time cumulative error of adjacent geophones to a certain extent, and improves the inversion accuracy of the near-surface velocity field. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 is a schematic diagram of the difference idea of an embodiment of the present invention.

[0049] Figure 2 is a schematic diagram of the overall flow of an embodiment of the present invention.

[0050] Figure 3 is a true velocity model diagram of a specific embodiment of the present invention.

[0051] Figure 4 is an initial gradient model diagram of a specific embodiment of the present invention.

[0052] Figure 5 is a schematic diagram of the inversion result of the traditional method in a specific embodiment of the present invention.

[0053] Figure 6 is a schematic diagram of the inversion result of the method of the present invention in a specific embodiment of the present invention.

[0054] Figure 7 is a comparison diagram of the velocity curves extracted at x = 8000 m from the inversion results of the present invention and the traditional method in a specific embodiment of the present invention.

[0055] Figure 8 is a ray density diagram of tomographic inversion of the traditional method in a specific embodiment of the present invention.

[0056] Figure 9 is a ray density diagram of tomographic inversion of the method of the present invention in a specific embodiment of the present invention.

[0057] Figure 10It is the true velocity model diagram of another specific embodiment of the present invention.

[0058] Figure 11 It is the initial gradient model diagram of another specific embodiment of the present invention.

[0059] Figure 12 It is the schematic diagram of the inversion result of the traditional method in another specific embodiment of the present invention.

[0060] Figure 13 It is the schematic diagram of the inversion result of the method in another specific embodiment of the present invention.

[0061] Figure 14 It is the comparison diagram of the velocity curves extracted at x = 6000 m from the inversion results of the present invention and the traditional method in another specific embodiment of the present invention.

[0062] Figure 15 It is the comparison diagram of the velocity curves extracted at z = 600 m from the inversion results of the present invention and the traditional method in another specific embodiment of the present invention. Specific Embodiment

[0063] Next, the technical solution of the present invention will be clearly and completely described in conjunction with the accompanying drawings. Obviously, the described embodiments are part of the embodiments of the present invention, rather than all of them. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative work shall fall within the protection scope of the present invention.

[0064] Aiming at the error caused by the conventional travel-time tomography inversion using the minimization of the absolute travel-time residual at the geophone as the criterion for velocity update, the present invention adds a new double-difference constraint condition to the calculation of the tomography inversion equations. This constraint minimizes the travel-time residuals of adjacent geophones, making the variation of the velocity field within a certain range more accurate, and can also eliminate the travel-time cumulative error of adjacent geophones to a certain extent, improving the inversion accuracy of the near-surface velocity field.

[0065] The following further describes the embodiments of the present invention in conjunction with the accompanying drawings:

[0066] As Figure 1 , Figure 2 shown, a first-arrival travel-time tomography inversion method based on double-difference constraint includes the following steps.

[0067] S1: Based on the high-frequency approximation theory and perturbation theory, establish the first-arrival ray tomography inversion equation based on travel time, and the specific steps are as follows.

[0068] Based on the high-frequency approximation theory, seismic waves can be regarded as propagating along an infinitely thin ray path, and the travel time required to propagate along the ray path is:

[0069] t = ∫L(s) s(x,z)dl (Formula 1)

[0070] In the formula, s(x,z) is the slowness field, which is the reciprocal of the velocity field v(x,z). x and z represent spatial coordinates, and L(s) is the ray path under the slowness field s(x,z).

[0071] When there is a perturbation relationship as shown in Formula 2 in the slowness field, a small perturbation δL will occur in the ray path. At this time, the travel time perturbation is δt. Specifically, as shown in Formulas 3 and 4,

[0072] s = s0 + δs + O(δs 2 ) (Formula 2)

[0073] In the formula, δs is the first-order perturbation quantity, and O(δs 2 ) is the high-order perturbation small quantity.

[0074]

[0075] t = t0 + δt (Formula 4)

[0076] Expanding Formula 3 gives the following Formula 5.

[0077]

[0078] The second integral on the right side of Formula 5 can be regarded as zero because Fermat's principle states that for fixed endpoints, the travel time along the ray path is stationary with respect to perturbations on the path, that is,

[0079] Combining Formula 4 and Formula 5 and omitting the high-order small quantity O(δs 2 ), the following Formula 6 is obtained.

[0080]

[0081] Under discrete conditions, there is a ray for each source and the corresponding geophone, and each ray has a corresponding tomography equation (Formula 6). Combining the tomography equations corresponding to all rays during the propagation process can form a large linear equation system, which is denoted in matrix form, thereby obtaining the first-arrival ray tomography inversion equation based on travel time. Specifically, the first-arrival ray tomography inversion equation based on travel time is as follows.

[0082] LΔs = Δt (Formula 7)

[0083] In the formula, L represents the ray path matrix, which is a large-scale sparse matrix. The number of rows is the number of rays, and the number of columns is the total number of discrete grids of the underground velocity (slowness) field. Each element in the matrix is the ray path length of the corresponding ray at the corresponding discrete grid.

[0084] By solving formula 7 using the LSQR algorithm, the slowness update Δs can be obtained, and then the velocity field can be updated.

[0085] S2: Introduce double-difference constraints when solving the first-arrival ray tomography inversion equation based on travel time.

[0086] The absolute travel time residual as the criterion for velocity update reduces the inversion accuracy and reliability. To solve the above problems and better constrain the velocity update amount of adjacent rays in the subsurface medium, the present invention proposes a new double-difference constraint condition and adds it to the calculation of the tomography inversion equation system. This constraint makes the change of the velocity field more accurate within a certain range by minimizing the travel time residual of adjacent geophones, and can also eliminate the travel time cumulative error of adjacent geophones to a certain extent, improving the inversion accuracy of the large-offset velocity field.

[0087] For convenience of description, the following discussion is carried out with n geophone points of the same seismic source. Usually, the first-arrival travel time tomography tries to minimize the absolute travel time residual at the geophone points. Ideally, the following conditions need to be satisfied:

[0088]

[0089] In the formula, t obs,i represents the observed travel time of the i-th channel, and t cal,i represents the calculated travel time of the i-th channel. As shown in the schematic diagram of the difference idea in Figure 1 , where the shot point represents the seismic source point, the black dot i and the gray dot i + 1 represent two adjacent channels of the same seismic source. The black and gray curves represent the ray paths corresponding to two adjacent channels of the same seismic source, and the arrows on them represent the propagation direction of the rays.

[0090] In this embodiment, for the adjacent i-th and (i + 1)-th channels corresponding to the same seismic source, the formula is as follows.

[0091]

[0092] In the formula, and are the observed travel time difference and the calculated travel time difference of the adjacent i-th and (i + 1)-th channels of the same seismic source respectively. When the velocity model is accurate, is a definite value. When updating the velocity model, if can be made very close to , then the relative travel time change of the ray at the adjacent geophones is more accurate. And and The calculation method can also reduce the cumulative calculation travel time error caused by the increase of the offset to a certain extent, thereby improving the inversion accuracy and reliability.

[0093] The double difference between is expressed in the following form:

[0094]

[0095] Ideally, the following conditions should be met:

[0096]

[0097] When using the LSQR method to solve the tomography equation, the single-difference constraint condition formula 8 makes the absolute travel-time residual approach 0, thus making the velocity field closer to the true velocity; while the double-difference constraint condition formula 11 makes the observed differential travel-time residual between adjacent geophones closer to the true differential travel-time residual, which can constrain the velocity change between adjacent rays, thereby reducing the error of the updated velocity field and reducing the ill-conditioning of the equation solution.

[0098] The double-difference constrained tomography inversion is to solve a constrained problem that satisfies the following conditions:

[0099]

[0100]

[0101] In the formula, L represents the ray path matrix of all rays, and L i+1 (L i ) represents the distribution of the ray path lengths corresponding to the geophones of the (i + 1)(i)-th channels in the underground discrete grid, and s cal represents the updated slowness.

[0102] When solving the above problem, the present invention does not add the double-difference constraint formula 13 as a penalty term to the original objective function. The reason for this treatment is that the introduction of the penalty term will make the entire objective function more complex, and it is quite difficult to determine the weight of the penalty factor. In view of the above considerations, the present invention adopts a method for dealing with the constrained problem without modifying the objective function. That is, in the normal iterative optimization process, after each iteration ends, Δs is "pulled back" to the constraint condition to ensure convergence while satisfying the constraint condition formula 13. In other words, when performing parameter update, the slowness updated each time is restricted so that the updated local slowness remains within a specific range. This method avoids the possible errors in the update process to a certain extent and improves the inversion accuracy of the absolute travel-time residual update speed. Since the double-difference travel time between adjacent 1 geophones is approximately 0, the present invention uses the trace distance of adjacent 5 channels to define the double-difference travel time in practical applications.

[0103] The following combines specific embodiments to specifically elaborate on the present invention:

[0104] Embodiment 1

[0105] First, the effectiveness of the method is verified using Example 1 - a simple multi - undulating layered velocity model as shown in Figure 3 . The size of the simple multi - layer velocity model is 801 * 201, with a grid spacing of 20m both horizontally and vertically. A total of 100 sources and 400 receivers are evenly distributed on the surface. The initial model is set as a constant - velocity gradient velocity model, as shown in Figure 4 . It can be seen that there is a relatively large difference between the initial gradient model and the true model. Then, the traditional tomography inversion method and the tomography inversion method of this method are used to update the initial velocity model. Among them, the traditional tomography inversion method refers to the traditional inversion method without adding double - difference constraints, and this method refers to the inversion method with double - difference constraints added. After completing the update of the initial velocity model, the accuracy of the velocity model is evaluated mainly in two aspects. Macroscopically, the structural form of the velocity field is mainly observed, and numerically, the numerical matching degree between the inversion velocity field and the true velocity field is extracted for a certain trace.

[0106] Figure 5 and Figure 6 respectively show the tomography inversion results using the traditional method and this method. By comparing the tomography inversion results of different methods, it can be qualitatively observed that both methods can basically invert the medium - and low - wavenumber information of the model. However, relatively speaking, the velocity structure of the tomography result obtained by this method is more accurate, with a smaller difference from the true velocity model and is closer to the true velocity model.

[0107] The velocity values at x = 8000m are extracted from the tomography inversion results of different methods and compared with the true velocity field, as shown in Figure 7 . By observing the tomography velocity curves obtained by different methods, since the initial velocity field is very different from the true velocity field, the traditional method has updated to a certain extent in the shallow layer, but there is still a large distance from the true model. The inversion result of this method shows a more significant improvement compared with the traditional method, being closer to the true value numerically and having a better inversion effect.

[0108] Figure 8 and Figure 9 are respectively the ray density of the tomography inversion by the traditional method and the ray density of the tomography inversion by this method. It can be seen that the deepest ray paths of both methods are at a depth of about 3000m, and the ray path distributions are relatively similar. However, the ray path distribution of this method fits the layered structure distribution of the actual velocity model better than that of the traditional method, and the ray density in the shallow layer is also more concentrated, resulting in a better inversion effect for the shallow layer.

[0109] Example 2

[0110] To further verify the effectiveness and processing effect of the method of the present invention, use as shown in Figure 10Example 2 shown - The Marmousi model is used to test this tomography inversion method. The size of the Marmousi model is 500 * 174, the grid spacing in the horizontal and vertical directions is 20m, and a total of 100 sources and 250 detectors are set. The initial velocity field is a constant velocity gradient velocity model, as Figure 11 shown.

[0111] The initial velocity field is tomographically inverted using the traditional method and this method, and the obtained results are as Figure 12 and Figure 13 shown. By qualitatively observing the tomographic inversion results of different methods, it can be seen that both methods can invert the general shape of the velocity in the area from the near-surface depth of 500m to the depth of 1500m of the Marmousi model. However, relatively speaking, 13 at the depth from 800m to 1000m and the horizontal position from 6000m to 9000m, the inversion of the horizon shape is more accurate and has a clearer inversion effect.

[0112] From the tomographic inversion results Figure 12 and Figure 13 of different methods, the velocity values at x = 6000m are extracted and compared with the true velocity field for analysis, as Figure 14 shown. By observing the tomographic velocity curves obtained by different methods, both methods generally invert the general trend of the velocity change, but this method is closer to the true model at the depth from 500m to 1000m. To further verify the effect of this method, the velocity curve at z = 600m is also extracted for comparison, and the results are as Figure 15 shown. At the depth of 600m, both methods invert the general trend of the velocity change, but relatively speaking, the tomographic inversion method with double-difference constraint is numerically better than the traditional method and is closer to the true model.

[0113] The advantages and positive effects of the present invention are:

[0114] The present invention adds a new double-difference constraint condition to the calculation of the tomography inversion equation set. This constraint makes the change of the velocity field more accurate within a certain range by minimizing the travel-time residuals of adjacent detectors, eliminates the travel-time cumulative error of adjacent detectors to a certain extent, and improves the inversion accuracy of the near-surface velocity field.

[0115] The above has described in detail an embodiment of the present invention, but the content described is only a preferred embodiment of the present invention and cannot be considered as limiting the scope of implementation of the present invention. All equivalent changes and improvements made according to the scope of the application of the present invention should still fall within the scope covered by the patent of the present invention.

Claims

1. A first arrival travel time tomographic inversion method based on double difference constraints, characterized in that: Including the following steps, S1: Based on the high-frequency approximation theory and the perturbation theory, establish the first-arrival ray tomography inversion equation based on travel time; S2: Introduce double-difference constraints when solving the first-arrival ray tomography inversion equation based on travel time.

2. The first arrival travel time tomographic inversion method based on double difference constraint according to claim 1, characterized in that: The S1 includes the following steps, S11: Based on the high-frequency approximation theory, seismic waves propagate along the ray path, and the required travel time is: t = ∫ L(s) s(x,z) dl (Formula 1) In the formula, s(x,z) is the slowness field, which is the reciprocal of the velocity field v(x,z), x and z represent spatial coordinates, and L(s) is the ray path under the slowness field s(x,z); S12: When there is a perturbation relationship as shown in Formula 2 in the slowness field, the ray path generates a small perturbation δL, and at this time the travel time perturbation is δt, as shown in Formula 3 and Formula 4, s = s0 + δs + O(δs 2 ) (Equation 2) where δs is the first-order perturbation quantity, and O(δs 2 ) is the small quantity of high-order perturbation, t = t0 + δt (Formula 4) Formula 3 is expanded to obtain the following Formula 5, S13: Combine Equation 4 and Equation 5, and neglect the higher-order small quantity O(δs 2 ), to obtain the following Equation 6, S14: Under discrete conditions, combine the tomography equations corresponding to all rays during the propagation process to obtain the first-arrival ray tomography inversion equation based on travel time.

3. A first arrival travel time tomographic inversion method based on double difference constraints according to claim 1 or 2, characterized in that: The first-arrival ray tomography inversion equation based on travel time is as follows, LΔs = Δt (Formula 7) In the formula, L represents the ray path matrix, which is a large-scale sparse matrix, the number of rows is the number of rays, the number of columns is the total number of discrete grids of the underground velocity (slowness) field, and each element in the matrix is the ray path length of the corresponding ray at the corresponding discrete grid.

4. A first arrival travel time tomography inversion method based on double difference constraint according to claim 3, characterized in that: The first-arrival ray tomography inversion equation based on travel time obtains the slowness update amount Δs through the LSQR algorithm.

5. A first arrival travel time tomography inversion method based on double difference constraint according to claim 1 or 2, characterized in that: The S2 includes the following steps, S21: When minimizing the absolute travel time residual at the geophone point in the first-arrival travel time tomography, the following conditions need to be met: where t obs,i represents the observed travel time of the i-th trace, and t cal,i represents the calculated travel time of the i-th trace; S22: Define for adjacent channels i and i + 1 corresponding to the same source, the formula is as follows, In the formula, and are respectively the observed travel time difference and the calculated travel time difference between adjacent i-th and (i + 1)-th channels of the same seismic source; S23: The double difference between is expressed in the following form: The following conditions need to be met: S24: Perform double-difference constrained tomography inversion.

6. A first arrival travel time tomographic inversion method based on double difference constraint according to claim 1 or 2, characterized in that: The double-difference constrained tomography inversion is to solve a constrained problem that meets the following conditions: where L represents the ray path matrix of all rays, L i+1 (L i ) represents the distribution of the ray path lengths corresponding to the i + 1 (i) geophone points in the underground discrete grid, s cal represents the updated slowness.

7. A first arrival travel time tomographic inversion method based on double difference constraint according to claim 1 or 2, characterized in that: In the S2, the double-difference travel time is defined by the channel spacing of adjacent 5 channels.

8. A device, characterized in that: Run the data processing method according to any one of claims 1 to 7.

9. A device, comprising a memory, a processor, and an algorithm stored in the memory and executable on the processor, characterized in that: When the processor executes the computer program, it implements the data processing method according to any one of claims 1 to 7.

10. A computer-readable storage medium storing a computer algorithm, characterized in that, When the computer algorithm is executed by the processor, it implements the data processing according to any one of claims 1 to 7.

Images