💍 🤳🏻 🌇 Geologists have their own minecraft: how to build what you don't know from what you know ⬅️ 🔚 🧓🏻

This is the beginning of a story about how mathematics first invaded geology, how then an IT specialist came and programmed everything, thereby creating a new profession of "digital geologist". This is a story about how stochastic modeling differs from kriging. It is also an attempt to show how you yourself can write your first geological software and, possibly, somehow transform the industry of geological and petroleum engineering.

Let's calculate how much oil is there

, , , , . , , . , (, ). , , , .

, . , $h$ , $\phi$ $s$ , «» $S$ , $h\cdot\phi\cdot s\cdot S$ ( -) . ( , )

— , , , , , , , .

$V$ , $u$ — $V$ . $\phi\left(u\right)$ $u$ , $s\left(u\right)$ — . $\intop_{V}s\left(u\right)\cdot\phi\left(u\right)du$ — . $\phi\left(u\right)$ , $s\left(u\right)$ ( , , ) $V$ . $V$ . , , ( , , ).

. , , , , . — , , . — , , , . , , .

. ; , , , .

, , . , , , , . , .

, ( , ). , . , , .

, , , . , , . , . . , , , - , .

, — , .

, , .. , , . , , , . — — «support». , , , , . .

, , (Krige, D.G. 1951. A statistical approach to some basic mine valuation problems on the Witwatersrand. Journal of the Chemical, Metallurgical and Mining Society of South Africa, December 1951. pp. 119–139.). , -, , , , . , , 5×5 , , , 1 , 50 ×50 ×1 — , , ( — upscaling).

. $z\left(u\right)$ , $u$ — . $z_{i}=z\left(u_{i}\right)$ .

\bar{z} (u) = F (u, z_{1}, . ., z_{n}, u_{1}, . . ., u_{n})

$\bar{z}\left(u\right)=F\left(u,z_{1},..,z_{n},u_{1},...,u_{n}\right)$

. ,

z_{k} \equiv F (u_{k}, z_{1}, . ., z_{n}, u_{1}, . . ., u_{n}), (1)

$z_{k}\equiv F\left(u_{k},z_{1},..,z_{n},u_{1},...,u_{n}\right), (1)$

\bar{z} = \sum_{i} z_{i} ν_{i},

$\bar{z}=\sum_{i}z_{i}\nu_{i},$

$\nu_{i}$ — , $u$ $u_{i}$ . , . , , , .

\bar{z} (u) = \sum_{i} z_{i} \frac{1}{{| u - u_{i} |}^{m}} \cdot {(\sum_{i} \frac{1}{{| u - u_{i} |}^{m}})}^{- 1}, (2)

$\bar{z}\left(u\right)=\sum_{i}z_{i}\frac{1}{\left|u-u_{i}\right|^{m}}\cdot\left(\sum_{i}\frac{1}{\left|u-u_{i}\right|^{m}}\right)^{-1},(2)$

$m$ — . $m=1$ (. . 1). $m=2$ , . , . $u=u_{k}$ , . , $u$ $u_{k}$ , $z_{k}$ 1, $z$ , $\bar{z}\left(u_{k}\right)$ $z_{k}$ .

Python .

Inverse distance interpolation

import numpy as np
import matplotlib.pyplot as pl

# num of data
N = 5

np.random.seed(0)

# source data
z = np.random.rand(N)
u = np.random.rand(N)

x = np.linspace(0, 1, 100)
y = np.zeros_like(x)

# norm weights
w = np.zeros_like(x)

# power
m = 2

# interpolation
for i in range(N):
    y += z[i] * 1 / np.abs(u[i] - x) ** m
    w += 1 / np.abs(u[i] - x) ** m

# normalization
y /= w

# add source data
x = np.concatenate((x, u))
y = np.concatenate((y, z))
order = x.argsort()
x = x[order]
y = y[order]

# draw graph
pl.figure()

pl.scatter(u, z)
pl.plot(x, y)

pl.show()

pl.close()

. , , - . , , . -, , , , ? , . , 1 . , , . , , (. . 2). .

, , , , . , , . , . — , , : , , , . , .

(2) :

\bar{z} (u) = \sum_{i} z_{i} \cdot f_{i} (u), (3)

$\bar{z}(u)=\sum_{i} z_i\cdot f_i(u), (3)$

, , . , , $i$ - $k$ - . .

, :

\bar{z} (u) = \sum_{i} λ_{i} \cdot c (∥ u - u_{i} ∥), (4)

$\bar{z}(u)=\sum_{i}\lambda_i\cdot c\left(∥u-u_i∥\right), (4)$

$c(h)$ - , $\lambda_i$ — . (4) — «», . $\lambda_i$ . , (1):

z_{k} = \sum_{i} λ_{i} \cdot c (∥ u_{k} - u_{i} ∥), (5)

$z_k=\sum_{i}\lambda_i\cdot c\left(∥u_k-u_i∥\right), (5)$

$\lambda$ . $Z=C\cdot\lambda$ , $Z$ — $z_k$ , $\lambda$ — , $C$ — «» $c(h)$ . , $\lambda=C^{-1}\cdot Z$ , (6),

\bar{z} (u) = \sum_{i} z_{i} \cdot g_{i} (u), (6)

$\bar{z}(u)=\sum_{i}z_i\cdot g_i\left(u\right), (6)$

g_{i} (u) = \sum_{k} C_{i, k}^{- 1} \cdot c (∥ u_{k} - u ∥) . (7)

$g_i(u)=\sum_{k}C^{-1}_{i, k}\cdot c\left(∥u_k-u∥\right). (7)$

(6) ( ). (4) , , dual kriging. (6) (3), , $g_i$ . $u$ $u_k$ , (7) , , $C$ , , $g_i(u_i)=1$ $g_i(u_k)=0$ $i≠k$ ( $f_i$ ).

(. . 3). , , .

, — Python:

Kriging interpolation

import numpy as np
import matplotlib.pyplot as pl

# num of data
N = 5

np.random.seed(0)

# source data
z = np.random.rand(N) - 0.5
u = np.random.rand(N)

x = np.linspace(0, 1, 100)
y = np.zeros_like(x)

# covariance function
def c(h):
    return np.exp(-np.abs(h ** 2 * 20.))

# covariance matrix
C = np.zeros((N, N))

for i in range(N):
    C[i, :] = c(u - u[i])

# dual kriging weights
lamda = np.linalg.solve(C, z)

# interpolation
for i in range(N):
    y += lamda[i] * c(u[i] - x)

# add source data
x = np.concatenate((x, u))
y = np.concatenate((y, z))
order = x.argsort()
x = x[order]
y = y[order]

# draw graph
pl.figure()

pl.scatter(u, z)
pl.plot(x, y)

pl.show()

pl.close()

, , . (5) (6) , .

, , . , , , .

. , (), , , , . , «», . , (6), , «», . , «» . : . , , , — . , , .

. , . . , , .

Libraries

from theory import probability
from numpy import linalg

, , .

, . () : u $z(u)$ , ( ). , $z_i=z(u_i)$ .

. . , ( , ).

, - — . - , . : $z(u)$ — ( ), $z(u_i)$ , $u_i$ , , . $z(u)$ — , .

( $Cov\left(z_1, z_2\right)=E\left(\left(z_1-E\left(z_1\right)\right)⋅\left(z_2-E\left(z_2\right)\right)\right)$ , ) . — . , u $z(u)$ , .

, , : $Cov\left(z\left(u\right),z\left(v\right)\right)=c\left(∥u-v∥\right)$ . -, $c(h)$ — «» (4), : $c(h)$ $h$ . , $c(500)/c(0) = 0.5$ , $500$ $50$ .., $c(1000) = 0$ , , .

, , , , . , , . , , . $Z$ $\zeta$ .

Z = A \cdot ζ .

$Z=A⋅ζ.$

$A$ , ( $Z$ ). :

C = E (Z \cdot Z^{T}) = E (A \cdot ζ \cdot ζ^{T} \cdot A^{T}) = A \cdot E (ζ \cdot ζ^{T}) \cdot A^{T} = A \cdot A^{T},

$C=E\left(Z\cdot Z^T\right)=E\left(A\cdot\zeta\cdot\zeta^T\cdot A^T\right)=A\cdot E\left(\zeta\cdot\zeta^T\right)\cdot A^T=A\cdot A^T,$

, $\zeta$ , , , $E\left(\zeta\cdot\zeta^T\right)$ . . $Z$ , , $C$ , $C=A\cdot A^T$ ( ) $Z$ $\zeta$ , , . . .

, ( , , — , ). , (. . 4).

Unconditional stochastic gaussian modeling

import numpy as np
import matplotlib.pyplot as pl

np.random.seed(0)

# source data
N = 100
x = np.linspace(0, 1, 100)

# covariance function
def c(h):
    return np.exp(-np.abs(h ** 2 * 250))

# covariance matrix
C = np.zeros((N, N))

for i in range(N):
    C[i, :] = c(x - x[i])

# eigen decomposition
w, v = np.linalg.eig(C)

A = v @ np.diag(w ** 0.5)

# you can check, that C == A @ A.T

# independent normal values
zeta = np.random.randn(N)

# dependent multinormal values
Z = A @ zeta

# draw graph
pl.figure()

pl.plot(x, Z)

pl.show()

pl.close()

, , ( , ). , 5.

Conditional stochastic gaussian simulation

import numpy as np
import matplotlib.pyplot as pl

np.random.seed(3)

# source data
M = 5

# coordinates of source data
u = np.random.rand(M)

# source data
z = np.random.randn(M)

# Modeling mesh
N = 100
x = np.linspace(0, 1, N)

# covariance function
def c(h):
return np.exp(−np.abs(h ∗∗ 2 ∗ 250))

# covariance matrix mesh−mesh
Cyy = np.zeros ((N, N))
for i in range (N):
    Cyy[ i , : ] = c(x − x[i])

# covariance matrix mesh−data
Cyz = np.zeros ((N, M))

# covariance matrix data−data
Czz = np.zeros ((M, M))
for j in range (M):
    Cyz [:, j] = c(x − u[j])
    Czz [:, j] = c(u − u[j])

# posterior covariance
Cpost = Cyy − Cyz @ np.linalg.inv(Czz) @ Cyz.T

# lets find the posterior mean, i.e. Kriging interpolation
lamda = np.linalg.solve (Czz, z)
y = np.zeros_like(x)

# interpolation
for i in range (M):
    y += lamda[i] ∗ c(u[i] − x)

# eigen decomposition
w, v = np.linalg.eig(Cpost)
A = v @ np.diag (w ∗∗ 0.5)

# you can check, that Cpost == A@A.T

# draw graph
pl.figure()
for k in range (5):
    # independent normal values
    zeta = np.random.randn(N)
    # dependent multinormal values
    Z = A @ zeta
    pl.plot(x, Z + y, color=[(5 − k) / 5] ∗ 3)
    pl.plot(x, Z + y, color=[(5 − k) / 5] ∗ 3, label=’Stochastic realizations’)

pl.plot(x, y, ’. ’, color=’blue’, alpha=0.4, label=’Expectation(Kriging)’)
pl.scatter(u, z, color=’red ’, label=’Source data’)
pl.legend()
pl.show()
pl.close()