First-arrival wave reverse time migration method and system based on excitation amplitude

Abstract

The present invention discloses a method and system for reverse time migration of first-arrival waves based on excitation amplitude, including obtaining seismic wavefield calculation parameters, forward extension of the source wavefield, identification of the first-arrival waves in the source wavefield, and preservation of the maximum amplitude and corresponding time of the first-arrival waves at each spatial grid node; reverse time extension of the detection point wavefield; and cross-correlating the maximum amplitude of the source wavefield with the detection point wavefield at the same time according to the excitation time at each spatial grid node, and superimposing the cross-correlation values at all times to obtain the reverse time migration result. By preserving the maximum amplitude and time of the first-arrival waves during the propagation of the source wavefield and then performing cross-correlation imaging with the detection point wavefield, the present invention can reduce the computational complexity and interference from multipath waves while ensuring the reverse time migration imaging capability, thereby improving the signal-to-noise ratio.

Description

Technical Field

[0001] The present invention relates to the technical field of earthquake source wave reverse time migration, in particular to a first arrival wave reverse time migration method based on excitation amplitude and a system thereof. Background Art

[0002] In the field of geophysics, the reverse time migration method is currently the most accurate and powerful migration method. Its implementation process is divided into three steps: (1) forward extension of the source wavefield, (2) reverse time extension of the receiver wavefield, and (3) using imaging conditions (such as cross-correlation imaging of the source wavefield and the receiver wavefield). This process requires the source wavefield value at each moment to be stored in the computer, which requires a large amount of computer memory and high computational cost. Compared with the cross-correlation imaging condition, the excitation amplitude imaging condition only needs to save the maximum amplitude value and corresponding time on each computational grid node during the forward extension of the source wavefield, thus greatly saving computer memory and improving computational efficiency. However, since the excitation amplitude imaging condition only saves one maximum amplitude and time on each grid node, it cannot image multipath waves, thereby limiting the imaging capability. This is considered a disadvantage of the excitation amplitude imaging condition. A solution strategy is to save several maximum amplitude values and their corresponding times, saving a certain amount of multipath waves for imaging. However, imaging multipath waves is a complex problem. Simply preserving the multipath waves and then using cross-correlation imaging is not enough. Targeted processing is required to achieve accurate imaging; otherwise, the resulting image will mostly be noise. In fact, the excitation amplitude imaging condition does not exclude multipath waves. In complex media, the amplitude of the multipath waves can sometimes be stronger than that of the first arrival. Therefore, when using the excitation amplitude imaging condition, the amplitudes and arrival times of the multipath waves are preserved and imaged.

[0003] In reverse time migration, whether using cross-correlation imaging or excitation amplitude imaging, the primary imaging method for subsurface media is the first arrival. For complex models, existing excitation amplitude reverse time migration methods preserve the excitation amplitude of the source wavefield, which includes both the first arrival and multipath wavefield amplitudes. For multipath wavefields, direct cross-correlation with the receiver wavefields not only fails to image the subsurface, but also generates migration noise, compromising imaging quality. Therefore, a first arrival reverse time migration method and system based on excitation amplitude are needed. Summary of the Invention

[0004] The present invention aims to provide a first-arrival wave reverse time migration method based on excitation amplitude and a system thereof.

[0005] The present invention specifically comprises the following steps:

[0006] Obtain the source wavefield calculation parameters, forward extend the source wavefield, identify the first arrival wave of the source wavefield, and retain the maximum amplitude and corresponding time of the first arrival wave at each spatial grid node;

[0007] Reverse time extension of the wave field at the detection point;

[0008] At each spatial grid node, the maximum amplitude of the source wavefield is cross-correlated with the wavefield of the receiver point at the same time according to the excitation time, and the cross-correlation values at all times are superimposed to obtain the reverse time migration result.

[0009] Furthermore, the calculation parameters of the earthquake source wave field are obtained, including the grid size, the number of grid nodes, the seismic wave propagation velocity and density parameters at each grid node, the spatial accuracy, time accuracy and time step parameters of the earthquake wave field numerical calculation, and the forward continuation of the earthquake source wave field is performed using the acoustic wave equation or the elastic wave equation. The finite difference method, the finite element method and the pseudo-spectral method can be used for numerical calculation. The earthquake source wave field is denoted as U S (x, t, s), where x represents the spatial position, t represents the time, and s represents the gun number.

[0010] Furthermore, the method for identifying the first arrival wave of the forward continuation of the source wave field includes the following calculation formula:

[0011]

[0012] Let F i is the first arrival wave identification factor, and its expression is:

[0013]

[0014] Where i and j are the time steps; nl is the selected time window, which is 2-3 main periods of seismic wavelets; w is the seismic wave field value, EW1 i EW2 is the sum of the energy values of all seismic wave fields from the i-nl+1th time step to the i-th time step; i is the sum of the energy values of all seismic wave fields from the 1st time step to the i-th time step; MCM i For EW1 i and EW2 i is the ratio of ; β is the stability factor, which is 0.2; th is the given threshold. When the F value corresponding to the wave field value at a certain moment is 1, it means that the wave field is the first arrival wave.

[0015] Furthermore, the method for identifying the first arrival wave adopts a method such as a long-short time average ratio method, a modified energy ratio method, a modified Coppen method or an Akaike information criterion method.

[0016] Furthermore, the maximum amplitude of the first arrival wave is found at each spatial node. When each spatial node is searched When the source wave field calculation is terminated, the acoustic wave equation or the elastic wave equation is used to calculate the reverse time extension of the detection point wave field. Numerical calculation methods such as the finite difference method, the finite element method and the pseudospectral method can be used.

[0017] Furthermore, the method of cross-correlating the maximum amplitude of the source wavefield with the wavefield of the receiver point at the same time according to the excitation time at each spatial grid node includes:

[0018] Using imaging conditions

[0019]

[0020] Where I(x,s) represents the reverse time migration result of the s-th shot, and the wave field of the detection point is recorded as U R (x, t, s), x represents the spatial position, t represents the time, s represents the gun number, t max represents the maximum propagation time of the maximum amplitude of the first arrival wave;

[0021] The migration results of all shot points are superimposed to obtain the reverse time migration results of the entire section:

[0022]

[0023] Where I(x) represents the superposition of the reverse time migration results of all shots, s max Indicates the maximum number of guns.

[0024] In another aspect, a first-break reverse time migration system based on excitation amplitude comprises

[0025] The calculation module is used to obtain the calculation parameters of the seismic wave field, perform forward extension of the source wave field, and perform reverse time extension of the detection point wave field;

[0026] Identification module, identifies the first arrival wave of the source wave field and retains the maximum amplitude and corresponding time of the first arrival wave at each spatial grid node;

[0027] The imaging module is used to cross-correlate the maximum amplitude of the source wave field with the wave field of the detection point at the same time according to the excitation time at each spatial grid node, and superimpose the cross-correlation values at all times to obtain the reverse time migration result.

[0028] Compared with the prior art, the present invention adopts the above technical solution, and its biggest feature is:

[0029] The present invention provides a method for reverse time migration using the first-arrival wave excitation amplitude. During the propagation of the source wave field, the maximum amplitude and time of the first-arrival wave are saved, and then cross-correlation imaging is performed with the detection point wave field. While maintaining the reverse time migration imaging capability, it can also reduce the amount of calculation, reduce the interference of multipath waves, and improve the signal-to-noise ratio. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 A schematic diagram of the first-arrival wave reverse time migration process of a first-arrival wave reverse time migration method based on excitation amplitude and a system thereof according to the present invention;

[0031] Figure 2 Schematic diagram of a three-layer horizontal layered medium model of a first-arrival wave reverse time migration method based on excitation amplitude and its system according to the present invention;

[0032] Figure 3 A comparison diagram of seismic wave field simulations of a first-arrival reverse time migration method based on excitation amplitude and its system described in the present invention;

[0033] in Figure 3 (a) A snapshot of the seismic wave field simulated by the three-layer model at time 0.75 s. Figure 3 (b) The first arrival wave field identified by this method, Figure 3 (c) Maximum amplitude of the first arrival wave field;

[0034] Figure 4 This is a comparison diagram of reverse time migration results of the first-arrival wave reverse time migration method based on excitation amplitude and the system thereof according to the present invention;

[0035] in Figure 4 (a) The reverse time migration results of the first arrival wave of the present invention, Figure 4 (b) Reverse time migration results of multipath waves;

[0036] Figure 5 Schematic diagram of the Sigsbee2a model of the first-arrival reverse time migration method based on excitation amplitude and its system according to the present invention;

[0037] Figure 6 This is the seismic record of the 100th shot using the first-arrival wave reverse time migration method based on excitation amplitude of the present invention;

[0038] Figure 7 Schematic diagram of the Sigsbee2a migration model of the first-arrival reverse time migration method based on excitation amplitude of the present invention;

[0039] Figure 8 A schematic cross-sectional diagram of excitation amplitude reverse time migration according to an embodiment of the excitation amplitude-based first-break reverse time migration method of the present invention;

[0040] Figure 9 Schematic cross-sectional view of reverse time migration of first-arrival wave excitation amplitude according to an embodiment of the first-arrival wave reverse time migration method based on excitation amplitude of the present invention;

[0041] Figure 10 In the embodiment of the first arrival wave reverse time migration method based on excitation amplitude of the present invention, Figure 8 A detailed enlarged view of the white box 1 in the reverse time migration profile of the excitation amplitude;

[0042] Figure 11 In the embodiment of the first arrival wave reverse time migration method based on excitation amplitude of the present invention, Figure 9 A magnified image of the area indicated by the white box 1 in the reverse time migration profile of the first arrival wave excitation amplitude;

[0043] Figure 12 In the embodiment of the first arrival wave reverse time migration method based on excitation amplitude of the present invention, Figure 8 and Figure 9 A magnified image of the reflection coefficient of the Sigsbee2a salt dome model in the area shown in Box 1;

[0044] Figure 13 The excitation amplitude reverse time migration profile in the embodiment of the first arrival wave reverse time migration method based on the excitation amplitude of the present invention is Figure 8 Enlarged image of the area indicated by white box 2;

[0045] Figure 14 In the embodiment of the first arrival wave reverse time migration method based on excitation amplitude of the present invention, the first arrival wave excitation amplitude reverse time migration profile is Figure 9 Enlarged image of the area indicated by white box 2;

[0046] Figure 15 In the embodiment of the first arrival wave reverse time migration method based on excitation amplitude of the present invention, Figure 8 and Figure 9 A magnified image of the reflection coefficient of the Sigsbee2a salt dome model in the area shown in Box 2; DETAILED DESCRIPTION

[0047] The technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0048] Flowchart as Figure 1 As shown, the present invention includes the following steps:

[0049] (1) Forward extension of the source wave field. In this process, the first arrival wave of the forward extension of the source wave field is identified, and the maximum amplitude and corresponding time of the first arrival wave are retained at each spatial grid node.

[0050] (2) Reverse time extension of the wave field at the detection point,

[0051] (3) At each spatial grid node, the maximum amplitude of the source wave field is cross-correlated with the wave field of the detection point at the same time according to the excitation time, and the cross-correlation values at all times are superimposed to obtain the reverse time migration result.

[0052] In this embodiment, the specific steps include:

[0053] Step 1: Given seismic wave field calculation parameters include grid size, number of grid nodes, physical parameters such as seismic wave propagation velocity and density at each grid node, spatial accuracy, time accuracy, and time step size of seismic wave field numerical calculation.

[0054] Step 2: Based on the parameters of step 1, the acoustic wave equation or elastic wave equation is used to forward extend the source wave field. Numerical calculation methods such as finite difference method, finite element method and pseudospectral method can be used for calculation. The source wave field is denoted as U S (x, t, s), where x represents the spatial position, t represents the time, and s represents the gun number.

[0055] Step 3, identify the first arrival wave of the source wave field. There are many methods for identifying the first arrival wave, such as the long-short time average ratio method (STA / LTA), the modified energy ratio method (MER), the modified Coppen method (MCM) and the Akaike information criterion method (AIC). These methods are all for identifying the first arrival wave of seismic records, and require an overall analysis of the wave field values at all recording times. To apply these methods in reverse time migration, a simple and direct way is to calculate and store the source wave field at all times, and then identify the first arrival wave, which undoubtedly requires a lot of computing time and storage space. In order to save computing power and improve computing efficiency, the present invention hopes to identify the first arrival wave during the calculation of the source wave field. By comparing with existing methods, the present invention selects the MCM method and improves it to meet the needs of computing efficiency. The traditional MCM method is:

[0056]

[0057] MCM i Molecule EW1 i When using equation 1 for calculation, nl additions need to be calculated at each spatial node at each time step. In order to reduce the amount of calculation, the present invention modifies equations (1) and (3) into the following forms:

[0058]

[0059] Where i and j are the time steps; nl is the selected time window, which is generally 2-3 main periods of the seismic wavelet; w is the seismic wave field value, EW1 i EW2 is the sum of the energy values of all seismic wave fields from the i-nl+1th time step to the i-th time step; i is the sum of the energy values of all seismic wave fields from the 1st time step to the i-th time step; MCM i For EW1 i and EW2 i is the ratio of ; β is a stabilizing factor to prevent the denominator from having a value of 0, and is generally taken as 0.2. i The numerator only needs to calculate one addition and one subtraction at each spatial node at each time step, which greatly reduces the amount of calculation.

[0060] Let F i is the first arrival wave identification factor, and its expression is:

[0061]

[0062] Where th is a given threshold. When the F value corresponding to the wave field value at a certain moment is 1, it means that the wave field is a first arrival wave.

[0063] Step 4: Find the maximum amplitude of the first arrival wave at each spatial node And save it in memory, when each spatial node is searched When , the calculation of the source wave field is terminated.

[0064] Step 5: Based on the parameters of step 1, the wave field of the detection point is reversely extended using the acoustic wave equation or the elastic wave equation. Numerical calculation methods such as the finite difference method, the finite element method, and the pseudospectral method can be used for calculation. The wave field of the detection point is denoted as U R (x,t,s);

[0065] Step 6: Use the imaging condition formula:

[0066]

[0067] Where x represents the spatial position, t represents the time, and s represents the shot number. I(x,s) represents the reverse time migration result of the sth shot.

[0068] Step 7: Superimpose the migration results of all shot points to obtain the reverse time migration results of the entire section:

[0069]

[0070] Where I(x) represents the superposition of the reverse time migration results of all shots.

[0071] In this embodiment, a three-layer geological model is used. Figure 2 The model is a three-layer horizontal stratified geological model with two stratum interfaces at depths of 1 km and 2 km respectively. The seismic wave field is simulated using the acoustic wave equation, as Figure 3 (a) The numerical calculation method adopts the finite difference method, the earthquake source is an explosion source, which is placed at a horizontal distance of 1 km from the surface, and the seismic wavelet is a Ricker wavelet with a main frequency of 20 Hz, a spatial step of 10 m, and a time step of 1 ms.

[0072] Figure 3 (a) is a snapshot of the seismic wave field at 0.75s. Figure 3 (b) is the first arrival wave field identified by the method of the present invention, and the MCM is calculated. i When nl is selected as the main period of 2 seismic wavelets, i.e. 100ms, the first arrival wave identification factor F is calculated. i The threshold is set to 0.9. Figure 3 (c) is the maximum amplitude of the first-arrival wave field, also known as the excitation amplitude.

[0073] Figure 4 (a) shows the reverse time migration image of the first arrival wave calculated using the method of the present invention. The total number of shots is 60, with shot locations ranging from 1 km to 4 km, and a shot spacing of 50 m. The image shows that the geological interfaces at 1 km and 2 km underground are accurately imaged. Figure 4 (b) shows the reverse time migration image of the multipath wave. The multipath wave reverse time migration method identifies the first arrival wave during the propagation of the source wavefield, and then considers the subsequent seismic wavefield as a multipath wave. The five maximum amplitudes of the multipath wave and their corresponding times are identified and cross-correlated with the wavefield at the detection point. It can be seen that although a relatively clear horizontal interface is present, its depth is inaccurate. This indicates that the migration result obtained by directly cross-correlating the multipath wave with the wavefield at the detection point is migration noise and cannot obtain an accurate image of the subsurface medium.

[0074] Effects of complex model applications

[0075] Figure 5 This is a seismic wave velocity distribution diagram of the Sigsbee2a salt dome model. When forward modeling earthquake records, the earthquake source is an explosive source and is placed on the surface, with a horizontal position of 0.7km to 21km, a shot distance of 98m, and a total of 215 shots. The seismic wavelet is the Ricker wavelet with a main frequency of 20Hz. The numerical calculation method uses the finite difference method, with a spatial step of 7m, a time step of 0.8ms, and a recording time of 7.4s. The simulated earthquake record is as follows Figure 6 shown.

[0076] The parameters of reverse time migration are the same as those of forward modeling. The migration velocity used is as follows: Figure 7 shown. Figure 8 and Figure 9 They are respectively the reverse time migration profile of the excitation amplitude and the reverse time migration profile obtained by the method of the present invention. Figure 8 The noise is stronger than Figure 9 .

[0077] Figure 10 and Figure 11 They are Figure 8 and Figure 9 , magnified image of the area indicated by white box 1. Figure 12 This is a magnified image of the reflectance coefficient of the Sigsbee2a salt dome model in the corresponding region, clearly indicating the location of the geological interface. The geological interface, indicated by the white arrow, is successfully imaged using the present method, whereas the excitation amplitude method failed to effectively image it. Furthermore, testing has shown that the signal-to-noise ratio of the present method is significantly higher than that of the excitation amplitude imaging method.

[0078] Figure 13 and Figure 14 They are Figure 8 and Figure 9 , an enlarged view of the area indicated by the white box 2. Figure 15 A magnified image of the reflectance coefficient of the Sigsbee2a salt dome model in the corresponding region. Similarly, the formation indicated by the white arrow was successfully imaged by the present invention's method, while the excitation amplitude method failed to effectively image it. Furthermore, the signal-to-noise ratio of the present invention's method was significantly higher than that of the excitation amplitude imaging method.

[0079] The reverse time migration profile obtained by the method of the present invention has a higher signal-to-noise ratio and a stronger imaging capability. The geological interface that the excitation amplitude method failed to image in the example was successfully imaged by the present invention. At the same time, this method only requires the amplitude of the first arrival wave in the source wave field, so there is no need to calculate the source wave field at all times, while the excitation amplitude method needs to search the source wave field at all times to obtain the maximum amplitude. As shown in the example, the earthquake record of Sigsbee2a is 7.4s. The excitation amplitude method needs to calculate the source wave field of 7.4s to obtain the maximum amplitude on all grid points in the calculation area of each shot. The method of the present invention only needs about 4s to terminate the calculation of the source wave field, so the calculation amount of the present invention is smaller and the efficiency is higher.

[0080] The above content is merely an example and explanation of the structure of the present invention. Those skilled in the art may make various modifications or additions to the described specific embodiments or replace them in a similar manner. As long as they do not deviate from the structure of the invention or exceed the scope defined by the claims, they should all fall within the scope of protection of the present invention.

Claims

1. A first-arrival reverse time migration method based on excitation amplitude, characterized in that: The specific steps include: Obtain the calculation parameters of the seismic wave field, perform forward extension on the source wave field, identify the first arrival wave of the source wave field, and retain the maximum amplitude and corresponding time of the first arrival wave at each spatial grid node; Reverse time extension of the wave field at the detection point; Cross-correlating the maximum amplitude of the source wavefield with the wavefield of the receiver point at the same time according to the excitation time at each spatial grid node, and superimposing the cross-correlation values at all times to obtain the reverse time migration result; The method for identifying the first arrival of the earthquake source wave field includes the following calculation formula: （2） （4） （5） make F i is the first arrival wave identification factor, and its expression is: （6） in i and j is the number of time steps; nl is the selected time window, which is 2-3 main periods of the earthquake wavelet; w is the seismic wave field value, It represents the square of the seismic wave field value at the jth time step. EW 1 i From the first time step to the i-nl The sum of the energy values of all seismic wave fields within a time step; EW 2 i From the first time step to the i The sum of the energy values of all seismic wave fields within a time step; MCM i For the i-nl+ 1 time step to the i The sum of the energy values of all seismic wave fields within a time step is EW 2 i The ratio of β is the stability factor, with a value of 0.

2. th It is a given threshold. When the wave field value at a certain moment corresponds to F When the value is 1, it means that the wave field is a first arrival wave.

2. The method for first-break reverse time migration based on excitation amplitude according to claim 1, characterized in that: Obtain the calculation parameters of the source wave field, including the grid size, the number of grid nodes, the seismic wave propagation velocity and density parameters on each grid node, the spatial accuracy, time accuracy and time step parameters of the seismic wave field numerical calculation, use the acoustic wave equation or the elastic wave equation to forward extend the source wave field, and use the finite difference method, finite element method and pseudo-spectral method numerical calculation method to calculate. The source wave field is recorded as , where x represents the spatial position, t Indicates time, s Indicates the gun number.

3. The method for first-break reverse time migration based on excitation amplitude according to claim 1, characterized in that: The method for identifying the first arrival wave adopts the long-short time average ratio method, the modified energy ratio method, the modified Coppen method or the Akaike information criterion method.

4. The method for first-break reverse time migration based on excitation amplitude according to claim 1, characterized in that: Find the maximum amplitude of the first arrival wave at each spatial node , when each spatial node is searched When the source wave field calculation is terminated, the acoustic wave equation or the elastic wave equation is used to perform reverse time extension on the detection point wave field, and the finite difference method, finite element method and pseudo-spectral numerical calculation method are used for calculation.

5. The method for first-break reverse time migration based on excitation amplitude according to claim 1, characterized in that: Methods for cross-correlating the maximum amplitude of the source wavefield with the wavefield of the receiver point at the same time according to the excitation time at each spatial grid node include: Using imaging conditions （7） in Indicates the s The reverse time migration result of the shot, the wave field of the detection point is recorded as , x represents the spatial position, t Indicates time, s Indicates the gun number. t max represents the maximum propagation time of the maximum amplitude of the first arrival wave; The migration results of all shot points are superimposed to obtain the reverse time migration results of the entire section: （8） in It represents the superposition of the reverse time migration results of all guns. s max Indicates the maximum number of guns.

6. A first-break reverse time migration system based on excitation amplitude, used to perform the first-break reverse time migration method based on excitation amplitude according to any one of claims 1 to 5, characterized in that: include The calculation module is used to obtain the calculation parameters of the seismic wave field, perform forward extension of the source wave field, and perform reverse time extension of the detection point wave field; Identification module, identifies the first arrival wave of the source wave field and retains the maximum amplitude and corresponding time of the first arrival wave at each spatial grid node; The imaging module is used to cross-correlate the maximum amplitude of the source wave field with the wave field of the detection point at the same time according to the excitation time at each spatial grid node, and superimpose the cross-correlation values at all times to obtain the reverse time migration result.

Images