TP3 - Analyse spectrale et modélisation

À rendre le 2026-04-12 - 14% du cours

Contexte
Objectifs
- Objectifs généraux
- Objectifs spécifiques
Données
Travail à faire
Tutoriel
Références

Contexte

Le TP3 porte sur la caractérisation gravimétrique d’un site de stockage géologique du CO$_2$, inspiré de l’étude de Andrés et al. (2016). Vous disposerez d’un levé gravimétrique, de trois contraintes de forage et d’un maillage de prismes rectangulaires déjà défini. Votre travail consistera à :

corriger les données gravimétriques,
calculer l’anomalie de Bouguer,
séparer la régionale et la résiduelle,
utiliser les dérivées spatiales pour mieux délimiter la structure,
construire un modèle 3D du socle à l’aide du maillage fourni,
estimer le volume du piège.

Le maillage 3D étant déjà fourni, vous n’avez pas à chercher la position latérale de chaque prisme. La géométrie en plan des cellules est fixée. L’inconnue principale est donc la profondeur du socle, ou l’épaisseur du piège, dans chaque cellule active.

L’objectif du TP n’est pas de faire une inversion gravimétrique avancée, mais de construire un modèle géologiquement plausible, cohérent avec les forages, compatible avec les dérivées spatiales et capable de reproduire raisonnablement l’anomalie résiduelle observée.

Objectifs

Objectifs généraux

Apprendre à filtrer des données géophysiques dans le domaine spectral.
Comprendre l’utilité et les limites des différents filtres géophysiques.
Relier une anomalie gravimétrique résiduelle à une géométrie 3D réaliste.

Objectifs spécifiques

Corriger les données gravimétriques fournies et obtenir l’anomalie de Bouguer.
Calculer l’anomalie régionale et l’anomalie résiduelle en effectuant un prolongement vers le haut de 100 m sur l’anomalie de Bouguer.
Calculer des dérivées spatiales utiles à la délimitation du piège.
Estimer la profondeur du socle dans un maillage 3D de prismes rectangulaires.
Estimer le volume du piège géologique et valider cette estimation par modélisation gravimétrique 3D.

Données

Veuillez télécharger les fichiers suivants :

TP3-GRAVITE-CO2.csv : levé gravimétrique.
TP3-GRAVITE-CO2-WELLS.csv : contraintes de forage.
TP3-GRAVITE-CO2-MESH.csv : maillage de modélisation 3D.

Toutes les données sont fournies dans les unités SI.

Le levé contient 10 201 stations espacées de 20 m sur une grille uniforme carrée.
La station de base est située aux coordonnées $(x, y)=(-1000, -1000)$.
La correction de dérive a déjà été effectuée.
Les corrections de latitude, d’altitude et de plateau n’ont pas été effectuées.
Les trois forages interceptent le toit du piège à 50 m de profondeur.
Le forage F1 est situé au centre d’une cellule du maillage aux coordonnées $(-208.33, 41.67)$ et confirme une épaisseur du piège de 75 m.
Le forage F2 est situé au centre d’une cellule du maillage aux coordonnées $(125.00, -208.33)$ et confirme une épaisseur du piège de 100 m.
Le forage F3 est situé au centre d’une cellule du maillage aux coordonnées $(41.67, 125.00)$ et confirme une épaisseur du piège de 125 m.
Le maillage de modélisation contient 144 cellules (12 x 12).
Certaines cellules sont inactives (Active = 0) et ne doivent pas être utilisées dans le modèle.
La profondeur du toit du piège (Top_depth_m) est imposée dans tout le maillage.
Le contraste de densité (Density_contrast_kg_m3) est imposé et identique dans toutes les cellules. Sa valeur est de -450 kg/m^3.

La topographie et l’emplacement des forages sur le site sont montrés à la Figure 1.

Carte topographique

Figure 1. Carte topographique du site à l’étude.

La gravité observée sur le site est montrée à la Figure 2.

Carte gravimétrique observée

Figure 2. Carte de la gravité observée sur le site à l’étude.

Travail à faire

Le TP doit être traité comme un tutoriel guidé. Vous devez reproduire les étapes suivantes dans l’ordre.

Étape 1 - Lire les trois fichiers

Chargez les trois fichiers dans pandas, puis vérifiez :

les colonnes disponibles,
le nombre de lignes,
les coordonnées minimales et maximales,
le nombre de cellules actives du maillage.

À la fin de cette étape, vous devriez être capables de répondre aux questions suivantes :

Combien y a-t-il de stations dans le levé ?
Combien y a-t-il de cellules au total dans le maillage ?
Combien de cellules sont actives ?
Dans quelles cellules tombent les trois forages ?

Étape 2 - Corriger les données

À partir de g_obs_SI, calculez :

la correction de latitude,
la correction d’altitude,
la correction de plateau,
l’anomalie de Bouguer.

À cette étape, enregistrez au minimum les colonnes suivantes dans votre DataFrame :

corr_lat
corr_alt
corr_pla
Bouguer

La carte finale attendue pour cette étape ressemble à la Figure 3.

Carte de Bouguer

Figure 3. Carte de l’anomalie de Bouguer attendue après les corrections gravimétriques.

Étape 3 - Isoler la résiduelle

À partir de l’anomalie de Bouguer :

calculez la régionale par prolongement vers le haut de 100 m,
calculez la résiduelle par différence.

À cette étape, enregistrez au minimum :

Regionale
Residuelle

Vos cartes devraient ressembler qualitativement aux Figure 4 et Figure 5.

Carte régionale

Figure 4. Exemple d’anomalie régionale obtenue par prolongement vers le haut de 100 m.

Carte résiduelle

Figure 5. Exemple d’anomalie résiduelle obtenue par soustraction de la régionale.

Étape 4 - Calculer les dérivées spatiales

À partir de l’anomalie résiduelle, calculez :

$\partial g / \partial x$,
$\partial g / \partial y$,
$\partial g / \partial z$,
$\partial^2 g / \partial z^2$,
le signal analytique.

À cette étape, enregistrez au minimum :

ddx
ddy
ddz
ddz2
SA

La Figure 6 montre un exemple de signal analytique qui peut vous aider à valider cette étape.

Signal analytique

Figure 6. Exemple de signal analytique calculé à partir de l’anomalie résiduelle.

Étape 5 - Interpréter le piège

À partir de la résiduelle et des dérivées :

délimitez la zone du piège,
vérifiez la cohérence avec la position des trois forages,
identifiez quelles cellules actives du maillage doivent contenir le piège.

Concrètement, à la fin de cette étape, vous devriez avoir :

une idée de la forme générale du piège;
une idée de son extension latérale;
une hypothèse géologique raisonnable sur les cellules les plus épaisses.

Étape 6 - Construire un premier modèle du socle

Attribuez une épaisseur à chaque cellule active :

imposez les épaisseurs connues dans les cellules des forages,
proposez une interpolation ou une distribution plausible pour les autres cellules,
convertissez cette épaisseur en profondeur du socle.

Le but de cette étape est d’obtenir un modèle initial, pas un modèle parfait. Ce premier modèle doit surtout :

être lisse;
respecter exactement les trois forages, mais de façon lissée dans l’espace;
rester compatible avec la forme de l’anomalie résiduelle.

Votre carte finale du socle devrait avoir une allure comparable à la Figure 7.

Carte du socle

Figure 7. Exemple de carte finale de profondeur du socle obtenue à partir du modèle 3D.

Étape 7 - Calculer la réponse gravimétrique du modèle

Calculez la réponse gravimétrique du maillage de prismes.

Le modèle de prismes représente déjà la source locale, mais l’anomalie résiduelle observée est un champ filtré. Pour comparer correctement le modèle et les données, vous devez donc appliquer au modèle le même traitement spectral qu’aux données. Autrement dit, la comparaison doit se faire entre deux quantités traitées de la même manière :

la résiduelle observée obtenue à partir de la Bouguer,
la résiduelle modélisée obtenue à partir de la gravité du maillage de prismes.

À la fin de cette étape, vous devez avoir :

une carte de la gravité modélisée;
une carte de la résiduelle modélisée;
une carte ou une mesure numérique du misfit.

Étape 8 - Ajuster le modèle

Ajustez les profondeurs du socle jusqu’à obtenir :

une géométrie plausible,
un accord raisonnable entre résiduelle observée et résiduelle modélisée,
un modèle compatible avec les forages.

Dans ce TP, l’« inversion » attendue doit être simple et explicite. Vous n’avez pas à coder une méthode d’optimisation avancée. La stratégie recommandée est la suivante :

construire un modèle initial du socle à partir des trois forages;
calculer la résiduelle modélisée de ce modèle initial;
effectuer une recherche exhaustive sur un facteur d’échelle global appliqué aux épaisseurs;
conserver le facteur qui minimise le misfit quadratique moyen entre résiduelle observée et résiduelle modélisée;
vérifier que les cellules des forages restent compatibles avec les contraintes géologiques;
produire la carte finale du socle et le volume.

Autrement dit, on vous demande ici une inversion très simple à un paramètre, appuyée sur un modèle initial géologiquement plausible.

Le point important est le suivant : vous ne devez pas tester des profondeurs complètement aléatoires cellule par cellule. La procédure attendue est :

construire une géométrie initiale plausible;
conserver cette géométrie générale;
ajuster son amplitude globale par recherche exhaustive;
choisir le meilleur facteur avec un critère quantitatif, par exemple le RMSE.

Étape 9 - Calculer le volume

Calculez le volume total du piège à partir des épaisseurs finales.

Tutoriel

Librairies

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

Vous pouvez aussi utiliser pathlib :

from pathlib import Path

root = Path(".")
df = pd.read_csv(root / "TP3-GRAVITE-CO2.csv")
wells = pd.read_csv(root / "TP3-GRAVITE-CO2-WELLS.csv")
mesh = pd.read_csv(root / "TP3-GRAVITE-CO2-MESH.csv")

Lecture des données

df = pd.read_csv("TP3-GRAVITE-CO2.csv")  # leve gravimetrique
wells = pd.read_csv("TP3-GRAVITE-CO2-WELLS.csv")  # contraintes de forage
mesh = pd.read_csv("TP3-GRAVITE-CO2-MESH.csv")  # maillage de modelisation

print(df.head())  # apercu du leve
print(wells)  # tableau des forages
print(mesh.head())  # apercu du maillage

Colonnes utiles :

Dans df : X, Y, Z, g_obs_SI
Dans wells : Well, X, Y, Trap_top_m, Trap_thickness_m, Socle_depth_m
Dans mesh : Cell_ID, X1, X2, Y1, Y2, X_center, Y_center, Active, Top_depth_m, Density_contrast_kg_m3

Préparer les grilles

Avant de faire des cartes ou des filtres, remettez les données en grille.

nx = df["X"].nunique()  # nombre de colonnes
ny = df["Y"].nunique()  # nombre de lignes

xx = df["X"].values.reshape(ny, nx)  # grille X
yy = df["Y"].values.reshape(ny, nx)  # grille Y
zz = df["Z"].values.reshape(ny, nx)  # grille topographique
gobs = df["g_obs_SI"].values.reshape(ny, nx)  # grille de gravite observee

Carte de couleur

Utilisez les conventions suivantes :

rainbow pour les cartes gravimétriques,
terrain pour la topographie,
Blues_r pour la profondeur du socle.

Exemple :

plt.figure(figsize=(8, 6))  # taille de la figure
plt.pcolormesh(xx, yy, gobs, cmap="rainbow")  # carte coloree
plt.gca().set_aspect("equal")  # meme echelle en x et y
plt.colorbar(label="$g_z$ (m/s$^2$)")  # barre de couleurs
plt.xlabel("$x$ (m)")  # axe x
plt.ylabel("$y$ (m)")  # axe y

Corrections gravimétriques

Les constantes à utiliser sont :

G = 6.6743e-11  # constante gravitationnelle
rho_plateau = 2400  # densite du plateau de Bouguer
grad_alt = 3.083e-6  # gradient de l'air libre
g_e = 9.7803267715  # gravite normale a l'equateur
c1 = 5.2790414e-3  # coefficient de la correction de latitude
R_e = 6.378e6  # rayon equatorial de la Terre
lat_hontomin_deg = 42.5716  # latitude de Hontomin

Calcul :

z0 = df["Z"].iloc[0]  # altitude de reference
y0 = df["Y"].iloc[0]  # coordonnee y de reference
lat0 = np.radians(lat_hontomin_deg)  # latitude de Hontomin en radians

df["corr_lat"] = -(g_e * c1 * np.sin(2 * lat0) / R_e) * (df["Y"] - y0)  # correction de latitude
df["corr_alt"] = -grad_alt * (df["Z"] - z0)  # correction d'altitude
df["corr_pla"] = 2 * np.pi * G * rho_plateau * (df["Z"] - z0)  # correction de plateau

df["Bouguer"] = (
    df["g_obs_SI"]  # gravite observee
    + grad_alt * (df["Z"] - z0)  # retrait de l'altitude
    - 2 * np.pi * G * rho_plateau * (df["Z"] - z0)  # retrait du plateau
    + (g_e * c1 * np.sin(2 * lat0) / R_e) * (df["Y"] - y0)  # retrait de la latitude
)  # anomalie de Bouguer

Ici, utilisez explicitement la formule du cours $\Delta g_\mathrm{L} = -\dfrac{g_\mathrm{e} c_1 \sin(2\phi)}{R_\mathrm{e}} (y - y_0)$ avec la latitude de Hontomín.

Fonction pour calculer l’anomalie régionale

def calcul_regionale(x, y, gz, h):
    nx = x.nunique()  # taille en x
    ny = y.nunique()  # taille en y
    grille = gz.values.reshape(ny, nx)  # remise en grille
    grille_pad = np.pad(grille, ((ny, ny), (nx, nx)), mode="reflect")  # ajout d'une bordure

    gz_hat = np.fft.fft2(grille_pad)  # transformee de Fourier

    dx = x.diff().dropna().max()  # pas en x
    dy = y.diff().dropna().max()  # pas en y
    kx = 2 * np.pi * np.fft.fftfreq(3 * nx, d=dx)  # nombres d'onde en x
    ky = 2 * np.pi * np.fft.fftfreq(3 * ny, d=dy)  # nombres d'onde en y
    KXX, KYY = np.meshgrid(kx, ky)  # grille spectrale

    filtre = np.exp(-h * np.sqrt(KXX**2 + KYY**2))  # filtre de prolongement vers le haut
    regionale_pad = np.fft.ifft2(gz_hat * filtre).real  # retour dans l'espace
    regionale = regionale_pad[ny:-ny, nx:-nx]  # retrait de la bordure
    return regionale.flatten()  # retour en vecteur

Utilisation :

df["Regionale"] = calcul_regionale(df["X"], df["Y"], df["Bouguer"], h=100)  # champ regional
df["Residuelle"] = df["Bouguer"] - df["Regionale"]  # champ local

Dans ce TP, la hauteur de prolongement demandée est 100 m.

Dérivées spatiales

Exemple de dérivée verticale :

def derivee_verticale(x, y, gz, ordre=1):
    nx = x.nunique()  # taille en x
    ny = y.nunique()  # taille en y
    grille = gz.values.reshape(ny, nx)  # remise en grille
    grille_pad = np.pad(grille, ((ny, ny), (nx, nx)), mode="reflect")  # ajout d'une bordure
    gz_hat = np.fft.fft2(grille_pad)  # transformee de Fourier

    dx = x.diff().dropna().max()  # pas en x
    dy = y.diff().dropna().max()  # pas en y
    kx = 2 * np.pi * np.fft.fftfreq(3 * nx, d=dx)  # nombres d'onde en x
    ky = 2 * np.pi * np.fft.fftfreq(3 * ny, d=dy)  # nombres d'onde en y
    KXX, KYY = np.meshgrid(kx, ky)  # grille spectrale

    filtre = np.sqrt(KXX**2 + KYY**2) ** ordre  # filtre de derivee verticale
    derivee_pad = np.fft.ifft2(gz_hat * filtre).real  # retour dans l'espace
    derivee = derivee_pad[ny:-ny, nx:-nx]  # retrait de la bordure
    return derivee.flatten()  # retour en vecteur

Exemple de dérivées horizontales :

def derivees_horizontales(x, y, gz):
    nx = x.nunique()  # taille en x
    ny = y.nunique()  # taille en y
    grille = gz.values.reshape(ny, nx)  # remise en grille
    grille_pad = np.pad(grille, ((ny, ny), (nx, nx)), mode="reflect")  # ajout d'une bordure
    gz_hat = np.fft.fft2(grille_pad)  # transformee de Fourier

    dx = x.diff().dropna().max()  # pas en x
    dy = y.diff().dropna().max()  # pas en y
    kx = 2 * np.pi * np.fft.fftfreq(3 * nx, d=dx)  # nombres d'onde en x
    ky = 2 * np.pi * np.fft.fftfreq(3 * ny, d=dy)  # nombres d'onde en y
    KXX, KYY = np.meshgrid(kx, ky)  # grille spectrale

    ddx_pad = np.fft.ifft2(gz_hat * (1j * KXX)).real  # derivee selon x
    ddy_pad = np.fft.ifft2(gz_hat * (1j * KYY)).real  # derivee selon y

    ddx = ddx_pad[ny:-ny, nx:-nx]  # retrait de la bordure
    ddy = ddy_pad[ny:-ny, nx:-nx]  # retrait de la bordure
    return ddx.flatten(), ddy.flatten()  # retour en vecteurs

Puis :

df["ddz"] = derivee_verticale(df["X"], df["Y"], df["Residuelle"], ordre=1)  # 1re derivee verticale
df["ddz2"] = derivee_verticale(df["X"], df["Y"], df["Residuelle"], ordre=2)  # 2e derivee verticale
df["ddx"], df["ddy"] = derivees_horizontales(df["X"], df["Y"], df["Residuelle"])  # derivees horizontales
df["SA"] = np.sqrt(df["ddx"]**2 + df["ddy"]**2 + df["ddz"]**2)  # signal analytique

Remarque importante :

Les dérivées servent surtout à interpréter la géométrie du piège.
Elles ne servent pas directement à calculer les profondeurs cellule par cellule.
Le signal analytique est particulièrement utile pour visualiser les limites du corps.

Gravité d’un prisme rectangulaire

La composante verticale de la gravité mesurée aux coordonnées d’observation $x, y, z$ et causée par un prisme rectangulaire de contraste de densité $\Delta\rho$ dont les côtés sont aux coordonnées $x_i, y_j, z_k$ est (Li & Chouteau, 1998)

$\Delta g_z = G \Delta\rho \sum_{i=1}^2 \sum_{j=1}^2 \sum_{k=1}^2 \mu_{ijk} \left[ \Delta x_i \ln{}(\Delta y_j + r_{ijk}) + \Delta y_j \ln{}(\Delta x_i + r_{ijk}) - \Delta z_k \arctan{\frac{\Delta x_i \Delta y_j}{\Delta z_k r_{ijk}}} \right]$,

où $G$ est la constante gravitationnelle,

$ \Delta x_i = x - x_i $,

$ \Delta y_j = y - y_j $,

$ \Delta z_k = z - z_k $,

$ r_{ijk} = \sqrt{\Delta x_i^2 + \Delta y_j^2 + \Delta z_k^2} $

$ \mu_{ijk} = (-1)^i (-1)^j (-1)^k $.

Gabarit Python :

def calculer_g_prisme(rho, x12, y12, z12, x, y, z):
    G = 6.6743e-11  # constante gravitationnelle
    gz = np.zeros_like(x, dtype=float)  # accumulation de la gravite

    for i in range(2):
        for j in range(2):
            for k in range(2):
                dx = x - x12[i]  # distance en x
                dy = y - y12[j]  # distance en y
                dz = z - z12[k]  # distance en z
                r = np.sqrt(dx**2 + dy**2 + dz**2)  # distance totale
                mu = (-1)**i * (-1)**j * (-1)**k  # signe du coin
                gz += mu * (
                    dx * np.log(dy + r)  # terme logarithmique 1
                    + dy * np.log(dx + r)  # terme logarithmique 2
                    - dz * np.arctan2(dx * dy, dz * r)  # terme arctangente
                )

    return -G * rho * gz  # gravite verticale du prisme

Construire le modèle 3D du socle

Le maillage fixe la géométrie latérale. Il vous reste à attribuer une épaisseur à chaque cellule active.

Commencez par :

garder les cellules actives;
créer une colonne d’épaisseur;
construire une première interpolation lisse à partir des trois forages;
imposer ensuite les forages par une correction gaussienne locale, plutôt qu’en remplaçant brutalement une seule cellule.

Exemple :

mesh_actif = mesh[mesh["Active"] == 1].copy()  # garder seulement les cellules actives
mesh_actif["Epaisseur_m"] = 0.0  # initialiser les epaisseurs

Ensuite, vous devez attribuer les autres épaisseurs. Une manière simple de démarrer consiste à interpoler depuis les forages :

for i, cellule in mesh_actif.iterrows():
    dx = wells["X"] - cellule["X_center"]  # distance aux forages en x
    dy = wells["Y"] - cellule["Y_center"]  # distance aux forages en y
    distance2 = dx**2 + dy**2  # distance au carre
    poids = np.exp(-distance2 / (2 * 220.0**2))  # poids gaussien

    mesh_actif.loc[i, "Epaisseur_m"] = np.sum(poids * wells["Trap_thickness_m"]) / np.sum(poids)  # moyenne ponderee

Pour respecter exactement les forages sans créer de pics trop brusques, appliquez ensuite une correction locale lissée :

def imposer_forages_en_douceur(mesh_modele, wells, sigma=150.0):
    mesh_modele = mesh_modele.copy()

    for _, well in wells.iterrows():
        masque = (
            np.isclose(mesh_modele["X_center"], well["X"])
            & np.isclose(mesh_modele["Y_center"], well["Y"])
        )
        epaisseur_actuelle = mesh_modele.loc[masque, "Epaisseur_m"].iloc[0]
        correction = well["Trap_thickness_m"] - epaisseur_actuelle

        dx = mesh_modele["X_center"] - well["X"]
        dy = mesh_modele["Y_center"] - well["Y"]
        distance2 = dx**2 + dy**2
        poids = np.exp(-distance2 / (2 * sigma**2))
        poids = poids / poids[masque].iloc[0]

        mesh_modele["Epaisseur_m"] += correction * poids

    return mesh_modele

mesh_actif = imposer_forages_en_douceur(mesh_actif, wells, sigma=150.0)

Ce n’est qu’un point de départ. Vous devrez ensuite ajuster le modèle en fonction :

de la position de l’anomalie résiduelle,
du signal analytique,
des dérivées spatiales,
de la géométrie plausible du piège.

Comment faire l’« inversion » de façon simple

Voici la procédure attendue pour ajuster votre modèle.

1. Construire un modèle initial

Le modèle initial doit :

respecter exactement les trois forages;
donner des épaisseurs nulles dans les cellules inactives;
être lisse et plausible dans les cellules actives.

2. Calculer la résiduelle modélisée du modèle initial

Calculez d’abord la gravité du maillage, puis appliquez le même prolongement vers le haut de 100 m qu’aux données pour obtenir une résiduelle modélisée comparable :

gz_modele = np.zeros(len(df))  # gravite totale du modele

for _, cellule in mesh_actif.iterrows():
    epaisseur = cellule["Epaisseur_m"]  # epaisseur de la cellule
    if epaisseur <= 0:  # ignorer les cellules sans piege
        continue

    x12 = [cellule["X1"], cellule["X2"]]  # limites en x
    y12 = [cellule["Y1"], cellule["Y2"]]  # limites en y
    z12 = [-(cellule["Top_depth_m"] + epaisseur), -cellule["Top_depth_m"]]  # base et toit du prisme

    gz_modele += calculer_g_prisme(
        cellule["Density_contrast_kg_m3"],  # contraste de densite
        x12,  # limites x
        y12,  # limites y
        z12,  # limites z
        df["X"].values,  # points d'observation x
        df["Y"].values,  # points d'observation y
        np.zeros(len(df)),  # points d'observation z
    )

gz_modele_series = pd.Series(gz_modele)  # conversion en serie pandas
regionale_modele = calcul_regionale(df["X"], df["Y"], gz_modele_series, h=100)  # meme filtrage que les donnees
residuelle_modele = gz_modele - regionale_modele  # champ local modele

3. Faire une recherche exhaustive sur un facteur d’échelle

L’idée est de garder la forme générale du modèle initial, mais de tester plusieurs amplitudes globales :

facteurs = np.linspace(0.2, 2.0, 181)  # facteurs testes

meilleur_facteur = None  # meilleur facteur trouve
meilleure_erreur = np.inf  # meilleure erreur trouvee

for facteur in facteurs:
    mesh_test = mesh_actif.copy()  # copie du modele initial
    mesh_test["Epaisseur_m"] = facteur * mesh_actif["Epaisseur_m"]  # mise a l'echelle globale
    mesh_test = imposer_forages_en_douceur(mesh_test, wells, sigma=150.0)  # forages imposes en douceur

    gz_test = np.zeros(len(df))  # gravite du modele teste

    for _, cellule in mesh_test.iterrows():
        epaisseur = cellule["Epaisseur_m"]  # epaisseur testee
        if epaisseur <= 0:  # ignorer les cellules vides
            continue

        x12 = [cellule["X1"], cellule["X2"]]  # limites x
        y12 = [cellule["Y1"], cellule["Y2"]]  # limites y
        z12 = [-(cellule["Top_depth_m"] + epaisseur), -cellule["Top_depth_m"]]  # limites z

        gz_test += calculer_g_prisme(
            cellule["Density_contrast_kg_m3"],
            x12,
            y12,
            z12,
            df["X"].values,
            df["Y"].values,
            np.zeros(len(df)),
        )

    residuelle_test = gz_test - calcul_regionale(df["X"], df["Y"], pd.Series(gz_test), h=100)  # residuelle du modele
    erreur = np.sqrt(np.mean((df["Residuelle"].values - residuelle_test) ** 2))  # RMSE

    if erreur < meilleure_erreur:
        meilleure_erreur = erreur  # mise a jour de la meilleure erreur
        meilleur_facteur = facteur  # mise a jour du meilleur facteur

print(meilleur_facteur, meilleure_erreur)  # facteur retenu et RMSE

Interprétation :

si facteur < 1, votre modèle initial était globalement trop épais;
si facteur > 1, votre modèle initial était globalement trop mince.

4. Construire le modèle final

Une fois le meilleur facteur trouvé :

mesh_final = mesh_actif.copy()  # copie du modele initial
mesh_final["Epaisseur_m"] = meilleur_facteur * mesh_actif["Epaisseur_m"]  # appliquer le meilleur facteur
mesh_final = imposer_forages_en_douceur(mesh_final, wells, sigma=150.0)  # forages imposes en douceur

mesh_final["Socle_m"] = mesh_final["Top_depth_m"] + mesh_final["Epaisseur_m"]  # profondeur finale du socle

5. Vérifier le résultat

Le but n’est pas d’obtenir un misfit nul, mais un modèle cohérent et défendable. Pour vous guider, posez-vous ces questions :

Les plus fortes épaisseurs sont-elles situées dans une zone cohérente avec les forages ?
Le socle varie-t-il de manière relativement lisse ?
Le modèle final reproduit-il la géométrie principale de la résiduelle ?
Le volume final est-il raisonnable compte tenu de l’étendue du piège ?

Calculer la gravité du maillage de prismes

gz_modele = np.zeros(len(df))  # gravite totale du maillage

for _, cellule in mesh_actif.iterrows():
    epaisseur = cellule["Epaisseur_m"]  # epaisseur locale
    if epaisseur <= 0:  # ignorer les cellules sans piege
        continue

    x12 = [cellule["X1"], cellule["X2"]]  # limites x
    y12 = [cellule["Y1"], cellule["Y2"]]  # limites y
    z12 = [-(cellule["Top_depth_m"] + epaisseur), -cellule["Top_depth_m"]]  # limites z

    gz_modele += calculer_g_prisme(
        cellule["Density_contrast_kg_m3"],
        x12,
        y12,
        z12,
        df["X"].values,
        df["Y"].values,
        np.zeros(len(df)),
    )

Comme l’anomalie résiduelle observée est filtrée, il faut appliquer exactement le même traitement spectral au modèle :

gz_modele_series = pd.Series(gz_modele)  # conversion en serie
regionale_modele = calcul_regionale(df["X"], df["Y"], gz_modele_series, h=100)  # meme filtrage
residuelle_modele = gz_modele - regionale_modele  # residuelle modelee

Le misfit peut ensuite être calculé :

erreur = df["Residuelle"].values - residuelle_modele  # ecart point par point
rmse = np.sqrt(np.mean(erreur**2))  # erreur quadratique moyenne
print(rmse)  # valeur du misfit

Vous pouvez aussi enregistrer l’erreur pour chaque facteur testé dans une liste et tracer la courbe RMSE en fonction du facteur. Cela permet de montrer clairement comment vous avez choisi le meilleur modèle.

Calcul du volume

dx = mesh["X2"].iloc[0] - mesh["X1"].iloc[0]  # largeur d'une cellule
dy = mesh["Y2"].iloc[0] - mesh["Y1"].iloc[0]  # hauteur d'une cellule
aire = dx * dy  # aire d'une cellule

volume = np.sum(mesh_actif["Epaisseur_m"] * aire)  # somme des volumes des cellules
print(volume)  # volume total du piege

Références

Andrés, J., Alcalde, J., Ayarza, P., Saura, E., Marzán, I., Martí, D., Martínez Catalán, J. R., Carbonell, R., Pérez-Estaún, A., García-Lobón, J. L., & Rubio, F. M. (2016). Basement structure of the Hontomín CO2 storage site (Spain) determined by integration of microgravity and 3-D seismic data. Solid Earth, 7(3), 827‑841. https://doi.org/10.5194/se-7-827-2016
Li, X., & Chouteau, M. (1998). Three-Dimensional Gravity Modeling In All Space. Surveys in Geophysics, 19(4), 339‑368. https://doi.org/10.1023/A:1006554408567