🎶 😄 🙏🏿 地質学者は独自のマインクラフトを持っています：あなたが知っていることからあなたが知らないことを構築する方法 🌼 🍲 👩🏾‍🤝‍👩🏻

これは、数学が最初に地質学に侵入した方法、その後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()