A seismic wavelet inversion method based on traveltime correction

Abstract

The invention relates to a seismic wavelet inversion method based on travel time correction, which belongs to the field of seismic inversion and specifically comprises the following steps: first, forward simulation is performed using an initial seismic wavelet and an initial velocity model to obtain a simulated seismic record, and the direct waves of the simulated seismic record and the actual seismic record obtained in the field are intercepted respectively, and then the intercepted direct wave of the actual seismic record and the direct wave of the simulated seismic record are subtracted to obtain a residual record, and the residual record is used as a disturbance for reverse time extension to obtain a companion wave field, and in the next step, each seismic record in the companion wave field is delayed in time using the initial velocity model information and the detection point position information and averaged, and finally, the average value is used as the wavelet update amount to update the wavelet. Model experiments show that the seismic wavelet obtained by inversion after time correction has a higher degree of agreement with the real seismic wavelet, and a velocity model with higher accuracy can be obtained by using it for velocity model inversion.

Description

Technical Field

[0001] The invention belongs to the field of seismic inversion, and in particular relates to a seismic wavelet inversion method based on travel time correction. Background Art

[0002] The full waveform inversion method is a pre-stack data imaging and velocity modeling algorithm. It estimates underground physical parameters by minimizing the difference between simulated data and observed data through optimization iteration. It can make full use of the travel time, amplitude, phase and other information in the seismic data to reconstruct the velocity, density and other parameter models of the underground medium. Since full waveform inversion has the potential to characterize the velocity model of complex underground structures, this method has attracted widespread attention in the industry and has become a research hotspot in the field of geophysics.

[0003] Full waveform inversion updates the model by matching the waveform. The quality of data waveform matching is extremely dependent on the quality of the given initial seismic wavelet. The closer the initial seismic wavelet is to the real seismic wavelet, the better the inversion effect. However, in conventional waveform inversion, the real seismic wavelet is unknown. Therefore, the theoretical wavelet given based on experience is often very different from the real wavelet, which seriously affects the subsequent inversion accuracy. Some scholars invert seismic wavelets through waveform inversion, but they do not consider the interference during travel. The inverted wavelet and the real wavelet often have a certain gap, which affects the accuracy of the velocity model in the subsequent full-wave inversion and restricts the application of this method in actual data processing. Summary of the invention

[0004] The present invention provides a seismic wavelet inversion method based on travel time correction in view of the above technical problems, belonging to the field of seismic inversion. The method uses the initial seismic wavelet and the initial velocity model for forward simulation to obtain a simulated seismic record, and respectively intercepts the direct wave of the simulated seismic record and the actual seismic record obtained in the field, and then the intercepted direct wave of the actual seismic record and the direct wave of the simulated seismic record are subtracted to obtain a residual record, and the residual record is used as a disturbance for reverse time extension to obtain a companion wave field, and the next step is to use the initial velocity model information and the detection point position information to time delay each seismic record in the companion wave field and average it, and finally, the average value is used as the wavelet update amount to update the wavelet. Model experiments show that the seismic wavelet obtained by inversion after time correction has a higher degree of consistency with the real seismic wavelet, and using it to perform velocity model inversion can obtain a more accurate velocity model.

[0005] The present invention adopts the following technical solutions:

[0006] A seismic wavelet inversion method based on traveltime correction, the method specifically comprising the following steps:

[0007] (1) Given the initial seismic wavelet W(t) and the initial velocity model V(x, z), where t represents time, x and z represent spatial coordinates, x=1, 2, 3, …, Nx, z=1, 2, 3, …, Nz, and Nx and Nz represent the total number of horizontal and vertical grid points in the model, respectively; based on the initial seismic wavelet and the initial velocity model, the second-order scalar acoustic wave equation is used for forward modeling to obtain the simulated seismic record R cal (x, z, t);

[0008] (2) Intercept the simulated earthquake records R cal (x, z, t) and the actual earthquake record R obtained in the field obs The direct wave of (x, z, t) is then subtracted from the intercepted direct wave of the actual earthquake record and the simulated direct wave of the earthquake record to obtain the residual record R cc (x, z, t) and then use the residual record as a disturbance to perform reverse time extension to obtain the adjoint record R bs (x, z, t);

[0009] (3) Use formula (1) to obtain the accompanying record R bs The time delay t required for each seismic record in (x, z, t) del (x, z, t), where x s , z s are the horizontal and vertical coordinates of the earthquake source, respectively, r , z r are the horizontal and vertical coordinates of the detection point respectively, and the formula (1) is: ;

[0010] (4) Using formula (2), the accompanying record R bs (x, z, t) is delayed and averaged to get the wavelet update wg k (t), where k is the number of iterations, R num is the number of detection points, the formula (2) is: ;

[0011] (5) Update the wavelet using formula (3), where α k is the iteration step length, obtained using formula (4); the formula (3) is , where w k+1 and w k+1 are the seismic wavelets obtained at the k+1th and kth iterations, k=1, 2, 3, …, N k , where N k is the total number of iterations, when k=1, w k is the initial seismic wavelet; the formula (4) is , α k is the k-th iteration step, T represents the cal Perform a transpose operation;

[0012] The above five steps need to be repeated in each iterative inversion of the seismic wavelet based on travel time correction.

[0013] The beneficial effects of the present invention compared with the prior art: The seismic wavelet obtained by inversion after time correction by the method of the present invention is more consistent with the real seismic wavelet. Using it for velocity model inversion can improve the accuracy of the inversion model, thereby accelerating the inversion convergence speed, saving time, and promoting the industrialization process of full waveform inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Figure 1 This is a flow chart of the conventional seismic wavelet inversion method;

[0015] Figure 2 It is a flow chart of a seismic wavelet inversion method based on travel time correction;

[0016] Figure 3 For the Marmousi model;

[0017] Figure 4 is the smoothed velocity model;

[0018] Figure 5 is the real earthquake wavelet;

[0019] Figure 6 For actual earthquake records;

[0020] Figure 7 is the initial seismic wavelet;

[0021] Figure 8 To simulate earthquake records;

[0022] Fig. 9 For accompanying records;

[0023] Fig.10 is the wavelet update amount;

[0024] Fig.11 is the seismic wavelet obtained by conventional inversion method;

[0025] Fig.12 is the seismic wavelet obtained based on the travel time correction method;

[0026] Fig.13 This is a comparison diagram of seismic wavelets;

[0027] Fig.14 This is the model obtained by conventional seismic wavelet inversion;

[0028] Fig.15 This is the model obtained by seismic wavelet inversion using the travel time correction method. DETAILED DESCRIPTION

[0029] The technical solution of the present invention is further explained below through embodiments in combination with the accompanying drawings, but the protection scope of the present invention is not limited in any form by the embodiments.

[0030] In the conventional seismic wavelet inversion, the present invention introduces travel time information to correct each seismic record in the accompanying wave field, and inverts the seismic wavelet based on the corrected seismic record (see the specific process). Figure 2 The experiment shows that the seismic wavelet obtained by inversion after time correction is more consistent with the real seismic wavelet. By using it for velocity model inversion, a more accurate velocity model can be obtained.

[0031] Example 1

[0032] The present invention takes the Marmousi model as an example (such as Figure 3 The model is illustrated in Figure 2. The horizontal length of the model is 4610m, the vertical depth is 1730m, and the grid step size in the x and z directions of the model is 10m.

[0033] Based on this model, a total of 461 blasts were conducted, with the blast points evenly distributed at intervals of 10 m directly above the model. 461 detection points were set, with the detection points also evenly distributed at intervals of 10 m directly above the model. The depths of the blast points and detection points were both 0 m.

[0034] In the following embodiments, an actual earthquake record is first obtained by simulation. The actual earthquake record is simulated by using a 10 Hz Ricker wavelet.

[0035] The specific implementation mode of the present invention is described in detail below. Figure 2 .

[0036] (1) The horizontal layered model with a horizontal length of 4610 m and a vertical depth of 1730 m is meshed. The mesh step size in the x and z directions of the model is 10 m. The real velocity model with a mesh size of 461 × 173 is obtained. The initial velocity model V(x, z) is obtained by smoothing it (e.g. Figure 4 As shown in the figure), where x and z represent the grid position coordinates (x=1, 2, 3, …, 461, z=1, 2,3, …, 173); given the real seismic wavelet W real (t), which is a 10 Hz Ricker wavelet (such as Figure 5 As shown in the figure, where t represents time, based on the real seismic wavelet and the real velocity model, the second-order scalar acoustic wave equation is used for forward simulation to obtain the actual seismic record Robs (x, z, t) (e.g. Figure 6 As shown); Given the initial seismic wavelet W(t) which is a 5 Hz Ricker wavelet (as shown Figure 7 As shown in the figure, where t represents time; based on the initial seismic wavelet and initial velocity model, the second-order scalar acoustic wave equation is used for forward simulation to obtain the simulated seismic record R cal (x, z, t) (e.g. Figure 8 shown);

[0037] (2) Intercept the simulated earthquake records R cal (x, z, t) and the actual earthquake record R obs The direct wave of (x, z, t) is then subtracted from the intercepted direct wave of the actual earthquake record and the simulated direct wave of the earthquake record to obtain the residual record R cc (x, z, t) and then use the residual record as a disturbance to perform reverse time extension to obtain the adjoint record R bs (x, z, t) (e.g. Fig. 9 shown);

[0038] (3) Use formula (1) to obtain the accompanying record R bs The time delay t required for each seismic record in (x, z, t) del (x, z, t), where x s , z s are the horizontal and vertical coordinates of the earthquake source, respectively, r , z r are the horizontal and vertical coordinates of the detection point respectively, and the formula (1) is: ;

[0039] (4) Using formula (2), the accompanying record R bs (x, z, t) is delayed and averaged to get the wavelet update wg k (t) (e.g. Fig.10 As shown), where k is the number of iterations, R num is the number of detection points, the formula (2) is: ;

[0040] (5) Update the wavelet using formula (3), where α k is the iteration step length, obtained using formula (4), the formula (3) is , the formula (4) is ;

[0041] The above five steps need to be repeated in each iterative inversion of the seismic wavelet based on travel time correction.

[0042] In order to illustrate the effectiveness of the method described in the present invention, a comparison is made here with a conventional seismic wavelet inversion method. The process of the conventional seismic wavelet inversion method is as follows: Figure 1 shown. Fig.11 is the seismic wavelet obtained by iterating 10 times the conventional inversion method. Fig.12 The seismic wavelets obtained by iterating 10 times based on the inversion method of travel time correction are compared with the real seismic wavelet ( Figure 5 ) are drawn into the same image, such as Fig.13 As shown in the figure, it can be seen that the seismic wavelet obtained by the inversion method of this paper completely coincides with the real wavelet, while the seismic wavelet obtained by the traditional method still has a certain gap with the real wavelet. The seismic wavelets obtained by the conventional method and the seismic wavelet obtained by the method of this paper are used for full wave inversion, and the inversion models obtained by iterative inversion 200 times are shown below. Fig.14 , Fig.15 As shown in the figure, it can be seen that the model structure is more clearly depicted, the event axis is more continuous, and the accuracy is higher when the seismic wavelet obtained in this paper is used for full waveform inversion. Model experiments show that this method can obtain more accurate seismic wavelets, which can be used for velocity model inversion to improve the accuracy of the inversion model, thereby accelerating the inversion convergence speed, saving time, and promoting the industrialization process of full waveform inversion.

Claims

1. A seismic wavelet inversion method based on travel time correction, characterized in that: The method specifically comprises the following steps: (1) Given the initial seismic wavelet W(t) and the initial velocity model V(x,z), where t represents time, x and z represent spatial position coordinates, x = 1, 2, 3, ..., Nx, z = 1, 2, 3, ..., Nz, and Nx and Nz represent the total number of horizontal and vertical grid points in the model, respectively; based on the initial seismic wavelet and the initial velocity model, the second-order scalar acoustic wave equation is used for forward simulation to obtain the simulated seismic record R cal (x,z,t); (2) Intercept the simulated earthquake records R cal (x,z,t) and the actual earthquake record R obtained in the field obs The direct wave of (x, z, t) is then subtracted from the intercepted direct wave of the actual earthquake record and the simulated direct wave of the earthquake record to obtain the residual record R cc (x,z,t), and use the residual record as a disturbance to perform reverse time extension to obtain the adjoint record R bs (x,z,t); (3) Using formula (1) to obtain the accompanying record R bs The time delay t required for each seismic record in (x,z,t) Del (x,z,t), where x s , z s are the horizontal and vertical coordinates of the earthquake source, respectively, r , z r are the horizontal and vertical coordinates of the detection point respectively, and the formula (1) is: (4) Using formula (2), the accompanying record R bs (x,z,t) is averaged after time delay to get the wavelet update wg k (t), where k is the number of iterations, R num is the number of detection points, the formula (2) is: The wavelet is updated using formula (3), where α k is the iteration step length, obtained using formula (4); the formula (3) is w k+1 (t) = w k (t)-α k *wg k (t), where w k+1 and w k are the seismic wavelets obtained at the k+1th and kth iterations, k=1,2,3,…,N k , where N k is the total number of iterations, when k=1, w k is the initial seismic wavelet; the formula (4) is α k is the k-th iteration step, T represents the cal (x,z,t) performs the transpose operation.

Images