WO2016083861A1

WO2016083861A1 - Mapping a geological parameter on an unstructured grid

Info

Publication number: WO2016083861A1
Application number: PCT/IB2014/002908
Authority: WO
Inventors: Victor ZAYTSEV; Pierre Biver; Hans Wackernagel; Denis Allard
Original assignee: Total Sa; Association Pour La Recherche Et Le Developpement Des Methodes Et Processus Industriels
Priority date: 2014-11-27
Filing date: 2014-11-27
Publication date: 2016-06-02

Abstract

The present invention relates to a method for a determination of a geological parameter in an unstructured grid representing a geological model. The method comprises receiving a description of a random field (Z(x), 101 ) of a petrophysical/geological variable and deriving a description of an auxiliary (102) Gaussian random field (Y(x)) through a transformation function. Then it is possible to determine (103) a decomposition in an Hilbertian basis of the transformation function, and for each block v_p of the unstructured grid, a change of support coefficient r_p by solving (a) being an input covariance function of the Gaussian random field or (b) being an input covariance function of the received random field. Using this block-specific change of support coefficients, it is possible to perform a mapping of the average values of the geological parameter over the blocks of the unstructured grid, honoring the desired structure of spatial distribution for the said parameter.

Description

MAPPING A GEOLOGICAL PARAMETER ON AN UNSTRUCTURED GRID

BACKGROUND OF THE INVENTION

The present invention relates to the representation of real geological subsoil and in particular for determination of a petrophysical/geological variable on an unstructured grid.

The approaches described in this section could be pursued, but are not necessarily approaches that have been previously conceived or pursued. Therefore, unless otherwise indicated herein, the approaches described in this section are not prior art to the claims in this application and are not admitted to be prior art by inclusion in this section. Furthermore, all embodiments are not necessarily intended to solve all or even any of the problems brought forward in this section.

Most of the geological models are represented with a regular grid. Indeed, the tools used for the determination of petrophysical/geological variables are not efficient on unstructured grids.

Nevertheless, the unstructured grids have the ability to adapt the size and the shape of the blocks depending of the interest of the various zones of the model: for instance, coarse blocks (convex polygons) may be used on zones that are of less importance for the hydrocarbons determination and fine blocks may be used on zones that are close to the region of interest. Unstructured grids have a more flexible geometry.

When prescribing petrophysical/geological variables on a geological model, the usual approach to populate unstructured grids with petrophysical variables is to simulate the petrophysical/geological variables on a fine regular grid and further to up-scale the simulated property to the target unstructured grid.

Nevertheless, as the auxiliary grid has to have very fine blocks, such an approach require many computations. There is thus a need for a method to determine petrophysical/geological variables directly on an unstructured grid and to limit the number of computations needed.

SUMMARY OF THE INVENTION

The invention relates to a method for mapping/determining a geological parameter in an unstructured grid representing a geological model.

The method comprises:

- receiving the unstructured grid, the unstructured grid having a plurality of blocks v_p ;

- receiving a distribution of a random field Z(x) of a petrophysical/geological variable;

- determining/deriving a transformation function to convert the random field Z(x) to a Gaussian random field Y(x), the transformation function being defined by the received distribution of the random field Z(x) and by a distribution of the Gaussian random field;

- determining a decomposition in an Hilbertian basis of the transformation function, {ψι, ί = 1 ... N_∞} being the coordinates of the transformation function in said Hilbertian basis;

- receiving:

- an input covariance function C(h) of the received random field; or

- an input covariance function p(h) of the Gaussian random field;

- for each block of the unstructured grid, determining the geological parameter based on the Gaussian random field, the change of support coefficient associated with said block and the coordinates of the transformation function in said Hilbertian basis.

Indeed, when determining a geological parameter on an unstructured grid, it is advantageous for modeling petrophysics to modify dynamically the distribution for each block. A key parameter for modifying the distribution is the change of support coefficient r_p for each block, said change of support coefficient may be different for each block.

The distribution of a random field Z(x) may be characterized by a marginal/univariate distribution of said random field and the covariance function of the said random field. The distribution of the Gaussian random field Y(x) may also be characterized by its marginal distribution and its covariance function. Determining/deriving a transformation function may also comprise deriving a covariance function p i) of the Gaussian random field Y(x). Advantageously, the determined geological parameter is a geological parameter that may be linearly averaged such as porosity. For the geological parameters such as permeability, which are not averaged linearly over the block, one who is experienced in the subject can use a transformation to a parameter which is linearly averaged over the block (classical transformation being power transformation). The conversion of the random field into a Gaussian space is performed as most of the common geostatistical tools are in Gaussian space (such as Kriging for instance).

The use of normalized Hermite polynomials as a Hilbertian basis is advantageous as they are compatible with Gaussian distributions.

An "unstructured grid" is a grid representing a geological model with a non-regular grid. Thus, the size and the shape of the blocks (or cells) in the grid are not always the same and do not follow a construction rule (except that the blocks are convex).

The decomposition in Hilbertian basis is not always an exact decomposition. Indeed, as only a limited number of coordinates are determined (the Hilbertian basis may be truncated and, thus, has only a limited number of members), there is no guarantee that the decomposition is exact. Thus, in the general case, the derived decomposition defines an approximation to the original function.

As Y(x) is a multivariate Gaussian random variable (N(0,1)), the multi-Gaussian assumption for Z(x) indicates that Z(x) = φ(Ύ(χ)) with φ Gaussian anamorphosis and φ^_1 the Gaussian transform function.

In addition, determining the geological parameter may comprise the computation of:

wherein the Hilbertian basis is noted { χι, ί = 1 ... N_∞}.

In one possible embodiment, determining the geological parameter may further comprise:

- for each couple of blocks v_p and v_q in the plurality of blocks respectively associated with r_p and r_q as change of support coefficients, computing a block-to-block covariance R_pq of an auxiliary Gaussian random vector by solving :

R_va = ——;—;— ί ί p (x. x')dxdx' ; or <

rpi- r_p p ₀ q¾ pq .

Furthermore, determining the geological parameter may further comprise: - for each current block of the unstructured grid, determining a realization of a

Gaussian random variable, said Gaussian random variable being converted from a random variable representing an average value of the received random field in said current block. The method may also comprise:

- receiving a set of conditioning data for points of the geological model;

- transforming the received conditioning data into Gaussian conditioning data by applying the inverse φ^_1 of the transformation function or the inverse of the decomposition in the Hilbertian basis of the transformation function φ;

- for each current point of the geological model where a conditioning data is received, determining a random value function of :

Y^NCS(v_a) + T^r_a or y — yNCS , / 1 _ _r2 T

I a — ^la ^Iv_a ^T V 'a ^l a wherein T_a are mutually independent standard Gaussian variables, v_a is the block containing said current point,r_a is the change of support coefficient determined for said block v_a, Y^NCS{v_a) is an average value of the Gaussian random variable Y(x) over block v_a and Y^^s is a Gaussian random variable (or an unconditional simulation of a Gaussian random variable) being converted from a random variable representing an average value of the received random field in the block v_a.

A second aspect of the invention relates to a computer program product comprising a computer readable medium, having thereon a computer program comprising program instructions. The computer program is loadable into a data- processing unit and adapted to cause the data-processing unit to carry out the method described above when the computer program is run by the data-processing unit.

Other features and advantages of the method and apparatus disclosed herein will become apparent from the following description of non-limiting embodiments, with reference to the appended drawings. BRIEF DESCRIPTION OF THE DRAWINGS

The present invention is illustrated by way of example, and not by way of limitation, in the figures of the accompanying drawings, in which like reference numerals refer to similar elements and in which:

- Figure 1 is a flow chart describing a possible embodiment of the present invention;

- Figure 2 is a possible embodiment for a device that enables the present invention.

DESCRIPTION OF PREFERRED EMBODIMENTS

Figure 1 is a flow chart describing a possible embodiment of the present invention.

The space coordinates of the geological model (106) being denoted as (x,y, z), x or x (for simplified notations), each convex block v <= R³ of an unstructured grid of said model may be analyzed and be represented (step 107) through a union of linear inequalities v = u_{i=1 n(V)} l^v) with:

d\ + e\x + fly≤ z≤ d\ + elx + f[y where ... , fl, f{ are constant parameters to be determined.

This representation enables computing the block-to block and block-to-point covariance with any target accuracy as well as to verify if a given point is contained in a given block.

In particular geometrical cases (such as parallelogram blocks, the representation through inequalities can be derived in an optimized manner).

In the general case, it is possible to use a general algorithm for representing a convex body as a set of non-intersecting systems of inequalities (Kohenblit & Shmerling, 2006) and thus to determine α-, α[, ... , fi , f[.

Any other representations of blocks are possible. Further, the unstructured grid may be represented as an optimized k-d tree (Friedman H., et al., 1 977) using the block centroids as tree terminal nodes. Using the k-d tree representation enables one to perform neighborhood searches in proven logarithmic time.

It is possible to receive a description of the target point-support petrophysical/geological variable Z(x) (a distribution of Z(x) , 1 01 and a covariance, 1 08) with a marginal distribution F (F(z) = E(Z(x) < z) , F is strictly monotonic and defines finite mean and variance and is considered as an input), x ε D c R³.

It is possible to describe (step 102) Y(x) (i G fl c R³), a distribution of a Gaussian field corresponding to the field Z(x) (i.e. to describe Y(x) which corresponds to Z(x) , see Geostatistics: Modeling Spatial Uncertainty Second Edition, J-P Chiles, P.Delfiner, 201 2, p.391 ). Thus, Y(x) is considered to have a multivariate Gaussian distribution with standard Gaussian marginal distribution with mean 0 and variance 1 . The covariance function p(h) of Y x) and the covariance function C(/i) of Z(x) are linked through the decomposition of the transformation function φ(. ) through the formula

This formula enables to derive p(h) in case if C(h) is given as input, and C(h) in case if p(h) is given as input.

The fields Y(x) and Z(x) are related through a transform function φ{. ) Z(x) = φ(Υ(χ)) which is fully determined by the marginal distribution F of Z(x) .

Transformation function <p(y) is equal to F o G^_1(y), where G(y) is the cumulative density function (CDF) of the Gaussian variable with 0 mean and 1 variance (and G^_1 is the inverse function of G), G(y) = -= §^y e^~dx.

In other words, <p(y) is the transformation function from the Gaussian space to the 'real' space.

In particular, it is possible to approximate (step 1 03) the transform function φ with a basis of normalized Hermitian polynomials te(y), i = 0 ...∞} φ(γ) = ∑?=_Q (PiXi(y), where {<p_it i = 1 ...∞} are the coordinates of φ in an Hilbertian basis to(y), i = 0 ...∞}.

For computational purpose, it may be difficult/impossible to determine an infinite number of coordinates, thus, it is possible to only compute the first N_∞ coordinates (with N_∞ a predetermined number of coordinates):

<p( ) = ^ <M(y)

£=0

If an input set of conditioning data {z_a, a = 1 ... n} (104) is received at fixed points {x_a, a = l ... n}, it is possible to transform the conditioning data {z_a, a = l ... n} into Gaussian data {y_a, a = 1 ... n} with the normal score transform function φ^'1 (step 105) :

When receiving a covariance function (108), it is possible to determine a change of support coefficient r_p for every block v_p (step 109).

The covariance function may for instance be derived from seismic data or well cores data. This change of support coefficient r_p may be computed (step 109) using the formula below, if p(h) is received (p(/i) being the covariance function of the random variable Y(x), defined on point support):

If C(h) (the covariance function of the variable Z(x) is received, defined on point support), it is possible to compute instead the variance of Z_Vp (the random variable representing the average value of Z(x) in the block v_p) :

Var 2i

The formula above is a polynomial equation for r_p, which can be solved (step 109) (i.e. inversed) with standard methods, e.g. Newton's method.

Once the change of support coefficient r_p is computed for every block v_p, it is possible to derive (step 1 10) the block transform function <p_Vp for every block by computing:

Given an input covariance function p(h) or C(h) and an unstructured grid of N blocks {Vp c D, p = 1 ...N}, a block-to-block covariance matrix of size N x N may be constructed, as described below.

In order to minimize the memory consumption of the algorithm the sparse storage format for a symmetric triangular matrix may advantageously be used and only the non-zero elements of the matrix are stored.

The construction of the matrix described below is only a possible embodiment. In particular, it is possible to use any other representation or to compute the values of the matrix whenever needed.

The matrix construction (step 1 1 1 ) may thus proceed as follows:

- represent the input grid {v_{i t} i = 1 ... N} as a k-d tree (or any other data structure capable for performing neighborhood search in an optimized way)

- for each block v_t of the grid, use the above-mentioned tree to determine the block neighborhood;

- For every identified block in the neighborhood, compute the block-to-block covariance R_pq for block v_p and v_q (see below) using the blocks representations through inequalities described above.

This approach of using the k-d tree for matrix construction may lead to significant reduction of the processing time. In order to further decrease the processing time, the set of grid blocks may be subdivided into small subsets which may be processed on separate CPUs.

If p(h) is received, the block-to-block covariance R_pq for block v_p and v_q may be defined as:

1 I f f

Rpq = — Tj— I I I p (x, x )dxdx

rP^rct \ ^Vp v \\ \\^Vq i \\ ^J *vv_vp ^J ->vv_aq where r_p and r_q are the change of support coefficients for block v_p and v_q .

If, C(h) is received, the block-to-block covariance R_pq for block v_p and v_q may be defined as: φί ^rp^lrqRpq

The formula above is a polynomial equation for R_pq, which can be solved (i.e. inversed) with standard methods, e.g. Newton's method.

Then, it is possible to compute (step 1 12) unconditional simulations of the Gaussian block variables {Yy. ^cs, i = 1 ... N} given the block-to-block covariances R_pq (retrieved, for instance, from the block-to-block covariance table).

Yv_p is a Gaussian random variable, underlying the random variable Z_Vp (average value of Z(x) in the block v_p), related with Z_Vp through the function cp_Vp , which is generally different for every block: Z_Vp = <p_Vp (Y_Vp) .

Yv_p ^cs is an unconditional simulation of the Gaussian random variable Y_Vp .

For reasons of computational optimization, it is possible to simulate {Y^NCS (Vi), i = 1 ... N} {Y( >p) being the average value of the Gaussian random variable Y(x) over block v_p, Y^NCS(v_p) being then unconditional simulation of Y(v_p)) and to use relation Y^NCS{y_p) = r_pY _p ^cs to switch between these variables.

Any simulation algorithm can be used to generate {Y^A = 1 ... N}, for instance, a discrete spectral decomposition method, or a sequential Gaussian simulation method. Indeed, given a covariance matrix between components = 1 ... N,j = 1 ... N}, it is possible to generate a multivariate Gaussian vector {Y^NCS(Vi), i = 1 ... N} (with zero mean). A classical Sequential Gaussian Simulation method can be used to generate the target field. When working with unstructured grids, it is possible to address properly the neighborhood search problem using non-optimized k-d trees (Bentley, 1975) for neighborhood search in Sequential Gaussian Simulation in the following manner:

- an empty k-d tree is created ; - at each step of the Sequential Gaussian Simulation, the k-d tree is used to find a fixed number n of closest neighbors to the currently simulated node.

- after every iteration of the Sequential Gaussian Simulation, the center of the new simulated block is added to the k-d tree.

In case of working with huge reservoir grids, the k-d tree can be optimized during the Sequential Gaussian Simulation procedure at time steps when the search performance drops down. The proposed methodology may provide a significant gain in performance comparing to direct search through all of the previously simulated neighbors or using an optimized k-d tree with deactivated terminal nodes.

For each point x_a {a = 1 ... n) having a received conditioning data, it is possible to simulate (step 1 13) a random value Y_a:

or, alternatively:

where T_a are mutually independent standard Gaussian variables, v_a stands for the block containing x_a and r_a stands for the change of support coefficient of v_a. T_a has a meaning of an error when predicting Y_a given Y^NCS{v_a) or, equivivalently Yv ^s■ It can be simulated using a standard generator of pseudo-random numbers with standard Gaussian distribution.

For the purposes of simplicity of explanation, the formulas above are given for the case when only one conditioning sample is contained within a block. The generalization for the case when multiple samples are contained within a block is straightforward.

This computation is performed to prepare the data for conditioning. Generated values {Y_a} are considered to be the values at points {x_a} of the point-support phenomenon, underlying the generated unconditional block-support phenomenon

(Y^NCS(v_a) or

Indeed, it is possible to condition the simulated field with the formula (part of step 1 14):

where λ^^κ(ν_ρ) is the simple kriging coefficients for the estimate of block v_p given {y_a, = 1 ... n} with the covariance function p(h).

The coefficients { ^s _a ^K(v_p), a = 1 ... n} for every block v_p may be determined with the below-described procedure.

For a given block v_p, find a subset of closest data points {x_ai, i = 1 ... n(v_p)}, where the number n(v_p) is either an input to the algorithm or is determined by a given neighbor-search rule. Then if a <t

i = 1 ... n( _p)} then the value of (v_p) is zero. Otherwise, the other coefficients i, {λ^(ν_ρ), ί = l ... n(v_p)} are determined from the following block-kriging system:

where p(x_ai, v_p) denotes the point-to-block covariance between point x_a. and block Given linear system can be solved by any standard method, i.e. Gaussian elimination method.

Then, it is possible to determine the conditionally simulated value Y r„CS _ = ¾——) for v^rv every block v_v (part of step 1 14).

If no conditioning data are received, it is possible to simply take Y$^s = Y

Thus, it is possible to back transform the conditionally simulated value Y in Gaussian space using the block-specific transform function <p_Vp (i.e. Ζξ_ρ = <p_v (Y ) (step 1 1 5).

This algorithm enables performing direct conditional simulations of blocks (1 16). Conditioning may be done by simple or ordinary kriging in the Gaussian space (for instance thanks to the G.Matheron mehod ). The block-to-block covariance in the Gaussian space may be also preserved.

Alternatively, it is possible to back-transform (step 1 1 7) the simulated variables {Yv ^CS} and {Y_a, a = 1... n} to the real space using the corresponding transform functions <p_Vp and φ respectively. The results of this back-transformation is called

{Zv_p ^cs} and {Z_a, a = 1 ... n} in the following description. The conditioning is then performed in the real space. To condition the simulated field, the following formula may be computed (step 1 1 8):

λί^κ(ν_ρ ^~) is computed as presented above but with the conditioning data {z_a, a = 1 ... ri) and with the covariance function C(7i). The coefficients {λα^Κ(ν_ρ), a = l ... n} for every block v_p are determined with the below-described procedure.

For a given block v_p, find a subset of closest datum points {χ_α;, ί = l ... n(v_p)}, where the number n(v_p) is either an input to the algorithm or is determined by a given neighbor-search rule. Then, if a i {a_t, i = 1 ... n(v_p)} then λ^^κ(ν_ρ) is set to zero. For the other coefficients in {λ^(ν_ρ), ί = l ... n(v_p)}, they are determined from the following block-kriging system:

where C(x_a., v_p) denotes the point-to-block covariance between point x_a. and block v„.

Given above linear system can be solved by any standard method, i.e. Gaussian elimination method.

An alternative solution for conditioning in the Gaussian space exists:

where ^s _a ^K(v_p) are simple kriging coefficients for the estimating Y_v given {y_a, a = 1 ... n} that are a obtained for every block v_p with the below-described procedure. For a given block v_p, find a subset of closest data points {x_a., i = 1 ... n(v_p)}, where the number n(v_p) is either an input to the algorithm or is determined by a given neighbor-search rule. Then if a <t

i = 1 ... n(y_p)} then ^s _a ^K(v_p) is set to zero. The other coefficients in [λ^ (v_p), i = 1 ... n(y_p)} are determined from the following kriging system: cov(Y_ai, Y_a2) -

where cov (γ_α., r_a.) = r_a.r_a.R_a._a. and cov (γ_α., r_Vp) = r_a.R_a._p.

Given above linear system can be solved by any standard method, i.e. Gaussian elimination method. The result of the conditional simulation in that case is Z^cs(v_p) =

Figure 2 is a possible embodiment for a device that enables the present invention.

In this embodiment, the device 200 comprise a computer, this computer comprising a memory 205 to store program instructions loadable into a circuit and adapted to cause circuit 204 to carry out the steps of the present invention when the program instructions are run by the circuit 204.

The memory 205 may also store data and useful information for carrying the steps of the present invention as described above.

The circuit 204 may be for instance:

- a processor or a processing unit adapted to interpret instructions in a computer language, the processor or the processing unit may comprise, may be associated with or be attached to a memory comprising the instructions, or

- the association of a processor / processing unit and a memory, the processor or the processing unit adapted to interpret instructions in a computer language, the memory comprising said instructions, or

- an electronic card wherein the steps of the invention are described within silicon, or

- a programmable electronic chip such as a FPGA chip (for « Field- Programmable Gate Array >>). This computer comprises an input interface 203 for the reception of data used for the above method according to the invention and an output interface 206 for providing all information needed for simulations on the unstructured grid. To ease the interaction with the computer, a screen 201 and a keyboard 202 may be provided and connected to the computer circuit 204.

Expressions such as "comprise", "include", "incorporate", "contain", "is" and "have" are to be construed in a non-exclusive manner when interpreting the description and its associated claims, namely construed to allow for other items or components which are not explicitly defined also to be present. Reference to the singular is also to be construed in be a reference to the plural and vice versa.

A person skilled in the art will readily appreciate that various parameters disclosed in the description may be modified and that various embodiments disclosed may be combined without departing from the scope of the invention.

For instance, when a covariance is mentioned in the above description and in the claims, the person skilled in the art would understand that any other function related to the covariance may be acceptable such as a variogram.

Claims

1 . A method for mapping a geological parameter in an unstructured grid representing a geological model, wherein the method comprises:

- receiving the unstructured grid (106), the unstructured grid having a plurality of blocks v_p ;

- receiving a distribution of a random field ( (x), 101 ) of a petrophysical/geological variable;

- determining (102) a transformation function to convert the random field (Z(x)) to a Gaussian random field ( (x)), the transformation function being defined by the received distribution of the random field and by a distribution of the Gaussian random field;

- determining (103) a decomposition in an Hilbertian basis of the transformation function, {φ₀ ί = 1 ... N_∞} being the coordinates of the transformation function in said Hilbertian basis;

- receiving:

- an input covariance function (108) C(h) of the received random field; or

- an input covariance function (108) p(h) of the Gaussian random field;

- determining, for each block v_p of the unstructured grid, a change of support coefficient r_p by solving:

r_p ² = p(x, x')dxdx', or

^ Sv ^{c {x}' ^x' ^dxdx' ^{= ∑}£ ¾^{2 i}'

2. A method according to claim 1 , wherein determining the geological parameter comprises the computation (1 10) of:

No, i=0 wherein the Hilbertian basis is noted { χ ϊ = 1 ... N_∞}.

3. A method according to one of the preceding claims, wherein determining the geological parameter further comprises:

- for each couple of blocks v_p and v_q in the plurality of blocks respectively associated with r_p and r_q as change of support coefficients, computing (1 1 1 ) a block-to-block covariance R_pq of an auxiliary Gaussian random vector by solving :

R_va = -^— -—— - f f p (x, x')dxdx' ; or

4. A method according to one of the preceding claims, wherein determining the geological parameter further comprises:

- for each current block of the unstructured grid, determining (1 12) a realization of a Gaussian random variable, said Gaussian random variable being converted from a random variable representing an average value of the received random field in said current block.

5. A method according to one of the preceding claims, wherein the method further comprises: - receiving a set of conditioning data (104) for points of the geological model;

- transforming (105) the received conditioning data into Gaussian conditioning data by applying the inverse of the transformation function or the inverse of the decomposition in the Hilbertian basis of the transformation function; - for each current point of the geological model where a conditioning data is received, determining (1 13) a random value function of :

Y^NCS(v_a) + T^r_a or

wherein T_a are mutually independent standard Gaussian variables, v_a is the block containing said current point, r_a is the change of support coefficient determined for said block v_a, Y^NCS{v_a) is an average value of the Gaussian random variable Y(x) over block v_a and Y^^s is a Gaussian random variable being converted from a random variable representing an average value of the received random field in the block v_a.

6. A non-transitory computer readable storage medium, having stored thereon a computer program comprising program instructions, the computer program being loadable into a data-processing unit and adapted to cause the data-processing unit to carry out the steps of any of claims 1 to 5 when the computer program is run by the data-processing device.