Skip to article frontmatterSkip to article content

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)

Loading...