Stable Q compensation reverse time migration method for seismic exploration data, medium and equipment

Abstract

The invention relates to a stable Q compensation reverse time migration method for seismic exploration data, a medium and equipment, and belongs to the technical field of geophysical exploration data processing. According to the method, based on a fractional order Laplace viscous acoustic wave equation, a physically guided time-varying regularization compensation operator is constructed in a wavenumber domain, the energy amplitude of a previous time step length is introduced as a reference, an energy growth monitoring mechanism is established, and the energy growth is monitored. And designing a three-level adaptive constraint system: local energy growth constraint, high-frequency adaptive attenuation and boundary region stabilization processing, and realizing adjustment of excessive compensation and adaptive suppression of noise in the wave field continuation process. According to the method, amplitude attenuation and phase frequency dispersion caused by stratum absorption are effectively recovered, meanwhile, the stability of numerical calculation can be guaranteed, the precision and resolution of deep structure imaging are remarkably improved, and good numerical stability and anti-noise capacity are achieved.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of exploration geophysics, and in particular to a stable Q-compensated reverse-time migration method, medium and device for seismic exploration data. BACKGROUND

[0002] When seismic waves propagate in the subsurface medium, they will be affected by the viscoelastic absorption and scattering of the medium, exhibiting a dual effect of amplitude attenuation and phase dispersion, which leads to reduced resolution and waveform distortion of deep seismic records. In reverse-time migration (RTM), if this attenuation effect is not effectively compensated, it will seriously restrict the accuracy and reliability of deep structure imaging. Therefore, building a reasonable viscoelastic compensation model is of great significance to improve the imaging quality of seismic data and the ability to identify reservoirs.

[0003] Reverse-time migration technology is based on the two-way wave equation as the theoretical basis, and can adapt to the precise imaging needs of complex structures. However, the traditional reverse-time migration method is usually based on the assumption of a completely elastic medium, ignoring the influence of stratum absorption on seismic waves. In practical applications, this neglect will lead to reduced imaging resolution, especially in deep and strongly attenuated areas. For this reason, viscoelastic reverse-time migration methods have emerged, which introduce a damping compensation term into the wave equation to restore the amplitude and phase characteristics of seismic waves.

[0004] The widely used viscoelastic compensation methods currently include inverse Q filtering, attenuation-compensated wave equation, etc. These methods can theoretically restore the absorption attenuation of the stratum, but in practical applications, they generally have problems such as high-frequency noise amplification and numerical instability. In particular, compensation methods based on fractional Laplace operators can separate amplitude attenuation and phase dispersion, but the high-frequency components grow exponentially during the compensation process, which easily leads to wave field divergence and affects the final imaging quality.

[0005] To suppress the high-frequency divergence during the compensation process, the traditional solution mainly uses window function filtering methods, such as Turkey window, Butterworth low-pass filter, etc. These methods maintain numerical stability by hard-cutting or smoothly transitioning high-frequency components. However, the window function method usually uses fixed parameters, which is difficult to adapt to the spatial variation characteristics of the subsurface medium, and the hard-cutting operation will introduce Gibbs phenomenon, causing loss of effective high-frequency signals and reducing imaging resolution.

[0006] With the development of viscoelastic compensation theory, researchers have begun to explore stable compensation strategies based on regularization ideas. This type of method attempts to balance compensation attenuation and noise suppression by adding regularization terms to the compensation equation or using a prediction-correction framework. However, existing regularization methods mostly use fixed regularization parameters, which cannot be adaptively adjusted according to local medium characteristics and wave field propagation states, and still have problems of insufficient compensation or excessive smoothing under complex geological conditions.

[0007] Adaptive regularization method represents the latest research direction in the field of viscoacoustic compensation. This method integrates medium physical properties (such as quality factor Q, velocity) into the design of regularization operator, so that the compensation process can dynamically adjust according to the local attenuation characteristics. Introducing adaptive regularization into the viscoacoustic reverse-time migration imaging process can effectively overcome the limitations of the traditional window function method, while ensuring numerical stability and preserving the effective signal frequency components to the greatest extent, significantly improving the deep imaging quality.

[0008] However, the existing adaptive regularization method still faces challenges in parameter selection, computational efficiency and adaptability. How to construct an energy growth control mechanism with clear physical meaning, strong adaptability and high computational efficiency, and realize accurate matching of compensation strength and medium attenuation characteristics and wave field evolution trend, is a key problem in the field of viscoacoustic reverse-time migration imaging that needs to be solved. Therefore, it is necessary to develop a compensation method that can dynamically monitor energy growth, adaptively adjust the constraint strength according to local medium properties and wave number characteristics, and has good numerical stability and physical self-consistency, so as to better meet the needs of high-precision viscoacoustic reverse-time migration imaging under complex geological conditions. SUMMARY

[0009] In view of the shortcomings of the prior art, the present application provides a new wave number domain energy growth control stable Q compensation reverse-time migration method. This method combines fractional Laplace compensation equation, adaptive constraint mechanism based on energy growth control and three-level linkage physical guidance strategy, and constructs a stable, efficient and adaptable viscoacoustic compensation reverse-time migration process.

[0010] The present application is realized by the following technical solutions: A stable Q compensation reverse-time migration method for seismic exploration data, the method comprising the following steps: S1: obtaining the velocity model of the target work area, the quality factor Q model representing the absorption and attenuation characteristics of the stratum, and the source wavelet; based on the main frequency of the source wavelet and the velocity model, calculating the reference wave number that varies with space ; S2: taking the fractional Laplace viscoacoustic wave equation as the physical model, performing forward time extension on the seismic wave field representing the source excitation, and storing wave field snapshots at time steps during the extension process; during the extension process, to compensate for the amplitude attenuation of the seismic wave caused by the absorption of the stratum, an adaptive regularization weight function is introduced in the wave number domain for the amplitude compensation term The constraint is performed to suppress the exponential growth of numerical instability and high-frequency noise caused by compensation. The improved physical model is shown in equation (4); wherein the constraint strength of the weight function is dynamically adjusted according to the energy growth rate g(k, t) of the seismic wave field in the wave number domain, and the quality factor Q, the wave number k and the time t; ; wherein, is the wave field, is the time, is the propagation speed of the wave in the medium, is the Laplace operator, is the fractional order, and correspond to the coefficients of the phase distortion term and the amplitude attenuation term respectively, and F and represent the Fourier transform and the inverse transform respectively; S3: using the same physical model and constraint mechanism as step S2, the receiver wave field received by the ground is reversely time-extended; S4: in the process of reverse extension of the receiver wave field in step S3, it is cross-correlated with the source wave field snapshot at the corresponding time stored in step S2 to obtain the migration imaging result reflecting the underground structure; S5: post-processing the migration imaging result to output the final migration profile for geological interpretation.

[0011] Further, the energy growth rate in step S2 is: ; wherein, to prevent the denominator from being zero, a term is added, and ; Further, in step S2, the wave number domain representation of the adaptive regularization weight function is: ; For the amplitude compensation term, the weight function is applied: ; wherein, is the wave number domain representation of the amplitude compensation term after regularization constraint.

[0012] Further, the adaptive regularization weight function contains a first constraint term ; for energy growth constraint, when is greater than a preset threshold , an exponential constraint is applied: ; wherein is the regularized constraint term of energy growth restriction, is the constraint strength; The constraint strength is determined by three factors, and the specific formula is: ; wherein, is the basic strength, and its expression is: ; is the time factor, and its expression is: ; is the wave number factor, and its expression is: ; In the above formula, is the reference medium factor, generally taking the average value of Q; is an intensity factor, which can be taken as 0.001-0.002, is the wave number, is the reference wave number; For components much higher than the reference wave number ( >5 ) to apply additional attenuation: ; wherein, is the high-frequency regularized constraint term, is the recording time length; At the same time, an exponential decay is introduced near the wave number domain boundary (normalized wave number > or > ): ; wherein, is the regularized constraint term of the boundary area, controls the steepness of the decay, and the threshold is usually taken as 0.8-0.9; is the wave number component in the x direction, is the wave number component in the z direction; the normalized wave number is defined as: , ; wherein: , .

[0013] Further, the cross-correlation calculation in step S4: ; wherein, is the migration result of the offset, is a source wave field, is a receiver wave field.

[0014] The application further provides a computer readable storage medium, which stores a computer program, and the computer program is suitable for being loaded and executed by a processor to perform the stable Q compensation reverse time migration method for seismic exploration data.

[0015] The application further provides a computer device, which comprises a memory and a processor, the memory stores a computer program, and the computer program is executed by the processor to make the processor perform the stable Q compensation reverse time migration method for seismic exploration data.

[0016] Compared with the prior art, the application has the following beneficial effects: The application develops a stable Q compensation method for viscoacoustic reverse time migration based on energy growth control, which can maintain numerical stability in long-time wave field continuation and effectively restore the energy of deep attenuation signals.

[0017] The application innovatively introduces an energy growth monitoring mechanism into the wave number domain compensation process, calculates the amplitude ratio of each wave number component at the front and rear time as an energy growth rate in real time, and sets a dynamic trigger threshold. For the compensation operator exceeding the threshold, the application designs a hierarchical response strategy: an adaptive attenuation coefficient is calculated according to the exceeding degree, the wave number height and the local Q value size, the more the exceeding degree, the higher the frequency and the smaller the Q value, the greater the attenuation intensity; the normal growth component is basically retained. This mechanism effectively suppresses high-frequency noise without affecting effective signals.

[0018] In addition, the application introduces a spatial adaptive adjustment mechanism based on local medium parameters in the wave field continuation process. The traditional compensation method uses a global uniform constraint parameter, and when processing transversely non-uniform media, it often has to choose one of the two - low Q strong attenuation area needs strong compensation, but too strong constraint will lead to deep signal cannot be effectively restored; the high Q weak attenuation area is insufficiently constrained and will amplify noise. The method directly couples the quality factor, velocity and other local parameters into the constraint strength function, so that the low Q area automatically reduces the constraint to ensure the compensation strength, and the high Q area automatically enhances the constraint to prevent excessive amplification, thereby realizing the compensation effect of considering both shallow and deep layers and stable whole profile. BRIEF DESCRIPTION OF DRAWINGS

[0019] Figure 1 is a workflow diagram of the stable Q compensation reverse time migration method based on wave number domain energy growth control; Figure 2 is a velocity model diagram; Figure 3 is a Q model diagram; Figure 4is a viscous acoustic single-shot record diagram being played forward; Figure 5 is a compensated viscous acoustic seismic record diagram; Figure 6 is a wave field snapshot diagram at 800 ms before compensation; Figure 7 is a wave field snapshot diagram at 800 ms after compensation; Figure 8 is a viscous acoustic data uncompensated viscous acoustic reverse time migration result; Figure 9 is a viscous acoustic data compensated viscous acoustic reverse time migration result; Figure 10 is a comparison diagram of amplitude curves before and after compensation at a horizontal direction of 2 km. DETAILED DESCRIPTION

[0020] In order to make the objects, technical solutions and advantages of the present application clearer, the present application will be further described in detail below with examples. It should be understood that the specific examples described herein are only used to explain the present application, and are not used to limit the present application.

[0021] The present application provides a stable Q compensation reverse time migration method based on wave number domain energy growth control, which successfully realizes high-resolution and high-fidelity imaging of seismic data in attenuating medium. The method is based on fractional order Laplacian viscous acoustic wave equation as a theoretical basis, realizes attenuation compensation by changing the sign of amplitude attenuation term, and innovatively constructs an adaptive regularization weight function related to medium properties (Q value, velocity), wave number and time in the wave number domain, and constructs a physically guided stable compensation equation. The weight function realizes intelligent suppression of high-frequency noise and maximum reservation of effective signals by real-time monitoring of wave field energy growth rate and dynamically triggering a three-level adaptive constraint mechanism (energy growth constraint, high-frequency attenuation, boundary stability). Subsequently, the forward and reverse continuation of the source wave field and the receiver wave field are carried out by using the regularization compensation wave equation, and the adaptive weight function is applied in the wave number domain to realize accurate regulation and control of the compensation process, and finally the high-resolution migration profile is obtained through the cross-correlation imaging condition, so that stable and efficient reverse time migration imaging in complex viscous medium is realized. The present application realizes high-precision reverse time migration imaging for complex velocity model, which will be described in detail below with reference to the accompanying drawings.

[0022] As shown in Figure 1 , the specific steps of the stable Q compensation reverse time migration method based on wave number domain energy growth control proposed by the present application are as follows: Step S10: calculating the reference wave number based on the source main frequency and the velocity. Input the velocity model, and the Q model is as shown in Figure 2 、 Figure 3As shown, the time step was 1 ms, the spatial step was 10 m, the model size was 1610 m × 3980 m, the wavelet used was the Ricker wavelet with a dominant frequency of 20 Hz, and the recording duration was 3.0 s. A total of 40 shots were uniformly arranged at a depth of 10 m, with offsets starting at 20 m and firing one shot every 100 m. The spatially varying reference wavenumber was: ; in, The dominant frequency of the earthquake source, For speed, It is the reference wavenumber.

[0023] Step S20: Construction of the adaptive regularized weighting function based on energy growth control and forward continuation of source wavefield attenuation compensation. When the source wavefield is forward continuated (from t=0 to T, where T is the recording duration) using the fractional-order Laplace viscous acoustic wave equation, the fractional-order Laplace viscous acoustic wave equation is: ; The first term on the right side of formula (2) is the phase distortion term, and the second term is the amplitude attenuation term. Wherein, For wave field, It is time. It is the speed at which a wave propagates in a medium. For the Laplace operator, It is the fractional order. and The coefficients corresponding to the phase distortion term and the amplitude attenuation term are expressed as follows: ; in, This is a reference speed. It is a reference frequency.

[0024] To compensate for the attenuation effect of formation absorption, it is only necessary to change the sign of the amplitude attenuation term in formula (2) while keeping the sign of the phase distortion term unchanged to achieve attenuation compensation. During the compensation process, the wave field is unstable; therefore, a term related to the medium properties and wavenumber is usually introduced into the compensation term. and time The relevant weighting function W yields the regularized viscous acoustic compensation wave equation: ; Among them, F and These represent the Fourier transform and the inverse transform, respectively.

[0025] To address the shortcomings of existing regularization methods that often employ fixed regularization parameters, this invention proposes an adaptive regularization weight function based on energy growth control, which better solves the problem of overcompensation under complex geological conditions. The specific construction process is as follows: At each time step, the amplitude compensation term is transformed to the wavenumber domain, and the wavenumber domain amplitude is calculated: ; in, It is the real part of the complex wave field. It is the imaginary part of the complex wave field. It is the wavenumber domain representation of the second term on the right side of formula (4).

[0026] Energy growth rate for: ; To prevent the denominator from being zero, a was added. Item, and Let be the wavenumber domain representation of the amplitude compensation term at time t. Let be the wavenumber domain representation of the amplitude compensation term at time t.

[0027] Regarding the energy growth constraint, when hour, In this embodiment, a set threshold is used. The value is 1.2, applying an exponential constraint: ; in It is a regularization constraint term of the energy growth constraint. It refers to the constraint strength.

[0028] The constraint strength is determined by three factors, and the specific formula is as follows: ; in, It is the basic strength, and its expression is: ; It is a time factor, and its expression is: ; It is the wavenumber factor, and its expression is: .

[0029] In the above formula, As a reference medium factor, the average value of Q is generally taken; This is an intensity factor, which can be set between 0.001 and 0.002. For wave number, The reference wavenumber is used.

[0030] For components much higher than the reference wavenumber ( >5 Apply additional attenuation: ; in, It is a high-frequency regularization constraint term. It records the duration.

[0031] Meanwhile, near the boundary of the wavenumber domain (normalized wavenumber) > or > Introducing exponential decay: ; in, It is a regularization constraint term for the boundary region. Control the steepness of decay, threshold The value is usually set between 0.8 and 0.9. It is the wavenumber component in the x-direction. This refers to the wavenumber component in the z-direction. Define the normalized wavenumber: , ; in: , ; Therefore, the wavenumber domain representation of the adaptive regularized weighting function proposed in this invention is: ; Applying a weighting function to the amplitude compensation term: ; in, It is the wavenumber domain representation of the amplitude compensation term after regularization constraints.

[0032] The wavefield is inversely transformed back to the spatial domain, and a snapshot of the wavefield is stored for subsequent imaging.

[0033] The viscous acoustic wave single-shot record currently in operation is as follows: Figure 4 As shown, the compensated viscous acoustic seismic record is as follows: Figure 5 As shown in the figure, comparing the two images clearly reveals that the uncompensated viscous acoustic wave forward modeling record exhibits weak deep signals and low resolution. After attenuation compensation processing, the amplitude and frequency components of the seismic waves are effectively recovered, the reflection energy in the mid-deep region is enhanced, the continuity of the phase axis is improved, and the seismic time resolution is significantly improved, laying a higher quality data foundation for subsequent high-precision imaging and inversion. In addition, Figure 6 This shows a snapshot of the wave field 800ms before compensation. Figure 7A snapshot of the wavefield at 800ms after compensation is shown. Before compensation, the wavefield exhibits rapid energy decay, wavefront blurring and dispersion, and insufficient deep illumination; while after attenuation compensation, the wavefield shows that the energy is effectively recovered spatially, and the illumination of deep and shadow areas is significantly enhanced.

[0034] Step S30: The wave field at the receiver point is extended in reverse (from t=T to 0) using the same regularized compensation wave equation. The same adaptive regularized weight function is applied in the wavenumber domain to consider the reverse time evolution characteristics and ensure the numerical stability of the compensation process.

[0035] Step S40: During the reverse extension of the receiver wavefield, at each time step, read the stored snapshot of the source wavefield at the corresponding time and perform cross-correlation calculation: ; in, For the offset imaging results, For the source wave field, The wave field at the detector point.

[0036] Step S50: Post-processing of imaging results. After all calculations are completed, the final offset imaging result is output and post-processed, such as Laplacian filtering to remove low-frequency noise and amplitude equalization, to improve the quality of the imaging result.

[0037] picture Figure 8 and Figure 9 The results of viscous acoustic reverse time migration (RTM) imaging before and after viscous acoustic data compensation are presented. Comparing the viscous acoustic RTM results, it can be seen that uncompensated viscous acoustic RTM imaging is limited by attenuation effects, exhibiting weak deep energy, low resolution, and structural blurring, severely restricting the potential for deep exploration. However, after viscous attenuation compensation, the RTM imaging quality is systematically improved: deep amplitude is recovered, profile energy is more balanced, vertical resolution is significantly improved, and deep structures and stratigraphic details are clearly distinguishable. Figure 10 The amplitude comparison curves at a horizontal distance of 2 km are shown. The results demonstrate that, compared to the uncompensated imaging, the amplitude energy of the viscosity-acoustic compensated reverse time migration imaging is significantly recovered, effectively improving the resolution and amplitude preservation of the imaging profile. Therefore, the stable Q-compensated reverse time migration method based on wavenumber domain energy growth control proposed in this invention can more precisely characterize the underground velocity structure, possessing high resolution and practical application value.

[0038] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any modifications, equivalent substitutions, and improvements made by those skilled in the art within the scope of the technology disclosed in the present invention, and within the spirit and principles of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A stable Q-compensated reverse-time migration method for seismic exploration data, characterized in that, The method includes the following steps: S1: Obtain the velocity model of the target work area, the quality factor Q model characterizing the formation absorption and attenuation properties, and the source wavelet; based on the dominant frequency of the source wavelet and the velocity model, calculate the spatially varying reference wavenumber. ; S2: Using the fractional-order Laplace viscous wave equation as the physical model, a positive time extension is performed on the seismic wavefield characterizing the source-generated wavefield, and a time-step snapshot of the wavefield is stored during the extension process. To compensate for the amplitude attenuation of the seismic waves caused by formation absorption during the extension process, an adaptive regularization weighting function is introduced in the wavenumber domain for the amplitude compensation term. Constraints are applied to suppress numerical instability caused by compensation and the exponential growth of high-frequency noise. The improved physical model is shown in Equation (4). The constraint strength of the weight function is adaptively and dynamically adjusted according to the energy growth rate g(k,t) of the seismic wave field in the wavenumber domain, as well as the quality factor Q, wavenumber k and time t. ; in, For wave field, It is time. It is the speed at which a wave propagates in a medium. For the Laplace operator, It is the fractional order. and The coefficients F and F, corresponding to the phase distortion term and amplitude attenuation term, respectively. These represent the Fourier transform and the inverse transform, respectively. S3: Using the same physical model and constraint mechanism as in step S2, perform reverse time extension on the receiver wavefield received from the Earth's surface; S4: During the reverse extension of the receiver wavefield described in step S3, cross-correlation calculation is performed between it and the corresponding moment source wavefield snapshot stored in step S2 to obtain the migration imaging result reflecting the subsurface structure. S5: Post-process the migration imaging results to output the final migration profile for geological interpretation.

2. The method for stable Q-compensated reverse-time migration of seismic exploration data according to claim 1, characterized in that, The energy growth rate in step S2 for: ； in, , For the wavenumber domain representation of the amplitude compensation term, Let be the wavenumber domain representation of the amplitude compensation term at time t.

3. The method for stable Q-compensated reverse-time migration of seismic exploration data according to claim 2, characterized in that, In step S2, the wavenumber domain representation of the adaptive regularization weight function is as follows: ； The amplitude compensation term applies a weighting function: ； in, It is the wavenumber domain representation of the amplitude compensation term after regularization constraints.

4. The method for stable Q-compensated reverse-time migration of seismic exploration data according to claim 3, characterized in that, The adaptive regularization weight function Include Regarding the energy growth constraint, when Greater than the preset threshold When applying exponential constraints: ； in, It is a regularization constraint term of the energy growth constraint. It is the constraint strength; The constraint strength is determined by three factors, and the specific formula is as follows: ； in, It is the basic strength, and its expression is: ; It is a time factor, and its expression is: ; It is the wavenumber factor, and its expression is: ; In the above formula, As a reference medium factor, the average value of Q is taken; This is an intensity factor, with a value ranging from 0.001 to 0.

002. For wave number, Reference wavenumber; for >5 Apply additional attenuation to the components: ； in, It is a high-frequency regularization constraint term. It records the duration; Meanwhile, near the boundary of the wavenumber domain, i.e., the normalized wavenumber > or > When exponential decay is introduced: ； in, It is a regularization constraint term for the boundary region. Control the steepness of decay, threshold Take 0.8-0.9; It is the wavenumber component in the x-direction. It is the wavenumber component in the z-direction; define the normalized wavenumber: ， ； in: , 。 5. A stable Q-compensated reverse-time migration method for seismic exploration data according to claim 4, characterized in that, The cross-correlation calculation described in step S4: ； in, For the offset imaging results, For the source wave field, The wave field at the detector point.

6. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program adapted to be loaded by a processor and executed by the stable Q-compensated reverse time migration method for seismic exploration data as described in any one of claims 1-5.

7. A computer device, the device comprising a memory and a processor, the memory storing a computer program, characterized in that, When the computer program is executed by the processor, the processor performs the stable Q-compensated reverse time migration method for seismic exploration data as described in any one of claims 1-5.

Images