3.5 Gradient spatial et gradiométrie
- Définition intuitive
- Calcul dans le domaine spatial
- Calcul dans le domaine spectral
- Exemple 1 : discrimination des sources superficielles
- Exemple 2 : localisation de contacts géologiques
- Conclusion
- Autres ressources
- Références
Définition intuitive
Le calcul des dérivées spatiales est une technique communément utilisée en géophysique pour faciliter l’interprétation des données gravimétriques et magnétométriques. La visualisation du gradient géophysique apporte une représentation supplémentaire des données qui tend à faire ressortir les bordures des sources d’anomalies. Les dérivées spatiales agissent comme des filtres passe-haut : elles amplifient les signaux à hautes fréquences par rapport aux basses fréquences dans le but de mettre en évidence les variations rapides dans les propriétés physiques du sol. Le gradient est donc utile pour des applications de délinéation des sources.
Les trois composantes du gradient spatial sont décrites ci-bas et la Figure 1 montre un exemple des dérivées horizontales d’une image. Remarquez bien l’effet de chacun des filtres selon l’orientation de la dérivée.
Figure 1. Calcul des gradients de premier ordre selon les deux dimensions de l’image.
Dérivée est-ouest (${\partial}/{\partial x}$)
Filtre passe-haut qui fait ressortir les variations locales dans la direction horizontale est-ouest. Pente du signal dans la direction $x$.
Dérivée nord-sud (${\partial}/{\partial y}$)
Filtre passe-haut qui fait ressortir les variations locales dans la direction horizontale nord-sud. Pente du signal dans la direction $y$.
Dérivée vertical (${\partial}/{\partial z}$)
Filtre passe-haut qui fait ressortir les variations locales. Corresponds approximativement à la différence entre le signal mesuré à deux hauteurs différentes divisée par cette hauteur :
\[\frac{\partial s}{\partial z} \approx \frac{s(z) - s(z + \Delta z)}{\Delta z} .\]Signal analytique
Le signal analytique $\vert A(x, y)\vert$ est défini comme la norme des gradients dans les trois directions, tel que
\[\vert A(x, y)\vert = \sqrt{\left(\frac{\partial s}{\partial x}\right)^2 + \left(\frac{\partial s}{\partial y}\right)^2 + \left(\frac{\partial s}{\partial z}\right)^2} .\]On utilise parfois aussi le gradient horizontal total $\vert H(x, y)\vert$ :
\[\vert H(x, y)\vert = \sqrt{\left(\frac{\partial s}{\partial x}\right)^2 + \left(\frac{\partial s}{\partial y}\right)^2} .\]Unités
Les unités du gradient spatial d’un champ de potentiel sont celles du champ par unité de distance. Dans la littérature géophysique le gradient de gravité est souvent exprimé en Eötvös.
- SI
- 1 Eötvös = $10^{-9}$ s$^{-2}$
- CGS
- 1 Eötvös = $10^{-9}$ Gal/cm$^2$
Calcul dans le domaine spatial
Soit $s(m, n)$ une carte géophysique obtenue suite à l’échantillonnage d’un signal géophysique sur une grille uniforme. Les indices $m$ et $n$ font respectivement référence aux indices des stations du levé dans la direction $x$ (p. ex. est-ouest) et $y$ (p. ex. nord-sud). De façon très générale, $s(m, n)$ est simplement une image de dimensions $n\times m$ pixels.
Dérivées horizontales
Les gradients se calculent dans le domaine spatial d’après leurs définitions intuitives. Prenons par exemple le gradient horizontal dans la direction $x$. Numériquement, le gradient correspond à la différence entre deux points consécutifs dans la direction $x$. On veut appliquer cette transformation pour toutes les stations ($m, n$).
Dans le domaine spatial, on peut écrire les gradients horizontaux comme des filtres qu’on applique à un signal $s$ avec un produit de convolution, dénoté $*$, tel que
\[\frac{\partial s}{\partial x} = \begin{bmatrix} -1 & 0 & 1 \end{bmatrix} * s ,\] \[\frac{\partial s}{\partial y} = \begin{bmatrix} -1\\ 0\\ 1 \end{bmatrix} * s .\]Les langages de programmation modernes ont souvent des librairies qui rendent facile l’application de ces opérations, par exemple en Python on utiliserait la démarche suivante pour calculer les gradients horizontaux.
import requests
from io import BytesIO
import numpy as np
from scipy import ndimage
def load_from_url(url):
img = Image.open(BytesIO(requests.get(url).content))
return np.array(img.convert('L'), dtype=float)
url = 'http://sipi.usc.edu/database/download.php?vol=misc&img=5.2.10'
data = load_from_url(url)
# Définir les filtres
ddx = np.array([[-1, 0, 1]]) # d/dx
ddy = np.array([[-1, 0, 1]]).T #d/dy
gx = ndimage.convolve(data, ddx) # produit de convolution
Essayez par vous-même. Vous pouvez remplacer le paramètre url pour pointer vers n’importe quelle image sur Internet.
Dérivée verticale
On ne peut pas obtenir le gradient vertical directement dans le domaine spatial. D’après l’équation (1), il faudrait avoir deux mesures du signal géophysique à chaque station : un à une hauteur de référence et une à une hauteur de $\Delta z$ par rapport à la première. Un levé qui procède à cette deuxième mesure simultanée pour obtenir le gradient vertical est un levé de gradiométrie.
Calcul dans le domaine spectral
La transformation de Fourier ($\mathcal{F}$) d’une fonction $f$ respecte la propriété de différentiation :
\[\mathcal{F} \left\{ \frac{\partial^n f(\boldsymbol{x})}{\partial x_j^n} \right\} = (ik_j)^n \ \hat{f}(\boldsymbol{k}) ,\]où $\boldsymbol{x}$ est un vecteur de position dans le domaine spatial et $\boldsymbol{k}$ est le vecteur des fréquences dans le domaine spectral. Ça rappelle aussi le théorème de convolution qu’on a vu précédemment, qui veut que la transformation de Fourier d’une convolution de deux signaux soit égale au produit de leurs transformations de Fourier respectives.
En effet, il est souvent plus facile de calculer les gradients spatiaux d’un signal en passant par sa transformation de Fourier dans le domaine spectral. Le protocole est le suivant :
- Transformer le signal d’origine vers le domaine spectral.
- Multiplier le spectre par la représentation spectrale d’un filtre.
- Calculer la transformation inverse du spectre filtré.
Dérivées horizontales
Sur une carte géophysique, le signal $s$ est échantillonné dans les directions horizontales. Il est donc assez facile de calculer les gradients horizontaux d’ordre $n$ avec la transformation de Fourier discrète $\hat{s}$ :
\[\frac{\partial^n s}{\partial x^n} = \mathcal{F}^{-1}\left\{ (ik_x)^n \ \hat{s}\right\}\]et
\[\frac{\partial^n s}{\partial y^n} = \mathcal{F}^{-1}\left\{ (ik_y)^n \ \hat{s}\right\} .\]Dérivée verticale
On a vu qu’on ne pouvait pas obtenir le gradient vertical directement dans le domaine spatial parce que sur une carte géophysique le signal est habituellement échantillonné dans les directions horizontales. Cependant, on sait que les potentiels gravitationnel et magnétique satisfont l’équation de Laplace et que :
\[\frac{\partial^2 s}{\partial x^2} + \frac{\partial^2 s}{\partial y^2} + \frac{\partial^2 s}{\partial z^2} = 0 .\]Les dérivées d’ordre 2 en $x$ et $y$ sont obtenues avec les équations (3) et (4). On obtient donc que
\[\frac{\partial^2 s}{\partial z^2} = - \left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}\right) s .\]On prend la transformation de Fourier de chaque côté :
\[\mathcal{F} \left\{ \frac{\partial^2 s}{\partial z^2} \right\} = - \mathcal{F} \left\{ \left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}\right) s \right\} ,\]donc
\[\mathcal{F} \left\{ \frac{\partial^2 s}{\partial z^2} \right\} = \left(k_x^2 + k_y^2 \right) \hat{s} .\]La dérivée verticale d’ordre 2 est finalement
\[\frac{\partial^2 s}{\partial z^2} = \mathcal{F}^{-1}\left\{ \left(k_x^2 + k_y^2 \right) \hat{s} \right\} .\]De façon plus générale on peut écrire la dérivée verticale d’ordre $n$ avec
\[\frac{\partial^n s}{\partial z^n} = \mathcal{F}^{-1}\left\{ \left(k_x^2 + k_y^2 \right)^{\frac{n}{2}} \hat{s} \right\} .\]Exemple 1 : discrimination des sources superficielles
Nous ferons maintenant un exemple de calcul des gradients verticaux d’ordre 1 et 2 avec la gravité produite par un modèle 3D de distribution de la densité.
Modèle pétrophysique
La Figure 2 montre deux sections différentes d’un modèle 3D de la densité du sous-sol. Une plaque de contraste de densité négatif -500 kg/m$^3$ est située à $z = -75$ m. Une sphère de contraste de densité positif +500 kg/m$^3$ est située au-dessus de la plaque à $z = -30$ m.
Figure 2. Vues en section verticale (gauche) et horizontale (droite) d’un modèle 3D de distribution de la densité.
Modélisation directe de la gravité
À cause du principe de superposition, la réponse gravimétrique mesurée en surface sera la somme des réponses individuelles de la plaque et de la sphère. La gravité mesurée à la surface ($z = 0$ m) est présentée à la Figure 3.
Figure 3. Réponse gravimétrique du modèle de densité contenant la plaque et la sphère.
On remarque que l’anomalie de Bouguer est légèrement asymétrique à cause de la sphère située au-dessus du côté droit de la plaque. Cependant, l’anomalie de Bouguer ne nous permet pas de faire une interprétation concluante quant à la géométrie exacte de la sphère et de la plaque. De plus, il n’est pas évident de savoir si la sphère se situe sous la plaque ou au-dessus de celle-ci à partir de ce signal.
Calcul du gradient vertical
On va maintenant calculer les gradients verticaux d’ordres 1 et 2 du signal (Figure 4).
Figure 4. Dérivée verticale d’ordre 1 (gauche) et d’ordre 2 (droite) du signal géophysique mesuré au-dessus du modèle 3D de la plaque et de la sphère.
Comme prévu, les gradients verticaux font ressortir les limites des sources qui se trouvent près de la surface par rapport aux corps plus profonds. La dérivée d’ordre 1 confirme que la sphère se situe au-dessus de la plaque et le gradient d’ordre 2 nous permet d’estimer que la sphère fait environ 20 m de diamètre, ce qui est exact d’après le modèle de densité utilisé pour générer le signal (Figure 2).
Exemple 2 : localisation de contacts géologiques
Dans cet exemple on utilise le gradient horizontal dans la direction $\boldsymbol{\hat{x}}$ pour déterminer l’emplacement exact d’un contact géologique (faille verticale, par exemple). La Figure 5 montre le gradient horizontal dans la direction $x$ de l’anomalie de Bouguer d’un contact géologique. Le gradient est calculé avec la fonction np.gradient(gz, x), où gz est un vecteur contenant les données de l’anomalie de Bouguer et x est le vecteur des positions des stations.
Figure 5. Anomalie de Bouguer et gradient horizontal de la gravité au-dessus d’un contact géologique.
Remarquez comment le point d’inflexion du champ gravitationnel correspondant au minimum de son gradient horizontal coïncide avec l’emplacement du contact.
Conclusion
En géophysique, le calcul des dérivées spatiales est utile pour :
- Faire ressortir les contrastes dans des directions spécifiques
- Mettre en évidence les anomalies près de la surface par rapport aux anomalies profondes
- Accentuer et/ou localiser les flancs d’une anomalie
On peut calculer le gradient d’un signal géophysique en suivant la démarche suivante :
- Transformer le signal vers le domaine spectral
- Multiplier le spectre par le filtre désiré
- Faire la transformation de Fourier inverse du spectre filtré
⚠️ Attention, il faut se souvenir que la qualité du signal reconstruit après le filtrage est fortement dépendante du nombre d’échantillons recueillis lors du levé géophysique. Il faut s’assurer que la configuration du levé par rapport aux anomalies anticipées respecte toujours le théorème d’échantillonnage.
Autres ressources
- La librairie Python SimPEG a été utilisée pour créer le modèle et calculer sa réponse gravimétrique.
Références
- Blakely, R. (1995). Transformations. Dans Potential Theory in Gravity and Magnetic Applications (pp. 311-358). Cambridge: Cambridge University Press.
- Elkins, T. A. (1951). The second derivative method of gravity interpretation. Geophysics, 16(1), 29-50.
- Tsuboi, C., & Kato, M. (1952). The First and Second Vertical Derivatives of Gravity. Journal of Physics of the Earth, 1(2), 95-96.