A wave field forward modeling simulation method and system based on nodal discontinuity finite element method

Abstract

This invention discloses a wavefield forward modeling method and system based on the nodal discontinuity finite element method. The method involves generating an unstructured mesh for the target continuously changing medium model, deriving relevant mesh parameters, and associating these parameters. Local solutions within the elements are approximated using Lagrange interpolation polynomials based on anchor points within the elements. The wave equation is multiplied by the test function, and the result is integrated within the elements using integration by parts to obtain the weak form of the equation system. A coordinate transformation is performed on the weak form, and the local solution approximation is applied. An interpolation polynomial approximation is then applied to the parameter matrix to obtain the weak form of the equation in the reference coordinate system. Numerical flux is calculated and substituted into the weak form of the equation to obtain the final spatial discretization scheme. This spatial discretization scheme is then incorporated into the time discretization method to obtain the final seismic wave forward modeling results. The forward modeling results are compared with those obtained using the widely accepted SEM method, and numerical errors are calculated. This invention offers high computational efficiency, high accuracy, and high stability.

Description

Technical Field

[0001] This invention belongs to the field of exploration geophysics research technology, specifically relating to a wavefield forward modeling method and system based on the nodal discontinuity finite element method. Background Technology

[0002] Numerical simulation of seismic wave propagation in complex subsurface models is a crucial foundation for full-wavefield inversion and is playing an increasingly important role in geophysics. Meanwhile, smoothly varying materials are very common in geological exploration, such as sedimentary environments. However, due to the variability of physical parameters, obtaining high-precision numerical solutions in such media is difficult.

[0003] Over the past few decades, various numerical methods have been employed in the computation of synthetic seismic records, such as finite difference methods (FDMs), finite element methods (FEMs), and spectral element methods (SEMs). Finite difference methods are the most widely used numerical methods in geophysics due to their simple numerical scheme and high computational efficiency. However, FDMs suffer from poor geometric flexibility, and the accuracy of solutions obtained in media containing specific interfaces, including acoustic-elastic coupling interfaces, bending interfaces, and free surfaces, is limited. The finite element method uses unstructured meshes to partition the computational domain. These meshes can clearly and accurately represent various interfaces present in the medium, offering high geometric flexibility, and approximate the solutions within the elements using low-order polynomials. A major drawback of FEMs is the need to invert a large global mass matrix at each time step, which significantly reduces the computational efficiency. The spectral element method is an extension of the finite element method, approximating the solutions within the elements using high-order polynomials based on Gauss-Lobatto-Legendre (GLL) interpolation points. The use of GLL interpolation points transforms the global mass matrix into a strictly diagonal matrix, thereby reducing the computational difficulty of inverting the mass matrix and improving computational efficiency. However, due to the use of GLL points, this method tends to use quadrilateral or hexahedral elements for mesh generation, which presents certain difficulties in generating such meshes when simulating complex geological structures.

[0004] Discontinuous finite element (DG) methods, due to their localized solution process and the ability to allow discontinuities in solutions at element boundaries, can serve as a viable alternative to the methods described above. The DG method was initially proposed by Reed and Hill, and further mathematically developed by Hesthaven and Warburton. Similar to the finite volume method, this method introduces numerical flux to handle discontinuities at element boundaries. Kaser and Dumbser were the first to apply the DG method to seismic wavefield forward modeling, proposing an arbitrary higher-order DG method with Riemann flux (ADER-DG) based on unstructured meshes. This method has subsequently been applied to viscoelastic and anisotropic media. Zhang et al. extended the ADER-DG method to multi-physics coupled media by explicitly handling element interfaces using Godunov numerical flux. A drawback of conventional discontinuous finite element methods is that physical properties must be constant within each element. When the element size approaches the wavelength, continuously varying physical properties can lead to errors in the simulation results. The nodal discontinuous finite element method (NDG) builds upon the DG method by introducing interpolation anchor points within the mesh elements and using anchor-based Lagrange interpolation polynomials to approximate local solutions. In complex models where the internal physical properties of elements are not constant, the nodal discontinuous finite element method can yield more accurate results. Summary of the Invention

[0005] The technical problem to be solved by this invention is to provide a wavefield forward modeling method and system based on the nodal discontinuous finite element method to address the shortcomings of the prior art. The method approximates the parameter matrix based on the physical property parameters at the anchor points inside the element, realizes the calculation and study of seismic wave propagation in continuously changing media, and analyzes the accuracy and feasibility of the method. This invention is used to solve the technical problems of low accuracy and error susceptibility of the discontinuous finite element method when dealing with continuously changing media models.

[0006] The present invention adopts the following technical solution:

[0007] A wave field forward modeling method based on the nodal discontinuity finite element method includes the following steps:

[0008] S1. Perform unstructured mesh generation on the target continuously changing medium model, export the mesh information and associate it with relevant physical property parameters to obtain the target continuously changing medium model parameters;

[0009] S2. Determine the wave equation;

[0010] S3. Determine the nodal discontinuity finite element method, and determine the spatial discretization format of the wave equation obtained in step S2 based on the target continuously changing medium model parameters obtained in step S1.

[0011] S4. Determine the time discretization method and the time step calculation method;

[0012] S5. Substitute the spatial discretization scheme of the wave equation obtained in step S3 into the time discretization method obtained in step S4, and iterate according to the time step obtained in step S4 to calculate the solution at each time step within the test duration, thus forming the final wave field forward modeling result.

[0013] Preferably, in step S1, the mesh information includes the coordinates of the cell nodes and the cell positional relationships.

[0014] Preferably, in step S2, the wave equations are first-order forms of the acoustic wave equation and the elastic wave equation, as follows:

[0015]

[0016] Where u is the solution vector, g is the source term, matrices A1 and A2 are coefficient matrices, t is time, and x and y are two-dimensional spatial coordinates.

[0017] Preferably, step S3 specifically includes:

[0018] S301. After dividing the computational domain into non-overlapping triangular elements, the local solution within the element is approximated using an interpolation polynomial.

[0019] S302. Multiply the test function with the wave equation and integrate the result within the element to obtain the equation system to be solved. Obtain the weak form of the equation system by integration by parts.

[0020] S303, Define reference unit D r The parameter matrix is ​​approximated by Lagrange polynomials, and coordinate transformation is performed to obtain the weak form of the equation system in the reference coordinate system. The boundary integral is then substituted into the weak form of the equation system to obtain the final discrete scheme.

[0021] More preferably, in step S301, the local solution u(x,t) within each unit is:

[0022]

[0023] Among them, l i (x) is based on unit D k The Lagrange interpolation polynomial with internal interpolation anchors, N p The number of interpolation points.

[0024] More preferably, in step S302, the weak form of the equation system is as follows:

[0025]

[0026] Where f is the boundary The flux at point u is the solution vector, g is the source term, matrices A1 and A2 are coefficient matrices, t is time, x and y are two-dimensional spatial coordinates, and l i For the test function, D k Let be the integration region contained in the k-th triangular unit. This represents the boundary of the triangular unit.

[0027] More preferably, in step S303, the discrete format is:

[0028]

[0029] Among them, u j J is the solution vector at the anchor point within the element. k For unit D k The coordinate transformation Jacobian determinant, l i For the test function, l j Let A be the interpolation polynomial for the local solution within the element. ξ,m Let ζ = (ξ, η) be the reference coordinate system, and l m Let g be the interpolation polynomial of the coefficient matrix. j The source vector value, T is the Jacobian determinant at the cell boundary. e This is the coordinate transformation matrix. and The coefficient matrix, D is the solution vector at the element boundary anchor point. r,e This is the e-th boundary of the reference unit.

[0030] Preferably, in step S4, the time discretization method adopts the fourth-order Runge-Kutta method, and the equations y'=f(t,y) and y(t0)=y0 have the following forms:

[0031]

[0032] Among them, y n For t n The solutions to the equation at time t, k1, k2, k3, k4 are given by t. n With t n+1 The values ​​of the function f(t,y) at each node, where Δt is the time step.

[0033] Preferably, in step S5, the numerical error e is specifically:

[0034]

[0035] Where n is the total number of time steps in the calculation, s NDG The result calculated by the method of this invention is s. SEMThe results are calculated using the widely accepted SEM method.

[0036] Secondly, embodiments of the present invention provide a wave field forward modeling simulation system based on the nodal discontinuity finite element method, comprising:

[0037] The parameter module performs unstructured mesh generation on the target continuously changing medium model, exports the mesh information and associates it with relevant physical property parameters to obtain the parameters of the target continuously changing medium model.

[0038] The equation module determines the wave equation;

[0039] The discrete module determines the spatial discrete format of the wave equation obtained by the equation module based on the parameters of the target continuously changing medium model obtained by the parameter module.

[0040] The computation module determines the time discretization method and the time step calculation method;

[0041] The simulation module substitutes the spatial discretization scheme of the wave equation obtained from the discretization module into the time discretization method obtained from the calculation module, and analyzes the accuracy of the obtained numerical error method.

[0042] Compared with the prior art, the present invention has at least the following beneficial effects:

[0043] A wave field forward modeling method based on the nodal discontinuity finite element method is proposed, comprising the following steps: model establishment and mesh generation; determination of the wave equation; determination of the spatial discretization method; determination of the temporal discretization method; numerical solution and method verification. The method steps are clear and well-defined, with close relationships between each step, clear logic, and ease of implementation.

[0044] Furthermore, in step S2, the mesh information includes element node coordinates. The purpose of this is to calculate the specific coordinates of anchor points within an element using these coordinates. Based on the anchor point coordinates and the continuously changing medium model, the physical property parameters at the anchor point can be determined, providing a basis for the subsequent polynomial approximation of the coefficient matrix. The mesh information also includes element positional relationships, the purpose of which is to recover the mesh in subsequent system implementation using these relationships.

[0045] Furthermore, step S301 approximates the local solutions within the non-overlapping triangular elements using interpolation polynomials. The purpose is to derive a discrete format for the local solutions, used in the derivation of the weak form of the equations in step S302. Step S303 approximates the parameter matrix using polynomials. This aims to represent the parameter matrix in discrete form, allowing the integral sign to be removed, and further incorporating the physical properties at the anchor points into the calculation, enabling the method to handle continuously changing media. Step S303 defines a reference element and a reference coordinate system, and transforms the equation system from non-reference triangular elements to the reference element for unified calculation, facilitating parallel computation and improving computational efficiency.

[0046] Furthermore, in step S4, the fourth-order Runge-Kutta method is used as the time discretization method. Its advantage is that this method indirectly adopts the Taylor expansion method, and uses a linear combination of function values ​​at several nodes to replace the derivative of the function. This avoids the need for differentiation calculation while ensuring the high accuracy of the method. It is also simple in form and easy to implement.

[0047] It is understandable that the beneficial effects of the second aspect mentioned above can be found in the relevant descriptions in the first aspect mentioned above, and will not be repeated here.

[0048] In summary, the simulation method of this invention has high computational efficiency and the obtained forward modeling results have high accuracy.

[0049] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0050] Figure 1 The waveform diagrams for the overthrust model are shown below. (a) is a schematic diagram of the longitudinal wave velocity of the smoothed overthrust model, (b) is a schematic diagram of the transverse wave velocity, and (c) is a schematic diagram of the density.

[0051] Figure 2 The image shows a snapshot of the wave field of the overthrust model, where (a) is t=0.8s, (b) is t=1.0s, (c) is t=1.2s, (d) is t=1.4s, (e) is t=1.6s, and (f) is t=1.8s.

[0052] Figure 3 Comparison figures of the numerical solutions and SEM results of particle displacements in the x and z directions are shown, where (a) is a comparison figure of particle displacements in the x direction and (b) is a comparison figure of particle displacements in the y direction.

[0053] Figure 4The numerical solution and SEM results of the particle displacement at (11000m, -1250m) are shown in the figure. (a) is the comparison figure of particle displacement in the x direction, and (b) is the comparison figure of particle displacement in the y direction.

[0054] Figure 5 A schematic diagram of a computer device provided in an embodiment of the present invention;

[0055] Figure 6 This is a block diagram of a chip according to an embodiment of the present invention;

[0056] Figure 7 This is a schematic diagram of the process of the present invention. Detailed Implementation

[0057] 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, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0058] In the description of this invention, it should be understood that the terms "comprising" and "including" indicate the presence of the described features, integrals, steps, operations, elements and / or components, but do not exclude the presence or addition of one or more other features, integrals, steps, operations, elements, components and / or collections thereof.

[0059] It should also be understood that the terminology used in this specification is for the purpose of describing particular embodiments only and is not intended to limit the invention. As used in this specification and the appended claims, the singular forms “a,” “an,” and “the” are intended to include the plural forms unless the context clearly indicates otherwise.

[0060] It should also be further understood that the term "and / or" as used in this specification and the appended claims refers to any combination and all possible combinations of one or more of the associated listed items, and includes such combinations. For example, A and / or B can represent three cases: A alone, A and B simultaneously, and B alone. Additionally, the character " / " in this invention generally indicates that the preceding and following objects have an "or" relationship.

[0061] It should be understood that although terms such as first, second, third, etc., may be used in the embodiments of the present invention to describe the preset range, these preset ranges should not be limited to these terms. These terms are only used to distinguish the preset ranges from one another. For example, without departing from the scope of the embodiments of the present invention, the first preset range may also be referred to as the second preset range, and similarly, the second preset range may also be referred to as the first preset range.

[0062] Depending on the context, the word "if" as used here can be interpreted as "when," "when," "in response to determination," or "in response to detection." Similarly, depending on the context, the phrase "if determination" or "if detection (of the stated condition or event)" can be interpreted as "when determination," "in response to determination," "when detection (of the stated condition or event)," or "in response to detection (of the stated condition or event)."

[0063] The accompanying drawings illustrate various structural schematic diagrams according to embodiments disclosed in this invention. These drawings are not to scale, and some details have been enlarged for clarity, and some details may have been omitted. The shapes of the various regions and layers shown in the drawings, as well as their relative sizes and positional relationships, are merely exemplary and may deviate from reality due to manufacturing tolerances or technical limitations. Furthermore, those skilled in the art can design regions / layers with different shapes, sizes, and relative positions as needed.

[0064] This invention provides a wavefield forward modeling method based on the nodal discontinuity finite element method. Employing this method based on upwind schemes and numerical flux through nodal discontinuity finite element methods, it tests media with continuously varying physical properties. This offers a novel approach to studying seismic wave propagation in continuously varying media, providing seismic records (displacement, velocity, acceleration) at specific locations and a snapshot of the overall wavefield of the model. This has significant implications for practical applications.

[0065] Please see Figure 7 This invention discloses a wavefield forward modeling simulation method based on the nodal discontinuity finite element method, comprising the following steps:

[0066] S1. For the target continuously changing medium model, perform unstructured mesh generation on the target continuously changing medium model, export mesh information including element node coordinates and element position relationships, and obtain the target continuously changing medium model parameters after associating relevant physical property parameters.

[0067] The research focused on a smoothed overthrust model, which includes three parameters: P-wave velocity (P-velocity), S-wave velocity (S-velocity), and density (Density). A triangular mesh was used for mesh generation. The derived mesh information includes mesh node coordinates, inter-element positional relationships, intra-element anchor point coordinates, and whether the element is located on a free boundary or an absorbing boundary.

[0068] The mesh is a triangular mesh, and the mesh size is defined as the side length of the triangular element. The size h is calculated as follows:

[0069]

[0070] Where N is the order of the interpolation polynomial, w min This represents the minimum wavelength used in the calculation process.

[0071] S2. Determine the wave equation;

[0072] The propagation of waves in fluids can be described by the sound wave equation, while in solids it can be represented by the elastic wave equation.

[0073] Under small deformation conditions, the elastic wave equation can be written as a system of first-order equations as follows:

[0074]

[0075]

[0076]

[0077]

[0078]

[0079] Where λ and μ are Lamé constants, ρ is density, and v = (u, v) is particle velocity.

[0080] The wave equations are first-order forms of the sound wave equation and the elastic wave equation, and their expressions are shown below:

[0081]

[0082] For the elastic wave equation, the velocity-stress vector u = (σ xx ,σ yy ,σ xy ,u,v), where g is the source term.

[0083] Matrices A1 and A2 represent the physical properties of the elastic medium, defined as follows:

[0084]

[0085]

[0086] The propagation of seismic waves in fluids is described as follows:

[0087]

[0088]

[0089]

[0090] For the acoustic wave equation, u = (p, u, v), matrices A1 and A2 are defined as follows:

[0091]

[0092]

[0093] Where λ is Lamé's constant and ρ is density.

[0094] S3. Determine the nodal discontinuity finite element method, and determine the spatial discretization format of the wave equation obtained in step S2 based on the target continuously changing medium model parameters obtained in step S1.

[0095] S301. After dividing the computational domain into non-overlapping triangular elements, the local solution within the element is approximated using an interpolation polynomial.

[0096] The computational domain is divided into non-overlapping triangular units as follows:

[0097]

[0098] Among them, Ω n This is the result after discretizing the computational domain.

[0099] The local solution within each element is represented as:

[0100]

[0101] Among them, l i (x) is based on unit D k Lagrange interpolation polynomial with internal interpolation anchors and N interpolation points. p The interpolation polynomial order N satisfies

[0102] S302. Define a set of test functions whose expressions are similar to l i (x) are the same. Multiply the test function with the wave equation and integrate the result in the unit to obtain the equation system to be solved. Obtain the weak form of the equation system by integration by parts.

[0103] To obtain the values ​​of the solution vector at the nodes, this invention introduces a series of test functions, whose definitions are similar to those of the interpolation function l. i (x) is the same; by multiplying the test function and the wave equation and integrating the result in each unit, we obtain the following system of equations:

[0104]

[0105] We obtain the weak form of the equation system using integration by parts:

[0106]

[0107] Where f is the boundary The flux at the location, the flux f will be the numerical flux f of the welcoming style in subsequent calculations. * What it replaced.

[0108] S303, Define reference unit D r The parameter matrix is ​​approximated by Lagrange polynomials, and coordinate transformation is performed to obtain the weak form of the equation system in the reference coordinate system. The boundary integral is then substituted into the weak form of the equation system to obtain the final discrete scheme.

[0109] Reference Unit D r The definition is as follows:

[0110] D r ={ζ=(ξ,η)|ξ,η≥-1;ξ+η≤0}

[0111] Where ζ=(ξ,η) is the reference coordinate system.

[0112] The coordinate transformation between ordinary elements and reference elements is achieved using the Jacobian determinant, as follows:

[0113]

[0114] Since the physical properties do not remain constant within the element, the parameter matrices A1 and A2 cannot be simply extracted from the integral; then, for and By expanding upon this, we obtain:

[0115]

[0116] in, Apply Lagrange polynomials to matrix A ξ and A η Approximation:

[0117]

[0118]

[0119] By approximating the parameter matrix using Lagrange polynomials and performing coordinate transformations, we can obtain the weak form of the equation system in the reference coordinate system, as shown below:

[0120]

[0121] The following stiffness matrix is ​​extracted from it:

[0122]

[0123]

[0124]

[0125] Furthermore, numerical flux is Godunov flux, defined as follows:

[0126]

[0127] Among them, F l For, F r Let T be the coordinate rotation matrix, and u be the coordinate rotation matrix. k for, for.

[0128]

[0129]

[0130] Where A1 is, in order to solve |A1|, we first perform eigenvalue decomposition on matrix A1: A1 = RΛR -1 ;

[0131] Λ=diag(-c p ,-c s ,0,c p ,c s )

[0132]

[0133] Matrix |A1| is derived from |Λ|=diag(|-c p |,|-c s |,0,|c p |,|c s |) and R are obtained

[0134]

[0135] Similarly, for F l and F r After performing a polynomial approximation, the boundary integral is obtained. Substituting this into the weak form yields the final discrete scheme.

[0136] For matrix F l and F r Perform polynomial approximation

[0137]

[0138] The boundary integral is expressed as:

[0139]

[0140]

[0141] Among them, D r,e Let ζ be the e-th boundary of the reference element, and l be the reference coordinate system after coordinate transformation. i For the test function, l j The interpolation polynomial for local solutions. This is the pre-computed matrix at the cell boundary.

[0142] S4. Determine the time discretization method and the time step calculation method;

[0143] The time discretization method employs the fourth-order Runge-Kutta method, which indirectly utilizes the Taylor expansion method. It replaces the function's derivative with a linear combination of function values ​​at several nodes, and then determines its coefficients using the Taylor expansion. For ease of expression, the equation to be solved is replaced with y'=f(t,y), y(t0)=y0, as shown below:

[0144]

[0145] Among them, y n For t n The solution to the equation at time t is given by k1, k2, k3, and k4, where k1, k2, k3, and k4 are the function values ​​at each node, and Δt is the time step.

[0146] k1=f(t n ,y n )

[0147]

[0148]

[0149] k4=f(t n +Δt,y n +Δtk3)

[0150] The time step of the time discretization method needs to satisfy the following Courant-Friedrichs-Lewy condition, as follows:

[0151]

[0152] Among them, c max The maximum wave velocity during the calculation process is given by C, where Δx is the minimum distance between interpolation points, and C is the maximum wave velocity during the calculation process. CFL The proportionality coefficient is generally less than 0.4.

[0153] S5. Substitute the spatial discretization scheme of the wave equation obtained in step S3 into the time discretization method obtained in step S4, and iterate according to the time step calculated in step S4 to calculate the solution at each time step within the test duration, thereby forming the final wave field forward modeling result. By comparing the obtained forward modeling result with the analytical solution or the semi-analytical solution widely recognized in the industry, the numerical error of the method is calculated, and the overall accuracy and reliability of the method are verified.

[0154] The numerical errors used to ensure the accuracy and reliability of the analytical method are shown below:

[0155]

[0156] Where n is the total number of time steps in the calculation, s NDG The result calculated by the method of this invention is s. SEM The results are calculated using the widely accepted SEM method.

[0157] Numerical error measures the reliability of the calculation results of a method; the smaller the error, the higher the accuracy of the method.

[0158] In another embodiment of the present invention, a wave field forward modeling simulation system based on the nodal discontinuity finite element method is provided. This system can be used to implement the above-mentioned wave field forward modeling simulation method based on the nodal discontinuity finite element method. Specifically, the wave field forward modeling simulation system based on the nodal discontinuity finite element method includes a parameter module, an equation module, a discretization module, a calculation module, and a simulation module.

[0159] The parameter module performs unstructured mesh generation on the target continuously changing medium model, derives the mesh information and associates it with relevant physical property parameters to obtain the target continuously changing medium model parameters.

[0160] The equation module determines the wave equation;

[0161] The discrete module determines the spatial discrete format of the wave equation obtained by the equation module based on the parameters of the target continuously changing medium model obtained by the parameter module.

[0162] The computation module determines the time discretization method and the time step calculation method;

[0163] The simulation module substitutes the spatial discretization scheme of the wave equation obtained from the discretization module into the time discretization method obtained from the calculation module, and analyzes the accuracy of the obtained numerical error method.

[0164] In another embodiment of the present invention, a terminal device is provided, comprising 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 to achieve a corresponding method flow or corresponding function. The processor described in this embodiment can be used for the operation of a wave field forward modeling simulation method based on the nodal discontinuity finite element method, including:

[0165] Unstructured meshing is performed on the target continuously changing medium model, mesh information is exported and associated with relevant physical property parameters to obtain the target continuously changing medium model parameters; wave equation is determined; nodal discontinuity finite element method is determined, and spatial discretization scheme of wave equation is determined based on target continuously changing medium model parameters; time discretization method and time step calculation method are determined; the spatial discretization scheme of wave equation is substituted into time discretization method, and iterative calculation is performed according to time step to calculate the solution at each time step within the test duration, thus forming the final wave field forward modeling result.

[0166] Please see Figure 5 The terminal device is a computer device. In this embodiment, the computer device 60 includes a processor 61, a memory 62, and a computer program 63 stored in the memory 62 and executable on the processor 61. When executed by the processor 61, the computer program 63 implements the fluid composition calculation method in the reservoir stimulation wellbore of this embodiment. To avoid repetition, these details are not elaborated here. Alternatively, when executed by the processor 61, the computer program 63 implements the functions of each model / unit in the wavefield forward modeling simulation system based on the nodal discontinuity finite element method of this embodiment. To avoid repetition, these details are not elaborated here.

[0167] Computer device 60 can be a desktop computer, laptop, handheld computer, cloud server, or other computing device. Computer device 60 may include, but is not limited to, a processor 61 and a memory 62. Those skilled in the art will understand that... Figure 5 This is merely an example of computer device 60 and does not constitute a limitation on computer device 60. It may include more or fewer components than shown, or combine certain components, or different components. For example, computer device may also include input / output devices, network access devices, buses, etc.

[0168] The processor 61 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. A general-purpose processor may be a microprocessor or any conventional processor.

[0169] The memory 62 can be an internal storage unit of the computer device 60, such as a hard disk or RAM of the computer device 60. The memory 62 can also be an external storage device of the computer device 60, such as a plug-in hard disk, smart media card (SMC), secure digital (SD) card, flash card, etc. equipped on the computer device 60.

[0170] Furthermore, the memory 62 may include both internal storage units of the computer device 60 and external storage devices. The memory 62 is used to store computer programs and other programs and data required by the computer device. The memory 62 can also be used to temporarily store data that has been output or will be output.

[0171] Please see Figure 6 The terminal device is a chip. In this embodiment, the chip 600 includes a processor 622, which may be one or more, and a memory 632 for storing computer programs executable by the processor 622. The computer program stored in the memory 632 may include one or more modules, each corresponding to a set of instructions. Furthermore, the processor 622 may be configured to execute the computer program to perform the wave field forward modeling simulation method based on the nodal discontinuity finite element method described above.

[0172] Additionally, chip 600 may also include a power supply component 626 and a communication component 650. The power supply component 626 can be configured to perform power management of chip 600, and the communication component 650 can be configured to enable communication of chip 600, such as wired or wireless communication. Furthermore, chip 600 may also include an input / output (I / O) interface 658. Chip 600 can operate on an operating system stored in memory 632.

[0173] In another embodiment of the present invention, a storage medium is also provided, specifically a computer-readable storage medium (memory). This computer-readable storage medium is a memory device in a terminal 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 terminal device and extended storage media supported by the terminal 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.

[0174] One or more instructions stored in a computer-readable storage medium can be loaded and executed by a processor to implement the corresponding steps of the wave field forward modeling method based on the nodal discontinuity finite element method in the above embodiments; one or more instructions in the computer-readable storage medium are loaded and executed by the processor to perform the following steps:

[0175] Unstructured meshing is performed on the target continuously changing medium model, mesh information is exported and associated with relevant physical property parameters to obtain the target continuously changing medium model parameters; wave equation is determined; nodal discontinuity finite element method is determined, and spatial discretization scheme of wave equation is determined based on target continuously changing medium model parameters; time discretization method and time step calculation method are determined; the spatial discretization scheme of wave equation is substituted into time discretization method, and iterative calculation is performed according to time step to calculate the solution at each time step within the test duration, thus forming the final wave field forward modeling result.

[0176] 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, not all, of the embodiments of the present invention. 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. Therefore, the following detailed description of the embodiments of the present 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 present invention without inventive effort are within the scope of protection of the present invention.

[0177] Example

[0178] This example uses a smoothed overthrust model as the experimental target. The model size is 3.72km × 16km, and the model data includes P-wave velocity, S-wave velocity, and density.

[0179] Please see Figure 1 , Figure 1 (a) Figure 1 (b) Figure 1 (c) Schematic diagrams showing the transverse wave velocity, longitudinal wave velocity, and density of the smoothed overthrust model. In this example, the seismic source is located at (8000m, 0m), the source time function is a Ricker wavelet, the center frequency is 10Hz, and the amplitude is 1×10⁻⁶. 14 The time delay is 0.05s. The receivers are distributed on a straight line from (11000m, -500m) to (11000m, -1750m), with a receiver spacing of 50m.

[0180] Please see Figure 2 , Figure 2 These are snapshots of the wave field at different times for the velocities of horizontal particles, from which surface waves propagating along the free surface can be clearly observed. Due to the inhomogeneity of the model, scattering and deformation of seismic waves can be observed at t = 1.4s, 1.6s, and 1.8s.

[0181] Please see Figure 3 , Figure 3 The figure shows a comparison between the obtained numerical solution and the SEM results, which demonstrates that the results obtained by the method used are in good agreement with the SEM results.

[0182] Please see Figure 4 To more clearly illustrate the accuracy of the method, Figure 4The results of the receiver located at (11000m, -1250m) were shown separately, with RMS errors of 0.543 and 0.567 for tangential and radial displacement, respectively, demonstrating the accuracy of the method of the present invention in processing continuously changing media.

[0183] In summary, this invention presents a wave field forward modeling method and system based on the nodal discontinuous finite element method. By using anchor points within the element and performing polynomial approximation on the coefficient matrix in the equation, the physical property parameters at the anchor points are incorporated into the calculation. This solves the problem of low accuracy in the traditional discontinuous finite element method when dealing with continuously changing media. Furthermore, the use of the fourth-order Runge-Kutta method as the time discretization method offers advantages such as high computational efficiency, high computational accuracy, and ease of implementation.

[0184] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the above-described division of functional units and modules is merely an example. In practical applications, the above functions can be assigned to different functional units and modules as needed, that is, the internal structure of the device can be divided into different functional units or modules to complete all or part of the functions described above. The functional units and modules in the embodiments can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit. The integrated unit can be implemented in hardware or as a software functional unit. Furthermore, the specific names of the functional units and modules are only for easy differentiation and are not intended to limit the scope of protection of this application. The specific working process of the units and modules in the above system can be referred to the corresponding process in the foregoing method embodiments, and will not be repeated here.

[0185] In the above embodiments, the descriptions of each embodiment have different focuses. For parts that are not described in detail or recorded in a certain embodiment, please refer to the relevant descriptions of other embodiments.

[0186] Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed in this invention can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementations should not be considered beyond the scope of this invention.

[0187] In the embodiments provided by this invention, it should be understood that the disclosed devices / terminals and methods can be implemented in other ways. For example, the device / terminal embodiments described above are merely illustrative. For instance, the division of modules or units is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection shown or discussed may be through some interfaces; the indirect coupling or communication connection between devices or units may be electrical, mechanical, or other forms.

[0188] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

[0189] Furthermore, the functional units in the various embodiments of the present invention can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit. The integrated unit can be implemented in hardware or as a software functional unit.

[0190] If the integrated module / unit is implemented as a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, all or part of the processes in the methods of the above embodiments can also be implemented by a computer program instructing related hardware. The computer program can be stored in a computer-readable storage medium, and when executed by a processor, it can implement the steps of the various method embodiments described above. The computer program includes computer program code, which can be in the form of source code, object code, executable files, or certain intermediate forms. The computer-readable medium can include: any entity or device capable of carrying the computer program code, recording media, USB flash drives, portable hard drives, magnetic disks, optical disks, computer memory, read-only memory (ROM), random-access memory (RAM), electrical carrier signals, telecommunication signals, and software distribution media, etc. It should be noted that the content included in the computer-readable medium can be appropriately added or removed according to the requirements of legislation and patent practice in the jurisdiction. For example, in some jurisdictions, according to legislation and patent practice, computer-readable media do not include electrical carrier signals and telecommunication signals.

[0191] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. 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, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0192] 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, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0193] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0194] The above content is only for illustrating the technical concept of the present invention and should not be construed as limiting the scope of protection of the present invention. Any modifications made to the technical solution based on the technical concept proposed in this invention shall fall within the scope of protection of the claims of this invention.

Claims

1. A wave field forward modeling simulation method based on the nodal discontinuity finite element method, characterized in that, Includes the following steps: S1. Perform unstructured mesh generation on the target continuously changing medium model, export the mesh information and associate it with relevant physical property parameters to obtain the target continuously changing medium model parameters; S2. Determine the wave equation; S3. Determine the nodal discontinuity finite element method. Based on the target continuously changing medium model parameters obtained in step S1, determine the spatial discretization scheme of the wave equation obtained in step S2, specifically: S301. After dividing the computational domain into non-overlapping triangular elements, the local solution within the element is approximated using an interpolation polynomial. S302. Multiply the test function by the wave equation and integrate the result within the element to obtain the system of equations to be solved. Obtain the weak form of the system of equations by integration by parts. The weak form of the system of equations is as follows: in, For the boundary Circulation at the location For the solution vector, For the source term, the matrix and The coefficient matrix, For time, , Two-dimensional spatial coordinates, For the test function, Let be the integration region contained in the k-th triangular unit. The boundary of the triangular unit; S303, Define Reference Unit The parameter matrix is ​​approximated by Lagrange polynomials and coordinate transformation is performed to obtain the weak form of the equation system in the reference coordinate system. The boundary integral is then substituted into the weak form of the equation system to obtain the final discrete scheme. S4. Determine the time discretization method and the time step calculation method; S5. Substitute the spatial discretization scheme of the wave equation obtained in step S3 into the time discretization method obtained in step S4, and iterate according to the time step obtained in step S4 to calculate the solution at each time step within the test duration, thus forming the final wave field forward modeling result.

2. The wave field forward modeling simulation method based on the nodal discontinuity finite element method according to claim 1, characterized in that, In step S1, the mesh information includes the coordinates of the element nodes and the positional relationships of the elements.

3. The wave field forward modeling simulation method based on the nodal discontinuity finite element method according to claim 1, characterized in that, In step S2, the wave equations are the first-order forms of the acoustic wave equation and the elastic wave equation, as follows: in, For the solution vector, For the source term, the matrix and The coefficient matrix, For time, , These are two-dimensional spatial coordinates.

4. The wave field forward modeling simulation method based on the nodal discontinuity finite element method according to claim 1, characterized in that, In step S301, the local solution within each element for: in, Based on unit Lagrange interpolation polynomial with internal interpolation anchors. The number of interpolation points.

5. The wave field forward modeling simulation method based on the nodal discontinuity finite element method according to claim 1, characterized in that, In step S303, the discrete format is: in, The solution vector at the anchor point within the cell. For unit The coordinate transformation Jacobian determinant, For the test function, The interpolation polynomial for the local solution within the element. for, As a reference coordinate system, The interpolation polynomial of the coefficient matrix, The source vector value, Let be the Jacobian determinant at the cell boundary. This is the coordinate transformation matrix. and The coefficient matrix, This is the solution vector at the unit boundary anchor point. This is the e-th boundary of the reference unit.

6. The wave field forward modeling simulation method based on the nodal discontinuity finite element method according to claim 1, characterized in that, In step S4, the time discretization method adopts the fourth-order Runge-Kutta method, and the equations... It has the following forms: in, for The solution to the equation at time t. , , , for and Functions at each node The value, For time step.

7. The wave field forward modeling simulation method based on the nodal discontinuity finite element method according to claim 1, characterized in that, In step S5, numerical error Specifically: in, The total number of time steps in the calculation result. The calculation results are those of the method of this invention. The results are calculated using the widely accepted SEM method.

8. A wave field forward modeling simulation system based on the nodal discontinuity finite element method, characterized in that, include: The parameter module performs unstructured mesh generation on the target continuously changing medium model, exports the mesh information and associates it with relevant physical property parameters to obtain the parameters of the target continuously changing medium model. The equation module determines the wave equation; The discretization module determines the spatial discretization scheme of the wave equations obtained from the equation module based on the parameters of the target continuously changing medium model obtained from the parameter module. Specifically: After dividing the computational domain into non-overlapping triangular elements, the local solutions within the elements are approximated using interpolation polynomials. Multiply the test function by the wave equation and integrate the result within the element to obtain the system of equations to be solved. The weak form of the system of equations is obtained by integration by parts, and the weak form of the system of equations is as follows: in, For the boundary Circulation at the location For the solution vector, For the source term, the matrix and The coefficient matrix, For time, , Two-dimensional spatial coordinates, For the test function, Let be the integration region contained in the k-th triangular unit. The boundary of the triangular unit; Define reference unit The parameter matrix is ​​approximated by Lagrange polynomials and coordinate transformation is performed to obtain the weak form of the equation system in the reference coordinate system. The boundary integral is then substituted into the weak form of the equation system to obtain the final discrete scheme. The computation module determines the time discretization method and the time step calculation method; The simulation module substitutes the spatial discretization scheme of the wave equation obtained from the discretization module into the time discretization method obtained from the calculation module, and analyzes the accuracy of the obtained numerical error method.

Images