Seismic wave numerical simulation method based on k-space compensation and related device

Abstract

The invention discloses a seismic wave numerical simulation method based on k-space compensation and a related device, and belongs to the technical field of oil exploration. The method comprises the following steps: constructing a first-order fractional order Laplace viscous acoustic wave equation of a seismic wave field, and deducing a viscous acoustic wave equation only containing velocity field variables according to the first-order fractional order Laplace viscous acoustic wave equation; converting the viscous acoustic wave equation only containing the velocity field variable into a viscous acoustic wave equation only containing the displacement field variable; deducing a recursion formula of the displacement field according to the viscous acoustic wave equation only containing the displacement field variable, and further obtaining a k-space compensation-based wave field recursion formula according to the recursion formula of the displacement field; and obtaining seismic wave field parameters and underground medium parameters, and inputting the seismic wave field parameters and the underground medium parameters into a k-space compensation-based wave field recursion formula to obtain a seismic wave numerical simulation result. According to the method, the time accumulation error can be effectively suppressed, so that the method is suitable for a larger time step length; and excessive extra calculation amount is not introduced, so that the method is more efficient on the premise of ensuring the precision.

Description

A numerical simulation method and related apparatus for seismic waves based on k-space compensation Technical Field

[0001] This invention belongs to the field of petroleum exploration technology and relates to a numerical simulation method and related apparatus for seismic waves based on k-space compensation. Background Technology

[0002] Due to the inherent viscosity of the strata, seismic waves experience amplitude loss, phase distortion, and reduction in high-frequency components during propagation in real subsurface media. These attenuation effects are prevalent in laboratory measurements (Wuenschel, 1965; Wang et al., 2007) and actual exploration (McDononal et al., 1958). To accurately simulate wave propagation underground, attenuation and dissipation must be considered during wavefield extrapolation. Among various models incorporating attenuation effects (Aki and Richards, 1980; Day and Minster, 1984; Carcione et al., 1988; Blanch et al., 1995; Carcione et al., 2002; Zhu and Harris, 2014), the decoupled fractional Laplace wave equation (DFL equation) (Zhu and Carcione, 2014; Zhu and Harris, 2014) has received widespread attention. Compared to other viscous wave equations, the DFL equation has three main advantages. First, the DFL equation accurately characterizes the frequency-independent Q. Second, the decoupled Laplace operator separates phase dispersion and amplitude attenuation effects, which is beneficial for attenuation-compensated reverse time migration (RTM) (Zhu et al., 2014; Li et al., 2016; Wang et al., 2017; Zhu and Sun, 2017) and full waveform inversion (FWI) (Chen and Zhou, 2017; Xue et al., 2017). Third, unlike fractional-order time derivatives, fractional-order Laplace operators are calculated using Fast Fourier Transform (FFT), avoiding the need to store previous wavefields and improving computational efficiency (Chen and Holm, 2004; Carcione, 2010; Treeby and Cox, 2010).

[0003] However, there are still some problems in the numerical solution of this equation, which can be summarized in two main aspects. The first difficulty is that the fractional derivative it contains is a function of both wavenumber and spatial variation, and efficient computation of such mixed-domain operators has always been a mathematical challenge. Several strategies have been proposed to address the mixed-domain problem. The first category is to approximate the fractional Laplace operator using numerical methods, such as the low-rank decomposition method (Sun et al., 2015; Chen et al., 2016), the Hermitian approximation (Yao et al., 2017), the phase-shifting interpolation method (Ji et al., 2017), the matrix transformation method (Chen et al., 2019), and the domain decomposition method (Xue et al., 2018). The second method constructs constant fractional-order DFL viscous acoustic wave equations (Chen et al., 2016, 2019; Wang et al., 2018; Xing and Zhu, 2018; Yang and Zhu, 2018), which naturally avoids the mixed domain problem.

[0004] Secondly, while the pseudospectral method commonly used to solve this equation can guarantee sufficient spatial accuracy, the equation also contains time derivatives related to the attenuation term, resulting in a first-order time discretization accuracy for the entire wave equation. For long-term or long-distance simulations, the accumulated error will gradually amplify. Although using a small time sampling step can suppress time dispersion, it will significantly increase the computational cost. To address this, scholars have proposed numerous methods to suppress time accumulation errors, such as the Lax-Wendroff method (Lax and Wendroff, 1964; Dablain, 1986; Chen, 2011), the rapid expansion method (Pestana and Stoffa, 2010; Tessmer, 2011), the one-step method (Zhang and Zhang, 2009; Sun et al., 2015), the k-space method (Bojarski, 1982; Mast et al., 2001), the Fourier difference method (Song and Fomel, 2011; Song et al., 2013), and the pseudo-analysis method (Etgen and Brandsberg-Dahl, 2009; Crawley et al., 2010). Although these methods also introduce a certain degree of time accumulation error, their temporal accuracy is far superior to the traditional difference-pseudospectral method (FD-PS method).

[0005] The k-space method has been widely used in simulating seismic wave propagation in recent years due to its high simulation accuracy and applicability to mixed domain problems. Compared with the k-space compensation scheme for the second-order wave equation (Soubaras and Zhang, 2008; Liu et al., 2014; Zhou et al., 2021), the k-space compensation scheme for the first-order wave equation (Tabei et al., 2002; Firouzi et al., 2012; Song et al., 2012; Chen et al., 2016) is easier to combine with variable density cases and perfectly matched layer absorbing boundary conditions (Roden and Gedney, 2000; Berenger, 2007; Komatitsch and Martin, 2007). In addition, when using the staggered mesh method (Fornberg, 1987, 1990; Chen, 1996; Ozdenvar and McMechan, 1996), it can effectively alleviate the artifacts and stability problems caused by regular mesh (non-staggered mesh) methods. While the k-space compensation scheme for the first-order wave equation, using the staggered pseudospectral method, allows for larger time steps and achieves better accuracy than the standard pseudospectral method for homogeneous media, the computational cost becomes extremely high for highly heterogeneous media. This is because solving the mixing domain operators using the low-rank decomposition method requires additional Fourier transforms. Furthermore, the computational cost increases further when there are too many mixing domain operators in the first-order wave equation. Summary of the Invention

[0006] The purpose of this invention is to provide a numerical simulation method and related apparatus for seismic waves based on k-space compensation, so as to solve the technical problems of large computational load and low calculation accuracy in the existing technology for highly heterogeneous media.

[0007] To achieve the above objectives, the present invention employs the following technical solution:

[0008] In a first aspect, the present invention provides a numerical simulation method for seismic waves based on k-space compensation, comprising the following steps:

[0009] Construct the first-order fractional Laplace viscous acoustic wave equation for the seismic wave field.

[0010] Based on the first-order fractional Laplace viscous acoustic wave equation, derive the viscous acoustic wave equation containing only velocity field variables.

[0011] Transform the viscous acoustic wave equation containing only velocity field variables into a viscous acoustic wave equation containing only displacement field variables;

[0012] Based on the viscous acoustic wave equation containing only displacement field variables, the recursive formula of the displacement field is derived, and then the wave field recursive formula based on k-space compensation is obtained by transforming the recursive formula of the displacement field.

[0013] Seismic wavefield parameters and subsurface medium parameters are obtained and input into the wavefield recursive formula based on k-space compensation to obtain the numerical simulation results of seismic waves.

[0014] Furthermore, the expression for the first-order fractional Laplace viscous acoustic wave equation of the seismic wave field is as follows:

[0015]

[0016] in,

[0017]

[0018] In the formula, v x v represents the particle velocity component in the x-direction. z Let be the component of the particle velocity in the z-direction; ρ represent the density of the medium; and c represent the velocity of the medium. c0 represents the velocity of the medium at a reference frequency of ω0; γ is a variable, calculated using the following formula:

[0019]

[0020] In the formula, Q is the quality factor of the medium.

[0021] Furthermore, the step of deriving the viscous acoustic wave equation containing only velocity field variables based on the first-order fractional Laplace viscous acoustic wave equation specifically includes:

[0022] Transforming equation (1) to the wavenumber-frequency domain, we get:

[0023]

[0024] In the formula, and Corresponding to P and v respectively x and v z Wavenumber-frequency domain form; operator and Corresponding operators and Fourier transform form;

[0025] The wavenumber-frequency domain form of operator H is:

[0026]

[0027] In the formula, Indicates wave number;

[0028] By eliminating the pressure term in equation (4) through algebraic operations, the viscous acoustic wave equation containing only velocity field variables is obtained as follows:

[0029]

[0030] In the formula, ρ represents the density of the medium.

[0031] Furthermore, the step of transforming the viscous acoustic wave equation containing only velocity field variables into a viscous acoustic wave equation containing only displacement field variables specifically includes:

[0032] Based on the relationship between particle velocity and displacement, equation (6) can be written in matrix form:

[0033]

[0034] Applying the eigenvalue decomposition method, equation (8) can be rewritten as:

[0035]

[0036] For equation (9) to hold, the eigenvalues ​​of the coefficient matrix in equation (8) must be 0. Therefore:

[0037]

[0038] Combine equation (6) and Substituting into equation (10), we get:

[0039]

[0040] The solution is:

[0041]

[0042] in,

[0043]

[0044] Since equations (11) and (12) are equivalent, the viscous acoustic wave equation containing only displacement field variables is obtained from equation (9):

[0045]

[0046] Furthermore, the step of deriving the recursive formula for the displacement field from the viscous acoustic wave equation containing only displacement field variables specifically includes:

[0047] Based on the properties of matrix operations, equation (14) can be split into:

[0048]

[0049] The expression for equation (15) transformed into the wavenumber-time domain is:

[0050]

[0051] In the formula, and This represents the displacement of the particle in the wavenumber-time domain;

[0052] The analytical solution to equation (16) is obtained from the properties of differential equations:

[0053]

[0054] From equation (17), the recursive relationship of the particle's displacement function at times t+δt, t-δt, and t is:

[0055]

[0056] In the formula, and Represents the displacement of a particle in the wavenumber domain;

[0057] Furthermore, the expressions for matrices A and R are:

[0058]

[0059] Furthermore, the step of transforming the displacement field recursive formula to obtain the wave field recursive formula based on k-space compensation specifically includes:

[0060] The matrix A in the recursive formula of the displacement field is split into... in

[0061]

[0062] Then equation (18) can be written as:

[0063]

[0064] set up and Substituting this into equation (21), we obtain the first-order differential equation:

[0065]

[0066] and Let them contain 2 and 4 elements respectively.

[0067]

[0068] Substituting equations (20), (23), and matrix R into equation (22) and transforming it to the spatial domain, we finally obtain the wave field recursive formula based on k-space compensation.

[0069] Furthermore, the expression for the wavefield recursive formula based on k-space compensation is as follows:

[0070]

[0071] Among them, x1=(x,z), x2=(x+δx / 2,z), x3=(x,z+δz / 2);

[0072] Operator L 11 ,L 12 ,L 21 ,L 22 The expressions are as follows:

[0073]

[0074] In the formula, F is the Fourier transform operator; F -1 It is the inverse Fourier transform operator.

[0075] Secondly, the present invention provides a seismic wave numerical simulation system based on k-space compensation, comprising:

[0076] The equation construction module is used to construct the first-order fractional Laplace viscous acoustic wave equations for the seismic wavefield.

[0077] The derivation module is used to derive the viscous acoustic wave equation containing only velocity field variables based on the first-order fractional Laplace viscous acoustic wave equation.

[0078] The first transformation module is used to transform the viscous acoustic wave equation containing only velocity field variables into the viscous acoustic wave equation containing only displacement field variables.

[0079] The second transformation module is used to derive the recursive formula of the displacement field based on the viscous acoustic wave equation containing only displacement field variables, and then transform the recursive formula of the displacement field to obtain the wave field recursive formula based on k-space compensation.

[0080] The simulation module is used to acquire seismic wavefield parameters and subsurface medium parameters, and input them into the wavefield recursive formula based on k-space compensation to obtain the numerical simulation results of seismic waves.

[0081] Thirdly, the present invention provides a computer device including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps of the method described above.

[0082] Fourthly, the present invention provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the steps of the method described above.

[0083] Compared with the prior art, the present invention has the following beneficial effects:

[0084] This invention discloses a numerical simulation method and related apparatus for seismic waves based on k-space compensation. First, a first-order fractional Laplace viscous acoustic wave equation for the seismic wavefield is constructed. Then, a viscous acoustic wave equation containing only velocity field variables is derived. Next, this velocity-only equation is transformed into a displacement-only equation, and a recursive formula for the displacement field is derived. Finally, a wavefield recursive formula based on k-space compensation is obtained for seismic wave numerical simulation. The k-space compensation scheme for the first-order viscous acoustic wave equation of this invention offers better time accuracy and more relaxed stability conditions, allowing for high-precision wavefield extrapolation with larger time steps. Low-rank decomposition addresses the mixed domain problem introduced by the fractional Laplace operator, reducing the number of computationally time-consuming mixed domain operators and significantly decreasing the computational load, thereby improving overall computational efficiency. This provides strong technical support for fields such as seismic exploration and oil and gas resource exploration. Attached Figure Description

[0085] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation on the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0086] Figure 1 is a flowchart of the method of the present invention;

[0087] Figure 2 is a schematic diagram of the system of the present invention;

[0088] Figure 3 shows a comparison between the solutions obtained by the k-space method in the embodiment of the present invention and the analytical solutions. In Figure 3(a), Q = 10; in Figure 3(b), Q = 20; in Figure 3(c), Q = 50; and in Figure 3(d), Q = 100.

[0089] Figure 4 shows wave field snapshots of the FD-PS method and the k-space method in a homogeneous medium according to the embodiments of the present invention. Figure 4(a) shows the FD-PS method, and Figure 4(b) shows the k-space method.

[0090] Figure 5(a) shows the BP velocity model of an embodiment of the present invention, and Figure 5(b) shows the Q model;

[0091] Figure 6 shows the wavefield snapshots generated by the FD-PS method and the k-space method in the embodiments of the present invention. Figure 6(a) is the reference simulation result, Figure 6(b) is the k-space method simulation result with a time step of 2ms, Figure 6(c) is the FD-PS method simulation result with a time step of 1ms, and Figure 6(d) is the FD-PS method simulation result with a time step of 2ms.

[0092] Figure 7 shows the seismic records generated by the FD-PS method and the k-space method in the embodiments of the present invention and their differences with the reference solution. Figure 7(a) shows the seismic records and errors of the k-space method with a time step of 2ms, and Figure 7(b) shows the seismic records and errors of the FD-PS method with a time step of 1ms.

[0093] Figure 8 shows a comparison between the seismic traces generated by different numerical methods in the embodiments of the present invention and the reference solution. Figure 8(a) shows the seismic trace generated by the FD-PS method with a time step of 2ms, Figure 8(b) shows the seismic trace generated by the FD-PS method with a time step of 1ms, Figure 8(c) shows the seismic trace generated by the k-space method with a time step of 2ms, and Figure 8(d) shows the seismic trace generated by the k-space method with a time step of 1ms.

[0094] Figure 9 shows the CPU runtime of the FD-PS method and the k-space method in the embodiments of the present invention;

[0095] Figure 10 is a schematic diagram of the computer device structure of the present invention. Detailed Implementation

[0096] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, 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. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.

[0097] Therefore, the following detailed description of the embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0098] It should be noted that similar labels and letters in the following figures indicate similar items. Therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures.

[0099] In the description of the embodiments of the present invention, it should be noted that if terms such as "upper," "lower," "horizontal," or "inner" indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings, or the orientation or positional relationship commonly used when the product of the invention is in use, they are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation, and therefore should not be construed as a limitation of the present invention. Furthermore, terms such as "first" and "second" are only used to distinguish descriptions and should not be construed as indicating or implying relative importance.

[0100] Furthermore, the use of the term "horizontal" does not imply that the component must be absolutely horizontal, but rather that it can be slightly tilted. For example, "horizontal" simply means that its direction is more horizontal than "vertical," and does not mean that the structure must be completely horizontal, but can be slightly tilted.

[0101] In the description of the embodiments of the present invention, it should also be noted that, unless otherwise explicitly specified and limited, the terms "set," "install," "connect," and "link" should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral connection; 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; and they can refer to the internal connection of two components. Those skilled in the art can understand the specific meaning of the above terms in the present invention according to the specific circumstances.

[0102] The present invention will now be described in further detail with reference to the accompanying drawings:

[0103] Referring to Figure 1, this embodiment of the invention discloses a numerical simulation method for seismic waves based on k-space compensation. By obtaining seismic wavefield parameters and subsurface medium parameters and inputting them into the wavefield recursive formula based on k-space compensation, the numerical simulation results of seismic waves can be obtained. This embodiment proposes a new k-space compensation scheme for the first-order viscous acoustic wave equation. This method can not only effectively suppress time accumulation errors, thus being applicable to larger time steps, but also avoids introducing excessive additional computation, making it more efficient while ensuring accuracy. The derivation steps of the wavefield recursive formula based on k-space compensation are as follows:

[0104] First, a first-order fractional Laplace viscous acoustic wave equation for the seismic wave field is constructed.

[0105] The constant fractional-order 2D velocity-pressure viscous acoustic wave equation can be expressed as:

[0106]

[0107] in

[0108]

[0109] P represents the pressure of the wave field, v x and v z Let be the components of the particle velocity in the x and z directions, respectively, where ρ represents the density of the medium and c represents the velocity of the medium. c0 represents the velocity of the medium at a reference frequency ω0, and the variable γ is related to Q, with the expression being:

[0110]

[0111] Where Q is the quality factor of the medium.

[0112] Secondly, the matrix form of the viscous acoustic wave equation containing only velocity field variables is derived.

[0113] Transforming equation 1 to the wavenumber-frequency domain, we get:

[0114]

[0115] in and Corresponding to P and v respectively x and v z Wavenumber-frequency domain form; operator and Corresponding operators and Fourier transform form; operator

[0116]

[0117] This is the wavenumber-frequency domain form of operator H; Indicates wave number.

[0118] By performing algebraic operations, we can eliminate the pressure term in equation 4 and obtain... and relation

[0119]

[0120] Third, it is transformed into a viscous acoustic wave equation containing only displacement field variables.

[0121] Because the particle's velocity and displacement have the following relationship

[0122]

[0123] Therefore, Equation 6 can be written in the following matrix form.

[0124]

[0125] Applying the eigenvalue decomposition method, Equation 8 is rewritten as follows:

[0126]

[0127] To make equation 9 true, let the eigenvalues ​​of the coefficient matrix of equation 8 be 0, that is...

[0128]

[0129] Six styles and Substituting into equation 10, we get

[0130]

[0131] The solution is

[0132]

[0133] in

[0134]

[0135] Since Form 11 and Form 12 are equivalent, Form 9 can be written as

[0136]

[0137] Fourth, derive the recursive formula for the displacement field.

[0138] According to equation 14, it can be decomposed into...

[0139]

[0140] Transforming equation 15 to the wavenumber-time domain, its expression is:

[0141]

[0142] in and This represents the displacement of the particle in the wavenumber-time domain. By the properties of differential equations, the analytical solution to equation 16 is obtained as follows:

[0143]

[0144] From Equation 17, the recursive relationship between the particle's displacement function at times t+δt, t-δt, and t is:

[0145]

[0146] in and Represents the displacement of a particle in the wavenumber domain, and has

[0147]

[0148] Fifth, a recursive formula for the wave field based on k-space compensation is obtained.

[0149] Matrix A is split into its components through matrix decomposition. in

[0150]

[0151] At this point, Form 18 can be written as

[0152]

[0153] If we assume and Substituting these equations into equation 21, we obtain the first-order differential equation.

[0154]

[0155] As can be seen from the above definition, and They contain 2 and 4 elements respectively, therefore let

[0156]

[0157] It is important to note here that no approximations were made during the derivation of Equation 22, therefore Equation 22 is the exact solution of the viscous wave equation (Equation 1) in the wavenumber domain. Specifically, as Q → +∞, L... m =sinc(kc0δt / 2) and R=O, and equation 22 degenerates into an exact solution of the attenuationless acoustic equation in the wavenumber domain.

[0158] Substituting equations 20, 23, and matrix R into equation 22, and transforming to the spatial domain, we obtain the wavefield recursive formula for k-space compensation of the first-order viscous acoustic wave equation.

[0159]

[0160] Where x1 = (x, z), x2 = (x + δx / 2, z), x3 = (x, z + δz / 2), the operator L 11 ,L 12 ,L 21 ,L 22 The expression is:

[0161]

[0162] Referring to Figure 2, this embodiment of the invention discloses a seismic wave numerical simulation system based on k-space compensation, characterized in that it includes an equation construction module, a derivation module, a first transformation module, a second transformation module, and a simulation module.

[0163] It should be noted that the equation construction module is used to construct the first-order fractional Laplace viscous acoustic wave equation for the seismic wavefield; the derivation module is used to derive the viscous acoustic wave equation containing only velocity field variables based on the first-order fractional Laplace viscous acoustic wave equation; the first transformation module is used to transform the viscous acoustic wave equation containing only velocity field variables into the viscous acoustic wave equation containing only displacement field variables; the second transformation module is used to derive the recursive formula of the displacement field based on the viscous acoustic wave equation containing only displacement field variables, and then transform the recursive formula of the displacement field to obtain the wavefield recursive formula based on k-space compensation; the simulation module is used to obtain the seismic wavefield parameters and subsurface medium parameters, and input them into the wavefield recursive formula based on k-space compensation to obtain the seismic wave numerical simulation results.

[0164] Example:

[0165] First, the accuracy of the proposed method was verified by testing a homogeneous model. The homogeneous medium model used had a grid size of 300×300 points, a spatial sampling interval of 20m, and a temporal sampling rate of 3.0ms. The medium reference velocity was c0 = 2000m / s, and the dominant source frequency was 12.5Hz. The analytical solution was obtained by convolution of the source wavelet and the Green's function. Figure 3 shows the waveforms of the k-space viscous acoustic wave analytical solution under different Q values. It can be seen that the waveforms simulated by k-space under different attenuation levels are in good agreement with the analytical solution. Therefore, the accuracy of k-space is preliminarily verified.

[0166] Figure 4 shows the wave fields simulated in a homogeneous medium by the k-space and FD-PS methods (using second-order finite difference for time and pseudospectral method for space) with different time steps. It can be seen that the wave field of the FD-PS method exhibits increasingly severe forward phase perturbations as the time step increases, indicating that the simulation results are significantly affected by time dispersion. However, the wave field of the k-space method shows no obvious numerical dispersion. This further demonstrates that the k-space method can compensate for the cumulative error introduced by the discrete-time derivative in the viscous acoustic wave equation.

[0167] A complex BP model was tested to verify the effectiveness of the proposed method in compensating for strong attenuation effects. Figure 5 shows the velocity model and the Q model. The velocity model has a grid size of 634 × 451 points, and the Q model is calculated using the empirical formula Q = 3.516(c0 / 1000). 2.2 The transformation yielded a Rayleigh wavelet with a dominant frequency of 20Hz, positioned at (5000m, 4400m).

[0168] Figure 6 shows a wavefield snapshot with a propagation time of 1 s. The reference simulation result is a reference wavefield snapshot generated by the FD-PS method with a time sampling interval of 0.2 ms. For clearer comparison, the area within the box is magnified. In Figure 6d, some obvious wavefront perturbations are observed, indicating that the FD-PS method suffers severe temporal dispersion at a time step of 2 ms. When the time step is reduced to 1 ms, the wavefront perturbations are significantly reduced (see Figure 6c). Comparing Figures 6a-6c, the wavefield snapshots are almost identical, demonstrating that the proposed k-space method is still capable of accurate wavefield simulation even when its time sampling interval is twice that of the FD-PS method.

[0169] Figure 7 shows the differences between the seismic records generated by the FD-PS method and the k-space method and the reference record, with a propagation time of 4 s. Figure 7 shows that the difference between the seismic record generated by the FD-PS method with a time step of 1 ms and the reference record is still significantly greater than the error of the seismic record simulated by the k-space method with a large time step of 2 ms.

[0170] Figure 8 compares the seismic traces within 1.0–3.5 s at an offset of 1000 m. Observing the relative error values ​​reveals a similar conclusion to the above: the seismic traces generated by the FD-PS method with a time step of 2 ms exhibit significant numerical dispersion (see Figure 8a). When the time step is reduced to 1 ms, the relative error decreases significantly (see Figure 8b). Figure 8c shows the calculations obtained by the proposed k-space method with a time step of 2 ms. It can be seen that its accuracy is almost the same as that of the FD-PS method with a time step of 1 ms, and significantly lower than that of the FD-PS method with a time step of 2 ms. Furthermore, when the proposed k-space method uses a 1 ms time step for simulation, the seismic traces it generates (see Figure 8d) are almost completely identical to the reference traces.

[0171] Finally, Figure 9 compares the computational efficiency of the proposed method with that of the FD-PS method. It can be seen that when the time step is 2ms, the CPU runtime of the k-space method is slightly shorter than that of the FD-PS method. However, when the simulation accuracy is comparable, the CPU runtime of the k-space method (time step 2ms) is significantly shorter than that of the FD-PS method (time step 1ms). This is because the k-space method reduces the computationally time-consuming mixed-domain operator, thus improving its computational efficiency.

[0172] In one embodiment of the present invention, referring to FIG10, a computer device is provided. The computer device includes a processor and a memory. The memory stores a computer program, which includes program instructions. The processor executes the program instructions stored in the computer storage medium. The processor may be a Central Processing Unit (CPU), or other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. It is the computing and control core of the terminal, suitable for implementing one or more instructions, specifically suitable for loading and executing one or more instructions in the computer storage medium to achieve a corresponding method flow or corresponding function. The processor described in this embodiment of the present invention can be used for the operation of a seismic wave numerical simulation method based on k-space compensation.

[0173] This invention also provides a storage medium, specifically a computer-readable storage medium (Memory), which is a memory device in a computer device used to store programs and data. It is understood that the computer-readable storage medium here can include both the built-in storage medium in the computer device and extended storage media supported by the computer device. The computer-readable storage medium provides storage space that stores the terminal's operating system. Furthermore, this storage space also stores one or more instructions suitable for loading and execution by a processor. These instructions can be one or more computer programs (including program code). It should be noted that the computer-readable storage medium here can be high-speed RAM or non-volatile memory, such as at least one disk storage device. The processor can load and execute one or more instructions stored in the computer-readable storage medium to implement the corresponding steps of the k-space compensation-based seismic wave numerical simulation method in the above embodiments.

[0174] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0175] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions specified in one or more blocks of the flowchart illustrations and / or one or more blocks of the block diagrams.

[0176] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means that implement the functions specified in one or more flowcharts and / or one or more block diagrams.

[0177] These computer program instructions may also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process, such that the instructions, which execute on the computer or other programmable apparatus, provide steps for implementing the functions specified in one or more flowcharts and / or one or more block diagrams.

[0178] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the specific implementation of the present invention. Any modifications or equivalent substitutions that do not depart from the spirit and scope of the present invention should be covered within the scope of protection of the claims of the present invention.

Claims

1. A numerical simulation method for seismic waves based on k-space compensation, characterized in that, Includes the following steps: A first-order fractional Laplace viscous acoustic wave equation for the seismic wavefield is constructed. Based on this equation, a viscous acoustic wave equation containing only velocity field variables is derived. This equation is then transformed into one containing only displacement field variables. A recursive formula for the displacement field is derived from this equation, and further transformed using this formula to obtain a wavefield recursive formula based on k-space compensation. Seismic wavefield parameters and subsurface medium parameters are obtained and input into the k-space compensated wavefield recursive formula to obtain the numerical simulation results of the seismic waves.

2. The numerical simulation method for seismic waves based on k-space compensation according to claim 1, characterized in that, The expression for the first-order fractional Laplace viscous acoustic wave equation of the seismic wavefield is as follows: in, In the formula, v x v represents the particle velocity component in the x-direction. z Let be the component of the particle velocity in the z-direction; ρ represent the density of the medium; and c represent the velocity of the medium. c0 represents the velocity of the medium at a reference frequency of ω0; γ is a variable, calculated using the following formula: In the formula, Q is the quality factor of the medium.

3. The numerical simulation method for seismic waves based on k-space compensation according to claim 2, characterized in that, The steps for deriving the viscous acoustic wave equation containing only velocity field variables based on the first-order fractional Laplace viscous acoustic wave equation specifically include: transforming equation (1) to the wavenumber-frequency domain, resulting in: In the formula, and Corresponding to P and v respectively x and v z Wavenumber-frequency domain form; operator and Corresponding operators and The Fourier transform form; the wavenumber-frequency domain form of operator H is: In the formula, The wave number is represented by: Through algebraic operations, the pressure term in equation (4) is eliminated, resulting in the viscous acoustic wave equation containing only velocity field variables: In the formula, ρ represents the density of the medium.

4. The numerical simulation method for seismic waves based on k-space compensation according to claim 3, characterized in that, The step of transforming the viscous acoustic wave equation containing only velocity field variables into a viscous acoustic wave equation containing only displacement field variables specifically includes: rewriting equation (6) in matrix form based on the relationship between the particle's velocity and displacement: Applying the eigenvalue decomposition method, equation (8) can be rewritten as: For equation (9) to hold, the eigenvalues ​​of the coefficient matrix in equation (8) must be 0. Therefore: Combine equation (6) and Substituting into equation (10), we get: The solution is: in, Since equations (11) and (12) are equivalent, the viscous acoustic wave equation containing only displacement field variables is obtained from equation (9):

5. The numerical simulation method for seismic waves based on k-space compensation according to claim 4, characterized in that, The steps for deriving the recursive formula for the displacement field based on the viscous acoustic wave equation containing only displacement field variables specifically include: splitting equation (14) into the following parts according to the properties of matrix operations: The expression for equation (15) transformed into the wavenumber-time domain is: In the formula, and This represents the displacement of the particle in the wavenumber-time domain; the analytical solution of equation (16) is obtained from the properties of the differential equation: From equation (17), the recursive relationship of the particle's displacement function at times t+δt, t-δt, and t is: In the formula, and Let A represent the displacement of a particle in the wavenumber domain; and let the expressions for matrices A and R be:

6. The numerical simulation method for seismic waves based on k-space compensation according to claim 5, characterized in that, The step of transforming the recursive formula of the displacement field to obtain the wave field recursive formula based on k-space compensation specifically includes: splitting matrix A in the recursive formula of the displacement field into... in Then equation (18) can be written as: set up and Substituting this into equation (21), we obtain the first-order differential equation: and Let them contain 2 and 4 elements respectively. Substituting equations (20), (23), and matrix R into equation (22) and transforming it to the spatial domain, we finally obtain the wave field recursive formula based on k-space compensation.

7. The numerical simulation method for seismic waves based on k-space compensation according to claim 1, characterized in that, The expression for the wavefield recursive formula based on k-space compensation is as follows: Where x1 = (x, z), x2 = (x + δx / 2, z), x3 = (x, z + δz / 2); operators The expressions are as follows: In the formula, F is the Fourier transform operator; F -1 It is the inverse Fourier transform operator.

8. A numerical simulation system for seismic waves based on k-space compensation, characterized in that, include: The system comprises the following modules: an equation construction module for constructing the first-order fractional Laplace viscous acoustic wave equations for the seismic wavefield; a derivation module for deriving viscous acoustic wave equations containing only velocity field variables from the first-order fractional Laplace viscous acoustic wave equations; a first-order transformation module for transforming the viscous acoustic wave equations containing only velocity field variables into viscous acoustic wave equations containing only displacement field variables; a second-order transformation module for deriving the recursive formula for the displacement field from the viscous acoustic wave equations containing only displacement field variables, and then transforming the recursive formula for the displacement field to obtain the wavefield recursive formula based on k-space compensation; and a simulation module for acquiring seismic wavefield parameters and subsurface medium parameters, inputting them into the wavefield recursive formula based on k-space compensation, and obtaining the numerical simulation results of the seismic waves.

9. A computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the computer program, it implements the steps of the method as described in any one of claims 1-7.

10. A computer-readable storage medium storing a computer program, 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