Simultaneous Source Acquisition and Separation on General Related

Sampling Grids

Field of the invention

[0001] The present invention relates to methods for acquiring and separating contributions from two or more different simultaneously emitting sources in a common set of measured signals representing a wavefield where the sources traverse general, related sampling grids. In particular, the present invention relates to acquiring and separating contributions from two or more different simultaneously emitting seismic sources traversing general, related sampling grids and

separating sets of recorded and/or processed seismic signals.

Description of related art

[0002] A common and long-standing problem in physical

wavefield experimentation is how to separate recorded signals from two or more simultaneously emitting sources. In

particular, for more than a decade, the simultaneous source problem has (arguably) been the most pertinent problem to solve to efficiently acquire data for 3D reflection seismic imaging of complex Earth subsurface structures.

[0003] Simultaneously emitting sources, such that their signals overlap in the (seismic) record, is also known in the industry as "blending". Conversely, separating signals from two or more simultaneously emitting sources is also known as "deblending" . And the data from such acquisitions as "blended data". Halliday et al . (2014) describe methods for

simultaneous source acquisition relying on highly-controlled marine vibroseis records.

[0004] Modern digital data processing of wavefields (or signals) uses a discretized version of the original wavefield, say g, that is obtained by sampling g on a discrete set. The Nyquist-Shannon sampling theorem shows how g can be recovered from its samples; for an infinite number of equidistant samples and given sample rate k _s , perfect reconstruction is guaranteed provided that the underlying signal was bandlimited to \k\≤k _N = k _s /2 (Shannon, 1949; Nyquist, 1928), where k _N is the so-called Nyquist wavenumber. The Nyquist-Shannon sampling theorem is equally relevant both to signals generated from a single source being recorded on multiple receivers (receiver- side sampling) as well as signals generated from multiple sources and recorded at a single receiver (source-side

sampling) .

[0005] Assume that the wavefield g is measured at a specific recording location for a source that is excited at different source positions along a straight line. The sampling theorem then dictates how the source locations must be sampled for a given frequency of the source and phase velocity of the wavefield. One aspect of the sampling problem is as follows. Consider that instead of using one source, one wants to use two (or more) sources to for instance increase the rate at which data can be acquired. The second source is triggered simultaneously or close in time with the first source while moving along another arbitrarily oriented line to excite the wavefield h. At the recording location the wavefields

interfere and the sum of the two wavefields, f = g + h, is measured. There is no known published exact solution to

perfectly separate the wavefields g and h that were produced from each source from the combined measurement / (e.g., see Ikelle, 2010; Abma et al., 2015; Kumar et al, 2015; Mueller et al. , 2015) .

[0006] It may therefore be seen as an object of the invention, to present new and/or improved methods for acquiring and separating simultaneous-source data, particularly, including methods that accommodate and can exploit sources acquired on general but related sampling grids to increase the efficiency and/or the quality with which sources can be simultaneously acquired and separated. Brief summary of the invention

[0007] Methods for separating or deblending wavefields generated by two or more sources contributing to a common set of measured or recorded signals are provided, where the sources have traversed general, related sampling grids, suited for seismic applications and other purposes, substantially as shown in and/or described in connection with at least one of the figures, and as set forth more completely in the claims.

[0008] Advantages, aspects and novel features of the present invention, as well as details of an illustrated embodiment thereof, may be more fully understood from the following description and drawings .

Brief Description of the Drawings

[0009] In the following description reference is made to the attached figures, in which:

Figs. 1A, B illustrates how in a conventional marine seismic survey all signal energy of two source typically sits inside a "signal cone" (horizontally striped) bounded by the propagation velocity of the recording medium and how this energy can be split in a transform domain by applying a modulation to the second source ;

Fig. 2 is the function from Eq. (3) for a set of apparition blended synthetic data with noise. The left panel shows horizontally striped the region described by Eq. (1), while the horizontally striped part in the center of the right panel shows the diamond shaped conditions Eq. (4), and the aliasing/mixing effects from Eqs . (3) and (5) is highlighted by speckles/dots.

Fig. 3 shows the effect of the source ghost when acquiring a data with varying source depths can be thought of as ghost sources (open white stars) excited at mirror locations across the sea surface with

opposite polarity compared to the true sources (solid black stars), where the responses due to the open white stars represents the deghosted the data on the source side ;

Fig. 4 Top panel contains blended data in the Fourier domain; the lower panel the recovered four functions.

Fig. 5 Left panel shows the shot location grids for four simultaneous sources and the location of 64 point diffractors in a subsurface model. Right panel shows a zoom-in on the shot location grids around the origin. Two of the sources acquire their shot locations with one grid spacing while the other two sources acquire their sources with a different grid spacing.

Fig. 6 Shows the encoding time delays for the shot locations in the shot location grids of each of the four simultaneous sources of Fig. 5.

Fig. 7 Shows the regions of the inline-crossline (Kx- Ky) horizontal wavenumber domain that each of the four sources of Fig. 5 map into after encoding the shots with the the encoding time delays of Fig. 6 and using the different (smaller) scales (grid spacings) for sources 2 and 3 compared to the scales (grid spacings) used for sources 1 and 4.

Fig. 8 Shows three slices (top, middle, bottom) through the 3D volume of synthetic simultaneous source data modeled for the four sources and their shot locations in Fig. 5, encoded using the time delays of Fig. 6, simultaneously shooting into a single receiver at the origin above a model with 64 point scatterers (left column) . In the middle column, the same three slices of the separated, remodeled data are shown. In the right column, the residual (i.e., the unmodeled part of the total simultaneous source wavefield) is shown. Fig. 9 Shows three slices (top, middle, bottom) through the 3D volume of the frequency-inline-crossline (F-Kx- Ky) wavenumber amplitude spectrum of the synthetic simultaneous source data modeled for the four sources and their shot locations in Fig. 5, encoded using the time delays of Fig. 6, simultaneously shooting into a single receiver at the origin above a model with 64 point scatterers (left column) . In the middle column, the same slices are shown of the separated, remodeled data. In the right column, the residual (i.e., the unmodeled part of the total simultaneous source

wavefield) is shown.

Fig. 10 Shows three slices (top, middle, bottom) through the 3D volume of the frequency-inline-crossline (F-Kx-Ky) wavenumber amplitude spectrum of the

synthetic simultaneous source data modeled for the fourth source and its shot locations in Fig. 5, activated independently and without encoding time delays, shooting into a single receiver at the origin above a model with 64 point scatterers (left column) . In the middle column, the same slices are shown for the separated data. In the right column, the difference between the modeled (reference) and separated data for the fourth source is shown.

Fig. 11 Shows three slices (top, middle, bottom) through the 3D volume of the frequency-inline-crossline (F-Kx-Ky) wavenumber amplitude spectrum of the

synthetic simultaneous source data modeled for the third source and its shot locations in Fig. 5,

activated independently and without encoding time delays, shooting into a single receiver at the origin above a model with 64 point scatterers (left column) . In the middle column, the same slices are shown for the separated data. In the right column, the difference between the modeled (reference) and separated data for the third source is shown. Fig. 12 Shows three slices (top, middle, bottom) through the 3D volume of synthetic simultaneous source data modeled for the fourth source and its shot

locations in Fig. 5, activated independently and without encoding time delays, shooting into a single receiver at the origin above a model with 64 point scatterers (left column) . In the middle column, the same slices are shown for the separated data. In the right column, the difference between the modeled

(reference) and separated data for the fourth source is shown .

Fig. 13 Shows three slices (top, middle, bottom) through the 3D volume of synthetic simultaneous source data modeled for the third source and its shot

locations in Fig. 5, activated independently and without encoding time delays, shooting into a single receiver at the origin above a model with 64 point scatterers (left column) . In the middle column, the same slices are shown for the separated data. In the right column, the difference between the modeled

(reference) and separated data for the third source is shown .

Fig. 14 Shows synthetic ghosted simultaneous source data modeled for a simple three-reflector model, with the sources shooting into a single receiver and doubly encoded and the results after decoding the first encoding (2 ^nd and 3 ^rd panel from the left) , resorting (4 ^th panel from the left) and the deghosted result after decoding the second encoding (5 ^th panel from the left) .

Fig 15. Discloses methods for towing of source arrays at different tow depths but offset to each other.

Fig 16. Discloses various configurations for towing one or more source arrays at varying or alternating tow depths . Fig 17. Discloses sub-arrays of a seismic source array that are towed at different depth levels encoded using seismic apparition principles.

Detailed Description

[0010] The following examples may be better understood using a theoretical overview as presented below.

[0011] The slowest observable (apparent) velocity of a signal along a line of recordings in any kind of wave experimentation is identical to the slowest physical propagation velocity in the medium where the recordings are made. As a result, after a spatial and temporal Fourier transform, large parts of the frequency-wavenumber ( a)k ) spectrum inside the Nyquist

frequency and wavenumber tend to be empty. In particular, for marine reflection seismic data (Robertsson et al . , 2015), the slowest observable velocity of arrivals corresponds to the propagation velocity in water (around 1500m/s) .

[0012] Fig. 1(A) illustrates how all signal energy when represented in or transformed into the frequency-wavenumber ( a)k ) domain sits inside a "signal cone" centered at k = 0 and bounded by the propagation velocity of the recording medium.

[0013] It is well known, for example, that due to the

"uncertainty principle", a function and its Fourier transform cannot both have bounded support. As (seismic) data are necessarily acquired over a finite spatial (and temporal) extent, the terms "bounded support" and "limited support" herein are used not in the strict mathematical sense, but rather to describe an "effective numerical support", that can be characterised, e.g., by the (amplitude) spectrum being larger than a certain value. For instance, larger than a certain noise threshold, or larger than the quantization error of the analog-to-digital converters used in the measurement equipment. Further, it is understood that by explicitly windowing space and/or space-time domain data, the support of a function may be spread over a larger region of, e.g., the wavenumber-frequency domain and in such cases the term

"bounded support" and "limited support" will also be

understood as "effective numerical support" as it will still be possible to apply the methods described herein.

[0014] Furthermore, the terms "cone" and "cone-shaped" used herein are used to indicate the shape of the "bounded" or "effective numerical" support of the data of interest (e.g., the data that would be recorded firing the sources

individually [i.e. non-simultaneously] ) in the frequency- wavenumber domain. In many cases, it will still be possible to apply the methods described herein if the actual support is approximately conic or approximately cone-shaped. For example, at certain frequencies or across certain frequency ranges the support could be locally wider or less wide than strictly defined by a cone. Such variations are contemplated and within the scope of the appended claims. That is, the terms "cone" and "cone-shaped" should be understood to include

approximately conic and approximately cone-shaped. In

addition, in some cases we use the terms "bounded support" or "limited support" and "effective numerical support" to refer to data with "conic support" or "cone-shaped support" even though in the strict mathematical sense a "cone" is not bounded (as it extends to infinite temporal frequency) . In such cases, the "boundedness" should be understood to refer to the support of the data along the wavenumber axis/axes, whereas "conic" refers to the overall shape of the support in the frequency-wavenumber domain.

[0015] Note that the term "cone-shaped support" or similar refers to the shape of the support of e.g. the data of

interest (in the frequency-wavenumber domain), if it were regularly sampled along a linear trajectory in 2D or Cartesian grid in 3D. That is, it refers only to the existence of such a support and not to the actual observed support of the data of interest in the simultaneous source input data or of the separated data of interest sampled as desired. The support of both of these depends on the chosen regularly or irregularly sampled straight or curved input (activation) and output

(separation) lines or grids. Such variations are within the scope of the appended claims.

[0016] For example consider a case where the input data are acquired using simultaneous curved shot lines. In this case, the methods described herein can either be applied directly to the input data, provided the curvature has not widened the support of the data interest such that it significantly overlaps with itself. In this case, the support used in the methods described herein can be different from cone-shaped. Alternatively, the methods described herein are used to reconstruct the data of interest in a transform domain which corresponds to, e.g., best-fitting regularly sampled and/or straight activation lines or Cartesian grids, followed by computing the separated data of interest in the non- transformed domain at desired regular or irregularly sampled locations .

[0017] In a wavefield experiment it may be that a source is excited sequentially for multiple source locations along a line while recording the reflected wavefield on at least one receiver. The source may be characterized by its temporal signature. In the conventional way of acquiring signals representing a wavefield the source may be excited using the same signature from source location to source location, denoted by integer n. Next, consider the alternative way of acquiring such a line of data using a periodic sequence of source signatures: every second source may have a constant signature and every other second source may have a signature which can for example be a scaled or filtered function of the first source signature. Let this scaling or convolution filter be denoted by a(t), with frequency-domain transform Α(ω) .

Analyzed in the frequency domain, using for example a receiver gather (one receiver station measuring the response from a sequence of sources) recorded in this way, can be constructed from the following modulating function m(n) applied to a conventionally sampled and recorded set of wavefield signals:

which can also be written as

[0018] By applying the function m in Eq. 0.1 as a modulating function to data fn) before taking a discrete Fourier

transform in space (over the following result

can be obtained:

which follows from a standard Fourier transform result

(wavenumber shift) (Bracewell, 1999) .

[0019] Eq. 0.2 shows that the recorded data / will be scaled and replicated into two places in the spectral domain as illustrated in Fig. 1(B) and as quantified in Tab. I for different choices of

TAB. I. Mapping of signal to cone centered at and cone centered at for different choices of Α for signal separation or signal apparition in Eq. (0.2) . [0020] Part of the data will remain at the signal cone

centered around k = 0 (denoted by H ₊ in Fig. 1(b)) and part of the data will be scaled and replicated to a signal cone centered around k _N (denoted by H_) . It can be observed that by only knowing one of these parts of the data it is possible to predict the other.

[0021] This process may be referred to as "wavefield

apparition" or "signal apparition" in the meaning of "the act of becoming visible". In the spectral domain, the wavefield caused by the periodic source sequence is nearly "ghostly apparent" and isolated.

[0022] A particular application of interest that can be solved by using the result in Eq. (0.2) is that of simultaneous source separation. Assume that a first source with constant signature is moved along an essentially straight line with uniform sampling of the source locations where it generates the wavefield g . Along another essentially straight line a second source is also moved with uniform sampling. Its signature is varied for every second source location according to the simple deterministic modulating sequence m(n) ,

generating the wavefield h. The summed, interfering data /= g + h are recorded at a receiver location.

[0023] In the frequency-wavenumber domain, where the recorded data are denoted by is partitioned into two components H ₊ and H_ with H = H ₊ + H_ where the //.-component is nearly "ghostly apparent" and isolated around the Nyquist- wavenumber [Fig. 1(B)], whereas G and H ₊ are overlapping wavefields around Furthermore,

H_ is a known, scaled function of H. The scaling depends on the chosen Α(ω) function (Tab. I), and can be deterministically removed, thereby producing the full appearance of the transformed wavefield H. When H is found, then G = F— H yielding the separate wavefields g and h in the time-space domain. [0024] Although the above description has focused on acquisition along essentially straight lines, the methodology applies equally well to curved trajectories such as coil- shaped trajectories, circles, or other smoothly varying trajectories or sequences of source activations.

[0025] The concept may be extended to the simultaneous acquisition of more than two source lines by choosing

different modulation functions for each source.

[0026] Acquiring a source line where the first two source locations have the same signature, followed by two again with the same signature but modified from the previous two by the function Α(ω) and then repeating the pattern again until the full source line has been acquired, will generate additional signal cones centered around

[0027] Fig. 1(B) also illustrates a possible limitation of signal apparition. The and parts are separated within

the respective lozenge-shaped regions in Fig. 1(B) . In the triangle-shaped parts they interfere and may no longer be separately predicted without further assumptions. In the example shown in Fig. 1(B), it can therefore be noted that the maximum non-aliased frequency for a certain spatial sampling is reduced by a factor of two after applying signal

apparition. Assuming that data are adequately sampled, the method nevertheless enables full separation of data recorded in wavefield experimentation where two source lines are acquired simultaneously.

[0028] As is evident from Tab. I, the special case A = l corresponds to regular acquisition and thus produces no signal apparition. Obviously, it is advantageous to choose A

significantly different from unity so that signal apparition becomes significant and above noise levels. The case where A = —1 (acquisition of data where the source signature flips polarity between source locations) may appear to be the optimal choice as it fully shifts all energy from k = 0 to k _N (Bracewell, 1999) . Although this is a valid choice for modeling, it is not practical for many applications (e.g., for marine air gun sources, see Robertsson et al . , 2015 as it requires the ability to flip polarity of the source signal. The case where A = 0 (source excited every second time only) may be a straightforward way to acquire simultaneous source data but has the limitation of reduced sub-surface

illumination. A particularly attractive choice of for wave experimentation seems to let every second source be excited a time shift T later compared to neighbouring

recordings, that is, select

[0029] The above description assumes a simple modulating sequence m(n), and thus generating the wavefield h . In practice it is difficult to obtain perfectly periodic time shifts from a measurement setup. It is for example common practice for seismic vessels to shoot or trigger their sources at

predetermined (essentially equidistant) positions, and due to practical variations (vessel velocity etc.) it will be

difficult to realize shots at both predetermined locations and times .

[0030] Deviations from perfectly periodic acquisition can be termed non-periodic and grouped into non-periodic controlled (or intentional) and non-periodic uncontrolled cases (such as caused by currents, rough seas, etc., which are beyond

influence by the acquisition crew) . Furthermore, non-periodic acquisition can be composed of a periodic part, overlain by a non-periodic part. In all these cases, the signal cone will be scaled and replicated additional times along the wavenumber axis and the effects can be dealt with by various methods, including modelling and inverting such scaled replications using cyclic convolution functions as described in more detail later .

[0031] Note that periodic or aperiodic variations in source locations can similarly be used to apparate the wavefield signals. This can be understood for example by noting that a variation in the source location produces (angle-dependent) time shifts and therefore can be used to encode the data using the apparition methods described above.

[0032] For a sub-horizontally layered Earth, the recorded reflections from the interfaces between the strata lie

(approximately) on hyperbolic trajectories in the space-time domain. The change in two-way traveltime of such reflections as a function of the source-receiver distance (or offset) is known as the normal moveout (NMO) and depends on the zero- offset two-way traveltime and the corresponding average sound speed in the overlying strata.

[0033] Correction of the normal moveout (NMO correction) is a standard procedure in seismic data processing which aims to remove the offset dependent part of the traveltime and align the reflected arrivals according to their zero-offset

traveltime such that they can be summed yielding an initial "stack image" of the subsurface with increased signal-to-noise ratio .

[0034] NMO correction is a very efficient way to reduce the maximum time-dip in the recorded data. On the other hand NMO correction tends to stretch the data at large offsets and at early times, effectively changing (lowering) the frequency content in a space- and time-dependent manner. Let us consider the effect of NMO correction on simultaneous source data that have been acquired using e.g. seismic apparition, or similar, principles .

[0035] Because of the stretch, it follows that the NMO correction also modifies the apparition encoding filters a(t) in an offset- and time-dependent manner. However, note that such effects can be accurately predicted or modelled either from theory and first principles and/or numerical experiments. For example, if the encoding filters used were pure time delays, then the time delay after NMO correction can be predicted accurately by multiplying with an expression for NMO stretch due to Barnes (1992) :

where

[0036] Alternatively, the space-time dependent effect of the NMO correction on encoding filters can be considered by evaluating the effect of NMO correction at t ₀ on a discrete delta function and on , respectively, and

computing, e.g., the ratio of the resulting responses in the frequency domain. This yields a time- and offset-dependent frequency filter which can be used to predict the effective modulation function (also time- and offset dependent in general) after NMO correction.

[0037] Thus, an effective modulation function takes into account, e.g., the space-time dependent effects of the NMO correction, or any other coordinate transform, on the encoding filters .

[0038] The well-known convolution theorem states that

convolution in the time or space domain corresponds to

multiplication in the temporal frequency or spatial frequency domain, respectively. The lesser-known dual of the convolution theorem states that multiplication in the space domain of d with a so-called modulation operator corresponds to

cyclic convolution of the (discrete) Fourier transform of the data, with the (discrete) Fourier transform of the

modulation operator , followed by inverse

(discrete) Fourier transform. Further, we note that cyclic convolution can be implemented conveniently as a matrix multiplication or using computationally fast methods such as the fast Fourier transform (FFT) .

[0039] Thus, for general aperiodic modulation functions, the recorded simultaneous source data can be modelled in the frequency-wavenumber domain as the sum of the fk-domain wavefields due to the individual sources, multiplied by one or more corresponding cyclic convolution matrices. Then, the fk- domain wavefields due to the individual sources can be

obtained by inverting the data through the model. Note that in this context, here and elsewhere, setting up and solving a system of equations can equally be understood as modelling and inversion, respectively.

[0040] Note that the effect of a general aperiodic modulation as compared to a periodic modulation can thus be understood as introducing additional, scaled replications (beyond the replications at (multiples of) +/- Nyquist horizontal

wavenumber) of the individual signal cones of the sources, which describe the known to be compact support of the sources, along the wavenumber axis/axes. Both the position and the scaling factor of the replications then are exactly given by the (discrete) Fourier transform of the aperiodic modulation function .

[0041] Finally, the S-transform (Stockwell, 1996) decomposes a time-signal into a time-frequency representation, localizing signals in both time and frequency. It provides a frequency- dependent resolution in accordance with the uncertainty principle while maintaining a direct relationship with the Fourier spectrum.

[0042] It is possible, then, to use the cyclic convolution principle in conjunction with the S-transform (or similar time-frequency decomposition) and NMO correction to improve the separation of aliased simultaneous source data, acquired, e.g., using seismic apparition principles, in the following manner :

1. NMO correct the (partially) aliased simultaneous source recordings using an invertible NMO transform, effectively reducing the maximum time-dip in the data.

2. Compute the S-transform of each NMO corrected

simultaneous source trace. The resulting time-frequency representation implies a frequency spectrum for every time sample of every trace. 3. Iterate over all time samples of the S-transformed NMO- corrected sim source traces. For every time sample:

a. Compute the offset (and current time) dependent effect of NMO stretch on the encoding filters using, e.g., Eq. 0.3 from Barnes (1992), providing the frequency spectrum of the effective encoding filters .

b. Loop over all the (positive/negative) frequencies.

For every frequency do:

i. Compute the wavenumber transform of the S- transformed NMO-corrected simultaneous source data at the current time sample and frequency. ii. Compute the wavenumber transform of the

effective modulation function,

iii. Using the cyclic convolution principle and the computed wavenumber transform of the effective modulation function, compute a forward

modelling operator that relates the source one and source two contributions at the current time sample and frequency to the corresponding sim source data.

iv. Invert the sim source data at the current time sample and current frequency using the forward modelling operator,

v. Compute the inverse wavenumber transform of the inverted source one and source two contributions and store the results.

c. Optional (if all the frequencies are computed) :

Compute the negative frequencies from the positive frequencies using complex conjugate mirroring.

d. Store the separated frequency spectra for the current time-step building the S-transforms of the separated data one time sample by one.

4. Inverse S-transform the obtained S-transform data for

source one and source two, obtaining the NMO corrected separated source one and source two data.

5. Inverse NMO correct separated si and s2 data.

[0043] The list of steps is merely included for completeness of the description of a method which improves the separation of aliased simultaneous source data.

[0044] The NMO is a coordinate transformation that reduces the spatial bandwidth of the recorded data, and therefore limiting the effect of aliasing. We now proceed to discuss methods that use other coordinate transformations, and also how several coordinate transformations can be used simultaneously.

Moreover, we also discuss how to make reconstruction in two steps: First by making partial reconstructions, using only the non-aliased parts; and secondly to use these partial

reconstructions to regularize and solve the full

reconstruction problem by means of directionality estimates, that imply local coordinate transformations specifying

directions with reduced bandwidth, and hence, reduced aliasing effects .

[0045] Further, to provide a more completele summary of methods for dealing with aliased simultaneous source data, we review the notation and recapitulate the theory for regular seismic apparition. We use the notation

for the Fourier transform in one variable, and consequently for the Fourier transform of two dimensional function with a time (t) and spatial (x) dependence.

[0046] Suppose that are two function

with frequency support in two cones of the form

[0047] The constraint comes from assuming that the functions represent the recording of a wavefield at time t at a fixed receiver coordinate, source coordinate x, and fixed depth, where the recorded wave field at this depth is a solution to the homogeneous wave equation with a velocity c . The wavefields are generated at discrete values of x which we assume to be equally spaced, i.e. of the form

[0048] We now assume that the two sources are used

simultaneously, in such a way that their mixture takes the form

i.e., the recorded data is now modelled as a sum of the two functions, but where one of them has been subjected to a periodic time shift. In a more general version more general filtering operations than time shifts can be applied. Let a _k filter operators (acting on the time variable) where the k dependence is such it only depends on if k is odd or even, i.e., that

[0049] It can be shown that

[0050] Now, due to the assumption of conic support of it holds that if

then only the terms where k = 0 above contribute, and the following simplified relation holds

[0051] In a similar fashion it holds for

[0052] This implies that for each pair satisfying (4), the values of can be obtained by solving the linear system of equations

[0053] This provides information on how to recover the wavefields for frequencies either up to the limit

or more generally, satisfying the (diamond shaped)

condition (4) . The overlaps of the cones are illustrated in Figure 2, where the left panel shows horizontally striped the cone given by (1), and the right panel the diamond shaped region described by (4) . In this approach the deblending takes place by considering the data available in the central cone of Figure 2.

[0054] An alternative approach for reconstruction, is by noting that if either of the support constraints (1) or (4) are satisfied, then it holds that for the values of of

interest that (3) reduces for instance to

implying that can be recovered from

In a similar fashion, can be recovered from

In this way, the deblending can be achieved by direct

consideration of the data in the shifted cones illustrated in Figure 2.

[0055] From (6) it is possible to recover the functions and f ₂ partially. Let be a filter such that has support inside

the region described by (4) . It is then possible to recover

[0056] For values of outside the region described by it is not possible to determine uniquely

without imposing additional constraints. Typically, seismic data can locally be well described by sums of plane waves with different directions. The plane waves carry the imprint of the source wavelet, and according to ray theory the data from such a plane event should have the same directionality for the frequency range that covers the source wavelet. We can use this information to construct a directionality penalty that we can use for the separation of the two wavefields f _x and f ₂ from the blended data d . This directionality penalty is equivalent to, by means of local coordinate transformations, imposing a bandwidth limitation in undesired directions for the purpose of suppressing aliasing.

[0057] One way of estimating local directionality is by means of so-called structure tensors. For the two known wavefields and the corresponding structure tensors are defined as

and similarly for T ₂ and h ₂ . Above, the function K describes a smooth, localizing windows, for instance a Gaussian. The eigenvalues of T _x and T ₂ will point-wise describe the local energy in the direction of maximum and minimum variation, and the associated eigenvectors contain the corresponding

directions. The tensors are computed as elementwise

convolutions of the outer product of the gradient of the underlying function, and this directly defines the

generalization to higher dimensions. For the sake of

simplicity, we describe here the two-dimensional case.

[0058] Let sj(t, x) and s ₂ (t, x) be the eigenvalues of Τ^ί, χ) , and let e¾;(t, x) and e ₂ (t, x) denote the corresponding eigenvectors. If the wavefield f _x only has energy in one direction in the vicinity around (t, x) covered by K, then this implies that

which in turn means that

[0059] The eigenvectors define local

coordinate transformation that describe directions of large and small variations. Along the directions specified by

only low-frequency components are to be dominant, and by suppressing the bandwidth of the reconstructions in these directions is an efficient way of de-aliasing the separated sources .

[0060] This property (8) is clearly not always satisfied (although counterparts in higher dimension hold more

frequently with increased dimensionality) , however it is a property that can be used as a penalty from which the blended data can be deblended. Even if (8) is not satisfied, the relation can be used to minimize the energy of the deblended data in the directions carried from and h ₂ , respectively.

[0061] From (8) we have a condition on the gradient of f _x and f ₂ when one the eigenvectors vanishes. For the more general case, we would need to formulate a penalty function that can deal with the cases where the components change gradually, and at places where the eigenvectors are equal in size, and equal amount of penalty should be used for the two directions. One such choice is to define

with

These functions have the property th and

implying that (8) will be satisfied in the case where there is locally only energy in one direction, and where an equal amount of penalty will be applied in the case where there is the same amount of energy in both directions. Note that the local coordinate transformations are now implicitly given in the operator S.

[0062] This definition now allows for the generalization of (8) to penalty functionals

and

for the two wavefields. The expressions above describe the energy in the undesirable direction, given the knowledge of the bandpass filtered versions respectively. The de-

aliasing is now taken place by punishing high frequencies (by the derivatives) along the directions given by the local coordinate transformations specified by and

[0063] Before we use these expressions to define a

minimization problem that describes the deblending procedure, we incorporate the original cone condition (1) in the

formulation. To this end, we will now work with sampled representations of and In the following, redefining the

notation, we will also use and to denote these sampled

values .

[0064] We let denote the inverse Fourier operator that is restricted to functions supported on the cone defined by (1) . Recall the definition of the apparition operator T from (2) . The relationship (2) is then satisfied for (the non-unique) solutions to

with the additional constraint that and have support on

the cone defined by (1) . To obtain a unique approximate solution, we now add the directionality penalties and consider

with the same cone constraint. To find the minima of (9), we compute the Frechet derivatives of the objective function (9) with respect to the functions and and equate them to zero

as they should at a minimum. The first term in (9) is

straightforward to derive, and concerning the other two terms it is readily verified using partial integrations that their Frechet derivatives are described by the elliptic operators

[0065] To formulate the solution to (9), let

Furthermore, introduce

and

Equating the Frechet derivatives of (9) with respect to and f ₂ to zero then yield the linear relationship for the solution of (9) . This equation can be solved using an iterative solver for linear equations, for instance the conjugate gradient method. The operators in A

are realized using standard FFT, and the operators in are computed using a combination of Fourier transforms and differential schemes, that also may be implemented by using FFT. The operator describe the fit to data, while the operator describe the de-aliasing that takes places using the local coordinate transformations induced from

De-aliasing using coordinate transformations Mutual parabolic change of coordinates

[0066] As seen, for the direct deblending using seismic apparition to work it is required that the energy is contained in either a combination of the cone constraint (1) along with the condition or in the larger domain determined

by (4) . Given the characteristics of seismic data,

improvements to these constraints can sometime be obtained by changes of coordinates. The simplest case where this can be utilized is by using a parabolic change of coordinates. This would apply in cases where the two source locations are close to each other. Let describe the coordinate map

for some choice of parameter q. Now, let

and let

Clearly, it holds that

It then holds that

and the regular apparition technique can now be directly applied to to deblend and hence also and

Separate coordinate systems and NMO

[0067] Next, still summarizing methods that can improve the separation of aliased simultaneous source data, we consider the situation where aliasing effects can be reduced for f ₁ and f ₂ when exposed to two different sets of coordinate

transformation. Let

be continuous bijective coordinate transformations, and let

[0068] The purpose of the coordinate transformations is to make the effective support of be smaller than and or at least that the energy content that is contained in regions of the form (4) relative to the total energy is higher for the functions in comparison to that of

In this way the aliasing effects when making reconstructions with limited frequency support hopefully be reduced.

[0069] However, the reduction of the aliasing effect for, say the first source, may have a negative influence on the

aliasing of the second source when applying the coordinate transformations to blended data. In many cases, the data cannot be separated directly, but the effects on both the sources need to be taken into account. The interplay between the two sources and the coordinate transformations can be modelled by means of linear systems, and by solving these systems we can obtain de-aliased deblending.

[0070] Recall that in the regular apparition setup, the functions f _x and f ₂ could be recovered using either the central frequency cone or one of the aliased versions. In the case where the central cone was used, a ( ) system of equations has to be solved for the recovery of the point-wise values of while when considering the shifted cones, the values of were obtained directly by a

correction factor. In the presence of the coordinate

transformations, neither one of the solutions will be directly applicable. However, it will still be the case that most of the energy (excluding induced aliasing effects from the opposite coordinate map) in the shifted cones will be due to only one of the two sources. We will make use of this fact, in combination with the fact that the sum of the two

representations should equal the measured data, and in this way set up a linear system of equations for the solution of the de-aliased-deblending problem.

[0071] To simplify the presentation, let us assume in this section that the filter operators are identity operators. Let us introduce

and

[0072] By assuming sufficiently dense sampling in time, and d ₂ can be computed by resampling of the measured data d. We now seek representations of using the (non-aliased) Fourier representations These representations read

and

respectively. Suppose now that τ is sampled discretely at equally spaced points, specifically

with a sampling in the corresponding frequency variable

where be an equally spaced sampling

of the spatial frequency parameter

[0073] Let us introduce the coordinates

[0074] We seek discrete Fourier representations of and respectively. For these representations, we want to include only the frequencies that satisfy the condition (4) or something similar. To this end, we introduce the mapping

that is indexing the points i.e., such that all the points satisfy a condition similar to (4) .

[0075] The Fourier representations can then be expressed using the element-wise defined Fourier matrices

[0076] Again, redefining the notation, we now use

also for the discrete representations, and similarly for d and d ₂ . We then have that

and [0077] The next step is now to provide a formulation where the contribution from the two sources can be separated. As a first step, we will make use of the fact that the energy present in the shifted Fourier cones of corresponds primarily of

energy from and respectively. Defining element-wise

by

[0078] Recall that a frequency shift in representations using the discrete Fourier transform can be represented by a

multiplication by Hence, define

as the operator that is shifting information to the cones of interest. This implies that the is almost in the null-space

of the mapping

and similarly, is almost in the null-space of the mapping

or in other words, the mapping

is almost block-diagonal. The matrix describing the map above will be square, and hence an iterative method such as the generalized minimal residual method (GMRES) could be applied to this problem, given the right hand side

[0079] However, using this information by itself means some loss of control in the residual to the data fit. To counteract that effect, we may also include an error term for the fit to the measured data. An overdetermined system for this case then reads

[0080] Let us end by discussing how to choose the coordinate transformation when dealing with NMO for different

parameterizations . Let us define the Heaviside function by

and define Let us define the NMO mapping

by

This mapping is bijective with inverse

coordinate transformation is chosen such that

is continuous.

[0081] Note that the presented methods of using coordinate transformations could favourably be applied simultaneously, where expected known aliasing effects can be reduced using for instance the NMO mappings mentioned above, along with a data driven de-aliasing part that uses the local coordinate

transformations for the directionality penalties.

[0082] In yet another example of a method to improve the dealiasing of simultaneous source data, if we apply an NMO correction to acquired simultaneous source data acquired using seismic apparition principles, one of the sources, say source A, will only produce a signal cone at k = 0 (i.e., at zero horizontal wavenumber) but that now is more compactly

supported as a result of the NMO correction. However, the other source, say source B, will produce two more compactly supported signal cones centered at k = 0 and k = k _n . We can therefore isolate the signal cone at k = k _n (i.e., at Nyquist horizontal wavenumber) and apply inverse NMO correction to achieve approximate separation of the wavefield due to source A. Such approximate separation methods could yield acceptable results, particularly when combined with other methods such as least-squares adaptive subtraction (Dragoset, 1995) .

[0083] Alternatively, one could use an apparition function that does not involve a time shift (e.g., A=0.5) . In this case, the apparition encoding function will be "immune" to the stretch effect of NMO correction such that the described non time-variant separation methods can be used. Yet another option to continue using the described non time-variant separation methods would be to use stretch-free NMO operators which have been proposed previously (e.g., Trickett, 2003;

Zhang et al . , 2012) .

[0084] Another method to handle aliased data when using, e.g., time-shifts as apparition encoding functions is to use NMO and create a look-up table of many solutions for different time shifts. We then piece together solutions sample by sample by choosing the appropriate sample in time and space for the time shift for a specific sample taking into account both the original time shift and the distortion introduced by NMO correction .

[0085] We note that (source) ghost reflections may also be considered simultaneous sources and, e.g., the methods

summarized above for improved separation in the presence of aliasing also apply to the problem of (source) deghosting by simultaneous source deblending in the presence of aliasing.

[0086] To understand source deghosting by simultaneous source deblending, consider the following example: a line of ghost sources with alternating depths such as shown in Figure 3 (open white stars) can be thought of as being constructed from ghost sources at a uniform distance from the original sea surface location the level closest to the original sea surface location) through the convolution of a modulation function :

[0087] where the lower case variables denote expressions in the time space domain, n _x is the shot number (spaced uniformly along a line) in the x-direction and n _y is the number of the parallel shot-line and is the modulating function:

[0088] The function a is a redatuming operator which is both a temporal and spatial filter [hence the convolution in space in equation (13)], which can be written in the frequency- wavenumber domain as : where The action of a can then be obtained from the

relation

[0089] The above equation (15) may be considered as a general example of a possible mathematical description of a redatuming operator, here in the frequency-wavenumber domain. In the frequency-wavenumber domain we obtain the following expression for the ghost wavefield using capital letters for the

variables in the frequency-wavenumber domain:

which follows from a standard Fourier transform result (wavenumber shift; Bracewell, 1999) .

[0090] Equation (16) shows that data from (staggered) ghost sources will be scaled and replicated or partitioned

into two places in the spectral domain. Part of the data will remain in a cone centred around with the limits of the

cone defined by the slowest propagation velocity of the wavefield in the medium and part of the data will be scaled and replicated or partitioned into a signal cone centered around along the with denoting the Nyquist

wavenumber representing the sampling frequency of the data. In other words, part of the data will remain at the signal cone centered around (denoted by in Fig. 1(B)) and

part of the data will be scaled and replicated to a signal cone centered around (denoted by

[0091] The carpet of desired (non-ghost) sources (solid black stars in Fig. 12 shown along a profile in the x-direction) with alternating depths can also be constructed from desired sources at a uniform depth z ₁ (the level closest to the

original sea surface location) through the convolution of another modulation function: where is the modulating function:

[0092] Again, the function b is a redatuming operator in the opposite direction compared to a such that: [0093] In the frequency-wavenumber domain we therefore obtain the following expression:

[0094] Again, equation (20) shows that the data from the

(staggered) desired sources will be scaled and

replicated or partitioned into two places in the spectral domain. Part of the data will remain at the signal cone centred around and part of the data will be scaled and

replicated or partitioned to a signal cone centered around

along the In other words, part of the data will remain at the signal cone centered around

(denoted by in Fig. 1(B)) and part of the data will be

scaled and replicated to a signal cone centered around

(denoted by

[0095] One method of source deghosting uses equations (16) and (20) . Equations (16) and (20) tell what the mix is of desired and ghost sources that occur around k = 0 (in the following denoted by F° ) :

[0096] In addition, equations (16) and (20) also tell what the mix of desired and ghost sources is that occur around k _N along the is (in the following denoted by

[0097] Equations (21) and (22) can be combined to give an expression for the wavefield of interest emitted from the desired source (i.e., the source-side deghosted wavefield) :

[0098] Since the operator A is known from equation (15), the deghosted wavefield (i.e., the wavefield of interest) can be computed explicitly. [0099] Note that A can be determined accurately for a 3D acquisition geometry. If the cross-line spacing is coarse or if only a single line of data is acquired, it may be necessary to resort two a 2D approximation of A.

[00100] Note that in this method for source deghosting we have not included any information about the reflection coefficient of the surface.

[00101] We can also make use of one more equation that relates the ghost sources to the desired sources:

where C is a redatuming operator that depends on the depth of the shallow sources z ₁ but that is related to the above defined operators A and B and knowledge or assumptions concerning the reflection coefficient of the free surface reflector, e.g. the sea surface, and is -1 in equation (24) :

[00102] Equations (16), (20), (24) and (25) therefore allow us to isolate the wavefields due to the (virtual) ghost sources and or the desired (physical) sources separately without knowing the redatuming operators A, B or C. In fact, the system of equations also allows to solve for the redatuming operator itself. Knowledge of such redatuming operators is important in many applications in seismic data processing such as in imaging of seismic data or in removing multiples in the data for instance. In this case, we obtain the following expression for the mix of desired and ghost sources that occur around (again denoted by

[00103] The following expression for the mix of desired and ghost sources that occur around k _N along the /c _x -axis (again denoted by

[00104] Equations (26) and (27) can be solved explicitly or iteratively for arbitrary values of z ₁ and sources Δζ such that we no longer need any information about what the redatuming operator looks like and instead can solve for the desired wavefield due to the physical sources only, and or the

wavefield due to the (virtual) ghost sources and or the redatuming operator.

[00105] Also in the source deghosting methods summarized above, for the direct deghosting by deblending to work it is required that the energy is contained in either a combination of a cone constraint such as (1) along with a condition such as

or in a larger domain such as determined by (4) .

Again, improvements to these constraints can be obtained by various changes of coordinates and considering their effect on the deblending as shown above, and, in some cases, even by using NMO or parabolic moveout corrections without considering their effect on the deblending, but instead using the

unmodified deblending expressions. Note that the cone-shaped supports or constraints and the related diamond-shaped

constraints in the frequency-wavenumber domain take on other shapes and forms in other transform domains and it should be understood when we refer to the support or the constraints as being so shaped in the fk-domain that this includes using the corresponding other shapes and forms in the other well-known transform domains as generally used in seismic data

processing .

[00106] Finally, note that in the following we also refer to virtual ghost sources as effective simultaneous sources, as opposed to real simultaneous sources which we reserve for separate physical sources.

Synchronous shooting on general related sampling grids

[00107] Following this summary of methods enabling acquisition and processing of perturbed and/or aliased simultaneous source or ghosted data, we present then, the new and/or improved methods for acquiring and separating simultaneous-source data acquired with general, but related sampling grids which increase the efficiency and/or the quality with which sources can be simultaneously acquired and separated. Key observations shaping the present invention are that: (1) no information about the relative source locations is required for encoding and decoding using the apparition principles described above, (2) no information about the source spatial sampling interval is required, other than for defining the support of the wave fields due to the individual sources relative to the Nyquist wavenumber (effective Nyquist wavenumber, in the case of irregular sampling), (3) the order in which the shot positions are traversed and activated is immaterial as long as all the sources traverse and activate at the shot positions in the same order, i.e., they acquire the source locations

synchronously, and (4) that Fourier transforms for data sampled at regular locations are just one type of transform that possess shift (for periodic modulation) and cyclic convolution (for general aperiodic modulation) properties.

[00108] A Bravais grid in two dimensions is a set of points of the form

where A is a 2x2 matrix. Associated with this grid is the dual grid

It can be verified that

[00109] Important examples include the regular quadratic grid for which

and the hexagon (rhombic) rid generated from

[00110] We will describe how the apparition principle works for Bravais grids. First, note that from the standard rules of the Fourier transform it holds that

Recall the two-dimensional Poisson sum formula

[00111] Using the Poisson sum formula in combination with the Fourier relation above, we have that

Now, let Each such element can thus be indexed as and respectively. Moreover, choose

for are non- negative integers. The Poisson sum can then be written

[00112] We now introduce the notation

which allows us to rewrite the above relation as

[00113] Note that if has such compact support that

if

then the left hand sum above reduces to a regular two- dimensional discrete Fourier transform which can be rapidly evaluated using FFT .

[00114] We now turn our attention to the apparition setup. Let and let be equal to the

number of elements in We will also use the notation modp to denote the modulus operator with respect to the set

[00115] Suppose that measurements are made using the

functions Morever, let h _p filter functions, for

instance time shift operators. We now assume that combinations of the p functions are measured at grid points

At each such point, linear combinations of the form

are measured during the acquisition. A Fourier transformation with respect to the time variable thus gives

and an additional Fourier series expansion of yields

[00116] Assume now that all of for each fixed have support in some domain Ω such that the set

can be mapped to some domain such that

Specifically, this mean that there exists a mapping

such that

[00117] In this case it follows that for each it is possible to obtain linear equations for the N _P unknowns by a suitable (as above) remapping of the left hand expression of (28) with and varying over intervals of

length respectively.

Example

[00118] Hexagonal sampling: A simple example of apparition for hexagonal sampling is shown in Fig 4. Here four functions are used, with the regular shift operators and

top panel contains blended data in the Fourier domain; the lower panel the recovered four functions. The previous

formalism can be made easier in this case.

Generalizations

[00119] The setup can be generalized by incorporating

additional coordinate transformations. The simplest such class of transformations concerns individual affine maps for each of the functions. In this case the expression A in (28) is replaced by individual matrices A _p as well as modulation factors corresponding to shifts in the affine maps. The relations (28) now still relate the Fourier transform of the encoded data d on a Bravais grid, with the Fourier transforms of the underlying functions f _p on individual Bravais grids.

[00120] An even more general setup can be obtained by

considering arbitrary coordinate transformations

(diffeomorphisms ) R _p and consider samplings at grids of the type R _p (k) . Such grids could for instance include unintentional curvature caused in the acquisition process or for instance coil sampling. Instead of considering the functions f _p

directly, we can consider new functions of the form

[00121] f _p (R _p (k)) = g _p (Ak).

[00122] The relation (28) can now be derived for the functions g _p . If the maps R _p are mild in some sense (like in the cases mentioned above) , then properties of support hold

approximately (using essential support) for the remapped functions .

[00123] Thus, the term "synchronous shooting" or "synchronous activation" will be used in the context of: a way of acquiring data with multiple sources where the relative locations of the consecutive shots of each source trace out the same or an affinely transformed version of the spatial pattern (i.e., a grid) , including Bravais sampling grids, as a function of time and where the sources are activated near-simultaneously and one or more source activation parameters are varied

(a) periodically across the spatial pattern of shots.

[00124] Further, the term "multi-scale apparition" will be used in the context of: synchronous shooting but where the

(regular) spatial pattern traced out as a function of time by the relative locations of consecutive shots of each of the sources is a scaled version of the (rotated or skewed) spatial pattern of the other sources.

[00125] Note that in synchronous shooting it is the relative locations of consecutive shots (not sources) that form the same or an affinely transformed version of the spatial pattern (or grid) . Hence the shot positions of the individual sources will, in general, be different and the actual shot positions of one source can be related to the actual shot positions of another source by, in the more general case, affine

transformations .

[00126] Similarly, for the case of multi-scale apparition, the actual shot positions of one source can be related to the actual shot positions of another source by affine

transformations. Thus, multi-scale apparition in the following is only a special case of the more general case of synchronous acquisition of affinely related grids.

[00127] Some limited non-exhaustive list of examples of

(regular) spatial patterns that can be used for synchronous shooting are: (1) Equidistant sampling with the sampling interval in y equal to the sampling interval in x, (2)

Equidistant sampling with the sampling interval in y unequal to the sampling interval in x, (3) (intentionally) Irregular sampling, and (4) Hexagonal sampling.

[00128] For synchronous shooting, the offset vectors relating the synchronously acquired positions, or more generally the affine transformations relating the shot locations, determine the amount of overlap in the subsurface coverage ranging from:

(1) no offset, i.e. the simultaneous shots are on top of each other and the amount of overlap is 100%. With the exception of acquiring the grids using two or more different source depths, there is nothing to encode with the corresponding exception of the different source ghost responses as all source-receiver positions define the same 3D Green' s function, to (2) small offsets on the order of where Δχ and Δγ

denote the inline shot spacing and the crossline line spacing, respectively. The simultaneous shots are on grids/patterns that are staggered with respect to each other by a small amount such as, e.g., half the inline and crossline shot spacing. In this case, clearly the shifts of the simultaneous sources are small compared to the extent (e.g., line length) of the shot pattern and all sources essentially still see the same subsurface. Reconstruction methods can be applied that take advantage of the fact that the sources see the same subsurface, to

(3) large offsets on the order of Lx~Ly where Lx and Ly denote, e.g., the shot line length and the crossline extent of all shotlines, respectively. In this case, the simultaneous shots are completely separated in terms of their midpoints and all sources see a different part of the full 3D Green's function and illuminate a different part of the subsurface. No overlap can be exploited and the number of unknowns not reduced. On the other hand, this type of acquisition is the most efficient way to acquire, e.g., a 3D common receiver gather in part or in full.

Kx-Ky space-filling multi-scale apparition

[00129] One advantageous aspect of the generalised apparition methods and equations described above for the cases where the sampling grids are related by affine transformations and/or are Bravais grids, is now illustrated in more detail, namely, the ability better fill the available (empty) Kx-Ky space by choosing suitably scaled (and offset/rotated/skewed) versions of the related source sampling grids. This is illustrated in Figures 5-13, where synthetic simultaneous source data are generated for a model with 64 point diffractors embedded in a homogeneous 3D background model (propagation velocity c = 1500 m/s) . Four sources are acquired simultaneously over the diffractors with regular equispaced shot grids. [00130] Fig. 5 Left panel shows the shot location grids for four simultaneous sources and the location of 64 point

diffractors in a subsurface model. Right panel shows a zoom-in on the shot location grids around the origin. Two of the sources acquire their shot locations with one grid spacing while the other two sources acquire their sources with a different grid spacing.

[00131] The shot spacing of the shot grids associated with sources two and three is ( 1-sqrt (2 ) /2 ) / (sqrt (2 ) /2 ) ~0.4145 times smaller than the scale and shot spacing of the shot grids associated with sources one and four.

[00132] Fig. 6 shows the encoding time delays for the shot locations in the shot location grids of each of the four simultaneous sources of Fig. 5. The 2D modulation functions (specifically, the time delays applied to each of the sources once they simultaneously reach a particular shot point) are shown in Figure Y.

[00133] Note that without changing the scale (shot spacing) of the shot grids, and assuming that source wavefields all occupy the same bandwidth, then aliasing occurs when the apparated part of the signal cone of sources two and three starts to overlap the unapparated part of these sources and the

unapparated source one. On the other hand, at that frequency, it is clear that the apparated part of the signal cone of source 4, which has been apparated to the corners of the Kx-Ky domain, is not yet overlapping with the unapparated part of source 4, nor with the unapparated source one. This is because the diagonal is a factor sqrt (2) longer and overlapping of sources 4 and 1 therefore occurs at a proportionally higher frequency. Thus, it is advantageous if we can reduce the relative bandwidth occupied by the sources two and three. This is possible using multi-scale apparition and indeed is the reason for the different scales (shot spacings) in Fig. 6. By decreasing the spacing for shots, and hence the scale of the shotgrids, two and three, we can create a situation where the wavefields from these sources occupy exactly the remaining space in the Kx-Ky at the aliasing frequency for apparated source 4 and unapparated source 1. This is illustrated in detail in Fig . 7.

[00134] Fig. 7 Shows the regions of the inline-crossline (Kx- Ky) horizontal wavenumber domain that each of the four sources of Fig. 5 map into after encoding the shots with the the encoding time delays of Fig. 6 and using the different

(smaller) scales (grid spacings) for sources 2 and 3 compared to the scales (grid spacings) used for sources 1 and 4. Source one is not apparated. Source four is apparated to the corners. Sources two and three are apparated to the

respectively .

Kx-Ky space-filling multi-scale apparition - Results

[00135] In the Figures 8,12,13, repeatedly three slices from 3D volumes of t,x,y data are shown. In all cases, the top panel is a time slice at t=0.796s, whereas the middle panel is a t-y slice at the central x position of the corresponding shot grid, and the bottom panel is a t-x slice at the central y position of the corresponding shot grid. Note that for the simultaneous source (remodeled) input data and residual, the shot positions from source one are shown for reference.

[00136] In the Figures 9-11, repeatedly three slices from 3D volumes of F,Kx,Ky data are shown (amplitude spectra) . In all cases, the top panel is a frequency slice at f=35Hz, whereas the middle panel is a F-Ky slice at Kx=0/m, and the bottom panel is a F-Kx slice at Ky=0/m. Note that for the

simultaneous source (remodeled) input data and residual, the wavenumber axes from source one are shown for reference.

[00137] Fig. 8 Shows three slices (top, middle, bottom) through the 3D volume of synthetic simultaneous source data modeled for the four sources and their shot locations in Fig. 5, encoded using the time delays of Fig. 6, simultaneously shooting into a single receiver at the origin above a model with 64 point scatterers (left column) . In the middle column, the same three slices of the separated, remodeled data are shown. In the right column, the residual (i.e., the unmodeled part of the total simultaneous source wavefield) is shown.

Note the low residual, indicating that the separated source model the input data very closely.

[00138] Fig. 9 Shows three slices (top, middle, bottom) through the 3D volume of the frequency-inline-crossline (F-Kx-Ky) wavenumber amplitude spectrum of the synthetic simultaneous source data modeled for the four sources and their shot locations in Fig. 5, encoded using the time delays of Fig. 6, simultaneously shooting into a single receiver at the origin above a model with 64 point scatterers (left column) . In the middle column, the same slices are shown of the separated, remodeled data. In the right column, the residual (i.e., the unmodeled part of the total simultaneous source wavefield) is shown .

[00139] In the frequency slices of Fig. 9, one can see the effect of multi-scale apparition where the wavefields due to sources 2 and 3 are mapped, due to their smaller scale (shot spacing) into the remaining space in the Kx-Ky domain left after apparating source 4 and the unapparated source 1. One can also note the low values (close to zero) in the amplitude spectrum of the residual (right column) .

[00140] Fig. 10 Shows three slices (top, middle, bottom) through the 3D volume of the frequency-inline-crossline (F-Kx- Ky) wavenumber amplitude spectrum of the synthetic

simultaneous source data modeled for the fourth source and its shot locations in Fig. 5, activated independently and without encoding time delays, shooting into a single receiver at the origin above a model with 64 point scatterers (left column) . In the middle column, the same slices are shown for the separated data. In the right column, the difference between the modeled (reference) and separated data for the fourth source is shown. [00141] Fig. 11 Shows three slices (top, middle, bottom) through the 3D volume of the frequency-inline-crossline (F-Kx- Ky) wavenumber amplitude spectrum of the synthetic

simultaneous source data modeled for the third source and its shot locations in Fig. 5, activated independently and without encoding time delays, shooting into a single receiver at the origin above a model with 64 point scatterers (left column) . In the middle column, the same slices are shown for the separated data. In the right column, the difference between the modeled (reference) and separated data for the third source is shown.

[00142] Comparing Fig. 10 and 11, the relative bandwidth occupied by the different sources four and three is very clearly visible. The relative bandwidth occupied by source three is significantly lower due to the reduced scale (i.e. shot spacing) of the grid.

[00143] Fig. 12 Shows three slices (top, middle, bottom) through the 3D volume of synthetic simultaneous source data modeled for the fourth source and its shot locations in Fig. 5, activated independently and without encoding time delays, shooting into a single receiver at the origin above a model with 64 point scatterers (left column) . In the middle column, the same slices are shown for the separated data. In the right column, the difference between the modeled (reference) and separated data for the fourth source is shown.

[00144] Fig. 13 Shows three slices (top, middle, bottom) through the 3D volume of synthetic simultaneous source data modeled for the third source and its shot locations in Fig. 5, activated independently and without encoding time delays, shooting into a single receiver at the origin above a model with 64 point scatterers (left column) . In the middle column, the same slices are shown for the separated data. In the right column, the difference between the modeled (reference) and separated data for the third source is shown. [00145] Comparing Fig. 12 and 13, the different scales at which the wavefield has been sampled are clearly visible, especially in the time slices. The separations shown in Fig. 10-13 show low error indicating that the sources relative bandwidth of the sources has been modeled accurately and the separation problem, due to the more parsimonious filling of the Kx-Ky domain, has become well-posed for the entire F-Kx-Ky domain at least until 35 Hz.

Double encoding acquisition and double decoding processing

[00146] Novel methods have been described herein for

simultaneous source separation and for source deghosting. We now describe how both simultaneous source separation and source deghosting can be achieved on single dataset (e.g.

common receiver gather for a single shot line) by acquiring the data in a specific and advantageous manner. Consider two sources SI and S2 at different depths zl and z2 where S2 is encoded according to seismic apparition principles (for example: every first of two consecutive shots, SI and S2 are activated simultaneously whereas every second of two

consecutive shots S2 is activated with a variation as

described herein, e.g., with a time delay of 10ms with respect to the activation time of SI) . Using the methods described herein, the data from SI and S2 can be separated. Traces from the separated data, which still have free-surface reflections or ghosts in them, can then be interleaved in such a manner that the previously described source deghosting method can be applied. More specifically, as an example, one can interleave all the odd traces from the separated SI data and all the even traces from the S2 data, resulting in a dataset with varying source depth, on which the deghosting method can be applied.

[00147] An example is shown in Figure 14 for a simple 3 reflector model. Source SI is activated at 15m depth and source S2 is activated at 25m depth. Every first of two consecutive shots, S2 is activated simultaneously with SI. Every second of two consecutive shots, S2 is activated 16ms following SI. Every fifth trace of the simultaneous source input data is shown in the leftmost panel. In the 2nd and 3rd panels from the left, the separated SI and S2 data are shown, respectively. In the 4th panel from the left, the even traces are taken from the separated S2 data and interleaved with the odd traces from the separated SI data. Source deghosting according to the method described previously is applied to this data and the resulting deghosted SI data are shown in the 5th panel from the left.

[00148] This example demonstrates how seismic apparition principles can be combined in a single acquisition, with source deghosting principles resulting in broadband data. More importantly, it shows how data can be double encoded and decoded with a particular acquisition and processing scheme. Variations of these methods may be derived, as well as further benefits may be obtained by performing the sim-source

separation and the source-side deghosting jointly. Next, we disclose another acquisition and processing method where the data are double encoded and can be double decoded.

[00149] As another embodiment of double (or more) encoding and double (or more) decoding, consider an acquisition with three sources where the first source SI is positioned half the inline shot spacing in front of the second source S2, and the third source S3 has an arbitrary crossline offset from the second source. In this particular embodiment, S2 and S3 are always activated at the same time. On every first of two consecutive shots, S2 and S3 are (jointly) activated

simultaneously with SI and on every second of two consecutive shots, S2 and S3 are (jointly) activated, e.g., 12 ms after SI. As, S2 and S3 are jointly encoded relative to SI using apparition principles, they can be separated from SI using the techniques described earlier. Then, interleaving the traces from the separated SI and S2+S3 datasets, it is clear that SI and S2 form a shotline that is twice as densely sampled compared to the shot spacing with which the simultaneous source data was acquired. On the other hand, every second trace still consists of a superposition of S2 and S3 and it is advantageous to realize that this new dataset corresponds to another case of seismic apparition type simultaneous source acquisition, namely where the S3 source is amplitude encoded relative to the interleaved separated SI and S2 sources with encoding function A=0 on every first of two consecutive shots. In such a case, we will use the terminology that S3 is encoded implicitly. Thus, S3 can be separated from the interleaved separated SI and S2 sources, yielding another dataset that can be used for processing. Once again, further benefits may be derived from performing the two decoding steps jointly.

[00150] The above methods show how different levels of encoding and decoding using apparition principles may be applied in a layered or recursive fashion, where the result of one decoding is possible following another decoding.

Improving crossline sampling by inline apparition

[00151] In many cases, in marine seismic acquisition, the crossline source sampling is poor compared to the inline source sampling. In the following, we disclose two methods that take advantage of a denser inline shot spacing, to improve the crossline source sampling.

[00152] It is advantageous to select the inline shot spacing such that it at least more than two times denser than the crossline shot line spacing. In the following, we describe, without loss of generality the case when the inline shot spacing is around four times denser than the crossline line spacing (e.g., inline shot spacing of 12.5m, crossline shot line spacing of 50m) . Inline apparition can then be used to simultaneously acquire two parallel shot lines and perfectly separate these in processing up to a temporal frequency corresponding to half the aliasing frequency of the unencoded data. Thus, if the 12.5m spacing aliases a horizontally propagating wave at 60Hz, inline apparition decodes the data perfectly up to at least 30Hz. On the other hand, since the shot line spacing in the crossline direction was 50m, a horizontally propagating wave propagating purely in the crossline direction will alias already at 15Hz. An important realization is that the decoding of inline apparated data also decodes any energy propagating in the crossline direction perfectly between 15 and 30Hz. So although the crossline line spacing of 50m still aliases such crossline propagating energy between 15 and 30Hz, we can obtain using the inline apparition technique, at twice as many shot lines in the crossline direction and completely overcome the crossline aliasing between 15 and 30Hz.

Reconstructing a single underlying function

[00153] It is clear that the complete crossline dealiasing between 15 and 30Hz is achieved when using the inline

apparition technique to acquire shot lines exactly in between the shots lines spaced every 50m. On the other hand, this is not the only way to achieve this. The simultaneous shot lines might be offset a small amount in the crossline direction, resulting in effective crossline gradient measurements which could be used in a crossline multi-channel interpolation.

Alternatively, the simultaneous shot lines can, when all taken together, be understood as a filtered (shifted) and decimated version of an underlying unaliased signal. Multiple of such filtered (shifted) and decimated versions can be reconstructed using the well-known Generalized sampling expansion by

Papoulis (1977) . Finally, if the simultaneous shot lines are not regularly sampled in crossline, irregularly spaced Fourier transforms or spectral reconstruction techniques may be used to de-alias the sim-source separated data in the crossline direction .

[00154] It should be noted that the above description very much is adopting a worst case analysis and that by applying various of the de-aliasing techniques described herein, one can reconstruct signals to much higher frequencies than half the original aliasing frequencies. In such cases, it may also be advantageous to adopt data-driven de-aliasing techniques such as the Anti-Leakage Fourier Transform (Xu et al . , 1995) in the crossline direction to reduce the crossline aliasing by more than a factor of two.

[00155] One embodiment will now be described in more detail, including with a workflow. The idea is to use closely located sources and seismic apparition for estimating source

derivatives and sources at sparse samples. From the (e.g., low-pass filtered) derivatives we can estimate the directional interpolation vector field to interpolate the data. We need to de-alias the measured data as well. Up to the Nyquist sampling rate of the un-apparated data. This can be done in two steps: 1. Estimate directional interpolation vector field from low- frequency data. 2. Obtain local shifts (in both source and receiver coordinates) . 3. Partition the data in smaller domains (using a partition of unity) . On these domains, use constant (local) shift parameters in the source and receiver coordinates. 4. Include the shifts in the apparition setup to only recover one function. The derivative information is now captured in the shifts. 5. Reconstructions can then be made at the original sampling rate in the source coordinate. Following this we can then recover higher frequency content using the directional interpolation vector field in a similar way as is described for the directional de-aliasing. This procedure can also be iterated and/or combined with (Normal) Moveout

corrections for broadening the frequency ranges. These ideas can also be used in combination with conventional sim-source separation methods where incoherency exploited in the

separation instead of apparition. In addition, these methods can be applied both inline and crossline.

Multi-dimensionality and separation in higher dimensional spaces

[00156] The methods described herein have mainly been

illustrated using so-called common receiver gathers, i.e., all seismograms recorded at a single receiver. Note however, that these methods can be applied straightforwardly, after (local) transforms over one or more receiver coordinates, to

individual or multiple receiver-side wavenumbers . Processing in such multi-dimensional or higher-dimensional spaces (i.e., higher than, e.g., the F,Kx,Ky space described previously) is known to reduce aliasing or overlapping issues, as the empty space increases dramatically relative to the bandwidth

occupied by the seismic signals.

[00157] We note that further advantages may derive from applying the current invention to three-dimensional shot grids instead of two-dimensional shot grids, where beyond the x- and y-locations of the simultaneous sources, the shot grids also extend in the vertical (z or depth) direction. Furthermore, the methods described herein could be applied to different two-dimensional shot grids, such as shot grids in the x-z plane or y-z plane. The vertical wavenumber is limited by the dispersion relation and hence the encoding and decoding can be applied similarly to 2D or 3D shotgrids which involve the z

(depth) dimension, including by making typical assumptions in the dispersion relation, e.g., in the 2D case, such as Ky=0/m

(when acquiring, e.g., shot grids in the x-z plane) .

Decoding after/during imaging or migration

[00158] The above discussion on separation after (local) transforms over one more receiver coordinates also makes it clear that seismic apparition principles can be applied in conjunction with and/or during the imaging process: using oneway or two-way wavefield extrapolation methods one can

extrapolate the recorded receiver wavefields back into the subsurface and separation using the apparition principles described herein can be applied after the receiver

extrapolation. Alternatively, one could directly migrate the simultaneous source data (e.g., common receiver gathers) and the apparated part of the simultaneous sources will be

radiated, and subsequently extrapolated, along aliased

directions, which can be exploited for separation (e.g. by recording the wavefield not in a cone beneath the sources, but along the edges of the model) .

Source configurations and types [00159] The methods described herein facilitate more efficient data acquisition. On the source side, this includes

configurations that allow for acquiring data at different source depths, in parallel (simultaneously) or sequentially. Figures 15, 16 and 17 illustrate source acquisition

configurations, for any type of source, that provide further benefits as will be clear from the descriptions.

[00160] In Figure 15, methods for towing of source arrays at different tow depths but offset to each other are disclosed. Multiple source arrays can be at a single depth level,

potentially excited simultaneously. This source configuration or layout is different to existing methods as the simultaneous or near-simultaneous shots of the sources at different tow depths do not interfere with each other, while the encoding methods disclosed herein enable separating the source side ghost. The disclosed configurations constitute an efficient way of acquiring seismic data for some of the methods

described herein.

[00161] In Figure 16, various configurations for towing one or more source arrays at varying or alternating tow depths are disclosed, including: (1) Periodically varying, e.g. on a sine or cosine trajectory, (2) Aperiodically varying, and (3) shooting simultaneously or shooting sequentially. The

configurations with varying tow depths teach how some of the methods described herein could be acquired efficiently.

[00162] In Figure 17, sub-arrays of a seismic source array that are towed at different depth levels, including at a constant depth and at a periodically or aperiodically varying depth, are shown. These subarrays are excited individually or can be excited in concert with subarrays from other source arrays. This method is different to traditional over/under sources as it does not aim to overcome the source side ghost notch by shooting over/under source sub-arrays simultaneously, but offers an acquisition scheme that can be processed and encoded for separation using the methods described herein. Examples are V-shaped or A-shaped configurations of source sub-arrays. Receiver configurations and types

[00163] As mentioned previously, the methods described herein are also applicable to the receiver side. A specific

embodiment which enables deghosting data on the receiver side is now presented in more detail. This embodiment exploits the tides that change the physical distance between receivers on the seafloor, the source position and the sea surface. When acquiring seismic data, for instance as part of a time-lapse survey, such changes can be utilized in a similar fashion to described source side deghosting approach, using at least two vintages of seismic data, to deghost the seismic wavefield on the receiver side. This method can also be combined with source-side enabled methods, jointly or sequentially, to provide source-side and receiver-side deghosted seismic data.

[00164] Furthermore, periodic and aperiodic modulation

functions can directly be imposed on individual receivers through various means to, for instance, have the amplitude strengths varying. Such modulation function, potentially in combination with modulation functions on the source side, can be employed taking advantage of the methods described herein to for instance more efficiently de-alias the data as part of the acquisition and / or seismic processing.

[00165] While various embodiments of the present invention have been described above, it should be understood that they have been presented by way of example only, and not of limitation. For example it should be noted that where filtering steps are carried out in the frequency-wavenumber space, filters with the equivalent effect can also be implemented in many other domains such as tau-p or space-time.

[00166] In addition, note that although the embodiments have focused on periodic and non-periodic variations in at least one source or source acquisition parameter, the methods described apply equally well to cases where such variations are in the receiver acquisition parameters, as can be

understood by applying such principles as source-receiver reciprocity and/or invoking invariance for translation along the horizontal coordinate ( s ) . A particular application in point is that of receiver-side deghosting.

[00167] Also note that while some of the methods and

embodiments have been described and illustrated by way of two- dimensional theory, processing, data, and/or examples, they can be applied/apply equally to processing of three- dimensional data and, as can be appreciated by someone of ordinary skill in the art, can be generalised to operate with three-dimensional processing on three-dimensional data or even four- or five-dimensional processing by jointly considering simultaneous source data recorded at multiple receivers.

[00168] In addition, we note that it can be advantageous to process and separate local subsets of simultaneous source data acquired using the methods and principles described herein. Processing local subsets can further reduce the aliasing ambiguity and improve the separation.

[00169] As should be clear to one possessing ordinary skill in the art, the methods described herein apply to different types of wavefield signals recorded (simultaneously or non- simultaneously) using different types of sensors, including but not limited to; pressure and/or one or more components of the particle motion vector (where the motion can be:

displacement, velocity, or acceleration) associated with compressional waves propagating in acoustic media and/or shear waves in elastic media. When multiple types of wavefield signals are recorded simultaneously and are or can be assumed (or processed) to be substantially co-located, we speak of so- called "multi-component" measurements and we may refer to the measurements corresponding to each of the different types as a "component". Examples of multi-component measurements are the pressure and vertical component of particle velocity recorded by an ocean bottom cable or node-based seabed seismic sensor, the crossline and vertical component of particle acceleration recorded in a multi-sensor towed-marine seismic streamer, or the three component acceleration recorded by a microelectromechanical system (MEMS) sensor deployed e.g. in a land seismic survey.

[00170] The methods described herein can be applied to each of the measured components independently, or to two or more of the measured components jointly. Joint processing may involve processing vectorial or tensorial quantities representing or derived from the multi-component data and may be advantageous as additional features of the signals can be used in the separation. For example, it is well known in the art that particular combinations of types of measurements enable, by exploiting the physics of wave propagation, processing steps whereby e.g. the multi-component signal is separated into contributions propagating in different directions (e.g., wavefield separation) , certain spurious reflected waves are eliminated (e.g., deghosting) , or waves with a particular (non-linear) polarization are suppressed (e.g., polarization filtering) . Thus, the methods described herein may be applied in conjunction with, simultaneously with, or after such processing of two or more of the multiple components.

[00171] Furthermore, in case the obtained wavefield signals consist of / comprise one or more components, then it is possible to derive local directional information from one or more of the components and to use this directional information in the reduction of aliasing effects in the separation as described herein in detail.

[00172] Further, it should be understood that the various features, aspects and functionality described in one or more of the individual embodiments are not limited in their

applicability to the particular embodiment with which they are described, but instead can be applied, alone or in various combinations, to one or more of the other embodiments of the invention .

[00173] For example, it is understood that the techniques, methods and systems that are disclosed herein may be applied to all marine, seabed, borehole, land and transition zone seismic surveys, that includes planning, acquisition and processing. This includes for instance time-lapse seismic, permanent reservoir monitoring, VSP and reverse VSP, and instrumented borehole surveys (e.g. distributed acoustic sensing) . Moreover, the techniques, methods and systems disclosed herein may also apply to non-seismic surveys that are based on wavefield data to obtain an image of the subsurface .

List of cited References

[Halliday et al . , 2014] Halliday and Laws, Seismic acquisition using phase-shifted sweeps: US Patent application

US20140278119A1 (2014) .

[Shannon, 1949] C. E. Shannon, Proc. Inst, of Radio Eng. 37, 10 (1949) .

[Nyquist, 1928] H. Nyquist, Trans. AIEE. 47, 617 (1928) .

[Ikelle, 2010] L. T. Ikelle, Coding and Decoding: Seismic Data: The Concept of Multishooting. (Elsevier, 2010), Vol. 39.

[Abma et al . , 2015] R. Abma, D. Howe, M. Foster, I. Ahmed, M. Tanis, Q. Zhang, A. Arogunmati and G. Alexander, Geophysics. 80, WD37 (2015) .

[Kumar et al . , 2015] R. Kumar, H. Wason and F. J. Herrmann, Geophysics. 80, WD73 (2015) .

[Mueller et al . , 2015] M. B. Mueller, D. F. Halliday, D. J. van Manen and J. 0. A. Robertsson, Geophysics. 80, V133

(2015) .

[Robertsson, 2015] Robertsson, J. 0. A., R. M. Laws, and J. E. Kragh, 2015, Marine seismic methods, in Resources in the near- surface Earth (eds. L. Slater and D. Bercovici), Treatise on Geophysics, 2 ^nd edition (ed. G. Schubert), Elsevier-Pergamon, Oxford .

[Bracewell, 1999] R. Bracewell, The Fourier Transform & Its Applications (McGraw-Hill Science, 1999) .

[Barnes, 1992] A. E. Barnes, GEOPHYSICS, 57(5), 749-751

(1992) . [Stockwell, 1996] R. G. Stockwell, L. Mansinha, and R. P.

Lowe. Signal Processing, IEEE Transactions on 44(4), 998-1001 (1996) .

[Dragoset, 1995] W. H. Dragoset. Geophysical applications of adaptive-noise cancellation. 65 ^th SEG meeting, Houston, USA, Expanded Abstracts, 1389-1392.

[Trickett, 2003] S. R. Trickett, Stretch-free stacking. In 2003 SEG Annual Meeting. Society of Exploration Geophysicists.

[Zhang et al . , 2012] B. Zhang, K. Zhang, S. Guo, and K. J. Marfurt, Nonstretching NMO correction of prestack time- migrated gathers using a matching-pursuit algorithm.

Geophysics, 78(1), U9-U18 (2012) .

[Papoulis, 1977] A. Papoulis, Generalized sampling expansion. Circuits and Systems, IEEE Transactions on 24, no. 11, 652-654 (1977) .

[Xu et al . , 2005] S. Xu, Y. Zhang, D. Pham, and G. Lambare, Antileakage Fourier transform for seismic data regularization . Geophysics 70, no. 4, V87-V95 (2005) .