👩🏽‍🔧 🛀🏿 🎆 地质学家拥有他们自己的Minecraft：如何根据所学知识构建您所不知道的 🖨️ 👴🏽 🧖🏿

这是关于数学如何首先入侵地质学，然后IT专家如何来对所有事物进行编程并由此创建“数字地质学家”这一新职业的故事的开篇。这是一个关于随机建模与克里金法有何不同的故事。这也是尝试展示您自己如何编写第一个地质软件，并可能以某种方式改变地质和石油工程领域的尝试。

让我们计算一下那里有多少油

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

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