A seismic data source wavelet estimation method based on full waveform inversion

Abstract

The present application relates to a kind of seismic data source wave estimation methods based on full waveform inversion, this method is under the framework of direct wave full waveform inversion, by the interlaced of source wave inversion and velocity inversion, high-precision source wave estimation is realized.In the framework of full waveform inversion, the kinematics and dynamics characteristics of wave field are fully considered, the accuracy and resolution of inversion result are guaranteed;Using direct wave information, the nonlinearity of inversion is reduced, and the calculation efficiency of inversion algorithm is guaranteed;Interlaced inversion of source wave and velocity as the core link, effectively avoid the influence of velocity error on source wave estimation.The method of the present application can estimate source wave efficiently and accurately using direct wave full waveform inversion in the case of shallow velocity error.

Description

Technical Field

[0001] This invention belongs to the field of seismic exploration technology, specifically relating to an alternating iterative inversion method for source wavelet and subsurface medium velocity, which is a method for high-precision source wavelet estimation when the subsurface velocity model is inaccurate. Background Technology

[0002] Source wavelet estimation is a crucial step in seismic data inversion, imaging, and interpretation. Conventional source wavelet estimation methods are based on convolution models, which treat the seismic record as a convolution of the source wavelet and the sequence of subsurface reflection coefficients. Source wavelet estimates can be obtained through linear inversion or statistical methods. However, source wavelet estimation based on convolution models is only applicable to post-stack seismic data and neglects complex reflection information such as multiples. Therefore, the obtained source wavelet is not entirely suitable for wave equation-based inversion and imaging.

[0003] The source wavelet is the initial condition for forward modeling of the wave equation, and its accuracy directly affects the final results of inversion and imaging. Currently, in wave equation-based processing such as full-waveform inversion and reverse-time migration imaging, the source wavelet is usually estimated based on direct waves or reflection Green's functions. The main advantage of using direct waves to estimate the source wavelet is its low computational cost and good fit to near-channel data. However, when shallow velocities are inaccurate, the source wavelet estimation method based on direct wave full-waveform inversion will lead to significant errors. Therefore, there is an urgent need to develop a method that can still obtain high-precision source wavelet estimation even with velocity errors. Summary of the Invention

[0004] The purpose of this invention is to provide a method for high-precision source wavelet estimation using short-path direct wave data, thereby addressing the problem of achieving high-precision source wavelet estimation even with velocity errors. Using the fitting of the full waveform information of the short-path direct wave as the objective function, an iterative full waveform inversion update of the source wavelet and subsurface medium velocity parameters is performed to achieve high-precision source wavelet estimation even with errors in the velocity model.

[0005] The objective of this invention is achieved through the following technical solution:

[0006] First, the seismic data undergoes necessary preprocessing to determine the initial source wavelet and initial velocity model used for inversion. Then, forward modeling is performed using the current source wavelet and velocity model to obtain the simulated direct wave data for the current iteration and calculate the objective function for the current iteration. The associated source and associated wavefield of the source wavelet inversion are calculated, and the gradient of the source wavelet inversion is obtained. An appropriate step size is selected to update the source wavelet. The source wavelet update stopping condition is determined; if satisfied, the velocity is updated; otherwise, the wavelet update continues. Using the current source wavelet and current velocity model, simulated direct wave data is calculated, and the objective function value is calculated. The associated source and associated wavefield of the velocity inversion are calculated, and the gradient of the velocity inversion is calculated. An appropriate update step size is selected to update the velocity model. The velocity inversion stopping condition is determined; if satisfied, the next step is performed; otherwise, the velocity model is updated. Finally, an overall iteration stopping condition is determined; if satisfied, the source wavelet estimation result is output; otherwise, the source wavelet update continues.

[0007] A method for estimating seismic source wavelet based on full waveform inversion includes the following steps:

[0008] a. Perform preprocessing on the raw seismic data, such as regularization, denoising, and direct wave interception;

[0009] b. Determine the initial velocity model used for inversion, and select the zero vector as the initial source wavelet;

[0010] c. The objective function shown in formula (2) is the common objective function for source wavelet inversion and velocity inversion.

[0011]

[0012] Where σ is the objective function value, ns and nr are the total number of sources and receivers, respectively, and u and d represent the simulated direct wave data and the observed direct wave data, respectively;

[0013] d. Based on the current velocity model, use the current source wavelet to calculate the simulated direct wave data, and calculate the objective function of the current iteration according to formula (2);

[0014] e. Calculate the associated source of the source wavelet inversion according to formula (3), and calculate the gradient of the source wavelet inversion according to formula (4);

[0015] A = ud (3)

[0016]

[0017] Where A represents the adjoint source, s represents the source function, S represents the source matrix, the superscript T represents the transpose, L represents the forward modeling operator of the wave equation, and its inverse L -1 This represents the backpropagation operator of the wave equation;

[0018] f. Select an appropriate update step size to update the source wavelet;

[0019] g. Determine the stopping condition for the source wavelet update. If it is met, proceed to step h; otherwise, proceed to step d.

[0020] h. Using the current source wavelet, calculate the simulated direct wave data on the current velocity model, and calculate the objective function value of velocity inversion according to formula (1);

[0021] i. Calculate the associated source of velocity inversion according to formula (3), and calculate the gradient of velocity inversion according to formula (5);

[0022]

[0023] Where v represents velocity, v0 represents background velocity, U represents the propagating wave field of the source, and t represents time;

[0024] j. Select an appropriate update step size and update the velocity model;

[0025] k. Determine the stopping condition for velocity inversion. If the condition is met, proceed to step l; otherwise, proceed to step h.

[0026] 1. Determine the overall iteration stopping condition. If it is met, output the source wavelet estimation result; otherwise, execute step d.

[0027] Further, in step b, the expression for the initial velocity model can be written as equation (1).

[0028] v0(i)=v min +(i-1)*(v max -v min ) / (n-1), (1)

[0029] Where v0 is the initial velocity model value, v min and v max These represent the minimum and maximum velocities obtained from velocity analysis, where i is the longitudinal grid coordinate of the model, and n is the maximum number of grid points in the longitudinal direction of the model. The zero vector is chosen as the initial source wavelet.

[0030] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0031] This invention achieves iterative updates of the source wavelet and velocity model within the framework of direct-wave full-waveform inversion, solving the problem of high-precision source wavelet estimation under conditions of inaccurate velocity models. Direct-wave full-waveform inversion fully utilizes the kinematic and dynamic information of the seismic wavefield, ensuring high-resolution source wavelet estimation. Continuous velocity updates based on direct-wave full-waveform inversion during the inversion process effectively correct for the impact of inaccurate velocity on source wavelet estimation. Alternating full-waveform inversion of the source wavelet and velocity model enables high-precision source wavelet estimation even when the initial velocity model is inaccurate.

[0032] This invention has the following characteristics:

[0033] 1. This method only uses near-channel direct wave seismic data, requires a short recording time for inversion and a small model depth, and the algorithm as a whole has high computational efficiency.

[0034] 2. This method uses source wavelet inversion and velocity inversion alternately, which can achieve high-precision source wavelet estimation even when the initial velocity model is inaccurate;

[0035] 3. This method performs source wavelet estimation within the framework of full waveform inversion, making full use of the kinematic and dynamic information of the wave field and ensuring the accuracy and resolution of the inversion results;

[0036] 4. This method does not require prior information on the source wavelet. In the inversion, the source wavelet estimation, which is essentially propagated in reverse time, is used as the initial wavelet, which can robustly achieve high-precision source wavelet estimation. Attached Figure Description

[0037] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation on the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0038] Figure 1 The flowchart of the earthquake source wavelet estimation method based on full waveform inversion described in this invention;

[0039] Figure 2a Realistic speed model;

[0040] Figure 2b Initial velocity model.

[0041] Figure 3a True source wavelet;

[0042] Figure 3b Initial source wavelet;

[0043] Figure 4a Direct inversion results of the source wavelet;

[0044] Figure 4b Results of the first staggered iteration source wavelet inversion;

[0045] Figure 4c Results of the second staggered iteration source wavelet inversion;

[0046] Figure 4d Results of the third staggered iteration source wavelet inversion;

[0047] Figure 4e Results of source wavelet inversion in the 4th staggered iteration;

[0048] Figure 5a Direct inversion of source wavelet extraction shot normalization display;

[0049] Figure 5b Interleaved iterative source wavelet inversion shot extraction normalization display. Detailed Implementation

[0050] The present invention will be further described below with reference to embodiments:

[0051] The seismic source wavelet estimation method based on full waveform inversion described in this invention is implemented using the MATLAB platform. The method includes the following steps:

[0052] 1. Install the MATLAB software platform on a Windows 7 or Linux system, requiring MATLAB R2016a or later, and include the Parallel Computing Toolbox.

[0053] 2. Perform data preprocessing. Perform preprocessing on the raw background noise data, including regularization, denoising, and direct wave interception.

[0054] 3. Through velocity analysis, the maximum and minimum values ​​of the model velocity are obtained, and a linear gradient model that gradually increases with depth within the velocity range is generated as the initial velocity model for inversion. The expression of the initial velocity model can be written as equation (1).

[0055] v0(i)=v min +(i-1)*(v max -v min ) / (n-1), (1)

[0056] Where v0 is the initial velocity model value, v min and v maxThese represent the minimum and maximum velocities obtained from velocity analysis, where i is the longitudinal grid coordinate of the model, and n is the maximum number of grid points in the longitudinal direction of the model. The zero vector is chosen as the initial source wavelet.

[0057] 4. The objective function shown in formula (2) is the common objective function for source wavelet inversion and velocity inversion.

[0058]

[0059] Where σ is the objective function value, ns and nr are the total number of sources and receivers, respectively, and u and d represent the simulated direct wave data and the observed direct wave data, respectively.

[0060] 5. Based on the current velocity model, use the current source wavelet to calculate the simulated direct wave data, and calculate the objective function of the current iteration according to formula (2).

[0061] 6. Calculate the associated source of the source wavelet inversion according to formula (3), and calculate the gradient of the source wavelet inversion according to formula (4);

[0062] A = ud (3)

[0063]

[0064] Where A represents the adjoint source, s represents the source function, S represents the source matrix, the superscript T represents the transpose, L represents the forward modeling operator of the wave equation, and its inverse L -1 This represents the backpropagation operator of the wave equation.

[0065] 7. Select an appropriate update step size to update the source wavelet.

[0066] 8. Determine the stopping condition for the source wavelet update. If it is met, proceed to step 9; otherwise, proceed to step 5.

[0067] 9. Using the current source wavelet, calculate the simulated direct wave data on the current velocity model, and calculate the objective function value of velocity inversion according to formula (2).

[0068] 10. Calculate the associated source of velocity inversion according to formula (3); calculate the gradient of velocity inversion according to formula (5), that is, the velocity inversion gradient is the zero-delay cross-correlation between the forward propagation wave field of the source and the backward propagation wave field of the associated source.

[0069]

[0070] Where v represents velocity, v0 represents background velocity, U represents the propagating wave field of the source, and t represents time.

[0071] 11. Select an appropriate update step size to update the velocity model.

[0072] 12. Determine the stopping condition for velocity inversion. If it is met, proceed to step 13; otherwise, proceed to step 9.

[0073] 13. Determine the overall iteration stopping condition. If it is met, output the source wavelet estimation result; otherwise, proceed to step 5.

[0074] Example 1

[0075] The overall process of this invention is as follows: Figure 1 As shown.

[0076] The actual velocity model used is as follows: Figure 2a As shown, this is the Marmousi velocity model, simulating a land-based sampling scenario without a water layer. The initial velocity model used is a smoothed result of the real velocity model; the initial velocity in the shallow part of the model differs significantly from the real velocity, facilitating verification of the effectiveness of the method of this invention. The real source wavelet used is as follows: Figure 3a As shown, this is a 10Hz dominant frequency Ricker wavelet; the initial source wavelet used is as follows: Figure 3b As shown, this is the zero vector. A total of 48 seismic sources are uniformly distributed on the model surface, and each real seismic source has the same source function, as shown below. Figure 3a As shown.

[0077] First, the source wavelet is directly inverted using the direct wave full waveform inversion, and the result is as follows: Figure 4a As shown in the figure, there are significant lateral variations in the 48 source functions obtained through inversion. This lateral inconsistency is due to source wavelet estimation errors caused by velocity errors. Then, using the staggered iterative inversion of source wavelets and velocities proposed in this invention, the results of the 1st, 2nd, 3rd, and 4th staggered iterative inversions are shown in the figure. Figure 4b , Figure 4c , Figure 4d and Figure 4e As shown, the lateral consistency of the source wavelet inversion results improves with staggered iterations, indicating that the influence of velocity errors is gradually being corrected. The inversion results are then normalized by shot extraction, as shown in Figure 5a and... Figure 5b The images show the results of a single-shot inversion of the source wavelet using direct inversion and the results of a single-shot inversion using staggered iterative inversion. The blue and red lines represent the actual source wavelet and the estimated wavelet result, respectively. It is evident that the inversion method proposed in this invention can obtain source wavelet estimates with more accurate amplitude and phase. Therefore, the method of this invention can obtain high-precision source wavelet inversion results even when errors exist in the underground velocity model.

[0078] This invention discloses a source wavelet estimation method for seismic data based on full-waveform inversion. The core step is to perform staggered iterative inversion of the source wavelet and subsurface velocity within a direct-wave full-waveform inversion framework. Utilizing direct-wave full-waveform inversion reduces the nonlinearity of the inversion while ensuring the overall efficiency of the algorithm. Performing source wavelet inversion within the full-waveform inversion framework fully utilizes the kinematic and dynamic information of the wavefield, guaranteeing the accuracy and resolution of the inversion results. The staggered inversion of source wavelet and subsurface velocity corrects for the influence of velocity errors on source wavelet estimation during the inversion process, thus achieving high-precision source wavelet estimation even when the initial velocity model contains errors.

[0079] Note that the above description is merely a preferred embodiment of the present invention and the technical principles employed. Those skilled in the art will understand that the present invention is not limited to the specific embodiments described herein, and various obvious changes, readjustments, and substitutions can be made without departing from the scope of protection of the present invention. Therefore, although the present invention has been described in detail through the above embodiments, the present invention is not limited to the above embodiments, and may include many other equivalent embodiments without departing from the concept of the present invention, the scope of which is determined by the scope of the appended claims.

Claims

1. A method for estimating seismic source wavelet based on full waveform inversion, characterized in that, Includes the following steps: a. Perform regularization, denoising, and direct wave interception preprocessing on the raw seismic data; b. Determine the initial velocity model used for inversion, and select the zero vector as the initial source wavelet; c. The objective function shown in formula (2) is the common objective function for source wavelet inversion and velocity inversion. Where σ is the objective function value, ns and nr are the total number of sources and receivers, respectively, and u and d represent the simulated direct wave data and the observed direct wave data, respectively; d. Based on the current velocity model, use the current source wavelet to calculate the simulated direct wave data, and calculate the objective function of the current iteration according to formula (2); e. Calculate the associated source of the source wavelet inversion according to formula (3), and calculate the gradient of the source wavelet inversion according to formula (4); A = ud (3) Where A represents the adjoint source, s represents the source function, S represents the source matrix, the superscript T represents the transpose, L represents the forward modeling operator of the wave equation, and its inverse L -1 This represents the backpropagation operator of the wave equation; f. Select an appropriate update step size to update the source wavelet; g. Determine the stopping condition for the source wavelet update. If it is met, proceed to step h; otherwise, proceed to step d. h. Using the current source wavelet, calculate the simulated direct wave data on the current velocity model, and calculate the objective function value of velocity inversion according to formula (1); i. Calculate the associated source of velocity inversion according to formula (3), and calculate the gradient of velocity inversion according to formula (5); Where v represents velocity, v0 represents background velocity, U represents the propagating wave field of the source, and t represents time; j. Select an appropriate update step size and update the velocity model; k. Determine the stopping condition for velocity inversion. If the condition is met, proceed to step l; otherwise, proceed to step h.

1. Determine the overall iteration stopping condition. If it is met, output the source wavelet estimation result; otherwise, execute step d.

2. The method for estimating seismic source wavelet based on full waveform inversion according to claim 1, characterized in that: Step b, the expression for the initial velocity model can be written as equation (1). v0(i)=v min +(i-1)*(v max -v min ) / (n-1), (1) Where v0 is the initial velocity model value, v min and v max These are the minimum and maximum velocities obtained from velocity analysis, i is the longitudinal grid coordinate of the model, and n is the maximum number of grid points in the longitudinal direction of the model. The zero vector is selected as the initial source wavelet.

Images