Holographic numerical simulation method and system for acoustic scattered wave field in earth medium, and storage medium

Abstract

The present invention relates to the field of seismic exploration. Disclosed are a holographic numerical simulation method and system for an acoustic scattered wave field in an earth medium, and a storage medium. The method comprises: setting a time domain acoustic wave equation for 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 on the basis of separation of velocity fields and separation of background and scattered fields, and obtaining a time domain scattered wave field wave equation on the basis of the scattered wave field equation; determining a system of space-wavenumber mixed domain scattered wave equations in a frequency domain on the basis of the time domain scattered wave field wave equation; and on the basis of the system of space-wavenumber mixed domain scattered wave equations in the frequency domain, calculating a spatial domain scattered wave field distribution as a result of holographic numerical simulation for an acoustic scattered wave field in the earth medium. In this way, an acoustic wave numerical simulation model is established on the basis of the characteristics of the earth medium, which ensures the consistency between external boundary conditions of the model and real boundary conditions, thereby realizing accurate holographic numerical simulation for an acoustic scattered wave field in an earth medium.

Description

A method, system, and storage medium for holographic numerical simulation of acoustic wave scattering wavefield in Earth's medium. Technical Field

[0001] This invention belongs to the field of seismic exploration, and in particular relates to a method, system and storage medium for holographic numerical simulation of acoustic wave scattering wavefield in the Earth's medium. Background Technology

[0002] Seismic exploration based on reflected waves is a widely used method in geological exploration. It can achieve good results when the geological structure is relatively homogeneous. However, as the application fields of seismic exploration continue to expand, the geological conditions of oil and mineral deposits are becoming more and more complex, and the non-homogeneity of underground media is constantly increasing. Only scattered wave seismic exploration technology for non-homogeneous and complex media can be applied to more complex geological conditions.

[0003] Scattered wave seismic exploration is based on scattering theory. In a non-homogeneous medium, seismic waves will be scattered when they encounter such a non-homogeneous geological body. The reflected wave of a layered medium is the superposition of the scattered waves generated by each scattering point on the layered medium. Different scales and different non-homogeneities result in different forms of scattered waves. Therefore, the shape and physical properties of the geological body can be inferred from the shape of the scattered waves.

[0004] To obtain accurate and reliable scattered wave imaging results, precise and efficient numerical simulation methods for scattered waves are essential. However, there are still many technical problems to be solved in the current research on numerical simulation of scattered waves: First, existing numerical simulation techniques for scattered waves mostly use approximations for boundary conditions, which differ significantly from the actual situation; Second, existing numerical simulation methods for scattered waves have huge computational and memory requirements, making it difficult to perform large-scale, detailed numerical simulations; Third, it is difficult to achieve both efficient solution of the seismic scattered wave field in the Earth's medium and high-precision numerical simulation. Summary of the Invention

[0005] To overcome the deficiencies and defects mentioned in the background art, the present invention provides a method, system, and storage medium for holographic numerical simulation of acoustic wave scattering wavefield in the Earth's medium.

[0006] To solve the above-mentioned technical problems, the technical solution proposed by this invention is as follows:

[0007] In a first aspect, this application provides a holographic numerical simulation method for the acoustic wave scattering wavefield of the Earth's medium, comprising:

[0008] S1: Define the time-domain acoustic wave equations specific to the characteristics of the Earth's medium;

[0009] S2: Based on velocity field separation, background field and scattered field separation, the time-domain acoustic wave equation is divided into background wave field equation and scattered wave field equation; and based on the scattered wave field equation, the time-domain scattered wave field wave equation is obtained.

[0010] S3: Determine the set of spatial wavenumber mixing domain scattering wave equations in the frequency domain based on the time-domain scattering wave field wave equations;

[0011] S4: The spatial domain scattered wave field distribution is calculated based on the spatial wavenumber mixing domain scattering wave equations in the frequency domain as the result of the holographic numerical simulation of the acoustic wave scattered wave field in the Earth medium.

[0012] In a second aspect, this application provides a holographic numerical simulation system for the acoustic wave scattering wave field of the Earth's medium, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the steps of the method described in the first aspect.

[0013] Thirdly, this application provides a computer storage medium, including a memory, a processor, and a computer program stored on the memory and executable on the processor, wherein the processor executes the computer program to implement the steps of the method described in the first aspect.

[0014] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0015] 1. Based on the characteristics of the Earth's medium, this invention establishes a numerical simulation model for acoustic waves, ensuring the consistency between the model's external boundary and the actual boundary conditions, thereby achieving accurate holographic numerical simulation of the acoustic wave scattering field in the Earth's medium.

[0016] 2. The function product Fourier transform conjecture proposed in this invention can transform the problem of solving three-dimensional partial differential equations into a problem of solving a series of one-dimensional ordinary differential equations, which greatly reduces the amount of computation and memory requirements of the simulation.

[0017] 3. Based on iterative algorithms and holographic Fourier transform methods, this invention achieves accurate numerical simulation calculations on the basis of efficient solution of seismic wave simulation. Attached Figure Description

[0018] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0019] Figure 1 is a flowchart of a holographic numerical simulation method for acoustic wave scattering wavefield in the Earth medium provided by a preferred embodiment of this application;

[0020] Figure 2 is a schematic diagram of a two-dimensional computational model without air provided in a preferred embodiment of this application;

[0021] Figure 3 is a comparison diagram of the total field of the present invention and the finite element method provided in the preferred embodiment of this application;

[0022] Figure 4 is a comparison curve of the total field of the present invention and the finite element method at x=0m provided in the preferred embodiment of this application;

[0023] Figure 5 is a comparison diagram of the scattering field of the present invention and the finite element method according to the preferred embodiment of this application;

[0024] Figure 6 is a comparison curve of the scattering field of the present invention and the finite element method at x=0m provided by the preferred embodiment of this application;

[0025] Figure 7 is an iterative convergence curve provided by a preferred embodiment of this application;

[0026] Figure 8 is a schematic diagram of a three-dimensional computational model without air provided in a preferred embodiment of this application;

[0027] Figure 9 is a comparison diagram of the total field of the present invention and the finite element method when y=0m is provided in the preferred embodiment of this application;

[0028] Figure 10 is a comparison curve of the total field of the present invention and the finite element method at the y=0m and x=0m cutoff lines provided in the preferred embodiment of this application.

[0029] Figure 11 is a comparison diagram of the scattering field of the present invention and the finite element method at y=0m cross section provided by the preferred embodiment of this application;

[0030] Figure 12 is a comparison curve of the scattering field of the present invention and the finite element method at the y=0, x=0m cutoff line provided in the preferred embodiment of this application;

[0031] Figure 13 is an iterative convergence curve of the three-dimensional model provided in the preferred embodiment of this application;

[0032] Figure 14 is a schematic diagram comparing the efficiency and memory usage of the two-dimensional model provided in the preferred embodiment of this application with that of the conventional finite element method;

[0033] Figure 15 is a schematic diagram comparing the efficiency and memory usage of the three-dimensional model provided in the preferred embodiment of this application with that of the conventional finite element method;

[0034] Figure 16 is a schematic diagram of a two-dimensional air-bearing combined model provided in a preferred embodiment of this application;

[0035] Figure 17 shows the calculation results of the two-dimensional combined model provided in the preferred embodiment of this application;

[0036] Figure 18 is a schematic diagram of a three-dimensional air-integrated model provided in a preferred embodiment of this application;

[0037] Figure 19 shows the calculation results of the three-dimensional combined model provided in the preferred embodiment of this application;

[0038] Figure 20 is a wave field profile at y=0m provided in a preferred embodiment of this application;

[0039] Figure 21 is a wave field planar diagram of z=0m and z=250m provided in the preferred embodiment of this application. Detailed Implementation

[0040] To facilitate understanding of the present invention, the present invention will be described more fully and in detail below with reference to the accompanying drawings and preferred embodiments, but the scope of protection of the present invention is not limited to the following specific embodiments.

[0041] Unless otherwise defined, all technical terms used herein have the same meaning as commonly understood by those skilled in the art. The technical terms used herein are for the purpose of describing particular embodiments only and are not intended to limit the scope of the invention.

[0042] Unless otherwise defined, the technical or scientific terms used in this invention shall have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in this invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Similarly, the terms "an" or "a" and similar terms do not indicate a quantity limitation, but rather indicate the presence of at least one. The terms "connected" or "linked" and similar terms are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. "Up," "down," "left," "right," etc., are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship also changes accordingly.

[0043] It should be understood that achieving high-precision numerical simulation of the seismic scattered wave field of the Earth's medium while efficiently solving the problem is a major challenge in the existing technology. Based on this, this application provides a holographic numerical simulation method for the acoustic wave scattered wave field of the Earth's medium, which can be used for two-dimensional and three-dimensional holographic numerical simulation of the acoustic wave scattered wave field of the Earth's medium. In this method, by using 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, thus obtaining the true spectrum of the scattered wave field. On the other hand, all spectral information corresponding to all wave numbers is utilized, so the numerical simulation model of the scattered wave field is consistent with the seismic geophysical model of the Earth's medium. Therefore, the kinematic and dynamic characteristics of seismic waves can be accurately simulated, hence the name holographic numerical simulation.

[0044] Please refer to Figure 1. This application provides a holographic numerical simulation method for acoustic wave scattering wavefield in the Earth's medium, comprising:

[0045] S1: Define the time-domain acoustic wave equations specific to the characteristics of the Earth's medium;

[0046] S2: Based on velocity field separation, background field and scattered field separation, the time-domain acoustic wave equation is divided into background wave field equation and scattered wave field equation; and based on the scattered wave field equation, the time-domain scattered wave field wave equation is obtained.

[0047] S3: Determine the set of spatial wavenumber mixing domain scattering wave equations in the frequency domain based on the time-domain scattering wave field wave equations;

[0048] S4: The spatial domain scattered wave field distribution is calculated based on the spatial wavenumber mixing domain scattering wave equations in the frequency domain as the result of the holographic numerical simulation of the acoustic wave scattered wave field in the Earth medium.

[0049] The above-mentioned holographic numerical simulation method for acoustic wave scattering wave field in the Earth medium establishes an acoustic wave numerical simulation model based on the characteristics of the Earth medium, ensuring the consistency between the external boundary of the model and the actual boundary conditions, thereby achieving accurate holographic numerical simulation of acoustic wave scattering wave field in the Earth medium.

[0050] It is worth noting that in this application, a numerical simulation model of acoustic waves is established based on the characteristics of the Earth's medium. The area above the ground is an air layer, and the external boundary of the model, i.e. the natural boundary, is completely consistent with the actual geophysical boundary conditions. Therefore, the numerical simulation model of acoustic waves is consistent with the seismic geological physical model of the Earth's medium, and the kinematics and dynamics of seismic waves can be accurately simulated.

[0051] The steps of the above-mentioned holographic numerical simulation method for acoustic wave scattering wavefield in the Earth's medium are described in detail below with a complete example:

[0052] Step 1: Give the time-domain acoustic wave equation

[0053] Sound wave propagation in a medium satisfies the following wave equation:

[0054] Where ψ represents the wave field, c represents the sound wave propagation speed in m / s, f is the source function, t is time in s, and β is a viscosity term added to account for energy dissipation during wave propagation; this term is 0 when there is no attenuation.

[0055] Step 2: Velocity Field Separation

[0056] It should be noted that in the numerical simulation of the wave field, the total velocity field is considered as the superposition of the background velocity field and the anomalous velocity field, where the background velocity field is represented by c0 and the total velocity field is represented by c. The relationship between the two is proposed in this application as follows:

[0057] Where α is a physical quantity characterizing the anomalous velocity field, which is a function of spatial coordinates. When α = 0, the total velocity field is the same as the background velocity field. When α > 0, the total velocity field is smaller than the background velocity field. When 0 ≥ α > -1, the total velocity field is larger than the background velocity field.

[0058] Step 3: Separation of background field and scattered field

[0059] After the velocity field is separated, the sound wave field can also be regarded as the superposition of the background field and the scattered field, that is: ψ=ψ0+ψ a (3);

[0060] Where, ψ a Let ψ0 be the scattering field and ψ0 be the background field.

[0061] Replacing the total velocity field with the sum of the anomalous velocity and the background velocity, and the sound wave field with the sum of the background field and the scattered field, that is, substituting equations (2) and (3) into equation (1), we get:

[0062] In the formula, β is a viscous term that considers the energy attenuation of the wave during propagation. Represents the Laplace operator. This indicates finding the partial derivative, f s Here, t is the source function;

[0063] It can be broken down into the following two formulas:

[0064] Background wave field equation:

[0065] Scattered wave field equation:

[0066] make Given that ψ0, then f0 is known. Let γ = 2α + α 2 Equation (6) can be written as:

[0067] Equation (7) above is the wave equation satisfied by the time-domain scattered wave field.

[0068] Step 4: Background Field Calculation

[0069] The background velocity models are all in full space, half space, or a homogeneous layered medium, and the sources can be either plane waves or point source acoustic waves, all of which satisfy the background field equation (5). In the case of full space, the background field of plane waves and point sources can be given directly; in the case of layered media, the wave field of plane waves and point sources can be obtained directly using analytical methods.

[0070] Specifically, in the case of the entire space where the meaning of the result is directly given, plane waves and point sources have analytical formulas, and their distributions can be directly determined. The solution for a plane wave is: p(t) = p0e^(-t / t). jωt ;

[0071] In the formula, p0 represents the amplitude, i represents the imaginary number, and t represents time;

[0072] The solution for the point source is:

[0073] In the formula, k represents the wave number, k = 2πf / c, c refers to the wave velocity, and g(ω) represents the frequency domain Ricker wavelet, which satisfies the following relationship;

[0074] Where f represents the frequency of the Ricker wavelet, f m It is the main frequency of the Reck wavelet.

[0075] Step 5: Ordinary Differential Equation of Scattered Waves in Spatial Wavenumber Mixing Domain in Frequency Domain

[0076] Performing a Fourier transform on the time term of equation (7) and converting it to the frequency domain yields the equation

[0077] in, The wave is a frequency-domain scattered wave, where ω is the angular frequency. Let be the source function in the frequency domain, γ represent the intermediate variable, i is an imaginary number, and β is the viscosity term that considers the energy attenuation of the wave during propagation.

[0078] Perform a Fourier transform on equation (8) in the x-direction (two-dimensional case) or a two-dimensional Fourier transform in the x and y directions (three-dimensional case) to transform it to the wavenumber domain, where γ = 2α + α 2 Let be a function of (x, y, z), which is related to... Fourier transform of multiplication This indicates that the ordinary differential equation for the scattering wave in the spatial wavenumber mixing domain is now obtained in the frequency domain:

[0079] Where γ=2α+α 2 Let be a function of (x, y, z). Indicates γ and Fourier transform of multiplication The wave is a spatial wavenumber-domain scattered wave in the frequency domain, where z is the vertical coordinate. For a spatial wavenumber domain seismic source in the frequency domain, k x k is the wave number corresponding to the x-direction. y Let k be the wave number in the y-direction. In the two-dimensional case, k y =0.

[0080] Step Six: Transform the system of ordinary differential equations into a series of ordinary differential equations

[0081] Because there are many pairs of kx and ky, each pair of kx and ky corresponds to an equation, and all of them together form a system of equations. In equation (9), The Fourier transform according to the convolution theorem is: Spectra of scattered wavefields at different wavenumbers Because they are closely related through convolution, equation (9) is actually a set of ordinary differential equations corresponding to the wavenumber. Solving this set of ordinary differential equations directly yields the scattered wave field spectra corresponding to all wavenumbers. However, this is a holistic problem that is not easy to break down and therefore not easy to compute in parallel. The computation time and memory requirements are similar to those of conventional methods such as finite element and finite difference methods. For large-scale, detailed numerical simulation problems, it is difficult to implement due to the large computation time and memory requirements.

[0082] To address this problem, this application proposes the function product Fourier transform conjecture, which states that for velocity anomaly regions, for any wave numbers kx and ky, there always exists a wave function... This causes the anomalous velocity γ(x,y,z) and the scattered wave field to... The Fourier transform of the product is equal to the wave function. and scattered wave field spectrum The product, ignoring variables, has the following mathematical expression:

[0083] Substituting equation (10) into equation (9), equation (9) can be written as:

[0084] Equation (11) is the ordinary differential equation corresponding to wave numbers kx and ky. Different wave numbers are uncorrelated, thus realizing 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 have a small amount of computation, which can greatly reduce the amount of computation and memory requirements.

[0085] Step 7: Solving Ordinary Differential Equations.

[0086] make Equation (11) can be written as:

[0087] There are usually no anomalous velocities in the upper and lower boundary regions of the solution domain, therefore The general solution to this equation can be written as: Among them, Ae -kz For a downward wave, Be kz It is an upward wave. At the upper boundary Zmin, there is only an upward wave, and at the lower boundary Zmax, there is only a downward wave.

[0088] Therefore, the boundary conditions can be obtained as follows:

[0089] The boundary value problem is:

[0090] For the boundary value problem corresponding to a series of wave numbers kx and ky in equation (15), a five-diagonal linear equation system is obtained by using the one-dimensional finite element method with quadratic interpolation. The wavenumber domain scattered wave field can be obtained by quickly solving the five-diagonal linear equation system using the chasing method. The spatial domain scattered wave field distribution can be obtained by first determining the distribution of the wave field, and then using the inverse Fourier transform. This solution method can fully consider both computational accuracy and efficiency, thus enabling efficient and high-precision numerical simulation of large-scale complex models.

[0091] Step 8: Iterate to find the exact solution of the scattered wave field.

[0092] This embodiment establishes an iterative algorithm, through which the scattered wave field gradually approaches the exact solution through iterative calculation.

[0093] Step 1:

[0094] According to equation (9), the system of one-dimensional ordinary differential equations in the wavenumber domain of any frequency space is as follows:

[0095] To generate initial values, the iteration number i = 0. According to the Born approximation method, let γ(x,y,z) = 0. Equation (16) is transformed into an ordinary differential equation corresponding to the wave numbers kx and ky, and the Born approximation solution can be obtained.

[0096] Step 2:

[0097] Start iterating, i = i + 1;

[0098] Based on the conjecture proposed in this application It can be obtained

[0099] Step 3:

[0100] Based on the conjecture (10) proposed in this invention, and the wave function Substituting both into equation (16), we get Satisfied ordinary differential equations

[0101] Solving equation (17), we get

[0102] Step 4:

[0103] For all wave numbers kx and ky, determine whether the following formula holds true.

[0104] Where ε is a decimal close to 0, if equation (18) holds, If the equation is an exact solution, exit the iteration. The smaller ε 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, then let the corrected solution be the solution with the perturbation added to the current solution. That is, the corrected solution is Proceed to the next step.

[0105] Step 5:

[0106] Will Substituting into (16), we get

[0107] Subtracting (19) from (17) yields the solution for the perturbation. One-dimensional ordinary differential equation system

[0108] Step 6:

[0109] right The system of one-dimensional ordinary differential equations that satisfy the given conditions can be solved using the Born approximation. According to the Born approximation method, by setting γ(x,y,z) to 0 on the left-hand side of the equations, we can obtain the solution. Based on the conjecture (10) proposed in this invention, there exists a Make the following equation true

[0110] In the formula, The wave function representing the perturbation field;

[0111] Seek In equation (20) use By substitution of variables, we obtain a one-dimensional ordinary differential equation:

[0112] Solving for the given information The corrected solution is Return to Step 2 to begin the next iteration, until equation (18) is reached. Established.

[0113] This patent proposes a function product Fourier transform conjecture, which transforms the problem of solving three-dimensional partial differential equations into solving one-dimensional ordinary differential equations corresponding to a series of wave numbers kx and ky. The ordinary differential equations corresponding to different wave numbers are easily solved in parallel with minimal computational cost, thus significantly reducing computational load and memory requirements. By employing the holographic Fourier transform, on the one hand, it ensures that the external boundary conditions of the numerical simulation region are completely consistent with the real physical boundary conditions, thereby obtaining the true spectrum of the scattered wave field; on the other hand, it utilizes all spectral information corresponding to all wave numbers. Therefore, the holographic numerical simulation model of the scattered wave field is consistent with the seismic geophysical model of the Earth's medium, and thus the kinematic and dynamic characteristics of seismic waves can be accurately simulated, hence the name holographic numerical simulation.

[0114] The proposed holographic numerical simulation algorithm for acoustic wave scattering is then tested.

[0115] The test computer was configured with a 12th Gen Intel(R) Core(TM) i9-12900KS, 3.40GHz, and 128GB of RAM.

[0116] 1. Correctness verification

[0117] To verify the correctness of the algorithm proposed in this invention, the commercial software COMSOL Multiphysics was used to perform a numerical simulation of the sound wave, and the finite element numerical solution was used as a reference solution to verify the correctness of the algorithm.

[0118] A. Two-dimensional

[0119] The two-dimensional case was validated. The computational model is shown in Figure 2. The background velocity is 3000 m / s, α is 0.1 (i.e., the anomalous velocity is 2727 m / s), the background sound wave is a plane wave with a frequency of 15 Hz, the horizontal partitioning nodes Nx and the vertical partitioning nodes Nz are both 201, uniform sampling is used in the spatial domain, the number of horizontal wavenumbers Nkx is 401, the wavenumber range is -0.5 to 0.5, and uniform sampling is used in the wavenumber domain. The convergence condition was reached after four iterations, with a computation time of 1.5 s.

[0120] Figure 3 shows the total field calculated by this invention and the total field calculated by the finite element method, along with their errors. The relative root mean square error between the two is less than 1%, verifying the correctness and accuracy of the algorithm proposed in this invention. Figure 4 shows a comparison of the waveforms in the vertical z-direction when the calculation model is x = 0m; the two waveforms almost completely overlap.

[0121] Figure 5 shows the scattered field calculated by this invention and the total field calculated by the finite element method, along with their errors. Figure 6 shows the comparison curves of the scattered fields calculated by this invention and the finite element method at x = 0m in the calculation model.

[0122] The iterative convergence of this algorithm is characterized by relative error, as shown in the following equation:

[0123] The iterative curve is shown in Figure 7, which shows that the algorithm can converge stably.

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

[0125] B. Three-dimensional

[0126] The three-dimensional prism model is shown in Figure 8. The numerical simulation study area is 1000m × 1000m × 500m. The background sound wave is a plane wave with a velocity of 3000m / s and a calculation frequency of 15Hz. The parameter α, representing the anomalous velocity, is 0.1, resulting in an anomalous velocity of 2727m / s. The anomalous body is located at the center of the model and is a cube with a side length of 100m. The spatial domain mesh is divided into Nx, Ny, and Nz with 201 elements, using uniform sampling. The wavenumbers Nkx and Nky are 401, with a wavenumber range of -0.06 to 0.06, also using uniform sampling. The simulation reached convergence after 3 iterations, with a computation time of 1.5s.

[0127] Taking the section at y = 0m, the total field calculated by this invention and the total field calculated by the finite element method, along with their errors, are shown in Figure 9. The relative root mean square error between the two is less than 1%, verifying the correctness and accuracy of the algorithm proposed in this invention in the three-dimensional case. Figure 10 shows a comparison of the waveforms in the vertical z-direction when the calculation model is y = 0m and x = 0m; the two almost completely overlap.

[0128] At the y=0m section, the scattered field calculated by this invention and the total field calculated by the finite element method, along with their errors, are shown in Figure 11. It can be seen that the two are very close, verifying the correctness and accuracy of the algorithm proposed in this invention. Near the boundary, the scattered field obtained by the finite element method still exhibits boundary effects, while the result of this invention shows no boundary effects, consistent with the actual situation. Figure 12 shows a comparison of the waveforms in the vertical z-direction when the calculation model is y=0m and x=0m; the two are also basically identical.

[0129] Figure 13 shows the iterative convergence curve of the model in the three-dimensional case, which shows that the algorithm can also converge stably in the three-dimensional case.

[0130] In summary, the correctness and accuracy of the algorithm proposed in this invention can be verified in both two-dimensional and three-dimensional cases.

[0131] 2. Efficiency

[0132] The efficiency of the algorithm is explained by comparing it with the conventional finite element method.

[0133] Let ΔN be the factor by which the total number of grid nodes changes, and Δt be the factor by which the computation time changes. The relationship between them is Δt = ΔN. a , where the value of 'a' varies in different algorithms.

[0134] The algorithm efficiency statistics for the two-dimensional model are as follows:

[0135] Table 1. Iteration time and memory usage statistics of the proposed algorithm and the COMSOL finite element algorithm in the two-dimensional case.

[0136] The two-dimensional model and the Comsol finite element method were compared under the same mesh partitioning. The algorithm proposed in this invention has an a-to-score of 0.6, while the finite element method has an a-to-score of 1.1. This means that when the number of nodes increases fourfold, the computation time of the algorithm proposed in this patent increases by a factor of 2.3, which is less than fourfold, indicating a sublinear increase. In contrast, the finite element method increases by a factor of 4.2, which is greater than fourfold, indicating an exponential increase. This shows that the computation time of the algorithm proposed in this invention exhibits a sublinear change with the number of nodes, while the computation time of the finite element method increases exponentially with the number of nodes. Therefore, in the two-dimensional case, the more nodes there are, the more significant the algorithm's advantage becomes.

[0137] The efficiency statistics for the 3D model are as follows:

[0138] Table 2. Iteration time and memory usage statistics of the proposed algorithm and the COMSOL finite element algorithm in three dimensions.

[0139] Figure 14 shows a comparison of efficiency and memory usage between the two-dimensional model and the conventional finite element method, while Figure 15 shows a comparison of efficiency and memory usage between the three-dimensional model and the conventional finite element method.

[0140] The 3D model and the Comsol finite element method were compared under the same mesh partitioning. The algorithm proposed in this invention has an a-to-value ratio of 0.8, while the finite element method has an a-to-value ratio of 1.5. This means that when the number of nodes increases by 8 times, the computation time of the algorithm proposed in this patent increases by a factor of 5.27, which is less than 8 times and represents a sublinear increase. In contrast, the finite element method increases by approximately 22 times, which is greater than 8 times and represents an exponential increase. Therefore, the algorithm of this invention has a more significant advantage in 3D cases, with computational efficiency more than three orders of magnitude higher for small-scale models (81×81×41 mesh). This computational efficiency advantage of this invention increases rapidly as the scale of the numerical simulation increases.

[0141] Therefore, the algorithm proposed in this invention has a much lower memory usage and computation time than the traditional finite element method, whether in two-dimensional or three-dimensional cases, and the advantages of this invention become more obvious as the computation scale increases.

[0142] Model Examples

[0143] Two-dimensional:

[0144] The design incorporates an air layer as shown in Figure 16, with an air layer depth of 300m and a subsurface depth of 1000m. The air layer velocity is 340m / s, and the background field consists of plane waves with a background velocity of 3000m / s. Two anomalies are included: a low-velocity anomaly with an α value of 0.1 and a velocity of 2727m / s, and a high-velocity anomaly with an α value of -0.2 and a velocity of 3750m / s. The model has 501 nodes in both the x and z directions, and 501 wavenumbers are selected. Both the spatial and wavenumber domains are uniformly sampled. The computation frequency is 15Hz. The model reached convergence after three iterations.

[0145] Figure 17 below shows the calculation results of the two-dimensional combined model. 17-1 is the plane wave background field, 17-2 is the scattered field, the waveforms are clear, the anomaly differences are obvious, and the greater the difference in velocity from the background, the stronger the anomaly field response. 17-3 is the total field, showing that the wave field below the low-velocity body can pass through, while the wave field below the high-velocity body cannot, indicating that the high-velocity body has a shielding effect on the wave field. From the two-dimensional combined model with air, it can be seen that the present invention can perform calculations with air in a two-dimensional context, and the model is realistic.

[0146] 3D:

[0147] The design incorporates an air layer as shown in Figure 18, with an air layer depth of 200m and a subsurface depth of 500m. The air layer velocity is 340m / s, and the background field consists of plane waves with a velocity of 3000m / s. Two anomalies are included: a low-velocity anomaly with an α of 0.1 and a velocity of 2727m / s, and a high-velocity anomaly with an α of -0.2 and a velocity of 3750m / s. The low-velocity anomaly is located in the negative x-direction. The number of nodes in the x and y directions is 501 each, and the number of nodes in the z-direction is also 501. 501 wavenumbers are selected, and sampling is uniform in both the spatial and wavenumber domains. The computation frequency is 15Hz. The model reaches convergence after three iterations.

[0148] Figure 19 below shows the calculation results of the three-dimensional combined model. 19-1 is the plane wave background field, and 19-2 is the scattered field. The pattern is the same as in the two-dimensional case: the greater the difference between the anomalous 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 below the low-velocity body can pass through, while the wave field below the high-velocity body cannot, indicating that the high-velocity body has a shielding effect on the wave field. From the three-dimensional combined model with air, it can be seen that this invention can perform calculations with air in a three-dimensional context. The model is realistic and can adapt to large-scale, efficient, and high-precision numerical calculations.

[0149] Figure 20 below shows the wavefield profile at y=0m and its logarithmic representation. It demonstrates that the wavefield is continuous at the interface between air and the subsurface medium, simulating a realistic situation considering air. The logarithmic diagram also reveals anomalous wavefields within the air; without considering air, the wavefield situation in the air would be impossible to determine. Therefore, this further illustrates the necessity of this invention starting from a realistic model and considering the air layer.

[0150] Figure 21 below shows the wave field planar diagrams of the ground and the center of the anomalous body, which can still clearly reflect the occurrence of the anomalous body. The wave field amplitude is larger at the anomalous body and becomes weaker at the ground, which is in line with the physical law. This once again illustrates the correctness and rationality of the algorithm proposed in this invention.

[0151] As can be seen, this invention develops 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 is stable and convergent, and the solution efficiency of small-scale numerical simulations is 3 to 5 orders of magnitude higher than conventional methods. The memory usage is reduced by 2 to 4 orders of magnitude compared to conventional methods, and the computational advantages become more pronounced as the computational scale increases. This invention develops an efficient and high-precision numerical simulation method for acoustic wave scattering in the real Earth medium, applicable to large-scale arbitrarily complex conditions, namely, a holographic numerical simulation method for acoustic wave scattering.

[0152] This application also provides a holographic numerical simulation system for the acoustic wave scattering wavefield of the Earth's medium, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the steps of the above-described method. This holographic numerical simulation system for the acoustic wave scattering wavefield of the Earth's medium can implement various embodiments of the holographic numerical simulation method for the acoustic wave scattering wavefield of the Earth's medium and achieve the same beneficial effects, which will not be elaborated here.

[0153] This application also provides a computer storage medium, including a memory, a processor, and a computer program stored on the memory and executable on the processor. When the processor executes the computer program, it implements the steps of the above-described method. This computer storage medium can implement various embodiments of the holographic numerical simulation method for acoustic wave scattering wavefields in the Earth's medium and achieve the same beneficial effects, which will not be elaborated here.

[0154] The preferred embodiments of the present invention have been described in detail above. It should be understood that those skilled in the art can make numerous modifications and variations based on the concept of the present invention without creative effort. Therefore, all technical solutions that can be obtained by those skilled in the art based on the concept of the present invention through logical analysis, reasoning, or limited experimentation on the basis of existing technology should be within the scope of protection defined by the claims.

Claims

1. A holographic numerical simulation method for acoustic wave scattering wavefield in the Earth's medium, characterized in that, include: S1: Define the time-domain acoustic wave equations specific to the characteristics of the Earth's medium; S2: Based on velocity field separation, background field and scattered field separation, the time-domain acoustic wave equation is divided into background wave field equation and scattered wave field equation; and based on the scattered wave field equation, the time-domain scattered wave field wave equation is obtained. S3: Determine the set of spatial wavenumber mixing domain scattering wave equations in the frequency domain based on the time-domain scattering wave field wave equations; S4: The spatial domain scattered wave field distribution is calculated based on the spatial wavenumber mixing domain scattering wave equations in the frequency domain as the result of the holographic numerical simulation of the acoustic wave scattered wave field in the Earth medium.

2. The holographic numerical simulation method for acoustic wave scattering wavefield in the Earth's medium according to claim 1, characterized in that, The time-domain acoustic wave equation in S1 satisfies the following relationship: Where ψ represents the sound wave field, c represents the sound wave propagation speed in m / s, and f s Let be the source function, t be time (in seconds), and β be the viscous term considering energy attenuation during wave propagation. Represents the Laplace operator. This indicates the partial derivative.

3. The holographic numerical simulation method for acoustic wave scattering wavefield in the Earth's medium as described in claim 2, characterized in that, S2 include: S21: The total velocity of sound waves is separated into a background velocity field and an anomalous velocity field, satisfying the following relationship: Where c is the total velocity of sound wave propagation, c0 represents the background velocity field, and α is the spatial coordinate function representing the anomalous velocity field; S22: The acoustic wave field after velocity field separation is separated into a background field and a scattered field, satisfying the following relationship: ψ=ψ0+ψ a (3); Where ψ is the sound wave field, ψ a Let ψ0 be the scattering field and ψ0 be the background field. S23: Based on equations (2) and (3), equation (1) can be split into a background wave field equation and a scattered wave field equation. The background wave field equation is as follows: In the formula, β is a viscous term that considers the energy attenuation of the wave during propagation. Represents the Laplace operator. This indicates finding the partial derivative, f s Here, t is the source function; The equation for the scattered wave field is as follows: S24: Order Let γ = 2α + α 2 γ is an intermediate variable, and f0 represents an intermediate variable of the background field. Equation (5) is transformed into the time-domain scattering wave field wave equation:

4. The holographic numerical simulation method for acoustic wave scattering wavefield in the Earth's medium as described in claim 1, characterized in that, S3 includes: S31: Take a Fourier transform of the time term t in the time-domain scattered wave field wave equation to obtain the equation in the frequency domain, as follows: in, The wave is a frequency-domain scattered wave, where ω is the angular frequency. Here, γ represents the frequency domain source function, i is an imaginary number, and β is a viscosity term that considers the energy attenuation during wave propagation. S32: Performing a two-dimensional Fourier transform of equation (7) in the x and y directions to the wavenumber domain, we obtain the following equation for the spatial wavenumber mixed domain scattering wave: Where γ=2α+α 2 Let be a function of (x, y, z). Indicates γ and Fourier transform of multiplication The wave is a spatial wavenumber-domain scattered wave in the frequency domain, where z is the vertical coordinate. For a spatial wavenumber domain seismic source in the frequency domain, k x k is the wave number corresponding to the x-direction. y Let k be the wave number in the y-direction, and k is the wave number in the two-dimensional case. y =0.

5. The holographic numerical simulation method for acoustic wave scattering wavefield in the Earth's medium as described in claim 4, characterized in that, S4 includes: S41: Combine multiple sets of spatial wavenumber mixing domain scattering wave equations into a set of spatial wavenumber mixing domain scattering wave equations, and transform them into a series of ordinary differential equations; S42: The spatial domain scattered wave field distribution is obtained by solving the series of ordinary differential equations based on a set iterative method. The steps of the set iterative method include a set function product Fourier transform conjecture.

6. The holographic numerical simulation method for acoustic wave scattering wavefield in the Earth's medium according to claim 5, characterized in that, The series of ordinary differential equations in S41 satisfy the following relationship:

7. The holographic numerical simulation method for acoustic wave scattering wavefield in the Earth's medium according to claim 6, characterized in that, S42 includes: S421: According to equation (9), the system of one-dimensional ordinary differential equations in the wavenumber domain of any frequency space is as follows: To generate initial values ​​for equation (10), let i = 0. According to the Born approximation method, let γ(x,y,z) = 0, transform equation (9) into an ordinary differential equation corresponding to the wave numbers kx and ky, and solve for the Born approximation solution. S422: Let i = i + 1, start iterative calculation, and conjecture based on the set function product Fourier transform. Find the wave function S423: Based on the established conjecture of the Fourier transform of the product of functions, the wave function... Substituting both into equation (10), we get The satisfied ordinary differential equation is as follows: Solving equation (11) yields the solution for this problem. S524: For all wave numbers kx and ky, obtain And determine whether the following expression is true: Where ε is a specified decimal number close to 0; If equation (12) holds, then the solution obtained in this case... The iterative calculation ends when the exact solution of the equation is obtained, and the results of the holographic numerical simulation of the acoustic wave scattering wave field in the Earth's medium are obtained. If equation (12) does not hold, then let the corrected solution be the current solution plus the perturbation solution. Even if the corrected solution is Continue to the next calculation; S425: Will Substituting into equation (10), we get: Subtracting (13) from (11) yields the solution for the perturbation. A system of one-dimensional ordinary differential equations: S426: Solving equation (14) using the Born approximation yields... Based on the established conjecture of the Fourier transform of the product of functions, there exists a... (k x ,k y ,z), such that the following equation holds: In the formula, The wave function representing the perturbation field; Seeking In equation (15) use By substitution of variables, we obtain a one-dimensional ordinary differential equation: Obtain the perturbation solution Let the modified solution be Return to S522 for the next iteration calculation.

8. A holographic numerical simulation system for the acoustic wave scattering wavefield of the Earth's medium, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes a computer program, it implements the steps of any one of claims 1 to 7.

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 the processor executes a computer program, it implements the steps of any one of claims 1 to 7.

Images