A rock physics regularization inversion method based on geological continuity

Abstract

The present application relates to the oil and gas exploration technical field, specifically relates to a kind of rock physics regularization inversion method based on geological continuity.The method includes: the rock physics model of target area is established;Rock physics inversion optimization model is established based on the rock physics model;According to the rock physics inversion optimization model, the spatial independent rock physics inversion optimization objective function is established, and the initial geological model of reservoir physical property parameter is obtained;According to the directional derivative of reservoir physical property parameter geological model, the geological continuity index of physical property parameter geological model is established;Rock physics inversion optimization model with geological continuity regularization term is established based on the inversion optimization model and geological continuity index;The optimal solution of physical property parameter is solved using gradient descent method.The method of the present application can realize the efficient and high-precision prediction of multiple physical property parameters, and can effectively suppress the noise in the inversion result.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of oil and gas exploration, and particularly relates to a rock physics regularization inversion method based on geological continuity. BACKGROUND

[0002] It is a challenging work to predict reservoir properties based on seismic data. The popular method is to build a forward model to describe the probability distribution between seismic data and reservoir properties, and then use point estimation or probabilistic simulation algorithm to infer reservoir property parameters from seismic data. Among them, the process of inferring reservoir property parameters can be regarded as an inversion problem. Initially, multivariate statistical methods were the mainstream method to build forward models and realize property parameter prediction (Doyen, 1988; Fournier, 1989; Lucet and Mavko, 1991). However, multivariate statistical methods are simple statistical analysis of specific data, and often lack universality. Therefore, researchers combined mathematical physics methods to build rock physics forward models that can describe physical processes, and then invert property parameters. In 2001, Mukerji et al. proposed a statistical rock physics model and realized the inversion prediction of porosity based on seismic elastic attributes. In 2004, Eidsvik et al. proposed a reservoir facies prediction method based on rock physics model combined with pattern recognition method. In the same year, Nie et al. introduced genetic algorithm to realize rock physics inversion based on BISQ model. In 2006, Bachrach combined Gassmann model and equivalent fluid theory to build rock physics model, and jointly inverted porosity and saturation in actual gas field. In recent years, more accurate rock physics models and more advanced statistical calculation methods have been introduced into the field of rock physics inversion (Fang and Yang, 2015; Wen et al., 2022).

[0003] Under the framework of Bayesian theory, prior regularization methods are widely used in various inversion problems. In 2010, Bosch et al. summarized a series of prior regularization methods and their applications in geophysical field. Among them, regularization methods can be used to improve the accuracy of seismic inversion (Grana and Rossa, 2010; and Berkhout, 1992; Tarantola, 2005), alleviate the multi-solution problem of nonlinear rock physics inversion (Wen et al., 2022; Grana, 2018), add constraints such as geological continuity and lithofacies correlation (Escobar et al., 2006; de Figueiredo et al., 2018; Li Kun et al., 2020).

[0004] Although some regularization methods have achieved certain success in petrophysical inversion problems, the regularization method based on geological continuity needs to be further developed. The geological continuity refers to the change of the reservoir medium underground is of strong spatial correlation. The model commonly used to describe the geological continuity is Markov chain (Ulvmoen and Omre, 2010), but this model is limited to the longitudinal continuity of discrete parameters of the reservoir. SUMMARY

[0005] The main purpose of the present application is to provide a petrophysical inversion regularization method based on geological continuity. In order to describe the transverse and longitudinal continuity of continuous parameters such as porosity and saturation, the present application quantifies the geological continuity of continuous parameters by using directional derivative, and develops a petrophysical regularization inversion method based on the above geological continuity. The present application realizes efficient and high-precision prediction of physical parameters such as reservoir porosity.

[0006] To solve the above problems, the present application adopts the following technical scheme:

[0007] The present application provides a petrophysical regularization inversion method based on geological continuity, which comprises the following steps:

[0008] Step 1. Establishing a petrophysical model F(l) of a target area, wherein l represents a physical parameter to be inverted;

[0009] Step 2. Establishing a petrophysical inversion optimization model based on the petrophysical model F(l);

[0010] Step 3. According to the petrophysical inversion optimization model, establishing a spatially independent petrophysical inversion optimization target generating function to obtain an initial geological model of the reservoir physical parameter;

[0011] Step 4. According to the directional derivative of the geological model of the reservoir physical parameter, establishing a geological continuity index of the physical parameter geological model l(x);

[0012] Step 5. Based on the inversion optimization model in step 2 and the geological continuity index in step 4, establishing a petrophysical inversion optimization model with a geological continuity regularization term;

[0013] Step 6. Based on the initial geological model in step 3 and the regularization inversion optimization model in step 5, using gradient descent method to solve the optimal solution of the physical parameter.

[0014] Further, the petrophysical model F(l) in step 1 satisfies the following relationship with the seismic attribute:

[0015] d=F(l)+e,

[0016] Where F(l) is the rock physics model, l represents the physical property parameters to be inverted, d represents the seismic attribute data, and e represents the random observation error.

[0017] Furthermore, in step 2, based on the rock physics model F(l) and seismic attribute data from step 1, a rock physics inversion optimization model is established at a single grid point using the least squares method. The model is shown in the following formula:

[0018]

[0019] Where d represents seismic attribute data, and W represents the inversion weight matrix. This represents the inner product of vectors e1 and e2 under the weight matrix W.

[0020] Further, in step 3, x∈X represents the location in the geological model, X represents the location range of the geological model in the target area, and d(x) and l(x) represent the seismic attribute data and physical property parameters at grid location x, respectively; based on the inversion optimization model in step 2, the spatially independent rock physics inversion optimization objective function is established as follows:

[0021] T1(l)=∫ x∈X dF(l(x)),dF(l(x)) W dx,

[0022] Where T1(l) represents the spatially independent rock physics inversion optimization objective function;

[0023] The initial geological model for the obtained reservoir physical property parameters is shown in the following equation:

[0024]

[0025] Furthermore, the geological continuity index of the physical property parameter geological model l(x) established in step 4 is shown in the following formula:

[0026]

[0027] Among them, l (k) Let λ represent the k-th component of the physical property parameter l. (k) Represents the component l of the physical property parameter (k) Geological continuity regularization weight parameters.

[0028] Furthermore, in step 5, the established rock physics inversion optimization model with a geological continuity regularization term is shown in the following equation:

[0029]

[0030] Where E(l) represents the objective functional of the regularized inversion optimization model, and Nl denotes the type of the petrophysical parameter to be inverted;

[0031] The seismic profile is divided into J rows and I columns, x i,j denotes the grid position of the seismic profile j row i column, d i,j = d(x i,j ) and l i,j = l(x i,j ) respectively denote the observed data and the reservoir petrophysical parameter at grid position x i,j , where i = 1, 2,..., I and j = 1, 2,..., J respectively denote the spatial grid longitudinal and lateral indices, and the discretized formula of the geologically continuous regularization petrophysical inversion optimization model is:

[0032]

[0033] where E(L) denotes the discretized geologically continuous regularization petrophysical inversion optimization objective, L = {l i,j} denotes the discretized petrophysical parameter geological model.

[0034] Further, the specific steps of solving the optimal solution of the petrophysical parameter by using the gradient descent method in step 6 are:

[0035] (1) Solving the spatial independent petrophysical inversion optimization model point by point to establish the initial model of the petrophysical parameter

[0036]

[0037] (2) Based on the initial geological model of the reservoir petrophysical parameter described in step (1), the Jacobian matrix of the petrophysical model about the initial model of the petrophysical parameter is calculated:

[0038]

[0039] where F i,j (l ) denotes the Jacobian matrix of the petrophysical model F(l) about the initial model of the petrophysical parameter l at grid position x .

[0040] (3) Based on the initial geological model of the reservoir petrophysical parameter described in step (1), the observation error between the synthetic seismic attribute and the seismic attribute data is calculated:

[0041]

[0042] where F i,j (l ) denotes the synthetic seismic attribute based on the initial model of the petrophysical parameter l at grid position x and seismic attribute data d i,j between the observation error.

[0043] (4) Based on the obtained Jacobian matrix and observation error, the gradient descent coefficient A of the regularization inversion optimization model is calculated i,j :

[0044] (5) The optimal solution of the regularization inversion optimization model is iteratively calculated:

[0045]

[0046] wherein the functional gradient represents the gradient descent direction of the optimization objective on the position x i,j , Λ = diag (2λ (k) ) Nl×Nl represents the geological continuity regularization weight parameter matrix; represents the Laplacian operator; l i,j (n+1) represents the result of the physical property parameter at the n th iteration, γ n represents the iteration step length.

[0047] Further, the gradient descent coefficient of the regularization inversion optimization model is represented by the following formula:

[0048]

[0049] Further, the numerical approximation of

[0050]

[0051] h1 and h2 are distance parameters of grid division.

[0052] Compared with the prior art, the present application has the following advantages:

[0053] The method of the present application considers the geological continuity of the reservoir physical property parameters, and the obtained result is closer to the actual exploration situation. The results of practical application show that in the aspect of oil and gas exploration, the method of the present application can realize efficient and high-precision prediction of multiple physical property parameters, and can effectively suppress noise in the inversion result. BRIEF DESCRIPTION OF DRAWINGS

[0054] Figure 1 the flow chart for solving step 6 of the method described in the embodiments of the present application;

[0055] Figure 2 are the porosity inversion results: (a) porosity inversion based on rock physics model and (b) rock physics regularization inversion of the method described in the present application;

[0056] Figure 3 (a) saturation inversion based on rock physics and (b) rock physics regularized inversion of the method described herein. DETAILED DESCRIPTION

[0057] It should be noted that the following detailed description is exemplary in nature and is intended to provide further description of the application. Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs.

[0058] It is to be understood that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of example embodiments in accordance with the present application. As used herein, the singular forms "a", "an" and "the" are intended to include the plural forms as well, unless the context clearly indicates otherwise. It will be further understood that the terms "comprises" and / or "comprising," when used in this specification, specify the presence of stated features, steps, operations, elements, and / or components, but do not preclude the presence or addition of one or more other features, steps, operations, elements, components, and / or groups thereof.

[0059] In order to enable persons skilled in the art to more clearly understand the technical solutions of the present application, the technical solutions of the present application will be described in detail below with specific examples.

[0060] Example 1

[0061] The present application provides a rock physics regularized inversion method based on geological continuity, which comprises the following steps:

[0062] Step 1. Establish a rock physics model F(l) of the target area, where l represents the physical property parameter to be inverted; the model represents the basic relationship between the physical property parameter to be inverted and the seismic attribute, and satisfies

[0063] d = F(l) + e,

[0064] Where F(l) is the rock physics model, l represents the physical property parameter to be inverted, d represents the seismic attribute data, and e represents the random observation error.

[0065] Step 2. Based on the rock physics model F(l), a rock physics inversion optimization model is established;

[0066] Based on the rock physics model in step 1 and the seismic attribute data, a rock physics inversion optimization model using least squares is used at a single grid point;

[0067]

[0068] Where d represents the seismic attribute data, W represents the inversion weight matrix, represents the inner product of vectors e1 and e2 under the weight matrix W.

[0069] Step 3. According to the rock physics inversion optimization model, a spatially independent rock physics inversion optimization objective function is established to obtain an initial geology model of reservoir physical parameters;

[0070] Let x e X represent a position in the geology model, X represent the position range of the target area geology model, d(x) and l(x) represent the seismic attribute data and physical parameters at grid position x, respectively. At this time, according to the inversion optimization model in step 2, the spatially independent rock physics inversion optimization objective function is established as

[0071] T1(l) = ∫ x∈X d-F(l(x)), d-F(l(x)) W dx,

[0072] wherein T1(l) represents the spatially independent rock physics inversion optimization objective function. Further, an initial geology model of reservoir physical parameters is calculated as

[0073]

[0074] Step 4. According to the directional derivative of the geology model of reservoir physical parameters, a geology continuity index of the physical parameter geology model l(x) is established.

[0075] In order to describe the geology continuity of the physical parameters, the reservoir parameter geology model l(x) can be regarded as a multi-element vector function with respect to the spatial position x. In multi-element calculus, the directional derivative of the kth component function l (k) (x) in the multi-element vector function l(x) is denoted as The directional derivative describes the local change rate of the multi-element scalar function l (k) (x) along the v direction through the position x. In practical applications, the smaller the directional derivative is, the better the continuity along the given direction at the fixed position is. The maximum square value of all directional derivatives at the fixed position x describes the geology continuity of the reservoir parameter l (k) (x), that is,

[0076]

[0077] Further, the geology continuity index DC(l (k) (x) of the kth reservoir physical parameter geology model l (k) (x) is established as

[0078]

[0079] Finally, a weight λ (k) is introduced to establish the geology continuity index of the physical parameter geology model l(x) as

[0080]

[0081] where, l (k) represents the kth component of the property parameter l, λ (k) represents the property parameter component l (k) geological continuity regularization weight parameter.

[0082] Step 5. Based on the inversion optimization model in step 2 and the geological continuity index in step 4, a rock physics inversion optimization model with a geological continuity regularization term is established;

[0083] The rock physics inversion optimization model with a geological continuity regularization term is established as:

[0084]

[0085] where, E(l) represents the objective functional of the regularization inversion optimization model.

[0086] Further, the seismic profile is divided into J rows and I columns, x i,j represents the grid position of the jth row and ith column of the seismic profile, d i,j = d(x i,j ) and l i,j = l(x i,j ) respectively represent the observation data and the reservoir property parameter at the grid position x i,j , where i = 1, 2,..., I and j = 1, 2,..., J represent the longitudinal and transverse indices of the spatial grid respectively. At this time, the discretization formula of the rock physics inversion optimization model with geological continuity regularization is:

[0087]

[0088] where, E(L) represents the discretized rock physics inversion optimization objective with geological continuity regularization, represents the discretized property parameter geological model.

[0089] Step 6. Based on the initial geological model in step 3 and the regularization inversion optimization model in step 5, the gradient descent method is used to solve the optimal solution of the property parameter, as shown in Figure 1 . The specific steps are:

[0090] Point-by-point solving of the spatially independent rock physics inversion optimization model:

[0091]

[0092] Establishing the initial model of the property parameter

[0093] Further, for any property parameter l0, the rock physics model F(l) can be linearized as:

[0094] F(l) ≈ F(l0) + J(l0)(l - l0),

[0095] where J(l0) denotes the Jacobian matrix of the rock physics model F(l) at l0:

[0096]

[0097] where N d denotes the number of seismic attributes.

[0098] Further, at the grid location x i,j , the Jacobian matrix of the rock physics model F(l) at l0 can be computed based on the petrophysical parameter initial model

[0099]

[0100] Meanwhile, the observation error between the synthetic seismic attributes and the seismic attribute data can be computed based on the petrophysical parameter initial model

[0101]

[0102] At this time, the first term of the optimization model can be expanded as:

[0103]

[0104] where

[0105]

[0106]

[0107] Note that is a constant that can be ignored in the inversion process.

[0108] In the second term of the optimization model, the forward difference quotient is used to approximate the gradient

[0109]

[0110]

[0111] where hi and h2 are distance parameters of the grid partition, but not necessarily the true distance.

[0112] Further, the approximated rock physics regularized inversion problem is obtained as:

[0113] ​​​

[0114] where, Approximated discretized geologic continuity regularized petrophysical inversion optimization objective.

[0115] At this time, the objective functional is calculated according to the variational method The gradient of the L functional is as follows:

[0116]

[0117] where, denotes the objective functional The gradient descent direction at position x i,j , the weight matrix of the regularization function. denotes the numerical approximation of the Laplace operator:

[0118]

[0119] Finally, the geologic continuity regularized petrophysical inversion optimization model is solved by using the gradient descent method:

[0120]

[0121] where, the functional gradient denotes the gradient descent direction of the optimization objective at position x i,j , l i,j (n+1) denotes the result of the physical property parameter at the n th iteration, γ n denotes the iteration step length.

[0122] Example 2

[0123] Using the method described in Example 1, porosity is used for inversion, and the results are as shown in Figure 2 : Figure 2 (a) Using the physical property parameter inversion method based on the petrophysical model, an initial geologic model of the physical property parameter is obtained; further, as shown in Figure 2 (b), based on the petrophysical regularization inversion method, a final geologic model of the physical property parameter is obtained.

[0124] Example 3

[0125] Using the method described in Example 1, saturation is used for inversion, and the results are as shown in Figure 3 : Figure 3 (a) Using the physical property parameter inversion method based on the petrophysical model, an initial geologic model of the physical property parameter is obtained; further, as shown in Figure 3 (b), based on the petrophysical regularization inversion method, a final geologic model of the physical property parameter is obtained. ​

[0126] The above embodiments are the preferred embodiments of the present application, but the embodiments of the present application are not limited to the above embodiments, and any changes, modifications, substitutions, combinations, simplifications, etc. made without departing from the spirit and principles of the present application should be equivalent replacement manners and should be included in the protection scope of the present application.

Claims

1. A geologically consistent rock physics regularization inversion method, characterized in that, The method comprises the following steps: Step 1. Establishing a petrophysical model F(l) of a target area, wherein l represents a physical property parameter to be inverted; Step 2. Establishing a petrophysical inversion optimization model based on the petrophysical model F(l); Step 3. According to the petrophysical inversion optimization model, a spatially independent petrophysical inversion optimization objective function is established to obtain an initial geological model of the reservoir physical property parameter; Step 4. According to the directional derivative of the reservoir physical property geological model, a geological continuity index of the physical property geological model l(x) is established; Step 5. Based on the inversion optimization model in step 2 and the geological continuity index in step 4, a petrophysical inversion optimization model with a geological continuity regularization term is established; Step 6. Based on the initial geological model in step 3 and the regularization inversion optimization model in step 5, a gradient descent method is used to solve the optimal solution of the physical property parameter; The specific steps of solving the optimal solution of the physical property parameter by the gradient descent method in step 6 are as follows: (1) Point by point solution of space-independent rock physics inversion optimization model, to establish the initial model of physical parameters The seismic profile is divided into J rows and I columns, x i,j represents the grid position of the seismic profile j row i column, d i,j = d(x i,j ) and l i,j = l(x i,j ) represent the observation data and reservoir property parameters at the grid position x i,j , respectively; (2) Based on the initial geological model of the reservoir physical property parameter in step (1), the Jacobian matrix of the petrophysical model with respect to the initial model of the physical property parameter is calculated: wherein, represents the Jacobian matrix of the petrophysical model F(l) with respect to the initial model of the petrophysical properties i,j at the grid position x ; (3) Based on the initial geological model of the reservoir physical property parameter in step (1), the observation error between the synthetic seismic attribute and the observation data is calculated: wherein, represents the observed error between the synthetic seismic attributes i,j based on the initial model of physical properties at the grid locations x and the observed data d i,j ; (4) Based on the resulting Jacobian matrix and the observation errors, compute the gradient descent coefficients A of the regularization inversion optimization model i,j : (5) The optimal solution of the regularization inversion optimization model is iteratively calculated: where the functional gradient denotes the optimization objective at the n-th step denotes the gradient descent direction at position x i,j denotes the geological continuity regularization weight parameter matrix; denotes the Laplacian operator; l i,j (n+1) denotes the result of the property parameter at the n-th iteration, γ n denotes the iteration step size.​ 2. The geologically consistent rock physics regularization inversion method of claim 1, wherein, In step 1, the petrophysical model F(l) and the seismic attribute satisfy the following relationship: d = F(l) + e, Wherein F(l) is the petrophysical model, l represents the physical property parameter to be inverted, d represents the seismic attribute data, and e represents the random observation error.

3. The geologically consistent rock physics regularization inversion method of claim 2, wherein, In step 2, based on the petrophysical model F(l) and the seismic attribute data in step 1, a petrophysical inversion optimization model is established at a single grid point by using the least square method, and the model is shown in the following formula: where d denotes seismic attribute data, W denotes an inversion weight matrix, denotes the inner product of the vectors e1 and e2 under the weight matrix W.

4. The petrophysical regularization inversion method based on geological continuity according to claim 3, characterized in that, In step 3, x∈X represents the position of the geological model, X represents the position range of the geological model of the target area, d(x) and l(x) represent the seismic attribute data and the physical property parameter at the grid position x respectively; according to the inversion optimization model in step 2, a spatially independent petrophysical inversion optimization objective function is established as follows: T1(l) = ∫ x∈X < d-F(l(x)), d-F(l(x)) > W dx, Wherein T1(l) represents the spatially independent petrophysical inversion optimization objective function; The obtained initial geological model of the reservoir physical property parameter is shown in the following formula:

5. The petrophysical regularization inversion method based on geological continuity according to claim 4, characterized in that, In step 4, the geological continuity index of the physical property geological model l(x) is established as shown in the following formula: wherein, l (k) denotes the k-th component of the physical property parameter l, λ (k) denotes the physical property parameter component l (k) geological continuity regularization weight parameter.

6. The petrophysical regularization inversion method based on geological continuity according to claim 5, characterized in that, In step 5, the petrophysical inversion optimization model with the geological continuity regularization term is established as shown in the following formula: where E(l) represents the objective functional of the regularization inversion optimization model, N l denotes the type of the physical property to be inverted; The discretization formula of the geological continuity regularization petrophysical inversion optimization model is as follows: where E(L) represents the discretized geologic continuity regularized petrophysical inversion optimization objective, L = {l i,j} represents the discretized petrophysical parameter geologic model.

7. The method of geologically consistent rock physics regularization inversion according to claim 1, wherein, The gradient descent coefficient of the regularization inversion optimization model is represented by the following formula:

8. The geologically consistent rock physics regularization inversion method of claim 1, wherein, The numerical approximation is: Wherein h1 and h2 are distance parameters of grid division.

Images