A method and system for variational fractional viscoacoustic wave equation attenuation-compensated reverse time migration

Abstract

The application discloses a variable fractional order viscous acoustic wave equation attenuation compensation reverse time migration method and system, and the variable fractional order viscous acoustic wave equation attenuation compensation reverse time migration method comprises the following steps: a stable attenuation compensation viscous acoustic wave equation is constructed by introducing a regularization term; and based on the stable attenuation compensation viscous acoustic wave equation, a local cross-correlation imaging condition based on the Nyquist sampling theorem is used to carry out attenuation compensation reverse time migration of viscous acoustic waves. The application adopts the method of adding a regularization term, can effectively inhibit the exponential growth of high-frequency noise in seismic records in the compensation process, and ensures the numerical stability of the attenuation compensation process; the local cross-correlation imaging condition based on the Nyquist sampling theorem can effectively reduce the calculation and storage cost of Q compensation reverse time migration on the basis of adapting to complex structure imaging, and improve the practicability of the reverse time migration imaging of viscoelastic media.

Description

Technical Field

[0001] This invention belongs to the field of geological exploration technology, and specifically relates to a method and system for attenuation compensation of variable fractional-order viscous acoustic wave equations and reverse time migration. Background Technology

[0002] As oil and gas exploration targets become increasingly complex, refined, and deep, the requirements for seismic imaging technology are becoming more and more demanding. In the past decade, thanks to the rapid development of computer science computing, seismic data processing has shifted from post-stack to pre-stack, giving rise to seismic imaging theories and technologies represented by reverse time migration and full waveform inversion.

[0003] Reverse time migration (RTM), based on wave equations for wavefield extrapolation, is the most accurate imaging method currently available for seismic wave migration. Theoretically, it can achieve reflection repositioning and diffraction convergence for any type of wavefield, without being limited by the dip angle of the strata. Due to its unique imaging advantages for complex structures, it has become an indispensable key technology in the field of migration imaging.

[0004] In actual seismic exploration, seismic waves often encounter viscous strata, resulting in strong absorption effects, amplitude attenuation, and frequency scattering. This manifests in seismic records as weakened amplitude and waveform distortion, rendering traditional acoustic wave theory inadequate for describing the propagation of seismic waves in viscous media. To accurately describe the frequency-dependent absorption and attenuation properties of seismic wave propagation, researchers have proposed using fractional derivatives to describe the constitutive relations of viscoelastic media and combining this with the dispersion relation of the constant Q model to establish a novel fractional Laplace operator viscous acoustic wave equation. This equation is concise, and its characterization parameters include only seismic wave velocity, quality factor Q, and reference frequency. Forward modeling shows that this equation achieves high accuracy in simulating the constant Q model. Furthermore, the amplitude and phase terms in this equation are decoupled; wavefield compensation can be achieved simply by changing the sign of the amplitude term. Based on these advantages, the fractional Laplace operator viscous acoustic wave equation is currently widely used in seismic simulation, Q-compensated reverse time migration, and full waveform inversion.

[0005] Currently, absorption attenuation compensation methods can be divided into two main categories. The first category is attenuation compensation based on seismic traces, mainly including unsteady-state deconvolution and inverse Q-filtering. Although this type of method is simple and efficient to operate, it essentially ignores the fact that seismic wave attenuation occurs along the propagation path. The second category of attenuation compensation methods is based on pre-stack depth migration of the viscous wave equation, eliminating the influence of medium absorption during the wavefield extension process of seismic migration. This type of method implements Q-value compensation along the attenuation path of the seismic wave, which is theoretically more convincing than the first type of method. In the process of realizing attenuation compensation in the reverse time migration of the viscous wave equation, the inverse process of wavefield attenuation is used to compensate for the absorbed energy and correct the wavelet phase distortion caused by velocity dispersion. Since the formation absorption exhibits an exponential attenuation form, the direct compensation strategy is to exponentially amplify the attenuated wavefield. However, this will also amplify background noise such as random interference, leading to numerical instability in the compensation results. The main solutions currently include: filtering methods, stability factor methods, regularization methods, and least squares reverse time migration of the viscous acoustic equation. Among these methods, filtering inevitably causes high-frequency damage, and the selection of its cutoff frequency is often uncertain, leading to redundant calculations if chosen blindly. The stability factor method, which characterizes the attenuation compensation operator by including only the ratio of the dispersive wave field to the viscous wave field, does not involve any form of energy growth and is very effective in suppressing numerical instability in compensation. However, this method requires separate calculation of only the dispersive wave field, increasing the computational load of one forward modeling step. The least squares reverse time migration method indirectly avoids numerical instability by matching the forward modeling results with the observation records during the iteration process, but its enormous computational load is beyond the capacity of current hardware, especially for large-scale three-dimensional reverse time migration. Therefore, how to add regularization constraints to the unstable amplitude terms in the viscous acoustic wave equations during the compensation process has become a key research direction for many scholars.

[0006] Consistent with conventional reverse-time migration (RTM), Q-compensated RTM also consists of three parts: forward-extending of the source wavefield, reverse-time extension of the detector wavefield, and obtaining seismic imaging results using imaging conditions. Q-compensated RTM also faces challenges such as massive computation, huge storage requirements, and migration noise. Extensive research on migration noise has significantly improved this problem. However, computational and storage burdens remain a technical bottleneck restricting the large-scale application of RTM. To address this, some researchers have proposed excitation amplitude imaging conditions, using only certain feature points in the wavefield for imaging. This method offers high imaging resolution, avoids large storage requirements, and eliminates the need for wavefield reconstruction, making it highly practical. However, this imaging condition does not fully utilize wavefield information, and the imaging results are easily affected by the signal-to-noise ratio of the data. Conventional cross-correlation imaging conditions utilizing all information can achieve accurate imaging under complex geological conditions, but this comes with enormous computational and storage costs. Therefore, there is an urgent need to develop a variational fractional-order viscous acoustic equation attenuation-compensated RTM method. Summary of the Invention

[0007] To address the above problems, this invention discloses a method for attenuation compensation and reverse-time migration of a variational fractional-order viscous acoustic wave equation, comprising the following steps:

[0008] A stable attenuation-compensated viscous acoustic wave equation is constructed by introducing a regularization term;

[0009] Based on the viscous acoustic wave equation with stable attenuation compensation, the attenuation compensation and reverse time migration of viscous acoustic waves are performed using the local cross-correlation imaging condition based on the Nyquist sampling theorem.

[0010] Furthermore, the regularization term is determined by the following formula:

[0011]

[0012] Where ε is the regularization factor; c is the seismic wave velocity; τ is the control parameter for the amplitude term of the equation; p is the wave field; t is time; ▽ 2 γ is the Laplace operator; γ is the fractional parameter.

[0013] Furthermore, the viscous acoustic wave equation for stable attenuation compensation is determined by the following formula:

[0014]

[0015] Where η is the control parameter for the phase term of the equation.

[0016] Furthermore, the specific steps for establishing the local cross-correlation imaging conditions are as follows:

[0017] During the source wave field extension process, a time window is designed for each grid point. The source wave field within the set time window range is stored with the time corresponding to the maximum amplitude of the source wave field at each grid point as the center.

[0018] During the back propagation of the detector's wave field, the wave field values ​​corresponding to the time window range are extracted one by one.

[0019] Reverse time migration imaging is performed based on the source wavefield and detector wavefield values ​​within the time window range.

[0020] Furthermore, the local cross-correlation imaging conditions for amplitude compensation are determined by the following formula:

[0021]

[0022] Where x is the spatial location; I(x) is the image value at spatial location x; t amax is the imaging time; l is half the length of the time window; Source wave field for amplitude compensation; The detector wavefield is for amplitude compensation; the superscript * indicates compensation; s is the source; r is the detector.

[0023] Furthermore, the source normalization formula for the local cross-correlation imaging conditions is as follows:

[0024]

[0025] in, The source wave field is attenuated; the superscript · indicates attenuation.

[0026] A variable fractional-order viscous acoustic wave equation attenuation-compensated reverse-time migration system includes:

[0027] A building block is used to construct a stable attenuation-compensated viscous acoustic wave equation by introducing a regularization term;

[0028] The compensation unit is used to perform attenuation compensation and reverse time migration of viscous sound waves based on the viscous sound wave equation with local cross-correlation imaging conditions based on the Nyquist sampling theorem, according to the stable attenuation compensation.

[0029] Furthermore, the compensation unit is specifically used for:

[0030] During the source wave field extension process, a time window is designed for each grid point. The source wave field within the set time window range is stored with the time corresponding to the maximum amplitude of the source wave field at each grid point as the center.

[0031] During the back propagation of the detector's wave field, the wave field values ​​corresponding to the time window range are extracted one by one.

[0032] Reverse time migration imaging is performed based on the source wavefield and detector wavefield values ​​within the time window range.

[0033] Furthermore, the viscous acoustic wave equation for stable attenuation compensation is determined by the following formula:

[0034]

[0035] Where η is the control parameter for the phase term of the equation.

[0036] A computer-readable storage medium storing at least one computer-executable program, which, when executed by the computer, implements the variational fractional-order viscous acoustic wave equation attenuation compensation reverse time migration method as described in any of the preceding claims.

[0037] Compared with the prior art, the beneficial effects of the present invention are:

[0038] 1) By adding a regularization term, the exponential growth of high-frequency noise in seismic records during the compensation process can be effectively suppressed, ensuring the numerical stability of the attenuation compensation process.

[0039] 2) By adopting local cross-correlation imaging conditions based on the Nyquist sampling theorem, the computation and storage costs of Q-compensated reverse time migration can be effectively reduced while adapting to imaging of complex structures, thus improving the practicality of reverse time migration imaging of viscoelastic media.

[0040] Other features and advantages of the invention will be set forth in the description which follows, and will be apparent in part from the description, or may be learned by practicing the invention. The objects and other advantages of the invention may be realized and obtained by means of the methods pointed out in the description, claims and drawings. Attached Figure Description

[0041] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0042] Figure 1 A two-dimensional BP cloud model according to an embodiment of the present invention is shown;

[0043] Figure 2 A schematic diagram of a detector wavefield snapshot at t=0 during the wave equation absorption attenuation compensation process according to an embodiment of the present invention is shown.

[0044] Figure 3 The reverse time migration results of a two-dimensional BP gas cloud model according to an embodiment of the present invention are shown;

[0045] Figure 4 A schematic diagram of offset single-track comparison according to an embodiment of the present invention is shown;

[0046] Figure 5 A schematic diagram of the Marmousi model according to an embodiment of the present invention is shown;

[0047] Figure 6 A schematic diagram comparing the migration profiles obtained by reverse time migration under different imaging conditions and the extracted migration single channel using the Marmousi model according to an embodiment of the present invention is shown.

[0048] Figure 7 A schematic diagram comparing the offset profiles obtained with different search step sizes and the extracted offset single track is shown according to an embodiment of the present invention.

[0049] Figure 8A schematic diagram of different types of reverse time migration test results on a three-dimensional pushover model according to an embodiment of the present invention is shown. Detailed Implementation

[0050] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0051] This invention proposes a method for attenuation compensation and reverse-time migration of a variable fractional-order viscous acoustic wave equation, comprising the following steps:

[0052] S101: Construct a stable attenuation-compensated viscous acoustic wave equation by introducing a regularization term;

[0053] The specific expression of the second-order constant-density fractional Laplace operator viscous acoustic wave equation reported in existing literature is as follows:

[0054]

[0055]

[0056] Where c is the seismic wave velocity; c0 is the reference seismic wave velocity defined at the reference angular frequency; p is the wave field; t is time; ▽ 2 η is the Laplace operator; γ is the fractional parameter, which takes values ​​between 0 and 1 / 2 for any quality factor; η is the control parameter for the phase term of the equation; τ is the control parameter for the amplitude term of the equation; ω0 is the reference angular frequency; and Q is the quality factor.

[0057] In formula (1) p controls the velocity dispersion of seismic waves; Controlling the amplitude attenuation of seismic waves is advantageous for developing reverse-time migration methods that can simultaneously compensate for amplitude attenuation and velocity dispersion. In the decoupled viscous wave equation, amplitude compensation and phase correction are achieved by changing the sign of the amplitude term while keeping the sign of the phase term unchanged. Velocity dispersion refers to the phenomenon where the propagation speed of different frequency components of seismic waves varies due to the absorption and attenuation effect of the subsurface medium during propagation.

[0058] Since formula (1) is decoupled, the counter-time propagation of the detector wavefield during the compensation process is physically equivalent to phase correction. Therefore, velocity dispersion can be ignored when adding the regularization term. The compensation form of the viscous acoustic wave equation, considering only amplitude attenuation, can then be expressed as:

[0059]

[0060] Directly solving equation (3) leads to instability in the numerical simulation because compensation is the inverse process of attenuation, which involves exponential energy growth. Noise in the seismic data will exacerbate the exponential amplification effect over time, resulting in numerical instability. In fact, the unstable wavefield p during compensation can be represented by introducing an intermediate wavefield q and a time-varying exponential term, i.e.:

[0061] p = e αt q (4)

[0062] Where e is the natural constant and α is the growth factor.

[0063] The growth factor α is determined by the following formula:

[0064]

[0065] It should be noted that α is a positive number because τ is a real number less than 0. Substituting the wave field p with formula (4) into formula (3), we get:

[0066]

[0067] Therefore, formula (3) can be rewritten as:

[0068]

[0069] In the first equation of formula (7) Without the first-order time derivative representing decay, the equation can be stably simulated. The intermediate wave field q obtained from this can be regarded as a stable solution of equation (3), while the actual solution p grows exponentially with time and space wavenumber.

[0070] To ensure the stability of the anti-propagating wavefield during the compensation process, a constraint term εα is introduced into the exponential term of formula (7). 2 Then the form of formula (4) changes to:

[0071]

[0072] The regularization factor ε is a very small positive number. Substituting equation (8) into equation (6) and discarding higher-order terms, we obtain the viscous acoustic wave equation with no velocity dispersion and stable amplitude attenuation compensation after adding the regularization term, which is expressed as follows:

[0073]

[0074] Now, taking the velocity dispersion term into account, the final viscous acoustic wave equation for stable attenuation compensation is:

[0075]

[0076] The last term in formula (10) p is the constructed regularization term. By selecting a suitable regularization factor ε, the stability of the wavefield backpropagation process can be guaranteed, and stable compensation of Q-compensated reverse time migration can be achieved.

[0077] The effectiveness of the stability compensation method proposed in this invention, which employs a regularization term, is verified using a two-dimensional BP gas cloud model. Figure 1 As shown, the two-dimensional BP cloud model has a grid size of 161×398 and a grid spacing of 10m. Figure 1 (a) shows the velocity model, and (b) shows the quality factor model. A 25Hz Ricker wavelet was used as the vertical source, placed at grid number (z,x) (1, 200), with a reference frequency of 200Hz and a time sampling interval of 1ms, for a total recording time of 2.0s. Figure 2 As shown, using the above parameters, a forward modeling of formula (1) is performed to obtain a single-shot seismic record with absorption attenuation. This record is then used as the source backpropagation data to obtain a compensated detector wavefield snapshot at time t=0. Wherein, Figure 2 (a) shows a detector wave field snapshot with direct compensation and without regularization strategy. It can be seen that the wave field snapshot with direct compensation has shown instability to varying degrees. Figure 2 (b) shows a snapshot of the detector wavefield obtained using a regularization strategy, where the regularization factor is set to ε = 3.5 × 10⁻⁶. -7 There are no stability issues, which verifies the effectiveness of the regularization method proposed in this invention.

[0078] Further reverse-time migration was performed on the two-dimensional BP gas cloud model, with a total of 60 shots, a shot spacing of 60m, and the first shot being 150m from the left edge of the model. Figure 3 As shown, where, Figure 3 (a) is the reference profile obtained by using a fully elastic acoustic reverse time migration algorithm with attenuated seismic data. Figure 3 (b) An uncompensated migration profile obtained using a fully elastic acoustic reverse-time migration algorithm with attenuated seismic data. Figure 3 (c) Imaging results obtained by applying attenuation-compensated reverse-time migration based on a regularization strategy to attenuated seismic data. A comparison shows that the migration profile without attenuation compensation... Figure 3 (b) has a relatively weak amplitude. Due to the presence of a strongly attenuated gas reservoir, the imaging results of the lower part of the reservoir are poor, with significant damage to the migration energy and resolution, and most reflection interfaces are blurred. The attenuation-compensated profile... Figure 3 (c) can better display the amplitude information that should be there. The fault under the air layer shows more effective information and is better illuminated.

[0079] like Figure 4 As shown, to display more detailed information, in Figure 3 Three offset tracks were extracted from the horizontal positions x=1km, x=2km, and x=3km in the offset results for comparison, and the results are as follows: Figure 4 As shown in (a), (b), and (c), the solid black line represents the reference offset channel, the dashed black line represents the offset channel without attenuation compensation, and the hollow black circle corresponds to the result of attenuated reverse-time offset. By comparison, it is easy to see that the uncompensated curve has weaker energy, a wider waveform, and significantly reduced resolution compared to the reference channel, especially... Figure 4 (b) The single channel shown passes through the gas cloud region, which severely reduces the reliability of the migration imaging results. In contrast, in the compensated curve, the attenuation amplitude and the distorted phase are well recovered, the agreement with the reference migration channel is high, and the location of the subsurface reflector is correctly corrected.

[0080] By adding a regularization term, the exponential growth of high-frequency noise in seismic records during the compensation process can be effectively suppressed, ensuring the numerical stability of the compensation process and significantly improving the imaging quality of reverse time migration of seismic data containing attenuation.

[0081] S102: Based on the viscous acoustic wave equation with stable attenuation compensation, the attenuation compensation and reverse time migration of the viscous acoustic wave are performed using the local cross-correlation imaging condition based on the Nyquist sampling theorem.

[0082] Consistent with conventional reverse-time migration, the traditional cross-correlation imaging conditions also apply to Q-compensated reverse-time migration. The corresponding Q-compensation form is shown in equation (11):

[0083]

[0084] Where x is the spatial location; I(x) is the image value at spatial location x, and t max Indicates the maximum recording duration; Source wave field for amplitude compensation; The detector wave field is for amplitude compensation; the superscript * indicates compensation; s is the source; r is the detector.

[0085] Although the above-mentioned cross-correlation imaging condition formula (11) can accurately image the target and is unconditionally stable, the obtained imaging values ​​lack physical meaning. Therefore, some scholars have proposed a normalized cross-correlation imaging condition, and the corresponding Q-compensation formula is as follows:

[0086]

[0087] in, The superscript · represents the attenuated source wavefield. Applying formula (12) not only yields accurate reflection coefficients but also provides illumination compensation for deep imaging, while eliminating the influence of the source wavelet on the imaging results. Since the source wavefield and the detector wavefield are extended in opposite directions in the time direction, using formula (11) or (12) for Q-compensated reverse time migration will result in massive storage and computational demands, especially for three-dimensional Q-compensated reverse time migration. To address this, some scholars have proposed excitation amplitude imaging conditions, considering that the form of attenuation compensation can be expressed as:

[0088]

[0089] Among them, t amax The time at which the maximum amplitude value is recorded at each grid point in the positive continuation of the source wavefield is the imaging time. In fact, formula (13) is derived from formula (12) by retaining only the time equal to t. amax This method, which uses only excitation amplitude points for imaging, significantly reduces the storage burden of reverse time migration and eliminates the need for source wavefield reconstruction. However, because it does not fully utilize wavefield information, the imaging results are susceptible to the signal-to-noise ratio of the data.

[0090] Weighing the advantages and disadvantages of the two imaging conditions mentioned above, this invention proposes a local cross-correlation imaging condition, the specific steps of which are as follows:

[0091] First, during the source wavefield extension process, a time window is designed for each grid point. Centered on the time corresponding to the detection of the maximum amplitude of the source wavefield at each grid point, the corresponding source wavefield is stored within the set time window range. For example, the time window size is usually set to 1.5 times the apparent period corresponding to the dominant frequency of the seismic wavelet. Each grid point will have a time window of the same length, but the data stored in each time window is different. Local cross-correlation means that only local data within a given range is retained in the data at all times, and the range of data to be stored is defined by the time window.

[0092] Secondly, during the back propagation of the detector's wave field, the wave field values ​​corresponding to the above-mentioned time window range are extracted one by one.

[0093] Finally, based on the source wavefield and detector wavefield values ​​within the aforementioned time window, the following attenuation compensation formula based on local cross-correlation imaging conditions is applied to perform reverse time migration imaging:

[0094] The formula for the local cross-correlation imaging condition with amplitude compensation is as follows:

[0095]

[0096] Where l is half the length of the local cross-correlation calculation window, and the corresponding source normalization formula is:

[0097]

[0098] It is worth noting that in this invention, the time interval corresponding to the search for the maximum amplitude of the wavefield at each grid point during the source wavefield extension process is defined as the search step size. Generally, this search step size is the same as the time sampling interval. However, in practical data processing, such as high-precision three-dimensional reverse-time migration, a smaller spatial grid is often required. To meet stringent stability conditions, the time sampling interval for wavefield extension is often smaller. In this case, oversampling may occur in the local time window. If the source wavefield is stored directly using the simulated time step, the storage capacity of the source wavefield will increase significantly, affecting the practicality of the local cross-correlation imaging conditions. Therefore, this invention combines the Nyquist sampling theorem with the local cross-correlation imaging conditions, setting the search step size t for the maximum amplitude at each grid point during the search for the source wavefield extension process. s The formula is:

[0099]

[0100] Among them, f max f is the maximum frequency of the seismic signal; min Assuming f is the minimum frequency of the seismic signal, this invention assumes f min It is 0Hz.

[0101] The effectiveness of the local cross-correlation imaging conditions of this invention was verified using the Marmosui model. For simplicity, all tests on the local cross-correlation algorithm were conducted in acoustic media, and the conclusions obtained are also applicable to viscous elastic media. It should be noted that all velocity models used for reverse-time migration imaging were obtained from the true velocity models using Gaussian smoothing with a smoothing factor of 10. Figure 5 As shown, the Marmousi model has a grid of 663×234 points with a grid spacing of 10m, and the model velocity ranges from 2000m / s to 4500m / s. A Ricker wavelet with a dominant frequency of 25Hz was used as the vertical source. The first shot was 100m from the left edge of the model, and a total of 60 shots were fired with a shot spacing of 100m. The seismic signal was received by 663 geophones uniformly distributed on the surface, with a time sampling interval of 1ms and a recording duration of 2.5s. The time window length was set to 101 search steps, with each search step being 1ms.

[0102] like Figure 6 As shown, Figure 6 (a) is the offset profile obtained using global cross-correlation imaging conditions, which serves as the reference profile; Figure 6 (b) shows the offset profile using local cross-correlation imaging conditions. As can be seen, the results of the two imaging conditions are not significantly different. Figure 6 (c) Shows the normalized seismic traces located at a horizontal position of 3600m. Figure 6 (d) Shows the normalized seismic traces located at a horizontal position of 5300m. To show more detailed information, from... Figure 6 Seismic traces at horizontal locations of 3600m and 5300m were extracted from (a) and (b) and displayed respectively. Figure 6 (c) and Figure 6 In (d), after comparison, even for the complex Marmousi model, the shift results of the local cross-correlation imaging conditions proposed in this invention can be compared with those of the Marmousi model. Figure 6 The good seismic trace matching in (a) proves the effectiveness of the local cross-correlation imaging conditions.

[0103] Next, the Marmosui model is used again to verify the practicality of the local cross-correlation imaging conditions based on the Nyquist sampling theorem in this invention. Given the dominant frequency, the empirical formula f... max =2.5f d The maximum frequency f of the seismic record signal can be calculated. d The frequency of the seismic wavelet is given. According to formula (16), the maximum search step size is 8 ms.

[0104] like Figure 7 As shown, reverse time migration imaging was performed using local cross-correlation imaging conditions. In the tests, a search step size of 101 time windows and a search step size of 1 ms were used. The imaging results are as follows. Figure 7 As shown in (a); a search window of 25 steps was used, with each search step lasting 4 ms. The imaging results are as follows. Figure 7 (b); Using a time window of 13 search steps, with each search step being 8ms, the imaging results are as follows. Figure 7 As shown in (c). It is difficult to discern the differences between the above images directly from the migration results. To more clearly demonstrate the details, normalized seismic traces at a horizontal position of 4500m were extracted from the migration results and compared. Figure 7(d) shows that, under the Nyquist sampling theorem, the imaging accuracy is almost identical for different search step sizes. Table 1 compares the storage overhead and runtime of three imaging conditions during the reverse time migration of the 60-shot Marmousi model. In Table 1, "local" represents the local cross-correlation imaging condition, "local + sampling" represents the local cross-correlation imaging condition based on the Nyquist sampling theorem, and "global" represents the traditional cross-correlation imaging condition. It can be seen that as the search step size increases, the storage overhead and computation time decrease significantly. When the search step size is 8 ms, the storage overhead of the local cross-correlation imaging condition based on the Nyquist sampling theorem is only 12.5% ​​of that of the local cross-correlation imaging condition and only 0.52% of that of the global cross-correlation imaging condition. Moreover, its computation time is only 55.2% of that of the local cross-correlation imaging condition and 17.86% of that of the global cross-correlation imaging condition. This demonstrates the practicality of the local cross-correlation imaging condition based on the Nyquist sampling theorem proposed in this invention, which significantly reduces the storage and computation costs of reverse time migration imaging while ensuring imaging accuracy.

[0105] Table 1 Comparison of storage overhead and runtime under different imaging conditions

[0106]

[0107] The calculation results mentioned in Table 1 above were all tested on a PC with a main frequency of 3.2GHz and a graphics card model of RTX 1080.

[0108] Finally, as Figure 8 As shown, the feasibility of applying the regularization strategy and the local cross-correlation imaging condition based on the Nyquist sampling theorem proposed in this invention to Q-compensated reverse time migration is verified using a three-dimensional pushover model. The time window length of the selected local cross-correlation imaging condition based on the Nyquist sampling theorem is 13 search steps, and the search step size is 8ms. Figure 8 (a) is a three-dimensional velocity model of the thrust body, and the quality factor is obtained from an empirical formula:

[0109]

[0110] The model has a total grid size of 100×250×250, a grid spacing of 20m, a time sampling interval of 1ms, and a recording duration of 1.5s. The seismic source is a Ricker wavelet with a dominant frequency of 25Hz and a reference frequency of 200Hz. The regularization parameter ε = 3.5×10⁻⁶ is selected during the compensation process. -7 The observation system was set up as follows: 20 survey lines were evenly distributed on the model surface, with 20 shots on each line, for a total of 400 shots. The first shot of each survey line was 50m from the left edge of the model, and a total of 62,500 geophones were evenly distributed on the surface. Figure 8(b) shows the results of fully elastic acoustic reverse time migration using unattenuated seismic data as a reference profile for absorption attenuation compensation. Figure 8 (c) The results of fully elastic acoustic reverse time migration of attenuated seismic data show that the amplitude and phase of the absorbed attenuated seismic data are distorted. Conventional reverse time migration results in weak imaging energy, misalignment of the phase axis, and lower resolution than the reference profile. Figure 8 (d) shows the results of Q-compensated reverse time migration of attenuated seismic data. It can be clearly seen that the compensated profile matches the reference profile in both energy and phase, which confirms the feasibility of the stable and efficient viscous acoustic equation absorption attenuation compensation reverse time migration method proposed in this invention.

[0111] By employing local cross-correlation imaging conditions based on the Nyquist sampling theorem, the computational and storage costs of Q-compensated reverse time migration can be effectively reduced while adapting to imaging of complex structures, thus improving the practicality of reverse time migration imaging of viscoelastic media.

[0112] Based on the above-mentioned method for attenuation compensation and reverse time migration using the variable fractional-order viscous acoustic wave equation, this invention proposes a variable fractional-order viscous acoustic wave equation attenuation compensation and reverse time migration system, comprising:

[0113] A building block is used to construct a stable attenuation-compensated viscous acoustic wave equation by introducing a regularization term;

[0114] The compensation unit is used to perform attenuation compensation and reverse time migration of viscous sound waves based on the viscous sound wave equation with local cross-correlation imaging conditions based on the Nyquist sampling theorem, according to the stable attenuation compensation.

[0115] The compensation unit is specifically used for:

[0116] During the source wave field extension process, a time window is designed for each grid point. The source wave field within the set time window range is stored with the time corresponding to the maximum amplitude of the source wave field at each grid point as the center.

[0117] During the back propagation of the detector's wave field, the wave field values ​​corresponding to the time window range are extracted one by one.

[0118] Reverse time migration imaging is performed based on the source wavefield and detector wavefield values ​​within the time window range.

[0119] The viscous acoustic wave equation for stable attenuation compensation is determined by the following formula:

[0120]

[0121] Where η is the control parameter for the phase term of the equation.

[0122] The present invention provides a computer-readable storage medium storing at least one computer-executable program, which, when executed by the computer, implements the variational fractional-order viscous acoustic wave equation attenuation compensation reverse time migration method as described above.

[0123] Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for attenuation compensation and reverse-time migration of a variable fractional-order viscous acoustic wave equation, characterized in that, Includes the following steps: A stable attenuation-compensated viscous acoustic wave equation is constructed by introducing a regularization term; Based on the viscous acoustic wave equation with stable attenuation compensation, the attenuation compensation and reverse time migration of viscous acoustic waves are performed using the local cross-correlation imaging condition based on the Nyquist sampling theorem. The regularization term is determined by the following formula: in, As a regularization factor; This refers to the velocity of seismic waves; These are the control parameters for the amplitude term of the equation; For wave field; For time; For the Laplace operator; For fractional-order parameters.

2. The method for attenuation compensation and reverse-time migration of the variational fractional-order viscous acoustic wave equation according to claim 1, characterized in that, The viscous acoustic wave equation for stable attenuation compensation is determined by the following formula: in, These are the control parameters for the phase term of the equation.

3. The method for attenuation compensation and reverse-time migration of the variational fractional-order viscous acoustic wave equation according to claim 1, characterized in that, The specific steps for establishing the local cross-correlation imaging conditions are as follows: During the source wave field extension process, a time window is designed for each grid point. The source wave field within the set time window range is stored with the time corresponding to the maximum amplitude of the source wave field at each grid point as the center. During the back propagation of the detector's wave field, the wave field values ​​corresponding to the time window range are extracted one by one. Reverse time migration imaging is performed based on the source wavefield and detector wavefield values ​​within the time window range.

4. The method for attenuation compensation and reverse-time migration of the variational fractional-order viscous acoustic wave equation according to claim 1 or 3, characterized in that, The amplitude-compensated local cross-correlation imaging conditions are determined by the following formula: in, Spatial location; For spatial location The image value at that location; The imaging time; l It is half the length of the time window; Source wave field for amplitude compensation; The detector wave field is for amplitude compensation; the superscript * indicates compensation. s The epicenter; r For detectors.

5. The method for attenuation compensation and reverse-time migration of the variational fractional-order viscous acoustic wave equation according to claim 4, characterized in that, The source normalization formula for the local cross-correlation imaging conditions is as follows: in, For the attenuated source wave field, superscript This indicates attenuation.

6. A variable fractional-order viscous acoustic wave equation attenuation-compensated reverse-time migration system, characterized in that, include: A building block is used to construct a stable attenuation-compensated viscous acoustic wave equation by introducing a regularization term; The compensation unit is used to perform attenuation compensation and reverse time migration of viscous sound waves based on the viscous sound wave equation with local cross-correlation imaging conditions based on the Nyquist sampling theorem, based on the stable attenuation compensation. The regularization term is determined by the following formula: in, As a regularization factor; This refers to the velocity of seismic waves; These are the control parameters for the amplitude term of the equation; For wave field; For time; For the Laplace operator; For fractional-order parameters.

7. The variational fractional-order viscous acoustic wave equation attenuation compensation reverse-time migration system according to claim 6, characterized in that, The compensation unit is specifically used for: During the source wave field extension process, a time window is designed for each grid point. The source wave field within the set time window range is stored with the time corresponding to the maximum amplitude of the source wave field at each grid point as the center. During the back propagation of the detector's wave field, the wave field values ​​corresponding to the time window range are extracted one by one. Reverse time migration imaging is performed based on the source wavefield and detector wavefield values ​​within the time window range.

8. The variational fractional-order viscous acoustic wave equation attenuation compensation reverse-time migration system according to claim 6, characterized in that, The viscous acoustic wave equation for stable attenuation compensation is determined by the following formula: in, These are the control parameters for the phase term of the equation.

9. A computer-readable storage medium storing at least one computer-executable program, characterized in that, When the at least one program is executed by the computer, it implements the variable fractional-order viscous acoustic wave equation attenuation compensation reverse time migration method as described in any one of claims 1-5.

Images