Rock physical regularization inversion method based on geological continuity

Abstract

The invention relates to the technical field of oil-gas exploration, in particular to a rock physical regularization inversion method based on geological continuity. The method comprises the following steps: establishing a rock physical model of a target area; based on the rock physical model, establishing a rock physical inversion optimization model; according to the rock physics inversion optimization model, establishing a spatially independent rock physics inversion optimization target generic function, and obtaining an initial geologic model of reservoir physical property parameters; according to the directional derivative of the reservoir physical property parameter geologic model, establishing a physical property parameter geologic model geologic continuity index; based on the inversion optimization model and the geological continuity index, establishing a rock physical inversion optimization model with a geological continuity regular term; and solving the optimal solution of the physical property parameters by adopting a gradient descent method. According to the method, efficient and high-precision prediction of various physical property parameters can be realized, and noise in an inversion result can be effectively suppressed.

Description

Technical Field

[0001] The invention relates to the technical field of oil and gas exploration, and in particular to a rock physics regularization inversion method based on geological continuity. Background Art

[0002] Predicting reservoir properties based on seismic data is a challenging task. The current popular method is to establish a forward model to describe the probability distribution between seismic data and reservoir properties, and then use point estimation or probability simulation algorithms to infer reservoir physical parameters from seismic data. Among them, the process of inferring reservoir physical parameters can be regarded as an inversion problem. At first, multivariate statistical methods were the mainstream methods for establishing forward models and realizing physical 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. To this end, researchers combined mathematical physics methods to establish a rock physics forward model that can characterize physical processes, and then invert physical 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 phase prediction method based on rock physics models and combined with pattern recognition methods. In the same year, Nie et al. introduced genetic algorithms to realize rock physics inversion based on the BISQ model. In 2006, Bachrach combined the Gassmann model and equivalent fluid theory to establish a rock physics model, and jointly inverted porosity and saturation in actual gas fields. 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] In the framework of Bayesian theory, a priori regularization methods are widely used in various inversion problems. In 2010, Bosch et al. summarized a series of a priori regularization methods and their applications in the field of geophysics. 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), and 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 some success in rock physics inversion problems, regularization methods based on geological continuity need to be further developed. Geological continuity refers to the fact that the changes in reservoir media in the underground are strongly spatially correlated. The model commonly used to describe geological continuity is the Markov chain (Ulvmoen and Omre, 2010), but this model is limited to the longitudinal continuity of discrete parameters of the reservoir. Summary of the invention

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

[0006] To solve the above problems, the present invention adopts the following technical solutions:

[0007] The present invention provides a rock physics regularization inversion method based on geological continuity, the method comprising the following steps:

[0008] Step 1. Establish a rock physics model F(l) of the target area, where l represents the physical property parameters to be inverted;

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

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

[0011] Step 4. Establishing geological continuity index of physical property parameter geological model l(x) according to directional derivative of reservoir physical property parameter geological model;

[0012] 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;

[0013] Step 6. Based on the initial geological model in step 3 and the regularized inversion optimization model in step 5, the gradient descent method is used to solve the optimal solution of the physical property parameters.

[0014] Furthermore, the rock physics model F(l) in step 1 and the seismic attributes satisfy the following relationship:

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

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

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

[0018]

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

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

[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 universal function;

[0023] The initial geological model of reservoir physical property parameters is shown as follows:

[0024]

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

[0026]

[0027] Among them, l (k) represents the kth component of the physical property parameter l, λ (k) Represents the physical property parameter component l (k) Geological continuity regularization weight parameter.

[0028] Furthermore, in step 5, the established rock physics inversion optimization model with geological continuity regularization term is as follows:

[0029]

[0030] Where E(l) represents the objective functional of the regularized inverse optimization model, Nl Indicates the type of physical property parameters to be inverted;

[0031] The seismic profile is divided into J rows and I columns, x i,j represents the grid position of row j and column i in the seismic profile, d i,j =d(x i,j ) and l i,j = l(x i,j ) represent the grid positions x i,j The above observation data and reservoir physical parameters, where i = 1, 2, ..., I, j = 1, 2, ..., J represent the vertical and horizontal indicators of the spatial grid respectively. The discretization formula of the rock physics inversion optimization model with geological continuity regularization is:

[0032]

[0033] Where E(L) represents the discretized geological continuity regularized rock physics inversion optimization objective, L = {l i,j} represents the discretized physical parameter geological model.

[0034] Furthermore, the specific steps of using the gradient descent method to solve the optimal solution of the physical property parameters in step 6 are:

[0035] (1) Solve the spatially independent rock physics inversion optimization model point by point to establish the initial model of physical property parameters

[0036]

[0037] (2) Based on the initial geological model of reservoir physical property parameters in step (1), the Jacobian matrix of the rock physics model with respect to the initial model of physical property parameters is calculated:

[0038]

[0039] in, Indicates that at grid position x i,j The initial model of the rock physics model F(l) is about the physical parameters The Jacobian matrix of .

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

[0041]

[0042] in, Indicates that at grid position x i,j Initial model based on physical property parameters Synthetic seismic attributes of and earthquake attribute data d i,j The observation error between .

[0043] (4) Based on the obtained Jacobian matrix and observation error, calculate the gradient descent coefficient A of the regularized inverse optimization model i,j :

[0044] (5) Iteratively calculate the optimal solution of the regularized inverse optimization model:

[0045]

[0046] Among them, the functional gradient Indicates the optimization target at step n At position x i,j The gradient descent direction on the (k) ) Nl×Nl represents the geological continuity regularization weight parameter matrix; represents the Laplace operator; l i,j (n+1) represents the result of the physical property parameter in the nth iteration, γ n Represents the iteration step size.

[0047] Furthermore, the gradient descent coefficient of the regularized inverse optimization model is expressed by the following formula:

[0048]

[0049] Furthermore, Numerical approximation of :

[0050]

[0051] h1 and h2 are used as distance parameters for meshing.

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

[0053] The method of the present invention takes into account the geological continuity of reservoir physical parameters, and the results obtained are closer to the actual exploration situation. The results of practical application show that in oil and gas exploration, the method of this patent can achieve efficient and high-precision prediction of multiple physical parameters and can effectively suppress the noise in the inversion results. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 A flowchart for solving step 6 of the method according to an embodiment of the present invention;

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

[0056] Figure 3 Saturation inversion results: (a) saturation inversion based on rock physics and (b) rock physics regularized inversion using the method of the present invention. DETAILED DESCRIPTION

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

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

[0059] In order to enable those skilled in the art to more clearly understand the technical solution of the present invention, the technical solution of the present invention will be described in detail below in conjunction with specific embodiments.

[0060] Example 1

[0061] The present invention provides a rock physics regularization inversion method based on geological continuity, the method comprising the following steps:

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

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

[0064] Among them, F(l) is the rock physics model, l is the physical property parameter to be inverted, d is the seismic attribute data, and e is 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 and seismic attribute data in step 1, a rock physics inversion optimization model is optimized using least squares at a single grid point;

[0067]

[0068] Where d represents 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 universal function is established to obtain an initial geological model of reservoir physical property parameters;

[0070] Let x∈X represent the position in the geological model, X represents the position 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 the grid position x, respectively. At this time, according to the inversion optimization model in step 2, the spatially independent rock physics inversion optimization objective functional is established as

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

[0072] Where T1(l) represents the spatially independent rock physics inversion optimization objective functional. Then, the initial geological model of reservoir physical parameters is calculated:

[0073]

[0074] Step 4. Establishing geological continuity index of physical property parameter geological model l(x) according to directional derivative of reservoir physical property parameter geological model;

[0075] In order to describe the geological continuity of physical parameters, the reservoir parameter geological model l(x) can be regarded as a multivariate vector function about the spatial position x. In multivariate calculus, the kth component function l in the multidirectional vector function l(x) is (k) The directional derivative of (x) is denoted by Directional derivative Describes the multivariate scalar function l (k) (x) is the local rate of change along the v direction through the position x. In practical applications, the smaller the directional derivative, the better the continuity along a given direction at a fixed position. The maximum square value of all directional derivatives at a fixed position x describes the reservoir parameter l (k) (x) geological continuity, i.e.

[0076]

[0077] Then, the geological model l of the kth reservoir physical property parameters is established (k) (x) geological continuity index DC(l (k) ,X):

[0078]

[0079] Finally, the weight λ is introduced (k) , establish geological model of physical parameters l(x) geological continuity index:

[0080]

[0081] Among them, l (k) represents the kth component of the physical property parameter l, λ (k) Represents the physical 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] Establish a rock physics inversion optimization model with geological continuity regularization term:

[0084]

[0085] Where E(l) represents the objective functional of the regularized inverse optimization model.

[0086] Further, the seismic profile is divided into J rows and I columns, x i,j represents the grid position of row j and column i in the seismic profile, d i,j =d(x i,j ) and l i,j = l(x i,j ) represent the grid positions x i,j The above observation data and reservoir physical parameters, where i = 1, 2, ..., I, j = 1, 2, ..., J represent the vertical and horizontal indicators 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 rock physics inversion optimization objective of the discretized geological continuity regularization, Represents a geological model of discretized physical parameters.

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

[0090] Solve the spatially independent rock physics inversion optimization model point by point:

[0091]

[0092] Establishing the initial model of physical property parameters

[0093] Furthermore, for any physical 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) represents the Jacobian matrix of the rock physics model F(l) at l0:

[0096]

[0097] Among them, N d Indicates the type of earthquake attribute.

[0098] Furthermore, at grid position x i,j Based on the initial model of physical property parameters The rock physics model F(l) can be calculated in The Jacobian matrix at :

[0099]

[0100] At the same time, based on the initial model of physical property parameters Compute the observed error between synthetic seismic attributes and seismic attribute data:

[0101]

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

[0103]

[0104] in,

[0105]

[0106]

[0107] Notice 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 Approximate Gradient

[0109]

[0110]

[0111] Among them, h1 and h2 are used as distance parameters for mesh generation, but they are not necessarily the actual distances.

[0112] Furthermore, the approximate rock physics regularized inversion problem is obtained as follows:

[0113]

[0114] in, Approximate discretization of geological continuity regularization for rock physics inversion optimization objectives.

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

[0116]

[0117] in, represents the target functional At position x i,j The gradient descent direction on , Represents the weight matrix of the regularization function. This represents a numerical approximation to the Laplace operator:

[0118]

[0119] Finally, the gradient descent method is used to solve the rock physics inversion optimization model with geological continuity regularization:

[0120]

[0121] Among them, the functional gradient Indicates the optimization target at step n At position x i,j The gradient descent direction on l i,j (n+1) represents the result of the physical property parameter in the nth iteration, γ n Represents the iteration step size.

[0122] Example 2

[0123] The method described in Example 1 was used to perform inversion using porosity, and the results were as follows: Figure 2 As shown: Figure 2 (a) Using the inversion method of physical property parameters based on rock physics model, the initial geological model of physical property parameters is obtained; then, Figure 2 As shown in (b), based on the rock physics regularization inversion method, the final physical property parameter geological model is obtained.

[0124] Example 3

[0125] Using the method described in Example 1, saturation is used for inversion, and the results are as follows: Figure 3 As shown: Figure 3 (a) Using the inversion method of physical property parameters based on rock physics model, the initial geological model of physical property parameters is obtained; then, Figure 3 As shown in (b), based on the rock physics regularization inversion method, the final physical property parameter geological model is obtained.

[0126] The above embodiments are preferred implementation modes of the present invention, but the implementation modes of the present invention are not limited to the above embodiments. Any other changes, modifications, substitutions, combinations, and simplifications that do not deviate from the spirit and principles of the present invention should be equivalent replacement methods and are included in the protection scope of the present invention.

Claims

1. A rock physics regularization inversion method based on geological continuity, characterized in that: The following steps are involved: Step 1. Establish a rock physics model F(l) of the target area, where l represents the physical property parameters to be inverted; Step 2. Based on the rock physics model F(l), a rock physics inversion optimization model is established; Step 3. According to the rock physics inversion optimization model, a spatially independent rock physics inversion optimization objective universal function is established to obtain an initial geological model of reservoir physical property parameters; Step 4. Establishing geological continuity index of physical property parameter geological model l(x) according to directional derivative of reservoir physical property parameter geological model; 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; Step 6. Based on the initial geological model in step 3 and the regularized inversion optimization model in step 5, the gradient descent method is used to solve the optimal solution of the physical property parameters.

2. The rock physics regularized inversion method based on geological continuity according to claim 1 is characterized in that: In step 1, the rock physics model F(l) and the seismic attributes satisfy the following relationship: d=F(l)+e, Among them, F(l) is the rock physics model, l is the physical property parameter to be inverted, d is the seismic attribute data, and e is the random observation error.

3. The rock physics regularized inversion method based on geological continuity according to claim 1 is characterized in that: In step 2, based on the rock physics model F(l) and seismic attribute data in step 1, a rock physics inversion optimization model is established using the least squares method at a single grid point. The model is shown in the following formula: Where d represents seismic attribute data, W represents the inversion weight matrix, Represents the inner product of vectors e1 and e2 under the weight matrix W.

4. The rock physics regularized inversion method based on geological continuity according to claim 1 is characterized in that: In step 3, x∈X represents the position in the geological model, X represents the position range of the geological model in the target area, d(x) and l(x) represent the seismic attribute data and physical property parameters at the grid position x, respectively; according to the inversion optimization model in step 2, the spatially independent rock physics inversion optimization objective function is established as: T1(l)=∫ x∈X <d-F(l(x)),d-F(l(x))> W right, Where T1(l) represents the spatially independent rock physics inversion optimization objective universal function; The initial geological model of reservoir physical property parameters is shown as follows:

5. The rock physics regularized inversion method based on geological continuity according to claim 1, characterized in that: The geological continuity index of the physical parameter geological model l(x) established in step 4 is as shown in the following formula: Among them, l (k) represents the kth component of the physical property parameter l, λ (k) Represents the physical property parameter component l (k) Geological continuity regularization weight parameter.

6. The rock physics regularized inversion method based on geological continuity according to claim 1, characterized in that: In step 5, the established rock physics inversion optimization model with geological continuity regularization term is as follows: Where E(l) represents the objective functional of the regularized inverse optimization model, N l Indicates the type of physical property parameters to be inverted; The seismic profile is divided into J rows and I columns, x i,j represents the grid position of row j and column i in the seismic profile, d i,j =d(x i,j ) and l i,j = l(x i,j ) represent the grid positions x i,j The above observation data and reservoir physical parameters, where i = 1, 2, ..., I, j = 1, 2, ..., J represent the vertical and horizontal indicators of the spatial grid respectively. The discretization formula of the rock physics inversion optimization model with geological continuity regularization is: Where E(L) represents the discretized geological continuity regularized rock physics inversion optimization objective, L = {l i,j } represents the discretized physical parameter geological model.

7. The rock physics regularized inversion method based on geological continuity according to claim 1, characterized in that: The specific steps of using the gradient descent method to solve the optimal solution of the physical property parameters in step 6 are: (1) Solve the spatially independent rock physics inversion optimization model point by point to establish the initial model of physical property parameters (2) Based on the initial geological model of reservoir physical property parameters in step (1), the Jacobian matrix of the rock physics model with respect to the initial model of physical property parameters is calculated: in, Indicates that at grid position x i,j The initial model of the rock physics model F(l) is about the physical parameters The Jacobian matrix of . (3) Based on the initial geological model of reservoir physical property parameters in step (1), the observation error between the synthetic seismic attribute and the seismic attribute data is calculated: in, Indicates that at grid position x i,j Initial model based on physical property parameters Synthetic seismic attributes of and earthquake attribute data d i,j The observation error between . (4) Based on the obtained Jacobian matrix and observation error, calculate the gradient descent coefficient A of the regularized inverse optimization model i,j : (5) Iteratively calculate the optimal solution of the regularized inverse optimization model: Among them, the functional gradient Indicates the optimization target at step n At position x i,j The gradient descent direction on , represents the geological continuity regularization weight parameter matrix; represents the Laplace operator; l i,j (n+1) represents the result of the physical property parameter in the nth iteration, γ n Represents the iteration step size.

8. The rock physics regularized inversion method based on geological continuity according to claim 7 is characterized in that: The gradient descent coefficient of the regularized inverse optimization model is expressed by the following formula:

9. The rock physics regularized inversion method based on geological continuity according to claim 7, characterized in that: Numerical approximation of : h1 and h2 are used as distance parameters for meshing.

Images