Multi-domain finite difference-pseudo spectrum hybrid calculation method for elastic wave simulation

Abstract

The invention discloses a multi-domain finite difference-pseudo spectrum hybrid calculation method for elastic wave simulation. For each sub-region obtained by dividing the target computational domain, dispersing the sub-region in the horizontal direction through equidistant grid points, and dispersing the sub-region in the vertical direction through Chebyshev points; for each sub-discrete region, calculating a spatial derivative in the elastic wave control equation in the horizontal direction by using a high-order finite difference method, and calculating a spatial derivative in the elastic wave control equation in the vertical direction by using a Chebyshev point pseudo-spectral method; performing time integration according to the two spatial derivatives to obtain wave field data at the next moment; and applying a characteristic value boundary condition under a local coordinate system at the boundary of each discrete sub-region to carry out wave field correction, thereby realizing transmission of wave field information between the sub-regions and high-precision free surface boundary treatment. The invention provides a novel mixing method for simulating an elastic wave field, and solves the problem of high-precision and high-efficiency simulation of elastic waves under the condition of rugged topography.

Description

Technical Field

[0001] The present invention relates to the technical field of elastic wave simulation, and particularly to a multi-domain finite difference-pseudo-spectral hybrid calculation method for elastic wave simulation. Background Art

[0002] Numerical simulation of elastic waves is extremely important in seismology and exploration geophysics. Accurately simulating the wave field characteristics and propagation process of the real earth under the condition of undulating terrain is crucial for studying the propagation law of seismic waves, predicting strong ground motion, imaging the internal structure of the earth, and underground resource exploration. With the rapid progress of modern seismology and exploration geophysics, the demand for computing resources in elastic wave simulation has increased significantly. Therefore, efficient elastic wave simulation has always been the core issue of research. In addition to optimizing and accelerating programs and promoting the progress of computing hardware, developing efficient elastic wave simulation algorithms is also indispensable.

[0003] Currently used elastic wave simulation algorithms mainly include: finite difference method, high-order finite difference method, Chebyshev pseudo-spectral method, multi-domain Chebyshev pseudo-spectral method, etc. However, these algorithms all have certain defects. For example, the finite difference method has low computational accuracy and requires dense grid division to improve computational accuracy; the high-order finite difference method has high computational accuracy, but lacks a matching boundary condition processing method; the Chebyshev pseudo-spectral method is difficult to parallelize and has a high computational cost; although the multi-domain Chebyshev pseudo-spectral method solves the problem of difficult parallelization, it lacks a high-precision and stable method to process the corner points between sub-regions. In short, the existing elastic wave simulation algorithms are difficult to solve the problem of high-precision and efficient simulation of elastic waves in the case of undulating terrain.

[0004] Therefore, the prior art still needs to be improved and developed. Summary of the Invention

[0005] The technical problem to be solved by the present invention is to provide a multi-domain finite difference-pseudo-spectral hybrid calculation method for elastic wave simulation in view of the above-mentioned defects of the prior art, aiming to solve the problem that the existing elastic wave simulation algorithms are difficult to achieve high-precision and efficient simulation of elastic waves in the case of undulating terrain.

[0006] The technical solution adopted by the present invention to solve the problem is as follows:

[0007] In a first aspect, an embodiment of the present invention provides a multi-domain finite difference-pseudo-spectral hybrid calculation method for elastic wave simulation, the method comprising:

[0008] Dividing a target calculation domain into a plurality of independent sub-regions;

[0009] For each of the said sub-regions, grid discretization is performed by equally spaced grid points in the horizontal direction and by Chebyshev points in the vertical direction to obtain a number of sub-discrete regions;

[0010] According to the elastic wave control equation constructed in the curvilinear coordinate system, for each of the said sub-discrete regions, the high-order finite difference method is used to calculate the spatial derivative in the horizontal direction in the elastic wave control equation, and the pseudo-spectral method based on Chebyshev points is used to calculate the spatial derivative in the vertical direction in the elastic wave control equation; time integration is performed based on the spatial derivative in the horizontal direction and the spatial derivative in the vertical direction to obtain the wave field data at the next moment;

[0011] At the local coordinate system at the boundary of each of the said sub-discrete regions, the eigenvalue boundary conditions of the corresponding boundary type are applied to obtain the corrected wave field data at the next moment; wherein, the boundary types include interface boundary, free surface boundary and absorbing boundary.

[0012] In one implementation, the target calculation domain includes undulating terrain.

[0013] In one implementation, the target calculation domain is divided into a number of independent sub-regions, including:

[0014] Along the direction perpendicular to the free surface, the target calculation domain is divided into a number of sub-regions with arbitrary heights and horizontal layers; wherein, the edges of the said sub-regions only overlap.

[0015] In one implementation, the number of Chebyshev points is less than or equal to 15.

[0016] In one implementation, the elastic wave control equation is constructed using the first-order velocity-stress equation.

[0017] In one implementation, the elastic wave control equation constructed in the curvilinear coordinate system is:

[0018]

[0019] Wherein, U = (v x v z τ xx τ zz τ xz ) T represents the velocity and stress vectors of the wave field; v x and v z represent the components of the velocity in the x direction and the z direction; τ xx and τ zz represent the normal stress components; τ xzdenote the shear stress components; A and B denote the elastic coefficient matrices in the curvilinear coordinate system; λ and μ denote the Lame constants; ρ denotes the density; (x, z) denote the coordinates in the physical plane, and (ξ, ζ) denote the coordinates in the computational plane; denote partial derivatives.

[0020] In one implementation, at the local coordinate system at the boundary of each of the sub-discrete regions, eigenvalue boundary conditions of the corresponding boundary type are applied to obtain the wave field data at the next moment after correction, including:

[0021] For each of the sub-discrete regions, the wave field data at the next moment is transformed from the physical coordinate system to the wave field data in the local coordinate system through a coordinate transformation matrix; wherein, the local coordinate system is established according to the normal direction perpendicular to the current boundary and pointing to the outside of the current sub-discrete region, and the tangent direction;

[0022] Apply eigenvalue boundary conditions of the boundary type corresponding to the current boundary to the wave field data at the next moment in the local coordinate system to obtain the wave field data at the next moment that satisfies the eigenvalue boundary conditions;

[0023] Through the coordinate transformation matrix, transform the wave field data at the next moment that satisfies the eigenvalue boundary conditions from the local coordinate system to the physical coordinate system to obtain the wave field data at the next moment after correction.

[0024] In one implementation, the interface boundary condition is:

[0025]

[0026] where Z p = ρv p and Z s = ρv s denote the impedances of the longitudinal wave and the transverse wave respectively, and denote the velocities of the longitudinal wave and the transverse wave respectively, λ, μ denote the Lame constants, ρ denotes the density; W = (v n v t τ nn τ tt τ nt ) T denotes the wave field data in the local coordinate system (n, t), denotes the wave field data after applying the eigenvalue boundary conditions; denote the impedances of the longitudinal wave and the transverse wave in the adjacent sub-regions respectively; denote the wave field data of the adjacent sub-regions in the local coordinate system (n, t) respectively;

[0027] The free surface boundary condition is:

[0028]

[0029] The absorbing boundary conditions are as follows:

[0030]

[0031] Wherein, Z p = ρv p and Z s = ρv s represent the impedances of the longitudinal wave and the transverse wave respectively, and represent the velocities of the longitudinal wave and the transverse wave respectively, λ and μ represent Lame constants, and ρ represents density; W = (v n v t τ nn τ tt τ nt ) T represents the wave field data in the local coordinate system (n, t), represents the wave field data after applying the eigenvalue boundary conditions; represent the impedances of the longitudinal wave and the transverse wave in adjacent sub-regions respectively; represent the wave field data of adjacent sub-regions in the local coordinate system (n, t) respectively.

[0032] Advantages of the present invention: The embodiments of the present invention combine the advantages of the multi-domain method, the finite difference method, and the pseudo-spectral method to simulate the elastic wave field, and solve the problem of high-precision and efficient simulation of elastic waves in the case of undulating terrain. BRIEF DESCRIPTION OF THE DRAWINGS

[0033] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the drawings in the following description are only some embodiments recorded in the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0034] Figure 1 is a schematic flowchart of a multi-domain finite difference-pseudo-spectral hybrid calculation method for elastic wave simulation provided by an embodiment of the present invention.

[0035] Figure 2 is a schematic diagram of the mapping relationship between an arbitrary curve grid and a Cartesian grid provided by an embodiment of the present invention.

[0036] Figure 3 is a schematic diagram of model division provided by an embodiment of the present invention.

[0037] Figure 4 It is a schematic diagram of the global model provided by an embodiment of the present invention.

[0038] Figure 5 It is a schematic diagram of the foothill model provided by an embodiment of the present invention.

[0039] Figure 6 It is a waveform comparison diagram of the horizontal velocity components at all geophones in the foothill model provided by an embodiment of the present invention.

[0040] Figure 7 It is a waveform comparison diagram of the vertical velocity components at all geophones in the foothill model provided by an embodiment of the present invention. Detailed implementation manners

[0041] The present invention discloses a multi-domain finite difference-pseudo spectral hybrid calculation method for elastic wave simulation. To make the objectives, technical solutions and effects of the present invention clearer and more definite, the following further describes the present invention in detail with reference to the accompanying drawings and by way of examples. It should be understood that the specific examples described herein are only used to explain the present invention and are not used to limit the present invention.

[0042] Those skilled in the art of the present technology can understand that unless specifically stated otherwise, the singular forms "a", "an", "the" and "said" used herein may also include the plural forms. It should be further understood that the term "comprising" used in the specification of the present invention means the presence of the described features, integers, steps, operations, elements and / or components, but does not exclude the presence or addition of one or more other features, integers, steps, operations, elements, components and / or their groups. It should be understood that when we say that an element is "connected" or "coupled" to another element, it can be directly connected or coupled to other elements, or there may also be intermediate elements. In addition, the "connection" or "coupling" used herein may include wireless connection or wireless coupling. The phrase "and / or" used herein includes all or any unit and all combinations of one or more related listed items.

[0043] Those skilled in the art of the present technology can understand that unless otherwise defined, all terms (including technical terms and scientific terms) used herein have the same meaning as the general understanding of those of ordinary skill in the field to which the present invention belongs. It should also be understood that terms such as those defined in a general dictionary should be understood to have a meaning consistent with the meaning in the context of the prior art, and will not be interpreted with an idealized or overly formal meaning unless specifically defined as here.

[0044] In view of the above defects of the prior art, the present invention provides a multi-domain finite difference-pseudo spectral hybrid calculation method for elastic wave simulation. The method includes: dividing a target calculation domain into a number of independent sub-regions; for each of the sub-regions, performing grid discretization through equally spaced grid points in the horizontal direction and through Chebyshev points in the vertical direction to obtain a number of sub-discretized regions; according to the elastic wave control equation constructed in the curvilinear coordinate system, for each of the sub-discretized regions, using the high-order finite difference method to calculate the spatial derivative in the horizontal direction of the elastic wave control equation, and using the pseudo spectral method based on Chebyshev points to calculate the spatial derivative in the vertical direction of the elastic wave control equation; performing time integration based on the spatial derivative in the horizontal direction and the spatial derivative in the vertical direction to obtain the wave field data at the next moment; applying the eigenvalue boundary conditions of the corresponding boundary type in the local coordinate system at the boundary of each of the sub-discretized regions to obtain the corrected wave field data at the next moment; where the boundary types include interface boundary, free surface boundary and absorbing boundary. The present invention combines the advantages of multi-domain method, finite difference method and pseudo spectral method to simulate the elastic wave field, and solves the problem of high-precision and efficient simulation of elastic waves in the case of undulating terrain.

[0045] As Figure 1 shown, the method specifically includes:

[0046] Step S100: Divide the target calculation domain into a number of independent sub-regions.

[0047] Specifically, the target calculation domain in the present application can represent the calculation area where elastic wave simulation is required. The target calculation domain can include undulating terrain, so as to simulate various propagation phenomena of elastic waves in the case of undulating terrain. Due to the existence of surface undulations, it is necessary to consider the irregularity of the calculation domain during the numerical simulation of elastic wave propagation. Therefore, in this embodiment, the complex target calculation domain is decomposed into smaller and more easily simulated calculation parts, so as to obtain a number of independent sub-regions.

[0048] In one implementation, the dividing the target calculation domain into a number of independent sub-regions includes:

[0049] Along the direction perpendicular to the free surface, divide the target calculation domain into a number of sub-regions with arbitrary heights and horizontal layers; where each of the sub-regions only overlaps at the edges.

[0050] Specifically, for the computational domain to be divided that contains undulating terrain, along the vertical direction downward along the ground surface (i.e., the direction perpendicular to the free surface), it is divided into a series of independent sub-regions, and only the edges of the sub-regions are allowed to overlap. This series of sub-regions are horizontal layers of arbitrary height, that is, the heights of adjacent sub-regions can be different. This flexible height division method can significantly enhance the flexibility of modeling. In practical applications, the vertical height within the sub-region and the number of Chebyshev grid points in the vertical direction of the sub-region can be appropriately adjusted according to the medium parameters to avoid spatial oversampling, thereby improving the calculation efficiency.

[0051] Step S200: For each of the sub-regions, the grid is discretized by equally spaced grid points in the horizontal direction and by Chebyshev points in the vertical direction to obtain a number of sub-discretized regions.

[0052] Specifically, this embodiment adopts different discretization strategies in the horizontal and vertical directions. Taking a sub-region as an example, in the horizontal direction, the model is discretized by equally spaced grid points used in the finite difference method, that is, the grid points are equally spaced, and the horizontal distance between each grid point is the same. In the vertical direction, the grid is discretized by Chebyshev points (Chebyshev-Gauss-Lobatto, CGL) used in the Chebyshev pseudo-spectral method, that is, the grid points are distributed according to Chebyshev points.

[0053] In one implementation, the number of Chebyshev points is less than or equal to 15.

[0054] Specifically, this embodiment adopts a Chebyshev point configuration scheme of 15 points or less, which can effectively control the height of the sub-region and enhance the flexibility of modeling. And it can also avoid the over-dense phenomenon of Chebyshev points at both ends of the region, reducing the problem of low calculation efficiency caused by too small time step in numerical simulation, thereby significantly improving the calculation efficiency.

[0055] Step S300: According to the elastic wave control equation constructed in the curvilinear coordinate system, for each of the sub-discretized regions, the spatial derivative in the horizontal direction of the elastic wave control equation is calculated using the high-order finite difference method, and the spatial derivative in the vertical direction of the elastic wave control equation is calculated using the pseudo-spectral method based on Chebyshev points; time integration is performed based on the spatial derivative in the horizontal direction and the spatial derivative in the vertical direction to obtain the wave field data at the next moment.

[0056] Specifically, to achieve elastic wave simulation in the case of undulating terrain, this embodiment needs to construct an elastic wave control equation in a curvilinear coordinate system to describe the propagation process of elastic waves in the medium. In practical applications, the elastic wave control equation is solved in each discrete sub-region. Among them, the high-order finite difference method is used to calculate the spatial derivative in the horizontal direction, which has high accuracy and good computational efficiency; the pseudo-spectral method based on Chebyshev points is used to calculate the spatial derivative in the vertical direction, and the pseudo-spectral method has exponential convergence at Chebyshev points and is suitable for high-precision calculations. Combining the spatial derivatives in the horizontal and vertical directions, the wave field data at the next moment is calculated through a time integration method to achieve elastic wave simulation calculations, and finally the wave field data describing the state of elastic waves at various time and space positions is obtained. Further, a visualization tool can be used to present the wave field data as a graph or animation to more intuitively understand the propagation behavior of elastic waves.

[0057] In one implementation, the elastic wave control equation is constructed using the first-order velocity-stress equation.

[0058] Specifically, there are various types of elastic wave control equations. This embodiment comprehensively considers the wave field calculation methods of subsequent boundaries and finally selects the first-order velocity-stress equation to establish the elastic wave control equation.

[0059] In one implementation, the elastic wave control equation constructed in the curvilinear coordinate system is:

[0060]

[0061]

[0062] where U = (v x v z τ xx τ zz τ xz ) T represents the velocity and stress vectors of the wave field; v x and v z represent the components of the velocity in the x and z directions; τ xx and τ zz represent the normal stress components; τ xz represents the shear stress component; A and B represent the elastic coefficient matrices in the curvilinear coordinate system; λ and μ represent the Lamé constants; ρ represents the density; (x, z) represents the coordinates of the physical plane, and (ξ, ζ) represents the coordinates of the computational plane; represents the partial derivative.

[0063] Figure 2Shows the mapping relationship between an arbitrary curvilinear grid and a Cartesian grid. On the left, the physical space (x, z) is shown, and on the right, the computational space (ξ, ζ) is shown. For the elastic wave simulation in the case of undulating terrain, this embodiment will map the physical plane of the surface to another computational domain through a coordinate transformation, so that the grid lines in the computational domain after the coordinate transformation are all straight lines. In other words, this embodiment will map the irregular grid in the physical space to the regular rectangular grid in the computational space. In addition, Figure 2 It can also verify the aforementioned grid discretization method, that is, the model is discretized with equally spaced grid points in the horizontal direction, and the grid is discretized according to Chebyshev points in the vertical direction. Specifically, the selected numerical scheme is applied in each sub-discretization region to solve the elastic wave control equation. Among them, the application process of the numerical scheme in each sub-discretization region is as follows: the high-order finite difference method is used to calculate the spatial derivative in the ξ direction (i.e., the horizontal direction), and the pseudo-spectral method based on Chebyshev points is used to calculate the spatial derivative in the ζ direction (i.e., the vertical direction). This method is convenient for subsequent application of eigenvalue boundary conditions, thereby significantly improving the calculation accuracy when dealing with strong velocity interfaces and free surface boundaries.

[0064] For example, as Figure 2 shown, in the right figure of the computational space, the template for calculating the spatial derivative using the high-order finite difference method in the ξ direction and the template for calculating the spatial derivative using the Chebyshev pseudo-spectral method in the ζ direction are shown.

[0065] In one implementation, the finite difference format corresponding to the high-order finite difference method uses a 12th-order accurate central co-located grid, and this finite difference format can remove non-physical high-frequency numerical oscillations through filtering.

[0066] In another implementation, the finite difference format can also use other co-located grid finite difference formats other than the 12th-order accurate central co-located grid.

[0067] Step S400: Apply the eigenvalue boundary conditions of the corresponding boundary type in the local coordinate system at the boundaries of each of the sub-discretization regions to obtain the corrected wave field data at the next moment; where the boundary types include interface boundaries, free surface boundaries, and absorption boundaries.

[0068] Specifically, in the numerical simulation of elastic wave propagation, due to the limited computational domain, the existence of boundaries has a great influence on the propagation of elastic waves. Therefore, the key issue in elastic wave simulation under undulating terrain lies in the implementation of boundary conditions. To obtain accurate simulation results, this embodiment pre-sets eigenvalue boundary conditions of different boundary types to correctly handle various types of boundary situations. In practical applications, eigenvalue boundary conditions of the corresponding boundary types, such as interface boundary conditions, free surface boundary conditions, and absorbing boundary conditions, are applied in the local coordinate system of the boundary of each sub-discrete region to achieve wave field correction, thereby obtaining more realistic and accurate wave field data.

[0069] For example, Figure 3 shows a schematic diagram of model division based on the multi-domain finite difference-pseudo-spectral hybrid method, which shows three sub-regions. The topmost boundary is the free surface boundary, the boundary between sub-regions is the interface boundary, and the bottom boundary is the absorbing boundary. Eigenvalue boundary conditions are applied in the local coordinate system of the sub-region boundary. Here, n represents the normal direction perpendicular to the boundary and pointing outside the current sub-domain, and t represents the tangent direction. Figure 4 shows a schematic diagram of the global model. The gray shaded areas on the left, right, and bottom of the model represent the waveform attenuation areas.

[0070] At the boundary of each of the sub-discrete regions, in the local coordinate system, eigenvalue boundary conditions of the corresponding boundary types are applied to obtain the corrected updated wave field data, including:

[0071] For each of the sub-discrete regions, the wave field data at the next moment is transformed from the physical coordinate system to the local coordinate system through a coordinate transformation matrix; where the local coordinate system is established according to the normal direction perpendicular to the current boundary and pointing outside the current sub-discrete region, and the tangent direction.

[0072] Apply the eigenvalue boundary conditions of the boundary type corresponding to the current boundary to the wave field data at the next moment in the local coordinate system to obtain the wave field data at the next moment that satisfies the eigenvalue boundary conditions.

[0073] Through the coordinate transformation matrix, the wave field data at the next moment that satisfies the eigenvalue boundary conditions is transformed from the local coordinate system to the physical coordinate system to obtain the corrected wave field data at the next moment.

[0074] Specifically, three different eigenvalue boundary conditions are preset in this embodiment. For different boundaries, the corresponding eigenvalue boundary conditions need to be satisfied. Information transfer is achieved between different subdomains through the eigenvalue boundary conditions, so as to ensure the high-precision processing of strong velocity interfaces and free surface boundary conditions. Taking a sub-region as an example, the wave field in the local coordinate system is obtained through coordinate transformation, and the eigenvalue boundary condition is applied in the local coordinate system of the sub-region boundary to obtain a new wave field. Among them, the local coordinate system is established based on n and t, where n represents the normal direction perpendicular to the boundary and pointing to the outside of the current sub-region, and t represents the tangent direction. The new wave field is equivalent to the wave field that corrects the boundary, and further optimizes the elastic wave simulation scheme based on partitioning.

[0075] For example, the wave field W = (v n v t τ nn τ tt τ nt ) T in the local coordinate system (n, t) is obtained by applying a coordinate transformation to the wave field U in the physical coordinate system (x, z).

[0076] The formula is: W = TU,

[0077] where T is the coordinate transformation matrix, and the T matrix is:

[0078]

[0079] where, represents the unit normal vector of the interface pointing to the outside of the current subdomain.

[0080] Applying the characteristic boundary condition in the local coordinate system, a numerical value of the wave field W( new ) that satisfies the boundary condition is obtained.

[0081] Using the coordinate transformation equation, the new wave field W( new ) is rotated back to the original coordinate system (i.e., the physical coordinate system) to obtain U (new) .

[0082] The formula is: U (new) = T -1 W( new ).

[0083] Furthermore, in one implementation, the interface boundary condition is:

[0084]

[0085] where, Z p = ρv p and Z s = ρv srespectively represent the impedances of the longitudinal wave and the transverse wave, and respectively represent the velocities of the longitudinal wave and the transverse wave, λ, μ represent the Lame constants, and ρ represents the density; W = (v n v t τ nn τ tt τ nt ) T represents the wave field data in the local coordinate system (n, t), represents the wave field data after applying the eigenvalue boundary condition; respectively represent the impedances of the longitudinal wave and the transverse wave in adjacent sub-regions; respectively represent the wave field data of adjacent sub-regions in the local coordinate system (n, t).

[0086] The free surface boundary condition is:

[0087]

[0088] The absorbing boundary condition is:

[0089]

[0090] where, Z p = ρv p and Z s = ρv s respectively represent the impedances of the longitudinal wave and the transverse wave, and respectively represent the velocities of the longitudinal wave and the transverse wave, λ, μ represent the Lame constants, and ρ represents the density; W = (v n v t τ nn τ tt τ nt ) T represents the wave field data in the local coordinate system (n, t), represents the wave field data after applying the eigenvalue boundary condition; respectively represent the impedances of the longitudinal wave and the transverse wave in adjacent sub-regions; respectively represent the wave field data of adjacent sub-regions in the local coordinate system (n, t).

[0091] To prove the technical effect of the present invention, the simulation effects of the multi-domain finite difference-pseudospectral hybrid method (multidomain FDM / PSM) of the present invention and the traditional finite difference method (FDM) are compared through experiments. Figures 5 - 7 The comparison results of the two methods are shown. Among them, Figure 5For the foothill model: The multi-domain finite-difference - pseudo-spectral hybrid method uses 80 sub-regions, with a total of 80×9×1201 grid points; while the traditional finite-difference method uses 1281×2401 grid points. Figure 5 The black inverted triangles in Figure 5 represent the receivers located at the free surface, labeled R1 to R19, and the positions of the explosion sources are marked by the red stars. Figure 6 For the waveform comparison of the horizontal velocity component (Vx) at all receivers in the foothill model: The solid line represents the waveform obtained by the finite-difference method, and the dashed line represents the waveform obtained by the multi-domain finite-difference - pseudo-spectral hybrid method. Figure 7 For the waveform comparison of the vertical velocity component (Vz) at all receivers in the foothill model: The solid line represents the waveform obtained by the finite-difference method, and the dashed line represents the waveform obtained by the multi-domain finite-difference - pseudo-spectral hybrid method.

[0092] In summary, the present invention discloses a multi-domain finite-difference - pseudo-spectral hybrid calculation method for elastic wave simulation. The method includes: dividing the target calculation domain into a number of independent sub-regions; for each of the sub-regions, discretizing the grid in the horizontal direction by equally spaced grid points and in the vertical direction by Chebyshev points to obtain a number of sub-discretized regions; according to the elastic wave control equation constructed in the curvilinear coordinate system, for each of the sub-discretized regions, using the high-order finite-difference method to calculate the spatial derivative in the horizontal direction of the elastic wave control equation and using the pseudo-spectral method based on Chebyshev points to calculate the spatial derivative in the vertical direction of the elastic wave control equation; performing time integration based on the spatial derivative in the horizontal direction and the spatial derivative in the vertical direction to obtain the wave field data at the next moment; applying the eigenvalue boundary conditions of the corresponding boundary type in the local coordinate system at the boundaries of the sub-discretized regions to obtain the corrected wave field data at the next moment; where the boundary types include interface boundaries, free surface boundaries, and absorption boundaries. The present invention combines the advantages of the multi-domain method, the finite-difference method, and the pseudo-spectral method to simulate the elastic wave field, can better handle the calculation of undulating terrain, and thus effectively solves the problem of high-precision and efficient simulation of elastic waves in the case of undulating terrain.

[0093] It should be understood that the application of the present invention is not limited to the above examples. For those of ordinary skill in the art, improvements or transformations can be made according to the above description, and all such improvements and transformations should fall within the protection scope of the appended claims of the present invention.

Claims

1. A multi-domain finite difference-pseudo spectral hybrid calculation method for elastic wave simulation, characterized in that The method includes: Dividing a target computational domain into a number of independent sub-regions; For each of the sub-regions, performing grid discretization by equally spaced grid points in the horizontal direction and by Chebyshev points in the vertical direction to obtain a number of sub-discretized regions; According to the elastic wave control equation constructed in a curvilinear coordinate system, for each of the sub-discretized regions, using the high-order finite difference method to calculate the spatial derivative in the horizontal direction in the elastic wave control equation, and using the pseudo-spectral method based on Chebyshev points to calculate the spatial derivative in the vertical direction in the elastic wave control equation; performing time integration based on the spatial derivative in the horizontal direction and the spatial derivative in the vertical direction to obtain the wave field data at the next moment; Applying the eigenvalue boundary conditions of the corresponding boundary types in the local coordinate system at the boundaries of the sub-discretized regions to obtain the corrected wave field data at the next moment; wherein the boundary types include interface boundaries, free surface boundaries, and absorbing boundaries.

2. The multi-domain finite difference-pseudo-spectral hybrid calculation method for elastic wave simulation according to claim 1, characterized in that The target computational domain includes undulating terrain.

3. The multi-domain finite difference-pseudo-spectral hybrid calculation method for elastic wave simulation according to claim 2, characterized in that Dividing the target computational domain into a number of independent sub-regions includes: Along the direction perpendicular to the free surface, dividing the target computational domain into a number of sub-regions with arbitrary heights and horizontal layers; wherein, only the edges of the sub-regions overlap.

4. The multi-domain finite difference-pseudo spectral hybrid calculation method for elastic wave simulation according to claim 1, wherein The number of Chebyshev points is less than or equal to 15.

5. The multi-domain finite difference-pseudo-spectral hybrid calculation method for elastic wave simulation according to claim 1, characterized in that The elastic wave control equation is constructed using the first-order velocity-stress equation.

6. The multi-domain finite difference-pseudo spectral hybrid calculation method for elastic wave simulation according to claim 5, characterized in that The elastic wave control equation constructed in the curvilinear coordinate system is: where U = (v x v z τ xx τ zz τ xz ) T represents the velocity and stress vectors of the wave field; v x and v z represent the components of the velocity in the x - direction and z - direction; τ xx and τ zz represent the normal stress components; τ xz represents the shear stress component; A and B represent the elastic coefficient matrices in the curvilinear coordinate system; λ and μ represent the Lame constants; ρ represents the density; (x, z) represents the coordinates of the physical plane, and (ξ, ζ) represents the coordinates of the computational plane; represents the partial derivative.

7. The multi-domain finite difference-pseudo spectral hybrid calculation method for elastic wave simulation according to claim 1, characterized in that Applying the eigenvalue boundary conditions of the corresponding boundary types in the local coordinate system at the boundaries of the sub-discretized regions to obtain the corrected wave field data at the next moment includes: For each of the sub-discretized regions, transforming the wave field data at the next moment from the physical coordinate system to the local coordinate system through a coordinate transformation matrix; wherein, the local coordinate system is established according to the normal direction perpendicular to the current boundary and pointing to the outside of the current sub-discretized region, and the tangent direction; Applying the eigenvalue boundary conditions of the boundary type corresponding to the current boundary to the wave field data at the next moment in the local coordinate system to obtain the wave field data at the next moment that satisfies the eigenvalue boundary conditions; Through the coordinate transformation matrix, transforming the wave field data at the next moment that satisfies the eigenvalue boundary conditions from the local coordinate system to the physical coordinate system to obtain the corrected wave field data at the next moment.

8. The multi-domain finite difference-pseudo spectral hybrid calculation method for elastic wave simulation according to claim 1, characterized in that The interface boundary condition is: Among them, Z p = ρv p and Z s = ρv s represent the impedances of the longitudinal wave and the transverse wave respectively, and represent the velocities of the longitudinal wave and the transverse wave respectively, λ, μ represent the Lame constants, and ρ represents the density; W = (v n v t τ nn τ tt τ nt ) T represents the wave field data in the local coordinate system (n, t), represents the wave field data after applying the eigenvalue boundary condition; represent the impedances of the longitudinal wave and the transverse wave in adjacent sub-regions respectively; represent the wave field data of adjacent sub-regions in the local coordinate system (n, t) respectively; The free surface boundary condition is: Among them, Z p = ρv p and Z s = ρv s represent the impedances of the longitudinal wave and the transverse wave respectively, and represent the velocities of the longitudinal wave and the transverse wave respectively, λ, μ represent Lame constants, ρ represents density; W = (v n v t τ nn τ tt τ nt ) T represents the wave field data in the local coordinate system (n, t), represents the wave field data after applying the eigenvalue boundary condition; The absorbing boundary condition is: where Z p = ρv p and Z s = ρv s represent the impedances of the longitudinal wave and the shear wave respectively, and represent the velocities of the longitudinal wave and the shear wave respectively, λ, μ represent the Lame constants, and ρ represents the density; W = (v n v t τ nn τ tt τ nt ) T represents the wave field data in the local coordinate system (n, t), represents the wave field data after applying the eigenvalue boundary condition; represent the impedances of the longitudinal wave and the shear wave in adjacent sub-regions respectively; represent the wave field data of adjacent sub-regions in the local coordinate system (n, t) respectively.

Images