Inversion-assisted heterogeneous model carbonate rock physical modeling method

Abstract

The invention provides an inversion-assisted heterogeneous model carbonate rock physical modeling method. The inversion-assisted heterogeneous model carbonate rock physical modeling method comprises the following steps: forming a Biot-Rayleigh rock physical modeling process by combining a dual-pore heterogeneous Biot-Rayleigh equation, a differential equivalent medium model and the like; in consideration of complex pore types in carbonate rocks, the content of soft pores is introduced to describe the relative content of soft pores and hard pores, and the relative content of soft pores and hard pores is set as a to-be-inverted parameter to be estimated through nonlinear inversion; the inversion result of the soft pore content can be used for describing the change condition of the relative content of different pore types along with the depth so as to correct the established Biot-Rayleigh rock physical model, and in order to improve the reliability of soft pore content inversion, an initial model, an initial particle and an initial search interval are estimated by adopting clay content experience. Actual well data tests show that the method provided by the invention has higher prediction precision than a conventional Gassmann modeling method.

Description

Technical Field

[0001] The invention relates to the field of unconventional reservoir seismic exploration, and in particular to an inversion-assisted heterogeneous model carbonate rock physics modeling method. Background Art

[0002] The seismic response of complex oil and gas reservoirs such as carbonate rocks is affected by the comprehensive influence of various physical property information inside them, such as lithology, facies, porosity, and saturated fluid related parameters. Rock physics is used to describe the mechanical and elastic properties of rocks, and establishes a quantitative relationship between rock elasticity and physical properties. This quantitative relationship can be an empirical formula obtained by fitting rock physics experimental data, but it relies on phenomenology and cannot be closely linked to objective physical principles. Theoretical rock physics models overcome this problem and physically describe the seismic / elastic response of rock properties. For example, equivalent medium theory, including inclusion models (such as self-consistent and differential equivalent medium models, etc.) and contact models (such as Hashin-Strikman and Voigt and Reuss limits, etc.). Another type is the two-phase medium theory that describes the problem of seismic wave propagation in fluid-containing porous rocks, such as the Gassmann equation, Biot theory, Biot / squirt (BISQ) model, and Biot-Rayleigth model. In order to facilitate the establishment of forward models and solve the corresponding inverse problems, current seismic rock physics inversion mostly uses rock physics models with uniform medium assumptions. Among them, the Gassmann equation is often used in seismic rock physics inversion to simulate the elastic properties of saturated rocks at the seismic scale. For example, in the prior art, when studying the problem of carbonate rock physics modeling, the influence of different pores, including hard pores, fractures, etc., was considered; in the modeling of seismic rock physics inversion, a modeling method that considers different pore types was also used to improve the accuracy of modeling; in order to facilitate the derivation of linear rock physics forward operators, the differential equivalent model and the Gassmann equation were combined to establish a rock physics model of carbonate reservoirs. The biggest problem with these methods is that they all use the Gassmann equation and its extended model with a homogeneous assumption. Although the complexity of pores in carbonate rocks is considered, it is impossible to describe the heterogeneous characteristics of carbonate rocks, and it is impossible to reasonably describe the rock property-elasticity relationship of complex carbonate reservoirs.

[0003] In summary, rock physics model is a bridge to establish rock elasticity and physical properties, and is an important basis for directly mining reservoir physical property information using seismic data. However, the current rock physics modeling methods for complex reservoirs such as carbonate rocks mostly use rock physics theory based on the assumption of uniform medium, which cannot describe the strong heterogeneity in carbonate rocks due to complex pores and other reasons. Summary of the invention

[0004] In order to overcome the problem that the rock physics modeling methods of complex reservoirs such as carbonate rocks in the prior art mostly adopt rock physics theories based on the assumption of uniform media, which cannot describe the strong heterogeneity in carbonate rocks due to complex pores and other reasons, the present invention provides an inversion-assisted heterogeneous model carbonate rock rock physics modeling method, which better describes the complex pore conditions in rocks, takes the soft pore content as a parameter to measure the soft and hard pores and as a parameter to be inverted, adopts a particle swarm-very fast simulated annealing algorithm to obtain the soft pore content inversion result, and corrects the established rock physics model to improve the model prediction accuracy.

[0005] The present invention provides the following technical solution: an inversion-assisted heterogeneous model carbonate rock physics modeling method, comprising the following steps:

[0006] S1. According to the volume fraction and modulus information of the mineral components in the well, the elastic modulus of the rock matrix is ​​obtained, and the initial values ​​of the soft pore content, inclusion content, soft pore porosity, and hard pore porosity are preset respectively;

[0007] S2, based on the rock matrix modulus, soft pores are added to the spherical inclusions and hard pores are added to the background main phase medium;

[0008] S3, adding spherical inclusions containing soft pores into the background main phase medium containing hard pores to form a dual porosity model rock dry skeleton;

[0009] S4. Based on the moduli of different fluids and the fluid saturation, the modulus of the mixed fluid is calculated, and the mixed fluid is added to the dry skeleton of the dual-porosity rock using the Biot-Rayleigh equation to calculate the modulus of the saturated fluid dual-porosity model rock;

[0010] S5, using the rock modulus of the dual porosity model to calculate the predicted P-wave velocity and S-wave velocity, and establish an inversion objective function;

[0011] S6, setting a search range for soft pore content, and generating initial particles of soft pore content based on an initial value of the soft pore content;

[0012] S7, using the particle swarm-fast simulated annealing optimization algorithm, the nonlinear inversion result of the soft pore content is calculated;

[0013] S8. Using the soft pore content inversion result as correction information, re-perform the rock physics modeling process, and return to execute steps S1 to S4 to obtain the final fluid-saturated dual porosity model rock modulus.

[0014] Furthermore, the aforementioned step S1 includes the following sub-steps:

[0015] S1.1. Based on the volume fraction and modulus information of the mineral components in the well, the elastic modulus of the rock matrix is ​​calculated using the Voigt-Reuss-Hill average, as follows:

[0016]

[0017] Among them, M mat is the elastic modulus of the rock matrix, M V represents the rock matrix elastic modulus obtained using the Voigt modulus, M R represents the rock matrix elastic modulus obtained using the Reuss model;

[0018] S1.2. Set the initial values ​​of soft pore content, inclusion content, soft pore porosity, and hard pore porosity as follows:

[0019]

[0020] in, is the initial value of soft pore content, is the initial value of inclusion content, is the initial value of soft pore porosity, is the initial value of hard pore porosity, α is the adjustment factor of the initial value of soft pore content, v sha is the clay content measured in the well.

[0021] Furthermore, in the aforementioned step S1.1, the rock matrix elastic modulus M R and M V As follows:

[0022]

[0023] M V =v qua M qua +v fel M fel +v sha M sha

[0024] Among them, v qua and M qua Represent the volume fraction and modulus of quartz, v fel and M fel represent the volume fraction and modulus of feldspar, v sha and M sha represent the clay content and modulus respectively.

[0025] Furthermore, the aforementioned step S2 includes the following sub-steps:

[0026] S2.1. Based on the rock matrix modulus, soft pores are added to the spherical inclusions as follows:

[0027]

[0028]

[0029]

[0030] in, and is the initial bulk modulus and initial shear modulus of the spherical inclusion; K sp and μ sp are the bulk modulus and shear modulus of the gradually added soft pores; is the initial value of soft pore porosity, K inc* and μ inc* is the equivalent bulk modulus and equivalent shear modulus of the spherical inclusion containing soft pores; P (*sp) and Q (*sp) It refers to the geometric factor of soft pores in the equivalent medium;

[0031] S2.1. Based on the rock matrix modulus, hard pores are added to the main phase of the background medium, as shown in the following expression:

[0032]

[0033]

[0034]

[0035] in, and is the initial bulk modulus and initial shear modulus of the background main phase medium; K hp and μ hp are the bulk modulus and shear modulus of the gradually added hard voids; is the initial value of hard pore porosity, K back* and μ back* is the equivalent bulk modulus and equivalent shear modulus of the background main phase medium containing hard pores; P (*h p ) and Q (*hp) Refers to the geometric factor of hard pores in the equivalent medium.

[0036] Furthermore, in the aforementioned step S3, the modulus of the dual-porosity rock skeleton is specifically expressed as follows:

[0037]

[0038]

[0039]

[0040] in, and is the initial bulk modulus and initial shear modulus of the main phase medium with hard pore background; K inc and μ inc are the bulk modulus and shear modulus of the soft pore inclusions added gradually; is the initial value of soft pore content, which is equal to the inclusion content; K dry* and μ dry* is the equivalent bulk modulus and equivalent shear modulus of the dry skeleton of the dual-porosity rock formed after the inclusions are added to the background main phase medium; P (*inc) and Q (*inc) It refers to the geometric factor of spherical inclusions containing soft pores in an equivalent medium.

[0041] Furthermore, the aforementioned step S4 includes the following sub-steps:

[0042] S4.1. Estimate the modulus of the mixed fluid based on the modulus and fluid saturation of different fluids. Assume that there is a mixture of two fluids in the rock pores. The bulk modulus and density expressions of the mixed fluid are as follows:

[0043]

[0044] ρ f =s1ρ1+(1-s1)ρ2

[0045] Among them, K f is the bulk modulus of the mixed fluid, ρ f is the density of the mixed fluid, K1 and K2 are the bulk moduli of the first and second fluids, ρ1 and ρ2 are the densities of the first and second fluids, and s1 is the saturation of the first fluid;

[0046] S4.2. The wave propagation equation of the dual porosity model is described using the Biot-Rayleigh equation. The Christoffel matrix of the Biot-Rayleigh equation is as follows:

[0047]

[0048] a 11 =A+2N+i(Q2φ1-Q1φ2)x1,a 12 =Q1+i(Q2φ1-Q1φ2)x2

[0049] a 13 =Q2+i(Q2φ1-Q1φ2)x3,a 21 =Q2-iR1φ2x3

[0050] a 22= R1 - iR1φ2x2, a 23 = -iR1φ2x3

[0051] a 31 = Q2 + iR2φ1x1, a 32 = iR2φ1x2, a 33 = R2 + iR2φ1x3

[0052] b 11 = iω(b1 + b2) - ρ 11 ω 2 ,b 12 = -iωb1 - ρ 11 ω 2

[0053] b 13 = -iωb2 - ρ 13 ω 2 ,b 21 = -iωb1 - ρ 12 ω 2

[0054] b 22 = iωb1 - ρ 22 ω 2 ,b 23 = 0, b 32 = 0

[0055] b 31 = -iωb2 - ρ 13 ω 2 ,b 33 = -ρ 33 ω 2 + iωb2

[0056] where k is the wave number, a 11 、a 12 、a 13 、a 21 、a 22 、a 23 、a 31 、a 32 、a 33 、b 11 、b 12 、b 13 、b 21 、b 22 、b 23 、b 31 、b 32 、b 33 are the coefficients of the wave equation, ρ 11 、ρ 12 、ρ 13 are the density parameters of the main phase medium, ρ 22and ρ 33 is the inclusion density parameter, φ1 is the main phase medium porosity, φ2 is the inclusion porosity, ω is the angular frequency, i is the imaginary number sign, x1, x2, x3 are attenuation factor parameters, A, N, Q1, Q2, R1, R2 are Biot stiffness coefficients

[0057] S4.3. Solve the above Christoffel matrix to obtain the P-wave and S-wave velocities of the dual-porosity rock saturated with fluid, as follows:

[0058]

[0059] in, and are the predicted P-wave velocity and S-wave velocity of fluid-saturated dual-porosity rock, μ dry is the shear modulus of the doubly porous rock, ρ pred represents the predicted density of fluid-saturated dual-porosity rock;

[0060] S4.4. Based on the fact that the P-wave velocity and S-wave velocity of fluid-saturated dual-porosity rock are functions of the soft pore content, a forward operator from rock physical parameters and soft pore content to saturated rock velocity information is established as follows:

[0061]

[0062] Where G represents the forward operator formed by Biot-Rayleigh dual-porosity rock physics modeling.

[0063] Furthermore, the aforementioned step S5 establishes the inversion objective function as follows:

[0064]

[0065] Among them, v sp represents the soft pore content to be inverted, J represents the inversion objective function, and denote the P-wave velocity and S-wave velocity predicted by dual-porosity rock physics modeling, and denote the measured P-wave velocity and S-wave velocity in the well, σ1 and σ2 denote the weight parameters of the P-wave velocity data item and S-wave velocity data item, λ is the regularization parameter, u vSP The expected soft pore content.

[0066] Furthermore, the aforementioned step S6 specifically includes: in the first inversion iteration, the search range of the soft pore content is set based on the initial value of the clay content, and the search range is generated by the following expression:

[0067]

[0068] After the second inversion iteration, the search range of the current iteration is set based on the soft hole content update result obtained in the previous iteration, which has the following expression:

[0069]

[0070] Where k represents the kth inversion iteration, N is the preset maximum number of iterations, and They represent the minimum and maximum values ​​of the soft pore content search range in the kth iteration, a is the adjustment factor of the initial value of the soft pore content, b is the adjustment factor of the soft pore content search range, and v sha is the clay content measured in the well, For the clay content expectation, is the estimated result of soft pore content in the k-1th iteration, is the expectation of the soft pore content estimation result in the k-1th iteration;

[0071] Based on the clay content in the well, generate initial particles with soft pore content; set M initial particles, particle group It is expressed as follows:

[0072]

[0073] in, The initial particle vector representing the soft hole content, represents the ith particle, and the soft pore content of the ith particle to be inverted is obtained by the following formula:

[0074]

[0075] Where rand(·) is a random number between [0,1], and c is the fine-tuning scale.

[0076] Furthermore, the aforementioned step S7 includes the following sub-steps:

[0077] S7.1. Substituting all particles into the objective function can obtain the fitness of different particles, and selecting the particle with the smallest fitness as the individual optimal particle:

[0078]

[0079]

[0080] Where k is the number of iterations, i is the i-th particle, is the soft pore content of the i-th initial particle, represents the individual optimal position of the i-th particle in the first iteration, represents the individual optimal position of the i-th particle in the k-th iteration, v sp,i,k represents the soft hole content of the i-th particle in the k-th iteration. J(v sp,i,k ) and J(v sp,i,k-1 ) represent the fitness of the ith particle in the kth and k-1th iterations respectively, which is calculated by substituting the particle into the inversion objective function; S7.2, the reception probability is calculated according to the fitness of different particles. After the reception probability of all particles is obtained, the roulette method is used to determine the optimal position of the group and the optimal particle of the group; the reception probability of the particle is calculated as follows:

[0081]

[0082] Among them, P accept (v sp,i,k ) represents the receiving probability of the i-th particle in the k-th iteration, represents the optimal position of the group in the kth iteration, represents the fitness of the optimal position of the group in the kth iteration, T k is the temperature value of the kth iteration, τ is a constant, and Σ is the summation symbol;

[0083] S7.3. Calculate the update rate based on the individual optimal position and the group optimal position. The expression is as follows:

[0084]

[0085] Among them, v sp,i,k represents the i-th particle in the k-th iteration, represents the individual optimal position of the i-th particle in the k-th iteration, represents the optimal position of the group at the kth iteration. p,i,k+1 and r sp,i,k Respectively represent the update rate of the i-th particle in the k+1th and kth iterations. rand1(·) and rand2(·) represent two random numbers between [0,1]. κ1 and κ2 are the update step sizes S7.4. Based on the update rate, update different particles in the particle swarm as follows:

[0086] v sp,i,k+1 =v sp,i,k +κ3·r sp,i,k+1

[0087] Among them, κ3 is the update step size, v sp,i,k+1 represents the i-th particle in the k+1-th iteration;

[0088] S7.5, judging whether the number of iterations reaches the preset maximum number of iterations, if yes, taking the updated soft hole content as the inversion result, otherwise executing step S7.6;

[0089] S7.6, substitute the updated soft pore content into step S1, return to execute steps S1 to S7, continue the inversion iteration until the preset maximum number of iterations is reached, and take the soft pore content output at the maximum number of iterations as the inversion result.

[0090] Furthermore, the aforementioned step S8 specifically includes: taking the inversion result of the soft pore content as known data, inputting it into step S1 to replace the initial value of the soft pore content, repeating steps S1 to S4, and predicting the longitudinal wave velocity, shear wave velocity and density, as shown in the following formula:

[0091]

[0092] in, is the inversion result of soft pore content, G represents the forward operator formed by the Biot-Rayleigh dual pore rock physics modeling from step S1 to step S4, To estimate the P-wave velocity using the soft pore content inversion results, To estimate the shear wave velocity using the inversion results of soft pore content, The density is estimated using the inversion results of soft pore content.

[0093] Compared with the prior art, the beneficial technical effects of the above technical solutions adopted by the present invention are as follows: the previous reservoir physical property parameter inversion adopts the rock physics theory of uniform medium assumption, such as Gassmann model, which cannot describe the complex pore structure in complex reservoirs or the heterogeneity caused by patchy fluids, which will lead to limited accuracy of reservoir prediction. The present invention introduces the Biot-Rayleigh model of heterogeneous media into carbonate rock physics modeling to improve the simulation accuracy of rock physics modeling. In addition, in order to consider the complex pore types in carbonate rocks, the soft pore content is introduced to describe the relative content of soft pores and hard pores, and it is set as the parameter to be inverted and estimated by nonlinear inversion; the inversion result of the soft pore content can describe the change of the relative content of different pore types with depth, so as to correct the established Biot-Rayleigh rock physics model and improve the prediction accuracy. Compared with conventional homogeneous medium rock physics modeling, the present invention has higher prediction accuracy and is more suitable for describing the rock elastic properties of complex heterogeneous reservoirs. BRIEF DESCRIPTION OF THE DRAWINGS

[0094] Figure 1 The figure is a flow chart of the method of the present invention.

[0095] Figure 2 Flow chart of rock physics modeling of heterogeneous carbonate rocks based on the Biot-Rayleigh model.

[0096] Figure 3Schematic diagram of the measured well data used in the test, in which (a) is a schematic diagram of the longitudinal wave velocity, (b) is a schematic diagram of the shear wave velocity, (c) is a schematic diagram of the density, (d) is a schematic diagram of the porosity, (e) is a schematic diagram of the water saturation, and (f) is a schematic diagram of the clay content.

[0097] Figure 4 These are the results of conventional Gassmann rock physics modeling predictions. In the figure, (a) is a schematic diagram of the P-wave velocity prediction results, (b) is a schematic diagram of the S-wave velocity prediction results, and (c) is a schematic diagram of the density prediction results.

[0098] Figure 5 Figure 2 shows the prediction results of Biot-Rayleigh rock physics modeling when the soft pore content is fixed. (a) is a schematic diagram of the P-wave velocity prediction results, (b) is a schematic diagram of the S-wave velocity prediction results, and (c) is a schematic diagram of the density prediction results.

[0099] Figure 6 Figure 1 shows the prediction results of Biot-Rayleigh rock physics modeling when inverting the soft pore content. In the figure, (a) is a schematic diagram of the P-wave velocity prediction results, (b) is a schematic diagram of the S-wave velocity prediction results, (c) is a schematic diagram of the density prediction results, and (d) is a schematic diagram of the soft pore content prediction results. DETAILED DESCRIPTION

[0100] In order to better understand the technical content of the present invention, specific embodiments are given and described as follows in conjunction with the accompanying drawings.

[0101] Various aspects of the invention are described herein with reference to the accompanying drawings, in which many illustrative embodiments are shown. The embodiments of the invention are not limited to those described in the accompanying drawings. It should be understood that the invention is implemented by any of the various concepts and embodiments described above, as well as the concepts and embodiments described in detail below, because the concepts and embodiments disclosed in the invention are not limited to any implementation. In addition, some aspects disclosed in the invention may be used alone or in any appropriate combination with other aspects disclosed in the invention.

[0102] refer to Figure 1 The present invention provides an inversion-assisted heterogeneous model carbonate rock physics modeling method, comprising the following steps:

[0103] S1. According to the volume fraction and modulus information of the mineral components in the well, the elastic modulus of the rock matrix is ​​obtained, and the initial values ​​of the soft pore content, inclusion content, soft pore porosity, and hard pore porosity are preset respectively;

[0104] S2, based on the rock matrix modulus, soft pores are added to the spherical inclusions and hard pores are added to the background main phase medium;

[0105] S3, adding spherical inclusions containing soft pores into the background main phase medium containing hard pores to form a dual porosity model rock dry skeleton;

[0106] S4. Based on the moduli of different fluids and the fluid saturation, the modulus of the mixed fluid is calculated, and the mixed fluid is added to the dry skeleton of the dual-porosity rock using the Biot-Rayleigh equation to calculate the modulus of the dual-porosity model rock saturated with fluid;

[0107] S5, using the rock modulus of the dual porosity model to calculate the predicted P-wave velocity and S-wave velocity, and establish an inversion objective function;

[0108] S6, setting a search range for soft pore content, and generating initial particles of soft pore content based on an initial value of the soft pore content;

[0109] S7, using the particle swarm-fast simulated annealing optimization algorithm, the nonlinear inversion result of the soft pore content is calculated;

[0110] S8. Using the soft pore content inversion result as correction information, re-perform the rock physics modeling process, and return to execute steps S1 to S4 to obtain the final fluid-saturated dual porosity model rock modulus.

[0111] Figure 2 The heterogeneous carbonate rock physical modeling process based on the Biot-Rayleigh model adopted in the present invention is shown, corresponding to steps S1 to S4.

[0112] As a preferred embodiment of the present invention, step S1 first uses Voigt-Reuss-Hill (VRH) to averagely calculate the elastic modulus of the mineral mixture. Based on the volume fraction and modulus information of the mineral components in the well, the elastic modulus of the rock matrix is ​​averaged using Voigt-Reuss-Hill (VRH), as follows:

[0113]

[0114] Among them, M mat is the elastic modulus of the rock matrix, M V represents the rock matrix elastic modulus obtained using the Voigt modulus, M R It represents the rock matrix elastic modulus obtained using the Reuss model. The specific expression is as follows:

[0115]

[0116] M V =v qua M qua +v fel M fel +v sha Msha

[0117] Among them, v qua and M qua Represent the volume fraction and modulus of quartz, v fel and M fel represent the volume fraction and modulus of feldspar, v sha and M sha represent the clay content and modulus respectively.

[0118] Secondly, the initial values ​​of soft pore content and soft pore porosity in the well are formed. Since there is no measured soft pore content in the well, a reasonable initial value of the soft pore content needs to be set as an important auxiliary parameter. Since soft pores are mostly developed in clay, it can be considered that there is a certain correlation between the soft pore content and the clay content. The initial values ​​of soft pore content, inclusion content, soft pore porosity, and hard pore porosity are set as follows:

[0119]

[0120] in, is the initial value of soft pore content, is the initial value of inclusion content, is the initial value of soft pore porosity, is the initial value of hard pore porosity, and α is the adjustment factor of the initial value of soft pore content.

[0121] As a preferred embodiment of the present invention, step S2 specifically uses an isotropic differential effective medium (DEM) model to add soft pores and hard pores to the inclusions and the main phase medium, and estimates the elastic modulus of the inclusions containing soft pores and the main phase medium containing hard pores, respectively. It includes the following sub-steps:

[0122] S2.1. Based on the rock matrix modulus, soft pores are added to the spherical inclusions as follows:

[0123]

[0124]

[0125]

[0126] in, and is the initial bulk modulus and initial shear modulus of the spherical inclusion; K sp and μ sp are the bulk modulus and shear modulus of the gradually added soft pores; is the initial value of soft pore porosity, K inc* and μ inc* is the equivalent bulk modulus and equivalent shear modulus of the spherical inclusion containing soft pores; P(*sp) and Q (*sp) Refers to the geometric factor of soft pores in the equivalent medium.

[0127] S2.1. Based on the rock matrix modulus, hard pores are added to the main phase of the background medium, as shown in the following expression:

[0128]

[0129]

[0130]

[0131] in, and is the initial bulk modulus and initial shear modulus of the background main phase medium; K hp and μ hp are the bulk modulus and shear modulus of the gradually added hard voids; is the initial value of hard pore porosity, K back* and μ back* is the equivalent bulk modulus and equivalent shear modulus of the background main phase medium containing hard pores; P (*hp) and Q (*hp) Refers to the geometric factor of hard pores in the equivalent medium.

[0132] As a preferred embodiment of the present invention, in step S3, inclusions are added to the main phase medium using an isotropic DEM model to obtain a dual-porosity carbonate rock dry skeleton. Based on the modulus of the spherical inclusions containing soft pores and the background main phase medium containing hard pores, the inclusions are added to the background medium main phase to form the modulus of the dual-porosity rock dry skeleton, as shown in the following expression:

[0133]

[0134]

[0135]

[0136] in, and is the initial bulk modulus and initial shear modulus of the main phase medium with hard pore background; K inc and μ inc are the bulk modulus and shear modulus of the soft pore inclusions added gradually; is the initial value of soft pore content, which is equal to the inclusion content; K dry* and μ dry* is the equivalent bulk modulus and equivalent shear modulus of the dry skeleton of the dual-porosity rock formed after the inclusions are added to the background main phase medium; P (*inc) and Q (*inc)It refers to the geometric factor of spherical inclusions containing soft pores in an equivalent medium.

[0137] As a preferred embodiment of the present invention, step S4 includes the following sub-steps:

[0138] S4.1. Estimate the modulus of the mixed fluid based on the modulus and fluid saturation of different fluids. Assume that there is a mixture of two fluids in the rock pores. The bulk modulus and density of the mixed fluid are expressed as follows:

[0139]

[0140] ρ f =s1ρ1+(1-s1)ρ2

[0141] Among them, K f is the bulk modulus of the mixed fluid, ρ f is the density of the mixed fluid, K1 and K2 are the bulk moduli of the first and second fluids respectively, ρ1 and ρ2 are the densities of the first and second fluids respectively, and s1 is the saturation of the first fluid.

[0142] S4.2. Use the Biot-Rayleigh equation to add the mixed fluid into the dry skeleton of the dual-porosity rock to estimate the rock modulus of the dual-porosity model saturated with fluid. The Biot-Rayleigh equation is a wave propagation equation that describes the dual-porosity model, and its Christoffel equation is:

[0143]

[0144] a 11 =A+2N+i(Q2φ1-Q1φ2)x1,a 12 =Q1+i(Q2φ1-Q1φ2)x2

[0145] a 13 =Q2+i(Q2φ1-Q1φ2)x3,a 21 =Q2-iR1φ2x3

[0146] a 22 =R1-iR1φ2x2, a 23 =-iR1φ2x3

[0147] a 31 =Q2+iR2φ1x1,a 32 =iR2φ1x2,a 33 =R2+iR2φ1x3

[0148] b 11 =iω(b1+b2)-ρ 11 ω 2 , b12 =-iωb1-ρ 11 ω 2

[0149] b 13 =-iωb2-ρ 13 ω 2 , b 21 =-iωb1-ρ 12 ω 2

[0150] b 22 =iωb1-ρ 22 ω 2 , b 23 =0, b 32 =0

[0151] b 31 =-iωb2-ρ 13 ω 2 , b 33 =-ρ 33 ω 2 +iωb2

[0152] Where k is the wave number, a 11 、a 12 、a 13 、a 21 、a 22 、a 23 、a 31 、a 32 、a 33 、b 11 、b 12 、b 13 、b 21 、b 22 、b 23 、b 31 、b 32 、b 33 is the wave equation coefficient, ρ 11 , 12 , 13 is the density parameter of the main phase medium, ρ 22 and ρ 33 is the inclusion density parameter, φ1 is the main phase medium porosity, φ2 is the inclusion porosity, ω is the angular frequency, i is the imaginary number sign, x1, x2, x3 are attenuation factor parameters, A, N, Q1, Q2, R1, R2 are Biot stiffness coefficients

[0153] S4.3. Solve the above Christoffel matrix to obtain the P-wave and S-wave velocities of the dual-porosity rock saturated with fluid, as follows:

[0154]

[0155] in, and represent the P-wave velocity and S-wave velocity of fluid-saturated dual-porosity rock, μ dry is the shear modulus of the doubly porous rock, ρ pred represents the predicted density of a fluid-saturated dual-porosity rock.

[0156] S4.4. Based on the function of soft pore content based on the P-wave velocity and S-wave velocity of dual-porosity rock saturated with fluid, a forward operator from rock physical parameters and soft pore content to saturated rock velocity information is established as follows:

[0157]

[0158] where G represents the forward operator formed by the Biot-Rayleigh dual-porosity rock physics modeling described in steps S1 to S4.

[0159] As a preferred embodiment of the present invention, in step S5, the predicted P- and S-wave velocities of fluid-saturated dual-porosity rocks are obtained by the Biot-Rayleigh equation, the soft pore content is set as the parameter to be inverted, and the inversion objective function is established as follows:

[0160]

[0161] Among them, v sp represents the soft pore content to be inverted, J represents the inversion objective function, and denote the P-wave velocity and S-wave velocity predicted by dual-porosity rock physics modeling, and denote the measured P-wave velocity and S-wave velocity in the well, σ1 and σ2 denote the weight parameters of the P-wave velocity data item and S-wave velocity data item, λ is the regularization parameter, u vSP The expected soft pore content.

[0162] As a preferred embodiment of the present invention, step S6 is specifically: in the first inversion iteration, the search range of the soft pore content is set based on the initial value of the clay content, and the search range is generated by the following expression:

[0163]

[0164] After the second inversion iteration, the search range of the current iteration is set based on the soft hole content update result obtained in the previous iteration, which has the following expression:

[0165]

[0166] Where n represents the nth inversion iteration, N is the preset maximum number of iterations, and They represent the minimum and maximum values ​​of the soft pore content search range in the nth iteration, a is the adjustment factor of the initial value of the soft pore content, b is the adjustment factor of the soft pore content search range, and v sha is the clay content, For the clay content expectation, is the estimated result of soft pore content in the n-1th iteration, is the expectation of the soft pore content estimation result in the n-1th iteration;

[0167] Based on the clay content in the well, generate initial particles with soft pore content; set M initial particles, particle group It is expressed as follows:

[0168]

[0169] in, The initial particle vector representing the soft hole content, The i-th particle representing the soft pore content, the i-th particle of the soft pore content parameter to be inverted is obtained by the following formula:

[0170]

[0171] Where rand(·) is a random number between [0,1], and c is the fine-tuning scale.

[0172] As a preferred embodiment of the present invention, step S7 includes the following sub-steps:

[0173] S7.1. Substituting all particles into the objective function can obtain the fitness of different particles, and selecting the particle with the smallest fitness as the individual optimal particle:

[0174]

[0175]

[0176] Where k is the number of iterations, i is the i-th particle, is the soft pore content of the ith initial particle, represents the individual optimal position of the i-th particle in the first iteration, represents the individual optimal position of the i-th particle in the k-th iteration, v sp,i,k represents the i-th particle in the k-th iteration. J(v sp,i,k ) and J(v sp,i,k-1 ) represent the fitness of the ith particle in the kth and k-1th iterations, respectively, and are calculated by substituting the particle into the inversion objective function;

[0177] S7.2. Calculate the receiving probability according to the fitness of different particles. After obtaining the receiving probability of all particles, use the roulette method to determine the optimal position of the group and the optimal particle of the group. The receiving probability of the particle is calculated as follows:

[0178]

[0179] Among them, P accept (v sp,i,k ) represents the receiving probability of the i-th particle in the k-th iteration, represents the optimal position of the group in the kth iteration, represents the fitness of the optimal position of the group in the kth iteration, T k is the temperature value of the kth iteration, τ is a constant, and Σ is the summation symbol;

[0180] S7.3. Calculate the update rate based on the individual optimal position and the group optimal position. The expression is as follows:

[0181]

[0182] Among them, v sp,i,k represents the i-th particle in the k-th iteration, represents the individual optimal position of the i-th particle in the k-th iteration, represents the optimal position of the group at the kth iteration. p,i,k+1 and r sp,i,k Respectively represent the update rate of the i-th particle in the k+1th and kth iterations. rand1(·) and rand2(·) represent two random numbers between [0,1]. κ1 and κ2 are the update step sizes S7.4. Based on the update rate, update different particles in the particle swarm as follows:

[0183] v sp,i,k+1 =v sp,i,k +κ3·r i,k+1

[0184] Among them, κ3 is the update step size, v sp,i,k+1 represents the i-th particle in the k+1-th iteration;

[0185] S7.5, judging whether the number of iterations reaches the preset maximum number of iterations, if yes, taking the updated soft hole content as the inversion result, otherwise executing step S7.6;

[0186] S7.6, substitute the updated soft pore content into step S1, return to execute steps S1 to S7, continue the inversion iteration until the preset maximum number of iterations is reached, and take the soft pore content output at the maximum number of iterations as the inversion result.

[0187] Further, as a preferred embodiment of the present invention, step S8 specifically comprises: taking the inversion result of the soft pore content as known data, inputting it into step S1 to replace the initial value of the soft pore content, repeating steps S1 to S4, and predicting the longitudinal wave velocity and the shear wave velocity, as shown in the following formula:

[0188]

[0189] in, is the inversion result of soft pore content, G represents the forward operator formed by the Biot-Rayleigh dual pore rock physics modeling from step S1 to step S4, To estimate the P-wave velocity using the soft pore content inversion results, To estimate the shear wave velocity using the inversion results of soft pore content, The density is estimated using the inversion results of soft pore content.

[0190] like Figure 3 As shown, Figure 3 The well data curves used for testing are shown, in which (a) shows the P-wave velocity, (b) shows the S-wave velocity, (c) shows the density, (d) shows the porosity, (e) shows the water saturation, and (f) shows the clay content. Figure 4 The prediction results of the conventional rock physics modeling method based on the Gassmann equation are shown. In the figure, (a) shows the prediction results of the P-wave velocity, (b) shows the prediction results of the S-wave velocity, and (c) shows the density prediction results. The black curve is the measured data, and the blue dotted line is the prediction result. It can be seen that the prediction accuracy of the conventional modeling method is low, and the prediction results of the P-wave velocity and S-wave velocity are not ideal. Figure 5 The prediction results of the heterogeneous carbonate rock physics modeling method based on the Biot-Rayleigh model are shown. In the figure, (a) shows the prediction results of the P-wave velocity, (b) shows the prediction results of the S-wave velocity, and (c) shows the density prediction results. The soft pore content is estimated empirically, that is, it is estimated using the clay content; it can be seen that compared with the conventional Gassmann method, the prediction effect of the Biot-Rayleigh rock physics modeling of the empirical soft pore content is significantly improved, and the correlation between the predicted P-wave velocity and S-wave velocity and the real data is significantly improved. Figure 6The prediction results of the inversion-assisted Biot-Rayleigh model rock physics modeling method are shown. In the figure, (a) shows the prediction results of P-wave velocity, (b) shows the prediction results of S-wave velocity, (c) shows the prediction results of density, and (d) shows the prediction results of soft pore content. The soft pore content is obtained by inversion; the black curve is the measured data, the blue dotted line is the prediction result, and the red curve is the inverted soft pore content; it can be seen that among the three methods, the inversion-assisted Biot-Rayleigh method has the highest prediction accuracy, and the predicted P-wave velocity and S-wave velocity have the highest correlation with the real data, which shows the effectiveness of the method.

[0191] Although the present invention has been described above with preferred embodiments, it is not intended to limit the present invention. A person skilled in the art of the present invention may make various modifications and improvements without departing from the spirit and scope of the present invention. Therefore, the scope of protection of the present invention shall be determined by the definition of the claims.

Claims

1. An inversion-assisted heterogeneous model carbonate rock physics modeling method, characterized in that: The following steps are involved: S1. According to the volume fraction and modulus information of the mineral components in the well, the elastic modulus of the rock matrix is ​​obtained, and the initial values ​​of the soft pore content, inclusion content, soft pore porosity, and hard pore porosity are preset respectively; S2, based on the rock matrix modulus, soft pores are added to the spherical inclusions and hard pores are added to the background main phase medium; S3, adding spherical inclusions containing soft pores into the background main phase medium containing hard pores to form a dual porosity model rock dry skeleton; S4. Based on the moduli of different fluids and the fluid saturation, the modulus of the mixed fluid is calculated, and the mixed fluid is added to the dry skeleton of the dual-porosity rock using the Biot-Rayleigh equation to calculate the modulus of the saturated fluid dual-porosity model rock; S5, using the rock modulus of the dual porosity model to calculate the predicted P-wave velocity and S-wave velocity, and establish an inversion objective function; S6, setting a search range for soft pore content, and generating initial particles of soft pore content based on an initial value of the soft pore content; S7, using the particle swarm-fast simulated annealing optimization algorithm, the nonlinear inversion result of the soft pore content is calculated; S8. Using the soft pore content inversion result as correction information, re-perform the rock physics modeling process, and return to execute steps S1 to S4 to obtain the final fluid-saturated dual porosity model rock modulus.

2. The inversion-assisted heterogeneous model carbonate rock physics modeling method according to claim 1, characterized in that: Step S1 includes the following sub-steps: S1.

1. Based on the volume fraction and modulus information of the mineral components in the well, the elastic modulus of the rock matrix is ​​calculated using the Voigt-Reuss-Hill average, as follows: Among them, M mat is the elastic modulus of the rock matrix, M V represents the rock matrix elastic modulus obtained using the Voigt modulus, M R represents the rock matrix elastic modulus obtained using the Reuss model; S1.

2. Set the initial values ​​of soft pore content, inclusion content, soft pore porosity, and hard pore porosity as follows: in, is the initial value of soft pore content, is the initial value of inclusion content, is the initial value of soft pore porosity, is the initial value of hard pore porosity, α is the adjustment factor of the initial value of soft pore content, v sha is the clay content measured in the well.

3. The inversion-assisted heterogeneous model carbonate rock physics modeling method according to claim 2 is characterized in that: In step S1.1, the rock matrix elastic modulus M R and M V As follows: M V =in qua M qua +v fel M fel +v sha M sha Among them, v qua and M qua Represent the volume fraction and modulus of quartz, v fel and M fel denote the volume fraction and modulus of feldspar, v sha and M sha represent the clay content and modulus respectively.

4. The inversion-assisted heterogeneous model carbonate rock physics modeling method according to claim 2, characterized in that: Step S2 includes the following sub-steps: S2.

1. Based on the rock matrix modulus, soft pores are added to the spherical inclusions as follows: in, and is the initial bulk modulus and initial shear modulus of the spherical inclusion; K sp and μ sp are the bulk modulus and shear modulus of the gradually added soft pores; is the initial value of soft pore porosity, K inc* and μ inc* is the equivalent bulk modulus and equivalent shear modulus of the spherical inclusion containing soft pores; P (*sp) and Q (*sp) It refers to the geometric factor of soft pores in the equivalent medium; S2.

1. Based on the rock matrix modulus, hard pores are added to the main phase of the background medium, as shown in the following expression: in, and is the initial bulk modulus and initial shear modulus of the background main phase medium; K hp and μ hp are the bulk modulus and shear modulus of the gradually added hard voids; is the initial value of hard pore porosity, K back* and μ back* is the equivalent bulk modulus and equivalent shear modulus of the background main phase medium containing hard pores; P (*hp) and Q (*hp) Refers to the geometric factor of hard pores in the equivalent medium.

5. The inversion-assisted heterogeneous model carbonate rock physics modeling method according to claim 1, characterized in that: In step S3, the modulus of the dry skeleton of the dual-porosity rock is specifically expressed as follows: in, and is the initial bulk modulus and initial shear modulus of the main phase medium with hard pore background; K inc and μ inc are the bulk modulus and shear modulus of the soft pore inclusions added gradually; is the initial value of soft pore content, which is equal to the inclusion content; K dry* and μ dry* is the equivalent bulk modulus and equivalent shear modulus of the dry skeleton of the dual-porosity rock formed after the inclusions are added to the background main phase medium; P (*inc) and Q (*inc) It refers to the geometric factor of spherical inclusions containing soft pores in an equivalent medium.

6. The inversion-assisted heterogeneous model carbonate rock physics modeling method according to claim 1, characterized in that: Step S4 includes the following sub-steps: S4.

1. Estimate the modulus of the mixed fluid based on the modulus and fluid saturation of different fluids. Assume that there is a mixture of two fluids in the rock pores. The bulk modulus and density expressions of the mixed fluid are as follows: r f =s1ρ1+(1-s1)ρ2 Among them, K f is the bulk modulus of the mixed fluid, ρ f is the density of the mixed fluid, K1 and K2 are the bulk moduli of the first and second fluids, ρ1 and ρ2 are the densities of the first and second fluids, and s1 is the saturation of the first fluid; S4.

2. The wave propagation equation of the dual porosity model is described using the Biot-Rayleigh equation. The Christoffel matrix of the Biot-Rayleigh equation is as follows: a 11 =A+2N+i(Q2φ1-Q1φ2)x1,a 12 =Q1+i(Q2φ1-Q1φ2)x2 a 13 =Q2+i(Q2φ1-Q1φ2)x3,a 21 =Q2-iR1φ2x3 a 22 =R1-iR1φ2x2,a 23 =-iR1φ2x3 a 31 =Q2+iR2φ1x1,a 32 =iR2φ1x2,a 33 =R2+iR2φ1x3 b 11 =iω(b1+b2)-ρ 11 oh 2 ,b 12 =-iωb1-ρ 11 oh 2 b 13 =-iωb2-ρ 13 oh 2 ,b 21 =-iωb1-ρ 12 oh 2 b 22 =iωb1-ρ 22 oh 2 ,b 23 =0,b 32 =0 b 31 =-iωb2-ρ 13 oh 2 ,b 33 =-ρ 33 oh 2 +iωb2 Where k is the wave number, a 11 、a 12 、a 13 、a 21 、a 22 、a 23 、a 31 、a 32 、a 33 、b 11 、b 12 、b 13 、b 21 、b 22 、b 23 、b 31 、b 32 、b 33 is the wave equation coefficient, ρ 11 , 12 , 13 is the density parameter of the main phase medium, ρ 22 and ρ 33 is the inclusion density parameter, φ1 is the main phase medium porosity, φ2 is the inclusion porosity, ω is the angular frequency, i is the imaginary sign, x1, x2, x3 are the attenuation factor parameters, A, N, Q1, Q2, R1, R2 are the Biot stiffness coefficients; S4.

3. Solve the above Christoffel matrix to obtain the P-wave and S-wave velocities of the dual-porosity rock saturated with fluid, as follows: in, and are the predicted P-wave velocity and S-wave velocity of fluid-saturated dual-porosity rock, μ dry is the shear modulus of the doubly porous rock, ρ pred represents the predicted density of fluid-saturated dual-porosity rock; S4.

4. Based on the fact that the P-wave velocity and S-wave velocity of fluid-saturated dual-porosity rock are functions of the soft pore content, a forward operator from rock physical parameters and soft pore content to saturated rock velocity information is established as follows: Where G represents the forward operator formed by Biot-Rayleigh dual-porosity rock physics modeling.

7. The inversion-assisted heterogeneous model carbonate rock physics modeling method according to claim 1, characterized in that: Step S5 establishes the inversion objective function as follows: Among them, v sp represents the soft pore content to be inverted, J represents the inversion objective function, and denote the P-wave velocity and S-wave velocity predicted by dual-porosity rock physics modeling, and denote the measured P-wave velocity and S-wave velocity in the well, σ1 and σ2 denote the weight parameters of the P-wave velocity data item and S-wave velocity data item, λ is the regularization parameter, u vSP The expected soft pore content.

8. The inversion-assisted heterogeneous model carbonate rock physics modeling method according to claim 1, characterized in that: Step S6 is specifically: in the first inversion iteration, the search range of the soft pore content is set based on the initial value of the clay content, and the search range is generated by the following expression: After the second inversion iteration, the search range of the current iteration is set based on the soft hole content update result obtained in the previous iteration, which has the following expression: Where k represents the kth inversion iteration, N is the preset maximum number of iterations, and They represent the minimum and maximum values ​​of the soft pore content search range in the kth iteration, a is the adjustment factor of the initial value of the soft pore content, b is the adjustment factor of the soft pore content search range, and v sha is the clay content measured in the well, For the clay content expectation, is the estimated result of soft pore content in the k-1th iteration, is the expectation of the soft pore content estimation result in the k-1th iteration; Based on the clay content in the well, generate initial particles with soft pore content; set M initial particles, particle group It is expressed as follows: in, The initial particle vector representing the soft hole content, represents the ith particle, and the soft pore content of the ith particle to be inverted is obtained by the following formula: Where rand(·) is a random number between [0,1], and c is the fine-tuning scale.

9. The inversion-assisted heterogeneous model carbonate rock physics modeling method according to claim 1, characterized in that: Step S7 includes the following sub-steps: S7.

1. Substituting all particles into the objective function can obtain the fitness of different particles, and selecting the particle with the smallest fitness as the individual optimal particle: When k = 1 When k>1 Where k is the number of iterations, i is the i-th particle, is the soft pore content of the i-th initial particle, represents the individual optimal position of the i-th particle in the first iteration, represents the individual optimal position of the i-th particle in the k-th iteration, v sp,i,k represents the soft hole content of the i-th particle in the k-th iteration. J(v sp,i,k ) and J(v sp,i,k-1 ) represent the fitness of the ith particle in the kth and k-1th iterations, respectively, and are calculated by substituting the particle into the inversion objective function; S7.

2. Calculate the receiving probability according to the fitness of different particles. After obtaining the receiving probability of all particles, use the roulette method to determine the optimal position of the group and the optimal particle of the group. The receiving probability of the particle is calculated as follows: Among them, P accept (v sp,i,k ) represents the receiving probability of the i-th particle in the k-th iteration, represents the optimal position of the group in the kth iteration, represents the fitness of the optimal position of the group in the kth iteration, T k is the temperature value of the kth iteration, τ is a constant, and Σ is the summation symbol; S7.

3. Calculate the update rate based on the individual optimal position and the group optimal position. The expression is as follows: Among them, v sp,i,k represents the i-th particle in the k-th iteration, represents the individual optimal position of the i-th particle in the k-th iteration, represents the optimal position of the group at the kth iteration. p,i,k+1 and r sp,i,k Respectively represent the update rate of the i-th particle in the k+1th and kth iterations. rand1(·) and rand2(·) represent two random numbers between [0,1]. κ1 and κ2 are the update step sizes S7.

4. Based on the update rate, update different particles in the particle swarm as follows: v sp,i,k+1 =v sp,i,k +κ3·r sp,i,k+1 Among them, κ3 is the update step size, v sp,i,k+1 represents the i-th particle in the k+1-th iteration; S7.5, judging whether the number of iterations reaches the preset maximum number of iterations, if yes, taking the updated soft pore content as the inversion result, otherwise executing step S7.6; S7.6, substitute the updated soft pore content into step S1, return to execute steps S1 to S7, continue the inversion iteration until the preset maximum number of iterations is reached, and take the soft pore content output at the maximum number of iterations as the inversion result.

10. The inversion-assisted heterogeneous model carbonate rock physics modeling method according to claim 1, characterized in that: Step S8 specifically includes: taking the inversion result of the soft pore content as known data, inputting it into step S1 to replace the initial value of the soft pore content, repeating steps S1 to S4, and predicting the longitudinal wave velocity, shear wave velocity and density, as shown in the following formula: in, is the inversion result of soft pore content, G represents the forward operator formed by the Biot-Rayleigh dual pore rock physics modeling from step S1 to step S4, To estimate the P-wave velocity using the soft pore content inversion results, To estimate the shear wave velocity using the inversion results of soft pore content, The density is estimated using the inversion results of soft pore content.

Images