Well logging lithofacies classification method integrating rock physics model and Bayesian inversion

Abstract

The present invention provides a well logging lithofacies classification method that integrates a rock physics model and Bayesian inversion. This method includes: Step 1: Detailed analysis of local geological and physical rock characteristics to identify and classify lithologies; Step 2: Performing lithofacies petrophysical analysis to obtain a statistically representative training database; Step 3: Performing Markov Chain Monte Carlo random inversion based on Bayesian inversion; and Step 4: Synthesizing the well logging lithofacies data with actual lithofacies data. This method overcomes the traditional difficulty in solving high-dimensional integrals by using simulation to computationally solve high-dimensional integrals, reducing the difficulty and increasing the practicality of the method.

Description

Technical Field

[0001] The present invention relates to the field of well logging lithofacies classification, and in particular to a well logging lithofacies classification method integrating a rock physics model and Bayesian inversion. Background Art

[0002] In the field of oil and gas exploration, geomorphological analysis and classification have become crucial tasks for petroleum geologists. Lithofacies have a significant influence on reservoir geometry and porosity distribution. Therefore, classifying lithofacies using rock physics models can enhance the ability to interpret sedimentary systems and predict associated reservoir sand bodies using ultra-fine seismic amplitude information. Regardless of the specific method used, a key requirement for classification is the establishment of a good, statistically representative training database. Careful rock physics and geological analysis can identify the primary logging lithofacies classes. Core data, thin section descriptions, and estimates of fluid and wellbore invasion effects all assist in this process. Once a representative training dataset is obtained, it can be used as input for various classification techniques. This allows statistical prediction methods to be integrated with geological interpretation and rock physics analysis, providing better control over reservoir characterization results. Furthermore, seismic lithofacies identified from well logs can be corrected based on seismic amplitude, thereby improving the ability to characterize oil and gas reservoirs from seismic data.

[0003] Seismic lithofacies are sedimentary units on the seismic scale that are seismically resolvable if the thickness of the sedimentary unit is greater than one-quarter of the seismic wavelength.

[0004] Geostatistical inversion plays a crucial role in reservoir prediction. It consists of two parts: modeling and stochastic inversion. Bayesian inversion, a novel uncertainty inversion method based on Bayesian theory and combined with a stochastic inversion algorithm, has demonstrated strong practical application in well logging lithofacies prediction. Random sampling is performed within a hypothetical distribution of reservoir parameters to construct a Markov chain. The stability properties of the Markov chain are exploited to stabilize the posterior probability, resulting in more accurate classification results.

[0005] Chinese patent application number CN201510094046.9 describes a method for simultaneously implementing seismic lithofacies identification and its quantitative uncertainty assessment. The method includes: determining lithofacies type and performing well logging lithofacies definition; establishing a rock physical response relationship between well logging physical properties and elasticity; establishing a probabilistic statistical relationship between lithofacies and well logging attributes; establishing a probabilistic statistical relationship between well-seismic elastic parameters and constructing a statistical relationship between lithofacies and seismic elastic parameters; inverting the probability distribution information of target layer elastic parameters; combining the inverted target layer elastic parameter probability information with the statistical relationship between lithofacies and seismic elastic parameters to obtain lithofacies probability information for the target layer; obtaining a maximum a posteriori solution for the lithofacies distribution based on this lithofacies probability information, and outputting the final model parameters. This method quantitatively characterizes the uncertainty of each step in lithofacies identification and its propagation and accumulation characteristics within the lithofacies identification process. It enables uncertainty analysis of seismic lithofacies identification, reduces reservoir evaluation risks, and has a wide range of applications.

[0006] In the Chinese patent application with application number CN201310455938.8, a method for seismic inversion of carbonate rock physical parameters is disclosed, which belongs to the field of petroleum geophysical exploration. The method comprises: (1) performing pre-stack AVO three-parameter inversion based on pre-stack angle gathers to obtain formation elastic parameters M, wherein M includes P-wave velocity, S-wave velocity, and density; (2) establishing a statistical rock physics model based on well logging data; (3) combining the rock physics model obtained in step (2) with stochastic simulation of reservoir physical conditions to obtain stochastic simulation results; (4) performing Bayesian classification simulation on the stochastic simulation results obtained in step (3) to obtain a posterior probability distribution, using M obtained in step (1) as the inversion input, and finding R with the maximum posterior probability distribution is the final inversion result.

[0007] In the Chinese patent application with application number CN201510633989.4, a phase-controlled porosity inversion method based on Bayesian classification is disclosed, which belongs to the field of petroleum geological exploration technology. The phase-controlled inversion method includes the following steps: (1) performing intersection analysis on the well logging curve of rock porosity and the well logging curves of multiple rock elastic parameters obtained through laboratory tests to obtain at least one rock elastic parameter data sensitive to rock porosity; (2) obtaining the screened porosity and the screened elastic parameter data; (3) establishing the initial well seismic model and obtaining the prior probability distribution of the screened porosity; (4) obtaining the resampled porosity and obtaining the fitted screened elastic parameter data and the fitted resampled elastic parameter data; (5) obtaining the resampling error; (6) obtaining the reconstructed elastic parameter data; (7) obtaining the prior sample sequence and performing Bayesian classification to obtain the prior probability and posterior probability distribution; (8) converting the known elastic parameter data volume into the porosity data volume through Bayesian inversion.

[0008] The above existing technologies are significantly different from the present invention and fail to solve the technical problem we want to solve. Therefore, we have invented a new logging lithofacies classification method that integrates rock physics model and Bayesian inversion. Summary of the Invention

[0009] The purpose of the present invention is to provide a dynamic simulation that realizes that the sampling distribution changes as the simulation progresses, making up for the defect that traditional Monte Carlo integration can only be a static simulation. It can combine multiple information and has a higher accuracy than traditional inversion methods by integrating rock physics models and Bayesian inversion logging lithofacies classification methods.

[0010] The object of the present invention can be achieved by the following technical measures: a well logging lithofacies classification method that integrates a rock physics model and Bayesian inversion, the well logging lithofacies classification method that integrates a rock physics model and Bayesian inversion includes:

[0011] Step 1: Analyze local geological and physical rock characteristics in detail to identify and classify lithology;

[0012] Step 2: Perform petrographic rock physics analysis to obtain a statistically representative training database;

[0013] Step 3, perform Markov chain Monte Carlo random inversion based on Bayesian inversion;

[0014] Step 4: Synthesize the logging lithofacies data and the actual lithofacies data.

[0015] The purpose of the present invention can also be achieved by the following technical measures:

[0016] In step 1, the logging data are divided into seismic rock layer groups and sedimentary unit groups with shear-like structural characteristics. Through detailed analysis of local geological and physical rock characteristics, the understanding of the variability of logging data is improved, the optimization of logging parameters is guided, and lithology identification and classification are carried out.

[0017] In step 1, well logging lithofacies are classified by lithology including sand, silt, and clay, bedding structure including massive, interbedded, or disordered, lithofacies including grain size, clay location, and cementation, and seismic properties including compressional wave velocity, shear wave velocity, and density.

[0018] In step 1, the relationship between velocity and porosity is analyzed from the existing logging data. The core data is comprehensively analyzed based on the following knowledge of physical and rock properties to classify the lithofacies:

[0019] (1) The cementation of rock phases produces a stiffening effect, resulting in relatively high seismic velocities and impedance;

[0020] (2) Cemented sand bodies, layered sand bodies, and interlayer sandstones and shales also have relatively high velocities;

[0021] (3) As the clay content in the rock increases, the density of the rock will increase;

[0022] (4) Silty shale is denser than pure shale;

[0023] (5) For particle-supported sediments, an increase in clay content reduces porosity, increases density, and hardens the rock;

[0024] (6) For clay minerals, due to their inherent porosity, the porosity will increase with the increase of clay mineral content, and the rock skeleton will be weakened. Therefore, when the clay content of the rock is 40%, its wave propagation speed will reach its peak value;

[0025] (7) Since clay particles are flaky, the shear strength of shale is relatively weak, so the ratio of P-wave propagation velocity to S-wave propagation velocity in shale is higher than that in sandstone;

[0026] Based on the existing logging core data, the lithofacies can be divided into six types: gravel and conglomerate, thick sandstone, interlayer sandstone and shale, shale containing a large amount of silt (content is about greater than 30%), relatively pure shale, and disordered sedimentation.

[0027] In step 2, the gamma ray and P-wave velocities of the logging data are analyzed, and by combining physical and geological knowledge, key logging parameters for formation classification are selected to obtain the lithofacies classification data of the logging data and obtain a representative training data set.

[0028] In step 2, the logging data is subjected to median filtering to remove high-frequency noise and smoothed so that the thickness of the individual phases classified from the logging data is better correlated with the seismic significance layer, thereby obtaining better clustering results.

[0029] In step 2, a type of well is selected from the logging data to identify seismic lithofacies. Since gamma rays are a good indicator of clay in quartz-rich sediments, gamma rays and P-wave velocity are combined to determine different lithofacies and used together as classification parameters to classify the lithofacies.

[0030] In step 3, assuming that the reservoir parameters obey Gaussian distribution, the existing Gibbs sampling and Metropolis-Hastings sampling algorithms (i.e., MH algorithms) are used to extract a Markov chain whose own transition probability is only related to the previous state when the state transitions, thus completing the Markov process.

[0031] In step 3, the MH sampling method is used to randomly search the target space, and the obtained points are selected using the judgment criteria. When the generated random number is less than the acceptance probability, the obtained point is recorded. Otherwise, re-sampling is performed and the points are re-taken and recorded, and the parameter candidate points at the next position are comprehensively counted. Repeat the above steps. When all parameter judgments are completed or the iteration termination conditions are reached, the inversion result value is output.

[0032] In step 3, for the facies data, a Markov variable is obtained. As the state changes, the transition probability of the variable is only related to the previous state and has nothing to do with other states.

[0033] P(X t+1 =S j |X0=S0,X1=S1,…,X t =S t )=P(X t+1 =S j |X t =S i ) (1)

[0034] X represents a random variable, X t is the state of variable X at a certain time t, S0, S1, ..., S i ,S j ∈Ω is all possible cases of X, P(X t+1 =S j |X t =S i ) indicates that X t =S i Status changes to X t+1 =S j Probability of state;

[0035] The expression of state transition probability is:

[0036] P(i→j):=P i,j =P(X t+1 =S j |X t =S i ) (2)

[0037] Set the variable X to S at a certain time t k The probability of Then the value of the variable at time t+1 is S i The probability of

[0038]

[0039] Assume that there are n cases of X, then:

[0040]

[0041] Eventually converges to a stationary distribution π (*) ,Right now: π (*) Satisfy π (*) P=π (*) ;

[0042] Due to the constraints of the actual situation and the influence of other factors in the algorithm solution process, π(x) and P do not satisfy π(i)P i,j =π(j)P j,i , so that the final Markov chain cannot be stabilized, so an acceptance probability α is added i,j , so that the equation holds true, that is: π(i)α i,j P i,j =π(j)α j,i P j,i

[0043] Acceptance probability α i,j The formula is:

[0044]

[0045] The acceptance probability is used to determine whether to perform the next update. Through continuous state updates and judgments, the target distribution eventually reaches a stable condition.

[0046] In step 3, the specific steps of MH sampling are:

[0047] (1) Given the state transition matrix P and the assumed stationary distribution π(x), set the number of state transitions n1 and the required number of sample points n2;

[0048] (2) Given an initial distribution, sample from the distribution to obtain the initial state x0;

[0049] (3) Assuming that the time variation parameter (iterative loop parameter) is t, the maximum number of iterations is n1+n2-1, and t starts the loop at 0, perform the following steps:

[0050] i. According to the conditional probability distribution P(x|x1), take the sample point x;

[0051] ii. Sampling from the [0,1] distribution; i.e., position u;

[0052] iii. If satisfied Then accept this sampling and update the Markov chain state. At this time, x t+1 =x;

[0053] iv. If the above conditions are not met, the update will not be accepted. t+1 =xt .

[0054] In step 3, the processing method for logging data is as follows: ① Based on the experience of geostatistics, an initial model of logging lithofacies can be given as a prior distribution model; ② The lithofacies division in step 2 is performed according to the existing logging data, and the mean and variance of each lithofacies reservoir parameter are counted respectively; ③ In the target solution space, the MH sampling method is used to randomly search the docking space to obtain candidate seed points; ④ The MH judgment criterion is used to select the above seed points to determine whether to accept the candidate point, so as to count the parameter candidate points at the next position; ⑤ Repeat the above steps ③-④; when all parameter judgments are completed or the loop termination condition is reached, the inversion lithofacies result is output.

[0055] In step 4, by comparing the resolution of the inversion results with the lithofacies classification results, the accuracy of the algorithm is verified, its final effect is analyzed, and a conclusion is drawn on whether the algorithm is effective.

[0056] The well logging lithofacies classification method, which integrates a rock physics model and Bayesian inversion, utilizes known well logging data to establish a rock physics model, determine training data, and then classify well logging lithofacies using Bayesian inversion techniques based on Markov Chain Monte Carlo. First, lithofacies are defined using geological knowledge, and existing well logging data are classified. Rock physics analysis of these lithofacies is then performed to generate a robust, statistically representative training database for subsequent inversion. MH sampling is then performed, and a stable Markov chain is generated through iteration, ultimately yielding high-resolution well logging lithofacies inversion results.

[0057] This method classifies well logging data into seismic rock groups and sedimentary units with shear-like structural features. Research has shown that detailed analysis of local geological and physical rock characteristics can improve understanding of logging data variability, guide optimization of logging parameters, and perform lithologic identification and classification. Because density logging is occasionally corrupted by erosion and rough borehole surfaces in real-world logging situations, S-wave velocity is of limited use in classification. Therefore, gamma-ray and P-wave velocities are used as parameters for lithofacies classification. Petrophysical analysis of the logging data is performed to draw lithofacies conclusions and extract a feature training set.

[0058] Based on the concept of Bayesian inversion, the Markov Chain Monte Carlo method is used to invert and classify well logging lithofacies. By introducing the Markov process into the Monte Carlo simulation through a computer, a dynamic simulation is achieved in which the sampling distribution changes as the simulation progresses. This overcomes the limitation of traditional Monte Carlo integration, which is limited to static simulation. It can incorporate multiple information and achieve higher accuracy than traditional inversion methods. In practice, the Markov Chain Monte Carlo method has overcome the traditional difficulty in solving high-dimensional integrals. By using simulation to solve high-dimensional integrals, the method can be computationally solved, reducing the difficulty and increasing its practicality. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 A flowchart for the implementation of the well logging lithofacies classification method that integrates rock physics modeling and Bayesian inversion;

[0060] Figure 2 Schematic diagram showing the accuracy of the inversion results obtained by the classification method using only P-wave velocity and the lithofacies classification method using Gamma rays and P-wave velocity as features in this paper;

[0061] Figure 3 Schematic diagram of the classification success rates of the six lithofacies using only Gamma rays (GR), only P-wave velocity, and both GR and P-wave velocity (the method adopted in the present invention);

[0062] Figure 4 From left to right are the lithofacies of the test well, the lithofacies obtained by the traditional complete probability distribution function (PDF) inversion method, and the lithofacies classification method of the well logging lithofacies classification using the fusion rock physics model and Bayesian inversion used in the present invention;

[0063] Figure 5 Schematic diagram comparing the inversion average classification accuracy obtained by using the traditional complete probability distribution function (PDF) inversion method, the traditional Markov chain Monte Carlo inversion method based on the Bayesian framework, and the logging lithofacies classification method that integrates the rock physics model and Bayesian inversion adopted in the present invention. DETAILED DESCRIPTION

[0064] It should be noted that the following detailed descriptions are exemplary and intended to provide further explanation of the present invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by those skilled in the art to which the present invention belongs.

[0065] It should be noted that the terms used herein are only for describing specific embodiments and are not intended to limit the exemplary embodiments according to the present invention. As used herein, unless the context clearly indicates otherwise, the singular form is intended to include the plural form. In addition, it should be understood that when the terms "comprise" and / or "include" are used in this specification, they indicate the presence of features, steps, operations and / or combinations thereof.

[0066] like Figure 1 As shown, Figure 1 The flowchart of a well logging lithofacies classification method that integrates rock physics model and Bayesian inversion in an embodiment of the present invention is as follows. The well logging lithofacies classification method that integrates rock physics model and Bayesian inversion includes the following steps:

[0067] In step 101, the logging data are divided into seismic rock layer groups and sedimentary unit groups with shear-like structural characteristics. Through detailed analysis of local geological and physical rock characteristics, the understanding of the variability of logging data can be improved, and the optimization of logging parameters can be guided to perform lithology identification and classification.

[0068] Well logging lithofacies are classified based on lithology (sand, silt, and clay), bedding structure (massive, interbedded, or disordered), lithofacies (grain size, clay location, and cementation), and seismic properties (P-wave velocity, S-wave velocity, and density). The relationship between velocity and porosity is analyzed using existing well logging data. Lithofacies classification is then performed through a comprehensive analysis of core data based on the following knowledge of physical rock properties.

[0069] (1) The cementation of rock phases produces a stiffening effect, and the seismic velocity and impedance will be relatively high.

[0070] (2) Cemented sand bodies, layered sand bodies, and interlayer sandstones and shales also have relatively high velocities.

[0071] (3) As the clay content in the rock increases, the density of the rock will increase.

[0072] (4) Silty shale is denser than pure shale.

[0073] (5) For particle-supported sediments, an increase in clay content reduces porosity, increases density, and hardens the rock.

[0074] (6) For clay minerals, due to their inherent porosity, the porosity will increase with the increase of clay mineral content, and the rock skeleton will be weakened. Therefore, when the clay content of the rock is about 40%, its wave propagation velocity will reach its peak.

[0075] (7) Since clay particles are flaky, the shear strength of shale is relatively weak, so the ratio of the P-wave propagation velocity to the S-wave propagation velocity in shale is higher than that in sandstone.

[0076] The existing logging core data can be divided into 6 types of lithofacies.

[0077] In step 102, the gamma ray and P wave velocities of the logging data are analyzed, and combined with good physical and geological knowledge, "key" logging parameters for formation classification are selected to obtain the lithofacies classification data of the logging and obtain a representative training data set.

[0078] The well logging data is median filtered to remove high-frequency noise and smoothed. This improves the correlation between the thickness of the individual facies classified from the logging data and the seismic layers of significance, resulting in better clustering results. A well type is selected from the logging data to identify seismic lithofacies. Because gamma rays (GR) are a good indicator of clay in quartz-rich sediments, gamma rays and P-wave velocity are combined to identify different lithofacies. These two parameters are used together as classification parameters to classify the lithofacies.

[0079] In step 103, assuming that the reservoir parameters follow a Gaussian distribution, the existing Gibbs sampling and Metropolis-Hastings algorithms are used to extract a Markov chain whose transition probability is only related to the previous state during state transition, completing the Markov process. Taking MH sampling as an example, the MH sampling method is used to randomly search the target space, and the resulting points are selected using a judgment criterion. When the generated random number is less than the acceptance probability, the resulting point is recorded. Otherwise, resampling is performed, and the points are recorded. A comprehensive statistical analysis is then performed to determine the parameter candidate point for the next position. The above steps are repeated. When all parameter judgments are completed or the iteration termination condition is met, the inversion result value is output.

[0080] For lithofacies data, a Markov variable is obtained. As the state changes, the transition probability of the variable is only related to the previous state and has nothing to do with other states.

[0081] P(X t+1 =S j |X0=S0,X1=S1,…,X t =S t )=P(X t+1 =S j |X t =S i ) (1)

[0082] S0,S1,...,S i ,S j ∈Ω is all possible cases of X.

[0083] The expression of state transition probability is:

[0084] P(i→j):=P i，j =P(X t+1 =S j |X t =S i ) (2)

[0085] Set the variable X to S at a certain time t k The probability of Then the value of the variable at time t+1 is Si The probability of

[0086]

[0087] Assume that there are n cases of X, then:

[0088]

[0089] Eventually converges to a stationary distribution π (*) ,Right now: π (*) Satisfy π (*) P=π (*) .

[0090] Due to the constraints of the actual situation and the influence of other factors in the algorithm solution process, π(x) and P do not satisfy π(i)P i,j =π(j)P j,i , so that the final Markov chain cannot be stabilized, so an acceptance probability α is added i,j , so that the equation holds true, that is: π(i)α i,j P i,j =π(j)α j,i P j,i

[0091] Acceptance probability α i,j The formula is:

[0092]

[0093] The acceptance probability is used to determine whether to perform the next update. Through continuous state updates and judgments, the target distribution eventually reaches a stable condition.

[0094] The specific steps of MH sampling are:

[0095] (1) Given the state transition matrix P and the assumed stationary distribution π(x), set the number of state transitions n1 and the required number of sample points n2;

[0096] (2) Given an initial distribution, sample from the distribution to obtain the initial state x0;

[0097] (3) Assuming that the time variation parameter (iterative loop parameter) is t, the maximum number of iterations is n1+n2-1, and t starts the loop at 0, perform the following steps:

[0098] i. According to the conditional probability distribution P(x|x1), take the sample point x;

[0099] ii. Sampling from the [0,1] distribution. That is, position u;

[0100] iii. If satisfied Then accept this sampling and update the Markov chain state. At this time, x t+1 =x;

[0101] iv. If the above conditions are not met, the update will not be accepted. t+1 =x t .

[0102] In summary, the following methods are used to process well logging data: ① Based on geostatistical experience, an initial model of well logging lithofacies can be developed as a prior distribution model; ② Based on the existing well logging data, the lithofacies division in step 102 is performed, and the mean and variance of each lithofacies reservoir parameter are calculated; ③ Within the target solution space, a random search is performed using the MH sampling method to obtain candidate seed points; ④ The MH judgment criterion is used to select these seed points, determine whether to accept the candidate point, and then calculate the parameter candidate point at the next location; ⑤ Repeat steps ③-④. When all parameter judgments are completed or the loop termination condition is met, the inversion lithofacies results are output.

[0103] In step 104, the accuracy of the algorithm is verified by comparing the resolution of the inversion result with the lithofacies classification result, and its final effect is analyzed to determine whether the algorithm is effective. The process ends.

[0104] The following are several specific embodiments of the present invention:

[0105] Example 1

[0106] In a specific embodiment 1 of the present invention, Figure 2 As shown in the figure, the application of this method was verified using actual well logging data. A representative training dataset was obtained by combining the existing lithofacies classification method using only P-wave velocity as a feature with the method used in this invention, which combines gamma rays with P-wave velocity. Inversion was performed using the Bayesian framework-based Markov Chain Monte Carlo random inversion method, and the results obtained from the two methods were compared. It can be seen that the inversion accuracy of both methods increases with the increase in the number of MCMC samplings, and the accuracy of this method is significantly superior to that of the traditional algorithm.

[0107] Example 2

[0108] In a specific embodiment 2 of the present invention, Figure 3 In this study, existing well logging data was classified into six lithofacies using this method, a method using only P-wave velocity, and a method using only gamma rays. The classification accuracy for each of the six lithofacies was calculated. As can be seen from the figure, the inversion results obtained using only gamma rays are the least accurate, and the overall performance of this method is significantly superior to the other two methods.

[0109] Example 3

[0110] In a specific embodiment 3 of the present invention, Figure 4 Maintaining the same conditions for combining gamma rays with P-wave velocities for classification as described in this method, we compared the MCMC method based on the Bayesian inversion framework described in this method with the traditional classification method using the complete probability distribution function (PDF). From left to right, we show the measured values, the traditional PDF method, and this method. The inversion result obtained by this method on the far right is the closest to the lithofacies derived from the measured data.

[0111] Example 4

[0112] In a specific embodiment 4 of the present invention, Figure 5 A comprehensive comparison was made of the average accuracy of the inversion results obtained by the traditional PDF method, the traditional MCMC method, and the logging lithofacies classification method that integrates the rock physics model and Bayesian inversion. As can be seen from the image, the PDF method has the worst accuracy, and the method of the present invention is the best.

[0113] Finally, it should be noted that the above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art may modify the technical solutions described in the aforementioned embodiments or substitute equivalents for some of the technical features therein. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention shall be included within the scope of protection of the present invention.

[0114] Except for the technical features described in the specification, all other technical features are known technologies to those skilled in the art.

Claims

1. A well logging lithofacies classification method integrating rock physics model and Bayesian inversion is characterized by: The well logging lithofacies classification method that integrates rock physics model and Bayesian inversion includes: Step 1: Analyze local geological and physical rock characteristics in detail to identify and classify lithology; Step 2: Perform petrographic rock physics analysis to obtain a statistically representative training database; Step 3, perform Markov chain Monte Carlo random inversion based on Bayesian inversion; Step 4, synthesizing the logging lithofacies data and the actual lithofacies data; In step 1, the well logging data are divided into seismic rock layer groups and sedimentary unit groups with shear-like structural characteristics. Through detailed analysis of local geological and physical rock characteristics, the understanding of the variability of the well logging data is improved, which guides the optimization of logging parameters and performs lithology identification and classification. In step 1, the well logging lithofacies are classified by lithology including sand, silt and clay, bedding structure including massive, interbedded or disordered, lithofacies including grain size, clay location and cementation, and seismic properties including compressional wave velocity, shear wave velocity and density; In step 2, the gamma ray and P-wave velocities of the logging data are analyzed. Combining physical and geological knowledge, key logging parameters for formation classification are selected to obtain lithofacies classification data for the logging, thus obtaining a representative training data set. In step 2, the well logging data is subjected to median filtering to remove high-frequency noise and smoothing, so that the thickness of the individual phases classified from the well logging data is more closely correlated with the seismic significance layer, thereby obtaining better clustering results. In step 2, a type of well is selected from the well logging data to identify seismic lithofacies. Since gamma rays are a good indicator of clay in quartz-rich sediments, gamma rays and P-wave velocity are combined to determine different lithofacies. These two parameters are used together as classification parameters to classify the lithofacies. In step 3, assuming that the reservoir parameters follow a Gaussian distribution, the existing Gibbs sampling and Metropolis-Hastings sampling algorithms (MH algorithms) are used to extract a Markov chain whose transition probability is only related to the previous state during state transition, thus completing the Markov process. In step 3, the MH sampling method is used to randomly search the target space, and the points obtained are selected using the judgment criteria. When the generated random number is less than the acceptance probability, the obtained point is recorded. Otherwise, re-sampling is performed, and the points are recorded. The parameter candidate points at the next position are comprehensively counted. Repeat the above steps. When all parameter judgments are completed or the iteration termination condition is reached, the inversion result value is output. In step 3, the specific steps of MH sampling are: (1) Given the state transition matrix P and the assumed stationary distribution π(x), set the number of state transitions n1 and the required number of sample points n2; (2) Given an initial distribution, sample from the distribution to obtain the initial state x0; (3) Assuming the iterative loop parameter is t, the maximum number of iterations is n1+n2-1, and t starts the loop at 0, perform the following steps: i. According to the conditional probability distribution P(x|x t ), from which sample point x is taken; ii. Sampling from the [0,1] distribution; i.e., position u; iii. If satisfied Then accept this sampling and update the Markov chain state. At this time, x t+1 =x; iv. If the above conditions are not met, the update will not be accepted. t+1 =x t ; In step 4, by comparing the resolution of the inversion results with the lithofacies classification results, the accuracy of the algorithm is verified, its final effect is analyzed, and a conclusion is drawn on whether the algorithm is effective.

2. The well logging lithofacies classification method integrating rock physics model and Bayesian inversion according to claim 1 is characterized in that: In step 1, the relationship between velocity and porosity is analyzed from the existing logging data. The core data is comprehensively analyzed based on the following knowledge of physical and rock properties to classify the lithofacies: (1) The cementation of rock phases produces a stiffening effect, resulting in relatively high seismic velocities and impedance; (2) Cemented sand bodies, layered sand bodies, and interlayer sandstones and shales also have relatively high velocities; (3) As the clay content in the rock increases, the density of the rock will increase; (4) Silty shale is denser than pure shale; (5) For particle-supported sediments, an increase in clay content reduces porosity, increases density, and hardens the rock; (6) For clay minerals, due to their inherent porosity, the porosity will increase with the increase of clay mineral content, and the rock skeleton will be weakened. Therefore, when the clay content of the rock is 40%, its wave propagation speed will reach its peak value; (7) Since clay particles are flaky, the shear strength of shale is relatively weak, so the ratio of P-wave propagation velocity to S-wave propagation velocity in shale is higher than that in sandstone; Based on the existing logging core data, the lithofacies can be divided into six types: gravel and conglomerate, thick sandstone, interlayer sandstone and shale, shale containing a large amount of silt, relatively pure shale, and disordered sedimentation.

3. The well logging lithofacies classification method integrating rock physics model and Bayesian inversion according to claim 1 is characterized in that: In step 3, for the facies data, a Markov variable is obtained. As the state changes, the transition probability of the variable is only related to the previous state and has nothing to do with other states. P(X t+1 =S j |X0=S0,X1=S1,…,X t =S t )=P(X t+1 =S j |X t =S i ) (1) X represents a random variable, X t is the state of variable X at a certain time t, S0, S1, ..., S i , S j ∈Ω is all possible cases of X, P(X t+1 =X j |X t =S i ) indicates that X t =S i Status changes to X t+1 =S j Probability of state; The expression of state transition probability is: P(i→j):=P i,j =P(X t+1 =S j |X t =S i ) (2) Set the variable X to S at a certain time t k The probability of Then the value of the variable at time t+1 is S i The probability of Assume that there are n cases of X, then: Eventually converges to a stationary distribution π (*) , that is: lim n→∞ π 0 P t =π (*) ,π (*) Satisfy π (*) P=π (*) ; Due to the constraints of the actual situation and the influence of other factors in the algorithm solution process, π(x) and P do not satisfy π(i)P i，j =π(j)P j，i , so that the final Markov chain cannot be stabilized, so an acceptance probability α is added i，j , so that the equation holds true, that is: π(i)α i，j P i，j =π(j)α j，i P j，i Acceptance probability α i，j The formula is: The acceptance probability is used to determine whether to perform the next update. Through continuous state updates and judgments, the target distribution eventually reaches a stable condition.

4. The well logging lithofacies classification method integrating rock physics model and Bayesian inversion according to claim 1, characterized in that: In step 3, the processing method for logging data is as follows: ① Based on the experience of geostatistics, an initial model of logging lithofacies can be given as a prior distribution model; ② The lithofacies division in step 2 is performed according to the existing logging data, and the mean and variance of each lithofacies reservoir parameter are counted respectively; ③ In the target solution space, the MH sampling method is used to randomly search the docking space to obtain candidate seed points; ④ The MH judgment criterion is used to select the above seed points to determine whether to accept the candidate point, so as to count the parameter candidate points at the next position; ⑤ Repeat the above steps ③-④; when all parameter judgments are completed or the loop termination condition is reached, the inversion lithofacies result is output.

Images