🦊 🙅🏻 🏥 Geologen haben ihre eigene Minecraft: Wie man aus dem, was man weiß, baut, was man nicht weiß 🙍 🚴🏼 🥠

Dies ist der Beginn einer Geschichte darüber, wie Mathematik zuerst in die Geologie eindrang, wie dann ein IT-Spezialist kam und alles programmierte, wodurch ein neuer Beruf des „digitalen Geologen“ geschaffen wurde. Dies ist eine Geschichte darüber, wie sich die stochastische Modellierung vom Kriging unterscheidet. Es ist auch ein Versuch zu zeigen, wie Sie selbst Ihre erste geologische Software schreiben und möglicherweise die Industrie der Geologie- und Erdöltechnik irgendwie verändern können.

Berechnen wir, wie viel Öl dort ist

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

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