emphasizing the influence of the inhomogeneous subsurface
Proefschrift
ter verkrijging van de graad van doctor
aan de Technische Universiteit Delft,
op gezag van de Rector Magnificus prof. dr. ir. J.T. Fokkema,
voorzitter van het College voor Promoties,
in het openbaar te verdedigen
op dinsdag 19 september 2006 om 15:00 uur
door
Edith Jacoba VAN VELDHUIZEN
Samenstelling promotiecommissie:
Rector Magnificus, voorzitter
Prof. dr. ir. A.J. Berkhout, Technische Universiteit Delft, promotor Prof. dr. ir. A. Gisolf, Technische Universiteit Delft
Prof. dr. ir. C.P.A. Wapenaar, Technische Universiteit Delft Prof. dr. ir. R.J. Arts, Technische Universiteit Delft Prof. dr. D.G. Simons, Technische Universiteit Delft
Dr. ir. G. Blacqui`ere, TNO Defensie en Veiligheid, Den Haag G. Hampson, B.Sc. Chevron Energy Technology Company, USA
ISBN 90-9020966-2
Copyright c2006, by E.J. van Veldhuizen, Delft University of Technology, Delft, The Netherlands.
All rights reserved. No part of this publication may be reproduced, stored in a retrieval system or transmitted in any form or by any means, electronic, mechanical, photocopying, recording or otherwise, without the prior written permission of the author.
SUPPORT
The research for this thesis was financially supported by the Delphi consortium and by the Netherlands Organisation for Applied Scientific Research TNO.
Typesetting system: LATEX.
1 Introduction 1
1.1 Seismic exploration . . . 1
1.2 Acquisition geometry design . . . 2
1.3 Integrated approach to acquisition geometry analysis . . . 4
1.4 Migration in terms of double focusing . . . 6
1.5 Focal beam analysis as an integrated acquisition geometry analysis method . . . 7
1.6 Objectives and outline of this thesis . . . 8
1.6.1 Objectives . . . 8
1.6.2 Outline . . . 10
2 Focal beam analysis 11 2.1 Operator notation . . . 12
2.2 Imaging by double focusing . . . 13
2.3 Focal beams . . . 16
2.3.1 Focal beams for dipping reflectors . . . 18
2.3.2 Stationary and non-stationary acquisition geometries . . . 20
2.4 Resolution at the target . . . 21
2.5 Angle-dependent amplitude imprint at the target . . . 23
2.6 Focal beam analysis for converted wave data . . . 25
2.7 Conclusions from this chapter . . . 29
3 Implementation and practical aspects 31 3.1 Computing focal beams and focal functions . . . 32
3.2 Overview of 3-D seismic modeling methods . . . 35
3.2.1 Physical modeling and numerical modeling . . . 36
3.2.2 Numerical modeling: an abundance of algorithms . . . 37
3.3 Recursive depth extrapolation in the space-frequency domain . . . . 40
3.3.1 Features of the method . . . 43
3.3.2 Implementation for converted wave acquisition . . . 44
3.3.3 Modeled versus measured point source response . . . 45
3.4 Role of velocity model and migration algorithm: extra options . . . . 50
3.5 Conclusions from this chapter . . . 52
4 Linking survey design to dynamic reservoir characterization 53 4.1 Forward models for AVP data . . . 54
4.1.1 P-P data . . . 54
4.1.2 P-S data . . . 56
4.1.3 Combined P-P P-S data . . . 56
4.2 Inversion for local medium properties: the role of the AVP imprint . 57 4.3 Time-lapse inversion: the link to hydro-carbon saturation . . . 59
4.4 Illustrations . . . 64
4.4.1 P-P and P-S inversion for an OBC acquisition geometry . . . 65
4.4.2 Ziggy model: ‘full 3-D’ versus marine geometry . . . 76
4.5 Conclusions from this chapter . . . 82
5 Observations from examples 83 5.1 P-S OBC acquisition in a model with dipping boundaries . . . 85
5.1.1 Purpose of this example . . . 85
5.1.2 Macro model and analysis set-up . . . 85
5.1.3 Results . . . 86
5.1.4 Conclusions from this example . . . 90
5.2 Model with a weak low velocity lens . . . 92
5.2.1 Purpose of this example . . . 92
5.2.2 Macro model, acquisition geometry, and target points . . . . 92
5.2.3 Results . . . 94
5.2.4 Conclusions from this example . . . 105
5.3 Distinction of different factors of influence . . . 107
5.3.1 Purpose of this example . . . 107
5.3.2 Macro model and analysis set-up . . . 107
5.3.3 Results . . . 109
5.3.4 Conclusions from this example . . . 116
5.4 The Ziggy model: subsalt image quality . . . 117
5.4.1 Purpose of this example . . . 117
5.4.2 Macro model . . . 117
5.4.4 Conclusions from this example . . . 132
5.5 Extended visualization of resolution functions . . . 133
5.5.1 Use of RMS along the time axis . . . 134
5.5.2 From point diffractor to target boundary . . . 135
5.6 Conclusions from this chapter . . . 139
6 Conclusions and recommendations 141 6.1 Conclusions . . . 142 6.1.1 Fundamentals . . . 142 6.1.2 Strategy . . . 142 6.1.3 Method . . . 143 6.1.4 Implementation . . . 143 6.1.5 Specifics . . . 143 6.2 Recommendations . . . 144
A The data matrix in 2-D and 3-D space 145
B Fourier transform in time and space 151
C Focal beams: different formulations 153
Symbols and abbreviations 159
Bibliography 163
Summary 173
Samenvatting 177
Curriculum vitae 181
Introduction
1.1
Seismic exploration
Geophysical methods for the exploration of the earth’s hydrocarbon resources have been used and developed since the early decades of the 20th century. The topic of this thesis is the seismic reflection method. In this method, elastic waves are sent into the subsurface, and subsequently the energy is recorded that arrives back at the surface. This recorded energy is due to reflection, diffraction and refraction at subsurface boundaries. These boundaries are interfaces between layers of the earth that have different acoustic and elastic properties. The seismic data are combined with any available geological information and well data – data acquired by measurements within boreholes – to determine as accurately as possible the subsurface structure and the material properties of the subsurface layers. The seismic method, involving the application of various advanced data processing techniques, has developed into one of the most valuable tools for finding accumulations of oil and gas. Furthermore, when hydrocarbon reservoirs have been located and are being produced, time-lapse seismic surveys are used to monitor the progress of the production.
1.2
Acquisition geometry design
Ideally, in the seismic method downgoing wavefields are generated by seismic sources on a large fine grid at the surface, and the upcoming wavefields are measured by detectors on a fine grid covering a large area of the surface. In practice there are limits to the amount of data that is gathered. Time is limited, budget is limited, the amount of equipment is limited, data storage capacity is limited, data processing capacity is limited, and so on. This means that compromises have to be made. Fur-thermore, all necessary information should be obtainable by some optimum – but limited – amount of measurements; large redundancy should be avoided. The goal is to obtain the best data quality, enabling optimal data processing and interpreta-tion, at the lowest possible cost. This is where acquisition geometry design comes in. Typical parameters are the required minimum and maximum distance (offset ) between source and detector and the optimum spacing between neighboring sources or detectors. Parameters like these determine the range of depths and geological dips that can be imaged, and they also influence the feasibility of velocity estimation and of attenuation of multiple reflections. Many such aspects are treated by Cordsen et al. (2000), Vermeer (2002), and Galbraith (2004).
Conventionally, acquisition geometry analysis and design is based on common mid-point (CMP) analysis – for an overview see Evans (1997). In this approach, the attributes that are used to judge an acquisition geometry are indirect measures like common mid-point fold, offset distribution and azimuth distribution. It is assumed that the earth is horizontally layered and that post-stack migration1will be applied. In many cases the conventional approach to geometry analysis is not sufficient. In areas with a complex subsurface structure, survey design can no longer be based on mid-point analysis, since reflection points are not situated at the mid-point locations (see for example Campbell et al., 2002). The assumption of a horizontally layered earth needs to be replaced by the introduction of a subsurface model. Awareness of this requirement has led to model-based illumination analysis approaches, the model being a macro model as used in seismic migration. For this kind of analysis, often ray tracing is applied to a macro subsurface model, and attributes for survey design are obtained from the ray path information: common reflection point (CRP) analysis is replacing the common mid-point analysis. Counting the number of rays – the hit count – in each CRP bin yields attributes of fold, offset and azimuth that are related to the CRP’s at the target instead of to the CMP’s at the surface. Examples of this approach can be found in Slawson et al. (1994), Muerdter and Ratcliff (2001a), Campbell et al. (2002) and Chang et al. (2002). Dong et al. (2005)
1For an explanation of exploration geophysics terminology, and for an explanation of common
use 2-D wave-equation modeling besides ray tracing to additionally compute the energy distribution at the target.
Furthermore, in industry a shift has occurred from post-stack processing to pre-stack processing, which means that an analysis method based on post-stack processing is no longer of interest. What we are looking for is a method that is in accordance with pre-stack processing; a method that yields direct measures for image quality at the target as they would be obtained by pre-stack processing. So there must be a shift from indirect attributes computed at the surface, to direct measures for the image quality at the target depth.
A laborious way to obtain direct measures for image quality at the target is to carry out a full simulation of the seismic experiment and to migrate the obtained synthetic seismic data. The image quality is visible in the migration result. However, there are more efficient ways to obtain the same information, without this sequence of model-ing and migration. The effect of the acquisition geometry on the spatial resolution of a seismic image can be evaluated by determination of the range of wavenumbers at the target (Beylkin, 1985). For a homogeneous subsurface, this has been investigated and applied by e.g. Vermeer (1999). For an inhomogeneous subsurface, the range of wavenumbers can be determined either by making use of ray tracing (Gibson and Tzimeas, 2002; Lecomte and Pochon-Guerin, 2005) or by using wave-equation based propagators (Toxopeus et al., 2003; Xie et al., 2005). These methods do not require a full forward simulation and a full migration. Moreover, often a
target-oriented approach is followed, i.e. the sought-after attributes are computed for one
or more target boundaries instead of for the full subsurface volume of the macro model. Additionally, simulated Kirchhoff migration amplitudes can be calculated, as described by Laurain et al. (2002). Application of this method is demonstrated by Moldoveanu et al. (2003) and Brink et al. (2004). Since computational power is increasing, it is expected that methods based on the full sequence of 3-D modeling and migration will gain interest (e.g. Jurick et al., 2003). Nevertheless, in practice this laborious sequence is not likely to be used for an acquisition design procedure, since that requires an iterative application. Furthermore, it is not clear how changes in the acquisition geometry should be determined during the iterative procedure of modeling and migration.
struc-wave−equation
hit count use angle information target−oriented full simulation
combined analysis for sources and detectors
separate analysis for sources and detectors acquisition geometry analysis
ray tracing
inhomogeneous subsurface homogeneous subsurface
Figure 1.1: Classification of acquisition geometry analysis methods. The gray boxes lead to the method that is the topic of this thesis.
ture, enabling a better estimation of the attainable image quality at the target. A number of illumination studies is reviewed by Laurain et al. (2004).
1.3
Integrated approach to acquisition geometry analysis
The above described developments in acquisition design can be considered in a broader sense, from the perspective of the seismic value chain (Berkhout, 2004), shown in Figure 1.2. The seismic value chain shows the procedure of seismic explo-ration in a conceptual way, making it easier to see the big picture. The three nodes in the chain, coupled by a double loop, are seismic acquisition, structural imaging, and reservoir characterization. Insight is provided into the interactions between these nodes: the arrows in Figure 1.2 indicate the interactions that take place or that should take place between the different nodes. The arrows that point from left to right indicate ‘influence on’, and the arrows pointing from right to left indicate ‘imposing requirements on’.
reservoir structural imaging characterization seismic acquisition t = T1, T2, . . . Figure 1.2: The seismic value chain.
In addition, the imaging process should be organized such that imaging artifacts are avoided as much as possible, given the available data.
The arrow pointing from ‘structural imaging’ to ‘seismic acquisition’ indicates that structural imaging imposes requirements on the acquisition geometry. The above-mentioned effects, like limited resolution, can be estimated. Thus, requirements can be determined that the acquisition geometry should meet to obtain sufficient structural image quality. This is what is done in most acquisition geometry design procedures.
The third node in the chain is reservoir characterization. Between structural imag-ing and reservoir characterization there is a similar interaction as between seismic acquisition and structural imaging. Imaging artifacts will influence characterization, and a proper characterization therefore requires a sufficient image quality. Whereas structural imaging has a global character, reservoir characterization goes into detail on a local scale. Structural imaging yields an image of the structure of the geology, involving angle-averaged reflection amplitudes. The main goal of reservoir character-ization is to determine local rock and pore parameters from angle-dependent reflec-tion amplitudes. The requirements on the structural imaging process for obtaining sufficiently accurate angle-dependent amplitudes are stricter than the requirements for obtaining angle-averaged amplitudes.
Combining the above-described interactions, links can be made between the begin-ning of the chain and the end of the chain. If the goal is to find rock and pore properties of the target, the acquisition geometry and data pre-processing should be designed such that those properties can be estimated with sufficient accuracy. Gen-erally this will yield stricter design criteria than for structural imaging alone. Since angle-dependent reflection amplitudes are needed in characterization, the acquisition geometry should ensure the conservation of angle-dependent amplitude information at the target.
• the forward path tells that acquisition geometry may introduce imaging
arti-facts, that have to be taken into account in imaging and characterization – the acquisition geometry determines the accuracy of the imaging and characteri-zation result;
• on the other hand, imaging and characterization should impose requirements
on acquisition, to realize the desired accuracy – the required imaging and characterization accuracy determine the acquisition geometry.
In an integrated approach to acquisition geometry design, all three stages of acqui-sition, imaging, and characterization are taken into account.
One aspect of this approach is that the effectiveness of acquisition geometries should be assessed with the specific goal(s) in mind, as different goals require different means. The requirements of the acquisition geometry are determined by the com-plexity of the overburden. For an approximately horizontally layered earth, being a one-dimensional situation, simpler techniques can be applied than for an area of complex geology.
Bear in mind that the link between seismic acquisition and reservoir characterization is made via structural imaging. The description of seismic imaging – or migration – as a double focusing process naturally leads to an acquisition geometry analysis method in which the source geometry and detector geometry are analyzed separately. This is the topic of the next section.
1.4
Migration in terms of double focusing
has propagated from the source(s) to a subsurface location, which we call ‘focusing in emission’. The second-mentioned summation leads to focusing of the wavefield that has propagated from that subsurface location to the detectors, called ‘focusing in detection’. In other words, focusing is achieved by a weighted summation, and double focusing is achieved by two weighted summations.
This physical description of all migration algorithms allows an explicit separation of the two summation steps: focusing of the sources and focusing of the detectors. Considering the two focusing results separately, leads to the introduction of focal
beams, being the intermediate, ‘half-migrated’ results.
1.5
Focal beam analysis as an integrated acquisition geometry
analysis method
In focal beam analysis (Berkhout et al., 2001; Volker et al., 2001), migrated seismic experiments are simulated for specific target points in the subsurface. As mentioned, an important feature is that this is done separately for the source geometry and for the detector geometry.
The forward wave propagation that occurs during a seismic experiment is simulated by modeling a Green’s function for the target point under consideration. Figure 1.3 shows an example of forward modeling of wave propagation between a subsurface target point and the acquisition surface. It shows the energy that travels upward when this point is excited by an incoming wave and acts as a point diffractor, or it can be seen as an “exploding point diffractor”. This simulation alone already shows that the energy from the point diffractor that arrives at the surface clearly varies in strength for different surface locations. This means that at surface locations where little or no energy arrives, the target point cannot be properly detected by a potential detector. It also means that a source placed at such a surface location cannot properly illuminate this target point.
After the forward modeling step, for both the source geometry and the detector geometry a focusing step is done. This yields the focal beams. In short, the
fo-cal source beam represents the contribution of the source geometry to the image
quality for one target point, including focusing artifacts and effects due to the sub-surface structure. Similarly, the focal detector beam shows the contribution of the detector geometry. This gives the opportunity to separately judge and adjust the configuration of the sources and the configuration of the detectors.
angle-dependent amplitude imprint (the focal functions). These attributes are helpful in structural imaging (the second node of the seismic value chain) and in dynamic reservoir characterization (the third node of the seismic value chain), respectively.
x [km] z [km] 14 16 18 20 22 24 26 28 30 32 0 2 4 6
Figure 1.3: ‘Exploding point diffractor’. A 2-D cross-section of the energy of the upward propagating 3-D wavefield is shown. The response at the surface is a Green’s function for the location of the point diffractor. The figure shows that this location cannot be illuminated and detected equally well at all surface locations.
In Figure 1.4 a flow chart for acquisition design is shown. By an analysis procedure, in our case focal beam analysis, the chosen acquisition geometry is evaluated. The subsurface model and the migration operators are taken into account. The output of the analysis, in our case the focal beams and focal functions, should be judged using quality criteria, dictated by characterization as explained in section 1.3. In the following chapters the focal beam analysis method will be thoroughly ex-plained and illustrated.
1.6
Objectives and outline of this thesis
1.6.1 Objectives
The objective of this research is to implement the integrated approach to acquisition geometry analysis using focal beam analysis, explicitly dealing with all three nodes of the seismic value chain.
no yes
final geometry parameters focal beams, focal functions
Focal beam analysis (macro description)
subsurface model
(Green’s functions) migration operators
(sources and detectors) initial geometry parameters
OK? parameters analysis update geometry imaging / characterization requirements
Figure 1.4: Focal beam analysis in an integrated approach to acquisition design (signal part).
• dealing with a 3-D inhomogeneous subsurface; • separate evaluation of source and detector geometry;
• dealing with the requirements imposed by imaging and reservoir
characteriza-tion.
1.6.2 Outline
In Chapter 2 the method of focal beam analysis is explained. It contains the main theoretical aspects, as proposed by Berkhout et al. (2001), and several related aspects, with special emphasis on the applicability of the method to converted wave acquisition.
Chapter 3 deals with the practical implementation of focal beam analysis. The seismic modeling method used in this thesis for modeling of the Green’s functions is recursive depth extrapolation in the space-frequency domain. The method operates in the seismic frequency band. This method is also used to perform the focusing steps. The features of the method are discussed, and the accuracy is investigated by a comparison with physically modeled data. Furthermore, the role of the velocity model and the role of the migration algorithm are discussed.
In Chapter 4 the link is made between the acquisition geometry and the AVP (amplitude versus ray parameter) imprint in dynamic reservoir characterization. The estimated AVP imprint provides a means to analyze the influence of this im-print on inversion of AVP data. This procedure is applicable to both static and dynamic reservoir characterization. The concept is illustrated with synthetic ex-amples, considering an OBC (ocean-bottom cable) experiment and considering a sub-salt situation.
In Chapter 5 the method is illustrated further, by several examples featuring differ-ent subsurface situations and differdiffer-ent types of acquisition geometries. The influence is shown of spatial aliasing, limited acquisition aperture, migration operator limita-tions, and a structurally complex overburden.
Chapter 6 contains the conclusions of this research, and recommendations for future developments.
Focal beam analysis
Focal beam analysis, as briefly explained in 1.5, is a method for assessment of ac-quisition geometries that is directly linked to pre-stack migration. The formulation of migration as a double focusing procedure enables a separation of the assessment of image quality imperfections caused by the source geometry – shown by the focal
source beam – and image quality imperfections caused by the detector geometry –
shown by the focal detector beam. The use of a matrix formulation particularly suits this concept, since this leads to the separation of the ‘source side’ and ‘detector side’ of the seismic image in a natural fashion.
The combined influence of the source geometry and detector geometry on seismic image quality is assessed by joining the focal source beam and focal detector beam. The thus obtained resolution function and AVP imprint – where AVP stands for amplitude versus ray parameter – show the acquisition-related deviation from an ideal structural image and the acquisition-related angle-dependent amplitude im-print, respectively.
2.1
Operator notation
In this thesis, an operator notation is used. In the continuous situation the opera-tors represent integrals. In practice however, the data of a seismic survey is always a discrete spatial sampling of the wavefield. Therefore, seismic data traces can be conveniently arranged in a so-called data matrix (see also Appendix A). The rows and columns of the matrix correspond to spatial locations. Each element of the data matrix corresponds to one frequency component of a seismic data trace. Af-ter removing waves that have traveled along the surface, the data matrix contains frequency-domain signals that can be expressed in terms of wavefield operators de-scribing propagation and reflection in the subsurface. In the notation of the data matrix and the operator matrices, the following conventions are used:
• Matrices are indicated with capital bold symbols, for example P.
• A row vector or column vector that is part of a matrix is indicated by adding
a subscript to the capital bold symbol. The subscript indicates which row or column is taken. The distinction between a row vector or column vector is made by using a dagger superscript (†) when the vector is a row vector:
– Pj is the jth column of matrix P.
– P†j is the jth row of matrix P.
• An element of the matrix is indicated by two subscripts, one that indicates the
row and one that indicates the column: Pij is the element on the ith row and
jth column of matrix P.
Furthermore, the depth levels z that the matrices are related to, are indicated. Examples, which will be discussed more elaborately in the next sections, are:
• S(z0) indicates seismic sources at depth level z0.
• W(zm, z0) is a propagation operator that describes propagation from level z0 to level zm.
• Consequently, S(zm, z0) = W(zm, z0)S(z0) indicates a wavefield at level zm
that results from propagating the source wavefield at level z0 to level zm.
• R(zm, zm) is a reflection operator that gives the relation between a downgoing
wavefield at level zmand an upgoing wavefield at level zm.
2.2
Imaging by double focusing
The focal beam analysis method originates from the ‘imaging by double focusing’ concept (Berkhout, 1982, 1997). This can be explained using the WRW model, describing seismic data in terms of matrices in the frequency domain (Berkhout, 1982). In this model, each monochromatic component P(z0, z0) of the wavefield that is recorded at the surface z0, can be decribed in the space-frequency domain as
P(z0, z0) = D(z0)
K
k=1
[W(z0, zk)R(zk, zk)W(zk, z0)]S(z0). (2.2.1) The depth levels zk are the levels where the waves are reflected. The designation
‘WRW model’ stems from the two matrix operators W and the matrix operator R, which are explained below. Figure 2.1 shows a schematic representation of the model. This equation is valid for stationary acquisition geometries and stationary parts of non-stationary acquisition geometries, as will be explained in section 2.3.2. In that section also the situation for non-stationary acquisition geometries is discussed. The symbols in equation 2.2.1 have the following meaning:
• S(z0): source matrix, containing the source wavelet and the locations of the source arrays. One column represents one source array.
• W(zk, z0): forward wave field propagation matrix. Each column contains a discretized Green’s function describing wave propagation from one point (one lateral location) at the surface z0 to many points at depth level zk.
• R(zk, zk): reflectivity matrix, describing the conversion of an incident wave
field into a reflected wave field. One column of this matrix gives the reflec-tivity relation between the incident wave field at one lateral location and the reflected wave field at all lateral locations; one row of this matrix gives the re-lation between the incident wave field at all lateral locations and the reflected wave field at one lateral location. An extensive explanation of the reflectivity function is given by de Bruin et al. (1990).
• W(z0, zk): forward wave field propagation matrix. Each column contains a
discretized Green’s function describing wave propagation from one point at the reflection level zk to many points at the surface z0.
• D(z0): detector matrix, containing the detector wavelet and the locations of the detector arrays. One row represents one detector array.
detectors. Therefore, one column of this data matrix represents a common source gather (shot record) and one row represents a common receiver gather. Other data gathers, such as CMP gathers or common offset gathers, can also be identified in the data matrix. The data matrix is explained in more detail in Appendix A, for the two-dimensional as well as the three-dimensional situation. The relations between the different data domains are given in Appendix B.
Note that this representation contains only primary reflections. Surface multiples can be included in the model by creating a feed-back loop. Internal multiples can be included by defining the W operators according to the generalized primary rep-resentation as given by Wapenaar (1996). However, since noise is not considered in this thesis, multiples are left aside here.
z0 z1 z2 zk S(z0) D(z0) P(z0, z0) W(z1, z0) W(z2, z0) W(zk, z0) W(z0, z1) W(z0, z2) W(z0, zk) R(z1, z1) R(z2, z2) R(zk, zk)
Figure 2.1: The WRW model, after the formulation of equation 2.2.1. For each depth level
z a WRW-term is added to the total expression for the measured pressure P(z0, z0).
Double focusing is the process of simulating focusing source arrays and simulating
removing the propagation effects. The purpose is retrieval of the reflectivity function, i.e. migration. In the matrix formulation of the WRW model, double focusing can be represented as the application of two focusing operators F to the recorded seismic data:
Pij(zm, zm) = F†i(zm, z0)P(z0, z0)Fj(z0, zm)
= [F†i(zm, z0)D(z0)]W(z0, zk)R(zk, zk)W(zk, z0)[S(z0)Fj(z0, zm)].
(2.2.2) In this expression one reflection level (zk) is considered, and focusing is applied for
depth level zm. The subscript i indicates that the ith row of the matrix is selected. A
row vector is marked by a dagger symbol (†). The subscript j indicates that the jth column of the matrix is selected. Focusing operator F†i(zm, z0) represents focusing from the acquisition level z0 to a subsurface grid point at depth zm with lateral
location (x, y)i, and similarly focusing operator Fj(z0, zm) represents focusing from
the acquisition level z0 to a subsurface grid point at depth zm with lateral location
(x, y)j. Note that [F†i(zm, z0)D(z0)] represents the simulated focusing detector array, and [S(z0)Fj(z0, zm)] represents the simulated focusing source array.
In the confocal version of double focusing, focusing of the sources and focusing of the detectors is applied for the same subsurface grid point (i = j). In the bifocal version of double focusing, the two focusing steps are applied for subsurface grid points with different lateral locations (i = j). If the focusing level is the same as the reflection level (zm=zk), ideally the focusing operators remove all influence of
propagation, W(z0, zm) and W(zm, z0), and acquisition geometry, D(z0) and S(z0). In the confocal case, this yields
Pjj(zm, zm) = Rjj(zm, zm) + εjj(z= zm), (2.2.3)
where Rjj(zm, zm) represents one diagonal element of the reflectivity matrix. εjj(z=
zm) contains the responses from different depth levels. Those responses are out of
focus, but they are part of the full response. An imaging step selecting t = 0 yields a diagonal element of the reflectivity function, being the angle-averaged reflection coefficient for the grid point under consideration. By application of this imaging condition the responses from different depth levels are discarded. The bifocal result of double focusing to depth level zmis
Pij(zm, zm) = Rij(zm, zm) + εij(z= zm), (2.2.4)
where Rij(zm, zm) represents one off-diagonal element of the reflectivity matrix (i=
j). After Radon transformation of the bifocal results for a range of grid points, the
angle-dependent reflection coefficients can be found at zero-intercept time (Berkhout, 1997).
2.3
Focal beams
In expression 2.2.2, all terms on the right-hand side of the reflectivity matrix R(zk, zk)
concern the source side of the seismic experiment. These terms represent the down-ward propagation from the focusing source array. This part of the double focusing expression we call the focal source beam Sj (Berkhout et al., 2001):
Sj(zk, zm) = W(zk, z0)S(z0)Fj(z0, zm), (2.3.5)
with k = 1, 2, ..., M . Figure 2.1 shows how the depth levels in the medium are defined. All terms in expression 2.2.2 that are on the left-hand side of the reflec-tivity matrix R(zk, zk) concern the detector side of the seismic experiment: the
upward propagation to the focusing detector array. This part of the double focusing expression we call the focal detector beam D†j:
D†j(zm, zk) = F†j(zm, z0)D(z0)W(z0, zk), (2.3.6)
with k = 1, 2, ..., M . Alternative – computationally useful – formulations of the focal beams are given in Appendix C. For both the focal detector beam and the focal source beam, choosing k = m implies that the depth slice of the focal beam at the depth of the focal point is considered. In an ideal situation
Sj(zm, zm) = Ij(zm), (2.3.7)
and
D†j(zm, zm) = I†j(zm). (2.3.8)
I(zm) is an identity matrix. Consequently, Ij(zm) is a column vector where the
jth element is the only one non-zero element, and I†
j(zm) is a row vector where the
jth element is the only one non-zero element. Expressions 2.3.7 and 2.3.8 indicate that all influence of acquisition geometry and propagation is perfectly removed by focusing.
the angles of incidence at the target point as a result of the chosen source geometry, and the focal detector beam shows the take-off angles at the target point that can be detected by the detectors. In the Radon domain the focal beams can be presented as a function of the lateral ray parameter components px and py. The linear Radon
transformation of the focal detector beam, with respect to target point (x, y)j, is
achieved as follows: Dj(p, zm, τ ) =∆x∆y yn xn Dj(xj− xn, zm, ω)ejω[τ +p·(xj−xn)], (2.3.9)
where τ is the intercept time, ω is the radial frequency, x = (x, y), and p = (px, py).
∆x and∆y are the discrete sampling intervals of x and y. Note that the ray parameter vector p here contains the lateral ray parameter components. The tilde symbol (∼) indicates the Radon domain. Equation 2.3.9 gives the Radon-transformed focal beam for lateral location xj, indicated by the subscript j. Note that the Radon transform
is shifted such that location xj is the lateral origin of the coordinate system. Since
we are only interested in zero intercept time, the Radon transform simplifies to Dj(p, zm, τ = 0) =∆x∆y yn xn Dj(xj− xn, zm, ω)ejωp·(xj−xn). (2.3.10)
The same result can be obtained by a Fourier transform of Dj(x, zm, ω), yielding
Dj(k, zm, ω), with k = (kx, ky), followed by mapping of the lateral wavenumbers kx
and ky to lateral ray parameters px and py. In matrix notation, expression 2.3.10
reads
D†j(zm, zm) = Dj†(zm, zm)ΓT, (2.3.11)
where each element of Radon transformation operator Γ has the form
Γkn=∆x∆y ejωpk·(xj−xn) (2.3.12)
The length of row vector D†j(zm, zm) is NxNy, where Nxis the number of x-samples
and Ny is the number of y-samples. The length of row vector D†j(zm, zm) is KxKy,
where Kxis the number of px-samples and Kyis the number of py-samples. The size
of operator Γ is KxKy by NxNy. Bear in mind that the two lateral coordinates x
and y are arranged along one matrix dimension, see figures A.4 and A.5 in Appendix A. The coordinates in the Radon domain, pxand py, are organized similarly.
The Radon transformation of the focal source beam is analogous. The matrix for-mulation of the Radon transform of the focal source beam reads
Sj(zm, zm) = ΓSj(zm, zm), (2.3.13)
plane wave, and its lateral components is
px=|p| sin φ cos θ = |p| cos φ tan α, (2.3.14)
py =|p| sin φ sin θ = |p| cos φ tan β, (2.3.15)
where θ is the azimuth, and φ is the angle with the normal to the horizontal plane. The angles α and β are the angles with the normal to the horizontal plane as observed in the (px, pz)-plane and (py, pz)-plane, respectively. This is illustrated in Figure 2.2.
α β φ θ px py pz p
Figure 2.2: Angles involved with ray parameter p = (px, py, pz), which is denoted by the vector. The relations between the angles are given by equations 2.3.14 and 2.3.15. The angles are measured counterclockwise.
The focal beams have been defined with respect to a single target point. However, for a global view of the illumination effects across a target boundary, focal beams can be created for a distribution of target points on the target bounday.
2.3.1 Focal beams for dipping reflectors
In the acquisition geometry analysis an angle-independent reflector at target depth is considered instead of the true target reflectivity. This is because we are only interested in angle-dependent effects due to the acquisition geometry, and not in the effects due to the target reflectivity. Furthermore, the reflector is assumed to be horizontal, as the focal beams are computed for a constant depth level zm. However,
detector beam. This effect can be accounted for by shifts of the angles α and β of expressions 2.3.14 and 2.3.15. The angles are adjusted according to
α→ α − αd, (2.3.16)
β → β − βd, (2.3.17)
where the subscript d stands for ‘dip’, and where αd and βd are the dip angles of
the reflector in x-direction and y-direction, respectively. This is illustrated in Figure 2.3, for the (px, pz)-plane only; the situation for the (py, pz)-plane is similar.
(a) (b) reflector reflector plane wavefront x x z z αd αd α α αeff αeff px px,d pz pz,d
Figure 2.3: When computing focal beams, initially a horizontal reflector is assumed (a). To obtain the focal beam for a dipping reflector with local dip α (b), in the Radon domain the
horizontal ray parameter components pxand pyare adjusted according to the dip angle, such
that px= py= 0 always corresponds to a plane wave perpendicular to the reflector plane.
The effective angle αeff is the difference between the angle α for a horizontal reflector and
the dip angle αd. For simplicity, only the (x, z)-plane and corresponding (px, pz)-plane are
shown; the situation for the (y, z)-plane is similar. The subscript d denotes ‘dip’. The angles are measured counterclockwise.
azimuth and angle:
px,d=|p|
sin φ cos θ− cos φ tan αd
1 + tan φ cos φ tan αd
, (2.3.18)
py,d =|p|
sin φ sin θ− cos φ tan βd
1 + tan φ sin φ tan βd
. (2.3.19)
These expressions reduce to expressions 2.3.14 and 2.3.15 for αd = 0 and βd =
0, respectively. This adjustment of the ray parameter components ensures that
px= py= 0 always corresponds to a plane wave perpendicular to the reflector plane.
Note that a dipping reflector will affect both azimuth and angle. As a result, for different dips the focal detector beam and the focal source beam occupy different areas in the Radon space.
The incidence angles and take-off angles as shown by the focal beams can be com-bined to specular reflection angles. This combination of the focal detector beam and the focal source beam shows the angle-dependent amplitude imprint on target reflectivity as caused by the overburden, acquisition geometry, and migration opera-tors. The result, the AVP imprint, will be discussed more elaborately in section 2.5 and Chapter 4. Since the reflector dip can be included in the focal beams with little effort, the angle-dependent amplitude imprint at the target can easily be evaluated for ranges of different reflector dips.
2.3.2 Stationary and non-stationary acquisition geometries
The focal beams can be computed for any configuration of sources or detectors. However, for assessment of the combined effect of the sources and detectors, the focal beams should be computed for stationary parts of the acquisition geometry. The focal source beam concerns the wave path between the sources and the target point, and the focal detector beam concerns the wave path between the target point and the detectors. When combining the focal source beam and focal detector beam, the wave paths are joined. As a result, wave propagation between the sources and the detectors is considered. Therefore, only sources and detectors should be included between which there is wave propagation. In practice this means that all sources can be included that correspond to the same active detector configuration. These sources and the corresponding detectors form one stationary part of the survey. In marine surveys the location of the hydrophones will be different for each shot, so one shot at a time will be considered. The WRW model for one shot reads
Pl(z0, z0) = D[l](z0)
K
k=1
Note that the detector matrix depends on the location of the source array. It contains the detector arrays corresponding to source array l.
The focal beams are formulated accordingly:
S[l]j (zk, zm) = W(zk, z0)Sl(z0)Flj(z0, zm), (2.3.21)
and
D†[l]j (zm, zk) = F†j(zm, z0)D[l](z0)W(z0, zk). (2.3.22)
To get a qualitative impression of the ‘illumination capability’ of all sources included in the survey, the focal source beams for stationary parts can be summed:
S[L]j (zk, zm) =
l
S[l]j (zk, zm), (2.3.23)
where l denotes one stationary part of the survey and L denotes the full survey. Similarly, to see how well the target point can be detected potentially by all detectors included in the survey, the focal detector beams for stationary parts should be summed:
D†[L]j (zm, zk) =
l
D†[l]j (zm, zk). (2.3.24)
However, in the computation of the focal functions – the resolution function and AVP imprint, see sections 2.4 and 2.5 – these ‘full survey focal beams’ are not used. For non-stationary acquisition geometries, focal functions are determined for each stationary part of the survey. Subsequently, the focal functions for the complete survey are obtained by a summation of the focal functions for all stationary parts.
2.4
Resolution at the target
The resolution function is the confocal imaging result for a point diffractor at target point j at depth zm:
δjPih(zm, zm) = F†i(zm, z0)δjP(z0, z0)Fh(z0, zm), (2.4.25)
with i = h varying around focal point j. δjP(z0, z0) is the seismic response of a unit point diffractor:
δjP(z0, z0) = D(z0)W(z0, zm)δjR(zm, zm)W(zm, z0)S(z0). (2.4.26)
δjR(zm, zm) represents a unit point diffractor at grid point j:
δjRhi= 1 i = h = j
δjRhi= 0 i= j or h = j,
δjR = 0 0 . . . . . . 0 0 1 0 . . . 0 .. . 0 0 . . . 0 .. . ... ... . .. ... 0 0 . . . . . . 0 j j h-i ?
Each value of i (=h) in equation 2.4.25 gives one element of the total resolution function. Considering the formulation of the focal beams as given in equations 2.3.6 and 2.3.5, it appears that each element of the resolution function is obtained by using a different pair of a focal detector beam and a focal source beam. In practice the computation of the resolution function can be simplified by using the time-reversed versions of the focal beams, equations C.4 and C.3. In that case the resolution function can be obtained by an element-by-element multiplication of one focal source beam and one focal detector beam. This can be seen by substituting equation 2.4.26 and the reflectivity matrix δjR(zm, zm) (expressions 2.4.27) in equation 2.4.25:
δjPih(zm, zm) = F†i(zm, z0)D(z0)W(z0, zm)δjR(zm, zm) × W(zm, z0)S(z0)Fh(z0, zm) = F†i(zm, z0)D(z0)Wj(z0, zm) × W†j(zm, z0)S(z0)Fh(z0, zm) = Dij(zm, zm)Sjh(zm, zm), (2.4.28)
for h = i =−N . . . N, where 2N + 1 is the number of points in the target area that is included in the focal beams. In this expression, Dij(zm, zm) and Sjh(zm, zm) are
single elements that correspond to the elements of the focal beams of equations C.4 and C.3, respectively. So indeed the resolution function is an element-by-element multiplication of the focal detector beam and the focal source beam in the space-frequency domain.
For a non-stationary acquisition geometry, the resolution functions for all stationary parts need to be summed (see section 2.3.2):
δjP [L] ih(zm, zm) = l δjP [l] ih(zm, zm) = l D[l]ij(zm, zm)S [l] jh(zm, zm) (2.4.29)
Note that in equation 2.4.25 the resolution function is formulated for a single fre-quency component. If a seismic frefre-quency band is considered, the resolution function is obtained by applying an imaging condition, i.e. by summing all monochromatic resolution functions. The resolution function is the seismic image of the point dif-fractor at the target point, at depth zm. The obtained result shows the influence of
the source geometry, the detector geometry, the migration operators, and the struc-ture of the overburden on the lateral image resolution at depth zm. The influence
on the vertical image resolution, and the artifacts occurring at different depth levels, can be assessed by computing the resolution functions for a range of depth levels.
2.5
Angle-dependent amplitude imprint at the target
We will now consider a plane reflecting target boundary. As mentioned in section 2.2, in bifocal imaging the two focusing steps are applied for subsurface grid points with different lateral locations, ideally yielding off-diagonal elements of the reflectivity matrix (equations 2.2.2 and 2.2.4). If the grid point j, where the sources are focused, is kept constant, while grid point i, where the detectors are focused, is varied, a column of the reflectivity matrix R(zk, zk) is obtained. This column vector contains
the angle-dependent reflectivity operator for grid point j. The process of obtaining angle-dependent reflectivity by pre-stack migration is described by de Bruin et al. (1990). As mentioned, in the Radon-transformed bifocal imaging result the angle-dependent reflection coefficients can be found at zero intercept time.
If acquisition and focusing are not perfect, however, the true angle-dependent re-flection coefficients are not obtained, but the angle-dependent rere-flection coefficients with an imprint that is due to the imperfect acquisition and imperfect focusing. In focal beam analysis, the goal is not to find the angle-dependent reflectivity of the target, but to find this angle-dependent amplitude imprint or AVP imprint, where AVP stands for amplitude versus ray parameter. Therefore, in focal beam analysis the AVP imprint is modeled using an angle-independent reflector at target depth. The AVP imprint is the Radon transformed bifocal imaging result for a range of grid points on that reflector. The bifocal imaging result for one pair of grid points is
∆Pij(zm, zm) = F†i(zm, z0)∆P(z0, z0)Fj(z0, zm), (2.5.30)
where the i and j indicate lateral coordinate pairs (x, y)i and (x, y)j at the
tar-get depth level zm. ∆P(z0, z0) is the seismic response of a reflector with angle-independent (and unit) reflectivity – in that case the reflectivity matrix R is an identity matrix:
For a range of grid points i, the bifocal imaging result is
∆Pj(zm, zm) = F(zm, z0)∆P(z0, z0)Fj(z0, zm), (2.5.32)
where each row of F(zm, z0) contains a focusing operator for one grid point i. The bifocal imaging expression can be written as the product of a focal source beam and a focal detector beam – both are vectors – by combining equations 2.5.30, 2.5.31, 2.3.6, and 2.3.5. For one grid point i, this yields
∆Pij(zm, zm) = D†i(zm, zm)Sj(zm, zm), (2.5.33)
and for a range of grid points i
∆Pj(zm, zm) = D(zm, zm)Sj(zm, zm), (2.5.34)
where every row of D(zm, zm) contains the focal beam for a grid point i. The AVP
imprint for grid point j can be obtained from the bifocal imaging results for a range of grid points i. Therefore, it appears that the AVP imprint is obtained by using one focal source beam, for grid point j, and many focal detector beams, for all grid points i. In practice the computation of the AVP imprint is simplified by assuming that the focal detector beam does not vary as a function of location i around the target point j, and by using relation C.18. In that case for each point i a shifted version of focal detector beam D†j(zm, zm) can be used, and the bifocal imaging
result can be written as
∆Pj(zm, zm) = D(zm, zm)Sj(zm, zm)
= DH(zm, zm)Sj(zm, zm)
= [DT]∗(zm, zm)Sj(zm, zm).
(2.5.35)
This equation reveals that each element of ∆Pj(zm, zm) is obtained by a
space-invariant convolution: ∆Pij(zm, zm) = N v=−N D∗j(i−v)(zm, zm)Svj(zm, zm), (2.5.36)
which can be written in its functional form as ∆Pj(xi, zm) =
N
v=−N
D∗j(xi− xv, zm, ω)Sj(xv, zm, ω), (2.5.37)
where x = (x, y). This convolution is a multiplication in the Radon domain. Fur-thermore, complex conjugation in the space-frequency domain implies reversal of the ray parameter in the Radon domain. Consequently, the AVP imprint is
where p = (px, py). For a non-stationary acquisition geometry, the AVP imprints
for all stationary parts need to be summed (see section 2.3.2): ∆ Pj[L](p, zm) = l ∆ Pj[l](p, zm) = l D[l]j (−p, zm) S [l] j (p, zm), (2.5.39)
where l denotes a single stationary part of the survey and L denotes the full survey. In short, the AVP imprint can be approximated by an element-by-element multipli-cation of the Radon-transformed focal detector beam and focal source beam with focal point location j. The broad-band AVP imprint is obtained by summing the monochromatic AVP imprints along lines of constant ray parameter1 (de Bruin et al., 1990). The result represents the angle-dependent amplitude imprint on target reflectivity. This imprint is caused by the source geometry, the detector geometry, and the migration operators, where the influence of the structure of the overburden is taken into account. In Chapter 4 the significance of the AVP imprint will be discussed in the light of reservoir characterization.
2.6
Focal beam analysis for converted wave data
In industry, there is a growing interest in converted wave data, although the added value of converted wave acquisition in practical situations is still uncertain at the time of writing. In this section I will briefly describe some aspects and issues of converted wave acquisition and processing, and explain how these relate to focal beam analysis.
I will limit this discussion to P-S waves: converted waves where the downgoing wave is a P-wave, conversion occurs during reflection at a subsurface target boundary, and the upgoing wave is an S-wave. Multiply converted waves that, for example, travel through a salt body as an S-wave and through the surrounding sediments as a P-wave do have useful properties (Al´a’i and Verschuur, 2006), but the signal strength of those wave types generally is low compared to P-P waves and P-S waves (Stewart et al., 2002). The cause of this is that for transmission of waves through an interface, most of the wavefield continues as the same wave type and only a small portion of the energy is converted (Aki and Richards, 2002). However, the reflection coefficients for reflection with or without conversion have the same order of magnitude – although there is no conversion for normal incidence. Therefore, the wave types most likely to be measured at the surface are the ones that have converted only at reflection and
1In this procedure, only positive frequencies are taken into account, in accordance with the
not at transmission, and the ones that have not converted. In practice this comes down to P-S and P-P, respectively.
The main distinguishing aspects of P-S wave propagation compared to P-P wave propagation are (Caldwell, 1999; Stewart et al., 2002)
• the lower propagation velocity – and therefore smaller wavelengths – of the
S-waves, causing asymmetrical wave paths with respect to the reflection point,
• the different angle-dependent amplitude behavior for reflection and
transmis-sion at a subsurface interface, and
• the difference in how P-wave velocities and S-wave velocities are affected by
rock and pore properties.
These aspects account for the potential added value of converted wave data in seismic exploration. Subsurface areas that are not illuminated by P-P waves, might in some cases be illuminated by P-S waves because of the different wave paths. Furthermore, the smaller wavelength of the S-waves can result in a better image resolution.2 With respect to the second aspect, S-wave properties of rocks may show greater variation than P-wave properties, providing a significant contrast and reflection (Stewart et al., 2002). Apart from the greater variation in rock properties, P-S reflectivity as a function of angle has a different nature than P-P reflectivity. This can yield better results for the inversion of angle-dependent P-S amplitudes than for the inversion of P-P amplitudes (see sections 4.1 and 4.4.1). The third aspect can result in a better P-S imaging than P-P imaging in the presence of gas clouds, because S-waves are less affected by pore properties than P-waves (Stewart et al., 2003; Knapp et al., 2002).
The difference in illumination and the increased image resolution can be assessed by focal beam analysis. Moreover, focal beam analysis can assess the results of inversion of P-P AVP data and P-S AVP data as differently affected by the AVP imprint. This will be explained in Chapter 4.
Focal beam analysis for converted wave acquisition is not very different from focal beam analysis for P-wave acquisition. It makes use of the generalization of the WRW model for multi-component data as described by Berkhout and Verschuur (2000), which was applied earlier by Winthaegen and Verschuur (2002) and Cox (2004). For converted wave data, the double focusing expression of equation 2.2.2 is written as
PSP,ij(zm, zm) =
[F†SS,i(zm, z0)DS(z0)]WSS(z0, zk)RSP(zk, zk)WPP(zk, z0)[SP(z0)FPP,j(z0, zm)].
(2.6.40)
2In practice, the improved resolution for P-S acquisition is not observed, due to the different
F†SS,i(zm, z0) represents S-wave focusing without mode conversions, DS(z0) is the de-tector matrix for detection of S-waves, WSS(z0, zk) represents S-wave propagation
without mode conversions, RSP(zk, zk) is the reflectivity matrix for reflection with
conversion of P-waves to S-waves, WPP(zk, z0) represents P-wave propagation with-out mode conversions, SP(z0) is the P-wave source matrix, and FPP,j(z0, zm)
rep-resents P-wave focusing without mode conversions. Consequently, the focal source beam for converted wave acquisition is given by
SPP,j(zk, zm) = WPP(zk, z0)SP(z0)FPP,j(z0, zm), (2.6.41)
with k = 1, 2, ..., M , and the focal detector beam for converted wave acquisition is given by
D†SS,j(zm, zk) = F†SS,j(zm, z0)DS(z0)WSS(z0, zk), (2.6.42)
with k = 1, 2, ..., M . The focal beams are evaluated for a target point j at depth level zm, where j indicates a set of lateral coordinates (x, y)j. For converted wave
acquisition this target point j is the conversion point at the reflecting target bound-ary. Equations 2.6.41 en 2.6.42 show that in focal beam analysis for converted wave data, the P-waves and the S-waves are treated separately.
The creation of the resolution function for converted wave data is straightforward. The resolution function represents the imaged response of a point diffractor, and by combining the focal source beam – for P-wave propagation – and the focal detector beam – for S-wave propagation – in the same way as would be done for the P-P case (equation 2.4.28), implicitly a P-S diffraction surface is migrated instead of a P-P diffraction surface.
The computation of the AVP imprint for converted wave data is straightforward as well. For P-S acquisition, the Radon-transformed focal source beam shows the an-gles of P-wave illumination, and the Radon-transformed focal detector beam shows S-wave detection angles. As explained in section 2.5, the AVP imprint is obtained by an element-by-element multiplication of the Radon-transformed focal detector beam and focal source beam. The elements correspond to ray parameter values. As the ray parameters are conserved at reflection – independently of wave type – this multiplication holds for both P-P and P-S acquisition. The AVP imprint involves factors of influence above the target, and not the angle-dependent reflectivity behav-ior of the target boundary itself. Therefore, for computation of the AVP imprint, the nature of the reflection at the target is irrelevant – naturally it is of importance when assessing the consequences of the AVP imprint.
Besides the possible advantages of P-S exploration there are also various issues re-lated to converted wave (marine) seismics, some of which are listed here:
the ocean bottom to enable recording of S-waves. An accurate positioning can be challenging in deep-water areas.
• The detectors should have sufficient vector fidelity to enable an accurate
iden-tification of the different wave modes.
• Large offsets are needed to record P-S signals with a sufficient amplitude,
since there is no conversion at normal incidence and the reflection amplitudes increase with angle.3
• Because of the smaller wavelengths, P-S waves need a denser sampling of
detectors than P-P waves to avoid spatial aliasing.
• In processing, sensitivity to velocity errors leads to mis-positioning of the
con-version point and of the image.
Focal beam analysis relates to these matters in different ways. Signal quality issues related to decomposition and vector fidelity are not taken into account, and position-ing is assumed to be precise. On the other hand, the requirement of large offsets and dense sampling can be analyzed and demonstrated by focal beam analysis. These are the typical acquisition geometry characteristics that can be assessed by focal beam analysis. Determination of the location of the conversion point is not an issue, since the conversion point location is the chosen target point location. The image artifacts caused by errors in the velocity model are not addressed in this thesis. However, the role of the velocity model in focal beam analysis will be discussed in Chapter 3 (section 3.4), where we see that velocity sensitivity analysis is an extra – but yet unexplored – option of focal beam analysis.
As a final remark I should mention that for P-S acquisition, it is important to con-sider anisotropy (Thomsen, 1999), since S-waves generally have a stronger anisotropy than P-waves. This is beyond the scope of this work, however.
2.7
Conclusions from this chapter
Focal beam analysis is a method for assessment of subsurface image quality as af-fected by the seismic acquisition geometry. The distinguishing features are:
• It has a direct link to pre-stack migration: migration – in terms of double
focusing – is an intrinsic part of the method. This important feature ensures that the analysis result is directly related to the seismic image quality as it would be obtained by pre-stack migration.
• It is based on a separate treatment of the source geometry and the detector
geometry. This enables a separate assessment of the effects on image qual-ity caused by the source geometry and the detector geometry. Consequently, options for improvement of the geometry can be identified separately for the detector side and the source side of the acquisition. Moreover, knowing the separate results provides more insight in the combined result – showing the combined effects from both sides – since it is known which effects are due to the source geometry and which are due to the detector geometry.
Implementation and practical
aspects
In the previous chapter, the theory of focal beam analysis was explained. No infor-mation was given how the focal beams and focal functions are computed. This will be done in this chapter.
In Chapter 1, several seismic modeling methods have been mentioned that are used in illumination studies and in resolution analysis. The method of choice depends on the goals and on the available computation power and time. For homogeneous media and simple inhomogeneous media the focal beams can be computed efficiently by a Fourier transform (Volker et al., 2000). However, the emphasis on complex subsurface structures leads to the choice for a recursive depth extrapolation method based on the one-way wave equation.
3.1
Computing focal beams and focal functions
To use the concept of focal beam analysis in practice, a method has to be chosen to compute the focal beams. This means that the expressions 2.3.5 and 2.3.6 (see section 2.3) need to be evaluated.
Both the expression for the focal detector beam and the expression for the focal source beam contain a forward propagation term (W) and a focusing term (F). For the simulation of forward wave propagation, a method should be chosen that models the true wave propagation as accurately as possible. For focusing, a method should be chosen that corresponds to the migration algorithm that is applied in practice. From a computational point of view it is convenient to use the alternative definitions of the focal beams instead of the definitions given in section 2.3. This means that expressions C.3 and C.4 are evaluated instead of expressions 2.3.5 and 2.3.6. As a consequence, the time-reversed focal beams are computed. To obtain the focal beams as given in equations 2.3.5 and 2.3.6, the time axis of the computed focal beams should be reversed. By choosing this approach, the computation of focal beams consists of the following steps:
• Define a point diffractor at the target point.
• Simulate forward wave propagation from the target point to the acquisition
level (W†j(zm, z0) and Wj(z0, zm)). This yields a Green’s function.
• Apply a sampling function. For computation of the focal source beam, the
traces are selected that correspond to the locations of the sources. For com-putation of the focal detector beam, the traces are selected that correspond to the locations of the receivers.
• Apply focusing (F(z0, zk) and F(zk, z0)) from the acquisition level to the target level, for an area around the target point.
This area around the target point should have a sufficient extent to include the locations where the main imaging artifacts occur (e.g. spatial aliasing), leading to a range of lateral locations for which focusing is applied.
The application of the sampling function and the focusing step are repeated for all stationary parts of the acquisition geometry, as explained in section 2.3.2.
This relation between the focal beams and the focal functions is shown diagrammat-ically in Figure 3.1.
Focal detector beam Focal source beam Focal functions
Spatial domain (x, y, ω) Radon domain (px, py, ω) multiply multiply Radon transform Radon transform Dj(zk, zm) Sj(zm, zk) Dj(zk, zm) Sj(zm, zk) δjPih(zk, zk) i = h = 0 . . . N ∆ Pj(zk, zk) x x x y y y px px px py py py k = m k = m
Figure 3.1: The relations between the focal detector beam, the focal source beam and the focal functions (resolution function at the top, and AVP imprint at the bottom), at focusing
level zk. px and py are the horizontal ray parameter components. The multiplications are
element-by-element multiplications. Note that the AVP imprint is not the Radon trans-formed resolution function.
As explained in section 2.3.1, the local target dip can easily be included in the computation by applying a ray parameter shift to the focal detector beam and the focal source beam in the Radon domain.
3.1.1 The ideal focal beam and ideal focal functions
