🧠 👨🏽‍🎓 📴 Los geólogos tienen su propio Minecraft: cómo construir lo que no sabes a partir de lo que sabes 🔽 🤸🏽 ⛹🏿

Este es el comienzo de una historia sobre cómo las matemáticas invadieron por primera vez la geología, cómo luego vino un especialista en TI y programó todo, creando así una nueva profesión de “geólogo digital”. Esta es una historia sobre cómo el modelado estocástico se diferencia del kriging. También es un intento de mostrar cómo usted mismo puede escribir su primer software geológico y, posiblemente, transformar de alguna manera el campo de la ingeniería geológica y petrolera.

Calculemos cuánto aceite hay

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

, . , $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()