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

Abstract

The invention discloses a first-motion wave reverse time migration method and system based on excitation amplitude. The method comprises the steps of obtaining seismic wave field calculation parameters, performing forward continuation on a seismic source wave field, identifying a first-motion wave of the seismic source wave field, and reserving the maximum amplitude and the corresponding time of the first-motion wave on each space grid node. Carrying out reverse time continuation on a wave detection point wave field; and on each space grid node, according to the excitation time, carrying out cross-correlation on the maximum amplitude of the seismic source wave field and the detection point wave field at the same moment, and superposing cross-correlation values at all moments to obtain a reverse time migration result. According to the method, the maximum amplitude and the time of the first-motion wave are stored in the propagation process of the seismic source wave field, and then cross-correlation imaging is performed on the seismic source wave field and the detection point wave field, so that the calculated amount and the interference of multipath waves can be reduced while the reverse time migration imaging capability is ensured, and the signal-to-noise ratio is improved.

Description

Technical Field

[0001] The 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 wave field, (2) reverse time extension of the receiver wave field, and (3) using imaging conditions (such as cross-correlation imaging of the source wave field and the receiver wave field). This process requires the source wave field 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 the corresponding time on each computational grid node during the forward extension of the source wave field, 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 to be a disadvantage of the excitation amplitude imaging condition. A solution strategy is to save several maximum amplitude values ​​and their corresponding times, so as to save a certain amount of multipath waves for imaging. However, imaging of multipath waves is a complex problem. It is not enough to simply save the multipath waves and then use cross-correlation imaging. Targeted processing methods are needed to achieve the purpose of accurate imaging. Otherwise, most of the noise generated in the imaging is noise. In fact, the excitation amplitude imaging condition does not exclude multipath waves. In complex media, sometimes the amplitude of the multipath wave is stronger than the amplitude of the first arrival wave. In this way, when using the excitation amplitude imaging condition, the amplitude of the multipath and its arrival time are saved and imaged.

[0003] In reverse time migration, whether it is the cross-correlation imaging condition or the excitation amplitude imaging condition, the first arrival wave is mainly used to image the underground medium. For complex models, the excitation amplitude of the source wave field saved by the existing excitation amplitude reverse time migration method includes the amplitude of the first arrival wave field and the multipath wave field. For the multipath wave field, direct cross-correlation with the detection point wave field not only cannot image the underground medium, but will generate migration noise, affecting the imaging quality. Therefore, a first arrival wave reverse time migration method based on excitation amplitude and its system 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 calculation parameters of the source wave field, forward extend 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 on 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 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.

[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 on each grid node, the spatial accuracy, time accuracy and time step parameters of the numerical calculation of the earthquake wave field, and the forward extension of the earthquake source wave field using the acoustic wave equation or the elastic wave equation can be calculated using the finite difference method, the finite element method and the pseudo-spectral method numerical calculation method. 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 extension 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 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 ith 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, for example, the long-short time average ratio method, the modified energy ratio method, the modified Coppen method or the 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 calculation of the source wave field 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 using numerical calculation methods such as the finite difference method, the finite element method, and the pseudo-spectral method.

[0017] Furthermore, the method of cross-correlating the maximum amplitude of the source wave field with the wave field 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 sth 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, to carry out forward extension of the source wave field, and to carry out reverse time extension of the detection point wave field;

[0026] An 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 the 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. This method can reduce the amount of calculation, reduce the interference of multipath waves, and improve the signal-to-noise ratio while ensuring the reverse time migration imaging capability. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 A schematic diagram of a 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 A schematic diagram of a three-layer horizontal layered medium model of a first-arrival wave reverse time migration method based on excitation amplitude and a system thereof according to the present invention;

[0032] Figure 3 A comparison diagram of seismic wave field simulation of a first-arrival wave reverse time migration method based on excitation amplitude and a system thereof according to 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 first arrival wave reverse time migration result of the present invention, Figure 4 (b) Reverse time migration results of multipath waves;

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

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

[0038] Figure 7 A schematic diagram of a Sigsbee2a migration model of the first arrival wave 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 of an embodiment of the first-arrival wave reverse time migration method based on excitation amplitude of the present invention;

[0040] Fig. 9 It is a cross-sectional schematic diagram 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] Fig.10 In an 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] Fig.11 In an embodiment of the first arrival wave reverse time migration method based on excitation amplitude of the present invention Fig. 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] Fig.12 In an embodiment of the first arrival wave reverse time migration method based on excitation amplitude of the present invention Figure 8 and Fig. 9 A magnified image of the reflection coefficient of the Sigsbee2a salt dome model in the area shown in Box 1;

[0044] Fig.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] Fig.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 Fig. 9 Enlarged image of the area indicated by white box 2;

[0046] Fig.15 In an embodiment of the first arrival wave reverse time migration method based on excitation amplitude of the present invention Figure 8 and Fig. 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 comprises 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 implementation example, the specific detailed steps include:

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

[0054] Step 2: Based on the parameters in 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 pseudo-spectral 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. Applying 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 calculation amount and improve calculation efficiency, the present invention hopes to identify the first arrival wave during the calculation of the source wave field. Through comparison with existing methods, the present invention selects the MCM method and improves it to meet the needs of calculation efficiency. The traditional MCM method is:

[0056]

[0057] MCM i EW1 i When using equation 1 for calculation, it is necessary to calculate nl additions 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 time steps; nl is the selected time window, which is generally 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 ith 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] Wherein, th is a given threshold. When the F value corresponding to the wave field value at a certain moment is 1, it indicates 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 source wave field is terminated.

[0064] Step 5: According to 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. The wave field of the detection point can be calculated using numerical calculation methods such as the finite difference method, the finite element method, and the pseudo-spectral method. 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 result of the entire section:

[0069]

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

[0071] In this implementation example, a three-layer geological model is used. Figure 2 It is a three-layer horizontal layered geological model, and the depths of the two stratum interfaces are 1km and 2km underground respectively. The acoustic wave equation is used to simulate the seismic wave field, such as Figure 3 As shown in (a), the numerical calculation method adopts the finite difference method, the earthquake source adopts the explosion source, which is placed at a horizontal distance of 1 km from the surface, the seismic wavelet is the Ricker wavelet, the main frequency is 20 Hz, the spatial step is 10 m, and the time step is 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, that is, the excitation amplitude.

[0073] Figure 4 (a) is the first arrival wave reverse time migration image calculated by the method of the present invention. The total number of shots is 60, the shot positions are from 1km to 4km, and the shot distance is 50m. From the figure, it can be seen that the geological interface at 1km and 2km underground is accurately imaged. Figure 4 (b) is the reverse time migration image of the multipath wave. The method of reverse time migration of multipath waves is to identify the first arrival wave during the propagation of the source wave field, and the subsequent seismic wave field is considered to be a multipath wave. The five maximum amplitudes of the multipath wave and their corresponding times are identified and cross-correlated with the detection point wave field. It can be seen that although there is a relatively obvious horizontal interface, its depth is inaccurate. This shows that the migration result obtained by directly cross-correlating the multipath wave with the detection point wave field is a migration noise, and an accurate underground medium image cannot be obtained.

[0074] Complex model application effect

[0075] Figure 5 This is a seismic wave velocity distribution diagram of the Sigsbee2a salt dome model. When simulating earthquake records in the forward model, the earthquake source is an explosion 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 Figure 7 shown. Figure 8 and Fig. 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 Fig. 9 .

[0077] Fig.10 and Fig.11 They are Figure 8 and Fig. 9 , an enlarged view of the area indicated by the white box 1. Fig.12 The enlarged image of the reflection coefficient of the Sigsbee2a salt dome model in the corresponding area can clearly indicate the location of the geological interface. It can be seen that the geological interface indicated by the white arrow is successfully imaged by the method of the present invention, while the excitation amplitude method is not effectively imaged. At the same time, after testing, the signal-to-noise ratio of the method of the present invention is significantly higher than that of the excitation amplitude imaging method.

[0078] Fig.13 and Fig.14 They are Figure 8 and Fig. 9 , an enlarged view of the area indicated by the white box 2. Fig.15 The figure is an enlarged view of the reflection coefficient of the Sigsbee2a salt dome model in the corresponding area. Similarly, the stratum indicated by the white arrow is successfully imaged by the method of the present invention, while the excitation amplitude method does not effectively image it. At the same time, the signal-to-noise ratio of the method of the present invention is 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, the present method only needs 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 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 takes 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 contents are merely examples and explanations of the structure of the present invention. The technicians in this technical field may make various modifications or additions to the specific embodiments described 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 protection scope of the present invention.

Claims

1. A method for first-arrival reverse time migration 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 on each spatial grid node; Reverse time extension of the wave field at the detection point; 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.

2. The method for first-arrival reverse time migration based on excitation amplitude according to claim 1, characterized in that: Obtain the calculation parameters of the earthquake 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 numerical calculation of the earthquake wave field, and use the acoustic wave equation or the elastic wave equation to forward extend the earthquake source wave field. The finite difference method, finite element method and pseudo-spectral numerical calculation method can be used for calculation. The earthquake source wave field is recorded as U S (x, t, s), where x represents the spatial position, t represents the time, and s represents the gun number.

3. The method for first-arrival reverse time migration based on excitation amplitude according to claim 1, characterized in that: The method for identifying the first arrival wave of the source wave field includes the following calculation formula: Let F i is the first arrival wave identification factor, and its expression is: Where i and j are 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 ith time step; MCM i For EW1 i and EW2 i is the ratio of ; β is the stability factor, which takes a value of 0.

2. 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.

4. The method for first-arrival 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.

5. The method of first-arrival 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 calculation of the source wave field 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.

6. The method of first-arrival reverse time migration based on excitation amplitude according to claim 1, characterized in that: The method of cross-correlating the maximum amplitude of the source wave field with the wave field of the receiver point at the same time according to the excitation time at each spatial grid node includes: Using imaging conditions Where I(x,s) represents the reverse time migration result of the sth 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; The migration results of all shot points are superimposed to obtain the reverse time migration results of the entire section: Where I(x) represents the superposition of the reverse time migration results of all shots, s max Indicates the maximum number of guns.

7. A first-break reverse time migration system based on excitation amplitude, used to execute the first-break reverse time migration method based on excitation amplitude according to any one of claims 1 to 6, characterized in that: include The calculation module is used to obtain the calculation parameters of the seismic wave field, to carry out forward extension of the source wave field, and to carry out reverse time extension of the detection point wave field; An 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