US 2021/0263176 A1 · MULTI-POINT GEOSTATISTICAL PRESTACK INVERSION METHOD BASED ON RENEWAL PROBABILITY RATIO CONSTANT THEORY

Abstract

A multi-point geostatistical prestack inversion method based on a renewal probability ratio constant theory, comprising the following steps: sorting information, work area gridding and distributing well data, assigning an initial attribute value to a simulation work area, selecting a data template with an appropriate size, inversion and judging iteration termination; the method overcomes the defects that in the prior art, oil and gas exploration and development are more and more difficult, exploration degree is higher and higher, new oil and gas reservoirs are more and more difficult to find, reservoir parameters are inaccurate, and exploration uncertainty is increased, and the method has the advantages that prior information is obtained through a multi-point geostatistics method and then screening is carried out through a minimum objective function such that the complexity of seismic inversion is reduced.

A multi-point geostatistical prestack inversion method based on a renewal probability ratio constant theory, comprising following specific steps:

{circle around (1)} sorting out information including:
checking whether original seismic information and well information is complete, with the original seismic information being prestack seismic data, and the well information comprising drilling core and logging data;
according to well information interpretation result, it being possible that lithofacies corresponding to different well depths and elastic parameters of different lithofacies can be determined, with the elastic parameters comprising density, longitudinal wave velocity, and shear wave velocity;
establishing facies proportions of different lithofacies according to lithofacies data, and establishing cumulative distribution functions of elastic parameters of different lithofacies according to elastic parameter data;
and
establishing a training image conforming to reservoir characteristics of a work area;
{circle around (2)} performing work area gridding and distributing well data comprising:
selecting an appropriate grid size according to an actual work area range;
performing grid division on the work and area;
establishing a grid model;
gridding core data and logging data according to a plane position and well depth of each well;
with the core data being lithofacies, interpreting the logging data including density, longitudinal wave velocity, and shear wave velocity;
and
distributing the core data and the logging data as hard data to a divided nearest neighbor grid node;
wherein:
information contained in actual drilling within a work area being called hard data or condition data, with the hard data comprising lithofacies, density, longitudinal wave data, and shear wave velocity;
a grid distributed with the hard data of the lithofacies being called a known grid or a known point;
a grid without lithofacies data being called to-be-evaluated grid or also called to-be-evaluated point;
after performing work area gridding and distributing well data, the distributing well data being drilling data and logging data;
the known grid being a grid distributed with lithofacies data, the known grid being a grid through which drilling passes, and the to-be-evaluated grid being a grid through which no drilling passes;
the known grid of the work area containing lithofacies, density, longitudinal wave velocity, and shear wave velocity, and the to-be-evaluated grid not having any data;
{circle around (3)} assigning an initial attribute value to a simulation work area
according to a statistical cumulative distribution function of a rock elastic parameter of the work area, an initial attribute value being assigned to a to-be-evaluated point of the work area, wherein the initial attribute value refers to elastic parameter;
{circle around (4)} selecting a data template with an appropriate size comprising:
determining a shape and size of a data template according to morphological characteristics of sedimentary facies;
it being possible that the shape of the data template can be determined according to heterogeneity of the sedimentary facies, it being possible that a two-dimensional ellipse or a three-dimensional ellipsoid can be adopted when the heterogeneity is strong, and it being possible that a two-dimensional rectangle or a three-dimensional cuboid can be adopted when the heterogeneity is weak;
the size of the data template adopting a size of 5×7 or 5×5×7;
{circle around (5)} obtaining inversion
wherein the inversion is obtained by following three sub-steps:
sub-step 5.1:
ordering access of simulation nodes;
wherein:
firstly, pseudo-random numbers are randomly generated for all grids in the work area with numerical values of the pseudo-random numbers being between 0 and 1, and then the size of a data template of the to-be-evaluated grid is taken to search a number of condition points around;
the numerical value of each to-be-evaluated point is equal to a sum of the pseudo-random number and the number of condition points;
all to-be-evaluated points are sorted from large to small;
and
a grid with a large value is simulated in priority to obtain a simulation path of an entire work area;
sub-step 5.2:
obtaining suggested data pattern;
if primary inversion is carried out and the to-be-evaluated grid is processed,
firstly the number of condition points, a position relative to the to-be-evaluated point and a lithofacies type in a data template range with the to-be-evaluated point as a center are determined to form a data event presented in a vector form with the to-be-evaluated point as the center;
the data event is utilized to randomly scan a training image established in step {circle around (1)} to obtain a first completely matched data pattern therefrom as a suggested data pattern;
and
the suggested data pattern represents a sedimentary mode, which is a lithofacies data structure with the size of the data template and where the lithofacies data structure of the condition points in the data template can be fully embodied in the suggested data pattern;
if iterative inversion is carried out and the to-be-evaluated grid is processed,
firstly, the data event is formed according to the to-be-evaluated point and a condition point of a data template range of the to-be-evaluated point, and the training image established in step {circle around (1)} is scanned randomly through the data event to obtain a plurality of completely matched data patterns from training images as candidate data patterns;
the candidate data patterns are sorted according to the lithofacies type of a central point to make cumulative probability distribution;
the candidate data pattern can contain up to a first 50 data patterns scanned from first-to-last training images, and when less than 50 candidate data patterns, all the data patterns of the training images scanned with the data event are used as candidate data patterns, and proportions of different lithofacies are calculated according to differences of central lithofacies of the candidate data patterns and denoted as ax1:
ax2:
. . . :
axm and a probability of elastic parameters of different facies is calculated according to a last elastic parameter according to following formula:
b
=
1
(
8
⁢
π
3
⁢
σ
ρ
i
·
σ
Vp
i
·
σ
V
Si
)
n
2
·
exp
⁡
[
∑
j
=
1
n
⁢
(
(
ρ
i
-
μ
ρ
)
2
2
⁢
σ
ρ
i
2
+
(
V
p
i
-
μ
V
⁢
p
)
2
2
⁢
σ
V
⁢
⁢
ρ
i
2
+
(
V
S
i
-
μ
V
⁢
s
)
2
2
⁢
σ
Vs
i
2
)
]
wherein in the formula,
b represents the probability of the lithofacies elastic parameter,
represents a number of grids in a range of the data template of the to-be-evaluated point,
ρi represents a density value of an ith grid in the range of data template surrounding the to-be-evaluated point,
Vpi represents the longitudinal wave velocity of the ith grid in the range of the data template surrounding the to-be-evaluated point,
VSi represents the shear wave velocity of the ith grid in the range of the data template surrounding the to-be-evaluated point,
μρ represents a mean value of the lithofacies density,
μVp represents the mean value of a lithofacies longitudinal wave velocity,
μVs represents the mean value of the lithofacies shear wave velocity,
σρi represents a variance of the lithofacies density,
σVpi represents the variance of the lithofacies longitudinal wave velocity, and
σVSi represents the variance of the lithofacies shear wave velocity;
a proportion of elastic parameters of different lithofacies is denoted as:
bx1:
bx2:
. . . :
bxx;
finally, a joint probability of different facies is calculated with (ax1·bx1):(ax2·bx2) :
. . . :
(axm·bxm), and denoted as Px1:
Px2:
Pxm;
the joint probability is normalized to obtain P′x1:
P′x2:
. . . :
P′xm, and the joint probability of different lithofacies is made into a cumulative distribution function;
then the candidate data patterns are sorted according to the lithofacies of a central grid, followed by randomly sampling to obtain a random number between 0 and 1, and the suggested data pattern is determined according to a position of a cumulative distribution function formed by the random number in the joint probability and the cumulative probability distribution of candidate modes;
sub-step 5.3:
obtaining suggested elastic parameter;
wherein according to the cumulative distribution function of the elastic parameters of different rock facies counted in step {circle around (1)}, the elastic parameters of the rock facies at each position are sampled in the suggested data pattern;
the elastic parameter comprises a density value, a longitudinal wave velocity value and a shear wave velocity value;
sub-step 5.4:
calculating a synthetic record of a suggested elastic parameter in the suggested data pattern;
pattern is calculated;
wherein:
a reflection coefficient is calculated based on following Connolly formula by using the suggested elastic parameter and a preset incident angle;
EI(m, θ)=Vp(1+tan2 θ)·VS(−8K sin2 θ)·ρ(1-4K sin2 θ)
in the Connolly formula,
K
=
V
s
2
V
p
2
,
EI (m, 0) represents the reflection coefficient when the preset incident angle is θ, θ represents the preset incident angle, m represents the elastic parameters, comprising Vp, Vs and ρ, Vp represents the longitudinal wave velocity, VS represents the shear wave velocity, and ρ represents density;
based on the reflection coefficient and seismic ground wave, initial forward modeling is calculated by adopting following convolution formula:
g(m)=w(θ)*EI(m,θ)
in the convolution formula, g(m) represents initial forward modeling, m represents the elastic parameter, and w(θ) represents a seismic wavelet when the preset incident angle is θ;
and
sub-step 5.5:
selecting the elastic parameter with a highest matching rate with an original prestack data;
wherein:
the synthetic record of a last elastic parameter and the suggested elastic parameter are calculated respectively by sub-step 5.3 and sub-step 5.4, and the last elastic parameter is denoted as g(ml), and the suggested elastic parameter is denoted as g(ms);
the original prestack elastic parameter is denoted as go;
herein the last elastic parameter can be an initial elastic parameter or a suggested elastic parameter in a previous iteration process;
a plurality of synthetic seismic traces is compared with an actual seismic record by adopting following seismic comparison formula, and a stratum model parameter with a minimum objective function and corresponding suggested data pattern are selected as an optimal inversion result:
m
=
{
m
l
,

g
⁡
(
m
l
)
-
g
o

<

g
⁡
(
m
s
)
-
g
o

⁢
m
s
,

g
⁡
(
m
l
)
-
g
o

≥

g
⁡
(
m
s
)
-
g
o

in the seismic comparison formula, m represents a selected elastic parameter, m1 represents the last elastic parameter, and ms represents the suggested elastic parameter in the suggested data model;
and
{circle around (6)} outputting a simulation result after iteration is judged to be terminated
herein sub-steps 5.1 to 5.5 being repeated until all grids are traversed, that is, simulation at all to-be-evaluated points being completed is that one-time inversion realization is achieved;
wherein:
when iterative inversion times are less than 7, step {circle around (5)} being repeated and next iterative inversion being carried out;
and
when the iterative inversion is completed at a seventh time, a whole inversion process being completed, a multi-point geostatistical seismic inversion model being obtained, and a simulation result being output.

260 more elements

Multi-point geostatistical prestack inversion method based on renewal probability ratio constant theory

Life of this patent

Multi-point geostatistical prestack inversion method based on renewal probability ratio constant theory

Life of this patent

Abstract

Claims8