These Audrey NEAU 2009 PDF
These Audrey NEAU 2009 PDF
These Audrey NEAU 2009 PDF
THÈSE
Présentée à
Audrey NEAU
Docteur
Spécialité : Exploration géophysique
Madame Frédérique Fournier et Monsieur Jean-Paul Chilès ont accepté d'être les rap-
porteurs de cette thèse, et je les en remercie. Ils ont également contribué par leurs nom-
breuses remarques et suggestions à améliorer la qualité de ce mémoire, et je leur en suis
très reconnaissante.
Messieurs Philippe Doyen, Hervé Perroud et Mrinal Sen m'ont fait l'honneur de parti-
ciper au Jury de soutenance ; je les en remercie profondément. Remerciements particuliers
à Monsieur Hervé Perroud pour avoir accepté le rôle de président du jury de cette thèse,
la première issue de son projet de Master.
Je ne saurais jamais assez remercier Issam Tarrass pour sa bonne humeur quotidienne,
pour son soutien moral permanent et pour sa gentillesse. J'ai eu la chance de partager son
iv Remerciements
bureau pendant trois ans... un véritable ami à présent. Un grand merci à toi. Et un gros
M.... pour la n de ta thèse !
Je remercie tous ceux sans qui cette thèse ne serait pas ce qu'elle est : aussi bien par
les discussions que j'ai eu la chance d'avoir avec eux, leurs suggestions ou contributions.
Je pense ici en particulier à Messieurs Marc Elias et Eric Maocec pour le champ Alpha,
à Lorette Anquelle, Yoann Guilloux et Dimitri Modin pour le champ Bêta, et à Cécile
Pabian, Jean Luc Broyer, Alain Jegoux et Matthieu Pellerin pour le champ Gamma. Merci
aux collègues de CGGVeritas, Rémi, Raphael, Robert pour leur aide sur GeoSI. Je remercie
également Messieurs Jean Luc Piazza, Emmanuel Chavanne, Pierre Biver, et Léon Barens
pour les discussions et suggestions en géomodélisation. Un grand merci à Benoit Paternos-
ter pour ces discussions où l'on vient avec une question et dont on ressort la tête pleine de
nouveaux horizons...
Je passe ensuite une dédicace spéciale à tous les jeunes gens que j'ai eu le plaisir de
côtoyer durant ces quelques années à Pau, je pense notamment à Alain-Christophe, Ma-
thieu, Livinus, Maelle, Khalil, Ravi, Paritosh, Désirée, Céline, et tous les étudiants des
promotions 2008 et 2009 du Master Génie Pétrolier de Pau, à qui j'ai dispensés quelques
cours de caractérisation réservoir.
Mes remerciements vont également vers François Levêque pour m'avoir lancée sur cette
voie, pour avoir cru en ma motivation et pour s'être décarcassé pour que je puisse assister
à tous les camps de terrain. Je ne saurai jamais assez remercier Guy Sénéchal pour ses
encouragements précieux qui m'ont permis de garder la tête hors de l'eau quand les choses
étaient diciles en Licence et pour avoir partagé le projet Rustrel avec moi deux années
de suite.
Enn, mais pas des moindres, un grand merci à mes amis, mes trois petits rayons de
soleil, Florence, Aurélie et Julia. Et aussi à Vicky pour m'avoir supporté au quotidien, dans
les moments diciles. Un grand merci à Anne et Bernard pour m'avoir "adoptée" et pour
avoir été présents à tout instant. Petit clin d'oeil à Annick et Nicole pour avoir suivi mon
parcours et pour votre présence le jour J.
Ces remerciements ne sauraient être complets si ma famille n'était citée. J'ai gardé le
meilleur pour la n : un grand merci à toute ma famille qui m'a toujours fait conance
et soutenue dans mes choix. Tout au long de mes études, vous avez su m'encourager et
me remonter le moral dans les moments diciles, et m'avez permis de réaliser mes projets...
A tous, merci.
Remerciements v
Résumé :
La caractérisation sismique des réservoirs pétroliers nécessite l'intégration de plusieurs
techniques telles que la lithosismique, la géomodélisation, la géostatistique, l'utilisation
des algorithmes évolutionnaires et la pétrophysique. L'information sismique est d'abord
utilisée pour la description de l'architecture externe des réservoirs car son utilisation pour
la description des faciès ne se fait pas sans dicultés. L'objectif de cette thèse est d'appor-
ter des outils nouveaux pour aider à l'utilisation de l'information sismique pour caractériser
les réservoirs.
Un premier travail a consisté à évaluer l'impact des incertitudes structurales sur les inver-
sions pétroélastiques et les conséquences en terme de classication de faciès. Ensuite, nous
considérons la modélisation sismique comme aide à l'évaluation du modèle réservoir. Cette
modélisation permettra de faire le lien entre les simulateurs réservoir ou les géomodeleurs
et la réponse sismique du réservoir.
Nous développons ensuite deux approches alternatives aux méthodes traditionnelles en in-
version pétroélastique et pétrophysique. La première utilise la méthode géostatistique des
déformations graduelles pour créer des réalisations de propriétés réservoirs. Elle permet
de créer des propriétés à l'échelle réservoir, conditionnées aux puits, tout en respectant
une fonction coût basée sur la comparaison des données sismiques réelles et issues de ces
réalisations.
La seconde méthode repose sur le principe de la classication supervisée et utilise des ré-
seaux de neurones pour analyser la forme des traces sismiques. Une première étape consiste
à générer un volume d'apprentissage contenant tous les modèles pétrophysiques envisa-
geables pour un champ donné. Ces modèles sont analysés par les réseaux de neurones. Les
neurones ainsi identiés sont appliqués aux données réelles, pour identier des relations
pétrophysique/sismique identiques aux données d'apprentissage.
Toutes les méthodologies sont validées sur plusieurs réservoirs choisis pour leurs particula-
rités géologiques (complexité structurale, lithologie du réservoir).
Abstract :
Seismic characterization of hydrocarbon reservoirs is based on various techniques : litho-
seismic, geomodeling, geostatistics, evolutionary algorithms, and petrophysics. Seismic in-
formation is rst used for the description of the reservoir structure, but then its relationship
with facies description is a dicult task. The aim of this thesis is to develop new tools for
seismic reservoir characterization.
A rst work has consisted in evaluating the impact of structural uncertainties on petroe-
lastic inversion and its consequences in terms of facies classication. Then, we consider
seismic modeling as an aid to reservoir model evaluation. This modeling step will make the
connection between the reservoir simulators (or geomodelers) and the seismic response of
the reservoir.
Then we develop two alternative approaches for petroelastic and petrophysical inversion.
The rst one uses the gradual deformation method to generate reservoir property realiza-
tions. This method generates properties at the reservoir scale, conditioned by the wells,
while respecting a cost function based on the comparison of actual and synthetic seismic
data.
The second method is based on supervised classication principle and uses neural networks
to analyze the waveform of seismic traces. A rst step is to generate a volume containing
all possible petrophysical models for the concerned eld. These models are analyzed by the
neural networks. The neurons identied are applied on the actual data to recognize similar
petrophysical/seismic relationships.
All methods are tested and validated on actual reservoirs, chosen for their specic features
(structural complexity, reservoir lithology).
Remerciements iii
Résumé vi
Abstract vii
Liste des gures xix
Liste des tableaux xx
1 Introduction générale 1
1.1 Problématique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Contexte scientique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.1 L'ingénierie du réservoir pétrolier . . . . . . . . . . . . . . . . . . . . 4
1.2.2 Les propriétés réservoirs . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 Organisation de la thèse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.2 Calcul des données sismiques synthétiques . . . . . . . . . . . . . . . . . . . 44
3.2.1 Le modèle convolutionnel 1D . . . . . . . . . . . . . . . . . . . . . . 44
3.2.2 Le tracé de rayons . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.2.3 Modélisation numérique de l'équation des ondes . . . . . . . . . . . . 45
3.3 Transformation en grille régulière . . . . . . . . . . . . . . . . . . . . . . . . 46
3.4 Evaluation des modélisations sur le champ Alpha . . . . . . . . . . . . . . . 48
3.4.1 Eets du changement de support Sismique-Réservoir-Voxet . . . . . 48
3.4.2 Sismogrammes synthétiques 1D . . . . . . . . . . . . . . . . . . . . . 50
3.4.3 Modélisation sismique à partir d'une grille réservoir . . . . . . . . . . 51
3.5 Lissage des propriétés du voxet . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.6 Comparaison des sismiques réelles et synthétiques . . . . . . . . . . . . . . . 55
3.7 Conclusions et perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
Annexes 173
Annexe 1 : SEG Annual Conference (San Antonio 2007) 173
Annexe 2 : EAGE Annual Conference (Rome 2008) 179
Annexe 3 : SEG Annual Conference (Las Vegas 2008) 185
Bibliographie 193
Table des gures
1.2.2 Construction du modèle réservoir. (a) Modèle structural ; (b) grille réser-
voir ; (c) remplissage du réservoir avec des propriétés (ici facies : en rouge,
les argiles ; autres couleurs, les sables avec diérents niveaux de cimenta-
tions). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2.1 La propagation d'une onde P sur une interface entre deux milieux induit
la génération d'une onde P rééchie Pr , d'une onde P transmise Pt , d'une
onde S rééchie Sr et d'une onde S transmise St . . . . . . . . . . . . . . . 16
2.2.2 Les diérentes approximations de l'équation de Zoeppritz. L'approxima-
tion d'Aki et Richards, en pointillé vert, est retenue pour sa validité dans
l'intervalle des angles considérés dans le domaine pétrolier. Tiré de Egre-
teau 2005. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.3.1 Localisation géographique du champ Alpha. . . . . . . . . . . . . . . . . . 21
2.3.2 Evolution stratigraphique des bassins oshore angolais. De Browneld and
Charpentier 2006b. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.3.3 Interprétation géologique régionale du champ Alpha. . . . . . . . . . . . . 23
2.3.4 Vue 3D de la grille réservoir R-Alpha, avec une section sismique et l'em-
placement du puits. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.3.5 Deux couches de la grille réservoir R-Alpha. La propriété illustrée repré-
sente les faciès lithologiques : argiles en rouge, sables grossiers (en vert) à
ns (en bleu). La limite des données sismiques étudiées a été ajoutée sur
l'image. Le point blanc correspond à la position du puits W-Alpha. . . . . 25
2.3.6 Logs pétrophysiques, pétroélastiques et lithologiques de W-Alpha. De gauche
à droite : Gamma-Ray, lithologie, densité, vitesses P et S, porosité eective
et saturation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.3.7 Diagramme croisé des propriétés Impédances relatives versus coecient de
Poisson relatif, coloré par les facies au puits. . . . . . . . . . . . . . . . . . 26
2.3.8 Courbes d'enfouissement estimées pour le puits W-Alpha, pour les proprié-
tés densité, vitesses P et S. Les points rouges correspondent aux données
du puits W-Alpha, les points verts aux données du puits voisin. . . . . . . 27
2.3.9 Section sismique Near suivant la direction inline (à gauche) et crossline (à
droite). Le réservoir est individualisé par les horizons bleus et la trajectoire
du puits est représentée en jaune. . . . . . . . . . . . . . . . . . . . . . . . 28
2.3.10 Ondelettes obtenues par calibration au puits sur le stack complet des don-
nées sismiques (à gauche), sur les stacks proches et moyens osets (à droite). 28
2.3.11 Sections en impedances P (à gauche) et S (à droite) pour les inversions
deterministes (en haut) et stochastiques (en bas). L'inversion stochastique
est réalisée uniquement sur le réservoir R-Alpha. . . . . . . . . . . . . . . 29
2.4.1 Localisation géographique du champ Beta. . . . . . . . . . . . . . . . . . . 30
xi
xii Table des gures
3.4.3 Eet du changement d'échelle sur une trace sismique. (a) l'impédance P
de l'inversion stochastique à l'échelle sismique, (b) l'impédance P dans le
voxet, (c) trace sismique réelle, (d) synthétique obtenu à partir de a, (e)
synthétique du voxet obtenu à partir de b et (f) diérence entre les deux
synthétiques. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.4.4 Eet du changement d'échelle sur les sections sismiques. En haut, la sec-
tion issue de la modélisation sur le volume à l'échelle sismique (quelques
cellules vides dans le volume d'impédance sismique apparaissent comme
des points sur l'image sismique) ; en bas, la section issue de la modéli-
sation sur le volume à l'échelle réservoir. Toutes ces modélisations sont
réalisées par convolution 1D sur chacune des traces du volume 3D. Les
traces synthétiques de la gure 3.4.3 ont été replacées sur la gure. . . . . 52
3.5.1 Application d'un ltre moyen sur un point donné et ses p traces de part
et d'autre, en x et y. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.5.2 Application d'un ltre moyen sur la section sismique issue des propriétés
du PEM du champ Alpha. En haut, la section brute ; en bas, la section
ltrée. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.5.3 Application d'un ltre moyen sur la section sismique issue des propriétés
du PEM du champ Bêta. En haut, la section brute ; en bas, la section ltrée. 55
3.6.1 Comparaison des données sismiques synthétiques et réelles du champ Bêta,
basée sur les cartes de Kohonen, avec 20 neurones. A gauche, la classica-
tion des données synthétiques issues des résultats d'inversion portés dans
la grille réservoir ; à droite, la classication des données réelles. . . . . . . 56
3.6.2 Comparaison des données sismiques synthétiques et réelles, basée sur les
cartes de Kohonen, avec 20 neurones. A gauche, la classication des don-
nées synthétiques issues du PEM ; à droite, la classication des données
réelles ; en bas, la concordance des classes : bleu foncé si la classe est la
même sur les deux volumes, bleu clair sinon. . . . . . . . . . . . . . . . . . 58
4.5.1 Slices of the F1 Probability cube along the main channel of R-Alpha and
the F1 most probable value in the initial reservoir grid. . . . . . . . . . . . 82
4.5.2 Facies repartition in the initial reservoir model, on layer 21, with dierent
shale thresholds in percent. . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.5.3 Facies repartition in the initial reservoir model, on layer 21, with the shale
threshold of 87%, determined from the well data calibration. . . . . . . . . 84
4.5.4 Facies repartition in the 21st layer of the structural models. The rst line
shows the geological interpretation, the picture on the left of the second
line shows the facies repartition in the initial reservoir grid. . . . . . . . . 85
4.5.5 Facies repartition in the 105th layer of the structural models. The rst line
shows the geological interpretation, the picture on the left of the second
line shows the facies repartition in the initial reservoir grid. . . . . . . . . 86
6.2.1 Workow for the supervised classication with neural networks. . . . . . . 119
Table des gures xvii
6.3.1 Example of learning phase for an Articial Neural Network. An input vec-
tor is given to the network. The winning neuron (in red) is determined,
and then updated so as to be closer to the input data. . . . . . . . . . . . 121
6.3.2 Time varying input patterns such as seismic traces can not be evaluated
by Euclidian distance. Here the repetition of the same trace with time shift
will be identied as dierent signals if the Euclidian distance is used. If we
use the correlation between the trace and the neuron, this series will be
seen as the same sequence. . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
6.4.1 A seismic section for the Beta eld, showing the reservoir of interest deli-
mited by the green horizons. . . . . . . . . . . . . . . . . . . . . . . . . . . 124
6.4.2 Geological correlation for the wells from Beta eld. We can see the dierent
geological bodies and how some markers appear or disappear laterally. . . 125
6.4.3 Correlation sheet for the wells from Beta eld, with only some of the mar-
kers shown, for clarity. The property shown here is Vp . . . . . . . . . . . . 125
6.4.4 Well log blocking results for well W-Beta3. This gure shows the original
logs in black and the blocked logs in green for (a) the density ρ, (b) the P
velocity Vp and (c) the S velocity Vs . The near seismic trace is compared
to (d) the synthetic of the blocked logs, (e) the synthetic of the original
logs. The well-to-seismic tie with the original logs is of poor quality. . . . . 126
6.4.5 Well log blocking optimization procedure for well W-Beta3. This gure
shows the original logs in black, the blocked logs in green and the optimized
logs in red for (a) the density ρ, (b) the P velocity Vp and (c) the S velocity
Vs . The near seismic trace is compared to (d) the synthetic of the blocked
logs, (e) the synthetic of the best logs and (f) the synthetic of the original
logs. The well-to-seismic tie has been greatly improved. . . . . . . . . . . . 127
6.4.6 The normalized tness function of the blocking optimization. . . . . . . . 127
6.4.7 At each location U, for each parameter and each layer, the SGS method
draws values Y from the CCDF provided by the parameter statistics ex-
tracted from the wells. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
6.4.8 The Massive Modeling. The sections are respectively ρ, Vp , Vs and the
synthetic seismograms. The horizontal scale is the trace number. . . . . . 129
6.4.9 Non supervised Kohonen maps realized with an interval thickness of 50ms
and 180ms, with 20 neurons. . . . . . . . . . . . . . . . . . . . . . . . . . . 130
6.4.10 Non supervised Kohonen maps with respectively 5, 20 and 50 neurons.
Next to each map, we show the resulting neurons (for the 50 neurons map,
we show only the second half). Some traces from the 50 neurons map
represent noise. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
6.4.11 Non supervised Kohonen maps realized with neighbourhood radius of 1
and 10 neurons. The rst line shows the neuron repartition, the second
line the neurons. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
6.4.12 Fitness map for the non supervised Kohonen maps realized with neighbou-
rhood radius of 1 and 10 neurons. . . . . . . . . . . . . . . . . . . . . . . . 133
6.5.1 A small part of the Massive Modeling volume, with the neuron repartition
(top line) and the tness (bottom line) for each trace. . . . . . . . . . . . . 134
6.5.2 A small part of the color-coded visualization of the classication. On the
top-row, the traces are shown as they are in the dataset. On the bottom-
row, the traces are ordered by class and tness. Since we show about 300
traces, once ordered we see only a small portion of two classes on the
bottom row. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
xviii Table des gures
6.5.3 A graph showing the percentage of model traces exceeding the tness on
the abscissa. For example, 50% of the model traces have a tness greater
than 90%, which is fully satisfactory. . . . . . . . . . . . . . . . . . . . . . 135
6.5.4 A part of the table for the classication analysis. We obtain the parame-
ter average and standard deviation for each neuron. Here are shown the
thickness, ρ and Vp for the two rst layers and the three rst neurons. . . 135
6.5.5 Likelihood of each model traces with the synthetics of the ve wells for
Beta eld. A likelihood of 1 is the best correlation value between a model
trace and a synthetic trace. . . . . . . . . . . . . . . . . . . . . . . . . . . 136
6.5.6 Predictive modeling : the model traces from the training phase have been
applied to the actual seismic traces. The maps show the classication and
the tness. The white circle on the tness map gives the location of well
W-Beta1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
6.5.7 Likelihood of each model traces with the actual seismic trace of the ve
wells for Beta eld. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
6.5.8 Descrictive modeling : the training phase has been performed on the actual
seismic traces. The maps show the classication and the tness. . . . . . . 139
6.5.9 Likelihood of each model traces with the actual seismic trace at the ve
wells for Beta eld. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
6.5.10 Graph showing the percentage of model traces exceeding the tness on the
abscissa. For example, 60% of the model traces have a tness greater than
50%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
6.5.11 A small part of the color-coded visualization of the classication. The colors
displayed on the top row are quite dull, this means that the tness for the
traces presented here is low. See gure 6.5.2 . . . . . . . . . . . . . . . . . 141
6.5.12 Likelihood of each model traces with the synthetic traces of the ve wells
for Beta eld. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
6.5.13 The classication map from the petrophysical training with geological in-
terpretation on top, with structural interpretation on the bottom picture. 142
6.5.14 The classication map from the seismic training with geological interpre-
tation on the top, with structural interpretation on the bottom. . . . . . . 143
6.6.1 Well log blocking results. This gure shows the original logs in black and
the blocked logs in green for (a) the density ρ, (b) the P velocity Vp and (c)
the S velocity Vs . The near seismic trace is compared to (d) the synthetic of
the blocked logs, (e) the synthetic of the original logs. The well-to-seismic
tie with the original logs is of poor quality. . . . . . . . . . . . . . . . . . . 146
6.6.2 Well log blocking optimization procedure. This gure shows the original
logs in black, the blocked logs in green and the optimized logs in red for
(a) the density ρ, (b) the P velocity Vp and (c) the S velocity Vs . The near
seismic trace is compared to (d) the synthetic of the blocked logs, (e) the
synthetic of the best logs and (f) the synthetic of the original logs. The
well-to-seismic tie has been greatly improved. . . . . . . . . . . . . . . . . 146
6.6.3 The Massive Modeling. The sections are respectively ρ, Vp , Vs and the
synthetic seismograms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
6.6.4 Non supervised Kohonen maps realized with an interval thickness of 50ms
and 100ms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
Table des gures xix
Introduction générale
1.1 Problématique
Les modèles de réservoirs pétroliers jouent un rôle de plus en plus prépondérant dans
l'industrie pétrolière. Ils sont utilisés tout le long de la vie d'un gisement pour planier les
études complémentaires à eectuer, pour optimiser l'implantation de nouveaux puits, mais
aussi et surtout, estimer les réserves d'hydrocarbures en place et simuler l'exploitation du
prospect réel.
Pour faire tout cela, les spécialistes ont besoin de connaître les propriétés-clés qui per-
mettront d'expliquer (et donc de reproduire) les phénomènes physiques qui gouvernent
le réservoir. Ils ont à leur disposition un certain nombre de données qui leur permettent
d'obtenir ces propriétés, soit de manière directe, soit de manière indirecte. Les travaux de
caractérisation réservoir s'appuient essentiellement sur deux types de données : les diagra-
phies et les données sismiques.
Les diagraphies : Une diagraphie ("well log") consiste à mesurer, pendant (diagraphies
instantanées) ou après (diagraphies diérées) un forage, les caractéristiques des roches
traversées, à l'aide de diérentes sondes. La diagraphie représente en général tout
enregistrement d'une caractéristique d'une formation géologique en fonction de la
profondeur. Les diagraphies instantanées, enregistrées pendant le forage ("LWD :
Logging While Drilling"), permettent de récupérer, lorsqu'elles traversent un milieu :
la teneur en hydrocarbures et/ou eau ;
la porosité ;
la perméabilité.
Toutes ces propriétés sont essentielles pour simuler la vie du réservoir et les spé-
1
2 Chapitre 1. Introduction générale
cialistes cherchent donc à les intégrer dans les modèles. Cependant, les diagraphies
fournissent des informations très locales, car limitées à l'emplacement des puits, et
ces derniers ne sont présents qu'en petit nombre sur toute la surface d'un prospect.
Les modèles de réservoirs calibrés à partir de ces données diagraphiques sont donc
mal conditionnés dès que l'on s'éloigne des puits.
Il est donc évident que les diagraphies ne susent pas pour contraindre les modèles
de réservoir.
Les données sismiques : La sismique est une technique de mesure indirecte qui consiste
à enregistrer en surface des échos issus de la propagation dans le sous-sol d'une onde
sismique provoquée. Ces échos sont générés par les hétérogénéités du sous-sol et se
manifesteront par la présence de réecteurs sur les enregistrements. Le temps d'arri-
vée de l'écho permet de situer dans l'espace la position d'un réecteur et l'amplitude
apporte des informations indirectes sur certains paramètres physiques.
Après un traitement adapté, les données sismiques nous donnent une image de la
structure du sous-sol, ainsi que des informations sur sa nature.
Les données sismiques représentent des mesures avec une bonne couverture spatiale
sur l'ensemble du prospect. Leur utilisation est couramment mise en avant lors de
la construction de l'enveloppe géométrique d'un modèle réservoir. En revanche, les
informations qu'elles contiennent en terme de paramètres physiques importants pour
la description du réservoir sont largement sous-utilisées.
L'intégration quantitative des données sismiques dans les modèles de réservoir, comme nou-
velle contrainte des propriétés loin des puits, est une approche de plus en plus étudiée de
nos jours. Le travail de cette thèse s'inscrit dans cette optique, et a pour but d'améliorer
la cohérence entre les modèles de réservoirs et les données sismiques.
La validation du réservoir passe par une phase de simulation de l'écoulement des uides
pour créer une courbe de production qui doit être comparable à la courbe réelle si le réser-
voir est correctement déterminé. On eectue pour cela un changement d'échelle du modèle
géologique au modèle réservoir de maillage plus grossier. Si l'on parvient à reconstituer sa
courbe de production, cela signie que le modèle est l'une des solutions envisageables pour
le réservoir réel. En revanche, si les courbes de production réelles et simulées dièrent, il
faut revenir au modèle et le modier. Cette validation dynamique est dicile à mettre en
oeuvre car la courbe de production n'indique pas où les erreurs se trouvent dans le modèle
réservoir. De plus, elle est coûteuse puisqu'il s'agit d'entrer le modèle de réservoir, très vo-
lumineux et très complexe, dans un simulateur d'écoulement qui reproduit le mouvement
des uides (comportement dynamique) à travers le réservoir pendant son exploitation, ainsi
que la production du champ étudié.
Dans le cadre de cette thèse, nous proposons une approche complémentaire qui consiste
à travailler directement sur le modèle réservoir dynamique. Nous présentons donc dans
ce manuscrit des méthodes de caractérisation réservoir qui visent à intégrer l'information
dans le modèle réservoir, sans passer par le modèle géologique, ce qui évite le changement
d'échelle entre ces deux supports.
Pour valider le réservoir, nous choisissons d'évaluer la qualité du modèle réservoir par
la reconstitution des données sismiques. Il s'agit donc d'une validation statique, car cette
validation est faite avec un volume unique de données sismiques. Cette approche a l'avan-
tage de rendre le modèle réservoir cohérent avec la sismique.
Cela signie qu'à partir de la géométrie et des propriétés du réservoir, un volume sismique
synthétique est créé et comparé aux données sismiques réelles. Les propriétés statiques du
modèle réservoir peuvent donc être validées de cette manière, ce qui réduit leur niveau
d'incertitudes pour la phase de simulation dynamique. En revanche, la comparaison entre
sismique réelle et synthétique n'est pas triviale.
En eet, la création d'un volume sismique synthétique suppose de connaître la fonction
de transfert entre les propriétés du réservoir et la sismique (de la physique des roches aux
données sismiques) et cela est loin d'être le cas. En revanche, nous pouvons faire appel à
des méthodes existantes pour transformer les propriétés du réservoir en propriétés signi-
catives pour la sismique (propriétés élastiques) et ensuite créer une sismique synthétique
idéalisée. Cela implique donc que les données sismiques synthétiques et réelles ne sont pas
directement comparables.
Un autre problème concerne la diérence des échelles d'observation pour les diérents
supports d'information. En eet, entre les données de puits, les données sismiques et les
modèles de réservoir, les informations ne sont pas stockées au même niveau de détail. Il est
donc important de se pencher sur ce problème pour assurer la compatibilité des diérentes
données entre elles, sachant que le support d'information nal est le modèle réservoir. Les
méthodes que nous proposons abordent d'une manière ou d'une autre cette thématique de
changement d'échelle.
qui nous permettent de simuler les propriétés réservoirs sont brièvement décrits.
Le modèle de réservoir représente donc la somme des connaissances disponibles, à savoir les
données sismiques, de puits, leurs interprétations, etc . . . Plus cette connaissance s'enrichit
et plus les incertitudes liées au réservoir diminuent. Nous résumons ci-après l'état de l'art
sur le réservoir pétrolier.
La création du premier modèle réservoir permet de faire les premières études de sensi-
bilité de la sismique aux paramètres pétrophysiques et donne une estimation des réserves
d'hydrocarbures. Si le gisement est économique, il sera mis en production. Cela consiste
tout d'abord à installer un certain nombre de puits (dits puits de production) à des endroits
stratégiques (dôme du réservoir par exemple). La pression interne du réservoir engendre
une déplétion naturelle du réservoir, de courte durée. D'autres puits (puits d'injection)
peuvent ensuite être installés pour assister la récupération des hydrocarbures.
Pour le reste de la thèse, nous nous focaliserons sur le modèle statique, et toutes références
ultérieures aux modèles de réservoirs concerneront la partie statique.
La construction du modèle d'un réservoir se fait généralement en trois étapes (voir gure
1.2.2) :
La construction du modèle structural :
La trame géométrique du réservoir est obtenue par l'interprétation des données sis-
miques, en terme d'histoire dépositionnelle et tectonique, croisée avec l'interprétation
des puits, sur le positionnement des unités stratigraphiques. Le modèle structural
nous donne l'enveloppe, l'extension du réservoir, les objets géologiques présents, ainsi
que les extensions des failles et les contacts entre les uides (eau, huile, gaz).
Pour l'obtenir, on procède au pointé des horizons et des failles sur un volume sis-
mique. On obtient alors des nuages de points dénissant une surface tridimension-
nelle dans l'espace. Ces données surfaciques nous donne une première visualisation
du réservoir (image du haut sur la gure 1.2.2).
Quelle que soit l'approche utilisée, le modèle réservoir doit être continuellement remis en
cause, pour chaque nouvelle donnée acquise à intégrer au modèle. Cette mise à jour est
dicile à mettre en place, conséquence de la complexité des phénomènes physiques qui
s'appliquent aux réservoirs.
Le rééchantillonnage des données de puits Les puits nous donnent accès à l'infor-
mation directe mais ponctuelle des propriétés réservoir à un endroit donné. Comme nous
1.2. Contexte scientique 7
Fig. 1.2.2 Construction du modèle réservoir. (a) Modèle structural ; (b) grille réservoir ;
(c) remplissage du réservoir avec des propriétés (ici facies : en rouge, les argiles ; autres
couleurs, les sables avec diérents niveaux de cimentations).
8 Chapitre 1. Introduction générale
cherchons à récupérer ces propriétés sur l'intégralité du modèle de réservoir, les données de
puits doivent être intégrées dans la grille. Cette intégration n'est possible que si l'informa-
tion puits est correctement positionnée dans l'espace. Les puits doivent tout d'abord être
recalés sur la sismique, ce qui est une thématique activement explorée car dicile (White
and Simm 2003 et Walden and White 1984).
Quand la calibration est correcte, cette intégration est réalisée en identiant les cellules qui
intersectent la trajectoire du puits et en faisant une moyenne des échantillons se retrouvant
au sein d'une même cellule. La méthode utilisée pour moyenner les données dans une cellule
dépendra de la propriété considérée. Une simple moyenne arithmétique pourra être utilisée
pour le volume d'argile par exemple.
Le processus de rééchantillonnage ne doit pas déformer la morphologie géologique ou la
géométrie des vitesses de soniques, car elles sont utilisées pour récupérer les vitesses des
ondes dans le réservoir et pour caler les données de puits aux données sismiques. D'autres
méthodes devront être considérées pour ces propriétés dites "non additives" (Narasimhan
1983).
En ce qui concerne les variables discrètes (comme les faciès), elles sont "moyennées" en
faisant le compte d'échantillon pour chaque faciès et en aectant la valeur du faciès "do-
minant" dans chaque cellule. La validation du rééchantillonnage est généralement réalisée
en vériant que le log original et le log rééchantillonné sont cohérents et que le nouveau
log préserve susamment les détails géologiques et les hétérogénéités verticales.
Un autre contrôle est de vérier les statistiques des logs originaux et rééchantillonnés
(moyenne ; écart-type, variance, corrélation entre plusieurs propriétés, ...). L'eet de sup-
port (log versus grille réservoir) est ainsi facilement visible : l'écart-type et la variance
diminuent avec l'élargissement du support de l'information.
On cherche les paramètres d'un modèle mathématique qui nous donnerait une ré-
ponse aussi proche que possible des données observées.
Il existe de nombreuses relations quantitatives dans la littérature, reliant les propriétés
élastiques des roches aux propriétés pétrophysiques. Nous présentons ci-après les relations
les plus communes.
1
κ= (Vcl +φ) (1−Vcl −φ)
(1.2)
κshale + κquartz
1
µ= (Vcl +φ) (1−Vcl −φ)
(1.3)
µshale + µquartz
Avec ρ la masse volumique (ρshale pour l'argile, ρquartz pour le quartz), Vcl le volume
d'argile, φ la porosité, κ le module d'incompressibilité et µ le module de cisaillement.
Pour une roche clastique sèche, la densité, le module d'incompressibilité et le module de
cisaillement sont estimés à partir des mesures eectuées sur carottes et aux puits. Ainsi,
ρshale et ρquartz sont des mesures faites sur carottes en laboratoire et Vcl et φ sont des
mesures de puits. Ces trois équations nous permettent de recréer les vitesses des ondes
compressives Vp et des ondes de cisaillement Vs par les équations de l'élastodynamisme et
d'étudier la réponse sismique de ces paramètres :
10 Chapitre 1. Introduction générale
r
V = κ+ 43 µ
p
q
ρ
(1.4)
Vs =
µ
ρ
Il existe dans la littérature une multitude de modèles (Avseth et al. 2005), mais aucun
ne donne une vision absolue des relations entre propriétés des roches et des uides et les
paramètres sismiques mesurables. Wang 2001 propose une excellente introduction à la phy-
sique des roches avec application aux données sismiques. Castagna et al. 1993 ore une vue
d'ensemble très complète des diérents types de modèles, leurs applications et implications.
de Gassmann sont ainsi fréquemment utilisées pour estimer les modules des roches.
La conversion opérée par le PEM lie le domaine réservoir au modèle élastique dans lequel
on peut décrire la propagation des ondes. En général, un PEM est spécique au champ sur
lequel il a été calibré.
Il existe également une notion de PEM inverse, à partir de laquelle nous cherchons à
résoudre un problème inverse pour obtenir les distributions des propriétés statiques du ré-
servoir (porosité, teneur en argile et la saturation) à partir des attributs sismiques observés.
Coleou et al. 2006 et Bornard et al. 2005 montre un schéma d'inversion particulier, basé
sur l'utilisation d'un PEM pour récupérer, directement pendant l'inversion, les propriétés
pétrophysiques. En pratique, la conversion opérée par le PEM est non unique et incertaine.
Le PEM est donc combiné aux géostatistiques pour évaluer les incertitudes associées aux
propriétés dérivées des données sismiques.
Nous terminons enn par le septième chapitre, où nous concluons sur l'ensemble des
travaux et nous proposons les perspectives envisagées pour chaque méthode.
2.1 Introduction
Les méthodologies développées dans la suite de la thèse ont été testées et développées
sur plusieurs réservoirs, choisis pour leur contexte géologique et les données disponibles. Ces
données étant condentielles, les trois cas étudiés sont dénommés Alpha, Bêta et Gamma.
Les diérences entre ces modèles sont de plusieurs ordres :
La géologie : Les champs Alpha et Bêta sont des réservoirs turbiditiques situés en eaux
profondes ; le champ Gamma est un réservoir carbonaté avec une faible épaisseur d'eau.
Les données sismiques : Les champs Alpha et Bêta ont une résolution sismique assez
ne qui permet de délimiter la stratigraphie et les corps géologiques, ainsi que leur contenu.
Pour le champ Gamma, les réecteurs sismiques sont basses fréquences et peuvent masquer
l'information sous-jacente.
Les contextes géologiques des champs présentés ici sont inspirés de publications de l'United
States Geological Survey. On se reportera à Browneld and Charpentier 2006b pour le
champ Alpha, à Browneld and Charpentier 2006a pour le champ Bêta et à Pollastro 2003
pour la champ Gamma.
Dans le domaine pétrolier, les méthodes d'imagerie sismique consistent à générer des ébran-
lements dans le sous-sol. Ces ébranlements se propagent sous la forme de vibrations qui
réagissent lorsqu'elles traversent les milieux géologiques. Cela engendre des phénomènes de
réexion et de transmission sur les limites de couches géologiques. L'étude de ces vibrations
mécaniques (les ondes sismiques) est basée sur la théorie de la propagation des ondes.
Les ondes rééchies remontent en surface et sont enregistrées sous la forme d'une énergie
en fonction du temps. Elles sont traitées puis interprétées en terme de caractéristiques de
couches du sous-sol.
15
16 Chapitre 2. Présentation des données utilisées
à l'interface en quatre types d'ondes : une onde P rééchie Pr , une onde P transmise Pt , une
onde S rééchie Sr et une onde S transmise St (voir gure 2.2.1). La loi de Snell-Descartes
nous permet de relier chacune de ces ondes par :
sinθPr sinθPt sinθSr sinθSt
= = = (2.1)
Vp1 Vp2 Vs1 Vs2
Fig. 2.2.1 La propagation d'une onde P sur une interface entre deux milieux induit la
génération d'une onde P rééchie Pr , d'une onde P transmise Pt , d'une onde S rééchie Sr
et d'une onde S transmise St .
Pendant le traitement sismique, les échos primaires (ou réexions) sont les signaux les
plus utilisés pour déduire les caractéristiques et les propriétés de la subsurface. Le rapport
de l'amplitude de l'onde rééchie sur l'amplitude de l'onde incidente est nommé coecient
de réexion. A l'incidence normale (θPr = 0), les ondes converties n'interviennent pas et le
coecient de réexion Rpp est déni par :
Vp2 ρ2 − Vp1 ρ1 Ip2 − Ip1
Rpp = = (2.2)
Vp2 ρ2 + Vp1 ρ1 Ip2 + Ip1
d'amplitude avec l'angle dit AVA (ou l'oset correspondant, AVO). L'équation 2.2 n'est
alors plus utilisable et il faut nous tourner vers les équations dites de " Knott-Zoeppritz ".
Les équations de Zoeppritz supposent deux milieux plans, semi-innis à extension in-
nie.
Aki and Richards 1980 propose une formulation simple de l'équation complète de Zoep-
pritz, dans une forme matricielle normalisée par l'amplitude de l'onde incidente A0 = 1 :
RP P cos θ1
TP P − sin θ1
(2.4)
P ∗ =
RP S − cos 2φ1
TP S − sin 2θ1
avec
Cette équation, qui reste très complexe, n'est pas utilisable de manière directe pour
analyser la variation des coecients de réexion en fonction de ces paramètres. Il est plus
simple de considérer les approximations des équations de Zoeppritz, qui sont nombreuses
dans la littérature. Notamment, l'approximation d'Aki et Richards dénie par :
2
β2
1 ∆α β ∆β 1 ∆
R(θ) = 2
(1 + tan θ) 2
− 4 2 sin θ 2
+ (1 − 4 2 sin θ (2.6)
2 α α β 2 α
Nous justions la sélection de cette approximation par le fait que dans le domaine pétrolier,
l'information pertinente pour le réservoir est enregistrée pour des angles d'incidence allant
jusqu'à 30 voire 40◦ . La gure 2.2.2 montre le domaine de validité angulaire des diérentes
approximations de l'équation de Zoeppritz, dont l'approximation d'Aki et Richards, face à
l'équation originale.
18 Chapitre 2. Présentation des données utilisées
Le traitement
Le traitement sismique est conventionnellement ciblé sur l'interprétation structurale
d'un réservoir. Si l'on souhaite extraire des informations stratigraphiques et lithologiques
des données sismiques pour caractériser le réservoir, il est nécessaire d'appliquer un traite-
ment qui conservera l'énergie des ondes à haute fréquence et les amplitudes relatives du si-
gnal (ou amplitudes préservées). Il faut éviter de détruire les informations utiles contenues
dans les variations relatives de l'amplitude sismique. Nous présentons ci-après certaines
étapes du traitement qui peuvent aecter l'étude des amplitudes.
pas complètement satisfaisante. Pour corriger ces deux eets, on utilisera plutôt une cor-
rection des amplitudes sismiques en fonction du temps et du déport. Pour cela, l'une des
méthodes fréquemment utilisées fait appel aux prols sismiques verticaux, qui fournissent
la trace sismique idéale à l'endroit du puits pour corréler avec les données sismiques (Camp-
bell et al. 2005). Les méthodes utilisées peuvent avoir des eets importants sur la réponse
AVO, comme par exemple des gradients d'amplitude majorés ou des diminutions de l'in-
tercept comme le montre Cambois 2000. Pour éviter cela, la réponse AVO avant égalisation
doit être modélisée pour servir de référence à la procédure.
Déconvolution La déconvolution est un ltre inverse utilisé pour récupérer les hautes
fréquences, égaliser les amplitudes, récupérer une ondelette à phase nulle directement liée à
la réectivité du sous-sol, ou encore pour tout autre processus aectant la forme d'onde. La
déconvolution peut améliorer des données sismiques aectées par le ltrage qui se produit
naturellement lors de la propagation dans les couches géologiques.
La principale contribution de l'ondelette qui est intégrée dans une trace sismique brute est
généralement la source de la signature, mais il y a aussi des contributions de nombreuses
autres parties du champ d'onde de la source au récepteur, telles que les structures super-
cielles autour de la source et du récepteur, les fantômes de la source et du récepteur, la
réponse du récepteur, la réponse du système d'enregistrement et les eets de l'atténuation
le long de la trajectoire des rais. Idéalement, la déconvolution supprime l'eet de toutes
ces contributions de la trace sismique avant interprétation.
L'atténuation des multiples L'atténuation des multiples joue un rôle majeur sur le
bon fonctionnement de l'analyse AVA. Ce problème est particulièrement important pour
les données sismiques marines. La superposition, dans une même fenêtre de temps, des
signaux primaires et de multiples, va modier de manière générale la réponse sismique en
terme d'amplitude, de phase, et de contenu fréquentiel.
Il existe de nombreuses méthodes pour atténuer les multiples (Yilmaz 2001), notamment
l'utilisation d'un ltre de Radon parabolique sur les collections CMP puisque dans ce
domaine les primaires s'aplatissent, mais les multiples restent hyperboliques.
avec Θ l'angle d'incidence à l'oset x ; t0 le temps d'arrivée à oset nul sur la surface
considérée ; VRM S la vitesse RMS jusqu'à la surface considérée ; VIN T la vitesse d'intervalle
20 Chapitre 2. Présentation des données utilisées
Le champ Alpha se situe au large des côtes angolaises dans le bassin du Congo inférieur,
sur la marge continentale ouest africaine. On se trouve dans un contexte en eaux profondes
(tranche d'eau de plus de 1500 mètres).
L'évolution des bassins côtiers de l'Angola provient des événements de la tectonique des
plaques qui a causé le rifting du supercontinent Gondwana et la séparation de l'Afrique
de l'Amérique du sud. Le rifting a commencé au Crétacé Inférieur et un élargissement
progressif de l'Océan Atlantique Sud s'est amorcé. Après le début de cet élargissement, les
bassins de l'Angola, dont le bassin du Congo, ont continué à subsider, permettant le dépôt
d'épais sédiments post-rift. L'évolution tectonique et sédimentaire de la marge angolaise
peut être divisée en huit phases séquentielles. La gure 2.3.2 résume ces huit phases.
Le champ Alpha s'est développé pendant les phases VI, VII et VIII, du Tertiaire Inférieur
au Miocène Moyen. Ces phases sont constituées d'une grande proportion d'argiles et de
quelques carbonates, reétant l'ouverture continue de l'Atlantique. L'Oligocène a été une
période de larges dépôts de turbidites. Les dépôts oligocènes sont essentiellement turbidi-
tiques, avec des chenaux sableux intercalés. La formation Malembo, d'âge Oligocène, est
22 Chapitre 2. Présentation des données utilisées
une unité lithostratigraphique transgressive qui contient des strates formées proche littoral
et en contexte de talus. Cette formation inclut certains des plus importants réservoirs dans
les eaux profondes de l'Angola.
La structure du champ Alpha, présenté en contexte régionale sur la gure 2.3.3, est un
anticlinal pincé en contexte compressif, initié à la fois par les mouvements salifères et par
une déformation de glissements gravitaires. Cette structure, inclue dans la formation Ma-
lembo, se forme pendant le déplacement gravitaire oligocène et se stabilise vers la n du
Miocène. L'unique puits dévié (W-Alpha) de la structure a été foré sur le anc sud-ouest
de l'anticlinal et traverse le at spot des réservoirs oligocènes. Cette structure est aectée
par des failles d'extrados qui ne brisent pas la couverture.
2.3. Le champ Alpha - modèle turbiditique 23
On retrouve dans ces réservoirs des chenaux de types érosifs (dans l'unité E) de dimension
kilométrique et assez épais, à des chenaux de types dépositionnels (dans l'unité I), très ns.
Le chenal de l'intervalle E est particulièrement visible, de plusieurs centaines de mètres de
large, tandis que les chenaux des autres intervalles sont très ns, et à peine visible sur les
données sismiques disponibles. Les illustrations des travaux sur le champ Alpha montreront
donc essentiellement le chenal de l'intervalle E.
La grille réservoir du champ Alpha (vue 3D sur la gure 2.3.4), que nous appellerons
R-Alpha, s'étend sur une surface de 7km x 7km et englobe les cinq unités réservoirs E, F,
G, H et I. Ce modèle, déni en profondeur, a été construit à partir d'une interprétation
structurale régionale (horizons majeurs et failles de l'anticlinal) et des données de l'unique
puits disponible. La grille est dénie par 114 ∗ 159 ∗ 113 cellules, soient 2.048.238 cellules
au total, de dimension moyenne de 60m x 70m x 10m (ou 5ms d'épaisseur). Comme le
montre la gure 2.3.5, la grille réservoir est plus large que les données sismiques qui nous
intéressent.
Fig.2.3.4 Vue 3D de la grille réservoir R-Alpha, avec une section sismique et l'emplace-
ment du puits.
2.3. Le champ Alpha - modèle turbiditique 25
Fig. 2.3.5 Deux couches de la grille réservoir R-Alpha. La propriété illustrée représente
les faciès lithologiques : argiles en rouge, sables grossiers (en vert) à ns (en bleu). La limite
des données sismiques étudiées a été ajoutée sur l'image. Le point blanc correspond à la
position du puits W-Alpha.
Un seul puits a été foré sur la structure anticlinale du champ. Il s'agit d'un puits
fortement dévié, avec des déviations allant jusqu'à 45◦ au fond du puits. Parmi les logs
disponibles, nous avons les logs pétrophysiques du Gamma-Ray, perméabilité, porosité ef-
fective et totale, saturation, volume d'argile, et les logs pétroélastiques, densité, vitesses
P et S, coecient de Poisson. Nous disposons également d'une interprétation en faciès,
montrée sur la gure 2.3.6 avec les autres logs d'intérêt. Sept faciès ont été interprétés à
partir de ces logs, et sont légendés ainsi sur la gure 2.3.7 :
1. Les argiles
2. Les sables laminés à huile
3. Les sables moyens à huile
4. Les sables grossiers à huile
5. Les sables cimentés à huile
6. Les débris ows
7. Les sables à eau (toutes granulométries confondues)
La discrimination des faciès en fonction des propriétés est plutôt de bonne qualité, comme
illustré sur la gure 2.3.7. Le chevauchement du nuage de points ne dépasse pas les 33%
(ce pourcentage est donné par la proportion de faciès extérieurs à un faciès dominant dans
une cellule élémentaire du diagramme criosé).
On remarquera que les faciès argiles, débris ows et sables à eau sont confondus sur le
diagramme et n'ont pas de discrimination pétrophysique. Ces trois faciès peuvent donc
éventuellement être regroupés ensemble dans le code faciès 1.
La bonne discrimination entre le faciès 1 et les sables à huile est une observation promet-
teuse pour les analyses AVO. Cela garantit que plus les écarts seront importants entre
argiles, sables à eau et sables à huile, plus les anomalies AVO seront fortes.
Fig. 2.3.7 Diagramme croisé des propriétés Impédances relatives versus coecient de
Poisson relatif, coloré par les facies au puits.
2.3. Le champ Alpha - modèle turbiditique 27
Fig. 2.3.8 Courbes d'enfouissement estimées pour le puits W-Alpha, pour les propriétés
densité, vitesses P et S. Les points rouges correspondent aux données du puits W-Alpha,
les points verts aux données du puits voisin.
Il s'agit d'un bloc 3D de données sismiques marines acquises avec une tranche d'eau
d'environ 1500m. Les caractéristiques principales de la zone étudiée sont les suivantes :
Echantillonnage spatial : 12.50m x 12.50m
Dimension : 151 Inlines x 661 crosslines x 1668 échantillons (tmin = 1500ms)
Pas d'échantillonnage en temps : 3ms
Fréquence dominante : 30 Hz (résolution de 21m)
Bande passante : 5 - 100 Hz
Trois substacks (Near stack 5◦ -18◦ , mid stack 18◦ -32◦ et far stack 32◦ -47◦ ) ont été créés à
partir de cette acquisition. La gure 2.3.9 montre une section de la sismique Near suivant
les directions inline et crossline. Le réservoir est individualisé par les horizons bleus.
28 Chapitre 2. Présentation des données utilisées
Fig. 2.3.9 Section sismique Near suivant la direction inline (à gauche) et crossline (à
droite). Le réservoir est individualisé par les horizons bleus et la trajectoire du puits est
représentée en jaune.
Deux types d'ondelette (gure 2.3.10) ont été obtenus par calage de la sismique et
du puits : une ondelette sur le stack complet, une ondelette par stack angulaire. De par
la complexité de la géologie à l'endroit du puits (fort pendage, failles), ce calage est de
qualité moyenne (corrélation autour de 60%). Ces ondelettes représentent tout de même
de manière correcte le signal et son contenu fréquentiel.
Fig. 2.3.10 Ondelettes obtenues par calibration au puits sur le stack complet des données
sismiques (à gauche), sur les stacks proches et moyens osets (à droite).
Le champ Bêta est situé dans la zone oshore au large du Nigeria, dans le Golf de
Guinée. On se trouve dans un contexte en eaux profondes (tranche d'eau de 1400 mètres).
Le Golfe de Guinée est bordé par des zones costales et oshores de la Côte d'Ivoire,
du Ghana, du Togo, du Bénin et de la partie ouest du Nigeria. Ces bassins partagent des
caractéristiques stratigraphiques et structurales communes, car ils contiennent des roches
de l'Ordovicien à l'Holocène. La limite Est correspond à la province du Delta du Niger et
la limite Ouest à la côte ouest africaine. L'histoire géologique du Golfe de Guinée se di-
vise en trois grandes périodes : la période intracratonique du Protérozoïque au Jurassique,
la période de rifting du Jurassique au Crétacé, la période de dérive du Crétacé à l'Holocène.
Ces trois périodes représentent le développement des bassins du Golfe de Guinée. Les
roches de la période intracratonique dans le bassin du Bénin, au large du Delta du Niger,
sont représentés par des roches de type grés, argiles et conglomérats, déposés dans des en-
vironnements deltaïques et uviaux. Pour la période rifting, on retrouve le même type de
roches dans des environnements deltaïques lacustres et uviaux. Les dépôts se retrouvent
dans des séries de grabens et semi-grabens. La n de la période rifting est marquée par
une discordance majeure qui sépare ces sédiments des roches marines de la période de
dérive. Les sédiments de la période de dérive dans le Golfe de Guinée sont de manière
prédominante des grés, des argiles et quelques carbonates, déposés par des successions de
régressions et de transgressions. Cette géologie régionale est représentée sur la gure 2.4.2
(Corredor et al. 2005).
2.4. Le champ Bêta - modèle turbiditique 31
Fig. 2.4.2 Interprétation géologique régionale du champ Bêta, d'après Corredor et al.
2005.
Le bassin sous-jacent a été rempli par des séquences de sédiments provenant du Delta
du Niger, qui se sont progressivement accumulés à la sortie du système de rivières Niger-
Bénoué. C'est l'un des plus grands deltas du monde, étendu sur plus de 300 kilomètres, et
l'épaisseur des sédiments clastiques qui le constituent est estimée à 12 kilomètres.
Fig. 2.4.3 Interprétation géologique du réservoir R-Bêta avec une section temps des
données sismiques (à gauche).
Le remplissage initial a été réalisé à partir des propriétés types de uides et faciès,
interprétés dans la grille réservoir par le scénario géologique le plus probable. Les propriétés
pétrophysiques ont ensuite été analysées aux puits, leurs intervalles de valeur récupérés par
2.4. Le champ Bêta - modèle turbiditique 33
faciès et types de uides. Les propriétés sont ensuite générées dans la grille par simulation
géostatistique.
Fig. 2.4.5 Deux couches de la grille réservoir R-Bêta. La propriété illustrée représente
le Net-To-Gross. La limite des données sismiques étudiées a été ajoutée sur l'image. Les
points blancs correspondent aux positions des puits W-Bêta1 à W-Bêta5.
Fig. 2.4.6 Implantation des puits du champ Bêta dans les intervalles chenalisés du réser-
voir.
34 Chapitre 2. Présentation des données utilisées
Parmi les logs disponibles, nous avons les logs pétrophysiques du Gamma-Ray, perméa-
bilité, porosité eective et totale, saturation, volume d'argile, et les logs pétroélastiques,
densité, vitesses P et S, coecient de Poisson. Nous disposons également d'une interpréta-
tion en faciès, montrée sur la gure 2.4.7 avec les autres logs d'intérêt. Les faciès interprétés
sont légendés ainsi :
1. Argiles silteuses
2. Argiles turbiditiques avec laminations sableuses très nes
3. Silts et sables ns en bancs alternés centimétriques
4. Silts et sables ns en bancs alternés décimétriques
5. Sables moyens
6. Sables grossiers à graveleux
7. Matrice argileuse à quartz
8. Matrice sablo-argileuse à quartz
9. Sables ns en bancs massifs décimétriques
10. Sables moyens en bancs massifs plurimétriques
11. Niveau conglomératique de base de chenal ou Lag
12. Sills ou dykes sableux métriques
La discrimination des faciès en fonction des propriétés est assez dicile, avec un chevau-
chement entre 41% et 52% (en fonction du puits), comme illustré sur la gure 2.4.8. Ceci
est dû au fait que la réponse d'un faciès varie légèrement de puits à puits et d'intervalle en
intervalle.
2.4. Le champ Bêta - modèle turbiditique 35
Un modèle pétroélastique a été établi à partir des cinq puits W-Beta, pour reconstituer
les propriétés pétroélastiques ρ, Vp et Vs à partir des propriétés pétrophysiques. Ce PEM
consiste en une série de fonctions déterministes, transformant l'espace volume d'argile -
porosité eective en propriétés pétroélastiques ρ, Vp , Vs et Ip, Is, P R. La comparaison
entre logs réels et logs estimés par le PEM est très satisfaisante à l'échelle du champ (voir
gure 2.4.9).
Fig. 2.4.9 Comparaison des logs réels (en noir) et simulés par le PEM (en rouge) pour
Vp , pour tous les puits du champ Bêta.
36 Chapitre 2. Présentation des données utilisées
Nous disposons de cinq substacks pour ces données (near stack 2◦ -12◦ , mid stack1 12◦ -22◦ ,
mid stack2 22◦ -32◦ , far stack 32◦ -42◦ and ultra-far stack 42◦ -50◦ ). La gure 2.4.10 montre
une section de la sismique réexion Near suivant la direction inline et crossline. Le réservoir
est individualisé par les horizons rouge et orange.
Fig. 2.4.10 Section sismique Near suivant la direction inline (à gauche) et crossline (à
droite). Le réservoir est individualisé par les horizons rouge et orange et la trajectoire des
puits est représentée en rouge.
Nous disposons d'une ondelette pour chacun des substacks (gure 2.4.11). Elles ont
été calées sur les sismiques à partir des logs de puits. Le calage est de bonne qualité pour
les puits 1 à 3 (corrélation de plus de 85%), et de qualité moyenne pour les puits 4 et 5
(corrélation de 70% à 65%).
2.4. Le champ Bêta - modèle turbiditique 37
Fig. 2.4.11 Ondelettes obtenues par calibration aux puits sur les diérents substacks. En
mauve : ondelette near ; en rouge : ondelette mid1 ; en vert : ondelette mid2 ; en jaune :
ondelette far.
Nous disposons de 5 puits sur ce champ, dont les données sont très bien corrélées aux
données sismiques, et très bien représentées par le PEM. Ce champ est donc bien contraint
au niveau pétrophysique.
2.5. Le champ Gamma - modèle carbonaté 39
Les données étant condentielles, seules des généralités sur le champ sont présentées ici,
ainsi que quelques particularités des données pertinentes en regard des travaux de ce ma-
nuscrit.
Géologiquement, le golfe Persique est une énorme plate-forme, un socle cristallin recou-
vert de sédiments, dont l'épaisseur dépasse 6 km dans le golfe lui-même, et où toutes les
ères géologiques sont représentées.
Cinq phases d'évolution tectonique sont associées à la plaque arabique :
1. Une phase précambrienne de compression et d'accrétion continentale qui a conduit à
l'assemblage de la plaque arabique.
2. Une phase dévonienne d'extension intracratonique liée à l'ouverture de la Paléo-Téthys.
3. Une phase d'âge permien d'alternances de compressions et d'extensions correspondant
à l'orogène hercynien.
4. Une phase mésozoique d'extension et de rifting continental. Les dépôts sont dominés par
une sédimentation carbonatée et évaporitique.
5. La dernière phase, du Mésozoique Supérieur à l'Actuel, durant laquelle la plaque ara-
bique se retrouve dans un environnement de marge active.
Fig. 2.5.1 Principes de la dolomitisation hydrothermale dans les niveaux non perméables
du réservoir Gamma.
Les données sismiques : Il s'agit de données sismiques marines acquises en eaux peu
profondes (70 mètres maximum). On observe une importante diérence de contenu fréquen-
tiel et résolution entre les diérents substacks. La sismique proche oset est très bruitée
et aectée par des multiples internes, notamment sur l'intervalle juste au-dessus du réser-
voir d'intérêt. La sismique lointain oset, très basse fréquence, n'a plus les multiples. La
sismique intermédiaire est un compromis entre altération par les multiples et qualité du
signal.
La qualité moyenne de la calibration puits-sismique (75% au mieux) est un frein pour l'in-
tégration de l'information sismique dans le réservoir. Il faudra donc rester prudent sur les
résultats des travaux qui feront intervenir les ondelettes.
Partie I
41
Chapitre 3
3.1 Introduction
Dans le domaine pétrolier, un soin particulier est apporté à la réalisation du modèle
géologique tout d'abord, puis à la réalisation du modèle réservoir. Ce dernier est le sup-
port principal de l'information puisqu'il est ensuite utilisé comme donnée d'entrée dans
des simulateurs sophistiqués, an de simuler et prédire le comportement du réservoir lors
de son exploitation (Enge et al. 2007). Ces simulateurs faisant appel à des méthodes nu-
mériques très complexes et étant très coûteux tant en mémoire qu'en temps de calcul, une
approche essai-erreur sur les propriétés du réservoir est préjudiciable. En eet, les erreurs
peuvent provenir des propriétés statiques et/ou des propriétés dynamiques. Pour permettre
une meilleure intégration, nous proposons une approche complémentaire : la modélisation
sismique comme méthode de validation du modèle réservoir. Ainsi, nous pouvons valider
les propriétés statiques en regard de la sismique, et donc diminuer leurs incertitudes. Lors
de la phase de simulation dynamique du réservoir, les erreurs seront donc essentiellement
provoquées par les incertitudes sur les propriétés dynamiques.
La modélisation sismique est un processus de calcul par lequel un modèle géologique est
transformé en un enregistrement sismique synthétique. Si l'on considère l'utilité de la modé-
lisation sismique, seule, sans faire référence à la grille réservoir, les enregistrements synthé-
tiques générés par ce processus ont plusieurs vocations et utilités au cours de l'exploration
d'un champ pétrolier et sont souvent générés avant et après l'acquisition des données de
sismique de réexion. Les sismiques synthétiques générés avant acquisition sont typique-
ment utilisées pour ajuster les paramètres de l'acquisition (géométrie, nombre de sources
et de récepteurs, . . . ), an d'assurer que l'objectif géologique sera imagé et interprétable
sur les données après traitement.
Quand ces volumes synthétiques sont générés après l'acquisition, l'interprétateur cherche
généralement à déterminer les caractéristiques des roches et des uides, qui sont à l'origine
des évènements observés sur les données sismiques réelles.
La modélisation sismique devient alors une aide à l'interprétation (Breton et al. 2007).
La diérence entre données réelles et synthétiques permet de mettre en évidence les er-
reurs du modèle réservoir, qu'il faut alors corriger, tout en maintenant la cohérence avec
les autres types de données (puits, carottes, . . . ) qui servent de contraintes au modèle.
Toutefois, ces comparaisons sont réalisées avec plus ou moins de pertinence. Pour que le
volume synthétique soit aussi proche des données réelles que possible, il faut connaître les
structures géologiques, leurs propriétés, mais aussi le signal sismique. Le choix de l'onde-
43
44 Chapitre 3. Modélisation sismique à partir de la grille réservoir
lette, qui représente la forme du signal à modéliser, est essentiel pour reproduire le contenu
fréquentiel et la phase du signal sismique.
Le support du modèle réservoir est une grille irrégulière en espace. On peut simuler di-
rectement un volume sismique régulièrement échantillonné à partir d'une grille réservoir,
mais cela passe par l'utilisation de méthodes complexes et diciles à mettre en oeuvre.
Pour que la modélisation sismique sur grille réservoir soit facile et rapide à mettre en
oeuvre, nous créons un support intermédiaire entre données sismiques et modèle réservoir.
dans une petite fenêtre de temps, où l'on entend par ondelette la source après sa modi-
cation par l'atténuation et les multiples au-dessus de la fenêtre.
Dans la modélisation par tracé de rayons (Virieux and Farra 1991), on fait une approxi-
mation haute fréquence de la solution exacte de l'équation d'onde. Cette approximation
repose sur plusieurs hypothèses :
Le champ de vitesse du réservoir varie à plus basse fréquence que l'onde sismique
(sa longueur d'onde est plus grande que celle de l'onde).
L'onde est supposée plane.
La direction de propagation suit la loi de Snell-Descartes.
Ces simplications font du tracé de rayons un bon candidat pour simuler rapidement un vo-
lume sismique synthétique à partir d'un modèle "basse fréquence". Ces méthodes sont très
utilisées pour l'interprétation structurale mais peu indiquées pour faire de la caractérisation
réservoir.
Les bords du modèle sont traités comme des frontières physiques, ce qui provoque
des artefacts.
Ces inconvénients peuvent être résolus par la mise en place de schémas numériques adaptés.
En revanche, ce type de modélisation génère des données sismiques brutes, auxquelles il
nous faut appliquer toutes les étapes du traitement. Ces méthodes ne sont pas applicables
à notre démarche. Puisqu'il s'agit de faire de la modélisation sismique à partir d'un modèle
réservoir 3D dans le but de valider ce dernier, il nous faut considérer des méthodes rapides
et simples à mettre en oeuvre.
C'est pourquoi nous ne considérerons à partir de maintenant que la méthode de la convo-
lution multi-1D pour créer un volume synthétique tridimensionnel. Cette méthode de mo-
délisation ne prend pas en compte les multiples et les artefacts de traitement. On simule
donc une sismique dans un monde "parfait".
Comme mentionnée en 1.2.1.2, la grille réservoir est basée sur le maillage structuré
irrégulier (voir gure 3.3.1). Dans la classication des maillages donnée par Souche 2005
et Moyen 2005, "irrégulier" signie que les noeuds de la grille ne sont pas positionnés de
manière prévisible dans l'espace, contrairement à un maillage régulier pour lequel il sut
des coordonnées à l'origine et d'une fonction mathématique pour passer d'un noeud à un
autre. Les coordonnées de chaque noeud devront donc être stockées de manière explicite.
Généralement, sur la grille réservoir, cet aspect est particulièrement visible par l'hétérogé-
néité de la taille des cellules (voir gure 3.3.1). "Structuré" signie que les connexions entre
chaque noeud sont prévisibles. De cette manière, un noeud ayant les coordonnées (u, v, w)
aura pour voisins les noeuds :
(u − 1, v, w),
(u + 1, v, w),
(u, v − 1, w),
(3.4)
(u, v + 1, w),
(u, v, w − 1),
(u, v, w + 1)
Ces maillages sont particulièrement adaptés en modélisation géologique. Cependant,
comme nous cherchons à créer un volume sismique synthétique comparable aux données
réelles, il est évident que la grille réservoir n'est pas un support adéquat en regard des
méthodes de modélisation sismique choisies. Il est donc nécessaire de passer par un support
intermédiaire, dont le maillage permettra de simuler des amplitudes sismiques pour chaque
échantillon du cube sismique original.
De tels maillages, dénis "réguliers structurés", sont appelés voxets (contraction de
voxel set ) ou plus simplement grilles cartésiennes. Un voxet est un objet tridimensionnel,
régulièrement échantillonné dans toutes les directions de l'espace, et fréquemment utilisé
en traitement sismique (Cognot et al. 1994). Dans le cas d'un volume sismique, le voxet est
déni par des distances entre chaque noeud correspondant au bin sismique suivant les axes
3.3. Transformation en grille régulière 47
Fig. 3.3.1 Cellules des maillages "CPG" (irrégulier structuré) et cartésiens (régulier
structuré).
Pour éviter d'avoir à gérer cela, nous utiliserons un "voxet déformé". Globalement, ce
voxet déformé peut être perçu comme un ensemble de puits virtuels à chaque point du
quadrillage déni en x, y. Verticalement, l'échantillonnage est basé sur l'intersection de
chaque pseudo-puits avec la grille réservoir, c'est-à-dire que pour chaque couche du mo-
dèle, les valeurs des propriétés et l'épaisseur des couches seront respectées (voir les gures
3.3.2 et 3.3.3). Ainsi, la complexité structurale du réservoir est préservée, tout en évitant le
changement d'échelle suivant l'axe vertical. Le produit de convolution prendra en charge le
changement d'échelle verticale. Dans la suite du manuscrit, nous ferons référence au voxet
déformé sous le nom de voxet.
Fig. 3.3.2 Transformation d'une cellule réservoir en une série de voxels (cellules du voxet).
48 Chapitre 3. Modélisation sismique à partir de la grille réservoir
Fig. 3.3.3 L'intersection de la grille réservoir (exemple du champ Alpha) avec un pseudo
puits (en rouge) permet de récupérer les propriétés du réservoir sans changement d'échelle.
La gure 3.4.2 montre le nombre d'échantillons sismiques moyennés par cellule du ré-
servoir, et l'écart-type de ces échantillons, sous la forme d'un pourcentage de variation par
rapport à la moyenne. Le fait que l'écart-type des échantillons dans le modèle réservoir
reste faible, suggère que le changement d'échelle a été réalisé de manière satisfaisante.
3.4. Evaluation des modélisations sur le champ Alpha 49
Fig.3.4.1 Eet du changement d'échelle entre les diérents supports d'information sur le
champ Alpha. En haut, impédances P à l'échelle sismique ; au centre, les mêmes impédances
P portées dans la grille réservoir ; en bas, les impédances portées de la grille vers le voxet.
50 Chapitre 3. Modélisation sismique à partir de la grille réservoir
Fig. 3.4.2 Statistiques du changement d'échelle entre les diérents supports d'informa-
tion sur le champ Alpha. En haut, nombre d'échantillons moyennés par cellule du réservoir ;
la zone bleue correspond à des cellules de plus grande taille ; en bas, écart-type des échan-
tillons, en pourcentage.
Cette comparaison est d'abord faite sur une trace, au niveau du puits. Les impédances
à l'échelle sismique et à l'échelle réservoir sont transformées en coecients de réexion,
puis convoluées avec une ondelette. Elles sont montrées sur la gure 3.4.3.
3.4. Evaluation des modélisations sur le champ Alpha 51
Fig. 3.4.3 Eet du changement d'échelle sur une trace sismique. (a) l'impédance P de
l'inversion stochastique à l'échelle sismique, (b) l'impédance P dans le voxet, (c) trace
sismique réelle, (d) synthétique obtenu à partir de a, (e) synthétique du voxet obtenu à
partir de b et (f) diérence entre les deux synthétiques.
L'analyse illustrée ici sur le modèle disponible pour le champ Alpha permet de mettre
en évidence les zones où l'évaluation des propriétés du réservoir sera imprécise, voire er-
ronée. On peut voir l'inuence du changement d'échelle sur les propriétés, comme sur les
synthétiques.
Fig. 3.4.4 Eet du changement d'échelle sur les sections sismiques. En haut, la section
issue de la modélisation sur le volume à l'échelle sismique (quelques cellules vides dans
le volume d'impédance sismique apparaissent comme des points sur l'image sismique) ;
en bas, la section issue de la modélisation sur le volume à l'échelle réservoir. Toutes ces
modélisations sont réalisées par convolution 1D sur chacune des traces du volume 3D. Les
traces synthétiques de la gure 3.4.3 ont été replacées sur la gure.
La comparaison qualitative montre que les résultats sont tout à fait acceptables. Les
sections de la gure 3.4.4 sont assez similaires. Les principaux réecteurs sont observés sur
les synthétiques.
Cependant, dans la modélisation dans la grille réservoir, nous n'avons pas d'information
sur les épontes, c'est-à-dire sur les intervalles au dessus et en dessous de la grille. Nous
avons donc un contraste très important sur les synthétiques, pour le premier et le dernier
3.5. Lissage des propriétés du voxet 53
réecteurs.
En raison de la grande taille des cellules de la grille réservoir, nous pouvons observer un
eet de pixellisation sur la section synthétique produite, ce qui peut rendre la comparaison
dicile. Cela est dû aux cellules de la grille réservoir, qui restent visibles dans le processus
de modélisation.
Fig. 3.5.1 Application d'un ltre moyen sur un point donné et ses p traces de part et
d'autre, en x et y.
54 Chapitre 3. Modélisation sismique à partir de la grille réservoir
Lors de la création du voxet, on crée une propriété qui conserve la position et l'index
des couches de la grille réservoir. Le ltre est ensuite appliqué couche après couche. Pour
un ltre de dimension p, en un point donné a(x, y), on calcule la moyenne des points de p
traces de part et d'autre de a, tel que illustré sur le schéma 3.5.1.
En appliquant ce ltre sur les impédances dans les grilles réservoir des champs Alpha
et Bêta, on obtient l'image du bas des gures 3.5.2 et 3.5.3, qui ont plus l'apparence
d'une image sismique, et permettent de faire une comparaison de meilleure qualité avec les
données réelles.
3.5.2 Application d'un ltre moyen sur la section sismique issue des propriétés du
Fig.
PEM du champ Alpha. En haut, la section brute ; en bas, la section ltrée.
3.6. Comparaison des sismiques réelles et synthétiques 55
3.5.3 Application d'un ltre moyen sur la section sismique issue des propriétés du
Fig.
PEM du champ Bêta. En haut, la section brute ; en bas, la section ltrée.
Une approche pour réaliser cette comparaison de manière plus quantitative serait d'utili-
ser les réseaux de neurones comme outil pour classer les traces sismiques synthétiques et
56 Chapitre 3. Modélisation sismique à partir de la grille réservoir
réelles. Nous utilisons pour cela les réseaux de neurones dits "cartes de Kohonen" qui sont
expliqués dans le chapitre 6.
Les cartes de Kohonen sont un type de réseaux de neurones non supervisé qui vont nous
permettre de réaliser une classication sans supervision des amplitudes sismiques. L'appren-
tissage est réalisé sur les données synthétiques, et les neurones résultants sont appliqués
aux données réelles.
La gure 3.6.1 nous montre la comparaison de la répartition des classes du réseau de neu-
rones pour les données réelles et les données synthétiques issues des résultats d'inversion
portés dans la grille réservoir. Cette gure permet de valider cette approche : la répartition
des classes est quasi-identique. La seule diérence visible provient du changement d'échelle
entre les deux volumes sismiques.
Fig. 3.6.1 Comparaison des données sismiques synthétiques et réelles du champ Bêta,
basée sur les cartes de Kohonen, avec 20 neurones. A gauche, la classication des données
synthétiques issues des résultats d'inversion portés dans la grille réservoir ; à droite, la
classication des données réelles.
3.6. Comparaison des sismiques réelles et synthétiques 57
La gure 3.6.2 montre l'application de cette approche sur les données synthétiques
issues des propriétés pétrophysiques obtenues avec le PEM : l'image du bas montre les
zones où un neurone se trouve au même endroit sur les sismiques réelles et synthétiques,
donc les zones où les traces synthétiques sont bien corrélées avec les traces sismiques.
La couleur foncée montre la concordance des classes entre les deux volumes de données.
L'image étant principalement claire, on voit que les données synthétiques générées par le
PEM ne correspondent pas aux données réelles. Dans ces zones, les propriétés du réservoir
doivent donc être modiées.
Cette approche permet de vérier si la sismique synthétique est proche des données réelles
ou pas, sur une carte 2D qui identie les zones cohérentes. La seule diculté est de trouver
un moyen de tirer avantage de ces cartes.
58 Chapitre 3. Modélisation sismique à partir de la grille réservoir
Fig. 3.6.2 Comparaison des données sismiques synthétiques et réelles, basée sur les cartes
de Kohonen, avec 20 neurones. A gauche, la classication des données synthétiques issues
du PEM ; à droite, la classication des données réelles ; en bas, la concordance des classes :
bleu foncé si la classe est la même sur les deux volumes, bleu clair sinon.
3.7. Conclusions et perspectives 59
Les études en caractérisation réservoir utilisent les données disponibles pour construire
un modèle réservoir et déterminer les propriétés qui le caractérisent. Le modèle utilisé
correspond à la vision la plus probable du réservoir, alors qu'il est possible de construire
plusieurs modèles réservoirs à partir des mêmes données. La gestion des incertitudes ap-
paraît à chaque étape de la vie du réservoir. En caractérisation réservoir, ces incertitudes
englobent les paramètres géométriques (taille du réservoir, emplacement des horizons, des
failles, etc.) et les paramètres pétrophysiques (perméabilité, porosités, etc.). Mais le plus
souvent, les incertitudes que l'on regarde en caractérisation réservoir ne correspondent
qu'aux incertitudes sur les propriétés.
Nous cherchons à remplir la grille réservoir de propriétés pétrophysiques dans le but de
simuler son comportement pendant sa production. Mais sans la prise en compte des incer-
titudes liées aux paramètres géométriques, ou incertitudes structurales, nous ignorons si
les résultats des diverses études sont correctement positionnés dans l'espace.
Les incertitudes structurales sont les sources d'erreurs les plus diciles à évaluer lors de la
construction de la grille réservoir. Pour illustrer l'eet de ces incertitudes, nous décrivons
ici la cas de la grille réservoir du champ Alpha, qui montre un décalage entre les horizons
pointés sur la sismique et les horizons correspondant de la grille réservoir. Un simple re-
port en temps ne sut pas à corriger ces diérences. Dans le cas présent, les diérences
observées sont dues à une confusion dans le modèle de vitesse à utiliser. Il résulte de ces
diérences que l'interprétation géologique portée dans la grille montre un chenal près de
20ms au dessus de sa position réelle.
Dans le cas du champ Alpha, l'interprétation géologique était déjà connue, ce qui a per-
mis d'identier le problème rapidement. Mais que se passe-t-il lorsque l'on commence les
travaux de caractérisation réservoir, et que l'interprétation n'est pas encore nalisée ? Com-
ment ces incertitudes vont-elles aecter notre vision du réservoir ?
Les incertitudes structurales doivent être prises en compte dans toutes études de caracté-
risation réservoir. Elles doivent au moins être évaluées avant de choisir un unique modèle
de réservoir. Pour montrer l'eet que peuvent avoir ces incertitudes sur le remplissage du
réservoir, nous proposons d'utiliser le champ Alpha pour réaliser une étude de la répartition
des faciès dans la grille réservoir par classication des attributs sismiques, à laquelle nous
ajoutons une étape de réalisation stochastique de la géométrie du modèle réservoir.
61
62 Chapitre 4. L'impact des incertitudes structurales
4.1 Introduction
Reservoir characterization studies use all available data to build a reservoir model and
determine its properties. The model used is the more probable vision of the reservoir, even
though it is possible to build several reservoirs models from these data. Uncertainty ma-
nagement appears at each stage of the life of the reservoir. In reservoir characterization,
these uncertainties include geometric parameters (size of reservoir, location of horizons,
faults, etc.) and petrophysical parameters (permeability, porosity, etc.). But in most cases,
the uncertainties that we use in reservoir characterization reect only the petrophysical
properties.
We want to inll the reservoir grid with petrophysical properties in order to simulate its
dynamic behaviour during the production. But if we do not take into account the uncer-
tainties associated to geometric parameters (or structural uncertainties), we do not know
if the results of various studies are properly positioned in space.
One question we need to ask is : How will these uncertainties aect our view of the reser-
voir ?
In the case of the Alpha eld presented hereafter, the geological interpretation is already
known. We therefore propose to use the Alpha eld to study the impact of uncertainties on
the structural interpretation of the reservoir (Neau et al. 2007). The dierences of position
observed on this dataset will be used to calibrate the uncertainties for this process.
In the example of Alpha elds, the uncertainties are related to the velocity model for the
time/depth conversion. Other types of uncertainty can lead to misinterpretation (Thore
et al. 2002). The most important ones, in addition to uncertainties related to the time/depth
conversion, are as follows :
Uncertainties related to seismic picking :
The character of a seismic reector can vary laterally, induced by changes in static
or dynamic properties, by variations in facies or types of uids, or by the presence
of multiples.
Uncertainties related to fault zones :
A reector can be degraded by the presence of a shadow zone, or may be lost in a
grinding zone.
All these uncertainties are the most dicult sources of errors to assess during the
construction of the reservoir grid. To take into account these uncertainties, we propose a
typical study of facies repartition in the reservoir grid by classication of seismic attributes,
to which we add a step of stochastic simulation of reservoir models (see workow 4.1.1).
Fig. 4.1.1 Workow to evaluate the structural uncertainty eect on reservoir characte-
rization. Seismic and well data are used to create a base case model. Multiple structural
realizations are simulated using interpretation uncertainties. The rest of the study consists
in inverting, then classifying the attributes to recover facies in the dierent reservoir grids.
equiprobable simulations of geometry of horizons and faults from their reference positions.
The P-Fields method enables to accommodate the use of geostatistical methods and the
available data to characterize the structural uncertainty in a reservoir model.
The theoretical formulation for 3D P-Fields is the following (see gure 4.2.1) : consider
a grid at the nodes of which an attribute z must be simulated. At the grid node x, the
range of possible values is given by the local conditional cumulative distribution function
(CCDF) :
F (x, z) = P (Z(x)) ≤ z|N (x)) (4.1)
Where N (x) ⊂ {z1 , . . . , zN } are the data available over the entire grid and N (x) are the
local conditioning information. Once the attribute has been generated, each simulated value
64 Chapitre 4. L'impact des incertitudes structurales
zsim at the location x may be obtained by computing the inverse of the local CCDF :
The values of p are considered as a realization of a uniform random function P (x), simu-
lated with a given variogram.
Fig. 4.2.1 Operating principle for the P-Fields simulation, applied on a surface. An
uncertainty value zsim (u) is added up to the initial surface at each node u. This value is
drawn in the local CDF. From Thore et al. 2002.
Fig. 4.2.2 Inuence of the correlation length on the simulated horizons. The dotted red
lines represent the uncertainty interval, the black line the reference horizon. The blue and
green lines are the simulated horizons respectively for a small and a long correlation length.
The simulation is constrained by two wells.
Note on well integration At the well location, the exact depth of a marker is known,
and the uncertainty at this point is null. Around the well location, the nearby points (at
one correlation length away) are sensitive to the constraining points and should have lesser
uncertainties. This introduces a bias in the original structural uncertainties of these points.
Note on model coherency As cited previously, the generation of new surfaces should
be done by maintaining the relationships between horizons and faults. A correlation factor
between horizons is introduced in the process. It denes how structural uncertainties from
the reservoir top will aect the other horizons of the model. Indeed, when a simulation
is done on the reservoir top, the transformations are transmitted to the other surfaces
depending on this correlation factor. Then, small perturbations are applied individually
on each horizon, so that these perturbations, added up to the uncertainties from the top
horizon, are within the local uncertainty range. For example, a correlation factor of 80%
means that the uncertainty applied on a horizon results from 64% of deformations on the
reservoir top and from 36% of local perturbations.
Fig. 4.2.3 Misplaced grid layers versus references horizons, overlaid on Alpha Field seismic
data. Reference horizons are blue, corresponding grid layers are grey. Dierences on the
far right-hand side should not be taken into account since the reservoir grid is bigger than
the geological formation.
Horizon Names Grid Layers Raw Uncertainty Final Uncertainty Correlation Factor
(in ms) (in ms) (in %)
E Top 1 23 34,5 100
F Top 20 28 42 80
F Base 51 30 45 80
G Base 77 28 42 80
H Base 101 40 60 70
I Base 112 41 61,5 70
Tab. 4.1 Reference horizons and corresponding grid layers, with their uncertainties and
correlation factors.
We also add the well, in the higher part of the anticline, to constrain the simulations
at the marker locations. We dene the correlation length as a variogram controlling the
spatial variability of the new horizons. The reference horizons picked on seismic data are
quite monotonous so we can take a long correlation length. The variogram used in the
study is derived from attribute analysis on seismic data :
Lx = 1800 m
Ly = 750 m
Orientation = 75◦ from North.
To characterize the simulated models, the rock volume is estimated between the reser-
voir top and bottom horizons. By creating many models, we can see the density functions
4.2. Multiple reservoir grid realizations 67
Fig. 4.2.4 Computed dierence maps between reference horizons and corresponding grid
layers, with their respective histograms. Dierences on the far right-hand side should not
be taken into account since the reservoir grid is bigger than the geological formation.
68 Chapitre 4. L'impact des incertitudes structurales
for the rock volume behaviour ; this approach enables the use of statistical analysis for the
models (Samson et al. 1996).
Note on rock volume computation The rock volume computation impose that the
reservoir model is dened in depth. The rest of the study consists in an inversion process
and a supervised classication, all operating in the time domain. Generating the structural
models in depth, then converting them in time would imply that we add up the structural
uncertainties from the time/depth conversion. To avoid this, the rock volume computation
was adapted so that the unit is not m3 , but m2 ∗ms. This way, we can generate the dierent
structural models directly in time.
Fig. 4.2.5 Statistical analysis from the structural simulation. The top panel gives for
each realization the rock volume. The bottom left panel shows the triangular PDF and the
right panel shows the CDF. The red marker on each panel corresponds to the base case
model.
We also take the initial structural model (percentile 81%). All those models are shown
on gure 4.2.6. No aberration can be seen on the simulated horizons of each model. The
horizons have a monotonous behaviour and respect the constraining markers at the well
location. All those models are good candidates for the rest of the study.
Fig. 4.2.6 The eight structural models on which the rest of the study will take place.
Respectively, the initial model, percentiles 0%, 5%, 10% (srt line), 50%, 90%, 95% and
100% (second line) are represented.
Moreover, this technique gives multiple equiprobable attribute realizations, reecting the
non unicity of the solution.
The Gaussian posterior PDF P (m|s) is constructed by combining the Gaussian prior PDF
P (m) and equation 4.4. Escobar et al. 2006 have developed an SGS-type algorithm to
sample the posterior PDF : it performs a trace-by-trace decomposition of the global PDF
into a number of approximate local PDFs, which are conditioned by previously visited
traces. Multiple realizations can therefore be eciently generated by sampling these local
Gaussian PDFs.
Seismic modeling for the acoustic case is performed with a convolutional model :
s(t) = w(t) ∗ r(t) = Gθ m (4.5)
with w the seismic wavelet and r the reectivity series in time, sampled irregularly because
of the grid layout. The reectivity for a given interface is based on impedance contrasts :
∆z ∆ln(z)
r= ≈ valid for small reectivity variations (4.6)
2z 2
Seismic modeling for the elastic case is performed by computing the angle-variant re-
ectivity with Fatti approximation of Zoeppritz equations :
1 Is
rθ ≈ 2
∆ln(Ip) − 4 sin2 θ∆ln(Is) (4.7)
2cos θ Ip
For the prior impedances of the inversion, we need to inll the models with geological
information. These prior impedances are computed by ltering the well logs (density, P
and S velocities) to keep the low frequency content, then 3D kriging is performed with
these ltered logs (see gure 4.3.1).
Fig. 4.3.1 Original well logs in dark blue and ltered logs with a gaussian lter in green,
respectively for P and S velocities and density.
As previously said, each impedance column is simulated from local PDF with standard
deviation depending on assumed uncertainties for each kind of data (seismic, well, prior
model).
For the prior model, a constant uncertainty of 15% is applied on seismic impedances. This
means that for each cell of the stratigraphic grid, impedances to be computed will have
values between plus and minus 15% around the prior model values in logarithmic domain,
following equation 4.6.
A correlation coecient c = 0, 7, observed on the well logs has been applied as constraint
between P and S impedances. A degree of uncertainty is also applied on the well logs, to
account for the calibration of well data with the seismic. The better the well-to seismic
calibration, the lower this degree of uncertainty will be. Here, we choose 5%, due to the
good quality calibration.
The uncertainty for the seismic data is integrated by making assumptions on signal quality,
meaning the signal to noise ratio. A constant value of NS = 6 to dene the noise statistics
for each trace (the component CsI,θ of equation 4.4) during the inversion process.
Then we dene variograms to control the behaviour of impedances well beyond the seismic
72 Chapitre 4. L'impact des incertitudes structurales
resolution. The vertical component of these variograms is dened by the statistical analysis
of well logs. We obtain a range of 18 ms, meaning 18 layers from the stratigraphic grids.
The lateral component is determined from an exponential variogram with a range of 1000
m, computed from seismic time-slice analysis.
Fig. 4.3.2 Inversion results at well location for P (left) and S (right) impedances. Dark
blue curves are the original P and S impedances.
The rst criterion for quality control of GeoSI inversion is to verify the coherency bet-
ween actual and simulated impedances at well locations. Figure 4.3.2 shows the P and
S mean impedances from the inversion process versus actual impedances from well logs.
Here the correlation between inverted and actual logs is very good, with the inverted P
impedances having a slightly higher variance than the actual P impedance log.
The second criterion is the synthetic seismic obtained from inverted impedances and their
respective residuals with the actual seismic data. Figure 4.3.3 shows an example of this
on the initial model, with the actual and synthetic seismic volumes and the residuals. We
see that the inverted impedances are in good agreement with the actual data, with overall
RMS error of 4%.
The gure 4.3.4 shows an independent realization and the mean from all realizations.
4.3. Stochastic inversion 73
Fig. 4.3.3 Inversion results for the initial model. Top : Actual seismic data, center :
synthetic seismic data and bottom : residuals.
74 Chapitre 4. L'impact des incertitudes structurales
Fig. 4.3.4 Inversion results for the initial model : an independent realization of P impe-
dances versus the mean. The actual seismic is represented in the background, as well as
the P impedance log.
space property and lithologies from a training set. Within each lithologic group,
the property is assumed to have a Gaussian distribution with conditional PDF.
The centroid and covariance matrix of the Gaussian density are inferred from a
training population of samples that should be available for each lithotype. Bayesian
framework is used to predict the group membership of a sample.
Discriminant analysis (Gnadesikan et al. 1989) :
Discriminant analysis looks for functions of the properties that optimally (in a least-
squares sense) separate the groups. These functions are linear functions of the proper-
ties and are calculated from a training data set of samples with known membership.
The purpose of the method is to select a combination of linear coecients that maxi-
mize the dispersion of group centroids yet minimize the dispersion of cases within
groups. This is illustrated on gure 4.4.1.
This method assigns to each trace a membership probability for a given lithology,
depending on the attributes. As it is in Nivlet et al. 2007, a trace x(t) = (x1 , . . . , xp )
with p attributes and {C1 , . . . , CN } the N probable facies, is assigned a set of N
probabilities p(Ci /x). The goal being to assign to this seismic trace one facies from
{C1 , . . . , CN }, this amounts to :
p(x|Ci )p(Ci )
p(Ci |x) = PN (4.8)
j=1 p(x|Cj )p(Cj )
with p(Ci ) the a priori probability for facies proportions and p(x|Ci ) the training
PDF for a given facies.
Fig. 4.4.1 Distribution of samples from two lithotypes in a space of two physical pro-
perties z1 and z2 . The probability densities for each group are represented near to the
axes, showing signicant overlap between both groups. The discrimination is improved
with linear combinations of the original variables. From Bosch et al. 2002.
76 Chapitre 4. L'impact des incertitudes structurales
We are going to use this last method, the discriminant analysis, since the results with
this method have a higher success rate than with other methods (Bosch et al. 2002).
Before performing the facies classication, we need to address the dierences between well
and seismic data. The easiest way to integrate geological information is to constrain the
classication with the facies interpretation of well-log data. The attributes from the well
data are used directly as a training set to calibrate the classication. Then the operator is
applied on seismic attributes to discriminate the dierent seismic facies. However, we must
be very careful with this approach since well and seismic data are not dened at the same
scale and resolution.
To address these issues, we should test the discrimination at well-log scale and then at
seismic scale.
Fig. 4.4.2 Raw crossplot of Ip versus P R coloured by the 7 facies interpreted on the logs.
Figure 4.4.2, which displays the crossplot P impedances log (Ip) versus Poisson Ra-
tio log (P R), coloured with the 7 geological facies, shows a good discrimination between
each facies. The overall overlap (proportion of facies dierent from the dominant one in an
elementary cell of the crossplot) for this crossplot is around 33%. Facies 1, 6 and 7 (respec-
tively shales, debris ows and water-sands) can not be discriminated in the crossplot. To
improve the overall discrimination in the crossplot, we group them together in the facies
code F1. The other facies codes are : F2. Laminated oil sands ; F3. Fine oil sands ; F4.
Coarse oil sands ; F5. Cimented oil sands.
With this modication, the overlap is down to 15% (see gure 4.4.3).
4.4. Facies classication 77
One of the eects we can see on this crossplot (especially on F1) is the compaction
eect, which results in a global increase of impedances with depth. In the Alpha eld, this
eect is remarkable due to the anticline structure of the eld and to the complex history
of deposition, which results in multiple reservoir intervals.
To attenuate this eect, Ip and P R have been detrended by computing trend curves Ip =
f (z) and P R = f (z). Figure 4.4.4 displays the new discrimination between the facies, with
an overlap of 11%. Quantitatively, the overlap score increases to 27%, which is essentially
due to a lesser dissemination of F1 points. The nal overlap score after detrending is
satisfactory for going on with the analysis.
Fig. 4.4.4 Crossplot of Ip versus P R coloured by facies, after removing the compaction
eect from both properties.
78 Chapitre 4. L'impact des incertitudes structurales
The clusters identied for each facies are not composed of 100% of this facies. If we
analyze the composition of each cluster, we obtain the pie-charts shown in gure 4.4.5.
When we will compute the lithocubes, each sample will be represented by a pie-chart
giving the proportion of all facies.
Fig. 4.4.5 Crossplot of Ip and P R detrended, with polygons delimiting the clusters for
each facies. The composition of each cluster is shown in the pie-charts (from left to right,
Facies 1 to Facies 5).
Seismic data are dened at a coarser scale than well logs. The discrimination observed
at log scale should be validated at seismic scale to ensure that the discrimination between
facies is still viable. To do so, the log data are ltered in the seismic bandwidth with a
band-pass lter.
The quality control of the nal discrimination consists in comparing the lithofacies log
at log scale and the simulated seismic facies log built from the attributes. Figure 4.4.6
compares initial facies log with facies log predicted from the attributes at seismic scale.
The two logs are quite similar, meaning that the main reservoir heterogeneities can be
properly identied from detrended P impedances and Poisson Ratio at the seismic scale.
4.4. Facies classication 79
Fig. 4.4.6 Comparison between initial facies interpretation of well log data (left), with
the facies interpretation of detrended properties at the seismic scale (right).
At well position, even though mists between well-logs and inversion results are ac-
ceptable, they are sucient to create inconsistencies in the seismic facies model. We have
therefore recalibrated the discrimination between seismic facies and detrended properties
based on the following data : a minicube composed of a 1x1 corridor traces extracted
around well positions from seismic attributes and the seismic facies being obtained from
the interpretation of well-log data at the seismic scale. For the adjoining traces, we suppose
that the facies associated to the sample ti at well location stays the same for the samples
ti of neighbouring traces in the minicube.
Figure 4.4.7 compares the cluster repartition between logs, ltered logs and seismic attri-
butes. This gure shows a global good agreement between the three analyses, meaning that
this discrimination can be applied to the seismic attribute cubes.
Fig. 4.4.8 PDF operator in the crossplot of detrended Ip versus detrended P R coloured
by facies : left : original well logs, center : ltered logs in the seismic bandwidth, right :
minicube from seismic attributes.
Figure 4.4.9 shows a slice of the actual seismic and the ve facies probability cubes
computed from the PDF operators, through the main channel of reservoir R-Alpha. Pro-
babilities associated to this interpretation are high ; they are slightly lower in the channel
sequences, where geology changes rapidly. In conclusion, 3D seismic facies analysis identies
qualitatively the main geological features of Alpha eld.
Fig. 4.4.9 Probability cubes computed from the PDF operators. The rst image shows
a slice of the actual seismic block through the main channel of reservoir R-Alpha, along a
horizon from the initial structural model. The other images show a slice of the probability
cubes for each facies. Facies F1 range is between 0 and 100% ; for the other facies codes,
the range is between 0 and 25%.
82 Chapitre 4. L'impact des incertitudes structurales
variable. Indeed, by upscaling the facies cube in the reservoir grid, we would have the pro-
blem of choosing which facies should be integrated in a cell, as it is shown on the example
described before. Direct upscaling would result in heterogeneity or facies connectivity loss
(channel continuity for example) (see Nivlet et al. 2007).
Hadj-Kacem and Pivot 2005 and Duplantier et al. 2006 describe a new methodology to
associate facies in a reservoir model. They propose to use the most probable facies cube
to integrate facies probabilites in the reservoir model, and thus, to avoid any bias in the
probabilities due to abnormal values. The remaining values are averaged using classical pro-
cedures. Here is an example : a reservoir grid cell containing height samples from the facies
probability cube at the seismic scale, which values are : R = {25; 32; 25; 51; 37; 25; 48; 25}.
This gives : arithmetic mean Ra = 33, 5, geometric mean Rg = 32, 12 and harmonic mean
Rh = 30, 92. From a statistical point of view, we should have a facies probability greater
than 30%. However, if we use the most probable facies method, we have four occurences
of facies probability 25% on eight, meaning a 50% probability, so we should have a 25%
probability for the considered facies in this reservoir cell. By applying this methodology on
each reservoir model, we upscale facies or probability cubes at the reservoir scale. Figure
4.5.1 shows the F1 probability cube upscaled in the reservoir grid, with the most probable
method.
Fig.4.5.1 Slices of the F1 Probability cube along the main channel of R-Alpha and the
F1 most probable value in the initial reservoir grid.
We will use this methodology on the probability cubes, meaning for each facies proba-
bility cube, we will take the most representative value inside the reservoir cell. Then we
create a facies cube, considered as the best scenario possible for facies repartition in the
4.5. Facies repartition in the reservoir grid 83
reservoir grid. For that, we use the dominant facies method. This method for building the
facies property should account for cases where there are dominant facies proportion. This
is the case for the reservoir R-Alpha, where shale proportions are very important, as shown
by the F1 probability cube (the average proportion is greater than 70%). We can dene a
threshold for these facies proportions. In the event the proportion of shale falls below the
threshold, one of the sand facies will take over in the cell. This threshold is determined by
tying the interpreted facies log and the seismic facies log derived from the seismic attributes
(shown in gure 4.4.6.
This calibration step is important : if the shale threshold is too high, the sand facies will
be over-represented ; if the threshold is too low, the sand facies will be under-represented.
This is illustrated on gure 4.5.2. For example, in the main channel of reservoir R-Alpha,
the shale threshold is around 87%. This means that for the interval of this channel, if the
F1 probability exceeds 87%, we will found shales ; otherwise, we will nd one of the sand
facies (see gure 4.5.3).
However, since the reservoir R-Aplha is in fact constituted in 5 reservoir intervals, the shale
proportion diers from one reservoir interval to the next (main channel threshold : 87% ;
last reservoir interval : 75%). Therefore, the calibration is done with dierent threshold
values on the whole reservoir.
We do these upscaling and dominant facies generation steps for each structural model.
Fig. 4.5.2 Facies repartition in the initial reservoir model, on layer 21, with dierent
shale thresholds in percent.
84 Chapitre 4. L'impact des incertitudes structurales
Fig. 4.5.3 Facies repartition in the initial reservoir model, on layer 21, with the shale
threshold of 87%, determined from the well data calibration.
The layers around the 21st and 105th layers were also investigated in order to see if the
expected geological bodies were recovered on neighbouring layers, but this was not the case
4.5. Facies repartition in the reservoir grid 85
Fig. 4.5.4 Facies repartition in the 21st layer of the structural models. The rst line
shows the geological interpretation, the picture on the left of the second line shows the
facies repartition in the initial reservoir grid.
86 Chapitre 4. L'impact des incertitudes structurales
Fig. 4.5.5 Facies repartition in the 105th layer of the structural models. The rst line
shows the geological interpretation, the picture on the left of the second line shows the
facies repartition in the initial reservoir grid.
4.6. Conclusion 87
4.6 Conclusion
An important issue in understanding heterogeneous reservoirs is determining the rela-
tive importance of structural and petrophysical variability. Structural uncertainty is not
typically incorporated into the construction of a reservoir model, unlike petrophysical un-
certainty.
We have presented a simple methodology to evaluate their eects in a reservoir model.
This methodology is based on using 3D reservoir grids. Multiple realizations of structural
models have been simulated and used in a standard reservoir modeling workow. For each
reservoir grid, we obtain the facies repartition and systematically investigate it. Each model
shows dierent facies layouts with diverse possible geological interpretations.
The results further reveal that evaluating the structural uncertainty eects in reservoir
grids is an important step before choosing the denitive model to be used in future studies.
By visualizing the structural uncertainty, we can gain a greater appreciation of the re-
liability of any analysis performed in reservoir characterization studies. Structural uncer-
tainties are quantiable and can be carried through a reservoir characterization workow
by building several reservoir models. The example shown here with the Alpha eld de-
monstrates how structural uncertainty can lead to dierent facies repartition and to wrong
interpretation in the reservoir models.
Partie II
Alternatives à l'approche
directe
89
Chapitre 5
5.1 Introduction
Une étude réservoir a pour objectif nal la création d'un modèle nous permettant de
prédire dans le futur le comportement du réservoir au cours de sa production. Pour prédire
les propriétés pétrophysiques qui nous permettront de faire cette simulation, nous cher-
chons à intégrer autant d'information que possible dans le réservoir.
Les méthodes traditionnelles de simulation réservoir créent un modèle de lithotype (envi-
ronnement + faciès géologiques + propriétés pétrophysiques) :
par des techniques "objet" : processus booléens, modélisation d'objets longs, pro-
cessus génétiques ;
par des techniques "pixel" : simulations gaussiennes, simulation d'indicatrice, sta-
tistiques multipoints.
Les modèles réservoirs résultants de ces processus sont conditionnés aux puits, mais laté-
ralement, nous n'avons aucune certitude sur leur validité.
Les méthodes traditionnelles de caractérisation réservoir font appels aux données sismiques
pour contraindre le réservoir. Les données sismiques sont la seule source d'information qui
couvre l'intégralité du réservoir. Elles sont l'expression des contrastes des propriétés pétroé-
lastiques aux limites des couches géologiques. Les valeurs réelles des propriétés ne sont donc
pas connues dans l'absolue. Un second problème est que les attributs dérivés des données
sismiques sont à une échelle diérente de celle du réservoir, et il faut donc soit eectuer
un changement d'échelle de ces attributs (upscaling ), avec tous les problèmes associés à ce
changement d'échelle (Ringrose 2007), soit simuler les propriétés réservoirs en se servant
des attributs comme variables secondaires.
Nous proposons une nouvelle méthode d'inversion à mi-chemin entre les disciplines simula-
tion réservoir et caractérisation sismique du réservoir. Ce processus d'inversion permet de
générer des attributs sismiques directement à l'échelle réservoir. Il est basé sur le couplage
de deux techniques :
les déformations graduelles, pour la génération de réalisations géostatistiques d'at-
tributs ;
la modélisation sismique, pour l'optimisation à partir des données sismiques.
91
92 Chapitre 5. Inversion par déformation graduelle
Les déformations graduelles ont donc été développées pour dépasser les problèmes de temps
d'exécution d'une part, et la dimension du réservoir d'autre part. Les déformations gra-
duelles sont donc d'un outil capable de modier les réalisations de propriétés réservoirs,
tout en maintenant les contraintes aux données géologiques et en reproduisant un com-
portement connu (les données dynamiques pour l'utilisation courante des déformations
graduelles, les données sismiques dans l'application décrite ici), tout en ayant un temps
d'exécution raisonnable, et s'appliquant à tous types de modèles, quelle que soit sa taille.
En utilisant un mécanisme simple de perturbation des réalisations, la méthode des défor-
mations graduelles cherche dans l'espace du modèle a priori les solutions qui honorent les
données que l'on veut inverser. Ce mécanisme assure le respect des statistiques du modèle
(Caers 2007).
Les méthodes statistiques sont largement utilisées pour modéliser des propriétés réservoirs
distribuées dans l'espace. Dans de nombreux cas, il n'y a pas susamment de données
pour en déduire avec précision les paramètres structuraux d'un modèle stochastique. Ces
paramètres, qui dénissent la géométrie des corps, sont souvent exprimés par des distri-
butions a priori et des variogrammes. La simulation doit donc inclure un paramètre sur
lequel repose tout cet aspect structural, découplé de la simulation de la distribution de la
propriété. Prenons L l'opérateur covariance et X un bruit blanc gaussien standard. On a
alors :
Yn = L [Xn ] (5.3)
94 Chapitre 5. Inversion par déformation graduelle
Ceci nous donne une réalisation gaussienne Y dont la covariance impose la structure
spatiale de la réalisation. Appliqué à la formulation des déformations graduelles, nous
avons :
Yn (L, t) = L[Xn (t)] = L [Yn−1 .cos(t) + Un .sin(t)] (5.4)
Sous cette forme, on peut donc appliquer les déformations graduelles sur les réalisations
tout en modiant la fonction de covariance (voir gure 5.2.1). Cette formulation correspond
aux déformations graduelles structurales.
Dans d'autres cas, il peut être utile de ne modier qu'une partie du modèle réservoir.
En eet, pour un réservoir déjà étudié, des changements locaux, liés à des données récem-
ment acquises, peuvent s'avérer nécessaires. Reprendre l'étude du réservoir depuis le début
serait long et fastidieux, et le fait de réajuster une zone par le processus global pourrait
dégrader la qualité d'une autre zone. Il existe une formulation des déformations graduelles
qui permet d'appliquer des modications locales dans le modèle réservoir. Ces modica-
tions ne peuvent s'appliquer que si la réalisation secondaire prend elle-même en compte
l'aspect local de la modication.
Ravalec-Dupin et al. 2000 ont donc proposé une formulation particulière de simulation
FFT, appelée FFTMA (pour Fast Fourier Transform Moving Average). Cette formulation
permet de considérer X le bruit blanc de tout le modèle réservoir par {X1 , . . . , Xn } la
partition de X en n zones indépendantes. Chaque zone peut donc être manipulée indi-
viduellement. Il se pose alors un problème de continuité spatiale entre les zones. Mais les
déformations graduelles locales s'appliquent uniquement sur le champ de bruit blanc. Ainsi,
en appliquant la fonction de covariance, sur le bruit blanc modié, la continuité spatiale
de la propriété est préservée (comme montré sur la gure 5.2.1).
Fig. 5.2.1 Les diérentes formulations des déformations graduelles. Les images du haut
montrent une chaîne de réalisation réservoir à partir de l'équation de base des déformations
graduelles. Au milieu, la formulation des déformations graduelles structurales (la covariance
change, mais pas la distribution de la propriété). Les images du bas montrent la formulation
locale (seule une petite partie du réservoir est réévaluée). De Hu et al. 1999.
z est un bruit blanc gaussien et g vient de la covariance par C = g ⊗ g̃ avec g̃(x) = g(−x).
g étant une fonction analytique complexe, nous passons dans le domaine spectral pour que
le produit de convolution se simplie en produit de multiplication.
C'est le cas de la méthode développée par Ravalec-Dupin et al. 2000 et appelée Trans-
formée de Fourier (TF) Rapide à Moyenne Mobile (ou Fast Fourier Transform - Moving
Average : FFT-MA). Cette méthode fait le calcul y = µ + g ⊗ z dans le domaine spec-
tral, ce qui a comme premier avantage de transformer le produit de convolution en simple
multiplication, soit : T F (g ⊗ z) = G.Z . La forme générale de la matrice de covariance est :
Z +∞
C(x) = S(f ).e2iπfx dx (5.6)
−∞
Et sa transformée inverse :
Z +∞
S(f ) = C(x).e−2iπfx df (5.7)
−∞
de la variance de Y est :
¯ ¯ ¯ ¯
h i h i
E Y Ý = E G.Ǵ.Z.Ź = G.Ǵ (5.8)
L'espérance mathématique d'un bruit blanc est zéro. De plus, le théorème d'orthogonalité
permet de dire :
¯
h i
E Y Ý = 0 si f 6= f´ (5.9)
¯ ¯
h i
E Y Ý = G.Ǵ si f = f´ (5.10)
Pour générer la réalisation gaussienne, il faut procéder suivant les étapes décrites dans
le tableau 5.1.
1 Construire la fonction de covariance C (en se servant du variogramme)
2 Générer le bruit blanc z
3 Calculer les TF de C et z pour obtenir S et Z
4 Dériver G de S
5 Faire le produit G.Z
6 Faire la TF inverse de G.Z , pour obtenir g ⊗ z
7 Dériver y par l'équation 5.5
Tab. 5.1 Etapes permettant la construction d'une réalisation gaussienne par la méthode
de simulation FFT-MA.
graduelles, Y est conditionnée aux données par la relation suivante (Chiles and Delner
1999) :
Yc (u) = Ydk (u) + [Y (u) − Yk (u)] (5.14)
avec u le point objectif, Yc la réalisation conditionnée, Ydk l'estimation par krigeage des
données connues et Yk l'estimation par krigeage de la réalisation non conditionnée. Les
estimations par krigeage Xk sont réalisées par :
n
(5.15)
X
Xk (u) = λi (u)X(ui )
i=1
avec λi le poids associés à chaque point de données ui , que l'on dérive de Cλ = B , où C est
l'opérateur covariance (le même que celui utilisé pour les déformations graduelles) et B est
un vecteur de corrélation spatiale entre la valeur à estimer en u et l'observation aux point
ui . Cette formulation nous permet de ramener la réalisation non conditionnelle à zéro et
d'ajouter ensuite les données connues, tout en respectant la structure spatiale (voir gure
5.2.2).
Fig. 5.2.2 Principe du krigeage des réalisations non conditionnées dans le processus
des déformations graduelles. La ligne 1 représente la simulation non-conditionnée avec le
point de données et l'estimation par krigeage du point ; la ligne 2 montre la simulation
conditionnée.
hyperellipse de dimension N .
Pour illustrer ce phénomène, prenons un exemple de modèle réservoir à deux cellules, avec
Y1 la réalisation initiale et Y2 la réalisation secondaire, de moyenne Y0 et covariance iden-
tiques :
a1 a2
Y1 = , Y2 = (5.16)
b1 b2
Le principe de déformation graduelle nous donne :
tel que cela a été décrit dans le chapitre 3 page 43. La grille réservoir est tout d'abord ré-
échantillonnée latéralement à l'échelle sismique, en créant un voxet, c'est-à-dire un ensemble
de pseudo-puits à chaque noeud du maillage sismique. Il n'est nécessaire de réaliser cette
opération qu'une seule fois pour tout le processus d'inversion par déformations graduelles
puisque la structure du réservoir n'est pas remise en cause. Les propriétés du réservoir sont
ensuite récupérées le long de chaque pseudo-puits, sans changement d'échelle. Ceci sera
répété pour chaque itération, après combinaison et conditionnement des réalisations. Les
coecients de réexion sont alors calculés, puis convolués avec une ondelette pour obtenir
les traces sismiques. Le volume ainsi obtenu peut maintenant être comparé aux données
réelles. Nous utilisons pour cela une fonction coût basée sur une norme L2 quantiant les
écarts entre données réelles et synthétiques pour chaque échantillon.
1 h sim i2
z(t) = × d (t) − dobs (5.20)
2
Cette fonction coût est minimisée en respect du paramètre de combinaison t par une
méthode de recherche directe. Une grille réservoir est constituée de plusieurs milliers (voir
millions) de cellules. L'optimisation d'un problème d'inversion sur la grille réservoir devrait
donc dépendre de la performance de chaque cellule. Le fait d'optimiser uniquement sur t
nous permet de réduire un problème à N dimensions à un problème à une dimension.
5.3.3 Itération
Comme vu sur le schéma 5.1.1, il existe deux plans d'itération dans le processus d'in-
version par déformation graduelle.
L'itération sur t
Le paramètre t est déni sur [0 ; 2π], car la loi de déformation est périodique : pour
t = 0, la réalisation optimale est égale à la réalisation initiale ; pour t = π2 , la réalisation
optimale est égale à la réalisation secondaire. Il est impossible de simuler cet intervalle
en continu dans les déformations graduelles. Il faut donc discrétiser l'intervalle tel que
t = c × 2π avec c = [0 ; 1]. Pour une réalisation initiale et une réalisation secondaire
données, on peut donc obtenir plusieurs réalisations optimales diérentes, donc plusieurs
valeurs de la fonction coût dont le comportement-type est montré sur le graphe 5.3.2. Ils
existent donc des minima locaux et un minimum global. Tant que ce dernier n'est pas
identié, il faut changer la valeur de t dans la combinaison des réalisations. À chaque
itération, les déformations graduelles améliorent ou au moins maintiennent l'ajustement
aux données réelles.
L'itération globale
Lorsque le minimum global de la fonction coût est trouvé, la réalisation optimale obte-
nue devient la réalisation initiale de l'itération suivante. Une nouvelle réalisation secondaire
est générée puis par combinaison la chaîne de réalisation est créée en fonction de t. Et le
processus global recommence.
Fig. 5.3.2 Comportement type de la fonction coût pendant une itération donnée, en
fonction de c. En c = 0, nous avons la réalisation initiale, en c = 1 la réalisation secondaire.
En c = 0.05, on note un minima local de la fonction coût et en c = 0.95, le minima global
(en rouge) qui correspond à la réalisation optimale de l'itération.
Fig. 5.4.3 Vitesse P du modèle initial sur trois couches diérentes dans le modèle réser-
voir.
5.4. Application sur le champ Alpha 103
Fig. 5.4.2 Application du PEM à l'échelle log : pour chaque faciès (colonne 1) et pour
chaque log, on extrait les lois d'enfouissement. Les logs calculés, à l'échelle réservoir, re-
produisent les logs initiaux.
Tab. 5.2 Relations V p − ρ de Gardner par lithologies pour l'équation ρ = dVpf . Les unités
sont km/s et g/cm3
Les propriétés obtenues sont montrés en gure 5.4.5 et leurs histogrammes sur la gure
5.4.6. À partir de ces propriétés "réelles", nous simulons le volume sismique qui servira de
contrainte pour le processus d'inversion.
Fig. 5.4.5 Densité (à gauche) et vitesse P (à droite) réelles sur la couche 51 du modèle
réservoir.
L'anamorphose n'altère pas les caractéristiques spatiales des données, telles que les
relations d'enfouissement. Cependant, si l'on simule des valeurs sortant de l'intervalle de
départ, la transformation inverse ne peut se faire directement car nous ne connaissons pas
la relation correcte entre les CDF. Ce problème est généralement résolu en faisant des hy-
pothèses sur le comportement de la propriété aux limites de la distribution, telles que des
extrapolations linéaires.
Dans notre procédure, seule la réalisation initiale doit être transformée en données gaus-
siennes. Après combinaison avec la réalisation secondaire (qui elle est gaussienne par
construction avec les FFT), on retourne dans le domaine expérimental, puis on conditionne
la réalisation optimale aux données de puits.
106 Chapitre 5. Inversion par déformation graduelle
Fig. 5.4.7 Processus de l'anamorphose. A gauche : Porosité brute (en haut, histogramme ;
en bas, la fréquence cumulée). A droite : Porosité transformée par anamorphose (en haut,
histogramme gaussien ; en bas, la fréquence cumulée gaussienne).
Les données sismiques synthétiques sont simulées comme expliqué dans le chapitre 3. Un
voxet déformé est d'abord créé à partir de la grille réservoir, puis les propriétés sont trans-
férées d'un support à l'autre. La modélisation est ensuite réalisée avec l'ondelette du champ
5.4. Application sur le champ Alpha 107
La grille du réservoir Alpha est dénie en temps, il n'y a donc pas besoin de réaliser
une conversion temps/profondeur. Si la grille avait été dénie en profondeur, il aurait fallu
convertir la grille en temps avec un modèle de vitesse avant d'eectuer la convolution.
Fig. 5.4.8 Fonction coût pour une inversion par déformations graduelles. Ici, la fonction
coût est calculée sur 300 itérations. La diminution de l'erreur entre réalisations initiales et
nales est de 69%.
D'autres exécutions du processus ont montré des diminutions de la fonction coût d'en-
viron 50% (de 49 à 59%) pour seulement 100 itérations. Les résultats sont assez similaires
à ceux montrés ci-après.
Réalisations réservoirs
La gure 5.4.9 en vitesse P (première ligne) et densité (deuxième ligne) et les gures
5.4.10 à 5.4.13 en impédance montrent les réalisations réservoirs pour le modèle initial
(première colonne), le modèle nal (deuxième colonne) et le modèle réel (dernière colonne)
pour plusieurs couches du modèle réservoir. La gure 5.4.14 montre une section verticale à
travers le réservoir, avec la section initiale (première colonne), la section nale (deuxième
colonne) et la section réelle (troisième colonne).
108 Chapitre 5. Inversion par déformation graduelle
Pour chacune de ces images, on peut voir que les réalisations nales s'approchent d'une
version lissée de la solution. Les grandes structures sont récupérées mais pas les détails les
plus ns.
Fig. 5.4.9 Vitesse P (première ligne) et densité (deuxième ligne) pour les réalisations
initiales (première colonne), nales (deuxième colonne) et réelles sur la couche 51 du modèle
réservoir.
Fig. 5.4.10 Impédances pour les réalisations initiales (première colonne), nales
(deuxième colonne) et réelles sur la couche 20 du modèle réservoir.
5.4. Application sur le champ Alpha 109
Fig. 5.4.11 Impédances pour les réalisations initiales (première colonne), nales
(deuxième colonne) et réelles sur la couche 40 du modèle réservoir.
Fig. 5.4.12 Impédances pour les réalisations initiales (première colonne), nales
(deuxième colonne) et réelles sur la couche 50 du modèle réservoir.
Fig. 5.4.13 Impédances pour les réalisations initiales (première colonne), nales
(deuxième colonne) et réelles sur la couche 80 du modèle réservoir.
110 Chapitre 5. Inversion par déformation graduelle
Fig. 5.4.14 Section en vitesse P pour les réalisations initiales (première colonne), nales
(deuxième colonne) et réelles.
Fig. 5.4.15 Comparaison des histogrammes initiaux, naux et réels pour la vitesse (pre-
mière ligne) et la densité (deuxième ligne).
La gure 5.4.15 montre les histogrammes des propriétés pour les réalisations initiales,
nales et réelles. On peut voir que l'on se rapproche de la distribution réelle des propriétés,
avec une proportion plus élevée des valeurs proches de la moyenne.
Volumes sismiques
Maintenant nous regardons les résultats de l'inversion par déformations graduelles après
transformation des réalisations en volumes sismiques. La gure 5.4.16 montre les réalisa-
tions initiales (première colonne), nales (deuxième colonne) et réelle (dernière colonne),
tandis que la gure 5.4.17 montre les résidus initiaux et naux. Nous voyons que l'inver-
sion par déformations graduelles a largement permis de restituer la sismique réelle, même
si les résidus naux ne sont pas nuls. Les amplitudes sismiques restantes correspondent aux
structures nes du modèle que la méthode n'a pas su récupérer.
Fig. 5.4.16 Sections sismiques avant et après le processus d'inversion par déformations
graduelles, et la sismique "réelle".
Fig. 5.4.17 Résidus avant et après le processus d'inversion par déformations graduelles.
(LG ), petit axe (Lg ) et orientation (en degré) du variogramme (en nombre de cellule de la
grille réservoir). La diminution de la fonction coût est également indiquée.
A l'itération 150 du processus décrit dans les sections précédentes, la fonction coût
montrait une amélioration de 60%. Pour les trois processus successifs, à l'itération 150,
l'amélioration de la fonction coût est de 58,54%. On obtient donc des résultats très sem-
blables en terme de réduction de l'erreur.
Les gures 5.4.18 à 5.4.21 montrent les réalisations obtenues avec ces trois processus suc-
cessifs et les comparent aux réalisations obtenues précédemment.
Ces résultats montrent que des structures de plus petites tailles sont détectées par l'in-
version, comme par exemple sur la gure 5.4.21 où l'on voit une forme ne triangulaire
(entourée en rouge), alors qu'avec un seul variogramme, la réalisation ne retrouvent pas
cette forme.
Fig. 5.4.18 Réalisations en impédances pour les réalisations initiales (première colonne),
nales (deuxième colonne) et réelles sur la couche 20 du modèle réservoir. La réalisation
nale avec plusieurs variogrammes (en bas) est comparée à la réalisation nale à un vario-
gramme (en haut).
5.4. Application sur le champ Alpha 113
Fig. 5.4.19 Réalisations en impédances pour les réalisations initiales (première colonne),
nales (deuxième colonne) et réelles sur la couche 40 du modèle réservoir. La réalisation
nale avec plusieurs variogrammes (en bas) est comparée à la réalisation nale à un vario-
gramme (en haut).
Fig. 5.4.20 Réalisations en impédances pour les réalisations initiales (première colonne),
nales (deuxième colonne) et réelles sur la couche 50 du modèle réservoir. La réalisation
nale avec plusieurs variogrammes (en bas) est comparée à la réalisation nale à un vario-
gramme (en haut).
114 Chapitre 5. Inversion par déformation graduelle
Fig. 5.4.21 Réalisations en impédances pour les réalisations initiales, nales et réelles sur
la couche 80 du modèle réservoir. La réalisation nale avec plusieurs variogrammes (en bas)
est comparée à la réalisation nale à un variogramme (en haut). Les carrés rouges montrent
une structure ne correctement recrée avec l'utilisation de plusieurs variogrammes.
Cela montre bien la nécessité de dénir plusieurs niveaux d'hétérogénéités dans les
simulations, ce qui n'est possible qu'en eectuant plusieurs inversions successives avec dif-
férents variogrammes. Comme la simulation FFT et le krigeage utilisent le même vario-
gramme pour dénir le comportement spatial global des propriétés, une alternative serait
de dénir des variogrammes à structures emboîtées. Le processus travaillerait d'abord à
large échelle, puis à échelle de plus en plus ne
sismiques identiques (en fonction de la taille des cellules) à comparer à des traces sismiques
réelles diéréntes.
6.1 Introduction
The goal of seismic reservoir characterization is to derive reservoir properties from 3D
seismic data with constraints by well logs and a priori geological knowledge. The most im-
portant objective is to nd a relationship between seismic traces and reservoir properties.
Three approaches can be considered.
In the rst one, those relationships are derived through the use of geostatistics, neural
networks, or multi-linear regression. For example, Fournier and Derain 1992 and deGroot
et al. 1993 propose to extract these relationships through statistical calibration of seismic
attributes and reservoir parameters at well location. Chandra et al. 2003 use neural net-
works to perform a non-supervised classication of seismic data. A seismic facies map is
created using unsupervised neural networks, and then the porosity trends by facies are ap-
plied on the map to derive a porosity volume. Their results clearly show the eectiveness of
neural network to identify the stratigraphic features, whereas conventional analysis would
have given less satisfactory results. The use of Neural Networks to link seismic parameters
with well information and predict reservoir properties between wells has been successful in
several instances. But in most cases, the reservoir parameters are the results of a pattern
recognition based on seismic trace shapes. These approaches do not rely on petrophysics
to derive the reservoir properties.
The second approach performs the inversion of the available data and tries to derive the
reservoir properties with calibrated relationships on the seismic attributes, such as Coleou
et al. 2006 or such as the work shown in section 4.4 from chapter 4.
The third approach uses seismic modeling to visualize the seismic responses due to re-
servoir property changes. Julien et al. 2002 have created an application called "Massive
Modeling", which creates a large amount of models with a realistic range of structural and
petrophysical variations in the reservoir.
117
118 Chapitre 6. Petrophysical inversion by neural supervised classication
In unsupervised neural network, seismic data are classied without a priori information.
The results respect logical or mathematical criteria that cannot be evaluated in terms of
petrophysical properties. In supervised neural nets, the classication is constrained by the
seismic response of a priori models. This approach is limited by the amount of data avai-
lable.
We use the advantages of both neural network approaches (unsupervised and supervised
neural nets) : the logical classication of data in the unsupervised classication, and the
constraints from known relationships in the supervised classication. We use a non super-
vised neural network that we apply rstly on a training set, and then we impose the output
neurons on the actual data. The originality of our methodology lies in the generation of
the training set (Massive Modeling and geostatistical simulation) and the petrophysical
discrimination analysis.
The methodology is tested on two dierent datasets : a clastic case (Beta eld) and a
carbonate case (Gamma eld).
The workow for the global methodology is illustrated on gure 6.2.1 and consists in
the following steps :
1. Well log data preparation.
Well log blocking.
Optimization of blocked well logs.
Massive Modeling data generation.
2. Neural network training.
Network learning phase on the training dataset
Validation of petrophysical meaning of the classication
3. Neural network application.
Network validation phase on the test dataset.
Computation of petrophysical property volumes.
The purpose of the well log preparation step is to decrease the number of parameters
one will have to handle, and to generate all possible 1D petrophysical models. For each
well, the logs are blocked, then optimized in regards to the seismic trace at the same lo-
cation. Once available, these "best" logs are input to the Massive Modeling scheme. New
6.2. The methodology 119
petrophysical models are generated through a SGS algorithm, based on the parameter and
thickness variations.
In the neural network training step, the Massive Modeling dataset is analyzed by the neural
network, to obtain model traces which explain the trace shape diversity. We verify that
this classication has a petrophysical meaning, then go to the next step.
The goal of the neural network application step is to classify the actual seismic data with
the model traces obtained from the petrophysical training. The classication is validated
with the seismic/petrophysical discrimination, the explanation rate of the seismic data and
by the coherency with the geological and/or the facies interpretation.
There is an alternative to this workow. The basic workow trains the network on the
synthetic training set, then calibrate the resulting model traces on the seismic data (this
approach will be called "petrophysical training"). The other way is to train the network
on the seismic data, and calibrate the representative traces on the petrophysical models
("seismic training"). Both are worth investigating and we show results for both of them.
Fig. 6.2.1 Workow for the supervised classication with neural networks.
120 Chapitre 6. Petrophysical inversion by neural supervised classication
Neural networks are constituted of several layers of interconnected neurons in input and
data output. Following a learning phase, neurons are linked to the data by minimizing the
error between simulated signal and known signal. Optimization by neural networks requires
no a priori information on input/output relationships, but the quality of optimization hea-
vily depends on the learning phase. Indeed, if the signal is complex, the learning phase will
have to be important to stabilize the neural connections to the output data.
Neural Networks are used since the early 1980s in geophysics and since the 1990s in
reservoir characterization (Poulton 2002). We want to perform a trace shape classica-
tion using a robust pattern-recognition process such as the Kohonen Self Organizing Map
(KSOM) method (Kohonen 1990) ; it is well suited to seismic data due to its ability to
lter out noise and identify representative traces in the reservoir intervals.
The KSOM is a type of Articial Neural Network that is trained using unsupervised lear-
ning. It describes a mapping from a higher dimensional input space to a lower dimensional
map space (one dimensional in this case). KSOMs are used to study the repartition of data
in a N-dimensional space by simply inspecting the data for regularities and characteristics.
This kind of networks organizes itself in such a way as to create an ordered description of
the data.
Where α(τ ) is a monotonically decreasing learning coecient and Xi (t) is the input vector.
The neighbourhood function φ(i, t) depends on the distance between the winning neuron
and neuron i. Regardless of its original value, the neighbourhood function shrinks with
time. At the beginning when the neighbourhood is broad, the self-organizing takes place
on the global scale. This will move each unit in the neighbourhood closer to the input
pattern.
As time progresses the learning rate and the neighbourhood size are reduced. When the
number of iteration is reached, each neuron is independent. If the parameters are well
6.3. Neural Networks : Kohonen Self-Organizing Maps 121
chosen, the nal network should capture the natural clusters in the input data. Figure
6.3.1 shows an illustration of this process.
Fig. 6.3.1 Example of learning phase for an Articial Neural Network. An input vector
is given to the network. The winning neuron (in red) is determined, and then updated so
as to be closer to the input data.
Fig. 6.3.2 Time varying input patterns such as seismic traces can not be evaluated by
Euclidian distance. Here the repetition of the same trace with time shift will be identied
as dierent signals if the Euclidian distance is used. If we use the correlation between the
trace and the neuron, this series will be seen as the same sequence.
At the end of this learning phase, we obtain the model traces that best classify the
input space. Each model trace is assigned a color and a number. Then every trace within
the actual seismic interval is correlated to all the model traces and is classied to belong to
the neuron with the highest correlation. The trace is assigned the same number and color
of the model trace and the correlation. At the end, a repartition map and a tness map
are created.
A neural network calibration is based on the selection of random actual traces at the
beginning of the process : two successive processes of KSOM will not obtain the same
output and/or order of output. The nal maps (classication and correlation) will change
from one process to the next.
The number of neurons constituting the Map The more neurons you will have,
the more precise the map will be. But, depending on the number of neurons dened for
the KSOM, the network can be overtted or undertted. An undertted network can not
represent correctly the signal in a dataset, and this can lead to non identied information.
An overtted network is perfectly adapted to the dataset but will try to explain the noise
in these data.
A trial-error approach should be done to estimate the correct number of neurons. We need
to nd some compromise depending on the problem we want to study. No method exists
to solve this particular question.
6.3. Neural Networks : Kohonen Self-Organizing Maps 123
The size of the neighbourhood The neighbourhood radius gives the size of the active
environment for a neuron. This means that the neurons inside the radius of the winning
neuron will be modied regarding to the learning law. With the number of iteration, this
radius decreases until all neurons are independent.
This parameter is very important since it controls the repartition of the neurons in the
map. For example, if we constrain the learning process with a big neighbourhood radius,
the algorithm will have some trouble to nd specic zones. Inversely, with a small radius,
we will map too many specic zones, which have no real signication.
The learning rate (α) The learning rate is the maximum amplitude of the modications
to be applied on each neuron at each step (or iteration) :
τ
α(τ ) = α0 (1 − ) (6.2)
τ∞
With α0 the learning rate at τ = 0 and τ∞ the number of iteration.
This parameter controls the convergence rate. When the value of α is high, the algorithm
converges quickly toward a representation of the data space, but we can have some insta-
bilities in the process, leading to the loss of the learning. For example, if α = 1, the output
neuron will be the last observation to be presented, regardless of all observations analyzed
before that. Generally, this learning rate is quite weak (α = 0.01 typically). This way, the
modications on one neuron will be gradual during the learning process.
The number of iteration τ The number of iteration τ has an inuence on the neigh-
bourhood radius and the learning rate. At the beginning of the learning phase, if we have
a neighbourhood radius of 2 for example, this means that the winning neuron and the four
neurons around it (two on one side, two on the other side of the map) are updated when
an input trace is presented. Also, the learning rate is strong, so the model traces undergo
an important change. At the end of the learning phase, each neuron is independent (only
the winning neuron is updated) and there is little change on the shape of the trace since
almost all input data has been presented and the model traces are stable/stationary.
124 Chapitre 6. Petrophysical inversion by neural supervised classication
Fig. 6.4.1 A seismic section for the Beta eld, showing the reservoir of interest delimited
by the green horizons.
The geological interpretation from Beta eld provides several stratigraphic markers for
each well, but the number diers between the wells. Basically, those markers follow the
correlation sheet shown in gure 6.4.2.
We need to create a few new markers so as to obtain the same number of macro-layers
for each well (here 18 macro-layers were dened). An example of results for this blocking
step is shown on well W-Beta3 in gure 6.4.4, where the synthetic of the blocked logs is
6.4. Data preparation for the Beta eld 125
Fig. 6.4.2 Geological correlation for the wells from Beta eld. We can see the dierent
geological bodies and how some markers appear or disappear laterally.
Fig.6.4.3 Correlation sheet for the wells from Beta eld, with only some of the markers
shown, for clarity. The property shown here is Vp .
126 Chapitre 6. Petrophysical inversion by neural supervised classication
Fig. 6.4.4 Well log blocking results for well W-Beta3. This gure shows the original logs
in black and the blocked logs in green for (a) the density ρ, (b) the P velocity Vp and (c)
the S velocity Vs . The near seismic trace is compared to (d) the synthetic of the blocked
logs, (e) the synthetic of the original logs. The well-to-seismic tie with the original logs is
of poor quality.
Now that the optimization function is dened, we can initialize the procedure with ini-
tial blocked logs. If the logs were petrophysical properties, we need a petroelastic model
(PEM) to compute the petroelastic logs (P and S velocities and density). We use the wa-
6.4. Data preparation for the Beta eld 127
velet extracted from the seismic data to compute the synthetic seismograms.
We also dene percentages of acceptable variations for the properties and the layer thick-
ness. The initial population is drawn from these variations. The genetic algorithm needs a
large population size, since we have to perturb 4 parameters (ρ, Vp , Vs and the thickness)
for each of the 20 layers, that is, 80 parameters.
An example of the blocking optimization is shown in gure 6.4.5 for Well W-Beta3. The
variations for the petroelastic properties were dened by : ∆ρ = 5%, ∆Vp = 20% and
∆Vs = 15%. A variation for the thickness of the macro-layers was also dened : ∆z = 20%.
The tness function is shown hereafter. This step is repeated for each well of the eld.
Fig. 6.4.5 Well log blocking optimization procedure for well W-Beta3. This gure shows
the original logs in black, the blocked logs in green and the optimized logs in red for (a)
the density ρ, (b) the P velocity Vp and (c) the S velocity Vs . The near seismic trace is
compared to (d) the synthetic of the blocked logs, (e) the synthetic of the best logs and
(f) the synthetic of the original logs. The well-to-seismic tie has been greatly improved.
Fig. 6.4.7 At each location U, for each parameter and each layer, the SGS method draws
values Y from the CCDF provided by the parameter statistics extracted from the wells.
SGS is the most used algorithm to generate conditional realizations. For details on
the SGS algorithm, one can consult Chiles and Delner 1999 and Goovaerts 1997. In the
gaussian domain, we draw a continuous line from the conditional cumulative distribution
function (CCDF) provided by the parameter statistics extracted from the wells, with dif-
ferent values for each location of the simulated dataset (see gure 6.4.7). When we go back
to the empirical domain, we obtain the values for each property in one layer. By repeating
this for each layer of the logs, we obtain a new simulated log for which we compute the
synthetic. Then this new log is considered as a hard data (a constraint) when we proceed
to the next location. Figure 6.4.8 shows the properties and synthetics of a small part of
the Massive Modeling volume. Each 50 traces, an actual well has been placed to constrain
the simulation.
6.4. Data preparation for the Beta eld 129
Fig. 6.4.8 The Massive Modeling. The sections are respectively ρ, Vp , Vs and the synthetic
seismograms. The horizontal scale is the trace number.
The thickness of the seismic interval is an important parameter since it will aect the sta-
bility of the network. This eect will depend upon the geological environment. In a eld
where amplitude variations are very wide on short distance, such as clastic environment,
it is recommended to classify a thicker sequence, so as to obtain more stable model traces.
The tests we performed lead us to an optimal thickness of 100ms.
The gure 6.4.9 shows the neuron repartition map for interval of 50ms and 180ms thick,
centered on the reservoir of interest. We can see that on the 50ms map, the neuron repar-
tition denes very small zones, whereas the 180ms map shows larger zones represented by
one neuron. Both maps were realized with 20 neurons.
For the rest of the Beta eld study, the interval thickness to be considered will be t = 100ms
(the size of the reservoir plus buers so as to avoid the high frequency aspect of the neuron
map).
Fig.6.4.9 Non supervised Kohonen maps realized with an interval thickness of 50ms and
180ms, with 20 neurons.
Number of neurons We need to determine the size of the Kohonen Map, i.e. the num-
ber of neurons that will best represent the data. This number of neurons should avoid
overtting. Figure 6.4.10 shows three maps with dierent sizes and their respective neu-
rons. A stated before, this eld has a high petrophysical variability. We can expect a lot
of trace shapes for the classication. The petrophysical variability is such that we could
obtain more than 100 neurons. As it is, such a classication would be dicult to interpret.
The map with 5 neurons on gure 6.4.10 is clearly unadapted. If we explore the seismic
traces characterized by the rst neuron, for example, we see a lot of traces which are really
dierent from the model trace for this neuron. The map with 50 neurons is also unadapted
since some model traces are pure noise (overtting). The tests we performed lead us to an
optimal number of 35 neurons.
6.4. Data preparation for the Beta eld 131
Fig. 6.4.10 Non supervised Kohonen maps with respectively 5, 20 and 50 neurons. Next
to each map, we show the resulting neurons (for the 50 neurons map, we show only the
second half). Some traces from the 50 neurons map represent noise.
Neighborhood radius The neighbourhood radius controls how many neurons will be
activated and modied along with the winning neuron during the learning phase. This
means that we end up with neurons with more or less progressive changes between them.
A small neighbourhood radius will mean a lot of specic neurons ; a big neighbourhood
radius will mean a lot of continuity between the neurons. The petrophysical variability
of the eld implies a lot of specic trace shape. The geological setting implies grouping
132 Chapitre 6. Petrophysical inversion by neural supervised classication
some of them. For example, with a facies analogy, in a channel, we can nd coarse and
ne deposits, that will have specic trace shape. But we would like to have those two
shapes close together : a "macro" neuron channel could be divided in two "micro" neurons
coarse and ne. This is illustrated on gure 6.4.11 with neighbourhood radius of 1 and 10
neurons. For the map with radius 10, we can see that the light blue and green neurons
are always close. Also, we can see that the tness map for the map with radius 1 has a
slightly more elevated average value. This is because the network is less constrained by the
neighbourhood radius.
Fig. 6.4.11 Non supervised Kohonen maps realized with neighbourhood radius of 1 and
10 neurons. The rst line shows the neuron repartition, the second line the neurons.
6.5. Neural Network classication applied to Beta eld 133
Fig. 6.4.12 Fitness map for the non supervised Kohonen maps realized with neighbou-
rhood radius of 1 and 10 neurons.
Fig. 6.5.1 A small part of the Massive Modeling volume, with the neuron repartition
(top line) and the tness (bottom line) for each trace.
Visualization of the classication We can easily visualize the result of the classica-
tion, by representing each trace with color-coding for the class number and the tness (see
gure 6.5.2). The background color indicates the class number to which the record belongs,
and the saturation of the color indicates the tness. For example, all red traces belong to
class 1. Dark red means a very high tness ; light red means a weak tness. At the top of
the window, the traces are not sorted and are in the normal order of the seismic block. At
the bottom of the window, they are sorted by class number and tness. This display shows
the general organization of the model traces on the dataset.
Fig. 6.5.2 A small part of the color-coded visualization of the classication. On the top-
row, the traces are shown as they are in the dataset. On the bottom-row, the traces are
ordered by class and tness. Since we show about 300 traces, once ordered we see only a
small portion of two classes on the bottom row.
6.5. Neural Network classication applied to Beta eld 135
Classication QC A graph is created using the tness map from the Kohonen network,
representing the percentage of traces with tness greater than or equal to the tness abs-
cissa. This way we can see if the output neurons from the Kohonen classication are well
suited with the traces of the dataset.
Fig. 6.5.3 A graph showing the percentage of model traces exceeding the tness on the
abscissa. For example, 50% of the model traces have a tness greater than 90%, which is
fully satisfactory.
The graph 6.5.3, applied on the Massive Modeling volume, shows that the model traces
explain the petrophysical models.
Classication analysis The analysis of classication provides a link between the trace
characteristics determined by the classication and the input parameters of the seismic
modeling.
For each class, the algorithm retrieves all traces belonging to the class and with tness above
the threshold. It then calculates the average and standard deviation on all these traces for
each input parameter of the model and for each layer. These results are presented in the
table 6.5.4. Such tables enable to follow the behaviour of the petrophysical models in a
statistical sense.
Fig. 6.5.4 A part of the table for the classication analysis. We obtain the parameter
average and standard deviation for each neuron. Here are shown the thickness, ρ and Vp
for the two rst layers and the three rst neurons.
136 Chapitre 6. Petrophysical inversion by neural supervised classication
Trace likelihood To ensure that the model traces reect the petrophysical discrimina-
tion in the seismic response, we calculate the likelihood of each model trace to explain the
synthetic traces at the well location. This is illustrated on gure 6.5.5, where we see that
distinct groups of model traces explain the synthetics of the 5 wells.
Fig. 6.5.5 Likelihood of each model traces with the synthetics of the ve wells for Beta
eld. A likelihood of 1 is the best correlation value between a model trace and a synthetic
trace.
Then we classify with the neural network the actual seismic dataset with the model
traces derived from the descriptive step. We obtain the maps shown on gure 6.5.6, where
one can visualize the geological and structural content of the seismic interval. The map on
the top, showing the repartition of the neurons, displays the shapes of two channels, one
going from Southwest to Northeast, the other one going from West to East. This map also
displays some irregular lineaments, which are in fact faults.
The map on the bottom shows the tness of the winning model traces with the concerned
seismic trace. The average tness is around 85% ; in the channelized zones, the value is a bit
less. This means that the supervised classication was successful in linking the petrophysical
models with the actual seismic response.
6.5. Neural Network classication applied to Beta eld 137
Fig. 6.5.6 Predictive modeling : the model traces from the training phase have been
applied to the actual seismic traces. The maps show the classication and the tness. The
white circle on the tness map gives the location of well W-Beta1.
138 Chapitre 6. Petrophysical inversion by neural supervised classication
Fig. 6.5.7 Likelihood of each model traces with the actual seismic trace of the ve wells
for Beta eld.
We check the discrimination by calculating the likelihood of each model trace to explain
the actual seismic traces at the well location (see gure 6.5.7).
The likelihood graph shows that a distinct groups of model traces explain the seismic
traces on 4 wells (likelihood ≥ 0.8) while other can be linked to the rst well (0.6 ≤ likeli-
hood ≥ 0.7).
This means that there was no model quite comparable to the petrophysical models from
the well W-Beta1 in the seismic data. This is conrmed by the tness on the map 6.5.6 at
well W-Beta1, which is the lowest tness on the map.
We perform the training phase on the actual seismic data (see gure 6.5.8). The model
traces repartition shows several geological shapes : two channels (dened by green colors),
debris ow deposits (in yellow, orange and red) and shale deposits (dark blue and purple).
We also see some long lineaments which cut across the model trace repartition : the faults.
The tness map has a very high average value (greater than 90%), except for the channelized
zones, where it is down to 70-80%.
6.5. Neural Network classication applied to Beta eld 139
Fig. 6.5.8 Descrictive modeling : the training phase has been performed on the actual
seismic traces. The maps show the classication and the tness.
140 Chapitre 6. Petrophysical inversion by neural supervised classication
We calculate the model trace likelihood with the seismic traces at well location. At this
stage, since we do not know the petrophysical models which characterize the traces, this
graph (6.5.9) shows the model trace variability against the actual seismic.
Fig. 6.5.9 Likelihood of each model traces with the actual seismic trace at the ve wells
for Beta eld.
Fig. 6.5.10 Graph showing the percentage of model traces exceeding the tness on the
abscissa. For example, 60% of the model traces have a tness greater than 50%.
6.5. Neural Network classication applied to Beta eld 141
In a rst QC, we can check the explanation rate of the petrophysical model (see gure
6.5.10). 60% of the model traces explain 50% of the Massive Modeling dataset ; this is not
fully satisfactory.
Fig. 6.5.11 A small part of the color-coded visualization of the classication. The colors
displayed on the top row are quite dull, this means that the tness for the traces presented
here is low. See gure 6.5.2
Fig. 6.5.12 Likelihood of each model traces with the synthetic traces of the ve wells for
Beta eld.
Fig. 6.5.13 The classication map from the petrophysical training with geological inter-
pretation on top, with structural interpretation on the bottom picture.
6.5. Neural Network classication applied to Beta eld 143
Fig. 6.5.14 The classication map from the seismic training with geological interpretation
on the top, with structural interpretation on the bottom.
144 Chapitre 6. Petrophysical inversion by neural supervised classication
On the map 6.5.13, the channels interpreted on the data are very well dened by the
model traces (bluish and purple colors). The debris ow deposits inside the channels and
in the crescent-shaped area (top-middle in the gure) are fairly well dened by the red
to yellow model traces. A large area of the map is green, which correspond to the shale
deposits. The structural interpretation shows that some of the faults can be seen in the
classication map (such as the top-left and bottom-left area of the map).
On the map 6.5.14, the channels are fairly well dened by the model traces (green and
light blue colors). The debris ow deposits inside the channels and in the crescent-shaped
area are very well dened by the red and orange model traces. The shale deposits are
characterized by several model traces (from dark blue to purple colors) with very wide
coherent zones. The structural interpretation of this map shows that most of the faults can
be seen in the classication map.
The classication map, from the seismic training, actually shows some variability (seve-
ral distinct model traces) for the shale deposits, whereas the classication map from the
petrophysical training explains the shales with very similar model traces (green colors).
This was to be expected, since the ve wells from the eld were drilled in dierent channel
sequences. The petrophysical models for shales could not be derived from the geostatistical
variability given by the wells.
6.6. Data Preparation for the Gamma Field 145
The data being condential, they are not fully described in this chapter.
The interval of interest is 50ms thick. The geology shows a layer-cake model with three
main facies : dolomites, limestones and anhydrites with low petrophysical variability. The
objective of this study is to see if seismic facies can be discriminated by the methodology.
We create 19 new markers (i.e. 20 macro-layers) for each well. By creating these new
markers, we increase the probability of introducing variability for the Massive Modeling
step, since the well segmentation approach relies on log shape.
An example of results for this blocking step is shown on one of the wells in gure 6.6.1,
where the synthetic of the blocked logs is compared with the actual near seismic.
146 Chapitre 6. Petrophysical inversion by neural supervised classication
Fig. 6.6.1 Well log blocking results. This gure shows the original logs in black and the
blocked logs in green for (a) the density ρ, (b) the P velocity Vp and (c) the S velocity Vs .
The near seismic trace is compared to (d) the synthetic of the blocked logs, (e) the synthetic
of the original logs. The well-to-seismic tie with the original logs is of poor quality.
Fig. 6.6.2 Well log blocking optimization procedure. This gure shows the original logs
in black, the blocked logs in green and the optimized logs in red for (a) the density ρ, (b)
the P velocity Vp and (c) the S velocity Vs . The near seismic trace is compared to (d) the
synthetic of the blocked logs, (e) the synthetic of the best logs and (f) the synthetic of the
original logs. The well-to-seismic tie has been greatly improved.
6.6. Data Preparation for the Gamma Field 147
Fig. 6.6.3 The Massive Modeling. The sections are respectively ρ, Vp , Vs and the synthetic
seismograms.
148 Chapitre 6. Petrophysical inversion by neural supervised classication
Thickness inuence The gure 6.6.4 shows the neuron repartition map for interval
of 50ms and 100ms thick, centered on the reservoir of interest. We can see that the two
maps show similar features, even if the zones on the 50ms map are smaller. There are no
signicant high frequency variations on the 50ms map, as shown in the rst case study.
Because of this, the two maps can be interpreted. Both maps were realized with 30 neurons.
For the rest of the study, the interval thickness to be considered will be t = 70ms (the exact
size of the reservoir on the seismic data plus buers).
Fig.6.6.4 Non supervised Kohonen maps realized with an interval thickness of 50ms and
100ms.
Number of neurons Figure 6.6.5 shows three maps with dierent number of neurons
and their respective model traces. This eld shows low petrophysical variability on the well
data. We can expect a small number of specic shapes for the classication. The map with
5 neurons on gure 6.6.5 is clearly unadapted. If we explore the seismic traces characterized
by the rst neuron, for example, we see a lot of traces which are signicantly dierent from
the model trace for this neuron. The map with 50 neurons is overtted since some model
traces are pure noise. We should choose an intermediate number of neurons, as shown with
the 30 neuron map, where there are no noisy traces. 30 neuron maps seem to be the best
size.
6.6. Data Preparation for the Gamma Field 149
Fig. 6.6.5 Non supervised Kohonen maps with respectively 5, 30 and 50 neurons. Next
to each map, we show the resulting neurons (for the 50 neurons map, we show only the
second half). Some traces from the 50 neurons map represent noise.
Neighbourhood radius The petrophysical variability of the eld implies a small number
of specic trace shapes. But we would like to have as many shapes as possible to be able to
discriminate between facies. We should have a very small neighbourhood radius, this way,
the network can nd dierent shapes. This is illustrated on gure 6.6.6 with neighbourhood
radius of 1 and 10 neurons. We can see that between the two tests, the neurons (line 2) are
very similar. This means that even a radius of 1 is too much for the network. To optimize
150 Chapitre 6. Petrophysical inversion by neural supervised classication
the search for original patterns, we will use a neighbourhood radius of 0, and only the
winning neuron will be updated in the learning phase.
Fig. 6.6.6 Non supervised Kohonen maps realized with neighbourhood radius of 1 and
10 neurons. The rst line shows the neuron repartition, the second line the neurons. The
last line shows the tness.
6.7. Neural network classication applied to Gamma eld 151
Fig. 6.7.1 A small part of the Massive Modeling volume, with the neuron repartition
(top line) and the tness (bottom line) for each trace.
152 Chapitre 6. Petrophysical inversion by neural supervised classication
Fig. 6.7.2 A small part of the color-coded visualization of the classication. The colors
displayed on the top row are very bright, this means that the tness for the traces presented
here is high. See gure 6.5.2 for explanation.
In a rst QC, we can check the explanation rate of the petrophysical model (see gure
6.7.3). 70% of the model traces explain 80% of the Massive Modeling dataset, which is a
very good result. The parameter variations for each model traces can be analyzed (see table
6.7.4). For the P velocity for example, the statistic show a very small scattering around the
average value (the standard deviation is around 2% of the mean), for all the layers and all
model traces. This means that the network describes very specic petrophysical models in
the Massive Modeling volume.
Fig. 6.7.3 Graph showing the percentage of model traces exceeding the tness on the
abscissa. 70% of the model traces have a tness equal or greater than 80%, which is fully
satisfactory.
6.7. Neural network classication applied to Gamma eld 153
Fig. 6.7.4 Part of the classication analysis table, showing the parameter average and
standard deviation for each neuron. Here are shown the thickness, ρ and Vp for the two
rst layers and the four rst neurons.
Fig. 6.7.5 Likelihood of each model traces with the synthetic traces of the twelve wells
for Gamma eld.
The likelihood graph 6.7.5 shows the groups of model traces that explain the synthetic
traces of the wells. The explanation rate is fully satisfactory so we can go ahead with the
study.
The map on the bottom shows the tness of the winning model traces with the concerned
seismic traces. The average tness is around 80% and in the three zones described just
before, the tness is around 60%.
Fig. 6.7.6 Predictive modeling : the model traces from the training phase have been
applied to the actual seismic traces. The maps show the classication and the tness. The
location of the two well clusters has been removed for condentiality.
6.7. Neural network classication applied to Gamma eld 155
Fig. 6.7.7 Likelihood of each model traces with the actual seismic trace at the well
location for Gamma Field.
We check the discrimination by calculating the likelihood of each model traces with the
seismic traces at the well locations (see gure 6.7.7). This graph shows that a few distinct
groups of model traces explain the seismic traces, while the other have a low likelihood on
several neighbouring wells. The overall likelihood is quite low.
All these observations mean that the supervised classication was not quite successful
in linking the petrophysical models with the actual seismic traces.
We perform the training phase on the actual seismic data (see gure 6.7.8). The model
trace repartition shows essentially three units (orange to green), surrounded by blue and
purple zones. The largest zone (bottom-left of the gure) shows 2 purple lineaments which
are in fact collapse zones (permeable drains). The tness map has a very high average
value (greater than 90%), except for some zones which cannot be linked to any geological
interpretation. These zones have a low quality signal which might be explained by gas-
bearing deposits.
156 Chapitre 6. Petrophysical inversion by neural supervised classication
Fig. 6.7.8 Descriptive modeling : the training phase has been performed on the actual
seismic traces. The maps show the classication and the tness.
6.7. Neural network classication applied to Gamma eld 157
We calculate the model trace likelihood with the seismic traces at the well locations.
This graph (6.7.9) shows the model trace variability against the actual seismic. We see two
kinds of behavior on this graph : the model traces which have a high likelihood on several
neighbouring wells, and the model traces that have an overall low likelihood except for one
well.
Fig. 6.7.9 Likelihood of each model traces with the synthetic traces of the well location
for Gamma Field.
158 Chapitre 6. Petrophysical inversion by neural supervised classication
Then we classify with the neural network the Massive Modeling dataset with the model
traces derived from the descriptive step, to see if they nd a comparable seismic response
within the generated petrophysical models.
Fig. 6.7.10 Graph showing the percentage of model traces exceeding the tness on the
abscissa. For example, 70% of the model traces have a tness greater than 70%.
Fig. 6.7.11 A small part of the color-coded visualization of the classication. The colors
of some groups of traces displayed on the top row are quite dull, this means that the tness
for these traces is low. See gure 6.5.2 for explanation.
The explanation rate of the Petrophysical models, provided by the gure 6.7.10 is very
high : 70% of the model traces have a tness greater than 70%, which is fully satisfactory.
The gure 6.7.11 shows that in spite of the high explanation rate, some area have a very
low tness.
6.7. Neural network classication applied to Gamma eld 159
Fig. 6.7.12 Likelihood of each model traces with the synthetic traces of the twelve wells
for Gamma eld.
The nal step is to check the likelihood of the seismic model traces to explain the pe-
trophysical models from the Massive Modeling dataset (see gure 6.7.12).
We observe that the maximum likelihood value for all model traces is around 1. We observe
two kinds of behaviour on this gure : the model traces either have a very high likelihood
with the well synthetics, or they have a average tness for all the wells. For the latter,
this means that the seismic data contains more variability than the petrophysical models
from the Massive Modeling volume, discriminated by the predictive modeling. There are 2
clusters of wells in the Gamma eld. They actually represent two families of very similar
petrophysical properties. The seismic data contains variabilities which were not sampled
by the wells.
160 Chapitre 6. Petrophysical inversion by neural supervised classication
Fig. 6.7.13 The classication map from the petrophysical training with geological inter-
pretation on top, with structural interpretation on the bottom picture.
6.7. Neural network classication applied to Gamma eld 161
Fig. 6.7.14 The classication map from the seismic training with geological interpretation
on the top, with structural interpretation on the bottom.
162 Chapitre 6. Petrophysical inversion by neural supervised classication
On the map 6.7.14, three zones are identied by the red to green model traces. We think
that these zones correspond to dominant facies signatures which are interpreted as the
dolomitic (left side), anhydritic (top-right area) and limestone (bottom-right area) domi-
nant facies. The structural interpretation of this map shows that the faults can be seen in
the classication map, and the main collapse zones (permeable drains) are dened by the
purple model traces.
The classication map, from the seismic training, actually shows more variability (several
distinct model traces) on the whole map, whereas the classication map from the petro-
physical training shows only some of the model traces. This was to be expected, since only
specic petrophysical models have been generated with the two well clusters. The missing
petrophysical models could not be derived from the geostatistical variability given by the
wells.
6.8. Concluding remarks and perspectives 163
Rigorous study of this technique on two 3D datasets with dierent geological settings
has produced promising results to derive properties that could be used in a reservoir mo-
deling sense. On the clastic case, the geological bodies were identied by the methodology.
A seismic facies interpretation was also possible on these data. On the carbonate case, a
macro-facies interpretation has been done, and the collapse zones identied.
The choice of parameters for our methodology is highly data-dependent, but tools are
available to quantify their eect, and guide the user through the study.
We are currently implementing the last step of our methodology, where we get the pe-
trophysical models assigned to each seismic trace.
164 Chapitre 6. Petrophysical inversion by neural supervised classication
Partie III
Conclusions et perspectives
165
Concluding remarks and future research
The work of this thesis have resulted in algorithms and tools programmed in Java and
implemented in SISMAGE
r
, the seismic interpretation platform developed by TOTAL.
167
168 Concluding remarks and future research
taking into account dierent scales of observation ; these approaches form the essential
work of this thesis.
Future work
Our work was developed in pursuit of operational eciency. This led to inevitable
simplications. For example, the method for seismic modeling is deliberately simplistic
(multi-1D convolution) to be as interactive as possible.
The methods developed in this thesis have been tested and validated on three actual case
studies chosen for their geological diversity. It will be interesting to test these methods
on other cases to ensure their robustness. The approaches are complemented with tools
to discuss the validity of results. All the tools are operational and will remain in constant
evolution.
Possible improvements
The gradual deformation inversion remains a product which could be improved :
Optimization of property space exploration : we have used the basic formulation of
the gradual deformation. It may be replaced by a more sophisticated approach that
uses several secondary realizations and a gradient search (Hu and Ravalec-Dupin
2004).
Optimization based on actual seismic data : the processes were constrained by syn-
thetic data, calculated from previous inversion results providing "solutions" to com-
pare our property realizations. The next step will be to directly use the actual data.
169
Using a petroelastic model : the current processes recover the density and seismic
velocities. The implementation may include a PEM to recover the petrophysical
properties.
171
Annexe 1 : SEG Annual Conference (San
Antonio 2007
Structural uncertainty effects on Reservoir grid infilling
Audrey NEAU* 1,2, PhD. Student, [email protected], Pierre THORE* 1, Béatrice de VOOGD* 2
1: TOTAL Pau, FRANCE
2: Modeling and Imaging in Geosciences, University of Pau, UMR5212 UPPA-CNRS-TOTAL FRANCE
attribute for reservoir characterization is impedance Panel 3 grid is G3. G2 top and bottom horizons are globally
(product of density and velocity), either acoustic or elastic, below initial horizons, whereas G3 top and bottom horizons
which contains relationships between rock physics and are globally above initial horizons.
seismic (Mukerji et al, 2001). From the multiple
realizations of GeoSI, we use P and S impedances to A hundred P (IP) and S (IS) impedance cubes are generated
calibrate attributes. for several structural models, using the same layering and
the same prior model for each of them. The multiple
Relationships, between petroelastical (impedances) and realizations obtained permit to evaluate reservoir
petrophysical (facies) reservoir properties are analyzed at uncertainty. The standard deviation ranges from 5% to 10%
wells. They are reviewed at seismic scale by extracting around the mean of all realizations. Figure 4 shows the P
minicubes around well locations. mean impedance on the 26th layer for each model.
The cross-plot Ip-Is colored by facies (or other property) at Considering the mean impedance for each structural model
the well is divided into different parts. These parts must be leads to differences mainly due to structural uncertainties.
as homogeneous as possible in the domain delimited by the After validation of the inversion at the well, we are able to
attributes. The different zones revealed by the cross-plot go ahead and classify the inversion results.
are projected on the seismic cubes to generate geological
property cubes. The impedance attribute cubes are then used to estimate the
most likely facies by cross-plot analysis. Figure 5 shows
Example from a West African offshore reservoir the facies classification from the inverted traces compared
to the known facies at the well, for each structural model.
The reservoir studied is a complex of submarine channel For each cross-plot, facies are analyzed separately.
reservoirs. This is a thick Oligocene sequence of silty clays
that includes important sandstone reservoirs in deepwater Except for facies 1 (too much values) and facies 5 (not
areas. The deposits are mainly turbiditic, with interbedded enough values), we observe that G2 generally has lower IS
channel sands. values (and to a lesser extent, IP too), while G3 displays
globally higher IS values. This shows that probability cubes
Near, mid, and far offset stacks are available at a resolution drawn from these classifications are affected by structural
of 3 ms over a 4.8 Km x 1.8 Km area. Only one well in the uncertainties.
higher part of the structure is available, with raw and
interpreted logs. Five major lithofacies are defined from Conclusions
these logs and are used for facies classification (shales (1);
laminated sands (2); fine sands (3); coarse sands (4); These results show that structural uncertainty effects need
cimented sands (5)). to be assessed on property probability cubes before
performing reservoir grid infilling. This has been done by
Reservoir top and bottom surfaces picked from the seismic considering not one but a class of structural models, all
data are used to build a corner-point stratigraphic grid of consistent with the seismic and well data.
114 x 159 x 113 = 2 048 238 nodes. Panel 1 in Figure 3
shows the initial reservoir grid with well path (red track). The approach described here is a practical and functional
The horizontal resolution of the grid ranges from 45 to 70 mean to insert the concept of structural uncertainty in the
m. The vertical resolution ranges between 3 and 15 m. reservoir modeling workflow.
Top and bottom layers are extracted from the initial Acknowledgements
reservoir grid as reference horizons. An uncertainty map
(based on seismic picking derived uncertainties) is The authors would like to thank TOTAL for the release of
associated to each of them. this work for publication. The authors would also like to
acknowledge G.R.C. and C.G.G. collaborators for the use
With the reference structural model and the uncertainty of GeoSI software.
maps, one hundred equiprobable structural models are
generated; two of them are illustrated in Figure 3. The
initial grid (Panel 1, Figure 3) is G1, Panel 2 grid is G2 and
Figure 3: Structural models realized with ALEA: Panel 1: Initial grid; Panels 2 and 3: Deformed grid with uncertainty map.
Figure 4: Layer 26 of P mean impedance on 100 realizations from GeoSI respectively for G1, G2 and G3 grids.
Figure 5: Supervised facies classification, respectively for G1, G2 and G3 grids. Facies correspond to their
legend in the text.
REFERENCES
Fournier, F., P. Y. Dequirez, C. G. Macrides, and M. Rademakers, 2002, Quantitative lithostratigraphic interpretation of seismic
data for characterization of the Unayzah Formation in central Saudi Arabia: Geophysics, 67, 1372–1381
Mukerji, T., P. Avseth, I. Takahashi, and E. F. Gonzalez, 2001, Statistical rock physics: Combining rock physics, information
theory, and geostatistics to reduce uncertainty in seismic reservoir characterization: The Leading Edge, 3, 313–319.
Thore, P., A. Shtuka, M. Lecour, T. Ait-Ettajer, and R. Cognot, 2002, Structural uncertainties: Determination, management, and
applications: Geophysics, 67, 840–852.
Williamson, P. R., A. J. Cherrett, and R. Bornard, 2007, Geostatistical stochastic elastic inversion - An efficient method for
integrating seismic and well data: Presented at the 69th Conference and Exhibition, EAGE.
SUMMARY
Blocking well logs is an initial and mandatory step for many processes. It is an upscaling procedure and
has no unique solution. In order to find a set of "optimal blocking" for geophysical applications we
propose to use a criterion based on seismic difference between synthetic seismogram computed on the fine
scale logs and synthetic seismogram computed on the blocked logs. The optimisation process is performed
with a genetic algorithm because genetic algorithms have the ability to solve non linear problem and to
provide a set of solutions which gives an idea of the uncertainty in the model space. Thanks to the limited
number of parameters, blocked well logs can be used lately in many applications such as Massive
Modelling.
Figure 1: Modelling scheme. If petrophysic logs are used the first step is to use a petroelastic
model to transform them into petroelastic logs, this first step is skipped otherwise.
Initialization: The process starts with an initial blocking based on the limits of geological
layers. For each log the blocked value in each layer has been computed using a standard
procedure (usually averaging). The logs can be neutron-density, shale-ratio, porosity… and a
petroelastic model that is used to transform these logs into petroelastic log i.e. compressional
velocity (Vp) shear velocity (Vs) and density (Ro) or directly the petroelastic logs themselves
(see figure 1). Associated with these initial values, we define a range of acceptable variations
for both the parameters and the thicknesses of the layers. The initial population used in the
GA is drawn randomly within these ranges of variation.
Parameterization and Evolution: All the parameters (ΔZi, Roi, Vpi, Vsi, i ε [1,N], where N is
the number of layers) are coded into a single chromosome .Typically the blocking is limited
to the reservoir i.e. a thickness of 50 to 300 meters and the number of layers usually does not
exceed 30. Even in these conditions the number of parameters is over 100 and requires a large
population size in order to keep enough diversity i.e. to avoid the problem of generic drift
(Mallick, 1995). In order to keep the population size reasonable (around 1000) and therefore
to diminish the CPU time, we use GA in a layer stripping way. Starting from the top the upper
layers (e.g. 10) are optimized during the first iterations (e.g. 50) then the next layers (with an
overlapping on the previous layers) during the next (50) iterations and so on… Figure 3
shows the evolution of the cost function for two separate runs. At the beginning, the upper
part of the log is optimized and the cost function rapidly diminishes then stabilizes and starts
to diminish again when the next part of the log is optimized.
Examples
Figure 2 shows an example of optimization on a segment of a well of 200 meters and with 30
layers in the blocked logs. The final match between the synthetics is not perfect particularly at
the top and the medium part of the synthetics. The residual discrepancy can be attributed to a)
the limited number of layers in the blocked well b) the constraint of the initial blocked logs c)
the small range of variation allowed for each parameter.
Figure 3 presents the cost function curves (inverse of fitness curve) for 2 runs of the GA.
During the first 50 iterations only the upper part of the log is optimized, then in the next 50
iterations, the medium part and finally the bottom part of log are optimized. The fitness is
changing very rapidly at the beginning of each step and then stabilizes when the optimal is
reached. This can also be seen on figure 4 which displays the best synthetic trace at each
iteration. The main variability is concentrated from iteration 1 to 15 and then from 51 to 75.
Note that the upper part of the traces remains unchanged after iteration 50 since the upper of
the log is only optimized during the 50 first iterations.
Figure 3: Evolution of the cost function of the best individual at each iteration (two runs
of the GA optimisation). The population size is 2000 and the number of layer is 30. In
the 50 first iterations, the upper part of the model is optimized, then the middle part in
the next 50 and the remaining part in the last 50.
Figure 4: Best synthetic trace at each iteration. In the 50 first iterations the upper part
of the model is optimized, then the middle part and finally the remaining part. The last
5 traces correspond to the synthetic trace computed with the original logs.
Example of Application:
Massive Seismic Modelling (Julien et al 2002) consists in computing the seismic response of
a set of pseudo wells. It is used to compute the tuning curves for pre or post stack attributes
Figure 5: Interpolation of blocked well logs. The synthetics traces are computed on
pseudo well logs generated by geostatistical interpolation of optimally blocked well
logs.
Conclusion:
A technique for optimizing blocking of well logs at the geological scale has been presented.
The criterion for optimality is the comparison of the synthetic seismograms generated with
original logs and the blocked logs. We use genetic algorithm for optimization. The blocked
logs after optimization can be used for any application which relates geology to seismic.
Bibliography:
Chouinard P.N., Paulson K.V., 1998, A Markov-Gauss algorithm for blocking well
logs: Geophysics 53, 1118-1121.
Julien, P., F. Pivot, A. Douillard, Y. El Ouair, and S. Toinet, 2002, The importance of
prestack massive seismic modeling for AVO calibration and seismic reservoir
characterization: 72nd Annual International Meeting, SEG, Expanded Abstracts , 1731-1734.
Mallick S., 1995, Model-base inversion of amplitude-variation-with-offset data using genetic
algorithm: Geophysics 60, 939-954.
Painter S., Beresford G., Paterson L., 1995, On the distribution of seismic reflection
coefficients and seismic amplitudes: Geophysics 60, 1187-1194.
Prüßmann J., 1994, Averaging sonic and impedance logs at variable intervals: 64th Annual
International Meeting, SEG, Expanded Abstracts , 81-85, Session: BG3.5.
Walden A.T., Hosken W.J., 1988, Choosing the averaging interval when calculating primary
reflection coefficients from well logs: Geophysical Prospecting 36, 799-824.
Combining the Gradual Deformation Method with Seismic Forward Modeling to constrain Reservoir
Models
Audrey NEAU* (1, 2), PhD. Student, Pierre THORE (1), Béatrice de VOOGD (2).
1: TOTAL Pau, France
2: Modeling and Imaging in Geosciences, UMR5212 UPPA-CNRS-TOTAL, University of Pau, France
INTRODUCTION
To achieve the best possible history match of production data, the The basic GDM formulation
reservoir model must synthesize as much information as possible. We Several versions of GDM exist, but they are all based on the same core
decrease the number of possible models by collecting all available ge- scheme. GDM requires Gaussian realizations and aims at generating
ological, geophysical as well as dynamic reservoir information, and gradual (meaning continuous) perturbations on an initial property re-
integrate them in the reservoir description. alization such as to better match the data.
Constraining reservoir properties simulations to seismic data leads to The cornerstone for GDM is that the sum of two Gaussian random
much more accurate reservoir models. Due to difference in scale, the functions is a Gaussian random function (Hu, 2000). The basic rela-
information derived from seismic data must be upscaled into the reser- tionship is of the form:
voir, leading to much inconsistencies and information losses.
Y (t) = Y0 .cos(t) +Y1 .sin(t). (1)
We propose a new methodology to directly recover seismic attributes
at reservoir scale, by combining the gradual deformation parameteri- Y0 and Y1 are two independent Gaussian random functions with iden-
zation for the geostatistical reservoir model, and an inversion process tical mean, variance and variogram. Y0 is considered as the initial
to match seismic data. The proposed methodology includes: realization and Y1 as a given complementary realization.
Whatever the gradual deformation coefficient t, Y (t) is a Gaussian ran-
1. Combination of the current and new random realizations un- dom function with the same two-order statistics than Y0 and Y1 . By
der GDM operating principles; changing t (from t = 0 to t = π2 ), we create a continuous chain of
2. Conditioning with the static hard data (well logs); realizations, connecting Y0 and Y1 .
3. Synthetic seismic generation with the new realization; Post-conditioning:
4. Minimization of the cost function on the GD coefficient t; The realizations, obtained through linear combination, are unlikely to
5. The optimized realization is set as the new current one; honor the static hard data. A post-processing step consists in using
6. Go back to step 1 until the objective function is satisfactory. Kriging techniques to modify the realizations and make them consis-
tent with these data. Conditioning the realization Y to the localized
Combining the Gradual Deformation Method with Seismic Forward Modeling to constrain Reservoir Models
observations is done through the relationship: To generate a seismic dataset, the reservoir grid is regularly resam-
pled in the XY direction at the seismic scale. By this means, a set of
Yc = YdK + [Y −YK ]. (2) pseudo-wells are computed, and their intersections with the reservoir
grid allows to recover the elastic parameters for each layers. We call
Y is the unconditional realization. YdK and YK are respectively the krig- this ensemble of pseudo-wells a voxet.
ing estimates of the known data and the Y values at the same location.
Yc is the conditional realization. Usually, a reservoir grid (and likewise the voxet) is defined in depth
This formalism ensures that the conditional realization honors the lo- domain, while the seismic is in time domain. A time reference (a
calized observations, while preserving the spatial variability. velocity model) is needed to convert the voxet from depth to time do-
main.
Anamorphosis: Then, the petroelastical properties are converted into reflection coeffi-
This method also extends to non-Gaussian realizations, provided that cients that are convolved with the seismic wavelet. Any method can be
we transfer the initial experimental properties into the Gaussian do- used to obtain the reflection coefficients, but here we choose a simple
main. To do so, an experimental variable x can be transformed into a 1D convolution algorithm. The resulting traces give a good approxi-
Gaussian variable V by Gaussian anamorphosis: mation of a zero-offset seismic.
V = φ −1 (F(x)), (3) A very important consideration on this whole process is the issue of
upscaling. Indeed, the construction of the regular grid usually involves
where φ denotes the standard normal cumulative distribution function vertical upscaling and lateral downscaling, depending on the size of
and F the sample cumulative distribution function. After combina- the reservoir grid cells. This transformation (reservoir grid into voxet)
tion, we back-transform the optimized realization in the experimental is non unique and must be done carefully to minimize the scale change
domain. effects (Thore, 2006).
The layer thickness, particularly, has an important impact on the verti-
cal sampling of elastic properties, and thus on the quality of generated
GEOSTATISTICAL SIMULATION ALGORITHM: synthetics. One solution is to ensure that the reservoir grid has a lay-
THE FFT MOVING AVERAGE METHOD ering thin enough to overcome this issue.
Otherwise, we need to interpolate the vertical samples to avoid arti-
Geostatistical simulation can provide us with independent realizations facts due to the reservoir grid discretization.
of reservoir properties, all consistent with the available informations
(well log data, spatial variograms). There are different existing algo-
rithms to create those new realizations. STOCHASTIC OPTIMIZATION APPROACH
The most efficient for the Gradual Deformation scheme is the Fast
Fourier Transform Moving Average method or FFTMA (Le Ravalec Objective function:
et al., 2000), which performs Gaussian random realizations by gener- Constraining stochastic reservoir models to actual data is an optimiza-
ating random numbers with the needed spatial variability. The FFTMA tion problem involving the definition of a minimization process be-
method allows the generation of a Gaussian realization y with: tween the actual data and the responses generated from the property
realizations. Any objective function can be defined in the process.
y = μ + g ⊗ z, (4) Here, we quantify this misfit by defining an objective function J with
the squared-difference between real (dobs ) and simulated (dsim ) data:
where μ is the mean, z a white Gaussian noise and g is a complex
analytic function, coming from the covariance matrix C = g ⊗ g̃, de- 1 Nobs i
pending on the spatial variability of the realization. To deal with the J= ∑ (dsim − dobs
2 i=1
i
)2 , (5)
use of the convolution product to determine C, then y, we need to
solve equation 4 into the spectral domain. With the Fourier Transform
with Nobs the number of observations.
(FT), we transform the convolution product into a simple multiplica-
A reservoir model is a high dimensional object. The number of cells
tion: FT (g ⊗ z) = G ∗ Z. If S = FT(C), then we can create y with:
amounts to N = 104 − 107 . It is, therefore, very irrelevant to optimize
the objective function with respect to the reservoir size.
• Build C with the variogram from the initial realization; The GD technique reduces the optimization to the single parameter t.
• Generate the white noise z; In this way, we reduce a N-dimensional problem to a one-dimensional
• Calculate the Fourier Transform of C and z to obtain S and Z, one. Starting with Y0 the initial realization and Y1 a new realization,
respectively; we create a chain of realizations Y (t) with equation 1 and we choose
• Derive G from S; the optimized realization which objective function is the lesser.
• Make the product G ∗ Z;
• Calculate the Inverse Fourier Transform of G ∗ Z to obtain Iterative approach:
g ⊗ z; The optimization scheme shown above might not reduce the objective
• Derive y from equation 4. function to a low-enough level in a one-step process. Consequently,
we must repeat the procedure by choosing Y (t) as the new initial re-
alization and generating a new complementary one, until the objective
Full mathematical details for the FFTMA generator can be seen in Le function is satisfactory. At iteration n, equation 1 can be rewritten:
Ravalec et al. (2000).
Yn (t) = Yn−1 .cos(t) + un .sin(t). (6)
SYNTHETIC SEISMIC GENERATION where Yn−1 is the optimized realization at iteration n − 1, and un are
the complementary realizations at each iteration. Yn is the new opti-
Seismic modeling is an essential part of seismic inversion algorithms. mized realization which improves (or at least maintains) the objective
To evaluate the quality of GD inversion, we need to generate the syn- function from iteration n − 1.
thetic seismic from the optimized realizations.
A reservoir grid is neither regularly sampled in space, nor in depth.
Combining the Gradual Deformation Method with Seismic Forward Modeling to constrain Reservoir Models
Synthetic model description: The Gradual Deformation Method is a parameterization technique that
We have a 3D reservoir model consisting of 26*26*113 = 76 388 cells, allows to continuously perturb reservoir property realizations to solve
defined in the time domain (no time to depth conversion is required), an inverse problem. It takes into account scale change in the inversion
with the following properties: process by simulating new realizations in a reservoir grid, while min-
imizing an objective function between the actual seismic data and the
• ”True” P-velocity Vpreal and density ρreal , from which we synthetic seismic datasets obtained from these realizations.
create the ”actual” seismic dataset; It has been tested on a 3D synthetic reservoir model and shows im-
• Initial P-density Vpinit and density ρinit as input realizations provements of the reservoir realizations. The algorithm is not fully
for the process. Here, they were created from a supposed ge- optimized but the run time is reasonable enough to consider using the
ological knowledge (repartition of facies in the grid); method on actual data.
• P-velocity Vplog and density ρlog logs upscaled into the reser- This feasibility test on seismic inversion directly in a reservoir model
voir grid to constrain the new realizations to the real proper- was done with the basic formulation for GDM. One of the first things
ties. to consider to optimize this algorithm is the improved GDM formula-
tion (Hu et al, 2004), which consists in combining the initial realiza-
tion with not one, but several complementary realizations. They also
The process needs correlation lengths (given in number of cells) to defined a gradient-oriented search to determine the best complemen-
define the variogram of the complementary realizations. One for the tary realizations for the next iteration and the associated deformation
maximum horizontal range for the x-direction of the reservoir model, coefficient t. Moreover, we can consider the use of parallel processing,
Lx , and one for the y-direction, Ly . There is also a vertical range Lz , by dividing the reservoir model in groups of layers with some layers
defined implicitly in the process, to assure that the realizations by lay- in common between one group and the next. Our current researches
ers do not differ too much from each other. Here, we define Lx = 10 include these improvements.
cells, Ly = 20 cells and Lz = 2 layers.
We also provide a wavelet for the generation of both actual and syn-
thetic seismic datasets. The optimized realizations in the reservoir grid ACKNOWLEDGMENTS
are exported in a Voxet. The ”actual” and synthetic datasets are mod-
eled using a simple 1D convolution technique. Figure 2 shows the The authors would like to thank TOTAL for the release of this work
different seismic datasets. for publication.
Figure 2: From left to right, seismic sections from initial and final
realizations respectively, and ”actual” seismic section.
Combining the Gradual Deformation Method with Seismic Forward Modeling to constrain Reservoir Models
EDITED REFERENCES
Note: This reference list is a copy-edited version of the reference list submitted by the author. Reference lists for the 2008
SEG Technical Program Expanded Abstracts have been copy edited so that references provided with the online metadata for
each paper will achieve a high degree of linking to cited sources that appear on the Web.
REFERENCES
Hu, L.-Y., 2000, Gradual deformation and iterative calibration of Gaussian-related stochastic models: Mathematical Geology, 32,
87–108.
Hu, L.-Y., and M. Le Ravalec, 2004, An improved gradual deformation method for reconciling random and gradient searches in
stochastic optimizations: Mathematical Geology, 36, 703–719.
Le Ravalec, M., B. Noetinger, and L.-Y. Hu, 2000, The FFT moving average (FFTMA) generator: An efficient numerical method
for generating and conditioning Gaussian simulations: Mathematical Geology, 32, 701–723.
Thore, P., 2006, Accuracy and limitations in seismic modeling of reservoir: 76th Annual International Meeting, SEG, Expanded
Abtracts, 1674–1676.
191
Bibliographie
deGroot, P. F. M., Bril, A. H., Florist, F. J., and Campbell, E. A. (1996). Monte carlo
simulation of wells. Geophysics , 61(3), 631638.
Doyen, P. M. (2007). Seismic reservoir characterization : an Earth modelling perspective .
EAGE Publications BV.
Dubrule, O. (2003). Geostatistics for seismic data integration in Earth Models . European
Association of Geoscientists and Engineers.
Dumay, J. and Fournier, F. (1988). Multivariate statistical analyses applied to seismic
facies recognition. Geophysics , 53(9), 11511159.
Duplantier, O., Hadj-Kacem, N., and Vittori, J. (2006). New strategies for seismic facies
upscaling to the reservoir grid scale. SPE Annual Technical Conference and Exhibition,
San Antonio, Texas , (SPE101945).
Egreteau, A. (2005). Etude des variations de l'amplitude de la réectivité du sous-sol après
imagerie sismique en profondeur . Ph.D. thesis, Ecole des Mines de Paris.
Enge, H. D., Buckley, S. J., Rotevant, A., and Howell, J. A. (2007). From outcrop to
reservoir simulation model : Workow and procedures. Geosphere , 3(6), 469490.
Escobar, I., Williamson, P., Cherrett, A., Doyen, P. M., Bornard, R., Moyen, R., and
Crozat, T. (2006). Fast geostatistical stochastic inversion in a stratigraphic grid. SEG
Technical Program Expanded Abstracts , 25(1), 20672071.
Fang, J. H., Karr, C. L., and Stanley, D. A. (1992). Genetic algorithm and its application
to petrophysics. (SPE26208).
Fournier, F. and Derain, J. F. (1992). Seismic data integration in reservoir simulations
through a multivariate statistical calibration approach. SEG Technical Program Expan-
ded Abstracts , 11(1), 9598.
Fournier, F. and Derain, J.-F. (1995). A statistical methodology for deriving reservoir
properties from seismic data. Geophysics , 60(5), 14371450.
Francis, A. M. (2002). Deterministic inversion : Overdue for retirement ? PETEX 2002
Conference and Exhibition, London, UK .
Francis, A. M. (2005). Limitations of deterministic and advantages of stochastic seismic
inversion. CSEG Recorder , pages 511.
Froidevaux, R. (1993). Probability eld simulation. In A. Soares, editor, Geostatistics
Troia 1992 , volume 1, pages 7384. Kluwer, New York.
Gardner, G. H. F., Gardner, L. W., and Gregory, A. R. (1974). Formation velocity and
density - the diagnostic basics for stratigraphic traps. Geophysics , 39, 770780.
Gassmann, F. (1951). Über die elastizität poröser medien. Vierteljahrsschrift der na-
turforschenden geselishaft in Zürich , 96, 123. English traduction available at http:
//sepwww.stanford.edu/sep/berryman/PS/gassmann.pdf.
Gill, D., Shomrony, A., and Fligelman, H. (1993). Numerical zonation of log suites and log
facies recognition by multivariate clustering. AAPG Bulletin , 77, 17811791.
196 Bibliographie
Gnadesikan, R., Blasheld, R. K., Breiman, L., Dunn, O. J., Friedman, J. H., Fu, K.,
Hartigan, J., Kttenring, J. R., Lachenbruch, P. A., Olshen, R. A., and Rohlf, J. (1989).
Discriminant analysis and clustering. Statistical Science , 4, 3469.
Gomez, J. and Srivastava, R. (1990). Isim3d : an ansi-c tree-dimensional multiple indicator
for conditional simulation. Computer and Geosciences , 16, 395410.
Goovaerts, P. (1997). Geostatistics for natural resources evaluation . Applied Geostatistics
Series, New York, USA.
Goovaerts, P. (2002). Geostatistical modeling of spatial uncertainty using p-eld simulation
with conditional probability elds. International Journal of Geographical Information
Science , 16, 167178.
Hadj-Kacem, N. and Pivot, F. (2005). "crafty" techniques for populating géomodèles with
seismic attributes data. EAGE 67th Conference and Exhibition - Madrid, Spain, 13 - 16
June 2005 .
HadjKacem, N. (2006). Intégration des données sismiques pour une modélisation statique
et dynamique plus réaliste des réservoirs . Ph.D. thesis, Ecole des Mines de Paris.
Himberg, J., Korpiaho, K., Mannila, H., Tikanmaki, J., and Toivonen, H. (2001). Time
series segmentation for context recognition in mobile devices. International Conference
on Data Mining (ICDM 2001) - San Jose, USA, 29 November - 2 December 2001 , pages
203210.
Holland, J. H. (1975). Adaptation in Natural and Articial Systems : An Introductory
Analysis with Applications to Biology, Control, and Articial Intelligence . The University
of Michigan.
Hu, L.-Y. (2000). Gradual deformation and iterative calibration of gaussian-related sto-
chastic models. Mathematical Geology , 31(1), 87108.
Hu, L.-Y. and Ravalec-Dupin, M. L. (2004). An improved gradual deformation method
for reconciling random and gradient searches in stochastic optimizations. Mathematical
Geology , 36(6), 703719.
Hu, L.-Y., Ravalec-Dupin, M. L., Blanc, G., Roggero, F., and Noetinger, B. (1999). Redu-
cing uncertainties in production forecasts by constraining geological modeling to dynamic
data. Annual Technical Conference of the SPE, Huston, USA, (SPE56703).
Insalaco, E., Virgone, A., Courme, B., Gaillot, J., Kamali, M., Moallemi, A., Lotfpour, M.,
and Monibi, S. (2006). Upper dalan member and kangan formation between the zagros
mountains and oshore fars, iran : depositional system, biostratigraphy and stratigraphic
architecture. GeoArabia , 11(2), 75176.
Journel, A. G. and Ying, Z. (2001). The theoretical links between sequential gaussian
simulation, gaussian truncated simulation and probability eld simulation. Mathematical
Geology , 33(1), 3140.
Julien, P., Pivot, F., Douillard, A., El-Ouair, Y., and Toinet, S. (2002). The importance
of pre-stack massive seismic modeling for avo calibration and seismic reservoir characte-
rization. SEG Technical Program Expanded Abstracts , 21(1), 17311734.
Bibliographie 197
Kohonen, T. (1990). The self-organizing map. Proceedings of the IEEE , 78(9), 14651480.
Krebes, E. (2004). Seismic forward modeling. CSEG Recorder , pages 2839.
Lecerf, D., Navion, S., and Boelle, J.-L. (2009). Azimuthal residual velocity analysis in
oset vector for waz imaging. 71st EAGE Conference and Exhibition - Amsterdam, The
Netherlands, 8 - 11 June 2009 .
Lecomte, I. and Pochon-Guerin, L. (2005). Simulated 2d/3d psdm images with a fast,
robust, and exible t-based ltering approach. SEG Technical Program Expanded Abs-
tracts , 24(1), 18101813.
Li, Y. and Kareem, A. (1991). Simulation of multivariate non stationary random processes
by t. Journal of Engineering Mechanics , 117(5), 10371058.
Linari, V. (2004). A practical approach to well-seismic data calibration. The Leading Edge ,
23(8), 774775.
Mavko, G., Mukerji, T., and Dvorkin, J. (2003). The Rock Physics Handbook : Tools for
Seismic Analysis of Porous Media . Cambridge University Press.
Menezes, C. and Gosselin, O. (2006). From logs scale to reservoir scale : Upscaling of the
petroelastic model. Annual Conference and Exhibition of the EAGE, Vienna, Austria ,
(SPE100233).
Moyen, R. (2005). Paramétrisation 3D de l'espace en géologie sédimentaire : Le modèle
GEOCHRON . Ph.D. thesis, Institut National Polytechnique de Lorraine, Ecole Natio-
nale Supérieure de Géologie.
Mukerji, T., Jorstad, A., Avseth, P., Mavko, G., and Granli, J. R. (2001). Mapping li-
thofacies and pore-uid probabilities in a north sea reservoir : Seismic inversions and
statistical rock physics. Geophysics , 66(4), 9881001.
Narasimhan, T. N. (1983). A note on volume-averaging. In G. F. Pinder, editor, Flow
through porous media , pages 4650. Computational Mechanics Publication.
Neau, A., Thore, P., and de Voogd, B. (2007). Structural uncertainty eects on reservoir
grid inlling. SEG Technical Program Expanded Abstracts , 26(1), 13971401.
Neau, A., Thore, P., and de Voogd, B. (2008). Combining the gradual deformation method
with seismic forward modeling to constrain reservoir models. SEG Technical Program
Expanded Abstracts , 27(1), 19101914.
Nivlet, P., Lefebvre, F., and Piazza, J. (2007). 3d seismic constraint denition in deep-
oshore turbidite reservoir. Oil and Gas Science and Technology , 62(2), 249264.
Pennington, W. D. (1997). Seismic petrophysics : An applied science for reservoir geophy-
sics. The Leading Edge , 16, 241244.
Pollastro, R. M. (2003). Total petroleum systems of the paleozoic and jurassic, greater
ghawar uplift and adjoining provinces of central saudi arabia and northern arabian-
persian gulf. U.S. Geological Survey Bulletin , 2202-H, 107p. Available at http://
pubs.usgs.gov/bul/b2202-h/b2202-h.pdf.
198 Bibliographie
Résumé :
La caractérisation sismique des réservoirs pétroliers nécessite l'intégration de plusieurs techniques
(lithosismique, géomodélisation, géostatistique et pétrophysique). L'information sismique est d'abord
utilisée pour décrire la structure des réservoirs car son utilisation pour la description des faciès ne
se fait pas sans dicultés. L'objectif de cette thèse est d'apporter des outils nouveaux basés sur
l'utilisation de l'information sismique pour caractériser les réservoirs.
Un premier travail a consisté à évaluer l'impact des incertitudes structurales et leurs conséquences
en terme de classication de faciès. Ensuite, nous considérons la modélisation sismique comme aide
à l'évaluation du modèle réservoir pour faire le lien entre les simulateurs et la réponse sismique du
réservoir. Nous développons ensuite deux approches alternatives aux méthodes traditionnelles en
caractérisation réservoir.
La première utilise les déformations graduelles pour créer des réalisations de propriétés réservoirs.
Ces propriétés sont à l'échelle réservoir, conditionnées aux puits, et respectent une fonction coût
basée sur la comparaison des données sismiques réelles et synthétiques.
La seconde méthode repose sur la classication supervisée par réseaux de neurones pour analyser la
forme des traces sismiques. Une première étape consiste à générer tous les modèles pétrophysiques
envisageables et à les analyser par réseaux de neurones. Les neurones ainsi identiés sont appliqués
aux données réelles, pour trouver des relations pétrophysiques/sismiques identiques aux données
d'apprentissage.
Tous les travaux sont validés sur des champs réels, choisis pour leurs particularités géologiques
(lithologie du réservoir, complexité structurale).
Abstract :
Seismic characterization of hydrocarbon reservoirs is based on various techniques : lithoseismic,
geomodeling, geostatistics, evolutionary algorithms, and petrophysics. Seismic information is rst
used for the description of the reservoir structure, but then its relationship with facies description is
a dicult task. The aim of this thesis is to develop new tools for seismic reservoir characterization.
A rst work has consisted in evaluating the impact of structural uncertainties on petroelastic in-
version and its consequences in terms of facies classication. Then, we consider seismic modeling
as an aid to reservoir model evaluation. This modeling step will make the connection between the
reservoir simulators (or geomodelers) and the seismic response of the reservoir.
Then we develop two alternative approaches for petroelastic and petrophysical inversion. The rst
one uses the gradual deformation method to generate reservoir property realizations. This method
generates properties at the reservoir scale, conditioned by the wells, while respecting a cost function
based on the comparison of actual and synthetic seismic data.
The second method is based on supervised classication principle and uses neural networks to
analyze the waveform of seismic traces. A rst step is to generate a volume containing all possible
petrophysical models for the concerned eld. These models are analyzed by the neural networks.
The neurons identied are applied on the actual data to recognize similar petrophysical/seismic
relationships.
All methods are tested and validated on actual reservoirs, chosen for their specic features (struc-
tural complexity, reservoir lithology).