A Viscoelastic Wave Reverse Time Migration Imaging Method Based on Finite Difference Algorithm

Abstract

This invention discloses a viscoelastic wave reverse time migration imaging method based on the finite difference algorithm. The method includes: S1, deriving a viscoelastic wave attenuation equation suitable for the finite difference algorithm by introducing a weighting function approximation and combining the complex modulus of the constant Q model; S2, establishing a viscoelastic wave compensation equation by introducing a compensation term to compensate for the energy attenuation of seismic waves propagating in a viscoelastic medium, and combining the wavefield compensation term with low-pass filtering to achieve wavefield stability compensation; S3, performing P-wave and S-wave field separation on the source wavefield and the detector wavefield respectively; S4, applying the source normalized cross-correlation imaging condition to the separated wavefield to obtain attenuation-compensated reverse time migration imaging results; S5, setting up a seismic wave excitation source and a signal receiving device for wavefield data acquisition. The method of this invention can effectively improve the accuracy and resolution of seismic imaging in viscous media.

Description

Technical Field

[0001] This invention relates to the field of seismic wave migration imaging technology, and in particular to a viscoelastic wave reverse time migration imaging method based on a finite difference algorithm. Background Technology

[0002] Reverse-time migration (RTM) is an imaging method based on wavefield extrapolation, widely used for precise imaging of complex structures (such as steeply tilted reflectors and subsalt regions). Due to the universal absorption and attenuation properties of media, seismic waves experience amplitude attenuation and velocity dispersion during propagation in subsurface media. Traditional elastic wave RTM lacks corresponding compensation, severely impacting the accurate imaging and interpretation of seismic data. Therefore, developing viscoelastic wave attenuation-compensated RTM is of great significance for high-precision seismic exploration.

[0003] In the field of exploration seismology, it is generally assumed that the seismic quality factor Q does not change with frequency, i.e., the constant Q model. Based on the constant Q assumption, methods for simulating the propagation of seismic waves in viscoelastic media are broadly classified into two categories: standard linear body models and fractional derivative models. The former mainly describes the stress-strain relationship by introducing stress, strain relaxation time, and relaxation function, and can approximate the constant Q model well when using multiple standard linear bodies. By using memory variables, this type of equation can be solved well using the finite difference method, resulting in high computational efficiency. However, the traditional standard linear body model uses stress and strain relaxation time to represent the parameter Q in the equation (a physical quantity used to measure the viscosity of the medium; the smaller the Q value, the stronger the viscosity of the medium), which is not conducive to Q-value modeling and inversion. When the Q value changes, the stress and strain relaxation time need to be recalculated, increasing the computational complexity. In addition, since amplitude attenuation and velocity dispersion in the standard linear body equation are coupled and difficult to separate, it is difficult to use the standard linear body equation to independently compensate for amplitude attenuation or phase distortion, which severely limits the application of this method in attenuation-compensated reverse time migration imaging.

[0004] Another approach, the fractional derivative model, is directly derived from the constant Q model and can effectively describe the constant Q effect. Because amplitude attenuation and velocity dispersion are decoupled in this equation, and the parameter Q is retained, facilitating the derivation of the gradient formula, it is widely used in numerical simulations of viscous medium wave fields, Q-migration imaging, and inversion. However, the modulus expression derived from the fractional derivative model contains two variable fractional Laplace operators, requiring multiple Fourier transforms for solution, resulting in low computational efficiency. This problem is even more severe in viscoelastic wave equations, significantly limiting the application of this method in viscoelastic wave attenuation-compensated reverse-time migration imaging.

[0005] Wavefield separation is an unavoidable topic in elastic wave migration imaging. Traditional wavefield separation methods mainly include direct Helmholtz decomposition, wavenumber domain separation, and Helmholtz decomposition based on the Poisson equation. Although direct Helmholtz decomposition can achieve separation of P-waves and S-waves, the amplitude and phase of the separated wavefields change, making it difficult to meet the requirements of seismic imaging. Wavenumber domain separation methods require multiple Fourier transforms, resulting in low computational efficiency and susceptibility to noise interference, thus limiting their practical application. Summary of the Invention

[0006] This invention aims to at least partially solve one of the technical problems in related technologies. To this end, the first objective of this invention is to propose a viscoelastic wave reverse time migration imaging method based on a finite difference algorithm. Building upon the original method, the divergence and curl of the wave field are first calculated. Then, the separated P-wave field is obtained by solving the Poisson equation corresponding to the divergence. Finally, the separated S-wave field is obtained by subtracting the P-wave field from the original wave field, thereby effectively improving the computational efficiency of wave field separation.

[0007] To achieve the above objectives, a first aspect of the present invention proposes a viscoelastic wave reverse time migration imaging method based on a finite difference algorithm, the method comprising:

[0008] S1, combining the complex modulus of the constant Q model, and by introducing a weight function approximation, derives the viscoelastic wave attenuation equation suitable for solving the finite difference algorithm;

[0009] S2, by introducing a compensation term, establishes a viscoelastic wave compensation equation to compensate for the energy attenuation of seismic waves propagating in a viscoelastic medium, and combines the wavefield compensation term with low-pass filtering to achieve wavefield stability compensation;

[0010] S3 separates the P-wave and S-wave fields of the source wavefield and the detector wavefield, respectively.

[0011] S4. Apply source-normalized cross-correlation imaging conditions to the separated wavefield to obtain attenuation-compensated reverse time migration imaging results.

[0012] S5, a seismic wave excitation source and signal receiving device are set up for wavefield data acquisition.

[0013] Furthermore, the viscoelastic wave reverse time migration imaging method based on the finite difference algorithm according to the above embodiments of the present invention may also have the following additional technical features:

[0014] According to an embodiment of the present invention, step S1 includes:

[0015] S11, According to the approximate constant Q model, the complex modulus M(ω) of the viscous medium is expressed as:

[0016]

[0017] in, Here, ρ is the reference modulus, v0 is the reference velocity, Q is the quality factor, ω is the angular frequency, ω0 is the reference angular frequency, and i is the imaginary unit.

[0018] S12, drawing on the complex modulus expression of the generalized standard linear volume model, defines the following complex weighting function:

[0019]

[0020] Where, τ σl and τ εl W represents the stress relaxation time and strain relaxation time, which are independent of the Q value. R and W I Let L be the real and imaginary parts of the complex weighting function, and L be the number of weights.

[0021] S13, using W(ω)-W R (ω0) approximates the part in parentheses in equation (1) to obtain:

[0022]

[0023] Among them, the stress relaxation time τ σl strain relaxation time τ εl This is obtained by solving an optimization problem;

[0024] S14, W(ω)-W R (ω0) can be rewritten in the following form:

[0025]

[0026] in,

[0027]

[0028] Therefore, the complex modulus M(ω) in equation (1) can be rewritten as:

[0029]

[0030] S15, in the frequency domain, the two-dimensional viscoelastic wave equation is expressed as:

[0031]

[0032] Among them, (σ xx ,σ zz ,σ xz ) is the stress vector, (ε) xx ,ε zz ,ε xz M is the strain vector.P0 M is the P-wave reference modulus. S0 Q is the S-wave reference modulus. P For P-wave quality factor, Q S For S-wave quality factor, Let be the first spatial partial derivative in the x-direction. Let g be the first spatial partial derivative in the z direction, and the expressions for g and h are given in equation (5);

[0033] S16, define a set of auxiliary variables as follows and

[0034]

[0035] S17, multiply both sides of equation (8) by (1-iωτ) σl And sorted out:

[0036]

[0037] S18, Substituting equation (8) into equation (7), and transforming equations (7) and (9) back to the time-space domain, we obtain the following viscoelastic wave equation:

[0038]

[0039] Where l = 1, 2, ..., L, l is the number of weights in the weighting function W(ω), and the other coefficients are:

[0040]

[0041] Equation (10) is solved directly using the staggered grid finite difference algorithm.

[0042] According to an embodiment of the present invention, step S2 includes:

[0043] S21, starting from the complex modulus, i.e., equation (1), its imaginary part -iM0Q -1 Related to amplitude decay, therefore, the complex modulus M′ with compensation effect P,S The expression for (ω) is as follows:

[0044]

[0045] S22, combining the weighted function approximation and the relation ω≈kv P0,S0 Where ω is the angular frequency, k is the wave number, and v P0,S0 As the reference velocity for the P-wave or S-wave, equation (12) is rewritten as:

[0046]

[0047] The expressions for g and h are given in equation (5);

[0048] S23, Substituting the complex modulus in equation (13) into equation (7), we derive the viscoelastic wave compensation equation in the following form:

[0049]

[0050]

[0051] in, and Represents the positive and inverse Fourier transforms, l = 1, 2, ..., L, and has

[0052]

[0053] The first-order partial derivatives in the formula are all calculated using the finite difference method.

[0054] According to an embodiment of the present invention, step S3 includes:

[0055] S31, according to Helmholtz decomposition, has

[0056] ▽ 2 w = u, u P =▽(▽·w),u S =-▽×▽×w (16)

[0057] Where u is the original elastic wave field vector, w is the auxiliary wave field vector, and u P The separated P-wave field vector, u S Let be the separated S-wave field vector, ▽ be the gradient, ▽· be the divergence, and ▽× be the curl operator;

[0058] S32, define a scalar potential function respectively. And a vector potential function ψ, as detailed below:

[0059]

[0060] S33, since the P-wave is an irrotational field and the S-wave is a divergence-free field, then:

[0061]

[0062] ▽×u=▽×u S =▽×(▽×ψ)=▽ 2 ψ(18b)

[0063] S34. Based on equations (17) and (18), first solve for the divergence of the original wave field u, and then obtain the scalar potential function by solving the Poisson equation. Finally, the scalar potential function The gradient is used to obtain the separated P-wave field. The S-wave field can be obtained by subtracting the P-wave field from the original wave field u.

[0064] According to an embodiment of the present invention, step S4 includes:

[0065] S41, the source wave field is positively extended using the viscoelastic wave attenuation equation, and the source wave field is separated into P-waves and S-waves.

[0066] S42, the detector wavefield is extended in reverse using the viscoelastic wave compensation equation, and P-wave and S-wave separation is performed on the detector wavefield; to prevent wavefield instability, the amplitude compensation term is low-pass filtered, where the low-pass filtering is combined with the inverse Fourier transform in the compensation equation for calculation; taking equation (14d) as an example, assuming χ is the low-pass filtering factor, the combined equation (14d) is rewritten as follows:

[0067]

[0068] S43. Applying source-normalized cross-correlation imaging conditions to the source wavefield and detector wavefield, we obtain the imaging results of attenuation-compensated reverse-time migration (PP) waves and PS waves. The expressions for the imaging results of PP waves and PS waves are as follows:

[0069]

[0070] Among them, I PP (x) represents the PP wave imaging result, I PS (x) represents the PS wave imaging result, where x is the spatial location coordinate and t is the wavefield time. For the P-wave field of the earthquake source, For the P-wave field at the detector end, This represents the S-wave field at the detector end.

[0071] According to an embodiment of the present invention, step S5 includes: selecting a source wavelet type and frequency that matches the actual situation according to actual needs, setting the detector to receive seismic signals, setting the receiving recording time length and sampling interval according to the model size and spatial interval, and then starting the reverse time migration calculation to finally obtain the subsurface medium migration imaging result.

[0072] This invention combines the complex modulus of the constant Q model with an approximation using a weighting function to develop a viscoelastic wave attenuation equation suitable for finite difference algorithms. By introducing a compensation term, a viscoelastic wave compensation equation is further constructed, and the wavefield compensation term is combined with low-pass filtering calculations, efficiently achieving wavefield stability compensation. After obtaining the source wavefield and detector wavefield, the P-wave and S-wave fields are separated separately. Then, by applying the source normalized cross-correlation imaging condition, attenuation-compensated reverse time migration (Q-RTM) imaging is finally obtained, effectively improving the accuracy and resolution of seismic imaging in viscous media. The advantages of this invention are:

[0073] (1) In terms of wave field simulation accuracy: The viscoelastic wave equation based on the weight function approximation proposed in this invention can better simulate the propagation law of seismic waves in viscous media and obtain calculation results similar to those of the traditional fractional viscoelastic wave equation; the constructed viscoelastic wave compensation equation can effectively compensate for the absorption and attenuation effect of the viscoelastic medium on seismic waves.

[0074] (2) In terms of the accuracy of P-wave and S-wave separation: The Helmholtz decomposition method based on the Poisson equation proposed in this invention can accurately separate the P-wave and S-wave fields, and only requires solving the Poisson equation once, which has high computational efficiency.

[0075] (3) Imaging results: The viscoelastic wave attenuation compensation reverse time migration imaging method based on weight function approximation proposed in this invention can effectively improve the resolution of reverse time migration imaging results, obtain calculation results similar to those of traditional fractional equations, and have good agreement with the reference solution.

[0076] (4) In terms of calculation method and efficiency: The traditional fractional-order viscoelastic wave equation requires multiple Fourier transforms to solve, which has low calculation efficiency. This invention develops a viscoelastic wave equation suitable for finite difference solution. In Q-RTM, its calculation efficiency is improved by about 54% compared with the traditional fractional-order equation, and it has a wider range of industrial application value.

[0077] Additional aspects and advantages of the invention will be set forth in part in the description which follows, and in part will be obvious from the description, or may be learned by practice of the invention. Attached Figure Description

[0078] Figure 1 This is a flowchart of a viscoelastic wave reverse time migration imaging method based on the finite difference algorithm according to an embodiment of the present invention;

[0079] Figure 2 A parameter diagram of a gas chimney model according to an embodiment of the present invention;

[0080] Figure 3 A 0.72s wave field snapshot of a gas chimney model according to an embodiment of the present invention;

[0081] Figure 4 A comparison diagram of longitudinal and transverse wave field separation for a chimney model according to an embodiment of the present invention;

[0082] Figure 5 This is a diagram showing the reverse-time migration result of a gas chimney model according to an embodiment of the present invention;

[0083] Figure 6 This is a comparison diagram of a gas chimney model with reverse time offset from one channel according to an embodiment of the present invention. Detailed Implementation

[0084] Embodiments of the present invention are described in detail below, examples of which are illustrated in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and intended to explain the present invention, and should not be construed as limiting the present invention.

[0085] The following description, with reference to the accompanying drawings, describes a viscoelastic wave reverse time migration imaging method based on a finite difference algorithm, according to an embodiment of the present invention.

[0086] like Figure 1 As shown, the viscoelastic wave reverse time migration imaging method based on the finite difference algorithm of this invention may include the following steps:

[0087] S1, combining the complex modulus of the constant Q model, and by introducing a weighting function approximation, derives a viscoelastic wave attenuation equation suitable for solving by the finite difference algorithm.

[0088] According to an embodiment of the present invention, step S1 includes:

[0089] S11, According to the approximate constant Q model, the complex modulus M(ω) of the viscous medium is expressed as: ：

[0090]

[0091] in, Here, ρ is the reference modulus, v0 is the reference velocity, Q is the quality factor, ω is the angular frequency, ω0 is the reference angular frequency, and i is the imaginary unit.

[0092] S12, drawing on the complex modulus expression of the generalized standard linear volume model, defines the following complex weighting function:

[0093]

[0094] Where, τ σl and τ εl W represents the stress relaxation time and strain relaxation time, which are independent of the Q value. R and WI Let L be the real and imaginary parts of the complex weighting function, and L be the number of weights. In one embodiment of the invention, L = 3.

[0095] S13, using W(ω)-W R (ω0) approximates the part in parentheses in equation (1) to obtain:

[0096]

[0097] Among them, the stress relaxation time τ σl strain relaxation time τ εl The results are obtained by solving an optimization problem. For example, when L = 3, the stress relaxation time and strain relaxation time in the complex weighted function W(ω) are shown in Table 1.

[0098] Table 1

[0099] l <![CDATA[τ σl (s)]]> <![CDATA[τ εl (s)]]> 1 0.96434144E-1 2.64139184E-1 2 1.1208990E-2 2.4713628E-2 3 1.5780938E-3 4.0529653E-3

[0100] S14, W(ω)-W R (ω0) can be rewritten in the following form:

[0101]

[0102] in,

[0103]

[0104] Therefore, the complex modulus M(ω) in equation (1) can be rewritten as:

[0105]

[0106] S15, in the frequency domain, the two-dimensional viscoelastic wave equation is expressed as:

[0107]

[0108] Among them, (σ xx ,σ zz ,σ xz ) is the stress vector, (ε) xx ,ε zz ,ε xz M is the strain vector. P0 M is the P-wave reference modulus. S0 Q is the S-wave reference modulus. P For P-wave quality factor, Q S For S-wave quality factor, Let be the first spatial partial derivative in the x-direction. Let g be the first spatial partial derivative in the z direction, and the expressions for g and h are given in equation (5);

[0109] S16, define a set of auxiliary variables as follows and

[0110]

[0111] S17, multiply both sides of equation (8) by (1-iωτ) σl And sorted out:

[0112]

[0113] S18, Substituting equation (8) into equation (7), and transforming equations (7) and (9) back to the time-space domain, we obtain the following viscoelastic wave equation:

[0114]

[0115]

[0116] Where l = 1, 2, ..., L, l is the number of weights in the weighting function W(ω), and the other coefficients are:

[0117]

[0118] Equation (10) can be solved directly using the staggered grid finite difference algorithm, which has high computational efficiency.

[0119] S2 establishes a viscoelastic wave compensation equation by introducing a compensation term to compensate for the energy attenuation of seismic waves propagating in a viscoelastic medium. The wavefield compensation term is then combined with low-pass filtering to achieve wavefield stability compensation.

[0120] According to an embodiment of the present invention, step S2 includes:

[0121] S21, because the amplitude and phase effects are coupled in the viscoelastic wave attenuation equation based on the weighting function approximation, they are difficult to separate directly. Starting from the complex modulus, i.e., equation (1), its imaginary part -iM0Q -1 Related to amplitude decay, therefore, the complex modulus M′ with compensation effect P,S The expression for (ω) is as follows:

[0122]

[0123] S22, combining the weighted function approximation and the relation ω≈kv P0,S0 Where ω is the angular frequency, k is the wave number, and v P0,S0 As the reference velocity for the P-wave or S-wave, equation (12) is rewritten as:

[0124]

[0125] The expressions for g and h are given in equation (5);

[0126] S23, Substituting the complex modulus in equation (13) into equation (7), we derive the viscoelastic wave compensation equation in the following form:

[0127]

[0128] in, and Represents the positive and inverse Fourier transforms, l = 1, 2, ..., L, and has

[0129]

[0130] The first-order partial derivatives in the formula are all calculated using the finite difference method.

[0131] S3 separates the longitudinal and transverse wave fields of the source wave field and the detector wave field, respectively.

[0132] It is important to understand that the separation of longitudinal and transverse wave fields is an essential step in the reverse-time migration of viscoelastic waves, and its separation accuracy and efficiency have a significant impact on the results of the reverse-time migration of viscoelastic waves.

[0133] According to an embodiment of the present invention, step S3 includes:

[0134] S31, according to Helmholtz decomposition, has

[0135] ▽ 2 w = u, u P =▽(▽·w),u S =-▽×▽×w (16)

[0136] Where u is the original elastic wave field vector, w is the auxiliary wave field vector, and u P The separated P-wave field vector, u S Let be the separated S-wave field vector, ▽ be the gradient, ▽· be the divergence, and ▽× be the curl operator;

[0137] S32, define a scalar potential function respectively. And a vector potential function ψ, as detailed below:

[0138]

[0139] S33, since the P-wave is an irrotational field and the S-wave is a divergence-free field, then:

[0140]

[0141] ▽×u=▽×u S =▽×(▽×ψ)=▽ 2ψ(18b)

[0142] S34. Based on equations (17) and (18), first solve for the divergence of the original wave field u, and then obtain the scalar potential function by solving the Poisson equation. Finally, the scalar potential function The gradient is used to obtain the separated P-wave field, and the S-wave field is obtained by subtracting the P-wave field from the original wave field u. This method only requires solving the Poisson equation once.

[0143] S4. Apply source-normalized cross-correlation imaging conditions to the separated wavefield to obtain attenuation-compensated reverse time migration imaging results.

[0144] It is important to understand that, due to the viscoelasticity of the medium, seismic waves are prone to amplitude attenuation and velocity dispersion when propagating underground. Conventional reverse time migration methods are prone to problems such as low resolution and inaccurate migration positions. However, the attenuation-compensated reverse time migration (Q-RTM) method compensates for the seismic wave field during wavefield extension, which can theoretically effectively eliminate the adverse effects of the viscoelasticity of the medium on seismic wave propagation and imaging.

[0145] According to an embodiment of the present invention, step S4 includes:

[0146] S41, the source wave field is positively extended using the viscoelastic wave attenuation equation, and the source wave field is separated into P-waves and S-waves.

[0147] S42, the detector wavefield is extended in reverse using the viscoelastic wave compensation equation, and P-wave and S-wave separation is performed on the detector wavefield; to prevent wavefield instability, the amplitude compensation term is low-pass filtered, where the low-pass filtering is combined with the inverse Fourier transform in the compensation equation for calculation; taking equation (14d) as an example, assuming χ is the low-pass filtering factor, the combined equation (14d) is rewritten as follows:

[0148]

[0149] S43. Applying source-normalized cross-correlation imaging conditions to the source wavefield and detector wavefield, we obtain the imaging results of attenuation-compensated reverse-time migration (PP) waves and PS waves. The expressions for the imaging results of PP waves and PS waves are as follows:

[0150]

[0151] Among them, I PP (x) represents the PP wave imaging result, I PS (x) represents the PS wave imaging result, where x is the spatial location coordinate and t is the wavefield time. For the P-wave field of the earthquake source, For the P-wave field at the detector end, S5 represents the S-wave field at the detector end. A seismic wave excitation source and signal receiving device are installed at S5 for wavefield data acquisition.

[0152] According to an embodiment of the present invention, step S5 includes: selecting a source wavelet type and frequency that matches the actual situation according to actual needs, setting the detector to receive seismic signals, setting the receiving recording time length and sampling interval according to the model size and spatial interval, and then starting the reverse time migration calculation to finally obtain the subsurface medium migration imaging result.

[0153] To make the objectives, technical solutions, and advantages of this invention clearer, based on the technical solution flow of this invention, a complex BP-gas chimney model is used as an example to test the viscoelastic wave attenuation compensation reverse time migration imaging method. The BP-gas model has a size of 398×161 computational grids, with a grid spacing of 8 meters both horizontally and vertically. The specific steps are as follows:

[0154] (1) As Figure 2 The diagram shows the parameters of the chimney model. Figure 2 'a' represents the P-wave reference velocity for the complex BP-gas chimney model. Figure 2 b is the quality factor of the complex BP-gas chimney model, and its transverse wave velocity is the longitudinal wave velocity divided by 1.6, i.e., v S0 =v P0 / 1.6, the shear wave quality factor is the longitudinal wave quality factor divided by 1.3, i.e., Q S =Q P / 1.3, the grid spacing for both the horizontal and vertical directions of the model is 8 meters. Based on these parameters, the inverse time migration parameters were set. Specifically: the time sampling interval is 0.6 milliseconds, the recording length is 2.4 seconds, with a total of 4000 sampling points, and the reference angular frequency is ω0 = 40π rad / s. The seismic source uses a Ricker wavelet with a dominant frequency of 25 Hz, placed at a depth of 8 meters to excite the seismic signal. Horizontally, a total of 57 shots are used, with a shot spacing of 56 meters (7 grid spacings), evenly distributed among the grid points. 398 geophones are placed on the ground to receive the seismic signal, with a spacing of 8 meters between the geophones and a low-pass filter cutoff frequency of 135 Hz.

[0155] (2) Absorbing Boundary Conditions. Due to the limited simulation area, absorbing boundaries need to be added when seismic waves propagate in the model. In this example, a 25-layer mixed absorbing boundary condition is used to suppress reflections generated by the model boundaries.

[0156] (3) After setting the relevant parameters, according to the Q-compensated reverse time migration imaging process in step S4, firstly calculate the source forward propagation wave field according to the viscoelastic wave attenuation equation derived in step S1 of the technical solution, and perform P-wave and S-wave separation on the source wave field according to step S3 of the technical solution; then calculate the detector reverse propagation wave field using the viscoelastic wave compensation equation derived in step S2 of the technical solution and the seismic record, and perform P-wave and S-wave separation on the detector wave field; finally, apply the source normalized cross-correlation imaging condition to the forward propagation wave field and the reverse propagation wave field according to step S4 of the technical solution.

[0157] (4) During the reverse time-shift imaging process, Figure 3 This shows a snapshot of the forward-modeled wave field of the 29th shot, in which... Figure 3 a and Figure 3 d represents the value of v calculated from the fractional-order viscoelastic wave equation. x and v z result, Figure 3 b and Figure 3 e represents the result calculated using the new equation. Figure 3 c and Figure 3 f represents the difference between the calculation results of the traditional fractional equation and the new equation, which shows that the difference between the two is small.

[0158] (5) The effect of longitudinal and transverse wave separation greatly affects the accuracy of viscoelastic wave attenuation compensation reverse time migration imaging results. Figure 4 The results of P-wave and S-wave separation of the source wavefield of the 29th shot are shown, among which... Figure 4 a and Figure 4 d represents the separated longitudinal wave component. Figure 4 b and Figure 4 e represents the separated transverse wave component. Figure 4 c and Figure 4 f represents the difference between the original wavefield and the separated P-wave and S-wave wavefields (see original wavefield). Figure 3 b and Figure 3 e). By Figure 3 c and Figure 3 It can be seen that the longitudinal and transverse wave separation method proposed in this invention has good separation accuracy.

[0159] (6) Figure 5 The figures show the RTM results for the gas chimney model obtained using different methods. Figure 5 a- Figure 5 b represents the reference solution calculated using unattenuated recording + elastic wave RTM. Compared to the reference solution, the elastic wave equation attenuation recording method has no compensation, and the calculated results ( Figure 5 c- Figure 5d) The energy is weak, the resolution is low, and some wavefields cannot be correctly repositioned. In contrast, the migration profiles calculated using several attenuation-compensated Q-RTM methods have clearer interfaces, stronger energy, significantly improved resolution, and better consistency with the reference solution.

[0160] (7) To specifically compare waveform differences, from Figure 5 Several single tracks were sampled at x=400 and 800m for comparison (see...) Figure 6 This allows for a clearer observation of the aforementioned phenomena. Therefore, the viscoelastic wave attenuation compensation reverse time migration Q-RTM method proposed in this invention can obtain calculation results similar to the reference solution and the traditional fractional-order viscoelastic wave equation Q-RTM. However, the fractional-order viscoelastic wave equation takes 10.177 hours, while the new viscoelastic wave equation Q-RTM provided by this invention takes only 4.704 hours (see Table 2 for details). Therefore, the new method described in this invention improves computational efficiency by approximately 54% compared to the traditional fractional-order viscoelastic wave equation Q-RTM, and has broader application prospects.

[0161] Table 2

[0162]

[0163] In the description of this specification, references to terms such as "one embodiment," "some embodiments," "example," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0164] Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Thus, a feature defined as "first" or "second" may explicitly or implicitly include at least one of that feature. In the description of this invention, "a plurality of" means at least two, such as two, three, etc., unless otherwise explicitly specified.

[0165] In this invention, unless otherwise explicitly specified and limited, the terms "installation," "connection," "linking," and "fixing," etc., should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral part; they can refer to a mechanical connection or an electrical connection; they can refer to a direct connection or an indirect connection through an intermediate medium; they can refer to the internal communication of two components or the interaction between two components, unless otherwise explicitly limited. Those skilled in the art can understand the specific meaning of the above terms in this invention according to the specific circumstances.

[0166] Although embodiments of the present invention have been shown and described above, it is understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments within the scope of the present invention.

Claims

1. A viscoelastic wave reverse-time migration imaging method based on a finite-difference algorithm, characterized in that, The method comprises: S1, in combination with the complex modulus of the Q model, an attenuation equation of viscoelastic wave suitable for solving by a finite difference algorithm is derived by introducing a weight function approximation; wherein, step S1 comprises: S11, complex modulus of the viscous medium according to the approximate constant Q model is represented as: （1） wherein is the reference modulus, p is the density, v0 is the reference velocity, Q is the quality factor, ω is the angular frequency, ω0 is the reference angular frequency, and i is the imaginary unit. S12, referring to the complex modulus expression of the generalized standard linear solid model, a complex weight function is defined as follows: （2） where τ σl and τ εl are the stress and strain relaxation times, respectively, independent of Q, W R and W I are the real and imaginary parts of a complex weighting function, and L is the number of weights. S13, using The time-space domain viscoelastic wave equation is obtained by approximating the bracketed part of equation (1), and is solved using the staggered-grid finite difference algorithm. S2, a viscoelastic wave compensation equation is established by introducing a compensation term, which is used to compensate the energy attenuation of seismic waves in viscoelastic media during propagation, and the wave field compensation term is combined with low-pass filtering for calculation to realize stable compensation of the wave field; wherein, step S2 comprises: S21, starting from the complex modulus, i.e. equation (1), its imaginary part is related to the amplitude attenuation, therefore, the complex modulus with compensation effect is expressed as follows: (12) S22, combining the weighted function approximation and the relationship where, ω is the angular frequency, k is the wave number, is the reference velocity of P or S wave, equation (12) is rewritten as: (13)； S23, the complex modulus in equation (13) is brought into the two-dimensional viscoelastic wave equation, and a derived viscoelastic wave compensation equation is obtained; S3, longitudinal and transverse wave field separation is respectively performed on the source wave field and the receiver wave field, and the wave field separation is realized based on the Helmholtz decomposition method of the Poisson equation; S4, the source normalization cross-correlation imaging condition is applied to the separated wave field to obtain the attenuation compensation reverse time migration imaging result; wherein, step S4 comprises: forward propagating the source wave field using the viscoelastic wave attenuation equation; backward propagating the receiver wave field using the viscoelastic wave compensation equation, wherein the low-pass filtering is combined with the inverse Fourier transform in the compensation equation for calculation; S5, a seismic wave excitation source and a signal receiving device are set for wave field data acquisition.

2. The viscoelastic wave reverse time migration imaging method based on the finite difference algorithm of claim 1, wherein, Step S1 further comprises: S13, utilizing Approximating the bracketed part of equation (1), we get: （3） where the stress relaxation time τ σl , the strain relaxation time τ εl is obtained by solving an optimization problem; S14, to is rewritten as follows: （4） Wherein, （5） Thus, the complex modulus in equation (1) is rewritten as: （6）； S15, in the frequency domain, the two-dimensional viscoelastic wave equation is expressed as: (7a) (7b) (7c) wherein is a stress vector, is a strain vector, M P0 is a P-wave reference modulus, M S0 is a S-wave reference modulus, Q P is a P-wave quality factor, Q S is a S-wave quality factor, is a first order spatial derivative in the x-direction, is a first order spatial derivative in the z-direction, the expressions for g and h are given in equation (5); S16, define a set of auxiliary variables as follows , and : (8a) (8b) (8c)； S17, multiply both sides of equation (8) by and rearranged to be: (9a) (9b) (9c)； S18, equation (8) is brought into equation (7), and equation (7) and (9) are transformed back to the time-space domain to obtain the following viscoelastic wave equation: (10a) (10b) (10c) (10d) (10e) (10f) (10g) wherein, l = 1, 2, …L, l is the number of weighting functions the number of weighting in the equation, and other coefficients are: (11)； Equation (10) is directly solved using the staggered grid finite difference algorithm.

3. The viscoelastic wave reverse time migration imaging method based on the finite difference algorithm of claim 2, wherein, Step S2 further comprises: S23, the complex modulus in equation (13) is brought into equation (7) to derive the following form of the viscoelastic wave compensation equation: (14a) (14b) (14c) (14d) (14e) (14f) (14g) wherein and represent the forward and inverse Fourier transform, l = 1, 2, …, L and have (15) In the formula, the first-order partial derivatives are calculated using the finite difference method.

4. The viscoelastic wave reverse time migration imaging method based on the finite difference algorithm of claim 3, wherein, Step S3 comprises: S31, according to the Helmholtz decomposition, there is (16) where u is the original elastic wavefield vector, w is the auxiliary wavefield vector, is the separated P-wavefield vector, is the separated S-wavefield vector, is the gradient, is the divergence and is the curl operator; S32, define a scalar potential function and a vector potential function respectively as follows: (17) S33, since the P wave is a non-rotational field and the S wave is a non-divergence field, there is: (18a) (18b) S34, according to equations (17) and (18), first solve the divergence of the original wave field u, and then solve the Poisson equation to get the scalar potential function , and finally take the gradient of the scalar potential function to get the separated P wave field. The S wave field can be obtained by subtracting the P wave field from the original wave field u.

5. The viscoelastic wave reverse time migration imaging method based on the finite difference algorithm of claim 4, wherein, Step S4 comprises: S41, the source wave field is forward propagated using the viscoelastic wave attenuation equation, and longitudinal and transverse wave separation is performed on the source wave field; S42, the wave field of the detector is back-propagated using the viscoelastic wave compensation equation, and the wave field of the detector is separated into P and S waves; in order to prevent wave field compensation instability, the amplitude compensation term is low-pass filtered, wherein the low-pass filtering is combined with the inverse Fourier transform in the compensation equation for calculation; taking equation (14d) as an example, assuming χ is a low-pass filter factor, then the combined equation (14d) is rewritten in the following form: (19) S43, the source normalization cross-correlation imaging condition is applied to the source wave field and the receiver wave field to obtain the imaging results of the attenuation compensation reverse time migration PP wave and PS wave, wherein, the expressions of the imaging results of the PP wave and the PS wave are: (20a) (20b) wherein, is the PP wave imaging result, is the PS wave imaging result, x is the spatial position coordinate, t is the wavefield time, is the source P wavefield, is the receiver end P wavefield, is the receiver end S wavefield.

6. The viscoelastic wave reverse time migration imaging method based on finite difference algorithm of claim 1, wherein, Step S5 comprises: According to actual needs, the source wave type and frequency that are consistent with the actual situation are selected, the receiver receives the seismic signal, the time length of the received record and the sampling interval are set according to the model size and the spatial interval, then the reverse time migration calculation is started, and finally the underground medium migration imaging result is obtained.

Images