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

Abstract

The invention relates to the technical field of seismic inversion, and particularly discloses a variational reasoning full-waveform inversion method and device based on geological prior information, and the method comprises the steps: constructing a plurality of filters; acting the plurality of filters obtained by construction on a seismic migration section to calculate a correlation coefficient matrix; according to a logging information fitting theory variation function, performing low-rank approximation on the calculated correlation coefficient matrix to obtain a correlation coefficient matrix after low-rank approximation; randomly generating a plurality of speed samples based on the correlation coefficient matrix after low-rank approximation; calculating the gradient of the objective function according to the plurality of speed samples; according to the gradient of the objective function, calculating the gradient of posteriori distribution by using a Bayesian formula; and based on the gradient of the full waveform inversion posteriori distribution, utilizing a Stein variational gradient descent algorithm to carry out multiple times of iterative updating on the plurality of speed samples, and outputting a target speed model obtained by inversion of the plurality of speed samples. According to the scheme, the multiplicity of inversion is reduced, and the inversion efficiency is improved.

Description

Technical Field

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

[0002] High-precision seismic imaging requires a high-precision depth-domain velocity model. Full waveform inversion (FWI) comprehensively utilizes information such as seismic wave travel time, amplitude, and phase, and can effectively make up for the missing medium and low frequency components in conventional velocity modeling, playing a key bridging role between travel time tomography and migration imaging. Thanks to the significant improvement in computer performance, FWI has received much attention in recent years. However, limited by the observation system, data bandwidth, approximation degree of the forward operator, noise interference, and the strong non-linearity of the wave equation, FWI is actually a strongly non-linear inverse problem under severely insufficient information conditions, with serious multi-solution problems. Currently, the vast majority of FWI methods adopt local optimization algorithms, and their inversion effects strongly depend on the initial model. When the initial model is inaccurate, FWI will fall into local minima. For the multi-solution inverse problem with insufficient information, compared with simply improving the algorithm to obtain the "optimal solution" in the sense of data matching, it is more reasonable to solve the inverse problem in the Bayesian framework, obtain the posterior probability density distribution of the model, and evaluate the uncertainty of the FWI solution.

[0003] FWI under the Bayesian framework has at least two advantages compared with traditional deterministic FWI. One is that it can conveniently use various prior information (such as geological understanding, etc.) to make up for the insufficient seismic data information and narrow the FWI solution space; the other is to provide an uncertainty evaluation of the velocity model, which is beneficial to assessing exploration risks and subsequent geological interpretation. With the significant improvement in computer computing performance, the research on FWI methods under the Bayesian framework has gradually developed. The Bayesian FWI method can make full use of various prior information to narrow the solution space and provide an uncertainty evaluation of the solution, which is an effective means to solve the multi-solution problem of FWI. However, the current Bayesian inversion has a relatively high multi-solution degree and low computational efficiency.

[0004] In view of the above problems, no effective solution has been proposed yet. Summary of the Invention

[0005] The embodiments of this specification provide a variational inference full waveform inversion method and device based on geological prior information to improve the computational efficiency of full waveform inversion and reduce the multi-solution degree of inversion.

[0006] The embodiments of this specification provide a variational inference full waveform inversion method based on geological prior information, including:

[0007] Constructing a variety of filters; applying the constructed variety of filters to the seismic migration profile to calculate the correlation coefficient matrix;

[0008] Perform low-rank approximation on the calculated correlation coefficient matrix according to the theoretical variogram fitted from well logging data to obtain the correlation coefficient matrix after low-rank approximation;

[0009] Based on the correlation coefficient matrix after low-rank approximation, randomly generate multiple velocity samples;

[0010] Calculate the gradient of the full waveform inversion objective function according to the multiple velocity samples; calculate the gradient of the posterior distribution of full waveform inversion using Bayes' formula according to the gradient of the objective function; based on the gradient of the posterior distribution of full waveform inversion, use the Stein variational gradient descent algorithm to perform multiple iterative updates on the multiple velocity samples, and output the target velocity model obtained by inverting the multiple velocity samples.

[0011] In one embodiment, after using the Stein variational gradient descent algorithm to perform multiple iterative updates on the multiple velocity samples based on the gradient of the posterior distribution of full waveform inversion and outputting the target velocity model obtained by inverting the multiple velocity samples, it further includes:

[0012] Determine the posterior distribution of the target velocity model to evaluate the uncertainty of the inversion result.

[0013] In one embodiment, construct multiple filters, including:

[0014] For a two-dimensional migration profile, set six filters; the six filters include: vertical mean filter, horizontal mean filter, vertical gradient filter, horizontal gradient filter, vertical curvature filter, and horizontal curvature filter.

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

[0016] Through the multiple filters, convert the geological patterns captured from the seismic migration profile into score values;

[0017] Define the multiple score values corresponding to the geological patterns around each node in the seismic migration profile as the score vector corresponding to each node;

[0018] According to the score vectors corresponding to any two nodes in the seismic migration profile, calculate the similarity of the geological patterns around the any two nodes to obtain the multi-point pattern correlation between the any two nodes;

[0019] Generate a correlation coefficient matrix based on the multi-point pattern correlation between any two nodes.

[0020] In one embodiment, using the Stein variational gradient descent algorithm to perform multiple iterative updates on the multiple velocity samples includes:

[0021] Update multiple velocity samples with a given step size and using the Stein variational gradient descent algorithm;

[0022] 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 the Bayesian formula based on the gradient of the objective function;

[0023] Repeat 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 be stable.

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

[0025] The embodiments of this specification also provide a variational inference full waveform inversion device based on geological prior information, including:

[0026] A calculation module, configured to construct multiple filters; and further configured to apply the multiple constructed filters to a seismic migration profile to calculate a correlation coefficient matrix;

[0027] An approximation module, configured to perform low-rank approximation on the calculated correlation coefficient matrix according to a theoretical variogram fitted from well logging data to obtain a low-rank approximated correlation coefficient matrix;

[0028] A generation module, configured to randomly generate multiple velocity samples based on the low-rank approximated correlation coefficient matrix;

[0029] An inversion module, configured to calculate the gradient of the full waveform inversion objective function according to the multiple velocity samples; further configured to calculate the gradient of the full waveform inversion posterior distribution using the Bayesian formula based on the gradient of the objective function; further configured to perform multiple iterative updates on 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 configured 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; calculate the similarity of the geological patterns around any two nodes according to the score vectors corresponding to any two nodes in the seismic migration profile to obtain the multi-point pattern correlation between any two nodes; generate a correlation coefficient matrix based on the multi-point pattern correlation between any two nodes.

[0031] An embodiment of this specification also provides a computer device, including a processor and a memory for storing instructions executable by the processor. When the processor executes the instructions, the steps of the variational inference full waveform inversion method based on geological prior information described in any of the above embodiments are implemented.

[0032] An embodiment of this specification also provides a computer-readable storage medium, on which computer instructions are stored. When the instructions are executed, the steps of the variational inference full waveform inversion method based on geological prior information described in any of the above embodiments are implemented.

[0033] In an embodiment of this specification, a variational inference full waveform inversion method based on geological prior information is provided. Multiple filters can be constructed, and the multiple constructed filters are applied to a seismic migration profile to calculate a correlation coefficient matrix. The calculated correlation coefficient matrix is low-rank approximated according to a theoretical variogram fitted from well logging data, and a low-rank approximated correlation coefficient matrix is obtained. Based on the low-rank approximated correlation coefficient matrix, multiple velocity samples are randomly generated. According to the multiple velocity samples, the gradient of the full waveform inversion objective function is calculated. According to the gradient of the objective function, the gradient of the full waveform inversion posterior distribution is calculated using Bayes' formula. Based on the gradient of the full waveform inversion posterior distribution, the multiple velocity samples are iteratively updated multiple times using the Stein variational gradient descent algorithm, and the target velocity model inverted from the multiple velocity samples is output. In the above solution, a correlation coefficient matrix with a geostatistical pattern of low-rank approximation is proposed, which removes false correlations and reduces memory and computational consumption. Random initial samples with geological significance are randomly generated based on the low-rank approximated correlation coefficient matrix, making the prior information more accurate. Then, based on the geostatistical prior information, the Stein variational gradient descent algorithm is used for efficient variational inference full waveform inversion, greatly reducing the non-uniqueness of the inversion and improving the inversion efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] The drawings described herein are used to provide a further understanding of this specification, form a part of this specification, and do not limit this specification. In the drawings:

[0035] Figure 1 Shows a flowchart of the variational inference full waveform inversion method based on geological prior information in an embodiment of this specification;

[0036] Figure 2 Shows the Marmousi model in an embodiment of this specification;

[0037] Figure 3 Shows the reverse time migration profile of the Marmousi model in an embodiment of this specification;

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

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

[0040] Figure 6 A variogram in an embodiment of this specification is shown;

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

[0042] Figure 8 Building a stochastic model using geostatistical priors in an embodiment of this specification is shown;

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

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

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

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

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

[0048] Figure 14 A schematic diagram of a computer device in an embodiment of this specification is shown. Detailed implementation manners

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

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

[0051] Bayesian inversion uses the prior probability density distribution of the model to estimate the posterior probability density distribution of the model, and is an effective means to evaluate the uncertainty of the FWI solution. Bayesian inversion methods mainly include sampling methods and variational inference methods. The sampling method is represented by the Markov-Monte Carlo method. The disadvantage of this type of method is the low sample acceptance rate and large computational amount. Another method to solve the Bayesian problem is the variational inference method. By minimizing a certain objective functional, the prior probability density distribution is continuously updated to make it sufficiently close to the posterior probability density distribution to be solved. The variational inference method can use the gradient information of the objective function, so the computational efficiency is much higher than that of the sampling method.

[0052] Among them, the Stein variational gradient descent method (SVGD) is relatively commonly used, which can be used to estimate complex multi-modal distribution problems, is easy to parallelize, and has high computational efficiency. This algorithm generates a series of samples according to the prior distribution, and correlates the gradient information of different samples through a kernel function, which can prevent the current sample from falling into local extrema. The samples after multiple iterations of update 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 result.

[0053] In Bayesian inversion, prior information is crucial for the posterior inversion result. Currently, extremely broad prior information is adopted in Bayesian full waveform inversion, that is, a uniform distribution with the same variance is satisfied at each sampling point, which will greatly increase the computational amount. Therefore, the method in the embodiments of this specification incorporates geostatistical prior information into Bayesian FWI. By applying the filter in multi-point geostatistics to the seismic migration profile, the correlation coefficient matrix of the geological patterns around any two points can be obtained. And a low-rank approximation is made to the correlation coefficient matrix, only retaining the correlation coefficients when the distance between two points is less than the range, thereby removing false correlations and greatly reducing the memory and computational consumption. Multiplying the above matrix by a random vector can obtain random samples with geological significance.

[0054] The method in the embodiments of this specification uses geostatistical prior information to generate multiple random samples with geological significance, and uses the SVGD algorithm for efficient variational inference full waveform inversion, greatly reducing the solution search space, reducing the non-uniqueness of the inversion, and can also determine the uncertainty estimation of the solution.

[0055] Based on this, an embodiment of this specification provides a variational inference full waveform inversion method based on geological prior information. Figure 1 The flowchart of the variational inference full waveform inversion method based on geological prior information in an embodiment of this specification is shown. Although this specification provides method operation steps or device structures as shown in the following embodiments or drawings, more or fewer operation steps or module units may be included in the method or device based on routine or non-creative labor. In steps or structures where there is no necessary causal relationship logically, the execution order of these steps or the module structure of the device is not limited to the execution order or module structure described in the embodiments of this specification and shown in the drawings. When the method or module structure is applied to an actual device or terminal product, it can be executed sequentially or in parallel according to the method or module structure shown in the embodiments or drawings (for example, in an environment of parallel processors or multi-threaded processing, or even a distributed processing environment).

[0056] Specifically, as Figure 1 shown, a variational inference full waveform inversion method based on geological prior information provided by an 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 a computer device. Multiple filters can be constructed. The filters are used to extract key information related to geological patterns from the seismic migration profile. The key information extracted by different types of filters is different. By constructing multiple filters, multiple key information related to geological patterns can be extracted from the seismic migration profile.

[0059] When applying the constructed multiple filters to the seismic migration profile, it is actually using these filters to analyze and process the geological patterns in the profile. Through the operations of the filters, key information related to geological patterns can be captured from the seismic profile, and then based on this key information, the correlation coefficient matrix can be calculated. The correlation coefficient matrix is used to characterize the degree of geological pattern association between different nodes in the seismic migration profile, providing important basic information for subsequent full waveform inversion.

[0060] Step S102, perform low-rank approximation on the calculated correlation coefficient matrix according to the well logging data fitting theoretical variogram to obtain the low-rank approximated correlation coefficient matrix.

[0061] Well logging data contains rich geological information. By analyzing and processing well logging data, a theoretical variogram can be fitted. The theoretical variogram is a function used to describe the variation characteristics of a variable in space. In geostatistics, it represents the semi-variance of the increment 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 geological pattern association between different nodes in the seismic migration profile, but there may be some problems with the original correlation coefficient matrix. Since it indicates 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, consuming a large amount of computing and storage resources. Using the theoretical variogram to perform low-rank approximation on the correlation coefficient matrix is to process the correlation coefficient matrix according to 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 usually considered not to be correlated) are retained, otherwise the correlation coefficients are set to zero. The essence of this is to simplify the matrix and remove those unreasonable and false correlations. It can make the correlation coefficient matrix after low-rank approximation more accurately reflect the true correlation of the subsurface geological pattern, so that the entire inversion process is based on a more realistic basis.

[0062] After low-rank approximation, the correlation coefficients between a large number of points with distances greater than the range in the matrix are set to zero, and the scale of the matrix is significantly reduced, thus greatly reducing the memory occupancy and the amount of computation during the calculation process, improving the computational efficiency, making full waveform inversion more feasible in practical applications, especially when dealing with large-scale data, the advantages are more obvious.

[0063] At the same time, a more accurate correlation coefficient matrix provides a more reliable basis for subsequent inversion steps. The velocity model generated by perturbing the correlation coefficient matrix after low-rank approximation is more in line with geological reality. In the processes of calculating the gradient of the full waveform inversion objective function and calculating the gradient of the posterior distribution using Bayes' formula, it can more accurately reflect the relationship between the velocity model and the observed data, helping to guide the inversion towards a direction closer to the true velocity model, thereby improving the accuracy of the inversion results and reducing the non-uniqueness of the inversion.

[0064] After low-rank approximation, the resulting correlation coefficient matrix after low-rank approximation is more in line with the actual situation of underground geology, while greatly reducing memory and computational consumption. This correlation coefficient matrix after low-rank approximation will be used in subsequent steps, such as randomly generating multiple velocity samples, etc., to provide more accurate and efficient basic data support for full waveform inversion. When using the Stein variational gradient descent algorithm to iteratively update the velocity samples, the samples generated based on a more reasonable correlation coefficient matrix can converge to a stable posterior distribution faster, reducing unnecessary iteration times, accelerating the convergence speed of the inversion, and saving computational time.

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

[0066] After obtaining the correlation coefficient matrix after low-rank approximation, multiple velocity samples can be randomly generated in a specific way. Generally, some random factors will be combined to achieve the perturbation process.

[0067] In the field of stochastic geological modeling, the following method can be used to perturb and generate random samples:

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

[0069] where m represents the column vector formed by arranging the two-dimensional perturbed random velocity model in columns, 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 form the first column of the two-dimensional velocity model, and so on. nx and nz are the number of columns and rows of the two-dimensional model respectively, m smooth is the vector composed of the smooth background model, R is a random vector, following a uniform distribution in the interval (0, R max ), R max is the maximum perturbation amount, and Γ′ is the correlation coefficient matrix after low-rank approximation.

[0070] By the above method, the correlation coefficient matrix after low-rank approximation is combined with the random vector, so that the generated velocity samples contain both the information of geological pattern correlation (reflected by the correlation coefficient matrix) and introduce a certain degree of randomness (brought by the random vector). These multiple velocity samples have geological significance, and they can more realistically simulate various possible situations of the underground velocity model, providing rich and reasonable initial conditions for the subsequent full waveform inversion. Subsequently, based on these generated velocity samples, operations such as further calculating the gradient of the full waveform inversion objective function will be carried out to promote the inversion process to 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' formula according to the gradient of the objective function; based on the gradient of the full waveform inversion posterior distribution, use the Stein variational gradient descent algorithm to perform multiple iterative updates on the multiple velocity samples, and output the target velocity model obtained by inverting the multiple velocity samples.

[0072] Specifically, the gradient of the full waveform inversion objective function can be calculated based on the multiple velocity samples. The goal of full waveform inversion is to find an optimal velocity model such that the synthetic seismic data forward modeled based on this model matches the actual observed data as much as possible. The objective function is used to measure this degree of matching and can usually be expressed as the error between the synthetic data and the observed data, such as in the form of the mean square error between 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 each of the multiple velocity samples, the variation of the objective function with respect to the velocity model can be determined.

[0073] In one embodiment, based on each velocity sample, the acoustic wave equation and its adjoint equation can be solved using the finite difference method, and the gradient of the full waveform inversion objective function can be calculated using the zero-delay cross-correlation of the source wavefield and the adjoint wavefield. The acoustic wave equation describes the propagation process of seismic waves in a medium, which includes physical quantities such as the source wavelet, density, seismic wave propagation velocity, pressure field, horizontal and vertical components of particle vibration velocity, etc. The finite difference method is a numerical calculation method that discretizes continuous space and time and approximates the solution of the acoustic wave equation by calculating the physical quantities at discrete points. Through this method, the propagation of seismic waves under a given initial velocity sample can be simulated to obtain the source wavefield. The adjoint equation is another equation related to the acoustic wave equation and plays a key role in the process of calculating the gradient of the objective function. Similarly, the adjoint equation is solved using the finite difference method to obtain the adjoint wavefield. The calculation of the adjoint wavefield is in the opposite direction of the propagation of the source wavefield, calculating from the maximum time to zero time.

[0074] The zero-delay cross-correlation of the source wavefield and the adjoint wavefield means that the gradient is the zero-delay cross-correlation of the source wavefield propagating in the positive time direction and the adjoint wavefield propagating in the negative time direction. Through this calculation method, the gradient of the full waveform inversion objective function with respect to the velocity model can be obtained, providing an important basis for subsequent updating of the velocity model using this gradient to gradually approach the optimal velocity model.

[0075] The gradient of the full waveform inversion posterior distribution can be calculated using Bayes' formula according to the gradient of the objective function. 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 ocbs), this process can be described by the following formula:

[0076]

[0077] where d obs is the observed data, p(d obs ) is the normalization constant, and p(d obs |m) is the likelihood function.

[0078] Combined with the gradient of the objective function, the Bayesian formula 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 the prior information (such as the prior probability density distribution obtained by the geostatistical method previously) and the observed data information, by calculating the gradient of the posterior distribution, a more comprehensive understanding of the uncertainty and possibility of the velocity model can be obtained, providing more accurate guidance for subsequent iterative updates.

[0079] Based on the gradient of the full waveform inversion posterior distribution, the Stein variational gradient descent algorithm can be used to perform multiple iterative updates on the multiple velocity samples, and the target velocity model inverted from the multiple velocity samples is output. The Stein variational gradient descent algorithm is an algorithm for approximating complex posterior distributions, with the advantages of being able to estimate complex multimodal distribution problems, being easy to parallelize, and having high computational efficiency. By performing iterative updates on multiple velocity samples, these samples gradually approach the posterior probability density distribution.

[0080] In the above embodiments, a correlation coefficient matrix of a geostatistical pattern with low-rank approximation is proposed, which removes spurious correlations and reduces memory and computational consumption. Random initial samples with geological significance are generated based on the low-rank approximation correlation coefficient matrix to make the prior information more accurate. Then, based on the geostatistical prior information, the Stein variational gradient descent algorithm is used for efficient variational inference full waveform inversion, greatly reducing the non-uniqueness of the inversion and improving the inversion efficiency.

[0081] In some embodiments of this specification, constructing multiple filters may include: for a two-dimensional migration profile, setting six filters; 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 the two-dimensional migration profile can detect geological patterns in the seismic migration profile from different angles. Different filters have their own unique detection functions and can comprehensively capture the characteristics of geological patterns in different aspects. The vertical mean filter and horizontal mean filter can detect the central 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 the changes in curvature. By combining 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 information loss that may occur when detecting from a single angle.

[0083] Since the underground geological structure is often very complex, a single filter is difficult to accurately characterize its features. The six filters in this embodiment detect geological patterns from different angles and can better adapt to the complex and changeable geological structure. Whether it is the central position, edge features, or curvature changes of geological structures, effective detection and analysis can be carried out, thus providing valuable information for full-waveform inversion 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 truly reflect the actual relationship of geological patterns around different nodes, providing a reliable basis for constructing an accurate correlation coefficient matrix subsequently.

[0085] An accurate correlation coefficient matrix is crucial for full-waveform inversion. The comprehensive geological pattern information obtained through these six filters can construct a correlation coefficient matrix that better reflects the actual geological situation, thereby providing more accurate prior information for inversion. In full-waveform inversion, the accuracy of prior information helps to narrow the search space of solutions, reduce the number of possible solutions, reduce the multi-solution nature of inversion, and make the inversion result more likely to approach the true velocity model.

[0086] It can be understood that other numbers of filters can also be set. For example, two filters, three filters, four filters, or eight filters can be set, etc. The types of filters can be selected according to actual needs.

[0087] In some embodiments of this specification, applying multiple constructed 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 through the multiple filters; defining a plurality of score values corresponding to the geological patterns around each node in the seismic migration profile as a score vector corresponding to each node; calculating 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 generating a correlation coefficient matrix based on the multi-point pattern correlation between any two nodes.

[0088] Specifically, in full waveform inversion based on geological prior information, multiple constructed filters (such as vertical mean filter, horizontal mean filter, vertical gradient filter, horizontal gradient filter, vertical curvature filter, horizontal curvature filter) have different functions and can detect geological patterns in a seismic migration profile from different perspectives. When these filters are applied to the seismic migration profile, each filter will perform specific operations and analyses on the geological patterns in the profile, thereby converting complex geological patterns into specific score values. For example, the vertical mean filter may give a score based on the mean characteristics in the vertical direction around a node, and the horizontal gradient filter will give a corresponding score based on the gradient characteristics in the horizontal direction. These score values reflect the characteristics of geological patterns from different aspects.

[0089] After obtaining the score values of the geological patterns around each node under each filter, combining these multiple score values is defined as the score vector corresponding to that node. Each score vector can be regarded as a quantitative representation of the characteristics of the geological pattern around that node. Through the score vector, the comprehensive characteristics of the geological pattern around the node can be described more comprehensively, rather than just single-dimensional information.

[0090] According to the score vectors corresponding to any two nodes in the seismic migration profile, a suitable method 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 degree of the two score vectors. The higher the similarity, the more similar the geological patterns around these 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 more importantly, represents the degree of association of the geological patterns around them.

[0091] Based on the multi-point pattern correlations between any two nodes, arranging these correlation values in matrix form according to certain rules generates a correlation coefficient matrix. The elements of this matrix represent the strength of the geological pattern correlations between any two nodes in the seismic migration profile. The correlation coefficient matrix provides important prior information for subsequent full waveform inversion. It can help us better understand the associations of geological patterns at different positions in the seismic migration profile, thereby more reasonably considering the continuity and correlations of the geological structure during the inversion process and improving the accuracy and reliability of the inversion results.

[0092] Through the above method, the geological pattern information in the seismic migration profile can be quantified and sorted to obtain a correlation coefficient matrix that can reflect the geological pattern correlations. This matrix plays an important constraining role in full waveform inversion, guiding the inversion process towards a direction more consistent with the geological reality, reducing the non-uniqueness of the inversion, and improving the inversion efficiency and accuracy.

[0093] In some embodiments of this specification, using the Stein variational gradient descent algorithm to iteratively update the multiple velocity samples multiple times may include: given a step size and using the Stein variational gradient descent algorithm to update the multiple velocity samples; calculating the gradient of the full waveform inversion objective function according to the multiple velocity samples; calculating the gradient of the full waveform inversion posterior distribution using Bayes' formula according to 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 be stable.

[0094] A step size can be given, and the multiple velocity samples are updated according to 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 then the velocity samples are updated again. Repeat this process until the 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 be stable. After multiple iterative updates of the finally output multiple velocity samples, the inverted target velocity model is obtained.

[0095] According to the updated multiple velocity samples, the gradient of the full waveform inversion objective function is calculated through corresponding mathematical methods. The objective function is used to measure the matching degree between the velocity model and the actual observed data, and the gradient represents the change direction and rate of the objective function in the velocity model space. For example, in full waveform inversion, numerical simulation methods based on the acoustic wave equation may be used, combined with techniques such as the finite difference method, to calculate the gradient of the objective function with respect to the velocity model. By calculating the gradient, the change trend of the current velocity sample on the objective function can be understood, providing a direction for subsequent updates.

[0096] According to the gradient of the objective function calculated previously, the Bayesian formula is used to calculate the gradient of the full-waveform inversion posterior distribution. In the Bayesian framework, the posterior distribution combines prior information and observational data information, reflecting the probability distribution of the velocity model given the observational data. By calculating the gradient of the posterior distribution, the variation of the velocity model in the probability space can be further understood, and the uncertainty of the velocity model can be considered more comprehensively. For example, using the Bayesian formula, the gradient of the objective function is combined with the prior probability density distribution and the likelihood function to derive and calculate the expression of the gradient of the posterior distribution.

[0097] Repeat the above steps, that is, use the Stein variational gradient descent algorithm to update multiple velocity samples again, and then calculate the gradient of the objective function and the gradient of the posterior distribution. Continuously perform such an iterative process until the set stopping condition is met. The stopping condition can be that the number of iterations reaches the pre-set maximum number of iterations, which is a simple and direct way to stop, ensuring that the algorithm ends within a certain range of computing resources; it can also be that the full-waveform inversion posterior distribution tends to be stable, that is, the change in the posterior distribution is very small, indicating that the algorithm has converged to a relatively stable state. At this time, it is considered that the obtained velocity samples are close to the optimal solution and the iteration can be stopped.

[0098] In the above way, the Stein variational gradient descent algorithm can continuously adjust the velocity samples to gradually approach the optimal solution of the full-waveform inversion, while considering the uncertainty of the velocity model. This method has important application value in full-waveform inversion, can improve the accuracy and reliability of the inversion results, and provide a more accurate velocity model 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 the Stein variational gradient descent algorithm (SVGD), the kernel function plays a key role. The kernel function is used to map the data in the low-dimensional space to the high-dimensional space. Through this mapping, more complex operations and analyses can be performed in the high-dimensional space without explicitly calculating the coordinates of the high-dimensional space. The radial basis function k(m, m′) is as follows:

[0101]

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

[0103] In the SVGD algorithm, the radial basis function is used as the kernel function to correlate the gradient information of different samples. Specifically, when calculating the update direction of the samples, the similarity between different samples is calculated through the kernel function, and the gradient information of samples with high similarity is fused. As a result, when the algorithm updates the samples, it can take into account the mutual relationship between the samples and avoid the current sample falling into a local extremum.

[0104] The radial basis function can map the data into a high-dimensional space, making the data that is linearly inseparable in the low-dimensional space become linearly separable or easier to process in the high-dimensional space. In full waveform inversion, the posterior distribution of the velocity model is usually a complex multimodal distribution. The SVGD algorithm using the radial basis function as the kernel function can better approximate this complex posterior distribution, thereby improving the accuracy of the inversion. Compared with some other complex kernel functions, the calculation of the radial basis function is relatively simple. It only involves the calculation of the distance between sample points and exponential operations. In the case of a large number of samples, it can ensure the calculation efficiency to a certain extent. This is very important for processing a large number of velocity samples in full waveform inversion because full waveform inversion itself has a large amount of calculation and requires an efficient algorithm to support. The scale constant in the radial basis function 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, determining the value according to the median of the relative distances of the samples is an effective way to adjust the parameters according to the data characteristics, which can make the algorithm better adapt to the specific full waveform inversion problem.

[0105] In some embodiments of this specification, after using the Stein variational gradient descent algorithm to perform multiple iterative updates on the multiple velocity samples based on the gradient of the posterior distribution of the full waveform inversion and outputting the target velocity model inverted from the multiple velocity samples, it may further include: determining the posterior distribution of the target velocity model to evaluate the uncertainty of the inversion result.

[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 the observed data and prior information. The determination of the posterior distribution can be achieved by analyzing the distribution of multiple velocity samples obtained during the statistical iteration process. For example, calculating statistical quantities such as the mean and variance of the samples, or constructing a probability density function to describe the shape of the posterior distribution. Based on the posterior distribution of the determined target velocity model, the uncertainty of the inversion result can be evaluated. Characteristics such as the width and variance of the posterior distribution can reflect the degree of uncertainty of the velocity model. A wider posterior distribution indicates that the values of the velocity model have greater uncertainty, and there may be multiple reasonable velocity models; while a narrower posterior distribution indicates that the values of the velocity model are relatively more certain. By evaluating the uncertainty of the inversion result, it helps to better understand the reliability of the inversion result. In subsequent geological interpretation and exploration decision-making, the uncertainty factors of the velocity model can be more reasonably considered to reduce exploration risks.

[0107] Each embodiment in this specification is described in a progressive manner. For the same or similar parts among the embodiments, reference can be made to each other. Each embodiment focuses on the differences from other embodiments. Specifically, reference can be made to the descriptions of the relevant processing-related embodiments mentioned above, and details will not be repeated here.

[0108] The above describes 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 can be performed in a different order than in the embodiments and still achieve the desired results. Additionally, the processes depicted in the figures do not necessarily require the specific order or continuous order shown to achieve the desired results. In certain embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0109] The above method will be described below in conjunction with a specific embodiment. However, it should be noted that this specific embodiment is only for better explaining this specification and does not constitute an improper limitation of this specification.

[0110] In this specific embodiment, a variational inference full waveform inversion method based on geological prior information is provided. In this specific embodiment, the implementation of the variational inference full waveform inversion method based on geological prior information mainly includes the following: 1) setting a full waveform inversion uncertainty evaluation method based on the Stein variational gradient descent algorithm; 2) setting geostatistical prior information; 3) statistically analyzing the posterior distribution of the inversion result and evaluating its uncertainty.

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

[0112] 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 obs ), and this process can be described by the following formula:

[0113]

[0114] where d obs is the observed data, p(d obs ) is the normalization constant, and p(d obs |m) is the likelihood function.

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

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

[0117] where E q denotes the expectation with respect to q. By optimizing the objective function (2), the probability distribution q * (m) obtained when the KL divergence reaches its minimum is the best approximation of the posterior probability density distribution p(m|d ocbs ).

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

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

[0120]

[0121] where m ∈ [m1, m2, … m n are n velocity models, n is the number of samples, can be regarded as the iterative update direction vector corresponding to n velocity models. ε is a small constant and can be regarded as the step size. Let the probability density function after the transformation be q T (m). According to mathematical transformation, we know that:

[0122]

[0123] Among them, trace represents the trace of the matrix, which is called the Stein operator and satisfies:

[0124]

[0125] where Φ(m) is the vector obtained after the function Φ acts on the samples. According to the description of Equation (2), the KL divergence reflects the closeness of two probability density functions, and Equation (4) shows that: when there exists Φ(m) such that the expectation on the right side of Equation (4) reaches the 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. Along this direction, the probability density function q T (m) is closer to the posterior probability density function p(m|d obs ) to be obtained than q(m) before the update. Therefore, as long as we find the one that maximizes the following objective function:

[0126]

[0127] where F represents the functional space, transform the above optimization problem to the Hilbert space for solution, and we get:

[0128]

[0129] where m j represents the j-th velocity model (j = 1, 2,..., n) in m1, m2,... m n , and k(m, m′) is the kernel function, usually taking the radial basis function:

[0130]

[0131] where h > 0 is the proportionality constant, and this patent takes where is the median of the relative distances of all samples. Substitute Equation (7) into Equation (4), and take in Equation (3) as the negative gradient direction of the KL divergence, then the update of the initial samples can be realized. The SVGD algorithm is summarized as follows:

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

[0133] 2) Iteratively update the samples, and the iterative formula is:

[0134]

[0135] Among them, the superscript l represents the number of iterations, the subscript i represents the i-th sample (i = 1, 2, …, n), and ε l represents the step size of the l-th iteration, and respectively represent the current and the i-th sample after iterative update. The initial sample corresponds to l = 0.

[0136]

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

[0138]

[0139] where d syn (m) is the synthetic seismic data obtained by forward modeling based on m, and Σ d is the covariance matrix of the data, which is usually assumed to be a diagonal matrix. At this time, obtaining the gradient of the likelihood function can be transformed into obtaining the gradient of the conventional L2 objective function in full waveform inversion:

[0140]

[0141] The above formula 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, ρ = ρ(x) is the density with the unit of kg / m 3 , v = v(x) is the propagation velocity of the seismic wave, is the pressure field, and are the horizontal and vertical components of the particle vibration velocity respectively, where the superscript → represents the forward wave field along the direction of increasing time, that is, the wave field is calculated from the zero time towards the increasing time direction. x = (x, z) is the two-dimensional spatial coordinate, and x s =(x s , z s ) represents the spatial position where the source wavelet is located, and δ is the Dirac-Delta function. The adjoint wave equation satisfied by the adjoint wave field obtained using the adjoint state method is:

[0144]

[0145] where r represents the adjoint source, and are adjoint wave field variables, and the superscript ← represents that the wave field propagation is along the direction of decreasing time, that is, from the maximum time towards the zero time.

[0146] The gradient formula for updating the velocity in full waveform inversion can be obtained by solving the derivative of the objective function with respect to the velocity model v as follows:

[0147]

[0148] where T max is the maximum recording duration of the seismic data, in seconds. This equation shows that the gradient is the zero-delay cross-correlation of the source wavefield propagating forward in time and the adjoint wavefield propagating backward in time .

[0149] The SVGD algorithm generates a series of samples according to the prior distribution, and then updates the samples based on the gradient of the KL divergence objective function. During the iterative update process, the gradient of the data term and the prior information are introduced through the Bayesian formula. At the same time, the information of different samples is correlated through the kernel function k(m, m′). This correlation is equivalent to a kind of constraint. When a sample falls into a local extremum, its adjacent samples can help it jump out of the local extremum. The samples after multiple iterative updates constitute an approximate estimate of the posterior probability density distribution. In this specific embodiment, it is used for FWI velocity 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 carried out independently. Therefore, in this specific embodiment, multi-GPU parallel computing for samples is realized in the program design, thereby improving the computational efficiency of variational inference FWI.

[0150] (2) Set the geostatistical prior information.

[0151] In Bayesian inversion, the prior information is crucial for the posterior inversion results. Currently, in Bayesian full waveform inversion, extremely broad prior information is adopted, that is, each sampling point satisfies a uniform distribution with the same variance, which will greatly increase the computational amount. Therefore, generating prior samples with geological significance is beneficial to narrowing the search space of the solution and accelerating the convergence speed of the inversion.

[0152] In the field of stochastic geological modeling, the following method can be used to perturb and generate random samples:

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

[0154] where m represents the column vector formed by arranging the two-dimensional perturbation random velocity model in columns, m = [m1 m2 … m nx T , m1 (i = 1, 2, …, nx) represents the i-th column vector in the two-dimensional model (that is, the first nz elements of m form the first column of the two-dimensional velocity model, and so on. nx and nz are the number of columns and rows of the two-dimensional model respectively), m smooth ​is a vector composed of a smooth background model, and R is a random vector that follows a uniform distribution in the interval (0, R max ), where R max is the maximum perturbation amount, and Γ is a matrix composed of the correlation coefficients between any two points of 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, a square filter with a radius of l is set. A geological pattern p obtained by using this filter around the sampling point S(i, j) in the migration profile can be expressed as:

[0157]

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

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

[0160]

[0161] Among them, 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 central position of the image, f3 and f4 can detect the edge pattern of the image, and f5 and f6 can detect the change in curvature.

[0162] Through these six filters, the geological patterns captured from the seismic profile can be converted into score values. The six score values are: E1(i, h), E2(i, h), E3(i, j), E4(i, h), E5(i, h), and E6(i, h). The expressions of the six score values are:

[0163]

[0164] In the formula, k = 1, 2,..., 6, i ∈ [m + 1, nz - m], j ∈ [m + 1, nx - m]. All geological patterns in the migration profile will be converted into score values, and six score values will be obtained at each point. Here, the six score values corresponding to the geological pattern around a node S(i, j) are defined as a score vector, and its expression 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 section, a score vector can be obtained according to the above method. The similarity between two patterns can be judged based on the score vectors corresponding to two nodes. This similarity between patterns not only represents the relationship between the two nodes, but more importantly, represents the relationship between the geological patterns around the two nodes. 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′), and Corr is the correlation coefficient operation.

[0169] However, according to the correlation coefficients in formula (20) to form the correlation coefficient matrix Γ in formula (15), it is extremely large and requires a large amount of computing and storage resources. Moreover, this matrix indicates that there is a correlation between all points in the space, which does not conform to the underground geological law. Therefore, the variogram in two-point geostatistics is introduced. Define the variogram r(h) of the variable Z(x) between the point x and x + h in the x direction, denoted as:

[0170]

[0171] where Var represents variance. The variogram represents the incremental semi-variance of the variable between two points with a distance of h. Generally, the experimental variogram is statistically obtained from well logging data, and then the theoretical variogram is further fitted. It can be seen from the above formula that the variogram indicates that there is no correlation between sampling points in the space when the distance is greater than a certain distance, and this distance is called the range.

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

[0173]

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

[0175] (3) Statistically invert the posterior distribution of the result and evaluate its uncertainty.

[0176] Using the initial sample set with geological significance generated in the previous step, perform full waveform inversion based on the SVGD algorithm, as shown in equations (9)-(15). When certain termination conditions are met, the inversion ends. The termination conditions can be simply set as the number of iterations reaching the set maximum number of iterations, or by continuously statistically analyzing the posterior distribution of the inversion results. When the posterior distribution gradually stabilizes, the inversion ends.

[0177] Finally, calculate the mean and standard deviation of the inversion results of all samples, and select the posterior probability density distribution at the characteristic positions of the model to evaluate the uncertainty of its 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 under certain conditions, and these solutions are all possible results.

[0179] The mean can represent the average situation of multiple inversion solutions, and the standard deviation can illustrate the dispersion of these solutions. All velocities corresponding to non-zero probabilities (frequencies) 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 evaluation. These uncertainty evaluations provide experience for well location deployment and geological risk assessment in subsequent oil and gas exploration.

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

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

[0183] Step 2) Fit the theoretical variogram according to the logging data and determine the range. Only retain the correlation coefficients when the distance between two points is less than the range to obtain the low-rank approximated correlation coefficient matrix Γ′.

[0184] Step 3) Given the number of samples and the maximum perturbation amount of the model, and according to equation (16), perturb to generate a series of initial samples with geological significance.

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

[0186] Step 5) According to the gradient of the objective function, use Bayes' formulas (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 Eqs. (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, at which point the inversion ends and the final inverted velocity model of all samples is output.

[0188] Step 7) Statistically analyze the posterior distribution of the inversion results of all samples and conduct uncertainty evaluation.

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

[0190] Next, the variational inference full waveform inversion method based on geological prior information in the embodiments of this specification will be applied to an actual scenario. Specifically, it includes the following content.

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

[0192] Figure 2 The Marmousi model is shown. Figure 2 (a) is the accurate model, Figure 2 (b) is the random model sampled from a uniform prior distribution, Figure 2 (c) is the velocity prior distribution. To verify that the random samples based on geostatistical prior information proposed in this embodiment have geological significance and thus reduce the inversion complexity, the Figure 2 Marmousi velocity model shown in (a) is used for testing. The velocity model has a grid size of 200×120, a grid spacing of 20 m, and the maximum and minimum velocities are 4670 m / s and 1500 m / s, respectively. Figure 3 The reverse time migration profile of the Marmousi model is shown.

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

[0194] Applying the above six filters to the migrated section gives six filtered score sections as shown in Figure 5 the figure. Figure 5 The figure shows the filtered score sections obtained after filtering the migrated section. Figure 5 (a) is a vertical mean filtered section, Figure 5 (b) is a horizontal mean filtered section, Figure 5 (c) is a vertical gradient filtered section, Figure 5 (d) is a horizontal gradient filtered section, Figure 5 (e) is a vertical curvature filtered section, Figure 5 (f) is a horizontal curvature filtered section.

[0195] Figure 7 The figure shows the multi-point mode correlation matrix. Figure 7 (a) is the full-rank correlation coefficient matrix, Figure 7 (b) is the correlation coefficient matrix after low-rank approximation. The multi-point mode correlation coefficient matrix is obtained using Eqs. (20) and (21), as shown in Figure 7 (a) of the figure. Taking the accurate Marmousi model velocity at x = 2 km as the pseudo-well data to statistically experiment the variogram, and then fitting the theoretical variogram, the range is determined to be 180 m, as shown in Figure 6 the figure, which shows the curve of the variogram. Therefore, only consider that any point in the model is only correlated with the sampling points within the range of 180 m × 180 m around it, that is, only retain the correlation coefficients when the distance between two points is less than the range, so as to obtain the correlation coefficient matrix after low-rank approximation, removing the false correlations and reducing the memory and computational consumption, as shown in Figure 7 (b) of the figure.

[0196] Figure 8 The figure shows the establishment of a stochastic model using geostatistical priors. Figure 8 (a) is a smooth background model, Figure 8 (b) is a stochastic 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] Set the maximum perturbation in Eq. (16) to 800 m / s. Multiplying the correlation coefficient matrix after low-rank approximation by the random vector can obtain the stochastic perturbation based on geostatistical priors, as shown inFigure 8 As shown in (b) of, add it to Figure 8 the smooth velocity model shown in (a) of to obtain a stochastic model based on geostatistical prior, as Figure 8 shown in (c) of. Compared with the stochastic model obtained by sampling from the uniform prior distribution shown in (c) of Figure 2 , the stochastic model generated based on geostatistical prior has geological significance, thus narrowing the search space of the solution and improving the computational efficiency. Figure 2 As shown in (b) of

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

[0199] Use the above uniform prior information and geostatistical prior information to generate 100 random samples respectively as the initial models of variational inference FWI, and test the influence of different prior information on the full waveform inversion results and velocity uncertainty evaluation. Arrange 8 shot points at a depth of 360 m, with a lateral spacing of 500 m, and the first shot point is located at 200 m. The geophones are set on the surface, excited in the middle, and received by full arrays on both sides. The seismic record duration is 5 s, the sampling interval ΔT = 1 ms, and the source is a Ricker wavelet. When inverting, filter the data to the main frequencies of 4 Hz and 10 Hz respectively, and conduct two independent inversions: (1) First invert the data with a main frequency of 4 Hz, (2) Then use the low-frequency inversion results of all samples as their respective initial models, and continue to invert the data with a main frequency of 10 Hz.

[0200] Set the step size of the SVGD algorithm as the empirical step size of full waveform inversion. The inversion based on the two prior information iterates 100 times each in the inversions of 4 Hz and 10 Hz respectively. In addition, the inversion based on the uniform prior is separately set to iterate 600 times each in the inversions of 4 Hz and 10 Hz. Statistically calculate the mean and standard deviation of the inversion results of all samples as Figures 9 to 11 shown, and plot the marginal distributions at four positions of (0.6 km, 2 km), (1.2 km, 2 km), (1.8 km, 2 km), and (2.4 km, 2 km) as histograms, as Figure 11 shown, and the black dashed line is the accurate velocity at each position.

[0201] Figure 9 is the variational inference full waveform inversion result based on uniform prior information, with each frequency band iterating 100 times. The mean of the low-frequency inversion result (as Figure 9 shown in (a) of) does not recover the structural form of the true model and there is low-frequency noise; while the mean of the high-frequency inversion result (as Figure 9 shown in (c) of) can recover the structural form of the velocity model, but the inversion in the marginal and deep regions is poor and even falls into local extrema, which may be caused by insufficient illumination. Figure 9(b) shows the standard deviation of the low-frequency inversion results. The standard deviation of the high-frequency inversion results (as shown in Figure 9 (d)) exhibits structural features similar to the mean. From the marginal distributions of the final inversion results at the four positions shown in Figure 11 (a) to Figure 11 (d), it can be seen that the above inversion did not converge, the inversion results among samples varied greatly, and most deviated from the accurate velocity.

[0202] Therefore, the number of iterations for each frequency band was adjusted to 600 times, and the mean and standard deviation of the inversion results were statistically analyzed, as shown in Figure 10 . The mean of the low-frequency inversion results (as shown in Figure 10 (a)) has recovered the structural form of the velocity model; while the mean of the high-frequency inversion results (as shown in Figure 10 (c)) has improved resolution, however, the deep region is still trapped in local extrema. The standard deviation of the high-frequency inversion results (as shown in Figure 10 (d)) is lower than the standard deviation of the low-frequency inversion results (as shown in Figure 10 (b)). And the standard deviation statistically analyzed after 600 iterations (as shown in Figure 10 (d)) is generally slightly lower than the standard deviation statistically analyzed after 100 iterations (as shown in Figure 9 (d)), indicating that with the increase in the number of iterations, the uncertainty of the inversion has decreased and the inversion accuracy has improved. From the marginal distributions of the final inversion results at the four positions shown in Figure 12 (e) to Figure 12 (h), it can be seen that the inversion results of all samples at the position of (0.6 km, 2 km) converge to the accurate velocity, and the inversion results at the other three positions do not converge, indicating that when using uniform prior for variational inference full waveform inversion, the solution search space is large, and the inversion results of most samples differ greatly from the accurate model.

[0203] Figure 9 Using the prior information sampling in Figure 2 (c) to generate 100 initial samples, perform full waveform inversion based on the SVGD algorithm, statistically analyze the (a) mean and (b) standard deviation after 100 iterations with a dominant frequency of 4 Hz, use the inversion results at 4 Hz as the initial samples to continue the inversion at 10 Hz, and statistically analyze the (c) mean and (d) standard deviation after 100 iterations

[0204] Figure 10 Using the prior information sampling in Figure 2 (c) to generate 100 initial samples, perform full waveform inversion based on the SVGD algorithm, statistically analyze the (a) mean and (b) standard deviation after 600 iterations with a dominant frequency of 4 Hz, use the inversion results at 4 Hz as the initial samples to continue the inversion at 10 Hz, and statistically analyze the mean (as shown inFigure 10 as shown in (c) of Figure 10 and the standard deviation (as shown in (d) of

[0205] Figure 11 Using geostatistical prior information sampling to generate 100 initial samples, perform full waveform inversion based on the SVGD algorithm. After 100 iterations with a main frequency of 4 Hz, the mean value in (a) and the standard deviation in (b) are statistically analyzed. The inversion result at 4 Hz is used as the initial sample to continue the inversion at 10 Hz. After 100 iterations, the mean value (as shown in (c) of Figure 11 and the standard deviation (as shown in (d) of Figure 11 are statistically analyzed.

[0206] The results of variational inference full waveform inversion based on geostatistical prior information are as shown in Figure 12 . The mean value of the low-frequency inversion result (as shown in (a) of Figure 12 ) can better recover the structural form of the velocity model; while the mean value of the high-frequency inversion result (as shown in (c) of Figure 12 ) has further improved resolution, and the inversion results in the deep and edge regions are much better than the inversion mean based on the uniform prior distribution (as shown in (c) of Figure 9 , Figure 10 ). The standard deviation of the high-frequency inversion result (as shown in (d) of Figure 12 ) is lower than the standard deviation of the low-frequency inversion result (as shown in (b) of Figure 10 ), and is much lower than the inversion mean based on the uniform prior distribution (as shown in (d) of Figure 9 , Figure 10 ). From the marginal distributions of the final inversion results at four positions shown in (i) of Figure 12 to (l) of Figure 12 , it can be seen that the inversion results of all samples at the four positions have a small difference from the accurate velocity or converge well to the accurate velocity. This shows that in variational inference full waveform inversion, introducing geostatistical prior information can narrow the search space of the solution, and a small number of iteration times can make the inversion results of all samples good.

[0207] Figure 12 Shows the results of variational inference full waveform inversion based on geostatistical prior information. Specifically, Figure 12 (a) to (d) of Figure 9 show the marginal distributions of the final inversion results of sample 10 Hz at different positions. Figure 12 The position corresponding to (a) of Figure 12 is (0.6 km, 2 km), Figure 12 the position corresponding to (b) of Figure 12The position corresponding to (d) is (2.4 km, 2 km). (e) to (h) of 12 show the calculation Figure 10 The marginal distribution of the final inversion result of sample 10 Hz at different positions. Figure 12 The position corresponding to (e) is (0.6 km, 2 km), Figure 12 The position corresponding to (f) is (1.2 km, 2 km), Figure 12 The position corresponding to (g) is (1.8 km, 2 km), Figure 12 The position corresponding to (h) is (2.4 km, 2 km). Figure 12 (i) to (l) of show the calculation Figure 11 The comparison diagram of the marginal distribution of the final inversion result of sample 10 Hz at different positions. Figure 12 The position corresponding to (i) is (0.6 km, 2 km), Figure 12 The position corresponding to (j) is (1.2 km, 2 km), Figure 12 The position corresponding to (k) is (1.8 km, 2 km), Figure 12 The comparison diagram of the marginal distribution at the position corresponding to (l) which is (2.4 km, 2 km). The velocity shown by the black dotted line is the accurate velocity at the corresponding position.

[0208] In this specific embodiment, 8 NVIDIA Tesla P40 graphics cards are used to test the inversion time-consuming for the above different iteration times. The calculation time-consuming for 100 samples iterated 100 times is 14 h, and the calculation time-consuming for iterating 600 times is 84 h. Therefore, from the perspective of velocity uncertainty evaluation, compared with the uniform prior information, the initial samples established based on the geostatistical prior have geological significance, reduce the search space of the solution, lower the inversion difficulty, greatly improve the calculation efficiency, and improve the inversion accuracy and reduce the non-uniqueness of the inversion.

[0209] Based on the same inventive concept, an apparatus for variational inference full waveform inversion based on geological prior information is also provided in the embodiments of this specification, as described in the following embodiments. Since the principle of the apparatus for variational inference full waveform inversion based on geological prior information to solve problems is similar to that of the method for variational inference full waveform inversion based on geological prior information, the implementation of the apparatus for variational inference full waveform inversion based on geological prior information can refer to the implementation of the method for variational inference full waveform inversion based on geological prior information, and the repeated parts will not be elaborated. Hereinafter, the term "unit" or "module" can be a combination of software and / or hardware that can achieve a predetermined function. Although the apparatuses described in the following embodiments are preferably implemented in software, implementation in hardware, or a combination of software and hardware is also possible and contemplated. Figure 13 is a structural block diagram of the apparatus for variational inference full waveform inversion based on geological prior information in the embodiments of this specification, asFigure 13 As shown, it includes: a calculation module 1301, an approximation module 1302, a generation module 1303, and an inversion module 1304. The following is an explanation of this structure.

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

[0211] The approximation module 1302 is used to perform low-rank approximation on the calculated correlation coefficient matrix according to the theoretical variogram fitted from well logging data 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 low-rank approximated correlation coefficient matrix;

[0213] The inversion module 1304 is used to calculate the gradient of the full waveform inversion objective function according to the multiple velocity samples; it is also used to calculate the gradient of the full waveform inversion posterior distribution using Bayes' formula according to the gradient of the objective function; it is also used to perform multiple iterative updates on 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 inverted from the multiple velocity samples. In some embodiments of this specification,

[0214] In some embodiments of this specification, the device further includes a determination module, and the determination module is used to: determine the posterior distribution of the target velocity model to evaluate the uncertainty of the inversion result.

[0215] In some embodiments of this specification, the calculation module is specifically used to: for a two-dimensional migration profile, set six filters; 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 the multiple score values corresponding to the geological patterns around each node in the seismic migration profile as the score vector corresponding to each node; calculate the similarity of the geological patterns around any two nodes according to the score vectors corresponding to any two nodes in the seismic migration profile to obtain the multi-point pattern correlation between any two nodes; 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 multiple velocity samples for multiple times, including: given a step size and using the Stein variational gradient descent algorithm to update the multiple velocity samples; calculating the gradient of the full waveform inversion objective function according to the multiple velocity samples; calculating the gradient of the posterior distribution of the full waveform inversion using the Bayesian formula according to 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 posterior distribution of the full waveform inversion tends to be stable.

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

[0219] From the above description, it can be seen that the embodiments of this specification achieve the following technical effects: a correlation coefficient matrix of a geostatistical model with low-rank approximation is proposed, false correlations are removed, and memory and computational consumption are reduced. Random initial samples with geological significance are generated based on the low-rank approximation correlation coefficient matrix to make the prior information more accurate. Then, based on the geostatistical prior information, the Stein variational gradient descent algorithm is used for efficient variational inference full waveform inversion, greatly reducing the non-uniqueness of the inversion and improving the inversion efficiency.

[0220] The embodiments of this specification also provide a computer device, which can specifically refer to Figure 14 the schematic diagram of the composition structure of the computer device for the variational inference full waveform inversion method based on geological prior information provided by the embodiments of this specification. The computer device can specifically include an input device 141, a processor 142, and a memory 143. Among them, the memory 143 is used to store instructions executable by the processor. 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 may specifically be one of the main devices for information exchange between a user and a computer system. The input device may include a keyboard, a mouse, a camera, a scanner, a light pen, a handwriting input board, a 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 may also acquire and receive data transmitted from other modules, units, and devices. The processor may be implemented in any suitable manner. For example, the processor may take the form of, for example, a microprocessor or a processor and a computer-readable medium storing computer-readable program code (such as software or firmware) executable by the (micro)processor, logic gates, switches, an application specific integrated circuit (ASIC), a programmable logic controller, and an embedded microcontroller, and so on. The memory may specifically be a memory device for storing information in modern information technology. The memory may include multiple levels. In a digital system, anything that can store binary data can be a memory; in an integrated circuit, a circuit with a storage function without a physical form is also called a memory, such as a RAM, a FIFO, etc.; in a system, a storage device with a physical form is also called a memory, such as a memory module, a TF card, etc.

[0222] In this embodiment, the functions and effects specifically implemented by this computer device may be explained by comparison with other embodiments and will not be elaborated here.

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

[0224] In this embodiment, the above storage medium includes but is not limited to a random access memory (RAM), a read-only memory (ROM), a cache, a hard disk drive (HDD), or a memory card. The memory may be used to store computer program instructions. The network communication unit may be set according to the standards specified by the communication protocol and is an interface for network connection communication.

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

[0226] Obviously, those skilled in the art should understand that the various modules or steps of the embodiments of the present specification described above can be implemented by a general-purpose computing device. They can be concentrated on a single computing device or distributed on a network composed of multiple computing devices. Optionally, they can be implemented by program codes executable by the computing device. Thus, they can be stored in a storage device and executed by the computing device. And in some cases, the steps shown or described can be executed in a sequence different from that here, or they can be separately fabricated into individual integrated circuit modules, or multiple modules or steps among them can be fabricated into a single integrated circuit module for implementation. In this way, the embodiments of the present specification are not limited to any specific combination of hardware and software.

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

[0228] The above are only the preferred embodiments of this specification and are not used to limit this specification. For those skilled in the art, various changes and modifications can be made to the embodiments of this specification. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of this specification shall be included within the protection scope of this specification.

Claims

1. A variational inference full waveform inversion method based on geological prior information, characterized in that Including: Constructing multiple filters; Applying the constructed multiple filters to a seismic migration profile to calculate a correlation coefficient matrix; Performing low-rank approximation on the calculated correlation coefficient matrix according to a theoretical variogram fitted from well logging data to obtain a low-rank approximated correlation coefficient matrix; Randomly generating multiple velocity samples based on the low-rank approximated correlation coefficient matrix; Calculating the gradient of the full waveform inversion objective function according to the multiple velocity samples; calculating the gradient of the full waveform inversion posterior distribution using Bayes' formula according to the gradient of the objective function; based on the gradient of the full waveform inversion posterior distribution, using the Stein variational gradient descent algorithm to perform multiple iterative updates on the multiple velocity samples, and outputting the target velocity model inverted from the multiple velocity samples.

2. The variational inference full waveform inversion method based on geological prior information according to claim 1, wherein After performing multiple iterative updates on 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 inverted from the multiple velocity samples, it further includes: Determining the posterior distribution of the target velocity model to evaluate the uncertainty of the inversion result.

3. The variational inference full waveform inversion method based on geological prior information according to claim 1, wherein Constructing multiple filters, including: For a two-dimensional migration profile, setting six filters; 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.

4. The variational inference full waveform inversion method based on geological prior information according to claim 1, wherein Applying the constructed multiple filters to a seismic migration profile to calculate a correlation coefficient matrix, including: Converting the geological patterns captured from the seismic migration profile into score values through the multiple filters; Defining the multiple score values corresponding to the geological patterns around each node in the seismic migration profile as the score vector corresponding to each node; Calculating the similarity of the geological patterns around any two nodes according to the score vectors corresponding to any two nodes in the seismic migration profile to obtain the multi-point pattern correlation between any two nodes; Generating a correlation coefficient matrix based on the multi-point pattern correlation between any two nodes.

5. The variational inference full waveform inversion method based on geological prior information according to claim 1, characterized in that Performing multiple iterative updates on the multiple velocity samples using the Stein variational gradient descent algorithm, including: Given a step size and using the Stein variational gradient descent algorithm to update the multiple velocity samples; Calculating the gradient of the full waveform inversion objective function according to the multiple velocity samples; calculating the gradient of the full waveform inversion posterior distribution using Bayes' formula according to 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 be stable.

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

7. A variational inference full waveform inversion device based on geological prior information, characterized in that, Including: A calculation module, configured to construct multiple filters; and further configured to apply the constructed multiple filters to a seismic migration profile to calculate a correlation coefficient matrix; An approximation module, configured to perform low-rank approximation on the calculated correlation coefficient matrix according to a theoretical variogram fitted from well logging data to obtain a low-rank approximated correlation coefficient matrix; A generation module, configured to randomly generate multiple velocity samples based on the low-rank approximated correlation coefficient matrix; An inversion module, configured to calculate the gradient of the full-waveform inversion objective function according to the multiple velocity samples; It is further configured to calculate the gradient of the full-waveform inversion posterior distribution by using the Bayesian formula according to the gradient of the objective function; it is further configured to perform multiple iterative updates on the multiple velocity samples by 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.

8. The variational inference full waveform inversion device based on geological prior information according to claim 7, characterized in that The calculation module is specifically configured to: convert the geological patterns captured from the seismic migration profile into score values through the multiple filters; define the multiple score values corresponding to the geological patterns around each node in the seismic migration profile as the score vector corresponding to each node; calculate the similarity of the geological patterns around any two nodes according to the score vectors corresponding to any two nodes in the seismic migration profile, and obtain the multi-point pattern correlation between any two nodes; Generate a correlation coefficient matrix based on the multi-point pattern correlation between any two nodes.

9. A computer device, characterized in that, It includes a processor and a memory for storing processor-executable instructions, and when the processor executes the instructions, it implements the steps of the method according to any one of claims 1 to 6.

10. A computer-readable storage medium having computer instructions stored thereon, characterized in that, When the instructions are executed by the processor, it implements the steps of the method according to any one of claims 1 to 6.

Images