A viscoelastic scalar p-wave equation construction and wave field numerical simulation method

Abstract

The application discloses a kind of viscoelastic scalar P-wave equation construction and wave field numerical simulation method, belong to deep sea oil and gas field exploration development field.For the problems of insufficient physical completeness of existing scalar acoustic equation, high cost of multi-component elastic wave equation calculation, and difficulty of viscoelastic equation of constant Q model to fully reflect the wave field dynamics characteristics, the application establishes the dispersion relation based on constant Q model, constructs the viscoelastic multi-component wave equation containing fractional order space differential operator and time derivative coupling term, then introduces the dual projection operator to realize P-wave and S-wave mode decoupling, obtains the viscoelastic scalar equation describing only P-wave, and uses pseudo-spectral method and finite difference method for numerical simulation.The method can more accurately depict the propagation dynamics of seismic wave in strong attenuation medium, effectively decouples amplitude attenuation and phase dispersion effect, and provides reliable theory and calculation tool for high-precision seismic data processing and imaging.

Description

Technical Field

[0001] This invention belongs to the field of deep-sea oil and gas field exploration and development technology, and particularly relates to a method for constructing viscoelastic scalar P-wave equations and numerically simulating wave fields. Background Technology

[0002] In the field of deep-sea oil and gas field exploration and development, the wave equation is the core mathematical tool for characterizing the propagation mechanism of seismic waves in subsurface media, and its accuracy directly determines the upper limit of the accuracy of subsequent seismic data processing and interpretation. For a long time, the scalar acoustic equation, by simplifying the approximation of the elastic wave field with a scalar field, has significantly reduced the computational burden in forward modeling and migration imaging of complex geological structures. However, as exploration targets shift towards more challenging strongly attenuating media and observation systems develop towards multi-component systems, this simple acoustic approximation shows significant shortcomings in physical completeness, making it difficult to truly reflect the elastic effects of wave propagation in multi-component data acquisition and amplitude-preserving processing scenarios.

[0003] To address the need for physical interpretation of complex media, a theory of multi-component elastic wave equations has been developed. While this theory provides a more comprehensive physical foundation for wavefield description, it commonly faces challenges in practical wavefield simulations, such as crosstalk between P-wave and S-wave modes, polarity reversal, and difficulty in decoupling multi-component wavefields. Since shear wave velocities in subsurface media are typically lower than P-wave velocities, and shear wave wavelengths are shorter at the same dominant frequency, this imposes more stringent dispersion conditions and stability constraints on mesh generation and time steps in numerical simulations. Using the traditional finite difference method for numerical solutions presents significant challenges in terms of computational cost and time consumption, making it difficult to meet the efficiency requirements of actual production.

[0004] Furthermore, to describe the propagation of seismic waves in attenuating media, a forward modeling technique based on the constant Q model of viscoelastic wave equations has been developed. While this technique can simulate the energy loss process of waves to some extent, its limitations in theoretical foundations and numerical implementation methods prevent it from effectively and consistently characterizing key dynamic features such as amplitude attenuation, phase changes, and dispersion during seismic wave propagation. This discrepancy between the simulation results and the actual wavefield dynamics significantly hinders subsequent high-precision seismic data processing, structural interpretation, and reservoir parameter inversion, leading to uncertainties in the assessment of subsurface structures and fluid distribution. Therefore, how to construct a high-fidelity modeling technique that accurately characterizes the dynamic behavior of seismic waves, especially P-waves, in strongly attenuating media while ensuring computational efficiency has become a pressing technical challenge in this field. Summary of the Invention

[0005] To address the aforementioned technical problems, this invention proposes a method for constructing viscoelastic scalar P-wave equations and numerically simulating wave fields, thereby resolving the issues present in the prior art.

[0006] Firstly, to achieve the above objectives, this invention provides a method for constructing viscoelastic scalar P-wave equations and numerically simulating wave fields, comprising the following steps:

[0007] Establish the dispersion relationship of viscoelastic media based on the constant Q model;

[0008] Based on the aforementioned dispersion relation of the viscoelastic medium, a viscoelastic multi-component wave equation containing a coupling term between a fractional-order spatial differential operator and a time derivative is constructed.

[0009] By introducing a dual projection operator to decouple the P-wave and S-wave modes of the viscoelastic multi-component wave equation, a viscoelastic scalar P-wave equation describing only the dynamic characteristics of the P-wave is obtained.

[0010] The viscoelastic scalar P-wave equations were numerically discretized using the pseudospectral method and the finite difference method to achieve forward modeling of the seismic wave field.

[0011] Optionally, the process of establishing the dispersion relation of viscoelastic media based on the constant Q model includes:

[0012] By treating the quality factor of the medium as a constant within the seismic exploration frequency band, we derive the complex wave number expression characterizing dispersion and attenuation.

[0013] The complex wavenumber expression is truncated and Taylorized within the target frequency band, transforming the frequency-wavenumber relationship into an expression composed of integer powers of the wavenumber and powers of the frequency.

[0014] By utilizing Fourier duality, the combinatorial expression is mapped into a constant-order partial differential equation in the time-space domain.

[0015] Optionally, the process of constructing the viscoelastic multicomponent wave equation includes:

[0016] Complex moduli related to frequency and wavenumber are introduced for P-waves and S-waves, respectively, and the complex moduli are determined by the dispersion relation of the viscoelastic medium.

[0017] Based on the differential properties of the Fourier transform, the complex modulus is transformed into the time-space domain to obtain a mechanical modulus expression that includes the fractional Laplace operator and the time partial derivative.

[0018] Based on the stress-strain constitutive relation, strain-displacement relation, and momentum conservation equation, the viscoelastic multi-component wave equation with displacement field as unknown is obtained by combining the equations.

[0019] Optionally, the process of decoupling P-wave and S-wave modes includes:

[0020] Construct a normalized gradient operator and its corresponding divergence operator, wherein the normalized gradient operator and the divergence operator satisfy the relationship that they are inverse operators of each other;

[0021] The displacement field is projected using the divergence operator to define a scalar P-wave field.

[0022] The dual projection formed by the normalized gradient operator and the divergence operator is applied to both sides of the viscoelastic multicomponent wave equation to eliminate the modes related to the S-wave and derive the viscoelastic scalar P-wave equation.

[0023] Optionally, the process of numerical discretization using the pseudospectral method and the finite difference method includes:

[0024] In the wavenumber domain, the spatial derivative and fractional Laplace operator in the viscoelastic scalar P-wave equation are calculated using Fourier transform and inverse transform.

[0025] In the time domain, the time derivative term in the viscoelastic scalar P-wave equation is discretized and approximated using a central difference scheme.

[0026] The wave field is iteratively propagated by alternately performing spatial derivative calculations in the wavenumber domain and time-step updates in the time domain.

[0027] Optionally, the method further includes:

[0028] A source term containing a volume force source term or a moment tensor source term is introduced into the viscoelastic scalar P-wave equation;

[0029] Based on the mathematical form of the source term in the wavenumber-frequency domain, the control law of different source loading methods on the wavefield radiation pattern and amplitude distribution is analyzed.

[0030] Optionally, the mechanical modulus is:

[0031] ;

[0032] in, For mechanical modulus, For density, A coefficient vector related to the medium velocity and Q value. For the Laplace operator, This is the time partial derivative.

[0033] Optionally, the definitions of the normalized gradient operator and the divergence operator are as follows:

[0034] ;

[0035] ;

[0036] in, For the normalized gradient operator, For space directional partial derivative, For space directional partial derivative, It is a divergence operator.

[0037] Secondly, the present invention also provides a computer terminal device, comprising:

[0038] One or more processors;

[0039] A memory, coupled to the processor, for storing one or more programs;

[0040] When the one or more programs are executed by the one or more processors, the one or more processors implement the steps of the viscoelastic scalar P-wave equation construction and wave field numerical simulation method in the first aspect above.

[0041] Thirdly, the present invention also provides a computer-readable storage medium having a computer program stored thereon, wherein when the computer program is executed by a processor, it implements the steps of the viscoelastic scalar P-wave equation construction and wave field numerical simulation method in the first aspect described above.

[0042] Compared with the prior art, the present invention has the following advantages and technical effects:

[0043] This invention provides a viscoelastic scalar P-wave equation construction and wavefield numerical simulation method. By considering the viscous characteristics of the subsurface medium and introducing five key physical parameters, it can more precisely characterize the seismic wave propagation process in non-uniform, strongly attenuated media. The wavefield constructed by this method maintains a high degree of consistency with the synthetic P-wave obtained from the viscoelastic full-wave equation in terms of kinematic and dynamic characteristics while preserving amplitude. Furthermore, the derived equations successfully decouple amplitude attenuation and phase dispersion effects at the operator level. Compared to existing technologies, this invention retains the key physical influence of viscoelasticity on P-wave amplitude and phase while significantly reducing computational redundancy and complexity by avoiding explicit solutions for S-waves, thus achieving more accurate dynamic simulation capabilities. This scheme provides a physical modeling foundation for high-fidelity seismic imaging in strongly attenuated media and can be effectively integrated with quality factor compensated reverse-time migration procedures. Attached Figure Description

[0044] The accompanying drawings, which form part of this invention, are used to provide a further understanding of the invention. The illustrative embodiments of the invention and their descriptions are used to explain the invention and do not constitute an undue limitation of the invention. In the drawings:

[0045] Figure 1 The diagram shows the scalar P-wave equation for viscoelastic media and the wave field diagram of synthesized P-waves in viscoelastic media according to an embodiment of the present invention, wherein (a) is the forward modeling simulation diagram of the scalar P-wave and (b) is the forward modeling simulation diagram of the synthesized P-wave.

[0046] Figure 2This is a scalar P-wave combined wavefield diagram according to an embodiment of the present invention;

[0047] Figure 3 This is a single-channel diagram of the combined wavefield according to an embodiment of the present invention;

[0048] Figure 4 This is a schematic diagram of the simulation results of the synthetic P-wave and scalar P-wave fields, where the transverse wave velocity is reduced to Vs=500m / s while keeping the other modeling parameters unchanged. (a) is the simulation result of the synthetic P-wave, and (b) is the simulation result of the scalar P-wave field.

[0049] Figure 5 These are wavefield diagrams corresponding to different source loading methods in embodiments of the present invention, where (a) represents the wavefield diagram with only vertical body force applied. (b) is the simultaneous application of vertical body forces. With longitudinal moment components (c) is the loading of only the longitudinal moment component. ;

[0050] Figure 6 These are single-channel wavefield diagrams corresponding to different source loading methods in embodiments of the present invention, wherein (a) represents the wavefield diagram with only vertical body force applied. (b) is the simultaneous application of vertical body forces. With longitudinal moment components (c) is the loading of only the longitudinal moment component. ;

[0051] Figure 7 These are radiation pattern diagrams for different source loading methods according to embodiments of the present invention, where (a) represents loading only vertical body force. (b) is the simultaneous application of vertical body forces. With longitudinal moment components (c) is the loading of only the longitudinal moment component. . Detailed Implementation

[0052] It should be noted that, unless otherwise specified, the embodiments and features described in the present invention can be combined with each other. The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0053] It should be noted that the steps shown in the flowchart in the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although a logical order is shown in the flowchart, in some cases the steps shown or described may be executed in a different order than that shown here.

[0054] Example 1

[0055] like Figure 1As shown, this embodiment provides a method for constructing viscoelastic scalar P-wave equations and numerically simulating wave fields, including:

[0056] Establish the dispersion relationship of viscoelastic media based on the constant Q model;

[0057] Based on the aforementioned dispersion relation of the viscoelastic medium, a viscoelastic multi-component wave equation containing a coupling term between a fractional-order spatial differential operator and a time derivative is constructed.

[0058] By introducing a dual projection operator to decouple the P-wave and S-wave modes of the viscoelastic multi-component wave equation, a viscoelastic scalar P-wave equation describing only the dynamic characteristics of the P-wave is obtained.

[0059] The viscoelastic scalar P-wave equations were numerically discretized using the pseudospectral method and the finite difference method to achieve forward modeling of the seismic wave field.

[0060] Furthermore, the process of establishing the dispersion relation of viscoelastic media based on the constant Q model includes:

[0061] By treating the quality factor of the medium as a constant within the seismic exploration frequency band, we derive the complex wave number expression characterizing dispersion and attenuation.

[0062] The complex wavenumber expression is truncated and Taylorized within the target frequency band, transforming the frequency-wavenumber relationship into an expression composed of integer powers of the wavenumber and powers of the frequency.

[0063] By utilizing Fourier duality, the combinatorial expression is mapped into a constant-order partial differential equation in the time-space domain.

[0064] Specifically, the implementation process of this embodiment includes:

[0065] Within the seismic exploration frequency band, the quality factor Q changes slowly with frequency and can be approximated as a constant. Its complex modulus is:

[0066] (1);

[0067] In the formula, It is the complex modulus. At the reference frequency Reference modulus at that location It is an imaginary unit, representing attenuation intensity. can be derived from formula We find that Q is the quality factor, which determines the strength of medium dissipation. From this, we obtain the complex wave number of the plane wave dispersion relation:

[0068] (2);

[0069] In the formula, the complex wave number medium propagation speed , At the reference frequency The reference phase velocity is given by $e$, where $e$ represents an exponential function. Equation (2) shows that waves of different frequencies propagate at different speeds in the medium, which is the dispersion effect. The complex exponential term is related to energy attenuation. When $Q$ varies with space, the operator corresponding to equation (2) is a space-wavenumber hybrid domain operator, which has low accuracy in direct numerical calculation. Therefore, a truncated Taylor approximation is performed on equation (2) in the target frequency band, and the ω-k relationship is written as a combination:

[0070] (3);

[0071] Then, using Fourier duality, constant-order partial differential equations in the time-space domain are obtained. It's frequency. Determines the attenuation intensity. It's the speed of transmission. It is a reference frequency. It is the complex wave number. Since the power exponent does not depend on spatial coordinates, spatial discretization can be directly performed using the pseudospectral method without the need for global averaging of Q, thereby improving the accuracy of wave field simulation in non-uniform viscoelastic media.

[0072] Furthermore, the process of constructing the viscoelastic multi-component wave equation includes:

[0073] Complex moduli related to frequency and wavenumber are introduced for P-waves and S-waves, respectively, and the complex moduli are determined by the dispersion relation of the viscoelastic medium.

[0074] Based on the differential properties of the Fourier transform, the complex modulus is transformed into the time-space domain to obtain a mechanical modulus expression that includes the fractional Laplace operator and the time partial derivative.

[0075] Based on the stress-strain constitutive relation, strain-displacement relation, and momentum conservation equation, the viscoelastic multi-component wave equation with displacement field as unknown is obtained by combining the equations.

[0076] Specifically, the implementation process of this embodiment includes:

[0077] To simultaneously characterize the dispersion and absorption of P / S waves in viscoelastic media, a method is introduced... Frequency-dependent complex modulus , used to describe in Viscoelastic response under the mechanism.

[0078] (4);

[0079] In the formula, For density, , representing P-wave and S-wave respectively. Substituting equation (3) into equation (4) yields:

[0080] (5);

[0081] In the formula, , is a coefficient vector containing multiple , Mechanism-related constants. Equation (5) "encodes" the dispersion and attenuation terms into a set of arbitrary-order constants. and The exponentiation facilitates mapping to spatial-temporal differential operators. This construction is consistent with the fractional Laplace's approach in the viscous acoustic equation, achieving operator decoupling and stable numerical implementation of attenuation and dispersion.

[0082] According to Fourier duality, as shown in equation (6), where Representing wave number The power of, where It is the first Wave number corresponding to each mechanism: Represents the Laplace operator power, when When the number is not even, this is a fractional Laplace operator. Represents the imaginary unit Multiply by angular frequency , Indicates the first The angular frequency corresponding to each mechanism. Let be the time partial derivative, and represent the first-order time derivative in the time domain.

[0083] , (6);

[0084] Transform equation (5) into its spatiotemporal form:

[0085] (7);

[0086] in, For mechanical modulus, For density, A coefficient vector related to the medium velocity and Q value. For the Laplace operator, This is the time partial derivative.

[0087] Under the assumption of isotropic deformation, the stress-strain constitutive model is written as:

[0088] (8);

[0089] All of these are stresses, representing three different directions; For strain in three different directions These represent the mechanical modulus under P-wave and S-wave mechanisms, respectively.

[0090] The strain-displacement relationship is as follows:

[0091] (9);

[0092] In the formula, For displacements in the x and z directions, For space directional partial derivative, For space Directional partial derivative.

[0093] The momentum equation is:

[0094] (10);

[0095] It represents the second-order time partial derivative in the time domain.

[0096] Substituting (8) and (9) into (10) and eliminating stress and strain, we obtain the viscoelastic displacement equation:

[0097] (11);

[0098] Equation (11) is passed A fractional-order Laplace-time derivative coupling term is explicitly introduced, thereby controlling phase dispersion and amplitude decay simultaneously at the operator level.

[0099] Furthermore, the process of decoupling the P-wave and S-wave modes includes:

[0100] Construct a normalized gradient operator and its corresponding divergence operator, wherein the normalized gradient operator and the divergence operator satisfy the relationship that they are inverse operators of each other;

[0101] The displacement field is projected using the divergence operator to define a scalar P-wave field.

[0102] The dual projection formed by the normalized gradient operator and the divergence operator is applied to both sides of the viscoelastic multicomponent wave equation to eliminate the modes related to the S-wave and derive the viscoelastic scalar P-wave equation.

[0103] Specifically, the implementation process of this embodiment includes:

[0104] To obtain the wave equation containing only the scalar P mode in an isotropic viscoelastic medium, a normalized gradient and divergence operator pair is introduced as follows:

[0105] (12);

[0106] (13);

[0107] in, For the normalized gradient operator, For space directional partial derivative, For space directional partial derivative, Here, T represents the divergence operator, and T represents the matrix transpose operator.

[0108] Its key properties are:

[0109] (14);

[0110] I represents the unit operator, therefore and They are inverse operators:

[0111] (15);

[0112] To extract the scalar P mode, first define the projected scalar field of the P mode:

[0113] (16);

[0114] and utilize Construct an orthogonal projection of the vector field onto the spinless subspace, Insert to the right end of (11), then multiply both sides of the equation on the left. Thus, the viscoelastic scalar P-wave equation is obtained:

[0115] (17);

[0116] That is, at the operator level by and The constructed dual projection achieves orthogonal decoupling of the P / S modes. Because... , Retaining the dispersion characteristics of constant Q, the characterization of phase dispersion and amplitude attenuation in equation (17) is consistent with the projection of the full vector equation onto the P mode.

[0117] Furthermore, the process of numerical discretization using the pseudospectral method and the finite difference method includes:

[0118] In the wavenumber domain, the spatial derivative and fractional Laplace operator in the viscoelastic scalar P-wave equation are calculated using Fourier transform and inverse transform.

[0119] In the time domain, the time derivative term in the viscoelastic scalar P-wave equation is discretized and approximated using a central difference scheme.

[0120] The wave field is iteratively propagated by alternately performing spatial derivative calculations in the wavenumber domain and time-step updates in the time domain.

[0121] Specifically, the implementation process of this embodiment includes:

[0122] To numerically realize the fractional differential term, a pseudospectral method is used to transform the wave field from the time domain to the wavenumber domain, extending the spatial differential operator from integer to fractional order. For the two-dimensional wave equation, arbitrary components... The Fourier transform is:

[0123] (18);

[0124] F represents the positive Fourier transform. Using the differential properties of the Fourier transform, we obtain the derivative form of this arbitrary component in the wavenumber domain. Taking the x-direction as an example, the principle of finding the spatial derivative of the Fourier transform is as follows:

[0125] (19);

[0126] This represents performing an inverse Fourier transform. As shown above, to obtain the derivative of a variable with respect to the horizontal or vertical direction, one should first perform a Fourier transform on that variable at that point, and then compare it with the wavenumber and imaginary unit of the corresponding direction. Multiply them, then perform an inverse Fourier transform.

[0127] Taking the approximate solution of the stress term in equation system (8) as an example, Taking the Fourier transform of the right side of the equation, we get:

[0128] (20);

[0129] The stress and strain components are the frequency-wavenumber domain Fourier transform results corresponding to the spatiotemporal domain.

[0130] To calculate the first and second-order partial derivatives of the wave equation with respect to time, the wave field... Perform a Taylor expansion around the time coordinate t, and Abbreviation as ,exist:

[0131] (twenty one);

[0132] (twenty two);

[0133] In the formula, It is the step size in the time direction. , , These are the first, second, and third order partial derivatives of the wave field p in the time domain, respectively. Subtracting equations (21) and (22) above, we get:

[0134] (twenty three);

[0135] After sorting, we get:

[0136] (twenty four);

[0137] In equation (24), the last term on the right-hand side represents a very small perturbation around coordinate t. The cube and higher powers of redundant terms. When the disturbance When it approaches 0, the wave field The higher-order derivative with respect to time t is a quantity that approaches 0 and can be ignored, resulting in the following equation:

[0138] (25);

[0139] Tron The approximate formula for the second-order partial differential of time t is:

[0140] (26);

[0141] In the formula: It is the step size in the time direction. It is the truncation error of the second-order difference.

[0142] For the left side of equation (11), the traditional pseudospectral method uses a second-order central difference to solve for the time derivative, and its solution formula can be written from equation (26):

[0143] (27);

[0144] Will Recorded as , Recorded as , Recorded as .

[0145] Then equation (11) changes as follows:

[0146] (28);

[0147] Similarly, the calculation formula for equation (17) is:

[0148]

[0149] (29).

[0150] Furthermore, the method also includes:

[0151] A source term containing a volume force source term or a moment tensor source term is introduced into the viscoelastic scalar P-wave equation;

[0152] Based on the mathematical form of the source term in the wavenumber-frequency domain, the control law of different source loading methods on the wavefield radiation pattern and amplitude distribution is analyzed.

[0153] Specifically, the implementation process of this embodiment includes:

[0154] In the viscoelastic scalar P equation (17), the source term is introduced as follows:

[0155] (30);

[0156] in, As a source of physical strength, For moment density tensor, An operator consisting of medium parameters and differential operators:

[0157] (31);

[0158] (32);

[0159] The frequency domain of the normalization operator is:

[0160] , (33);

[0161] In the formula, the unit vector Indicates the direction of propagation. Therefore:

[0162] (34);

[0163] Right now Only depend With medium parameters, with azimuth angle This is irrelevant, indicating that angular anisotropy is not caused by... Introduction.

[0164] Taking the spacetime Fourier transform of equation (30) yields:

[0165] (35);

[0166] and then:

[0167] (36);

[0168] in, for - Modular Green's function. Equation (30) shows that the angular factor is entirely determined by:

[0169] (37) Confirmed.

[0170] The three loading methods exhibit consistent kinematic wavefronts under the scalar P framework, with the main differences being in angular amplitude and phase symmetry.

[0171] Example 2

[0172] In this embodiment, a computer terminal device is provided, including:

[0173] One or more processors;

[0174] A memory, coupled to the processor, for storing one or more programs;

[0175] When the one or more programs are executed by the one or more processors, the one or more processors implement the steps of the above-described viscoelastic scalar P-wave equation construction and wave field numerical simulation method.

[0176] In this embodiment, a computer-readable storage medium is also provided, on which a computer program is stored. When the computer program is executed by a processor, it implements the steps of the above-described method for constructing the viscoelastic scalar P-wave equation and numerically simulating the wave field.

[0177] Figure 1 This invention presents the scalar P-wave equation for viscoelastic media and a schematic diagram of the wave field for synthesized P-waves in viscoelastic media. Figure 1 (a) is the forward modeling diagram of the scalar P-wave. Figure 1 (b) is the forward modeling diagram of the synthesized P-wave. The computational domain is a 300×300 regular grid with a grid spacing of 10m. The medium parameters are Vp=2000m / s, Vs=1176.5m / s, Qp=20, Qs=20, and the reference frequency is... =10Hz. The source was placed at the center of the model (1.5,1.5)km. The vector wave clearly showed the coexistence and coupled propagation of P / S modes. The corresponding synthetic P wave and scalar P wave were highly consistent in travel time and phase.

[0178] Figure 2 The wave fields are represented from top left to bottom right as representing four scenarios: elastic, amplitude attenuation only, phase distortion, and viscoelastic. Forward modeling is performed for each of these four scenarios, and the z-components of the four vector wave fields are spliced ​​together by quadrant at the same time to form a combined wave field.

[0179] Figure 3To extract a single trace from the same location in the combined wavefield, a comparison reveals that: closing the reciprocal time coupling term and the fractional Laplace term yields an elastic response. Retaining only the reciprocal time coupling term shows significant amplitude decay while the phase remains essentially unchanged. Retaining only the fractional Laplace term exhibits significant dispersion with slight energy change. Retaining both the reciprocal time coupling term and the fractional Laplace term together displays typical viscoelastic characteristics. Therefore, it can be determined that the reciprocal time coupling term dominates the attenuation effect, and the fractional Laplace term dominates the dispersion effect, indicating that the newly constructed equation can decouple the dispersion and attenuation effects.

[0180] Figure 4 To reduce the shear wave velocity to Vs = 500 m / s while keeping other modeling parameters constant, the simulation results of the synthetic P-wave and scalar P-wave fields were compared. Figure 4 (a) shows the simulation results of the synthetic P-wave. Figure 4 (b) shows the simulation results of the scalar P-wave field. It can be seen that the synthesized P-wave exhibits significant numerical artifacts (local ringing and speckled energy), while the scalar P-wave field remains stable in terms of phase continuity and amplitude purity. The scalar P equation does not explicitly evolve the S-mode, avoiding numerical contamination introduced by the short wavelength of the S-wave. Therefore, when the Vp / Vs ratio is large and Vs is very low, scalar P modeling can maintain higher numerical robustness and wavefield purity, providing more reliable input for subsequent amplitude-preserving imaging and parameter inversion.

[0181] Figure 5 These are wavefield diagrams corresponding to three different source loading methods, among which... Figure 5 (a) is when only vertical body forces are applied. , Figure 5 (b) is the simultaneous application of vertical body forces. With longitudinal moment components , Figure 5 (c) is for loading only the longitudinal moment component. To eliminate the influence of differences in medium and spectrum, the same model and wavelet were used uniformly in the experiment. and The time function only loads the vertical direction. At that time, the energy of the upper and lower wave fields is not uniform, and only the longitudinal moment component is loaded. At that time, the energy of the upper and lower wave fields is uniform.

[0182] Figure 6 These are single-channel wavefields corresponding to three different source loading methods, among which Figure 6 (a) is when only vertical body forces are applied. , Figure 6 (b) is the simultaneous application of vertical body forces. With longitudinal moment components , Figure 6 (c) is for loading only the longitudinal moment component. This shows the magnitude of the upper and lower wave fields.

[0183] Figure 7 The dominant effect of different source loading methods on wavefield directivity, among which Figure 7 (a) is when only vertical body forces are applied. hour, The main lobe is strongest vertically, with the horizontal direction being the node, and the upper and lower sides having similar amplitudes but opposite phases. Figure 7 (b) is the simultaneous application of vertical body forces. With longitudinal moment components hour, . and The coherent superposition of components leads to asymmetry, with one side contributing to the growth while the other side cancels it out. For amplitude weighting. Figure 7 (c) is for loading only the longitudinal moment component. hour, It exhibits an even function distribution, with vertical in-phase enhancement and horizontal minimum, and is symmetrical in the ring zone.

[0184] This invention provides a viscoelastic scalar P-wave equation construction and wavefield numerical simulation method. By considering the viscous characteristics of the subsurface medium and introducing five key physical parameters, it can more precisely characterize the seismic wave propagation process in non-uniform, strongly attenuated media. The wavefield constructed by this method maintains a high degree of consistency with the synthetic P-wave obtained from the viscoelastic full-wave equation in terms of kinematic and dynamic characteristics while preserving amplitude. Furthermore, the derived equations successfully decouple amplitude attenuation and phase dispersion effects at the operator level. Compared to existing technologies, this invention retains the key physical influence of viscoelasticity on P-wave amplitude and phase while significantly reducing computational redundancy and complexity by avoiding explicit solutions for S-waves, thus achieving more accurate dynamic simulation capabilities. This scheme provides a physical modeling foundation for high-fidelity seismic imaging in strongly attenuated media and can be effectively integrated with quality factor compensated reverse-time migration procedures.

[0185] The above are merely preferred embodiments of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for constructing viscoelastic scalar P-wave equations and numerically simulating wave fields, characterized in that, include: Establish the dispersion relationship of viscoelastic media based on the constant Q model; Based on the aforementioned dispersion relation of the viscoelastic medium, a viscoelastic multi-component wave equation containing a coupling term between a fractional-order spatial differential operator and a time derivative is constructed. By introducing a dual projection operator to decouple the P-wave and S-wave modes of the viscoelastic multi-component wave equation, a viscoelastic scalar P-wave equation describing only the dynamic characteristics of the P-wave is obtained. The process of decoupling the P-wave and S-wave modes includes: Construct a normalized gradient operator and its corresponding divergence operator, wherein the normalized gradient operator and the divergence operator satisfy the relationship that they are inverse operators of each other; The displacement field is projected using the divergence operator to define a scalar P-wave field. The dual projection formed by the normalized gradient operator and the divergence operator is applied to both sides of the viscoelastic multi-component wave equation to eliminate the modes related to the S-wave and derive the viscoelastic scalar P-wave equation. The viscoelastic scalar P-wave equations were numerically discretized using the pseudospectral method and the finite difference method to achieve forward modeling of the seismic wave field.

2. The method according to claim 1, characterized in that, The process of establishing the dispersion relation of viscoelastic media based on the constant Q model includes: By treating the quality factor of the medium as a constant within the seismic exploration frequency band, we derive the complex wave number expression characterizing dispersion and attenuation. The complex wavenumber expression is truncated and Taylorized within the target frequency band, transforming the frequency-wavenumber relationship into an expression composed of integer powers of the wavenumber and powers of the frequency. By utilizing Fourier duality, the combinatorial expression is mapped into a constant-order partial differential equation in the time-space domain.

3. The method according to claim 2, characterized in that, The process of constructing the viscoelastic multi-component wave equation includes: Complex moduli related to frequency and wavenumber are introduced for P-waves and S-waves, respectively, and the complex moduli are determined by the dispersion relation of the viscoelastic medium. Based on the differential properties of the Fourier transform, the complex modulus is transformed into the time-space domain to obtain a mechanical modulus expression that includes the fractional Laplace operator and the time partial derivative. Based on the stress-strain constitutive relation, strain-displacement relation, and momentum conservation equation, the viscoelastic multi-component wave equation with displacement field as unknown is obtained by combining the equations.

4. The method according to claim 1, characterized in that, The process of numerical discretization using the pseudospectral method and the finite difference method includes: In the wavenumber domain, the spatial derivative and fractional Laplace operator in the viscoelastic scalar P-wave equation are calculated using Fourier transform and inverse transform. In the time domain, the time derivative term in the viscoelastic scalar P-wave equation is discretized and approximated using a central difference scheme. The wave field is iteratively propagated by alternately performing spatial derivative calculations in the wavenumber domain and time-step updates in the time domain.

5. The method according to claim 1, characterized in that, The method further includes: A source term containing a volume force source term or a moment tensor source term is introduced into the viscoelastic scalar P-wave equation; Based on the mathematical form of the source term in the wavenumber-frequency domain, the control law of different source loading methods on the wavefield radiation pattern and amplitude distribution is analyzed.

6. The method according to claim 3, characterized in that, The mechanical modulus is: ； in, For mechanical modulus, For density, A coefficient vector related to the medium velocity and Q value. For the Laplace operator, This is the time partial derivative.

7. The method according to claim 1, characterized in that, The definitions of the normalized gradient operator and the divergence operator are as follows: ； ； in, For the normalized gradient operator, For space directional partial derivative, Let z be the partial derivative in the z-direction of space. It is a divergence operator.

8. A computer terminal device, characterized in that, include: One or more processors; A memory, coupled to the processor, for storing one or more programs; When the one or more programs are executed by the one or more processors, the one or more processors perform the steps of the method as described in any one of claims 1-7.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the steps of the method as described in any one of claims 1-7.

Images