Seismic inversion reservoir porosity modeling method driven by phase-controlled stochastic simulation

Abstract

The present invention discloses a seismic inversion reservoir porosity modeling method driven by phase-controlled random simulation, the core of which is the lithology-elastic parameter calibration relationship, phase-separation probability density function, phase-controlled geostatistical random simulation and co-simulation, and seismic inversion. The wide coverage of seismic data is fully utilized to obtain a reservoir porosity model that is faithful to hard-constrained data such as drilling and logging, and faithful to the lateral distribution characteristics of seismic data. The method introduces phase-controlled geostatistical random simulation into seismic inversion reservoir porosity modeling. In the process of geostatistical random simulation, the phase-separation probability density function model is used to overcome the problems of low modeling accuracy and strong uncertainty caused by strong heterogeneity of underground strata, improve the accuracy of reservoir porosity modeling, and reduce the uncertainty of modeling.

Description

Technical Field

[0001] The present invention belongs to the field of petroleum geophysical exploration, and the specific technology is a seismic inversion reservoir porosity modeling method driven by phase-controlled random simulation. The method introduces phase-controlled geostatistical random simulation into seismic inversion reservoir porosity modeling. In the process of geostatistical random simulation, the phase probability density function is used to overcome the problems of low modeling accuracy and strong uncertainty caused by strong heterogeneity of underground strata, improve the accuracy of reservoir porosity modeling, and reduce the uncertainty of modeling. It provides stable and reliable data support for multiple links of the oil and gas industry, such as sweet spot prediction of oil and gas reservoirs, well site deployment, hydraulic fracturing design, and dynamic simulation of oil and gas reservoirs. Background Art

[0002] Oil and gas reservoir porosity modeling is an important part of oil and gas exploration and development. Its results can be used in many important links in oil and gas exploration and development, such as oil and gas sweet spot prediction, well site deployment, hydraulic fracturing design, and dynamic simulation of oil and gas reservoirs. The traditional reservoir porosity modeling method uses drilling, logging, and core data as hard constraint data, and seismic data as soft constraint data. It uses geostatistical estimation or random simulation algorithms to extrapolate the inter-well reservoir porosity, and obtains the distribution model of the reservoir porosity of the entire work area in the underground three-dimensional space, thus realizing reservoir porosity modeling.

[0003] With the development of the petroleum industry, the exploration of conventional oil and gas reservoirs has gradually deepened, the potential for oil and gas production value-added has become smaller and smaller, and the focus of exploration and development has gradually shifted from shallow and medium-layer oil and gas reservoirs to deep and ultra-deep oil and gas reservoirs. With the advancement of theory and technology, global oil and gas exploration has a trend of advancing into deep layers (4500-6000m) and ultra-deep layers (>6000m), and more and more oil and gas reservoirs have been discovered in deep layers. However, due to the characteristics of deep and ultra-deep reservoirs such as deep burial depth, complex formations, and strong heterogeneity, exploration and development are difficult. Traditional reservoir porosity modeling methods have problems such as low accuracy and large uncertainty, making it difficult to achieve satisfactory results.

[0004] For deep and ultra-deep reservoirs, the seismic inversion reservoir porosity modeling method driven by random simulation is the most appropriate. This method uses geostatistical random simulation to iteratively update elastic parameters, regenerate synthetic seismic data and match them with actual seismic data until the inversion objective function converges; then, based on the statistical relationship between elastic parameters and porosity, the elastic parameters obtained by seismic inversion are converted to porosity to achieve reservoir porosity modeling. However, this modeling method of first inverting elastic parameters and then obtaining reservoir porosity still has serious shortcomings. First, there may not be a consistent correlation between reservoir porosity and elastic parameters. The relationship between reservoir porosity and elastic parameters of different lithofacies may not be consistent. If elastic parameters (for example, seismic wave impedance) are directly converted to porosity, there may be extremely strong multiple solutions. In particular, for deep and ultra-deep reservoirs with strong heterogeneity, the relationship between elastic parameters and reservoir porosity is more complicated, which further limits the actual effect of porosity modeling. Secondly, the porosity model obtained by indirect conversion is not directly derived from the original seismic data. Therefore, its distribution characteristics may not be faithful to the original seismic data, and the accuracy and reliability of the reservoir porosity distributed between wells are not high. In addition, the existing seismic inversion reservoir porosity modeling methods have taken reservoir heterogeneity into account to a certain extent. For example, when stochastically simulating elastic parameters in geostatistics, the probability density function of the layered segment or the self-updating probability density function is used, but artificial stratification of the inverted geological segments still cannot reasonably characterize the strong heterogeneity of deep and ultra-deep reservoirs. Summary of the invention

[0005] The purpose of the present invention is to provide a phase-controlled random simulation driven seismic inversion reservoir porosity modeling method, which makes full use of the wide coverage of seismic data to obtain a reservoir porosity model that is faithful to hard constraint data such as drilling and logging, and is faithful to the lateral distribution characteristics of seismic data. The method introduces phase-controlled geostatistical random simulation into seismic inversion reservoir porosity modeling. In the process of geostatistical random simulation, a phase-separated probability density function is used to overcome the problems of low modeling accuracy and strong uncertainty caused by strong heterogeneity of underground strata, improve the accuracy of reservoir porosity modeling, and reduce the uncertainty of modeling.

[0006] In order to achieve the above technical objectives, the present invention provides the following technical solutions.

[0007] The core of the seismic inversion reservoir porosity modeling method driven by phase-controlled random simulation is the lithology-elastic parameter calibration relationship, phase separation probability density function, phase-controlled geostatistical random simulation and co-simulation, and seismic inversion. The method includes the following steps in sequence:

[0008] (1) Based on logging elastic parameters and lithofacies interpretation, establish the lithofacies-elastic parameter calibration relationship;

[0009] (2) Based on the logging porosity, elastic parameters and lithofacies interpretation, the phase-separated porosity probability density function and the phase-separated porosity-elastic parameter joint probability density function are established;

[0010] (3) Inputting initial elastic parameters, and obtaining initial reservoir lithofacies distribution based on the lithofacies-elastic parameter calibration relationship;

[0011] (4) Based on the well logging porosity, reservoir lithofacies distribution, and phase-separated porosity probability density function, multiple realizations of the reservoir porosity model are obtained by phase-controlled geostatistical stochastic simulation;

[0012] (5) Using the reservoir porosity model as an auxiliary variable, multiple realizations of the elastic parameters are obtained by phase-controlled geostatistical stochastic co-simulation based on the well logging elastic parameters, reservoir lithofacies distribution, and the joint probability density function of phase-separated porosity and elastic parameters;

[0013] (6) Perform seismic forward modeling on multiple realizations of elastic parameters to obtain corresponding synthetic seismic data;

[0014] (7) Inputting actual seismic data, calculating the seismic inversion objective function value, and selecting the best matching elastic parameters and the corresponding porosity model;

[0015] (8) Based on the lithofacies-elastic parameter calibration relationship, the updated reservoir lithofacies distribution is obtained by the best matching elastic parameters;

[0016] (9) Iterating the above processes (4) to (8) until the seismic inversion objective function value reaches a preset tolerance requirement;

[0017] (10) The best matching porosity model obtained in the last iteration is the final modeling result.

[0018] Specific:

[0019] Step (1) specifically includes the following sub-steps:

[0020] Firstly, the distribution histograms of wave impedance are statistically analyzed for different lithofacies, that is, the distribution histograms of wave impedance are statistically analyzed for limestone facies and dolomite facies respectively.

[0021] Then, the probability density function of the wave impedance distribution of different lithofacies is fitted with one-dimensional Gaussian distribution as the theoretical probability distribution, namely, the phase-separated wave impedance probability density function.

[0022] Next, the wave impedance value corresponding to the intersection point of the phase-separation wave impedance probability density function is used as the threshold for dividing different lithofacies, and the lithofacies-wave impedance calibration relationship shown below is obtained:

[0023]

[0024] Among them, c1 With c 2 represents dolomite facies and limestone facies, z is the wave impedance, and z* is the threshold for dividing different facies.

[0025] Step (2) specifically includes the following sub-steps:

[0026] First, the probability density function of phase-separated porosity is established. The distribution histogram of porosity of different lithofacies is statistically analyzed, and then the probability density function of porosity distribution of different lithofacies is fitted using one-dimensional Gaussian distribution as the theoretical probability distribution, which is the probability density function of phase-separated porosity.

[0027] Then, the joint probability density function of phase-separated porosity and wave impedance is established. The two-dimensional distribution histogram of porosity and wave impedance of different lithofacies is statistically analyzed, and then the two-dimensional Gaussian distribution is used as the theoretical joint probability distribution to fit the joint probability density function of porosity and wave impedance of different lithofacies, which is the joint probability density function of phase-separated porosity-wave impedance.

[0028] In step (3), the wave impedance inversion result obtained by deterministic seismic wave impedance inversion is used as the initial wave impedance.

[0029] Step (4) specifically includes the following sub-steps:

[0030] First, a random path is established in the entire grid to be simulated (i.e., in the three-dimensional space to be modeled). Based on the random path and the logging porosity data, the grid points to be simulated are selected, and the local mean and local variance of the porosity of the grid points to be simulated are obtained through Kriging estimation.

[0031] Then, based on the reservoir lithofacies distribution corresponding to the grid points to be simulated, the phase porosity probability density function of the corresponding lithofacies is selected. The following Gaussian probability density function is defined by the phase probability density function:

[0032]

[0033] Among them, G(a,b) is a Gaussian probability density function with a as mean and b as variance, z sk is the local mean of porosity estimated by Kriging, is the local variance of porosity estimated by Kriging, is the phase probability density function and The conversion process is:

[0034]

[0035] Among them, μ and σ are the mean and standard deviation of the phase separation probability density function respectively.

[0036] Afterwards, Random sampling is performed to obtain a random simulated value of porosity. The specific method of random sampling is to generate a random number p from a uniform distribution in the interval (0, 1), and then obtain a random number p from the Gaussian probability density function. Generate a random number y so that its corresponding probability value is p, that is:

[0037]

[0038] Then execute the random number y The reverse process

[0039]

[0040] That is, a random simulation value of porosity is obtained

[0041] Finally, the next grid point is simulated along the random path until the simulation of the entire grid is completed, and a random simulation realization of the porosity model is obtained. Multiple different random paths are established in the entire grid to be simulated, and multiple realizations of the porosity model can be obtained.

[0042] Step (5) specifically includes the following sub-steps:

[0043] First, a random path is established in the entire grid to be simulated. Based on the random path and the logging wave impedance, the grid points to be simulated are selected, and the reservoir porosity model is used as an auxiliary variable. The local mean and local variance of the wave impedance of the grid points to be simulated are obtained through cokriging estimation.

[0044] Then, based on the reservoir lithofacies distribution corresponding to the grid point to be simulated, the phase porosity-wave impedance joint probability density function of the corresponding lithofacies is selected. Then, the phase condition probability density function of the wave impedance of the grid point f(z 1 |z 2 ,c i ):

[0045]

[0046] Among them, c i is the corresponding lithofacies distribution of the simulated grid points, f(z 2 |c i ) is the probability density function of phase-separated porosity.

[0047] Then, the following Gaussian probability density function is defined by the wave impedance phase condition probability density function:

[0048]

[0049] G(a,b) is a Gaussian probability density function with a as mean and b as variance, but in this case, z csk is the local mean of the wave impedance estimated by cokriging, is the local variance of the wave impedance estimated by cokriging, The probability density function of the wave impedance phase separation condition f(z 1 |z 2 ,c i )and The conversion process,

[0050]

[0051] At this time, μ and σ are the mean and standard deviation of the probability density function of the wave impedance phase separation condition respectively.

[0052] Next, Perform random sampling to obtain a random simulated value of porosity. The specific method of random sampling is the same as step (4), that is, generate a random number y such that:

[0053]

[0054] Then execute the random number y The reverse process

[0055]

[0056] That is, a random co-simulation value of the wave impedance is obtained.

[0057] Finally, the next grid point is simulated along the random path until the simulation of the entire grid is completed, and a co-simulation realization of the wave impedance is obtained. Based on multiple different random paths established in the entire grid to be simulated, multiple realizations of the wave impedance can be obtained.

[0058] Step (6) specifically includes the following sub-steps:

[0059] First, the corresponding reflection coefficient sequence is calculated from the wave impedance:

[0060]

[0061] Among them, r(i) and z(i) represent the values ​​of the reflection coefficient sequence and the wave impedance at the i-th sampling point respectively.

[0062] Then, the synthetic seismic data is obtained by the seismic data convolution model:

[0063] d=w*r (12)

[0064] Among them, d is the synthetic seismic data, w is the seismic wavelet, and r is the reflection coefficient sequence.

[0065] Step (7) specifically includes the following sub-steps:

[0066] First, the value of the seismic inversion objective function is calculated based on the actual seismic data and the synthetic seismic data:

[0067] F=||sd|| 2 +λ|||zz prior || 2 (13)

[0068] Among them, F is the seismic inversion objective function, s is the actual seismic data, z is the wave impedance, and z prior is the prior wave impedance, λ is the prior wave impedance regularization constraint coefficient, which can usually be taken as 0.1, and ||.|| is the L2 norm. The prior wave impedance should contain effective low-frequency information and have lateral stability. As with the given initial wave impedance, the present invention uses the wave impedance inversion result obtained by deterministic seismic wave impedance inversion as the prior wave impedance.

[0069] Then, among the seismic inversion objective function values ​​obtained by multiple realizations of wave impedance, the wave impedance corresponding to the minimum seismic inversion objective function value is selected, which is the best matching wave impedance. The porosity model corresponding to the best matching wave impedance is the best porosity model.

[0070] Step (9) specifically includes the following sub-steps:

[0071] If the seismic inversion objective function value satisfies the following iteration tolerance condition, the seismic inversion iteration is stopped:

[0072]

[0073] Among them, l is the current iteration number, and ε is the iteration tolerance parameter. BRIEF DESCRIPTION OF THE DRAWINGS

[0074] Figure 1 is a flow chart of the present invention;

[0075] Figure 2 is the calibration relationship between lithofacies and wave impedance in the embodiment;

[0076] Figure 3 is the joint probability density function of phase porosity and wave impedance in the embodiment. DETAILED DESCRIPTION

[0077] The following is a further detailed description of the process of the present invention by taking the porosity modeling of carbonate reservoirs using seismic wave impedance inversion (ie, the elastic parameter is wave impedance, and the lithofacies are divided into two categories: limestone phase and dolomite phase) as an example.

[0078] The flowchart corresponding to the above steps is as follows Figure 1 shown.

[0079] Step (1): Input the well logging wave impedance data and lithofacies interpretation data to establish the lithofacies-wave impedance calibration relationship. The specific implementation process is as follows.

[0080] Firstly, the distribution histograms of wave impedance are statistically analyzed for different lithofacies, that is, the distribution histograms of wave impedance are statistically analyzed for limestone facies and dolomite facies respectively.

[0081] Then, the probability density function of the wave impedance distribution of different lithofacies is fitted with one-dimensional Gaussian distribution as the theoretical probability distribution, namely, the phase-separated wave impedance probability density function.

[0082] Next, the wave impedance value corresponding to the intersection point of the phase-separation wave impedance probability density function is used as the threshold for dividing different lithofacies, and the lithofacies-wave impedance calibration relationship shown below is obtained:

[0083]

[0084] Among them, c 1 With c 2 represents dolomite facies and limestone facies, z is the wave impedance, and z* is the threshold for dividing different facies.

[0085] Figure 2 The calibration relationship between lithofacies and wave impedance is shown. f(z|c 1 ) and f(z|c 2 ) is the fitted probability density function of the phase-separated wave impedance.

[0086] Step (2): Input the logging porosity, elastic parameters and lithofacies interpretation to construct the phase-separated porosity probability density function and the phase-separated porosity-wave impedance joint probability density function. The specific implementation process is as follows:

[0087] First, establish the phase-separated porosity probability density function. Similar to step (1), statistically analyze the porosity distribution histograms of different lithofacies, and then use the one-dimensional Gaussian distribution as the theoretical probability distribution to fit the probability density function of the porosity distribution of different lithofacies, which is the phase-separated porosity probability density function.

[0088] Then, the joint probability density function of phase-separated porosity and wave impedance is established. The two-dimensional distribution histogram of porosity and wave impedance of different lithofacies is statistically analyzed, and then the two-dimensional Gaussian distribution is used as the theoretical joint probability distribution to fit the joint probability density function of porosity and wave impedance of different lithofacies, which is the joint probability density function of phase-separated porosity-wave impedance.

[0089] Down Figure 3 The phase-separated porosity-wave impedance joint probability density function f(z|c1 ) and f(z|c 2 ) is the phase-separated porosity-wave impedance joint probability density function, where z = (z 1 ,z 2 ), at this time z 1 With z 2 represent wave impedance and porosity respectively.

[0090] Step (3): Input the initial wave impedance, and obtain the initial reservoir lithofacies distribution from the initial wave impedance based on the lithofacies-wave impedance calibration relationship. The present invention uses the wave impedance inversion result obtained by deterministic seismic wave impedance inversion as the initial wave impedance, and the deterministic seismic wave impedance inversion can adopt any mature inversion method, such as sparse pulse inversion, model inversion, color inversion, etc.

[0091] Step (4): Based on the well logging porosity, reservoir lithofacies distribution, and phase-separated porosity probability density function, multiple realizations of the reservoir porosity model are obtained by phase-controlled geostatistical stochastic simulation. The specific implementation process of phase-controlled geostatistical stochastic simulation is as follows:

[0092] First, a random path is established in the entire grid to be simulated (i.e., in the three-dimensional space to be modeled). Based on the random path and the logging porosity data, the grid points to be simulated are selected, and the local mean and local variance of the porosity of the grid points to be simulated are obtained through Kriging estimation.

[0093] Then, based on the reservoir lithofacies distribution corresponding to the grid points to be simulated, the phase porosity probability density function of the corresponding lithofacies is selected. The following Gaussian probability density function is defined by the phase probability density function:

[0094]

[0095] Among them, G(a,b) is a Gaussian probability density function with a as mean and b as variance, z sk is the local mean of porosity estimated by Kriging, is the local variance of porosity estimated by Kriging, is the phase probability density function and The conversion process is:

[0096]

[0097] Among them, μ and σ are the mean and standard deviation of the phase separation probability density function respectively.

[0098] Afterwards, Random sampling is performed to obtain a random simulated value of porosity. The specific method of random sampling is to generate a random number p from a uniform distribution in the interval (0, 1), and then obtain a random number p from the Gaussian probability density function. Generate a random number y so that its corresponding probability value is p, that is:

[0099]

[0100] Then execute the random number y The reverse process

[0101]

[0102] That is, a random simulation value of porosity is obtained

[0103] Finally, the next grid point is simulated along the random path until the simulation of the entire grid is completed, and a random simulation realization of the porosity model is obtained. Multiple different random paths are established in the entire grid to be simulated, and multiple realizations of the porosity model can be obtained.

[0104] Step (5): Using multiple realizations of the reservoir porosity model as auxiliary variables, based on the well logging impedance, reservoir lithofacies distribution, and phase-separated porosity-impedance joint probability density function, multiple realizations of the impedance are obtained by phase-controlled geostatistical stochastic co-simulation. Similar to phase-controlled geostatistical stochastic simulation, the specific implementation process of phase-controlled geostatistical stochastic co-simulation is as follows.

[0105] First, a random path is established in the entire grid to be simulated. Based on the random path and the logging wave impedance, the grid points to be simulated are selected, and the reservoir porosity model is used as an auxiliary variable. The local mean and local variance of the wave impedance of the grid points to be simulated are obtained through cokriging estimation.

[0106] Then, based on the reservoir lithofacies distribution corresponding to the grid point to be simulated, the phase porosity-wave impedance joint probability density function of the corresponding lithofacies is selected. Then, the phase condition probability density function of the wave impedance of the grid point f(z 1 |z 2 ,c i ):

[0107]

[0108] Among them, c i is the corresponding lithofacies distribution of the simulated grid points, f(z 2 |c i ) is the probability density function of phase-separated porosity.

[0109] Then, the following Gaussian probability density function is defined by the wave impedance phase condition probability density function:

[0110]

[0111] Similar to step (4), G(a,b) is a Gaussian probability density function with a as mean and b as variance, but in this case, z csk is the local mean of the wave impedance estimated by cokriging, is the local variance of the wave impedance estimated by cokriging, The probability density function of the wave impedance phase separation condition f(z 1 |z 2 ,c i )and The conversion process,

[0112]

[0113] At this time, μ and σ are the mean and standard deviation of the probability density function of the wave impedance phase separation condition respectively.

[0114] Next, Perform random sampling to obtain a random simulated value of porosity. The specific method of random sampling is the same as step (4), that is, generate a random number y such that:

[0115]

[0116] Then execute the random number y The reverse process

[0117]

[0118] That is, a random co-simulation value of the wave impedance is obtained.

[0119] Finally, the next grid point is simulated along the random path until the simulation of the entire grid is completed, and a co-simulation realization of the wave impedance is obtained. Based on multiple different random paths established in the entire grid to be simulated, multiple realizations of the wave impedance can be obtained.

[0120] Step (6): Perform seismic forward modeling on multiple wave impedance realizations to obtain corresponding synthetic seismic data. The specific implementation process is as follows.

[0121] First, the corresponding reflection coefficient sequence is calculated from the wave impedance:

[0122]

[0123] Among them, r(i) and z(i) represent the values ​​of the reflection coefficient sequence and the wave impedance at the i-th sampling point respectively.

[0124] Then, the synthetic seismic data is obtained by the seismic data convolution model:

[0125] d=w*r (12)

[0126] Among them, d is the synthetic seismic data, w is the seismic wavelet, and r is the reflection coefficient sequence.

[0127] Step (7): Input actual seismic data, calculate the seismic inversion objective function value, and select the best matching wave impedance and the corresponding porosity model. The specific implementation process is as follows.

[0128] First, the value of the seismic inversion objective function is calculated based on the actual seismic data and the synthetic seismic data:

[0129] F=||sd|| 2 +λ|||zz prior || 2 (13)

[0130] Among them, F is the seismic inversion objective function, s is the actual seismic data, z is the wave impedance, and z prior is the prior wave impedance, λ is the prior wave impedance regularization constraint coefficient, which can usually be taken as 0.1, and ||.|| is the L2 norm. The prior wave impedance should contain effective low-frequency information and have lateral stability. As with the given initial wave impedance, the present invention uses the wave impedance inversion result obtained by deterministic seismic wave impedance inversion as the prior wave impedance.

[0131] Then, among the seismic inversion objective function values ​​obtained by multiple realizations of wave impedance, the wave impedance corresponding to the minimum seismic inversion objective function value is selected, which is the best matching wave impedance. The porosity model corresponding to the best matching wave impedance is the best porosity model.

[0132] Step (8): Based on the lithofacies-wave impedance calibration relationship, the updated reservoir lithofacies distribution is obtained by the best matching wave impedance.

[0133] Step (9): Iterate the above steps (4) to (8) until the seismic inversion objective function value reaches the preset tolerance requirement. The specific implementation process is as follows.

[0134] If the seismic inversion objective function value satisfies the following iteration tolerance condition, the seismic inversion iteration is stopped:

[0135]

[0136] Among them, l is the current iteration number, and ε is the iteration tolerance parameter.

[0137] Step (10): Output the best matching porosity model obtained in the last iteration, i.e., the final modeling result.

Claims

1. Seismic inversion reservoir porosity modeling method driven by phase-controlled stochastic simulation, It is characterized in that The following steps are involved: (1) Based on logging elastic parameters and lithofacies interpretation, establish the lithofacies-elastic parameter calibration relationship; (2) Based on logging porosity, elastic parameters and lithofacies interpretation, the probability density function of phase-separated porosity and the joint probability density function of phase-separated porosity and elastic parameters are established; (3) Inputting initial elastic parameters and obtaining initial reservoir lithofacies distribution based on the lithofacies-elastic parameter calibration relationship; (4) Based on the well logging porosity, reservoir lithofacies distribution, and phase-separated porosity probability density function, multiple realizations of the reservoir porosity model are obtained by phase-controlled geostatistical stochastic simulation; (5) Using the reservoir porosity model as an auxiliary variable, multiple realizations of the elastic parameters are obtained by phase-controlled geostatistical stochastic co-simulation based on the well logging elastic parameters, reservoir lithofacies distribution, and the joint probability density function of phase-separated porosity and elastic parameters; (6) Perform seismic forward modeling on multiple realizations of elastic parameters to obtain corresponding synthetic seismic data; (7) Input actual seismic data, calculate the seismic inversion objective function value, and select the best matching elastic parameters and corresponding porosity model; The value of the seismic inversion objective function: F=||s-d|| 2 +λ||z-z prior || 2 (13) Among them, F is the seismic inversion objective function, s is the actual seismic data, d is the synthetic seismic data, z is the wave impedance, z prior is the prior wave impedance, λ is the prior wave impedance regularization constraint coefficient, ||.|| is the L2 norm; (8) Based on the lithofacies-elastic parameter calibration relationship, the updated reservoir lithofacies distribution is obtained by the best matching elastic parameters; (9) Iterating the above processes (4) to (8) until the seismic inversion objective function value reaches a preset tolerance requirement; (10) The best matching hole obtained in the last iteration.

2. The phase-controlled random simulation driven seismic inversion reservoir porosity modeling method according to claim 1, It is characterized in that Step (1) specifically includes the following sub-steps: Firstly, the wave impedance distribution histograms are statistically analyzed for different lithofacies, that is, the wave impedance distribution histograms are statistically analyzed for limestone phase and dolomite phase respectively; Then, the probability density function of the wave impedance distribution of different lithofacies is fitted using one-dimensional Gaussian distribution as the theoretical probability distribution, namely, the phase-separated wave impedance probability density function; Next, the wave impedance value corresponding to the intersection point of the phase-separation wave impedance probability density function is used as the threshold for dividing different lithofacies, and the lithofacies-wave impedance calibration relationship shown below is obtained: Among them, c 1 With c 2 represents the dolomite phase and limestone phase, z is the wave impedance, z * It is the threshold for dividing different lithofacies.

3. The phase-controlled random simulation driven seismic inversion reservoir porosity modeling method according to claim 2, It is characterized in that Step (2) specifically includes the following sub-steps: Firstly, the probability density function of porosity by phase is established; the distribution histogram of porosity by different lithofacies is statistically analyzed, and then the probability density function of porosity distribution by different lithofacies is fitted using one-dimensional Gaussian distribution as the theoretical probability distribution, which is the probability density function of porosity by phase; Then, the joint probability density function of phase-separated porosity and wave impedance is established; the two-dimensional distribution histograms of porosity and wave impedance are statistically analyzed for different lithofacies, and then the two-dimensional Gaussian distribution is used as the theoretical joint probability distribution to fit the joint probability density function of porosity and wave impedance of different lithofacies, which is the joint probability density function of phase-separated porosity-wave impedance.

4. The phase-controlled random simulation driven seismic inversion reservoir porosity modeling method according to claim 3, It is characterized in that In step (3), the wave impedance inversion result obtained by deterministic seismic wave impedance inversion is used as the initial wave impedance.

5. The phase-controlled random simulation driven seismic inversion reservoir porosity modeling method according to claim 4, It is characterized in that Step (4) specifically includes the following sub-steps: First, a random path is established in the entire grid to be simulated; based on the random path and the logging porosity data, the grid points to be simulated are selected, and the local mean and local variance of the porosity of the grid points to be simulated are obtained through Kriging estimation; Then, based on the reservoir lithofacies distribution corresponding to the grid points to be simulated, the phase porosity probability density function of the corresponding lithofacies is selected; and the following Gaussian probability density function is defined by the phase probability density function: Among them, G(a,b) is a Gaussian probability density function with a as mean and b as variance, z sk is the local mean of porosity estimated by Kriging, is the local variance of porosity estimated by Kriging, is the phase probability density function and The conversion process is: Among them, μ and σ are the mean and standard deviation of the phase separation probability density function respectively; Afterwards, Random sampling is performed to obtain a random simulation value of porosity; the specific method of random sampling is to generate a random number p from a uniform distribution in the interval (0, 1), and then obtain a random number p from a Gaussian probability density function. Generate a random number y so that its corresponding probability value is p, that is: Then execute the random number y The reverse process That is, a random simulation value of porosity is obtained Finally, the next grid point is simulated along the random path until the simulation of the entire grid is completed, that is, a random simulation realization of the porosity model is obtained; multiple different random paths are established in the entire grid to be simulated, and multiple realizations of the porosity model can be obtained.

6. The phase-controlled random simulation driven seismic inversion reservoir porosity modeling method according to claim 5, It is characterized in that Step (5) specifically includes the following sub-steps: First, a random path is established in the entire grid to be simulated; based on the random path and the logging wave impedance, the grid points to be simulated are selected, and the reservoir porosity model is used as an auxiliary variable. The local mean and local variance of the wave impedance of the grid points to be simulated are obtained through cokriging estimation; Then, based on the reservoir lithofacies distribution corresponding to the grid point to be simulated, the phase-separated porosity-wave impedance joint probability density function of the corresponding lithofacies is selected; then, the phase-separated porosity-wave impedance joint probability density function, the phase-separated porosity probability density function and the porosity value of the grid point to be simulated are used to obtain the phase-separated conditional probability density function of the wave impedance of the grid point f(z 1 |z 2 ,c i ): Among them, c i is the corresponding lithofacies distribution of the simulated grid points, f(z 2 |c i ) is the probability density function of phase-separated porosity; Then, the following Gaussian probability density function is defined by the wave impedance phase condition probability density function: G(a,b) is a Gaussian probability density function with a as mean and b as variance, but in this case, z csk is the local mean of the wave impedance estimated by cokriging, is the local variance of the wave impedance estimated by cokriging, The probability density function of the wave impedance phase separation condition f(z 1 |z 2 ,c i )and The conversion process, At this time, μ and σ are the mean and standard deviation of the probability density function of the wave impedance phase separation condition respectively; Next, Perform random sampling to obtain a random simulated value of porosity; random sampling generates a random number y such that: Then execute the random number y The reverse process That is, a random co-simulation value of the wave impedance is obtained; Finally, the next grid point is simulated along the random path until the simulation of the entire grid is completed, that is, a co-simulation realization of the wave impedance is obtained; based on multiple different random paths established in the entire grid to be simulated, multiple realizations of the wave impedance are obtained.

7. The phase-controlled random simulation driven seismic inversion reservoir porosity modeling method according to claim 6, It is characterized in that Step (6) specifically includes the following sub-steps: First, the corresponding reflection coefficient sequence is calculated from the wave impedance: Among them, r(i) and z(i) represent the values ​​of the reflection coefficient sequence and the wave impedance at the i-th sampling point respectively; Then, the synthetic seismic data is obtained by the seismic data convolution model: d=w*r (12) Among them, d is the synthetic seismic data, w is the seismic wavelet, and r is the reflection coefficient sequence.

8. The phase-controlled random simulation driven seismic inversion reservoir porosity modeling method according to claim 7, It is characterized in that Step (7) specifically includes the following sub-steps: First, the value of the seismic inversion objective function is calculated from the actual seismic data and the synthetic seismic data; the prior wave impedance should contain effective low-frequency information and have lateral stability; the wave impedance inversion result obtained by deterministic seismic wave impedance inversion is used as the prior wave impedance; Then, among the seismic inversion objective function values ​​obtained by multiple realizations of wave impedance, the wave impedance corresponding to the minimum seismic inversion objective function value is selected, which is the best matching wave impedance; and the porosity model corresponding to the best matching wave impedance is co-simulated, which is the best porosity model.

9. The phase-controlled random simulation driven seismic inversion reservoir porosity modeling method according to claim 8, It is characterized in that Step (9) specifically includes the following sub-steps: If the seismic inversion objective function value satisfies the following iteration tolerance condition, the seismic inversion iteration is stopped: Among them, l is the current iteration number, and ε is the iteration tolerance parameter.

Images