Numerical simulation method and device for wave equation of viscous-acoustic medium on undulating surface

Abstract

The invention discloses an undulating surface viscoacoustic medium wave equation numerical simulation method and device, and the method comprises the steps: building a viscoacoustic medium wave equation according to the position of a seismic source, stratum longitudinal wave speed data, stratum density data, stratum quality factor data, a seismic source wavelet function and a surface elevation; establishing a mapping relation expression between the physical domain data and the computational domain data; determining a viscoacoustic medium wave equation in a computational domain based on the mapping relation expression and the viscoacoustic medium wave equation; determining a wave field of the viscous-acoustic medium wave equation in the computational domain by using a fully staggered grid high-order finite difference method; and mapping the wave field of the computational domain to obtain a wave field of a physical domain.

Description

Technical Field

[0001] This article relates to the field of geophysical exploration technology, and particularly to a numerical simulation method and device for wave equations in a visco-acoustic medium with a rugged surface. Background Art

[0002] Numerical simulation of wave equations is a seismic wave numerical simulation method that directly solves wave equations on a discrete grid using numerical algorithms for seismic wave simulation. Among them, the finite difference method has the advantage of low computational cost and is widely used in forward modeling, reverse time migration, and full waveform inversion. However, in exploration areas with severe surface undulations and where seismic wave attenuation needs to be considered, traditional finite difference methods cannot adapt to the rugged surface and cannot efficiently and accurately simulate the propagation of wave fields in visco-acoustic media, making it difficult to apply in actual production. Therefore, there is an urgent need to propose a numerical simulation method that is suitable for a rugged surface and can efficiently and accurately simulate the propagation of wave fields in visco-acoustic media to solve the problems in actual production. Summary of the Invention

[0003] The inventors found that:

[0004] The numerical simulation method of wave equations in visco-acoustic media can be divided into frequency domain methods and time domain methods. Frequency domain methods can be directly solved, but they require Fourier transforms, resulting in high computational costs. The commonly used classical standard linear solid equation can be directly solved in the time domain, but in the equation, Q needs to be converted into multiple stress-strain relaxation times, which requires a high amount of memory and has an unclear physical meaning. To solve the problem of a rugged surface, there are many methods, and among them, the coordinate transformation method is relatively effective. Compared with using regularly encrypted grids in the Cartesian coordinate system to fit the rugged surface, the coordinate transformation method transforms the grid of the physical domain with undulations into a horizontal computational domain, with lower computational costs. And the coordinate transformation method can easily incorporate free surface conditions. In numerical simulation algorithms, to balance computational accuracy and efficiency, time domain staggered grid finite difference is a better choice. However, due to the need for vertical stretching in the coordinate transformation method, the conventional staggered grid finite difference method cannot accurately calculate partial derivatives.

[0005] To address the above problems, a wave equation for visco-acoustic media based on an explicit expression of the quality factor is provided, which has a low memory requirement and a clear physical meaning for the equation; for the problem of a rugged surface, a conversion method for a rugged surface based on an analytical elevation mapping relationship is implemented; for the fact that the conventional staggered grid finite difference method cannot accurately calculate partial derivatives, a full staggered grid is used for high-order finite difference forward modeling.

[0006] In a first aspect, the present application provides a numerical simulation method for the viscoacoustic medium wave equation on a rough surface, including: establishing a viscoacoustic medium wave equation according to the position of the seismic source, formation P-wave velocity data, formation density data, formation quality factor data, seismic source wavelet function, and surface elevation; establishing a mapping relationship expression between physical domain data and computational domain data; determining the viscoacoustic medium wave equation in the computational domain based on the mapping relationship expression and the viscoacoustic medium wave equation; using the fully staggered grid high-order finite difference method to determine the wave field of the viscoacoustic medium wave equation in the computational domain in the computational domain; and mapping the wave field of the computational domain to obtain the wave field of the physical domain.

[0007] In a second aspect, an embodiment of the present invention further provides a device for numerical simulation of the viscoacoustic medium wave equation on a rough surface. The device includes: a memory and a processor; the memory is used to store a program for numerical simulation of the viscoacoustic medium wave equation on a rough surface, and the processor is used to read and execute the program for numerical simulation of the viscoacoustic medium wave equation on a rough surface, and execute the method described in any one of the above embodiments.

[0008] In a third aspect, an embodiment of the present invention further provides a computer-readable storage medium, on which a data processing program is stored, and the data processing program is executed by a processor to perform the method for numerical simulation of the viscoacoustic medium wave equation on a rough surface described in any one of the above embodiments.

[0009] Compared with the related art, the present application provides a numerical simulation method and device for the viscoacoustic medium wave equation on a rough surface. The method includes: establishing a viscoacoustic medium wave equation according to the position of the seismic source, formation P-wave velocity data, formation density data, formation quality factor data, seismic source wavelet function, and surface elevation; establishing a mapping relationship expression between physical domain data and computational domain data; determining the viscoacoustic medium wave equation in the computational domain based on the mapping relationship expression and the viscoacoustic medium wave equation; using the fully staggered grid high-order finite difference method to determine the wave field of the viscoacoustic medium wave equation in the computational domain in the computational domain; and mapping the wave field of the computational domain to obtain the wave field of the physical domain. In response to the problem of surface undulation, the present application uses a conversion method for a rough surface based on an analytical elevation mapping relationship to convert the viscoacoustic medium wave equation into the viscoacoustic medium wave equation in the computational domain, and performs forward modeling using the fully staggered grid high-order finite difference method. This method reduces the computational cost and obtains more accurate numerical simulation results of the viscoacoustic medium.

[0010] Other features and advantages of the present application will be described in the subsequent specification, and some of them will become obvious from the specification, or be understood by implementing the present application. Other advantages of the present application can be realized and obtained through the solutions described in the specification and the accompanying drawings. Description of the Drawings

[0011] The accompanying drawings are used to provide an understanding of the technical solutions of the present application and form a part of the specification. Together with the embodiments of the present application, they are used to explain the technical solutions of the present application and do not constitute a limitation to the technical solutions of the present application.

[0012] Figure 1 It is a flowchart of the numerical simulation method for the visco-acoustic medium wave equation on the undulating ground surface in the embodiment of the present application;

[0013] Figure 2 It is a schematic diagram of the numerical simulation device for the visco-acoustic medium wave equation on the undulating ground surface in the embodiment of the present application;

[0014] Figure 3 It is a schematic diagram of the velocity field and Q field of the 3D loess tableland model in some exemplary embodiments;

[0015] Figure 4 It is a schematic diagram of the wave field snapshot of the acoustic wave simulation of the 3D loess tableland model in some exemplary embodiments;

[0016] Figure 5 It is a schematic diagram of the shot gather record of the acoustic wave simulation of the 3D loess tableland model in some exemplary embodiments;

[0017] Figure 6 It is a schematic diagram of the wave field snapshot of the visco-acoustic wave simulation of the 3D loess tableland model in some exemplary embodiments;

[0018] Figure 7 It is a schematic diagram of the shot gather record of the visco-acoustic wave simulation of the 3D loess tableland model in some exemplary embodiments. Detailed implementation manners

[0019] The present application describes multiple embodiments, but the description is exemplary rather than restrictive, and it is obvious to those of ordinary skill in the art that there can be more embodiments and implementation solutions within the scope covered by the embodiments described in the present application. Although many possible feature combinations are shown in the accompanying drawings and discussed in the detailed implementation manners, many other combination ways of the disclosed features are also possible. Unless specifically restricted, any feature or element of any embodiment can be combined with any other feature or element in any other embodiment, or can replace any other feature or element in any other embodiment.

[0020] The present application includes and contemplates combinations of features and elements known to those of ordinary skill in the art. The embodiments, features and elements disclosed in the present application may also be combined with any conventional features or elements to form a unique invention scheme defined by the claims. Any features or elements of any embodiment may also be combined with features or elements from other invention schemes to form another unique invention scheme defined by the claims. Therefore, it should be understood that any feature shown and / or discussed in the present application may be implemented individually or in any appropriate combination. Therefore, except for the limitations made according to the attached claims and their equivalents, the embodiments are not subject to other restrictions. In addition, various modifications and changes may be made within the scope of protection of the attached claims.

[0021] In addition, when describing representative embodiments, the specification may have presented the method and / or process as a specific sequence of steps. However, to the extent that the method or process does not rely on the specific order of the steps described herein, the method or process should not be limited to the steps of the specific order described. As will be understood by those of ordinary skill in the art, other sequences of steps are also possible. Therefore, the specific sequence of the steps set forth in the specification should not be interpreted as a limitation to the claims. In addition, the claims for the method and / or process should not be limited to the steps of performing them in the order written, and those skilled in the art can easily understand that these sequences can be changed and still remain within the spirit and scope of the embodiments of the present application.

[0022] The embodiment of the present invention provides a method for numerically simulating a wave equation of a viscoacoustic medium on an undulating surface. Figure 1 As shown, the method includes steps S100-S140:

[0023] S100: Establish a viscoacoustic medium wave equation according to the location of the earthquake source, the formation P-wave velocity data v(x), the formation density data ρ(x), the formation quality factor data Q(x), the earthquake source wavelet function f(t) and the surface elevation Γ(x, y); x represents the physical domain data x=(x, y, z);

[0024] S110: Establishing a mapping relationship expression between physical domain data and computational domain data;

[0025] S120: Determine the viscoacoustic medium wave equation in the calculation domain based on the mapping relationship expression and the viscoacoustic medium wave equation;

[0026] S130: Determine the wave field of the viscoacoustic medium wave equation in the calculation domain by using a fully staggered grid high-order finite difference method;

[0027] S140: Mapping the wave field in the computational domain to obtain the wave field in the physical domain.

[0028] In this embodiment, the data in the physical domain includes: the position of the seismic source x = (x, y, z), the formation P-wave velocity data v(x), the formation density data ρ(x), the formation quality factor data Q(x), the seismic source wavelet function f(t), and the surface elevation Γ(x, y).

[0029] In an exemplary embodiment, the mapping relationship between the physical domain data and the computational domain data is expressed as:

[0030]

[0031] In the above relationship expression, the physical domain data x = (x, y, z), the computational domain data α = (α, β, γ), Γ(x, y) is the surface elevation, z max is the maximum depth of the computational region, γ max = z max + Γ max where Γ max is the maximum elevation.

[0032] Among them, the physical domain data and the computational domain data are in a one-to-one correspondence relationship.

[0033] In an exemplary embodiment, based on the above mapping relationship expression, converting the input data in the physical domain to the computational domain and determining the viscoacoustic medium wave equation in the computational domain includes:

[0034] Step 1. Convert the physical domain data to the computational domain data based on the mapping relationship expression and calculate the partial derivatives in the corresponding coordinate system;

[0035] Step 2. Obtain the viscoacoustic medium wave equation in the computational domain based on the viscoacoustic medium wave equation and the partial derivatives.

[0036] In an exemplary embodiment, the viscoacoustic medium wave equation with explicit expression of the quality factor in the physical domain is:

[0037]

[0038] In the above formula, ρ(x) is the formation density, v(x, t) is the particle velocity, u(x, t) is the stress, Q(x) is the formation quality factor, v u (x) is the formation P-wave velocity, s represents the sum of weights, ∧ (l) (x, t) is the memory variable, x s is the position of the seismic source, H(t) is the step function, f(t) is the seismic source time function; M (l)The weight representing the l-th relaxation mechanism, N represents the number of groups of relaxation mechanisms, l represents the relaxation mechanism number, τ (l) represents the relaxation time of the l-th relaxation mechanism.

[0039] In an exemplary embodiment, the acoustic wave equation in the viscous acoustic medium in the computational domain is:

[0040]

[0041] where is the generalized partial derivative operator, and α s is the position of the seismic source in the computational domain.

[0042] In an exemplary embodiment, the generalized partial derivative operator is:

[0043]

[0044]

[0045] In an exemplary embodiment, the process of using the fully staggered grid high-order finite difference method to determine the wave field of the acoustic wave equation in the viscous acoustic medium in the computational domain is as follows:

[0046] Step 1. Use the fully staggered grid high-order finite difference method to perform differential discretization on the acoustic wave equation in the viscous acoustic medium in the computational domain to obtain second-order in time and high-order in space differences;

[0047] Step 2. Set the time step, spatial interval, and source function that meet the conditions, and perform iteration along the time direction to determine the simulated wave field in the computational domain.

[0048] In an exemplary embodiment, the fully staggered grid finite difference of the acoustic wave equation in the viscous acoustic medium is:

[0049]

[0050] where u n+1 represents the stress at the calculation point (i, j, k) at the n + 1 moment, and u n represents the stress at the calculation point (i, j, k) at the n moment; Δt represents the time step, and S represents the total weight; M α represents the first-order derivative of the finite difference operator along the α direction, and M β represents the first-order derivative of the finite difference operator along the β direction, and M γ represents the first-order derivative of the finite difference operator along the γ direction; is the x component of the particle velocity data at the n + 1 / 2 moment, is the y component of the particle velocity data at the n + 1 / 2 moment, The z - component of the particle velocity data at the time \(n + 1 / 2\), The x - component of the particle velocity data at the time \(n+3 / 2\), The y - component of the particle velocity data at the time \(n + 3 / 2\), The z - component of the particle velocity data at the time \(n + 3 / 2\), ∧ (l)，n The memory variable at the time \(n\); ∧ (l)，n+1 The memory variable at the time \(n + 1\).

[0051] In an exemplary embodiment, the condition satisfied by the time step and the spatial interval is:

[0052]

[0053] where \(c\) is the sum of the difference coefficients, \(v\) max is the maximum velocity, \(Q\) min is the minimum value of the quality factor, \(c\) x，max is the maximum value of the surface coordinate stretching coefficient \(c\) x of \(c\) y，max is the maximum value of the surface coordinate stretching coefficient \(c\) y of \(c\) z，max is the maximum value of the surface coordinate stretching coefficient \(c\) z of \(c\).

[0054] In a second aspect, an embodiment of the present invention further provides a device for numerical simulation of the visco - acoustic medium wave equation on a rugged surface, as Figure 2 shown. The device includes: a memory 200 and a processor 210; the device includes: a memory and a processor; the memory is used to save the program for numerical simulation of the visco - acoustic medium wave equation on a rugged surface, and the processor is used to read and execute the program for numerical simulation of the visco - acoustic medium wave equation on a rugged surface, and execute the method described in any one of the above embodiments.

[0055] In a third aspect, an embodiment of the present invention further provides a computer - readable storage medium, on which a data - processing program is stored, and the data - processing program is executed by a processor to perform the method for numerical simulation of the visco - acoustic medium wave equation on a rugged surface described in any one of the above embodiments.

[0056] The embodiments of the present application have the following technical effects:

[0057] (1) Based on the visco - acoustic medium wave equation with explicit expression of the quality factor, the calculation cost is reduced;

[0058] (2) Based on the conversion method of the rugged surface based on the analytical elevation mapping relationship, the problems brought by the rugged surface are solved;

[0059] (3) Using a full - staggered grid for high - order finite - difference forward simulation to obtain accurate numerical simulation results of the visco - acoustic medium.

[0060] Example 1

[0061] Step 1. Determine the model parameters of the physical domain:

[0062] In this step, the model parameters of the physical domain are the formation longitudinal wave velocity data v(x), the formation density data ρ(x), the formation quality factor data Q(x), the source wavelet function f(t), and the surface elevation Γ(x, y); x represents the physical domain data x = (x, y, z).

[0063] Figure 3 For the velocity field and Q field of the 3D loess tableland model, the model size is 475 * 1301 * 800 grids, and the grid size is 10 m.

[0064] Step 2. Establish the mapping relationship expression between the physical domain data and the computational domain data;

[0065] Map the irregular physical domain x = (x, y, z) to the regular computational domain α = (α, β, γ):

[0066]

[0067] where Γ(x, y) is the surface elevation, and z max is the maximum depth of the computational region;

[0068] γ max = z max + Γ max Γ max is the maximum surface elevation.

[0069] After the above coordinate transformation, the undulating surface in the physical domain becomes a horizontal surface in the computational domain, and the free surface boundary condition can be realized by using a simple mirror method.

[0070] Step 3. Establish the viscoacoustic medium wave equation

[0071] Step 31. In the viscoacoustic medium, establish the constitutive relation equation satisfied by the stress u(x, t) and the velocity v(x, t):

[0072]

[0073] where K(x, t) is the bulk modulus varying with time; x represents the irregular physical domain data x = (x, y, z), t and t′ represent two different moments respectively, and v(x, t′) represents the velocity at the moment t′.

[0074] Step 32. Determine the function bulk modulus

[0075] The function bulk modulus K(x, t) can be represented by N sets of relaxation mechanisms:

[0076]

[0077] Among them, is the relaxed bulk modulus, Q(x) is the quality factor, and H(t) is the step function; v u (x) is the formation P-wave velocity, M (l) represents the weight of the l-th relaxation mechanism, N represents the number of groups of relaxation mechanisms, l represents the relaxation mechanism number, and τ (l) represents the relaxation time of the l-th relaxation mechanism.

[0078] Step 33. Taking the derivative of K(x, t), we can obtain:

[0079]

[0080] Among them, δ(t) is the impulse function, s is the total weight,

[0081] Step 34. Substituting the equation in Step 33 into the constitutive relation equation in Step 31:

[0082]

[0083] Among them,

[0084] Step 35. Combining the equation in Step 34 with Newton's second law, we can obtain the viscoacoustic wave motion equation expressed by Q:

[0085]

[0086] Among them, x s represents the source position, and f(t) is the source time function.

[0087] Step 4. Based on the mapping relationship expression between the physical domain data and the computational domain data, the partial derivatives in the physical domain can be transformed into:

[0088]

[0089] Step 5. Based on the mapping relationship expression, the partial derivatives in Step 4, and the viscoacoustic medium wave equation, determine the viscoacoustic medium wave equation in the computational domain:

[0090]

[0091] Among them, it contains a generalized partial derivative operator:

[0092]

[0093] Among them, c x 、cy and c z is the surface coordinate stretching coefficient.

[0094] Step 6: Set numerical simulation parameters

[0095] Step 61. To ensure the stability of the solution, the time step and the spatial interval need to satisfy the following conditions:

[0096]

[0097] where c is the sum of the difference coefficients, v max is the maximum velocity, Q min is the minimum value of the quality factor, c x，max is the maximum value of the surface coordinate stretching coefficient c x is the maximum value of the surface coordinate stretching coefficient c y，max is the maximum value of the surface coordinate stretching coefficient c y is the maximum value of the surface coordinate stretching coefficient c z，max is the maximum value of the surface coordinate stretching coefficient c z is the maximum value of the surface coordinate stretching coefficient c, Δt is the time step, and Δh is the spatial interval.

[0098] Step 62. Calculation of the surface coordinate stretching coefficients c x , c y and c z :

[0099]

[0100] where z max is the maximum depth of the calculation area, γ max = z max + Γ max and Γ max is the maximum elevation.

[0101] Step 7: Use a fully staggered grid to perform a difference discretization of the viscoacoustic wave equation in Step 5 to obtain its second-order-in-time and high-order-in-space finite-difference format, as follows:

[0102]

[0103] where M i (i = α, β, γ) represents the first derivative of the finite-difference operator; the stress, velocity, and stored wave fields are all calculated iteratively according to the formula.

[0104] In a fully staggered grid, two types of grids at different positions are coupled together. If the source loading is not appropriate, pseudo-waves will be generated. In this embodiment, in addition to loading the source at the source position (i, j, k), a source 0.25 times as strong is also loaded at the auxiliary points (i±1 / 2, j±1 / 2, k), (i±1 / 2, j, k±1 / 2), and (i / 2, j±1 / 2, k±1 / 2). The final seismic record is obtained by weighted superposition of the seismic records of the main source and the auxiliary sources. Taking the calculation of u at the point (i, j, k) as an example n+1 the x-component of the velocity at (i±m / 2, j, k) and (i, j, k±m / 2) is used for calculation and and the y-component of the velocity at (i, j±m / 2, k) and (i, j, k±m / 2) is used for calculation and and the z-component of the velocity at (i, j, k±m / 2) is used for calculation where m = (1, 2...M), and M is half of the length of the finite difference operator. In the embodiment of this article, M = 4 is set, achieving spatial eighth-order accuracy, and the staggered grid difference coefficients are calculated using Taylor expansion

[0105] By performing iterative stepping in the time direction using the above difference format, wave field data such as stress and velocity at each calculation point at each moment can be obtained

[0106] Step 8: Use the linear interpolation algorithm to map the simulated wave field in the computational domain back to the real physical domain to obtain the final simulated wave field and shot gather record

[0107] Applying the method provided by the present invention to the established 3D loess tableland model has obtained good numerical simulation results Figure 3 show the velocity field and Q field of the 3D loess tableland model. The model size is 475*1301*800 grids, and the grid size is 10m Figure 4 and Figure 5 are the wave field snapshots and shot gather records of acoustic wave simulation of the 3D loess tableland model respectively. The time interval of numerical simulation is 0.5ms Figure 5 and Figure 6 are the wave field snapshots and shot gather records of visco-acoustic wave simulation of the 3D loess tableland model obtained in the embodiment of the present invention respectively. The time interval of numerical simulation is 0.5ms. Comparing with the wave field snapshots and shot gather records of acoustic wave simulation of the loess tableland model, it can be seen that the wave field attenuation effect is very obvious; and the overall characteristics of the shot gather of the visco-acoustic medium are very consistent with the real data, and typical "black triangle" noise is generated below the source position. In addition, the calculation cost of this method is relatively low, and it has the potential to be applied in actual production

[0108] In this example, by adopting the above method, the problems of inaccurate simulation and high computational cost in the numerical simulation of the viscoacoustic wave equation under the condition of undulating surface are solved. The numerical simulation of the viscoacoustic wave equation can be carried out stably, accurately and efficiently, and a wave field consistent with the actual collected data can be obtained.

[0109] Those of ordinary skill in the art can understand that all or some of the steps in the methods disclosed above, and the functional modules / units in the systems and devices, can be implemented as software, firmware, hardware, and appropriate combinations thereof. In the hardware implementation, the division of the functional modules / units mentioned above does not necessarily correspond to the division of physical components; for example, one physical component may have multiple functions, or one function or step may be executed by several physical components in cooperation. Some or all of the components may be implemented as software executed by a processor, such as a digital signal processor or a microprocessor, or may be implemented as hardware, or may be implemented as an integrated circuit, such as an application-specific integrated circuit. Such software can be distributed on a computer-readable medium, which may include a computer storage medium (or non-transitory medium) and a communication medium (or transitory medium). As is well known to those of ordinary skill in the art, the term computer storage medium includes volatile and non-volatile, removable and non-removable media implemented in any method or technology for storing information, such as computer-readable instructions, data structures, program modules, or other data. Computer storage media include, but are not limited to, RAM, ROM, EEPROM, flash memory or other memory technologies, CD-ROM, digital versatile disc (DVD) or other optical disc storage, magnetic cassette, tape, magnetic disk storage or other magnetic storage devices, or any other medium that can be used to store the desired information and can be accessed by a computer. In addition, as is well known to those of ordinary skill in the art, communication media typically contain computer-readable instructions, data structures, program modules, or other data in a modulated data signal such as a carrier wave or other transmission mechanism, and may include any information delivery medium.

Claims

1. A method for numerical simulation of wave equations of viscoacoustic media on undulating surfaces, characterized in that: The method includes: Establishing a viscoacoustic medium wave equation based on the position of the seismic source, formation P-wave velocity data, formation density data, formation quality factor data, seismic source wavelet function, and surface elevation; Establishing a mapping relationship expression between physical domain data and computational domain data; Determining the viscoacoustic medium wave equation in the computational domain based on the mapping relationship expression and the viscoacoustic medium wave equation; Using the fully staggered grid high-order finite difference method to determine the wave field in the computational domain of the viscoacoustic medium wave equation in the computational domain; Mapping the wave field in the computational domain to obtain the wave field in the physical domain.

2. The method for numerical simulation of wave equations of viscoacoustic media on undulating surfaces according to claim 1, characterized in that: The mapping relationship expression between the physical domain data and the computational domain data is: In the above relationship expression, the physical domain data x = (x, y, z), the computational domain data α = (α, β, γ), Γ(x, y) is the surface elevation, and z max is the maximum depth of the computational area, and γ max = z max + Γ max , and Γ max is the maximum surface elevation.

3. The numerical simulation method for the wave equation of the visco-acoustic medium with undulating surface according to claim 1, characterized in that, The physical domain data and the computational domain data are in a one-to-one correspondence.

4. The numerical simulation method of the fluctuating surface viscoacoustic medium wave equation according to claim 1, wherein The determining the viscoacoustic medium wave equation in the computational domain based on the mapping relationship expression and the viscoacoustic medium wave equation includes: Converting the physical domain data into the computational domain data based on the mapping relationship expression and calculating the partial derivatives in the corresponding coordinate system; Obtaining the viscoacoustic medium wave equation in the computational domain based on the viscoacoustic medium wave equation and the partial derivatives.

5. The numerical simulation method of the acoustic wave equation in a viscoelastic medium with undulating surface according to claim 4, characterized in that The viscoacoustic medium wave equation is: In the above formula, ρ(x) is the formation density, v(x, t) is the particle velocity, u(x, t) is the stress, Q(x) is the formation quality factor, v u (x) is the formation P-wave velocity, s is the total weight, ∧ (l) (x, t) is the memory variable, x s is the position of the seismic source, H(t) is the step function, f(t) is the source time function; M (l) represents the weight of the l-th relaxation mechanism, N represents the number of groups of relaxation mechanisms, l represents the relaxation mechanism number, and τ (l) represents the relaxation time of the l-th relaxation mechanism, and t and t′ represent two different moments respectively.

6. The numerical simulation method for the wave equation of the viscoacoustic medium with undulating surface according to claim 4, wherein The viscoacoustic medium wave equation in the computational domain is: Among them, is the generalized partial derivative operator, and α s is the position of the seismic source in the computational domain.

7. The method for numerical simulation of wave equations of viscoacoustic media on undulating surfaces according to claim 6, characterized in that: The generalized partial derivative operator is as follows:

8. The numerical simulation method for the wave equation of the viscoacoustic medium with undulating surface according to claim 1, characterized in that The process of using the fully staggered grid high-order finite difference method to determine the wave field in the computational domain of the viscoacoustic medium wave equation in the computational domain is: Performing difference discretization on the viscoacoustic medium wave equation in the computational domain using the fully staggered grid high-order finite difference method to obtain its second-order time and high-order space difference format; Setting a time step, spatial interval, and seismic source function that meet predetermined conditions and performing iteration in the time direction to determine the simulated wave field in the computational domain.

9. The numerical simulation method of the wave equation for a visco-acoustic medium with a rough surface according to claim 8, characterized in that The fully staggered grid finite difference format of the viscoacoustic medium wave equation is: Among them, u n+1 represents the stress at the calculation point (i, j, k) at the (n + 1)th moment, u n represents the stress at the calculation point (i, j, k) at the nth moment, Δt represents the time step, M α represents the first-order derivative of the finite difference operator along the α direction, M β represents the first-order derivative of the finite difference operator along the β direction, M γ represents the first-order derivative of the finite difference operator along the γ direction; is the x-component of the particle velocity data at the (n + 1 / 2)th moment, is the y-component of the particle velocity data at the (n + 1 / 2)th moment, is the z-component of the particle velocity data at the (n + 1 / 2)th moment, is the x-component of the particle velocity data at the (n + 3 / 2)th moment, the y-component of the particle velocity data at the (n + 3 / 2)th moment, the z-component of the particle velocity data at the (n + 3 / 2)th moment, ∧ (l)，n is the memory variable at the nth moment, ∧ (l)，n+1 the memory variable at the (n + 1)th moment.

10. The method for numerical simulation of wave equations of viscoacoustic medium on undulating surface according to claim 8, characterized in that: The predetermined conditions are: where c is the sum of the difference coefficients, v max is the maximum speed, Q min is the minimum value of the quality factor, c x，max is the maximum value of the surface coordinate stretching coefficient c x is the maximum value of the surface coordinate stretching coefficient c y，max is the maximum value of the surface coordinate stretching coefficient c y is the maximum value of the surface coordinate stretching coefficient c z，max is the maximum value of the surface coordinate stretching coefficient c z is the maximum value.

11. An apparatus for numerical simulation of wave equation in a visco-acoustic medium with undulating surface, characterized in that, The device includes: a memory and a processor; the memory is used to store a program for performing numerical simulation of the viscoacoustic medium wave equation on a rugged surface, and the processor is used to read and execute the program for performing numerical simulation of the viscoacoustic medium wave equation on a rugged surface and execute the method according to any one of claims 1-10.

12. A computer-readable storage medium, on which a data processing program is stored, and the data processing program is executed by a processor to perform the method for numerical simulation of the viscoacoustic medium wave equation on a rugged surface according to any one of claims 1-10.

Images