WO2018087648A2

WO2018087648A2 - Seismic acquisition and processing method

Info

Publication number: WO2018087648A2
Application number: PCT/IB2017/056926
Authority: WO
Inventors: Kurt Eggenberger; Johan Robertsson; Dirk-Jan Van Manen; Fredrik Andersson
Original assignee: Seismic Apparition Gmbh
Priority date: 2016-11-10
Filing date: 2017-11-06
Publication date: 2018-05-17
Also published as: GB2555821A; WO2018087648A3

Abstract

Methods are described for optimally preconditioning and separating the unknown contributions of two or more sources from a commonly acquired set of wavefield signals while varying at least one acquisition parameter between the sources.

Description

Seismic Acquisition and Processing Method

Field of the invention

[0001] The present invention relates to methods for separating contributions from two or more different sources in a common set of measured signals, representing a wavefield,

particularly of seismic sources and of 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 vibrator 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 are robust with respect to perturbations of any kind and/or aliasing.

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, 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:

Fig. 1A illustrates how in a conventional marine seismic survey all signal energy of two sources typically sits inside a "signal cone" (horizontally striped) bounded by the propagation velocity of the recording medium.

Fig. IB illustrates how this energy can be split in a transform domain by applying a modulation to the second source. The aliasing/mixing effects are highlighted by speckles/dots in the area of overlap.

Fig. 1C illustrates the impact slow traveling, wavefields of low-frequency content (e.g. surface waves) can have on the aliasing/mixing effects.

Fig. ID illustrates in addition to Fig 1C the impact that perturbations can have on the shape/sharpness of energy-bounded areas .

Fig. 2 illustrates a master workflow that includes the methods (also called techniques) being discussed herein and that are attributed to preserving/enhancing signal apparition effects and/or to address data aliasing/data ambiguity. The displayed arrangement represents one of many possible alternative embodiments. Individual methods can be applied in arbitrary order and can be included or omitted from the master workflow. Furthermore, any method can inform another method.

Fig. 3 illustrates a method (technique 1) of optimal preservation and use of periodicity in presence of noise .

Fig. 4 illustrates a method (technique 2) of

determining and compensating for data perturbations attributed to the source side.

Fig. 5 illustrates a method (technique 3) of

determining and compensating for data perturbations originating heterogeneities of any kind between source and receiver station.

Fig. 6 illustrates a method (technique 4) of

determining region of effective numerical support for parameters and/or workflow selection.

Fig. 7 illustrates a method (technique 5) of wavefield separation to narrow the region of effective numerical support .

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. [0012] 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) . However, for land or seabed applications and with the receivers being coupled to a non-acoustic medium, the slowest observable velocity can be significantly less. For certain subsets of marine seismic data, namely so-called common offset gathers, and depending on the particular

subsurface geology of the area being surveyed, the slowest observable velocity can also be significantly less than

1500m/s .

[0013] Wavefield energy can be split into body waves and surface waves. In general surface waves typically represent the slowest observable velocities.

[0014] For exploration-seismic applications, the body waves are of primary interest as they illuminate the target interval in the subsurface. Therefore, the common practice is to remove surface waves prior to imaging. However, it is reported that surface waves can also be used to precondition and/or to support wavefield imaging (e.g. see Strobbia et al . , 2011)

[0015] Surface waves generated by a particular excitement of a source and for a particular point in time, tend to have a much lower frequency content than corresponding body waves, as a result of stronger attenuation in the near surface layer.

[0016] Body waves can undergo mode conversions where a

pressure wave is converted into a shear wave or vice-versa. By far the majority of seismic sources are intended to generate compressional waves. Shear wave modes travel slower in the subsurface and generally have less high frequency content than compressional waves for a given elastic medium. Compressional- to shear-wave velocity (Vp/Vs) ratios around 2 are common for deeper strata/sediments. However, in the near surface this ratio can be much higher, possibly exceeding values of 10. In case of receivers being placed in an acoustic medium, mode- converted waves shear waves or shear waves more in general are not recorded, whereas compressional waves are. This is the case for towed marine seismic acquisition. With the receivers coupled to an elastic or inelastic medium, shear waves and pressure waves are recorded. This for example corresponds to seabed acquisition and land seismic acquisition.

[0017] Furthermore, the effects of subsurface heterogeneities will be different on surface waves than on body waves. This is linked to the different travel paths and hence, to the

different illumination and sampling of the subsurface by the different wave types.

[0018] 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.

[0019] 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.

[0020] 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.

[0021] 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.

[0022] 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 .

[0023] 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, Robertsson et al (2016) 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 (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: m(n)=±[l + (-l)»]+±i4[l-(-l)»],

which can also be written as m(n) =-[l + eⁱⁿⁿ] +-A[l-eⁱⁿⁿ]. (0.1) [0024] 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 n) , (/c) = ^F( (n)), the following result can be obtained:

T(fn)m(ri)) = ^ (fe) + ^F(/c - k_N) , (0.2) which follows from a standard Fourier transform result

(wavenumber shift) (Bracewell, 1999) .

[0025] 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 k = 0 (H+) and cone centered at k = k_N (H_) for different choices of Α(ω) for signal separation or signal apparition in Eq. (0.2) .

[0026] 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. [0027] This process may be referred to as "wavefield apparition" or "signal apparition" in the meaning of "the act of becoming visible" and is described in Robertsson et al .

(2016) . In the spectral domain, the wavefield caused by the periodic source sequence is nearly "ghostly apparent" and isolated .

[0028] 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 deterministic modulating sequence m(n) , generating the wavefield h. The summed, interfering data f = g + h are

recorded at a receiver location.

[0029] In the frequency-wavenumber domain, where the recorded data are denoted by F = G + H, the H-part 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 k = 0. 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.

[0030] Although the above description has focused on

acquisition along essentially straight lines, the methodology applies equally well to curved trajectories such as coil- shaped tra ectories, circles, or other smoothly varying trajectories or sequences of source activations.

[0031] The concept may be extended to the simultaneous

acquisition of more than two source lines by choosing different modulation functions for each source and it can be applied to higher dimensional source sampling in space (van Manen et al . , 2016a)

[0032] 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 +k_N/2.

[0033] Fig. 1(B) also illustrates a possible limitation of signal apparition. The H₊ and H_ 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 and van Manen et al . (2016b) describe methods how to address this. 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.

[0034] Fig 1 (C) illustrates an example where the recorded wavefield is composed of faster traveling body waves

comprising a broad-band frequency content including relatively higher frequencies and slower traveling surface waves (and/or mode converted waves) that mainly comprise lower frequencies. It can be noted that the maximum non-aliased frequency for a certain spatial sampling is further reduced, compared to Fig 1 (B) .

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

traveltime such that they can be summed resulting in an initial "stack image" of the subsurface with increased signal- to-noise ratio. However, the mix of surface waves and body waves (possibly together with mode conversions) can render this procedure ineffective.

[0036] Seismic apparition, as described by Robertsson et al .

(2016) relies on a periodic modulating sequence m(n) to

generate the wavefield h . In practice it is difficult though 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.

[0037] 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 .

[0038] To some extent such undesired perturbations can be compensated for by using the method by van Manen et al .

(2016c) . However, in this invention we disclose other

techniques to complement seismic apparition simultaneous source separation and to narrow the effective numerical support .

[0039] In this invention we refer to the process of narrowing the effective numerical support of the data in the frequency- wavenumber domain as a process that, for example, removes noise such as ground-roll prior to simultaneous source

separation. This can for instance be done in a common shot gather where the ground roll from multiple simultaneously emitting sources appears coherent. The process of narrowing the effective numerical support of the data in the frequency- wavenumber domain can also refer to removing certain

perturbations such as for instance source and/or receiver static variations for instance induced by near-surface

variations. If simultaneous source separation is carried out in for instance the common offset domain, the latter (i.e., receiver static) perturbations will also lead to a broadening of the effective numerical support of the data in the

frequency-wavenumber domain. It is sometimes the case during wavefield acquisition, that local, geology-driven effects such as anisotropy, topography, subsurface heterogeneities, etc., introduce additional wavefield perturbations that can be much larger than the wavefield "perturbations" introduced in the signal apparition acquisition method (i.e., the perfectly periodic variations in the source modulation functions) .

[0040] Fig 1(D) shows the impact of such perturbations (e.g., source-side statics) on the wavefield acquisition example shown in Fig 1 (C) , where the individual signal cones start to blur, resulting in a non-controlled replication along the wavenumber axis/axes, and smearing of the energy across the frequency-wavenumber space. This effectively constitutes a broadening of the numerical support of the wavefields that are encoded. Similarly, such broadening effects can also be caused by acquisition driven noise (e.g. harmonics) . To summarize, the lower the apparent velocity of an arrival, the broader the effective numerical support of the data (e.g., Fig 1C) . This in turns causes aliasing to occur at lower frequencies during apparition. In land seismic data acquisition, ground-roll is an example of a noise type that arrives at the recording station (and/or propagates away from the source) with a particularly low apparent velocity. Since it is often

considered to be noise, it is therefore of interest to remove the ground-roll in order to narrow the effective numerical support of the data that are being encoded. Another effect that also leads to a broadening of the effective numerical support of the data are perturbations in the data such as receiver statics for instance (e.g., Fig ID) . Such

perturbations introduce irregularities in the data that may mask the periodicity of the modulation function. As a result, the data will appear increasingly incoherent in the frequency- wavenumber domain introducing ambiguity as the broadened numerical support of the data being encoded mapped to multiple different regions of the frequency-wavenumber domain

increasingly overlap. Herein we refer to both

such mechanisms or effect as examples of aliasing. It is therefore a primary objective of this invention to describe methods that reduce aliasing such that the effective numerical support of the data being encoded is narrowed in one or more regions of the frequency-wavenumber domain thereby enabling more effective simultaneous source separation by means of signal apparition.

[0041] The effects sketched in Fig 1 (C) and Fig 1 (D) can be dealt with either in acquisition and/or seismic processing, separately or jointly, and by various methods and their combinations which are described in more detail subsequently.

Exemplary master workflow

[0042] The methods discussed herein pre-condition wavefield data in a manner that allows exploiting the signal apparition effect introduced by modulation functions that periodically change with shot position. Hence, any pre-processing flow has to honour and/or restore the periodic nature of the modulation functions. This in contrast to any conventional processing flow where such elements do not usually require extensive consideration. In addition, the use of bespoke modulation functions on the source side also can offer opportunities to address long-standing data processing challenges. This may both involve careful design of the modulation functions themselves as well as the pattern of how they are applied across the shots. Fig 2. discloses a workflow that comprises various key elements to be considered when

recovering/preserving/optimizing the effects of signal

apparition, which will be covered in detail subsequently. The complexity of such a workflow depends on the type of data (e.g., marine, land, seabed, borehole, etc.) . Hence, the arrangement and the number of workflow elements can vary.

Furthermore, it is understood that the techniques discussed herein apply in various data processing domains.

Technique 1 - optimal preservation and use of periodicity in the presence of noise

[0043] The periodicity introduced in signal-apparition based methods by a systematic variation of the modulation functions from shot to shot has to be preserved/enhanced in the pre^¬ processing up to the processing step where the signal

apparition effect is exploited.

[0044] For instance, it is understood that a certain, data- dependent level of pre-processing is required to condition the data for simultaneous source separation.

[0045] Typically, such pre-processing is geared towards enhancement of the signal-to-noise ratio. This can be achieved by removing noise, by enhancing signal, or by a combination of both. Noise can have an impact on the effective numerical support as shown in Figures 1C and ID. Considering the land seismic case, group forming is an example of such a pre^¬ processing technique commonly employed in seismic processing to attenuate noise and to enhance the signal.

[0046] Whereas some of these techniques are applied in the common shot domain, more recent techniques also explore the use of super groups where neighboring common shot gathers are summed or stacked, including in a frequency-dependent weighted manner (Neklyudov et al . , 2015) . Without proper consideration both in the survey design as well as in the actual processing stage such and other non-linear processing steps can

compromise or destroy the coding relations in the data, making apparition decoding/separation difficult or even impossible.

[0047] A second example of a process that can distort the coding relations in the data, that is the periodic shot-to- shot variation of the modulation function for one or more simultaneous shots, and which can potentially create

decoding/separation ambiguity is the mitigation and

cancellation of harmonic energy introduced by the source. So- called harmonics represent another long-standing engineering and processing challenge (Seriff and Kim, 1970) .

[0048] Method 1 or technique 1, as illustrated in Fig. 3, comprises a new workflow with the prime objective to ensure optimal preservation/enhancement of the encoding periodicity in the presence of noise and noise attenuation processes. In addition, and as a secondary objective only, any remaining degrees of freedom in the choice of the periodicity of the modulation functions can be exploited to attenuate certain noise types. Since the concept of signal apparition allows utilizing a wide range of modulation functions that can be superimposed and combined arbitrarily, careful selection of these modulation functions can allow mitigation and/or

cancelling acquisition related noises while at the same time preserving/enhancing the characteristics of the encoded data that are to be exploited for signal apparition

decoding/separation. One embodiment of method 1 that addresses group forming is shown in Fig.3. Another embodiment addresses the attenuation of harmonics.

[0049] Referring to the example of group forming (Fig. 3) an optimal acquisition parameter set is numerically or

analytically established in a first step by using, e.g., optimization routines explore the multi-dimensional space of possible periodic arrangements and choises of modulation functions that fit the survey requirements and preserve the periodicity when propagated through the group forming process.

[0050] In case of forming super groups (e.g. summing up neighboring shotpoints and/or receiver positions in a common offset) and when acquiring data with least two simultaneously shooting sources, one way of preserving periodicity would be to stack an uneven number of shot points as summing data from an even number of sources with different consecutive amplitudes averages out the encoding.

[0051] An optimal modulation function would also depend on the acquisition geometry envisaged. However, the problem can also be considered the other way around, where the acquisition depends on the modulation function envisaged that is

potentially in interplay with a multitude of other sources. A first set of acquisition parameters then can be tested in the field while the level of apparition and/or the quality of signal separation is examined. Examination can come in various guises and consist of qualitative and quantitative measures. Quantitative examples would determine how strongly and how well the apparated energy appears. This can be achieved by minimization/maximization of energy within a certain region in the transform domain using a cost function, and/or by the use of coherency attributes.

[0052] Similarly, it can also be investigated how well certain noise elements are attenuated. With group-forming mainly targeting surface waves (i.e. ground roll), this could imply choosing a low-frequency, high-wavenumber window for analysis. However, different analysis windows and in other domains are possible too. It is important though to choose an analysis window which reflects the characteristics of a particular noise manifestation and in which it is prominently represented (e.g. strong amplitudes compared to other possible events in the window) .

[0053] Within a bandlimited analysis window containing a particular noise embodiment, the level of noise attenuation can be assessed qualitatively and quantitatively. The measure of root mean square (RMS) is a common criterion to determine the quality of noise attenuation quantitatively. Other

examples of measures used are power spectral densities.

[0054] During an acquisition, a feedback loop enables

deterministically changing the acquisition parameters as a function of shot point position. [0055] Having considered sequential approaches until now, a joint process can also be envisaged. For our example this would imply simultaneous group-forming and wavefield decoding using signal apparition.

[0056] Should one decide to bypass group forming on the receiver side and or source side, one might trade in the challenge of maintaining encoding periodicity in the data through the group forming process with the challenge or need of dealing with ground roll. Ground roll widens the region of effective numerical support and therefore requires another technique to narrow it deal with it. This is discussed more in detail later in Method 5 (Fig. 8)

[0057] Alternatively, a hybrid approach can be envisaged where only a subset of the data is going through noise attenuation (e.g. group forming) or where part of the recorded wavefield is decoded prior to a noise attenuation process, whereas the remaining part will go through noise attenuation first and then be decoded. The periodicity and nature of the modulation functions can also provide opportunities to more efficiently tackle noise (e.g. harmonics), as a secondary objective on the back of preserving periodicity. By carefully choosing source sweep functions, activation times and locations of different vibroseis sources, it is to a large degree possible to

mitigate interference between sources in a slip-sweep

operation. With seismic apparition it is thus possible to excite multiple sources (or points sources, e.g a single vibrator truck) at the same time with different modulation functions that allow to separate them again subsequently. This for instance gives rise to a "simultaneous sweep" where distance separation is not requirement anymore. Furthermore, such a simultaneous sweep can be combined with a slip sweep. Instead of having a single source (or source array, a group of sources that defines a single, averaged shotpoint) operating in slip-sweep mode, multiple, simultaneously shooting sources, potentially next to each other, can operate in a slip seep mode . [0058] Hence, when acquiring data in such a simultaneous sweep mode, this can also be used for group forming on the source side by careful selection of modulation function to attenuate noise (e.g. acquisition noise) .

[0059] It is, however, also possible through careful design of modulation functions, to maintain the point-source nature and simply use the positive and negative interference of point- source data to tackle specific noise. Assume a simultaneous sweep acquisition with 3 vibrator trucks arranged in a

triangle equidistant from each other. This corresponds to a triple source simultaneous shooting configuration, acquiring point sources. In can be envisaged that the sources are simultaneous emitting with amplitudes of A and B on

consecutive shot points: Source 1 = A, B, B, A, B, B, A, B, B, ... and Source 2 = B, A, B, B, A, B, B, A, B, ... and Source 3 = B, B, A, B, B, A, B, B, A, .... This produces a composite wavefield that allows to attenuate certain wavefield

components (e.g. noise) and to reinforce others (e.g. signal) . These modulation function can be combined with small time shifts to further address specific noise types like harmonic noise .

Technique 2&3 - Determining and compensating for data

perturbations

[0060] Perturbations in the recorded seismic data (e.g. source statics) will spread energy mostly along the wavenumber axis/axes in the frequency-wavenumber domain as sketched in Fig. ID. The mitigation of such perturbations can be addressed at various stages within the data processing sequence.

[0061] When employing simultaneous source shooting, source- side perturbations cannot be easily compensated for on a composite shot record. Therefore, it would be advantageous if source-side perturbations are predicted prior the excitement of individual shots, e.g., by the use of any kind of legacy data available. [0062] A multitude of processes and a combination of them can be used to predict these shifts and can incorporate, for example, ray tracing using a near surface model, remote sensing information, LiDAR information, or any topographic, geologic and geophysical information.

[0063] An example of such source-side perturbations are source-field statics that will have to be anticipated in the modulation function prior to firing a particular pair of two or more simultaneous shot to maintain the apparition

periodicity recorded on the receiver side of the simultaneous shot records. Since such source-side perturbations in the most general case affect all simultaneous sources, the modulation function for a particular simultaneous source needs to both anticipate the static time delay/advance associated with that source, as well as the static time delay/advance associated with a designated master source (i.e., typically the source which is unencoded) . Thus, in the case of, e.g., two sources where the second source is apparition-encoded using a simple time delay on every second shot, the timing proposed of the second source with respect to the first, which anticipates all the source statics, consists of the nominal encoding time- delay (i.e., zero on odd shots and the encoding timeshift on the even shots) less the known or estimated static under the second source plus the known or estimated static under the first or reference source. This description refers to method 2, further aspects of which are illustrated in Fig. 4 and with an example given in Fig. 5 A-D.

[0064] Such compensation can also be determined in an

adaptive, data-driven fashion where data is analyzed in real time, during the ongoing acquisition, potentially updating a model and to directly feed back into acquisition by changing one or a multitude of acquisition parameters on the fly.

[0065] When referring to static correction the terminology used herein also includes dynamic corrections and can also encompass lateral compensations. [0066] In addition to perturbations on the source side, perturbations in travel times can also arise from

heterogeneities in the subsurface, in between source and receiver positions. The complexity of these perturbations, which require an exact understanding of the subsurface (often the reason why a seismic survey is being conducted) , render it very difficult to include a compensation of such a

perturbation to the acquisition parameters (e.g., firing interval) as in method 2.

[0067] Hence, asks for a data-driven methodology. From the encoding of the various sources when capitalizing on signal apparition, the nominal positions of the apparated energy is known. Traces are then shifted, scaled, stretched, squeezed, in a time variant of frequency variant manner to enhance the coherence of the energy sitting in the predicted areas or to maximize the energy sitting in individual cones. Cross- correlation between neighboring traces may provide a starting model. The obtained perturbations corrections can be inverted for near-surface models. One embodiment of the method is shown in Fig 6.

[0068] In many cases, in particular due to near-surface variations, static shifts and amplitude variations will be introduced from shot location to shot location. Such a

variation will have a detrimental effect on the intended modulation function to be emitted. In this invention, we refer to the combined effect of the source static and amplitude variations from shot-point to shot-point and the intended (desired) modulation function as the effective emitted

modulation function.

Technique 4 - Determining region of effective numerical support for selection of parameters and/or workflows

[0069] Velocities of surfaces waves can vary significantly within a survey area. These variations will impact the level of aliasing. Furthermore, heterogeneities in the subsurface, often more prominent in the near surface, will introduce perturbations and therefore contribute to data ambiguity.

[0070] Furthermore, the level of aliasing and/or data

ambiguity can have an impact on the workflow used to decode the simultaneous sources data set. The workflow comprises elements discussed herein, but can also comprise a particular decoding algorithm that is dependent on the data.

[0071] A high level of aliasing and/or data ambiguity can require the use of highly sophisticated decoding algorithms which are computationally expensive. Hence, being able to determine the level of efficient numerical support can result in considerable computational savings by streamlining the decoding solution to the data at hand.

[0072] Furthermore, knowledge of the level of numerical support can also loop back to an ongoing acquisition so that acquisition parameters can be modified so that other workflows discussed herein can take benefit.

[0073] For example, the group-forming workflow can be adapted, as part of an ongoing acquisition, to ensure an optimal signal-to-noise ratio output by preserving the periodicity inherent to the signal apparated data.

[0074] The method or technique for determining the level of numerical support can comprise of a data driven and/or a model based approach.

[0075] The model based approach can include real-time model updates that refine the model based on which the level of numerical support is determined.

[0076] The level of numerical support can be determined on a subset of data, for example only using the lower frequencies where the aliasing is minimized or even entirely absent.

However, other band limited data subsets are possible too for instance defining windows with reduced data ambiguity. [0077] Therefore, a data driven adaptive scheme is proposed that determines the data complexity and selects workflow elements accordingly. One embodiment of this method is shown in Fig . 7.

Technique 5 - Wavefield separation to narrow the effective numerical support in one or more regions

[0078] The variability in wave propagation velocities and, consequently, of the corresponding minimum apparent

velocities, can be a challenge for decoding of sources as the level of aliasing/data ambiguity is increased. Therefore, it can be beneficial to separate the wavefield into individual wavefield components, prior to performing certain processing steps like source decoding. Wavefield separation can also encompass separating wavefield components into up- and down- going wavefields.

[0079] The different wave types also sample different levels of subsurface heterogeneity and/or anisotropy and therefore are affected by different levels of perturbations. Separating wavefields allows better and/or more adequately compensating for such perturbations, allowing to preserve and/or enhance the periodicity of the modulation functions to better decode the simultaneous source data.

[0080] After separation and subsequent decoding, individual data sets can be re-combined again.

[0081] Alternatively, and in some embodiments, such a wavefield separation can also be seen as noise attenuation where one of the separated wavefield parts is discarded after potentially having extracted a set of attributes first to inform the processing of the main data set.

[0082] In another embodiment, wavefield separation can be jointly performed with simultaneous source decoding. [0083] However, source decoding can be performed on separated wavefields individually.

[0084] Information on related to the wavefield propagation, for instance polarization information derived from the

wavefield records and can be calculated in 2D, 3D, etc., and can be used to further enhance the ability to decode signal apparated data. The wave propagation attributes provide an additional dimension in which to separate wavefields and hence, also to enhance source decoding.

[0085] Separated wave fields and/or corresponding source decoded data can also be used to serve as an input data set to the other methods discussed h> :rein .

[0086] In summary, wavefield separation is a relevant tool for signal apparition as it allows to narrow the area of effective numerical support to facilitate data encoding. One embodiment of the method is shown in Fig. 8.

Quality control

[0087] Parameter selection can be challenging with a lack of reference data available. To account for this a workflow is introduced that uses a test line or test shots for quality control, for example of the decoding.

[0088] The test line or test shots can be acquired as part of an ongoing survey. However, they can also stem from a legacy seismic survey.

[0089] The size of the subset to be compared with can vary and it can even be a subset of a data gather (e.g. truncated in offset and/or time) .

[0090] The quality control is performed on decoded data against an actual measured data set. Multi-component applications and extensions

[0091] 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 component of the particle motion vector (where the motion can be:

displacement, velocity, or acceleration) associated with compressional waves propagating in acoustic 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.

[0092] 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 signal 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 multi-components.

[0093] 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 and in detail elsewhere.

General

[0094] All methods (or techniques) discussed herein apply to all embodiments of signal apparition. Whereas the discussion takes the example of simultaneous source separation, it is understood that the methodologies discussed equally applies to source and/or receiver deghosting, residual shot noise

attenuation, seismic interference removal etc. and

combinations of them.

[0095] All methods discussed herein can also be executed in a data driven, adaptive manner whilst acquiring the survey.

[0096] All the methods discussed represent techniques that can be combined in any possible manner and they can inform each other. The combination of all of these / some of these methods into a workflow as shown in Fig. 2 further adds to the

uniqueness described herein

[0097] It is further 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] J. 0. A. Robertsson, 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 .

[Strobbia et al . , 2011] C. Strobbia, A. Laake, P. Vermeer, A. Glushchenko, Surface waves: use them then lose them: Surface- wave analysis, inversion and attenuation in land reflection seismic surveying, Near Surface Geophysics, 9(6), pp: 503-514.

[Robertsson et al . , 2016] J. 0. A. Robertsson, L. Amundsen, A. S. Pedersen, 2016, Express Letter: Signal apparition for simultaneous source wavefield separation: Geophys . J. Int., 206(2), 1301-1305.

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

[van Manen et al . , 2016a] van Manen, D. J., Andersson, F., Robertsson, J. 0. A., and Eggenberger, K., 2016, Simultaneous source acquisition and separation on general related sampling grids: GB Patent application No. 1608297.6 filed on 12 May 2016.

[van Manen et al . , 2016b] D. J. van Manen, F. Andersson, J. 0. A. Robertsson, K. Eggenberger, 2016, De-aliased source separation method: GB patent application No. 1605161.7 filed on 28 March 2016.

[van Manen et al . , 2016c] van Manen, D. J., Andersson, F., Robertsson, J. 0. A., and Eggenberger, K., 2016, Source separation method: GB Patent application No. 1603742.6 filed on 4 March 2016.

[Seriff and Kim, 1970] A. J. Seriff and W. H. Kim, Effect of harmonic distortion in the use of vibratory surface sources, Geophysics, 35(2), pp:234-246.

[Neklyudov et al . , 2015] D. Neklyudov, A. Bakulin, M.

Dmitriev, P. Golikov, Enhancing 3D broadband land seismic data with smart super groups for processing and FWI, EAGE expanded abstract, Th P7 06.

Claims

A method for simultaneous source separation based on signal apparition wherein at least one process is applied to recorded wavefield data such that the effective numerical support of a representation of the wavefield data is narrowed in at least one region in the frequency- wavenumber domain.

The method of claim 1, wherein the recorded wavefield data are land seismic recordings and a ground roll attenuation process is used to narrow the effective numerical support of the data in at least one region i the frequency-wavenumber domain.

3. The method of claim 2, wherein the ground roll

attenuation process comprises the application of noise attenuation filters after data acquisition.

The method of claim 2, wherein the ground roll

attenuation process comprises group forming.

The method of claim 1, wherein a harmonics attenuation process is used to narrow the effective numerical suppo of the data in least one region of in the frequency- wavenumber domain.

The method of any of the preceding claims, wherein the process of narrowing the effective numerical support o the data in one or more regions in the frequency- wavenumber domain preserves the characteristics of an intended source modulation function.

7. The method of any of the preceding claims, wherein the process of narrowing the effective numerical support of the data in one or more regions in the frequency- wavenumber domain is achieved by modifying the emitted source signatures at different shot-points such that the characteristics of an effectively emitted modulation function are similar to the intended modulation function,

The method of claim 7, wherein the characteristics of the effectively emitted modulation function are derived by modelling perturbations.

The method of claim 8, wherein determining the effective emitted modulation function involves exploring,

minimizing, maximizing or solving for a cost function at planned surveying locations, while varying at least one parameter .

10. The method of claim 9, wherein cross-correlation is used to establish a starting model.

11. The methods of claims 8 to 10, wherein the

characteristics of the effective emitted modulation function are used to update a subsurface model.

12. The method of claim 11, wherein a subsurface model created or updated such that the updated model allows adjusting acquisition parameters of an ongoing survey in that area.

The method of any of the preceding claims, wherein the effective numerical support of the data in the frequency-wavenumber is narrowed in one or more regions by wavefield separation.

14. The method of claim 13, wherein separated wavefield quantities are recombined after simultaneous source separation .

The method of claims 13 and 14, wherein

directionality information of the wavefield i used to separate the wavefield.

16. The method of any of the preceding claims wherein signal apparition includes the steps of

(a) Obtaining wavefield recordings based on the activation of at least two sources along one or more activation lines varying at least one parameter between the sources from one activation to the following selected from one or more of a group consisting of source signal amplitude, source signal spectrum, source activation time, source location at activation time and source depth, such that the varying causes one or more scaled replications of a wavefield with effective numerical support in the frequency-wavenumber domain of at least one of the sources along the wavenumber axis or axes with the scaled replications partially overlapping a wavefield with effective numerical support of one or more of the other sources with the wavefield with effective numerical support the one or more of the other sources being not affected or also replicated and scaled by the varying;

(b) Separating a contribution of at least one of the at least two sources to the obtained wavefield recordings as generated by the at least two sources individually in the absence of the other sources.

17. The method of claim 16, wherein the obtained

wavefield recordings consist of or comprise multiple components .

18. The method of claim 17, wherein one or more of the multiple components have the source contributions

separated independently.

19. The method of claim 17, wherein the combination of the one or more of the multiple components consists of one or more of: a wavefield separation step, a deghosting step, a redatuming step, a polarisation filtering step, and a multi-channel processing step.

20. The methods of any of the preceding claims applied to all applications of signal apparition for acquisition and processing of seismic data, including simultaneous source separation, wavefield deghosting, residual shot noise attenuation and seismic interference removal.

21. The method of any of the preceding claims applied to land seismic data, marine seismic data, seabed seismic data, permanent monitoring seismic data, time-lapse seismic data, transition zone seismic data or borehole seismic data with (near) surface or downhole placed receivers and/or sources such as VSP, 3D VSP, or

distributed acoustic sensing seismic data.

22. The method of claim 17, wherein two or more of the multiple components have the source contributions

separated jointly.

23. The method of claim 18, wherein one or more of the multiple components have the source contributions

separated using information derived from one or more of the other multiple components.

24. The method of claim 17, wherein one or more of the multiple components are combined before one or more products of the combination have the contributions separated .

25. The method of claim 23, wherein the information

derived from one or more of the other multiple components comprises local directionality information.

26. The method of claim 25, wherein local directionality information is determined jointly from two or more of the multiple components.