A stable q-compensated reverse time migration method, medium and device for seismic exploration data

Abstract

The present application relates to a kind of stable Q compensation reverse time migration method, medium and equipment for seismic exploration data, belong to geophysical exploration data processing technical field.The method is based on fractional order Laplace viscous wave wave equation, constructs physical guide time-varying regularization compensation operator in wave number domain, by introducing the energy amplitude of previous time step as reference benchmark, establishes energy growth monitoring mechanism, and designs three-level adaptive constraint system: local energy growth constraint, high-frequency adaptive attenuation and boundary region stabilization processing, in wave field continuation process, realize the adjustment of overcompensation and the adaptive suppression of noise.The method can effectively restore the amplitude attenuation and phase dispersion caused by stratum absorption while ensuring the stability of numerical calculation, significantly improve the accuracy and resolution of deep structure imaging, with good numerical stability and noise resistance.

Description

Technical Field

[0001] This invention relates to the field of exploration geophysics, and more particularly to a stable Q-compensated reverse time migration method, medium, and apparatus for seismic exploration data. Background Technology

[0002] When seismic waves propagate through subsurface media, they are subjected to viscoelastic absorption and scattering, resulting in a dual effect of amplitude attenuation and phase dispersion. This leads to reduced resolution and waveform distortion in deep seismic records. In reverse time migration (RTM) imaging, if this attenuation effect is not effectively compensated, it will severely limit the accuracy and reliability of deep structural imaging. Therefore, constructing a reasonable viscoelastic-acoustic compensation model is of great significance for improving the imaging quality of seismic data and reservoir identification capabilities.

[0003] Reverse-time migration (RTM) is based on the two-way wave equation and can meet the precise imaging requirements of complex structures. However, traditional RTM methods are usually based on the assumption of a perfectly elastic medium, neglecting the influence of formation absorption on seismic waves. In practical applications, this neglect leads to reduced imaging resolution, especially in deep layers and regions with strong attenuation. To address this, the viscous acoustic RTM method was developed, which aims to recover the amplitude and phase characteristics of seismic waves by introducing attenuation compensation terms into the wave equation.

[0004] Currently widely used methods for viscosity-acoustic compensation include inverse Q filtering and attenuation compensation wave equation methods. While these methods can theoretically recover the absorption attenuation of the formation, they generally suffer from problems such as high-frequency noise amplification and numerical instability in practical applications. In particular, compensation methods based on the fractional-order Laplace operator, although capable of separating amplitude attenuation from phase dispersion, exhibit exponential growth of high-frequency components during the compensation process, easily leading to wavefield divergence and affecting the final imaging quality.

[0005] To suppress high-frequency divergence during the compensation process, traditional solutions mainly employ window function filtering methods, such as the Turkey window and the Butterworth low-pass filter. These methods maintain numerical stability by hard-truncating or smoothing the transition of high-frequency components. However, window function methods typically use fixed parameters, making it difficult to adapt to the spatial variations of the subsurface medium. Furthermore, hard truncation introduces the Gibbs phenomenon, resulting in the 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. These methods attempt to strike a balance between compensation attenuation and noise suppression by adding regularization terms to the compensation equation or employing a predictive-correction framework. However, existing regularization methods mostly use fixed regularization parameters, which cannot adaptively adjust according to local medium characteristics and wavefield propagation states. Under complex geological conditions, they still suffer from insufficient compensation or excessive smoothing.

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

[0008] However, existing adaptive regularization methods still face challenges in parameter selection, computational efficiency, and adaptability. How to construct a physically meaningful, highly adaptive, and computationally efficient energy growth control mechanism to achieve precise matching between compensation intensity and medium attenuation characteristics and wavefield evolution trends is a key problem urgently needing to be solved in the field of viscous acoustic reverse time migration imaging. Therefore, it is necessary to develop a compensation method that can dynamically monitor energy growth, adaptively adjust constraint intensity according to local medium properties and wavenumber characteristics, and possess good numerical stability and physical self-consistency, making it more suitable for the needs of high-precision viscous acoustic reverse time migration imaging under complex geological conditions. Summary of the Invention

[0009] To address the shortcomings of existing technologies, this invention provides a novel stable Q-compensated reverse-time migration method based on wavenumber domain energy growth control. This method combines a fractional-order Laplace compensation equation, an adaptive constraint mechanism based on energy growth control, and a three-level linkage physical guidance strategy to construct a stable, efficient, and highly adaptable viscous acoustic compensation reverse-time migration process.

[0010] This invention is achieved through the following technical solution:

[0011] A stable Q-compensated reverse-time migration method for seismic exploration data, the method comprising the following steps:

[0012] S1: Obtain the velocity model of the target work area, the quality factor Q model characterizing the formation absorption and attenuation characteristics, and the source wavelet; based on the dominant frequency of the source wavelet and the velocity model, calculate the spatially varying reference wavenumber. ;

[0013] S2: Using the fractional-order Laplace viscous wave equation as the physical model, a positive time extension is performed on the seismic wave field characterizing the source-generated wave field, and a time-step snapshot of the wave field 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 weighting function is adaptively and dynamically adjusted according to the energy growth rate g(k, t) of the seismic wavefield in the wavenumber domain, as well as the quality factor Q, wavenumber k, and time t. ;

[0014] 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.

[0015] 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;

[0016] 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.

[0017] S5: Post-process the migration imaging results to output the final migration profile for geological interpretation.

[0018] Furthermore, the energy growth rate in step S2 for:

[0019] ;

[0020] To prevent the denominator from being zero, a was added. Item, and ;

[0021] Furthermore, in step S2, the wavenumber domain representation of the adaptive regularization weight function is:

[0022] Applying a weighting function to the amplitude compensation term:

[0023] ;

[0024] in, It is the wavenumber domain representation of the amplitude compensation term after regularization constraints.

[0025] Furthermore, the adaptive regularization weight function Includes the first constraint term Regarding the energy growth constraint, when When applying exponential constraints:

[0026] ;

[0027] in It is a regularization constraint term of the energy growth constraint. It is the constraint strength;

[0028] The constraint strength is determined by three factors, and the specific formula is as follows:

[0029] ;

[0030] in, It is the basic strength, and its expression is: ;

[0031] It is a time factor, and its expression is: ;

[0032] It is the wavenumber factor, and its expression is: ;

[0033] 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, Reference wavenumber;

[0034] For components much higher than the reference wavenumber ( >5 Apply additional attenuation:

[0035] ;

[0036] in, It is a high-frequency regularization constraint term. It records the duration;

[0037] Meanwhile, near the boundary of the wavenumber domain (normalized wavenumber) > or > Introducing exponential decay:

[0038] ;

[0039] in, Control the steepness of decay, threshold The value is usually taken as 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:

[0040] , ;in:

[0041] , .

[0042] Furthermore, the cross-correlation calculation described in step S4:

[0043] ;

[0044] in, For the offset imaging results, For the source wave field, The wave field at the detector point.

[0045] The present invention also provides a computer-readable storage medium storing a computer program adapted to be loaded and executed by a processor, the computer program being used for the stable Q-compensated reverse time migration method for seismic exploration data.

[0046] The present invention also provides a computer device, the device including a memory and a processor, the memory storing a computer program, which, when executed by the processor, causes the processor to perform the stable Q-compensated reverse time migration method for seismic exploration data.

[0047] The beneficial effects of this invention compared to the prior art are as follows:

[0048] This invention develops a viscoacoustic reverse time-shift stable Q-compensation method based on energy growth control, which can maintain numerical stability in long-term wave field extension and effectively recover the energy of deep attenuated signals.

[0049] This invention innovatively introduces an energy growth monitoring mechanism into the wavenumber domain compensation process. It calculates the amplitude ratio of each wavenumber component before and after a given time step as the energy growth rate and sets a dynamic trigger threshold. For compensation operators exceeding the threshold, this invention designs a hierarchical response strategy: an adaptive attenuation coefficient is calculated based on the degree of exceedance, wavenumber level, and local Q-value. The greater the exceedance, the higher the frequency, and the smaller the Q-value, the greater the attenuation intensity; normally growing components are essentially preserved. This mechanism effectively suppresses high-frequency noise while leaving the effective signal unaffected.

[0050] Furthermore, this invention introduces a spatially adaptive adjustment mechanism based on local medium parameters during wavefield extension. Traditional compensation methods use globally uniform constraint parameters, which often suffer from trade-offs when dealing with transversely inhomogeneous media—strong compensation is needed in low-Q, high-attenuation regions, but excessive constraints can lead to ineffective recovery of deep signals; insufficient constraints in high-Q, weak-attenuation regions can amplify noise. This method directly couples local parameters such as quality factor and velocity into the constraint strength function, allowing the constraints to automatically decrease in low-Q regions to ensure compensation strength, and the constraints to automatically increase in high-Q regions to prevent excessive amplification, thereby achieving a compensation effect that considers both shallow and deep layers and is stable across the entire profile. Attached Figure Description

[0051] Figure 1 This is a flowchart of the stable Q-compensated reverse-time migration method based on wavenumber domain energy growth control;

[0052] Figure 2 It is a velocity model diagram;

[0053] Figure 3 It is a Q-model diagram;

[0054] Figure 4 This is a single-shot recording of viscous acoustic waves during a live performance;

[0055] Figure 5 This is a compensated viscous acoustic seismic record;

[0056] Figure 6 This is a snapshot of the wave field 800ms before compensation;

[0057] Figure 7 It is a snapshot of the wave field at 800ms after compensation;

[0058] Figure 8 The results of the adhesive acoustic wave data without compensation for the reverse time shift of the adhesive acoustic wave;

[0059] Figure 9 It is the result of compensating for the reverse time migration of the acoustic wave data;

[0060] Figure 10 This is a comparison chart of amplitude curves before and after compensation at a horizontal distance of 2km. Detailed Implementation

[0061] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to examples. It should be understood that the specific examples described herein are merely illustrative and not intended to limit the invention.

[0062] This invention provides a stable Q-compensated reverse-time migration method based on wavenumber domain energy growth control, successfully achieving high-resolution, high-fidelity imaging of seismic data in attenuated media. The method is based on the fractional-order Laplace operator viscous acoustic wave equation, achieving attenuation compensation by changing the sign of the amplitude attenuation term. It innovatively constructs an adaptive regularized weighting function in the wavenumber domain, related to medium properties (Q-value, velocity), wavenumber, and time, thus building a physically guided stable compensation equation. This weighting function dynamically triggers a three-level adaptive constraint mechanism (energy growth constraint, high-frequency attenuation, and boundary stabilization) by real-time monitoring of the wavefield energy growth rate, achieving intelligent suppression of high-frequency noise and maximizing the preservation of effective signals. Subsequently, the regularized compensation wave equation is used to perform forward and reverse extensions of the source wavefield and receiver wavefield, respectively. The adaptive weighting function is applied in the wavenumber domain to achieve precise control of the compensation process. Finally, a high-resolution migration profile is obtained through cross-correlation imaging conditions, thereby achieving stable and efficient reverse-time migration imaging in complex viscous acoustic media. This invention achieves high-precision reverse-time migration imaging for complex velocity models. The invention will be described in detail below with reference to the accompanying drawings.

[0063] like Figure 1 As shown, the stable Q-compensated reverse-time migration method based on wavenumber domain energy growth control proposed in this invention has the following specific steps:

[0064] Step S10: Calculate the reference wavenumber based on the dominant frequency and velocity of the seismic source. Input the velocity model, Q model as follows: Figure 2 , Figure 3 As 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:

[0065] ;

[0066] in, The dominant frequency of the earthquake source, For speed, It is the reference wavenumber.

[0067] 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:

[0068] ;

[0069] 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:

[0070] ;

[0071] in, This is a reference speed. It is a reference frequency.

[0072] 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:

[0073] ;

[0074] Among them, F and These represent the Fourier transform and the inverse transform, respectively.

[0075] 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:

[0076] At each time step, the amplitude compensation term is transformed to the wavenumber domain, and the wavenumber domain amplitude is calculated:

[0077] ;

[0078] 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).

[0079] Energy growth rate for:

[0080] ;

[0081] 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.

[0082] Regarding the energy growth constraint, when hour, In this embodiment, a set threshold is used. The value is 1.2, applying an exponential constraint:

[0083] ;

[0084] in It is a regularization constraint term of the energy growth constraint. It refers to the constraint strength.

[0085] The constraint strength is determined by three factors, and the specific formula is as follows:

[0086] ;

[0087] in, It is the basic strength, and its expression is: ;

[0088] It is a time factor, and its expression is: ;

[0089] It is the wavenumber factor, and its expression is: .

[0090] 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.

[0091] For components much higher than the reference wavenumber ( >5 Apply additional attenuation:

[0092] ;

[0093] in, It is a high-frequency regularization constraint term. It records the duration.

[0094] Meanwhile, near the boundary of the wavenumber domain (normalized wavenumber) > or > Introducing exponential decay:

[0095] ;

[0096] in, 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:

[0097] , ;in:

[0098] , ;

[0099] Therefore, the wavenumber domain representation of the adaptive regularized weighting function proposed in this invention is:

[0100] ;

[0101] Applying a weighting function to the amplitude compensation term:

[0102] ;

[0103] in, It is the wavenumber domain representation of the amplitude compensation term after regularization constraints.

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

[0105] 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 7 A 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.

[0106] 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.

[0107] 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:

[0108] ;

[0109] in, For the offset imaging results, For the source wave field, The wave field at the detector point.

[0110] 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.

[0111] 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.

[0112] 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 characteristics, 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 weighting function is adaptively and dynamically adjusted according to the energy growth rate g(k, t) of the seismic wavefield 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: 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 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