Efficient viscoacoustic wave attenuation compensation reverse time migration method based on weight function approximation

Abstract

This invention discloses an efficient viscosonic wave attenuation compensation reverse-time migration method based on a weighted function approximation. The method includes: S1, introducing a weighted function approximation to derive a viscosonic wave attenuation equation applicable to finite difference solutions; S2, based on the viscosonic wave attenuation equation, establishing a viscosonic wave compensation equation to compensate for energy attenuation during seismic wave propagation in a viscous medium; S3, implementing an attenuation-compensated reverse-time migration imaging process, wherein the wavefield compensation term is combined with low-pass filtering calculations to achieve wavefield stability compensation; S4, configuring the source and receiver settings to perform reverse-time migration imaging. The method of this invention efficiently achieves wavefield stability compensation and attenuation-compensated reverse-time migration imaging, effectively improving the accuracy and resolution of seismic imaging in viscosonic media.

Description

Technical Field

[0001] This invention relates to the field of seismic wave migration imaging technology, and in particular to a highly efficient method for retrograde migration with viscosonic wave attenuation compensation based on weighting function approximation. Background Technology

[0002] Viscoelasticity is a fundamental property of media. Due to the viscoelasticity of the medium, seismic waves typically experience amplitude reduction and phase distortion during propagation in subsurface media, leading to changes in the kinematic and dynamic characteristics of seismic waves and severely impacting the processing and interpretation of seismic data. Therefore, accurately and efficiently simulating the propagation patterns of seismic waves in viscous media and performing targeted compensation is of significant research importance for seismic wave migration imaging.

[0003] The seismic quality factor Q is a commonly used parameter to describe the viscosity of a medium; the smaller the Q value, the stronger the viscosity of the medium, and the stronger its absorption and attenuation effect on seismic waves. In the field of exploration seismology, it is usually assumed that the seismic quality factor Q does not change with frequency, i.e., the constant Q model. Based on the constant Q assumption, seismic wave attenuation simulation methods can be roughly divided into two categories: traditional mechanical models and fractional derivative models. The former simulates the approximate constant Q effect by simulating multiple sets of springs (elastic) and dampers (viscous) in series or parallel, with the generalized standard linear body model being the most representative. This method introduces multiple sets of memory variables, allowing direct solution of the equations using the finite difference method, resulting in high computational efficiency. However, this equation no longer contains the quality factor Q, but instead uses multiple sets of stress and strain relaxation times. When the Q value changes, the relaxation time needs to be recalculated; since the equation does not contain Q parameters, the sensitive kernel for Q cannot be directly derived, which is not conducive to Q parameter inversion modeling. In addition, the amplitude attenuation and phase distortion effects in this equation are coupled, making it difficult to simultaneously achieve accurate compensation for amplitude attenuation and phase information. Although some scholars have proposed polynomial correction strategies for this type of equation, the process is cumbersome and the effect is not good, so it is rarely used in actual attenuation-compensated reverse time migration (Q-RTM).

[0004] In contrast, the fractional derivative model, directly derived from the constant Q model, can be directly characterized using parameters such as reference velocity, reference angular frequency, and quality factor. Its equation form is relatively simple and it is widely used in numerical simulations of wavefields in viscous acoustic media, Q-migration imaging, and inversion. However, the fractional viscous acoustic equation based on the constant Q model contains two variable fractional Laplace operators, requiring multiple Fourier transforms to solve, resulting in low computational efficiency and severely limiting its application in seismic data processing.

[0005] To address the aforementioned issues, some scholars have combined the generalized standard linear body model and fractional-order constant Q theory, constructing weighting functions to approximate the dispersion and attenuation-related components of the constant Q complex modulus. This resulted in the derivation of a set of viscous acoustic wave equations solvable using finite-difference algorithms. Unlike traditional generalized standard linear body equations, the new equations retain the parameter Q, facilitating Q-inversion. Furthermore, the stress and strain relaxation times introduced in the equations are independent of the Q value; therefore, when the Q value changes, there is no need to repeatedly calculate the relaxation time, making it better suited for numerical simulations of complex media with Q anisotropy.

[0006] However, similar to the generalized standard linear body equation, the amplitude and phase changes in the new equation are coupled and difficult to separate, which severely limits the application of this method in Q-RTM imaging. Summary of the Invention

[0007] 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 an efficient viscous acoustic wave attenuation compensation reverse-time migration method based on a weighted function approximation. This method efficiently achieves wavefield stability compensation and attenuation compensation reverse-time migration imaging, effectively improving the accuracy and resolution of seismic imaging in viscous acoustic media.

[0008] To achieve the above objectives, a first aspect of the present invention proposes an efficient method for compensating for viscoacoustic wave attenuation based on a weighting function approximation, the method comprising:

[0009] S1. By introducing a weighting function approximation, the viscous acoustic wave attenuation equation applicable to finite difference solutions is derived.

[0010] S2, Based on the aforementioned viscous acoustic wave attenuation equation, establish a viscous acoustic wave compensation equation to compensate for the energy attenuation of seismic waves propagating in a viscous medium.

[0011] S3, implement the attenuation-compensated reverse time migration imaging process, in which the wavefield compensation term is combined with the low-pass filter calculation to achieve wavefield stabilization compensation;

[0012] S4, configure source and receiver settings to perform reverse time migration imaging.

[0013] Furthermore, the efficient viscosonic wave attenuation compensation reverse-time migration method based on weighted function approximation 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, ρ is the reference modulus, v0 is the reference velocity, Q is the quality factor, ω and ω0 are the angular frequency and the reference angular frequency, respectively, 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 These are the stress relaxation time and strain relaxation time, which are independent of the Q value, respectively; W R and W I These are the real and imaginary parts of the complex weighting function, respectively, and L is 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 viscous acoustic wave equation is expressed as:

[0031]

[0032] in, For the frequency domain wave field, ρ is the density, and ω is the angular frequency;

[0033] S16, Substitute the complex modulus in equation (6) into equation (7), and define the following auxiliary variables.

[0034] Equations (7) and (8) can be rewritten in the following form:

[0035]

[0036]

[0037] S17, Transforming equation (9) back to the time domain, we obtain the viscous acoustic wave attenuation equation:

[0038]

[0039] in,

[0040]

[0041] Equation (10) can be solved directly using finite difference.

[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 Since it is related to amplitude decay, the complex modulus M′(ω) with compensation effect is expressed as follows:

[0044]

[0045] S22, combined with the compensation complex modulus M′(ω), the frequency wavenumber domain compensation viscous acoustic wave equation can be written in the following form:

[0046]

[0047] in Let ρ be the frequency domain wave field, ρ be the density, and k be the wave number.

[0048] S23, Substitute the complex modulus in equation (12) into equation (13):

[0049]

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

[0051] S24, combined with the approximate expression ω≈kv0, where k is the wave number and v0 is the reference velocity, equation (14) is rewritten as:

[0052]

[0053] S25, transforming equation (15) back to the frequency-space domain, we get:

[0054]

[0055] S26, define the following auxiliary variables:

[0056]

[0057] S27, combining auxiliary variables and transforming equation (17) back to the time-space domain, we get:

[0058]

[0059] Among them, F -1 For the inverse Fourier transform, the other coefficients are as follows:

[0060]

[0061] The first and second time partial derivatives in the formula are calculated using the finite difference method.

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

[0063] S31, using the viscoacoustic attenuation equation to positively extend the source wave field;

[0064] S32 uses the viscous acoustic wave compensation equation to reverse extend the detector wave field. In order to prevent the wave field compensation from being unstable, the amplitude compensation term is subjected to low-pass filtering. The low-pass filtering is combined with the inverse Fourier transform in the compensation equation for calculation.

[0065] S33, applying source-normalized cross-correlation imaging conditions to the source wavefield and detector wavefield, obtaining Q-RTM imaging results:

[0066]

[0067] Where I(x,z) represents the imaging result, x and z are the spatial coordinates, t is the wavefield time, T is the recording length, and σ s (x,z,t) represents the source wave field, σ r (x,z,t) represents the detector wave field, and ε is a decimal to avoid dividing the denominator by zero.

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

[0069] Based on actual needs, select the source wavelet type and frequency that match the actual situation, set the detector to receive the seismic signal, set the receiving recording time length and sampling interval according to the model size and spatial interval, and then start the reverse time migration calculation to finally obtain the subsurface medium migration imaging results.

[0070] This invention combines the complex modulus of the constant Q model with an approximation using a weighting function to develop a viscous acoustic wave attenuation equation suitable for solving using the finite difference algorithm. By introducing a compensation term, a viscous acoustic wave compensation equation is further constructed, and the wavefield compensation term is combined with low-pass filtering calculations, efficiently achieving wavefield stability compensation and attenuation compensation in reverse time migration (Q-RTM) imaging, effectively improving the accuracy and resolution of seismic imaging in viscous acoustic media. The advantages of this invention are as follows:

[0071] (1) Imaging results: The viscous acoustic wave equation based on the weight function approximation proposed in this invention can effectively compensate for the absorption and attenuation effect of viscous media on seismic waves, significantly improve the resolution of reverse time migration imaging results, and obtain calculation results similar to those of traditional fractional equations.

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

[0073] 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

[0074] Figure 1 This is a flowchart of an efficient viscous acoustic wave attenuation compensation reverse time migration method based on weight function approximation according to an embodiment of the present invention;

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

[0076] Figure 3 The wave field diagrams for the forward and reverse propagation of the 29th shot in a gas chimney model according to an embodiment of the present invention are shown.

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

[0078] Figure 5 The results and amplitude spectrum of a gas chimney model according to an embodiment of the present invention are shown. Detailed Implementation

[0079] 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.

[0080] The following description, with reference to the accompanying drawings, describes an efficient viscoacoustic attenuation compensation method for reverse time migration based on weighted function approximation proposed in an embodiment of the present invention.

[0081] like Figure 1 As shown, the efficient viscoacoustic wave attenuation compensation reverse-time migration method based on weighted function approximation of this invention may include the following steps:

[0082] S1. By introducing a weighting function approximation, we derive a viscous acoustic attenuation equation applicable to finite difference solutions.

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

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

[0085]

[0086] in, ρ is the reference modulus, v0 is the reference velocity, Q is the quality factor, ω and ω0 are the angular frequency and the reference angular frequency, respectively, and i is the imaginary unit;

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

[0088]

[0089] Where, τ σl and τ εl These are the stress relaxation time and strain relaxation time, which are independent of the Q value, respectively; W R and W I These are the real and imaginary parts of the complex weighting function, respectively, and L is the number of weights;

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

[0091]

[0092] Among them, the stress relaxation time τ σl strain relaxation time τ εl The results were obtained by solving an optimization problem. For example, the stress relaxation time and strain relaxation time corresponding to L=3 and L=5 are shown in Table 1 and Table 2.

[0093] Table 1

[0094] 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

[0095] Table 2

[0096] l <![CDATA[τ σl (s)]]> <![CDATA[τ εl (s)]]> 1 1.8230838E-1 4.5748839E-1 2 3.2947348E-2 6.3276617E-2 3 8.4325390E-3 15.4145588E-3 4 2.3560480E-3 4.2784094E-3 5 5.1033826E-4 12.3424456E-4

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

[0098]

[0099] in,

[0100]

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

[0102]

[0103] S15, in the frequency domain, the two-dimensional viscous acoustic wave equation is expressed as:

[0104]

[0105] in, For the frequency domain wave field, ρ is the density, and ω is the angular frequency;

[0106] S16, Substitute the complex modulus in equation (6) into equation (7), and define the following auxiliary variables.

[0107]

[0108] Equations (7) and (8) can be rewritten in the following form:

[0109]

[0110] S17, Transforming equation (9) back to the time domain, we obtain the viscous acoustic wave attenuation equation:

[0111]

[0112] in,

[0113]

[0114] Equation (10) can be solved directly using finite difference, which is computationally efficient.

[0115] S2. Based on the viscous acoustic wave attenuation equation, a viscous acoustic wave compensation equation is established to compensate for the energy attenuation of seismic waves propagating in viscous media.

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

[0117] S21, due to the coupling between amplitude and phase effects in the viscosity wave attenuation equation based on the weighting function approximation, it is difficult to separate them directly. Starting from the complex modulus, i.e., equation (1), its imaginary part -iM0Q -1 Since it is related to amplitude decay, the complex modulus M′(ω) with compensation effect is expressed as follows:

[0118]

[0119] S22, combined with the compensation complex modulus M′(ω), the frequency wavenumber domain compensation viscous acoustic wave equation can be written in the following form:

[0120]

[0121] in Let ρ be the frequency domain wave field, ρ be the density, and k be the wave number.

[0122] S23, Substitute the complex modulus in equation (12) into equation (13):

[0123]

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

[0125] S24, combined with the approximate expression ω≈kv0, where k is the wave number and v0 is the reference velocity, equation (14) is rewritten as:

[0126]

[0127] S25, transforming equation (15) back to the frequency-space domain, we get:

[0128]

[0129] S26, define the following auxiliary variables:

[0130]

[0131] S27, combining auxiliary variables and transforming equation (17) back to the time-space domain, we get:

[0132]

[0133] Among them, F -1 For the inverse Fourier transform, the other coefficients are as follows:

[0134]

[0135] The first and second time partial derivatives in the formula are calculated using the finite difference method.

[0136] S3 implements the attenuation-compensated reverse time migration imaging process, in which the wavefield compensation term is combined with the low-pass filter calculation to achieve wavefield stabilization compensation.

[0137] Specifically, due to the viscoelasticity of the medium, seismic waves are prone to amplitude attenuation and velocity dispersion when propagating underground, resulting in low resolution and inaccurate migration positions in conventional reverse time migration imaging profiles. The attenuation-compensated reverse time migration (Q-RTM) method, by appropriately compensating for the seismic wavefield during wavefield extension, can theoretically effectively eliminate the adverse effects of the viscoelasticity of the medium on seismic wave propagation and imaging.

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

[0139] S31, using the viscoacoustic attenuation equation to positively extend the source wave field;

[0140] S32 uses the viscous acoustic wave compensation equation to reverse extend the detector wave field. In order to prevent the wave field compensation from being unstable, the amplitude compensation term is subjected to low-pass filtering. The low-pass filtering is combined with the inverse Fourier transform in the compensation equation for calculation.

[0141] S33, applying source-normalized cross-correlation imaging conditions to the source wavefield and detector wavefield, obtaining Q-RTM imaging results:

[0142]

[0143] Where I(x,z) represents the imaging result, x and z are the spatial coordinates, t is the wavefield time, T is the recording length, and σ s (x,z,t) represents the source wave field, σ r (x,z,t) represents the detector wave field, and ε is a decimal to avoid dividing the denominator by zero.

[0144] S4, configure source and receiver settings to perform reverse time migration imaging.

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

[0146] Based on actual needs, select the source wavelet type and frequency that match the actual situation, set the detector to receive the seismic signal, set the receiving recording time length and sampling interval according to the model size and spatial interval, and then start the reverse time migration calculation to finally obtain the subsurface medium migration imaging results.

[0147] 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 viscosonic 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 13 meters both horizontally and vertically. The specific steps are as follows:

[0148] (1) Based on the reference velocity and quality factor of the complex BP-gas chimney model ( Figure 2 The reverse time migration parameters were set. Specifically, the time sampling interval was 1 millisecond, the recording length was 3 seconds, and the reference angular frequency was ω0 = 40π rad / s. The seismic source used a Ricker wavelet with a dominant frequency of 25 Hz, placed at a depth of 13 meters to excite the seismic signal. A total of 57 shots were used horizontally, with a shot spacing of 91 meters, evenly distributed among the grid points. 398 geophones were placed on the ground to receive the seismic signal, with a geophone spacing of 13 meters and a low-pass filter cutoff frequency of 120 Hz.

[0149] (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, 20 layers of mixed absorbing boundary conditions are used to suppress reflections generated by the model boundaries.

[0150] (3) After setting the relevant parameters, according to the process of attenuation compensation reverse time migration imaging in step S3, first calculate the source forward propagation wave field according to the viscosonic wave attenuation equation derived in step S1 of the technical solution, then use the viscosonic wave compensation equation derived in step S2 of the technical solution and the seismic record to calculate the detector reverse propagation wave field, and finally apply the source normalization cross-correlation imaging condition to the forward propagation wave field and the reverse propagation wave field. Figure 3 The results show a comparison of the source and detector wavefields of the 29th shot from the gas chimney model. Figure 3 'a' represents the result of the fractional equation calculation. Figure 3 b and Figure 3 c represents the calculated results for the parameters L = 3 and 5 in the new equation. Figure 3 d and Figure 3 e are respectively Figure 3 a and Figure 3 b, Figure 3 a and Figure 3 As shown in the figure, the difference in c indicates that the calculation results of the viscosonic wave equation based on the weight function approximation are in good agreement with the source wave field and detector wave field calculated by the traditional fractional-order viscosonic wave equation. Figure 4 The figures show the RTM results for the gas chimney model obtained using different methods. Figure 4 'a' represents the reference solution for calculating the acoustic wave equation using acoustic recording. Compared to the reference solution, the acoustic wave equation attenuation recording method does not offer compensation, and the calculated result ( Figure 4b) The energy is weaker and the resolution is lower; in contrast, the migration profiles calculated using several attenuation-compensated Q-RTM methods have clearer interfaces, stronger energy, and significantly improved resolution, showing good consistency with the reference solution. To specifically compare the waveform differences, from... Figure 4 Several single tracks were selected for comparison (see...) Figure 5 This allows for a clearer observation of the aforementioned phenomena. Therefore, the viscosonic wave Q-RTM method based on weighted function approximation proposed in this invention can obtain computational results similar to the reference solution and the traditional fractional-order viscosonic wave equation Q-RTM. However, the fractional-order viscosonic wave equation takes 12091.75 seconds, the method using a weighted function of L=3 takes 3724.70 seconds, and the method using a weighted function of L=5 takes 4053.81 seconds (see Table 3 for details). Therefore, the new method described in this invention improves computational efficiency by approximately 67% compared to the traditional fractional-order viscosonic wave equation, and has broader application prospects.

[0151] Table 3

[0152]

[0153] 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.

[0154] 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.

[0155] 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.

[0156] 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 highly efficient method for compensating for viscous acoustic wave attenuation based on weighted function approximation, characterized in that, The method includes: S1. By introducing a weighting function approximation, the viscous acoustic wave attenuation equation applicable to finite difference solutions is derived. S2, Based on the aforementioned viscous acoustic wave attenuation equation, establish a viscous acoustic wave compensation equation to compensate for the energy attenuation of seismic waves propagating in a viscous medium. S3, implement the attenuation-compensated reverse time migration imaging process, in which the wavefield compensation term is combined with the low-pass filter calculation to achieve wavefield stabilization compensation; S4, Configure the source and receiver settings to perform reverse time migration imaging; wherein, step S1 includes: S11, According to the approximate constant Q model, the complex modulus of the viscous medium Represented as: （1） in, For reference modulus, ρ For density, v 0 is the reference speed. Q For quality factors, ω and ω 0 represents the angular frequency and the reference angular frequency, respectively, and i is the imaginary unit; S12, drawing on the complex modulus expression of the generalized standard linear volume model, defines the following complex weighting function: （2） in, τ σl and τ εl These are the stress relaxation time and strain relaxation time, which are independent of the Q value, respectively. W R and W I These are the real and imaginary parts of the complex weighting function, respectively. L The weighted number; S13, use Approximating the part in parentheses in equation (1), we get: （3） Among them, stress relaxation time τ σl , Response relaxation time τ εl This is obtained by solving an optimization problem; S14, Rewrite it in the following form: （4） in, （5） Therefore, the complex modulus in equation (1) Rewritten as: （6）。 2. The efficient viscoacoustic wave attenuation compensation reverse-time migration method based on weighted function approximation according to claim 1, characterized in that, Step S1 also includes: S15, in the frequency domain, the two-dimensional viscous acoustic wave equation is expressed as: （7） in, For frequency domain wave fields, ρ For density, x and z Spatial location coordinates, ω Angular frequency; S16, Substitute the complex modulus in equation (6) into equation (7), and define the following auxiliary variables. : （8） Equations (7) and (8) can be rewritten in the following form: （9a） （9b） S17, Transforming equation (9) back to the time domain, we obtain the viscous acoustic wave attenuation equation: （10） Where t is the wave field time, （11） Equation (10) can be solved directly using finite difference.

3. The efficient viscoacoustic wave attenuation compensation reverse-time migration method based on weighted function approximation according to claim 2, characterized in that, Step S2 includes: S21, starting from the complex modulus, i.e., equation (1), its imaginary part Related to amplitude decay, therefore, the complex modulus with compensation effect The expression is as follows: （12） S22, combined with compensated complex modulus The frequency-wavenumber domain compensated viscous acoustic wave equation can be written in the following form: （13） in For frequency domain wave fields, ρ For density, k Wave number; S23, Substitute the complex modulus in equation (12) into equation (13): （14） in g and The expression is given in equation (5); S24, combined with approximate expression ,in, k For wave number, For the reference velocity, equation (14) is rewritten as: （15） S25, transforming equation (15) back to the frequency-space domain, we get: （16） S26, define the following auxiliary variables: （17） S27, combining auxiliary variables and transforming equation (17) back to the time-space domain, we get: （18） （19） in, F -1 For the inverse Fourier transform, the other coefficients are as follows: （20） The first and second time partial derivatives in the formula are calculated using the finite difference method.

4. The efficient viscoacoustic wave attenuation compensation reverse-time migration method based on weighted function approximation according to claim 3, characterized in that, Step S3 includes: S31, using the viscoacoustic attenuation equation to positively extend the source wave field; S32 uses the viscous acoustic wave compensation equation to reverse extend the detector wave field. In order to prevent the wave field compensation from being unstable, the amplitude compensation term is subjected to low-pass filtering. The low-pass filtering is combined with the inverse Fourier transform in the compensation equation for calculation. S33, applying source-normalized cross-correlation imaging conditions to the source wavefield and detector wavefield, obtaining Q-RTM imaging results: （21） in For the imaging results, T To record length, For the source wave field, For the detector wave field, ε To avoid dividing the denominator by a decimal.

5. The efficient viscoacoustic wave attenuation compensation reverse-time migration method based on weighted function approximation according to claim 1, characterized in that, Step S4 includes: Based on actual needs, select the source wavelet type and frequency that match the actual situation, set the detector to receive the seismic signal, set the receiving recording time length and sampling interval according to the model size and spatial interval, and then start the reverse time migration calculation to finally obtain the subsurface medium migration imaging results.

Images