Nouveaux algorithmes efficaces de modélisation 2D et 3D

Ce travail de thèse, réalisé dans le cadre d’une convention CIFRE, est le fruit d’une collaboration entre l’équipe de Géophysique du centre de Géosciences de Mines ParisTech et Magnitude, société spécialisée dans l’écoute et la surveillance de la microsismicité induite. Un séisme induit est un séisme déclenché directement ou indirectement par des activités humaines qui ont modifié le champ des contraintes géologiques et provoqué la rupture. Les séismes induits peuvent résulter de l’injection ou de l’extraction de fluides (gaz, pétrole, eau, C02,. . . ), de l’exploitation minière, de la construction de barrages, d’explosions nucléaires, de tirs en carrière, etc…

La surveillance de la sismicité induite s’est développée pendant les vingt dernières années comme discipline majeure dans le monitoring de risques liés à l’exploitation du sous sol. Elle consiste à écouter et localiser les microséismes dus à la fracturation du milieu. Les microséismes se caractérisent par une magnitude très faible (magnitude inférieure à 2). Dans le cadre de l’- exploitation de réservoirs, ils sont enregistrés par des sismomètres placés le plus souvent dans des puits pour être proches des séismes et s’affranchir du bruit de la surface mais un autre type de réseaux dit de surface est en train de se développer. L’écoute microsismique est utilisée généralement pour la surveillance des activités d’exploitation des hydrocarbures, des réservoirs géothermiques et des mines mais elle se développe aussi actuellement dans le cadre du stockage du CO2. Ce domaine a donc un très fort intérêt industriel mais tout développement est potentiellement applicable à plus grande échelle, dans un cadre plus académique de surveillance de la sismicité naturelle, comme celle enregistrée sur les volcans.

Le monitoring de la sismicité induite consiste à localiser les séismes dans la structure imagée afin de contraindre la géométrie des fractures. Lorsque c’est possible, les mécanismes au foyer sont aussi estimés pour contraindre le mode de rupture de la fracture. Le monitoring de la sismicité induite est actuellement la seule méthode permettant de suivre l’évolution de la fracturation au sein d’un réservoir. Bien que les sites soient de plus en plus instrumentés, la relation entre la fracturation et la sismicité n’est pas bien comprise. Afin de progresser dans la compréhension des phénomènes mis en jeu, il est primordial de localiser les événements sismiques et d’estimer les mécanismes au foyer de manière fiable.

La localisation d’un événement consiste à résoudre un problème inverse. Connaissant les temps d’arrivées des ondes P et S à un ensemble de capteurs, on cherche à déterminer la position de l’hypocentre. Si, dans les années 70, des méthodes linéarisées itératives étaient utilisées ([Lee & Lahr 1975], [Lahr 1999]), grâce à l’explosion des capacités informatiques, des méthodes non linéaires de type « grid search » sont maintenant préférées ([Tarantola & Valette 1982], [Lomax & Virieux 2000] et [Bardainne & Gaucher 2010]). La meilleure localisation correspond au point de l’espace et du temps qui minimise l’écart entre les temps observés aux capteurs et les temps synthétiques calculés pour un modèle de vitesse donné. Dans le contexte du monitoring microsismique, les capteurs sont le plus souvent alignés verticalement dans un puits, ce qui laisse une incertitude sur l’azimut du séisme (angle entre le Nord et la direction séisme-puits). Afin de lever cette incertitude azimutale, la polarisation de l’onde est mesurée au capteur. C’est dans ce but que les capteurs 3 composantes, pourtant plus coûteux que les capteurs mono composante, sont déployés dans les puits. La polarisation est alors utilisée soit directement pour donner l’azimut du séisme, soit incluse dans l’inversion, où l’écart entre les polarisations synthétiques et les polarisations observées sera minimisé conjointement aux temps (E. Gaucher, brevet, 2000). Le processus de localisation nécessite donc de modéliser les temps de trajet et les polarisations synthétiques entre une source et des récepteurs.

La détermination des mécanismes au foyer consiste également à résoudre un problème inverse. La méthode la plus basique consiste à utiliser la polarité du premier mouvement de l’onde P, qui porte une information sur le signe du lobe du diagramme de rayonnement d’où émerge le rai joignant le séisme au récepteur. Cette approche nécessite d’avoir préalablement localisé le séisme puis simplement de savoir modéliser l’angle (au niveau du séisme) que fait le rai joignant la source au récepteur. Un grand nombre d’observations de polarité doivent cependant être disponibles afin d’échantillonner suffisamment les différents lobes du diagramme de rayonnement pour contraindre la solution. Le monitoring étant souvent réalisé depuis un unique puits, cette méthode classique ne peut souvent pas être utilisée. L’amplitude du premier mouvement apporte alors une information supplémentaire ; en effet elle porte la polarité mais elle est aussi liée à la zone d’émergence du rai sur le lobe du diagramme de rayonnement. On peut ainsi réduire le nombre d’observations nécessaires à l’inversion du mécanisme au foyer mais cette approche nécessite alors de modéliser les amplitudes de la première arrivée.

Temps de premières arrivées : état de l’art 

Une des méthodes les plus répandues permettant le calcul des temps de trajet est la résolution de l’équation Eikonal par la méthode des différences finies. Plusieurs solveurs de l’Eikonal par différences finies de premier ordre ont été développés ([Vidale 1988, Vidale 1990], [Podvin & Lecomte 1991], [VanTrier & Symes 1991], [Hole & Zelt 1995] [Kim & R. Cook 1999], [Afnimar & Koketsu 2000], [Kim 2002], [Qian & Symes 2002], [Zhang et al. 2005], [Fomel et al. 2009]). D’autres schémas de différences finies d’ordres supérieurs on été proposés pour résoudre l’équation Eikonal ([Kao et al. 2004], [Zhang et al. 2006b], [Luo & Qian 2011], [Luo et al. 2012]). [Fomel et al. 2009] propose de résoudre l’équation Eikonal en perturbation au lieu des temps ce qui permet de mieux modéliser les fronts d’ondes à forte courbure ([Luo & Qian 2011], [Luo et al. 2012]).

Méthode de Vidale

Le temps de trajet de l’onde à été calculé par des méthodes variées. Pour les milieux simples comme le milieu homogène, le milieu à couches horizontales ou le milieu à gradient constant de la vitesse de propagation, des calculs analytiques peuvent donner des temps de propagation corrects.

En revanche, pour les milieux où la vitesse varie verticalement et latéralement, il est extrêmement complexe de calculer analytiquement les temps de trajet. Une des méthodes développées pour le calcul des temps de trajet est la résolution de l’équation Eikonal 2.115, qui relie les temps de propagation à la vitesse du milieu traversé. Dans cette section nous allons présenter une des premières méthodes de calcul des temps de trajet, basée sur l’équation Eikonal, développée par [Vidale 1988] et [Vidale 1990]. Cette formulation du problème se base sur la méthode des différences finies pour discrétiser les gradients des temps de trajet sur un modèle de vitesse échantillonné.

Temps de trajet en 2D

[Vidale 1988] a proposé une solution à l’équation Eikonal isotrope qui utilise une approximation des différences finies de premier ordre. Le schéma de calcul de temps de trajet est simple et efficace. La solution suppose que le champ de vitesse est divisé en une grille constituée de cellules carrées. Dans les cellules, les propriétés du milieu sont constantes et l’onde propagée est supposée plane. L’algorithme proposé peut être divisé en deux grandes parties : une partie de calcul local des temps de trajet et une partie de calcul global. La première partie consiste à calculer le temps de trajet en un point donné de la grille à partir des temps de trajet des points voisins supposés connus. La deuxième partie consiste à considérer le problème à grande échelle (échelle du modèle de vitesse). Il s’agit de propager les temps de trajets depuis la source vers tous les points du modèle en suivant un ordre précis.

Dans ce qui suit nous allons décrire les deux aspects de l’algorithme de calcul des temps de trajet en 2D tels que proposés par Vidale.

Approximation des ondes sphériques Près de la position de la source, l’approximation de l’onde plane présente des limites. En effet, au voisinage de la source, la courbure des fronts d’onde est très importante et contraste avec l’hypothèse des fronts d’onde plans faite plus haut. Vidale a développé un opérateur local de calcul des temps de trajets qui reproduit le comportement sphérique du front d’onde près de la source. Nous décrirons brièvement cet opérateur qui n’est utilisé pour le calcul des temps de trajet qu’au voisinage immédiat de la position de la source.

Considérons localement la maille décrite par la Figure 3.5, nous cherchons à connaître le temps d’arrivée TC au point C connaissant les temps TV , TH et TD des points respectifs V, H et D. Un front d’onde sphérique est défini sur un modèle 2D par trois paramètres (4 sur les modèles 3D) :

1- La coordonnée xS de la source virtuelle du front d’onde,
2- La coordonnée zS de la source virtuelle du front d’onde,
3- Le temps d’origine TS de la source virtuelle du front d’onde.

Table des matières

1 Introduction
2 Dynamique des milieux acousto-élastiques
2.1 Introduction
2.1.1 Équilibre des forces
2.1.2 Déformations
2.1.3 Contraintes
2.1.4 Relation Contraintes-Déformations : Loi de comportement
2.2 Milieu élastique
2.2.1 Équation de l’élastodynamique
2.2.2 Paramètres élastiques
2.2.3 Approximation des hautes fréquences en milieu élastique
2.2.4 Equation Eikonal : temps de trajet
2.2.5 Équation de transport : divergence géométrique
2.2.6 Conclusion préliminaire
2.3 Milieu acoustique
2.3.1 Équation de l’acoustodynamique
2.3.2 Paramètres acoustiques
2.3.3 Approximation des hautes fréquences en milieu acoustique
2.4 Bilan
3 Temps de premières arrivées : état de l’art
3.1 Méthode de Vidale
3.1.1 Temps de trajet en 2D
3.1.2 Temps de trajet en 3D
3.1.3 Discussion
3.2 Méthode de Podvin Lecomte
3.2.1 Temps de trajet en 2D
3.2.2 Temps de trajet en 3D
3.2.3 Discussion
3.3 Méthode de Fomel
3.3.1 Equation Eikonal factorisée
3.3.2 Schéma local
3.3.3 Calcul global
3.3.4 Discussion
3.4 Conclusion
4 Nouvelle méthode de résolution de l’Eikonal
4.1 Temps de trajet 2D
4.1.1 Opérateur de transmission 2D type Vidale
4.1.2 Opérateur de transmission 2D type Podvin
4.1.3 Opérateur de diffraction
4.1.4 Opérateurs de transmission 1D
4.2 Temps de trajet 3D
4.2.1 Les opérateurs de transmission 3D
4.2.2 Les opérateurs de transmission 2D
4.2.3 Les opérateurs de transmission 1D
4.2.4 Les opérateurs de diffraction
4.3 Résultats
4.3.1 Sensibilité aux dimensions de la grille
4.3.2 Comparaison sur un modèle homogène
4.3.3 Comparaison sur un modèle à couches
4.4 Conclusion
5 Calcul des angles à la source
5.1 Méthode de Vidale
5.2 Résolution de l’équation des Angles
5.2.1 Méthode directe
5.2.2 Méthode des perturbations
5.3 Résultats
5.3.1 Résultats 2D
5.3.2 Résultats 3D
5.4 Conclusion
6 Conclusion

Télécharger le rapport complet