Variational inference full waveform inversion method and device based on geological prior information

Abstract

The specification relates to the technical field of seismic inversion, and specifically discloses a variational inference full waveform inversion method and device based on geological prior information, wherein the method comprises: constructing multiple filters; applying the constructed multiple filters to a seismic migration profile to calculate a correlation coefficient matrix; fitting a theoretical variation function to the calculated correlation coefficient matrix according to well logging data to obtain a low-rank approximation correlation coefficient matrix; randomly generating multiple velocity samples based on the low-rank approximation correlation coefficient matrix; calculating the gradient of an objective function according to the multiple velocity samples; calculating the gradient of a posterior distribution by using a Bayesian formula according to the gradient of the objective function; and performing multiple iterative updates on the multiple velocity samples by using a Stein variational gradient descent algorithm based on the gradient of the full waveform inversion posterior distribution, and outputting a target velocity model obtained by inverting the multiple velocity samples. The above scheme reduces the multi-solution of inversion and improves the inversion efficiency.

Description

Technical Field

[0001] This specification relates to the field of seismic inversion technology, and in particular to a variational inference full waveform inversion method and apparatus based on geological prior information. Background Technology

[0002] High-precision seismic imaging requires high-precision depth-domain velocity models. Full waveform inversion (FWI) comprehensively utilizes information such as seismic wave travel time, amplitude, and phase, effectively compensating for the lack of mid- and low-frequency components in conventional velocity modeling, and playing a crucial bridging role between travel-time tomography and migration imaging. Thanks to significant improvements in computer performance, FWI has received considerable attention in recent years. However, limited by observation systems, data bandwidth, the approximation degree of forward modeling operators, noise interference, and the strong nonlinearity of the wave equation, FWI is essentially a strongly nonlinear inverse problem under conditions of severe information deficiency, exhibiting significant multiple solutions. Currently, most FWI methods employ local optimization algorithms, and their inversion results heavily depend on the initial model. When the initial model is inaccurate, FWI will get trapped in local minima. For the information-deficient and multi-solution inverse problem, rather than simply improving the algorithm to obtain the "optimal solution" in the sense of data matching, it is more reasonable to solve the inverse problem within a Bayesian framework, obtain the posterior probability density distribution of the model, and evaluate the uncertainty of the FWI solution.

[0003] Compared to traditional deterministic FWI, Bayesian FWI offers at least two advantages. First, it can easily utilize various prior information (such as geological knowledge) to compensate for insufficient information in seismic data, thus narrowing the FWI solution space. Second, it provides an assessment of the uncertainty of velocity models, which is beneficial for evaluating exploration risks and subsequent geological interpretation. With the significant improvement in computer computing performance, research on Bayesian FWI methods has gradually developed. Bayesian FWI methods can fully utilize various prior information to narrow the solution space and provide an assessment of solution uncertainty, making them an effective means of solving the problem of FWI ambiguity. However, current Bayesian inversion methods suffer from high ambiguity and low computational efficiency.

[0004] There is currently no effective solution to the above problems. Summary of the Invention

[0005] This specification provides a variational inference full waveform inversion method and apparatus based on geological prior information to improve the computational efficiency of full waveform inversion and reduce the ambiguity of inversion.

[0006] This specification provides an embodiment of a variational inference full waveform inversion method based on geological prior information, including:

[0007] Construct multiple filters; apply the constructed filters to the seismic migration profile to calculate the correlation coefficient matrix;

[0008] The correlation coefficient matrix is ​​approximated by a low-rank approximation of the theoretical variation function fitted to the well logging data to obtain the low-rank approximation correlation coefficient matrix.

[0009] Based on the correlation coefficient matrix after the low-rank approximation, multiple velocity samples are randomly generated.

[0010] Based on the multiple velocity samples, the gradient of the full waveform inversion objective function is calculated; based on the gradient of the objective function, the gradient of the full waveform inversion posterior distribution is calculated using Bayes' theorem; based on the gradient of the full waveform inversion posterior distribution, the Stein variational gradient descent algorithm is used to iteratively update the multiple velocity samples multiple times, and the target velocity model obtained by inverting the multiple velocity samples is output.

[0011] In one embodiment, after iteratively updating the plurality of velocity samples using the Stein variational gradient descent algorithm based on the gradient of the full waveform inversion posterior distribution and outputting the target velocity model obtained from the inversion of the plurality of velocity samples, the method further includes:

[0012] The posterior distribution of the target velocity model is determined to evaluate the uncertainty of the inversion results.

[0013] In one embodiment, multiple filters are constructed, including:

[0014] For the two-dimensional offset profile, six filters are set: vertical mean filter, horizontal mean filter, vertical gradient filter, horizontal gradient filter, vertical curvature filter, and horizontal curvature filter.

[0015] In one embodiment, the constructed multiple filters are applied to the seismic migration profile to calculate the correlation coefficient matrix, including:

[0016] The geological patterns captured from the seismic migration profiles are converted into score values ​​using the various filters described above.

[0017] The multiple score values ​​corresponding to the geological patterns around each node in the seismic migration profile are defined as the score vector corresponding to each node.

[0018] Based on the score vectors corresponding to any two nodes in the seismic migration profile, the similarity of the geological patterns around any two nodes is calculated to obtain the multi-point pattern correlation between any two nodes.

[0019] A correlation coefficient matrix is ​​generated based on the multi-point pattern correlation between any two nodes.

[0020] In one embodiment, the Stein variational gradient descent algorithm is used to iteratively update the plurality of velocity samples multiple times, including:

[0021] Given a step size, update multiple velocity samples using the Stein variational gradient descent algorithm;

[0022] Based on the multiple velocity samples, calculate the gradient of the full waveform inversion objective function; based on the gradient of the objective function, calculate the gradient of the full waveform inversion posterior distribution using Bayes' theorem;

[0023] Repeat the above steps until the number of iterations reaches the set maximum number of iterations or until the posterior distribution of the full waveform inversion tends to stabilize.

[0024] In one embodiment, the kernel function of the Stein variational gradient descent algorithm is a radial basis function.

[0025] This specification also provides an embodiment of a variational inference full waveform inversion device based on geological prior information, including:

[0026] The calculation module is used to construct various filters; it is also used to apply the constructed filters to the seismic migration profile to calculate the correlation coefficient matrix.

[0027] The approximation module is used to perform a low-rank approximation on the calculated correlation coefficient matrix based on the fitting theoretical variogram function of the well logging data, so as to obtain the low-rank approximated correlation coefficient matrix.

[0028] The generation module is used to randomly generate multiple velocity samples based on the correlation coefficient matrix after the low-rank approximation.

[0029] The inversion module is used to calculate the gradient of the full waveform inversion objective function based on the multiple velocity samples; it is also used to calculate the gradient of the full waveform inversion posterior distribution using Bayes' theorem based on the gradient of the objective function; and it is also used to perform multiple iterations of updating the multiple velocity samples using the Stein variational gradient descent algorithm based on the gradient of the full waveform inversion posterior distribution, and output the target velocity model obtained by inverting the multiple velocity samples.

[0030] In one embodiment, the calculation module is specifically used to: convert the geological patterns captured from the seismic migration profile into score values ​​using the multiple filters; define multiple score values ​​corresponding to the geological patterns around each node in the seismic migration profile as score vectors corresponding to each node; calculate the similarity of the geological patterns around any two nodes based on the score vectors corresponding to any two nodes in the seismic migration profile to obtain the multi-point pattern correlation between any two nodes; and generate a correlation coefficient matrix based on the multi-point pattern correlation between any two nodes.

[0031] This specification also provides a computer device, including a processor and a memory for storing processor-executable instructions, wherein the processor executes the instructions to implement the steps of the variational inference full waveform inversion method based on geological prior information described in any of the above embodiments.

[0032] This specification also provides a computer-readable storage medium storing computer instructions that, when executed, implement the steps of the variational inference full waveform inversion method based on geological prior information described in any of the above embodiments.

[0033] This specification provides a variational inference full waveform inversion method based on geological prior information. It constructs multiple filters and applies them to seismic migration profiles to calculate a correlation coefficient matrix. A low-rank approximation is performed on the calculated correlation coefficient matrix using a theoretical variogram fitting to well logging data, resulting in a low-rank approximated correlation coefficient matrix. Based on this low-rank approximated correlation coefficient matrix, multiple velocity samples are randomly generated. The gradient of the full waveform inversion objective function is calculated based on these velocity samples. The gradient of the objective function is then used to calculate the gradient of the full waveform inversion posterior distribution using Bayes' theorem. Based on the gradient of the full waveform inversion posterior distribution, the Stein variational gradient descent algorithm is used to iteratively update the multiple velocity samples, outputting the target velocity model obtained from the inversion of these multiple velocity samples. The above scheme proposes a correlation coefficient matrix of geostatistical models with low-rank approximation, which removes spurious correlations and reduces memory and computational consumption. Based on the low-rank approximation correlation coefficient matrix, geologically significant random initial samples are randomly generated to make the prior information more accurate. Then, based on the geostatistical prior information, the Stein variational gradient descent algorithm is used to perform efficient variational inference full waveform inversion, which greatly reduces the ambiguity of the inversion and improves the inversion efficiency. Attached Figure Description

[0034] The accompanying drawings, which are included to provide a further understanding of this specification and form part of it, do not constitute a limitation thereof. In the drawings:

[0035] Figure 1 A flowchart of a variational inference full waveform inversion method based on geological prior information in one embodiment of this specification is shown;

[0036] Figure 2 The Marmousi model is shown in one embodiment of this specification;

[0037] Figure 3 A reverse-time offset profile of the Marmousi model in one embodiment of this specification is shown;

[0038] Figure 4 Six filters are shown in one embodiment of this specification;

[0039] Figure 5 Six filter score profiles are shown in one embodiment of this specification;

[0040] Figure 6 The variation function in one embodiment of this specification is shown;

[0041] Figure 7 A multi-point pattern correlation matrix is ​​shown in one embodiment of this specification;

[0042] Figure 8 This specification illustrates an embodiment of establishing a stochastic model using geostatistical priors;

[0043] Figure 9 The full waveform inversion result of variational inference based on uniform prior information in one embodiment of this specification is shown;

[0044] Figure 10 The mean and standard deviation of the statistical inversion results in one embodiment of this specification are shown;

[0045] Figure 11 The mean and standard deviation of the statistical inversion results in one embodiment of this specification are shown;

[0046] Figure 12 The full waveform inversion result of variational inference based on geostatistical prior information is shown in one embodiment of this specification;

[0047] Figure 13 A schematic diagram of a variational inference full waveform inversion device based on geological prior information is shown in one embodiment of this specification;

[0048] Figure 14 A schematic diagram of a computer device according to one embodiment of this specification is shown. Detailed Implementation

[0049] The principles and spirit of this specification will now be described with reference to several exemplary embodiments. It should be understood that these embodiments are given merely to enable those skilled in the art to better understand and implement this specification, and are not intended to limit the scope of this specification in any way. Rather, these embodiments are provided to make this disclosure more thorough and complete, and to fully convey the scope of this disclosure to those skilled in the art.

[0050] Those skilled in the art will recognize that the embodiments described in this specification can be implemented as a system, apparatus, method, or computer program product. Therefore, the disclosure of this specification can be specifically implemented in the following forms: entirely hardware, entirely software (including firmware, resident software, microcode, etc.), or a combination of hardware and software.

[0051] Bayesian inversion estimates the posterior probability density distribution of a model using its prior probability density distribution, making it an effective means of evaluating the uncertainty of FWI solutions. Bayesian inversion methods mainly include sampling methods and variational inference methods. Sampling methods, represented by Markov-Monte Carlo methods, suffer from low sample acceptance rates and high computational costs. Another approach to solving Bayesian problems is variational inference, which continuously updates the prior probability density distribution by minimizing a certain objective functional, making it sufficiently close to the desired posterior probability density distribution. Variational inference can utilize the gradient information of the objective function, thus its computational efficiency is much higher than that of sampling methods.

[0052] Stein's variational gradient descent (SVGD) method is commonly used for estimating complex multimodal distribution problems. It is easily parallelized and computationally efficient. This algorithm generates a series of samples based on the prior distribution and correlates the gradient information of different samples using a kernel function, preventing the current sample from getting trapped in local extrema. The samples, updated through multiple iterations, constitute an approximate estimate of the posterior probability density distribution. The variational inference full waveform inversion method based on geological prior information in the embodiments of this specification can apply the SVGD algorithm to FWI and can also evaluate the uncertainty of the inversion results.

[0053] In Bayesian inversion, prior information is crucial to the posterior inversion results. Currently, Bayesian full-waveform inversion employs extremely broad prior information, meaning each sampling point satisfies a uniform distribution with the same variance, which significantly increases computational complexity. Therefore, the method in this specification incorporates geostatistical prior information into the Bayesian FWI. Applying a multi-point geostatistical filter to the seismic migration profile yields the correlation coefficient matrix of geological models surrounding any two points. A low-rank approximation is then applied to the correlation coefficient matrix, retaining only the correlation coefficients when the distance between two points is less than the range, thus removing spurious correlations and significantly reducing memory and computational costs. Multiplying this matrix by a random vector yields geologically meaningful random samples.

[0054] The method in the embodiments of this specification utilizes prior information from geostatistics to generate multiple geologically significant random samples, and uses the SVGD algorithm to perform efficient variational inference full waveform inversion, which greatly reduces the search space for solutions, reduces the multiple solutions in the inversion, and can also determine the uncertainty estimate of the solution.

[0055] Based on this, the embodiments of this specification provide a variational inference full waveform inversion method based on geological prior information. Figure 1 A flowchart of a variational inference full-waveform inversion method based on geological prior information according to an embodiment of this specification is shown. While this specification provides method operation steps or apparatus structures as shown in the following embodiments or figures, more or fewer operation steps or module units may be included in the method or apparatus based on conventional or non-inventive effort. In steps or structures where there is no logically necessary causal relationship, the execution order of these steps or the module structure of the apparatus is not limited to the execution order or module structure described in the embodiments and figures of this specification. When the method or module structure is applied in actual devices or end products, it can be executed sequentially or in parallel according to the method or module structure shown in the embodiments or figures (e.g., in a parallel processor or multi-threaded processing environment, or even a distributed processing environment).

[0056] Specifically, such as Figure 1 As shown, the variational inference full waveform inversion method based on geological prior information provided in one embodiment of this specification may include the following steps.

[0057] Step S101: Construct multiple filters; apply the constructed multiple filters to the seismic migration profile to calculate the correlation coefficient matrix.

[0058] The method in this embodiment can be applied to computer equipment. Various filters can be constructed. These filters are used to extract key information related to geological models from seismic migration profiles. Different types of filters extract different key information. By constructing multiple filters, various key information related to geological models can be extracted from seismic migration profiles.

[0059] When the constructed filters are applied to seismic migration profiles, they are essentially used to analyze and process the geological patterns within the profiles. Through filter operations, key information related to the geological patterns can be captured from the seismic profiles, and based on this key information, a correlation coefficient matrix can be calculated. The correlation coefficient matrix characterizes the degree of correlation between geological patterns at different nodes in the seismic migration profile, providing crucial foundational information for subsequent full-waveform inversion.

[0060] Step S102: Based on the fitting theoretical variation function of the well logging data, perform a low-rank approximation on the calculated correlation coefficient matrix to obtain the low-rank approximation correlation coefficient matrix.

[0061] Well logging data contains rich geological information. By analyzing and processing this data, a theoretical variogram can be fitted. The theoretical variogram is a function used to describe the spatial variation characteristics of a variable; in geostatistics, it represents the incremental semivariogram of a variable between two points at a certain distance. After obtaining the theoretical variogram, it is applied to the previously calculated correlation coefficient matrix. The correlation coefficient matrix reflects the degree of correlation between geological patterns between different nodes in a seismic migration profile, but the original correlation coefficient matrix may have some problems. Because it represents that all points in space are correlated, this does not conform to the actual laws of subsurface geology, and such a matrix is ​​usually very large, requiring a lot of computational and storage resources. Using the theoretical variogram to perform a low-rank approximation of the correlation coefficient matrix involves processing the correlation coefficient matrix based on the spatial correlation law reflected by the variogram. Specifically, only the correlation coefficients when the distance between two points is less than the range (a distance value determined by the variogram; when the distance between two points is greater than this value, they are generally considered not correlated), are retained; otherwise, the correlation coefficients are set to zero. This approach essentially simplifies the matrix, removing unreasonable and spurious correlations. This allows the correlation coefficient matrix after the low-rank approximation to more accurately reflect the true correlations of subsurface geological models, thus grounding the entire inversion process in a more realistic foundation.

[0062] After performing the low-rank approximation, the correlation coefficients between a large number of points in the matrix whose distances are greater than the range are set to zero, and the size of the matrix is ​​significantly reduced. This greatly reduces the memory usage and the amount of computation during the calculation process, improves the computational efficiency, and makes full waveform inversion more feasible in practical applications, especially when processing large-scale data.

[0063] Meanwhile, a more accurate correlation coefficient matrix provides a more reliable foundation for subsequent inversion steps. The velocity model generated by perturbing the correlation coefficient matrix after the low-rank approximation is more consistent with geological reality. In the process of calculating the gradient of the objective function for the full waveform inversion and calculating the gradient of the posterior distribution using Bayes' theorem, it can more accurately reflect the relationship between the velocity model and the observation data. This helps guide the inversion towards a direction that is closer to the true velocity model, thereby improving the accuracy of the inversion results and reducing the ambiguity of the inversion.

[0064] After low-rank approximation, the resulting correlation coefficient matrix better reflects the actual underground geology, while significantly reducing memory and computational costs. This low-rank approximated correlation coefficient matrix will be used in subsequent steps, such as randomly generating multiple velocity samples, providing more accurate and efficient basic data support for full waveform inversion. When iteratively updating the velocity samples using the Stein variational gradient descent algorithm, samples generated based on a more reasonable correlation coefficient matrix can converge to a stable posterior distribution faster, reducing unnecessary iterations, accelerating the convergence speed of the inversion, and saving computation time.

[0065] Step S103: Based on the correlation coefficient matrix after the low-rank approximation, randomly generate multiple velocity samples.

[0066] After obtaining the correlation coefficient matrix after the low-rank approximation, multiple velocity samples can be randomly generated using specific methods. Generally, some random factors are incorporated to implement the perturbation process.

[0067] In the field of stochastic geological modeling, random samples can be generated by perturbation in the following way:

[0068] m = m smooth +Γ′·R

[0069] Where m represents a column vector composed of the two-dimensional perturbation stochastic velocity model arranged in column order, m = [m1 m2 … m nx ] T m i (i = 1, 2, ..., nx) represents the i-th column vector in the two-dimensional model (i.e., the first nz elements of m constitute the first column of the two-dimensional velocity model, and so on, where nx and nz are the number of columns and rows of the two-dimensional model, respectively, and m... smooth Let R be a vector composed of a smooth background model, and let R be a random vector following the property (0, R). max Uniform distribution in the interval, R max Let Γ′ be the maximum perturbation, and Γ′ be the correlation coefficient matrix after the low-rank approximation.

[0070] By combining the correlation coefficient matrix obtained from the low-rank approximation with random vectors in the above manner, the generated velocity samples contain information about the correlation between geological models (reflected by the correlation coefficient matrix) while also introducing a degree of randomness (provided by the random vectors). These multiple velocity samples are geologically significant, enabling a more realistic simulation of various possible scenarios in the subsurface velocity model and providing rich and reasonable initial conditions for subsequent full-waveform inversion. Following this, further operations such as calculating the gradient of the full-waveform inversion objective function will be performed based on these generated velocity samples to advance the inversion process and obtain a more accurate target velocity model.

[0071] Step S104: Calculate the gradient of the full waveform inversion objective function based on the multiple velocity samples; calculate the gradient of the full waveform inversion posterior distribution using Bayes' theorem based on the gradient of the objective function; and iteratively update the multiple velocity samples using the Stein variational gradient descent algorithm based on the gradient of the full waveform inversion posterior distribution to output the target velocity model obtained by inverting the multiple velocity samples.

[0072] Specifically, the gradient of the objective function for full waveform inversion can be calculated based on the multiple velocity samples. The goal of full waveform inversion is to find an optimal velocity model that makes the synthetic seismic data obtained from forward modeling based on this model match the actual observation data as closely as possible. The objective function is used to measure this degree of matching and can usually be expressed as the error between the synthetic and observation data, such as using the mean square error of the two. The gradient represents the rate of change of the objective function in each direction. By calculating the gradient of the objective function for multiple velocity samples, it is possible to determine how the objective function changes with the velocity model.

[0073] In one embodiment, the acoustic wave equation and its adjoint equation can be solved using the finite difference method for each velocity sample. The gradient of the objective function for full waveform inversion can then be calculated using the zero-delay cross-correlation between the source wavefield and the adjoint wavefield. The acoustic wave equation describes the propagation process of seismic waves in a medium, including physical quantities such as the source wavelet, density, seismic wave propagation velocity, pressure field, and the horizontal and vertical components of particle vibration velocity. The finite difference method is a numerical calculation method that discretizes continuous space and time, approximating the acoustic wave equation by calculating the physical quantities at discrete points. This method can simulate the propagation of seismic waves under a given initial velocity sample, obtaining the source wavefield. The adjoint equation is another equation related to the acoustic wave equation and plays a crucial role in calculating the gradient of the objective function. Solving the adjoint equation using the finite difference method also yields the adjoint wavefield. The calculation of the adjoint wavefield is opposite to the propagation direction of the source wavefield, proceeding from the maximum time to the zero time.

[0074] The zero-delay cross-correlation between the source wavefield and the accompanying wavefield refers to the zero-delay cross-correlation between the source wavefield propagating forward in time and the accompanying wavefield propagating backward in time. This calculation method yields the gradient of the full-waveform inversion objective function with respect to the velocity model, providing an important basis for subsequently updating the velocity model using this gradient to gradually approximate the optimal velocity model.

[0075] The gradient of the posterior distribution of the full waveform inversion can be calculated using Bayes' theorem based on the gradient of the objective function. Bayesian inversion estimates the posterior probability density distribution p(m|d) of the model using the prior probability density distribution p(m) of the velocity model m. ocbsThis process can be described by the following formula:

[0076]

[0077] Where d obs For observation data, p(d) obs p(d) is the normalization constant. obs |m) is the likelihood function.

[0078] By combining the gradient of the objective function, Bayes' theorem is used to calculate the gradient of the posterior distribution. The gradient of the posterior distribution reflects the variation of the posterior probability density distribution in the velocity model space. Considering prior information (such as the prior probability density distribution obtained through geostatistical methods) and observational data, calculating the gradient of the posterior distribution provides a more comprehensive understanding of the uncertainties and probabilities of the velocity model, offering more accurate guidance for subsequent iterative updates.

[0079] Based on the gradient of the posterior distribution obtained from the full waveform inversion, the Stein variational gradient descent algorithm can be used to iteratively update the multiple velocity samples, outputting the target velocity model obtained from the inversion of the multiple velocity samples. The Stein variational gradient descent algorithm is an algorithm for approximating complex posterior distributions, possessing advantages such as the ability to estimate complex multimodal distribution problems, ease of parallel computation, and high computational efficiency. By iteratively updating multiple velocity samples, these samples gradually approximate the posterior probability density distribution.

[0080] In the above embodiments, a correlation coefficient matrix of geostatistical models with low-rank approximation is proposed, which removes spurious correlations and reduces memory and computational consumption. Based on the low-rank approximation correlation coefficient matrix, geologically significant random initial samples are randomly generated to make the prior information more accurate. Then, based on the geostatistical prior information, the Stein variational gradient descent algorithm is used to perform efficient variational inference full waveform inversion, which greatly reduces the ambiguity of the inversion and improves the inversion efficiency.

[0081] In some embodiments of this specification, multiple filters may be constructed, including: setting six filters for a two-dimensional offset profile; the six filters include: a vertical mean filter, a horizontal mean filter, a vertical gradient filter, a horizontal gradient filter, a vertical curvature filter, and a horizontal curvature filter.

[0082] In this embodiment, the vertical mean filter, horizontal mean filter, vertical gradient filter, horizontal gradient filter, vertical curvature filter, and horizontal curvature filter set for two-dimensional migration profiles can detect geological patterns in seismic migration profiles from different angles. Each filter has its own unique detection function, comprehensively capturing the characteristics of geological patterns in different aspects. The vertical mean filter and horizontal mean filter can detect the center position of the image, the vertical gradient filter and horizontal gradient filter can detect the edge patterns of the image, and the vertical curvature filter and horizontal curvature filter are used to detect changes in curvature. Through the combination of these filters with different functions, information about geological patterns can be obtained from multiple dimensions, making the understanding of geological patterns more comprehensive and in-depth, and avoiding the information loss that may result from detecting from only a single angle.

[0083] Because underground geological structures are often highly complex, a single filter cannot accurately characterize their features. The six filters in this embodiment probe geological patterns from different angles, enabling them to better adapt to complex and varied geological structures. They can effectively detect and analyze aspects such as the center location, edge features, and curvature changes of geological structures, thus providing valuable information for full waveform inversion even under complex geological conditions.

[0084] Since the six filters can detect geological patterns from multiple angles, the key information obtained contains richer geological information and can more realistically reflect the actual relationship between geological patterns around different nodes, providing a reliable foundation for the subsequent construction of an accurate correlation coefficient matrix.

[0085] An accurate correlation coefficient matrix is ​​crucial for full waveform inversion. The correlation coefficient matrix constructed using comprehensive geological model information obtained through these six filters better reflects the actual geological conditions, thus providing more accurate prior information for the inversion. In full waveform inversion, the accuracy of prior information helps to narrow the solution search space, reduce the number of possible solutions, decrease the ambiguity of the inversion, and make the inversion results more likely to approximate the true velocity model.

[0086] Understandably, other numbers of filters can also be set, such as two, three, four, or eight filters. The type of filter can be selected according to actual needs.

[0087] In some embodiments of this specification, applying the constructed multiple filters to a seismic migration profile to calculate a correlation coefficient matrix may include: converting geological patterns captured from the seismic migration profile into score values ​​using the multiple filters; defining multiple score values ​​corresponding to geological patterns around each node in the seismic migration profile as score vectors corresponding to each node; calculating the similarity of geological patterns around any two nodes based on the score vectors corresponding to any two nodes in the seismic migration profile to obtain multi-point pattern correlation between any two nodes; and generating a correlation coefficient matrix based on the multi-point pattern correlation between any two nodes.

[0088] Specifically, in the full waveform inversion based on prior geological information, various filters (such as vertical mean filter, horizontal mean filter, vertical gradient filter, horizontal gradient filter, vertical curvature filter, and horizontal curvature filter) are constructed with different functions, enabling the detection of geological patterns in seismic migration profiles from different angles. When these filters are applied to the seismic migration profile, each filter performs specific calculations and analyses on the geological patterns in the profile, thereby transforming complex geological patterns into specific score values. For example, the vertical mean filter might give a score based on the mean characteristics in the vertical direction around a node, while the horizontal gradient filter would give a corresponding score based on the gradient characteristics in the horizontal direction. These score values ​​reflect the characteristics of the geological patterns from different perspectives.

[0089] After obtaining the scores of the geological models surrounding each node under various filters, these multiple scores are combined to define the score vector corresponding to that node. Each score vector can be seen as a quantitative representation of the characteristics of the geological models surrounding that node. Through the score vector, the comprehensive characteristics of the geological models surrounding the node can be described more comprehensively, rather than just information from a single dimension.

[0090] Based on the score vectors corresponding to any two nodes in a seismic migration profile, appropriate methods can be used to calculate their similarity. For example, distance metrics between vectors (such as Euclidean distance, cosine similarity, etc.) can be used to measure the similarity between two score vectors. The higher the similarity, the more similar the geological patterns around the two nodes are. By calculating the similarity between any two nodes, the multi-point pattern correlation between any two nodes is obtained. This correlation not only reflects the relationship between the two nodes themselves, but also represents the degree of association between the geological patterns around them.

[0091] Based on the multi-point model correlation between any two nodes, these correlation values ​​are arranged into a matrix according to certain rules, thus generating the correlation coefficient matrix. The elements of this matrix represent the strength of the correlation between geological models between any two nodes in the seismic migration profile. The correlation coefficient matrix provides important prior information for subsequent full-waveform inversion, helping us better understand the association between geological models at different locations in the seismic migration profile. This allows for a more reasonable consideration of the continuity and correlation of geological structures during the inversion process, improving the accuracy and reliability of the inversion results.

[0092] Using the methods described above, geological model information from seismic migration profiles can be quantified and organized to obtain a correlation coefficient matrix that reflects the interrelationships of geological models. This matrix plays a crucial constraining role in full-waveform inversion, guiding the inversion process towards a direction that better reflects actual geological conditions, reducing the ambiguity of inversion, and improving inversion efficiency and accuracy.

[0093] In some embodiments of this specification, the Stein variational gradient descent algorithm is used to iteratively update the plurality of velocity samples, which may include: given a step size and updating the plurality of velocity samples using the Stein variational gradient descent algorithm; calculating the gradient of the full waveform inversion objective function based on the plurality of velocity samples; calculating the gradient of the full waveform inversion posterior distribution using Bayes' theorem based on the gradient of the objective function; repeating the above steps until the number of iterations reaches the set maximum number of iterations or until the full waveform inversion posterior distribution tends to stabilize.

[0094] A given step size can be used to update multiple velocity samples based on the gradient of the objective function and the gradient of the posterior distribution. After each update, the gradient of the full waveform inversion objective function and the gradient of the posterior distribution are recalculated, and the velocity samples are updated again. This process is repeated until a stopping condition is met, such as the number of iterations reaching the set maximum number of iterations or the full waveform inversion posterior distribution tending to stabilize. Finally, after multiple iterations, the output multiple velocity samples yield the inverted target velocity model.

[0095] Based on the updated velocity samples, the gradient of the objective function for full waveform inversion is calculated using appropriate mathematical methods. The objective function measures the degree of fit between the velocity model and the actual observation data, while the gradient represents the direction and rate of change of the objective function in the velocity model space. For example, in full waveform inversion, numerical simulation methods based on the acoustic wave equations, combined with techniques such as the finite difference method, might be used to calculate the gradient of the objective function with respect to the velocity model. By calculating the gradient, the changing trend of the current velocity samples on the objective function can be understood, providing guidance for subsequent updates.

[0096] Based on the gradient of the objective function calculated earlier, the gradient of the posterior distribution obtained from the full waveform inversion is calculated using Bayes' theorem. Within the Bayesian framework, the posterior distribution integrates prior information and observed data, reflecting the probability distribution of the velocity model given the observed data. By calculating the gradient of the posterior distribution, we can further understand the variation of the velocity model in the probability space and more comprehensively consider the uncertainties of the velocity model. For example, by combining the gradient of the objective function with the prior probability density distribution and the likelihood function using Bayes' theorem, the expression for the gradient of the posterior distribution is derived and calculated.

[0097] Repeat the aforementioned steps, that is, update multiple velocity samples again using the Stein variational gradient descent algorithm, and then calculate the gradient of the objective function and the gradient of the posterior distribution. Continue this iterative process until a predetermined stopping condition is met. The stopping condition can be that the number of iterations reaches a pre-set maximum number of iterations, a simple and direct way to stop, ensuring the algorithm terminates within a certain computational resource limit; or it can be that the posterior distribution of the full waveform inversion tends to stabilize, i.e., the change in the posterior distribution is very small, indicating that the algorithm has converged to a relatively stable state. At this point, the obtained velocity samples are considered close to the optimal solution, and iteration can stop.

[0098] Through the above method, Stein's variational gradient descent algorithm can continuously adjust the velocity samples, gradually approaching the optimal solution of the full waveform inversion, while taking into account the uncertainties of the velocity model. This method has significant application value in full waveform inversion, improving the accuracy and reliability of the inversion results and providing more accurate velocity models for fields such as geological exploration.

[0099] In some embodiments of this specification, the kernel function of the Stein variational gradient descent algorithm is a radial basis function.

[0100] Specifically, in Stein's Variational Gradient Descent (SVGD) algorithm, the kernel function plays a crucial role. The kernel function maps data from a low-dimensional space to a high-dimensional space. This mapping allows for more complex calculations and analyses in the high-dimensional space without explicitly calculating the coordinates. The radial basis function k(m,m′) is shown below:

[0101]

[0102] Where m and m′ are two sample points in the sample space, h>0, and is a scaling constant, which is taken in this embodiment. in The median of the relative distances among all samples, where n is the number of samples. It represents the square of the Euclidean distance between two points.

[0103] In the SVGD algorithm, radial basis functions are used as kernel functions to correlate gradient information from different samples. Specifically, when calculating the sample update direction, the similarity between different samples is calculated using the kernel function, and the gradient information of samples with high similarity is fused. This allows the algorithm to consider the interrelationships between samples when updating samples, avoiding the current sample from getting trapped in local maxima.

[0104] Radial basis functions (RBFs) can map data to a high-dimensional space, making linearly inseparable data in low-dimensional space linearly separable or easier to process in high-dimensional space. In full waveform inversion, the posterior distribution of the velocity model is usually a complex multimodal distribution. The SVGD algorithm, which uses RBFs as its kernel function, can better approximate this complex posterior distribution, thereby improving the accuracy of the inversion. Compared to some other complex kernel functions, the calculation of RBFs is relatively simple. It only involves distance calculation and exponential operation between sample points, which can guarantee computational efficiency to a certain extent in the case of large-scale samples. This is very important for handling a large number of velocity samples in full waveform inversion, because full waveform inversion itself is computationally intensive and requires efficient algorithms. The proportionality constant in the RBF can be adjusted according to the characteristics of the data. By reasonably adjusting the value, the width of the kernel function can be controlled, thereby affecting the similarity calculation between samples and the convergence speed of the algorithm. In this embodiment, the value determined based on the median of the relative distance between samples is an effective way to adjust parameters according to data characteristics, enabling the algorithm to better adapt to specific full waveform inversion problems.

[0105] In some embodiments of this specification, after iteratively updating the plurality of velocity samples using the Stein variational gradient descent algorithm based on the gradient of the posterior distribution of the full waveform inversion and outputting the target velocity model obtained by inverting the plurality of velocity samples, the method may further include: determining the posterior distribution of the target velocity model to evaluate the uncertainty of the inversion results.

[0106] Specifically, after obtaining the target velocity model, its corresponding posterior distribution is further determined. The posterior distribution describes the probability distribution of different values ​​of the velocity model after considering all observation data and prior information. The posterior distribution can be determined by analyzing the distribution of multiple velocity samples obtained during statistical iteration, such as calculating the mean, variance, and other statistics of the samples, or by constructing a probability density function to describe the shape of the posterior distribution. Based on the determined posterior distribution of the target velocity model, the uncertainty of the inversion results can be evaluated. The width and variance of the posterior distribution reflect the degree of uncertainty of the velocity model. A wider posterior distribution indicates greater uncertainty in the values ​​of the velocity model, suggesting the existence of multiple reasonable velocity models; while a narrower posterior distribution indicates that the values ​​of the velocity model are relatively certain. Evaluating the uncertainty of the inversion results helps to better understand the reliability of the inversion results, allowing for more reasonable consideration of the uncertainty factors of the velocity model in subsequent geological interpretation and exploration decisions, thereby reducing exploration risks.

[0107] The various embodiments in this specification are described in a progressive manner. Similar or identical parts between embodiments can be referred to interchangeably. Each embodiment focuses on describing the differences from other embodiments. For details, please refer to the foregoing descriptions of the relevant processing embodiments; they will not be repeated here.

[0108] The foregoing has described specific embodiments of this specification. Other embodiments are within the scope of the appended claims. In some cases, the actions or steps recited in the claims may be performed in a different order than that shown in the embodiments and may still achieve the desired result. Furthermore, the processes depicted in the drawings do not necessarily require the specific or sequential order shown to achieve the desired result. In some embodiments, multitasking and parallel processing are possible or may be advantageous.

[0109] The above method will be described below with reference to a specific embodiment. However, it is worth noting that this specific embodiment is only for better illustration of this specification and does not constitute an improper limitation of this specification.

[0110] This specific embodiment provides a variational inference full waveform inversion method based on geological prior information. The implementation of this variational inference full waveform inversion method based on geological prior information mainly includes the following: 1) setting up a full waveform inversion uncertainty evaluation method based on the Stein variational gradient descent algorithm; 2) setting up geostatistical prior information; 3) statistically analyzing the posterior distribution of the inversion results and evaluating their uncertainty.

[0111] (1) Full waveform inversion method based on SVGD algorithm.

[0112] Bayesian inversion uses the prior probability density distribution p(m) of the velocity model m to estimate the posterior probability density distribution p(m|d) of the model. obs This process can be described by the following formula:

[0113]

[0114] Where, d obs For observation data, p(d) obs p(d) is the normalization constant. obs |m) is the likelihood function.

[0115] Variational inference uses KL divergence to describe the current distribution q(m) and the posterior probability density distribution p(m|d). obs Error of )

[0116] KL[q(m)||p(m|d obs )]=E q [log q(m)]-E q [log p(m|d ols (2)

[0117] Among them, E q Let q represent the expectation of q. The probability distribution q is obtained by optimizing the objective function (2) to minimize the KL divergence. * (m) represents the posterior probability density distribution p(m|d). ocbs The best approximation of ).

[0118] Most algorithms in variational inference assume that the posterior distribution of the model follows a Gaussian distribution; however, the velocity in actual underground rock strata does not satisfy this assumption. The Stein variational gradient descent algorithm samples an initial set of samples from the prior distribution, applies a smoothing transformation to the initial distribution, and minimizes the KL divergence through a local optimization algorithm, thereby approximating a multimodal posterior distribution.

[0119] Assuming the initial probability density distribution of the velocity model m is q(m), the above smoothing transformation is:

[0120]

[0121] Where m∈[m1,m2,…m n [There are n velocity models, where n is the number of samples.] This can be viewed as the iterative update direction vectors corresponding to n velocity models. ε is a small constant, which can be considered as the step size. Let the transformed probability density function be q. T (m), according to mathematical transformation, we know:

[0122]

[0123] Where, trace represents the trace of the matrix. It is called the Stein operator and satisfies:

[0124]

[0125] Where Φ(m) is the vector obtained after the function Φ is applied to the sample. According to the explanation of equation (2), the KL divergence reflects the closeness of the two probability density functions, and equation (4) shows that when there exists Φ(m) that makes the expectation on the right side of equation (4) reach its maximum value, and when ε is close to 0, the left side of equation (4) actually gives the fastest descent direction of the KL objective function, and the probability density function q after updating the sample m along this direction. T (m) is closer to the desired posterior probability density function p(m|d) than the original q(m). obs Therefore, we only need to find the objective function that maximizes the following:

[0126]

[0127] Where F represents the functional space, the above optimization problem is transformed into a Hilbert space for solution, yielding:

[0128]

[0129] Where, m j This represents m1, m2, ... m n In the j-th velocity model (j=1,2,…,n), k(m,m′) is the kernel function, usually taken as the radial basis function:

[0130]

[0131] Where h>0, is a proportionality constant, which this patent takes as... in The median of the relative distances among all samples. Substitute equation (7) into equation (4), and change the expression in equation (3) to... By taking the negative gradient direction of the KL divergence, the initial samples can be updated. The SVGD algorithm can be summarized as follows:

[0132] 1) Generate an initial sample set by sampling based on the prior probability density p(m) of the model.

[0133] 2) Iteratively update the samples, using the following formula:

[0134]

[0135] Where the superscript l represents the iteration number, the subscript i represents the i-th sample (i = 1, 2, ..., n), and ε l This represents the step size of the l-th iteration. and These represent the current and the i-th sample after the iteration update, respectively, with the initial sample corresponding to l=0.

[0136]

[0137] The gradient of the posterior distribution in the above equation is transformed into the gradient of the likelihood function using Bayes' theorem (1). In Bayesian full waveform inversion, it is usually assumed that the likelihood function follows a Gaussian distribution:

[0138]

[0139] Where d syn (m) represents the synthetic seismic data obtained based on m-forward modeling, Σ d Let be the covariance matrix of the data, usually assumed to be a diagonal matrix. In this case, finding the gradient of the likelihood function can be transformed into finding the gradient of the conventional L2 objective function in full waveform inversion:

[0140]

[0141] The above equation can be solved using the adjoint state method. Taking the first-order wave equation as an example:

[0142]

[0143] In the formula, s(t) is the source wavelet, and ρ = ρ(x) is the density, with units of kg / m³. 3 v = v(x) is the propagation speed of the seismic wave. For pressure field, and These represent the horizontal and vertical components of the particle's vibration velocity, respectively. The superscript → indicates the forward propagation wave field along the direction of increasing time, i.e., the wave field is calculated from time zero in the direction of increasing time. x=(x,z) are two-dimensional spatial coordinates, x s =(x s ,z s () represents the spatial location of the source wavelet, and δ is the Dirac-Delta function. The associated wave equation satisfied by the associated wave field obtained using the associated state method is:

[0144]

[0145] Where r represents the accompanying source. and For the accompanying wave field variable, the superscript ← indicates that the wave field propagation is along the direction of decreasing time, that is, calculated from the maximum time to the zero time.

[0146] Solving for the derivative of the objective function with respect to the velocity model v yields the gradient formula used for updating the velocity in the full waveform inversion:

[0147]

[0148] Where T max The maximum recording duration of the earthquake data is expressed in seconds. This formula indicates that the gradient is the source wavefield propagating in the positive time direction. Accompanying wave field propagating in reverse time Zero-latency cross-correlation.

[0149] The SVGD algorithm generates a series of samples based on the prior distribution and then updates the samples based on the gradient of the KL divergence objective function. During the iteration update process, the gradient and prior information of the data items are introduced through Bayes' theorem, and the information of different samples is associated through the kernel function k(m,m′). This association is equivalent to a constraint. When a sample gets stuck in a local extremum, its neighboring samples can help it escape the local extremum. The samples after multiple iterations of updates constitute an approximate estimate of the posterior probability density distribution. This specific embodiment uses it for FWI speed modeling and evaluates the uncertainty of the inversion results. In addition, it can be seen from Equation (10) that the gradient calculation for each sample is performed independently. Therefore, this specific embodiment implements multi-GPU parallel computing for samples from the program design perspective, thereby improving the computational efficiency of variational inference FWI.

[0150] (2) Set up prior information for geostatistics.

[0151] In Bayesian inversion, prior information is crucial to the posterior inversion results. Currently, Bayesian full-waveform inversion employs extremely broad prior information, meaning each sampling point follows a uniform distribution with the same variance, which significantly increases computational complexity. Therefore, generating geologically meaningful prior samples helps narrow the solution search space and accelerates the inversion convergence speed.

[0152] In the field of stochastic geological modeling, random samples can be generated by perturbation in the following way:

[0153] m = m smooth +Γ·R (16)

[0154] Where m represents a column vector composed of two-dimensional perturbation stochastic velocity models arranged in column order, m = [m1m2…m ... nx ] T m1 (i = 1, 2, ..., nx) represents the i-th column vector in the two-dimensional model (i.e., the first nz elements of m constitute the first column of the two-dimensional velocity model, and so on, where nx and nz are the column number and row number of the two-dimensional model, respectively). smoothLet R be a vector composed of a smooth background model, and let R be a random vector following the property (0, R). max Uniform distribution in the interval, R max Let Γ be the maximum disturbance, and let Γ be the matrix of correlation coefficients between any two points in the model.

[0155] However, the correlation between two points is difficult to characterize complex geological structures. In this embodiment, geostatistics and seismic migration profiles are used to efficiently characterize the pattern correlation (i.e., multi-point correlation) between different sampling points.

[0156] Given a seismic migration profile S, and a square filter with radius l, a geological model p obtained around the sampling point S(i,j) in the migration profile using this filter can be expressed as:

[0157]

[0158] in i and j represent the vertical and horizontal coordinates of the sampling point, respectively.

[0159] For the two-dimensional offset profile, a total of six filters are set:

[0160]

[0161] Here, f1, f2, f3, f4, f5, and f6 are the vertical mean, horizontal mean, vertical gradient, horizontal gradient, vertical curvature, and horizontal curvature filters, respectively. f1 and f2 can detect the center position of the image, f3 and f4 can detect the edge patterns of the image, and f5 and f6 can detect changes in curvature.

[0162] These six filters convert the geological patterns captured from seismic profiles into scores: E1(i,h), E2(i,h), E3(i,j), E4(i,h), E5(i,h), and E6(i,h). The expressions for these six scores are:

[0163]

[0164] In the formula, k = 1, 2, ..., 6, i ∈ [m+1, nz-m], j ∈ [m+1, nx-m]. All geological models in the offset profile are converted into score values, and six score values ​​are obtained at each point. Here, the six score values ​​corresponding to the geological models around a node S(i,j) are defined as a score vector, the expression of which is:

[0165] e(i,j)=[E1(i,j),E2(i,j),E3(i,j),E4(i,j),E5(i,j),E6(i,j)] (20)

[0166] Therefore, for any point in the offset profile, a score vector can be obtained using the above method. The similarity between two patterns can be determined based on the score vectors corresponding to two nodes. This similarity represents not only the relationship between the two nodes but also the relationship between the geological patterns surrounding them. Here, this similarity is defined as multi-point pattern correlation, expressed as:

[0167] γ=Corr[e(i,j),e(i′,j′)] (21)

[0168] This formula represents the magnitude of the pattern correlation between S(i,j) and S(i′,j′), where Corr is the correlation coefficient calculation.

[0169] However, the correlation coefficient matrix Γ in equation (15) formed by the correlation coefficients in equation (20) is extremely large, requiring massive computational and storage resources. Furthermore, this matrix represents the correlation between all points in space, which does not conform to the laws of underground geology. Therefore, we introduce the variogram function from two-point geostatistics, defining the variogram function r(h) of variable Z(x) between point x and x+h, denoted as:

[0170]

[0171] Where Var represents variance, and the variogram represents the incremental semivariogram of a variable between two points at a distance h. The experimental variogram is typically derived from well logging data, and then the theoretical variogram is fitted. As can be seen from the formula above, the variogram indicates that there is no correlation between sampling points in space when the distance exceeds a certain threshold; this distance is called the range.

[0172] Therefore, in this embodiment, only the correlation coefficient when the distance between two points is less than the range is retained; otherwise, it is set to zero, thereby removing spurious correlations and significantly reducing memory and computational consumption. The essence of this process is a low-rank approximation of the matrix. Assuming that any point in two-dimensional space is only correlated with the four points adjacent to it (up, down, left, and right), the correlation coefficient matrix Γ′ after the low-rank approximation is expressed as:

[0173]

[0174] Where n = nx × nz, and the element γ in the matrix i,j This indicates that the element is located on the diagonal of the i-th row from the main diagonal, and is the j-th element of that diagonal. Substituting the correlation coefficient matrix Γ′ after the low-rank approximation into formula (16), multiple geologically significant prior samples can be obtained, thereby accelerating the convergence speed of the full waveform inversion based on the SVGD algorithm and efficiently evaluating the uncertainty of the solution.

[0175] (3) Statistically determine the posterior distribution of the inversion results and evaluate its uncertainty.

[0176] Using the geologically significant initial sample set generated in the previous step, a full waveform inversion based on the SVGD algorithm is performed, as shown in equations (9)-(15). The inversion ends when a certain termination condition is met. The termination condition can be simply set as the number of iterations reaching the set maximum number of iterations, or it can be achieved by continuously statistically analyzing the posterior distribution of the inversion results. When the posterior distribution gradually stabilizes, the inversion ends.

[0177] Finally, the mean and standard deviation of all sample inversion results are calculated, and the posterior probability density distribution of the model feature locations is statistically analyzed to evaluate the uncertainty of the velocity inversion results.

[0178] Conventional full waveform inversion yields a definite inversion result, while the variational inference full waveform inversion in this patent yields multiple inversion solutions that satisfy certain conditions, and these solutions are all possible results.

[0179] The mean represents the average of multiple inversion solutions, while the standard deviation describes the discreteness of these solutions. All velocities corresponding to non-zero probabilities (frequency values) in the posterior probability density distribution are possible solutions.

[0180] The process of analyzing the mean, standard deviation, and posterior probability distribution is called uncertainty assessment. These uncertainty assessments provide experience for subsequent well location deployment and geological risk assessment in oil and gas exploration.

[0181] In summary, the variational inference full waveform inversion and its uncertainty evaluation method based on geostatistical prior information specifically includes the following steps.

[0182] Step 1) Apply the six filters shown in Equation (18) to the seismic migration profile to obtain six geological model score profiles. Then, calculate the multi-point model correlation between any two sampling points of the model based on the score vector, and thus form the correlation coefficient matrix Γ.

[0183] Step 2) Fit the theoretical variation function based on the well logging data and determine the range. Only retain the correlation coefficient when the distance between two points is less than the range to obtain the correlation coefficient matrix Γ′ after low-rank approximation.

[0184] Step 3) Given the number of samples and the maximum perturbation of the model, and according to Equation (16), a series of geologically significant initial samples are generated by perturbation.

[0185] Step 4) Based on each initial velocity sample, solve the acoustic wave equation (13) and its adjoint equation (14) using the finite difference method, and use the zero-delay cross-correlation between the source wave field and the adjoint wave field, i.e., equation (15), to calculate the gradient of the full waveform inversion objective function.

[0186] Step 5) Based on the gradient of the objective function, use Bayes' formula (1) and (11)-(12) to obtain the gradient of the posterior distribution.

[0187] Step 6) Given a small step size ε, update all samples using the SVGD algorithm as shown in equations (9)-(10), and return to step 4) until the number of iterations reaches the set maximum number of iterations, or the posterior distribution gradually stabilizes, the inversion ends, and the final inverted velocity model of all samples is output.

[0188] Step 7) Calculate the posterior distribution of all sample inversion results and evaluate the uncertainty.

[0189] This specific embodiment proposes a variational inference full waveform inversion and velocity uncertainty evaluation method based on geostatistical prior information. Compared with the traditional Bayesian full waveform inversion method, 1) this scheme proposes a geostatistical model correlation coefficient matrix with low-rank approximation, which removes spurious correlations and reduces memory and computational consumption; 2) multiplying the low-rank approximation model correlation coefficient matrix with a random vector can obtain geologically meaningful random initial samples, making the prior information more accurate; 3) based on geostatistical prior information, this scheme uses the SVGD algorithm for efficient variational inference full waveform inversion, which greatly reduces the ambiguity of the inversion and obtains an uncertainty evaluation of the solution.

[0190] The variational inference full waveform inversion method based on geological prior information described in the embodiments of this specification will be applied to a real-world scenario. Specifically, the following content will be included.

[0191] Step 1. Generate random samples based on prior information from geostatistics.

[0192] Figure 2 The Marmousi model is shown. Figure 2 (a) is the accurate model. Figure 2 (b) is a random model obtained by sampling from a uniform prior distribution. Figure 2 (c) represents the velocity prior distribution. To verify that the random samples based on geostatistical prior information proposed in this embodiment have geological significance, thereby reducing the inversion complexity, the following method is used: Figure 2 The Marmousi velocity model shown in (a) was tested. The velocity model has a grid number of 200×120 and a grid spacing of 20m. The maximum and minimum velocities are 4670m / s and 1500m / s, respectively. Figure 3 The reverse time-shifted profile of the Marmousi model is shown.

[0193] Figure 4 Six types of filters are shown. Figure 4 (a) is a vertical mean filter. Figure 4(b) is a water-average filter. Figure 4 (c) is the vertical gradient filter. Figure 4 (d) is the horizontal gradient filter. Figure 4 (e) is the vertical curvature filter. Figure 4 f is a horizontal curvature filter.

[0194] Applying the above six filters to the offset profile yields the following result: Figure 5 The six filter score profiles are shown. Figure 5 The filtered score profile obtained after filtering the offset profile is shown. Figure 5 (a) is the vertical mean filtering profile. Figure 5 (b) is the water-average filtered profile. Figure 5 (c) is the vertical gradient filter profile. Figure 5 (d) represents the horizontal gradient filter profile. Figure 5 (e) is the vertical curvature filtering profile. Figure 5 (f) is the horizontal curvature filter profile.

[0195] Figure 7 The multi-point pattern correlation matrix is ​​shown. Figure 7 (a) is the full-rank correlation coefficient matrix. Figure 7 (b) is the correlation coefficient matrix after low-rank approximation. Using equations (20) and (21), the multi-point model correlation coefficient matrix is ​​obtained, as follows: Figure 7 As shown in (a). The accurate Marmousi model velocity at x = 2 km was used as the statistical variogram function of the pseudo-well data, and then the theoretical variogram function was fitted to determine the range as 180 m, as shown in (a). Figure 6 As shown, the curve of the variation function is illustrated. Therefore, it is assumed that any point in the model is only correlated with the sampling points within a 180m × 180m radius around it, that is, only the correlation coefficient when the distance between two points is less than the range is retained. This yields the correlation coefficient matrix after the low-rank approximation, eliminating spurious correlations and reducing memory and computational consumption. Figure 7 As shown in (b).

[0196] Figure 8 This demonstrates the use of geostatistical priors to build stochastic models. Figure 8 (a) is a smooth background model. Figure 8 (b) represents a random perturbation based on geostatistical priors, with a maximum perturbation of 800 m / s. Figure 8 (c) is a stochastic model based on geostatistical priors.

[0197] Setting the maximum perturbation in equation (16) to 800 m / s, multiplying the correlation coefficient matrix after the low-rank approximation with the random vector yields the random perturbation based on geostatistical priors, such as... Figure 8 As shown in (b), compare it with Figure 8 The smooth velocity models shown in (a) are added together to obtain a stochastic model based on geostatistical priors, such as Figure 8 As shown in (c). With the use of Figure 2 The uniform prior distribution shown in (c) is obtained by sampling. Figure 2 Compared to the stochastic model shown in (b), the stochastic model generated based on geostatistical priors has geological significance, thereby reducing the search space for understanding and improving computational efficiency.

[0198] Step 2. Variational inference full waveform inversion and uncertainty evaluation.

[0199] Using the above-mentioned uniform prior information and geostatistical prior information, 100 random samples were generated as the initial model for variational inference FWI, and the influence of different prior information on the full waveform inversion results and velocity uncertainty evaluation were compared. Eight shot points were arranged at a depth of 360m with a lateral spacing of 500m, and the first shot was located at 200m. The geophones were set on the surface, with the middle excitation and the two sides receiving. The seismic record duration was 5s, the sampling interval ΔT=1ms, and the source was a Ricker wavelet. During the inversion, the data with the main frequency of 4Hz and 10Hz were filtered respectively, and two independent inversions were performed: (1) the data with the main frequency of 4Hz was inverted first, and (2) the low-frequency inversion results of all samples were used as their respective initial models to continue inverting the data with the main frequency of 10Hz.

[0200] The step size of the SVGD algorithm was set to the empirical step size for full waveform inversion. Inversions based on both priors were iterated 100 times each at 4Hz and 10Hz. Additionally, an inversion based on a uniform prior was iterated 600 times each at 4Hz and 10Hz. The mean and standard deviation of all sample inversion results were calculated as follows: Figures 9 to 11 As shown, histograms were plotted showing the edge distribution at four locations: (0.6km, 2km), (1.2km, 2km), (1.8km, 2km), and (2.4km, 2km). Figure 11 As shown, the black dashed lines represent the precise speeds at their respective positions.

[0201] Figure 9 For variational inference full waveform inversion results based on uniform prior information, each frequency band is iterated 100 times. The mean of the low-frequency inversion results (e.g.) Figure 9 As shown in (a), the actual model's structural form was not recovered, and low-frequency noise was present; while the mean of the high-frequency inversion results (as shown in (a)) was lower. Figure 9 As shown in (c), the structure of the velocity model can be recovered, but the inversion of the edge and deep regions is poor, and even falls into local extrema, which may be due to insufficient illumination. Figure 9(b) shows the standard deviation of the low-frequency inversion results. The standard deviation of the high-frequency inversion results (e.g.) Figure 9 As shown in (d), it exhibits structural features similar to the mean. Through Figure 11 (a) to Figure 11 As can be seen from the marginal distribution of the final inversion results at the four positions shown in (d), the above inversion did not converge, the inversion results between samples differed greatly, and most of them deviated from the accurate speed.

[0202] Therefore, the number of iterations for each frequency band was adjusted to 600, and the mean and standard deviation of the inversion results were statistically analyzed, such as... Figure 10 As shown. The mean of the low-frequency inversion results (e.g. Figure 10 As shown in (a), the construction of the velocity model has been recovered; while the mean of the high-frequency inversion results (as shown in (a)) has been recovered. Figure 10 (As shown in (c)) the resolution has improved, but the deep regions are still trapped in local extrema. The standard deviation of the high-frequency inversion results (e.g.) Figure 10 (d) shows the standard deviation of the lower frequency inversion results (as shown in the figure). Figure 10 As shown in (b), the value decreased. Furthermore, the standard deviation after 600 iterations (as shown in [reference]) was also lower. Figure 10 (d) shows that the overall standard deviation is slightly lower than the standard deviation of the statistics after 100 iterations (as shown in the figure). Figure 9 As shown in (d), this illustrates that as the number of iterations increases, the uncertainty of the inversion decreases and the inversion accuracy improves. Through... Figure 12 (e) to Figure 12 As shown in (h), the marginal distribution of the final inversion results at the four locations shows that the inversion results of all samples at the (0.6km, 2km) location converge to the accurate velocity, while the inversion results at the other three locations do not converge. This indicates that when using uniform priors for variational inference full waveform inversion, the search space for the solution is large, and the inversion results of most samples differ significantly from the accurate model.

[0203] Figure 9 use Figure 2 100 initial samples were generated from the prior information in (c), and full waveform inversion based on the SVGD algorithm was performed. After 100 iterations at a main frequency of 4Hz, the mean (a) and standard deviation (b) were calculated. The inversion results at 4Hz were used as initial samples to continue the inversion at 10Hz. After 100 iterations, the mean (c) and standard deviation (d) were calculated.

[0204] Figure 10 use Figure 2 100 initial samples were generated from the prior information of (c), and full waveform inversion based on the SVGD algorithm was performed. After 600 iterations at a main frequency of 4Hz, the mean of (a) and the standard deviation of (b) were calculated. The inversion result of 4Hz was used as the initial sample to continue the inversion at 10Hz. After 600 iterations, the mean was calculated (e.g. Figure 10 (as shown in (c)) and standard deviation (as shown in...) Figure 10 (as shown in (d)).

[0205] Figure 11 100 initial samples were generated using prior information from geostatistics. Full waveform inversion based on the SVGD algorithm was then performed. After 100 iterations at a dominant frequency of 4Hz, the mean (a) and standard deviation (b) were calculated. The 4Hz inversion results were used as initial samples for further inversion at 10Hz. After 100 iterations, the mean (e.g.) was calculated. Figure 11 (as shown in (c)) and standard deviation (as shown in...) Figure 11 (as shown in (d)).

[0206] The full waveform inversion result based on variational inference using prior information from geostatistics is as follows: Figure 12 As shown. The mean of the low-frequency inversion results (e.g. Figure 12 As shown in (a), the structure of the velocity model can be recovered relatively well; while the mean of the high-frequency inversion results (as shown in (a)) can be recovered relatively well. Figure 12 (c) shows that the resolution is further improved, and the inversion results for deep and edge regions are far better than the inversion mean based on a uniform prior distribution (e.g. Figure 9 (c) Figure 10 (as shown in (c)). Standard deviation of high-frequency inversion results (e.g.) Figure 12 (d) shows the standard deviation of the lower frequency inversion results (as shown in the figure). Figure 10 (b) shows a decrease. And it is much lower than the inversion mean based on a uniform prior distribution (as shown in [reference]). Figure 9 (d) Figure 10 (as shown in (d)). Figure 12 (i) to Figure 12 As shown in (l), the marginal distribution of the final inversion results at the four locations indicates that the inversion results of all samples at the four locations are only slightly different from the accurate velocity, or converge well to the accurate velocity. This suggests that in variational inference full waveform inversion, introducing geostatistical prior information can reduce the search space of the solution, and a small number of iterations are sufficient to achieve good inversion results for all samples.

[0207] Figure 12 The results of variational inference full waveform inversion based on prior information from geostatistics are shown. Specifically, Figure 12 (a) through (d) show the calculations Figure 9 The edge distribution of the final inversion result of sample 10Hz at different locations. Figure 12 The location corresponding to (a) is (0.6km, 2km). Figure 12 The location corresponding to (b) is (1.2km, 2km). Figure 12 The location corresponding to (c) is (1.8km, 2km). Figure 12The location corresponding to (d) is (2.4km, 2km). (e) to (h) of 12 show the calculation... Figure 10 The edge distribution of the final inversion result of sample 10Hz at different locations. Figure 12 The location corresponding to (e) is (0.6km, 2km). Figure 12 The location corresponding to (f) is (1.2km, 2km). Figure 12 The location corresponding to (g) is (1.8km, 2km). Figure 12 The location corresponding to (h) is (2.4km, 2km). Figure 12 (i) to (l) show the calculation Figure 11 Comparison of edge distribution at different locations in the final inversion result of sample 10Hz. Figure 12 The location corresponding to (i) is (0.6km, 2km). Figure 12 The location corresponding to (j) is (1.2km, 2km). Figure 12 The location corresponding to (k) is (1.8km, 2km). Figure 12 The location (l) corresponds to the edge distribution comparison map at (2.4km, 2km). The speed shown by the black dashed line is the accurate speed at the corresponding location.

[0208] In this specific embodiment, the inversion time was tested using eight NVIDIA Tesla P40 graphics cards for different iteration numbers. The computation time for 100 iterations with 100 samples was 14 hours, and the computation time for 600 iterations was 84 hours. Therefore, from the perspective of speed uncertainty evaluation, compared to uniform prior information, the initial samples established based on geostatistical priors have geological significance, narrow the search space for solutions, reduce the difficulty of inversion, greatly improve computational efficiency, and also improve the accuracy of inversion while reducing the ambiguity of inversion solutions.

[0209] Based on the same inventive concept, this specification also provides a variational inference full waveform inversion device based on geological prior information, as described in the following embodiments. Since the principle of the variational inference full waveform inversion device based on geological prior information is similar to that of the variational inference full waveform inversion method based on geological prior information, the implementation of the variational inference full waveform inversion device based on geological prior information can refer to the implementation of the variational inference full waveform inversion method based on geological prior information, and repeated details will not be elaborated further. As used below, the terms "unit" or "module" can refer to a combination of software and / or hardware that performs a predetermined function. Although the device described in the following embodiments is preferably implemented in software, hardware implementation, or a combination of software and hardware, is also possible and contemplated. Figure 13 This is a structural block diagram of a variational inference full waveform inversion device based on geological prior information, as described in this specification. Figure 13 As shown, it includes: calculation module 1301, approximation module 1302, generation module 1303 and inversion module 1304. The structure is described below.

[0210] The calculation module 1301 is used to construct various filters; it is also used to apply the constructed filters to the seismic migration profile to calculate the correlation coefficient matrix.

[0211] The approximation module 1302 is used to perform a low-rank approximation on the calculated correlation coefficient matrix based on the fitting theoretical variation function of the well logging data, so as to obtain the low-rank approximated correlation coefficient matrix.

[0212] The generation module 1303 is used to randomly generate multiple velocity samples based on the correlation coefficient matrix after the low-rank approximation;

[0213] The inversion module 1304 is used to calculate the gradient of the full waveform inversion objective function based on the plurality of velocity samples; it is also used to calculate the gradient of the full waveform inversion posterior distribution using Bayes' theorem based on the gradient of the objective function; and it is further used to perform multiple iterative updates on the plurality of velocity samples using the Stein variational gradient descent algorithm based on the gradient of the full waveform inversion posterior distribution, outputting the target velocity model obtained by inverting the plurality of velocity samples. In some embodiments of this specification,

[0214] In some embodiments of this specification, the apparatus further includes a determination module, which is used to: determine the posterior distribution of the target velocity model in order to evaluate the uncertainty of the inversion results.

[0215] In some embodiments of this specification, the calculation module is specifically used to: set six filters for a two-dimensional offset profile; the six filters include: a vertical mean filter, a horizontal mean filter, a vertical gradient filter, a horizontal gradient filter, a vertical curvature filter, and a horizontal curvature filter.

[0216] In some embodiments of this specification, the calculation module is specifically used to: convert the geological patterns captured from the seismic migration profile into score values ​​through the multiple filters; define multiple score values ​​corresponding to the geological patterns around each node in the seismic migration profile as score vectors corresponding to each node; calculate the similarity of the geological patterns around any two nodes based on the score vectors corresponding to any two nodes in the seismic migration profile to obtain the multi-point pattern correlation between any two nodes; and generate a correlation coefficient matrix based on the multi-point pattern correlation between any two nodes.

[0217] In some embodiments of this specification, the Stein variational gradient descent algorithm is used to iteratively update the plurality of velocity samples, including: given a step size and updating the plurality of velocity samples using the Stein variational gradient descent algorithm; calculating the gradient of the full waveform inversion objective function based on the plurality of velocity samples; calculating the gradient of the full waveform inversion posterior distribution using Bayes' theorem based on the gradient of the objective function; repeating the above steps until the number of iterations reaches the set maximum number of iterations or until the full waveform inversion posterior distribution tends to stabilize.

[0218] In some embodiments of this specification, the kernel function of the Stein variational gradient descent algorithm is a radial basis function.

[0219] As can be seen from the above description, the embodiments of this specification achieve the following technical effects: a correlation coefficient matrix of geostatistical models with low-rank approximation is proposed, which removes spurious correlations and reduces memory and computational consumption. Based on the low-rank approximation correlation coefficient matrix, random initial samples with geological significance are randomly generated, making the prior information more accurate. Then, based on the geostatistical prior information, the Stein variational gradient descent algorithm is used to perform efficient variational inference full waveform inversion, which greatly reduces the ambiguity of the inversion and improves the inversion efficiency.

[0220] This specification also provides a computer device, which can be found in the following description. Figure 14 The diagram shown illustrates the computer device structure for the variational inference full waveform inversion method based on geological prior information provided in the embodiments of this specification. Specifically, the computer device may include an input device 141, a processor 142, and a memory 143. The memory 143 stores processor-executable instructions. When the processor 142 executes the instructions, it implements the steps of the variational inference full waveform inversion method based on geological prior information described in any of the above embodiments.

[0221] In this embodiment, the input device can specifically be one of the main devices for information exchange between the user and the computer system. The input device may include a keyboard, mouse, camera, scanner, light pen, handwriting input tablet, voice input device, etc.; the input device is used to input raw data and programs for processing these data into the computer. The input device can also receive data transmitted from other modules, units, and devices. The processor can be implemented in any suitable manner. For example, the processor can take the form of a microprocessor or processor and a computer-readable medium storing computer-readable program code (e.g., software or firmware) executable by the (micro)processor, logic gates, switches, application-specific integrated circuits (ASICs), programmable logic controllers, and embedded microcontrollers, etc. The memory can specifically be a memory device used to store information in modern information technology. The memory can include multiple layers; in digital systems, anything that can store binary data can be considered memory; in integrated circuits, a circuit without physical form but with storage function is also called memory, such as RAM, FIFO, etc.; in a system, a storage device with physical form is also called memory, such as a memory stick, TF card, etc.

[0222] In this embodiment, the specific functions and effects implemented by the computer device can be explained in comparison with other embodiments, and will not be repeated here.

[0223] This specification also provides a computer storage medium for a variational inference full waveform inversion method based on geological prior information, wherein the computer storage medium stores computer program instructions that, when executed, implement the steps of the variational inference full waveform inversion method based on geological prior information described in any of the above embodiments.

[0224] In this embodiment, the storage medium includes, but is not limited to, Random Access Memory (RAM), Read-Only Memory (ROM), cache, hard disk drive (HDD), or memory card. The memory can be used to store computer program instructions. The network communication unit can be an interface configured according to standards specified in the communication protocol for network connection communication.

[0225] In this embodiment, the specific functions and effects implemented by the program instructions stored in the computer storage medium can be explained by comparison with other embodiments, and will not be repeated here.

[0226] Obviously, those skilled in the art will understand that the modules or steps of the embodiments described above can be implemented using general-purpose computing devices. They can be centralized on a single computing device or distributed across a network of multiple computing devices. Optionally, they can be implemented using computer-executable program code, thereby storing them in a storage device for execution by a computing device. In some cases, the steps shown or described can be performed in a different order than those presented herein, or they can be fabricated as separate integrated circuit modules, or multiple modules or steps can be fabricated as a single integrated circuit module. Thus, the embodiments of this specification are not limited to any particular combination of hardware and software.

[0227] It should be understood that the above description is for illustrative purposes and not for limitation. Many embodiments and applications beyond the provided examples will be apparent to those skilled in the art upon reading the above description. Therefore, the scope of this specification should not be determined by reference to the above description, but rather by reference to the foregoing claims and the full scope of their equivalents.

[0228] The above description is merely a preferred embodiment of this specification and is not intended to limit this specification. Various modifications and variations can be made to the embodiments described herein by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this specification should be included within the scope of protection of this specification.

Claims

1. A variational inference full-waveform inversion method based on geological prior information, characterized in that, include: Construct multiple filters; The constructed filters are applied to the seismic migration profile to calculate the correlation coefficient matrix; Based on the fitting theoretical variogram function of well logging data, the calculated correlation coefficient matrix is ​​approximated by a low-rank approximation. Only the correlation coefficients between two points with a distance less than the range are retained, and the correlation coefficients with a distance not less than the range are set to zero, thus obtaining the correlation coefficient matrix after low-rank approximation. Based on the correlation coefficient matrix after the low-rank approximation, multiple velocity samples are randomly generated. Based on the multiple velocity samples, the gradient of the full waveform inversion objective function is calculated; based on the gradient of the objective function, the gradient of the full waveform inversion posterior distribution is calculated using Bayes' theorem; based on the gradient of the full waveform inversion posterior distribution, the Stein variational gradient descent algorithm is used to iteratively update the multiple velocity samples multiple times, and the target velocity model obtained by inverting the multiple velocity samples is output. The process involves applying various filters to a seismic migration profile to calculate a correlation coefficient matrix. This includes: converting geological patterns captured from the seismic migration profile into score values ​​using the filters; defining multiple score values ​​corresponding to the geological patterns around each node in the seismic migration profile as score vectors for each node; calculating the similarity of geological patterns around any two nodes based on their score vectors to obtain multi-point pattern correlation between them; and generating a correlation coefficient matrix based on the multi-point pattern correlation between any two nodes.

2. The variational inference full waveform inversion method based on geological prior information according to claim 1, characterized in that, After iteratively updating the multiple velocity samples using the Stein variational gradient descent algorithm based on the gradient of the full waveform inversion posterior distribution, and outputting the target velocity model obtained from the inversion of the multiple velocity samples, the method further includes: The posterior distribution of the target velocity model is determined to evaluate the uncertainty of the inversion results.

3. The variational inference full-waveform inversion method based on geological prior information according to claim 1, characterized in that, Construct various filters, including: For the two-dimensional offset profile, six filters are set: vertical mean filter, horizontal mean filter, vertical gradient filter, horizontal gradient filter, vertical curvature filter, and horizontal curvature filter.

4. The variational inference full-waveform inversion method based on geological prior information according to claim 1, characterized in that, The Stein variational gradient descent algorithm is used to iteratively update the multiple velocity samples, including: Given a step size, update multiple velocity samples using the Stein variational gradient descent algorithm; Based on the multiple velocity samples, calculate the gradient of the full waveform inversion objective function; based on the gradient of the objective function, calculate the gradient of the full waveform inversion posterior distribution using Bayes' theorem; Repeat the above steps until the number of iterations reaches the set maximum number of iterations or until the posterior distribution of the full waveform inversion tends to stabilize.

5. The variational inference full-waveform inversion method based on geological prior information according to claim 1, characterized in that, The kernel function of the Stein variational gradient descent algorithm is a radial basis function.

6. A variational inference full-waveform inversion device based on geological prior information, characterized in that, include: The calculation module is used to construct various filters; it is also used to apply the constructed filters to the seismic migration profile to calculate the correlation coefficient matrix. The approximation module is used to perform a low-rank approximation on the calculated correlation coefficient matrix based on the fitting theoretical variation function of the well logging data. It retains only the correlation coefficients where the distance between two points is less than the range, and sets the correlation coefficients where the distance between two points is not less than the range to zero, thus obtaining the low-rank approximation correlation coefficient matrix. The generation module is used to randomly generate multiple velocity samples based on the correlation coefficient matrix after the low-rank approximation. The inversion module is used to calculate the gradient of the full waveform inversion objective function based on the multiple velocity samples; It is also used to calculate the gradient of the full waveform inversion posterior distribution using Bayes' theorem based on the gradient of the objective function; and to perform multiple iterations of updating the multiple velocity samples using the Stein variational gradient descent algorithm based on the gradient of the full waveform inversion posterior distribution, and output the target velocity model obtained by inverting the multiple velocity samples. Specifically, the calculation module is used to: convert the geological patterns captured from the seismic migration profile into score values ​​through the multiple filters; define multiple score values ​​corresponding to the geological patterns around each node in the seismic migration profile as the score vector corresponding to each node; and calculate the similarity of the geological patterns around any two nodes based on the score vectors corresponding to any two nodes in the seismic migration profile to obtain the multi-point pattern correlation between any two nodes. A correlation coefficient matrix is ​​generated based on the multi-point pattern correlation between any two nodes.

7. A computer device, characterized in that, It includes a processor and a memory for storing processor-executable instructions, wherein the processor, when executing the instructions, implements the steps of the method according to any one of claims 1 to 5.

8. A computer-readable storage medium storing computer instructions thereon, characterized in that, When the instructions are executed by the processor, they implement the steps of the method according to any one of claims 1 to 5.

Images