FR 3047338 A1 20170804 FR 1650864 A 20160203 1 MODELISATION EN PUITS UTILISANT LA PRESSION ET LA CONTRAINTE REGIONALE INVERSEES État antérieur de la technique [01] Une faille peut être considérée comme une discontinuité de surface tridimensionnelle complexe finie dans un volume de terre ou de roche. Les fractures, tels que les joints, les veines, les dykes, les couches de dissolution sous pression à stylolithes etc., peuvent être propagées intentionnellement, pour augmenter la perméabilité dans les formations telles que de schiste, dans lesquelles l'optimisation du nombre, du placement et de la taille des fractures dans la formation augmente le rendement des ressources comme le gaz de schiste. [02] Dans la mécanique des milieux continus, une contrainte peut être considérée comme une mesure des forces internes qui agissent au sein d'un volume. Une telle contrainte peut être définie comme une mesure de la force moyenne par unité de surface au niveau d'une surface dans le volume sur lequel agissent les forces internes. Les forces internes sont parfois produites entre les particules dans le volume en réaction aux forces externes appliquées sur le volume. [03] On peut comprendre l'origine et l'évolution des failles ainsi que l'historique tectonique des régions faillées en associant l'orientation de la faille, le sens des glissements, les données géologiques et géodésiques à l'état de contrainte dans la croûte terrestre. Dans les problèmes inverses conventionnels, les directions des contraintes principales distantes et un rapport de leurs grandeurs sont contraints en analysant des données de champ sur les orientations des failles et les directions des glissements à partir d'artéfacts tels que les striations sur les surfaces de faille exposées. 3047338 2 [04] Les procédés conventionnels d'inversion des contraintes, utilisant les striations et/ou le rejet mesurés sur les failles, reposent essentiellement sur les hypothèses que le champ de contraintes est uniforme dans 5 la masse rocheuse renfermant les failles (en supposant l'absence de champ de contraintes perturbées), et que le rejet sur les failles a la même direction et le même sens que la contrainte en champ éloigné résolue sur le plan de faille. Toutefois, il a été montré que les directions des 10 rejets sont affectées par : anisotropie de conformité de faille causée par une géométrie de ligne de rejet nul irrégulière, anisotropie de frottement de faille (ondulations de surface) ; hétérogénéité de la rigidité de la roche hôte ; et perturbation du champ de 15 contraintes locales essentiellement en raison d'interactions mécaniques de failles adjacentes. Les interactions mécaniques dues à une géométrie de failles complexes dans les milieux hétérogènes peuvent être prises en compte lors de l'exécution de l'inversion des 20 contraintes. Dans ce cas, la détermination des paramètres d'une telle paléocontrainte (et de la pression de fluide à l'intérieur des surfaces de faille) en présence de multiples failles interagissantes est effectuée en exécutant de nombreuses simulations, lesquelles demandent 25 une quantité formidable de temps de calcul pour ajuster les données observées. L'espace conventionnel de paramètres doit être exploré pour chaque possibilité, et pour chaque simulation, le modèle et les post-traitements sont recalculés. 30 [05] Avec les procédés conventionnels, les équations de mouvement sont rarement utilisées, et les perturbations du champ de contraintes locales par un rejet de faille sont ignorées. Le rôle mécanique joué par les failles dans la déformation tectonique n'est pas explicitement inclus dans de telles analyses. Ceci étant, 3047338 3 un traitement mécanique relativement complet est appliqué pour l'inversion de paléocontrainte conventionnelle. Toutefois, les résultats pourraient être grandement améliorés si l'on pouvait utiliser de multiples types de 5 données pour mieux contraindre l'inversion. Résumé [6] Un procédé, un support lisible par ordinateur et un système pour prédire un attribut de contraintes 10 d'un volume souterrain sont divulgués. Le procédé comporte la simulation d'un modèle de contrainte en champ éloigné linéairement indépendant, un modèle de pression de discontinuité, et un modèle de dérive de pression de discontinuité du volume souterrain. Une valeur de 15 contrainte, une dérive de pression de discontinuité au niveau d'une surface de la Terre (p. ex. à z=0), une valeur de déformation, une valeur de déplacement, ou une combinaison de celles-ci est/sont calculées pour des points de données dans le volume souterrain sur la base 20 d'une superposition du modèle de contrainte en champ éloigné linéairement indépendant, du modèle de pression de discontinuité, et desdits un ou plusieurs modèles de dérive de pression de discontinuité. L'attribut de contrainte du volume souterrain est prédit, sur la base 25 de la valeur de contrainte calculée, la dérive de pression de discontinuité calculée à la surface de la Terre, la valeur de déformation calculée, la valeur de déplacement calculée, ou leur combinaison. [7] Le résumé ci-dessus est fourni afin de 30 présenter une sélection de concepts qui sont décrits davantage ci-dessous dans la description détaillée. Ce résumé ne vise pas à identifier des caractéristiques clés ou essentielles de l'invention revendiquée, tout comme il n'est pas destiné à être utilisé comme aide à la 35 limitation de l'étendue de l'invention revendiquée. 3047338 4 Brève description des dessins [8] Les dessins annexés, qui sont incorporés ici et constituent une partie de la présente spécification, 5 illustrent des modes de réalisation des présents enseignements et avec la description, servent à expliquer les principes des présents enseignements. Dans les figures : [9] La figure 1 illustre un diagramme d'un système 10 de modélisation de contraintes et de fractures illustratif, conformément à un ou plusieurs modes de réalisation divulgués. [10] La figure 2 illustre un schéma de principe d'un environnement de calcul illustratif pour effectuer la 15 modélisation de contraintes et de fractures en utilisant le principe de la superposition, conformément à un ou plusieurs modes de réalisation divulgués. [11] La figure 3 illustre un schéma de principe d'un moteur de modélisation de contraintes et de fractures 20 illustratif, conformément à un ou plusieurs modes de réalisation divulgués. [12] Les figures 4A-4C illustrent des schémas de principe qui comparent des techniques de récupération de paléocontrainte, conformément à un ou plusieurs modes de 25 réalisation divulgués. Plus particulièrement, la figure 4A illustre une technique qui utilise des distributions de rejet sur les failles comme entrée de données, mais qui produit des résultats instables. La figure 4B illustre une méthode de Monte-Carlo (2D), mais sans 30 l'optimisation qui utilise le principe de la superposition. La figure 4C illustre le procédé décrit ici, y compris des techniques qui utilisent le principe de la superposition et qui réduisent formidablement la complexité du modèle. 3047338 5 [13] Les figures 5A-5E illustrent des diagrammes d'un procédé illustratif appliqué à des plans de fractures et de failles conjuguées utilisant des jeux de données sans informations de grandeur, conformément à un 5 ou plusieurs modes de réalisation divulgués. Plus particulièrement, la figure 5A montre une orientation de 03 relative à une fracture ouverte (joints, veines, dykes) donnée par sa normale el en 3D. La figure 5B montre la même chose que la figure 5A, mais pour une orientation de 10 03 relative à un joint donné par sa normale projetée fi (p. ex. tracé sur affleurement). La figure 5C et la figure 5D montrent la même chose que la figure 5A et la figure 5B, mais pour une stylolithe. La figure 5E montre l'orientation de 02 et al relativement à des plans de 15 failles conjuguées donnés par l'un de la normale n en 3D et de l'angle de frottement interne O. [14] Les figures 6A et 6B illustrent des diagrammes qui comparent l'inversion de contraintes illustratives pour des considérations de faille normale et de faille 20 chevauchante, conformément à un ou plusieurs modes de réalisation divulgués. Plus particulièrement, l'inclinaison des deux plans de failles conjuguées est présentée pour une configuration de faille normale (figure 6A) et une configuration de faille chevauchante 25 (figure 6B). [15] Les figures 7A et 7B illustrent des diagrammes de fonctions de coût des inversions de contraintes de la figure 6, conformément à un ou plusieurs modes de réalisation divulgués. Plus particulièrement, la figure 30 7A montre la fonction de coût de la faille normale, et la figure 7B montre la fonction de coût de la faille chevauchante. [16] La figure 8 illustre un diagramme d'un modèle de configuration illustratif montrant des points de 3047338 6 données InSAR et une surface de faille, conformément à un ou plusieurs modes de réalisation divulgués. [17] La figure 9 illustre un diagramme qui compare les franges d'une grille InSAR d'origine et d'une grille 5 InSAR récupérée par un exemple de procédé de modélisation de contraintes et de fractures utilisant le principe de la superposition, conformément à un ou plusieurs modes de réalisation divulgués. [18] La figure 10 illustre un diagramme montrant une 10 vue de dessus et une vue en perspective frontale d'un tracé de la surface de coût pour un exemple de procédé appliqué au cas InSAR des figures 8 et 9, en indiquant la solution de coût, conformément à un ou plusieurs modes de réalisation divulgués. 15 [019] Les figures 11A-11D illustrent des diagrammes montrant les résultats d'un procédé illustratif appliqué à un scénario d'horizon aplani, conformément à un ou plusieurs modes de réalisation. Plus particulièrement, la figure 11A présente une configuration de modèle montrant 20 l'horizon et la surface de faille. La figure 11B montre une comparaison du rejet incliné d'origine (gauche) et du rejet incliné récupéré (droit). La figure 11C montre une comparaison du rejet horizontal d'origine (gauche) et du rejet horizontal récupéré (droit). La figure 11D montre 25 le déplacement vertical d'origine (à gauche) par rapport à l'horizon aplani (à gauche) et le déplacement vertical récupéré (à droite) par rapport à l'horizon aplani. [20] La figure 12 illustre un organigramme d'un procédé illustratif de modélisation de contraintes et de 30 fractures utilisant le principe de superposition, conformément à un ou plusieurs modes de réalisation exposés. [21] La figure 13 illustre un organigramme d'un procédé illustratif de modélisation de contraintes et de 35 fractures utilisant le principe de la superposition et 3047338 7 une fonction de coût, conformément à un ou plusieurs modes de réalisation divulgués. [22] La figure 14 illustre un modèle montrant les orientations des dykes à partir de l'inversion de la 5 pression en utilisant un bouchon vertical, conformément à un ou plusieurs modes de réalisation divulgués. [23] La figure 15 illustre un organigramme d'un procédé de modélisation de contraintes et de fractures utilisant le principe de la superposition, conformément à 10 un ou plusieurs modes de réalisation divulgués. [24] La figure 16 illustre un système informatique pour exécuter au moins une partie des procédés montrés aux figures 12, 13 et 15, conformément à un mode de réalisation. 15 Description détaillée [25] La description détaillée suivante se réfère aux dessins annexés. Dans la mesure du possible, les mêmes numéros de référence sont utilisés dans les dessins et la 20 description suivante pour se référer à des parties identiques ou semblables. Bien que plusieurs modes de réalisation et caractéristique de la présente invention soient décrits ici, des modifications, adaptations et autres mises en oeuvre sont possibles, sans s'écarter de 25 l'esprit et de l'étendue de la présente invention. [26] La présente invention décrit la modélisation de contraintes et de fractures en utilisant le principe de la superposition. D'après diverses données d'entrée, telles que la géométrie des failles, et des jeux de 30 données et mesure de données sélectionnables et facultatifs, dont le rejet de faille, les directions du rejet incliné ou des stries de faille, les mesures de contraintes, les données de fracture, les orientations des plans de failles secondaires, des données de système 35 de positionnement global (GPS), des données de radar à 3047338 8 ouverture synthétique interférométrique (InSAR), des données géodésiques provenant de clinomètres de surface, la télémétrie laser, etc., les systèmes et procédés décrits ici peuvent rapidement générer et/ou récupérer de 5 nombreux types de résultats. Les systèmes et procédés décrits ici appliquent le principe de la superposition aux surfaces de failles de géométrie complexe en 2D et/ou 3D, et les failles sont, de nature, d'une dimension finie et non infinie ou semi-infinie. Les résultats sont 10 souvent rendus en temps réel et peuvent inclure, par exemple, des paramètres de contraintes, déformations et/ou déplacements en temps réel en réponse à une interrogation d'utilisateur ou à un paramètre actualisé ; des états de contraintes distantes de multiples 15 événements tectoniques ; une prédiction de fracturation future prévue ; la différenciation des fractures préexistantes des fractures induites ; et ainsi de suite. Les diverses données d'entrée peuvent être dérivées de données de puits de forage, d'une interprétation 20 sismique, d'observations des champs, etc. [27] Les systèmes et procédés décrits ci-dessous sont applicables à de nombreuses opérations de réservoir et de subsurfaces différentes, y compris les opérations d'exploration et de production de gaz naturel et d'autres 25 hydrocarbures, le stockage de gaz naturel, la fracturation hydraulique et la stimulation matricielle pour augmenter la production d'un réservoir, la gestion des ressources d'eau y compris le développement et la protection environnementale des aquifères et autres 30 ressources d'eau, la capture et le stockage souterrain de gaz carbonique (CO2), et ainsi de suite. [28] Le système peut appliquer une technique d'éléments de frontière tridimensionnelle (3D) utilisant le principe de la superposition qui s'applique à 35 l'élasticité linéaire de milieux d'espaces complets ou de 3047338 9 demi-espaces hétérogènes isotropes. D'après des valeurs précalculées, le système peut évaluer une fonction de coût pour générer des résultats en temps réel, tels que des paramètres de contrainte, déformation et déplacement 5 pour n'importe quelle point dans un volume souterrain quand l'utilisateur modifie la valeur de contrainte en champ éloigné et la pression de faille. Dans une mise en oeuvre, le système utilise la géométrie des failles et les données de puits de forage, y compris, p. ex. 10 l'orientation de la fracture, des données de plans de failles secondaires, et/ou des mesures de contraintes in situ par fracture hydraulique, pour récupérer un ou plusieurs événements tectoniques et pressions de failles, ou un ou plusieurs tenseurs de contraintes représentés 15 par un rapport de grandeurs principales et l'orientation et les pressions de failles associées. Le système peut utiliser de nombreux types différents de données géologiques tirées de l'interprétation sismique, des lectures de puits de forage et de l'observation des 20 champs pour fournir une variété de résultats, tels qu'une propagation de fracture prédite d'après le champ de contraintes perturbées. [029] La figure 1 montre un système de modélisation de contraintes et de fractures illustratif 100, 25 conformément à un ou plusieurs modes de réalisation divulgués. La géométrie des failles est souvent connue (et facultativement un rejet de faille mesuré et des contraintes d'inégalité imposées telles que faille normale, faille chevauchante, etc., peuvent être connus). 30 L'utilisateur peut avoir accès à des données de puits de forage (orientation des fractures, mesures de contraintes in situ, plans de failles secondaires), à des données géodésiques (p. ex. InSAR, GPS, et clinomètre), ainsi qu'à des horizons interprétés. Un moteur de modélisation 35 de contraintes et de fractures 102 peut récupérer l'état 3047338 10 de contraintes distantes et le régime tectonique d'événements tectoniques pertinents et la pression de failles, ainsi que la discontinuité de déplacement sur les failles, et estimer le déplacement et les champs de 5 déformations et contraintes perturbées n'importe où dans le système. [030] En utilisant le principe de la superposition, le système 100 ou le moteur 102 de modélisation de contraintes et de fractures peut exécuter chacune de 10 trois simulations linéairement indépendantes de modèles de tenseur de contraintes et d'un modèle de pression de faille en temps constant quelle que soit la complexité de chaque modèle sous-jacent. Chaque modèle n'a pas à être recalculé. Ainsi, comme indiqué ci-dessus, les 15 applications du système 100 peuvent inclure l'interpolation des contraintes et la modélisation des fractures, la récupération d'événements tectoniques, le contrôle de qualité de failles interprétées, le calcul en temps réel de champ de contraintes et de déplacements 20 perturbés quand l'utilisateur exécute une estimation de paramètres, une prédiction de la propagation des fractures, une différenciation des fractures préexistantes des fractures induites, et de nombreuses autres applications. 25 [031] La figure 2 montre le système 100 de la figure 1 dans le contexte d'un environnement informatique, dans lequel la modélisation des contraintes et fractures en utilisant le principe de la superposition peut être exécuté, conformément à un ou plusieurs modes de 30 réalisation divulgués. [032] Un dispositif informatique 200 met en oeuvre un composant, tel que le moteur de modélisation de contraintes et de fractures 102. Le moteur de modélisation de contraintes et de fractures 102 est 35 illustré sous forme de logiciel, mais peut être mis en 3047338 11 oeuvre sous forme de matériel ou de combinaison de matériel et d'instructions logicielles. [33] Le dispositif informatique 200 peut être couplé de manière à communiquer avec des dispositifs de 5 détection (et dispositifs de commande) dans un environnement réel, par exemple, un véritable volume souterrain 202, un réservoir 204, un bassin sédimentaire, un fond marin, etc., et des puits associés 206 pour la production d'une ressource pétrolière, la gestion d'une 10 ressource d'eau ou des services de stockage de dioxyde de carbone, et ainsi de suite. [34] Le dispositif informatique 200 peut être un ordinateur, un réseau informatique, ou un autre dispositif comportant un processeur 208, une mémoire 210, 15 un système de stockage de données 212, et un autre matériel associé tel qu'une interface réseau 214 et un lecteur de support 216 pour lire et écrire un support de mémorisation amovible 218. Le support de mémorisation amovible 218 peut être, par exemple, un disque compact 20 (CD) ; un disque versatile numérique/disque vidéo numérique (DVD) ; un lecteur flash, etc. Le support de mémorisation amovible 218 contient des instructions, lesquelles, à leur exécution par le dispositif informatique 200, amènent le dispositif informatique 200 25 à exécuter un ou plusieurs procédés décrits ici. Ainsi, le support de mémorisation amovible 218 peut inclure des instructions pour mettre en oeuvre et exécuter le moteur de modélisation de contraintes et de fractures 102. Au moins certaines parties du moteur de modélisation de 30 contraintes et de fractures 102 peuvent être mémorisées sous forme d'instructions sur une instance donnée du support de mémorisation amovible 218, de dispositif amovible ou de dispositif de stockage de données local 212, à charger dans la mémoire 210 en vue de leur 35 exécution par le processeur 208. 3047338 12 [35] Bien que le moteur de modélisation de contraintes et de fractures 102 soit illustré sous forme de programme résidant dans la mémoire 210, dans d'autres modes de réalisation, le moteur de modélisation de 5 contraintes et de fractures 102 peut être mis en oeuvre sous forme de matériel, tel qu'un circuit intégré spécifique à l'application (ASIC) ou de combinaisons de matériel et de logiciel. [36] Le dispositif informatique 200 peut recevoir 10 des données entrantes 220, telles que la géométrie des failles et de nombreux autres types de données, depuis de multiples sources, telles que mesures de puits de forage 222, observations de champs 224, et interprétation sismique 226. Le dispositif informatique 200 peut 15 recevoir de nombreux types de jeux de données 220 par l'intermédiaire de l'interface réseau 214, laquelle peut aussi recevoir des données de l'Internet 228, telles que des données GPS et des données InSAR. [37] Le dispositif informatique 200 peut calculer et 20 compiler des résultats de modélisation, des résultats de simulation, et des résultats de commande, et un contrôleur d'affichage 230 peut produire des images de modèles géologiques et des images et données de simulation sur un afficheur 232. Les images peuvent être 25 une simulation 2D ou 3D 234 des résultats de contraintes et de fractures utilisant le principe de la superposition. Le moteur de modélisation de contraintes et de fractures 102 peut aussi générer une ou plusieurs interfaces utilisateur (UI) visuelles pour l'entrée et 30 l'affichage des données. [38] Le moteur de modélisation de contraintes et de fractures 102 peut aussi générer et/ou produire des signaux de commande destinés à être utilisés par l'intermédiaire de dispositifs de commande, p. ex. un 35 équipement de forage et d'exploration, ou des injecteurs 3047338 13 et vannes de commande de puits, dans une commande réelle du réservoir 204, d'un réseau de transport et de livraison, d'une installation de surface, et ainsi de suite. 5 [039] Ainsi, le système 100 peut inclure un dispositif informatique 200 et une unité d'affichage graphique interactive 232. Dans son ensemble, le système 100 peut constituer des simulateurs, modèles et l'exemple de moteur de modélisation de contraintes et de fractures 10 102. [40] La figure 3 montre le moteur de modélisation de contraintes et de fractures 102 de manière plus détaillée, conformément à un ou plusieurs modes de réalisation divulgués. Le moteur de modélisation de 15 contraintes et de fractures 102 illustré à la figure 3 comporte un tampon de jeux de données 302 ou un accès aux jeux de données 302, un moteur d'initialisation 304, des simulateurs de modèle de contraintes 306 ou un accès aux simulateurs de modèle de contraintes 306, un 20 sélectionneur de paramètres d'optimisation 308, un moteur d'évaluation de coût 310, et un tampon ou dispositif de sortie des résultats 312. D'autres composants, ou d'autres agencements des composants, peuvent aussi être utilisés pour permettre diverses mises en oeuvre du moteur 25 de modélisation de contraintes et de fractures 102. La fonctionnalité du moteur de modélisation de contraintes et de fractures 102 est décrite ci-après. [41] La figure 4 montre divers procédés de récupération de la paléocontrainte, conformément à un ou 30 plusieurs modes de réalisation divulgués. La figure 4A montre une technique qui utilise les distributions du rejet sur les failles comme entrée de données mais qui produit des résultats instables. La technique de la figure 4A ne peut pas augmenter les distributions du 35 rejet avec d'autres formes d'entrée de données, telles 3047338 14 que les données GPS, et autres. La figure 4B montre une méthode de Monte-Carlo (2D) mais sans l'optimisation utilisant le principe de la superposition. La figure 4C montre le procédé décrit ici, comportant des techniques 5 qui utilisent le principe de la superposition et qui réduisent dramatiquement la complexité du modèle. Le procédé montré à la figure 4C peut être mis en oeuvre par le moteur de modélisation de contraintes et de fractures 102. Par exemple, dans une mise en oeuvre, le moteur 10 d'initialisation 304, par le biais des simulateurs de modèle de contraintes 306, génère trois modèles précalculés de la contrainte en champ éloigné associée à un volume souterrain 202. Pour chacun des trois modèles, le moteur d'initialisation 304 précalcule, par exemple, 15 des valeurs de déplacement, déformation et/ou contrainte. Le sélectionneur de paramètres d'optimisation 308 met itérativement à l'échelle les valeurs de déplacement, déformation et/ou contrainte pour chaque modèle superposé afin de minimiser un coût au niveau du moteur 20 d'évaluation de coût 310. Les valeurs ainsi optimisées en temps réel sont utilisées pour générer des résultats particuliers 312. [042] Dans une mise en oeuvre, les trois paramètres de contrainte en champ éloigné linéairement indépendants 25 sont : (i) l'orientation vers le nord, et (ii) & (iii) les deux grandeurs principales. Ces paramètres de contrainte en champ éloigné sont modélisés et simulés pour générer un ensemble de variables pour chacun des trois modèles simulés. Quatre variables peuvent être 30 utilisées : le déplacement sur une faille, le champ de déplacement à n'importe quel point de données ou point d'observation, un tenseur des déformations à chaque point d'observation, et la contrainte tectonique. Le sélectionneur de paramètres d'optimisation 308 35 sélectionne un alpha pour chaque simulation, c.-à-d. un 3047338 15 ensemble d'Halphas" pour les quatre modèles de contrainte simulés destinés à faire office de paramètres d'optimisation modifiables pour converger itérativement sur les valeurs de ces variables afin de minimiser une ou 5 plusieurs fonctions de coût, qui sont décrites ci-dessous. Dans une mise en oeuvre, le sélectionneur de paramètres d'optimisation 308 sélectionne les paramètres d'optimisation de manière aléatoire pour commencer la convergence des paramètres de déformation, de contrainte 10 et de déplacement mis à l'échelle vers le coût le plus bas. Quand les paramètres mis à l'échelle (sensiblement optimisés) sont évalués pour avoir le coût le plus bas, les paramètres de déplacement, déformation et/ou contrainte mis à l'échelle peuvent être appliqués pour 15 prédire un résultat tel qu'une nouvelle contrainte tectonique. [43] Comme le procédé de la figure 4C utilise des valeurs précalculées superposées à partir de leurs simulations respectives, le moteur de modélisation de 20 contraintes et de fractures 102 peut produire des résultats rapidement, même en temps réel. Comme on l'a indiqué ci-dessus, le moteur de modélisation de contraintes et de fractures 102 peut rapidement récupérer de multiples événements tectoniques responsables des 25 conditions présentes du volume souterrain 202, peut différencier une fracturation induite d'une fracturation préexistante plus rapidement que les techniques conventionnelles, peut fournir une estimation de paramètres en temps réel au fur et à mesure que 30 l'utilisateur modifie un paramètre de contrainte, et/ou peut prédire rapidement une fracturation. [44] Bien que l'inversion de paléocontrainte conventionnelle puisse appliquer un scénario entièrement mécanique, le moteur de modélisation de contraintes et de 35 fractures 102 constitue une amélioration par rapport aux 3047338 16 techniques conventionnelles par le fait qu'il utilise de multiples types de données Les jeux de données 302 à utiliser sont généralement de deux types : ceux qui fournissent des informations d'orientation (tels que 5 fractures, plans de failles secondaires avec angle de frottement interne, et striations de faille, etc.), et ceux qui fournissent des informations de grandeur (tels que rejet de faille, données GPS, données InSAR, etc.). Conventionnellement, l'inversion de paléocontrainte a été 10 calculée en utilisant les mesures des rejets sur les plans de failles. [45] La technique montrée à la figure 4A exécute une opération qui se déroule en deux étapes : (i) la résolution du tenseur de contraintes distantes initial GR 15 sur les éléments de faille qui n'ont pas de données de déplacement relatif et la résolution des déplacements relatifs inconnus (b3 à la figure 4A) ; et (ii) l'utilisation des déplacements relatifs calculés et connus pour résoudre aR (figure 4A). Un solveur itératif 20 qui cycle entre les étapes (i) et (ii) jusqu'à la convergence est utilisé. [46] La technique montrée à la figure 4B repose sur un algorithme de Monte-Carlo. Il est toutefois impossible d'utiliser cette technique car son temps de calcul est 25 long, d'une complexité de O(n2+p), où n et p sont respectivement un nombre d'éléments triangulaires modélisant les failles et le nombre de points de données. Pour une simulation donnée, une contrainte en champ éloigné aléatoire aR est choisie, et la discontinuité de 30 déplacement correspondante û sur les failles est calculée. Ensuite, comme post-traitement aux points de données, et en fonction des types de mesures, les fonctions de coûts sont calculées en utilisant le champ de déplacement, de déformation ou de contrainte. Dans ce 35 cas de figure, pour des centaines de milliers de 3047338 17 simulations, le meilleur coût (proche de zéro) est retenu comme solution. [47] En revanche, le procédé illustré à la figure 4A, qui peut être exécuté par le moteur de modélisation 5 de contraintes et de fractures 102, étend l'invention à de nombreux autres types de données, et fournit un moteur de modélisation beaucoup plus rapide 102. Par exemple, une inversion de contrainte rapide et fiable résultante est décrite ici. Les différents types de données peuvent 10 être pondérés et combinés ensemble. Le moteur de modélisation de contraintes et de fractures 102 peut rapidement récupérer les événements tectoniques et la pression de faille ainsi que la discontinuité de déplacement sur les failles en utilisant divers jeux et 15 sources de données, puis obtenir une estimation du champ de déplacement et de déformation et contrainte perturbées n'importe où dans le milieu, en utilisant des données produites par l'interprétation sismique, les puits de forage et l'observation des champs. L'application du 20 principe de superposition permet à un utilisateur d'exécuter très rapidement l'estimation des paramètres. [48] Une technique numérique d'exécution des procédés est décrite ci-après. Un tenseur distant réduit utilisé pour la simulation est décrit par la suite, suivi 25 du principe de superposition proprement dit. Une estimation de la complexité est également décrite. [49] Dans une mise en oeuvre, une formulation appliquée par le moteur de modélisation de contraintes et de fractures 102 peut être exécutée en utilisant IBEM3D, 30 un successeur de POLY3D. POLY3D est décrit par F. Maerten, P. G. Resor, D. D. Pollard, et L. Maerten, Inverting for slip on three-dimensional fault surfaces using angular dislocations, Bulletin of the Seismological Society of America, 95 :1654-1665, 2005, et par A. L. 35 Thomas, Poly3D : a three-dimensional, polygonal element, 3047338 18 displacement discontinuity boundary element computer program with applications to fractures, faults, and cavities in the earth's crust, Master's thesis, Stanford University, 1995. IBEM3D est un code d'éléments de 5 frontière basé sur la solution analytique d'une dislocation angulaire dans un espace entier ou un demi-espace élastique homogène ou non homogène. Un solveur itératif peut être employé pour une plus grande vitesse ainsi que pour une parallélisation sur les architectures 10 à coeurs multiples. (Voir, par exemple, F. Maerten, L. Maerten, et M. Cooke, Solving 3d boundary element problems using constrained iterative approach, Computational Geosciences, 2009). Toutefois les contraintes d'inégalités ne peuvent pas être utilisées 15 car elles sont non linéaires et le principe de la superposition n'est pas applicable. Dans le code sélectionné, les failles sont représentées par des surfaces triangulées avec un déplacement discontinu. Les surfaces de failles tridimensionnelles peuvent de plus 20 prêt se rapprocher de surfaces ondulées et de lignes de rejet nul ondulées sans introduire de chevauchements ou d'espaces. [050] Des conditions de frontières mixtes peuvent être prescrites, et lorsque des conditions de frontières 25 à tractions sont spécifiées, le moteur d'initialisation 304 résout les composantes de Burger inconnues. Après la résolution du système, il est possible de calculer n'importe où, dans l'espace entier ou le demi-espace, le déplacement, la déformation ou la contrainte aux points 30 d'observation, dans le cadre d'un post-traitement. De manière spécifique, le champ de contrainte à n'importe quel point d'observation est donné par le champ de contrainte perturbée dû aux failles à glissement et à la contribution de la contrainte distante. En conséquence, 35 l'obtention du champ de contrainte perturbée dû au rejets 3047338 19 de failles ne suffit pas. De plus, l'estimation du rejet de faille à partir de l'interprétation sismique est donnée dans le sens du pendage. Rien n'est connu dans la direction azimutale, et un scénario mécanique complet est 5 utilisé pour récupérer les composantes inconnues du vecteur de glissement qui aura un effet sur le champ de contrainte perturbée. Le changement de la contrainte en champ éloignée imposée (orientation et/ou grandeurs relatives) modifie la distribution de rejet de faille et 10 en conséquence le champ de contrainte perturbée. En général, un code tel que IBEM3D convient bien aux calculs des vecteurs de déplacement complet sur les failles, et a été intensivement optimisé en utilisant une technique à matrice H. L'inconnu aux fins de la modélisation reste 15 l'estimation de la contrainte en champ éloigné qui doit être imposée comme condition de frontières. [051] Dans un mode de réalisation qui peut être mis en oeuvre par le moteur de modélisation de contraintes et de fractures 102, un modèle composé de multiples surfaces 20 de faille est soumis à un tenseur de contraintes en champ éloigné oR défini dans le système de coordonnées global par l'équation (1) : aR. a11 a12 a22 a23 ( 1 ) a33 En supposant une contrainte en champ éloigné sub25 horizontale (mais la présente méthodologie ne se limite pas à ce cas), l'équation (1) est simplifiée en l'équation (2) : a11 a12 (FR = a22 0 (2) a33 Comme l'addition d'une contrainte hydrostatique ne change 30 pas aR, le tenseur de contraintes en champ éloigné aR peut être écrit comme dans l'équation (3) : 3047338 20 trR = a22 a33 0 = 11 a12] à22 (3) a11 - a33 a12 r 0 5 [52] en champ savoir En utilisant En 10 [53] que utilisant l'équation éloigné o-R = fà Le où, défini conséquence, une définition d'une contrainte = 40129 avec trois inconnues est obtenue, à comme à tenseur de contraintes en champ éloigné, tel [0-1 01 à l'équation (3) peut être calculé en az à. - deux paramètres au lieu des trois fà 22) (112} la décomposition spectrale du oR réduit, (4) peut être obtenue : (4) l'équation (5) : (5) le tenseur de contraintes en champ éloigné aR est la matrice de valeurs principales, et à l'équation (6) : 15 n= [cos° -sin 01 no (6) sin 0 cos6 -I est la matrice de rotation autour de l'axe z global (puisqu'un tenseur de contraintes sub-horizontales est supposé). [54] En écrivant, à l'équation (7) : 20 Q2 = koi (7) l'équation (4) prend alors la forme de l'équation (8) : 6R(6, = cr1R79' [0 1 0 k] Ro (8) [55] En omettant le paramètre de mise à l'échelle ci dû à la propriété 1 discutée ci-dessous (quand a1=5 dans 25 la propriété 1), aR peut être exprimé en fonction de deux paramètres 5 et k, comme à l'équation (9) : i T r o 0-R(9,k) = Re 1_0 k Re (9) [56] Ces deux paramètres sont naturellement bornés par l'équation (10) : f-n-/2 19 n-/2 30(10) t. -10 10. Pour k=1, un tenseur de contraintes hydrostatiques est trouvé, lequel n'a aucun effet sur le modèle. Par ailleurs, un tenseur de 5 contraintes en champ éloigné lithostatiques (qui est donc fonction de la profondeur z) n'invalide pas la technique présentée, et l'équation (9) prend la forme de l'équation (11) : 6R(9,k, z) = zR70. [01 k°] R0 (11) 10 qui dépend linéairement de z. La définition de tenseur simplifiée donnée par l'équation (9) est utilisée dans les sections suivantes pour déterminer 6, k), ou de manière équivalente ..-- fà 11,à22,à12) conformément aux mesures. [057] Même quand un espace paramétrique 15 tridimensionnel est utilisé pour la simulation de Monte- Carlo en utilisant (0, k, p), où p est la pression de faille, trois composantes sont toujours utilisées pour la contrainte en champ éloigné, spécifiées par les paramètres (al, 0(2, a3, u4). Les conversions sont données 20 par l'expression (12) : (0, k) -) (cfoo, am, an) -' (ai, «2, «3) (12) où à l'équation (13) : o-00= cos2 0 + k sin2 9 { 601 = (k - 1)cos 0 sin 0 (13) un = cos2 0 + k sin2 0 et (al, a2, a3) sont donnés par l'équation (20) ci-après. 25 [058] Principe de superposition [059] Le moteur de modélisation de contraintes et de fractures 102 utilise le principe de superposition, un principe bien connu dans la physique de l'élasticité linéaire, pour récupérer le déplacement, la déformation 30 et la contrainte à n'importe quel point d'observation P en utilisant les valeurs spécifiques précalculées tirées des simulations linéairement indépendantes. Le principe de superposition stipule qu'une valeur donnée f peut être 3047338 22 déterminée par une combinaison linaire de solutions spécifiques. [060] Dans le moteur de modélisation de contraintes et de fractures 102, la récupération d'une contrainte en 5 champ éloigné implique la récupération des trois paramètres (al, a2, a3). En conséquence, le nombre de solutions linéairement indépendantes utilisé est trois. Autrement dit, dans l'équation (14) : f = F (a) = F(cr16(1) + a2o-(2) + a3a(3)) criF(6(1)) + a2F (cr(2)) + a3F (o-(3)) (14) =a1f1ry f f -2J2 - lY 3J3 10 où (al, a2, a3) sont des nombres réels et am (pour i=1 à 3) sont trois tenseurs de contraintes distantes linéairement indépendants. Si F est sélectionnée en tant que fonctions de Green de déformation, contrainte ou déplacement, les valeurs résultantes E, o, et u, à P 15 peuvent alors être exprimées comme une combinaison de trois solutions spécifiques, comme cela est montré ci-dessous. Ainsi, les champs de déformation, contrainte et déplacement d'une charge tectonique sont une combinaison linéaire des trois solutions spécifiques, et sont donnés 20 par l'équation (15) : (1) (2) (3) { Ep = alEp + a2Ep ± a3Ep (1) (2) (3) Op = aifFp + a2Cfp + a3Cfp (15) ( ( ( Up --= ailip1) ± a2Up2) ± a3Up3) [061] De la même façon, l'utilisation de (al, a2, a3) permet de récupérer les discontinuités de déplacement sur les failles, comme à l'équation (16) : 25 be = alibe(1) + a2be(2) + a3b,(3) (16) et n'importe quelle contrainte en champ éloigné est également donnée sous forme de combinaison des trois paramètres, comme à l'équation (17) : (i) _i_ (2) _i_ n, '(3) aR = cricri? -r- a2 aR -I- '.3 ''R 30 [062] Estimation de la complexité (17) 3047338 23 [63] Le modèle complet est souvent recalculé pour changer GR afin de déterminer les discontinuités de déplacement inconnues correspondantes. Ainsi, à n'importe quel point d'observation P, la contrainte est déterminée 5 comme une superposition de la contrainte en champ éloigné oR et du champ de contrainte perturbée dû aux éléments de glissement. [64] Pour un modèle composé de n éléments triangulaires discontinus, le calcul de l'état de 10 contrainte au point P implique de résoudre les discontinuités de déplacement inconnues sur les éléments triangulaires (dont la complexité est 0(n2)), puis d'exécuter approximativement 350n multiplications en utilisant le procédé standard. En revanche, en utilisant 15 le principe de superposition, les discontinuités de déplacement inconnues sur les éléments triangulaires n'ont pas à être recalculées et un nombre bien inférieur de multiplications (p. ex. 18) est exécuté par le moteur de modélisation de contraintes et de fractures 102. La 20 complexité est indépendante du nombre d'éléments triangulaires dans le modèle et est constante dans le temps. [65] Certaines applications directes des procédés vont maintenant être décrites, telles que l'évaluation en 25 temps réel de la déformation et du champ de contrainte perturbée quand un utilisateur varie un paramètre de contrainte en champ éloigné. L'estimation de la paléocontrainte en utilisant différents jeux de données 302 est également présentée davantage ci-dessous, tout 30 comme un procédé de récupération de multiples phases tectoniques, et une description de la manière dont l'exemple de procédé peut être utilisé pour le contrôle de qualité durant l'interprétation d'une faille. [66] Calcul en temps réel 3047338 24 [67] Avant de décrire le procédé d'inversion de paléocontrainte, nous décrivons tout d'abord un autre procédé pour calculer en temps réel la discontinuité de déplacement sur les failles, et les champs de 5 déplacement, déformation et contrainte aux points d'observation alors que l'orientation et/ou la grandeur de la contrainte en champ éloigné est modifiée. [68] Si la contrainte tectonique aR est donnée et trois solutions indépendantes sont connues, il existe un 10 triplet unique (ai, a2, 0(3) pour lequel l'équation (17) est satisfaite, et les équations (15) et (16) sont applicables. [69] Sous forme matricielle, l'équation (17) est écrite dans le format montré à l'équation (18) : (1) (2) (3) 600 600 600 a1 o-eo u(i) o-(2) (3) a = R 15 (18) (roi 2 - 19-01 (1)(2) (3) a3 011 0-11 0-11 °i1 ou, sous forme compacte, comme à l'équation (19) : Ao_a = o-R (19) [70] Comme les trois solutions particulières &il sont linéairement indépendantes, le système peut être 20 inversé, ce qui donne l'équation (20) : a = AQ16R (20) [71] A l'équation (20), AJ1- est précalculé par le moteur d'initialisation 304. Etant donné une contrainte distante sélectionnée par l'utilisateur, aR, le moteur de 25 modélisation de contraintes et de fractures 102 récupère les trois paramètres (al, a2, a3), puis le rejet de faille et les champs de déplacement, déformation et contrainte sont calculés en temps réel en utilisant les équations (16) et (15), respectivement. Pour cela, les trois 30 solutions particulières du déplacement, de la déformation et de la contrainte sont mémorisées à l'initialisation à chaque point d'observation, ainsi que la discontinuité de 3047338 25 déplacement sur les failles. Dans une mise en oeuvre, l'exemple de moteur de modélisation de contraintes et de fractures 102 permet à l'utilisateur de modifier l'orientation et la grandeur de GR, et d'afficher 5 interactivement le champ de déformation et de contrainte perturbée associé. [72] Inversion de la paléocontrainte au moyen des jeux de données [73] Comme on l'a vu plus haut, les principales 10 inconnues dans une modélisation directe de l'estimation de la distribution de rejet sur les failles, et en conséquence du champ de contrainte perturbée associé, sont l'orientation et les grandeurs relatives de la contrainte en champ éloigné GR. 15 [074] Si les mesures de champ sont connues à certains points d'observation donnés (p. ex. déplacement, déformation et/ou contrainte, orientation des fractures, plans de failles secondaires qui se sont formés au voisinage de failles principales, etc.), il est alors 20 possible de récupérer le triplet (al, a2, a3), et en conséquence de récupérer également la contrainte tectonique GR et le régime tectonique correspondant. La section suivante décrit la méthode de résolution et les fonctions de coût utilisées pour minimiser le coût de 25 différents types de jeux de données 302. [75] Méthode de résolution [76] L'application de la technique de Monte-Carlo permet de trouver les paramètres (al, 0(2, a3), ce qui minimise les fonctions de coût pour trois contraintes en 30 champ éloigné indépendantes données (voir l'équation 15). Toutefois, même si (al, a2, a3) impliquent un espace paramétrique tridimensionnel, cet espace peut être réduit à deux dimensions (à savoir, aux paramètres 0 et k), la conversion étant donnée par l'équation (20) et (6,k). 35 -*.(000, coi, Gli) -> (al, 0(2, a3) , où à l'équation (21) : 3047338 26 1 600 = cos2 0 + k sin2 0 obi = (k -- 1)cos 0 sin 0 (21) = cos2 0 + k sin2 0 an (pour une description détaillée, voir aussi la figure 4-3 et l'algorithme 1). En conséquence, le procédé de recherche (p. ex. la recherche de paramètres optimisés) 5 est accéléré en réduisant l'espace paramétrique par une dimension. [077] Un procédé d'échantillonnage simple peut être mis en oeuvre en considérant un domaine rectangulaire bidimensionnel dont les axes correspondent à 0 et k. Le 10 domaine 2D est échantillonné de manière aléatoire avec np points, et la fonction de coût associée (qui est définie dans les sections suivantes) est utilisée pour déterminer le point de coût minimum. Un affinement est alors créé autour du point sélectionné et la procédure est répétée 15 avec un domaine plus petit. L'algorithme (1) illustre une version simplifiée de la procédure, pour laquelle aucun affinement n'existe. L'exemple de procédé d'échantillonnage présenté ici peut être grandement optimisé par diverses techniques. 20 ALGORITHME (1) : Estimation de la paléocontrainte 25 Donnée d'entrée : géométrie des failles Donnée d'entrée : jeux de données Donnée produite : (TRI/ la paléocontrainte estimée Initialiser : calculer 3 simulations en utilisant 3 (i) (1 ((Focs, cfoi, --> (ai, az, a3) (12') 5 [087] En outre les équations (13) et (21) sont remplacées par : I600 = sin f3 - cos ie cosi 0 an = cos fi' sin 0 cos 0 (13' ) 600 = sin 13 - cos /3 sine 0 [88] f3 étant un angle défini par : /3= tan-1 (R' R-, 1) pour un régime de faille normale p = tan-1(R' - 1) pour un régime de faille à décrochement f3 = tan-1 (3 2- R') pour un régime de faille chevauchante [89] En outre, l'algorithme (1) peut être modifié en 10 l'algorithme (1') suivant : ALGORITHME (1') : Estimation de la paléocontrainte 15 Donnée d' entrée : géométrie des failles Donnée d' entrée : jeux de données Donnée produite : uR / / la paléocontrainte estimée 20 Initialiser : calculer 3 simulations en utilisant 3 (FM (1 en 3D. La figure 5B est identique 15 à la figure 5A, hormis que l'orientation de os est donnée par rapport à un joint donné par sa normale projetée h9. (p. ex. sur affleurement). La figure 5C et la figure 5D montrent la même chose que les figures 5A et 5B mais pour une stylolithe. La figure 5E montre l'orientation de 02 et 20 ai par rapport aux plans de failles conjuguées donnés par l'un de la normale e en 3D et de l'angle de frottement interne O. Le but est de déterminer le meilleur ajustement de la contrainte en champ éloigné aR, et en conséquence les paramètres (ai, CX2, a3), étant donné 25 certaines orientations de plans de fractures ouvertes pour lesquels les normales coïncident avec les directions de la contrainte la moins compressive a3P à P, ou de manière équivalente pour lesquels le plan de la fracture contient la contrainte la plus compressive (al), comme à 30 la figure 5A et la figure 5B. [96] En modifiant (al, a2, a3), l'état de contrainte à n'importe quel point d'observation P peut être calculé rapidement en utilisant the trois modèles précalculés. La fonction de coût à minimiser est donnée à l'équation 35 (24) : 3047338 32 =_'e.ffp)2±(er,2).ffp)2,1,2 ffrac(al, az, a3) 2mze où "." est le produit scalaire, n est la normale à un plan de fracture, et m est le nombre de points d'observation. La minimisation d'une fonction des trois 5 paramètres est exprimée par l'équation (25) : "frac = Minai,a2,a3t ffrac (al) a2, a3)) (25) [97] De la même façon, pour les couches de dissolution sous pression et stylolithes, la fonction de coût est définie comme à l'équation (22) en utilisant la 10 contrainte la moins compressive 03 comme à l'équation (26) (voir la figure 5-3 et la figure 5-4) : fsty/(ai, a2, a3) = 2m Ep ((of eP)2 + e/12)1/2 (26) [98] Utilisation des plans de failles secondaires [99] L'orientation des plans de failles secondaires 15 qui se développent au voisinage de failles actives plus importantes peut être estimée en utilisant un critère de défaillance de Coulomb, défini par l'équation (27) : tan(20) = 1/y (27) où e est l'angle des plans de défaillances par rapport à 20 la contrainte compressive principale maximum 01 et p est le coefficient de frottement interne. Deux plans de défaillances conjuguées s'intersectent le long de 02 et l'orientation de la faille est influencée par l'orientation des contraintes principales et la valeur du 25 frottement. [0100] La fonction de coût est en conséquence définie par l'équation (28) : 2 112 [0-2 np) 2 efiP + (Isin-lr.)I-u) ffau/t(al, a2,, C(3) = 2m 2,p (28) où 01 est la direction de la contrainte la plus 30 compressive et 02 est la direction de la contrainte principale intermédiaire. Le premier terme droit de l'équation (26) maintient une orthogonalité entre le 02 calculé et la normale au plan de faille, tandis que le (24) 3047338 33 second terme garantit que l'angle entre le o calculé et le plan de faille est proche de 0 (voir la figure 5-5). Exemple : Faille normale et faille chevauchante [0101] La figure 6 montre un exemple synthétique 5 utilisant une faille plane inclinée comme l'un de deux plans de failles conjuguées sélectionné aléatoirement, conformément à un ou plusieurs modes de réalisation divulgués. L'inclinaison des deux plans de failles conjuguées est présentée pour une configuration de faille 10 normale (la figure 6A) et une configuration de faille chevauchante (figure 6B). L'azimut de rejet et l'angle de rejet de chaque plan de failles conjuguées sont utilisés pour effectuer l'inversion et l'angle de frottement interne est 0=30. La faille active principale est 15 représentée par le plan rectangulaire incliné 602. [0102] Initialement, le modèle est contraint par une contrainte en champ éloigné à certains points d'observation 604, où les deux plans conjugués sont calculés en utilisant an angle de frottement interne de 20 30 degrés. Ensuite, pour chaque point d'observation 604, l'un des plans de failles conjuguées est choisi aléatoirement et utilisé comme données d'entrée pour l'inversion de contrainte. [0103] La figure 7 montre la fonction de coût de 25 l'exemple synthétique de la figure 6, conformément à un ou plusieurs modes de réalisation. La figure 7A montre la fonction de coût de la faille normale, et la figure 7B montre la fonction de coût de la faille chevauchante. Dans les deux cas, les tenseurs de contraintes 30 régionales, déplacement sur la faille, et plans de failles conjuguées prédits récupérés correspondent bien au modèle synthétique initial. [0104] Utilisation des striations de faille [0105] Dans le cas des striations de failles, la 35 fonction de coût est définie comme à l'équation (29) : 3047338 34 Ep(1 - (1 iie- c e J2)112 \112 (29) fstri(cti,a2,a3)=- où è4 et dé représentent respectivement le vecteur de gliCsement normalisé à partir d'une simulation et le vecteur de glissement mesuré. 5 [0106] Jeux de données contenant des informations de grandeur [0107] La grandeur des déplacements peut être utilisée pour déterminer l'orientation des contraintes et la grandeur du tenseur de contraintes distantes, plutôt que 10 simplement le rapport de contraintes principal. [0108] Pour cela, la procédure est semblable à celle décrite précédemment. Toutefois, étant donné les équations (15) et (16), il est évident qu'il existe un paramètre 5 pour lequel la discontinuité de déplacement 15 calculée sur les failles et les champs de déplacement, déformation et contrainte aux points d'observation sont gradués linéairement avec la contrainte en champ éloigné imposée. Autrement dit, comme à l'équation (30) : {80-R Sbe 60-R 61.1p &JR SEp 60-R Sap 20 [0109] Ceci conduit à la propriété suivante : [0110] Propriété 1 : La mise à l'échelle de la contrainte en champ éloigné par S E 9î met à l'échelle la discontinuité de déplacement sur les failles ainsi que les champs de déplacement, déformation et contrainte aux 25 points d'observation par 5. [0111] En utilisant cette propriété, les mesures aux points de données sont globalement normalisées avant tout calcul et le facteur de mise à l'échelle est noté (les simulations sont aussi normalisées, mais le facteur 30 d'échelle est sans importance). Après la résolution du système, les champs de contrainte en champ éloigné, (30) 3047338 35 déplacement et contrainte récupérés sont réduits par un facteur de 5J1. [0112] Utilisation des données GPS [0113] Dans le cas d'un jeu de données GPS, la 5 fonction de coût est définie à l'équation (31) : 2 ( 2 1/2 fgps(ai, a2, a3) = -Ep [(1 + 1 - (el; - 7;121 21n 1114 lillun) 11411 (31) où ie est l'élévation mesurée normalisée globalement changée au point P depuis l'horizon, et ûP est le 10 changement d'élévation calculé normalisé globalement pour un ensemble de paramètres (al, a2, a3) donné. Le premier terme droit de l'équation (31) représente une minimisation de l'angle entre les deux vecteurs de déplacement, tandis que le second e représente une 15 minimisation de la différence de la norme. [0114] Utilisation des données InSAR [0115] L'utilisation d'un jeu de données InSAR donne lieu à deux possibilités. Soit les vecteurs de déplacement globaux des mesures sont calculés en 20 utilisant le déplacement u le long de la direction de la ligne de visée du satellite g, auquel cas l'équation (32) est utilisée : uP = ilinsar - u. (32) 25 et la même procédure qui est utilisée pour les jeux de données GPS (ci-dessus) est appliquée avec le e4), calculé, soit les vecteurs de déplacement calculés sont calculés le long de la ligne de visée du satellite, auquel cas l'équation (33) est utilisée : 30 uf,=ii.g (33) où "." est the produit scalaire. La fonction de coût est en conséquence donnée par l'équation (34) : 2)\1/2 finsar uc a2, a3) = -Ep(1 - (34) up 3047338 36 [0116] Les figures 8, 9 et 10 présentent un exemple synthétique en utilisant un jeu de données InSAR. Plus particulièrement, la figure 8 montre une configuration de modèle montrant les points de données InSAR 802 ainsi que 5 la surface de faille 804, La figure 9 montre une comparaison des franges de la grille InSAR d'origine 902 et de la grille InSAR récupérée 904, et la figure 10 montre un tracé de la surface de coût, en n de 0 (axe x) et k (axe y), conformément à un ou plusieurs modes de 10 réalisation divulgués. A gauche de la figure est présentée une vue de dessus 1002 du tracé, et à droite de la figure est présentée une vue en perspective frontale 1004 du tracé. La solution de coût minimisée 1006 dans chaque vue est marquée par un petit cercle blanc (1006). 15 [0117] Pour utiliser un jeu de données InSAR, un modèle direct est exécuté en utilisant un plan de faille 804 et une grille d'observation (figure 8) au niveau de la surface du demi-espace (voir la surface renfermant les points de données 802). Une direction de satellite est 20 sélectionnée, et pour chaque point d'observation 802, le déplacement le long de la ligne de visée du satellite est calculé. Ensuite, le procédé de modélisation de contraintes et de fractures décrit ici est appliqué en utilisant la seconde forme de la fonction de coût InSAR 25 donnée à l'équation (34). La figure 9 compare l'interférogramme d'origine 902 (gauche) à l'interférogramme récupéré 904 (droit). La figure 10 montre combien la surface de coût peut être complexe, même pour un modèle synthétique simple. Dans une mise en 30 oeuvre, la surface de coût a été échantillonnée avec 500 000 points de données 802 (nombre de simulations), et a demandé 18 secondes sur un ordinateur portable moyen à processeur de 2 GHz et 2 GB de mémoire vive exécutant Linux Ubuntu version 8.10, 32 bits. 35 [0118] Utilisation d'un horizon aplani 3047338 37 [0119] En utilisant le plan moyen d'un horizon sismique donné (horizon aplani), le moteur de modélisation de contraintes et de fractures 102 calcule d'abord le changement d'élévation de chaque point 5 constituant l'horizon. Ensuite la fonction de coût GPS peut être utilisée, pour laquelle la composante uz est fournie, donnant l'équation (35) : ' (up fhorizon(al, az, a3) ' (35) `uni)2)1/2 Pll Si un pré- ou post-plissement de la zone est constaté, le 10 plan moyen ne peut plus être utilisé comme proxy. En conséquence, une surface d'ajustement lisse et continue qui supprime les déformations de faille tout en maintenant les plis doit être construite. Ensuite, la même procédure que celle pour le plan moyen peut être 15 utilisée pour estimer la paléocontrainte. Dans certaines circonstances et avant la définition de la surface d'ajustement continue, les bruits haute fréquence, tels que les ondulations et bosses de l'horizon d'entrée peuvent être filtrés tout en préservant les déformations 20 naturelles. [0120] La figure 11 montre les résultats obtenus lors de l'application d'un procédé de modélisation de contraintes et de fractures à un exemple synthétique en utilisant un horizon aplani, conformément à un ou 25 plusieurs modes de réalisation divulgués. La figure 11A présente une configuration de modèle montrant l'horizon 1102 et la surface de faille 1104. La figure 11B montre une comparaison du rejet incliné 1106 (gauche) et du rejet incliné récupéré 1108 (droit). La figure 11C montre 30 une comparaison du glissement horizontal d'origine 1110 (gauche) et du glissement horizontal récupéré 1112 (droit). La figure 11D montre le déplacement vertical d'origine 1114 (gauche) par rapport à l'horizon aplani 3047338 38 (gauche) et le déplacement vertical récupéré 1116 (droit) par rapport à l'horizon aplani. [0121] Comme le montre la figure 11, une faille conformée complexe est initialement contrainte par une 5 contrainte en champ éloigné, et en conséquence glisse pour accepter la contrainte distante. A chaque point d'un plan d'observation entrecoupant la faille, le moteur de modélisation de contraintes et de fractures 102 calcule le vecteur de déplacement résultant et déforme la grille 10 comme il convient. L'inversion se produit alors en utilisant la géométrie de la faille. Après l'aplanissement de la grille déformée, le changement d'élévation de chaque point est utilisé pour contraindre l'inversion et récupérer la contrainte en champ éloigné 15 précédemment imposée ainsi que le glissement de faille et le champ de déplacement. La comparaison des rejets incliné (figure 11B) et rejet horizontal (figure 11C) d'origine et récupéré montre qu'ils concordent bien (même échelle). Une bonne concordance est aussi constatée pour 20 le champ de déplacement au niveau de la grille d'observation (figure 11D). [0122] Utilisation d'informations de rejet incliné [0123] Quand des données de rejet incliné sont utilisées, la fonction de coût est définie comme à 25 l'équation (36) : 1 (_be) )1/2 2 fds a2, a3) = -Ep(1 - (36) bem où ben' est la grandeur de rejet incliné mesurée pour un élément triangulaire e, et bec est la grandeur de rejet incliné calculée. 30 [0124] Utilisation des informations disponibles [0125] Le moteur de modélisation de contraintes et de fractures 102 peut combiner les fonctions de coût précédemment décrites pour mieux contraindre l'inversion en utilisant les données disponibles (p. ex. données 3047338 39 d'orientation de plan de faille et de fracture, données GPS, données InSAR, données d'horizon aplani, mesures de rejet incliné tirées de la réflexion sismique, striations de faille, etc.). En outre, les données peuvent être 5 pondérées différemment, et chaque donnée peut aussi prendre en charge un poids pour chaque coordonnée. [0126] Multiples événements tectoniques [0127] Dans le cas de multiples événements tectoniques, il est possible d'en récupérer les 10 principaux, p. ex. ceux dont le régime tectonique et/ou l'orientation et/ou la grandeur sont notablement différents. L'algorithme 2 ci-dessous présente une façon de déterminer différents événements à partir de l'orientation des fractures (joints, stylolithes, plans 15 de failles conjuguées) mesurée le long des puits de forage. [0128] Après avoir effectué une première simulation, un coût est attaché à chaque point d'observation pour montrer la fiabilité de la contrainte tectonique 20 récupérée par rapport aux données attachées à ce point d'observation. Un coût de zéro implique une bonne fiabilité, tandis qu'une code de un signifie une fiabilité médiocre. Voir à la figure 7 un exemple de tracé du coût. En sélectionnant des points de données qui 25 sont en dessous d'une valeur limite donnée et en exécutant une autre simulation avec ces points, il est possible d'extraire une valeur de paléocontrainte plus précise. Les points de données restants au-dessus de la valeur limite peuvent alors être utilisés pour exécuter 30 une autre simulation avec l'état de paléocontrainte pour récupérer un autre événement tectonique. Si le graphique du nouveau coût montre des disparités, l'exemple de procédé ci-dessus est répété jusqu'à l'obtention de résultats satisfaisants. Durant la détermination des 35 phases tectoniques, les points d'observation sont classés 3047338 dans leur événement tectonique respectif. La chronologie des phases tectoniques reste toutefois indéterminée. 5 ALGORITHME (2) : Détection de multiples événements tectoniques Donnée d'entrée : E comme seuil d'utilisateur dans 10 [0, 1] Donnée d'entrée : S = toutes les fractures de tous les puits Soit : T=0 quand S 0 faire 15 Simulations : calculer le coût pour chaque fracture in S Si max (coût) REVENDICATIONS1. Procédé de prédiction d'un attribut de contrainte d'un volume souterrain, comprenant : une capture de données utilisant un capteur, le capteur faisant partie d'un outil de mesure pendant le forage ou d'un outil de diagraphie durant le forage dans un puits de forage; une simulation d'un ou de plusieurs modèles de contrainte en champ éloigné linéairement indépendants, d'un ou de plusieurs modèles de pression de discontinuité, et d'un ou de plusieurs modèles de dérive de pression de discontinuité pour le volume souterrain en utilisant les données capturées; un calcul, en utilisant un processeur, d'une valeur de contrainte, d'une dérive de pression de discontinuité au niveau d'une surface de la Terre, d'une valeur de déformation, d'une valeur de déplacement, ou d'une combinaison de celles-ci pour des points de données dans le volume souterrain sur la base d'une superposition desdits un ou plusieurs modèles de contrainte en champ éloigné linéairement indépendants, desdits un ou plusieurs modèles de pression de discontinuité, et desdits un ou plusieurs modèles de dérive de pression de discontinuité ; et une prédiction de l'attribut de contrainte du volume souterrain, sur la base de la valeur de contrainte calculée, la dérive de pression de discontinuité calculée à la surface de la Terre, la valeur de déformation calculée, la valeur de déplacement calculée, ou la combinaison de celles-ci. 2. Procédé selon la revendication 1, comprenant en outre une inversion, au moyen du processeur, desdits un ou plusieurs modèles de contrainte en champ éloigné linéairement indépendants. 3. Procédé selon la revendication 1, comprenant en outre 35 l'inversion, au moyen du processeur, desdits un ou plusieurs modèles de pression de discontinuité. 58 3047338 4. Procédé selon la revendication 3, dans lequel lesdits un ou plusieurs modèles de pression de discontinuité sont inversés par rapport à une différence entre une contrainte régionale inversée principale maximum et une contrainte régionale inversée 5 principale minimum. 5. Procédé selon la revendication 3, dans lequel un ou plusieurs modèles de pression de discontinuité sont inversés sur la base au moins en partie d'une dislocation angulaire dans un 10 espace complet élastique isotrope tridimensionnel ou un demi-espace élastique isotrope tridimensionnel. 6. Procédé selon la revendication 1, dans lequel lesdits un ou plusieurs modèles de pression de discontinuité comprennent 15 un modèle de pression d'une faille, d'un dyke, d'une chambre magmatique, d'un dôme de sel, ou d'une combinaison de ceux-ci. 7. Procédé selon la revendication 1, comprenant en outre une minimisation d'une fonction de coût pour déterminer un ou 20 plusieurs paramètres d'optimisation afin de prédire l'attribut de contrainte. 8. Procédé selon la revendication 7, dans lequel la prédiction de l'attribut de contrainte du volume souterrain comprend en 25 outre l'application desdits un ou plusieurs paramètres d'optimisation à la valeur de contrainte calculée, la dérive de pression de discontinuité calculée à la surface de la Terre, la valeur de déformation calculée, la valeur de déplacement calculée, ou la combinaison de celles-ci. 30 9. Procédé selon la revendication 1, dans lequel lesdits un ou plusieurs modèles de contrainte en champ éloigné linéairement indépendants comprennent trois modèles de contrainte en champ éloigné linéairement indépendants qui reposent sur des jeux de données différents, chaque jeu de données comprenant des données de géométrie de faille, des données d'orientation de fractures, des données d'orientation de stylolithes, des données de plans de failles secondaires, des données de rejet de faille, des 59 3047338 données de stries de faille, des données de système de positionnement global (GPS), des données de radar à ouverture synthétique interférométrique (InSAR), des données de télémétrie laser, des données de clinomètres, des données de déplacement 5 d'une faille géologique, des données de grandeur de contrainte de la faille géologique, ou une combinaison de celles-ci. 10. Procédé selon la revendication 1, dans lequel l'attribut de contrainte prédit comprend une inversion de contrainte, un 10 champ de contrainte, une valeur de contrainte en champ éloigné, une interpolation de contrainte dans un réservoir faillé complexe, un champ de contrainte perturbée, un rapport de contraintes et son orientation, un ou plusieurs événements tectoniques, une discontinuité de déplacement d'une faille, un 15 rejet de faille, un déplacement estimé, une déformation perturbée, une distribution de rejet sur les failles, un contrôle de qualité des failles interprétées, une prédiction de fracture, une prédiction de propagation de fracture en fonction du champ de contrainte perturbée, un calcul en temps réel des 20 champs de contrainte et déplacement perturbés tout en exécutant une estimation interactive de paramètres, ou une différenciation d'une fracture induite d'une fracture préexistante. 11. Support lisible par ordinateur sur lequel sont mémorisées 25 des instructions qui, à leur exécution par un processeur, sont configurées pour amener le processeur à exécuter des opérations, les opérations comprenant : une simulation d'un ou de plusieurs modèles de contrainte en champ éloigné linéairement indépendants, d'un ou de plusieurs 30 modèles de pression de discontinuité, et d'un ou de plusieurs modèles de dérive de pression de discontinuité pour le volume souterrain en utilisant des données capturées par un capteur faisant partie d'un outil de mesure pendant le forage ou d'un outil de diagraphie durant le forage dans un puits de forage; 35 un calcul, au moyen d'un processeur, d'une valeur de contrainte, d'une dérive de pression de discontinuité à z = 0, d'une valeur de déformation, d'une valeur de déplacement, ou d'une combinaison de celles-ci pour des points de données dans 60 3047338 le volume souterrain sur la base d'une superposition desdits un ou plusieurs modèles de contrainte en champ éloigné linéairement indépendants, desdits un ou plusieurs modèles de pression de discontinuité, et desdits un ou plusieurs modèles de dérive de 5 pression de discontinuité ; et une prédiction de l'attribut de contrainte du volume souterrain, sur la base de la valeur de contrainte calculée, la dérive de pression de discontinuité calculée à z = 0, la valeur de déformation calculée, la valeur de déplacement calculée, ou 10 la combinaison de celles-ci. 12. Support lisible par ordinateur selon la revendication 11, dans lequel les opérations comprennent en outre une inversion, au moyen du processeur, desdits un ou plusieurs modèles de 15 contrainte en champ éloigné linéairement indépendants et desdits un ou plusieurs modèles de pression de discontinuité. 13. Support lisible par ordinateur selon la revendication 12, dans lequel lesdits un ou plusieurs modèles de pression de 20 discontinuité sont inversés par rapport à une différence entre une contrainte régionale inversée principale maximum et une contrainte régionale inversée principale minimum. 14. Support lisible par ordinateur selon la revendication 12, 25 dans lequel un ou plusieurs modèles de pression de discontinuité sont inversés sur la base au moins en partie d'une dislocation angulaire dans un espace complet élastique isotrope tridimensionnel ou un demi-espace élastique isotrope tridimensionnel. 30 15. Support lisible par ordinateur selon la revendication 11, dans lequel lesdits un ou plusieurs modèles de pression de discontinuité comprennent un modèle de pression d'une faille, d'un dyke, d'une chambre magmatique, d'un dôme de sel, ou d'une 35 combinaison de ceux-ci. 16. A système informatique, comprenant : un processeur ; et 61 3047338 un système de mémoire comprenant un ou plusieurs supports lisibles par ordinateur non transitoires sur lequel sont mémorisées des instructions qui, à leur exécution par le processeur, sont configurées pour amener le système informatique 5 à exécuter des opérations, les opérations comprenant : une simulation d'un ou de plusieurs modèles de contrainte en champ éloigné linéairement indépendants, d'un ou de plusieurs modèles de pression de discontinuité, et d'un ou de plusieurs modèles de dérive de pression de 10 discontinuité pour le volume souterrain en utilisant des données capturées par un capteur faisant partie d'un outil de mesure pendant le forage ou d'un outil de diagraphie durant le forage dans un puits dé forage; un calcul, au moyen d'un processeur, d'une valeur de 15 contrainte, d'une dérive de pression de discontinuité à z = 0, d'une valeur de déformation, d'une valeur de déplacement, ou d'une combinaison de celles-ci pour des points de données dans le volume souterrain sur la base d'une superposition desdits un ou plusieurs modèles de 20 contrainte en champ éloigné linéairement indépendants, desdits un ou plusieurs modèles de pression de discontinuité, et desdits un ou plusieurs modèles de dérive de pression de discontinuité ; et une prédiction de l'attribut de contrainte du volume 25 souterrain, sur la base de la valeur de contrainte calculée, la dérive de pression de discontinuité calculée à z = 0, la valeur de déformation calculée, la valeur de déplacement calculée, ou la combinaison de celles-ci. 30 17.. Système informatique selon la revendication 16, dans lequel les opérations comprennent en outre l'inversion, au moyen du processeur, desdits un ou plusieurs modèles de contrainte en champ éloigné linéairement indépendants et desdits un ou plusieurs modèles de pression de discontinuité. 35 18. Système informatique selon la revendication 17, dans lequel lesdits un ou plusieurs modèles de pression de discontinuité sont inversés par rapport à une différence entre une contrainte 62 3047338 régionale inversée principale maximum et une contrainte régionale inversée principale minimum. 19. Système informatique selon la revendication 17, dans lequel 5 un ou plusieurs modèles de pression de discontinuité sont inversés sur la base au moins en partie d'une dislocation angulaire dans un espace complet élastique isotrope tridimensionnel ou un demi-espace élastique isotrope tridimensionnel. 10 20. Système informatique selon la revendication 16, dans lequel lesdits un ou plusieurs modèles de pression de discontinuité comprennent un modèle de pression d'une faille, d'un dyke, d'une chambre magmatique, d'un dôme de sel, ou d'une combinaison de 15 ceux-ci.