A numerical simulation method of power-law frequency-dependent Q effect decoupled from dispersion and attenuation

Abstract

The application discloses a power-law frequency-dependent Q effect numerical simulation method which is decoupled from frequency dispersion and attenuation, and comprises the following steps: S1, deducing complex modulus and viscoacoustic wave equation based on power-law frequency-dependent Q effect; S2, solving the viscoacoustic wave equation by using a numerical method; S3, stability condition analysis; and S4, setting a seismic source and a receiver. The method overcomes the defects of the prior art numerical simulation method of viscoacoustic waves, that is, the prior art numerical simulation method of viscoacoustic waves can only simulate constant Q effect and ignores the characteristics that Q changes with frequency. The method not only realizes numerical simulation of power-law frequency-dependent Q effect in a full frequency band, but also has the advantages of amplitude attenuation and phase dispersion being decoupled from each other, and the calculation amount is equivalent to that of a conventional constant Q method. The method can be further applied to Q compensation reverse time migration imaging and waveform inversion, and has important significance for accurate imaging and interpretation of underground structures.

Description

Technical Field

[0001] This invention relates to the field of seismic wave numerical simulation and imaging technology, specifically to a numerical simulation method for power-law frequency-varying Q-effects with decoupling of dispersion and attenuation. Background Technology

[0002] Due to the viscosity of the medium, seismic waves inevitably undergo absorption and attenuation as they propagate through the Earth's interior. This absorption and attenuation reduces the amplitude of seismic waves and distorts their phase. Accurately describing the absorption and attenuation effect of seismic waves is crucial for imaging and interpreting subsurface structures. Conventional numerical simulation methods for seismic wave absorption and attenuation (such as standard linear bodies and fractional derivative models) are mainly based on the assumption of constant Q (independent of frequency). This assumption is generally only valid under normal temperature and pressure conditions. However, in high-temperature and high-pressure environments or fluid-containing media, the quality factor Q changes with frequency and follows a power-law relationship. Therefore, further development of numerical simulation methods for wave equations that address the power-law frequency-varying Q effect is of great significance for studying the propagation laws of seismic waves in real media and improving the resolution of subsurface media imaging.

[0003] Viscosity is a fundamental property of media. When seismic waves propagate in viscous media, internal friction converts some mechanical energy into heat, resulting in energy loss and attenuation of seismic wave amplitude, which adversely affects seismic imaging and interpretation. The strength of absorption attenuation in a medium can be quantitatively characterized by defining a quality factor Q. To simulate the absorption attenuation effect of seismic waves during propagation, some scholars have proposed standard linear volume models and fractional derivative models, which have been widely used in seismic exploration. However, these methods are mainly based on the constant Q assumption, that is, the assumption that Q does not change with frequency. This assumption usually only matches the results observed in a limited frequency range under normal temperature and pressure conditions. In high-temperature and high-pressure environments or fluid-containing media, Q usually changes with frequency and follows a power-law relationship. For example, Berckhemer et al. (1979) observed creep in mantle peridotite and found that Q has a power-law relationship with frequency. Anderson and Minster (1979) proved that Chandler oscillation, tidal dissipation, and free oscillation data are consistent with the power-law dependence of Q on frequency. Ignoring the frequency variation of Q in seismic wave numerical simulations will inevitably lead to simulation errors, affecting seismic imaging, inversion, and interpretation.

[0004] To simulate the power-law frequency-varying Q-effect of seismic waves, a few scholars have proposed using the Q-value simulated by the relaxation mechanism of traditional mechanical models (standard linear bodies) to approximate the power-law frequency-varying Q-effect within the target frequency band. However, this method can only approximate the effect within a limited frequency band. Moreover, to obtain a better fit, multiple standard linear body models are usually superimposed, which complicates the equations and increases computational and memory consumption. On the other hand, the quality factor Q is not explicitly expressed in the wave equation, which is not conducive to parameter inversion. In addition, the amplitude attenuation and velocity dispersion effects simulated by this type of equation are coupled and difficult to separate, making it unsuitable for Q-compensated reverse time migration imaging. Therefore, there is an urgent need to develop a numerical simulation method for the power-law frequency-varying Q-effect that decouples dispersion and attenuation. Summary of the Invention

[0005] Existing numerical simulation methods for viscoacoustic waves can only simulate the constant Q effect, neglecting the frequency-dependent characteristics of Q. To overcome the shortcomings of existing technologies, this invention provides a numerical simulation method for the power-law frequency-varying Q effect with decoupled dispersion and attenuation. This method not only achieves numerical simulation of the power-law frequency-varying Q effect across the entire frequency band but also has the advantage of decoupling amplitude attenuation and phase dispersion, with computational complexity comparable to traditional constant Q methods. This invention can be further applied to Q-compensated reverse time migration imaging and waveform inversion, which is of great significance for accurate imaging and interpretation of subsurface structures. To achieve the above technical objectives, this invention adopts the following technical solution:

[0006] A numerical simulation method for power-law frequency-varying Q-effects with dispersion and attenuation decoupled includes the following steps:

[0007] S1: Derive the complex modulus and viscous acoustic wave equations based on the power-law frequency-varying Q effect;

[0008] S2: Solve the equation of viscous acoustic waves using numerical methods;

[0009] S3: Stability condition analysis;

[0010] S4: Source and receiver settings.

[0011] Preferably, step S1 includes:

[0012] S11: In the power-law frequency-varying Q model, the quality factor Q changes with frequency, expressed as:

[0013]

[0014] Where Q0 is the reference quality factor, ω and ω0 are the angular frequency and the reference angular frequency, respectively, and χ is the fractional exponent.

[0015] S12: From the Kramers-Kronig relation, we obtain an approximate expression for the complex modulus:

[0016]

[0017] Where M0 is the reference modulus, i is the imaginary unit, and the real and imaginary parts of the complex modulus represent phase dispersion and amplitude attenuation, respectively.

[0018] S13: Combining the approximate relationship between wave number k and angular frequency ω ≈ kc0, where k = |k| is the modulus of the wave number and c0 is the reference velocity, the complex modulus is rewritten in the following form:

[0019]

[0020] S14: In the two-dimensional model, the first-order momentum conservation equation is expressed as:

[0021]

[0022] Where v x v z Let σ be the velocity component, σ be the stress, and ρ be the density. The first-order time partial derivative, and It is the first-order spatial partial derivative;

[0023] S15: Relationship between strain ε and velocity The equation for viscous acoustic waves is derived as follows:

[0024]

[0025] in For the Laplace operator, the other parameters can be expressed as:

[0026]

[0027] S16: Ignoring amplitude attenuation and retaining phase dispersion, the viscous acoustic wave equation becomes:

[0028]

[0029] S17: Ignore phase dispersion, at this time The equation for viscous acoustic waves then becomes:

[0030]

[0031] Preferably, step S2 includes:

[0032] S21: For the variational order Laplace operator in the viscosonic wave equation (5) and First, we introduce a first-order Taylor series expansion to approximate the wavenumber domain operator corresponding to the above operator, that is:

[0033] k 1-χ ≈k(1-χlnk) (9)

[0034] k 2-χ ≈k 2 (1-χlnk) (10)

[0035] S22: Based on equations (9) and (10), the wavenumber-related operators are solved using the inverse Fourier transform, yielding the following results:

[0036]

[0037]

[0038] in, and These represent the forward and inverse Fourier transforms, respectively.

[0039] S23: Combining equations (11) and (12), the solution to the viscous acoustic wave equation (5) is as follows:

[0040]

[0041] S24: Solve the time partial derivatives in equation (13) using finite difference methods:

[0042]

[0043]

[0044] Where x = (x, z) is the spatial coordinate vector, t is the current time, and Δt is the time sampling interval.

[0045] Preferably, step S3 includes:

[0046] In the numerical simulation process, it is necessary to analyze the stability conditions of the equations in order to set the time sampling interval appropriately.

[0047] S31: Transforming equation (5) to the wavenumber domain and combining it with equations (14) and (15), we obtain:

[0048] σ n+1 =2σ n -σ n-1 -Δt 2 λk 2 σ n +Δt 2 ηk 2-χ σ n -Δtτk 1-χ (σ n -σ n-1 (16)

[0049] Where, σ n+1 =σ((n+1)Δt), σ n =σ(nΔt) and σ n-1 =σ((n-1)Δt) represents the stress wave field at times (n+1)Δt, nΔt, and (n-1)Δt, respectively;

[0050] S32: Rewrite equation (16) in matrix form:

[0051]

[0052] Where a=-λk 2 +ηk 2-χ b=-τk 1-χ ;

[0053] S33: Solving for the characteristic values ​​of equation (17) yields:

[0054]

[0055] S34: Setting the eigenvalue |ψ|≤1 yields the following stability condition:

[0056]

[0057] Based on stability conditions, and considering the model size and grid spacing, the time sampling interval Δt and the simulation record length are determined.

[0058] Preferably, step S4 includes: selecting the source wavelet and wavelet frequency according to actual needs, and setting the detector to receive the seismic signal. Once the above parameter settings are completed, numerical simulation can begin to obtain the seismic wave field and seismic record.

[0059] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0060] The method of this invention can achieve numerical simulation of the power-law frequency-varying Q-effect. Its advantages are as follows:

[0061] (1) Regarding the frequency-varying characteristics of Q: Traditional numerical simulation methods for viscous acoustic wave equations are based on the constant Q assumption and cannot simulate the frequency-varying characteristics of Q. This invention starts directly from the power-law frequency-varying Q model, derives a new viscous acoustic wave equation (Equation 5), and solves it using numerical methods (Equations 13–15), which can directly simulate the power-law frequency-varying Q effect across the entire frequency band.

[0062] (2) In terms of calculation method and efficiency: The traditional constant Q model method and the method described in this invention have similar forms of viscous acoustic wave equations, both containing two variable fractional Laplace operators. The numerical solution methods are similar (combining Taylor series expansion and Fourier transform), and the calculation efficiency is comparable.

[0063] (3) Decoupling of dispersion and attenuation: Similar to the traditional constant Q model method, this invention separates the real and imaginary parts of the complex modulus of the power-law frequency-varying Q model, deriving equations that include only phase dispersion and only amplitude attenuation effects, respectively. Using equations 7 and 8, seismic wave fields and seismic records with phase dispersion and amplitude attenuation effects can be simulated independently, respectively. Compared with the acoustic wave results, the amplitude calculated by the amplitude attenuation equation shows a large attenuation, while the phase remains basically the same; in contrast, the amplitude calculated by the phase dispersion equation is basically the same as the acoustic wave amplitude, while the phase shows a significant difference. Attached Figure Description

[0064] To more clearly illustrate the technical solutions of the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0065] Figure 1(a) shows the parameter-reference velocity c0 parameter diagram of the complex Marmousi model.

[0066] Figure 1(b) shows the parameters of the complex Marmousi model and the reference quality factor Q0.

[0067] Figure 1(c) shows the parameter-fractional χ parameter plot of the complex Marmousi model.

[0068] Figure 2 This is a snapshot of the wave field at 1.6s for the complex Marmousi model.

[0069] Figure 3 This is a seismic record diagram of a complex Marmousi model.

[0070] Figure 4 A comparison chart of single-track records for a complex Marmousi model. Detailed Implementation

[0071] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0072] A numerical simulation method for power-law frequency-varying Q-effects with dispersion and attenuation decoupled includes the following steps:

[0073] S1: Derive the complex modulus and viscous acoustic wave equations based on the power-law frequency-varying Q effect;

[0074] S11: In the power-law frequency-varying Q model, the quality factor Q changes with frequency, expressed as:

[0075]

[0076] Where Q0 is the reference quality factor, ω and ω0 are the angular frequency and the reference angular frequency, respectively, and χ is the fractional exponent.

[0077] S12: From the Kramers-Kronig relation, we obtain an approximate expression for the complex modulus:

[0078]

[0079] Where M0 is the reference modulus, i is the imaginary unit, and the real and imaginary parts of the complex modulus represent phase dispersion and amplitude attenuation, respectively.

[0080] S13: Combining the approximate relationship between wave number k and angular frequency ω ≈ kc0, where k = |k| is the modulus of the wave number and c0 is the reference velocity, the complex modulus is rewritten in the following form:

[0081]

[0082] S14: In the two-dimensional model, the first-order momentum conservation equation is expressed as:

[0083]

[0084] Where v x v z Let σ be the velocity component, σ be the stress, and ρ be the density. The first-order time partial derivative, and It is the first-order spatial partial derivative;

[0085] S15: Relationship between strain ε and velocity The equation for viscous acoustic waves is derived as follows:

[0086]

[0087] in For the Laplace operator, the other parameters can be expressed as:

[0088]

[0089] Since viscosity causes seismic wave amplitude attenuation and phase dispersion, in Q-compensated reverse time migration, in order to achieve amplitude compensation and phase correction, the two need to be decoupled. For this purpose, a decoupling equation is used to simulate the amplitude attenuation and phase dispersion effects separately. The first and second terms on the right side of equation (5) are related to phase dispersion, and the third term is related to amplitude attenuation.

[0090] S16: Ignoring amplitude attenuation and retaining phase dispersion, the viscous acoustic wave equation becomes:

[0091]

[0092] S17: Ignore phase dispersion, at this time The equation for viscous acoustic waves then becomes:

[0093]

[0094] S2: Solve the equation of viscous acoustic waves using numerical methods;

[0095] S21: For the variational order Laplace operator in the viscosonic wave equation (5) and First, we introduce a first-order Taylor series expansion to approximate the wavenumber domain operator corresponding to the above operator, that is:

[0096] k 1- χ≈k(1-χlnk)(9)

[0097] k 2- χ≈k 2 (1-χlnk)(10)

[0098] S22: Based on equations (9) and (10), the wavenumber-related operators are solved using the inverse Fourier transform, yielding the following results:

[0099]

[0100]

[0101] in, and These represent the forward and inverse Fourier transforms, respectively.

[0102] S23: Combining equations (11) and (12), the solution to the viscous acoustic wave equation (5) is as follows:

[0103]

[0104] S24: Solve the time partial derivatives in equation (13) using finite difference methods:

[0105]

[0106]

[0107] Where x = (x, z) is the spatial coordinate vector, t is the current time, and Δt is the time sampling interval.

[0108] S3: Stability condition analysis;

[0109] In the numerical simulation process, it is necessary to analyze the stability conditions of the equations in order to set the time sampling interval appropriately.

[0110] S31: Transforming equation (5) to the wavenumber domain and combining it with equations (14) and (15), we obtain:

[0111] σ n+1 =2σ n -σ n-1 -Δt 2 λk 2 σ n +Δt 2 ηk 2-χ σ n -Δtτk 1-χ (σ n -σ n-1 (16)

[0112] Where, σ n+1 =σ((n+1)Δt), σ n =σ(nΔt) and σ n-1 =σ((n-1)Δt) represents the stress wave field at times (n+1)Δt, nΔt, and (n-1)Δt, respectively;

[0113] S32: Rewrite equation (16) in matrix form:

[0114]

[0115] Where a=-λk 2 +ηk 2-χ b=-τk 1-χ ;

[0116] S33: Solving for the characteristic values ​​of equation (17) yields:

[0117]

[0118] S34: Setting the eigenvalue |ψ|≤1 yields the following stability condition:

[0119]

[0120] Based on stability conditions, and considering the model size and grid spacing, the time sampling interval Δt and the simulation record length are determined.

[0121] S4: Source and receiver settings: Select the source wavelet and wavelet frequency according to actual needs, and set the detector to receive the seismic signal. Once the above parameter settings are completed, the numerical simulation can be started to obtain the seismic wave field and seismic record.

[0122] Using the complex Marmousi model as an example, a numerical simulation test of the power-law frequency-varying Q-effect was conducted. The Marmousi model has a computational grid size of 501×353, with a grid spacing of 12 meters both horizontally and vertically. The specific steps are as follows:

[0123] As shown in Figure 1, based on the Marmousi model reference velocity c0, reference quality factor Q0, and fractional-order χ parameter, combined with the stability condition of equation (19), the time sampling interval is set to 1 millisecond, the recording length is 2 seconds, and the reference angular frequency ω0 is 2π radians / second; the seismic source adopts a Ricker wavelet with a main frequency of 20Hz, which is set at (3000, 12) meters to excite seismic waves, and the detector is set on the ground to receive seismic signals; according to equation (13), numerical simulation is performed and the seismic wave field and record are stored; then the numerical simulation results are compared and analyzed. Figure 2 A snapshot of the wave field at time 1.6s is shown. Figure 2 (a) shows the results calculated using the constant Q method (χ=0). Figure 2 (b) and (c) show the results of the numerical simulation method for the power-law frequency-varying Q-effect. Among them, Figure 2 (b) A constant fractional-order model (χ = 0.3) is used. Figure 2 (c) A true fractional-order χ² model is used. When the frequency variation effect of Q is ignored, Figure 2 (a) and Figure 2 There is a significant difference between (c) and its amplitude decay is more severe. Figure 3 The corresponding earthquake records are displayed, from which similar phenomena can be observed. Figure 2 This phenomenon. To compare the differences in amplitude and phase between different calculation results, we analyzed the following: Figure 3 Single tracks were sampled at lateral distances x = 1800, 2520, and 4200 meters for comparison. Figure 4 As shown, the curve calculated by the constant Q method (χ=0) has the smallest amplitude and the most advanced phase, indicating that the constant Q method suffers the most severe absorption attenuation. The method of this invention fully considers the frequency-varying effect of Q; its actual simulated Q value is larger than the reference Q value Q0, and the absorption attenuation is weaker than the constant Q result, thus resulting in a larger amplitude. It is evident that traditional methods based on the constant Q assumption ignore the frequency-varying characteristics of Q, leading to significant discrepancies between their simulation results and actual conditions. The method described in this invention fully considers the frequency-varying characteristics of Q, accurately simulating the propagation law and wavefield characteristics of seismic waves in the power-law frequency-varying Q effect model, which is of great significance for improving seismic imaging and interpretation.

Claims

1. A numerical simulation method for power-law frequency-varying Q-effects with dispersion and attenuation decoupled, comprising the following steps: S1: Derivation of the complex modulus and viscous acoustic wave equations based on the power-law frequency-varying Q effect, including: S11: In the power-law frequency-varying Q model, the quality factor Q changes with frequency, expressed as: （1） in For reference quality factor, and These are the angular frequency and the reference angular frequency, respectively. It is a fractional exponent; S12: From the Kramers-Kronig relation, we obtain an approximate expression for the complex modulus: （2） in The reference modulus is i, where i is the imaginary unit. The real and imaginary parts of the complex modulus represent phase dispersion and amplitude attenuation, respectively. S13: Combining wavenumber k and angular frequency Approximate relationship between ,in The modulus of the wave number, For reference speed, the complex modulus is rewritten in the following form: （3） S14: In the two-dimensional model, the first-order momentum conservation equation is expressed as: （4） in For velocity components, For stress, For density, The first-order time partial derivative, and It is the first-order spatial partial derivative; S15: Combined strain Relationship with speed The equation for viscous acoustic waves is derived as follows: （5） in For the Laplace operator, the other parameters are expressed as: （6） S16: Ignoring amplitude attenuation and retaining phase dispersion, the viscous acoustic wave equation becomes: （7） S17: Ignore phase dispersion, at this time The equation for viscous acoustic waves then becomes: （8）； S2: Solve the equation of viscous acoustic waves using numerical methods; S3: Stability condition analysis; S4: Source and receiver settings.

2. The numerical simulation method for power-law frequency-varying Q-effect with dispersion and attenuation decoupling according to claim 1, characterized in that, Step S2 includes: S21: For the variational order Laplace operator in the viscosonic wave equation (5) and First, we introduce a first-order Taylor series expansion to approximate the wavenumber domain operator corresponding to the above operator, that is: （9） （10） S22: Based on equations (9) and (10), the wavenumber-related operators are solved using the inverse Fourier transform, yielding the following results: （11） （12） in, and These represent the forward and inverse Fourier transforms, respectively. S23: Combining equations (11) and (12), the solution to the viscous acoustic wave equation (5) is as follows: （13） S24: Solve the time partial derivatives in equation (13) using finite difference methods: （14） （15） in For spatial coordinate vectors, For the current moment, The time sampling interval is denoted as .

3. The numerical simulation method for power-law frequency-varying Q-effect with dispersion and attenuation decoupling according to claim 1, characterized in that, Step S3 includes: In the numerical simulation process, it is necessary to analyze the stability conditions of the equations in order to set the time sampling interval appropriately. S31: Transforming equation (5) to the wavenumber domain and combining it with equations (14) and (15), we get: （16） in, , and Representing the first , and Stress wave field at time; S32: Rewrite equation (16) in matrix form: （17） in ; S33: Solving for the characteristic values ​​of equation (17) yields: （18） S34: Let the eigenvalues The following stability conditions can be obtained: （19） Based on stability conditions, and considering the model size and grid spacing, the time sampling interval is determined. And the simulated record length.

4. The numerical simulation method for power-law frequency-varying Q-effect with dispersion and attenuation decoupling according to claim 1, characterized in that, Step S4 includes: selecting the source wavelet and wavelet frequency according to actual needs, and setting the detector to receive the seismic signal. Once the above parameter settings are completed, numerical simulation can begin to obtain the seismic wavefield and seismic record.

Images