Exponential decay based time-frequency domain phase inversion method

Abstract

The application provides an exponential decay-based time-frequency domain phase inversion method for obtaining high-precision velocity modeling of the underground. The method can suppress the phase noise interference caused by the conventional time-frequency domain phase inversion method and improve the precision of velocity modeling. First, the seismic data is transformed into the time-frequency domain, a time-frequency domain phase inversion objective function is constructed, and a decay factor is introduced to achieve the effect of layered multi-scale inversion. Then, the gradient and the accompanying source of the time-frequency domain phase inversion objective function based on the exponential decay are derived. Finally, the local optimization algorithm is used to update the velocity model, and a better initial model is provided for the conventional full waveform inversion. The application uses the Marmousi velocity model for numerical test and verifies that the method can provide a better initial velocity modeling for the conventional FWI and obtain high-precision inversion results.

Description

Technical Field

[0001] This invention is a time-frequency domain phase inversion method based on exponential decay. Specifically, it involves performing short-time Fourier transform on seismic data to extract time-frequency domain phase information and adding an exponential decay term to model seismic wave velocity. This method can alleviate the main problem in full waveform inversion: period jumps caused by the lack of low-frequency information, and provide a high-precision velocity modeling result for subsequent conventional full waveform inversion. Background Technology

[0002] In the 1980s, Lailly and Tarantola reshaped Clearbout's migration imaging principle into a least-squares local optimization problem, pioneering the Full Waveform Inversion (FWI) theory. Conventional FWI primarily uses an L2-norm objective function to fit simulated and observed data, minimizing the residuals and updating the corresponding subsurface parameters to obtain high-precision inversion results. However, due to the computationally intensive nature of directly calculating the Fréchet derivative, Tarantola pointed out that the gradient of the velocity model can be directly obtained using the adjoint state method, i.e., the zero-delay cross-correlation between the forward-modeled wavefield and the residual backpropagated wavefield, effectively reducing the time consumption of time-domain FWI. However, due to the strong nonlinearity of conventional L2-norm-based FWI, its inversion accuracy largely depends on the low-frequency information of the seismic data, the selection of the initial model, and the accuracy of the source estimation. But in actual seismic exploration, it is difficult to obtain low-frequency data below 4Hz, and low-frequency data is easily contaminated by noise. This can cause conventional FWI gradients to be updated in the wrong direction, resulting in severe period jump problems and failing to obtain high-precision inversion results.

[0003] To overcome the strong nonlinearity of the objective function of conventional FWI, many methods have been proposed. Bunks proposed a multi-scale inversion strategy using seismic data from low to high frequencies in 1995. This frequency-divisional inversion strategy effectively alleviates the period jump problem in full-waveform inversion. Similarly, Pratt pioneered the application of FWI in the frequency domain. Frequency-domain full-waveform inversion has the natural advantage of frequency-divisional multi-scale inversion, requiring only discrete frequency values ​​to invert to a high-resolution velocity model. Luo et al. proposed a travel-time inversion method based on the wave equation; using a cross-correlation function to calculate the travel time difference between observed and simulated data avoids manually picking up the travel time of seismic events, effectively improving the accuracy of travel-time inversion compared to travel-time tomography. However, the high dependence of FWI on low-frequency data and the low quality of actual field acquisition data severely limit its application in practical exploration. Low-frequency reconstruction of seismic data is an effective way to alleviate the period jump problem in FWI. Some researchers proposed an envelope method for seismic data, pointing out that envelope signals can effectively reduce the nonlinearity of the FWI objective function and obtain large-scale structural information of the subsurface velocity model. Wu et al., in their work on digital signal theory, analyzed in detail how envelopes can demodulate missing low-frequency signals in seismic data, thereby inverting smooth large-scale subsurface structures, and proposed objective functions based on the L2 norm objective function, specifically envelope square and envelope logarithm. Xiong et al. proposed an improved envelope inversion method, analyzing in detail the problems of traditional envelope inversion, namely, the suppression of reflection information by the accompanying source, leading to the loss of unmatched reflected wave information in the accompanying source. The improved envelope inversion effectively solved this problem. Oh and Alkhalifah et al. combined signal envelopes and global cross-correlation objective functions, reducing amplitude errors and the strong nonlinearity of the conventional L2 norm objective function, while utilizing the macroscopic construction of the cross-correlation envelope inversion model, and successfully applied it to the North Sea OBC data.

[0004] Compared to amplitude information, phase information contains the kinematic characteristics of seismic wave propagation, revealing information such as the propagation path and subsurface medium. Sun and Schuster proposed a time-domain phase inversion objective function, ignoring variations in seismic data amplitude and focusing on utilizing seismic wave phase information for optimized inversion. Phase inversion can gradually correct and improve upon initial model errors through iteration and optimization to obtain relatively reliable inversion results. Kim extended the time-domain phase inversion method to the frequency domain, constructing a frequency-domain pure phase objective function. The phase inversion method based on the frequency-domain inversion framework can perform frequency-division inversion and obtain better initial models more efficiently. Luo et al. used Hilbert transform to extract instantaneous phase information of seismic signals and applied it to time-domain waveform inversion. Hu et al. developed a time-frequency domain weighted phase cross-correlation waveform inversion method, which utilizes short-time Fourier transform to fully consider the local characteristics of seismic data and obtains better inversion results, verifying the significant advantages of phase inversion.

[0005] Due to the inherent characteristics of signal phase information, the pure phase objective function is easily affected by noise, making it difficult to obtain the long-wavelength structure of the subsurface velocity model. Therefore, based on research in time-frequency domain phase inversion, this invention proposes a time-frequency domain phase inversion method based on exponential decay. This method can effectively suppress the noise influence introduced by time-frequency domain phase inversion, improving inversion accuracy. During the inversion process, the value of the attenuation factor is continuously adjusted to achieve a hierarchical, multi-scale inversion effect, suppressing noise interference, improving inversion accuracy, and providing a better initial velocity model for subsequent conventional FWI. Summary of the Invention

[0006] The main technical innovation of this invention is as follows: By extracting time-frequency domain phase information from seismic data through short-time Fourier transform, and adding an attenuation factor to construct a time-frequency domain phase inversion objective function based on exponential decay, high-precision inversion results of the subsurface velocity model are obtained. This method can suppress noise interference caused by the time-frequency domain phase inversion objective function. During the inversion process, the attenuation factor is adjusted, and the inversion proceeds from the shallow to the deep parts of the model, forming a layered multi-scale strategy to recover the low wavenumber components of the velocity model.

[0007] This invention is a time-frequency domain phase inversion method based on exponential decay, the main steps of which are as follows:

[0008] Step 1. Use the Crews toolkit to preprocess the acquired seismic data, and then input the processed seismic data as observation data into the time-frequency domain phase inversion method based on exponential decay.

[0009] Step 2. Based on the geological information, well logging information, and travel time information of the work area, set the initial velocity model for inversion as the initial model for the time-frequency domain phase inversion method based on exponential decay.

[0010] Step 3. Read the observation system and input the source wavelet. Perform forward modeling using the initial velocity model, and store the wavefield values ​​at each sampling time point for subsequent cross-correlation calculations to update the gradient of the velocity model.

[0011] To simplify the entire physical process and accelerate the inversion efficiency, this study is based on the constant density acoustic wave equation.

[0012]

[0013] Where u represents the wave field value, z and x are the longitudinal and transverse coordinates, t is time, and f(t) is the source wavelet.

[0014] Step 4. Time-frequency transformation of seismic data. The seismic data undergoes a short-time Fourier transform; in this study, the Gabor transform is used. The Gabor transform for time-domain data is defined as follows:

[0015]

[0016] u(τ) represents time-domain seismic data, h represents the Gaussian window function, and t and ω represent the time and frequency coordinates in the time spectrum. F represents the time-frequency domain representation of Gabor-fed seismic data. h [·] denotes the Gabor transform operator for the data. The inverse Gabor transform is defined as:

[0017]

[0018] in This represents the Gabor inverse transform operator for the data.

[0019] Step 5. Calculate the time-frequency domain phase information of the seismic signal. The time-frequency domain signal of the seismic data can be represented as:

[0020]

[0021] The exponential phase of time-frequency domain seismic data is then expressed as:

[0022]

[0023] in The time-frequency domain envelope of earthquake data.

[0024] Step 6. Construct a time-frequency domain phase inversion objective function based on exponential decay. Phase information contains the kinematic characteristics of seismic wave propagation, revealing information such as the propagation path of seismic waves and the subsurface medium. Time-frequency domain phase information focuses more on the local phase changes of the seismic signal. Furthermore, adding a decay factor to the pure phase objective function in the time-frequency domain can effectively alleviate the nonlinear characteristics of the objective function, providing a better initial velocity model for conventional full-waveform inversion and mitigating the period jump problem during the inversion process.

[0025] The objective function of a conventional FWI can be expressed as:

[0026]

[0027] Where u and d represent simulated seismic data and actual observation data, respectively. ns and nr represent the number of source shot points and the number of geophones, respectively, and v represents the velocity model parameters.

[0028] The objective function for time-frequency domain phase inversion is:

[0029]

[0030] and These represent the time-frequency phase of the simulated data and the observed data, respectively.

[0031] The objective function for time-frequency domain phase inversion based on exponential decay can be expressed as:

[0032]

[0033] Where a is the introduced attenuation factor

[0034] Step 7. Calculate the partial derivatives of the objective function with respect to the model velocity parameters and the associated seismic sources. First, the partial derivatives of the conventional FWI with respect to the velocity parameters can be expressed as:

[0035]

[0036] The corresponding accompanying epicenter is R. s =(ud), which is the residual record used for backpropagation to the model space. The partial derivative of the time-frequency domain phase inversion objective function with respect to the velocity parameters can be expressed as:

[0037]

[0038] Where Re[·] denotes taking the real part of the seismic data, and * denotes complex conjugate. Further simplification using the chain rule yields:

[0039]

[0040] in This can be further simplified to:

[0041]

[0042] in Let the Gabor inverse transform operator be represented, then the adjoint source of the objective function of conventional time-frequency domain phase inversion can be expressed as:

[0043]

[0044] Based on the above derivation, the partial derivative of the time-frequency domain phase inversion objective function based on exponential decay with respect to the velocity parameter can be expressed as:

[0045]

[0046] Where Re[·] denotes taking the real part of the seismic data, and * denotes complex conjugate. Further simplification using the chain rule yields:

[0047]

[0048] in This can be further simplified to:

[0049]

[0050] in The adjoint source of the Gabor inverse transform operator, based on the exponentially decaying time-frequency domain phase inversion objective function, can be expressed as:

[0051]

[0052] Step 8. Calculate the direction of velocity model update. The gradient of the conventional FWI objective function can be directly obtained by the zero-delay cross-correlation between the forward modeling wavefield and the residual backpropagated wavefield, expressed as:

[0053]

[0054] Where Nt represents the number of samples in the time domain. u represents the second-order partial derivative of the forward-modeled wavefield with respect to time. b This represents the residual backpropagation wave field.

[0055] Step 9. Update and iterate the velocity parameters using an optimization algorithm. This invention selects the steepest descent method, and the update and iteration formula can be expressed as:

[0056]

[0057] Where α k Indicates the update step size, v k+1 This indicates that the velocity model v from the previous step... k Update the obtained speed parameters.

[0058] This represents the update gradient of the velocity model.

[0059] Step 10. Output the inversion results. Determine whether the inversion results of each step meet the termination conditions set for the objective function. If they do, output the inversion results; otherwise, use the current inversion results as the initial velocity model for the next inversion and continue iterating until the termination conditions are met. The inversion process uses the objective function as a constraint, and the iteration process requires the objective function value to decrease continuously.

[0060] This invention proposes a time-frequency domain phase inversion method based on exponential decay for obtaining high-precision subsurface velocity modeling. This method can suppress phase noise interference from conventional time-frequency domain phase inversion methods, improving the accuracy of velocity modeling. First, seismic data is transformed to the time-frequency domain, a time-frequency domain phase inversion objective function is constructed, and an attenuation factor is introduced to achieve a layered, multi-scale inversion effect. Then, the gradient and associated seismic sources of the exponentially decaying time-frequency domain phase inversion objective function are derived. Finally, a local optimization algorithm is used to update the velocity model, providing a better initial model for conventional full-waveform inversion. This invention uses the Marmousi velocity model for numerical testing and verification, demonstrating that this method can provide a good initial velocity model for conventional FWI and obtain high-precision inversion results. Attached Figure Description

[0061] Figure 1 This is a flowchart of the time-frequency domain phase inversion method based on exponential decay;

[0062] Figure 2 It is the Marmousi velocity model, (a) the initial velocity model; (b) the actual velocity model;

[0063] Figure 3These are the results of conventional FWI inversion: (a) inversion results in the low-frequency band (4-6Hz); and (b) inversion results in the high-frequency band (4-10Hz).

[0064] Figure 4 The results are time-frequency domain phase inversion + conventional FWI inversion: (a) low-frequency band time-frequency domain phase inversion results (4-6Hz); (b) low-frequency band conventional FWI inversion results (4-6Hz) with (a) as the initial model; (c) high-frequency band conventional FWI inversion results (4-10Hz) with (b) as the initial model.

[0065] Figure 5 The results are based on exponential decay time-frequency domain phase inversion and conventional FWI inversion. (a) Low-frequency band time-frequency domain phase inversion results based on exponential decay (4-6Hz), with an attenuation factor of 4; (b) Low-frequency band time-frequency domain phase inversion results based on exponential decay (4-6Hz) with (a) as the initial model, with an attenuation factor of 3; (c) Time-frequency domain phase inversion results based on exponential decay (4-6Hz) with (b) as the initial model, with an attenuation factor of 1; (d) Low-frequency band conventional FWI inversion results (4-6Hz) with (c) as the initial model; (e) High-frequency band conventional FWI inversion results (4-10Hz) with (d) as the initial model.

[0066] Specific examples

[0067] Conventional full-waveform inversion obtains inversion results by optimizing the least-squares objective functional and updating subsurface physical parameters. However, the optimization process often gets stuck in local minima, failing to obtain inversion results from a high-precision subsurface velocity model. This is mainly because low-frequency data is difficult to obtain in actual seismic exploration, and the objective function of low-frequency seismic data has fewer local extrema, which further exacerbates the nonlinearity of the inversion, leading to periodic jumps in the matching process between observed and simulated data, resulting in erroneous inversion results. To address these issues and alleviate the periodic jump problem in waveform inversion, this invention proposes a time-frequency domain phase inversion method based on exponential decay. Phase information contains the kinematic characteristics of seismic wave propagation, revealing information such as the propagation path of seismic waves and the subsurface medium. However, time-frequency domain phase inversion introduces a lot of interference. By adding a decay factor and continuously adjusting its weight in the inversion process, a layered multi-scale inversion effect can be achieved, effectively mitigating the periodic jump problem in waveform inversion.

[0068] This invention utilizes finite-difference forward modeling for model testing. The experimental seismic source consisted of 40 shots, uniformly distributed across the surface of the velocity model. The source used a Ricker wavelet with a dominant frequency of 8 Hz, and 320 geophones were evenly distributed across the surface. The total recording time for the seismic data was 4 seconds, with a sampling interval of 2 ms. A local optimization algorithm was used to update and iterate the velocity model parameters, obtaining the inversion results for the low wavenumber components of the model.

[0069] This invention uses the Marmousi velocity model for numerical testing. First, raw seismic data is obtained through forward simulation using the finite difference wave equation, which is then used as observation data in the inversion process. Then, following the procedure outlined in the claims, an initial velocity model is set based on geological information, well logging information, and travel time information of the work area, serving as the initial model for the exponentially decaying time-frequency domain phase inversion method. An exponentially decaying time-frequency domain phase inversion objective function is constructed, and the partial derivatives of the objective function with respect to the velocity parameters are calculated to obtain the update amounts of the model parameters. A local optimization algorithm is then used to iterate the velocity model parameters, providing a better initial velocity model for conventional full-waveform inversion, obtaining higher-accuracy inversion results, and mitigating the period jump problem caused by the inversion process. The results of this invention are comprehensively compared with conventional FWI methods (…). Figure 3 Time-frequency domain phase inversion method ( Figure 4 ), based on the exponentially decaying time-frequency domain phase inversion method ( Figure 5 ).

[0070] Model parameters:

[0071] Table 1: Test parameters of the time-frequency domain phase inversion method based on exponential decay

[0072]

[0073] contrast Figure 3 a and Figure 4 It can be observed that the low-frequency time-domain phase inversion method yields better inversion results than the conventional FWI method. In the imaging results of the Marmousi model, the large-scale stratigraphic structure imaging results on the left side are more uniform. Further comparison... Figure 4 and Figure 5 By adding an attenuation factor to the time-frequency domain phase inversion, and adjusting the magnitude of the attenuation factor, the shallow part of the model can be inverted to the deep part of the model, achieving a layered inversion effect. This can provide a better and more uniform initial model for conventional FWI. Figure 5c) This method effectively alleviates the period jump problem in FWI. Phase information contains the kinematic characteristics of seismic wave propagation and can reveal information such as the propagation path of seismic waves and the subsurface medium. However, time-frequency domain phase inversion has certain noise effects during the inversion process. By introducing an attenuation factor, the period jump problem in conventional FWI is mitigated. Numerical tests verify that this method has certain advantages in subsurface velocity modeling.

Claims

1. A time-frequency domain phase inversion method based on exponential decay, comprising: Step 1: Use the Crews toolkit to preprocess the acquired seismic data as observation data; Step 2: Based on the geological information, well logging information and travel time information of the work area, set the initial velocity model for inversion; Step 3: Read the observation system, input the source wavelet, perform forward modeling using the initial velocity model, and store the wavefield values ​​at each sampling time point for subsequent cross-correlation calculation of the updated gradient of the velocity model. The study is based on the constant density acoustic wave equation. ; Where u represents the wave field value, z and x are the longitudinal and transverse coordinates respectively, t is time, and f(t) is the source wavelet; Step 4: Perform a short-time Fourier transform on the seismic data using the Gabor transform. The Gabor transform is defined as follows: ; u(τ) represents time-domain seismic data, h represents the Gaussian window function, and ω is the frequency. F represents the time-frequency domain representation of seismic data after Gabor transform. h [·] denotes the Gabor transform operator; The inverse Gabor transform is defined as: ; Where F h -1 [·] denotes the inverse Gabor transform operator; Step 5: The time-frequency domain signal of the seismic data is: ; The exponential phase of time-frequency domain seismic data is: ; in Represents the time-frequency domain envelope of seismic data; Step 6: Construct the time-frequency domain phase inversion objective function based on exponential decay: The objective function of a standard FWI is expressed as: ; Where u and d represent simulated data and actual observation data, respectively. The objective function for time-frequency domain phase inversion is: ; These represent the time-frequency phase of the simulated data and the observed data, respectively. The objective function for time-frequency domain phase inversion based on exponential decay is expressed as: ; Where a is the attenuation factor; Step 7: Calculate the partial derivatives of the objective function with respect to the model velocity parameters and the associated seismic source: The partial derivative of a conventional FWI with respect to the velocity parameter is: ; The corresponding accompanying epicenter is R s =(ud), the partial derivative of the time-frequency domain phase inversion objective function with respect to the velocity parameter is: ; Where Re[·] represents taking the real part of the seismic data, and * represents complex conjugate, which can be further simplified according to the chain rule to obtain: ; in Further simplification yields: ; The adjoint source of the objective function of conventional time-frequency domain phase inversion is expressed as: ; Based on the above derivation, the partial derivative of the time-frequency domain phase inversion objective function based on exponential decay with respect to the velocity parameter is expressed as: ; Further simplification using the chain rule yields: ; in Further simplification: ; The adjoint seismic source for the time-frequency domain phase inversion objective function based on exponential decay is: ; Step 8: Calculate the direction of velocity model updates: The gradient of the objective function in a conventional FWI model is directly obtained by the zero-delay cross-correlation between the forward-modeled wavefield and the residual backpropagated wavefield, expressed as: ; Where Nt represents the number of samples in the time domain. u represents the second-order partial derivative of the wavefield in the forward modeling with respect to time. b Represents the residual backpropagation wave field; Step 9: Update and iterate the velocity parameters using the steepest descent method. The update and iteration formula is as follows: ; Where α k Indicates the update step size, v k+1 This indicates that the velocity model v from the previous step... k Update the obtained speed parameters. This represents the update gradient of the velocity model; Step 10: Determine whether the inversion result of each step satisfies the termination condition set for the objective function. If it does, output the inversion result; if it does not, use the current inversion result as the initial velocity model for the next inversion and continue iterating until the termination condition is met. The inversion process uses the objective function as a constraint, and the iteration process needs to satisfy the continuous decrease of the objective function value.

Images