Surface-wave separation and its impact on seismic survey design

(1)

Surface-wave separation and its

impact on seismic survey design

(2)

(3)

Surface-wave separation and its

impact on seismic survey design

PROEFSCHRIFT

ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus Prof.ir. K.C.A.M. Luyben, voorzitter van het College voor Promoties,

in het openbaar te verdedigen op donderdag 9 april 2015 om 10.00 uur door Tomohide ISHIYAMA

Master of Earth Science geboren te Tokyo.

(4)

Prof.dr. W.A. Mulder

Copromotor:

Dr.ir. G. Blacqui`ere

Samenstelling promotiecommissie:

Rector Magnificus, Technische Universiteit Delft, voorzitter Prof.dr. W.A. Mulder, Technische Universiteit Delft, promotor Dr.ir. G. Blacqui`ere, Technische Universiteit Delft, copromotor Prof.dr.ir. A. Gisolf, Technische Universiteit Delft

Prof.dr. S.M. Luthi, Technische Universiteit Delft Dr.ir. D.J. Verschuur, Technische Universiteit Delft Prof. K.A. Berteussen, The Petroleum Institute Dr.ir. A.W.F. Volker, TNO

This research was a project of ADNOC R&D Oil Sub-Committee (R&D OSC) and the DELPHI consortium, and was financially supported by IN-PEX, ADNOC R&D OSC and the DELPHI consortium.

ISBN 978-94-6295-102-0

Copyright c⃝ 2015 by T. Ishiyama.

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.

Published by: Uitgeverij BOXPress

Printed by: Proefschriftmaken.nl, Uitgeverij BOXPress

(5)

(6)

(7)

Chapter 1 Introduction

1.1 Surface waves in seismic data

A seismic survey is an experiment designed to map subsurface structures and properties by recording seismic waves, often with the aim of locating reser-voirs and deposits of oil and gas (Sheriﬀ, 2002). Usually seismic waves are generated with artificial seismic sources. The recorded seismic data contain the arrivals of events such as surface waves, refractions, reflections, surface-related and internal multiples, etc.

Surface waves are often the most dominant events in the observed seismic data, in land and shallow-water environments. They distinguish themselves from the primaries by larger amplitudes, lower frequencies and lower appar-ent velocities. Traditionally, they were considered as coherappar-ent noise, masking primaries and hampering reservoir imaging and characterization. Nowadays, they are more and more regarded as useful signal for near-surface charac-terization. In either case, their separation from the primaries is of great importance. However, their complex properties, such as dispersion, multi-modality and spatial variability, make the surface-wave separation quite a challenge in processing.

A seismic survey is identified by the acquisition geometry, which in turn is specified by the survey parameters. Therefore, the observed seismic data are not only a result of the Earth’s response, described by its geophysical parameters, but also of the chosen survey parameters. The survey parameters should be chosen such that the acquired data have the quality required to achieve the objectives, e.g., for exploration, appraisal and development of oil and gas fields. Since the surface waves are often dominant, they impose certain requirements on the survey parameters, because these should allow for eﬀective surface-wave separation. Here, the word ‘separation’ rather than

(10)

Figure 1.1: Seismic data in a shallow-water environment, (a) and (d) with a spatial

sampling interval of 25 m, (b) and (e) with a spatial sampling interval of 50 m, (c) and (f) with blending, in the xy-t and the kxky-f domain with a vertical section at the top and a

horizontal time/frequency slice at the bottom. A dotted line in the section indicates the position of the slice, and vice versa. (to be continued)

‘removal’ is used to emphasize the potential application of surface waves for near-surface characterization.

Figure 1.1 shows typical shot gathers containing surface waves in the space-time (xy-t) and the wavenumber-frequency (kxky-f ) domain, just to

illustrate the appearance of surface waves in seismic data. The surface-wave characteristics and properties as well as their separation from the primaries will be discussed in detail in the following chapters, together with an over-view of earlier work. In this example, dispersive surface waves stand out clearly with their larger amplitudes and lower frequencies. Two modes, the fundamental mode and one higher mode, can be identified by their slow-nesses being larger than those of the primaries. Here, slowness is reciprocal of velocity. In addition, the surface waves are aliased in the case of a large spatial sampling interval in the acquisition (Figure 1.1(b) and 1.1(e)). Fur-thermore, in the case of blended acquisition, they are blended just as the primaries (Figure 1.1(c) and 1.1(f)). This example shows that surface waves have a diﬀerent appearance in various data sets for diﬀerent survey

(11)

paramet-1.2. Objectives and outline 3

Figure 1.1: (Continued) Red, magenta, green and blue arrows indicate the surface-wave

fundamental mode, a surface-wave higher mode, refractions and reflections, respectively. Filled arrows indicate the true energy, and the whitish versions of the arrows indicate the aliased energy. Notice that the surface waves are dispersive, multi-modal, aliased, and blended for (c) and (f). For the blended case, the arrows are for the various events from one source with a zero delay time in the xy-t domain, but from all three sources in the

kxky-f domain.

ers. Therefore, different data sets specified by different survey parameters may lead to differences in data quality in terms of surface-wave separation.

1.2 Objectives and outline

Objectives of this thesis

The objectives of this thesis, motivated by the above challenges and problems, are

• To establish a method of 3-D surface-wave estimation and separation,

which addresses the surface-wave properties and adopts recent advances in processing.

• To understand how the application of surface-wave separation aﬀects

(12)

data quality.

Outline of this thesis

The remainder of this thesis consists of six chapters, each self-contained with their own introduction and literature review.

In chapter 2, a forward model of surface waves in seismic data is de-scribed. First, we review the surface-wave properties in seismic data. Then, we introduce a forward model of primaries, which takes the survey specific-ation into account, and extend the representspecific-ation of the forward model to include surface waves. The definition of survey parameters and the paramet-erization of surface-waves introduced in this chapter have an important role in the thesis.

In chapter 3, the method of 3-D surface-wave estimation and separa-tion is established. There, we try to achieve the first objective of the thesis. First, we discuss the inverse problem for the above forward model of sur-face waves. Then, we introduce a method using a closed-loop approach to solve the inverse problem. We review the literature on conventional and re-cent approaches of 3-D surface-wave estimation and separation with their features and limitations. We also discuss recent advances in processing, such as a closed-loop or iterative approach. Finally, the established method is compared to a conventional method.

In chapter 4, we apply the established method to several data sets to demonstrate its virtues. These include 3-D geophone data sets in a land en-vironment simulated by a 3-D elastic finite-diﬀerence code, real 3-D geophone data sets acquired in a land environment, as well as real 3-D ocean-bottom-cable (OBC) hydrophone data sets acquired in a shallow-water environment. We also apply the method to under-sampled, asymmetrically sampled, irreg-ularly sampled and blended data sets to demonstrate the wide range of its applicability.

In chapter 5, we discuss the relationship between the survey parameters and the resulting data quality in order to find the essential types of survey parameters and their optimal values for a required data quality in the context of surface-wave separation. In this chapter, we try to achieve the second objective of the thesis. We carry out a case study and examine the impact on survey design in terms of surface-wave separation. The latter application should be included in survey design, in addition to the intended reflection imaging and amplitude-versus-oﬀset (AVO) applications.

In chapter 6, we further discuss the relationship between the survey parameters and the resulting data quality, again in the context of reflection

(13)

1.2. Objectives and outline 5

imaging and AVO applications. We carry out a case study for this purpose and further describe how the survey parameters aﬀect the resulting data quality.

Chapter 7 summarizes the main conclusions and includes

recommenda-tions for future research.

It should be noted that chapters 2 to 6 consist of a series of forthcoming papers, Ishiyama et al. (2014a), Ishiyama et al. (2014d) and Ishiyama & Blacqui`ere (2014).

(14)

(15)

Chapter 2 Forward model of surface waves in

seismic data

1

2.1 Introduction

Seismic data in land and shallow-water environments are often dominated by surface waves. Traditionally, they were thought of as coherent noise masking primaries and hampering reservoir imaging and characterization. Nowadays, they are more and more regarded as useful signal for near-surface character-ization. In either case, their separation from primaries is of importance.

Any events in discretely sampled seismic data, such as surface waves and primaries, are acquired by a seismic survey. The seismic survey is identified by the acquisition geometry, and the acquisition geometry is specified by the survey parameters. Therefore, the observed seismic data result not only from the Earth’s response grounded on its geophysical parameters but also from the survey specification described by its parameters. A forward model of seismic data should take this fact into account.

The goal of this chapter is to obtain a forward model of surface waves in seismic data. First, we review surface-wave properties in seismic data. Next, we discuss survey parameters and attributes. Then, a forward model of primaries is introduced, which takes the survey specification into account. Finally, the representation of the forward model is extended to include surface waves.

1

This chapter consists of parts of the papers submitted to and accepted by Geophys.

Prosp., Ishiyama et al. (2014a), Ishiyama et al. (2014d) and Ishiyama & Blacqui`ere (2014). The contents were re-arranged and made consistent with the other chapters of this thesis.

(16)

2.2 Surface-wave properties in seismic data

Surface waves often have larger amplitudes, lower frequencies and lower ap-parent velocities, compared to primaries. They are dispersive and multi-modal. Also, these properties vary spatially. In addition, surface waves are usually aliased due to a large spatial sampling interval in acquisition, e.g., a practical choice based on economical and operational constraints. In the case of blended acquisition, they are blended just as the primaries. These complex properties make surface-wave separation quite a challenge in processing.

Countless authors from several disciplines have discussed surface-wave characteristics since their existence was recognized. For a comprehensive overview of wave properties and a historical perspective of surface-wave applications, the papers by Socco & Strobbia (2004) and Socco et al. (2010) are referred to. In a land environment, surface waves are known as ground roll. Ground roll is commonly used as a synonym for Rayleigh waves because ground roll is in most cases caused by this type of wave. Rayleigh waves propagate parallel to the Earth’s surface without spreading their en-ergy into the Earth’s interior. They contain most of the enen-ergy within one wavelength from the surface and, therefore, their propagation is influenced only by near-surface properties. The propagation is characterized by dis-persion, i.e., different frequencies propagate with different phase velocities. Furthermore, the propagation is characterized by a multi-modal behaviour, i.e., each frequency propagates with several phase velocities at the same time, in other words, different modes exist simultaneously. In general, the funda-mental mode is primarily dominant over a wide frequency range. In land seismic data, the Rayleigh waves are stronger than the body waves, because seismic sources are located near the surface, and because they spread out with cylindrical divergence whereas the energy of the body waves decays more rapidly by spherical divergence. In a shallow-water environment, sur-face waves are referred to as mud roll or Scholte waves. Scholte waves have characteristics similar to those of Rayleigh waves. However, their propaga-tion velocity at the water bottom is slightly lower than that of Rayleigh waves at the free surface, due to the interaction with the overlying water. The dif-ference is related to the ratio between wavelength and water depth (Park

et al., 2005). For wavelengths considerably shorter than the water depth,

i.e., in deep-water conditions, the influence of the water layer becomes more significant. In contrast, for wavelengths longer than the water depth, i.e., in shallow-water conditions, the influence of the water layer is no longer sig-nificant and, in this case, Scholte waves can be regarded as Rayleigh waves (Bohlen et al., 2004). This is particularly true in a shallow-water

(17)

environ-2.2. Surface-wave properties in seismic data 9

Figure 2.1: Seismic data, (a) in the xy-t, (b) in the kxky-f , and (c) in the pxpy-f domain,

with a vertical section at the top and a horizontal time/frequency slice at the bottom. A dotted line in the section indicates the position of the slice, and vice versa. Red, magenta, green and blue arrows indicate the surface-wave fundamental mode, a surface-wave higher mode, refractions and reflections, respectively. Filled arrows indicate the true energy, and the whitish versions of the arrows indicate the aliased energy of the various types of waves. Notice that the surface waves are dispersive, multi-modal and aliased.

ment like the Gulf region in the Middle East where the water depth is about a few tens of meters. In OBC or ocean-bottom-node (OBN) seismic data, the Scholte waves are dominant both in hydrophone and geophone data since the receivers are deployed directly at the water bottom, whereas in streamer seismic data they are not so visible because the receivers are coupled to the water and not directly to the water bottom.

Figure 2.1 shows OBC seismic data in a shallow-water environment, 3-D OBC hydrophone data acquired oﬀshore Abu Dhabi, in the xy-t, the kxky

-f and the Radon (pxpy-f ) domain. This seismic data is the same as used

earlier; see also Figures 1.1(b) and 1.1(e). Surface waves are clearly seen, distinguishing themselves from primaries by larger amplitudes and lower fre-quencies. The fundamental mode and a higher mode can be identified. Each mode has a cone shape in the xy-t and the kxky-f domains and a hyperboloid

shape in the pxpy-f domain with a slowness larger than those of primaries.

(18)

shape in the xy-t domain, a wider cone shape in the kxky-f domain, and a

wider hyperboloid shape in the pxpy-f domain. Aliased energy is observed

especially in the kxky-f and the pxpy-f domains, wrapping around and

ob-liquely intersecting the true energy. Primaries are also found as the typical appearance for a flat-layered subsurface structure. Refractions are in a wider cone shape around a smaller two-way time range in the xy-t domain, in a narrower cone shape in the kxky-f domain, and in a more cylindrical shape in

the pxpy-f domain, compared to those of surface waves. Reflections appear

as flat events, i.e., around apexes of hyperboloids, with an infinitely small slowness over the whole two-way times in the xy-t domain, and around the origin (kx, ky) = (0, 0) and (px, py) = (0, 0) over their frequency range in the

kxky-f and the pxpy-f domain.

2.3 Seismic survey parameters and attributes

For 3-D seismic surveys, the relevant survey parameters are the four spatial sampling intervals and apertures of the template geometry (Vermeer, 2012). The four spatial sampling intervals are the receiver and source intervals, ∆xd

and ∆yd for receivers and ∆xs and ∆ys for sources, each in two sampling

directions that are usually orthogonal. The four spatial sampling apertures are the receiver and source apertures, Xd and Yd for the receivers and Xs

and Ys for the sources, oriented in the same way as the above four spatial

coordinates. Proper spatial sampling can be thought of as the ability to properly reconstruct seismic wavefields.

Figure 2.2 shows typical template geometries and their basic subsets. Here, the x-direction is considered as the in-line direction. For an areal geo-metry, the basic subset is a common-source gather or a common-receiver gather. For common-source gathers, receivers are deployed on a densely spaced grid (∆xd and ∆yd fine) while sources are on a sparsely spaced grid

(∆xsand ∆ys coarse). Receiver-spread widths (Xd and Yd) specify the

max-imum apertures. In contrast, for common-receiver gathers, sources are ar-ranged on a densely spaced grid (∆xs and ∆ys fine) while receivers are on a

sparsely spaced grid (∆xd and ∆yd coarse). Source-spread widths (Xs and

Ys) specify the maximum apertures. For an orthogonal geometry, the basic

subset is a cross-spread gather, where receiver-point and source-point inter-vals are quite fine (∆xd and ∆ys fine), whereas receiver-line and source-line

intervals are often coarse (∆yd and ∆xs coarse). Receiver-line and

source-line lengths specify the maximum apertures (Xd and Ys). It should be noted

(19)

2.3. Seismic survey parameters and attributes 11

Figure 2.2: Template geometries and their basic subsets, (a) and (b) for common-source

gathers, (c) and (d) for common-receiver gathers, (e) and (f) for cross-spread gathers. Red and blue indicate source and receiver locations, respectively.

can be directly calculated from these survey parameters, ∆xbin = ∆xd 2 , (2.1) ∆ybin= ∆ys 2 , (2.2) Mx= Xd 2∆xs xrep, (2.3) My = Ys 2∆yd yrep, (2.4) Mbin = MxMy, (2.5) ρbin= Mbin

∆xbin∆ybin

(20)

where ∆xbin and ∆ybin are the bin size; Mx, My and Mbin the in-line, the

cross-line and the nominal fold; xrep and yrep the template-repeat factors

resulting from rolling the template in the in-line and the cross-line directions while repeating a part of the template; ρbin the trace density.

Two of the four spatial coordinates, the set {∆xb, ∆yb, Xb, Yb}, specify

the spatial sampling of the basic subset, where the subscript b can be either

d or s, independently for each survey parameter but not in arbitrary

com-binations. Two other coordinates, the set {∆xB, ∆yB}, specify the

spa-tial redundancy of the basic subsets, i.e., the fold, where again the sub-script B can be d or s. Those maximum apertures, XB and YB, are

usu-ally the same as Xb and Yb, resulting from forming the template or from

rolling the patch with repeating receiver and/or source locations. For in-stance, the set {∆xd, ∆ys, Xd, Ys} specifies the spatial sampling of a

cross-spread gather, i.e., {∆xb, ∆yb, Xb, Yb} = {∆xd, ∆ys, Xd, Ys}, whereas the set

{∆xs, ∆yd} specifies the spatial redundancy of the cross-spread gather, i.e.,

{∆xB, ∆yB} = {∆xs, ∆yd}. Again note that traditional survey attributes,

maximum oﬀset and largest minimum oﬀset, can be directly obtained by

Xmax= Xb 2 , (2.7) Ymax= Yb 2 , (2.8) Xmin= ∆xB, (2.9) Ymin= ∆yB, (2.10)

where Xmax and Ymax are the maximum oﬀsets; Xmin and Ymin the largest

minimum oﬀsets. Furthermore, the survey eﬀort can be defined as a combined attribute of these survey parameters relative to a reference template,

Cxb = ∆xbref ∆xb Xb Xbref , (2.11) Cyb = ∆ybref ∆yb Yb Ybref , (2.12) CxB = ∆xB ref ∆xB xrep, (2.13) CyB = ∆yB ref ∆yB yrep, (2.14) Cx= CxbCxB, (2.15) Cy = CybCyB, (2.16)

(21)

2.4. Forward model of seismic data 13

Cb= CxbCyb, (2.17)

CB= CxBCyB, (2.18)

C = CxCy

= CbCB, (2.19)

where C_{·}is the survey eﬀort for each component; the subscript ‘ref’ denotes ‘reference’. The attributes in terms of symmetry can be also defined as

A∆xb = ∆xb ∆yb , (2.20) AXb = Yb Xb , (2.21) A∆xB = ∆xB ∆yB , (2.22)

where A∆xb and A∆xB are the aspect ratios of the spatial sampling intervals; AXb is the aspect ratio of the spatial sampling apertures. These survey parameters and attributes determine the specification of seismic data.

2.4 Forward model of seismic data

To describe 3-D seismic data with the survey parameters, the so-called WRW model of Berkhout (1982), a representation of primaries in the space-frequency domain, is introduced. Following the author, the word ‘detector’ is used for ‘receiver’. In this model, primaries can be described for each monochromatic component by operator matrices as

P (zd, zs) = D (zd)

∑

m

[W (zd, zm) R (zm, zm) W (zm, zs)] S (zs) . (2.23)

Each matrix multiplication represents a multi-dimensional spatial convolu-tion, and each element of every matrix is a complex number containing amp-litude and phase information. It should be noted that this equation is valid for stationary (parts of non-stationary) acquisition geometries. These matrices have the following meaning:

• P(zd, zs) is the data matrix of the primary wavefields recorded by

de-tectors at depth zd due to sources at depth zs, both depths being near

the surface z0 (Figure 2.3). A row and a column of the data

(22)

Figure 2.3: A common-source gather ⃗Ps(zd, zs) for a certain source s at ⃗xs with all

frequencies, (a) in the 2-D d-t domain, and (b) in the 3-D xy-t domain. A common-receiver gather ⃗P_d†(zd, zs) can be considered in the same way but for a certain detector d at ⃗xd.

Seismic data P(zd, zs) with all frequencies, (c) in the 3-D ds-t domain, and (d) in the 3-D ds-ω domain, are described, (e) and (f) by data matrix P(zd, zs) for each monochromatic

component.

and a certain source location s at (⃗xs, zs) = (xs, ys, zs), respectively.

(23)

2.4. Forward model of seismic data 15

Figure 2.4: Schematic reflectivity matrix R(zm, zm) for each monochromatic component. ∆xs= xs+1−xsand ∆ys= ys+1−ys; see e.g. Figures 2.3(a) and 2.3(b).

Furthermore, Xd = Max(xd)− Min(xd), Yd = Max(yd) − Min(yd),

Xs = Max(xs)− Min(xs) and Ys = Max(ys)− Min(ys), where Max(.)

and Min(.) are the maximum and the minimum values for each com-ponent. According to the row and column numbering, a row ⃗P_d†(zd, zs)

constitutes a common-receiver gather, and a column ⃗Ps(zd, zs)

consti-tutes a common-source gather. The dagger symbol † denotes a row vector. Common-oﬀset gathers and common-mid-point (CMP) gath-ers are identified in the data matrix as diagonals and anti-diagonals, respectively. An element Pds(zd, zs) constitutes one frequency

compon-ent of a single trace recorded by the dth detector and shot by the sth source, i.e., P (⃗xd, zd; ⃗xs, zs; ω).

• R(zm, zm) is the reflectivity matrix representing the conversion of the

incident wavefields into the reflected wavefields at depth zm (Figure

2.4). In the same manner as the data matrix, a row and a column of the reflectivity matrix correspond to a certain grid-point location i at (⃗xi, zm) = (xi, yi, zm) and a certain grid-point location j at (⃗xj, zm) = (xj,

yj, zm), respectively. An element Rij(zm, zm) represents one frequency

component of a reflected wavefield at the ith grid point generated by a delta-function source at the jth grid point, i.e., R(⃗xi, zm; ⃗xj, zm; ω).

• S(zs) is the source matrix, its columns ⃗Ss(zs) containing the source

properties at (⃗xs, zs), i.e., S(⃗xs, zs; ω).

• W(zm, zs) is the propagation matrix, its columns containing the

(24)

element Wjs(zm, zs) represents a direct wavefield at (⃗xj, zm) caused by

a delta-function source at (⃗xs, zs), i.e., W (⃗xj, zm; ⃗xs, zs; ω).

• W(zd, zm) is the propagation matrix, its columns containing the

up-going wavefields at zd generated by delta-function sources at zm. An

element Wdi(zd, zm) represents a direct wavefield at (⃗xd, zd) caused by

a delta-function source at (⃗xi, zm), i.e., W (⃗xd, zd; ⃗xi, zm; ω).

• D(zd) is the detector matrix, its rows ⃗D†_d(zd) containing the detector

properties at (⃗xd, zd), i.e., D(⃗xd, zd; ω).

Later on, the spatial coordinate z will be sometimes omitted for brevity. In the case that W represents one-way wave propagation, Equation (2.23) can be viewed as the Born approximation to the Lippmann-Schwinger equation with an one-way Green’s function W. It should be noted that Equation (2.23) is a general description of primaries. In the case of blended acquisition with the blending operator Γ, SΓ can replace S in this equation (Berkhout, 2008). For a full-wavefield model, including not only the primary wavefields but also the secondary wavefields generated by surface-related and internal multiples, a full-wavefield propagation matrix or Green’s function G can replace W in this equation (Kumar et al., 2014). Also note that with all frequencies the matrices in the frequency domain can be envisaged more naturally in the time domain. Since the discrete Fourier transform is invertible, either domain can be used depending on the purpose.

2.5 Forward model of surface waves

In the data matrix P, every column vector corresponds to a common-source gather. Using reciprocity, a common-receiver gather can be thought of as a common-source gather, and cross-spread gathers can be re-arranged into common-source gathers by sorting the cross-spread gathers into the corres-ponding vectors with lexicographical ordering. In this representation, 3-D seismic data generated by a single source containing both subsurface signals and surface waves can be described for each monochromatic component as

⃗

Pstot= ⃗Ps+ ⃗Ns, (2.24)

or with omitting the subscript s for brevity,

⃗

Ptot= ⃗P + ⃗N , (2.25)

where ⃗Ptot represents the common-source gather; ⃗P the subsurface signals; ⃗

(25)

2.5. Forward model of surface waves 17

signals’ refer to all events except the surface waves, i.e., ⃗P includes refractions,

reflections, surface-related and internal multiples, etc. The surface waves ⃗N can be described as

⃗ N = nn ∑ n ⃗ Nn, (2.26) ⃗ Nn= DHnS⃗n, (2.27)

or with assuming a unit-detector response, i.e., D = I,

⃗

Nn= HnS⃗n, (2.28)

where ⃗N and ⃗Nn are the total and the modal surface waves; ⃗Sn contains the

source properties for a surface-wave mode. The subscript n denotes a certain surface-wave mode; nnis the number of surface-wave modes. The matrix Hn

has the following meaning:

• Hn(zd, zs), with the spatial coordinates zd and zs being near the

sur-face z0, is the propagation matrix for a certain surface-wave mode n, its

columns containing the horizontally propagating wavefields at zd

gen-erated by delta-function sources at zs. An element (Hn)ds(zd, zs)

rep-resents a direct wavefield at (⃗xd, zd) caused by a delta-function source

at (⃗xs, zs), i.e., Hn(⃗xd, zd; ⃗xs, zs; ω).

As already stated, in the case of blended acquisition with the blending op-erator ⃗Γ, Sn⃗Γ can replace ⃗Sn in Equation (2.28). Also remember that with

all frequencies the vectors in the frequency domain can be envisaged more naturally in the time domain.

Following Aki & Richards (2002) and Socco & Strobbia (2004), each ele-ment of Hn can be described by

Hn(zd, zs) = ( e−⃗λn(ω)·⃗r √ |⃗r| ) e−iω ⃗pn(ω)·⃗r_, _(2.29)

where ⃗r is the receiver location relative to the source location, i.e., ⃗r =

(⃗xd, zd)− (⃗xs, zs) and |⃗r| is the oﬀset; ⃗λn(ω) is the frequency-dependent

modal attenuation factor; ⃗pn(ω) is the frequency-dependent modal slowness,

describing the dispersion surface. The factor √1

|⃗r| represents the cylindrical

spreading, and the factor e−⃗λn(ω)·⃗r_{represents the intrinsic attenuation of the} amplitude. In general, the eﬀect of the geometrical attenuation is much lar-ger than that of the intrinsic attenuation. The term ⃗pn(ω)· ⃗r represents the

(26)

Figure 2.5: The model parameters. (a) The dispersion surfaces of the surface waves with

a vertical section at the top and a horizontal frequency slice at the bottom. (b) The source wavelet of the fundamental mode. (c) The source wavelet of the higher mode. A dotted line in the section indicates the position of the slice, and vice versa. Red and magenta indicate the fundamental mode and the higher mode, respectively.

travel-time phase shift. Equations (2.26), (2.28) and (2.29) show that ⃗N is

parameterized by ⃗λn, ⃗pn and ⃗Sn. With the assumption that ⃗λn is zero, i.e.,

the intrinsic attenuation is negligible compared to the geometrical attenu-ation, ⃗N is forward-modelled with the model parameters ⃗pn and ⃗Sn.

Figures 2.5 to 2.8 show two examples of the forward model. A common-source gather is considered, consisting of a centered common-source and a receiver spread with a spatial sampling interval of 50 m and a width of 3200 m both in the x- and the y-direction, i.e., {∆xb, ∆yb, Xb, Yb} = {50 m, 50 m, 3200 m,

3200 m}, A∆xb = 1 and AXb = 1. The number of surface-wave modes is taken as two, the fundamental mode and one higher mode. For the first example, the model parameters, ⃗pn and ⃗Sn, are selected as exhibited in Figure 2.5.

Here, ⃗pn does not vary in azimuth. Figure 2.6 depicts the forward-modelled

surface waves in the xy-t, the kxky-f and the pxpy-f domain. The surface

waves are clearly seen; see Figure 2.1 for comparison with the observed surface waves. For the second example, ⃗Snis the same as in the first example, but ⃗pn

has azimuthal variability as exhibited in Figure 2.7. Figure 2.8 depicts the results of the forward model. The eﬀect of ⃗pnon phase shift can be identified

(27)

2.5. Forward model of surface waves 19

Figure 2.6: The forward-modelled surface waves, (a) in the xy-t, (b) in the kxky-f , and

(c) in the pxpy-f domain, with a vertical section at the top and a horizontal time/frequency

slice at the bottom. A dotted line in the section indicates the position of the slice, and vice versa.

by the change of the surface-wave shapes. Notice that the surface waves are more compressed and more sparsely distributed in the kxky-f and the pxpy-f

domain than in the xy-t domain.

Based on the above observations, it is also possible to express the modal surface waves ⃗Nn as

⃗˜

Nn= L ⃗Nn,

∴ ⃗Nn= LHN⃗˜n, (2.30)

or with normalization in the space-frequency (xy-f ) domain,

⃗˜

Nn= LB ⃗Nn,

∴ ⃗Nn= B−1LHN⃗˜n, (2.31)

where L and B are the Radon and the normalization operators (Appendix). The superscript H denotes ‘conjugate transposition’. The tilde symbol ˜· indicates the result of the Radon transform in the pxpy-f domain. The

(28)

Figure 2.7: The model parameters. The dispersion surfaces, of the surface waves with a

bird’s-eye view at the left, a vertical section at the upper right and a frequency slice at the lower right. A dotted line in the section indicates the position of the slice, and vice versa. Red and magenta indicate the fundamental mode and the higher mode, respectively.

(McMechan & Yedlin, 1981). The normalization in the xy-f domain followed by the Radon transform enhances the phase information in the pxpy-f

do-main (Park et al., 1998). Equations (2.26) and (2.31) show that ⃗N is also

explained by the surface-wave samples N⃗˜n in the pxpy-f domain. This is

possible because the surface waves are well distinguished and their samples can be selected with relative ease in this domain in terms of their amplitudes and areal separation from the subsurface signals.

2.6 Conclusions from this chapter

• For 3-D seismic surveys, the relevant survey parameters are the four

spatial sampling intervals and apertures of the template geometry. These survey parameters and the resulting attributes can determine the specification of the seismic data. The seismic experiment can be de-scribed with the survey parameters and, therefore, be forward-modelled with taking those into account.

(29)

2.6. Conclusions from this chapter 21

Figure 2.8: The forward-modelled surface waves, (a) in the xy-t, (b) in the kxky-f , and

(c) in the pxpy-f domain, with a vertical section at the top and a horizontal time/frequency

slice at the bottom. A dotted line in the section indicates the position of the slice, and vice versa.

• Surface waves can be parameterized by the frequency- and mode-dependent

dispersion surfaces and source properties and, therefore, forward-modelled with this parameterization. This addresses the surface-wave properties such as dispersion and multi-modes.

(30)

2.7 Appendix: Radon transform

Radon transform

The discrete linear Radon transform of seismic data can be described for each monochromatic component as

⃗˜

Ptot = L ⃗Ptot, (2.32)

where ⃗Ptotrepresents the common-source gather; the tilde symbol ˜· indicates

the result of the Radon transform in the Radon (pxpy-f ) domain. L is the

Radon operator, its rows representing a series of the slowness ⃗p, and its

columns representing a series of the receiver location relative to the source location ⃗r, i.e.,|⃗r| is the oﬀset. Each element of L can be provided by

L(⃗p, ⃗r, ω) = eiω ⃗p·⃗r. (2.33)

The number of rows in the matrix based on the sampling interval and aperture of the series of ⃗p, i.e., the number of traces in the pxpy-f domain, can be

chosen in accordance with the resolution and bandwidth required for seismic events to be imaged in the pxpy-f domain (Turner, 1990).

The inverse Radon transform can be described for each monochromatic component as

⃗

Ptot= LHP⃗˜tot, (2.34)

where LH is the adjoint of L, i.e., an approximate inverse of L. Since L is not an orthogonal matrix, i.e., LLH ̸= I, a least-squares method is approached to solve the inverse problem by minimizing the objective function

J (P⃗˜tot) = ∑ ω ⃗Ptot− LHP⃗˜tot 2 . (2.35)

To obtain a high resolution result, several approaches have been developed, such as Thorson & Claerbout (1985) with a stochastic sparsity constraint, Sacchi & Ulrych (1995) with a Chauchy-norm sparsity constraint, and Wang

et al. (2010) with a greedy sparsity approach. For a comprehensive overview

of sparse Radon transforms, the paper by Trad et al. (2003) is referred to.

Normalized Radon transform

The normalized Radon transform, the normalization in the xy-f domain followed by the Radon transform, can be described for each monochromatic

(31)

2.7. Appendix: Radon transform 23

component as

⃗˜

Ptot= LB ⃗Ptot, (2.36)

where B is the normalization operator. B is a diagonal matrix, where both dimensions represent receiver location. A non-zero diagonal element of B for a kth receiver location, Bkk, can be provided by

Bkk = |P k|

|Pk|2+ ε2

, (2.37)

where ε is a stabilization factor; Pk is an element of ⃗Ptot for the kth receiver

location. The normalized Radon transform enhances the phase information in the pxpy-f domain (Park et al., 1998).

The inverse normalized Radon transform can be described for each mono-chromatic component as

⃗

Ptot= B−1LHP⃗˜tot. (2.38)

Similar to the inverse Radon transform, a least-squares method is used to solve the inverse problem by minimizing the objective function

J (P⃗˜tot) = ∑ ω ⃗Ptot− B−1LHP⃗˜tot 2 . (2.39)

(32)

(33)

Chapter 3 Surface-wave estimation and

separation using closed loop

1

3.1 Introduction

To separate out surface waves from seismic data, many approaches have been developed in processing. Traditionally, areal filtering methods are used, in which the seismic data are transformed to a certain domain where the surface waves and the body waves are more easily separated in terms of amplitude, frequency and apparent velocity. Examples are the xy-f , the kxky-f and the

Radon (pxpy-τ ) domain (e.g. Yilmaz, 2001). However, the aliased energy of

the surface waves often overlaps with the body waves, making it still diﬃcult to separate them out even in one of these domains. Furthermore, a fixed set of filtering parameters is commonly used across the whole seismic data, ignoring the spatial variability of the surface waves. In addition to filtering methods, model-based methods are also used. In these methods, a set of model para-meters is required to obtain the forward-modelled surface waves. However, the model parameters are usually unknown. This requires labour-intensive testing for the model-parameter estimation. Furthermore, imperfect model parameters result in a substantial residual between the forward-modelled and the observed surface waves. This residual is seldom re-evaluated and usually left as is. Le Meur et al. (2008, 2010) proposed a data-driven, data-adaptive and model-based method. In their method, the model-parameter estimation is highly automated. However, the model parameters are estimated section by section in a 3-D cube, not in a true 3-D fashion, possibly resulting in model

1

This chapter is a part of the paper submitted to and accepted by Geophys. Prosp., Ishiyama et al. (2014a). The contents were made consistent with the other chapters of this thesis.

(34)

inconsistencies in the crossline direction. Furthermore, the residual between the forward-modelled and the observed surface waves is not re-evaluated and the model parameters are not estimated, although data-adaption may re-duce the residual and compensate for the imperfection of the model para-meters. Recently, near-surface model-based methods were introduced, using the so-called surface-wave method. Strobbia et al. (2010, 2011) proposed a near-surface model-based method. In their method, the model parameters are defined as elastic properties at all grid points of a near-surface model that has been estimated by an inversion scheme from surface-wave properties such as dispersion and multi-modes. In this near-surface model, the surface waves are then forward-modelled by a suitable algorithm, e.g., a 3-D elastic finite-difference code. However, the near-surface model obtained by the inversion scheme is significantly non-unique due to the uncertainties in the interpret-ation of the surface-wave properties as well as the complexity of the inverse problem (Ivanov et al., 2013). Furthermore, again, the residual between the forward-modelled and the observed surface waves is not re-evaluated and the near-surface model is not re-estimated based on this residual. It is difficult to feed this residual directly back to the near-surface model, because the near-surface model is estimated indirectly from the surface-wave properties. More recently, methods based on full-waveform inversion (FWI) have started to appear (e.g. Ernst, 2013). In this approach, the model parameters de-scribe a near-surface model that is estimated directly from the surface waves themselves, not indirectly via surface-wave properties such as dispersion and multi-modes. In this way, these surface-wave properties are not required in explicit form. Instead, a realistic model of the wavefield propagation in the near-surface is required to explain all surface-wave phenomena. This increases the complexity of the forward modelling significantly and makes the inverse problem more challenging. For the forward modelling, a 3-D elastic finite-difference algorithm is most popular. However, this approach is not flexible in terms of the choice of the parameterization. The model parameters should be defined as elastic properties at all grid points of a subsurface model. Other types of parameterizations could be better suited, although mapping them to grid-based elastic properties may not be feasible. In addition, this algorithm is computationally expensive, in particular if the wavefields need to be simulated over and over again in an inversion scheme. For these reasons, this approach is still demanding at present.

To address the challenges of surface waves and adopt recent advances in processing, a method of 3-D surface-wave estimation and separation using a closed-loop or iterative approach is considered. The closed-loop approach has been developed for several stages of processing (Berkhout, 2013). In this

(35)

3.2. Inverse problem 27

approach, at each processing stage, an optimal parameterization is used to describe the inverse problem. The closed loop contains a forward-modelling module, allowing the evaluation of the residual between forward-modelled and observed seismic data and, therefore, feedback between the model parameters and the observed seismic data. These parameters are iteratively estimated such that the residual is minimized. Examples are Kontakis & Verschuur (2014) for deblending, Lopez & Verschuur (2014) for multiple estimation and separation, Mulder & Plessix (2004); Davydenko & Verschuur (2014); Soni & Verschuur (2014) for full-wavefield migration, Staal et al. (2014) for joint migration inversion, etc. Furthermore, data reconstruction has been widely developed using this approach, based on decomposition of seismic data in various basis functions, such as Xu et al. (2005); Abma & Kabir (2006); Zwartjes & Gisolf (2006); Schonewille et al. (2009) for Fourier-based, Wang et al. (2010) for linear Radon-based, Kabir & Verschuur (1995); Trad

et al. (2003) for parabolic Radon-based, Herrmann & Hennenfent (2008)

for curvelet-based, and Berkhout & Verschuur (2006); Kutscha & Verschuur (2012) for focal-based. For each basis function, the corresponding seismic event is parameterized by the samples in the respective domain such as Four-ier, linear Radon, parabolic Radon, curvelet and focal domain. The model parameters are iteratively estimated under certain sparsity constraints such that the forward-modelled (inverse-transformed) seismic data match the ob-served seismic data.

The goal of this chapter is to establish a method of 3-D surface-wave estimation and separation using the closed-loop approach. First, we discuss the inverse problem for the forward model of surface waves. Then, the method is introduced and illustrated in a step-by-step manner. Furthermore, it is extended to include an iterative and a conventional filtering step. Finally, it is compared with a conventional filtering method.

3.2 Inverse problem

If the surface waves N are estimated with the model parameters ⃗⃗ˆ pn and

⃗

Sn, and subtracted from the observed seismic data ( ⃗P + ⃗N ), a

subsurface-signal estimate, P = ( ⃗⃗ˆ P + ⃗N )−N , is obtained. Here, the hat symbol ˆ⃗ˆ ·

indicates ‘estimated’. If N is not perfectly estimated, a non-zero residual,⃗ˆ

∆ ⃗N = ⃗N −N , is left. In this case,⃗ˆ P includes the term ∆ ⃗⃗ˆ N , i.e., P =⃗ˆ ⃗

P + ∆ ⃗N . To reduce this residual, an adaptive filter An for each

surface-wave mode is applied. The matrix An is a diagonal matrix, where both

(36)

model inaccuracies in terms of (intrinsic) attenuation and spatial variability. To further minimize the residual, an iterative inversion scheme is used. The surface waves N are obtained by minimizing the objective function⃗ˆ

J (⃗pn, ⃗Sn; An) = ∑ ω ⃗P + ∆ ⃗N 2 =∑ ω ( ⃗P + ⃗N )−N⃗ˆ 2 =∑ ω ( ⃗P + ⃗N )− nn ∑ n AnHnS⃗n 2 , (3.1)

see also Equation (2.28). The model parameters ⃗pn and ⃗Sn, as well as the

adaptive filter An, are determined in such a way that the residual ( ⃗P + ∆ ⃗N )

between the surface waves N and the observed seismic data ( ⃗⃗ˆ P + ⃗N ) is

minimized. Notice that this minimization scheme is supposed to act on N⃗ˆ

only, i.e., the term ∆ ⃗N should be minimized only. Therefore, ⃗P remains

untouched. To this end, a signal-protecting scheme is used as will be discussed later. Consequently, the resulting residual ( ⃗P + ∆ ⃗N ) closely corresponds to

the subsurface signals ⃗P . In this thesis, we call this approach ‘SWES’, which

stands for ‘surface-wave estimation and separation’.

3.3 Closed loop

To solve the inverse problem, an iterative closed-loop approach is used. Figure 3.1 shows the concept, consisting of three modules: parameter es-timation/selection, forward modelling, and adaptive subtraction. For each loop (i), the model parameters are estimated from the residual update, ( ⃗P + ∆ ⃗N )(i); the surface-wave estimate N⃗ˆ(i) is built by forward modelling; this estimate is adaptively subtracted from the observed seismic data ( ⃗P + ⃗N )

to obtain the subsurface-signal estimateP⃗ˆ(i). This estimate becomes the re-sidual update ( ⃗P +∆ ⃗N )(i+1)for the next loop. The procedure is iterated until it has reached a stopping criterion. In practice, each i-loop contains an inner

j-loop for each of the surface-wave modes. For each inner loop (i, j), the

para-meters are estimated from the residual update ( ⃗P +∆ ⃗N )(i,j), and the resulting modal surface-wave update ∆ ⃗N_j(i)is obtained. The parameter estimation and the surface-wave updating will be discussed later. This update is summed in terms of j to build the intermediate and the total surface-wave update,

(37)

3.3. Closed loop 29

Figure 3.1: The closed-loop concept for SWES.

respectively, i.e., ∆ ⃗N(i,j) ₌

nj ∑

j

∆ ⃗N_j(i) and ∆ ⃗N(i) _{= ∆ ⃗}_N(i,j=nn) ₌

nn ∑

j

∆ ⃗N_j(i)

where n_{.} is the corresponding upper limit for the index of summation. This update is subsequently summed in terms of i to build the total surface-wave estimate, i.e.,N⃗ˆ(i) =

ni ∑

i

∆ ⃗N(i). Notice that the loop is seamlessly closed once the parameters are estimated and the resulting modal surface-wave update ∆ ⃗N_j(i) is obtained. This algorithm can be expressed as

( ⃗P + ∆ ⃗N )(i) = ( ⃗P + ⃗N )−N⃗ˆ(i−1), (3.2) ( ⃗P + ∆ ⃗N )(i,j)= ( ⃗P + ⃗N )(i)− ∆ ⃗N(i,j−1), (3.3) ∆ ⃗N(i,j)= ∆ ⃗N(i,j−1)+ ∆ ⃗N_j(i), (3.4)

⃗ˆ

N(i)=N⃗ˆ(i−1)+ ∆ ⃗N(i), (3.5)

⃗ˆ

P(i)= ( ⃗P + ⃗N )−N⃗ˆ(i). (3.6)

Equations (3.2), (3.5) and (3.6) express the outer i-loop, and Equations (3.3) and (3.4) the inner j-loop.

(38)

Figure 3.2: The closed-loop kernel for SWES.

Now, the parameter estimation and the surface-wave updating are dis-cussed. For a loop (i, j = n), the parameters and/or updates, ⃗pn(i), ∆ ⃗Sn(i)

and A(i)n , are estimated from the residual update ( ⃗P + ∆ ⃗N )(i,n). Here,

⃗

pn(i)= ⃗pn(i−1)+∆⃗pn(i), and ∆ ⃗Sn(i)= ⃗Sn(i)−⃗Sn(i−1). First, at the parameter

estim-ation/selection module, the dispersion surface ⃗pn(i)is estimated by picking the

amplitude maxima of the residual update added to the modal surface-wave estimate in the normalized pxpy-f domain. As already stated, the

disper-sion surface is imaged as the amplitude maxima in the pxpy-f domain, and

the normalization enhances the phase information in the pxpy-f domain and

makes the picks accurate. Equation (2.29) then provides H(i)n , and the modal

surface-wave update H(i)n ∆ ⃗Sn(i−1) is subsequently obtained by forward

mod-elling. Second, at the adaptive subtraction module, the source update ∆ ⃗Sn(i)

(39)

3.3. Closed loop 31

the following misfit ∑

ω

( ⃗P + ∆ ⃗N )(i,n)− α(i)_n H(i)_n ∆ ⃗S_n(i−1) 2, (3.7) and the modal surface-wave update is obtained by

∆ ⃗S_n(i)= α(i)_n ∆ ⃗S_n(i−1), (3.8)

∆ ⃗N(i)_n = H(i)_n ∆ ⃗S_n(i). (3.9)

It should be noted that α(i)n , combined for all frequencies, defines a

convo-lution filter in the time domain. Optionally, the length of this filter in time can be restricted. Third, still in the adaptive subtraction module, the adapt-ive filter A(i)n is estimated by solving the local (area-by-area) Wiener filter,

obtained by minimizing the following misfit ∑

ω

( ⃗P + ∆ ⃗N )(i,n)− A(i)_n ∆ ⃗N(i)_n 2

, (3.10)

and the modal surface-wave update is finally obtained by

∆ ⃗N_n(i)= A(i)_n ∆ ⃗N(i)_n . (3.11)

Again, A(i)n , combined for all frequencies, defines a convolution filter in the

time domain, and the length of this filter in time can be optionally restricted. In practice, for Equation (3.10), a non-zero diagonal element of A(i)n for a kth

receiver location, (A(i)n )kk, is provided by

(A(i)_n )kk= Trk [ ( ⃗P + ∆ ⃗N )(i,n)_{(∆ ⃗}_N(i) n )H ] Trk [ ∆ ⃗N(i)_n (∆ ⃗N(i)_n )H_{+ ε}2_I] , (3.12)

where ε is a stabilization factor; Trk(.) denotes the trace of a matrix, i.e., the

sum of its diagonal elements, in a spatial window centered on the kth_receiver

location. In Equation (3.12), the numerator is the cross-correlation of ( ⃗P +

∆ ⃗N )(i,n) _{and ∆ ⃗}_N(i)

n , and the denominator is the auto-correlation of ∆ ⃗N

(i)

n .

The spatial window of the cross-correlation and the stabilization factor of the auto-correlation control the signal-protecting scheme. A wider spatial window and a larger stabilization factor lead to a Wiener filter protecting the signals ⃗P but probably predicting less energy for the surface waves ∆ ⃗N

(40)

lead to a Wiener filter predicting more energy for the surface waves ∆ ⃗N but

possibly attacking the signals ⃗P of ( ⃗P + ∆ ⃗N )(i,n)_{. This means that there is a}

trade-oﬀ between conservative signal protection and aggressive surface-wave estimation. The choice of the spatial window and the stabilization factor can be determined by testing.

Figure 3.2 shows the closed-loop kernel, summarizing the roles of the outer i-loop, the inner j-loop, the parameter estimation and the surface-wave updating. Notice that, in addition to the total surface-surface-wave estimate in Equation (3.5), the modal surface-wave estimate can be also obtained for a surface-wave mode, i.e., j = n, as

⃗ˆ N_n(i) =N⃗ˆ_n(i−1)+ ∆ ⃗N_n(i), (3.13) meaning N⃗ˆn(i)= ni ∑ i ∆ ⃗Nn(i).

3.4 Step-by-step illustration of the closed loop

To illustrate the closed-loop approach, we use 3-D OBC hydrophone data acquired oﬀshore Abu Dhabi in a shallow-water environment. A cross-spread gather is considered, consisting of a receiver line in the x-direction and a source line in the y-direction, each with a spatial sampling interval of 50 m and a length of 3200 m, i.e., {∆xb, ∆yb, Xb, Yb} = {50 m, 50 m, 3200 m, 3200 m},

A∆xb = 1 and AXb = 1. This fills a 3-D cube in the xy-t domain. Bad traces are edited. Figure 3.3 shows the seismic data in the xy-t, the kxky-f and

the pxpy-f domain. This seismic data is the same as used earlier to discuss

the surface-wave properties; see also Figure 2.1. Surface waves are clearly seen, distinguishing themselves from subsurface signals by larger amplitudes and lower frequencies. The fundamental mode and one higher mode can be identified, both with a slowness larger than those of subsurface signals.

The initial parameters for the closed loop are either based on a priori knowledge or determined from the seismic data. Figure 3.4(a) shows the initial dispersion surfaces ⃗pn(0)for the fundamental mode and the higher mode.

In this example, the initial dispersion surfaces were obtained by manually picking a 1-D dispersion curve, e.g., the amplitude maxima of a 2-D section in the pxpy-f domain, and by filling the picks along azimuth at each frequency.

Thus, the initial dispersion surfaces have a circular shape in a frequency slice. The initial source wavelets ∆ ⃗Sn(0) were prepared as a spike, minimum-phased

and band-limited in the frequency range. The initial adaptive filters A(0)n are

(41)

3.4. Step-by-step illustration of the closed loop 33

Figure 3.3: The seismic data, (a) in the xy-t, (b) in the kxky-f , and (c) in the pxpy-f

domain, with a vertical section at the top and a horizontal time/frequency slice at the bottom. A dotted line in the section indicates the position of the slice, and vice versa.

For the first loop and the fundamental mode (i, j) = (1, 1), the residual update ( ⃗P + ∆ ⃗N )(1) (to become also ( ⃗P + ∆ ⃗N )(1,1)) equals the seismic data ( ⃗P + ⃗N ); see Figure 3.3(a). Figures 3.5(a) and 3.5(b) show the initial

disper-sion surface ⃗p₁(0) and the initial source wavelet ∆ ⃗S₁(0), respectively. First, the dispersion surface ⃗p₁(1) is estimated. Figures 3.5(c) and 3.5(d) show the re-sidual update in the normalized pxpy-f domain, in which ⃗p1(1) is estimated by

automatically picking the amplitude maxima. The normalization enhances the phase information and makes the picks accurate; see Figure 3.3(c) for comparison. The fundamental mode is clearly seen as a hyperboloid shape with a slowness larger than those of the higher mode and subsurface signals. Aliased energy is observed, obliquely intersecting the true energy. To enhance the true energy and make the picks more accurate, optionally, a smoothing fil-ter and a windowing filfil-ter along ⃗p₁(0) can be applied. Figure 3.5(e) shows the estimated dispersion surface ⃗p₁(1) that explains the azimuthal variability; see Figure 3.5(a) for comparison. Figure 3.5(f) shows the resulting fundamental-mode update H(1)₁ ∆ ⃗S₁(0). Second, the source update ∆ ⃗S₁(1) is estimated by solving the global Wiener filter α(1)₁ and multiplying this filter with the

(42)

ini-Figure 3.4: The dispersion surfaces. (a) The initial dispersion surfaces, and (b) the estimated dispersion surfaces after three iterations, of the surface waves with a bird’s-eye view at the left, a vertical section at the upper right and a frequency slice at the lower right. A dotted line in the section indicates the position of the slice, and vice versa. Red and magenta indicate the fundamental mode and the higher mode, respectively.

(43)

tial source wavelet ∆ ⃗S₁(0). Figure 3.5(g) shows the estimated source update ∆S(1)₁ that explains the source bubble oscillation and ringing; see Figure 3.5(b) for comparison. Figure 3.5(h) shows the resulting fundamental-mode update ∆ ⃗N(1)₁ . Third, the adaptive filter A(1)₁ is estimated by solving the local Wiener filter. This filter compensates for model inaccuracies, e.g., re-lated to (intrinsic) attenuation and spatial variability. Figure 3.5(i) shows the resulting fundamental-mode update ∆ ⃗N₁(1). Finally, ∆ ⃗N₁(1)itself becomes the fundamental-modeN⃗ˆ₁(1) at this stage. For the higher mode (j = 2), the same procedure is followed. Figure 3.5(j) shows the resulting higher-mode update ∆ ⃗N₂(1). Figures 3.5(k) and 3.5(l) show the total surface waves N⃗ˆ(1) and the subsurface signals P⃗ˆ(1) after the first loop. P⃗ˆ(1) contains a surface-wave re-sidual at this stage, and becomes the rere-sidual update ( ⃗P + ∆ ⃗N )(2) for the second loop. For subsequent loops, the same procedure is iterated until it has reached a stopping criterion. Empirically, a few iterations are suﬃcient to converge to a stable solution.

Figure 3.6 shows the seismic data ( ⃗P + ⃗N ), the estimated surface waves ⃗ˆ

N(3) and the resulting estimated subsurface signals P⃗ˆ(3) after three itera-tions, in the xy-t and the kxky-f domain. Here, a moderate signal-protecting

scheme is used. As already stated, there is a trade-oﬀ between conservative signal protection and aggressive surface-wave estimation. Figures 3.7 and 3.8 show the results, each with a conservative and an aggressive signal protecting scheme. In the former case,P is protected, but less energy of⃗ˆ N is estimated.⃗ˆ

In the latter, more energy ofN is estimated, but⃗ˆ P is touched and distorted.⃗ˆ

This is observed especially in the kxky-f domain. The choice of the scheme

can be determined by the objective of surface-wave separation. IfN serves to⃗ˆ

obtain P for reservoir imaging and characterization, a conservative scheme⃗ˆ

should be used. On the contrary, an aggressive scheme would be preferred, ifN serves to use itself for near-surface characterization.⃗ˆ

(44)

Figure 3.5: Step-by-step illustration of the closed loop. (a) and (b) The initial dispersion

surface and the initial source wavelet of the fundamental mode; (c) and (d) the residual update in the pxpy-f domain without and with a smoothing/windowing filter; (e) and

(f) the estimated dispersion surface and the resulting fundamental mode update; (to be continued)

(45)

Figure 3.5: (Continued) (g) and (h) the estimated source update and the resulting

fundamental-mode update; (i) the fundamental-mode update after estimating the adaptive filter; (j) the higher-mode update after estimating the adaptive filter; (k) and (l) the estim-ated surface waves and the resulting subsurface signals after the first loop, with a vertical section at the top and a horizontal time/frequency slice at the bottom. A dotted line in the section indicates the position of the slice, and vice versa.

(46)

Figure 3.6: The results with a moderate signal protecting scheme. (a) and (d) The seismic

data; (b) and (e) the estimated surface waves; (c) and (f) the resulting subsurface signals, in the xy-t and the kxky-f domain with a vertical section at the top and a horizontal

time/frequency slice at the bottom. A dotted line in the section indicates the position of the slice, and vice versa.

(47)

Figure 3.7: The results with a conservative signal protecting scheme. (a) and (d) The

seismic data; (b) and (e) the estimated surface waves; (c) and (f) the resulting subsurface signals, in the xy-t and the kxky-f domain with a vertical section at the top and a horizontal

(48)

Figure 3.8: The results with an aggressive signal protecting scheme. (a) and (d) The

seismic data; (b) and (e) the estimated surface waves; (c) and (f) the resulting subsurface signals, in the xy-t and the kxky-f domain with a vertical section at the top and a horizontal

(49)

3.5. Post-SWES processing 41

3.5 Post-SWES processing

SWES estimates the surface-wave parameters, see e.g. Figure 3.4(b), and, consequently, estimates most of the surface waves including the aliased en-ergy. However, a residual may still exist in the form of unexplained noise. A part of the residual is present as remaining surface-wave samples ∆N in⃗˜

the pxpy-f domain. Optionally, to further minimize this residual, the surface

waves N can be further estimated by minimizing the objective function⃗ˆ J (∆N⃗˜n) = ∑ ω ⃗P + ∆ ⃗N 2 =∑ ω ( ⃗P + ⃗N )−N⃗ˆ 2 =∑ ω ( ⃗P + ⃗N )− nn ∑ n B−1LH∆N⃗˜n 2 , (3.14)

see also Equation (2.31). The model parameters ∆N⃗˜n for each surface-wave

mode are determined in such a way that the residual betweenN and ( ⃗⃗ˆ P + ⃗N )

(50)

is minimized. This minimization scheme is supposed to act on N only, by⃗ˆ

sparsity constraints implemented as amplitude thresholding and areal win-dowing to extract ∆N⃗˜n in the pxpy-f domain. Therefore, ⃗P remains

un-touched. Consequently, the resulting residual ( ⃗P + ∆ ⃗N ) approaches the

subsurface signals ⃗P . This approach can be viewed as an iterative filtering

in the pxpy-f domain. To solve the inverse problem, again, a closed-loop

approach is used. Figure 3.9 shows the closed-loop concept. Similar to the earlier closed-loop approach, each i-loop contains inner j-loop for each of the surface-wave modes. For each inner loop (i, j), the model parameters ∆N⃗˜_j(i) are estimated from the residual update ( ⃗P + ∆ ⃗N )(i,j), and the res-ulting modal surface-wave update ∆ ⃗N_j(i) is obtained. For a loop (i, j = n), ∆N⃗˜n(i) and ∆ ⃗Nn(i) are obtained by

∆N⃗˜_n(i)= C(i)_N

nLB( ⃗P + ∆ ⃗N )

(i,n)_, _(3.15)

∆ ⃗N_n(i)= B−1LH∆N⃗˜_n(i), (3.16)

where C(i)_N

n represents the sparsity constraints for each surface-wave mode. In a similar fashion as the earlier closed-loop approach, this update is summed in terms of j to build the intermediate and the total surface-wave update, re-spectively, i.e., ∆ ⃗N(i,j)=

nj ∑

j

γ_j(i)∆ ⃗N_j(i)and ∆ ⃗N(i)= ∆ ⃗N(i,j=nn)₌

nn ∑

j

γ_j(i)∆ ⃗N_j(i)

where γ_j(i) is the data-adaption scalar. This update is subsequently summed in terms of i to build the total surface-wave estimate, i.e., N⃗ˆ(i) =

ni ∑

i

∆ ⃗N(i), resulting in the surface-wave estimate. This algorithm can be expressed as

( ⃗P + ∆ ⃗N )(i,j)= ( ⃗P + ⃗N )−N⃗ˆ(i,j−1), (3.17)

⃗ˆ

N(i,j)=N⃗ˆ(i,j−1)+ γ_j(i)∆ ⃗N_j(i), (3.18)

⃗ˆ

P(i,j)= ( ⃗P + ⃗N )−N⃗ˆ(i,j), (3.19) where

⃗ˆ

N(i,0) =N⃗ˆ(i−1), (3.20) and the data-adaption scalar obtained by an algorithm to minimize the re-sidual, e.g., a steepest-descent code,

γ_j(i)= (( ⃗P + ∆ ⃗N ) (i,j)₎H_{∆ ⃗}_N(i) j 2(∆ ⃗N_j(i))H_{∆ ⃗}_N(i) j . (3.21)

Surface-wave separation and its impact on seismic survey design

Surface-wave separation and its

impact on seismic survey design

Surface-wave separation and its

impact on seismic survey design

Contents

Chapter 1

Introduction

1.1

Surface waves in seismic data

1.2

Objectives and outline

Chapter 2

Forward model of surface waves in

seismic data

1

2.1

Introduction

2.2

Surface-wave properties in seismic data

2.3

Seismic survey parameters and attributes

2.4

Forward model of seismic data

2.5

Forward model of surface waves

2.6

Conclusions from this chapter

2.7

Appendix: Radon transform

Chapter 3

Surface-wave estimation and

separation using closed loop

1

3.1

Introduction

3.2

Inverse problem

3.3

Closed loop

3.4

Step-by-step illustration of the closed loop

3.5

Post-SWES processing