A Viscous Medium Target Area Imaging Method Based on Virtual Detection Point Dyeing Algorithm

Abstract

The present invention belongs to the technical field of data processing, and particularly relates to a viscous medium target area imaging method based on a virtual geophone point coloring algorithm. 1. Set the coloring position; 2. Set the source, geophone, and seismic recording information; 3. Obtain the source seismic wavefield value; 4. Save the calculated true wavefield; 5. Calculate the true wavefield; 6. Determine whether the time loop has reached the maximum time step T. If not, repeat steps 3 to 5. If so, end the time loop; 7. Assign the wavefield value of the true wavefield reaching the coloring area to the coloring wavefield; 8. Calculate the true wavefield; 9. Determine whether the time loop has reached the maximum time step T. If not, repeat steps 7 to 8. If so, end the time loop; 10. Load the geophone wavefield data; 11. Solve the geophone seismic wavefield value; 12. Calculate the geophone seismic wavefield value and determine whether the time step 0 has been reached. If not, repeat steps 7 to 9. If so, end the loop; 13. Obtain the underground structure imaging result after processing using the imaging condition.

Description

Technical Field

[0001] The present invention belongs to the technical field of data processing, and particularly relates to a method for imaging a target area of a viscous medium based on a virtual geophone point coloring algorithm. Background Art

[0002] In oil exploration, understanding the subsurface medium conditions plays a crucial guiding role in oil development and production. Before the advent of seismic exploration, geologists inferred the location of oil and gas fields by studying surface outcrops to find oil and gas seeps and anticlinal structures. The emergence of seismic exploration technology from 1920 to 1930 shifted geological surveys from the surface to the subsurface. Seismic data interpreters searched for existing oil and gas reservoirs by looking for structural traps based on the interpretation results of seismic data. After decades of development, the oil and gas contained in conventional structural oil and gas reservoirs have been almost completely developed. The advent of the theory of subtle oil and gas reservoirs has provided a broader prospect for oil and gas exploration and development. During the research process of subtle lithologic oil and gas reservoirs, geophysicists found that large oil and gas reservoirs are likely to form in the subsalt area. However, the limited observation system is difficult to receive enough effective reflected waves. And when seismic waves propagate to the surface of a salt dome with a high reflectivity, the waveform will be greatly distorted, and the amplitude will also be greatly attenuated, making the deep energy complex and weak, resulting in an extremely low signal-to-noise ratio of seismic waves in the subsalt area. Therefore, developing precise imaging technology for the subsalt area has become the consensus in the current industry.

[0003] To solve this problem, geophysicists have proposed quite a number of solutions. For example, three-dimensional seismic logging technology that changes the observation system can receive seismic signals in multiple azimuths and multiple scales, greatly enriching the illumination of the underground structure. It is an effective imaging method for the subsalt area (Chen, 1992). A wide-azimuth three-dimensional acquisition system can collect rich azimuth information, which is beneficial for imaging deep layers, high-steep structures, and anisotropic rock masses (Tang, 2003). The acquisition aperture correction in the local angle domain can also significantly improve the amplitude of the imaging results (Cao, 2009). Weighting the weakly illuminated areas while retaining the strongly illuminated areas improves the imaging effect in the subsalt area (Gherasim, 2014). Since the reflection paths of multiple reflection waves are different from those of primary reflection waves, multiple waves can also be used to enhance the energy signal in the weakly illuminated subsalt area (Li Peng, 2006). Multiple scattering waves are also used to image areas lacking primary reflection wave illumination (Guitton, 2002). Taking multiple waves as an important reflection energy and using the imaging results of multiple waves to correct the results of traditional reverse time migration can enhance the illumination effect in the deep subsalt area (Liu, 2011). In addition, a relatively accurate velocity model is also the key to obtaining high-quality subsalt imaging results. Therefore, accurate velocity modeling technology also plays a crucial role in imaging the weakly illuminated subsalt area. Tang et al. (2011) extracted a subset of the total dataset during the velocity modeling process. However, this sub-dataset contains all the information required for velocity modeling, making this sub-dataset suitable for constructing a velocity model for a specific target structure. Inspired by this, Chen et al. (2014) proposed the elastic medium complex domain coloring algorithm. Its inspiration comes from the idea of the "fate map" in biology. The colored wave field related to the target area is regarded as "cells". By marking these "cells", their "life activities" - the reflected waves and transmitted waves generated by the colored wave field can be observed. Applying this algorithm can obtain a colored wave field related to the target area, and this wave field is synchronous with the real wave field during the propagation process. Applying the obtained colored wave field to the migration imaging condition can obtain an imaging result that erases everything except the target structure. Therefore, this algorithm can improve the signal-to-noise ratio of the target area. However, the colored wave field obtained by the complex domain coloring algorithm is affected by the complex velocity and has the disadvantages that the amplitude is much smaller than that of the real wave field and the waveform is distorted. Li et al. (2017) proposed the control equation of the elastic medium generalized coloring algorithm. The constructed colored wave field matches the real wave field much better than the complex domain coloring algorithm, and can achieve high-resolution imaging of the target area.

[0004] However, due to theoretical and technical limitations and considerations of production efficiency, in early seismic exploration research, the underground medium was simply simplified to a non-attenuating elastic medium. In fact, the underground medium usually has viscosity, and amplitude attenuation, phase distortion, and frequency band narrowing will occur during the propagation of seismic waves in the formation. If the formation viscosity is ignored, it will lead to more attenuation of seismic wave energy and waveform distortion in the subsalt area, seriously affecting the imaging accuracy of the underlying structure.

[0005] Currently, the generalized standard linear solid model (GSLS) that uses memory variables and viscoelastic parameters to characterize the wave field characteristics has been applied to describe the viscosity of the formation (Carcione, 1988). However, the memory variables used in this equation implicitly characterize the amplitude attenuation and phase dislocation phenomena. Although the seismic wave amplitude information can be restored in the compensated migration imaging, the phase information cannot be restored, resulting in inaccurate formation positions in the imaging results. And due to the excessive number of parameters characterized by the memory variables, it is not an ideal wave field continuation operator. In contrast, the commonly studied constant-Q model (Kjartansson, 1979) today has fewer characterized parameters and simpler calculations. According to the definition of this model, within the seismic exploration frequency band (<150HZ), the quality factor Q hardly changes with frequency, which is more in line with physical laws.

[0006] Carcion et al. (2002) based on the constant-Q model, first realized the application of this model in seismic wave propagation. Carcion et al. (2010) developed an approximate constant-Q wave equation for the large storage and calculation costs required for the fractional-order time partial derivative in the calculation equation. This equation contains a fractional-order Laplace operator and can be calculated by the fast Fourier transform, greatly simplifying the numerical calculation. On this basis, Zhu et al. (2014) proposed a viscoacoustic equation with an approximate constant-Q fractional-order Laplace operator. This equation contains two Laplace operators corresponding to the phase dispersion term and the amplitude attenuation term respectively, and its characterized parameters are only velocity and Q, which is a good method for realizing viscoacoustic reverse time migration. However, this method still does not solve the problems such as insufficient imaging illumination and low signal-to-noise ratio in the presence of salt dome overlying structures. Summary of the Invention

[0007] In order to solve the problems described in the background art, the present invention provides a viscosity medium target area imaging method based on a generalized coloring algorithm, which can improve the recognition accuracy of complex subsalt oil and gas reservoirs and provide core technical support for the rapid and efficient development of oil and gas resources.

[0008] The technical solution adopted by the present invention is: a viscosity medium target area imaging method based on a virtual geophone coloring algorithm, and the method includes the following steps:

[0009] Step 1. Set the dyeing position, and set N virtual detection points at the selected dyeing position to receive the wave field values propagated to the dyeing position;

[0010] Step 2. Set the source, detection points, and seismic record information;

[0011] Step 3. Use the fractional Laplacian operator constant-Q decoupled visco-acoustic wave equation and boundary conditions to obtain the source seismic wave field values at time step T, and save the wave field values;

[0012] Fractional Laplacian operator constant-Q decoupled visco-acoustic wave equation:

[0013]

[0014]

[0015] In the equation: c is the velocity, t is the time, p is the visco-acoustic wave field; Q is the quality factor; c 0 is the phase velocity defined at the reference frequency ω 0 ; s is the source; the first term on the right side of equation (1) mainly controls dispersion, and the second term mainly controls amplitude attenuation;

[0016] Step 4. Save the true wave field calculated by the fractional Laplacian operator constant-Q decoupled visco-acoustic wave equation at time step t

[0017] Step 5. Calculate the true wave field at time step t + 1;

[0018] Step 6. Determine whether the time loop has reached the maximum time step T. If not, repeat steps 3 to 5. If so, end the time loop;

[0019] Step 7. Assign the part of the wave field values of the saved true wave field reaching the dyeing area to the dyeing wave field;

[0020] Step 8. Calculate the true wave field at time step j + 1 according to the virtual detection point dyeing algorithm equation, and save the dyeing wave field

[0021] Generalized dyeing algorithm equation for viscous media:

[0022]

[0023]

[0024]

[0025] In the equation and represent the original wave field and the true source, and represents the colored wave field and the wave field value received by the virtual geophone point;

[0026] Step Nine: Determine whether the time loop has reached the maximum time step T. If not, repeat Steps Seven to Eight; if yes, end the time loop;

[0027] Step Ten: Load the wave field data of the geophone points with a time step of T + 1;

[0028] Step Eleven: Solve the seismic wave field value of the geophone points at time step T + 1 using the wave equation (1) and the boundary conditions;

[0029] Step Twelve: Calculate the seismic wave field value of the geophone points at time step T - t - 1, and determine whether it has reached time step 0. If not, repeat Steps Seven to Nine; if yes, end the loop;

[0030] Step Thirteen: Use the imaging condition to process the colored wave fields of all shots, the reverse propagated wave field information of the geophone points, and obtain the underground structure imaging result after denoising;

[0031] The imaging condition is as follows:

[0032]

[0033] where R is the reverse propagated wave field of the geophone point, is the colored wave field.

[0034] Furthermore, the finite difference method is used to solve the time derivative on the left side of the equation, and the pseudo - spectral method is used to obtain the spatial operator on the right side of the equation. Therefore, the time accuracy is second - order accuracy, and the spatial accuracy is spectral accuracy. Its calculation equation is:

[0035]

[0036]

[0037] where F and F -1 are the one - dimensional forward and inverse Fourier transforms respectively, and k is the discrete wave number.

[0038] Furthermore, N virtual geophone points are set at the coloring position in Step One, where N is greater than 80, and the number of coloring layers is not less than 5.

[0039] Furthermore, the setting of the coloring area conforms to the Kirchhoff integral equation.

[0040] When the real wave field propagates to the coloring area D, the Kirchhoff integral equation gives the solution of equation (2);

[0041]

[0042] where r‵ is the source position vector in space, r 0 is the vector from the origin of coordinates to any point on the integration surface D, r is the coordinate vector of any point outside the space, t 0 is the time point for recording the wave field, D = 1 represents the integration surface, which is also the boundary of the space, n is the normal vector of the boundary D, and G is the Green's function of the viscoacoustic equation;

[0043] where:

[0044]

[0045] is the true wave field obtained by the convolution of the source and the Green's function, then

[0046]

[0047] Equation (7) shows that the colored wave field obtained by the above method is the integral of the true wave field on the boundary D; that is, when the true wave field on the boundary D and its spatial derivative are known, the wave field at spatial positions on the side different from the source can be obtained, which conforms to the Huygens principle;

[0048] The colored wave field of the virtual geophone in the viscous medium conforms to the above conclusion, and constructing the colored wave field can be achieved by conforming to this conclusion:

[0049]

[0050] In the equation, is the true wave field, is the colored wave field, c is the velocity, t is the time, and Q is the quality factor; c 0 is the phase velocity defined at the reference frequency ω 0 ; s is the source,

[0051] D = 1 represents the coloring position.

[0052] Furthermore, the fractional-order viscoacoustic equation realizes the decoupling of the phase dispersion term and the amplitude attenuation term, and can change the sign of the amplitude attenuation term during the wave field propagation process to achieve compensation; during the compensation process, the compensation operator increases exponentially with the increase of the wave number, and the high-frequency components are severely amplified, resulting in numerical instability; to suppress the numerical instability phenomenon, by referring to the idea of inverse Q filtering, an adaptive stable propagation operator is used to stabilize the wave field, and its calculation equation is:

[0053]

[0054] where is the amplitude attenuation operator, and the rest is the stabilization operator, denoted as:

[0055]

[0056] Since the attenuation of seismic waves occurs during the entire wavefield propagation period, compensation should also be performed within each wavefield continuation step. For the first continuation moment, i.e., when l = 1, the stabilization coefficient is:

[0057]

[0058] For the l-th continuation moment, l ≥ 2, the stabilization coefficient is:

[0059]

[0060] σ in the equation 2 is the introduced stabilization factor, and its value is 2.5×10 -3 ;

[0061] The stabilization compensation operators of equations (11) and (12) are applied to the compensation wavefield update at each moment. Based on the constant-Q model fractional Laplacian operator viscoacoustic equation, the compensation form is:

[0062]

[0063]

[0064] In the equation, c is the velocity, t is the time, p is the viscoacoustic wavefield; Q is the quality factor; c 0 is the phase velocity defined at the reference frequency ω 0 ; s is the source.

[0065] Advantages of the present invention: A viscoelastic medium target area imaging method based on a generalized coloring algorithm is provided, which can improve the recognition accuracy of complex oil and gas reservoirs under salt, and provide core technical support for the rapid and efficient development of oil and gas resources. Its main advantages are as follows:

[0066] (1). Due to the existence of the minimum velocity in the conventional complex-domain coloring method, the amplitude value of the colored wavefield is much smaller than that of the real wavefield. The virtual geophone coloring algorithm used in this method can realize the amplitude-preserving propagation of the colored wavefield; The virtual geophone coloring algorithm used in this method can realize the amplitude-preserving propagation of the colored wavefield;

[0067] (2). In the conventional elastic wave coloring algorithm, the absorption and attenuation effect of the formation is ignored, and the propagation law of seismic waves in the subsurface cannot be accurately described. This method uses the fractional Laplacian operator constant-Q decoupled viscoacoustic equation to extend the coloring algorithm to viscoelastic media, which can more accurately depict the attenuation mechanism of seismic waves during subsurface propagation. By using the decoupling characteristics of this equation, the decoupling of amplitude and phase is beneficial to attenuation compensation reverse time migration;

[0068] (3) Introduce an adaptive stable propagation operator into the wave equation to alleviate the high-frequency outlier problem caused by direct compensation and effectively reduce the loss of effective wave information due to low-pass filtering;

[0069] (4) The reverse time migration imaging process of the viscoacoustic medium target area based on the virtual receiver coloring algorithm improves the signal-to-noise ratio and resolution of the imaging profile in the complex structure target area, which is conducive to accurately identifying underground structures and predicting reservoir intervals in seismic data interpretation;

[0070] (5) Use the imaging method of the viscoelastic medium target area based on the virtual receiver coloring algorithm to solve the problem of excessive computational cost in the conventional coloring algorithm due to the simultaneous calculation of two forward calculations. Description of the Drawings

[0071] Figure 1 is the imaging flow chart of the viscoacoustic medium target area based on the virtual receiver coloring algorithm;

[0072] Figure 2 is the wave field diagram of the virtual receiver coloring in different media;

[0073] Figure 3 is the comparison diagram of the original wave field and the colored wave field with different colored area thicknesses during the numerical solution of the pseudo-spectral method;

[0074] Figure 4 is the relative error diagram of the colored wave field corresponding to different colored area lengths during the numerical solution of the pseudo-spectral method;

[0075] Figure 5 is the horizontal layered velocity model diagram;

[0076] Figure 6 is the horizontal layered Q-value model diagram;

[0077] Figure 7 is the original wave field diagram;

[0078] Figure 8 is the colored wave field diagram obtained by the traditional coloring algorithm;

[0079] Figure 9 is the comparison diagram of the original wave field and the colored wave field obtained by the traditional coloring algorithm;

[0080] Figure 10 is the colored wave field diagram obtained by the viscoacoustic virtual receiver coloring algorithm;

[0081] Figure 11 is the comparison diagram of the original wave field and the colored wave field obtained by the viscoacoustic virtual receiver coloring algorithm;

[0082] Figure 12 is the salt dome velocity model diagram;

[0083] Figure 13 It is the Q - value model diagram of the salt dome;

[0084] Figure 14 It is the acoustic imaging diagram of the salt dome;

[0085] Figure 15 It is the attenuation imaging diagram of the salt dome;

[0086] Figure 16 It is the compensated imaging diagram of the salt dome;

[0087] Figure 17 It is the compensated conventional staining algorithm imaging diagram of the salt dome;

[0088] Figure 18 It is the compensated virtual geophone point staining algorithm imaging diagram of the salt dome;

[0089] Figure 19 It is the single - trace information diagram at the horizontal position of 2.4 km;

[0090] Figure 20 It is the double - salt - dome velocity model;

[0091] Figure 21 It is the Q - value model diagram of the double - salt dome;

[0092] Figure 22 It is the compensated imaging diagram of the double - salt dome;

[0093] Figure 23 It is the compensated virtual geophone point staining algorithm imaging diagram of the double - salt dome;

[0094] Figure 24 It is the single - trace information diagram at the horizontal position of 3.83 km;

[0095] Figure 25 It is the uniform velocity model diagram;

[0096] Figure 26 It is the sigbee model diagram. Specific implementation mode

[0097] Example 1

[0098] Referring to each figure,

[0099] A viscous medium target area imaging method based on the virtual geophone point staining algorithm, comprising the following steps:

[0100] Step 1: Set the staining position, and set N virtual geophone points at the selected staining position to receive the wave field values propagated to the staining position;

[0101] Step 2: Set the source, geophone, and seismic record information;

[0102] Step 3: Use the constant-Q decoupled visco-acoustic wave equation and boundary conditions of the fractional Laplacian operator to obtain the source seismic wave field values at time step T, and save the wave field values;

[0103] Constant-Q decoupled visco-acoustic wave equation of the fractional Laplacian operator:

[0104]

[0105]

[0106] In the equation: c is the velocity, t is the time, p is the visco-acoustic wave field; Q is the quality factor; c 0 is the phase velocity defined at the reference frequency ω 0 ; s is the source; the first term on the right side of equation (1) mainly controls dispersion, and the second term mainly controls amplitude attenuation;

[0107] Step 4: Save the true wave field calculated by the constant-Q decoupled visco-acoustic wave equation of the fractional Laplacian operator at time step t

[0108] Step 5: Calculate the true wave field at time step t + 1;

[0109] Step 6: Determine whether the time loop has reached the maximum time step T. If not, repeat Steps 3 to 5; if so, end the time loop;

[0110] Step 7: Assign the partial wave field values of the saved true wave field reaching the colored area to the colored wave field;

[0111] Step 8: Calculate the true wave field at time step j + 1 according to the virtual geophone coloring algorithm equation, and save the colored wave field

[0112] Generalized coloring algorithm equation for viscoelastic media:

[0113]

[0114]

[0115]

[0116] In the equation and represent the original wave field and the true source, and represent the colored wave field and the wave field values received by the virtual geophone;

[0117] Step 9: Determine whether the time loop has reached the maximum time step T. If not, repeat Steps 7 to 8; if so, end the time loop;

[0118] Step Ten: Load the detected wavefield data at time step T+1;

[0119] Step Eleven: Solve for the seismic wavefield values at the detection points at time step T+1 using the wave equation (1) and boundary conditions;

[0120] Step Twelve: Calculate the seismic wavefield values at the detection points at time step T-t-1, and determine whether the time step reaches 0; otherwise, repeat Steps Seven to Nine, and if so, end the loop;

[0121] Step Thirteen: Use the imaging condition to process the colored wavefields of all shots, the information of the reverse-propagated wavefields at the detection points, and after denoising, obtain the subsurface structure imaging results;

[0122] The imaging condition is as follows:

[0123]

[0124] where R is the reverse-propagated wavefield at the detection point, is the colored wavefield.

[0125] The finite difference method is used to solve the time derivative on the left side of the equal sign, and the pseudo-spectral method is used to obtain the spatial operator on the right side of the equal sign; therefore, the time accuracy is second-order accuracy, and the spatial accuracy is spectral accuracy. Its calculation equation is:

[0126]

[0127]

[0128] where F and F -1 are the one-dimensional forward and inverse Fourier transforms respectively, and k is the discrete wave number.

[0129] In Step One, N virtual detection points are set at the coloring positions, where N is greater than 80, and the number of coloring layers is not less than 5.

[0130] The setting of the coloring area conforms to the Kirchhoff integral equation,

[0131] When the true wavefield propagates to the coloring area D, the Kirchhoff integral equation gives the solution of equation (2);

[0132]

[0133] where r‵ is the source position vector in space, r 0 is the vector from the origin of coordinates to any point on the integration surface D, r is the coordinate vector of any point outside the space, t 0 is the time point for recording the wavefield, D = 1 represents the integration surface, which is also the boundary of the space, n is the normal vector of the boundary D, and G is the Green's function of the viscoacoustic equation;

[0134] wherein:

[0135]

[0136] is the true wave field obtained by convolving the seismic source with the Green's function, then

[0137]

[0138] Equation (7) shows that the colored wave field obtained by the above method is the integral of the true wave field on the boundary D; that is, when the true wave field on the boundary D is known and its spatial derivative the wave field at spatial positions on the side different from the seismic source can be obtained, which conforms to Huygens' principle;

[0139] The colored wave field of the virtual geophone in the viscous medium conforms to the above conclusion, and the construction of the colored wave field can be realized by conforming to this conclusion:

[0140]

[0141] In the equation, is the true wave field, is the colored wave field, c is the velocity, t is the time, and Q is the quality factor; c 0 is the phase velocity defined at the reference frequency ω 0 ; s is the seismic source,

[0142]

[0143] D = 1 represents the coloring position.

[0144] The fractional-order visco-acoustic equation realizes the decoupling of the phase dispersion term and the amplitude attenuation term, and can change the sign of the amplitude attenuation term to achieve compensation during the wave field propagation process; during the compensation process, the compensation operator increases exponentially with the increase of the wave number, and the high-frequency components are severely amplified, resulting in numerical instability; to suppress the numerical instability phenomenon, by referring to the idea of inverse Q filtering, an adaptive stable propagation operator is used to stabilize the wave field, and its calculation equation is:

[0145]

[0146] wherein is the amplitude attenuation operator, and the rest is the stabilization operator, denoted as:

[0147]

[0148] Since the attenuation of seismic waves occurs during the entire wave field propagation period, compensation should also be carried out within each wave field continuation step. For the first continuation moment, that is, when l = 1, the stabilization coefficient is:

[0149]

[0150] For the l-th extension moment, where l ≥ 2, the stabilization coefficient is:

[0151]

[0152] σ in the equation 2 is the introduced stabilization factor, and its value is 2.5×10 -3 ;

[0153] The stabilization compensation operators of equations (11) and (12) are applied to the compensation wavefield update at each moment. Based on the constant-Q model, the compensated form of the fractional Laplacian viscoacoustic equation is:

[0154]

[0155]

[0156] In the equation, c is the velocity, t is the time, p is the viscoacoustic wavefield; Q is the quality factor; c 0 is the phase velocity defined at the reference frequency ω 0 ; s is the source.

[0157] Table 1 - Table of Time Consumption and Memory Occupation of Different Models

[0158]

[0159] Figure 2 The four wavefield snapshots in

[0160] correspond to the acoustic wave coloring algorithm equation, the dispersion-dominated coloring algorithm equation, the attenuation-dominated coloring algorithm equation, and the viscoacoustic wave coloring algorithm equation respectively. It can be observed that, compared with the acoustic wave coloring algorithm, the dispersion-dominated coloring algorithm equation depicts the phase dislocation caused by the formation viscosity effect, the attenuation-dominated coloring algorithm equation more accurately describes the attenuation of the amplitude energy, and the viscoacoustic wave coloring algorithm equation has both effects. Figure 3 When numerically solving by the pseudospectral method, the relative error of the colored wavefield corresponding to different lengths of the coloring region. When the length of the coloring region is greater than 80 grid points, the relative error between the colored wavefield and the true wavefield will be controlled below 1%, which can effectively control the simulation error of the colored wavefield.

[0161] When calculating the spatial derivative of the wave field, the wave field values at n grid points along the wave field propagation direction are required. Assuming that a finite difference operator of order x in space is used during wave field forward modeling, x / 2 + 1 wave field values at grid points are needed to calculate this spatial derivative. And when the thickness of the colored region is greater than or equal to x / 2, the spatial derivative inside its wave field is determined and correct. After obtaining the determined and correct spatial derivative, the amplitude values and phases of the colored wave field and the true wave field can be made the same (Jia and Lu, 2017); however, in this paper, the pseudo-spectral method is used for solution. Among them, the Laplace operator is a global operator, and its spatial accuracy is spectrally accurate and approximately the same as that of a high-order finite difference in space. Its globality will inevitably introduce errors. Figure 4 The curve shown in it is the relative error of the colored wave field corresponding to colored regions with different thicknesses during the numerical solution by the pseudo-spectral method. When the number of grid points is greater than 5, the relative error can be controlled below 1%, and the simulation error of the colored wave field can be effectively controlled.

[0162] As mentioned above, there is a large difference in the amplitude between the colored wave field obtained by the conventional coloring algorithm for visco-acoustic media and the original wave field, and the phase of the colored wave field also changes. The virtual geophone coloring algorithm for visco-acoustic media can achieve the preservation of the amplitude of the colored wave field and the correction of the phase. Figure 5 , 6 are model parameters, where the dashed line is the coloring position. The wave field snapshot at a propagation time of 1 s is extracted as Figure 7 shown, where Figure 9 is the difference between the original wave field and the colored wave field obtained by the conventional coloring algorithm for visco-acoustic media. Due to the existence of the minimum velocity , the amplitude value of the conventional colored wave field is much smaller than that of the true wave field; Figure 11 is the difference between the original wave field and the colored wave field obtained by the virtual geophone coloring algorithm for visco-acoustic media. Their amplitudes and phases are the same.

[0163] In practice, if no compensation is carried out, as Figure 15 shown, the imaging effect of the deep layer under the salt will be very weak. Taking the acoustic data imaging as the reference solution, Figure 16 although the deep layer energy is compensated to a certain extent after the compensation treatment, the structure is still not clear. After applying the coloring algorithm, Figure 17 due to the existence of , the imaging amplitude and resolution of the colored wave field have a large gap with Figure 18 . Figure 18 Compared with Figure 17 , the imaging resolution of the target area is significantly improved. The structure of the target area is clear, the fault is crisp, and the break point is clear, which can better achieve the accurate imaging of the deep structure under the salt. To further illustrate the superiority of the coloring algorithm, Figure 14 , 16, 17, and 18 single-channel information at 2.4 km in the horizontal direction is extracted. Figure 19It is also possible to verify that the new staining algorithm and the staining algorithm proposed in this paper can both obtain imaging results with high imaging amplitude, high resolution, and correct structural migration.

[0164] From Figure 22 , a single-trace record is extracted from the 23-offset section, normalized, and then Figure 24 shown as Figure 24 in which the first layer position at 2 km is correctly migrated, and the imaging amplitude obtained by the staining algorithm is higher than that of the conventional reverse time migration imaging. Moreover, the artifacts at 2.16 km are suppressed. And due to the large dip angles of the two wings of the salt dome, it is usually very difficult to obtain seismic reflection signals from the two wings of the salt dome. Applying the staining algorithm can obtain an imaging result of the two wings of the salt dome with more focused energy than the conventional reverse time migration, as shown by the arrow in Figure 23 .

[0165] The reverse time migration time consumption and computational memory consumption of five velocity models are tested. The five velocity models are as shown in Figure 5 , 25, 12, 20, 26. Compared with the conventional reverse time migration, the staining algorithm needs to perform two synchronous forward wavefield propagations, one reverse propagation of receiver data, and two imaging conditions, which makes the staining algorithm require more computer memory and computational time. However, although the virtual receiver staining algorithm also needs to perform two forward wavefield propagations, since the wavefield at the boundary of the stained area is saved and only one forward wavefield propagation is required, compared with the conventional reverse time migration, the virtual receiver staining algorithm consumes the same computer memory as the conventional reverse time migration, which is less than the memory occupied by calculating the traditional staining algorithm. The required computational time is about 1.5 - 1.6 times that of the conventional reverse time migration imaging calculation time, which is less than 2 times that of the conventional reverse time migration imaging calculation time. The results are shown in Table 1.

Claims

1. A method for imaging the target area of viscous medium based on virtual geophone coloring algorithm, Characterized in that: The method for imaging the target area of viscous medium based on virtual geophone coloring algorithm includes the following steps: Step 1: Set the coloring position, and set N virtual geophones at the selected coloring position to receive the wave field values propagated to the coloring position; Step 2: Set the source, geophone, and seismic record information; Step 3: Use the fractional Laplacian operator with constant Q to decouple the visco-acoustic wave equation and boundary conditions to obtain the source seismic wave field values at time step T, and save the wave field values; Fractional Laplacian operator with constant Q decoupled visco-acoustic wave equation: In the equation: c is the velocity, t is the time, p is the viscoacoustic wave field; Q is the quality factor; c 0 is the phase velocity defined at the reference frequency ω 0 ; s is the source; Step 4: Save the true wave field calculated by decoupling the visco-acoustic wave equation with the constant Q of the fractional Laplacian operator at time step t Step 5: Calculate the real wave field at time step t + 1; Step 6: Determine whether the time loop reaches the maximum time step T. If not, repeat steps 3 to 5; if yes, end the time loop; Step 7: Assign the partial wave field values of the real wave field saved that reach the dyeing area to the dyeing wave field; Step VIII. Calculate the true wave field at time step j + 1 according to the virtual detection point coloring algorithm equation, and save the colored wave field Generalized coloring algorithm equation for viscous medium: In the equation and represent the original wavefield and the true source, and represent the colored wavefield and the wavefield values received by the virtual geophone; Step 9: Determine whether the time loop reaches the maximum time step T. If not, repeat steps 7 to 8; if yes, end the time loop; Step 10: Load the geophone wave field data at time step T + 1; Step 11: Use equation (1) wave equation and boundary conditions to solve the geophone seismic wave field values at time step T + 1; Step 12: Calculate the geophone seismic wave field values at time step T - t - 1, and determine whether it reaches time step 0; if not, repeat steps 7 to 9; if yes, end the loop; Step 13: Use the imaging condition to process the colored wave fields of all shots, the geophone back-propagated wave field information, and denoise to obtain the underground structure imaging result; The imaging condition is as follows: where R is the reverse wave field at the detection point, is the colored wave field.

2. A method for imaging the target area of viscous medium based on virtual geophone coloring algorithm according to claim 1, Characterized in that: The finite difference method is used to solve the time derivative on the left side of the equal sign, and the pseudo-spectral method is used to obtain the spatial operator on the right side of the equal sign; therefore, the time accuracy is second-order accuracy, and the spatial accuracy is spectral accuracy. Its calculation equation is: where F and F- 1 are one-dimensional forward and inverse Fourier transforms respectively, and k is the discrete wave number.

3. A method for imaging the target area of viscous medium based on virtual geophone coloring algorithm according to claim 1, Characterized in that: N virtual geophones are set at the coloring position in step 1, where N is greater than 80, and the number of coloring layers is not less than 5.

4. A method for imaging the target area of viscous medium based on virtual geophone coloring algorithm according to claim 1, Characterized in that: The setting of the coloring area conforms to the Kirchhoff integral equation, When the real wave field propagates to the coloring area D, the Kirchhoff integral equation gives the solution of equation (2); where r‵ is the source position vector in space, r 0 is the vector from the origin of coordinates to any point on the integration surface D, r is the coordinate vector of any point outside the space, t 0 is the time point for recording the wave field, D = 1 represents the integration surface, which is also the boundary of the space, n is the normal vector of the boundary D, and G is the Green's function of the viscoacoustic equation; Where: is the real wave field obtained by the convolution of the source and the Green's function, then Equation (7) indicates that the colored wave field obtained by the above method is the integral of the true wave field on the boundary D; that is, when the true wave field on the boundary D is known and its spatial derivative the wave field at spatial positions on the side different from the source can be obtained, which conforms to Huygens' principle; The colored wave field of the virtual geophone in the viscous medium conforms to the above conclusion, and the construction of the colored wave field can be realized by conforming to this conclusion: In the equation, is the true wave field, is the colored wave field, c is the velocity, t is the time, and Q is the quality factor; c 0 is the phase velocity defined at the reference frequency ω 0 ; s is the source, γ = arctan(1 / Q) / π c = c 0 cos(πγ / 2)D = 1 represents the colored position.

5. A method for imaging the target area of viscous medium based on virtual geophone coloring algorithm according to claim 1, Characterized in that: The fractional-order viscoacoustic equation decouples the phase dispersion term and the amplitude attenuation term, and can change the sign of the amplitude attenuation term to achieve compensation during the wavefield propagation; during the compensation process, the compensation operator increases exponentially with the increase of wavenumber, and the high-frequency components are severely amplified, resulting in numerical instability; to suppress the numerical instability phenomenon, by referring to the idea of inverse Q filtering, an adaptive stable propagation operator is used to stabilize the wavefield, and its calculation equation is: Among them, is the amplitude decay operator, and the remaining part is the stabilization operator, denoted as: Since the attenuation of seismic waves occurs during the entire wavefield propagation period, compensation should also be carried out within each wavefield continuation step. For the first continuation moment, i.e., when l = 1, the stabilization coefficient is: For the l-th continuation moment, l ≥ 2, the stabilization coefficient is: σ in the equation 2 is the introduced stabilization factor, and its value is 2.5×10 -3 ; The stabilization compensation operators in equations (11) and (12) are applied to the compensation wavefield update at each moment. Based on the constant Q model, the compensation form of the fractional-order Laplacian viscoacoustic equation is: In the equation, c is the velocity, t is the time, p is the visco - acoustic wave field; Q is the quality factor; c 0 is the phase velocity defined at the reference frequency ω 0 ; s is the seismic source.

Images