FR3031210B1

FR3031210B1 - Simulation d'onde elastique sismique pour milieu isotrope transversalement incline utilisant une grille decalee de lebedev

Info

Publication number: FR3031210B1
Application number: FR1561115A
Authority: FR
Inventors: Fan Jiang; Shengwen Jin
Original assignee: Landmark Graphics Corp
Current assignee: Landmark Graphics Corp
Priority date: 2014-12-31
Filing date: 2015-11-19
Publication date: 2019-08-02
Anticipated expiration: 2035-11-19
Also published as: AR102739A1; GB2548729A; WO2016108896A1; US10241222B2; CA2969101A1; FR3031210A1; CA2969101C; GB2548729B; US20170336522A1; GB201708516D0

Abstract

Le présent exposé divulgue des systèmes et des procédés permettant de simuler numériquement la propagation d'une onde sismique dans un milieu isotrope transversalement incliné (TTI), en utilisant une grille décalée adaptative de Lebedev. Dans divers modes de réalisation, la grille adaptative comprend de multiples zones horizontales (400, 402, 502, 510) ayant différents espacements de grille (dz1, dz2) associés, qui peuvent être déterminés en fonction d'un modèle de vitesse d'onde vertical. La simulation numérique peut comprendre la résolution itérative d'un jeu d'équations de différence finie comprenant des coefficients de différence finie qui varient spatialement en fonction de l'espacement de grille (dz1, dz2). Des modes de réalisation et des caractéristiques additionnels sont décrits.

Description

ARRIERE-PLAN

[1] Dans l'exploration du gaz et du pétrole, les formations rocheuses ou d'autres milieux sous la surface sont souvent caractérisés en se basant sur des relevés sismiques, c'est à dire, des mesures sismiques en association à la modélisation informatique du milieu et la simulation de la propagation de l'onde sismique provenant des sources sismiques vers des récepteurs sismiques. Par exemple, dans une géométrie type d'acquisition de terrain, une ou plusieurs sources sismiques artificielles, telles que des explosives, des vibrateurs montés sur un camion ou des marteaux, peuvent être développés sur la surface de la terre pour exciter les ondes sismiques dans la formation rocheuse sous la surface, et une pluralité de récepteurs sismiques, tels que des géophones, peut être distribuée sur la surface autour de la ou des sources afin de mesurer les ondes sismiques produites par l'excitation (telles que les réflexions sur des limites géologiques). Afin d'obtenir des informations sur la formation (telles que les épaisseurs et des matériaux dans les diverses couches, et l'emplacement des réservoirs de pétrole ou de gaz dans celle-ci) à partir des mesures sismiques, les mesures peuvent être comparées aux résultats d'une simulation informatique qui est basée sur un modèle informatique reflétant les suppositions concernant la formation sous la surface. Les écarts entre les mesures et la simulation suggèrent des imprécisions dans le modèle.

[2] Les formations rocheuses comprennent souvent de multiples couches horizontales composées de différents matériaux. Lorsqu'elles agissent comme un milieu de propagation pour les ondes sismiques, de telles formations démontrent une isotropie transversale verticale, c.-à-d., leurs propriétés matérielles sont indépendantes de la direction de la propagation de l'onde dans un plan horizontal (c.-à-d., un plan perpendiculaire à l'axe de symétrie verticale de l'isotropie). Un tel milieu isotrope transversal vertical (VTI) peut être adéquatement modélisé avec des approches existantes de modélisation élastique de différence finie, qui impliquent généralement la résolution numérique des équations d'onde élastique discrétisées sur un volume en 3D approprié de la formation, à l'aide d'une grille de discrétisation appropriée. En ce qui concerne le milieu isotrope transversalement incliné (TTI pour, en anglais, « Tilted Transversely Isotropie »), dans lequel une limite géologique et, par conséquent, le plan de l'isotropie fait un angle relativement au plan horizontal, cependant, les approches conventionnelles ne performent généralement pas de façon satisfaisante en ce que, afin d'éviter les résultats imprécis ou des artefacts de simulation, elles utilisent un espacement de grille tellement fin qu'elles peuvent très vite devenir intraitables en terme de calcul.

PRESENTATION

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, un procédé comprenant : la définition de positions d'au moins une source sismique et d'au moins un récepteur sismique par rapport à un milieu isotrope transversalement incliné (TTI) ; la définition d'une grille décalée adaptative de Lebedev sur au moins une partie du milieu TTI, la grille comprenant une pluralité de zones horizontales avec des espacements de grille associés, un espacement de grille associé à au moins une des zones étant différent d'un espacement de grille associé à une autre des zones ; et l'utilisation d'un processeur, calculant la propagation d'une onde sismique émise par l'au moins une source sismique à travers le milieu TTI et la réception de celle-ci au niveau du récepteur en résolvant un jeu d'équations d'onde élastique discrétisées sur la grille décalée adaptative de Lebedev.

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, le calcul de la propagation de l'onde sismique est basé au moins en partie sur un modèle de calcul du milieu TTI comprenant un ou plusieurs paramètres de champ.

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, le procédé comprend en outre : l'émission d'une onde sismique avec l'au moins une source sismique ; la mesure d'une onde sismique résultant de l'onde sismique émise avec l'au moins un récepteur sismique ; la comparaison de l'onde sismique mesurée avec l'onde sismique calculée au niveau du récepteur et, si une différence entre les deux dépasse un seuil défini, l'ajustement du modèle de calcul du milieu TTI en ajustant un paramètre de champ du modèle.

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, les paramètres de champ comprennent une vitesse d'onde de pression, une vitesse d'onde de cisaillement et les paramètres inélastiques de Thomsen.

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, les paramètres de champ comprennent un angle d'inclinaison et un angle d'azimut associés à une limite géologique inclinée du milieu TTI.

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, la définition de la grille décalée adaptative de Lebedev comprend la détermination des zones horizontales et des espacements de grille associés à celles-ci basé au moins en partie sur le modèle de vitesse d'onde de cisaillement spécifiant une vitesse d'onde de cisaillement verticalement variable.

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, les équations d'onde élastique discrétisées sur la grille comprennent des équations de différence finie d'au moins du deuxième ordre.

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, les équations de différence finie sont au moins du quatrième ordre.

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, les équations de différence finie est au moins du seizième ordre.

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, les équations de différence finie comprennent des coefficients de différence finie, le procédé comprenant en outre le calcul des coefficients de la différence finie sur la base des espacements de grille.

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, les coefficients de différence finie sont calculés en utilisant au moins une décomposition en valeurs singulières ou au moins une optimisation de moindre carré.

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, les coefficients de différence finie sont variables à l'intérieur d'une région de chevauchement comprenant des parties adjacentes de deux des zones horizontales, et dans lequel les coefficients de différence finie sont constants à l'intérieur de chacune des zones horizontales à l'extérieur de la région de chevauchement.

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, le calcul de la propagation de Tonde sismique comprend la mise à jour de variables de champ utilisant les équations de différence finie, et dans lequel la mise à jour des variables de champ dans une région de chevauchement qui comprend des parties adjacentes de deux des zones horizontales est basé sur des variables de champ dans chacune desdites deux zones horizontales.

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, la mise à jour de la variable de champ dans la région de chevauchement comprend la création de points de grille imaginaires en interpolant entre des points de grille réels.

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, la définition des positions de Tau moins une source et de Tau moins un récepteur comprend la spécification d'un type d'acquisition.

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, le calcul de la propagation de Tonde sismique comprend l'application de conditions de limite de couche à correspondance parfaite convolutionnelle.

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, le procédé comprend également : l'affichage, utilisant un écran ou un appareil d'impression, de la propagation de Tonde sismique émise par Tau moins une source sismique à travers le milieu TTI.

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, un système comprend : au moins une source sismique pour l'émission d'une onde sismique dans un milieu TTI; au moins un récepteur sismique configuré pour détecter une onde sismique se propageant à travers le milieu TTI ; et un dispositif de calcul configuré pour recevoir des informations sur les positions de l'au moins une source sismique et l'au moins un récepteur sismique par rapport à la formation rocheuse, définir une grille décalée adaptative de Lebedev sur au moins une partie de la formation rocheuse, la grille comprenant une pluralité de zones horizontales avec des espacements de grille associés, un espacement de grille associé à au moins l'une des zones différant d'un espacement de grille associé à une autre des zones, basé au moins en partie sur le modèle de calcul du milieu TTI, le calcul de la propagation de l'onde sismique émise par l'au moins une source sismique à travers la formation et la détection de celle-ci au niveau du récepteur en résolvant un jeu d'équations d'onde élastique discrétisées sur la gille décalée adaptative de Lebedev, et la comparaison de l'onde sismique détectée avec l'onde sismique calculée au niveau du récepteur et, lorsqu'une différence entre les ondes sismiques détectée et calculée dépasse un seuil défini, l'ajustement du modèle de calcul en ajustant au moins un paramètre de champ de celui-ci.

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, le dispositif de calcul est configuré en outre pour déterminer les zones horizontales et les espacements de grille associés à celles-ci basé au moins en partie sur un modèle de vitesse d'onde de cisaillement spécifiant une vitesse d'onde de cisaillement verticalement variable.

Selon un ou plusieurs mode(s) de réalisation de la présente divulgation, les équations d'onde élastique discrétisées comprennent des équations de différence finie d'au moins du seizième ordre.

BRÈVE DESCRIPTION DES FIGURES

[3] La FIG. 1 est un organigramme d'un procédé de relevé sismique conformément à divers modes de réalisation.

[4] Les FIG. 2A et 2B représentent une vue verticale en coupe et une vue du haut, respectivement, d'une formation TTI conformément à divers modes de réalisation.

[5] La FIG. 3 est un diagramme d'une cellule de grille décalée de Lebedev, conformément à divers modes de réalisation.

[6] Les FIG. 4A et 4B sont des diagrammes de deux plans de grille verticaux, séparés par une demi-cellule dans la direction horizontale, d'un exemple d'une grille décalée adaptative de Lebedev comprenant deux zones ayant des espacements de grille différents, conformément à divers modes de réalisation.

[7] La FIG. 5 est un diagramme schématique d'une grille adaptative, illustrant une interpolation de point de grille virtuelle, conformément à divers modes de réalisation.

[8] Les FIG. 6A et 6B sont des diagrammes de squelettes verticaux d'une grille décalée adaptative, illustrant des indices d'espacement de grille, conformément à divers modes de réalisation, pour des points de sous-grille à des emplacements verticaux de demi-entier ou d'entier, respectivement.

[9] La FIG. 7 est un organigramme d'un exemple de procédé permettant de simuler une propagation d'onde sismique dans un milieu TTI utilisant une grille décalée adaptative de Lebedev, conformément à divers modes de réalisation.

[10] La FIG. 8 est un schéma de blocs d'un exemple de système permettant d'implémenter le procédé de la FIG. 7, conformément à divers modes de réalisation.

DESCRIPTION DETAILLEE

[11] La présente divulgation décrit une approche permettant de simuler une propagation d'onde sismique qui est taillée sur mesure pour un milieu TTI, permettant une simulation d'onde sismique précise avec des espacements de grille qui sont suffisamment grands pour être praticables avec les ressources de traitement et de stockage de données actuellement disponibles. Dans divers modes de réalisation, cette approche utilise une grille de discrétisation qui combine les caractéristiques d'une grille décalée de Lebedev avec les caractéristiques d'une grille adaptative. Dans une grille décalée, le jeu de variables de champ discrétisés (par ex., les composants de la vitesse des particules et du tenseur de contrainte des équations d'onde élastique) est divisé en multiples sous-ensembles qui sont stockés à des positions différentes à l'intérieur d'une cellule de la grille (la cellule étant l'unité répétitive la plus petite de la grille). Dans une grille de Lebedev (expliquée en détail en référence à la FIG. 3), spécifiquement, chacun des composants de la vitesse de particule et de contrainte est divisé en quatre variables, donnant quatre groupes de sous-grille pour les composants de la vitesse de particule et quatre groupes de sous-grille pour les composants de la contrainte, qui sont placés au niveau de 8 coins d'une cellule de demi-grille, respectivement. Une grille adaptative est caractérisée par différents espacements de grille utilisés dans des zones respectives différentes du volume en 3D discrétisé. L'espacement de grille à l'intérieur de chaque zone est généralement choisi en se basant sur la ou les valeurs d'un ou de plusieurs paramètres de champ (tels que la vitesse d'onde) à l'intérieur de cette zone, avec une vue sur la recherche d'un compromis approprié entre une précision numérique élevée et un faible coût en calcul. Par conséquent, les opérations informatiques peuvent être améliorées en ce que l'usage de la mémoire est réduit, et le temps de calcul est diminué tout en maintenant des résultats de modélisation précis.

[12] La FIG. 1 est un organigramme d'un exemple de procédé de relevé sismique 100, qui comprend à la fois les simulations d'onde sismique conformément à divers modes de réalisation et correspondant aux mesures sismiques, aussi bien que des comparaisons entre les résultats des simulations et les mesures. Les simulations sont basées sur un modèle de calcul de la formation, qui doit être ajusté de façon itérative jusqu'à ce que les mesures physiques et les résultats du calcul de la simulation soient raisonnablement en accord (par ex., défini par un certain degré choisi d'accord, ou de convergence) l'un avec l'autre. À ce moment, le modèle ajusté de la formation peut être utilisé, par ex., pour également planifier une opération de forage, ou à d'autres fins.

[13] Le volet simulation du procédé implique, au niveau de l'opération 102, la définition de la géométrie d'acquisition par définition des positions de la ou des sources sismiques et du ou des récepteurs relativement à la formation devant être relevée, par ex., en termes de et en référence à en ligne, coupe transversale et l'élévation. Dans ce contexte, l'élévation correspond à la hauteur de l'emplacement géographique au-dessus de la surface de la terre (compris comme étant une surface gravitationnelle équipotentielle au niveau des océans). En ligne et en coupe transversale décrivent des directions d'acquisition (par ex., les « directions de prise » de la ou des sources). Par ex., la direction ouest-est peut être utilisée comme la direction en ligne, et la direction sud-nord peut être utilisée comme la direction en coupe transversale.

[14] Le fait de préciser les emplacements de la source et du récepteur détermine le type d'acquisition. Divers types d'acquisition sont bien connus et utilisés en routine dans les relevés sismiques par les hommes de métier : acquisition sur terre, par ex., une pluralité de récepteurs sont agencés autour d'une ou de plusieurs sources déployées au niveau de la surface (c.-à-d., à une élévation de zéro) ; dans l'acquisition en milieu marin, la ou les sources et les récepteurs sont placés au niveau de la surface de l'océan (c.-à-d., également à une élévation de zéro) ; dans l'acquisition au fond de l'océan, la source est généralement placée au niveau de la surface de l'océan, alors que les récepteurs sont placés le long du plancher océanique (correspondant à une élévation négative) ; pour une acquisition de profil sismique verticale, l'une ou plusieurs des sources et/ou des récepteurs sont localisés en-dessous de la surface, par ex., à diverses profondeurs à l'intérieur d'un trou de forage ; et pour une acquisition micro-sismique, les récepteurs sont placés au niveau de la surface de la terre, mais les sources sont déployées en dessous de la surface. Les procédés divulgués ici sont généralement applicables à tous ces procédés, et à d'autres. En raison des coûts élevés associés à l'excitation de multiples sources, une source unique est souvent utilisée en association avec de multiples récepteurs. Cependant, en principe, il est également possible d'utiliser de multiples sources en association à un récepteur unique. En outre, afin d'augmenter la quantité d'informations obtenue lors du relevé (par ex., dans le but d'obtenir une meilleure précision ou résolution spatiale, ou pour résoudre des ambiguïtés), certains relevés utilisent une pluralité de sources et une pluralité de récepteurs, par ex., agencées le long de deux jeux mutuellement perpendiculaires de lignes parallèles pour les sources et les récepteurs, respectivement.

[15] Comme le démontre la FIG. 1, le procédé 100 comprend également la mise à disposition, pour une utilisation dans des simulations ultérieures, d'un modèle de calcul de la formation sous la surface (voir l'opération 104). Le modèle peut comprendre des paramètres indicatifs des propriétés de la géométrie et/ou du matériau des diverses couches de la formation (ou d'autres sous-formations) (tels qu'une épaisseur de couche, des positions et des orientations, des densités du matériau, des coefficients élastiques, etc.) et/ou d'autres paramètres physiques dépendants de ou dérivés des propriétés de la géométrie et du matériau (telles que les vitesses de propagation des ondes). Ces paramètres sont généralement des paramètres de champ, c.-à-d., ils sont définis sous forme de fonction de et peuvent varier dépendamment de l'emplacement à l'intérieur de la formation. Le type et le nombre de paramètres appropriés pour modéliser une certaine formation dépendent généralement des propriétés de symétrie de la formation. Par ex., les milieux isotropes sont souvent caractérisés en termes de densité de matériaux p et de vitesse de propagation Vpo et Vso des ondes acoustiques et de cisaillement (qui sont également souvent appelées ondes « primaire » et « secondaire », respectivement) dans la direction verticale. Pour le milieu VTI aussi bien que pour le milieu isotrope transversal horizontal (HTI) (dans lequel l'axe de symétrie est horizontal et les diverses couches de matériaux se trouvent généralement dans des plans verticaux), trois paramètres additionnels caractérisant l'anisotropie — tels que les paramètres anisotropiques de Thomsen ε, δ et γ, qui sont bien connus dans le domaine — sont généralement précisés (au-delà de la densité p et des vitesses de l'onde acoustique et de l'onde de cisaillement

Vpo, Vs0).

[16] Le milieu TTI, qui présente un intérêt dans divers modes de réalisation de la présente divulgation, comprend généralement au moins une limite inclinée entre deux formations sous la surface, et il est également caractérisé par un angle d'inclinaison Θ et un angle d'azimut φ de cette limite inclinée. Comme il est illustré dans la FIG. 2A dans une vue verticale en coupe à travers la formation, l'angle d'inclinaison Θ est défini comme l'angle qu'une normale 200 de la limite 202 entoure avec l'axe vertical 204 (ou, de façon équivalente, l'angle entre le plan de, ou un plan tangentiel à, la limite 202 et un plan horizontal 206). Comme il est illustré dans la FIG. 2B dans une vue du haut de la formation, prise en association à la FIG. 2A, l'angle d'azimut φ de la limite inclinée 202 est l'angle auquel la projection 208 de la normale 200 jusque dans le plan horizontal entoure avec la direction de la coupe transversale 210. La limite 202 ne doit pas nécessairement être plane à travers toute l'intégralité de la région modélisée, mais peut avoir une courbure non-zéro ; ainsi, l'angle d'inclinaison θ et l'angle d'azimut φ peuvent varier en fonction de la position (en d'autres termes, ce sont des paramètres de champ). Le fait de considérer l'angle d'inclinaison θ et l'angle d'azimut φ comme des paramètres de champ permettent également de décrire de multiples limites géologiques qui pourraient être présentes à l'intérieur de la formation. En conformité avec divers modes de réalisation, les milieux TTI sont modélisés avec sept paramètres de champ : l'angle d'inclinaison Θ, l'angle d'azimut φ, la densité p, les vitesses de propagation de Tonde verticale VpQ et Vso. et les paramètres anisotropiques de Thomsen ε, δ et γ.

[17] De nouveaux en référence à la FIG. 1, lors de la préparation de la simulation d'une propagation d'ondes sismiques à travers un milieu modélisé, au niveau de l'opération 106, une grille de discrétisation est définie au-dessus d'une région d'intérêt (qui peut être coïncidante avec le milieu modélisé, ou comprendre seulement une partie de celui-ci). Dans divers modes de réalisation, la grille de discrétisation est une grille décalée de Lebedev, tel qu'il est expliqué en référence à la FIG. 3 ci-dessous. L'espacement de la grille peut être choisi en fonction d'un ou de plusieurs paramètres de champ du modèle. Par ex., pour une région de la grille avec un espacement uniforme de la grille, l'espacement de la grille dz peut être déterminé en se basant sur la vitesse de cisaillement minimal VSmin dans cette région, en association à la fréquence maximale fmax des ondes sismiques qui sont modélisées et le nombre de points de grille ngrid par longueur d'onde minimale, en conformité à la relation de dispersion du milieu isotrope :

[18] Dans une grille adaptative, telle qu'elle est utilisée dans divers modes de réalisation, l'espacement de la grille varie entre différentes zones de la région modélisée. Par ex., dans certains modes de réalisation, la région modélisée est divisée, le long d'une direction verticale, en de multiples zones horizontales c.-à-d., des zones séparées par des limites horizontales) basée sur une distribution verticale (c.-à-d., une variation dans la direction z) de la vitesse de Tonde de cisaillement Vs0 (x, y, z) (ou, si la vitesse de Tonde de cisaillement est égale à zéro, basée sur la distribution verticale de la vitesse de Tonde acoustique Vp0). Souvent, plus profonde est la position à l'intérieur de la formation, plus la vitesse est profonde. Dans certains modes de réalisation, si la vitesse est constante entre une profondeur Ζχ et une profondeur Z2 (par ex., à l'intérieur d'une couche d'eau), cette région constitue une zone de vitesse. En outre, si une valeur de vitesse anormale est trouvée, par ex., à l'intérieur d'un corps salé (qui tend à avoir une vitesse beaucoup plus

élevée et une densité plus faible que la roche environnante), la couche comprenant ce corps salé peut être considérée comme une zone. Les régions avec des vitesses anormalement faibles peuvent être fusionnées avec la couche qui se trouve au-dessus de celles-ci. Le nombre et la taille des zones horizontales peuvent être choisis pour équilibrer le compromis entre la précision et l'efficacité. Pour différents modèles de vitesse, la région modélisée peut être divisée, par ex., en 8 zones, ou seulement en 2 zones. Un nombre de zones particulièrement bénéfique pour la grille adaptative peut être obtenu si la vitesse augmente avec la linéarité de la profondeur, dans quel cas, la région modélisée peut être divisée en zone de profondeur égale. Une fois que l'espacement de la grille et l'étendue verticale des zones ont été déterminés, la taille de la grille et les tailles de la zone (c.-à-d., le nombre total de points de grille ou de point de grille par zone, dans chaque dimension) sont déterminées en se basant sur ceux-ci.

[19] Suivant la définition de la grille de discrétisation, les ondes sismiques peuvent être simulées, dans l'opération 108, en résolvant les équations des ondes élastiques discrétisées. Dans un milieu élastique isotrope, les équations d'ondes élastiques applicables (également appelées les équations de vitesse de contrainte) sont :

Dans laquelle, vÉsont des composants de la vitesse de particule et (Tjj sont les composants du tenseur de contrainte (i,j = x, y, z), qui constituent les variables de champ du modèle (c.-à-d., la variable pour laquelle les équations sont numériquement résolues) ;

rindiquent les dérivées de premier ordre des composants de la vitesse de particules et de la contrainte par rapport au temps ;

-indiquent les dérivées spatiales de premier ordre des composants de la vitesse de particules et de la contrainte dans la direction de k (k = x, y, z); Sÿ est le tenseur de Kronecker ; et λ et μ sont les constantes de Lame. La convention de somme est utilisée, c.-à-d., = Vxx + Vyy +

Vzz. Pour le milieu TTI, les constantes de Lame sont remplacées par le tenseur de raideur C/7, de sorte que la deuxième des équations susmentionnées devient :

Dans laquelle, Cÿ est calculée à partir du tenseur de raideurC^ pour le milieu VTI et la matrice de transformation de liaison R, qui est une fonction de l'angle d'inclinaison ά et l'angle d'azimut φ, selon Cÿ = RC^R^ (où R^est la matrice de transposition de R). L est la matrice du tenseur de raideur bien connue dans le milieu VTI :

Les composants individuels du tenseur de raideur C dépendent des, et peuvent être directement calculé à partir des, propriétés du matériel et des paramètres de champ apparenté du milieu modélisé, tels que la densité p, les vitesses de propagation de l'onde Vp0 et Vs0, et les paramètres anisotropiques de Thomsen ε, δ et γ.

[20] Les équations de contrainte-vélocité peuvent être discrétisées conformément à un schéma de différence finie dans lequel les dérivées spatiales des variables de champ sont exprimées, pour chaque point de grille (ou point de sous-grille, tel qu'il est applicable dans des grilles décalées et tel qu'il est expliqué ci-dessous), en termes d'une combinaison linéaire des valeurs de la variable de champ de multiples emplacements de grille entourante/point de sous-grille, avec des coefficients linéaires (dans ce contexte également appelé « coefficients de différence finie ») qui dépendent de l'espacement de la grille et sont, par conséquent, constants pour un espacement de grille uniforme, mais spatialement variables dans une grille adaptative possédant de multiples zones d'espacement de grille différents. Les coefficients de différence finie peuvent être calculés au moment de la définition de la grille, avant qu'il ne résolve de façon itérative les équations de contrainte-vélocité. Les grilles entourantes/points de sous-grille peuvent être choisis pour être distribués symétriquement autour de la grille/du point de sous-grille au niveau duquel la dérivée spatiale est évaluée. Le nombre de grilles entourantes/point de sous-grille utilisé correspond à Tordre du schéma de la différence finie ; plus élevé est Tordre, plus élevé sera à la fois la précision de la simulation et le coût de calcul associé. En conformité avec ceci, le schéma de différence finie est généralement de 2e ordre ou d'un ordre plus élevé ; dans certains modes de réalisation, un schéma de différence finie de 16e ordre (ou d'un ordre plus élevé) est utilisé.

[21] La résolution numérique des équations de vitesse-contrainte discrétisées implique généralement l'entrée dans une boucle temporelle dans laquelle les composants de vitesse de particule et de contrainte sont mises à jour de façon itérative (en utilisant des coefficients de différence finie). L'intervalle de temps (c.-à-d., l'inverse de la vitesse d'échantillonnage numérique) peut être déterminé en se basant sur l'espacement minimal de la grille dzmm, la vitesse maximale de Tonde de pression Vpmax, les coefficients de différence finie, et la dimensionnalité Z) du modèle (par ex., 2 pour la modélisation en 2D et 3D pour la modélisation en 3D) en conformité avec :

2|cî| représente la somme des valeurs absolues des coefficients de différence finie. Pour les grilles adaptatives, l'intervalle de temps peut être calculé pour chaque zone de vitesse, et ensuite l'intervalle de temps le plus petit parmi toutes les zones peut être choisi pour la simulation.

[22] Les opérations de modélisation et de simulation 102-108 permettent une quantification des ondes sismiques (par ex., en termes de composants de la vitesse de particule ou d'autres paramètres physiques qui sont directement calculés à partir de celui-ci) reçues au niveau des emplacements de récepteur comme une fonction du temps. Afin de tester la validité du modèle de la formation, ces résultats simulés peuvent être comparés à des mesures sismiques. Ainsi, le procédé de relevé sismique 100 comprend également un volet de mesure, qui implique, au niveau de l'opération 110, la définition physique d'une ou de plusieurs sources sismiques et d'un ou de plusieurs récepteurs sismiques au niveau des emplacements (par rapport à la formation) cohérents avec la géométrie d'acquisition mentionnée dans l'opération 102 à des fins de simulation. Ensuite, au niveau de l'opération 112, les ondes sismiques peuvent être excitées dans la formation en utilisant une ou des sources sismiques (par ex., par la détonation d'explosifs, un tir de pistolet à air ou en frappant le sol avec un marteau), et, au niveau de l'opération 114, les ondes sismiques ainsi obtenues au niveau des emplacements de récepteur peuvent être mesurées. Les résultats de la simulation et de la mesure peuvent être comparés dans l'opération 116. En cas de décalage dépassant un seuil prédéfini, le modèle de la formation peut être ajusté (par ex., en modifiant légèrement les paramètres de champ) (opération 118), et la simulation (opération 106, 108) peut être répétée avec le nouveau modèle. Une fois que la simulation et la mesure sont en accord (dans les limites de la tolérance définie par le seuil de décalage), le modèle de formation peut être jugé comme étant précis, et peut être utilisé (opération 120) par des géophysiciens, des géologues, des ingénieurs de forage, et d'autres, par ex., comme point de départ pour d'autres évaluations et/ou d'imagerie de la formation, pour motiver les décisions de forage et/ou pour guider les opérations de forage, etc.

[23] En se tournant maintenant vers les divers détails du volet simulation du procédé de relevé sismique 100, la figure 3 illustre une cellule de grille individuelle 300 d'une grille décalée de Lebedev telle qu'elle est utilisée, en conformité avec divers modes de réalisation, pour discrétiser les équations d'onde élastiques dans un milieu TTI. La cellule de grille 300, qui peut être cubique, est associée à un emplacement de grille (x, y, z), indiqué au niveau de 302, et se prolonge le long de 3 dimensions perpendiculaires de la grille par une longueur d'unité dans chaque direction positive (de sorte qu'elles soient définies par 8 coins à (x, y, z), (x+1, y, z), (x,y+l, z), (x, y, z+1), (x+1, y+1, z), (x+1, y, z+1), (x, y+1, z+1) et (x+1, y+1, z+1)). Les variables de champ associés au point de la grille (x, y, z) sont localisées au niveau des coins de la cellule de demi-grille 304 (c.-à-d., une cellule d'une longueur d'une demie unité se prolongeant du point de grille (x, y, z): (x, y, z), (x+1/2, y, z), (x,y+1/2, z), (x, y, z+1/2), (x+1/2, y+1/2, z), (x+1/2, y, z+1/2), (x, y+1/2, z+1/2) et (x+1/2, y+1/2, z+1/1). Tel qu'il est démontré, les composants de la vitesse de particules et de contraintes sont « décalés », c.-à-d., placé au niveau de différents emplacements de la cellule de demi-grille. En particulier, chaque composant de la vitesse de particule Vz (z = x, y, z) est divisé en 4 variables Vz/, Vzj, Vzjet Vz4 qui sont placées, respectivement, au niveau du point de la grille 302 (x, y, z) et des 3 coins diagonalement opposés des 3 faces de la cellule de demi-grille 304 se rejoignant au niveau du point de la grille 302 (illustrés sous forme de cercles). De la même façon, chaque composant de contrainte ^(z,7 = x,y,z) est divisé en 4 variables (Tjji, (7^2, <7ÿ3 et <T//4 ces variables sont localisés au niveau des 4 coins restant de la grille de demi-cellule 304 (illustrés sous forme de triangle). La grille décalée de Lebedev dans son intégralité (c.-à-d., comprenant tous les points de la grille) possède donc des composants de vitesse de particule localisés sur les points de la grille et au niveau des centres des faces des cellules de la grille, alors que les composants de contraintes sont localisés au niveau des centres des cellules de la grille et au niveau des centres des rebords des cellules de la grille. Dans la discussion qui suit, tous les points d'une grille décalée au niveau desquels des variables de champ sont placées sont appelés des « points de sous-grille », alors que le terme « point de la grille » est réservé aux points à valeur d'entier au niveau desquels les cellules de la grille sont ancrées. (Avec cette terminologie, les points de la grille forment un sous-ensemble de points de sous-grille).

[24] Avec la division des variables de champ en 4 groupes de sous-grille, les dérivées temporelles des composants de vitesse de particules s'apparentent aux dérivées spatiales des composants de contraintes comme suit :

avec (i, j, k) = (x, y, z), (y, z, x) ou (z, x, y).

De la même façon, les dérivées temporelles des composants de contrainte peuvent être divisées en 4 chacun :

Encore une fois, (i, j, k) = (x, y, z), (y, z, x) ou (z, x, y). Le Cjj [remplacé « xy » par « ij »] est spatialement dépendant étant donné qu'il est dérivé de la vitesse VpQ et Vsq, les paramètres anisotropes ε, δ et y, et l'angle d'inclinaison et l'angle d'azimut. En outre, étant donné que djj est défini au niveau d'une position de la grille, les valeurs relatives de Cn θ Ci6 sont différentes dans les 4 équations susmentionnées. De façon bénéfique, avec ces relations pour les 4 sous-grilles, la rotation d'un gradient et la divergence d'une rotation s'évaporent de façon inhérente. En outre, il n'y a pas lieu d'interpoler les dérivées spatiales. Suivant la solution numérique des équations d'onde élastique, les composants de la vitesse de particules et de contraintes au niveau de chaque point de la grille peuvent être calculés en additionnant les composants respectifs des 4 groupes de sous-grille

[25] Dans une grille décalée de Lebedev, les dérivées spatiales des variables de champ peuvent être exprimées légèrement différemment pour différents jeux de variables. Par ex., dans un schéma de différence finie de 4e ordre, les dérivées spatiales verticales (c.-à-d., dans la notation des coordonnées de la FIG. 3, la dérivée spatiale par rapport à z) des variables de champ

(où i,j,k = x.,y,z) (qui sont les variables de champ localisés au niveau des points de sous-grille avec des valeurs de demi entier de z) peuvent être exprimés sous forme de

et les dérivées spatiales verticales des variables de champ G — » @jk4 (°ù ijifc = x/y,z) (qui sont localisés au niveau des points de sous-grille avec des valeurs d'entier de z) peuvent être exprimées sous forme de

À noter que la sélection des emplacements de point de la grille dans ces exemples d'équations est telle que les variables de champ sont mises à jour en se basant sur d'autres variables de champ à des points de sous-grille agencés symétriquement autour de la variable devant être mise à jour. Par ex., la dérivée du temps de située à (x, y, z), dépend de la dérivée

, qui, à son tour est calculé à partir de (Jxzi (jV, J/, Z-2) (localisée, dans une grille décalée de Lebedev, à z-3/2), (7xzi (jV, J2, Z- 7) (localisé à z-1/2), (Jxz/ (x, y, z) (localisé à z+1/2) et (Txzi (x, y, Z+7) (localisé à z+3/2).) [26] Afin de réduire l'utilisation de la mémoire et économiser sur le temps de calcul, divers modes de réalisation utilisent une grille adaptative, c.-à-d., une grille comprenant de multiples zones avec différents espacements de grille. Par ex., la grille peut être divisée le long de la dimension verticale (z) en de multiples zones, les espacements de grille dépendant du profil vertical de la vitesse de l'onde de cisaillement. Les FIGS. 4A et 4B illustrent une grille décalée de Lebedev avec 2 zones 400, 402 séparées par une limite horizontale 404 juste au-dessus de z = k, montrant les plans x-z de la grille à des positions d'entier ou de demi-entier de y, respectivement.

[27] Lors de la résolution des équations discrétisées de Tonde élastique, les variables de champ dans la plupart des cellules de la grille peuvent être mises à jour en utilisant seulement les cellules environnantes de la grille à l'intérieur de la même zone. À proximité de la limite de la zone 404, cependant, certains variables de champ sont mises à jour en se basant sur des cellules de la grille des 2 zones 400, 402. Par conséquent, une région de chevauchement entourant la limite entre les 2 zones peut être définie comme la région dans laquelle la cellule de la grille comporte au moins une variable de champ associée qui est mise à jour en se basant en partie sur une cellule de la grille de l'autre zone. Comme il sera facilement compris par les hommes de métier, la taille de la zone de chevauchement dépend généralement de Tordre du schéma de la différence finie, et augmente avec l'accroissement de Tordre.

[28] En référence aux FIGS. 4A et 4B, pour un schéma de différence finie de 4e ordre, la région de chevauchement 406 comprend 4 couches de cellules de la grille, qui sont associés aux points de la grille au niveau de z = k-2, k-1, k et k+1. (À noter que la région de chevauchement 402 comprend son plan de limite supérieure au niveau de z = k-2, mais son plan de limite inférieure au niveau de z = k+2.) Par ex., les dérivées de temps des variables de champ placées au niveau de z = k-3/2 (qui comprennent associées avec des points de la grille au niveau de z = k-2), indiquées généralement au niveau de 408, sont évaluées en utilisant les valeurs des variables de champ placé au niveau de z = k-3, k-2 et k-1 (qui appartiennent à la zone 400 avec un espacement de grille plus petit) et des variables de champ placées au niveau de z = k (qui appartiennent à la zone 402 avec un espacement de grille plus grand) ; par conséquent, les cellules de la grille au niveau de z = k-2 appartiennent à la région de chevauchement. De la même façon, les dérivées de variables de champ placées au niveau de z = k+1 (qui comprennent Vjp Vjg, @jk2» ^/fe4 associées aux points de la grille au niveau de z = k+1), indiquées généralement au niveau de 410, sont évaluées en utilisant les valeurs des variables de champ placées au niveau de z = k-1/2 (qui appartiennent à la zone 400) et des variables de champ placées au niveau de z = k+1/2, k+3/2 et k=5/2 (qui appartiennent à la zone 402) ; par conséquent, les cellules de la grille au niveau de z = k+1 appartiennent à la région de chevauchement. Par contre, en utilisant le schéma de 4e ordre du schéma de la différence finie décrit ci-dessus, les variables de champ associées aux points de la grille au niveau de z = k-3 et plus bas peuvent être calculées à partir des variables de champ complètement à l'intérieur de la zone supérieure 400, et les variables de champ associées aux points de la grille au niveau de z = k+2 et plus peuvent être calculées à partir de la variable de champ complètement à l'intérieur de la zone inférieure 402, comme il sera facilement compris par un homme de métier. En étendant le schéma de la différence finie décrit ci-dessus à des ordres plus élevés, le nombre de couches de cellules compris dans la région de chevauchement peut être égal à l'ordre. Ainsi, pour un schéma de la différence finie d'un 16e ordre, conformément aux divers modes de réalisation, une région de chevauchement peut comprendre 16 couches des cellules de la grille, 8 dans chacune des 2 zones 400, 402.

[29] Étant donné que la grille et les points de la sous-grille de 2 zones adjacentes ayant des espacements de grille différents ne s'alignent généralement pas horizontalement, l'évaluation d'une dérivée spatiale verticale au niveau d'un point de la sous-grille dans la région de chevauchement peut dépendre d'un emplacement dans lequel aucun point de sous-grille n'existe. Ceci est illustrée dans la FIG. 5, dans laquelle la dérivée spatiale verticale au niveau du point de la sous-grille 500, localisé au niveau de z = k-1 à l'intérieur d'une première zone 502, dépend des points de la sous-grille 504, 505, 506 au niveau de z = k-5/2, k-3/2 et k-1/2 à l'intérieur de la même zone 502, et, de façon souhaitable, sur un point 508 au niveau de z = k +1/2, à l'intérieur d'une 2e zone 510, qui possède les mêmes coordonnées x-y que le point de la sous-grille 500. Le point 508 ne correspond cependant pas à un point de la sous-grille dans la 2e zone 510. Dans divers modes de réalisation, ce problème est résolu en créant un « point virtuel de la sous-grille » au niveau d'un emplacement (hors grille) 508, par ex., en interpolant entre des points avoisinants de la sous-grille 512, 513, 514, 515 au niveau de la même position verticale z = k+1/2.

[30] Dans divers modes de réalisation, les coefficients de la différence finie Cj (par ex., avec z = 1, 2, 3, 4 pour un schéma de différence finie de 4e ordre) utilisé à l'intérieur des expressions pour les dérivées spatiales

des variables de champ G sont calculées en fonction des distances entre les points de la sous-grille au niveau desquels la dérivée doit être évaluée (qui peut être le point de la grille (x, y, z) lui-même ou tout autre point à l'intérieur de la cellule de la grille associée, tels que le point au niveau de (x, y, z+1/2)) et les points de la sous-grille environnante utilisés dans l'expression ; ces distances sont également appelées les « indices d'espacement de la grille ». Dans une grille décalée adaptative, les indices d'espacement de la grille peuvent différer non seulement entre les diverses zones de différent espacement de grille, mais également entre diverses variables de champ localisées au niveau de différents points de la sous-grille à l'intérieur d'une cellule. Les FIGS. 6A et 6B illustrent les indices d'espacement de la grille dl, d2, d3 et d4 et d'un schéma de différence finie de 4e ordre pour les variables de champ au niveau d'un point de la sous-grille 600 avec une valeur de demi-entier de z et au niveau d'un point de la sous-grille 602 avec une valeur d'entier de z, respectivement, dans une cellule ancrée au niveau de z = k, juste en dessous de la limite de la zone. En fonction de l'emplacement du point de la sous-grille au niveau duquel la dérivée est évaluée par rapport à la limite de la zone, certains indices de l'espacement de la grille peuvent être symétriques ou non autour de ce point. Par ex., comme le démontre la FIG. 6A, les distances d3, d4 à partir du point de la sous-grille 600 par rapport à son voisin le plus proche des points de la sous-grille 604, 605 sont d'égale longueur, alors que les distances dl, d2 par rapport aux voisins les plus proches en 2e position 606, 607 sont différentes en raison des espacements de grille différents. Dans la FIG. 6B, les 4 indices d'espacement de grille sont différents.

[31] Les coefficients de différence finie peuvent être déterminés à partir d'un système linéaire d'équations qui est formé en rapprochant les exponentielles des dérivées spatiales des variables de champ avec une expansion de Taylor. Par ex., pour un schéma de différence finie de 4e ordre et utilisant la notation de l'indice d'espacement de la grille des FIG. 6A et 6B, les coefficients de différence finie Cz (z = 1, 2, 3, 4) peuvent être calculés en utilisant l'équation matricielle suivante :

Cette équation matricielle peut être résolue avec les techniques bien connues des hommes de métier, telles que la décomposition en valeurs singulières (pour lesquelles des résolutions commerciales sont facilement disponibles).

[32] Dans divers modes de réalisation, la résolution numérique des équations de l'onde élastique au-dessus d'une région modélisée d'intérêt implique l'application des conditions de limite spéciale aux limites de cette région afin d'éviter les artefacts numériques tels que la réflexion de la limite numérique (qui est la réflexion apparente des ondes sismiques réfléchies sur la limite (numérique), malgré l'absence de limites physiques coïncidante). Ces conditions de limite peuvent prendre la forme, par ex., de conditions de limite de couche convolutive à correspondance parfaite (C-PML), qui sont généralement connues des hommes de métier dans une forme appropriée pour être utilisées dans les grilles régulières (non-décalées). Afin d'incorporer les conditions de limite C-PML avec les grilles décalées de Lebedev, chaque condition de limite peut être divisée (tout comme les équations d'onde élastique elles-mêmes) en 4 équations, correspondant à 4 sous-grilles, afin d'obtenir des conditions de limite absorbantes pour chaque sous-grille. Par ex., une première condition de limite, C-PML1, peut être utilisée pour absorber vxi au niveau de la limite.

[33] La différence principale entre les conditions C-PML et les équations d'onde élastique applicables au gros de la région simulée réside dans leurs dérivées spatiales. Pour une grille non-décalée, la dérivée spatiale dx applicable à la limite est réécrite en format C-PML sous forme de :

Ici, &e est un réseau C-PML temporaire qui applique les coefficients absorbants à la dérivée spatiale, et bx et kx sont des facteurs atténuants. La dérivée spatiale (da^/dx) PML PeUt θίΓθ appliqué® dans '/θΠυ3ί'οη d'onde élastique au niveau de chaque limite pour absorber Tonde là-dessus afin d'éviter la réflexion sur la limite. (Des équations analogues s'appliquent à (Jyy et (JZz)· Pour appliquer les conditions C-PML à la grille décalée de Lebedev, les dérivées spatiales sont modifiées comme suit :

Ici, et Bxj_sont recalculés en se basant sur la position de la grille de Lebedev, et kx est un jeu de 1 (et donc rejeté de l'équation). Pour qui est un entier au niveau du point de la grille (x,y,z), et fi^sont les mêmes que les facteurs d'atténuation conventionnels CLX, bx. Pour qui est un demi entier au niveau du point de la grille (x+1/2, y+1/2, z+1/2), et bx^ sont définies pour absorber l'énergie provenant du point de demi-grille. Chaque facteur d'atténuation est spécifié pour chaque dérivée spatiale.

[34] La FIG. 7 illustre un exemple de procédé 700 permettant de simuler une propagation d'onde sismique dans un milieu TTI utilisant une grille décalée de Lebedev adaptative, conformément à divers modes de réalisation. Le point de début du procédé 700 est un modèle de calcul du milieu TTI (spécifiant, par ex., les angles d'inclinaison et d'azimut de la limite géologique inclinée, les vitesses d'onde acoustique et de cisaillement et les paramètres anisotropes de Thomsen), en association avec une spécification de la géométrie d'acquisition. Le procédé 700 implique la division de la région d'intérêt dans de zone horizontale différente (par ex., séparée par des plans de limite horizontaux) basée sur la distribution verticale de la vitesse de Tonde de cisaillement Vsq (et, si Vsq est zéro, de la vitesse de Tonde acoustique Vpo) (opération 702). Pour chaque sous-zone, un

espacement de grille approprié peut ensuite être calculé, par ex., en se basant sur la relation de dispersion susmentionnée, et basé sur l'espacement de la grille, la taille de la grille (c.-à-d., le nombre de points de la grille) à l'intérieur de chaque zone peut être déterminée (opération 704). À partir des espacements de la grille, les coefficients de différence finie peuvent ensuite être calculés (pour un schéma de différence finie donné d'un ordre spécifié) (opération 706). Ceci peut impliquer la construction d'une équation matricielle du coefficient de différence finie, par ex., basée sur une expansion de la série Taylor (opération 708), et la résolution de l'équation matricielle par, par ex., une décomposition en valeurs singulières ou un algorithme d'optimisation des moindres carrés (opération 709). Les coefficients de différence finie ainsi obtenus peuvent ensuite être stockés dans un ou plusieurs réseaux pour une récupération et une utilisation ultérieures (opération 710). Au niveau de l'opération 712, une grille décalée de Lebedev est implémentée en attribuant une mémoire à toutes les variables de champ dans toutes les cellules de la grille, les variables de champ comprenant 4 variables distinctes pour chaque trois composants de la vitesse de la particule et six composants du tenseur de contraintes, localisés à diverses positions dans la sous-grille.

[35] Les équations de vitesse-contrainte, discrétisées sur une grille décalée adaptative de Lebedev peuvent ensuite être résolues de façon itérative dans une boucle temporelle 714, dans laquelle, en alternance, les composants de la vitesse de particules et les composants de contraintes sont mis à jour dans les 4 sous-groupes de Lebedev (opération 716, 718). À l'extérieur de la ou des régions de chevauchement entre les différentes zones, les vitesses de particules peuvent être mises à jour à l'intérieur de chaque zone conformément à un schéma régulier de différence finie utilisant des coefficients constants de différence finie (calculées en se basant sur l'espacement de la grille à l'intérieur de cette zone) (opération 720). À l'intérieur de la région de chevauchement, les vitesses de particule peuvent être mises à jour en se basant sur les coefficients variables de différence finie et en utilisant des points imaginaires de la grille, s'il y a lieu, (opération 722), de la façon décrite ci-dessus. Les valeurs de la variable de champ d'un quelconque point imaginaire de la grille peuvent être adaptées par interpolation (par ex., entre les points avoisinants dans le plan horizontal, y compris le point imaginaire de la grille) (opération 724). (L'ordre des opérations 722 et 724 n'est généralement pas crucial étant donnés qu'ils n'affectent pas le résultat de la simulation de façon importante. Dans certains modes de réalisation, la mise à jour des composants de vitesse dans la région de chevauchement avant l'interpolation horizontale permet d'avoir certains bénéfices en matière d'efficacité de calcul). Les opérations 720, 722, 724 peuvent ensuite être répétées pour les composants de contrainte (opérations 726, 728, 730). La boucle temporelle 714 peut être traversée autant de fois que souhaité afin de recouvrir, avec la simulation, une période de temps d'intérêt. Par ex., afin de comparer la simulation aux mesures sismiques physiques, la simulation peut se prolonger sur une période de temps comparable à la durée de mesure totale à partir de l'excitation d'une onde sismique jusqu'à la dernière mesure de l'onde sismique avec un détecteur.

[36] La FIG. 8 est un schéma de blocs d'un exemple de systèmes permettant d'implémenter le procédé de la FIG. 1, conformément à divers modes de réalisation. Le système comprend généralement un système de simulation 800 permettant de réaliser des calculs de simulation de l'onde sismique, tel qu'il est décrit ici, et un système de mesure 802 permettant d'acquérir et de traiter des mesures physiques de l'onde sismique. Le système de mesure 802 comprend généralement une ou plusieurs sources sismiques 804 (telles qu'un explosif avec le circuit de déclenchement associé, un marteau, un fusil à air, etc.) et un ou plusieurs récepteurs sismiques 806 (tels que des géophones, des hydrophones, etc.), aussi bien qu'une unité de commande et de traitement 808 en communication avec la ou les sources 804 et le ou les récepteurs 806 pour contrôler leur fonctionnement (par ex., pour correctement minuter l'acquisition du signal par rapport à l'excitation des ondes sismiques) et le traitement des données acquises par le ou les récepteurs 806. L'unité de commande et de traitement 808 peut être implémentée avec une quelconque combinaison appropriée de matériel et/ou de logiciels, par ex., un ordinateur polyvalent exécutant des programmes logiciels appropriés, ou un ordinateur dédié (par ex., un processeur de signal numérique, un réseau de portes programmables par l'utilisateur, etc.).

[37] Le système de simulation 800 peut, de la même façon, être implémenté par une quelconque combinaison de matériel et de logiciel. Dans divers modes de réalisation, le système de simulation 800 comprend un ou plusieurs processeurs 810 (par ex., une unité centrale de traitement conventionnel (UC), une unité de traitement graphique, ou autre) configurés pour exécuter les programmes logiciels stockés dans la mémoire 812 (qui peuvent, par ex., être une mémoire RAM, une mémoire ROM ou une mémoire flash, etc.). En outre, le système de simulation 800 peut comprendre un écran 814, un ou plusieurs dispositifs de saisie par l'utilisateur 815 (tels que, par ex., un clavier, une souris, etc.), des dispositifs de stockage de données permanents 816 (tel que, par ex., un disque dur, un disque, etc.) et une interface réseau 817 qui facilite la communication avec l'unité de contrôle et de traitement 808 du système de mesure. Dans certains modes de réalisation, le système de simulation 800 reçoit des données provenant de (ou envoi des données vers) l'unité de contrôle et de traitement 808 à travers l'Internet, un réseau local ou un autre réseau. Dans certains modes de réalisation, les données provenant d'un système (par ex., les données de mesure provenant de l'unité de contrôle et de traitement 808) sont stockées sur un support lisible par ordinateur, et sont ensuite lues à partir de ce support par un autre système (par ex., le système de simulation 800). Par ailleurs, l'unité de contrôle et de traitement 808 et le système de simulation 800 peuvent être intégrés à un système informatique unique, par ex., sous forme de différents programmes logiciels s'exécutant sur le même ordinateur polyvalent. Les programmes informatiques peuvent être implémentés dans l'un quelconque des divers langages de programmation, par ex., sans limitation, C, C++, Object C, Pascal, Basic, Fortran, Matlab et Python.

[38] Les programmes logiciels du système de simulation 800 comprennent des instructions exécutables par le processeur implémentant les procédés de calcul décrits ici (par ex., le procédé de la figure 7), en se basant sur les entrées concernant la formation et la géométrie d'acquisition. Ces instructions peuvent être organisées sous forme de modules qui implémentent certaines fonctionnalités discrètes, telles que, par ex., : un module de définition de grille 820 qui détermine les zones horizontales en se basant sur une entrée du modèle de vitesse verticale et calcule l'espacement de la grille et la taille de la grille pour chaque zone ; un module de calcul du coefficient 822 qui calcule des coefficients de différence finie à partir de l'espacement de la grille et les stockent dans une mémoire 812 ; un module d'implémentation de la grille 824 qui attribue de la mémoire pour les variables stockées au niveau des points de la grille et de la sous-grille de la grille décalée de Lebedev ; et un module de simulation 826 qui implémente une boucle temporelle pour la mise à jour itérative des variables de champ. Les instructions peuvent aussi comprendre un module graphique 828 permettant d'afficher visuellement les ondes élastiques simulées, par ex., sur un écran 814 et/ou un module de comparaison 830 qui reçoit des entrées concernant les résultats de mesure provenant de l'unité de contrôle et de traitement 808, compare les mesures avec les résultats des simulations et présente les résultats de la comparaison à l'utilisateur (qui peut ensuite décider de réviser ou non le modèle de la formation), peut être sous la forme d'un affichage sur un écran 814 ou une impression d'une copie papier, et/ou la mise à jour automatique du modèle de formation conformément à un algorithme programmé. D'autres modules implémentant des fonctionnalités additionnelles peuvent être prévus. En outre, comme il sera facilement compris par des hommes de métier, la fonctionnalité globale permise par le système de simulation 800 peut être organisée et regroupée de plusieurs façons différentes (par ex., plus ou moins de modules ou des modules différents que ceux illustrés), ou implémentée dans son intégralité ou en partie avec des modules matériels dédiés à la place des modules logiciels. Ainsi, le mode de réalisation illustré n'est qu'un exemple illustratif parmi d'autres.

[39] En résumé, l'utilisation des modes de réalisation décrits ici peut entraîner une réduction substantielle des ressources informatiques, améliorant les opérations et le fonctionnement de l'ordinateur lui-même : beaucoup moins de mémoire peut être utilisée, et le temps de calcul pour des grands sous-ensembles de données peut être grandement réduit. Ainsi, la valeur des services fournis par une compagnie d'exploitation/d'exploration peut être augmentée d'un degré important.

[40] Même si des modes de réalisation spécifiques ont été illustrés et décrits ici, on doit apprécier que tout agencement configuré pour atteindre le même objectif peut être substitué par les modes de réalisation spécifiques illustrés. Cette divulgation est destinée à couvrir une quelconque et toutes les adaptations ou variations des divers modes de réalisation. Des combinaisons des modes de réalisation susmentionnés, et d'autres modes de réalisation qui n'ont pas été décrits ici, seront apparentes aux spécialistes du domaine lors de la lecture de la description précédente.

Claims

REVENDICATIONS

1. Procédé caractérisé en ce qu'il comprend : la définition de positions d'au moins une source sismique (804) et d'au moins un récepteur sismique (806) par rapport à un milieu isotrope transversalement incliné (TTI) ; la définition d'une grille décalée adaptative de Lebedev sur au moins une partie du milieu TTI, la grille comprenant une pluralité de zones horizontales (400, 402, 502, 510) avec des espacements de grille (dzl, dz2) associés, un espacement de grille associé à au moins une des zones étant différent d'un espacement de grille associé à une autre des zones ; et l'utilisation d'un processeur (810), calculant la propagation d'une onde sismique émise par l'au moins une source sismique (804) à travers le milieu TTI et la réception de celle-ci au niveau du récepteur (806) en résolvant un jeu d'équations d'onde élastique discrétisées sur la grille décalée adaptative de Lebedev.
2. Procédé selon la revendication 1, dans lequel le calcul de la propagation de l'onde sismique est basé au moins en partie sur un modèle de calcul du milieu TTI comprenant un ou plusieurs paramètres de champ.
3. Procédé selon la revendication 2, comprenant en outre : l'émission d'une onde sismique avec l'au moins une source sismique (804) ; la mesure d'une onde sismique résultant de l'onde sismique émise avec l'au moins un récepteur sismique (806) ; la comparaison de l'onde sismique mesurée avec l'onde sismique calculée au niveau du récepteur (806) et, si une différence entre les deux dépasse un seuil défini, l'ajustement du modèle de calcul du milieu TTI en ajustant un paramètre de champ du modèle.
4. Procédé selon la revendication 2 ou 3, dans lequel les paramètres de champ comprennent une vitesse d'onde de pression, une vitesse d'onde de cisaillement et les paramètres inélastiques de Thomsen.
5. Procédé selon l'une quelconque des revendications 2 à 4, dans lequel les paramètres de champ comprennent un angle d'inclinaison (Θ) et un angle d'azimut (φ) associés à une limite géologique inclinée (202) du milieu TTI.
6. Procédé selon l'une quelconque des revendications 1 à 5, dans lequel la définition de la grille décalée adaptative de Lebedev comprend la détermination des zones horizontales (400, 402, 502, 510) et des espacements de grille (dzl, dz2) associés à celles-ci basé au moins en partie sur le modèle de vitesse d'onde de cisaillement spécifiant une vitesse d'onde de cisaillement verticalement variable.
7. Procédé selon l'une quelconque des revendications 1 à 6, dans lequel les équations d'onde élastique discrétisées sur la grille comprennent des équations de différence finie d'au moins du deuxième ordre, par exemple d'au moins du quatrième ordre, en particulier d'au moins du seizième ordre.
8. Procédé selon la revendication 7, dans lequel les équations de différence finie comprennent des coefficients de différence finie, le procédé comprenant en outre le calcul des coefficients de la différence finie sur la base des espacements de grille (dzl, dz2).
9. Procédé selon la revendication 8, dans lequel les coefficients de différence finie sont calculés en utilisant au moins une décomposition en valeurs singulières ou au moins une optimisation de moindre carré.
10. Procédé selon la revendication 8 ou 9, dans lequel les coefficients de différence finie sont variables à l'intérieur d'une région de chevauchement (406) comprenant des parties adjacentes de deux des zones horizontales (400, 402, 502, 510), et dans lequel les coefficients de différence finie sont constants à l'intérieur de chacune des zones horizontales (400, 402, 502, 510) à l'extérieur de la région de chevauchement (406).
11. Procédé selon l'une quelconque des revendications 7 à 9, dans lequel le calcul de la propagation de l'onde sismique comprend la mise à jour de variables de champ utilisant les équations de différence finie, et dans lequel la mise à jour des variables de champ dans une région de chevauchement (406) qui comprend des parties adjacentes de deux des zones horizontales (400, 402, 502, 510) est basé sur des variables de champ dans chacune desdites deux zones horizontales.
12. Procédé selon la revendication 11, dans lequel la mise à jour de la variable de champ dans la région de chevauchement (406) comprend la création de points de grille imaginaires en interpolant entre des points de grille réels.
13. Procédé selon l'une quelconque des revendications 1 à 12, dans lequel la définition des positions de l'au moins une source (804) et de l'au moins un récepteur (806) comprend la spécification d'un type d'acquisition.
14. Procédé selon l'une quelconque des revendications 1 à 13, dans lequel le calcul de la propagation de l'onde sismique comprend l'application de conditions de limite de couche à correspondance parfaite convolutionnelle.
15. Procédé selon l'une quelconque des revendications 1 à 14, comprenant en outre : l'affichage, utilisant un écran (814) ou un appareil d'impression, de la propagation de l'onde sismique émise par l'au moins une source sismique (804) à travers le milieu TTI.
16. Système caractérisé en ce qu'il comprend : au moins une source sismique (804) pour l'émission d'une onde sismique dans un milieu TTI ; au moins un récepteur sismique (806) configuré pour détecter une onde sismique se propageant à travers le milieu TTI ; et un dispositif de calcul configuré pour recevoir des informations sur les positions de l'au moins une source sismique (804) et l'au moins un récepteur sismique (806) par rapport à la formation rocheuse, définir une grille décalée adaptative de Lebedev sur au moins une partie de la formation rocheuse, la grille comprenant une pluralité de zones horizontales (400, 402, 502, 510) avec des espacements de grille (dzl, dz2) associés, un espacement de grille associé à au moins l'une des zones différant d'un espacement de grille associé à une autre des zones, basé au moins en partie sur le modèle de calcul du milieu TTI, le calcul de la propagation de l'onde sismique émise par l'au moins une source sismique (804) à travers la formation et la détection de celle-ci au niveau du récepteur (806) en résolvant un jeu d'équations d'onde élastique discrétisées sur la gille décalée adaptative de Lebedev, et la comparaison de l'onde sismique détectée avec l'onde sismique calculée au niveau du récepteur (806) et, lorsqu'une différence entre les ondes sismiques détectée et calculée dépasse un seuil défini, l'ajustement du modèle de calcul en ajustant au moins un paramètre de champ de celui-ci.
17. Système selon la revendication 16, dans lequel le dispositif de calcul est configuré en outre pour déterminer les zones horizontales (400, 402, 502, 510) et les espacements de grille (dzl, dz2) associés à celles-ci basé au moins en partie sur un modèle de vitesse d'onde de cisaillement spécifiant une vitesse d'onde de cisaillement verticalement variable.
18. Système selon la revendication 16 ou 17, dans lequel les équations d'onde élastique discrétisées comprennent des équations de différence finie d'au moins du seizième ordre.