Atelier 1 - Calculateur - Note de cours et ateliers interactifs

Objectif¶

Ce calculateur vise à résoudre le krigeage ordinaire et simple avec un modèle de variogramme imbriqué anisotrope.

Exemple d’application du calculateur avec 4 observations et deux composantes dont une anisotrope.¶

Un gisement d’or est caractérisé par un variogramme sphérique présentant une anisotropie géométrique. La direction de meilleure continuité spatiale (azimut) est de 40° (sens horaire depuis le Nord). Le modèle variogrammique imbriqué comprend :

Une pépite (effet microstructurel) de variance (C_0 = 2) ppm²,
Un composant structuré sphérique de variance (C = 30) ppm²,
Une portée majeure (a_{40} = 60) m,
Une portée mineure (a_{130} = 40) m.

Données d’entrée

Les données échantillons autour du point d’intérêt ((100, 100)) sont :

Observation	Coordonnée X (m)	Coordonnée Y (m)	Teneur (ppm)
1	100	90	10
2	91	104	4
3	110	110	7
4	107	93	3

Instructions¶

Configurer les points connus : Activez/désactivez les observations et modifiez leurs coordonnées ou teneurs si nécessaire. Maximum de 7 observations.
Paramètres des variogrammes : Les paramètres du variogramme sont préconfigurés selon l’exemple ci-dessus, avec possibilité de les modifier.
Point de prédiction : Positionnez le point où vous souhaitez estimer la teneur (par défaut (x=100, y=100)).
Effectuer le calcul : Cliquez sur le bouton “Faire le calcul” pour obtenir :
- La matrice de covariance (K),
- Le vecteur de covariance (k),
- La matrice augmentée (K_{aug}) utilisée pour le krigeage ordinaire,
- Les poids de krigeage pour la méthode simple et ordinaire,
- Le multiplicateur de Lagrange (ordinaire),
- Les estimations et variances aux points de prédiction.

Note sur l’anisotropie¶

L’anisotropie géométrique est prise en compte en appliquant une transformation des coordonnées selon l’azimut et les portées majeures et mineures. L’angle est mesuré depuis le Nord, dans le sens horaire, ce qui correspond à la convention géologique classique.

Avertissements¶

Si la portée mineure est configurée plus grande que la portée majeure pour un composant, un message d’avertissement s’affichera.

Profitez de cet outil pour valider vos calculs des exercices.

Bon travail !

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import cdist
import ipywidgets as widgets
from IPython.display import display, HTML, clear_output

# --- Covariance ---
def spherical_cov(h, range_, sill):
    cov = np.zeros_like(h)
    hr = h / range_
    mask = h <= range_
    cov[mask] = sill * (1 - 1.5 * hr[mask] + 0.5 * hr[mask]**3)
    return cov

def exponential_cov(h, range_, sill): return sill * np.exp(-3*h / range_)
def gaussian_cov(h, range_, sill): return sill * np.exp(-3*(h / range_)**2)
def nugget_cov(h, sill): return sill * (h == 0)

# --- Anisotropie géométrique ---
def apply_anisotropy(coords, azimuth_deg, range_ag, range_ap):
    theta = np.deg2rad(-(90 - azimuth_deg))
    R = np.array([
        [np.cos(theta), -np.sin(theta)],
        [np.sin(theta),  np.cos(theta)]
    ])
    S = np.array([
        [1/range_ag, 0],
        [0, 1/range_ap]
    ])
    A = S @ R
    return (A @ coords.T).T

# --- Modèle imbriqué ---
def nested_covariance(h, params):
    cov_total = np.zeros_like(h)
    for p in params:
        if p['type'] == 'Sphérique': cov_total += spherical_cov(h, 1, p['sill'])
        elif p['type'] == 'Exponentiel': cov_total += exponential_cov(h, 1, p['sill'])
        elif p['type'] == 'Gaussien': cov_total += gaussian_cov(h, 1, p['sill'])
        elif p['type'] == 'Pépite': cov_total += nugget_cov(h, p['sill'])
    return cov_total

def compute_cov_matrix_anisotrope(coords, params, azimuths, ranges_ag, ranges_ap):
    n = len(coords)
    cov_total = np.zeros((n, n))
    for i, p in enumerate(params):
        coords_trans = apply_anisotropy(coords, azimuths[i], ranges_ag[i], ranges_ap[i])
        h = cdist(coords_trans, coords_trans)
        cov_total += nested_covariance(h, [p])
    return cov_total

def compute_cov_vector_anisotrope(coords, coords_pred, params, azimuths, ranges_ag, ranges_ap):
    cov_vec = np.zeros((len(coords_pred), len(coords)))
    for i, p in enumerate(params):
        coords_trans = apply_anisotropy(coords, azimuths[i], ranges_ag[i], ranges_ap[i])
        coords_pred_trans = apply_anisotropy(coords_pred, azimuths[i], ranges_ag[i], ranges_ap[i])
        h = cdist(coords_pred_trans, coords_trans)
        cov_vec += nested_covariance(h, [p])
    return cov_vec

def simple_kriging(xy_samples, values, xy_pred, params, azimuths, ranges_ag, ranges_ap, mean):
    C = compute_cov_matrix_anisotrope(xy_samples, params, azimuths, ranges_ag, ranges_ap)
    c0 = compute_cov_vector_anisotrope(xy_samples, xy_pred, params, azimuths, ranges_ag, ranges_ap)
    weights = np.linalg.solve(C, c0.T)
    y_est = mean + weights.T @ (np.array(values) - mean)
    sill_tot = sum(p['sill'] for p in params)
    sigma2 = sill_tot - np.sum(weights * c0.T, axis=0)
    return y_est, sigma2, weights, c0, C

def ordinary_kriging(xy_samples, values, xy_pred, params, azimuths, ranges_ag, ranges_ap):
    n = len(xy_samples)
    C = compute_cov_matrix_anisotrope(xy_samples, params, azimuths, ranges_ag, ranges_ap)
    C_aug = np.zeros((n+1,n+1))
    C_aug[:n,:n] = C
    C_aug[n,:n] = 1
    C_aug[:n,n] = 1
    c0 = compute_cov_vector_anisotrope(xy_samples, xy_pred, params, azimuths, ranges_ag, ranges_ap)
    y_est, sigma2, weights_list, lambdas = [], [], [], []
    sill_tot = sum(p['sill'] for p in params)
    c_aug_list = []
    for i in range(len(xy_pred)):
        c_aug = np.append(c0[i],1)
        c_aug_list.append(c_aug)   # sauvegarde du vecteur k_aug
        sol = np.linalg.solve(C_aug, c_aug)
        w, mu = sol[:-1], sol[-1]
        y_est.append(np.dot(w, values))
        sigma2.append(sill_tot - np.dot(w,c0[i]) - mu)
        weights_list.append(w)
        lambdas.append(mu)
    return np.array(y_est), np.array(sigma2), np.array(weights_list), np.array(lambdas), C_aug, np.array(c_aug_list)

def format_matrix_html(mat, title):
    html = f"<h4>{title}</h4><table border='1' style='border-collapse:collapse'>"
    for row in mat:
        html += "<tr>" + "".join(f"<td style='padding:5px;text-align:center'>{v:.3f}</td>" for v in row) + "</tr>"
    html += "</table><br>"
    return html

# --- Widgets ---

data_points = [
    (True, 100,  90, 10),
    (True,  91, 104,  4),
    (True, 110, 110,  7),
    (True, 107,  93,  3),
    (False, 0,    0,  0),
    (False, 0,    0,  0),
    (False, 0,    0,  0)
]

activate_points = []
pts_x = []
pts_y = []
vals = []

for i in range(7):
    active, x_val, y_val, v_val = data_points[i]
    activate_points.append(widgets.Checkbox(value=active, description=f'Point {i+1}', layout=widgets.Layout(width='120px')))
    pts_x.append(widgets.FloatText(value=x_val, description=f'x{i+1}', layout=widgets.Layout(width='180px')))
    pts_y.append(widgets.FloatText(value=y_val, description=f'y{i+1}', layout=widgets.Layout(width='180px')))
    vals.append(widgets.FloatText(value=v_val, description=f'val{i+1}', layout=widgets.Layout(width='180px')))

x0 = widgets.FloatText(value=100, description='x0', layout=widgets.Layout(width='180px'))
y0 = widgets.FloatText(value=100, description='y0', layout=widgets.Layout(width='180px'))
mean_ks = widgets.FloatText(value=0, description='Moyenne KS', layout=widgets.Layout(width='180px'))

options_cov = ['Aucun','Sphérique','Exponentiel','Gaussien','Pépite']
def create_cov_widgets(n):
    if n == 1:
        sill = 2.0
        range_ap = 1.0
        range_ag = 1.0
        angle = 40
        typ = 'Pépite'
    elif n == 2:
        sill = 30.0
        range_ap = 40.0
        range_ag = 60.0
        angle = 40
        typ = 'Sphérique'
    else:
        sill = 0.0
        range_ap = 1.0
        range_ag = 1.0
        angle = 0
        typ = 'Aucun'

    sill_w      = widgets.FloatText(value=sill, description=f'Sill {n}', layout=widgets.Layout(width='180px'))
    range_ap_w  = widgets.FloatText(value=range_ap, description=f'Portée ap {n}', layout=widgets.Layout(width='180px'))
    range_ag_w  = widgets.FloatText(value=range_ag, description=f'Portée ag {n}', layout=widgets.Layout(width='180px'))
    angle_w     = widgets.IntText(value=angle, description=f'Angle {n} (°)', layout=widgets.Layout(width='180px'))
    typ_w       = widgets.Dropdown(options=options_cov, value=typ, description=f'Type {n}', layout=widgets.Layout(width='180px'))
    return sill_w, range_ap_w, range_ag_w, angle_w, typ_w

cov_widgets = [create_cov_widgets(i) for i in range(1,4)]

output = widgets.Output()

def update(*args):
    with output:
        clear_output()
        xs, ys, vs = [], [], []
        warnings = []
        for i in range(7):
            if activate_points[i].value:
                try:
                    xval = float(pts_x[i].value)
                    yval = float(pts_y[i].value)
                    vval = float(vals[i].value)
                    xs.append(xval)
                    ys.append(yval)
                    vs.append(vval)
                except:
                    pass

        if len(xs) < 1:
            print("Au moins un point doit être activé et valide.")
            return

        xy_samples = np.column_stack([xs, ys])
        xy_pred = np.array([[float(x0.value), float(y0.value)]])

        params, azimuths, ranges_ag, ranges_ap = [], [], [], []
        for idx, (sill_w, range_ap_w, range_ag_w, angle_w, type_w) in enumerate(cov_widgets, start=1):
            sill = sill_w.value
            range_ap = range_ap_w.value
            range_ag = range_ag_w.value
            angle = angle_w.value
            typ = type_w.value
            if typ != 'Aucun' and sill > 0:
                if range_ap > range_ag:
                    print(f"⚠️ Composant {idx} a portée petite ({range_ap}) > portée grande ({range_ag}).")
                params.append({'sill': sill, 'type': typ})
                azimuths.append(angle)
                ranges_ag.append(range_ag)
                ranges_ap.append(range_ap)

        if len(params) == 0:
            print("Choisir au moins un composant variogramme valide.")
            return

        y_ks, var_ks, w_ks, c0_ks, C_ks = simple_kriging(xy_samples, vs, xy_pred, params, azimuths, ranges_ag, ranges_ap, mean_ks.value)
        y_ko, var_ko, w_ko, lambda_ko, C_aug, c_aug_ko = ordinary_kriging(xy_samples, vs, xy_pred, params, azimuths, ranges_ag, ranges_ap)

        display(HTML(format_matrix_html(C_ks, "Matrice de covariance $K_{KS}$")))
        display(HTML(format_matrix_html(c0_ks, "Vecteur covariance $k_{KS}$")))
        display(HTML(format_matrix_html(C_aug, "Matrice augmentée $K_{KO}$")))
        display(HTML(format_matrix_html(c_aug_ko, "Vecteur covariance $k_{KO}$")))

        display(HTML(format_matrix_html(np.array([w_ks[:,0]]), "Poids KS")))
        display(HTML(format_matrix_html(np.array([w_ko[0]]), "Poids KO")))
        display(HTML(format_matrix_html(np.array([[lambda_ko[0]]]), "Multiplicateur de Lagrange (KO)")))

        print(f"Estimation KS: {y_ks[0]:.3f} ; Variance KS: {var_ks[0]:.3f}")
        print(f"Estimation KO: {y_ko[0]:.3f} ; Variance KO: {var_ko[0]:.3f}")

box_pts = widgets.HBox([
    widgets.VBox(activate_points),
    widgets.VBox(pts_x),
    widgets.VBox(pts_y),
    widgets.VBox(vals)
])

box_cov = widgets.HBox([widgets.VBox(cw) for cw in cov_widgets])
box_krigeage = widgets.HBox([mean_ks, x0, y0])
calc_btn = widgets.Button(description="Faire le calcul", button_style='success')
calc_btn.on_click(update)

display(widgets.HTML("<h3>Points connus (max 7)</h3>"), box_pts)
display(widgets.HTML("<h3>Paramètres du variogramme imbriqué</h3>"), box_cov)
display(widgets.HTML("<h3>Paramètres de krigeage</h3>"), box_krigeage)
display(calc_btn, output)

Note de cours et ateliers interactifs

Chap8

Note de cours et ateliers interactifs

Atelier 2 - Krigeage Simple vs Ordinaire (1D)