A seismic petrophysical nonlinear inversion method for connected porosity content

Abstract

The application discloses a seismic rock physics nonlinear inversion method for connected porosity content, and comprises the following steps: calculating elastic parameters of rock skeleton, pore fluid and fluid-containing rock, and constructing a seismic rock physics model; obtaining total porosity and connected porosity content near a well based on logging data; constructing a target parameter sample set according to the total porosity and the connected porosity content near the well; estimating a prior probability density function based on the target parameter sample set, and constructing a posterior probability density function; updating a target parameter model, calculating a receiving probability, and judging whether the target parameter model is received; repeating the foregoing steps until a maximum iteration number is reached, outputting the target parameter model, calculating an expected mean value of a connected porosity content model, and taking the expected mean value as a final inversion result. The application can effectively improve the seismic quantitative prediction accuracy of unconventional reservoirs such as carbonate rock and tight sandstone reservoirs.

Description

Technical Field

[0001] This invention belongs to the field of seismic exploration technology for unconventional reservoirs, and specifically relates to a seismic rock physics nonlinear inversion method for the content of interconnected porosity. Background Technology

[0002] Seismic rock physics inversion combines rock physics theory with seismic inversion techniques, aiming to quantitatively predict reservoir physical parameters using seismic data. It is an important technical means for oil and gas exploration and development. With the continuous advancement of oil and gas exploration, conventional reservoirs are becoming increasingly depleted, and unconventional reservoirs such as carbonate rocks and tight sandstones have become the main exploration targets. However, these reservoirs generally have complex pore structures, and conventional seismic rock physics inversion methods do not take into account the complexity of the reservoir pore structure, making it difficult to quantitatively predict pore connectivity, which is closely related to permeability, thus limiting the seismic exploration and development of such reservoirs.

[0003] Conventional seismic rock physics inversion methods mainly use reservoir physical parameters (porosity, clay content, water saturation, etc.) as variables to be predicted, such as the existing technologies disclosed in references [1]-[4]. However, the above methods cannot quantitatively predict and evaluate reservoir pore structure information. Studies have found that unconventional reservoirs such as carbonate rocks and tight sandstones usually develop complex pore structures, and their pore spaces usually contain interconnected pores and isolated pores. The relevant pore parameters are of great significance for reservoir permeability evaluation and have guiding value for reservoir oil and gas development (see references [5]-[7] for details). In response, seismic rock physics inversion techniques that consider pore structure have received increasing attention. For example, Zhao et al. (2013) proposed a seismic prediction method for pore type in carbonate rocks (see reference [8]); Li et al. (2022) proposed a seismic inversion method for predicting the percentage of hard pores in tight sandstone (see reference [9]); and Li Hongbing et al. (2021) developed a seismic inversion method for predicting the aspect ratio of pores in complex reservoirs (see reference

[10] ). Although the above inversion techniques can quantitatively predict pore type or structural parameters, they do not consider pore connectivity in rock physics modeling and cannot quantitatively assess pore connectivity. In addition, in order to quantitatively predict the content of connected pores in reservoirs, it is necessary to establish a seismic rock physics model that combines connected pores and isolated pores. Its forward modeling operator has a high degree of nonlinearity, so the corresponding inversion problem has a strong multiple solution and analytical solutions are usually difficult to obtain. Conventional quasi-linear or deterministic inversion algorithms are difficult to apply (see references

[11] -

[12] ). It is necessary to use nonlinear optimization algorithms to solve the inverse problem and rely on reliable prior information to constrain the inversion process to improve the accuracy of the prediction results.

[0004] In summary, for carbonate and tight sandstone reservoirs with complex pore structures, it is necessary to consider pore connectivity during rock physics modeling and inversion, and to develop effective nonlinear inversion algorithms to improve the accuracy of inversion results. Therefore, developing seismic rock physics nonlinear inversion methods for connected pore content is of great significance for improving the accuracy of seismic quantitative prediction of such unconventional reservoirs.

[0005] The above-mentioned non-patent literature is cited as follows:

[0006] [1]Fjeldstad, T., Avseth, P., Omre, H., 2021, A one-stepBayesian inversionframework for 3D reservoir characterization based on a Gaussian mixturemodel-A Norwegian Sea demonstration[J].Geophysics, 86(2), R221-R236.

[0007] [2] Luo, C., Ba, J., Guo, Q., 2023, Probabilistic seismic petrophysical inversion with statistical double-porosity Biot-Rayleigh model[J]. Geophysics, 88(3), M 157-M171.

[0008] [3] Gui Jinyong, Gao Jianhu, Yong Xueshan, Li Shengjun, 2015, Reservoir parameter inversion method based on two-phase medium theory [J]. Chinese Journal of Geophysics, 9, 3424-3438.

[0009] [4] Li Kun, Yin Xingyao, Zong Zhaoyun. 2020. Rock physics driven phase-constrained pre-stack seismic probabilistic inversion method [J]. Science in China: Earth Science, 50, 124-146.

[0010] [5] Deng Jixin, Zhou Hao, Wang Huan, et al. 2015. Elastic wave dispersion response analysis based on microscopic pore structure characteristics of reservoir sandstone [J]. Chinese Journal of Geophysics, 58, 3389-3400.

[0011] [6] He Xilei, Lin Kai, Zhang Zuhao et al. 2018. Quantitative comparison of the effects of porosity and pore structure on reservoir characteristics [J]. Petroleum Geophysical Exploration, 57(2): 179-185.

[0012] [7] Lai, J., Wang, G., Wang, Z., et al., 2018, A review on pore structure characterization in tight Sandstones[J]. Earth-Science Reviews, 177, 436-457.

[0013] [8] Zhao, L., Nasser, M., Han, DH, 2013, Quantitative geophysical pore-type characterization and its geological implication in carbonate reservoirs[J]. Geophysical Prospecting, 61, 827-841.

[0014] [9]Li, K., Yin,

[0015]

[10] Li Hongbing, Zhang Jiajia, Pan Haojie, et al., 2021, Nonlinear synchronous inversion of pore structure and physical property parameters based on elastic impedance [J]. Science in China: Earth Sciences, 51, 1166-1180.

[0016]

[11] Liu, Q., Dong, N., Ji, Y., et al., 2018, Direct reservoir property estimation based on prestack seismic inversion[J]. Journal of Petroleum Science and Engineering, 171, 1475-1486.

[0017]

[12] Zhang Jiajia, Yin Xingyao, Zhang Guangzhi, et al. 2020. Prediction method of physical property parameters based on linearized rock physics inversion [J]. Petroleum Exploration and Development, 47(1), 57-64. Summary of the Invention

[0018] Technical Problems Solved: This invention proposes a seismic rock physics nonlinear inversion method for connected pore content, addressing the shortcomings of existing technologies. These include the prevalence of complex pore structures in carbonate and tight sandstone reservoirs, the failure of conventional seismic rock physics inversion methods to consider the complexity of reservoir pore structures, the use of reservoir physical parameters (porosity, clay content, water saturation, etc.) as predictable variables leading to insufficient quantitative prediction and evaluation of pore structure information, the lack of consideration for pore connectivity in rock physics modeling resulting in insufficient quantitative prediction and evaluation of pore connectivity, and the high nonlinearity of seismic rock physics models for both connected and isolated pores, leading to multiple solutions to the inversion problem, difficulty in obtaining analytical solutions, and the inapplicability of conventional quasi-linear or deterministic inversion algorithms. The aim is to improve the accuracy of seismic quantitative prediction for unconventional reservoirs such as carbonate and tight sandstone reservoirs.

[0019] Technical solution:

[0020] This invention discloses a seismic rock physics nonlinear inversion method for determining the interconnected porosity content. The seismic rock physics nonlinear inversion method includes the following steps:

[0021] Step A: Calculate the rock skeleton, pore fluid, and elastic parameters of fluid-containing rocks to construct an earthquake rock physics model;

[0022] Step B: Based on well logging data and seismic rock physics model, obtain the total porosity and connected pore content near the well.

[0023] Step C: Based on the total porosity and connected pore content near the well, construct a Monte Carlo simulation sample set for the target parameters;

[0024] Step D: Based on the Monte Carlo simulation sample set of the target parameters, estimate the prior probability density function and construct the posterior probability density function of the target parameters;

[0025] Step E involves iteratively updating the target parameters. In each iteration, an updated model for the target parameters is randomly generated based on the probability density function of the target parameters, and the model acceptance probability is calculated. The model acceptance probability is then used to determine whether to accept the updated model.

[0026] Step F is repeated, followed by step E, until the maximum number of iterations is reached. The target parameter model is output, and the expected mean of the connected porosity content model is calculated as the final inversion result.

[0027] As a preferred embodiment of the present invention, step A, which involves calculating the rock skeleton, pore fluid, and elastic parameters of fluid-containing rocks to construct an earthquake rock physics model, includes the following sub-steps:

[0028] Step A1: Set basic rock physical parameters, including the bulk modulus, shear modulus, and density of minerals, and the bulk modulus and density of fluids; determine the percentage content of rock minerals based on core X-ray spectroscopy analysis, and calculate the bulk modulus, shear modulus, and density of the rock matrix using the VRH model.

[0029]

[0030]

[0031]

[0032] In the formula, K m G m and ρ m These represent the bulk modulus, shear modulus, and density of the rock matrix, respectively, K. mi G mi ρ mi and v mi These represent the bulk modulus, shear modulus, density, and percentage content of the i-th mineral, respectively, where N is the total number of mineral types.

[0033] Step A2: Divide the rock pores into connected pores and isolated pores, and calculate the equivalent elastic modulus of the rock skeleton using a self-consistent model:

[0034] (1-φ)(K m -K * )P *m +φv cn (K cn -K * )P *cn +φv ds (K ds -K * )P *ds =0 (4);

[0035] (1-φ)(G m -G * )Q *m +φv cn (G cn -G * )Q *cn +φv ds (G ds -G * )Q *ds =0 (5);

[0036] In the formula, K * and G * These represent the equivalent bulk modulus and shear modulus of the rock skeleton, respectively; φ is the total porosity near the well; vcn and v ds These represent the content of connected pores and the content of isolated pores, respectively, satisfying v cn +v ds =1; K cn and G cn K represents the bulk modulus and shear modulus of the interconnected pore inclusions, respectively. ds and G ds These are the bulk modulus and shear modulus of an isolated porous inclusion, respectively; P *cn and Q *cn P represents the geometric factor of the connected pore inclusions. *ds and Q *ds P represents the geometric factor of an isolated porous inclusion. *m and Q *m These are the geometric factors of the rock matrix;

[0037] Step A3: If connected pores are defined as empty pores and isolated pores as fluid-containing pores, then:

[0038] K cn =G cn =G ds =0 (6);

[0039] K ds =K fl (7);

[0040] In the formula, K fl The bulk modulus of the fluid;

[0041] Step A4: Based on the rock skeleton containing isolated pores, the Biot-Gassmann model is used to calculate the bulk modulus and shear modulus of the reservoir rock for connected pores.

[0042]

[0043] G rk =G * (9);

[0044] In the formula, K rk and G rk These are the bulk modulus and shear modulus of the reservoir rock, respectively, K. sk Bulk modulus of a rock skeleton containing isolated pores:

[0045]

[0046] Step A5, the P-wave velocity, S-wave velocity, and density of the reservoir rock are expressed as follows:

[0047]

[0048]

[0049] ρ rk =(1-φ)ρ m +φρ fl (13);

[0050] In the formula, V P V S and ρ rk These represent the P-wave velocity, S-wave velocity, and density of the reservoir rock, respectively, ρ fl The density of the fluid.

[0051] As a preferred embodiment of the present invention, step B, which involves obtaining the total porosity and connected porosity content near the well based on logging data, includes the following sub-steps:

[0052] Step B1: Based on the P-wave impedance and S-wave impedance data from well logging observations, the total porosity and connected porosity content near the well are obtained by inverting the predicted data from the matching model with the well logging observation data, and an error function is constructed:

[0053]

[0054] in,

[0055]

[0056]

[0057] In the formula, v j =[φ j v cnj ] represents the target parameter, φ j v represents the total porosity near the well at the j-th logging sampling point. cnj The content of connected porosity at the j-th well logging sampling point; and These are the observed P-wave impedance and S-wave impedance at the j-th logging sampling point, respectively. and V represents the predicted P-wave impedance and S-wave impedance at the j-th logging sampling point, respectively. P V S and ρ rk These represent the P-wave velocity, S-wave velocity, and density of the reservoir rock, respectively, K. rk and G rk Let J represent the bulk modulus and shear modulus of the reservoir rock, respectively, j = 1, 2, ..., N. s N s κ represents the total number of logging sampling points, and κ is the weighting coefficient of the error function;

[0058] Step B2, to balance the magnitude difference between the longitudinal wave impedance and the transverse wave impedance, estimates the weighting coefficients of the error function in the following way:

[0059]

[0060] Step B3: Using the fast simulated annealing algorithm, estimate the target parameters that minimize the error function, and obtain the total porosity and connected porosity v at the j-th logging sampling point. j =[φ j v cnj ].

[0061] As a preferred embodiment of the present invention, step C, the process of constructing the target parameter sample set based on the total porosity and connected pore content near the well, includes the following sub-steps:

[0062] Step C1: Based on the total porosity and connected pore content near the well, a Gaussian mixture model of the target parameters is constructed. The expectation-maximization algorithm is used to estimate the parameters of the Gaussian mixture model. The probability distribution function of the target parameters is then:

[0063]

[0064] In the formula, v = [φ, v cn [] represents the target parameter, φ represents the total porosity near the well, and v cn For the content of interconnected pores, N k Let μ represent the k-th Gaussian distribution. v|k、 ∑ v|k and σ k Let N be the expectation, covariance, and weight coefficients of the k-th Gaussian distribution. c The total number of Gaussian distributions;

[0065] Step C2: Based on the probability distribution function of the target parameters, Monte Carlo random simulation is used to obtain the augmented target parameters. These augmented target parameters are then combined with the total porosity and connected porosity content near the well to construct a Monte Carlo simulation sample set v for the target parameters. mci* =[φ i* v cni* ], i* = 1, 2, 3, ..., N mc , where N mc This represents the total number of samples in Monte Carlo.

[0066] As a preferred embodiment of the present invention, step D, which involves estimating the prior probability density function and constructing the posterior probability density function based on the target parameter sample set, includes the following sub-steps:

[0067] Step D1: Based on the target parameter sample set, calculate the prior probability density function of the target parameters using the Epanechnikov kernel density estimation method.

[0068]

[0069] in,

[0070]

[0071] In the formula, P(v) is the prior probability density function of the target parameter, v = [φ, v cn [] represents the target parameter, φ represents the total porosity near the well, and v cn For the content of interconnected pores, N mc Let L be the total number of Monte Carlo samples, L be the dimension of the target parameters, f(x) be the Epanechnikov kernel function, and h be the number of samples. v The scaling factor for the kernel function;

[0072] Step D2, based on the Bayesian framework, the posterior probability density function of the target parameters is equivalent to:

[0073]

[0074] In the formula, Z represents the observed impedance data, including the observed P-wave impedance and S-wave impedance, which are obtained from pre-stack seismic inversion; P(v|Z) and P(Z|v) are the posterior probability density function and likelihood function of the target parameter, respectively; P(Z) is the marginal probability of the observed impedance data, which is a constant.

[0075] Step D3, assuming the observed impedance data error follows a Gaussian distribution with zero expectation, then the likelihood function of the target parameter is expressed as:

[0076]

[0077] In the formula, D is the dimension of the observed impedance data. Combining equations (19), (21), and (22), the posterior probability density function of the target parameter is expressed as:

[0078]

[0079] In the formula, Z pdt The predicted impedance data, including predicted P-wave impedance and predicted S-wave impedance, are obtained from the rock physics model constructed in step A; C nor Let ∑ be the normalization constant. z This is the covariance matrix of the error in the observed impedance data.

[0080] As a preferred embodiment of the present invention, step E, the process of iteratively updating the target parameters, includes the following sub-steps:

[0081] Based on the prior probability density function of the target parameters, an updated model of the target parameters is randomly generated. The model acceptance probability is calculated based on the Metropolis criterion. The model acceptance probability at the k-th iteration is expressed as:

[0082]

[0083] Among them, v * For the update model of the objective parameters, v (k) Let P be the target parameter before the k-th iteration update. acp Let P(v) be the model acceptance probability. * |Z) and P(v) (k) |Z) represents the posterior probability density function of the target parameter update model and the posterior probability density function of the target parameter before the k-th iteration update, respectively; Z represents the observed impedance data, including the observed P-wave impedance and S-wave impedance.

[0084] Generate a random number between 0 and 1. If the random number is greater than or equal to P... acp Then, the updated model that receives the target parameters is given, and the target parameters v before the (k+1)th iteration are set. (k+1) =v * Otherwise, the objective parameter model remains unchanged, and the objective parameter v before the (k+1)th iteration update is set to... (k+1) =v (k) .

[0085] As a preferred embodiment of the present invention, in step F, the maximum number of iterations k max It is a positive integer greater than or equal to 20000.

[0086] As a preferred embodiment of the present invention, step F, which involves outputting the target parameter model and calculating the expected mean value of the connected porosity content model as the final inversion result, includes the following sub-steps:

[0087] Select the target parameter model set v for the last M iterations. est =[v w v w+1 v w+2 , ..., v w+M-1 ], v w Let w be the target parameter before the w-th iteration update, where w = k. max -M-1, based on the objective parameter model set v est Estimate the expected mean of the connected pore content model as the final inversion result.

[0088] Beneficial effects:

[0089] First, the seismic rock physics nonlinear inversion method for connected porosity content of the present invention is used to invert the connected porosity content of complex porosity reservoirs. The inversion results can be used to evaluate reservoir permeability and improve the accuracy of seismic quantitative prediction of unconventional reservoirs such as tight sandstone and carbonate rocks.

[0090] Second, the seismic rock physics nonlinear inversion method for connected pore content of the present invention constructs a seismic rock physics model that considers both connected and isolated pores, which can be used to invert connected pore content; it develops a seismic nonlinear inversion method and combines it with a stable inversion process based on the prior probability density function of the target parameter. Compared with traditional seismic rock physics inversion methods, the present invention improves the accuracy of quantitative seismic prediction of complex porous reservoirs. Attached Figure Description

[0091] Figure 1 This is a flowchart of the seismic rock physics nonlinear inversion method for connected porosity content according to an embodiment of the present invention;

[0092] Figure 2 The P-wave impedance and S-wave impedance data obtained from well logging observations in this application embodiment;

[0093] Figure 3 The total porosity and connected pore content near the well are obtained based on well logging data in the embodiments of this application;

[0094] Figure 4 The following are the inversion results and impedance data prediction results using conventional methods in the embodiments of this application: where (a) is the total porosity inversion result, and (b) is the longitudinal wave impedance and transverse wave impedance data prediction results;

[0095] Figure 5 The following are the inversion results and impedance data prediction results using the method of the present invention in the embodiments of this application: (a) is the inversion result of total porosity and connected pore content, and (b) is the prediction result of longitudinal wave impedance and transverse wave impedance data. Detailed Implementation

[0096] The following embodiments are provided to enable those skilled in the art to more fully understand the present invention, but do not limit the invention in any way.

[0097] The technical solution of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0098] This invention provides a seismic rock physics nonlinear inversion method for determining the interconnected porosity content, the seismic rock physics nonlinear inversion method comprising the following steps:

[0099] Step A: Calculate the rock skeleton, pore fluid, and elastic parameters of fluid-containing rocks to construct an earthquake rock physics model;

[0100] Step B: Based on well logging data and seismic rock physics model, obtain the total porosity and connected pore content near the well.

[0101] Step C: Based on the total porosity and connected pore content near the well, construct a Monte Carlo simulation sample set for the target parameters;

[0102] Step D: Based on the Monte Carlo simulation sample set of the target parameters, estimate the prior probability density function and construct the posterior probability density function of the target parameters;

[0103] Step E involves iteratively updating the target parameters. In each iteration, an updated model for the target parameters is randomly generated based on the probability density function of the target parameters, and the model acceptance probability is calculated. The model acceptance probability is then used to determine whether to accept the updated model.

[0104] Step F is repeated, followed by step E, until the maximum number of iterations is reached. The target parameter model is output, and the expected mean of the connected porosity content model is calculated as the final inversion result.

[0105] like Figure 1 The seismic rock physics nonlinear inversion method described in this embodiment specifically includes the following steps:

[0106] I. Constructing a seismic rock physics model and obtaining total porosity and connected porosity content near the well based on well logging data:

[0107] Basic rock physical parameters are set, including the bulk modulus, shear modulus, and density of minerals, and the bulk modulus and density of fluids; based on core X-ray spectroscopy analysis, the percentage content of rock minerals is determined, and the bulk modulus, shear modulus, and density of the rock matrix are calculated using the VRH model.

[0108]

[0109]

[0110]

[0111] In the formula, K mi G m and ρ m These represent the bulk modulus, shear modulus, and density of the rock matrix, respectively, K. mi G mi ρ mi and v miThese represent the bulk modulus, shear modulus, density, and percentage content of the i-th mineral, respectively, while N represents the total number of mineral types.

[0112] The rock pores are divided into connected pores and isolated pores, and the equivalent elastic modulus of the rock skeleton is calculated using a self-consistent model:

[0113] (1-φ)(K m -K * )P *m +φv cn (K cn -K * )P *cn +φv ds (K ds -K * )P *ds =0 (4);

[0114] (1-φ)(G m -G * )Q *m +φv cn (G cn -G * )Q *cn +φv ds (G ds -G * )Q *ds =0 (5);

[0115] In the formula, K * and G * Here, φ represents the equivalent bulk modulus and shear modulus of the rock skeleton, respectively, φ represents the total porosity, and v represents the total porosity. cn and v ds These represent the content of connected pores and the content of isolated pores, respectively, satisfying v cn +v ds =1, P* and Q* are the geometric factors of the inclusion, and K is the geometric factor of the inclusion. cn and G cn K represents the bulk modulus and shear modulus of the interconnected pore inclusions, respectively. ds and G ds These are the bulk modulus and shear modulus of an isolated porous inclusion, respectively.

[0116] If we define connected pores as empty pores and isolated pores as fluid-containing pores, then:

[0117] K cn =G cn =G ds =0 (6);

[0118] K ds =K fl (7);

[0119] In the formula, K fl This is the bulk modulus of the fluid.

[0120] Based on a rock skeleton containing isolated pores, the Biot-Gassmann model is used to calculate the bulk modulus and shear modulus of the reservoir rock for connected pores:

[0121]

[0122] G rk =G * (9);

[0123] In the formula, K rk and G rk These are the bulk modulus and shear modulus of the reservoir rock, respectively, K. sk Bulk modulus of a rock skeleton containing isolated pores:

[0124]

[0125] The P-wave velocity, S-wave velocity, and density of the reservoir rock can be expressed as follows:

[0126]

[0127]

[0128] ρ rk =(1-φ)ρ m +φρ fl (13);

[0129] In the formula, V P V S and ρ rk These represent the P-wave velocity, S-wave velocity, and density of the reservoir rock, respectively, ρ fl Where K is the density of the fluid; in this embodiment, the rock physical parameters used are as follows: K m =67 GPa, μ m =43 GPa, K fl =0.05GPa, ρ m =2.66g / cm 3 , ρ fl =0.80g / cm 3 .

[0130] Based on P-wave and S-wave impedance data from well logging observations, the total porosity and connected porosity content near the well are obtained by inversion using a matched model prediction and well logging observation data, and an error function is constructed:

[0131]

[0132] in,

[0133]

[0134]

[0135] In the formula, v j =[φ j v cnj ] represents the target parameter, φ j v represents the total porosity near the well at the j-th logging sampling point. cnj The content of connected porosity at the j-th well logging sampling point; and These are the observed P-wave impedance and S-wave impedance at the j-th logging sampling point, respectively. and V represents the predicted P-wave impedance and S-wave impedance at the j-th logging sampling point, respectively. P V S and ρ rk These represent the P-wave velocity, S-wave velocity, and density of the reservoir rock, respectively, K. rk and G rk These are the bulk modulus and shear modulus of the reservoir rock, respectively, N. s κ represents the total number of logging sampling points, and κ is the weighting coefficient of the error function.

[0136] To balance the magnitude difference between longitudinal and transverse wave impedances, the weighting coefficients of the error function are estimated using the following method.

[0137]

[0138] A fast simulated annealing algorithm is used to estimate the target parameters that minimize the error function, and to obtain the total porosity and connected porosity v near the well. j =[φ j v cnj ], j = 1, 2, 3, ..., N s In the example, the weight coefficient of the error function is k = 1.35.

[0139] 2. Based on the total porosity and connected pore content near the well, construct a sample set of target parameters, then estimate the prior probability density function and construct the posterior probability density function.

[0140] Based on the total porosity and connected pore content near the well, a Gaussian mixture model of the target parameters is constructed. The expectation-maximization algorithm is used to estimate the parameters of the Gaussian mixture model. The probability distribution function of the target parameters is then:

[0141]

[0142] In the formula, v = [φ, v cn[ ] represents the target parameters (including total porosity and connected pore content), N k Let μ represent the k-th Gaussian distribution. v|k ,∑ v|k and σ k Let N be the expectation, covariance, and weight coefficients of the k-th Gaussian distribution. c This represents the total number of Gaussian distributions.

[0143] Based on the probability distribution function of the target parameters, Monte Carlo stochastic simulation is used to obtain the augmented target parameters. These augmented target parameters are then combined with the total porosity and connected porosity content near the well to construct a Monte Carlo simulation sample set v for the target parameters. mci* =[φ i* v cni* ], i* = 1, 2, 3, ..., N mc , where N mc This represents the total number of samples in Monte Carlo.

[0144] Based on the target parameter sample set, the Epanechnikov kernel density estimation method is used to calculate the prior probability density function of the target parameters:

[0145]

[0146] in,

[0147]

[0148] In the formula, P(v) is the prior probability density function of the target parameter, v = [φ, v cn [ ] represents the target parameters (including total porosity and connected pore content), N mc Let L be the total number of Monte Carlo samples, L be the dimension of the target parameters, f(x) be the Epanechnikov kernel function, and h be the number of samples. v The scaling factor for the kernel function is N; in this example, the total number of Gaussian distributions, the total number of Monte Carlo samples, and the scaling factor for the kernel function are all set to N. c =3, N mc =10,000, h ν =5.

[0149] Based on the Bayesian framework, the posterior probability density function of the objective parameter is equivalent to:

[0150]

[0151] In the formula, v = [φ, v cn] represents the target parameters (including total porosity and connected porosity content), Z represents the observed impedance data (including observed P-wave impedance and S-wave impedance), obtained from pre-stack seismic inversion, P(v|Z) and P(Z|v) are the posterior probability density function and likelihood function of the target parameters, respectively, and P(Z) is the marginal probability of the observed impedance data, which is a constant.

[0152] Assuming the observed impedance data error follows a Gaussian distribution with zero expectation, the likelihood function of the target parameter can be expressed as:

[0153]

[0154] In the formula, D is the dimension of the observed impedance data. Combining equations (19), (21), and (22), the posterior probability density function of the target parameter can be expressed as:

[0155]

[0156] In the formula, Z pdt The predicted impedance data (including predicted P-wave impedance and predicted S-wave impedance) are obtained from the rock physics model constructed in step one, C nor Let ∑ be the normalization constant. z This is the covariance matrix of the observed impedance data error;

[0157] 3. Update the target parameter model, determine acceptance, and repeat the update step until the maximum number of iterations is reached. Calculate the expected mean of the connected porosity content model:

[0158] Based on the prior probability density function of the target parameters, an updated model of the target parameters is randomly generated. The model acceptance probability is calculated based on the Metropolis criterion. Therefore, the model acceptance probability at the k-th iteration can be expressed as:

[0159]

[0160] Among them, v * For the update model of the objective parameters, v (k) Let P be the target parameter before the update in the k-th iteration. acp Let P(v|Z) be the model reception probability, P(v|Z) be the posterior probability density function of the target parameters, and Z be the observed impedance data (including observed P-wave impedance and S-wave impedance).

[0161] To determine whether to accept the updated model, generate a random number between 0 and 1. If the random number is greater than or equal to P... acp Then the updated model that receives the target parameters (i.e., in the (k+1)th iteration, v) (k+1) =v * Otherwise, the objective parameter model remains unchanged (i.e., in the (k+1)th iteration, v...). (k+1) =v(k) ).

[0162] Repeat step five until the number of iterations k reaches the maximum number of iterations k. max (k max (greater than or equal to 20000), select the target parameter model from the last M iterations (e.g., 1000), i.e., v est =[v w v w+1 v w+2 , ..., v w+M-1 ], w = k max -M+1, based on v est The expected mean of the interconnected porosity content is estimated to be the final inversion result. In this embodiment, the maximum number of iterations is set to k. max =30000.

[0163] Figure 2 These are the P-wave impedance and S-wave impedance data obtained from well logging observations in this embodiment. Figure 3 In this embodiment, the total porosity and connected pore content near the well are obtained based on well logging data. Figure 4 (a) in this embodiment shows the total porosity inversion result using conventional methods (i.e., without considering connected and isolated pores). Figure 4 (b) in the figure shows the predicted results of longitudinal and transverse wave impedance data using conventional methods in this embodiment. From... Figure 4 As can be seen in (b), the predicted results (dashed line) of the longitudinal wave impedance and transverse wave impedance data deviate significantly from the observed data (solid line). Figure 5 In this embodiment, (a) represents the inversion results of total porosity and connected pore content using the method of the present invention. Figure 5 (b) in the figure represents the predicted longitudinal and transverse impedance data using the method of the present invention in this embodiment. Figure 5 As shown in (b), the predicted results (dashed lines) of the longitudinal wave impedance and transverse wave impedance data are basically consistent with the observed data (solid lines). The calculated results... Figure 5 In (b) of the study, the correlation coefficients between the predicted P-wave impedance and S-wave impedance data and the observed data are 0.976 and 0.968, respectively. Figure 4 The correlation coefficients in (b) are 0.933 and 0.914, respectively, which shows that the method of the present invention can effectively improve the accuracy of reservoir seismic prediction.

[0164] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

[0165] The above are merely preferred embodiments of the present invention. The scope of protection of the present invention is not limited to the above embodiments. All technical solutions falling within the scope of the present invention's concept are within the scope of protection of the present invention. It should be noted that for those skilled in the art, any improvements and modifications made without departing from the principles of the present invention should be considered within the scope of protection of the present invention.

Claims

1. A seismic rock physics nonlinear inversion method for determining the interconnected porosity content, characterized in that, The seismic rock physics nonlinear inversion method includes the following steps: Step A: Calculate the rock skeleton, pore fluid, and elastic parameters of fluid-containing rocks to construct an earthquake rock physics model; Step B: Based on well logging data and seismic rock physics model, obtain the total porosity and connected pore content near the well. Step C: Based on the total porosity and connected pore content near the well, construct a Monte Carlo simulation sample set for the target parameters; Step D: Based on the Monte Carlo simulation sample set of the target parameters, estimate the prior probability density function and construct the posterior probability density function of the target parameters; Step E involves iteratively updating the target parameters. In each iteration, an updated model for the target parameters is randomly generated based on the probability density function of the target parameters, and the model acceptance probability is calculated. The updated model is then accepted based on the model acceptance probability. Step F, repeat step E until the maximum number of iterations is reached, output the target parameter model, calculate the expected mean of the connected porosity content model, and use it as the final inversion result; Step A, which involves calculating the rock skeleton, pore fluid, and elastic parameters of fluid-containing rocks to construct an earthquake rock physics model, includes the following sub-steps: Step A1: Set basic rock physical parameters, including the bulk modulus, shear modulus, and density of minerals, and the bulk modulus and density of fluids; determine the percentage content of rock minerals based on core X-ray spectroscopy analysis, and calculate the bulk modulus, shear modulus, and density of the rock matrix using the VRH model. (1)； (2)； (3)； In the formula, K m G m and ρ m These represent the bulk modulus, shear modulus, and density of the rock matrix, respectively, K. mi G mi ρ mi and v mi These represent the bulk modulus, shear modulus, density, and percentage content of the i-th mineral, respectively, where N is the total number of mineral types. Step A2: Divide the rock pores into connected pores and isolated pores, and calculate the equivalent elastic modulus of the rock skeleton using a self-consistent model: (4)； (5)； In the formula, K * and G * ϕ represents the equivalent bulk modulus and shear modulus of the rock skeleton, respectively; ϕ represents the total porosity near the well; v cn and v ds These represent the content of connected pores and the content of isolated pores, respectively, satisfying v cn + v ds = 1; K cn and G cn K represents the bulk modulus and shear modulus of the interconnected pore inclusions, respectively. ds and G ds These are the bulk modulus and shear modulus of an isolated porous inclusion, respectively. and These are the geometric factors of the interconnected pore inclusions. and These are the geometric factors of isolated porous inclusions. and These are the geometric factors of the rock matrix; Step A3: If connected pores are defined as empty pores and isolated pores as fluid-containing pores, then: (6)； (7)； In the formula, K fl The bulk modulus of the fluid; Step A4: Based on the rock skeleton containing isolated pores, the Biot-Gassmann model is used to calculate the bulk modulus and shear modulus of the reservoir rock for connected pores. (8)； (9)； In the formula, K rk and G rk These are the bulk modulus and shear modulus of the reservoir rock, respectively, K. sk Bulk modulus of a rock skeleton containing isolated pores: (10)； Step A5, the P-wave velocity, S-wave velocity, and density of the reservoir rock are expressed as follows: (11)； (12)； (13)； In the formula, V P V S and ρ rk These represent the P-wave velocity, S-wave velocity, and density of the reservoir rock, respectively, ρ fl The density of the fluid.

2. The seismic rock physics nonlinear inversion method for connected porosity content according to claim 1, characterized in that, Step B, the process of obtaining the total porosity and connected porosity content near the well based on well logging data, includes the following sub-steps: Step B1: Based on the P-wave impedance and S-wave impedance data from well logging observations, the total porosity and connected porosity content near the well are obtained by inverting the predicted data from the matching model with the well logging observation data, and an error function is constructed: (14)； in, (15)； (16)； In the formula, ν j =[ϕ j , v cnj ] represents the target parameter, ϕ j v represents the total porosity near the well at the j-th logging sampling point. cnj The content of connected porosity at the j-th well logging sampling point; and These are the observed P-wave impedance and S-wave impedance at the j-th logging sampling point, respectively. and V represents the predicted P-wave impedance and S-wave impedance at the j-th logging sampling point, respectively. P V S and ρ rk These represent the P-wave velocity, S-wave velocity, and density of the reservoir rock, respectively, K. rk and G rk Let J represent the bulk modulus and shear modulus of the reservoir rock, respectively, j=1,2,…,N. s N s κ represents the total number of logging sampling points, and κ is the weighting coefficient of the error function; Step B2, to balance the magnitude difference between the longitudinal wave impedance and the transverse wave impedance, estimates the weighting coefficients of the error function in the following way: (17)； Step B3: Using the fast simulated annealing algorithm, estimate the target parameters that minimize the error function, and obtain the total porosity and connected porosity ν at the j-th logging sampling point. j =[ϕ j , v cnj ].

3. The seismic rock physics nonlinear inversion method for connected porosity content according to claim 1, characterized in that, Step C, the process of constructing the target parameter sample set based on the total porosity and connected pore content near the well, includes the following sub-steps: Step C1: Based on the total porosity and connected pore content near the well, a Gaussian mixture model of the target parameters is constructed. The expectation-maximization algorithm is used to estimate the parameters of the Gaussian mixture model. The probability distribution function of the target parameters is then: (18)； In the formula, ν=[ϕ, v cn [] represents the target parameter, ϕ represents the total porosity near the well, and v cn For the content of interconnected pores, N k Let μ represent the k-th Gaussian distribution. ν|k Σ ν|k and σ k Let N be the expectation, covariance, and weight coefficients of the k-th Gaussian distribution. c The total number of Gaussian distributions; Step C2: Based on the probability distribution function of the target parameters, Monte Carlo random simulation is used to obtain the augmented target parameters. The augmented target parameters are then combined with the total porosity and connected porosity content near the well to construct a Monte Carlo simulation sample set ν of the target parameters. mci* = [ϕ i* , v cni* ], i*=1, 2, 3, ..., N mc , where N mc This represents the total number of samples in Monte Carlo.

4. The seismic rock physics nonlinear inversion method for connected porosity content according to claim 1, characterized in that, In step D, the process of estimating the prior probability density function and constructing the posterior probability density function based on the target parameter sample set includes the following sub-steps: Step D1: Based on the target parameter sample set, calculate the prior probability density function of the target parameters using the Epanechnikov kernel density estimation method. (19)； in, (20)； In the formula, P(ν) is the prior probability density function of the target parameter, ν=[ϕ, v cn [] represents the target parameter, ϕ represents the total porosity near the well, and v cn For the content of interconnected pores, N mc Let L be the total number of Monte Carlo samples, L be the dimension of the target parameters, f(x) be the Epanechnikov kernel function, and h be the number of samples. ν The scaling factor for the kernel function; Step D2, based on the Bayesian framework, the posterior probability density function of the target parameters is equivalent to: (21)； In the formula, Z represents the observed impedance data, including the observed P-wave impedance and S-wave impedance, which are obtained from pre-stack seismic inversion; P(ν|Z) and P(Z|ν) are the posterior probability density function and likelihood function of the target parameter, respectively; P(Z) is the marginal probability of the observed impedance data, which is a constant. Step D3, assuming the observed impedance data error follows a Gaussian distribution with zero expectation, then the likelihood function of the target parameter is expressed as: (22)； In the formula, D is the dimension of the observed impedance data. Combining equations (19), (21), and (22), the posterior probability density function of the target parameter is expressed as: (23)； In the formula, Z pdt The predicted impedance data, including predicted P-wave impedance and predicted S-wave impedance, are obtained from the rock physics model constructed in step A; C nor Σ is the normalization constant. z This is the covariance matrix of the error in the observed impedance data.

5. The seismic rock physics nonlinear inversion method for connected porosity content according to claim 1, characterized in that, Step E, the process of iteratively updating the target parameters, includes the following sub-steps: Based on the prior probability density function of the target parameters, an updated model of the target parameters is randomly generated. The model acceptance probability is calculated based on the Metropolis criterion. The model acceptance probability at the k-th iteration is expressed as: (24)； Where, ν * For the update model of the objective parameters, ν (k) Let P be the target parameter before the k-th iteration update. acp For the model's acceptance probability, and Z represents the posterior probability density function of the target parameter update model and the posterior probability density function of the target parameter before the k-th iteration update, respectively; Z represents the observed impedance data, including the observed P-wave impedance and S-wave impedance. Generate a random number between 0 and 1. If the random number is greater than or equal to P... acp Then, the updated model that receives the target parameters is set to let the target parameters ν before the (k+1)th iteration update. (k+1) =ν * Otherwise, the objective parameter model remains unchanged, and the objective parameter ν before the (k+1)th iteration update is set to... (k+1) =ν (k) .

6. The seismic rock physics nonlinear inversion method for connected porosity content according to claim 1, characterized in that, In step F, the maximum number of iterations k max It is a positive integer greater than or equal to 20000.

7. The seismic rock physics nonlinear inversion method for connected porosity content according to claim 1, characterized in that, In step F, the process of outputting the target parameter model and calculating the expected mean of the connected porosity content model as the final inversion result includes the following sub-steps: Select the target parameter model set ν for the last M iterations. est =[ν w, ν w+1, ν w+2, ..., ν w+M-1 ],ν w Let w be the target parameters before the w-th iteration update, where w=k max -M-1, based on the objective parameter model set ν est Estimate the expected mean of the connected pore content model as the final inversion result.

Images