FR 2854245 A1 20041029 FR 0305159 A 20030428 PROCEDE DE TRAITEMENT DE CUBES SISMIQUES CORRESPONDANT A DES ACQUISITIONS REALISEES POUR UNE MEME ZONE A DES MOMENTS DIFFERENTS. DOMAINE GENERAL ET ETAT DE LA TECHNIQUE La présente invention est relative au traitement de traces sismiques correspondant à des acquisitions réalisées pour une même zone à 10 moments différents. Il est commun, dans le domaine de la sismique, de mettre en évidence des évolutions de caractéristiques géologiques ou géophysiques d'une zone donnée de la croûte terrestre, en comparant des cubes sismiques correspondant à des acquisitions réalisées à deux moments 15 différents (sismique 4D ou Time-lapse). Par " cube sismique ", on entend ici et dans tout le présent texte, toute collection de traces sismiques dans laquelle chaque trace sismique est un signal d'amplitude en fonction du temps et est associée à une position donnée sur la zone au sol sur laquelle l'acquisition a été réalisée. Généralement, ainsi que l'illustre la figure 1, la mise en évidence des variations est faite en soustrayant deux cubes de données sismiques correspondant aux deux moments d'acquisition. C'est notamment cette méthode qui est utilisée pour mettre en évidence des évolutions dans les gisements d'hydrocarbures en production. 25 Toutefois, lorsque la différence en terme de signal sismique est de faible amplitude, les propriétés physiques que l'on cherche à mettre en évidence peuvent être facilement masquées par des artefacts comme par exemple le bruit dû à la géométrie imparfaite de l'acquisition des données sismiques. Le caractère aléatoire et non répétable du bruit de chaque acquisition fait que, comme l'illustre la figure 2, la simple soustraction des données sismiques augmente le niveau de bruit par rapport au signal recherché. Sur cette figure 2, le bruit est schématisé par des traits horizontaux (cube A) et des traits verticaux (cube B) qui se superposent sur le cube soustrait (cube A-B). De plus la superposition des bruits respectifs détruit leurs organisations spatiales statistiques (si elles existent) et donc rend plus difficile leurs filtrages. PRESENTATION DE L'INVENTION Un but de l'invention est de proposer une nouvelle technique qui permet de pallier les inconvénients des techniques antérieures. A la différence de la simple soustraction, l'invention propose de prendre comme cube de référence non pas un des cubes existant, c'est-àdire le cube enregistré initialement dans le temps pour la sismique 4D, mais la partie commune des deux cubes (ou plus) qui par définition est invariante. Ainsi, l'invention propose un procédé de traitement d'un ensemble de cubes sismiques correspondant à une même zone d'acquisition et à différents moments d'acquisition, caractérisé en ce qu'on détermine au moins un cube sismique qui est une estimation de la composante commune à au moins deux cubes sismiques correspondant à des moments 20 d'acquisition et on soustrait cette estimation de composante commune d'au moins l'un des deux cubes initiaux. De cette façon, on dispose de cubes différentiels qui ont été déterminés sans que les bruits des deux cubes initiaux ne s'y retrouvent superposés. Et ces cubes différentiels peuvent être facilement filtrés, puisque l'organisation spatiale statistique n' y a pas été détruite par la superposition de bruits. DESCRIPTION DES DESSINS - la figure 1 est une représentation schématique qui illustre la mise en évidence de l'évolution d'un volume d'hydrocarbure par soustraction des cubes sismiques enregistrés à différentes étapes de production; sur les cubes, la zone horizontale séparant les deux cubes symbolise le contact entre l'hydrocarbure (partie haute) et l'eau (partie basse); - la figure 2 est une représentation schématique qui illustre le fait que la soustraction des cubes bruités ne fait qu'augmenter le niveau de bruit par rapport aux variations que l'on met en évidence; - la figure 3 est une représentation schématique qui illustre la détermination d'un cube commun et la détermination, à partir de ce cube commun, de cubes de différences; - la figure 4 illustre un mode de mise en oeuvre possible pour la 10 détermination d'une estimation de la partie commune à deux cubes sismiques; - les figures 5 et 6 illustrent pour un instant d'échantillonnage dans le temps donné, des exemples de cubes commun et différentiels obtenus en mettant en oeuvre le procédé illustré sur les figures 3 et 15 4; - la figure 7 illustre un exemple de masque appliqué aux cartographies 2D; - la figure 8 illustre un autre mode de mise en oeuvre possible avec plusieurs itérations. - la figure 9 illustre un autre mode de mise en oeuvre possible pour calculer deux cubes avec une phase commune. DESCRIPTION D'UN OU PLUSIEURS EXEMPLES DE MISE EN OEUVRE Exemple général On considère deux cubes sismiques (données A et données B) correspondant à deux acquisitions réalisées, pour une même zone au sol, à deux moments différents. Les données peuvent être décomposées ainsi: données A = partie commune + bruit A données B = partie commune + bruit B + A o A représente les variations des données sismiques et o la partie commune se définit comme la partie cohérente ou invariante des données A et B. Dans une première étape, on détermine une estimation de cette partie commune. Différentes méthodes peuvent être utilisées à cet effet. Notamment, cette estimation peut être calculée par co-krigeage, ou encore par cross-corrélation ou par calcul de valeurs moyennes. Cette partie commune estimée représente la géologie qui n'a pas changée avec le temps. Elle a l'avantage de présenter un meilleur rapport signal à bruit que les données initiales. Dans un deuxième temps on détermine les deux différences partielles suivantes: Diff A = données A partie commune estimée 15 Diff B = partie commune estimée - données B La figure 3 illustre le cube de composante commune et les cubes de différences partielles (cubes différentiels) ainsi obtenus. Les deux différences partielles permettent de mettre en évidence les évolutions de la géologie. Elles vérifient: Diff A= bruit A - a A Diff B= -bruit B - (1 -c)A O a est le rapport de répartition du signal résiduel, lequel est souvent proche de 0.5 mais peut varier spatialement. Le rapport Signal/Bruit de ces deux différences partielles est alors ac.A / bruit A pour Diff A et (1-OE)A I bruit B pour Diff B. On met ensuite en oeuvre sur les différences partielles ainsi obtenues un traitement de filtrage déterministe, ce qui est possible grâce au fait que les organisations spatiales statistiques sont conservées et en 30 particulier n'ont pas été supprimées par la superposition de bruits et que le bruit de chaque jeux de données conserve ses propriétés intrinsèques initiales, comme sa cohérence spatiale statistique. Le filtrage indepedant pour chaque différence partielle peut etre mis en oeuvre par exemple par un filtre spatial deterministique dans le domain KxKy ou un filtre statistique par krigeage factoriel ou meme une fonction de lissage suivant un axe particulier permettant d'atténuer les bruits de l'acquisition. L'orientation 5 spatial du bruit est d'autant plus facile a déterminer par des mesures statistiques (variogrammes 2D/3D) que la partie commune contenant la majorite du signal sismique de la geologie a ete soustraite. Après filtrage, il est possible de reconstituer la différence totale à partir des deux différences partielles filtrées ainsi obtenues: données A - données B = Diff A + Diff B On comprend que la sommation des différences partielles après filtrage de celles-ci permet un meilleur rapport signal à bruit que ce qui est obtenu en calculant la différence totale par soustraction des données A et des données B. On notera qu'il n'y a pas de perte d'information, puisque la somme des différences partielles est égale à la différence totale. Egalement, on notera que, si le niveau de bruit est différent entre les données A et B initiales, une des deux différences partielles présente un meilleur rapport signal sur bruit que la différence totale. Dans certains cas, l'information sur l'évolution de la géologie ou géophysique du sous-sol sera plus apparente sur un cube correspondant à une différence partielle que sur le cube correspondant à la différence totale. Exemple de calcul de partie commune Principe La figure 4 illustre un exemple de calcul de la partie commune utilisant le Co-krigeage factoriel automatique spatial (ou analyse krigeante 30 multi-variables automatique) dans le domaine des fréquences. Les données temporelles initiales sont transformées par FFT en cubes de parties réelles et de parties imaginaires dans le domaine fréquentiel (domaine de Fourier). La partie commune de chaque plan fréquentiel réel et imaginaire est calculée par la méthode dite " co-krigeage factoriel automatique " avec un opérateur 2D dans le plan (x, y). Par le processus de transformée de Fourier inverse du cube 5 commun des parties réelles et du cube commun des parties imaginaires, on obtient alors le cube commun temporel des données sismiques. Exemple détaillé Un exemple plus détaillé de traitement est donné ci-dessous: 1) Dans une première étape, on lit les données qui correspondent aux deux collections de traces (cube A et cube B) à traiter. Ces données sont enregistrées dans un fichier o chaque trace est identifiée par un numéro de ligne et un numéro de colonne qui correspondent à une position au sol. 2) Dans une deuxième étape, on sélectionne une fenêtre de temps pour définir sur les traces sismiques l'intervalle de temps pris en considération pour le traitement.. On complète le cas échéant les portions de traces ainsi sélectionnées par des valeurs nulles pour obtenir pour chaque trace un nombre d'échantillons correspondant à une puissance de 20 2, 3, 5, etc. pour, selon les fonctions de FFT utilisées, permettre la FFT. 3) Dans une troisième étape, on applique à la succession d'échantillons correspondant à chacune des traces et à la fenêtre temporelle sélectionnée une transformée de Fourier rapide (FFT). On obtient alors au total quatre cubes en fréquence correspondant 25 à la partie réelle et la partie imaginaire de chaque sélection des cubes A et B. 4) Pour chaque plan de fréquence, on construit alors, dans le domaine réel d'une part et dans le domaine imaginaire d'autre part, deux cartes (2D) correspondant l'une aux données du cube A et l'autre aux 30 données du cube B. Puis on applique sur les cartes ainsi obtenues pour les cubes de données A et B un traitement de co-krigeage factoriel automatique afin d'en déduire une carte commune aux deux cartes initiales. Des rappels sur le cokrigeage factoriel sont donnés en Annexe I à la fin de la description. Le traitement de co-krigeage factoriel est réalisé de façon indépendante d'une part avec les partiesréelles et d'autre part avec les parties imaginaires. Ce traitement comporte par exemple, dans chacun de ces deux domaines, les étapes suivantes: - En fonction d'un rayon de variogramme donné, calcul d'un variogramme 2D pour la carte correspondant aux données A et pour la 10 carte correspondant aux données B; calcul d'un cross-variogramme, . Construction des matrices de covariances et cross-covariance à partir des 3 variogrammes ainsi obtenus. - Résolution du système de co-krigeage factoriel qui correspond à ces matrices, afin d'en déduire un opérateur spatial, c'est à dire fonction de 15 la position en ligne et en colonne de l'échantillon considéré sur la cartographie. Pour stabiliser l'inversion de la matrice du système, on rajoute sur la diagonale des matrices de covariance une constante qui correspond à un tres petit pourcentage de la racine carrée du produit des deux variances. - On convolue alors l'opérateur spatial ainsi obtenu avec les cartes initiales pour obtenir la partie commune des deux cartes. 5) Le traitement qui vient d'être décrit est répété pour chaque plan fréquence des parties réelles et imaginaires. On obtient deux cubes, un pour les parties réelles communes et l'autre pour les parties imaginaires communes. 6) Une fois tous les plans fréquences calculés, on met en oeuvre une transformée de Fourier inverse des cubes complexes communs pour obtenir un cube en temps de sismique commune. Exemple de résultats Les figures 5 et 6 illustrent des résultats obtenus avec ce traitement. Les données A et B ont été enregistrées à plusieurs années d'intervalle. La figure 5 montre que la partie commune apparaît avec un meilleur rapport signal sur bruit que les données initiales A et B. Elle représente la géologie invariante dans le temps, celle-ci est prise comme référence. La différence totale des deux jeux de données A et B est découplée en deux 5 différences partielles par rapport à leur partie commune. On observe que les variations de la réponse sismique des événements géologiques dans le temps (signature 4D) ressortent nettement plus sur la différence partielle des données A que sur la différence totale (Figure 6). Ceci est du au fait que le bruit visible sur la différence totale est reporté sur la différence 10 partielle des données B. Exemples de variantes de mise en oeuvre D'autres variantes de réalisation sont bien entendu envisageables. Masque de traitement ainsi que l'illustre la figure 7, le calcul de l'opérateur peut être limité à des régions précises correspondant à des données valides. Les variogrammes et cross-variogramme sont ainsi calculés sur les régions définies par le masque. Le masque peut par exemple cacher une zone sans 20 données due à une plate-forme de forage ou des données très bruitées par des effets de bords. L'utilisation d'un masque va améliorer la qualité de l'opérateur et donc du filtrage de bruit. Le masque est défini avec des valeurs valides en blanc (1) et des valeurs non valides en noir (0). Double itération il est parfois avantageux de faire le calcul du cube commun en deux itérations: (Figure 8) La première itération consiste à calculer un spectre commun d'amplitude. 30. A cet effet, pour chaque plan de fréquence, on détermine à partir des cartes de parties réelles et imaginaires, une cartographie d'amplitude, ainsi qu'une cartographie de phase. Puis, à partir des cartographies d'amplitude ainsi déterminées pour les données A et les données B pour chaque fréquence, on détermine une cartographie commune d'amplitude spectrale. Cette cartographie est de façon appropriée déterminée par cokrigeage factoriel automatique. La partie réelle et imaginaire de A et B est ensuite reformulée avec 5 l'amplitude commune ainsi obtenue et la phase initiale de A et B. Après un traitement de FFT inverse des cubes, on obtient deux cubes de sismique A' et B' ayant le même contenu fréquentiel. Ces deux cubes sont ensuite traités avec une seconde itération. La seconde itération calcule la partie commune réelle et imaginaire 10 de A' et B'. En fait, elle va fixer la phase commune des deux cubes sismiques A' et B' puisque le spectre d'amplitude est déjà commun après la première itération. On détermine ensuite la partie commune, de la même façon que ce qui a été précédemment décrit en référence à la figure 6. Une autre variante encore possible de mise en oeuvre est illustrée sur la figure 9. Dans cette variante de mise en oeuvre, on extrait des cubes d'amplitude associés respectivement au cube A et au cube B. En parallèle ou dans une deuxième étape, on sort des parties réelles 20 et imaginaires des cubes A et B, des cubes correspondant au cosinus et au sinus de la phase des données. On détermine ensuite des cubes correspondant l'un à la partie commune des cubes de cosinus, l'autre à la partie commune des cubes de sinus. Puis on détermine un cube dont les données sont le ratio des données du cube de sinus et du cube de cosinus ainsi obtenus: on dispose alors d'un cube donnant la partie commune de la phase. Ce cube de phase commune est combiné avec les deux cubes d'amplitude initialement extraits des cubes A et B pour en déduire des 30 cubes filtrés A' et B'. Ces deux cubes de données filtrées peuvent eux-mêmes faire l'objet de la seconde itération décrite en référence à la figure 5. En variante encore, il est possible de calculer le cube commun aux cubes de données A et B en combinant d'une part le cube de partie commune de phase et d'autre part le cube de partie commune d'amplitude. On notera que dans les exemples qui précèdent, on a calculé la 5 partie commune de chaque plan fréquentiel par une méthode géostatistique appelée co-krigeage factoriel automatique (ou analyse krigeante multivariables automatique). Autres traitements que le co-krigeage factoriel 2D D'autres traitements d'estimation de partie commune, et notamment de traitements mettant en oeuvre un co-krigeage factoriel, sont envisageables. On peut notamment déterminer une partie commune directement dans le domaine temps sur une succession de cartes (x,y) d'amplitude 15 sismique (ou d'autres données). Pour chaque instant considéré, on détermine alors la partie commune aux deux cartographies qui correspondent à cet instant en mettant en oeuvre un traitement de co-krigeage automatique factoriel. Egalement, le traitement peut utiliser des opérateurs 3D au lieu 20 d'opérateurs 2D. Le traitement de co-kriegeage automatique factoriel se généralise en effet avec une matrice de co-variance à trois dimensions. Il est à noter que l'utilisation d'opérateurs 2D (x,y) (ou 3D le cas échéant) rend la méthode très performante pour filtrer des bruits spatialement organisés. Bien entendu, également, il est possible de déterminer des parties communes à plus de deux cubes, en calculant par exemple des parties communes de parties communes. De cette façon, le procédé proposé peut être étendu à un nombre de cubes supérieurs à deux. Annexe I - Rappel sur le co-krigeage factoriel On décrit ci-dessous un exemple d'obtention d'une cartographie 5 commune par co-krigeage factoriel. On désigne par Z1 et Z2 les deux fonctions correspondant à ces deux cartographies. Dans un premier temps on calcule le variogramme croisé de ces deux fonctions, dont les valeurs sont: Y12(h)= - % (ZI (x)- ZI (x+h))(Z2(x) Z2(x+h)) N o x et x+h désignent les couples de points pris en considération selon la direction et pour la distance h pour lesquels on détermine la valeur du variogramme, et o N est le nombre de ces couples de points pour cette direction et cette distance. Connaissant ce variogramme croisé, on détermine ensuite une estimation de la fonction qui lui correspond, laquelle vérifie: Z 2(X) = E ZI + E 2 Z2 Z*12 a a fl2f2 o a et f3 sont deux indices muets désignant les points considérés autour du point x pour lequel on cherche à déterminer une estimation de ladite 20 fonction, Z' et Z2 étant la valeur audit point x, N étant le nombre de ces a /2 points et o l' et Sz sont des coefficients de pondération. Ces coefficients de pondération X et 2X se déterminent par inversion de l'équation de co-krigeage: C"ll1... C CN1211 C2211... C211 21X C,,1N... C, ,NN C,211... C,2NN _IAN C,2NX C2,11.. C2,N1 C2211... C22iN1 221 C,2X C211N... C21NN C2211... C22NN 22N C,2NX o les coefficients C12a3 et C21ac sont les valeurs de variance croisée des fonctions ZI et Z2 aux points correspondants aux indices a et f3 o les coefficients Cl lap et C22ap sont les valeurs de covariance 5 respectivement de la fonction Zi et de la fonction Z2 auxdits points. L'indice X correspond au point désigné par x précédemment. On notera que la matrice qui apparaît dans cette équation a l'avantage d'être inversible sous certaines conditions de calcul. De cette façon, en utilisant les covariances expérimentales, on 10 décompose de façon automatique les deux variables correspondant aux deux séries de données initiales en une composante commune et deux composantes résiduelles orthogonales. La régularité des données fait que la covariance expérimentale est connue pour toutes les distances utilisées. aucune interpolation n'est nécessaire et la matrice est définie positive La fonction alors obtenue est une estimation de la composante commune aux deux séries de données que constituent les deux cartographies. REVENDICATIONS 1. Procédé de traitement d'un ensemble de cubes sismiques correspondant à une même zone d'acquisition et à différents moments 5 d'acquisition, caractérisé en ce qu'on détermine au moins un cube sismique qui est une estimation de la composante commune à au moins deux cubes sismiques correspondant à des moments d'acquisition différents et on soustrait cette estimation de composante commune d'au moins l'un des deux cubes initiaux. 2. Procédé selon la revendication 1, caractérisé en ce qu'on met en oeuvre un traitement de filtrage sur un cube différentiel ainsi obtenu. 3. Procédé selon les revendications 1 ou 2, caractérisé en ce qu'on soustrait à l'un et à l'autre des deux cubes initiaux le cube qui correspond à l'estimation de leur composante commune, en ce qu'on met en oeuvre un 15 traitement de filtrage sur l'un et l'autre des deux cubes différentiels ainsi obtenus et en ce qu'on soustrait ou somme ces cubes différentiels filtrés pour reconstituer une estimation de la différence entre les deux cubes initiaux. 4. Procédé selon l'une des revendications précédentes, caractérisé 20 en ce que pour déterminer un cube qui est une estimation d'une composante commune à au moins deux cubes sismiques: - on convertit chacun de ces deux cubes sismiques dans le domaine fréquentiel pour obtenir pour chacun d'eux un cube réel et un cube imaginaire, - on calcule un cube sismique qui est une estimation de la composante commune aux cubes réels ainsi obtenus et un cube qui est une estimation de la composante commune aux deux cubes imaginaires ainsi obtenus, et - on convertit dans le domaine temporel les deux cubes, l'un réel, l'autre imaginaire, ainsi obtenus. 5. Procédé selon l'une des revendications précédentes, caractérisé en ce que, pour déterminer un cube qui est une estimation d'une composante commune à au moins deux cubes: - on convertit chacun de ces cubes dans le domaine fréquentiel pour en déduire pour chacun d'eux un cube d'amplitude et un cube de phase, - on calcule une estimation de la partie commune aux cubes d'amplitude ainsi obtenus. 6. Procédé selon la revendication 5, caractérisé en ce que - on convertit dans le domaine temporel le cube d'amplitude ainsi obtenu en l'associant respectivement aux différents 10 cubes de phase, - on convertit dans le domaine fréquentiel chacun des cubes ainsi obtenus dans le domaine temporel pour obtenir pour chacun de ces cubes un cube réel et un cube imaginaire, - on calcule un cube qui est une estimation de la composante 15 commune aux cubes réels et un cube qui est une estimation de la composante commune aux cubes imaginaires, et - on convertit dans le domaine temporel les deux cubes, l'un réel, l'autre imaginaire, ainsi obtenus. 7. Procédé selon l'une des revendications précédentes, caractérisé en ce que pour déterminer un cube qui est une estimation d'une composante commune à au moins deux cubes: - on convertit chacun de ces cubes dans le domaine fréquentiel 25 pour en déduire pour chacun d'eux un cube d'amplitude, ainsi qu'un cube correspondant à sa partie réelle et un cube correspondant à sa partie imaginaire, - on extrait des cubes correspondant à la partie réelle et des cubes correspondant à la partie imaginaire des cubes cosinus et 30 sinus, - on calcule une estimation de la partie commune aux cubes de cosinus et une estimation des cubes de sinus, - on déduit des cubes de parties communes ainsi calculés un cube de données de phase communes, 8. Procédé selon la revendication 7, caractérisé en ce que: - on convertit dans le domaine fréquentiel chacun des cubes ainsi obtenus dans le domaine temporel pour obtenir pour chacun de ces cubes un cube réel et un cube imaginaire, - on calcule un cube qui est une estimation de la composante commune aux cubes réels et un cube qui est une estimation 10 de la composante commune aux cubes imaginaires, et - on convertit dans le domaine temporel les deux cubes, l'un réel, l'autre imaginaire, ainsi obtenus 9. Procédé selon l'une des revendications 4 à 8, caractérisé en ce 15 que pour calculer un cube qui est une estimation d'une composante commune à au moins deux cubes sismiques, on met en oeuvre un traitement de co- krigeage factoriel. 10. Procédé selon la revendication 9, caractérisé en ce que pour calculer un cube qui est une estimation d'une composante commune à au 20 moins deux cubes sismiques: - on échantillonne chacune de ces cubes sissmiques pour en déduire pour chacune d'eux une succession de cartographie 2D - on détermine pour chaque couple de cartographies 2D ainsi obtenu des matrices de co-variance, correspondant aux variogrammes de l'une et 25 l'autre de ces deux cartographies, ainsi qu'à leur variogramme croisé, - on résout un système de co-krigeage factoriel qui correspond à ces matrices pour en déduire un opérateur spatial à appliquer aux deux cartographies 2D qui leur correspondent, afin d'en déduire une cartographie 2D qui correspond à leur composante commune, l'ensemble des 30 cartographies 2D ainsi obtenues permettant de reconstituer le cube sismique commun. 11. Procédé selon la revendication 10, caractérisé en ce que pour déterminer le système de co-krigeage à résoudre pour déduire l'opérateur spatial à appliquer à deux cartographies 2D, on applique aux données d'au moins l'une des deux cartographies, préalablement à la détermination des co-variances et des cross-variances, un masque destinée à annuler les données de ladite cartographie sur une certaine zone. 12. Procédé pour la mise en évidence d'évolutions de caractéristiques géologiques ou géophysiques sur des réservoirs d'hydrocarbure, dans lequel on met en oeuvre des acquisitions sismiques à plusieurs moments différents et on traite des cubes sismiques correspondant à ces acquisitions pour mettre en évidence les éventuelles 10 évolutions, caractérisé en ce qu'on met en oeuvre à cet effet un procédé de traitement selon l'une des revendications précédentes.