Optimisation numérique pour l’inversion EM — guide pas à pas (1D → 2D → 3D)

ien avant d’écrire la moindre ligne de code, posons une question simple : comment « voir » ce que l’on ne peut pas toucher ?

Imaginez que vous êtes sur un bateau dans un épais brouillard. Vous ne voyez pas le fond marin, mais vous avez un sonar. Vous envoyez un « ping » (une onde sonore) dans l’eau. C’est le problème direct : utiliser une physique connue pour prédire ce qui va se passer. Quelques secondes plus tard, un écho revient : c’est votre donnée. Votre cerveau fait immédiatement le calcul : « D’après le temps d’aller-retour et l’intensité, le fond est à 100 m et plutôt rocheux. » Cette opération mentale — convertir des données observées en un modèle invisible du monde — c’est l’essence de l’inversion : retrouver la cause (le modèle) à partir de l’effet (les mesures).

En effet, un médecin utilise un scanner (CT) pour la même raison : on ne peut pas voir l’intérieur du corps à l’œil nu, alors on fait passer des rayons X (physique), on mesure ce qui ressort (données), puis un ordinateur inverse ces mesures pour reconstruire un volume 3D de vos organes. L’inversion numérique est notre « scanner » pour la Terre : un processus mathématique et computationnel qui transforme des observations limitées et bruitées de surface en une image fiable — un modèle — du sous-sol invisible.

Plutôt que d’envoyer un « ping » sonar ou des rayons X, nous utilisons l’électromagnétisme (EM). Une inversion EM procède ainsi :

Ce n’est pas qu’un exercice académique. L’inversion EM est le moteur qui permet de trouver :

Dans tous ces cas, le forage est coûteux et risqué. L’inversion nous permet de « forer » des milliers de fois à l’ordinateur d’abord, afin de ne forer qu’une seule fois sur le terrain. Maintenant, construisons le moteur.

eci correspond au « ping sonar » de notre problème. En inversion géophysique, il nous faut d’abord un modèle direct — la physique qui prédit ce que nous devrions mesurer. En EM harmonique (domaine fréquentiel), cette physique est régie par l’équation d’onde vectorielle.

Cette équation est notre bête de somme : elle décrit comment le champ électrique $\mathbf{E}$ (notre « ping ») interagit avec le sous-sol. L’élément clé est le terme $k^2(\mathbf{x})$, qui encapsule les propriétés du milieu : $k^2(\mathbf{x}) = \omega^2\mu\epsilon(\mathbf{x}) - i\omega\mu\sigma(\mathbf{x})$. Cela signifie que si l’on modifie la conductivité $\sigma$ (notre modèle), l’équation produira une réponse prédite différente (nos données).

Toute cette simulation physique se résume par une fonction $\mathbf{F}(\mathbf{m})$ : vous lui fournissez un modèle $\mathbf{m}$ (la carte de conductivité) et elle renvoie les données simulées $\mathbf{d}_{\text{sim}} = \mathbf{F}(\mathbf{m})$. Voilà pour le problème direct. Pour le problème inverse, nous disposons déjà des mesures réelles $\mathbf{d}_{\text{obs}}$. Notre objectif est de trouver le modèle $\mathbf{m}$ qui soit :

On impose cet équilibre en définissant un « score » — la fonction objectif $\Phi$ — que l’on cherche à rendre aussi faible que possible. Elle comporte deux termes :

Ici, $\Phi_d$ est l’écart aux données (data misfit), qui mesure l’erreur entre prédiction et observation. $\Phi_m$ est la régularisation (ou « coût du modèle »), qui pénalise les modèles trop rugueux ou trop complexes. La valeur $\beta$ joue le rôle d’un curseur qui équilibre les deux : privilégie-t-on un ajustement parfait aux données, ou un modèle plus lisse ?

Le but de l’inversion devient alors simple : trouver le modèle $\mathbf{m}$ qui minimise $\Phi$. Les prochaines étapes présenteront les algorithmes précis (Gauss–Newton et Levenberg–Marquardt) utilisés pour atteindre ce minimum.

Soyons concrets. Notre $\mathbf{F}(\mathbf{m})$ n’est pas une idée en l’air : c’est une physique précise. Pour la plupart des levés EM, le cheval de bataille est la forme diffusive des équations de Maxwell — une simplification courante de l’équation d’onde complète de l’Éq. (1). Aux basses fréquences usuelles en géophysique, le courant de conduction ($i\omega\sigma\mathbf{E}$) domine, tandis que le courant de déplacement ($\omega^2\epsilon\mathbf{E}$) est si faible qu’on peut l’ignorer. On obtient alors l’EDP pratique suivante :

Cette équation est le cœur de $\mathbf{F}(\mathbf{m})$ : elle relie directement la propriété que l’on cherche — la conductivité $\sigma(\mathbf{x})$ (notre modèle $\mathbf{m}$) — au champ électrique $\mathbf{E}$ que l’on peut mesurer. Cette physique « diffusive » fournit aussi une règle d’or essentielle : la profondeur de peau, $\displaystyle \delta \approx \sqrt{\tfrac{2}{\mu\,\sigma\,\omega}}$.

La profondeur de peau $\delta$ indique jusqu’où notre « ping » EM peut pénétrer. Les structures bien plus petites que $\delta$ seront mal résolues (floutées), et celles bien plus profondes que $\delta$ ne renverront pratiquement aucun écho. Cette formule simple sert de premier reality check pour concevoir un levé et choisir si un modèle 1D ou 2D suffit, ou s’il faut passer à une simulation 3D complète. (Astuce débutant : plus la conductivité $\sigma$ ou la fréquence $\omega$ est grande, plus $\delta$ est petite ; on « voit » moins profond.)

n dispose d’un modèle direct (Éq. 3) qui prédit les champs électrique $\mathbf{E}$ et magnétique $\mathbf{H}$ partout dans le sous-sol. Mais nous ne pouvons pas tout mesurer : seulement quelques composantes en des points récepteurs à la surface. Les instruments enregistrent ces champs de surface et les convertissent en rapports robustes et interprétables. En MT/AMT/CSAMT, cela nous fournit deux observables clés :

D’autres méthodes (FDEM, TDEM) peuvent enregistrer des composantes de champ directement. Mais dans tous les cas, qu’il s’agisse de $\mathbf{Z}$, de $\mathbf{T}$ ou de champs bruts, ce sont simplement des « vues » différentes d’une même physique sous-jacente. Ce sont les données observées $\mathbf{d}_{\text{obs}}$ que l’inversion cherchera à reproduire.

Notre modèle direct $\mathbf{F}(\mathbf{m})$ résout l’EDP (Éq. 3). Son coût et sa complexité dépendent entièrement de la dimensionnalité — c’est-à-dire des hypothèses géométriques que l’on fait sur la Terre.

Point crucial : l’algorithme d’inversion que nous allons construire reste agnostique à ce choix. Seul change le « moteur physique » — le solveur direct 1D, 2D ou 3D que l’on y branche.

On inverse pour $\mathbf{m}=\log\sigma$ afin d’imposer la positivité (via l’exponentielle) et de rendre les pas d’optimisation multiplicatifs en $\sigma$. Les maillages numériques (couches 1D, cellules 2D/3D) sont des grilles de modèle ; les données du levé (impédances, champs) sont prédites aux emplacements récepteurs par interpolation et évaluation fonctionnelle. Les sensibilités sont prises par rapport à $\mathbf{m}$ : le jacobien $\mathbf{J}=\partial\mathbf{d}/\partial\mathbf{m}$ couple la physique à l’optimisation. (Astuce débutant : travailler en $\log\sigma$ évite les valeurs négatives de conductivité et stabilise souvent la recherche.)

L’inversion équilibre deux forces : ajuster les données fiables et rester géologiquement plausible. On code cela comme la somme d’un écart aux données pondéré et d’un régularisateur de modèle :

$\mathbf W_d$ met à l’échelle chaque donnée selon son incertitude (et un plancher d’erreur optionnel). Si l’écart-type $\sigma_i$ est connu, un choix classique est $(\mathbf W_d)_{ii}=1/\sigma_i$, de sorte qu’un « bon » ajustement correspond à un $\chi^2$ normalisé proche de 1.

$\mathbf W_m$ favorise la structure attendue dans le sous-sol. En régularisation de Tikhonov lisse, $\mathbf W_m$ empile des différences d’ordres 0/1/2, par ex. $\mathbf W_m=\operatorname{diag}(\alpha_s\,\mathbf I,\ \alpha_x\,\mathbf D_x,\ \alpha_y\,\mathbf D_y,\ \alpha_z\,\mathbf D_z)$. On peut basculer en variation totale (TV) pour des contacts plus nets, et inclure un modèle de référence $\mathbf m_{\mathrm{ref}}$ pour tirer la solution vers une géologie connue.

$\beta$ arbitre entre l’ajustement et la douceur. On commence souvent avec un $\beta$ élevé (stabilité d’abord), puis on le diminue (« refroidissement ») jusqu’à atteindre le niveau de bruit des données — c’est le principe de la discordance, visant $\Phi_d \approx N_{\text{data}}/2$.

Pour employer des solveurs à gradient, il nous faut la dérivée de $\Phi$ par rapport au modèle :

Ici $\mathbf J=\partial \mathbf F/\partial \mathbf m$ est la sensibilité (jacobien) et $\mathbf r$ le résidu. En 1D, on peut souvent former $\mathbf J$ explicitement ; en 2D/3D, on applique $\mathbf J$ et $\mathbf J^\top$ de façon implicite via des résolutions directe/adjointes — essentiel pour la mémoire et la vitesse.

L’écart quadratique est optimal pour un bruit gaussien. En présence de pics ou de décalages statiques, on remplace la norme au carré par une pénalité robuste $\Phi_d^{\rho}=\sum_i \rho\!\left((\mathbf W_d\,\mathbf r)_i\right)$, par ex. Huber ou Tukey. On conserve la sensibilité près de zéro tout en limitant l’influence des outliers, ce qui améliore la convergence et la stabilité du modèle.

• Définir des planchers d’erreur raisonnables dans $\mathbf W_d$ (p. ex., 5–10 % en amplitude, 1–2° en phase).
• Choisir $\mathbf W_m$ en accord avec la géologie (lisse vs. « blocs ») et fournir $\mathbf m_{\mathrm{ref}}$ lorsqu’une information a priori fiable est disponible.
• Programmer $\beta$ : commencer haut, puis réduire par un facteur 0.5–0.8 jusqu’à ce que $\Phi_d$ atteigne le niveau de bruit.
• Travailler en $\mathbf m=\log\sigma$ pour garantir la positivité ; appliquer des bornes par projection si nécessaire.

Nous allons maintenant transformer la fonction objectif en algorithme. L’idée centrale est simple : linéariser la carte directe autour du modèle courant, résoudre une mise à jour de modèle qui fait baisser l’objectif, et recommencer. Cela conduit au système Gauss–Newton / Levenberg–Marquardt (GN/LM) :

La matrice de gauche mêle trois ingrédients : la courbure « donnée » $\mathbf J^\top \mathbf W_d^\top \mathbf W_d\,\mathbf J$ (comment les variations du modèle affectent les données), la courbure de régularisation $\beta\,\mathbf W_m^\top \mathbf W_m$ (comment elles affectent la rugosité), et l’amortissement $\mu\,\mathbf I$ qui contrôle la taille du pas. Quand $\mu\to 0$ on retrouve le pur Gauss–Newton ; quand $\mu$ est grand, le pas se comporte comme une descente de gradient prudente et bien conditionnée. Après avoir résolu pour $\Delta\mathbf m$, on accepte un pas échelonné via une recherche linéaire ou une région de confiance, puis on ajuste $\beta$ et $\mu$ de manière adaptative. (Astuce débutant : pensez à $\mu$ comme à un bouton « prudence ». Plus il est élevé, plus on avance doucement.)

Dans la Étape 2, nous dériverons soigneusement ce système et nous l’implémenterons en pratique :

À la fin de l’Étape 2, vous aurez une boucle GN/LM opérationnelle qui atteint une misfit cible (principe de la discordance), respecte les bornes en $\mathbf m=\log\sigma$, et produit des modèles à la fois géologiquement plausibles et quantitativement justifiés.

ous transformons l’objectif en algorithme. Linéarisons la carte directe autour du modèle courant : $\mathbf F(\mathbf m+\Delta\mathbf m)\approx \mathbf F(\mathbf m)+\mathbf J\,\Delta\mathbf m$. Avec cette approximation quadratique, la minimisation de $\Phi$ conduit aux équations normales :

Le membre de droite est le gradient $\nabla\Phi(\mathbf m)=\mathbf J^\top\mathbf W_d^\top\mathbf W_d\,\mathbf r +\beta\,\mathbf W_m^\top\mathbf W_m(\mathbf m-\mathbf m_{\mathrm{ref}})$ avec le résidu $\mathbf r=\mathbf F(\mathbf m)-\mathbf d$. Les trois termes de gauche sont la courbure données, la courbure régularisation et l’amortissement $\mu$. Un petit $\mu$ donne un pas de Gauss–Newton ; un grand $\mu$ se comporte comme un pas de gradient prudent.

Résolvez pour $\Delta\mathbf m$ (généralement avec CG/LSQR, en matrix-free). Choisissez ensuite une longueur de pas $\alpha\in(0,1]$ via une recherche linéaire d’Armijo avec backtracking, en appliquant les bornes dans $\mathbf m=\log\sigma$ si elles existent. Les variantes LM en région de confiance adaptent $\mu$ au moyen d’un rapport de gain :

Si $\rho$ est élevé (le modèle quadratique prédit bien la baisse), diminuez $\mu$ et autorisez des pas plus audacieux ; si $\rho$ est faible, augmentez $\mu$ et réessayez avec un pas plus court. Cela stabilise l’itération même quand la linéarisation est imparfaite.

Démarrez avec un $\beta$ plus grand pour stabiliser les premières itérations, puis réduisez $\beta\leftarrow c\,\beta$ avec $c\in[0.5,0.8]$ une fois que la misfit baisse. La cible est le niveau de bruit : $\Phi_d \approx N_{\text{data}}/2$ pour un bon ajustement.

Cette boucle 1D est volontairement compacte : elle reflète exactement la logique que nous réutiliserons en 2D/3D, où des produits matrix-free remplacent les matrices explicites et où un CG préconditionné résout les équations normales. Ensuite, nous brancherons ce moteur sur un modèle direct 1D en couches pour retrouver un bassin conducteur simple — puis nous passerons à l’échelle en 2D et en 3D.

ettons à l’épreuve l’algorithme de l’Étape 2. Nous définissons un « vrai » modèle à trois couches : une couche résistive (100 $\Omega\cdot$m), puis un bassin conducteur (10 $\Omega\cdot$m), au-dessus d’un socle résistif (1000 $\Omega\cdot$m). Nous générons des données synthétiques avec notre solveur direct 1D, ajoutons 5 % de bruit gaussien pour simuler un levé réel, puis nous fournissons ces « observations » à notre fonction inv1d_gnlm. L’inversion démarre d’un demi-espace simple et non informé à 100 $\Omega\cdot$m.

Le solveur converge en quelques itérations. Il commence avec un $\beta$ élevé, ce qui fait émerger d’abord les grandes structures (on l’observe sur les modèles des premières itérations). À mesure que la misfit aux données $\Phi_d$ atteint sa cible (principe de la discordance), l’algorithme réduit automatiquement $\beta$ et permet d’ajuster les détails plus fins. Le modèle final résout clairement le bassin conducteur, montrant que notre boucle d’optimisation est stable et efficace. Astuce débutant : un $\beta$ trop faible trop tôt peut conduire à sur-ajuster le bruit ; commencez « stable », puis desserrez.

oici l’idée la plus puissante en inversion numérique : l’algorithme d’optimisation écrit pour le 1D ne change pas quand on passe en 2D ou 3D. Le pseudocode inv1d_gnlm.py, avec sa boucle while, sa recherche linéaire d’Armijo et son planning de $\beta$, est exactement le même moteur que nous utilisons en 3D avec des millions de paramètres.

La seule chose qui change, c’est la manière de calculer quatre composants clés :

ous avons construit un modèle complet de l’inversion EM. Nous sommes partis de la physique directe (Éq.1) avons défini une fonction objectif pondérée et stable (Éq.2), puis avons dérivé l’algorithme de travail Gauss–Newton / Levenberg–Marquardt (Éq.6) pour la minimiser.

L’idée clé est que la logique d’optimisation est indépendante de l’implémentation physique. En utilisant les adjoints pour le gradient et des solveurs matrix-free pour l’étape de mise à jour, le même algorithme s’étend d’un modèle 1D calculé en une seconde à une inversion 3D en dix heures. Cette modularité — séparer minimize(Objective) de Objective.compute(gradient, residual) — est au cœur de tous les cadres modernes d’inversion géophysique.