Earth medium sound wave scattering wave field holographic numerical simulation method and system and storage medium

Abstract

The invention relates to the field of seismic exploration, and discloses an earth medium sound wave scattering wave field holographic numerical simulation method and system and a storage medium, and the method comprises the steps: setting a time domain sound wave equation aiming at the characteristics of an earth medium; dividing the time domain acoustic wave equation into a background wave field equation and a scattered wave field equation based on velocity field separation and background field and scattered field separation; obtaining a time domain scattering wave field wave equation based on the scattering wave field equation; determining a spatial wave number mixed domain scattered wave equation set in the frequency domain based on the time domain scattered wave field wave equation; and calculating spatial domain scattering wave field distribution based on the spatial wave number mixed domain scattering wave equation set in the frequency domain as a holographic numerical simulation result of the earth medium sound wave scattering wave field. Therefore, the sound wave numerical simulation model is established according to the characteristics of the earth medium, and the consistency of the external boundary of the model and the real boundary condition is ensured, so that the accurate holographic numerical simulation of the earth medium sound wave scattering wave field is realized.

Description

Technical Field

[0001] The present invention belongs to the field of seismic exploration, and in particular relates to a holographic numerical simulation method, system and storage medium for earth medium acoustic wave scattering wave field. Background Art

[0002] Seismic exploration based on reflection waves is a widely used method in geological exploration. It can achieve good results when the geological structure is relatively uniform. However, with the continuous expansion of the application field of seismic exploration, the geological conditions of oil and minerals are becoming more and more complex, and the heterogeneity of underground media is increasing. Only scattered wave seismic exploration technology for heterogeneous and complex media can be applied to more complex geological conditions.

[0003] Scattered wave seismic exploration is based on scattering theory. In an inhomogeneous medium, seismic waves will scatter when encountering such an inhomogeneous geological body. The reflected wave of the layered medium is the superposition of scattered waves generated by each scattering point on the layered medium. The scales are different, the inhomogeneities are different, and the forms of scattered waves are also different. Therefore, the shape and physical properties of the geological body can be inferred from the form of the scattered waves.

[0004] In order to obtain accurate and reliable scattered wave imaging results, accurate and efficient scattered wave numerical simulation methods are indispensable. However, there are still many technical problems to be solved in the current research on scattered wave numerical simulation: First, the existing scattered wave numerical simulation technologies mostly use approximate processing for boundary conditions, which is significantly different from the actual situation; second, the existing scattered wave numerical simulation methods have huge requirements for computing performance and memory, making it difficult to carry out large-scale and precise numerical simulations; third, it is difficult to achieve both efficient solution and high-precision numerical simulation of the seismic scattering wave field of the earth's medium. Summary of the invention

[0005] In order to overcome the deficiencies and defects mentioned in the above background technology, the present invention provides a holographic numerical simulation method, system and storage medium for earth medium acoustic wave scattering wave field.

[0006] In order to solve the above technical problems, the technical solution proposed by the present invention is: In a first aspect, the present application provides a method for holographic numerical simulation of earth medium acoustic wave scattering wave field, comprising: S1: Setting the time domain acoustic wave equation for the characteristics of the earth medium; S2: Based on velocity field separation, background field and scattering field separation, the time domain acoustic wave equation is divided into a background wave field equation and a scattering wave field equation; and based on the scattering wave field equation, a time domain scattering wave field wave equation is obtained; S3: Determine a group of spatial wavenumber mixed domain scattering wave equations in the frequency domain based on the time domain scattering wave field wave equation; S4: Calculate the spatial domain scattering wave field distribution based on the spatial wave number mixed domain scattering wave equation group in the frequency domain as the result of holographic numerical simulation of the earth medium acoustic wave scattering wave field.

[0007] In a second aspect, the present application provides a holographic numerical simulation system for the earth medium acoustic wave scattering wave field, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor implements the steps of the method described in the first aspect when executing the computer program.

[0008] In a third aspect, the present application provides a computer storage medium, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor implements the steps of the method described in the first aspect when executing the computer program.

[0009] Compared with the prior art, the present invention has the following beneficial effects: 1. The present invention establishes an acoustic wave numerical simulation model according to the characteristics of the earth medium, ensuring the consistency between the external boundary of the model and the real boundary conditions, thereby realizing accurate holographic numerical simulation of the earth medium acoustic wave scattering wave field; 2. The function product Fourier transform conjecture proposed in the present invention can transform the problem of solving three-dimensional partial differential equations into the problem of solving a series of one-dimensional ordinary differential equations, greatly reducing the amount of simulation calculations and memory requirements; 3. The present invention is based on iterative algorithm and holographic Fourier transform method, and realizes accurate numerical simulation calculation on the basis of efficient solution of seismic wave simulation. BRIEF DESCRIPTION OF THE DRAWINGS

[0010] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0011] Figure 1 It is a flow chart of a holographic numerical simulation method of earth medium acoustic wave scattering wave field provided by a preferred embodiment of the present application; Figure 2 It is a schematic diagram of a two-dimensional calculation model without air provided in a preferred embodiment of the present application; Figure 3 It is a total field comparison diagram of the present invention and the finite element method provided in the preferred embodiment of the present application; Figure 4 It is a total field comparison curve of the present invention and the finite element method at x=0m provided in the preferred embodiment of the present application; Figure 5It is a comparison diagram of the present invention and the finite element scattered field provided by the preferred embodiment of the present application; Figure 6 It is a comparison curve of the scattered fields of the present invention and the finite element method at x=0m provided in the preferred embodiment of the present application; FIG7 is an iterative convergence curve provided by a preferred embodiment of the present application; Figure 8 It is a schematic diagram of a three-dimensional calculation model without air provided in a preferred embodiment of the present application; Fig. 9 It is a total field comparison diagram of the present invention and the finite element method at y=0m provided in the preferred embodiment of the present application; Fig.10 It is a total field comparison curve of the present invention and the finite element method at the y=0m, x=0m section provided in the preferred embodiment of the present application; Fig.11 It is a comparison diagram of the y=0m cross section of the present invention and the finite element scattered field provided by the preferred embodiment of the present application; Fig.12 It is a comparison curve of the scattered fields of the present invention and the finite element method at the y=0, x=0m section provided in the preferred embodiment of the present application; FIG13 is an iterative convergence curve diagram of a three-dimensional model provided in a preferred embodiment of the present application; FIG14 is a schematic diagram showing a comparison of the efficiency and memory usage of a two-dimensional model provided by a preferred embodiment of the present application and a conventional finite element method; FIG15 is a schematic diagram showing a comparison of the efficiency and memory usage of the three-dimensional model provided by the preferred embodiment of the present application and the conventional finite element method; Fig.16 It is a schematic diagram of a two-dimensional air combination model provided by a preferred embodiment of the present application; FIG17 is a calculation result of a two-dimensional combination model provided in a preferred embodiment of the present application; Fig.18 It is a schematic diagram of a three-dimensional air combination model provided by a preferred embodiment of the present application; FIG19 is a calculation result of a three-dimensional combined model provided in a preferred embodiment of the present application; Fig. 20 It is a y=0m wave field profile provided by the preferred embodiment of the present application; FIG. 21 is a wave field plane diagram at z=0m and z=250m provided in a preferred embodiment of the present application. DETAILED DESCRIPTION

[0012] In order to facilitate the understanding of the present invention, the present invention will be described more comprehensively and meticulously below in conjunction with the accompanying drawings and preferred embodiments of the present invention, but the protection scope of the present invention is not limited to the following specific embodiments.

[0013] Unless otherwise defined, all professional terms used below have the same meanings as those generally understood by those skilled in the art. The professional terms used herein are only for the purpose of describing specific embodiments and are not intended to limit the scope of protection of the present invention.

[0014] Unless otherwise defined, the technical terms or scientific terms used in the present invention shall have the usual meanings understood by persons with ordinary skills in the field to which the present invention belongs. The words "first", "second" and similar words used in the present invention do not indicate any order, quantity or importance, but are only used to distinguish different components. Similarly, words such as "one" or "one" do not indicate quantity restrictions, but indicate the existence of at least one. Words such as "connect" or "connected" are not limited to physical or mechanical connections, but may include electrical connections, whether direct or indirect. "Up", "down", "left", "right" and the like are only used to indicate relative positional relationships. When the absolute position of the object being described changes, the relative positional relationship also changes accordingly.

[0015] It should be understood that in the prior art, it is a major problem to achieve high-precision numerical simulation of the seismic scattering wave field of the earth medium on the basis of efficient solution. Based on this, the present application provides a holographic numerical simulation method for the earth medium acoustic wave scattering wave field, which can be used for holographic numerical simulation of two-dimensional and three-dimensional acoustic wave scattering wave fields of the earth medium; wherein, by means of holographic Fourier transform, on the one hand, it is ensured that the outer boundary conditions of the numerical simulation area are completely consistent with the real physical boundary conditions, so that the real spectrum of the scattering wave field can be obtained, and on the other hand, the spectral information corresponding to all wave numbers is fully utilized, so that the numerical simulation model of the scattering wave field is consistent with the seismic geological physical model of the earth medium, so that the kinematic characteristics and dynamic characteristics of the seismic wave can be accurately simulated, hence the name holographic numerical simulation.

[0016] See also Figure 1 The present application provides a holographic numerical simulation method of earth medium acoustic wave scattering wave field, comprising: S1: Setting the time domain acoustic wave equation for the characteristics of the earth medium; S2: Based on velocity field separation, background field and scattering field separation, the time domain acoustic wave equation is divided into a background wave field equation and a scattering wave field equation; and based on the scattering wave field equation, a time domain scattering wave field wave equation is obtained; S3: Determine a group of spatial wavenumber mixed domain scattering wave equations in the frequency domain based on the time domain scattering wave field wave equation; S4: Calculate the spatial domain scattering wave field distribution based on the spatial wave number mixed domain scattering wave equation group in the frequency domain as the result of holographic numerical simulation of the earth medium acoustic wave scattering wave field.

[0017] The above-mentioned holographic numerical simulation method of the earth medium acoustic wave scattering wave field establishes an acoustic wave numerical simulation model according to the characteristics of the earth medium, ensures the consistency of the external boundary of the model with the actual boundary conditions, and thus realizes accurate holographic numerical simulation of the earth medium acoustic wave scattering wave field.

[0018] It is worth pointing out that in the present application, a numerical simulation model of sound waves is established according to the characteristics of the earth's medium. There is an air layer above the ground, and the external boundary of the model, i.e. the natural boundary, is completely consistent with the real geophysical boundary conditions. Therefore, the numerical simulation model of sound waves is consistent with the seismic geophysical model of the earth's medium, and the kinematics and dynamics of seismic waves can be accurately simulated.

[0019] The following is a complete example to describe in detail the steps of the above-mentioned holographic numerical simulation method for the earth medium acoustic wave scattering wave field: Step 1: Give the time domain acoustic wave equation The propagation of sound waves in a medium satisfies the following wave equation: (1); in, Represents the wave field, It represents the speed of sound wave propagation in m / s. is the source function, is the time in seconds, The viscosity term is added to take into account the energy dissipation of the wave during propagation. This term is 0 when there is no attenuation.

[0020] Step 2: Velocity field separation It should be pointed out that in numerical simulation of wave fields, the total velocity field is regarded as the superposition of the background velocity field and the abnormal velocity field, where the velocity of the background velocity field is expressed as The total velocity field is expressed as Indicates that this application proposes the following relationship between the two: (2); in, is the physical quantity that characterizes the abnormal velocity field and is a function of the spatial coordinates. =0 The total velocity field is consistent with the background velocity field, The total velocity field is smaller than the background velocity field. The total velocity field is larger than the background velocity field.

[0021] Step 3: Separation of background field and scattered field After the velocity field is separated, the acoustic wave field can also be regarded as the superposition of the background field and the scattered field, that is: (3); in, is the scattered field, For the background field.

[0022] The total velocity field is replaced by the sum of the abnormal velocity and the background velocity, and the acoustic wave field is replaced by the sum of the background field and the scattered field. That is, equations (2) and (3) are substituted into equation (1), and we get: (4); In the formula, is the viscosity term related to energy attenuation during wave propagation. represents the Laplace operator, , represents partial derivative, is the source function, For time; Split into the following two formulas: Background wave field equation: (5); Scattering wave field equation: (6); make , known but Known, let , formula (6) can be written as: (7); The above equation (7) is the wave equation satisfied by the time domain scattering wave field.

[0023] Step 4: Background field calculation The background velocity models are all full space, half space or uniform layered medium, and the sources can be plane waves or point source acoustic waves, both of which satisfy the background field equation (5). In the full space case, the background field of plane waves and point sources can directly give the results; in layered media, the wave fields of plane waves and point sources can be directly obtained using analytical methods.

[0024] Specifically, in the full space case where the meaning of the result is given directly, there are analytical formulas for plane waves and point sources, and their distribution can be directly known. Among them, the solution for plane waves is: ; In the formula, represents amplitude, i represents imaginary number, and t represents time; The solution for a point source is: ; In the formula, k represents the wave number, k=2πf / c, c refers to the wave velocity, represents the frequency domain Ricker wavelet, which satisfies the following relationship; ; Where f represents the frequency of the Ricker wavelet, fm is the main frequency of the Ricker wavelet.

[0025] Step 5: Scattering wave ordinary differential equation in the frequency domain with spatial wavenumber mixing Perform Fourier transform on the time term of equation (7) and convert it to the frequency domain to obtain the equation (8); in, is the scattered wave in the frequency domain, ω is the angular frequency, is the frequency domain source function, represents the intermediate variable, is an imaginary number, It is a viscosity term related to energy attenuation during wave propagation.

[0026] For the above formula (8), The Fourier transform of the x and y directions (in the two-dimensional case) or the two-dimensional Fourier transform of the x and y directions (in the three-dimensional case) is converted to the wave number domain, where About The function of The Fourier transform of the multiplication is Represented. So far, the spatial wave number mixed domain scattered wave ordinary differential equation is obtained in the frequency domain: (9); in, About The function of express and The Fourier transform of the multiplication, is the scattered wave in the frequency domain, z is the vertical coordinate, is the spatial wavenumber domain source in the frequency domain, is the wave number corresponding to the x direction, is the wave number corresponding to the y direction. In the two-dimensional case, .

[0027] Step 6: Convert the ODE system into a series of ODEs Because there are many groups of kx and ky, one group of kx and ky corresponds to one equation, and all of them are added together to form a system of equations. In formula (9), According to the convolution theorem, the Fourier transform of , the spectrum of the scattered wave field at different wave numbers Therefore, Equation (9) is actually a set of ordinary differential equations corresponding to the wave number. By directly solving this set of ordinary differential equations, we can obtain the scattering wave field spectra corresponding to all wave numbers. However, this is an overall problem that is not easy to split and therefore difficult to parallelize. The computational time and memory requirements are similar to those of conventional methods such as finite element and finite difference methods. For large-scale and precise numerical simulation problems, it is difficult to implement due to the large computational time and memory requirements.

[0028] To solve this problem, this application proposes the function product Fourier transform conjecture, which is specifically: for the velocity anomaly area, for any wave number kx, ky, there is always a wave function , making the abnormal speed and scattered wave field The Fourier transform of the product is equal to the wave function and scattered wave field spectrum Product, ignoring variables, the mathematical expression of this conjecture is (10); Substituting formula (10) into formula (9), formula (9) can be written as: (11); Formula (11) is the ordinary differential equation corresponding to the wave numbers kx and ky. Different wave numbers are independent of each other. This realizes the transformation of the large-scale three-dimensional partial differential equation solution problem into the one-dimensional ordinary differential equation solution problem corresponding to a series of wave numbers kx and ky. The ordinary differential equations corresponding to different wave numbers are easy to solve in parallel and the amount of calculation is very small, which can greatly reduce the amount of calculation and memory requirements.

[0029] Step 7: Solution of ordinary differential equations.

[0030] make , formula (11) is written as: (12); There is usually no abnormal velocity in the upper and lower boundary areas of the solution region, so =0, the general solution of this equation can be written as ,in is a downward wave, At the upper boundary Zmin, there are only upward waves, and at the lower boundary Zmax, there are only downward waves, that is, (13); So the boundary conditions are: (14); The boundary value problem is: (15); For the boundary value problem corresponding to a series of wave numbers kx and ky in equation (15), the one-dimensional finite element method with quadratic interpolation is used to obtain a five-diagonal linear equation system. The five-diagonal linear equation system is quickly solved by the pursuit method to obtain the scattering wave field in the wave number domain: The distribution of the scattered wave field in the spatial domain can be obtained by inverse Fourier transform. This solution can fully consider the calculation accuracy and efficiency, thereby realizing efficient and high-precision numerical simulation of large-scale complex models.

[0031] Step 8: Iterate to obtain the exact solution of the scattered wave field.

[0032] This embodiment establishes a set of iterative algorithms, and the scattered wave field gradually approaches the exact solution through iterative calculation.

[0033] Step 1: According to formula (9), the one-dimensional ordinary differential equation system in the wave number domain of arbitrary frequency space is: (16); To generate the initial value, the number of iterations , according to the Born approximation, let , (16) is transformed into an ordinary differential equation corresponding to the wave numbers kx and ky, which can be solved to obtain the Born approximate solution . Step 2: Start iteration, ; Based on the hypothesis proposed in this application , can be obtained .

[0034] Step 3: According to the conjecture (10) proposed by the present invention, the wave function Substituting them into formula (16), we get Ordinary differential equations that satisfy (17); Solving equation (17), we obtain .

[0035] Step 4: For all wave numbers kx, ky, determine whether the following equation holds true: (18); in is a decimal close to 0. If (18) holds, For an exact solution to the equation, exit the iteration. The smaller it is, the closer the solution of equation (17) is to the solution of the original equation (16), and the higher the accuracy of the solution. If equation (18) does not hold, let the revised solution be the current solution plus the disturbance , that is, the corrected solution is , proceed to the next step.

[0036] Step 5: Will Substituting into (16), we get (19); Subtracting (19) from (17), we obtain the perturbation solution One-dimensional ordinary differential equations (20); Step 6: right The one-dimensional ordinary differential equations satisfying the Born approximation are solved. According to the Born approximation, let the left side of the equation is 0, and the solution is According to the conjecture (10) proposed by the present invention, there exists a , so that the following equation holds (twenty one); In the formula, represents the wave function of the disturbance field; Obtain , replace the formula (20) use Substituting variables, we get a one-dimensional ordinary differential equation: (twenty two); Solved , then the corrected solution is , return to Step 2 to start the next iteration until (18) Established.

[0037] This patent proposes the function product Fourier transform conjecture, which uses this conjecture to transform the problem of solving three-dimensional partial differential equations into the problem of solving one-dimensional ordinary differential equations corresponding to a series of wave numbers kx and ky. Ordinary differential equations corresponding to different wave numbers are easy to solve in parallel and the amount of calculation is very small, which can greatly reduce the amount of calculation and memory requirements. With the help of holographic Fourier transform, on the one hand, it ensures that the outer boundary conditions of the numerical simulation area are completely consistent with the real physical boundary conditions, so that the real spectrum of the scattered wave field can be obtained. On the other hand, the spectral information corresponding to all wave numbers is fully utilized, so that the holographic numerical simulation model of the scattered wave field is consistent with the seismic geophysical model of the earth medium, so the kinematic and dynamic characteristics of the seismic wave can be accurately simulated, hence the name holographic numerical simulation.

[0038] The acoustic scattering wave holographic numerical simulation algorithm proposed in the present invention is tested below.

[0039] The test computer configuration is 12th Gen Intel(R) Core(TM) i9-12900KS, 3.40 GH, and 128GB of memory.

[0040] 1. Correctness verification In order to verify the correctness of the algorithm proposed in the present invention, the commercial software COMSOL Multiphysics is used to perform numerical simulation of the sound waves, and the finite element numerical solution is used as a reference solution to verify the correctness of the algorithm of this patent.

[0041] A. Two-dimensional The calculation model is verified for the two-dimensional case. Figure 2 As shown, the background velocity is 3000m / s, α is 0.1, that is, the abnormal velocity is 2727m / s, the background sound wave is a plane wave, the frequency is 15Hz, the horizontal direction subdivision node Nx and the vertical direction subdivision node Nz are both 201, the spatial domain is uniformly sampled, the horizontal direction wave number Nkx is 401, the wave number selection range is -0.5~0.5, and the wave number domain is uniformly sampled. After four iterations, the convergence condition is reached, and the calculation time is 1.5s.

[0042] The total field calculated by the present invention and the total field calculated by the finite element method and their errors are as follows: Figure 3 As shown, the relative root mean square error of the two is less than 1%, which verifies the correctness and accuracy of the algorithm proposed in the present invention. Figure 4 This is a waveform comparison diagram in the vertical z direction when the calculation model is x=0m. The two are almost completely overlapped.

[0043] The scattered field calculated by the present invention and the total field calculated by the finite element method and their errors are as follows: Figure 5 shown. Figure 6 It is the comparison curve of the scattered field of the present invention and the finite element method at x=0m in the calculation model.

[0044] The iterative convergence calculated by this algorithm is characterized by relative error, as follows: ; The iteration curve is shown in Figure 7, which shows that the algorithm can converge stably.

[0045] In summary, the correctness and accuracy of the algorithm proposed in the present invention can be verified.

[0046] B. Three-dimensional Design a three-dimensional prism model such as Figure 8 As shown in the figure, the size of the numerical simulation study area is 1000m×1000m×500m, the background sound wave adopts plane wave, the background velocity is 3000m / s, and the calculation frequency is 15 Hz; the parameters characterizing the abnormal velocity If the value is 0.1, the velocity of the anomaly is 2727m / s. The anomaly is located in the center of the model and is a cube anomaly with a side length of 100m. The spatial domain grid division Nx, Ny and Nz are 201, uniform sampling; the wave number Nkx and Nky are 401, the wave number selection is -0.06~0.06, and uniform sampling. After 3 iterations, the convergence condition is reached, and the calculation time is 1.5s.

[0047] Take the section y=0m, the total field calculated by the present invention and the total field calculated by the finite element method and their errors are as follows: Fig. 9 As shown, the relative root mean square error of the two is less than 1%, which verifies the correctness and accuracy of the algorithm proposed in the present invention in three-dimensional conditions. Fig.10 This is a waveform comparison diagram in the vertical z direction when the calculation model is y=0m, x=0m. The two are almost completely overlapped.

[0048] At the section y=0m, the scattered field calculated by the present invention and the total field calculated by the finite element method and their errors are as follows: Fig.11 As shown, it can be seen that the two are very close, verifying the correctness and accuracy of the algorithm proposed in the present invention. Near the boundary, the scattered field obtained by the finite element method still has boundary effects, while the result of the present invention has no boundary effects, which is consistent with the actual situation. Fig.12 This is a waveform comparison diagram in the vertical z direction when the calculation model is y=0m, x=0m. The two are basically the same.

[0049] The iterative convergence curve of this model calculation in three-dimensional case is shown in Figure 13. It can be seen that the algorithm can also converge stably in three-dimensional case.

[0050] In summary, the correctness and accuracy of the algorithm proposed in the present invention can be verified in two-dimensional and three-dimensional situations.

[0051] 2. Efficiency The algorithm efficiency is illustrated by comparing it with the conventional finite element method.

[0052] The change factor of the total number of grid nodes is calculated as , the calculation time changes by a factor of , which satisfies the relationship , where a has different values ​​in different algorithms.

[0053] The algorithm efficiency statistics for the two-dimensional model are as follows: Table 1 Statistics of iteration time and memory usage of the proposed algorithm and COMSOL finite element algorithm in two-dimensional case

[0054] The two-dimensional model and the Comsol finite element method are compared under the same mesh generation conditions. , finite element method That is, when the number of nodes increases by 4 times, the calculation time of the algorithm proposed in this patent increases by 2.3 times, which is less than 4 times, which is a sublinear growth, while the time increased by the finite element method is 4.2 times, which is greater than 4 times, which is an exponential growth. It can be seen that the calculation time of the algorithm of the present invention changes sublinearly with the number of nodes, while the calculation time of the finite element method increases exponentially with the increase in the number of nodes. Therefore, in the two-dimensional case, the more nodes the present invention has, the more obvious the algorithm advantage.

[0055] The efficiency statistics of the 3D model are as follows: Table 2 Statistics of iteration time and memory usage of the proposed algorithm and COMSOL finite element algorithm in three-dimensional case

[0056] Among them, the comparison of efficiency and memory usage between the two-dimensional model and the conventional finite element method is shown in Figure 14, and the comparison of efficiency and memory usage between the three-dimensional model and the conventional finite element method is shown in Figure 15.

[0057] The three-dimensional model and the Comsol finite element method are compared under the same mesh generation conditions. , finite element method That is, when the number of nodes increases by 8 times, the calculation time of the algorithm proposed in this patent increases by 5.27 times, which is less than 8 times, which is a sublinear growth, while the time increased by the finite element method is about 22 times, which is greater than 8 times, which is an exponential growth. Therefore, the algorithm advantage of the present invention is more obvious in the three-dimensional case, and the calculation efficiency of the small scale (grid 81×81×41) is more than 3 orders of magnitude higher. As the scale of numerical simulation increases, the calculation efficiency advantage of the present invention increases rapidly.

[0058] Therefore, whether in two-dimensional or three-dimensional cases, the memory usage and calculation time of the algorithm proposed in the present invention are much lower than those of the traditional finite element method, and the larger the calculation scale, the more obvious the advantages of the present invention.

[0059] Model Example 2D: Design as Fig.16 The model with air shown has an air layer of 300m and an underground layer of 1000m. The air layer velocity is 340m / s, the background field is a plane wave, and the background velocity is 3000m / s; there are two abnormal bodies, a low-speed body with α of 0.1 and a velocity of 2727m / s, and a high-speed body with α of -0.2 and a velocity of 3750m / s; the number of nodes in the x-direction is 501, the number of nodes in the z-direction is 501, and the wave number is selected to be 501, and both the spatial domain and the wave number domain are uniformly sampled. The calculation frequency is 15Hz. The convergence standard is reached after three iterations.

[0060] Figure 17 below is the calculation result of the two-dimensional combined model. Figure 17-1 is the plane wave background field, and Figure 17-2 is the scattered field. The waveform is clear, the abnormal difference is obvious, and the greater the difference with the background speed, the stronger the abnormal field response; Figure 17-3 is the total field. It can be seen that the wave field under the low-speed body can pass through, and the wave field under the high-speed body cannot pass through. It can be seen that the high-speed body has shielding properties for the wave field. From the two-dimensional combined model with air, it can be seen that the present invention can be calculated with air in the two-dimensional case, and the model is real.

[0061] 3D: Design as Fig.18 The model with air shown has an air layer of 200m and an underground layer of 500m; the air layer velocity is 340m / s, the background field is a plane wave, and the background velocity is 3000m / s; there are two abnormal bodies, a low-speed body with α of 0.1 and a velocity of 2727m / s, and a high-speed body with α of -0.2 and a velocity of 3750m / s, and the negative x direction is a low-speed body. The number of nodes in the x and y directions is 501, the number of nodes in the z direction is 501, and the wave number is selected to be 501, and the spatial domain and the wave number domain are uniformly sampled. The calculation frequency is 15Hz. After three iterations, the convergence standard is reached.

[0062] Figure 19 below is the calculation result of the three-dimensional combined model, 19-1 is the plane wave background field, 19-2 is the scattered field, which is the same as the two-dimensional case. The greater the difference between the abnormal velocity and the background velocity, the stronger the response; 19-3 is the total field, which is also consistent with the two-dimensional case, that is, the wave field under the low-speed body can pass through, and the wave field under the high-speed body cannot pass through. It can be seen that the high-speed body has shielding properties for the wave field. From the three-dimensional combined model with air, it can be seen that the present invention can be calculated with air in the three-dimensional case, the model is real, and can adapt to large-scale, efficient and high-precision numerical calculations.

[0063] Down Fig. 20 The wave field profile at y=0m and its logarithmic display show that at the transition interface between air and underground medium, the wave field is continuous and can simulate the real situation of considering air. From the logarithmic diagram, it can be seen that there are also abnormal wave fields in the air. If the air is not considered, the wave field situation in the air cannot be known. Therefore, it is once again explained that the present invention starts from the real model and the necessity of considering the air layer.

[0064] The following Figure 21 is a wave field plane diagram of the ground and the center of the anomaly, which can still clearly reflect the existence of the anomaly. The wave field amplitude is large at the anomaly, and it becomes weaker when it reaches the ground, which is in line with the laws of physics. This once again illustrates the correctness and rationality of the algorithm proposed in the present invention.

[0065] It can be seen that the present invention has developed a numerical simulation method that can correctly and accurately simulate the acoustic wave field of the real earth medium. The method is correct, the algorithm converges stably, and the efficiency of small-scale numerical simulation is 3 to 5 orders of magnitude higher than that of conventional methods. The memory occupied is reduced by 2 to 4 orders of magnitude compared with conventional methods. The larger the calculation scale, the more obvious the calculation advantage. The present invention has developed a numerical simulation method for high-efficiency and high-precision real earth medium acoustic wave scattering waves suitable for large-scale arbitrary complex conditions, that is, a holographic numerical simulation method for acoustic wave scattering waves.

[0066] The present application also provides a holographic numerical simulation system for earth medium acoustic wave scattering wave field, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor implements the steps of the above method when executing the computer program. The holographic numerical simulation system for earth medium acoustic wave scattering wave field can implement various embodiments of the holographic numerical simulation method for earth medium acoustic wave scattering wave field, and can achieve the same beneficial effects, which will not be described in detail here.

[0067] The present application also provides a computer storage medium, including a memory, a processor, and a computer program stored in the memory and executable on the processor, and the processor implements the steps of the above method when executing the computer program. The computer storage medium can implement various embodiments of the holographic numerical simulation method of the earth medium acoustic wave scattering wave field, and can achieve the same beneficial effects, which will not be described in detail here.

[0068] The preferred specific embodiments of the present invention are described in detail above. It should be understood that a person skilled in the art can make many modifications and changes based on the concept of the present invention without creative work. Therefore, any technical solution that can be obtained by a person skilled in the art through logical analysis, reasoning or limited experiments based on the concept of the present invention on the basis of the prior art should be within the scope of protection determined by the claims.

Claims

1. A holographic numerical simulation method for the earth medium acoustic wave scattering wave field, characterized in that: include: S1: Setting the time domain acoustic wave equation for the characteristics of the earth medium; S2: Based on velocity field separation, background field and scattering field separation, the time domain acoustic wave equation is divided into a background wave field equation and a scattering wave field equation; and based on the scattering wave field equation, a time domain scattering wave field wave equation is obtained; S3: Determine a group of spatial wavenumber mixed domain scattering wave equations in the frequency domain based on the time domain scattering wave field wave equation; S4: Calculate the spatial domain scattering wave field distribution based on the spatial wave number mixed domain scattering wave equation group in the frequency domain as the result of holographic numerical simulation of the earth medium acoustic wave scattering wave field.

2. The method for holographic numerical simulation of earth medium acoustic wave scattering wave field according to claim 1, characterized in that: The time domain acoustic wave equation in S1 satisfies the following relationship: （1）； in, represents the acoustic wave field, It represents the speed of sound wave propagation in m / s. is the source function, is the time in seconds, is the viscosity term related to energy attenuation during wave propagation, where represents the Laplace operator, , It means partial derivative.

3. The method for holographic numerical simulation of earth medium acoustic wave scattering wave field according to claim 2, characterized in that S2 include: S21: Separate the total velocity of the acoustic wave into the background velocity field and the abnormal velocity field, satisfying the following relationship: （2）； in, is the total speed of sound wave propagation, represents the background velocity field, is the spatial coordinate function representing the abnormal velocity field; S22: Separate the acoustic wave field after velocity field separation into background field and scattered field, satisfying the following relationship: （3）； in, is the acoustic wave field, is the scattered field, For the background field; S23: According to equations (2) and (3), equation (1) is split into the background wave field equation and the scattered wave field equation. The background wave field equation is as follows: （4）； In the formula, is the viscosity term related to energy attenuation during wave propagation. represents the Laplace operator, , represents partial derivative, is the source function, For time; The scattered wave field equation is as follows: （5）； S24: Order ,make , is the intermediate variable, represents the intermediate variable of the background field, and transforms Equation (5) into the wave equation of the scattered wave field in the time domain: （6）。 4. The method for holographic numerical simulation of earth medium acoustic wave scattering wave field according to claim 1, characterized in that: The S3 includes: S31: Perform Fourier transform on the time term t in the time domain scattering wave field wave equation and convert it into the frequency domain equation as follows: （7）； in, is the scattered wave in the frequency domain, is the angular frequency, is the frequency domain source function, represents the intermediate variable, is an imaginary number, It is the viscosity term related to energy attenuation during wave propagation; S32: Perform a two-dimensional Fourier transform of equation (7) in the x and y directions to the wave number domain, and obtain the scattering wave equation in the spatial wave number mixing domain as follows: （8）； in, About The function of express and The Fourier transform of the multiplication, is the scattered wave in the frequency domain, z is the vertical coordinate, is the spatial wavenumber domain source in the frequency domain, is the wave number corresponding to the x direction, is the wave number corresponding to the y direction, in the two-dimensional case .

5. The holographic numerical simulation method of earth medium acoustic wave scattering wave field according to claim 4, characterized in that: The S4 includes: S41: Combine multiple sets of spatial wave number mixed domain scattering wave equations into a set of spatial wave number mixed domain scattering wave equations, and transform them into a series of ordinary differential equations; S42: Solving the series of ordinary differential equations based on a set iterative method to obtain the spatial domain scattering wave field distribution, wherein the step of setting the iterative method includes setting the function product Fourier transform conjecture.

6. The method for holographic numerical simulation of earth medium acoustic wave scattering wave field according to claim 5, characterized in that: The series of ordinary differential equations in S41 satisfy the following relationship: （9）。 7. The method for holographic numerical simulation of earth medium acoustic wave scattering wave field according to claim 6, characterized in that: The S42 includes: S421: According to equation (9), the one-dimensional ordinary differential equation system in the wave number domain of arbitrary frequency space is: （10）； Generate the initial value for equation (10), let , according to the Born approximation, let , transform equation (9) into an ordinary differential equation corresponding to the wave numbers kx, ky, and solve it to obtain the Born approximate solution ; S422: Order , start iterative calculation, according to the set function product Fourier transform conjecture , and obtain the wave function ; S423: According to the set function product Fourier transform conjecture, the wave function Substituting them into formula (10), we get The ordinary differential equations satisfied are as follows: （11）； Solve equation (11) and get the solution ; S524: For all wave numbers kx, ky, find , and determine whether the following formula is true: （12）； in, is a specified decimal number close to 0; If equation (12) holds true, then the solution The iterative calculation is ended to obtain the exact solution of the equation and the result of the holographic numerical simulation of the earth medium acoustic wave scattering wave field is obtained; If equation (12) does not hold, let the revised solution be the current solution plus the perturbation solution , that is, let the corrected solution be , continue to the next step of calculation; S425: Substituting into formula (10), we get: （13）； Subtracting (13) from (11) yields the perturbation solution The one-dimensional ordinary differential equation system: （14）； S426: Using Born approximation to solve equation (14) we get , according to the function product Fourier transform conjecture, there exists a , so that the following holds: （15）； In the formula, represents the wave function of the disturbance field; Obtain , replace the formula (15) use Substituting variables, we get a one-dimensional ordinary differential equation: （16）； Find the perturbation solution , let the corrected solution be , return to S522 to perform the next iterative calculation.

8. A holographic numerical simulation system for earth medium acoustic wave scattering wave field, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When a processor executes the computer program, the steps of any one of the methods of claims 1 to 7 are implemented.

9. A computer storage medium comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When a processor executes the computer program, the steps of any one of the methods of claims 1 to 7 are implemented.

Images