Self-adaptive attenuation compensation imaging method for bottoming region

Abstract

The invention relates to a self-adaptive attenuation compensation imaging method for a bottoming region. The method comprises the following steps: performing numerical simulation calculation along a time forward direction of simulation duration to obtain a forward wave field, storing part of the forward wave field, obtaining a reconstructed seismic source wave field by using the part of the forward wave field, performing numerical simulation calculation along a time reverse direction of the simulation duration to obtain a reverse wave field, and performing correlation imaging by using the reconstructed seismic source wave field and the reverse wave field to obtain a seismic source wave field. According to the method, numerical simulation can be carried out, meanwhile, the situation that a high-frequency numerical value is unstable and the problem of imaging phase dislocation can be effectively solved, the imaging quality of a bottoming high-attenuation area is improved, and the imaging requirements of subsequent reservoir prediction and hydrocarbon detection are met.

Description

Technical Field

[0001] The present invention relates to the technical field of oil and gas exploration and development, and more specifically, to an adaptive attenuation compensation imaging method for a diapir area. Background Art

[0002] The Yinggehai Basin has a high hydrocarbon generation intensity, good vertical migration and accumulation conditions, and multiple traps, with great exploration potential. However, due to the shielding effect of shallow natural gas and diapir activity, the imaging quality of deep-layer seismic data in the diapir fuzzy zone is poor, the signal-to-noise ratio is low, and the stratigraphic contact relationship is unclear. It is difficult to finely implement the deep-layer traps based on the existing imaging results, and it is difficult to accurately evaluate the effectiveness of the traps. In order to improve the imaging of the fuzzy zone and the sand body contact relationship, improve the reservoir prediction accuracy and release production capacity, a technology that can improve the imaging quality of the diapir fuzzy zone is urgently needed.

[0003] Reverse time migration imaging technology is one of the most accurate imaging methods currently. It can effectively improve the imaging quality of steep-angle structures and has been widely used in seismic imaging. Since the actual earth is a non-perfectly elastic medium, the kinetic energy of seismic waves is gradually converted into thermal energy during propagation, resulting in amplitude attenuation and phase dispersion problems. Conventional acoustic reverse time migration imaging is based on the acoustic wave assumption, ignoring the effect of formation viscosity on seismic waves, resulting in blurred imaging results, especially in areas with strong viscosity, which seriously affects the fine characterization of effective reservoirs. Therefore, in actual exploration, it is necessary to fully consider the absorption and attenuation properties of the formation and compensate for them.

[0004] Attenuation-compensated reverse time migration imaging can effectively improve the imaging quality in high-attenuation fuzzy areas. At present, attenuation-compensated migration imaging mainly includes two categories: the first category uses wave equations based on physical models, such as standard linear forward equations. However, this type of method usually requires a lot of computational cost, which is not conducive to actual three-dimensional exploration. In addition, the amplitude attenuation and phase dispersion operators in the equation are coupled with each other. When the amplitude is compensated by changing the sign of the amplitude attenuation term, it usually causes phase shift phenomenon, and the reflection axis cannot be correctly reset. In contrast, the wave equation based on the fractional-order Laplace operator has a relatively low computational cost, and the amplitude attenuation and phase dispersion operators can be independently decoupled and represented, so the phase shift problem is avoided when compensating the amplitude. In theory, changing the sign of the amplitude attenuation term in the wave equation can achieve amplitude compensation. However, directly changing the sign will cause the high-frequency divergence of its Green's function at large time points, thereby generating high-frequency noise in the numerical simulation process, affecting the final imaging accuracy. Introducing a low-pass filter in the numerical simulation process can suppress high-frequency noise, but the low-pass filter is a global filter operator in the time domain, which usually damages the effective signal and also reduces the reliability of the imaging results. Summary of the invention

[0005] In order to overcome the defect of insufficient imaging quality of the diapir fuzzy zone in the above-mentioned prior art, the present invention provides an adaptive attenuation compensation imaging method for the diapir zone, thereby improving the imaging fineness of the diapir fuzzy zone and improving the imaging quality.

[0006] In order to solve the above technical problems, the technical solution adopted by the present invention is:

[0007] An adaptive attenuation compensation imaging method for a diapir region comprises the following steps:

[0008] Step 1: Read the information of earthquake source, detection point, earthquake record, etc. from the earthquake data volume, and generate Ricker wavelet through earthquake source and detection point;

[0009] Step 2: Set the calculation simulation time;

[0010] Step 3: Taking the Ricker wavelet as the source, numerical simulation is performed along the positive direction of the simulation time to obtain the forward wave field, and part of the forward wave field is stored. The calculation process is from the starting point to the end point of the simulation time; taking the earthquake record as the source, numerical simulation is performed along the reverse direction of the simulation time to obtain the reverse wave field. The calculation process is from the end point to the starting point of the simulation time;

[0011] Step 4: Use the stored partial forward wave field to perform numerical simulation calculations to obtain the reconstructed source wave field;

[0012] Step 5: Perform cross-correlation imaging on the back propagation wave field obtained in step 3 and step 4 and the reconstructed source wave field. The process starts from the starting point to the end point of the simulation time and obtains the imaging result.

[0013] By performing numerical simulation calculations along the positive direction of the simulation duration to obtain the forward wave field, and storing part of the forward wave field, using this part of the forward wave field to obtain the reconstructed source wave field, and then performing numerical simulation calculations along the reverse direction of the simulation duration to obtain the reverse wave field, and using the reconstructed source wave field and the reverse wave field for correlation imaging, the imaging results can be obtained, which can effectively suppress the high-frequency numerical instability and imaging phase shift problems during numerical simulation, improve the imaging quality of the high attenuation area of ​​the diapir, and meet the imaging requirements of subsequent reservoir prediction and hydrocarbon detection. Specifically, after the source data are distributed along the time axis of the simulation duration, they are calculated in the forward and reverse directions of the time axis.

[0014] Preferably, in step three, an adaptive attenuation compensation forward operator is used to calculate the forward wavefield and the backward wavefield.

[0015] Preferably, in step 4, the stored partial forward wavelength is used through the adaptive attenuation compensation forward operator, and the reconstructed source wave field is calculated in the reverse direction of the simulation duration, and the calculation process starts from the end point of the simulation duration to the starting point.

[0016] Preferably, the formula of the adaptive attenuation compensation forward operator is: in, In the formula, represents the pressure wave field, c represents the longitudinal wave velocity of the medium, c0 is the phase velocity at the reference frequency ω0, t represents time, Q represents the quality factor, which represents the absorption and attenuation ability of the medium to seismic waves, is the amplitude compensation term, is the compensation stability term, β and σ are the stability compensation parameters, and γ are intermediate variables, and k represents the wave number.

[0017] The smaller the value of Q, the greater the attenuation intensity of seismic waves, and vice versa. and They are the amplitude compensation term and the compensating stability term, which can suppress high-frequency numerical instability and are the designed amplitude stability operators. β can retain more high-frequency information, and as σ increases, more high-wavenumber information is suppressed.

[0018] The adaptive stable compensation forward operator is used to accurately compensate the amplitude while avoiding the phase shift problem. In order to suppress the high-frequency numerical instability problem in numerical simulation, an amplitude stabilization operator is designed, which can adaptively suppress high-frequency noise according to the propagation time in numerical simulation, avoid numerical instability problems, and thus improve the final imaging accuracy.

[0019] Preferably, the fractional-order Laplace operator in the adaptive attenuation compensation forward operator is obtained using a fast Fourier transform based on an information transfer interface (MPI), and the transformation formula is: In the formula, F and F -1 They represent forward and inverse Fourier transform respectively.

[0020] Preferably, in step three, a hybrid absorption boundary condition based on a transmission boundary is used to suppress boundary reflection.

[0021] Since computers can only simulate limited space, boundary reflections are bound to occur. Therefore, a mixed absorption boundary condition based on a transmission boundary is used to suppress boundary reflections.

[0022] Preferably, in step 5, the back propagation wave field obtained in step 3 and step 4 and the reconstructed source wave field are cross-correlated imaged using an imaging formula, and the imaging formula is: In the formula, I p Represents imaging results, p s represents the reconstructed source wavefield, p r represents the back propagation wave field, X represents the spatial position, and Tmax is the maximum propagation time of the seismic wave.

[0023] Although reverse time migration imaging has been widely used in the field of seismic data processing, it still faces problems such as large amount of calculation and frequent disk I / O. In actual exploration, the problem of computational efficiency seriously limits the smooth promotion of migration imaging. To solve this problem, this method further applies multi-threaded multi-graphics processor (GPU) acceleration and MPI parallel acceleration technology, which can greatly improve the computational efficiency and meet the computational needs.

[0024] Preferably, the simulation duration is divided into sampling times according to the sampling rate of whole milliseconds. Further, the simulation duration can be set according to the focal depth, for example, a focal depth of 1 km corresponds to a simulation duration of 1 second.

[0025] Preferably, the unit sampling rate ranges from 1 to 4 ms.

[0026] Preferably, in step 1, a Ricker wavelet of 15-35 Hz is used as the seismic source.

[0027] The image accuracy calculated by seismic sources of different frequencies is different. The frequency of 20-25hz is conducive to high-precision imaging and has relatively high imaging accuracy.

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

[0029] 1. Imaging using this method can effectively suppress high-frequency numerical instability and imaging phase shift problems, improve the imaging quality of high-attenuation areas in the diapir, and meet the imaging requirements of subsequent reservoir prediction and hydrocarbon detection.

[0030] 2. Use the adaptive stable compensation forward operator to accurately compensate the amplitude while avoiding the phase shift problem. In order to suppress the high-frequency numerical instability problem in numerical simulation, an amplitude stability operator is designed. This operator can adaptively suppress high-frequency noise according to the propagation time in numerical simulation, avoid numerical instability problems, and thus improve the final imaging accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] Figure 1 It is a schematic flow chart of a method for adaptive attenuation compensation imaging of a diapir area according to the present invention;

[0032] Figure 2 is a model parameter diagram of Example 2 of an adaptive attenuation compensation imaging method for a diapir region of the present invention;

[0033] Figure 31 is a numerical simulation result diagram of Example 2 of an adaptive attenuation compensation imaging method for a diapir zone of the present invention, wherein (a) is the true solution; (b) is the traditional method; (c) is the method proposed by the present application; (d) is the difference between the true solution and the traditional method; (e) is the difference between the true solution and the proposed method;

[0034] Figure 4 It is a model parameter diagram of Example 3 of an adaptive attenuation compensation imaging method for a diapir zone of the present invention.

[0035] Figure 5 1 is an imaging result diagram of different methods of Example 3 of an adaptive attenuation compensation imaging method for a diapir zone of the present invention, wherein (a) is an acoustic wave reference solution; (b) is an imaging result without compensation; (c) is an imaging result with direct compensation; and (d) is a method proposed in the present application;

[0036] Figure 6 It is an amplitude energy plane diagram of Example 4 of an adaptive attenuation compensation imaging method for a diapir zone of the present invention, wherein (a) shallow layer; (b) middle layer; (c) deep layer;

[0037] Figure 7 1 is a diagram of actual data imaging results of Example 4 of an adaptive attenuation compensation imaging method for a diapir zone of the present invention, wherein (a) is a conventional method; (b) is a method proposed in the present application;

[0038] Figure 8 is a spectrum diagram of imaging results of Example 4 of an adaptive attenuation compensation imaging method for a diapir region of the present invention;

[0039] Fig. 9 It is another flow chart of the adaptive attenuation compensation imaging method of a diapir area of ​​the present invention. DETAILED DESCRIPTION

[0040] The drawings are only for illustrative purposes and cannot be construed as limiting the present invention. To better illustrate the present embodiment, some parts of the drawings may be omitted, enlarged, or reduced, and do not represent the size of the actual product. For those skilled in the art, it is understandable that some well-known structures and their descriptions may be omitted in the drawings. The positional relationships described in the drawings are only for illustrative purposes and cannot be construed as limiting the present invention.

[0041] The same or similar numbers in the drawings of the embodiments of the present invention correspond to the same or similar parts; in the description of the present invention, it should be understood that if the terms "upper", "lower", "left", "right", "long", "short" and the like indicate orientations or positional relationships based on the orientations or positional relationships shown in the drawings, they are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation. Therefore, the terms describing the positional relationship in the drawings are only used for illustrative purposes and cannot be understood as limitations on this patent. For ordinary technicians in this field, the specific meanings of the above terms can be understood according to specific circumstances.

[0042] The technical solution of the present invention is further described in detail below through specific embodiments and in conjunction with the accompanying drawings:

[0043] Example 1

[0044] like Figure 1 and 9 As shown, an adaptive attenuation compensation imaging method for a diapir region comprises the following steps:

[0045] Step 1: Read the information of earthquake source, detection point, earthquake record, etc. from the earthquake data volume, and generate Ricker wavelet through earthquake source and detection point;

[0046] Step 2: Set the calculation simulation time;

[0047] Step 3: Taking the Ricker wavelet as the source, the adaptive attenuation compensation forward operator is used to perform numerical simulation calculations along the positive direction of the simulation time to obtain the forward wave field, and store part of the forward wave field. The calculation process is from the starting point to the end point of the simulation time; taking the earthquake record as the source, the adaptive attenuation compensation forward operator is used to perform numerical simulation calculations along the reverse direction of the simulation time to obtain the reverse wave field. The calculation process is from the end point to the starting point of the simulation time. The formula of the adaptive attenuation compensation forward operator is: in, In the formula, represents the pressure wave field, c represents the longitudinal wave velocity of the medium, c0 is the phase velocity at the reference frequency ω0, t represents time, Q represents the quality factor, which represents the absorption and attenuation ability of the medium to seismic waves, is the amplitude compensation term, is the compensation stability term, β and σ are the stability compensation parameters, and γ are intermediate variables, and k represents the wave number;

[0048] Step 4: using the stored partial forward wavelength through the adaptive attenuation compensation forward operator, and calculating in the reverse direction of the simulation time to obtain the reconstructed source wave field, the calculation process is from the end point of the simulation time to the starting point;

[0049] Step 5: The back propagation wave field obtained in step 3 and step 4 and the reconstructed source wave field are cross-correlated imaged using the imaging formula. The process starts from the starting point to the end point of the simulation time and the imaging result is obtained. The imaging formula is: In the formula, I p Represents imaging results, p s represents the reconstructed source wavefield, p r represents the back propagation wave field, X represents the spatial position, and Tmax is the maximum propagation time of the seismic wave.

[0050] By performing numerical simulation calculations along the positive direction of the simulation duration to obtain the forward wave field, and storing part of the forward wave field, using this part of the forward wave field to obtain the reconstructed source wave field, and then performing numerical simulation calculations along the reverse direction of the simulation duration to obtain the reverse wave field, and using the reconstructed source wave field and the reverse wave field for correlation imaging, the imaging results can be obtained, which can effectively suppress the high-frequency numerical instability and imaging phase shift problems during numerical simulation, improve the imaging quality of the high attenuation area of ​​the diapir, and meet the imaging requirements of subsequent reservoir prediction and hydrocarbon detection. Specifically, after the source data are distributed along the time axis of the simulation duration, they are calculated in the forward and reverse directions of the time axis.

[0051] The smaller the value of Q, the greater the attenuation intensity of seismic waves, and vice versa. and They are the amplitude compensation term and the compensation stability term, which can suppress high-frequency numerical instability and are the designed amplitude stability operator. β can retain more high-frequency information, and as σ increases, more high-wave number information is suppressed. The adaptive stability compensation forward operator can accurately compensate for the amplitude while avoiding the phase shift problem. In order to suppress the high-frequency numerical instability problem in numerical simulation, an amplitude stability operator is designed. This operator can adaptively suppress high-frequency noise according to the propagation time in the numerical simulation, avoid numerical instability problems, and thus improve the final imaging accuracy.

[0052] Although reverse time migration imaging has been widely used in the field of seismic data processing, it still faces problems such as large amount of calculation and frequent disk I / O. In actual exploration, the problem of computational efficiency seriously limits the smooth promotion of migration imaging. To solve this problem, this method further applies multi-threaded multi-graphics processor (GPU) acceleration and MPI parallel acceleration technology, which can greatly improve the computational efficiency and meet the computational needs.

[0053] Furthermore, the fractional-order Laplace operator in the adaptive attenuation compensation forward operator is obtained using the fast Fourier transform based on the information passing interface (MPI), and the transformation formula is: In the formula, F and F -1They represent forward and inverse Fourier transform respectively.

[0054] Furthermore, in step 3, a hybrid absorbing boundary condition based on a transmission boundary is used to suppress boundary reflection. Since a computer can only simulate a limited space, boundary reflection is bound to occur, so a hybrid absorbing boundary condition based on a transmission boundary is used to suppress boundary reflection.

[0055] Furthermore, the simulation duration is divided into sampling times according to the sampling rate of whole milliseconds. Furthermore, the simulation duration can be set according to the corresponding setting of the focal depth, for example, a focal depth of 1 kilometer corresponds to a simulation duration of 1 second.

[0056] Furthermore, the unit sampling rate ranges from 1 to 4 ms.

[0057] Furthermore, in step 1, a Ricker wavelet of 15-35 Hz is used as the earthquake source.

[0058] The image accuracy calculated by seismic sources of different frequencies is different. The frequency of 20-25hz is conducive to high-precision imaging and has relatively high imaging accuracy.

[0059] Fig. 9 In the calculation, the maximum time and the minimum time refer to the forward and reverse calculations according to the time axis of the simulation duration until the calculation is completed, and the maximum number of shots refers to the p in the imaging calculation formula. s and p r Whether the number of inputs is reached, that is, whether all p s and p r Calculation completed.

[0060] The beneficial effects of this embodiment are as follows: a forward wave field is obtained by performing numerical simulation calculations along the positive direction of the simulation duration, and a part of the forward wave field is stored. The part of the forward wave field is used to obtain a reconstructed source wave field, and then numerical simulation calculations are performed in the reverse direction of the simulation duration to obtain a reverse wave field. The reconstructed source wave field and the reverse wave field are correlated for imaging to obtain imaging results. During the numerical simulation, high-frequency numerical instability and imaging phase shift problems can be effectively suppressed, thereby improving the imaging quality in the high attenuation zone of the diapir and meeting the imaging requirements of subsequent reservoir prediction and hydrocarbon detection.

[0061] Example 2

[0062] A two-layer model is used to test the accuracy of numerical simulation of the proposed high-precision adaptive attenuation compensation forward operator.

[0063] The verification model contains 400 × 400 grid points with a grid spacing of 10 m. The Ricker wavelet with a main frequency of 20 Hz is used as the excitation source. The time sampling rate is 1 ms. The velocity and Q parameters are as follows: Figure 2The numerical simulation results of three methods were tested, namely the accurate real solution, the traditional attenuation compensation method and the adaptive attenuation compensation method proposed in this patent. Figure 3 a-3c are the wave field simulation results of different methods at 700ms. Figure 3 d is the difference between the true solution and the traditional method, Figure 3 e is the difference between the true solution and the method proposed in this patent. Figure 3 In a-3c, the three numerical simulation methods can effectively simulate the propagation of seismic waves, but Figure 3 There is a significant simulation error in d; in contrast, Figure 3 Almost no residual is observed in e, which proves the accuracy and superiority of the proposed compensation equation in numerical simulation. Since the numerical simulation error is positively correlated with the simulation time, as the seismic wave gradually propagates, the traditional method will produce more cumulative errors, which in turn affects the accuracy of subsequent imaging.

[0064] The remaining features and working principles of this embodiment are consistent with those of Embodiment 1.

[0065] Example 3

[0066] In order to verify the suppression effect and effectiveness of the method proposed in this patent on high-frequency noise, a theoretical model is used for verification.

[0067] The verification model contains 400 and 150 grid points, with a grid spacing of 10 m. The Ricker wavelet with a main frequency of 20 Hz is used as the excitation source, the time sampling rate is 1 ms, 20 sources and 400 detectors are evenly arranged on the surface, and the velocity and Q parameters are as follows: Figure 4 As shown, Figure 4 a and 4b represent the velocity model and attenuation model respectively. A total of 4 methods were tested: Figure 5 aSound wave reference solution, Figure 5 b Imaging results without compensation, Figure 5 c Imaging results of direct compensation, Figure 5 d Imaging results proposed in this patent. In order to suppress boundary reflection, 20 layers of mixed absorption boundary conditions are introduced at the four boundaries. Affected by the absorption attenuation effect of the stratum, Figure 5 In b, we can observe obvious amplitude reduction, and the deep structure cannot be clearly presented; Figure 5 In c, the high-frequency noise caused by amplitude compensation seriously reduces the signal-to-noise ratio of the imaging result, making it difficult to observe the effective phase axis. In contrast, Figure 5 The amplitude attenuation effect in d is effectively compensated without generating obvious high-frequency noise. The imaging result is very close to the reference solution, which proves the effectiveness of the proposed method.

[0068] The rest of the working principle and working process of this embodiment are consistent with those of Embodiment 1.

[0069] Example 4

[0070] The seismic data of a diapir fuzzy zone are used to demonstrate the effectiveness of the proposed method in actual exploration.

[0071] Figure 6 is the influence range of the diapir fuzzy area on imaging at different depths, where Figure 6 a-6c represent the impact range of shallow, middle and deep layers respectively. The darker the color, the greater the amplitude energy. As shown in the figure, with the increase of depth, the impact range of the diapir fuzzy zone gradually increases. This is because the absorption and attenuation effect of the formation will become stronger with the increase of propagation time. Therefore, when compensating the amplitude, it is more reasonable to gradually increase the compensation intensity with the increase of extension time. The effective area of ​​this data is 400km 2 Taking a two-dimensional survey line in the work area as an example, the effectiveness and superiority of the proposed method in processing the actual work area are demonstrated. The horizontal and vertical lengths of the survey line are 8km and 2s respectively. Eight seismic sources are arranged on the sea surface, and 400 geophones are arranged on the seabed. The Ricker wavelet with a main frequency of 25Hz is determined as the seismic source by well seismic calibration. Figure 7 are the imaging results of different methods, among which, Figure 7 a is the imaging result of the traditional method, Figure 7 b is the high-precision adaptive attenuation compensation imaging method proposed in this patent. As shown in the figure, the traditional method has energy distortion in the shallow, middle and deep layers of the diapir core to varying degrees, which increases the uncertainty of reservoir prediction and hydrocarbon detection to varying degrees. In contrast, the proposed method effectively improves the signal-to-noise ratio of the seismic profile, and the lateral continuity of the phase axis is significantly enhanced, especially in the diapir core. Figure 8 The spectra of imaging results of different methods are shown in Figure 2, where the dotted line represents the proposed method and the solid line represents the traditional method. Figure 8 As shown in the figure, compared with the traditional method, the proposed method effectively compensates for the absorption attenuation effect of the formation and improves the amplitude intensity. Therefore, the method proposed in this patent can effectively improve the imaging quality of the diapir fuzzy area and the contact relationship of the sand body, as well as the problems of weak signals, poor imaging, and large lateral changes in velocity in the middle and deep layers, providing a reliable data basis for the implementation of the trap and improving the prediction accuracy of reservoir prediction and hydrocarbon detection.

[0072] The rest of the working principle and working process of this embodiment are consistent with those of Embodiment 1 or Embodiment 2.

[0073] Obviously, the above embodiments of the present invention are merely examples for clearly illustrating the present invention, and are not intended to limit the embodiments of the present invention. For those skilled in the art, other different forms of changes or modifications can be made based on the above description. It is not necessary and impossible to list all the embodiments here. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the claims of the present invention.

Claims

1. An adaptive attenuation compensation imaging method for a diapir region, characterized in that: The following steps are involved: Step 1: Read the information of earthquake source, detection point, earthquake record, etc. from the earthquake data volume, and generate Ricker wavelet through earthquake source and detection point; Step 2: Set the calculation simulation time; Step 3: Using the Ricker wavelet as the source, numerical simulation is performed along the positive direction of the simulation duration to obtain the forward wave field, and part of the forward wave field is stored. The calculation process starts from the starting point to the end point of the simulation duration. Taking the earthquake record as the source, numerical simulation is performed in the reverse direction of the simulation time to obtain the reverse wave field. The calculation process starts from the end point of the simulation time to the starting point. Step 4: Use the stored partial forward wave field to perform numerical simulation calculations to obtain the reconstructed source wave field; Step 5: Perform cross-correlation imaging on the back propagation wave field obtained in step 3 and step 4 and the reconstructed source wave field. The process starts from the starting point to the end point of the simulation time and obtains the imaging result.

2. The method for adaptive attenuation compensation imaging of a diapir region according to claim 1, characterized in that: In step 3, the forward wavefield and the reverse wavefield are calculated using the adaptive attenuation compensation forward operator.

3. The method for adaptive attenuation compensation imaging of a diapir region according to claim 2, characterized in that: In step 4, the stored partial forward wavelength is used through the adaptive attenuation compensation forward operator, and the reconstructed source wave field is calculated in the reverse direction of the simulation time, and the calculation process starts from the end point of the simulation time to the starting point.

4. The method for adaptive attenuation compensation imaging of a diapir region according to claim 2, characterized in that: The formula of the adaptive attenuation compensation forward operator is: in, In the formula, represents the pressure wave field, c represents the longitudinal wave velocity of the medium, c0 is the phase velocity at the reference frequency ω0, t represents time, Q represents the quality factor, which represents the absorption and attenuation ability of the medium to seismic waves, is the amplitude compensation term, is the compensation stability term, β and σ are the stability compensation parameters, and γ are intermediate variables, and k represents the wave number.

5. The method for adaptive attenuation compensation imaging of a diapir region according to claim 4, characterized in that: The fractional-order Laplace operator in the adaptive attenuation compensation forward operator is obtained using a fast Fourier transform based on the information passing interface (MPI), and the transformation formula is: In the formula, F and F -1 They represent forward and inverse Fourier transform respectively.

6. The method for adaptive attenuation compensation imaging of a diapir region according to claim 5, characterized in that: In step 3, a hybrid absorbing boundary condition based on a transmission boundary is used to suppress boundary reflection.

7. The method for adaptive attenuation compensation imaging of a diapir region according to claim 1, characterized in that: In step 5, the back propagation wave field obtained in step 3 and step 4 and the reconstructed source wave field are cross-correlated imaged using an imaging formula, and the imaging formula is: In the formula, I p Represents imaging results, p s represents the reconstructed source wavefield, p r represents the back propagation wave field, X represents the spatial position, and Tmax is the maximum propagation time of the seismic wave.

8. The method for adaptive attenuation compensation imaging of a diapir region according to claim 3, characterized in that: The simulation duration is divided into sampling times according to a sampling rate in units of whole milliseconds.

9. The method for adaptive attenuation compensation imaging of a diapir region according to claim 8, characterized in that: The unit sampling rate ranges from 1 to 4 ms.

10. The method for adaptive attenuation compensation imaging of a diapir region according to claim 1, characterized in that: In step 1, a 15-35 Hz Ricker wavelet is used as the source.

Images