Master Thesis in Geosciences
Imaging the sea surface from towed dual-sensor streamer data
A feasibility study
Orji Okwudili Chuks
Time varying sea surface image
Imaging the sea surface from towed dual-sensor streamer data A feasibility study
Orji Okwudili Chuks
Master Thesis in Geosciences Discipline: Geophysics Department of Geosciences
Faculty of Mathematics and Natural Sciences
UNIVERSITY OF OSLO
[
June 2009]
© Orji Okwudili Chuks, 2009
Tutor(s): Dr. Walter Söllner (PGS) and Prof. Leiv J. Gelius (UiO) This work is published digitally through DUO – Digitale Utgivelser ved UiO http://www.duo.uio.no
It is also catalogued in BIBSYS (http://www.bibsys.no/english)
All rights reserved. No part of this publication may be reproduced or transmitted, in any form or by any means, without permission.
I Acknowledgements
I thank God Almighty for his graces. My gratitude goes to my supervisors Dr. Walter Söllner and Prof. Leiv J. Gelius for their help and advices.
I am grateful to PGS and its staff, especially the R&D team.
I thank my wife and son for their patience and love throughout this work. I am immensely grateful to my parents and siblings with whose help I have come this far. I also like to thank my friends and classmates for their friendship.
Orji, Okwudili Chuks June 2009.
II
Abstract
Sea surface profile and reflection coefficient estimates are vital input parameters to various seismic data processing applications. The common assumption of a flat sea surface when processing seismic data can lead to misinterpretations and mislocations of events. A new method of imaging the sea surface from decomposed wavefields is presented. Wavefield separation is applied to the data acquired by a towed dual sensor streamer containing collocated hydrophones and geophones to obtain the up- and down-going wavefields of the related sensors. The up- and down-going wavefields of a given sensor are extrapolated to the sea surface where an imaging condition is applied in order to obtain the sea surface profile.
The image points values obtained at the sea surface are the reflection coefficient values of these points. Ray tracing and finite difference methods are used to generate different controlled data sets employed in this feasibility study to demonstrate the imaging principle and to test the image accuracy. Finally, a first field data example of a marginal weather line from the Norwegian North-Sea is presented.
III
Contents
Introduction 1
1 The marine seismic ghost 3
1.1 Introduction 3
1.2 Marine seismic receiver ghost 3
1.3 Relationship between hydrophone and geophone receiver ghosts 4 1.4 Detrimental effects of the marine seismic receiver ghost 7
1.5 Deghosting approaches 10
1.6 Dual-sensor streamer (Geostreamer) 11
2 Sea surface imaging technique 14
2.1 Wavefield separation 15
2.2 The low-frequency condition 19
2.3 Wavefield extrapolation 21
2.4 Imaging condition application 24
3 Synthetic data examples 28
3.1 Imaging based on test data from ray tracing 28
3.1.1 Definition of controlled model and data generation 28
IV
3.1.2 Feasibility study of wavefield decomposition and extrapolation 31
3.1.3 The imaged sea surface 43
3.1.4 Error analysis 44
3.1.5 Sensitivity/resolution of the imaging technique 47 3.2 Imaging based on test data from time domain finite difference (FDTD) calculations 48 3.2.1 Definition of controlled models and data generation 48
3.2.2 Quality control of the synthetic data 51
3.2.3 Flat sea surface (First case) 54
3.2.4 Simple sea surface (Second case) 56
3.2.5 Complex/rough sea surface (Third case) 58
4 Field data examples 61
4.1 Field data analysis 61
4.2 “Frozen” sea surface image 63
4.3 Time variant behaviour of the sea surface shape 65
5 Conclusions and future work 69
5.1 Conclusions 69
5.2 Future work 70
Appendices 72
Appendix A: Redatuming (Non-flat streamer case) 72
V
Appendix B: Method of processing generated synthetic and field data sets 73
Appendix C: A brief summary of dynamic ray tracing 75
Appendix D: A brief review of explicit time domain finite difference method 77
References 81
1
Introduction
Originally, oil exploration was generally a guessing game based on surface observations.
However, various remote sensing methods have been developed through the years that can be used to search for and evaluate subterranean formations associated with possible hydrocarbon accumulations. Seismic surveying represents the most important geophysical prospecting method, which involves sending signals into the earth and receiving them again at the surface. In a land-based seismic survey, a seismic signal is generated on or near the earth’s surface. These signals then travels downward into the subsurface of the earth while in a marine seismic survey the seismic signals also travel through a body of water overlying the earth’s subsurface. Seismic energy sources like vibrators and explosives for land-based case or airguns and marine vibrators etc., for marine case, generates the seismic signals which after propagating into the earth is partially reflected by the subsurface seismic reflectors (acoustic impedance boundaries) and is sensed/detected by seismic sensors/receivers (hydrophones or geophones) at or near the surface of the earth.
Seismic data sets are processed based on a combination of signal processing and wave- equation methods to give images (seismic sections) of the subsurface. The aim of this processing is to be able to extract as much information as possible regarding the subterranean formations. Seismic images or sections represent slices through the geological subsurface, which can be interpreted manually or by employing workstations to decide if an area is a possible hydrocarbon prospect. Interpretation of seismic sections to locate the possible positions of hydrocarbon accumulations requires well-processed seismic data. Wells are drilled to verify the presence of petroleum. Since wells are extremely expensive and time consuming, there is always a continued need to improve the processing and display of seismic images.
In marine seismic surveying sea surface reflections (ghosts) pose a major problem in processing seismic data. Though various techniques have been employed to tackle this problem none is very efficient. 3-D modeling of sea surface multiples needs the knowledge of the sea surface profile and reflection coefficient values. In the deghosting of seismic data
Introduction
2
using the wavefield separation principle, the sea surface profile is needed in building vertical particle velocity wavefields from pressure wavefields. Similarly, if the vertical particle velocity wavefield is acquired simultaneously with the pressure wavefield, the sea surface profile is needed in applying the low frequency condition. Seismic data acquired in rough weather conditions with significant wave heights (SWH) of above 4 m are normally rejected;
the sea surface profile during the acquisition if available can serve as a check. 4-D seismic interpretation will be improved if the sea surface profile during acquisition can be obtained.
Dual-sensor measurements have the potential to image both the sea surface variation and to estimate the reflection coefficients along the streamer.
The main thrust of this work is to combat the detrimental effects of this sea surface reflection on seismic sections and possibly eliminate it by imaging the sea surface. The sea surface profile and reflection coefficient estimates are obtained by extrapolating the wavefield to the sea surface and employing a normalised imaging condition. The sea surface profile and reflection coefficients obtained can be important input parameters for further processing of seismic data as opposed to the conventional assumption of a flat sea surface and a reflection coefficient value of -1 at this interface.
3
Chapter 1
The marine seismic ghost
1.1 Introduction
In marine seismic surveying the source and receiver arrays are placed at certain depths below the sea surface because of technical reasons. The water surface is responsible for distorting both the source waveform as well as the data recorded by the receiver. Ghosting represents reflection of up-going energy at the water surface caused by the large impedance contrast between water and air (Ghosh, 2000). Thus, we have both source and receiver ghosts. The effective source signal includes not only the direct pulse but also its ghost. A similar effect arises on the receiver side as well. Therefore, in marine seismic acquisition a receiver measures a total wavefield that is distorted due to both the source and receiver side ghosts.
1.2 Marine seismic receiver ghost
In traditional marine seismic explorations vessels towing arrays of airguns (source) and streamers (receivers) are employed. Streamers containing pressure-sensitive hydrophones are towed at certain depths under the sea surface. The fundamental reason for this is that the free surface of water is assumed to be pressure free, thus excess pressure cannot be measured there. Practically, this is also done to avoid noise generated at the water-air interface due to sea waves (Ghosh, 2000). In a marine seismic survey, seismic signals travel downward through a body of water overlying the subsurface before being partly reflected upwards from the various subsurface structures. In a calm sea, the water-air contact almost acts like a perfect reflector whose reflection coefficient is very close to -1.
Chapter 1: The marine seismic ghost
4
Consequently, a hydrophone detects not only an up-coming reflection but also the corresponding negative reflection associated with the water-air interface (cf. Fig. 1.1). The sea surface related ghost effect is an obstacle that is probably as old as marine seismic exploration. Many authors have discussed this effect among them are Giles and Johnson (1973), Brandsaeter et al. (1979), Lindsey (1960) and Waters (1987).
Figure 1.1: A sketch depicting the down-going wave or ghost (D), up-going wave (U) and the receiver (ZR).
In Figure 1.1 θ is the angle of incidence of the up-going wavefield at the receiver (ZR). The amount of energy reflected depends on the condition/state of the sea surface. Rough sea surface conditions gives lower reflectivity coefficient compared to a calm sea (Gelius, 2004).
1.3 Relationship between hydrophone and geophone receiver ghosts
A hydrophone is a scalar sensor (non-directional) and measures pressure data corresponding to a combined pulse s(t). Assuming a reflection coefficient of -1 at the sea surface and introducing a ghost filter
!
g(t) we can express the combined pulse s(t) as follows:
!
s t
( )
=so( )
t "so(
t"#)
!
=so
( )
t "g(t) ,!
" =2d
co (1.1)
5 Where
!
so
( )
t is the primary or the up-going wave,!
so
(
t"#)
is the down-going wave (ghost),!
" is the delay time of the ghost, d is the depth at which the streamer is towed,
!
co is the velocity of sound waves in water and * indicates convolution (Ghosh, 2000). Only for reasons of simplicity we assume vertically incident signal in this review (i.e. θ = 0 in Fig.
1.1). After Fourier transform we can represent the hydrophone ghost filter
!
Gh(") in the frequency domain as:
!
Gh(") = S(")
So(")=1#e#i"$ (1.2) Similarly, it can be shown that the geophone sensor, which is a directional sensor, is associated with a ghost filter
!
Gg(") that can be written as:
!
Gg(") =1+e#i"$, (1.3) where the positive z-axis is measured downwards. The amplitude spectrum of the hydrophone ghost filter can be written explicitly as:
!
Gh(") =2Sin "d
co
#
$ % &
' ( (1.4) Whereas the amplitude spectrum of the geophone ghost filter can be written as:
!
Gg(") =2Cos "d co
#
$ % &
' ( (1.5) The amplitude spectra of these ghost filters show that the geophone has periodic notches at angular frequencies, ωg (cf. Fig. 1.2), while the hydrophone’s notches fall at angular frequencies, ωh (cf. Fig. 1.3), where:
!
"h = k#co
d , k = 0, 1, 2, … , and
!
"g = k#co
d , k = 0.5, 1.5, 2.5, … (1.6)
If we assume that the primaries and the ghost pulses of both geophone and hydrophone can be represented by unity spikes; adding the measurements of these two sensors will give a down-going unity spike (ghost) multiplied by two with a flat amplitude spectrum (cf. Fig.
Chapter 1: The marine seismic ghost
6
1.4). Similarly subtracting the geophone measurement from that of the hydrophone gives an up-going unity spike (primary pulse) multiplied by two (cf. Fig. 1.5).
Figure 1.2: The up- and down-going unity spikes of the geophone (left) and the amplitude spectrum (right).
Figure 1.3: The up- and down-going spikes unity of the hydrophone (left) and the amplitude spectrum (right).
Figure 1.4:Down-going spike (unity x2) (left) and the amplitude spectrum (right).
7
Figure 1.5: Up-going spike (unity x2) (left) and the amplitude spectrum (right).
From the above we can immediately see that the hydrophone and geophone ghost filter amplitude spectra are complementary (cf. Fig. 1.6). This has serious implications, which we will exploit further in this thesis in order to find efficient techniques to eliminate ghost notches from a seismic data set.
Figure 1.6: A superimposition of receiver ghost amplitude spectrum of the hydrophone (blue) on that of the geophone (red) for a collocated dual-sensor at a 15 m depth. In both cases the wavefield is assumed to have vertical propagation (Carlson et al., 2007).
1.4 Detrimental effects of the marine seismic receiver ghost
The amplitude spectrum of the receiver ghost shows that the effective bandwidth is between zero and the first notch. Since
!
"=2#f,
!
f being the frequency of the signal, the depth (d) of the sensor is given as:
Chapter 1: The marine seismic ghost
8
!
d="co
#1
= co 2fi,
!
i=h,g (1.7) Where
!
fg and
!
fh are the frequencies of the first notches of the geophone and the hydrophone respectively.
These notches place a limit in depth to which a conventional streamer (streamer with hydrophones only) can be towed thus, the frequency spectrum is band limited. Therefore in order to obtain lower frequency content, streamers are towed deeper and higher frequencies are attenuated and vice versa (cf. Fig. 1.7). It then follows that receiver ghosts limit seismic resolution. Emphasis on high resolution requires a bandwidth of about 250 Hz (Pieuchot, 1984). Thus on one hand, the requirement of high frequencies for improving resolution and the consequent need to push the zeroes of the ghost filter beyond Nyquist range enforce a shallow towing depth for the streamers. On the other hand, to minimize sea wave-generated noise a reasonably large towing depth is required. These conflicting demands would be reconciled satisfactorily if we could retain a large depth for the streamers without compromising the requirement of a large bandwidth. An ideal situation is to tow the streamer above the water surface where we eliminate all ghosts and obtain a large bandwidth. This can of course not be realised in practice but a towed dual-sensor streamer can provide a fair solution.
Figure 1.7: A superposition of receiver ghost amplitude spectra for a pressure sensor towed at 8 m (black) and 15 m (blue) depth, vertical incidence case is assumed here (Carlson et al., 2007).
From Eqs. (1.2) and (1.4) it follows that the amplitude spectrum of the combined pulse in
9 case of a hydrophone sensor can be written as:
!
S(") = Gh(")• So(") (1.8) A sketch of the amplitude spectrum of the combined pulse S(ω) is shown in Figure 1.8. The main criterion is that the depth d should be chosen as
!
d"#co
$ , Ω being the maximum angular frequency of the original signal
!
So("). In such a case the spectrum of the combined pulse S(ω) is not distorted by the notches. Figure 1.8 shows an example of a choice of a too large d.
Figure 1.8: A sketch showing the multiplication of the amplitude spectrum of a primary signal,So(ω), with the amplitude spectrum of the ghost, G(ω), to obtain the amplitude spectrum of a combined signal , S(ω). This is an example of a too large towing depth d.
The deghost filter in the case of hydrophone data can be written formally as:
!
Fh(")= 1
Gh(") = 1
(1+re#i"$) (1.9)
Chapter 1: The marine seismic ghost
10
Where r is the reflection coefficient. In the various seismic deghosting techniques, the reflection coefficient r is often assumed to be -1. This assumption introduces errors in the interpretation of the resulting seismic image. The effect of rough sea surface condition on seismic images can be quite dramatic. An overlap between the desired data and its ghost can combine destructively to mask subsurface information and compound interpretation (Schneider et al, 1964). Sea surface reflection coefficients also vary along the streamer. If we employ the correct values of r in the processing of seismic data we will obtain better images of the subsurface. Thus, this work seeks to image the sea surface and possibly obtain reflection coefficient estimates along the streamer.
In 4-D seismic acquisition one factor affecting data repeatability is the influence of rough sea surface reflection, since the sea surface is assumed to be flat and stationary in processing conventional seismic data. The inclusion of rough sea impulses confuses interpretation of the time-lapse images by introducing an error into the signal, giving rise to spurious events in the difference image that ‘mimic’ the real structure in the section (Laws and Kragh, 2002a).
Based on the detrimental effects of ghost on seismic images, the need for a different method of signal processing cannot be over emphasized.
1.5 Deghosting approaches
Traditional deghosting of seismic data is not a trivial procedure because of the zeroes in the spectrum of the ghost filter (Ghosh, 2000). The actual sea surface profile is an important input parameter of a rough-sea deghosting solution. Kragh et al., (2002b) show how the sea surface wave height can be measured using very low-frequency pressure fluctuations (below 0.5 Hz) from the hydrophones. Amundsen et al., (2005) extended this approach, by using the derived sea surface depth variation to obtain the up-going pressure wavefield (deghosted seismic data) using the principle of wavefield separation. This approach has a limitation of employing a pressure gradient approximation. Again, the truncation involved in the calculation of the vertical particle velocity wavefields from the pressure wavefield measurements further limits the accuracy of this method. Kragh and Laws (2001a) applied a
11
statistical deconvolution method to remove the rough sea error. Important limitations of this approach are the assumptions of whitening and wavelet stationarity.
Ziolkowski (1971) shows how a number of airguns at different depths can be used to eliminate the zeroes in the amplitude spectrum of the source ghost making it possible for eventual deconvolution of the source signature. Ghosh (2000) used a similar approach to show how zeroes in the spectrum of recorded data caused by receiver-side ghosts can be eliminated. An important criterion pointed out by Ghosh (2000) is that the recordings should be carried out at two different depths chosen in such a way that the respective receiver ghost filters do not share common zeroes. A limitation of this approach is that vertical incidence is assumed and accurate knowledge of the reflection coefficients of the sea surface is still needed.
In this work we deghost marine seismic data by exploiting the relationship between hydrophone and geophone ghost filters (i.e. by collocating the two sensors). One of the sensors measures pressure wavefield (hydrophone) while the second sensor measures vertical particle velocity wavefield (geophone). By applying the principle of wavefield separation, the recorded seismic data is decomposed into up- and down-going particle velocity wavefields and up- and down-going pressure wavefields. Any of the up-going wavefields is an equivalent of a conventional deghosted seismic data, but generally with a better resolution (Carlson et al., 2007). Correspondingly, the down-going wavefields represent the ghost wavefields.
1.6 Dual sensor streamer (Geostreamer)
If the vertical component of the particle velocity wavefield is measured simultaneously with the pressure wavefield (collocated dual-sensor measurement), the combined seismic wavefield can be decomposed into up- and down-going pressure wavefields as well as up- and down-going vertical velocity wavefields (Schneider and Backus, 1964; Claerbout, 1976;
Barr and Sanders, 1989; Amundsen, 1993; Fokkema and van den Berg, 1993).
Chapter 1: The marine seismic ghost
12
Based on this principle a new generation of towed streamers comprising of collocated hydrophones and geophones (vertical particle velocity sensors), which respectively measure pressure wavefields and vertical particle velocity wavefields have been developed. This Geostreamer basically exploits the complementary relationship between these two sensors as demonstrated in Figs. 1.2-1.5. Figure 1.9 shows stacked seismic sections based on both geophone and hydrophone measurements. Notice the flat amplitude spectrum of the deghosted seismic data.
Figure 1.9: Unmigrated stack comparison. The images are from left to right, the pressure wavefield, the vertical velocity wavefield, and the up-going (deghosted) pressure wavefield derived from the pressure and velocity wavefields. Note the complementary amplitude spectra for pressure vs. particle velocity (Carlson et al., 2007).
A combination of these sensors allows a wider streamer towing depth range. Advantages of using dual sensor streamers in seismic acquisition are summarized by (Carlson et al., 2007).
A pertinent fact is that dual sensor streamer data generally contains a greater frequency bandwidth, which gives the related seismic images better resolution. With the dual sensor streamers we can tow deeper where weather and operational noises are minimal and the streamers are better behaved. The improved signal-to-noise ratio of the related seismic image is more reliable to interpret. Thus, deeper target amplitudes are stronger and more continuous. The down-going wavefield parts of the sensors (representing the sea surface ghosts) can relax the sea surface assumptions of SRME (sea surface related multiple elimination) if included in the prediction of multiples (Söllner et al., 2007, Frijlink et al.,
13
2009). Furthermore, the up-going and down-going wavefields can be extrapolated and summed to reconstruct the total wavefield at any desired recording depth to duplicate the parameters of any existing survey, allowing 4-D matching plus the benefits of improved image clarity (cf. Fig. 1.10). A limitation of the dual sensor data is that the geophone data is swamped by noise from 0-20 Hz. Any existing 4-D seismic section can be obtained from the dual sensor data in the following steps (Carlson et al., 2007):
Measure data at a given depth using dual sensor streamers
Apply wavefield separation at the acquisition depth
Build the geophone low-frequency (0-20 Hz) data from that of the hydrophone data (applying low-frequency conditioning)
Forward propagate the up-going wavefield to the desired depth (f.ex 8 m)
Backward propagate the down-going wavefield to the desired depth (f.ex 8 m)
Add the up- and down-going wavefields at the new depth
Figure 1.10: The images are from left to right, the hydrophone (pressure) stack acquired at 8 m depth (blue curve at the right), the equivalent total pressure stack at a stimulated depth of 8 m derived from the dual sensor data acquired at 15 m depth (red curve at right) and the comparative amplitude spectra (Carlson et al., 2007).
14
Chapter 2
Sea surface imaging technique
Seismic prospecting consists of generating very low–amplitude artificial earthquakes at predetermined times and positions. The induced seismic disturbances are recorded by receiver spreads making up the traces. Acquisition geometry is defined by the distribution of the source and receiver group. Seismic energy radiated by the source is split between body waves (compressional or shear waves) and surface waves (Raleigh and Love waves). This energy is usually partly reflected by subsurface reflectors (geologic layer boundaries/interfaces usually having different acoustic impedances). Reflection coefficient is the amplitude ratio of reflected energy to incident energy at material interfaces. In petroleum exploration we are mostly concerned with body waves. Seismic medium/layers can be considered as being elastic and obey Hooke’s law. The displacements that are observable at all points within the medium and, particularly, on the surface are solutions to wave equation.
In this thesis, we will consider this displacement vector as a 2-D scalar potential that gives rise to the P-wave equation. P-waves are compressional waves corresponding to longitudinal vibrations, which at every point in the medium have a particle motion parallel to the direction of propagation. The hydrophones employed here measure the pressure difference between the ambient motion and that caused by the source in the direction of propagation of the wave while the geophones measure the vertical particle velocity of the propagating wave.
The imaging technique as will be applied here consists of three main parts namely: wavefield separation, extrapolation and the imaging condition. Wavefield separation is applied to the dual sensor data at the acquisition depth to obtain four wavefields namely; up- and down- going pressure wavefields for hydrophone and up- and down-going vertical particle velocity wavefields for the geophone. The up- and down-going wavefields of any of the sensor can then be extrapolated to the sea surface where an imaging condition is applied.
15
2.1 Wavefield separation
The Geostreamer equipped with collocated dual sensors (geophones and hydrophones) represents a new generation of streamers developed by PGS. This section will demonstrate that such dual-sensor recording can be used to decompose the measurements into upward and downward travelling wavefields. The principle is fundamentally the same as that demonstrated earlier with unity spikes (cf. Figs. 1.2-1.5). In Figure 2.1, the yellow tube depicts the streamer (Geostreamer) with the collocated sensor mounted inside it; above this is a schematic of a rough sea surface. We can see that the geophone records the same polarity for both up- and down-going signals while the hydrophone records opposite polarity for up-going and down-going signals.
Figure 2.1: A schematic of the dual-sensor streamer (left), the up- and down-going wavefields of the hydrophone (blue) and that of the geophone (red). The amplitude measurements of the up- and down- going pulses as well as the total wavefields for the respective sensors are shown.
To demonstrate the working principle of a dual sensor streamer we first derive the relationship between pressure and particle velocity. Consider now the generalised Newton’s second law, which can be written as:
Orji.O.C
16
F = ma (2.1) Where F is force, m is mass and a is acceleration of the body. Consider a rectangular isotropic cube of fluid with zero viscosity, this cube has a base area A, height
!
z, mass
!
m, and density
!
" as shown in Figure 2.2.
Figure 2.2: A cube of non-viscous isotropic fluid with a net force ΔF (in the z-axis) acting on it in.
In order to compress this cube we need a net force ΔF. If we consider a net force in the z- axis only, the relationship between the acceleration,
!
dVz
dt , as a result of the pressure, ΔP, caused by this force is (Berkhout, 1982):
!
!
"#P
#z = "#F
#A#z =$dVz
dt (2.2) where t is time,
!
Vz is the vertical particle velocity and ΔP represents the net pressure. After Fourier transform of Eq. (2.2) and considering an infinitesimal volume element of the block we can relate the vertical particle velocity field and the pressure field as follows:
!
Vz = "1 i#$
dP dz ,
!
i" # d
dt (2.3)
17 If we now assume 2-D waves measured at a depth
!
z and introduce spatial Fourier transform of the horizontal coordinate so as to decompose the wavefield into monochromatic plane waves, the up-going plane wave can then be written formally as:
!
U=Aexp
[
i k(
xx+kzz+"t) ]
(2.4) while the down-going plane wave is:
!
D=Bexp
[
i k(
xx"kzz+#t) ]
(2.5)!
kx and
!
ky are horizontal wavenumbers along respectively x and y direction while
!
kz is the vertical wavenumber. The 2-D wavenumbers of a plane-wave fulfil the following dispersion relation:
!
kx2+kz2 = "
c
#
$ % &
' (
2
=k2 (2.6) Extension of the above formulation to 3-D is straightforward:
!
kx2+kz2+ky2 = "
c
#
$ % &
' (
2
=k2
In Eqs. (2.4) and (2.5), A and B are the amplitudes of the up-going and down-going plane waves respectively. Now, assuming a plane wave propagating with a given angle along the vertical direction (positive
!
z-axis is defined to be downwards), the hydrophone, P, will measure a total pressure wavefield of:
P =U+D (2.7) While the geophone will measure a total vertical velocity wavefield of:
!
Vz = "dp i#$dz= kz
$#
(
D"U)
=cos%(
D"U)
#c (2.8) where θ is the angle of incidence. Adding equations (2.7) and (2.8) will give:
2D
!
D=1
2 P+"#
kz Vz
$
% & '
( ) (2.9)
Orji.O.C
18 While subtracting (2.8) from (2.7) gives:
2U
!
U=1
2 P"#$
kz Vz
%
&
' (
) * (2.10) We can represent these equations in a matrix form as given below:
!
U D
"
# $ %
&
' =1 2
1 ()*
kz
1 )*
kz
"
#
$
$
$
$
%
&
' ' ' '
P Vz
"
# $ %
&
' (2.11)
Thus, by adding or subtracting the hydrophone and geophone measurements we can separate the up- and down-going wavefields. The principle of wavefield separation as derived above is illustrated schematically below (cf. Fig. 2.3). Note that U and D in the Eqs. (2.7-2.11) above denote respectively up- and down-going pressure wavefields. Divergent responses at small vertical wavenumbers (i.e. division by small values in Eq. (2.11)) cause problems for events with large incidence angles when the pressure and particle velocity data are combined in order to separate the wavefields as shown in Figure 2.4. This makes the expression unstable at high incidence angles (θ ≤ 90°). This is avoided by applying a dip filter. This filtering is carried out in ω-k domain prior to the wavefield separation. Wavefield separation constitutes the first step of imaging the sea surface.
-150 -100 -50 0 50 100 150
0 20 40 60 80 100
Time (msec)
Amplitude
-150 -100 -50 0 50 100 150
0 20 40 60 80 100
Time (msec)
Amplitude
- =
2 up -going P
pressure
-300 -200 -100 0 100 200 300
0 20 40 60 80 100
Time (msec)
Amplitude
-150 -100 -50 0 50 100 150
0 20 40 60 80 100
Time (msec)
Amplitude
-150 -100 -50 0 50 100 150
0 20 40 60 80 100
Time (msec)
Amplitude
+ =
2 down -going P
-300 -200 -100 0 100 200 300
0 20 40 60 80 100
Time (msec)
Amplitude
Total P Total G
€
=cosθρc Vz
Figure 2.3: Graphical illustration of wavefield separation principle. The “right” combination of the hydrophone (blue) and the geophone (red) amplitude measurements in time gives either the up-going or the down-going pressure wavefield.
19
Since the receiver ghost of the hydrophone has opposite polarity compared to that of the geophone receiver ghost (the geophone is directional), a combination of the two ghost functions suppresses receiver ghost as shown earlier. If the receiver ghost is suppressed then the streamer can be towed at a greater depth where it encounters reduced acquisition noise, therefore the seismic acquisition weather window is increased. An enhanced frequency bandwidth is obtained consequently increasing the resolution of the corresponding seismic image.
2.2 The low-frequency condition
In practice the velocity sensor will be swamped by strong mechanical noise over a typical frequency range of 0-20 Hz. For this frequency range, the velocity sensor data are predicted from the pressure sensor data by applying a low-frequency condition. The governing equation is derived below. The sea surface is associated with a reflection coefficient of r and the relationship between the amplitudes of the up- and down-going plane waves can be written as:
!
r= B
A, B = rA (2.12) The relationship between the scattered wavefields (i.e. muted direct wavefield) of vertical particle velocity sensor and pressure sensor follows from Eq. (2.8), i.e.
!
Vz = kz
"# rAe
$ikzz
$Aeikzz
( )
e(ikxx+i"t)!
= Akz
"# re
$ikzz
$e+ikzz
( )
e(ikxx+i"t) (2.13)From Eq. (2.7) we will also have:
!
P=
(
Ae+ikzz+rAe"ikzz)
e(ikxx+i#t)!
=A(eikzz+re"ikzz)e(ikxx+i#t) (2.14)
Combination of Eqs. (2.13) and (2.14) gives (2-D case):
!
Vz(kx,zR,")=
!
kzc
"
#
$ % &
' ( 1 )c
#
$ % &
' ( (re*ikzzR *eikzzR)
(re*ikzzR +eikzzR)P(kx,zR,") (2.15)
Orji.O.C
20
In case of a flat sea surface (i.e. r = -1), Eq. (2.15) takes the special form shown below, which is the low frequency condition:
!
Vz(kx,zR,")=
!
" kzc
#
$
% & ' ( ) 1
*c
$
% & '
( ) (e"ikzzR +eikzzR)
(eikzzR "e"ikzzR)P(kx,zR,#)=" kzc
#
$
% & ' ( ) 1
*c
$
% & '
( ) (1+e"i2kzzR)
(1"e"i2kzzR)P(kx,zR,#) (2.16) We define the following terms for Eq. (2.16):
the minus sign shows the opposite polarity relationship between the up-going wavefields of geophone and hydrophone
!
zR is the recording depth
!
"cis the acoustic impedance
!
1+e"i2kzzR is the geophone ghost function
!
1"e"i2kzzR is the hydrophone ghost function
!
P(kx,zR,") is the hydrophone pressure measurement
!
Vz(kx,zR,") is the geophone vertical velocity measurement
!
kz ="
c cos# is the vertical wavenumber (θ = angle of incidence)
The implication of the relationship between pressure and vertical velocity wavefields is that we can build the vertical velocity wavefield given the pressure wavefield data and vice versa.
As only the vertical component of particle velocity is reconstructed in Eq. (2.16), we multiply by cosine of the incidence angle given by (
!
kzc
" ) as shown in Figure 2.4. The low
frequency-conditioning equation, (Eq. (2.16)), implies a division of geophone ghost function by hydrophone ghost function and therefore is unstable at the notches (e.g at 0 Hz and 50 Hz for a receiver at 15 m depth). Hence, there is an upper limit in depth at which the streamer can be towed (Carlsen et. al., 2007).
21
Figure 2.4: A schematic representation of the relationship between the vertical wavenumber (
!
kz) and horizontal wavenumber (
!
kx).
2.3 Wavefield extrapolation
The second step in imaging the sea surface is to extrapolate in small depth increments the up-going pressure wavefields forward to the sea surface and the down-going pressure wavefields backward to the sea surface. We perform this extrapolation in the
!
" #kdomain for computational efficiency. Consider a flat straight streamer at a recording depth
!
z=zR below the sea surface (cf. Fig. 2.5). We introduce an arbitrary observation level,
!
zo, above this streamer for the purpose of extrapolation. Wavefield separation is performed at the streamer level and by extrapolating the resulting up- and down-going wavefields of a given sensor in a stepwise manner we move towards the sea surface. The related equations are shown below.
Figure 2.5: A sketch depicting a flat sea surface, an observation level (
!
zo) and a flat streamer at a recording depth (
!
zR) for the purpose of extrapolation.
Orji.O.C
22
In order to extrapolate the plane wave to the sea surface, we start with the scalar wave equation in 2-D. Let
!
p x,z,t
( )
be a pressure wavefield propagating with a velocity v in the!
x,z, plane and satisfies the scalar wave equation, Eq. (2.17).
!
"2p
"x2 +"2p
"z2 = 1 v2
"2p
"t2 (2.17) The wave equation above can be solved using the 2-D inverse Fourier transform
!
P k
(
x,z,")
of the
!
p x,z,t
( )
(Buttkus, 2000):
!
p x,
(
z,t)
="#
#
$
P k(
x,z,%)
"#
#
$
exp[
i k(
xx+%t) ]
dkxd% (2.18)The partial derivatives in Eq. (2.17) are then:
!
"2p
"x2
"2p
"t2
"2p
"z2
#
$
%
% %
&
%
%
%
=
'kx2P k
(
x,z,()
'(2P k
(
x,z,()
d2P k
(
x,z,()
dz2 )
*
%
% %
+
%
% %
#
$
%
% %
&
%
% %
',
-
, ',-
, exp[
i k(
xx+(t) ]
dkxd( (2.19)Thus the wave equation can be transformed into a normal differential equation for the Fourier transform
!
P(kx,z,"):
!
d2P dz2 + "2
v2 #kx2
$
% & '
( ) P=0 (2.20)
The general solution of the differential Equation above, (Eq. (2.20)), is:
!
P=Aexp
[
ikzz]
+Bexp[
"ikzz]
, (2.21)where the vertical wavenumber is given as
!
kz = "2 v2 #kx2
$
% & ' ( )
1 2
. The first term and second term on the right hand side of Eq. (2.21) describes up- and down-going wavefields respectively.
Now for a wavefield propagating upwards, only the first term applies thus, we set B = 0. It then follows from Eq. (2.18) that:
23
!
p x,
(
z,t)
="#
#
$
"##
$
Aexp[
i k(
xx+%t+kzz) ]
dkxd% (2.22)A comparison of Eq. (2.22) and Eq. (2.18) shows that for the 2-D Fourier transform of the up-going wavefield at a depth
!
z,
!
P k
(
x,z,")
=Aexp[
ikzz]
(2.23) Let now!
P k
(
x,zR,")
represent the scattered wavefield measured by the streamer at the acquisition depth!
z=zR, Eq. (2.23) can then be rewritten as:
!
P k
(
x,zR,")
=Aexp[
ikzzR]
(2.24) Similarly, at the observation or datum level!
z=zo we will obtain:
!
P k
(
x,zo,")
=Aexp[
ikzzo]
(2.25) It then follows that we can extrapolate the measured wavefield at!
z=zR to the observation level
!
z=zoby combining Eqs. (2.24) and (2.25) :
!
P k
(
x,zo,")
=P k(
x,zR,")
exp[
#ikz(
zR #zo) ]
(2.26) Equation (2.26) describes the forward extrapolation of the wavefield from the recording depth to a chosen observation level. The backward extrapolation can be obtained by assuming a down-going wavefield and setting A = 0 in Eq. (2.21). If we now apply our convention of denoting the up-going pressure wavefield as U and the down-going pressure wavefield as D; then the up-going wavefield can be forward propagated to the observation level,!
zo, by applying:
!
U(kx,zo,")=U(kx,zR,")exp
[
#ikz(zR #zo)]
(2.27) While the downward travelling wave is backward propagated to the observation level,!
zo, by applying:
!
D(kx,zo,")=D(kx,zR,")exp
[
#ikz(zo#zR)]
(2.28)Orji.O.C
24
Similar expressions can be written for the particle velocity sensor. Note that in the above equations we consider only scattered wavefields and assume that direct arrivals have been muted. We choose a flat sea surface and streamer here for simplicity. For a non-flat streamer, redatuming can be applied to condition the data as described in Appendix A. If the streamer is smoothly shaped but not flat and the data is acquired under a rough sea surface condition, we can choose the observation level to coincide with the datum level by employing Eq. (2A) (see Appendix A). Since the depth of the sensors varies with distance along the streamer because of the streamer weight and other factors, a variable depth wavefield separation that takes this effect into account is applied in reality. This is achieved by employing a processing software that reads the depth information for each group of sensors. The main purpose of extrapolation is to be able to image the sea surface and possibly recover the reflection coefficient at the water-air interface along the streamer. Hence, it is pertinent to choose an optimal extrapolation step-size in order to capture the sea surface wave height variations (typically 4 m SWH for a rough sea day (Kragh et al., 2002a)). SWH is the significant wave height defined in physical oceanography as the average wave height (trough to crest) of one-third of the largest waves.
2.4 Imaging condition application
Detrimental effects of sea surface reflection on seismic images have been demonstrated. The conventional assumption of a reflection coefficient of -1 at the sea surface can be overcome by employing the dual-sensor technique. In this study we seek to image the sea surface and show that it is indeed not flat and also varies with time. In the third step of the sea surface imaging we apply an imaging condition. After decomposing the dual-sensor data, the separated up- and down-going wavefields are extrapolated iteratively to the sea surface where the imaging condition is applied. The normalized imaging condition employed here is (Claerbout, 1971, 1976, Kaelin and Guitton, 2006, Schleicher et al, 2008):
I(x,z)
!
=
U x,z,"
( )
D x,z,"( )
"
#
U x,z,"
( )
U x,z,"( )
"
#
(2.29)In a chosen time window of traces/channels, cross-correlation of the extrapolated up-and down-going wavefields is applied (numerator of Eq. (2.29)). An autocorrelation of the up-
25
going wavefield is applied in the same manner (denominator of Eq. (2.29)). The correlations are carried out in the frequency domain; hence time-reversal is replaced by complex conjugating (denoted by the bar).
The zero-lag values of the cross-correlation are divided by the zero-lag values of the autocorrelation for each extrapolation step. The extrapolation depth corresponding to the maximum image value I(x,z) for a given channel corresponds to the sea surface depth (z) at the channel’s position (x). A plot of these extrapolation depths versus the channels’ positions gives a time stationary image of the sea surface along the streamer corresponding to the selected time window. For a flat sea surface the maximum image value will occur at zero- lag, which will always correspond to extrapolation to zero meter depth or the reference/datum. In the case of a rough sea surface the maximum image value will also occur at the zero-lag, which is not necessarily at the reference/datum. We expect that at the sea surface, the up- and down-going wavefields should be exactly the same except for a180 degrees phase shift.
At a given depth below the sea surface there is a time delay between the up- and the down- going wavefields. If we extrapolate towards the sea surface in a stepwise manner with a given step-size, this delay time gets smaller. At sea surface, the difference between the arrival times of up-going and down-going wavefields is zero (cf. Fig. 2.6). We can view this extrapolation as bringing the streamer systematically closer and closer to the sea surface.
Thus, when the streamer is exactly at the sea surface, the up- and down-going wavefields will be recorded exactly at the same time. In Figure 2.6 the sea surface is assumed flat for simplicity. In this flat case if we add the amplitudes of up- and down-going wavefields, they will cancel each other out. This is because the down-going wavefield is the same as the up- going wavefield, except for a phase shift of 180 degrees caused by the sea surface reflection.
Orji.O.C
26
Figure 2.6: A sketch showing the stepwise extrapolation of up-going and down-going wavefields to the sea surface.
If the sea-surface was time-invariant, then one imaging step would apply to the entire recording time. As the sea-surface is changing in time, one imaging condition applies solely in the time window of the extrapolated wavefields. From a sliding time window of the correlation analysis iteratively applied along the length of the traces, we can extract the time- variant behaviour of one image point. By a superposition of all the image points of a particular record, the sea surface state of the related acquisition window is obtained (PGS- 08-04US).
PGS proprietary processing software exists which can be used to image the sea surface in this way. In this work a flat sea surface will be assumed in generating synthetic seismic data based on dynamic ray tracing and another set of synthetic data will be generated based on the finite difference method assuming various sea surface conditions. The fundamental difference between these two methods is that ray tracing is based on a high frequency
27
approximation (i.e. assuming wavelengths far smaller than the velocity anomalies in the model). The described imaging technique will be applied to these synthetic data to obtain the sea surface image and possibly also the reflection coefficients along the streamer thus, demonstrating the validity of the procedure. One basic assumption in the synthetic data studies is that the sea surface is assumed to be stationary. Finally, we will apply this technique on a 2-D field data acquired from the Norwegian North Sea where the weather condition is often marginal. A description of the processing flow applied to both the synthetic and field data is given in Appendix B.
28
Chapter 3
Synthetic data examples
Both ray tracing and finite difference methods were employed to generate synthetic data sets used in investigating the imaging technique. In the ray tracing case because of the inherent high frequency assumption, a flat sea surface was assumed when generating the synthetic data. In the finite difference case, synthetic data sets were generated for three cases: flat sea surface, simple sea surface and rough sea surface conditions. The imaged sea surface conditions were compared with the modeled sea states.
3.1 Imaging based on test data from ray tracing
3.1.1 Definition of controlled model and data generation
A basic summary of the dynamic ray tracing (DRT) technique is given in Appendix C. The PGS proprietary software, Nucleus, was used to generate synthetic data for both hydrophones and geophones based on DRT. First, a vessel towing a single 3000 m long streamer at a distance of 250 m from the vessel in the inline direction (x-axis in this case) and a depth of 15 m is modeled. The vessel is also towing a source at a distance of 50 m and depth of 6 m. The streamer is a Geostreamer containing groups of 16 collocated pressure and vertical particle velocity sensors at 12.5 m group intervals forming the acquisition surface assumed to be flat in this case (cf. Fig. 3.1). The shot was fired at intervals of 25 m. A recording length of 4 seconds and a sampling interval of 2 milliseconds were used.
29
Figure 3.1: A plot of the vessel towing the streamer and the source.
The type of source used here is a single airgun with the time-function plot of the signature given in Figure 3.2. We can see that this signal is 400 milliseconds long. Though, this pulse is not an ideal source signature, it has a high primary to bubble ratio. We will show later that shorter pulses give better resolution of the reflectors.
Figure 3.2: A plot of the time-function of the source signature.
The corresponding earth model is horizontally layered containing 4 interfaces. The first interface at 0 m is the air-water interface (sea surface); the second interface at a depth of 1000 m is the sea floor while the remaining two interfaces respectively at depths of 2000 m
Orji.O.C
30
and 3000 m are geological reflectors. The last layer in Figure 3.3 is a half space. The layers were respectively assigned P-wave velocities and densities of: VP = 333 m/s, ρ = 0.01 g/cm3 (air), VP = 1500 m/s, ρ = 1 g/cm3 (water), VP = 2500 m/s, ρ = 1 g/cm3 (first sediment), VP = 3000 m/s, ρ = 1.5 g/cm3 (second sediment) and VP = 2500 m/s, ρ = 2 g/cm3 (half space). The total size of this 2-D model is 10,000 m in lateral direction and 3500 m in the vertical direction (cf. Fig. 3.3).
Figure 3.3: The 2-D horizontal plane-layer model showing the P-wave velocities.
A survey area which may be a point survey (if only a point is to be surveyed), a line survey (in which case the survey is carried out along a line) or an area survey (survey within an area) must be defined in order to acquire synthetic data set. In this case we defined a line survey that is 1000 m long in order to keep the model simple (cf. Fig. 3.4). The coordinates of the survey were defined based on the UTM coordinate system and events were traced based on the ray tracing method. The generated hydrophone and geophone synthetic data sets were in SEG-Y format and without any form of noise or NMO correction. The generated traces were displayed in trace (time-channel) plots.
