3.4 Séparation régionale–résiduelle

Sources du signal
- Exemples de sources régionales
- Exemples de sources locales
Prolongement des champs de potentiels
Exemple : cavité et contacts géologiques
Autres ressources
Références

Sources du signal

Nous avons vu que les contrastes dans les propriétés pétrophysiques sont la source des signaux géophysiques. Dans la nature, différents phénomènes peuvent se manifester par des contrastes de propriétés pétrophysiques à différentes échelles.

Par exemple, une suite de roches métamorphiques le long d’un gradient pression–température sur plusieurs km cause des variations très lentes de la densité du sous-sol. Cette source produit une anomalie gravimétrique qui est caractérisée par une grande longueur d’onde et une petite fréquence spatiale. À l’inverse, une cavité laissée dans la roche par des opérations minières antérieures cause un contraste de densité défini sur quelques mètres seulement. L’anomalie de Bouguer d’une telle cavité est caractérisée par une petite longueur d’onde et une fréquence spatiale bien plus élevée que celle du gradient métamorphique. Bien entendu, les notions de grande échelle (régionale) et petite échelle (locale) sont entièrement relatives aux dimensions de la grille de levé et à la fréquence d’échantillonnage.

Dans cette section, on verra comment séparer les sources régionales des sources locales. On appelle cette opération la séparation régionale–résiduelle. Cette opération se fait à l’aide d’un filtre qui peut être défini dans le domaine spectral : le filtre de prolongement vers le haut.

Exemples de sources régionales

Gradients métamorphiques
Contacts géologiques majeurs
Changements régionaux de l’épaisseur du mort-terrain
Proximité de batholithes ou plutons
Effet des marées (temporel mais affecte les données régionales)

Exemples de sources locales

Tube de lave
Aquifère
Galerie minière
Dykes (0.3 à 6 m de large en moyenne)
Formation de fer

Prolongement des champs de potentiels

Dans la pratique, on mesure les champs de potentiels sur la surface terrestre, soit au-dessus d’une distribution de propriétés pétrophysiques dans le sol. Il est possible de prolonger, ou d’extrapoler, la réponse de cette distribution à n’importe quel autre endroit, tant que le potentiel est harmonique dans la région où on fait le prolongement. Cette opération est possible parce que le champ gravitationnel obéit à l’équation de Laplace.

Comme on mesure déjà l’accélération gravitationnelle sur la surface terrestre dans les directions $\hat{x}$ et $\hat{y}$, il est particulièrement intéressant de faire des prolongements des champs de potentiels dans la direction verticale ($\hat{z}$). Comme on sait que l’intensité du champ gravitationnel varie en $1/r^2$, on peut se douter qu’il est possible d’atténuer ou d’amplifier le signal provenant des sources souterraines en effectuant des prolongements vers le haut ou vers le bas, respectivement.

Justification

Supposons que les contrastes de densité dans le sol peuvent être approximés comme une plaque dont la densité surfacique est $\sigma$, tel que montré à la Figure 1. Pour tenir compte des hétérogénéités dans les propriétés pétrophysiques du sol, la densité de la plaque varie selon les directions horizontales $x$ et $y$. On dénotera l’emplacement de chaque élément de la plaque avec les coordonnées $\xi$ et $\eta$. On s’intéresse à l’accélération gravitationnelle qu’on mesurerait au point $P$, situé en haut de la plaque.

svg

Figure 1. Schéma montrant l’approximation que la sous-surface est une plaque dont la densité surfacique est $\sigma$. Pour les besoins du développement mathématique $z$ est défini ici comme positif vers le bas.

La composante verticale de la gravité au point $P$ est

\[\Delta g_z(x, y, z_0-\Delta z) = G\Delta z\iint_{-\infty}^{\infty} \frac{\sigma(\xi, \eta)}{r^3} \,\mathrm{d}\xi\,\mathrm{d}\eta ,\]

où

\[r = \sqrt{(x-\xi)^2 + (y-\eta)^2 + \Delta z^2} .\]

Directement en haut de la plaque, pour toutes les sources situées à $z\lt z_0$, la composante verticale de l’accélération gravitationnelle est donnée par (pensez à la tranche de Bouguer) :

\[\Delta g_z(\xi, \eta) = 2\pi G \cdot \sigma(\xi, \eta) .\]

On peut donc réécrire la gravité au point $P$ comme :

\[\Delta g_z(x, y, z_0-\Delta z) = \frac{z}{2\pi}\iint_{-\infty}^{\infty} \frac{\Delta g_z(\xi, \eta)}{r^3} \,\mathrm{d}\xi\,\mathrm{d}\eta .\]

⚠️ L’équation prend la même forme pour une dérivation avec le champ magnétique.

L’équation (4) a une signification importante pour les méthodes géophysiques qui utilisent les champs de potentiels. C’est l’intégrale de prolongement. Elle implique qu’il est possible d’estimer le champ à n’importe quel point $P$, tant qu’il soit situé plus haut qu’un autre point où on a réellement mesuré $\Delta g_z$.

Cela veut dire que pour un levé géophysique réalisé au sol, on peut faire un prolongement vers le haut pour simuler quelle aurait été la réponse si on avait fait le levé du haut des airs. Inversement, on peut aussi simuler quelle aurait été la gravité mesurée en forage avec un prolongement vers le bas, quoique cette opération est très risquée comme nous le verrons dans ce qui suit.

Prolongement vers le haut

Pour alléger la notation dans ce qui suit on va dénoter le prolongement de la gravité $\Delta g_z(x, y, z_0-\Delta z)$ comme $\Delta g_P$. Reprenons à partir de l’équation (4) :

\[\Delta g_P = \frac{z}{2\pi}\iint_{-\infty}^{\infty} \frac{\Delta g_z(\xi, \eta)\,\mathrm{d}\xi\,\mathrm{d}\eta}{\left[(x-\xi)^2 + (y-\eta)^2 + \Delta z^2\right]^{3/2}} .\]

On reconnait bien évidemment la forme d’une opération de convolution ($*$) en 2D :

\[\Delta g_P = \iint_{-\infty}^{\infty} \Delta g_z(\xi, \eta)\ W(x-\xi, y-\eta)\,\mathrm{d}\xi\,\mathrm{d}\eta ,\]

où $W$ est la fonction de filtrage :

\[W(x-\xi, y-\eta) = \frac{z}{2\pi \left[(x-\xi)^2 + (y-\eta)^2 + \Delta z^2\right]^{3/2}} .\]

Le théorème de convolution dit que l’équation (6) est équivalente à

\[\mathcal{F}\left\{ \Delta g_P \right\} = \mathcal{F}\left\{ W \right\}\ \mathcal{F}\left\{ \Delta g_z(\xi, \eta) \right\} .\]

En prenant la transformation de Fourier 2D de chaque côté de l’équation (6) on obtient :

\[\iint_{-\infty}^{\infty} \Delta g_P e^{-i(k_x x + k_y y)}\,\mathrm{d}x\,\mathrm{d}y = \mathcal{F}\left\{ W \right\} \iint_{-\infty}^{\infty} \Delta g_z(\xi, \eta) e^{-i(k_x \xi + k_y \eta)}\,\mathrm{d}\xi\,\mathrm{d}\eta\]

et ainsi

\[\hat{g}_P(k_x, k_y) = \mathcal{F}\left\{ W \right\}\ \hat{g}_z(k_x, k_y) ,\]

où $\hat{g}_P$ et $\hat{g}_z$ sont respectivement les transformations de Fourier du prolongement de la gravité et de la gravité observée. Il ne reste qu’à calculer la transformation de Fourier du filtre $W$, qui donne

\[\mathcal{F}\left\{ W \right\} = \hat{W} = e^{-\Delta z \sqrt{k_x^2 + k_y^2}} .\]

On dit que $\hat{W}$ est la réponse spectrale du filtre de prolongement vers le haut. Grâce au théorème de convolution, le prolongement des champs de potentiel peut s’écrire comme une simple multiplication dans le domaine spectral.

Finalement, le prolongement du champ à une élévation $\Delta z$ par rapport au niveau de référence est

\[g_P(\Delta z) = \mathcal{F}^{-1} \left\{ e^{-\Delta z \sqrt{k_x^2 + k_y^2}}\ \hat{g}_z(k_x, k_y) \right\} .\]

Trois points importants découlent de cette équation :

Toutes les fréquences sont atténuées sauf les fréquences $(k_x,k_y) = (0,0)$.
Plus la fréquence d’une composante du signal est élevée, plus elle est atténuée.
Quand $\Delta z$ augmente l’atténuation augmente.

Le prolongement vers le haut est un filtre passe-bas qui atténue les hautes fréquences pour faire ressortir les anomalies régionales (celles qui ont des grandes longueurs d’ondes). Comme le prolongement vers le haut adoucit les données, cette opération est considérée comme stable : plus on s’éloigne des sources du champ plus le résultat du prolongement vers le haut tend vers la valeur moyenne des données.

Prolongement vers le bas

Jusqu’à maintenant on a seulement considéré un prolongement du champ dans une direction qui s’éloigne des sources. Qu’arrive-t-il si on tente de faire un prolongement vers le bas pour s’approche des sources? En d’autres mots, peut-on obtenir $\Delta g_z$ à partir de $\Delta g_P$? En principe cette opération est possible tant qu’on ne tente pas de faire un prolongement à l’intérieur de la région définie par la source. Prenons l’inverse de l’équation (10) :

\[\hat{g}_z(k_x, k_y) = \mathcal{F}^{-1}\left\{ W \right\}\ \hat{g}_P(k_x, k_y) ,\]

donc

\[\hat{g}_z(k_x, k_y) = e^{+\Delta z \sqrt{k_x^2 + k_y^2}}\ \hat{g}_P(k_x, k_y) .\]

Trois points importants :

Toutes les fréquences sont atténuées sauf les fréquences $(k_x,k_y) = (0,0)$.
Plus la fréquence d’une composante du signal est élevée, plus elle est amplifiée.
Quand $\Delta z$ augmente l’amplification augmente.

Le prolongement vers le bas est un filtre passe-haut qui amplifie les hautes fréquences et atténue les basses fréquences. Le prolongement vers le bas sous $z=0$ est possible si on arrive à démontrer que le potentiel gravitationnel est harmonique sous le sol, ce qui était une condition pour établir l’équation (4). C’est le cas partout en haut de la surface du sol, mais cette hypothèse ne tient rarement la route à cause des hétérogénéités dans le sol.

À cause de sa nature de filtre passe-haut, le prolongement vers le bas a pour effet de fortement amplifier le bruit dans les mesures. Il n’est donc pas toujours évident de l’utiliser en pratique. Le prolongement vers le bas n’est pas une opération stable parce que des petites variations dans la gravité à un point donné peuvent se traduire en d’énormes variations quand on se rapproche de la source. C’est à cause de l’exponentielle positive dans la représentation spectrale du filtre et de la dépendance en $1/r^2$ de $g$. Il faut souvent appliquer des filtres additionnels pour réduire le bruit généré par l’opération de prolongement vers le bas.

Séparation des sources

Il est possible de faire ressortir les anomalies locales en utilisant le principe de superposition et un prolongement vers le haut. Cette approche permet d’éviter l’amplification du bruit avec un prolongement vers le bas. En effet, puisque le prolongement vers le haut isole les sources régionales, on a qu’à le soustraire du signal original (l’anomalie de Bouguer) pour obtenir les sources locales. On appelle cette méthode la séparation régionale–résiduelle :

\[\Delta g_\mathrm{Résiduelle} = \Delta g_\mathrm{Bouguer} - \Delta g_\mathrm{Régionale} .\]

La séparation régionale–résiduelle est un filtre passe-haut.

Exemple : cavité et contacts géologiques

Les méthodes de prolongement ont une application évidente dans l’analyse des mesures géophysiques de champs de potentiels, comme en gravimétrie et en magnétométrie. Cette application est la possibilité de faire la séparation des sources régionales et locales du signal géophysique :

Les anomalies locales près de la surface sont caractérisées par des fréquences élevées qui sont atténuées par un prolongement vers le haut.
Les anomalies régionales et en profondeur sont caractérisées par des fréquences basses et sont atténuées par un prolongement vers le bas.

Nous allons démontrer cette application avec un exemple concret.

Distribution des propriétés pétrophysiques

Soit le modèle 3D de distribution de densité présenté à la Figure 2. Sous le mort-terrain d’une épaisseur de 20 m, la stratigraphie géologique est composée d’une série de trois groupes : un shale, un grès, et une roche volcanique (de gauche à droite). Une cavité de forme sphérique s’est formée dans le grès près de la base du mort terrain et on tente de caractériser sa géométrie avec la méthode de gravimétrie.

svg

Figure 2. Schéma d’un modèle 3D de densité. Une vue en section verticale et une vue en section horizontale est offerte. La stratigraphie géologique est une source de contraste de densité régionale et la cavité est une source locale.

Les propriétés pétrophysiques des sources sont résumés dans le tableau suivant :

Source	$\rho$ (kg/m$^3$)
Cavité	0
Mort-terrain	1000
Shale	2000
Grès	2400
Roche volcanique	2800

Modélisation directe de la gravité

L’accélération gravitationnelle mesurée à la surface terrestre ($z=0$ m) pour toutes les stations de mesures au-dessus du modèle géologique (sur le plan $xy$) a été calculée et est présentée à la Figure 3. Il s’agit de l’anomalie de Bouguer (on suppose que toutes les corrections ont été apportées).

svg

Figure 3. Anomalie de Bouguer mesurée au-dessus du modèle géologique.

La réponse gravimétrique du modèle géologique est relativement diffuse. On remarque d’abord qu’il existe un gradient dans l’orientation est-ouest du signal. Ce gradient est causé par les variations régionales dans la lithologie : la gravité est maximale au-dessus des roches volcaniques et diminue progressivement pour devenir minimale au-dessus du shale.

On aperçoit aussi un creux de forme circulaire au centre de la carte du signal gravimétrique, à $(x,y)=(0, 0)$. D’après nos connaissances a priori du modèle, il s’agit probablement de la réponse de la cavité. Cependant, cette anomalie est relativement mal définie par rapport au signal causé par la géologie régionale. Il serait souhaitable d’effectuer une séparation régionale–résiduelle pour mettre en évidence la réponse de la cavité.

Prolongement et séparation

Le protocole à suivre pour faire la séparation régionale–résiduelle est le suivant :

Transformer l’anomalie de Bouguer vers le domaine spectral
Multiplier le spectre par le filtre de prolongement vers le haut
Faire la transformation inverse du spectre pour reconstruire l’anomalie régionale
Soustraire l’anomalie régionale de l’anomalie de Bouguer pour obtenir l’anomalie résiduelle

La fonction Python suivante fait exactement ça. Prenez le temps de bien comprendre chaque étape du processus.

def separer_reg_res(gz, h, step):
    """Calcule le prolongement du champ `gz` à une hauteur `h`.
    Args:
        gz (2d array): Anomalie de Bouguer (m x n).
        h (int): Hauteur du prolongement en mètres.
        step (tuple): Pas (x, y) entre les stations en mètres.
    Returns:
        tuple: La régionale et la résiduelle.   
    """
    gz_hat = np.fft.fft2(gz)  # transformation Fourier
    k_x = 2 * np.pi * np.fft.fftfreq(n=len(gz), d=step[1])  # calcul kx
    k_y = 2 * np.pi * np.fft.fftfreq(n=len(gz), d=step[0])  # calcul ky
    KXX, KYY = np.meshgrid(k_x, k_y)  # mettre kx et ky en grille
    gz_hat *=  np.exp(-h * np.sqrt(KXX*KXX + KYY*KYY))  # filtrage
    regionale = np.fft.ifft2(gz_hat).real  # transformation inverse
    residuelle = gz - regionale  # séparation rés-rég
    return regionale, residuelle

La Figure 4 montre l’anomalie de Bouguer suite à un prolongement vers le haut de 100 m et le résultat de la séparation régionale–résiduelle.

svg

Figure 4. Séparation des anomalies régionale (A) et résiduelle (B) calculées à partir d’un prolongement vers le haut de 100 m de l’anomalie de Bouguer.

Wow! Comme prévu, le prolongement vers le haut de 100 m a atténué la réponse de la cavité pour mettre en évidence seulement les variations dans la stratigraphie (Figure 4A). La soustraction de l’anomalie de Bouguer et de l’anomalie régionale a donné l’anomalie résiduelle (Figure 4B), qui a correctement isolé la réponse gravimétrique de la cavité.

⚠️ Notez que l’anomalie résiduelle est environ 100 fois plus faible que la régionale.

Autres ressources

J’ai utilisé la librairie Python SimPEG pour créer le modèle pétrophysique en 3D et pour calculer sa réponse gravimétrique. Son utilisation est un peu trop avancée pour les travaux pratiques dans ce cours, mais c’est très utile pour faire des démonstrations.

Références

Blakely, R. (1995). Transformations. Dans Potential Theory in Gravity and Magnetic Applications (pp. 311-358). Cambridge: Cambridge University Press. doi:10.1017/CBO9780511549816.013.