3.3 Filtrage des signaux géophysiques
Types de filtres
En pratique, il existe plusieurs façons d’appliquer des filtres sur un signal, que ce soit au moyen d’opérations mathématiques sur un signal numérique ou avec des circuits électroniques pour un signal analogique. On distingue au moins quatre types de filtres en traitement de signal qui dépendent tous d’un ou de plusieurs paramètres reliés à une fréquence de coupure ($f_c$).
- Filtre passe-haut
- laisse passer les fréquences plus hautes que $f_c$.
- Filtre passe-bas
- laisse passer les fréquences plus basses que $f_c$.
- Filtre passe-bande
- laisse passer les fréquences entre $f_{c_1}$ et $f_{c_2}$, où $f_{c_1} \lt f_{c_2}$.
- Filtre coupe-bande
- empêche la propagation des fréquences entre $f_{c_1}$ et $f_{c_2}$, où $f_{c_1} \lt f_{c_2}$.
Des exemples d’images sur lesquelles on a appliqué des filtres passe-haut et passe-bas et les spectres de puissance des images résultantes sont montrés à la Figure 1. Notez comment le spectre de puissance de l’image est affecté par chaque opération.
Figure 1. Comparaison entre une image d’entrée et sa sortie après l’application d’un filtre passe-haut et d’un filtre passe-bas. Source du chat.
Les quatre types de filtres donnés ici sont des termes génériques qui désignent des familles d’algorithmes. Un filtre est considéré comme idéal quand il coupe complètement toutes les composantes du signal passé la fréquence de coupure. Les filtres appliqués à la Figure 1 ne sont pas idéaux : ils ont atténué et amplifié certaines fréquences sans faire une coupure nette.
Exemple : filtrage spectral
Dans cet exemple nous verrons comment filtrer un signal gravimétrique en passant par son domaine spectral.
Échantillonnage de la gravité
On procède à la collecte de mesures de l’accélération gravitationnelle le long d’un profil de 1 km pour améliorer notre connaissance de la géologie dans un secteur. Chaque station de mesure est espacée de 5 m ; la fréquence d’échantillonnage est donc $f_e = 0.2$ m$^{-1}$. Les échantillons recueillis sont présentés à la Figure 2. On a sélectionné la fréquence d’échantillonnage en tenant compte de la géométrie des anomalies géologiques attendues de façon à satisfaire le théorème d’échantillonnage.
Vous pouvez télécharger manuellement les données utilisées dans cet exemple ou les ouvrir directement à partir de mon site avec cette fonction :
import requests
from io import BytesIO
import numpy as np
def load_from_url(url):
address = BytesIO(requests.get(url).content)
return np.loadtxt(address)
url = 'https://clberube.org/glq2200/assets/data/profil-gravi-bruite.txt'
data = load_from_url(url)
x = data[:, 0] # la première colonne contient la position
gz = data[:, 1] # la deuxième la gravité
Figure 2. Échantillons de l’accélération gravitationnelle récoltés avec un espacement de 5 m sur la ligne 1.
On remarque que le signal est assez bruité parce que les valeurs de la gravité varient très rapidement pour des stations consécutives, ce qui produit un patron en dents de scie qui complique l’interprétation des anomalies principales dans le signal. On attribue ce bruit à un opérateur un peu trop pressé qui n’aurait pas pris toutes les précautions nécessaires pour stabiliser le gravimètre avant de faire chaque lecture.
Le bruit est caractérisé par des courtes longueurs d’ondes (hautes fréquences). Comme on s’intéresse surtout à des anomalies plus larges, comme celle identifiée sur la Figure 2, on va essayer de réduire le bruit en appliquant un filtre passe-bas.
Analyse du domaine spectral
Il faut d’abord calculer la transformation de Fourier discrète pour le signal échantillonné $\Delta g_z$. La transformation se calcule très facilement avec le module numpy.fft. Le résultat est un nombre complexe qui est décrit par une amplitude et un déphasage et les fréquences associées sont symétriques autour de l’origine (Figure 3).
# Calcul de la transformation directe
gz_tf = np.fft.fft(gz)
# Calcul des fréquences du spectre avec un pas de 5 m
freq = np.fft.fftfreq(len(gz), 5)
Figure 3. Spectres d’amplitude et de déphasage de la transformation de Fourier discrète du profil gravimétrique.
On remarque que la plupart des contributions au signal se font dans les fréquences inférieures à 0.025 m$^{-1}$ et qu’il n’y a aucune information dans les fréquences supérieures à la fréquence de Nyquist ($f_{N}$).
Coupure des hautes fréquences
La fréquence de coupure choisie ($f_c = 0.025$ m$^{-1}$) correspond à des longueurs d’ondes de l’ordre de 40 m. En coupant à cette fréquence, il faut comprendre qu’on va seulement retenir les anomalies dont la longueur d’onde est plus grande que 40 m. On effectue cette opération en ramenant à 0 toutes les composantes de la série de Fourier dont la fréquence est supérieure à $f_c$ :
freq_coupure = 0.025
gz_tf[abs(freq) > freq_coupure] = 0
En traçant le graphique des composantes qui restent on obtient le spectre présenté à la Figure 4.
Figure 4. Spectre d’amplitude après avoir coupé les fréquences supérieures à $f_c$.
Reconstruction du signal
Il ne reste plus qu’à reconstruire le signal à partir du spectre qui a été coupé à $f_c$ :
gz_reconstruit = np.fft.ifft(gz_tf)
Le spectre reconstruit est présenté à la Figure 5.
Figure 5. Signal gravimétrique reconstruit à partir du spectre coupé. Notez que le bruit a été éliminé.
On remarque que l’application du filtre passe-bas a été efficace pour filtrer le bruit à hautes fréquences qui était présent dans les échantillons. Comme le signal reconstruit ne contient plus aucune information dans les fréquences supérieures à 0.025 m$^{-1}$, le modèle géologique qui servira à interpréter les données devrait être discrétisé avec des blocs dont la largeur est comparable à la grille proposée dans la Figure 5. L’utilisation d’un modèle qui aurait des blocs plus étroits ne serait pas justifiée et le modèle proposé souffrirait de surinterprétation.
Exemple : filtrage spatial
Dans cet exemple nous verrons comment filtrer un signal gravimétrique dans le domaine spatial avec une opération de convolution très simple : la moyenne mobile.
Échantillonnage de la gravité
On s’intéresse maintenant à un profilage gravimétrique fait sur une deuxième ligne de mesure. Les paramètres du levé sont les mêmes que pour la ligne 1. Encore une fois, on constate que l’opérateur n’a pas été en mesure d’obtenir des données propres (Figure 6).
url = 'https://clberube.org/glq2200/assets/data/profil-gravi-bruite-2.txt'
data = load_from_url(url)
x = data[:, 0] # la première colonne contient la position
gz = data[:, 1] # la deuxième la gravité
Figure 6. Échantillons de l’accélération gravitationnelle récoltés avec un espacement de 5 m sur la ligne 2.
Opération de convolution
Un signal filtré ($s_\mathrm{f}$) avec une moyenne mobile du signal d’entrée $s$ est défini comme
\[s_\mathrm{f}[n] = \frac{1}{M} \sum_{i=0}^{M-1} s[n + i] ,\]où $M$ est la largeur du filtre. Cette opération dit :
Pour chaque échantillon $n$ du signal d’entrée : remplacer le point par la moyenne arithmétique du point $n$ et des $M-1$ points qui le suivent.
Le noyau de convolution du filtre de moyenne glissante est un vecteur de longueur $M$ dont tous les éléments valent $1/M$ :
\[[1/M, 1/M, ..., 1/M] .\]Prenons pour cet exemple $M=5$ et appliquons le filtre sur le signal avec la fonction numpy.convolve. Le résultat du filtrage est présenté à la Figure 7.
M = 5
filtre = np.ones(M)/M
gz_filtre = np.convolve(gz, filtre, mode='same')
Figure 7. Signal gravimétrique filtré par une opération de moyenne mobile avec $M=5$. Notez que le bruit à été en grande partie éliminé.
Beaucoup mieux! Le filtre de la moyenne mobile est un filtre passe-bas, car il atténue les hautes fréquences. On peut contrôler l’intensité de ce filtre en ajustant la largeur $M$ du noyau.