Impact de la teneur de coupure et de la corrélation spatiale sur la localisation des ressources
Tout au long de la session, vous apprendrez à définir une teneur de coupure, à caractériser la corrélation spatiale d’un gisement, et surtout à comprendre comment ces paramètres influencent directement l’estimation des ressources.
Mais pour l’instant… pourquoi ne pas explorer un peu ?
💡 Interagissez avec ce modèle 3D pour observer comment la localisation des ressources varie selon la teneur de coupure ou le degré de corrélation spatiale. Vous verrez que de petits changements peuvent avoir un grand impact. Un excellent moyen de développer votre intuition géologique, en jouant un peu les prospecteurs numériques !
⚠️ Cette visualisation interactive est encore en développement. Certaines fonctionnalités pourraient évoluer.
Note : L’anisotropie n’est pas bien implémentée
Source
import numpy as np
from numpy.fft import fftn, ifftn
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import ipywidgets as widgets
from IPython.display import display, clear_output
# IMPORTANT : Exécutez en premier dans une cellule séparée :
# %matplotlib widget
# --- Modèle de covariance sphérique 3D vectorisé et optimisé ---
def spherical_covariance_fft_3d(fft_shape, ax, ay, az, angle_x, angle_y, angle_z, sill=1.0):
"""
Calcule le modèle de covariance sphérique anisotrope 3D,
appliquant la rotation dans l'espace fréquentiel.
"""
nx, ny, nz = fft_shape
# Coordonnées fréquentielles normalisées (centrées)
x_lags = np.fft.fftfreq(nx) * nx
y_lags = np.fft.fftfreq(ny) * ny
z_lags = np.fft.fftfreq(nz) * nz
X, Y, Z = np.meshgrid(x_lags, y_lags, z_lags, indexing='ij')
# Matrices de rotation (ZYX convention)
rx, ry, rz = np.radians([angle_x, angle_y, angle_z])
Rx = np.array([[1,0,0],
[0,np.cos(rx), -np.sin(rx)],
[0,np.sin(rx), np.cos(rx)]])
Ry = np.array([[ np.cos(ry),0,np.sin(ry)],
[0,1,0],
[-np.sin(ry),0,np.cos(ry)]])
Rz = np.array([[np.cos(rz), -np.sin(rz), 0],
[np.sin(rz), np.cos(rz), 0],
[0,0,1]])
R = Rz @ Ry @ Rx
# Applique la rotation et la mise à l'échelle anisotrope en une étape vectorisée
coords = np.stack((X.ravel(), Y.ravel(), Z.ravel()), axis=1) # (N,3)
rotated = coords @ R.T # (N,3)
# Mise à l'échelle anisotrope
rotated[:,0] /= ax
rotated[:,1] /= ay
rotated[:,2] /= az
# Calcul des distances sphériques (norme)
h = np.linalg.norm(rotated, axis=1).reshape(fft_shape)
# Modèle sphérique (vectorisé)
cov = np.zeros_like(h)
mask = h < 1
h_masked = h[mask]
cov[mask] = sill * (1 - 1.5 * h_masked + 0.5 * h_masked**3)
return cov
# --- Simulation 3D FFT-MA ---
def fftma_3d(shape, ax, ay, az, angle_x, angle_y, angle_z, sill=1.0, seed=0):
np.random.seed(seed)
fft_shape = tuple(s*2 for s in shape)
cov = spherical_covariance_fft_3d(fft_shape, ax, ay, az, angle_x, angle_y, angle_z, sill)
white_noise = np.random.normal(size=fft_shape)
cov_fft = fftn(cov)
white_fft = fftn(white_noise)
z_fft = np.sqrt(np.abs(cov_fft)) * white_fft
field = np.real(ifftn(z_fft))
# Extraction de la partie centrale
slices = tuple(slice(s//2, s//2 + orig) for s, orig in zip(fft_shape, shape))
return field[slices]
# --- Conversion gaussien vers lognormal ---
def gaussian_to_lognormal(field, mean_ln, var_ln):
if mean_ln <= 0:
print("Moyenne lognormale doit être > 0. Retourne un champ nul.")
return np.zeros_like(field)
sigma_g2 = np.log(var_ln / (mean_ln**2) + 1)
sigma_g2 = max(sigma_g2, 0) # pas négatif
sigma_g = np.sqrt(sigma_g2)
mu_g = np.log(mean_ln) - sigma_g2 / 2
return np.exp(field * sigma_g + mu_g)
# --- Tracé 3D optimisé ---
def plot_matplotlib_3d(field, cutoff):
nx, ny, nz = field.shape
mask = field >= cutoff
coords = np.argwhere(mask)
count = coords.shape[0]
total = field.size
perc = 100*count/total if total > 0 else 0
if count == 0:
print(f"Aucun voxel >= {cutoff:.2f}. Pourcentage: {perc:.2f}%")
return None
x, y, z = coords.T
values = field[mask]
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111, projection='3d')
sc = ax.scatter(x, y, z, c=values, cmap='viridis', s=10, alpha=0.8,
vmin=0, vmax=0.95*field.max())
ax.set(xlabel='X', ylabel='Y', zlabel='Z')
ax.set_title(f'Simulation 3D (Teneur ≥ {cutoff:.2f})\n{perc:.2f}% voxels')
cbar = fig.colorbar(sc, ax=ax, pad=0.1)
cbar.set_label('Teneur')
ax.set_xlim(0, nx)
ax.set_ylim(0, ny)
ax.set_zlim(0, nz)
return fig
# --- Widgets ---
common_style = {'description_width': '180px'}
text_layout = widgets.Layout(width='250px')
slider_layout = widgets.Layout(width='350px')
mean_w = widgets.FloatText(value=1.0, description='Moyenne Lognormale (ppm)', style=common_style, layout=text_layout)
variance_w = widgets.FloatText(value=2.0, description='Variance Lognormale (ppm²)', style=common_style, layout=text_layout)
cutoff_w = widgets.FloatText(value=2.0, description='Teneur coupure (ppm)', style=common_style, layout=text_layout)
ax_w = widgets.FloatText(value=20.0, description='Portée Ax (m)', style=common_style, layout=text_layout)
ay_w = widgets.FloatText(value=20.0, description='Portée Ay (m)', style=common_style, layout=text_layout)
az_w = widgets.FloatText(value=20.0, description='Portée Az (m)', style=common_style, layout=text_layout)
angle_x_w = widgets.IntSlider(value=0, min=-90, max=90, step=5, description='Rotation axe X (°)', style=common_style, layout=slider_layout)
angle_y_w = widgets.IntSlider(value=0, min=-90, max=90, step=5, description='Rotation axe Y (°)', style=common_style, layout=slider_layout)
angle_z_w = widgets.IntSlider(value=0, min=-90, max=90, step=5, description='Rotation axe Z (°)', style=common_style, layout=slider_layout)
seed_w = widgets.IntText(value=42, description='Graine aléatoire', style=common_style, layout=text_layout)
button = widgets.Button(description='Exécuter la simulation', button_style='primary', icon='play', layout=widgets.Layout(width='90%'))
output = widgets.Output()
def on_button_clicked(b):
with output:
clear_output(wait=True)
print("Génération en cours... (cela peut prendre quelques secondes pour de grandes tailles)")
shape = (100, 100, 50)
try:
gaussian_field = fftma_3d(shape, ax_w.value, ay_w.value, az_w.value,
angle_x_w.value, angle_y_w.value, angle_z_w.value,
sill=1.0,
seed=seed_w.value)
lognormal_field = gaussian_to_lognormal(gaussian_field, mean_w.value, variance_w.value)
fig = plot_matplotlib_3d(lognormal_field, cutoff_w.value)
if fig:
display(fig)
except Exception as e:
print(f"Erreur lors de la simulation ou du tracé : {e}\nVérifiez les paramètres et la taille.")
button.on_click(on_button_clicked)
ui = widgets.VBox([
widgets.Label("Paramètres de Distribution et Continuité Spatiale:"),
mean_w,
variance_w,
ax_w,
ay_w,
az_w,
widgets.Label("Paramètres Économique:"),
cutoff_w,
widgets.Label("Paramètres de Rotation:"),
angle_x_w,
angle_y_w,
angle_z_w,
widgets.Label("Paramètre de Simulation:"),
seed_w,
button,
output
])
display(ui)