Seismic wave detection method and system, electronic equipment and storage medium

Abstract

The embodiment of the invention provides a seismic wave detection method and system, electronic equipment and a storage medium, and belongs to the field of seismic wave engineering detection. The method comprises the steps that seismic sources and detectors are arranged in advance according to actual working conditions, seismic waves are excited at the positions of the seismic sources, and seismic wave data are collected through the detectors; based on seismic wave data, solving a time fractional order viscoelastic acoustic wave equation by using a short memory fast algorithm so as to simulate the propagation behavior of seismic waves in a real stratum medium; and performing Bayesian full waveform inversion and reverse time migration imaging based on the seismic wave data to obtain a stratum structure and stratum characteristics. The memory effect of a time fractional order viscoacoustic wave equation is reasonably simplified through a short memory fast algorithm, so that the storage requirement for wave field historical data is reduced, efficient numerical solution is achieved, on the premise that precision is guaranteed, memory occupation is greatly reduced, and the cost is reduced. And the time fractional order viscoelastic acoustic wave equation has higher operability and calculation efficiency in engineering application.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic wave engineering detection, and in particular to a seismic wave detection method, system, electronic equipment and storage medium. Background Art

[0002] Traditional seismic wave exploration methods are based on the wave equation in the classical perfect elastic theory, which assumes that seismic waves do not attenuate and maintain phase when propagating in the medium. However, many scholars have found through experimental evidence that seismic waves will experience significant waveform attenuation and phase delay when propagating in actual formation media. This phenomenon is caused by the viscoelastic properties of the formation medium, which cannot be effectively described by the classical wave equation, making it difficult to accurately simulate the propagation behavior of seismic waves in complex formations.

[0003] To make up for this deficiency, the time fractional viscoacoustic wave equation was introduced into the study of seismic wave propagation. By introducing a time fractional derivative that is much smaller than 1, this equation can accurately describe the viscoelastic stress-strain relationship in coal mine media, making it have significant advantages in simulating the waveform attenuation and phase delay of seismic wave propagation. However, this improvement also brings new challenges - due to the existence of time fractional derivatives, in order to obtain an accurate numerical solution, it is necessary to record the complete historical information of the wave field during the calculation process, which leads to huge memory requirements and limits its application in engineering practice. Summary of the invention

[0004] The purpose of the embodiments of the present invention is to provide a seismic wave detection method, system, electronic device and storage medium, which are used to fully or at least partially solve the technical problem that in order to obtain an accurate numerical solution, it is necessary to record the complete historical information of the wave field during the calculation process, resulting in huge memory requirements and limiting its application in engineering practice.

[0005] In order to achieve the above object, an embodiment of the present invention provides a method for detecting seismic waves, comprising: Arrange seismic sources and geophones in advance according to actual working conditions, stimulate seismic waves at each seismic source location, and collect seismic wave data through geophones; Based on the seismic wave data, a short memory fast algorithm is used to solve the time fractional viscoelastic acoustic wave equation to simulate the propagation behavior of seismic waves in real stratum media; Bayesian full waveform inversion and reverse time migration imaging are performed based on the seismic wave data to obtain the formation structure and formation characteristics.

[0006] Alternatively, the time fractional viscoelastic acoustic wave equation is as follows:

[0007] In the formula, , Respectively represent the particles in , The velocity component in the axial direction, is the reference frequency, is the fractional order, represents the phase velocity, is the density of the medium, is the source term acting on the unit volume, , The order is The fractional derivative of Represents stress.

[0008] Optionally, solving the time fractional viscoelastic acoustic wave equation using a short memory fast algorithm includes: introducing an auxiliary function to convert the solution of the time fractional derivative of the time fractional viscoelastic acoustic wave equation into the calculation of a weak singular integral over a semi-infinite domain.

[0009] Optionally, the auxiliary function formula is as follows:

[0010] In the formula, Represents an auxiliary function, is a parameter introduced in the fractional derivative, which is related to the characteristics of the relaxation process and the memory effect of the material. Indicates that the auxiliary function is in t = the initial value when 0. is the strain function The first derivative with respect to time represents the rate of change of strain, is the strain function The second derivative with respect to time represents the acceleration of the strain, It represents the historical integral term, reflecting the memory effect of the strain function, that is, the current strain not only depends on the current state, but also is related to the past evolution process.

[0011] The optional, short memory fast algorithm formula is as follows:

[0012] In the formula, is the stress function, indicating the position x and time t The stress distribution under . represents the phase velocity, is the number of storage layers of the auxiliary function, represents the scaling factor, is the weight factor, which is related to the fractional derivative and represents the weight coefficient of a certain layer in the calculation. is the auxiliary variable of the fractional derivative, indicating thej The node positions or values ​​of the layer, is the auxiliary function, and the strain function Correlation describes the memory effect in the fractional-order model.

[0013] Optionally, Bayesian full waveform inversion and reverse time migration imaging are performed based on the seismic wave data to obtain the formation structure and formation characteristics, including: Obtain observation data through detectors arranged in advance according to actual working conditions; An initial geological model is constructed based on geological knowledge, historical data and prior information, and numerical simulation is performed using the initial geological model according to a short-memory fast algorithm to calculate theoretical seismic waveforms; Computers are used to simulate seismic wave propagation forward and short-memory fast algorithms are used to numerically solve the time fractional viscoelastic acoustic wave equation to predict the propagation behavior of seismic waves in complex underground structures. Calculate the error between the observed data and the theoretical data obtained by forward modeling, and determine whether the residual between the observed data and the forward modeling data meets the preset threshold requirements. If so, obtain the target geological model; if not, use the Bayesian statistical method to sample the geological model parameters; Based on the parameter distribution obtained by sampling, parameter estimation and uncertainty analysis are performed to obtain the probability distribution of geological structures.

[0014] Optionally, during the simulation of seismic wave propagation, the region is gridded in two-dimensional space, wherein the grid spacing in the horizontal direction and the vertical direction is set equal; On this basis of gridding, the seismic wave field is analyzed in the frequency domain through discrete Fourier transform to realize the calculation of spatial derivatives in the frequency domain, and the fast Fourier transform algorithm is used to improve the calculation efficiency.

[0015] On the other hand, the present invention also provides a seismic wave detection system, comprising: The acquisition unit is used to arrange the seismic sources and detectors in advance according to the actual working conditions, and to excite seismic waves at each seismic source position and collect seismic wave data through the detectors; A simulation unit, used for solving the time fractional viscoelastic acoustic wave equation based on the seismic wave data by using a short memory fast algorithm to simulate the propagation behavior of seismic waves in real stratum media; The inversion unit is used to perform Bayesian full waveform inversion and reverse time migration imaging based on the seismic wave data to obtain the formation structure and formation characteristics.

[0016] On the other hand, the present invention also provides an electronic device, 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-mentioned method for detecting seismic waves when executing the program.

[0017] On the other hand, the present invention also provides a non-transitory computer-readable storage medium having a computer program stored thereon, and when the computer program is executed by a processor, the steps of the above-mentioned method for detecting seismic waves are implemented.

[0018] Through the above technical solution, by using the short memory fast algorithm and reasonably simplifying the memory effect of the time fractional viscoacoustic wave equation, the storage demand for wave field historical data is reduced, and efficient numerical solution is achieved. Moreover, while ensuring accuracy, the memory usage is greatly reduced, making the time fractional viscoelastic acoustic wave equation more operable and computationally efficient in engineering applications.

[0019] Other features and advantages of the embodiments of the present invention will be described in detail in the subsequent detailed description. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] The accompanying drawings are used to provide a further understanding of the embodiments of the present invention and constitute a part of the specification. Together with the following specific implementations, they are used to explain the embodiments of the present invention, but do not constitute a limitation on the embodiments of the present invention. In the accompanying drawings: Figure 1 is an implementation flow chart of a seismic wave detection method provided by an embodiment of the present invention; Figure 2 is a flow chart for obtaining formation structure and formation characteristics provided by an embodiment of the present invention; Figure 3 Schematic diagram of a pseudo-spectral solution process provided by an embodiment of the present invention; Figure 4 is a schematic flow chart of a Bayesian full waveform inversion provided by an embodiment of the present invention; Figure 5 It is a structural schematic diagram of a seismic wave detection system provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0021] The specific implementation of the embodiment of the present invention is described in detail below in conjunction with the accompanying drawings. It should be understood that the specific implementation described here is only used to illustrate and explain the embodiment of the present invention, and is not used to limit the embodiment of the present invention.

[0022] See also Figure 1 As shown, it is a flowchart of an implementation method of a seismic wave detection method provided by an embodiment of the present invention, comprising the following execution steps: Step 100: Arrange seismic sources and geophones in advance according to actual working conditions, excite seismic waves at each seismic source location, and collect seismic wave data through geophones.

[0023] Specifically, the seismic source and the detector are designed and arranged according to the actual working conditions, seismic waves are excited at each source position, the seismic wave information is automatically collected and stored by the detector, and the collected data is subjected to processing such as denoising and file format conversion.

[0024] Exemplarily, the information of the seismic wave includes at least the amplitude, frequency and arrival time of the seismic wave.

[0025] Step 101: Based on the seismic wave data, a short memory fast algorithm is used to solve the time fractional viscoelastic acoustic wave equation to simulate the propagation behavior of seismic waves in real stratum media.

[0026] By adopting a short memory fast algorithm and reasonably simplifying the memory effect of the time fractional viscoacoustic wave equation, the storage requirements for wave field history data are reduced, and efficient numerical solution is achieved. This method significantly reduces memory usage while ensuring accuracy, making the time fractional viscoacoustic wave equation more operable and computationally efficient in engineering applications.

[0027] Step 102: Perform Bayesian full waveform inversion and reverse time migration imaging based on the seismic wave data to obtain the formation structure and formation characteristics.

[0028] Specifically, first, the medium parameter model of the target area is constructed based on the Bayesian full waveform inversion method through the collected seismic data. The Bayesian method can introduce prior information in the inversion process, so as to more accurately handle uncertainty factors and provide confidence intervals for the inversion results. In full waveform inversion, the stored information such as seismic wave amplitude, frequency, phase and arrival time is used to gradually optimize the formation model through the difference between forward simulation and observation data to match it with the actual seismic response, so as to obtain more accurate formation characteristics. Next, based on reverse time migration imaging technology, the information of reflected waves and scattered waves is used to further enhance the resolution of the model. The reverse time migration imaging process focuses the reflection event on its true underground position by back-propagating the seismic wave field, realizing fine structural imaging. Finally, through the combination of Bayesian full waveform inversion and reverse time migration imaging, not only can the high-precision distribution of formation medium parameters be obtained, but also the formation structure and characteristics can be restored more comprehensively and accurately through reverse time imaging.

[0029] In some embodiments, the process of designing a viscoelastic wave short memory fast forward algorithm to solve the viscoelastic seismic wave equation in step 101 includes: modeling in a numerical model according to the location of the area to be detected and the seismic excitation source and the detector, determining the relevant parameters in the viscoelastic wave equation, so as to complete the construction of the viscoelastic seismic wave equation group. The numerical model in this embodiment is based on the actual engineering application scenario, and by introducing the viscoelastic theory, it can more accurately simulate the absorption and velocity dispersion phenomena of seismic wave propagation in the actual stratum. The constructed viscoelastic wave equation group will be solved by a short memory fast algorithm. This embodiment aims to simplify the solution process of the time fractional viscoelastic wave equation into a solution task within a local time, thereby reducing the demand for historical wave field storage. Through the short memory algorithm (Efficient Short Memory Algorithm, ESMA), the viscoelastic wave equation is quickly solved, which is particularly suitable for engineering application scenarios with high computing loads.

[0030] Furthermore, the viscoelastic wave equation used to simulate the propagation of seismic waves in coal mine media is controlled by three sets of equations: The first group is the linear momentum conservation equation, which describes the motion state of mass elements in the medium under the action of force. Specifically, when seismic waves propagate in the medium, each medium unit is affected by stress and inertial force, resulting in displacement and velocity changes. The momentum conservation equation ensures that during the wave propagation process, the momentum change of each medium unit can match the external force it is subjected to, thereby satisfying the basic physical law of momentum conservation. The specific form of this equation is:

[0031] In the formula, is the component of the stress tensor, which represents the direction j Acting on the direction i The stress, and are the particles in , The displacement component in the axial direction, is the density of the medium, is the source term acting on the unit volume, yes About time The second-order partial derivative of represents the component of the particle acceleration, It is stress The spatial gradient of represents the change of stress with spatial position.

[0032] The second set of equations is used to describe the relationship between displacement and strain, which is called the geometric equation. Under the action of external force, the relative displacement between each particle will cause the shape of the object to change. The strain can be described by the uneven displacement between each particle. The geometric equation is:

[0033] In the formula, are the components of the strain tensor, representing the position and time Down, in the direction and The relative deformation degree, represents the displacement component of the particle in the direction of the coordinate axis, express right The first-order partial derivative of is the component describing the displacement gradient, express right The first-order partial derivative of , also represents another component of the displacement gradient.

[0034] The third group is the constitutive equations, which are used to describe the propagation characteristics of seismic waves in viscoelastic media (such as coal mine media). When seismic waves propagate in different media, their stress-strain relationship will be affected by the characteristics of the medium. Unlike elastic media, viscoelastic media will produce viscosity when subjected to force, which manifests as waveform attenuation and phase delay. This characteristic makes the viscoelastic theory more accurate in describing the propagation of seismic waves in actual strata, especially in coal mine media with viscoelastic characteristics, where the attenuation and hysteresis effects of the waves are more significant.

[0035] The stress-strain relationship in the constitutive equation provides a more realistic model for the seismic wave propagation process, so that the seismic wave propagation algorithm proposed in this application can more accurately reflect the true response of the stratum in engineering applications. This stress-strain relationship is not only suitable for simulating the attenuation characteristics of seismic waves, but also provides a basis for subsequent Bayesian full waveform inversion and reverse time migration imaging, so that waveform changes caused by the viscoelasticity of the medium can be better corrected during imaging, thereby improving the accuracy and reliability of stratum imaging. The specific expression of the constitutive equation is as follows:

[0036] in, The components of the stress tensor, is the convolution symbol, The order is The fractional derivatives of describe the nonlocal time dependence, is the order of the fractional order, is the bulk modulus, represents the components of the strain tensor, represents the time scaling factor, where is the reference time (usually the same as the reference frequency Related, satisfying = ),index It reflects the time dependence of the fractional-order model. For Kjartansson's The relaxation function in the model is expressed as:

[0037] In the formula, is the bulk modulus, is a step function, is the gamma function, is the reference time, is the reference frequency, is the fractional order, from the equation and as well as The time fractional viscoelastic acoustic wave equation can be derived as:

[0038] In the formula, , Respectively represent the particles in , The velocity component in the axial direction, is the reference frequency, is the fractional order, represents the phase velocity, is the density of the medium, is the source term acting on the unit volume, , The order is The fractional derivative of Represents stress.

[0039] In some embodiments, performing step 101 to solve the time fractional viscoelastic-acoustic wave equation using a short memory fast algorithm includes: introducing an auxiliary function to convert the solution of the time fractional derivative of the time fractional viscoelastic-acoustic wave equation into the calculation of a weak singular integral on a semi-infinite domain.

[0040] Specifically, the difficulty in solving the viscoelastic acoustic wave equation is that the equation is nonlocal and the singular kernel in its time fractional derivative is The change is very gradual. To obtain a numerical solution of corresponding accuracy, a large amount of historical data needs to be stored. When studying three-dimensional problems, Zhu used 500 layers of historical data to conduct numerical experiments on viscoelastic acoustic waves. For a three-dimensional case, The computer memory required is huge. The short memory efficient algorithm proposed in this application introduces an auxiliary function to transform the solution of time fractional derivative into the calculation of weak singular integral on semi-infinite domain, which effectively reduces the memory requirement. The following describes the process, first giving the definitions of Caputo fractional derivative and Gamma function: Definition 1: Let is a positive real number, , is a positive integer, the function Defined in the interval or On, called:

[0041] For function left Caputo fractional derivative, where is the symbol of fractional derivative, is the Gamma function, which represents the normalization factor of the fractional derivative. n represents the integer-order part of the fractional derivative, is the kernel function of the fractional derivative, which reflects the non-local memory effect. For function of n The derivative represents the higher-order derivative form when calculating fractional derivatives.

[0042] If the function Defined in the interval or On, called:

[0043] For function Right Caputo fractional derivative, where is the symbol of fractional derivative, is the Gamma function, which represents the normalization factor of the fractional derivative. n represents the integer-order part of the fractional derivative, is the kernel function of the fractional derivative, which reflects the non-local memory effect. For function of n The derivative represents the higher-order derivative form when calculating fractional derivatives.

[0044] Definition 2: The definition of the Gamma function is as follows:

[0045] In the formula t is the integration variable, Represents the region in the complex plane where real numbers are greater than 0.

[0046] From Definition 1, we can get The expression of the Caputo fractional derivative is:

[0047] In the formula is the strain function, which represents the physical quantity that changes with time and space. is the integral variable, representing the past moments in the time domain, is the kernel function, describing the memory effect of fractional derivatives, represents the gamma function.

[0048] From Definition 2, we can get The Gamma function of order is:

[0049] In the formula, is the Gamma function, which is a generalized function of factorial. α is the input parameter of the Gamma function, indicating the order, is the defining integral of the Gamma function, g is the integration variable, Indicates that the Gamma function is and The relationship between is an important property of the Gamma function.

[0050] Will Substitution The formula can be obtained:

[0051] in, is the fractional derivative, representing the strain The non-local time derivative reflects the memory effect of the material, g is the integral variable, and the second equation is obtained by the variable Replaced by, is the new integral variable defined after variable substitution, Current time t With past time The difference between , represents the time scale of the memory effect, describes the weighting of fractional derivatives with respect to time, The exponential decay term represents the rapid decay of the memory effect over time, It represents the derivative of the strain function with respect to time, which indicates the rate of change of strain with time.

[0052] In order to reduce the amount of data storage, the singular kernel in the fractional stress-strain relationship is localized, and the solution of the time fractional derivative is transformed into the calculation of the weak singular integral on the semi-infinite domain. Here, an auxiliary function containing the singular kernel is defined: :

[0053] In the formula, is an auxiliary function used to describe the non-local time dependence in fractional derivatives and to simplify the calculation of memory effects. is a control parameter related to variable substitution and used to adjust the weight of the calculation. The integral variable represents a point in time in the past, is the strain function, describing the position x and time The deformation, is the time derivative of the strain function, which indicates the rate of change of strain. For current time t With past time The time difference between , represents the time scale of the memory effect, is an exponential decay kernel, which controls the strength of the memory effect. As increases, this term decays rapidly, reflecting the local nature of temporal memory.

[0054] Will Substitution In the formula, the Caputo fractional derivative can be transformed into a weak singular integral form on a semi-infinite plane. This method effectively solves In the formula The singular point problem exists when , and it also gets rid of the fractional order The weak singular integral form is:

[0055] In the formula, represents the fractional derivative, represents the strain Non-local dependencies in time. is the fractional order, Normalization factor, is an auxiliary variable, controlling the integral weight, is a weight function, reflecting the weighted contribution of fractional derivatives to different time scales, is an auxiliary function used to represent the historical effect of strain.

[0056] right Derivative of both sides of the equation with respect to time t:

[0057] In the formula, is an auxiliary function used to represent the historical effect of strain. is the derivative of the auxiliary function with respect to time, indicating that at the current time t The rate of change, is the derivative of the strain function with respect to time, indicating the rate of change of strain at the current moment, is an auxiliary variable used to adjust the exponential decay rate of the auxiliary function. Represents the effect of the auxiliary function decaying over time. This formula is used to recursively derive the value of the auxiliary function in numerical calculations.

[0058] Then, by equation The expression of the auxiliary function can be obtained as:

[0059] In the formula, is an auxiliary function used to represent the historical effect of strain, including the influence of time decay and past state. is the initial condition term, indicating that t = 0, the initial value of the auxiliary function, is the cumulative memory effect term, indicating that at each time in the past For current time t contribution. is the parameter for controlling the exponential decay rate, is the derivative of the strain function with respect to time, reflecting the past moments The rate of change of strain.

[0060] For large enough , through integration by parts we can get:

[0061] In the formula, is an auxiliary function used to represent the historical effect of strain. is the initial condition term, indicating that t = 0 is the initial value of the auxiliary function. is the decay factor, which controls the time scale of the memory effect. is the strain function at the current time t The time derivative of represents the rate of change of strain at the current moment, is the strain change rate at the initial time t=0, It is a historical accumulation term, which represents the influence of the acceleration (second-order time derivative of strain) in the past time on the current time.

[0062] In the process of applying short memory efficient algorithms, The strain The value update method of is different from the traditional speed update, but follows the auxiliary function Update, the historical data generated during the update process is stored in the auxiliary function, which can achieve the purpose of reducing the amount of data storage.

[0063] According to the fractional stress-strain relationship, let the fractional First, by introducing a smooth truncation function , decompose it into a front asymptotic range and a back asymptotic range, use the front asymptotic range to capture the main contribution, truncate the infinitesimal contribution of the back asymptotic range, and further reduce the history storage length of the auxiliary function to reduce the amount of data storage. The expression of this process is: (17) In the formula, represents the weight, represents the scaling factor, , Equation The first term, the preasymptotic bound, can be computed using the generalized Laguerre-Gauss quadrature:

[0064] In the formula, The number of storage levels for auxiliary functions. represents a weighted integral that describes the nonlocal nature of fractional derivatives. represents the discrete weight, the weight factor associated with the fractional derivative, is the weight correction factor, including the scaling parameter and discrete node values , adjust the contribution of each layer, is the auxiliary function at discrete points The value of is the integral weight kernel function, which controls the contribution of each node to the total.

[0065] And the following relationship holds:

[0066] In the formula, The modified Gaussian integration node is represented by the scaling factor The adjusted integration point position, is the modified Gaussian integral weight, which is expressed by the scaling factor The adjusted weight value, is the Laguerre function The zero value of , then the expression of Gaussian integration point weight is:

[0067] In the formula is the weight value of the Gaussian integral, which is the zero of the Laguerre polynomial Related, Represents the Gamma function, generalized factorial, used to normalize the weight factor.

[0068] is the generalized Laguerre polynomial No. j zero points, indicating the locations of the integration nodes.

[0069] represents the generalized Laguerre polynomial at zero The value at Indicates the index range of the integration nodes, there are M + 1 node.

[0070] Equation The second term of the post-asymptotic range decays to:

[0071] In the formula, Indicates weight The integral of is a truncation function used to truncate the contribution of a smaller value range and to limit the calculation range to improve efficiency. Auxiliary functions to describe memory effects in fractional derivatives, Represents the upper bound term by truncating the parameter and scaling factor Determine the limit value, is an integer so that the second term has a sufficient decay rate and can be ignored.

[0072] Comprehensive equation , , , we get the stress update formula in the short memory efficient algorithm:

[0073] In the formula, is the stress function, indicating the position x and time t The stress distribution under represents the phase velocity, is the number of storage layers of the auxiliary function, represents the scaling factor, is the Gaussian integral weight, and is related to the node position Correspondingly, for weighted summation, represents the correction factor, combined with the scaling parameter and node value , adjust the weight contribution of each node. is an auxiliary function that describes the time memory effect and indicates that The nonlocal cumulative contribution on .

[0074] At this time, the calculation of the fractional stress-strain relationship is converted to using the Laguerre-Gauss scaling technique to solve a weak singular integral on a semi-infinite domain, and the storage length of the history variable is converted from recording the entire history wave field to storing Layer auxiliary functions effectively reduce the amount of data storage. Through the short memory fast algorithm, the problem of solving the viscoelastic wave equation is transformed into an equivalent problem that is local in time and non-local in space. The method of storing the entire historical wave field in the traditional numerical algorithm is transformed into using several Gaussian integral points to solve the integral problem on the semi-infinite plane. Only a limited number of auxiliary functions need to be stored to accurately solve the viscoelastic wave equation, which greatly improves the calculation efficiency.

[0075] In some embodiments, see Figure 2 As shown, when executing step 102, the following steps may be specifically performed: S200: Acquire observation data through detectors arranged in advance according to actual working conditions.

[0076] S201: construct an initial geological model based on geological knowledge, historical data and prior information, and use the initial geological model to perform numerical simulation and calculate theoretical seismic waveforms according to a short memory fast algorithm.

[0077] S202: Use computers to perform forward simulation of seismic wave propagation and adopt short-memory fast algorithm to numerically solve the time fractional viscoelastic acoustic wave equation to predict the propagation behavior of seismic waves in complex underground structures.

[0078] Specifically, in the process of seismic wave propagation simulation, the area is gridded in two-dimensional space, wherein the grid spacing in the horizontal and vertical directions is set equal; on the basis of this gridding, the seismic wave field is analyzed in the frequency domain by discrete Fourier transform to realize the calculation of spatial derivatives in the frequency domain, and the fast Fourier transform algorithm is used to improve the calculation efficiency.

[0079] In some embodiments, the wave equation is numerically solved by computer simulation using the short memory fast algorithm proposed in step 2 to predict the propagation behavior of seismic waves in underground complex structures. This process takes into account the heterogeneity of geological structures from a physical perspective, including the elastic and viscoelastic properties of different rocks, and the influence of layered structures on wave velocity, attenuation and phase, thereby more accurately reproducing the propagation process of seismic waves in actual underground media.

[0080] In this step, for the calculation of spatial derivatives, the present invention uses a pseudo-spectral method to calculate spatial derivatives. Pseudo-spectral method is a widely used numerical simulation method. Pseudo-spectral method can still maintain high accuracy when using large-pitch grids. The finite difference method requires the use of grids with smaller spacing to maintain accuracy, which will increase the amount of calculation. Moreover, compared with the finite difference method, the pseudo-spectral method can effectively suppress numerical dispersion and improve calculation accuracy. Therefore, this method is widely used in forward simulation of seismic exploration. The pseudo-spectral method uses the derivative theorem of Fourier transform to calculate spatial derivatives, and uses the difference method to calculate time derivatives.

[0081] In the simulation of seismic wave propagation, in order to facilitate calculation, we grid the area in two-dimensional space and set the grid spacing in the horizontal and vertical directions equal. Based on this gridding, the seismic wave field is analyzed in the frequency domain through discrete Fourier transform (DFT) to realize the calculation of spatial derivatives in the frequency domain. The fast Fourier transform (FFT) algorithm can improve the calculation efficiency.

[0082] In this embodiment, the pseudo-spectral method solution process is as follows: Figure 3 As shown, we first use the fast Fourier transform to transform the equation The spatial derivative of the stress in the image is transformed by Fourier transform to obtain:

[0083]

[0084] In the formula, is the Fourier transform operator, which is used to transform the original function from the spatial domain to the spectral domain. and are Fourier frequency variables, corresponding to x and z The wave number, is the Fourier transform of stress.

[0085] Equivalence and Perform inverse Fourier transform on both sides and get:

[0086]

[0087] In the formula, Represents the inverse Fourier transform, which is used to transform the function from the spectral domain to the spatial domain. and are Fourier frequency variables, corresponding to x and z The wave number, is the Fourier transform of stress. Using the above process, the spatial derivative of velocity can also be calculated.

[0088] From the above steps, it can be seen that in order to calculate the derivatives of stress and velocity in the horizontal or vertical direction, the stress value is first converted to the frequency domain using FFT. In the frequency domain, the data of each frequency point is multiplied by the spatial angular frequency and imaginary unit of the point to complete the calculation of the derivative. Subsequently, the derivative value in the corresponding direction can be obtained by returning to the time-space domain through the inverse Fourier transform. Using this feature of the Fourier transform, the complex spatial derivative calculation can be converted into a simple algebraic operation in the frequency domain, making the solution process of the seismic wave equation more concise and efficient.

[0089] S203: Calculate the error between the observed data and the theoretical data obtained by forward modeling, and determine whether the residual between the observed data and the forward modeling data meets the preset threshold requirements. If so, obtain the target geological model; if not, use the Bayesian statistical method to sample the geological model parameters.

[0090] It should be noted that the preset threshold can be set according to the specific application scenario and is not limited here.

[0091] In this embodiment, in order to update the model parameters more accurately, the Bayesian statistical method is used. Specifically, the Markov chain Monte Carlo (MCMC) method is used to sample the model parameters. Bayesian statistics obtains the posterior distribution of the parameters by combining the likelihood information of the prior distribution and the observed data, thereby effectively reflecting the uncertainty of the parameters. In the specific calculation process, the initial parameters of the seismic wave propagation model and its prior distribution are first set. Through the MCMC method, a possible set of parameter values ​​is sampled from the posterior distribution. The parameter values ​​generated by each sampling are used to calculate the propagation of seismic waves in the model and generate corresponding prediction data. Then, these predicted data are compared with the actual observed data, and the rationality of the sampling results is evaluated based on the degree of fit between the two (i.e., the likelihood function). During the iteration process, the MCMC algorithm will automatically adjust the sampling parameters to make them more inclined to the area that is highly consistent with the observed data, and finally obtain a set of stable parameter posterior distributions. Through this process, the probability distribution of each parameter value can be determined, and the optimal value of the model parameter and its uncertainty range can be estimated based on this distribution. The final result is a set of posterior distributions of all model parameters, which not only provides the most likely parameter estimates, but also reflects the confidence intervals of model parameters under different assumptions. This provides an important reference for the accuracy and robustness of seismic wave propagation models, making the prediction results more reliable.

[0092] S204: Based on the parameter distribution obtained by sampling, parameter estimation and uncertainty analysis are performed to obtain the probability distribution of the geological structure.

[0093] In some embodiments, see Figure 4 As shown, a schematic flow chart of a Bayesian full waveform inversion provided by an embodiment of the present invention includes observing seismic records, constructing a velocity model based on prior information, performing forward simulation, calculating the likelihood function, and judging whether it meets the requirements. If so, the velocity model is saved, and the posterior probability distribution is determined using the velocity model. If it does not meet the requirements, a new velocity model is generated to continue the forward simulation.

[0094] The memory effect of the time fractional viscoacoustic wave equation is reasonably simplified through a short-memory fast algorithm, thereby reducing the storage requirements for wave field historical data and achieving efficient numerical solution. Moreover, while ensuring accuracy, the memory usage is greatly reduced, making the time fractional viscoelastic acoustic wave equation more operable and computationally efficient in engineering applications.

[0095] See also Figure 5 FIG. 1 is a schematic diagram of a seismic wave detection system provided by an embodiment of the present invention, comprising: The acquisition unit 50 is used to arrange the seismic sources and detectors in advance according to the actual working conditions, and to excite seismic waves at each seismic source position and to collect seismic wave data through the detectors; A simulation unit 51 is used to solve the time fractional viscoelastic acoustic wave equation based on the seismic wave data using a short memory fast algorithm to simulate the propagation behavior of seismic waves in real stratum media; The inversion unit 52 is used to perform Bayesian full waveform inversion and reverse time migration imaging based on the seismic wave data to obtain the formation structure and formation characteristics.

[0096] The seismic wave detection system includes a processor and a memory. The acquisition unit, simulation unit, inversion unit, etc. are all stored in the memory as program units, and the processor executes the program units stored in the memory to realize corresponding functions.

[0097] The processor contains a kernel, which calls the corresponding program unit from the memory. One or more kernels can be set, and the calculation time can be reduced and the accuracy can be improved by adjusting the kernel parameters.

[0098] The memory may include non-permanent memory in a computer-readable medium, random access memory (RAM) and / or non-volatile memory, such as read-only memory (ROM) or flash RAM, and the memory includes at least one memory chip.

[0099] An embodiment of the present invention provides a storage medium on which a program is stored. When the program is executed by a processor, the seismic wave detection method is implemented.

[0100] An embodiment of the present invention provides an electronic device, the device comprising a processor, a memory, and a program stored in the memory and executable on the processor, wherein the processor implements the following steps when executing the program: arranging seismic sources and detectors in advance according to actual working conditions, exciting seismic waves at various seismic source positions, and collecting seismic wave data through the detectors; solving a time fractional-order viscoelastic acoustic wave equation based on the seismic wave data using a short memory fast algorithm to simulate the propagation behavior of seismic waves in real stratum media; Bayesian full waveform inversion and reverse time migration imaging are performed based on the seismic wave data to obtain the formation structure and formation characteristics. The device in this article can be a server, PC, PAD, mobile phone, etc.

[0101] Those skilled in the art will appreciate that the embodiments of the present application may be provided as methods, systems, or computer program products. Therefore, the present application may adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Moreover, the present application may adopt the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program codes.

[0102] The present application is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or box in the flowchart and / or block diagram, as well as the combination of the processes and / or boxes in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0103] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0104] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process. Figure 1 A process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0105] In a typical configuration, a computing device includes one or more processors (CPU), input / output interfaces, network interfaces, and memory.

[0106] The memory may include non-permanent memory in a computer-readable medium, random access memory (RAM) and / or non-volatile memory in the form of read-only memory (ROM) or flash RAM. The memory is an example of a computer-readable medium.

[0107] Computer readable media include permanent and non-permanent, removable and non-removable media that can be implemented by any method or technology to store information. Information can be computer readable instructions, data structures, program modules or other data. Examples of computer storage media include, but are not limited to, phase change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technology, compact disk read-only memory (CD-ROM), digital versatile disk (DVD) or other optical storage, magnetic cassettes, magnetic tape magnetic disk storage or other magnetic storage devices or any other non-transmission media that can be used to store information that can be accessed by a computing device. As defined herein, computer readable media does not include temporary computer readable media (transitory media), such as modulated data signals and carrier waves.

[0108] It should also be noted that the terms "include", "comprises" or any other variations thereof are intended to cover non-exclusive inclusion, so that a process, method, commodity or device including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, commodity or device. In the absence of more restrictions, the elements defined by the sentence "comprises a ..." do not exclude the existence of other identical elements in the process, method, commodity or device including the elements.

[0109] The above are only embodiments of the present application and are not intended to limit the present application. For those skilled in the art, the present application may have various changes and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application should be included within the scope of the claims of the present application.

Claims

1. A method for detecting seismic waves, characterized in that: include: Arrange seismic sources and geophones in advance according to actual working conditions, stimulate seismic waves at each seismic source location, and collect seismic wave data through geophones; Based on the seismic wave data, a short memory fast algorithm is used to solve the time fractional viscoelastic acoustic wave equation to simulate the propagation behavior of seismic waves in real stratum media; Bayesian full waveform inversion and reverse time migration imaging are performed based on the seismic wave data to obtain the formation structure and formation characteristics.

2. The method for detecting seismic waves according to claim 1, characterized in that: The time fractional viscoelastic acoustic wave equation is as follows: In the formula, , Respectively represent the particles in , The velocity component in the axial direction, is the reference frequency, is the fractional order, represents the phase velocity, is the density of the medium, is the source term acting on the unit volume, , The order is The fractional derivative of Represents stress.

3. The method for detecting seismic waves according to claim 1, characterized in that: Solving the time fractional viscoelastic acoustic wave equation using a short memory fast algorithm includes: introducing an auxiliary function to transform the solution of the time fractional derivative of the time fractional viscoelastic acoustic wave equation into the calculation of a weak singular integral on a semi-infinite domain.

4. The method for detecting seismic waves according to claim 3, characterized in that: The auxiliary function formula is as follows: In the formula is a parameter introduced in the fractional derivative, which is related to the time memory effect or relaxation characteristics. is the integration variable, which represents the historical effect in the time integration, Represents an auxiliary function, represents the first-order derivative of the strain function with respect to time, which represents the rate of change of strain with time. It represents the second-order derivative of the strain function with respect to time, indicating the acceleration or changing trend of the strain.

5. The method for detecting seismic waves according to claim 1, characterized in that: The short memory fast algorithm formula is as follows: In the formula, is the stress function, indicating the position x and time t The stress distribution under represents the phase velocity, is the number of storage layers of the auxiliary function, represents the scaling factor, is the weight factor, which is related to the fractional derivative and represents the weight coefficient of a certain layer in the calculation. is the auxiliary variable of the fractional derivative, indicating the j The node position or value of the layer, is the auxiliary function, and the strain function Correlation describes the memory effect in the fractional-order model.

6. The method for detecting seismic waves according to claim 1, characterized in that: Bayesian full waveform inversion and reverse time migration imaging are performed based on the seismic wave data to obtain the formation structure and formation characteristics, including: Obtain observation data through detectors arranged in advance according to actual working conditions; An initial geological model is constructed based on geological knowledge, historical data and prior information, and numerical simulation is performed using the initial geological model according to a short-memory fast algorithm to calculate theoretical seismic waveforms; Computers are used to simulate seismic wave propagation forward and short-memory fast algorithms are used to numerically solve the time fractional viscoelastic acoustic wave equation to predict the propagation behavior of seismic waves in complex underground structures. Calculate the error between the observed data and the theoretical data obtained by forward modeling, and determine whether the residual between the observed data and the forward modeling data meets the preset threshold requirements. If so, obtain the target geological model; if not, use the Bayesian statistical method to sample the geological model parameters; Based on the parameter distribution obtained by sampling, parameter estimation and uncertainty analysis are performed to obtain the probability distribution of geological structures.

7. The method for detecting seismic waves according to claim 6, characterized in that: During the simulation of seismic wave propagation, the region is gridded in two dimensions, where the grid spacing in the horizontal and vertical directions is set equal; On this basis of gridding, the seismic wave field is analyzed in the frequency domain through discrete Fourier transform to realize the calculation of spatial derivatives in the frequency domain, and the fast Fourier transform algorithm is used to improve the calculation efficiency.

8. A seismic wave detection system, characterized in that: include: The acquisition unit is used to arrange the seismic sources and detectors in advance according to the actual working conditions, and to excite seismic waves at each seismic source position and collect seismic wave data through the detectors; A simulation unit, used for solving the time fractional viscoelastic acoustic wave equation based on the seismic wave data by using a short memory fast algorithm to simulate the propagation behavior of seismic waves in real stratum media; The inversion unit is used to perform Bayesian full waveform inversion and reverse time migration imaging based on the seismic wave data to obtain the formation structure and formation characteristics.

9. An electronic device 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 the program, the steps of the seismic wave detection method as described in any one of claims 1-7 are implemented.

10. A non-transitory computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the seismic wave detection method as described in any one of claims 1 to 7 are implemented.

Images