High-precision fractional derivative seismic reverse time migration method and system

Abstract

The invention provides a high-precision fractional derivative seismic reverse time migration method and system, and belongs to the technical field of seismic imaging in oil-gas exploration and development. The method comprises the following steps of: 1, inputting seismic wavelets, and reversely transmitting seismic data; 2, obtaining a seismic source wave field by using the seismic wavelets; 3, obtaining a receiving wave field by using the back propagation seismic data; and 4, according to the seismic source wave field and the receiving wave field, obtaining a final migration image by using an imaging condition. According to the method, the operation efficiency and the seismic simulation precision are greatly improved, the imaging efficiency is improved, and the accuracy and the resolution of an imaging result are improved. Reverse time migration imaging is carried out on the attenuation medium complex structure, the method is suitable for high-precision imaging of complex geological conditions and attenuation salient regions, and powerful support is provided for follow-up migration velocity analysis and oil reservoir description work.

Description

Technical Field

[0001] The present invention belongs to the technical field of seismic imaging in oil and gas exploration and development, and in particular relates to a high-precision fractional-order derivative seismic reverse time migration method and system. Background Art

[0002] Traditional reverse time migration uses the cross-correlation of the source wave field obtained by forward propagation of seismic wavelets and the receiving wave field obtained by reverse propagation of seismic data to obtain the migration image. However, in viscoacoustic attenuation media, seismic data has obvious dispersion and attenuation effects. The results obtained by the traditional reverse time migration method have low signal-to-noise ratio, inaccurate layer positions, and poor quality. Therefore, compensation for the attenuation effect should be considered in reverse time migration to improve the resolution and reliability of seismic imaging.

[0003] The attenuation characteristics of the earth medium can be combined with the physical models of springs and viscoelastic pots, such as the common standard linear body (SLS) model, Kevin-Vogt model, etc., which are characterized by relatively simple calculations and can be solved on a large scale using finite differences. However, this type of model cannot accurately describe the constant Q characteristics of the earth's attenuation and does not match the measured results. Another type of attenuation model starts from the constant Q characteristics and derives a model based on fractional spatial derivatives, which can accurately describe the attenuation characteristics that are independent of earthquake frequency, but its phase velocity is inaccurate, especially in the low-frequency part. The generalized viscoacoustic wave equation combines the above two types of models and can accurately simulate the earthquake attenuation characteristics.

[0004] The generalized viscoacoustic wave equation also expresses the viscosity characteristics in the form of fractional time derivatives. However, the numerical calculation of fractional time derivatives is very difficult and requires a lot of storage space, which is not suitable for large-scale seismic simulation and imaging.

[0005] Therefore, an efficient solution is to convert the fractional time derivative into a spatial derivative for calculation. The fractional spatial derivative is efficiently calculated using the pseudospectral method. However, the pseudospectral method is not suitable for inhomogeneous media in theory. The usual practice is to use the average attenuation parameter instead of the global space-variant attenuation parameter for approximate calculation. This will lead to large simulation errors, especially in complex geological conditions such as oil and gas-rich areas. Some improvement schemes include increasing the number of fractional derivatives or using the integral method to improve accuracy, but this greatly increases the amount of calculation and reduces the efficiency of simulation and migration.

[0006] For attenuation compensation, the widely used method is to use inverse Q filtering to process seismic data, but it is based on the one-dimensional model assumption and cannot handle complex geological conditions. In reverse time migration, attenuation can be compensated during the forward and reverse propagation process. Conventional seismic reverse time migration cannot decouple dispersion and attenuation, and cannot correct the phase when compensating for amplitude, thus limiting the effect of attenuation compensation and the improvement of imaging resolution. Summary of the invention

[0007] The purpose of the present invention is to solve the difficulties existing in the above-mentioned prior art, and to provide a high-precision fractional-order derivative seismic reverse time migration method and system based on the generalized viscoacoustic wave equation, so as to solve the problems of low signal-to-noise ratio, phase difference and inaccurate layer position in seismic migration in attenuating media, significantly improve the efficiency and accuracy of viscoacoustic medium reverse time migration based on fractional-order derivatives, and improve the resolution and accuracy of the migration results.

[0008] The present invention is achieved through the following technical solutions:

[0009] A first aspect of the present invention provides a high-precision fractional-order derivative seismic reverse time migration method, the method comprising:

[0010] Step 1: Input seismic wavelet and back-transmitted seismic data;

[0011] Step 2: Use seismic wavelets to obtain the source wave field;

[0012] Step 3: Use the back-propagation seismic data to obtain the receiving wave field;

[0013] Step 4: Based on the source wavefield and the receiving wavefield, the imaging conditions are used to obtain the final migration image.

[0014] Preferably, the first step includes:

[0015] The collected seismic data are processed to obtain seismic wavelets and back-propagation seismic data.

[0016] Preferably, the second step includes:

[0017] The source wave field is obtained using the dispersion equation:

[0018]

[0019]

[0020] Among them, u S (x) represents the source wave field, c represents the viscoelastic velocity, t represents the time, ▽ 2 represents the Laplace operator, β represents the viscoelastic parameter related to the quality factor, ω0 is the reference frequency, c0 is the reference velocity, f(β(x)) represents the correction function, w(t) represents the seismic wavelet, x represents the spatial position, x0 represents the source position, and δ(x-x0) represents the Dirac-Delta function.

[0021] Preferably, the operation of the third step includes:

[0022] The received wave field is obtained using the attenuation compensation equation:

[0023]

[0024]

[0025]

[0026] Among them, u R (x) represents the received wave field, c represents the viscoelastic velocity, t represents the time, ▽ 2 represents the Laplace operator, β represents the viscoelastic parameter related to the quality factor, ω0 is the reference frequency, c0 is the reference velocity, f(β(x)) represents the correction function, and d(t,x) represents the back-propagation seismic data.

[0027] Preferably, the fourth step includes:

[0028] Imaging is performed according to the source cross-correlation imaging condition, and the final migration image I(x) is obtained using the following formula:

[0029]

[0030] where u S (x) and u R (x) represent the source wavefield and the receiving wavefield respectively, and t represents time.

[0031] A second aspect of the present invention provides a high-precision fractional-order derivative seismic reverse time migration system, the system comprising:

[0032] An input unit, used for inputting seismic wavelets and back-propagation seismic data;

[0033] A source wave field acquisition unit, connected to the input unit, for obtaining a source wave field using seismic wavelets;

[0034] A receiving wave field acquisition unit, connected to the input unit, for obtaining a receiving wave field using back-transmitted seismic data;

[0035] The migration imaging unit is connected to the source wave field acquisition unit and the receiving wave field acquisition unit respectively, and is used to obtain the final migration image according to the source wave field and the receiving wave field by using imaging conditions.

[0036] Preferably, the source wave field acquisition unit performs the following operations:

[0037] The source wave field is obtained using the dispersion equation:

[0038]

[0039]

[0040] Among them, u S(x) represents the source wave field, c represents the viscoelastic velocity, t represents the time, ▽ 2 represents the Laplace operator, β represents the viscoelastic parameter related to the quality factor, ω0 is the reference frequency, c0 is the reference velocity, f(β(x)) represents the correction function, w(t) represents the seismic wavelet, x represents the spatial position, x0 represents the source position, and δ(x-x0) represents the Dirac-Delta function.

[0041] Preferably, the receiving wave field acquisition unit performs the following operations:

[0042] The received wave field is obtained using the attenuation compensation equation:

[0043]

[0044]

[0045]

[0046] Among them, u R (x) represents the received wave field, c represents the viscoelastic velocity, t represents the time, ▽ 2 represents the Laplace operator, β represents the viscoelastic parameter related to the quality factor, ω0 is the reference frequency, c0 is the reference velocity, f(β(x)) represents the correction function, and d(t,x) represents the back-propagation seismic data.

[0047] Preferably, the offset imaging unit performs the following operations:

[0048] According to the source cross-correlation imaging condition, imaging is performed and the final migration image I(x) is obtained using the following formula:

[0049]

[0050] where u S (x) and u R (x) represent the source wavefield and the receiving wavefield respectively, and t represents time.

[0051] According to a third aspect of the present invention, a computer-readable storage medium is provided, wherein the computer-readable storage medium stores at least one computer-executable program, and when the at least one program is executed by the computer, the computer executes the steps of a high-precision fractional-order derivative seismic reverse time migration method of the present invention.

[0052] Compared with the prior art, the present invention has the following beneficial effects:

[0053] The present invention is based on the generalized viscoacoustic wave equation, converts the fractional time derivative into the fractional space derivative, greatly improves its computational efficiency, and proposes a frequency-independent correction function, thereby improving the accuracy of earthquake simulation.

[0054] The present invention uses the high-precision viscoacoustic equation for reverse time migration, adopts different equations to generate source wavefield and receiving wavefield, performs phase correction and energy compensation on the attenuation effect, thereby improving the efficiency of imaging and the accuracy and resolution of imaging results.

[0055] The present invention performs reverse time migration imaging of the complex structure of the attenuation medium, which is suitable for high-precision imaging of complex geological conditions and areas with significant attenuation, and provides strong support for subsequent migration velocity analysis and reservoir description work. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] Figure 1 A flowchart of the steps of the method of the present invention;

[0057] Figure 2 Comparison of analytical solutions of formula (1) and formula (5);

[0058] Figure 3a Results of conventional reverse time migration using unattenuated data (reference);

[0059] Figure 3b Traditional reverse time migration results using attenuation data;

[0060] Figure 3c Use the attenuation data to perform attenuation compensation on the reverse time migration results;

[0061] Figure 4 Image details of certain locations. DETAILED DESCRIPTION

[0062] The present invention is further described in detail below in conjunction with the accompanying drawings:

[0063] The invention provides a high-precision fractional-order derivative seismic reverse time migration method, which is mainly used in seismic migration in areas with significant attenuation to improve the resolution and reliability of reverse time migration imaging results.

[0064] The present invention transforms the generalized viscoacoustic wave equation into a fractional-order spatial derivative form, uses the pseudo-spectral method to greatly improve the calculation efficiency, and adds a frequency-independent space-varying correction function to correct the numerical error of the pseudo-spectral method to improve the simulation accuracy.

[0065] The equation to which the present invention is applied is as follows:

[0066] 1. Construct the fractional-order spatial derivative viscoacoustic wave equation:

[0067] The fractional time derivative viscoacoustic wave equation (FTD) can be written as:

[0068]

[0069] Where u is the scalar wave field, ω0 is the reference frequency, β is the viscoelastic parameter related to the quality factor, when β = 0, the above equation becomes a pure elastic wave equation, c(β) is the viscoelastic velocity, t represents time, ▽ 2 represents the Laplace operator.

[0070] Formula (1) is in frequency wave number, that is

[0071]

[0072] Where ω is the angular frequency, k is the circular wave number, is the complex wave field in the frequency wavenumber domain, c is the abbreviation of c(β). Considering that the attenuation characteristics of the earth medium are usually weak, ω / c0≈k is approximately valid, c0 is the reference velocity, and formula (2) can be approximated as:

[0073]

[0074] in,

[0075]

[0076] By using the inverse Fourier transform back to the time-space domain, we can get the fractional-order spatial derivative wave equation (FSD):

[0077]

[0078] By comparing the analytical solutions of formula (5) and formula (1), it can be obtained that the approximated formula (5) is accurate and its RMS error is very small.

[0079] 2. Construct the correction function. The above formula (5) can be quickly solved by pseudospectral method, but it can only be applied to homogeneous media. In non-homogeneous media, the average value of β is usually used. To replace the spatially varying β(x) (β(x) represents the spatially varying β). In order to improve the accuracy, a correction function f(β(x)) is introduced to satisfy:

[0080]

[0081] Considering the weak attenuation assumption, the correction function can be written as

[0082]

[0083] Among them, ω m is the center frequency of the frequency band. Based on the relatively narrow seismic frequency band and the limited spatial variation of β(x), the center frequency assumption (the last line of Formula 7) holds. This correction function is independent of the frequency variable. Substituting it into Formula (5) and converting it into the time-space domain equation, the corrected equation is:

[0084]

[0085] 3. Use the finite difference method and pseudospectral method to solve the equation (both formula (5) and formula (8) can be solved). The integer-order time and space derivatives can be directly calculated using the finite difference method, and the fractional-order derivatives can be calculated using the pseudospectral method, which is expressed as:

[0086]

[0087] in, and is the Fourier forward and inverse transform, and x is the spatial variable. Formula (9) is the solution formula for calculating the fractional derivatives in formulas (5) and (8).

[0088] 4. Establish the equation for attenuation compensation of reverse time migration attenuation. Formula (8) can be regarded as a decoupled viscoacoustic wave equation, in which the C1 term represents dispersion and the C2 term represents energy attenuation. Since seismic data contains both dispersion and energy attenuation, in the reverse time migration process, the same phase velocity is used for transmission and the same energy compensation is used to restore the seismic signal. Therefore, in the attenuation compensation equation, the dispersion term remains unchanged and the attenuation term is reversed. The attenuation compensation equation is as follows:

[0089]

[0090] Since the energy of high-frequency noise signals will also expand exponentially during the energy expansion process, resulting in numerical instability, a low-pass filter is used to ensure the stability of the algorithm during the compensation process.

[0091] 5: Calculate the source wave field and the receiving wave field. For the source wave field, only the dispersion is considered, that is, the energy absorption term C2 = 0, then the dispersion equation can be written as

[0092]

[0093] The dispersion equation only contains dispersion, which can ensure phase accuracy, but the equation does not include energy absorption. Compared with the original forward attenuation wave field, the energy has been exponentially expanded. This equation has only one fractional derivative term, which can improve the efficiency of numerical simulation. For the receiving wave field, both dispersion and energy must be considered, so the complete attenuation compensation equation (10) is used. u represents the calculated wave field.

[0094] The numerical solutions of formula (10) and formula (11) are realized by the method in step 3. Formula (5) and formula (8) in steps 1 and 2 can be used as the theoretical basis for subsequent development and can also be used to simulate the propagation of seismic waves in attenuating media.

[0095] Based on the decoupling characteristics of the above-mentioned fractional-order spatial derivative equation, the present invention uses the dispersion equation to obtain the source wave field, uses the attenuation compensation equation to obtain the receiving wave field, and finally uses the imaging condition to obtain the final migration image, thereby improving the efficiency of migration imaging and ensuring the imaging resolution.

[0096] Specifically, the method of the present invention comprises:

[0097] Step 1: Input seismic wavelet and back-transmitted seismic data;

[0098] Step 2: Use seismic wavelets to obtain the source wave field;

[0099] Step 3: Use the back-propagation seismic data to obtain the receiving wave field;

[0100] Step 4: Based on the source wavefield and the receiving wavefield, the imaging conditions are used to obtain the final migration image.

[0101] The embodiments of the method of the present invention are as follows:

[0102] Embodiment 1:

[0103] The first step includes:

[0104] The collected seismic data are processed using existing processing methods to obtain seismic wavelets and back-propagation seismic data.

[0105] Embodiment 2:

[0106] The second step includes:

[0107] The source wave field is obtained using the dispersion equation:

[0108]

[0109]

[0110] Among them, u S (x) represents the source wave field, c represents the viscoelastic velocity, t represents the time, ▽ 2 represents the Laplace operator, β represents the viscoelastic parameter related to the quality factor, ω0 is the reference frequency, c0 is the reference velocity, f(β(x)) represents the correction function, w(t) represents the seismic wavelet, x represents the spatial position, x0 represents the source position, and δ(x-x0) represents the Dirac-Delta function.

[0111] The parameters used in the above equation are mainly: seismic wavelet, velocity model and Q model. The velocity model (reference velocity) is regarded as the phase velocity of the reference frequency, and the viscoelastic velocity is the 0Hz phase velocity, which can be obtained according to the phase velocity formula; β can be obtained by the Q model, and the specific formula is:

[0112]

[0113] b i ={0.79788,0.31831,0.16787,-0.08260,-0.08730,-0.03774}.

[0114]

[0115] Embodiment three:

[0116] The third step includes:

[0117] The received wave field is obtained using the attenuation compensation equation:

[0118]

[0119]

[0120] Among them, u R (x) represents the received wave field, c represents the viscoelastic velocity, t represents the time, ▽ 2 represents the Laplace operator, β represents the viscoelastic parameter related to the quality factor, ω0 is the reference frequency, c0 is the reference velocity, f(β(x)) represents the correction function, and d(t,x) represents the back-propagation seismic data.

[0121] Embodiment 4:

[0122] The fourth step includes:

[0123] Imaging is performed according to the source cross-correlation imaging conditions, and the imaging cross-correlation conditions are as follows:

[0124]

[0125] where u S (x) and u R (x) represent the source wavefield and the receiving wavefield respectively, and t represents time.

[0126] According to the source wavefield and the receiving wavefield, the final migration image I(x) is obtained using formula (12).

[0127] Embodiment five:

[0128] A typical Marmousi model was selected for migration imaging. The model area is 7.5 km × 3 km, the spatial spacing is 10 m, and the seismic wavelet is a 20 Hz Ricker wavelet. Figure 3a to Figure 3c Results of the migration imaging are shown. Figure 3a This is the result of traditional reverse time migration using unattenuated data, which can be used as a reference result. Figure 3b and3c The result of reverse time migration using attenuation data, before compensation (such as Figure 3b As shown in the figure, the amplitude of the deep migration image is very small and the imaging quality of the underground structure is very poor. However, after attenuation compensation, the amplitude is restored and the overall resolution of the image is greatly improved, which is comparable to the reference result.

[0129] Figure 4 for Figure 3a , Figure 3b , Figure 3c The gray line (the curve on the left) is the reference line. Figure 4 It can be seen that when using attenuation data, the traditional reverse time migration (green line, i.e. the line on the left that is approximately a straight line at each distance) not only has a small amplitude, but also has obvious phase distortion and a large error in the imaging position; however, after amplitude compensation (red line, i.e. the curve on the right at each distance), the image not only has consistent energy with the reference, but also has its phase restored.

[0130] The present invention also provides a high-precision fractional-order derivative seismic reverse time migration system, the system comprising:

[0131] An input unit, used for inputting seismic wavelets and back-propagation seismic data;

[0132] A source wave field acquisition unit, connected to the input unit, for obtaining a source wave field using seismic wavelets;

[0133] A receiving wave field acquisition unit, connected to the input unit, for obtaining a receiving wave field using back-transmitted seismic data;

[0134] The migration imaging unit is connected to the source wave field acquisition unit and the receiving wave field acquisition unit respectively, and is used to obtain the final migration image according to the source wave field and the receiving wave field by using imaging conditions.

[0135] Preferably, the source wave field acquisition unit performs the following operations:

[0136] The source wave field is obtained using the dispersion equation:

[0137]

[0138]

[0139] Among them, u S (x) represents the source wave field, c represents the viscoelastic velocity, t represents the time, ▽ 2 represents the Laplace operator, β represents the viscoelastic parameter related to the quality factor, ω0 is the reference frequency, c0 is the reference velocity, f(β(x)) represents the correction function, and w(t) represents the seismic wavelet.

[0140] Preferably, the receiving wave field acquisition unit performs the following operations:

[0141] The received wave field is obtained using the attenuation compensation equation:

[0142]

[0143]

[0144]

[0145] Among them, u R (x) represents the received wave field, c represents the viscoelastic velocity, t represents the time, ▽ 2 represents the Laplace operator, β represents the viscoelastic parameter related to the quality factor, ω0 is the reference frequency, c0 is the reference velocity, f(β(x)) represents the correction function, and d(t,x) represents the back-propagation seismic data.

[0146] Preferably, the offset imaging unit performs the following operations:

[0147] According to the source cross-correlation imaging condition, imaging is performed and the final migration image I(x) is obtained using the following formula:

[0148]

[0149] where u S (x) and u R (x) represent the source wavefield and the receiving wavefield respectively, and t represents time.

[0150] The present invention also provides a computer-readable storage medium, which stores at least one computer-executable program. When the at least one program is executed by the computer, the computer executes the steps in a high-precision fractional-order derivative seismic reverse time migration method of the present invention.

[0151] The present invention belongs to the technical field of seismic data processing in oil and gas exploration and development, and is mainly used to improve the quality of reverse time migration imaging of seismic data in attenuating media, and can provide strong support for subsequent migration velocity analysis and reservoir description work.

[0152] The present invention has a broad application prospect in the exploration and development of oil and gas fields. First, in terms of technology, the viscoelastic medium seismic wave propagation theory and attenuation compensation method based on wave theory are still being improved and developed, and the present invention can provide theoretical support for related algorithm research; secondly, in terms of application, the present invention can be applied in actual production, and has a relatively good application prospect in complex exploration areas such as Northwest and North China of Sinopec, and can provide effective technical support for the migration velocity analysis and reservoir description work in relevant exploration areas of Sinopec.

[0153] Based on the generalized viscoacoustic wave equation, the present invention converts the fractional time derivative into the fractional space derivative, greatly improving its computational efficiency, and proposes a frequency-independent correction function to improve the accuracy of seismic simulation. The high-precision viscoacoustic equation is used for reverse time migration, and different equations are used to generate the source wave field and the receiving wave field, and the attenuation effect is phase corrected and energy compensated, thereby improving the efficiency of imaging and the accuracy and resolution of the imaging results. The present invention performs reverse time migration imaging on the complex structure of the attenuation medium, which is suitable for high-precision imaging of complex geological conditions and areas with significant attenuation, and provides strong support for subsequent migration velocity analysis and reservoir description.

[0154] The above technical solution is only one implementation mode of the present invention. For those skilled in the art, it is easy to make various types of improvements or modifications based on the principles disclosed in the present invention, and it is not limited to the technical solution described in the above specific embodiments of the present invention. Therefore, the above description is only preferred and does not have a restrictive meaning.

Claims

1. A high-precision fractional derivative seismic reverse time migration method, characterized by: The method comprises: Step 1: Input seismic wavelet and back-transmitted seismic data; Step 2: Use seismic wavelets to obtain the source wave field; Step 3: Use the back-propagation seismic data to obtain the receiving wave field; Step 4: Based on the source wavefield and the receiving wavefield, the imaging conditions are used to obtain the final migration image.

2. The high-precision fractional-order derivative seismic reverse time migration method according to claim 1, characterized in that: The first step includes: The collected seismic data are processed to obtain seismic wavelets and back-propagation seismic data.

3. The high-precision fractional-order derivative seismic reverse time migration method according to claim 1, characterized in that: The second step includes: The source wave field is obtained using the dispersion equation: Among them, u S (x) represents the source wave field, c represents the viscoelastic velocity, t represents the time, ▽ 2 represents the Laplace operator, β represents the viscoelastic parameter related to the quality factor, ω0 is the reference frequency, c0 is the reference velocity, f(β(x)) represents the correction function, w(t) represents the seismic wavelet, x represents the spatial position, x0 represents the source position, and δ(x-x0) represents the Dirac-Delta function.

4. The high-precision fractional-order derivative seismic reverse time migration method according to claim 1, characterized in that: The third step includes: The received wave field is obtained using the attenuation compensation equation: Among them, u R (x) represents the received wave field, c represents the viscoelastic velocity, t represents the time, ▽ 2 represents the Laplace operator, β represents the viscoelastic parameter related to the quality factor, ω0 is the reference frequency, c0 is the reference velocity, f(β(x)) represents the correction function, and d(t,x) represents the back-propagation seismic data.

5. The high-precision fractional-order derivative seismic reverse time migration method according to claim 1, characterized in that: The fourth step includes: According to the source cross-correlation imaging condition, imaging is performed and the final migration image I(x) is obtained using the following formula: where u S (x) and u R (x) represent the source wavefield and the receiving wavefield respectively, and t represents time.

6. A high-precision fractional-order derivative seismic reverse time migration system, characterized by: The system comprises: An input unit, used for inputting seismic wavelets and back-propagation seismic data; A source wave field acquisition unit, connected to the input unit, for obtaining a source wave field using seismic wavelets; A receiving wave field acquisition unit, connected to the input unit, for obtaining a receiving wave field using back-transmitted seismic data; The migration imaging unit is connected to the source wave field acquisition unit and the receiving wave field acquisition unit respectively, and is used to obtain the final migration image according to the source wave field and the receiving wave field by using imaging conditions.

7. The high-precision fractional-order derivative seismic reverse time migration system according to claim 6, characterized in that: The source wave field acquisition unit performs the following operations: The source wave field is obtained using the dispersion equation: Among them, u S (x) represents the source wave field, c represents the viscoelastic velocity, t represents the time, ▽ 2 represents the Laplace operator, β represents the viscoelastic parameter related to the quality factor, ω0 is the reference frequency, c0 is the reference velocity, f(β(x)) represents the correction function, w(t) represents the seismic wavelet, x represents the spatial position, x0 represents the source position, and δ(x-x0) represents the Dirac-Delta function.

8. The high-precision fractional-order derivative seismic reverse time migration system according to claim 6, characterized in that: The receiving wave field acquisition unit performs the following operations: The received wave field is obtained using the attenuation compensation equation: Among them, u R (x) represents the received wave field, c represents the viscoelastic velocity, t represents the time, ▽ 2 represents the Laplace operator, β represents the viscoelastic parameter related to the quality factor, ω0 is the reference frequency, c0 is the reference velocity, f(β(x)) represents the correction function, and d(t,x) represents the back-propagation seismic data.

9. The high-precision fractional-order derivative seismic reverse time migration system according to claim 6, characterized in that: The offset imaging unit performs the following operations: Imaging is performed according to the source cross-correlation imaging condition, and the final migration image I(x) is obtained using the following formula: where u S (x) and u R (x) represent the source wavefield and the receiving wavefield respectively, and t represents time.

10. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores at least one computer-executable program, and when the at least one program is executed by the computer, the computer executes the steps in the high-precision fractional-order derivative seismic reverse time migration method according to any one of claims 1 to 5.

Images