Field of the invention

[0001] The present invention relates to methods for acquiring and removing residual energy from a previous shot in a seismic survey .

Background

[0002] Seismic data can be acquired in land, marine, seabed, transition zone and boreholes for instance. Depending on in what environment the seismic survey is taken place the survey equipment and acquisition practices will vary.

[0003] In towed marine seismic data acquisition a vessel tows streamers that contain seismic sensors (hydrophones and sometimes particle motion sensors) . A seismic source usually but not necessarily towed by the same vessel excites acoustic energy in the water that reflects from the sub-surface and is recorded by the sensors in the streamers. The seismic source is typically an array of airguns customarily deployed as a set of sub-arrays, each of which includes a set of individual airguns. These are normally programmed to fire at the same instant, providing a close to instantaneous peak of energy followed by a longer, lower energy output as a result of oscillating air bubbles. A marine source can also be a marine vibrator for instance, which may be a single unit or a set of individual units composing an array. In either case, the intent is to provide a seismic source output which contains as far as possible a broad range of frequencies within the usable seismic frequency ranges, typically from 1-2 Hz up to around 500Hz. In modern marine seismic operations many streamers are towed behind the vessel (3D seismic data acquisition) . It is also common that several source and/or receiver vessels are involved in the same seismic survey in order to acquire data that is rich in offsets and azimuths between source and receiver locations.

[0004] In seabed seismic data acquisition, nodes or cables containing sensors (hydrophones and/or particle motion

sensors) are deployed on the seafloor. These sensors can also record the waves on and below the seabottom and in particular shear waves which are not transmitted into the water. Similar sources are used as in towed marine seismic data acquisition. The sources are towed by one or several source vessels.

[0005] In land seismic data acquisition, the sensors on the ground are typically geophones and the sources are commonly vibroseis trucks. Vibroseis trucks are usually operated in arrays with two or more vibroseis trucks emitting energy close to each other roughly corresponding to the same shot location. In this invention we refer to such source configurations as groups of sources.

[0006] Explosive sources may also be used onshore, which may be one large charge or a series of smaller ones.

[0007] Impulsive marine sources are traditionally formed from a combination of individual energy emitting source elements, typically being of the airgun type, by which a volume of compressed air is released into the water column to produce energy in the preferred frequency spectrum. Each airgun element is typically deployed a few metres below the surface, arranged into arrays of similar units.

[0008] There are various brand names and designs of such units, including but not limited to Sleeve Guns, GI Guns and Bolt Airguns and donut guns. All such units work in a similar way and will be referred to herein as "airgun" for the sake of convenience .

[0009] Each individual airgun unit has a specific volume of air, which can be configured by the user. As each unit is initiated, the air volume is ejected almost instantaneously into the water column, and the resulting bubble rises towards the surface, oscillating with a given periodicity with

decaying amplitude. This continues for up to a second or two. The periodicity is a function of the volume and pressure of the air.

[0010] Individual airgun elements are combined into sub-arrays in various configurations, consisting of airguns with a range of volumes such that the bubble periodicity is different for each airgun element. Airgun units are commonly combined together in such sub-arrays such that the overall output consists of a short, aligned initial output (referred to as the "peak") , followed by a period in which the various bubble periodicity times result in largely destructive interference, in order to make the overall radiating pressure wave, referred to as the sub-array signature, as close as possible to the idealized spike. Such a process is referred to as sub-array tuning, and the techniques involved in this are well

established practice and beyond the scope of this description.

[0011] Each airgun subarray is typically linear, though not universally so, and is usually deployed under some floatation device such that the in-line separation as well as the depth of the airgun elements is controlled and remains consistent at each shot point, resulting in as stable a signature as

possible between each shot.

[0012] The output from a single sub-array - which typically consists of a dozen or fewer individual airguns - is generally considered to be insufficient for mainstream seismic

exploration and reservoir management purposes. It is therefore common practice to use two or more sub-arrays, generally deployed laterally and/or in-line separated by a few

(generally twenty or fewer) metres apart. This separation is user-designed and is aimed at controlling the extent to which the sub-array elements interact with each other.

[0013] The overall result is an array, consisting of two or more sub-arrays, each consisting of multiple airguns, usually of varying volumes such that they form a tuned array. The sub- arrays may be at the same or different depths, depending on the geophysical objectives. For example, some recent

configurations may include a set of sub-arrays deployed at different depths, whose firing times may be staggered such that the down-going wavefront is uniform whilst the up-going wavefront exhibits destructive interference in order to reduce the so-called source ghost effect.

[0014] All of the units are generally (but not universally) excited such that the downgoing energy is created

simultaneously, resulting in a far field signature where the peaks are all aligned.

[0015] Deficiency of low frequencies is generally a concern for seismic sources. In addition to the technique just described, sources are sometimes towed at greater depth to attempt enhancing the lower frequency content. However, towing a source at greater depth will introduce ghost notches within the spectrum of interest for higher frequencies. A composite approach to combine a deeper towed source for lower frequency with a sources towed shallower (i.e. a broadband source) is therefore of great interest.

[0016] After a short period of time, since the source vessel is moving continuously, a subsequent shot is fired after a few seconds. This is generally between five and twenty seconds for mainstream seismic acquisition. The objective, quite apart from giving time for the source vessel to move, is also to allow the energy from each shot-point to decay before the next one is initiated. The energy that remains from the previous shot (residual energy) when starting to record the subsequent shot is referred to as residual shot noise. It typically has a low frequency content and therefore limits data quality in the lower part of the seismic data frequency band.

[0017] Some approaches for seismic data acquisition use shorter shot intervals (two or more seconds), often but not universally combined with some element of timing change on sequential shots in order to limit the impact of the insufficient decay time on sequential shot records. These approaches are referred to as "simultaneous source" and are discussed below. These approaches enable more source points per unit area, albeit at some compromise in terms of

interference or fold.

[0018] An alternative approach to conventional simultaneous source separation is referred to as "Signal Apparition"

(Robertsson et al . , 2016) and discussed in more detail below by which shot points include sequences of individual shots, typically very closely separated in time (for example, each shot point is separated within a few tens of milliseconds, rather than a few seconds) . Individual shots are then

separated using the signal apparition approach which in theory is exact at low frequencies (although for certain and the most common choices of so-called modulation sequences discussed below the separation suffers from poor signal-to-noise ratio at low frequencies) . The signal apparition approach is typically achieved with some variation of timing of shot sequences (but can also be achieved by other variations in shot sequences such as amplitude variations or source

signature variations) and also will benefit from the use of some type of reconstruction technique to mitigate or limit aliasing at higher frequencies. There are no theoretical limitations on the number of shots that can be separated in this way.

[0019] In the following we will refer to all methods for simultaneous source acquisition and separation as well as methods for signal apparition-based source acquisition and separation and quasi-simultaneous source acquisition and separation as methods for simultaneous source acquisition and separation. Note that the source elements will not be fired at exactly the same time as some form of encoding (usually in time) is necessary. The descriptor simultaneous refers to sources that are being excited during the record time of another source. [0020] Traditionally seismic data have been acquired sequentially: an impulsive source, typically formed of two or more airgun sub-arrays or vibroseis units is excited and data are recorded until the energy that comes back has diminished to an acceptable level and all reflections of interest have been captured after which a new shot at a different shot location is excited. Being able to acquire data from several sources at the same time or reducing the time interval from shot to shot is therefore clearly highly desirable. Not only would it allow to cut expensive acquisition time drastically but it could also better sample the wavefield on the source side which typically is much sparser sampled than the

distribution of receiver positions. It would also allow for better illumination of the target from a wide range of

azimuths as well as to better sample the wavefield in areas with surface obstructions. In addition, for some applications such as 3D VSP acquisition, or marine seismic surveying in environmentally sensitive areas, reducing the duration of the survey is critical to save costs external to the seismic acquisition itself (e.g., down-time of a producing well) or minimize the impact on marine life (e.g., avoiding mating or spawning seasons of fish species) .

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

[0022] Simultaneous source acquisition has a long history in land seismic acquisition dating back at least to the early 1980' s. Commonly used seismic sources in land acquisition are vibroseis sources which offer the possibility to design source signal sweeps such that it is possible to illuminate the sub ^¬ surface "sharing" the use of certain frequency bands to avoid simultaneous interference at a given time from different sources. 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. Such approaches are often referred to as slip sweep acquisition techniques. In marine seismic data contexts the term overlapping shooting times is often used for related practices. Moreover, it is also possible to design sweeps that are mutually orthogonal to each other (in time) such that the response from different sources can be isolated after acquisition through simple cross-correlation procedures with sweep signals from individual sources. We refer to all of these methods and related methods to as "time encoded simultaneous source acquisition" methods and "time encoded simultaneous source separation" methods.

[0023] The use of simultaneous source acquisition in marine seismic applications is more recent as marine seismic sources (i.e., airgun sources) do not appear to yield the same

benefits of providing orthogonal properties as land seismic vibroseis sources, at least not at a first glance. Western Geophysical was among the early proponents of simultaneous source marine seismic acquisition suggesting to carry out the separation as a pre-processing step by assuming that the reflections caused by the interfering sources have different characteristics. Beasley et al . (1998) exploited the fact that provided that the sub-surface structure is approximately layered, a simple simultaneous source separation scheme can be achieved for instance by having one source vessel behind the spread acquiring data simultaneously with the source towed by the streamer vessel in front of the spread. Simultaneous source data recorded in such a fashion is straightforward to separate after a frequency-wavenumber (ok) transform as the source in front of the spread generates data with positive wavenumbers only whereas the source behind the spread

generates data with negative wavenumbers only.

[0024] Another method for enabling or enhancing separability is to make the delay times between interfering sources incoherent (Lynn et al . , 1987) . Since the shot time is known for each source, they can be lined up coherently for a

specific source in for instance a common receiver gather or a common offset gather. In such a gather all arrivals from all other simultaneously firing sources will appear incoherent. To a first approximation it may be sufficient to just process the data for such a shot gather to final image relying on the processing chain to attenuate the random interference from the simultaneous sources (aka. passive separation) . However, it is of course possible to achieve better results for instance through random noise attenuation or more sophisticated methods to separate the coherent signal from the apparently incoherent signal (Stefani et al . , 2007; Ikelle 2010; Kumar et al . 2015) . In recent years, with elaborate acquisition schemes to for instance acquire wide azimuth data with multiple source and receiver vessels (Moldoveanu et al . , 2008), several methods for simultaneous source separation of such data have been described, for example methods that separate "random dithered sources" through inversion exploiting the sparse nature of seismic data in the time-domain (i.e., seismic traces can be thought of as a subset of discrete reflections with "quiet periods" in between; e.g., Akerberg et al . , 2008; Kumar et al . 2015) . A recent state-of-the-art land example of simultaneous source separation applied to reservoir characterization is presented by Shipilova et al . (2016) . Existing simultaneous source acquisition and separation methods based on similar principles include quasi random shooting times, and pseudo random shooting times. We refer to all of these methods and related methods to as "random dithered source acquisition" methods and "random dithered source separation" methods.

"Random dithered source acquisition" methods and "random dithered source separation" methods are examples of "space encoded simultaneous source acquisition" methods and "space encoded simultaneous source separation" methods.

[0025] A different approach to simultaneous source separation has been to modify the source signature emitted by airgun sources. Airgun sources comprise multiple (typically three) sub-arrays each comprised of several individual airguns or clusters of smaller airguns. Whereas in contrast to land vibroseis sources, it is not possible to design arbitrary source signatures for marine airgun sources, one in principle has the ability to choose firing time (and amplitude i.e., volume) of individual airgun elements within the array. In such a fashion it is possible to choose source signatures that are dispersed as opposed to focused in a single peak. Such approaches have been proposed to reduce the environmental impact in the past (Ziolkowski, 1987) but also for

simultaneous source shooting.

[0026] Abma et al . (2015) suggested to use a library of

"popcorn" source sequences to encode multiple airgun sources such that the responses can be separated after simultaneous source acquisition by correlation with the corresponding source signatures following a practice that is similar to land simultaneous source acquisition. The principle is based on the fact that the cross-correlation between two (infinite) random sequences is zero whereas the autocorrelation is a spike. It is also possible to choose binary encoding

sequences with better or optimal orthogonality properties such as Kasami sequences to encode marine airgun arrays (Robertsson et al . , 2012) . Mueller et al . (2015) propose to use a

combination of random dithers from shot to shot with

deterministically encoded source sequences at each shot point. Similar to the methods described above for land seismic acquisition we refer to all of these methods and related methods to as "time encoded simultaneous source acquisition" methods and "time encoded simultaneous source separation" methods .

[0027] Yet another approach is to fire a sequence of source arrays, one or more of which has a random time dither applied relative to the adjacent source points, but at a shorter time interval, for example, five seconds rather than the

conventional ten. This has the advantage of keeping the shallow part of each shot free from interference, whilst mitigating the drop in fold. For example, conventional

exploration seismic involves two identical source arrays, offset laterally from each other by, for example, 50m (source centre to source centre) . The firing cycle is Port - starboard - port - starboard, such that a source fires every ten seconds, into different sub-surface lines. This results in half-fold data relative to single source. Experiments with triple source using the same approach resulted in 1/3 fold data, considered insufficient. The partially overlapping approach in the above dual source example, would involve firing every 5 seconds, returning to full fold. Employing the same approach with three partially overlapping sources and a five second shot interval would result in limited fold drop and undisturbed shallow data. However, extrapolating this form three to four sources, for example (and temporarily ignoring the issues outlined above about overall sub-array capacity) would require, for example, a 2-3 second shot interval, resulting in limited undisturbed data lengths and loss of fold. Taking into consideration the practicalities, it has also been presented (for example, Hager, 2016), to arrange the firing sequence such that individual airgun sub-arrays may form part of more than one array, as noted above. However, the interference of adjacent shots (even mitigated by dither) and the loss of fold are unavoidable and their effects increase as attempts are made to increase the total number of arrays.

[0028] Recently there has been an interest in industry to explore the feasibility of marine vibrator sources as they would, for instance, appear to provide more degrees of freedom to optimize mutually orthogonal source functions beyond just binary orthogonal sequences that would allow for a step change in simultaneous source separation of marine seismic data.

Halliday et al . (2014) suggest to shift energy in ωλ-space using the well-known Fourier shift theorem in space to

separate the response from multiple marine vibrator sources. Such an approach is not possible with most other seismic source technology (e.g., marine airgun sources) which lack the ability to carefully control the phase of the source signature (e.g., flip polarity) .

[0029] The recent development of "signal apparition" suggests an alternative approach to deterministic simultaneous source acquisition that belongs in the family of "space encoded simultaneous source acquisition" methods and "space encoded simultaneous source separation" methods. Robertsson et al .

(2016) show that by using modulation functions from shot to shot (e.g., a short time delay or an amplitude variation from shot to shot) , the recorded data on a common receiver gather or a common offset gather will be deterministically mapped onto known parts of for instance the ω/c-space outside the conventional "signal cone" where conventional data is strictly located (Figure la) . The signal cone contains all propagating seismic energy with apparent velocities between water velocity

(straight lines with apparent slowness of +-1/1500 s/m in k - space) for the towed marine seismic case and infinite velocity

(i.e., vertically arriving events plotting on a vertical line with wavenumber 0) . The shot modulation generates multiple new signal cones that are offset along the wavenumber axis thereby populating the ωλ-space much better and enabling exact simultaneous source separation below a certain frequency

(Figure lb) . Robertsson et al . (2016) referred to the process as "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. A critical observation and insight in the "signal apparition" approach is that partially shifting energy along the ωλ-axis is sufficient as long as the source variations are known as the shifted energy fully predicts the energy that was left behind in the "conventional" signal cone. Following this methodology simultaneously emitting sources can be exactly separated using a modulation scheme where for instance

amplitudes and or firing times are varied deterministically from shot to shot in a periodic pattern.

[0030] Consider a seismic experiment where 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 (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 m(n) =-[l + e ^i7Tn ] +-A[l-e ^i7Tn ]. (0.1)

[0031] 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) = the following result can be obtained:

T(fn)m(ri)) = ^F(fc) + ^F(k-k _N ), (0.2) which follows from a standard Fourier transform result

(wavenumber shift) (Bracewell, 1999) .

[0032] 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 Α(ω) .

Α(ω) H_ = (l-A)/2 H ₊ = (l+A)/2

1 0 1 -1 1 0

0 1/2 1/2

½ 1/4 3/4

βίωΤ (l-e ^iwT )/2 (ΐ + β ^ίωΤ )/2

1 + β ^ίωΤ -β ^ίωΤ /2 1 + β ^ίωΤ /2

TAB. I. Mapping of signal to cone centered at k = 0 (H+) and cone centered at k = k _N ( _) for different choices of i4(ii>) for signal separation or signal apparition in Eq. (0.2) .

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

[0034] This process may be referred to as, "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.

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

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

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

[0038] The concept may be extended to the simultaneous

acquisition of more than two source lines by choosing

different modulation functions for each source.

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

[0040] Fig. 1(B) also illustrates a possible limitation of signal apparition. The H ₊ and H_ parts are separated within the respective diamond-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. [0041] 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 A = β ^ιωΤ .

[0042] In the prior art it has been suggested to combine different methods for simultaneous source acquisition. Miiller et al . (2015) outline a method based on seismic data

acquisition using airgun sources. By letting individual airguns within a source airgun array be actuated at different time a source signature can be designed that is orthogonal to another source signature generated in a similar fashion. By orthogonal, Miiller et al . (2015) refer to the fact that the source signatures have well-behaved spike-like autocorrelation properties as well as low cross-correlation properties with regard to the other source signatures used. On top of the encoding in time using orthogonal source signatures, Miiller et al . (2015) also employ conventional random dithering (Lynn et al . , 1987) . In this way, two different simultaneous source separation approaches are combined to result in an even better simultaneous source separation result.

[0043] Halliday et al . (2014) describe a method for

simultaneous source separation using marine vibrator sources that relies on excellent phase control in marine vibrator sources to fully shift energy along the wavenumber axis in the frequency-wavenumber plane.

[0044] For time-dithered simultaneous source acquisition Abma et al . (2012) and Jiang and Abma (2010), report that the very low frequencies are compromised when using short time dithers. To overcome this limitation, they suggest to resort to a large dither length of several hundred milliseconds and more

rendering particularly single-vessel multi-source acquisition impractical .

[0045] Signal apparition technology, however, does in

principle not suffer from this apparent drawback. However, there is a separate issue that needs to be considered with respect to the low frequency response for a particular choice of modulation function that is particularly attractive to use for practical reasons.

[0046] For many choices of the factor A governing the

modulation sequence, low frequencies can be separated just as well as higher frequencies (see Table I) . However, a

particularly attractive choice for the factor is A = β ^ιωΤ which amounts to a time shift. Typically a small time shift is chosen so that ωΤ < π for sufficiently high frequencies to avoid nulls (notches) within the frequency band of interest. However, for sufficiently low frequencies we then observe the same problem as well with a notch at DC where ωΤ is small (close to 0) . From inspecting Table I we see that in the limit no energy is shifted to the cone centered at the Nyquist wavenumber. Instead all energy for all sources remain

overlapping in the cone at zero wavenumber.

[0047] As described above, signal generated from previous shots is often still present when exciting a new shot. This remaining energy is referred to as residual shot noise.

[0048] Offshore seismic surveys, whether 2D or 3D in nature, share a requirement that the location of both source and receivers should be known as precisely as possible, in order to ensure that the best possible quality image is produced during data processing. Exact location of receiver positions is a function both of the ability to place them correctly (for example, lateral steering control on towed marine or accurate placement of sensors on the sea floor with OBS), and the ability to know where they are by a combination of acoustic and GPS connected devices. On the source side, the same level of control of geometry (lateral control, acoustic and GPS devices) is also applied. This has the effect of enabling a pre-plotted source firing point location (or sweep initiation for marine vibroseis) to be determined in advance. As the source vessel steams along the line, sources are excited precisely at the point where the pre-determined location is reached. This is referred to as "shooting on distance" or "shooting on position" and refers to a precise distance being travelled between each shot-point. The time taken to travel, however, always varies to some extent, simply because the relevant distance is over the sea floor, and any slight change vessel propulsion effort, waves, currents, through water speed or similar naturally occurring changes will also impact the time taken to move to the next pre-defined shot location. As most modern recording systems actually record continuously, it is therefore critical to know exactly the time stamp when each shot was actually initiated. Whilst the "notional" time between shots may be fairly uniform, the absolute time, which is recorded in the data trace headers, will thus vary. An alternative strategy is to fire at specific time intervals, which will result in more or less distance being covered in the time between this and the preceding shot. This will thus have the effect of regularizing the time sequencing and introducing a variation in the actual distance (related to the sub-surface) from one shot point location to the next. For the purposes of this invention, either approach will introduce some natural (and unavoidable) irregularity to the shot sequence/timing . Furthermore, as explained herein, such shot timings may be deliberately altered (either in some regular pattern, in a non-regular be pre-defined pattern, or via random dithers) from shot point to shot point in order to limit interference between subsequent shots.

Brief summary of the invention

[0049] Methods for acquiring and separating residual shot noise from seismic data exploiting induced or natural

variations of shot times and / or shot positions from shot point to shot point. Induced variations include, but are not limited to, methods based on the signal apparition approach mentioned above or a conventional method for simultaneous source separation such as methods based on random dithering 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.

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

[0051] In particular, the present invention may be regarded as relying on some kind of irregularity from shot point to shot point in terms of time or position. Such variations can be deliberately introduced through methods that often are used for simultaneous source acquisition, for instance encoded, but not limited to, by means of the method of signal apparition or methods using random dithering to encode sources. The

variations from shot point to shot point in time or position can also be natural due to for instance varying vessel speed when shooting on position or variations in position when shooting on time. The invention would apply equally to onshore and offshore seismic surveys, and for implosive, explosive or vibratory type sources. Finally, the invention applies to the removal of seismic interference noise which is a related problem that displays similar characteristics to residual shot noise. [0052] In the current invention we refer to the nominal shot point time interval as the approximate time that the vessel will take to move from one shot point to the next if the seismic data are shot on position.

[0053] In the case that data are shot on position, at each shot point a shot point interval is defined that corresponds to the time that a vessel carrying a source traverses from one shot point to the next where the source is actuated at a shot point at a time belonging to said shot point interval. Note that the shot point interval may not be identical to the nominal shot point time interval as variations in the shot point interval may occur from shot point to shot point either as a result of natural variations in for instance vessel speed or as a result of deliberately introduced time dithers.

[0054] Similarly we refer to the nominal shot point distance interval as the approximate distance that the vessel will traverse from one shot point to the next if the seismic data are shot on time.

[0055] In the case that data are shot on time, at each shot point a shot point interval is defined that in this instance corresponds to the distance that a vessel carrying a source traverses from one shot point to the next where the source is actuated at a shot point at a time belonging to said shot point interval. Note that the shot point interval may not be identical to the nominal shot point distance interval as variations in the shot point interval may occur from shot point to shot point either as a result of natural variations in for instance vessel speed or as a result of deliberately introduced spatial dithers.

[0056] In the current invention we refer to a residual signal as signal present in the signal acquired during a shot point interval and being caused by the source being actuated at a shot point at a time belonging to a previous shot point interval . [0057] In the current invention we refer to the average shot point distance interval as the average distance between shot points in an ensemble of signals recorded on at least one seismic sensor. This average can fro instance be calculated for instance as the arithmetic average of the distances between shot points in the ensemble. The average can also for instance be calculated using the median distance between shot points in the ensemble. These are just two examples of how the average shot point distance interval can be calculated and other suitable means will be known to those skilled in the art .

[0058] In the current invention we refer to the wave

propagation velocity of the recording medium as the velocity that waves would propagate with in an infinite medium with the acoustic or elastic properties of the medium where the seismic sensor is located. For seismic sensors being hydrophones in water, the wave propagation velocity of the recording medium will be that of water or approximately 1500 m/ s .

[0059] In the present description and claims for the sake of clarity and compactness, use is made of compact mathematical expressions, such as terms and equations. However, it is clear to a person skilled in the art that such terms and equations may be represented differently when implemented as numerical computational methods to perform the methods using computing devices. In the current invention we refer to numerical computational methods as computer based implementations of mathematical algorithms in finite-precision arithmetic. An example of a numerical computational method can be a computer- based implementation of so-called LU-decomposition to invert a square matrix. Hence such numerical representations of the terms and equations are meant to fall within the scope of the claims in the same way as the compact mathematical expressions used when they can be shown to represent such compact

expressions .

[0060] It may also be regarded as an advancement provided by the present invention to see the residual shot noise problem as a problem of simultaneous source separation albeit where the sources are not intentionally triggered simultaneously or near-simultaneously .

Brief Description of the Drawings

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

Figs. 1A, B illustrate how in a conventional marine seismic survey all signal energy of sources typically sits inside a "signal cone" 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 a second source;

Fig. 2 shows a synthetic seismic data example generated over a complex synthetic salt model using a finite- difference modeling method without adding residual shot noise;

Fig. 3 shows a synthetic seismic data example generated over a complex synthetic salt model using a finite- difference modeling method after adding residual shot noise with a random time shift between -0.1s and 0.1s applied to the residual shot noise at each trace location;

Fig. 4 shows the ω/t-spectrum of the synthetic data shown in Fig. 2;

Fig. 5 shows the difference in ω/t-spectra of the reference data without residual shot noise shown in Fig. 2 and Fig. 4 and the reconstruction from the data with residual shot noise shown in Fig. 3.

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

[0063] In mathematics, a well-posed problem has a unique solution. A linear system of equations is an example of a problem that may be well-posed or not.

[0064] A linear system of equations can be solved exactly if there are as many equations as there are unknowns and if the equations are linearly independent. In this case only one solution exists, and the problem is said to be well-posed. The coefficients of the linear system are often represented by means of a matrix. A system of equations is considered over- determined if there are more equations than unknowns. An over- determined system is often inconsistent (it has no solution) . However, an over-determined system will have solutions in some cases, for example if some equation occurs several times in the system, or if some equations are linear combinations of the others. Although it is not an exact solution, an over- determined system of equations can be solved ensuring a solution that is as consistent as possible with the available equations (in some mathematical sense, e.g., using a root mean square norm to quantify the mismatch between the solution and the system of equations) . This least squares solution can be obtained by solving the so-called normal equations, in which a new square matrix is constructed by matrix multiplication of the transpose of the original matrix times the original matrix along with a new right hand side that consists of the matrix- vector product of the transpose of the original matrix and the original right hand side. If the columns of the new square matrix are linearly independent then the normal equations are well-posed, and this defines the least squares solution to the over-determined system uniquely.

[0065] In contrast, in mathematics, a system of linear

equations is considered under-determined if there are fewer equations than unknowns. The terminology can be explained using the concept of constraint counting. Each unknown can be seen as an available degree of freedom. Each equation introduced into the system can be viewed as a constraint that restricts one degree of freedom.

[0066] Therefore, the critical case (between over-determined and under-determined) occurs when the number of equations and the number of free variables are equal. For every variable giving a degree of freedom, there exists a corresponding constraint removing a degree of freedom. The under-determined case, by contrast, occurs when the system has been under- constrained - that is, when the unknowns outnumber the

equations .

[0067] Although it is not widely appreciated, the problem of simultaneous source separation is always an over-determined or well-posed problem below a certain frequency. The method of signal apparition maximizes the region where the problem is over-determined. However, also methods for simultaneous source separation based on random dithers can be solved exactly below a certain frequency. The region where the solution can be exactly solved is half the size if a method using random dithers is used compared to if a method based on signal apparition is used (enabling exact separation within diamond shaped regions in ω/c-space as opposed to triangular shaped regions as illustrated in Fig IB; see Andersson et al . , 2017) .

[0068] In this invention we exploit the fact that the presence of the residual shot noise varies from shot to shot. This can be caused by a variation of source signal amplitude, source signal spectrum, source activation time, source location at activation time and source depth for instance.

[0069] Let the recorded data as a function of time t for shot number j be denoted by d(t,j). The data is a superposition of signal from shot number j denoted by /i(t,_/ ^' ) and residual shot noise that is slowly decaying signal from the previous shot number j— 1. The residual shot noise present at shot number j is denoted by ₂ (t— '- · this invention we exploit the fact that the presence of the residual shot noise varies from shot to shot. This can be caused by a variation of shot position or a variation time for instance. In the following

description we limit the description to a variation in time from shot point to shot point for instance caused by natural timing variations due to vessel speed / current variations in a marine seismic setting for instance. The parameter Tj represents this variation in shot time. It is typically known highly accurately as the precise firing time of shots is recorded .

[0070] Residual shot noise can also be present from shots that occurred earlier than just the previous shot. The current invention can also be generalized to remove the interference of such residual shot noise.

[0071] It will be apparent to those skilled in the art that it is also possible to vary the signal firing time from shot point to shot point such that the residual shot noise will occur at the same time from shot to shot. This is achieved using for instance the method of signal apparition. The following description can easily be generalized also to such scenarios as time shifts lining up shots or residual shot noise is a purely relative operation to each other.

[0072] The recorded data can therefore be written as: d(t,)) = A(t,)) + f ₂ (t - Tj,j) . (0.3)

[0073] Taking a Fourier transform in time of equation (0.3) results in:

T _t d(a> ) = Λ(ω,7) + / ₂ (ω,; - ^2πίτω , (0.4) where ω is the angular frequency and and f denote temporal Fourier transforms of f and f ₂ , respectively.

[0074] Next we consider a sequence of N regularly spaced shots and apply a discrete Fourier transform over the shot coordinate : ∑ _j ^N ₌₁ T _t d{t,i e- ² ^ ^k ' ^N = F^,k) + (/ ^? ₂ (6V) * w _w )(/c), (0.5) where k is the discrete wavenumber, F ₁ and F ₂ denote the signal and residual shot noise in the ω/t-domain, * denotes convolution in the wavenumber domain and νν _ω is given by:

[0075] Since propagating waves such as the signal /i(t,_/ ^' ) and the residual shot noise ₂ (t— ^τ ϊ'ί) always have a limited cone shaped support in the ω/c-domain bounded by the slowest

possible apparent velocity (1500m/s in marine seismic

recordings) we can solve equations (0.5) and (0.6) exactly up to a certain frequency where the number of unknowns (the F ₁ and F ₂ constituents along with the conic support constraint) are smaller than the number of equations.

[0076] Equations (0.5) and (0.6) can be solved frequency by frequency up to a maximum frequency where the number of unknowns is greater than the number of equations. That is, below that frequency, the problem is over-determined or well- posed.

[0077] Due to various loss mechanisms in real earth media, high frequencies decay rapidly with time in recorded

seismograms. Residual shot noise is therefore predominantly a problem in the low frequency part of the spectrum where the current invention is most effective.

[0078] Depending on the magnitude of the time shifts, the convolution filter described by equation (0.6) will shift more or less energy from the central signal cone. The larger the time shifts, the more energy will be shifted and the

estimation of the residual shot noise at low frequencies will be better conditioned, i.e., leading to a more accurate estimation of the residual shot noise at those frequencies.

[0079] However, there is another effect that counteracts the ill conditioned problem arising from smaller time shifts. At lower frequencies the cone shaped support of the propagating waves is narrower than at higher frequencies (in the limit it becomes a point at 0 Hz centered at zero wavenumber) . The narrower the cone shaped support is at low frequencies results in a better conditioned problem for estimating the residual shot noise.

[0080] The current invention can be applied to legacy data with natural dithers in shot firing time as long as the actual firing times are known. In cases where actual firing times are unknown it can be possible to estimate shot dithering times from the data themselves. For instance, by

crosscorrelating portions of neighboring traces where the residual shot noise is dominant (e.g., before the first arrival of signal in the seismic data) , it is possible to determine the time shifts from shot to shot.

[0081] Instead or in addition to removing residual shot noise, the present invention can also be used to remove seismic interference arrivals or repetitive rig noise for instance as long as the arrival time of the noise trains can be predicted in each shot record. The arrival times can for instance be determined using a similar cross-correlation procedure as described in the previous paragraph. The seismic interference noise can be caused by a seismic vessel conducting a different survey nearby. It can also be caused by a separate vessel in the same survey either attempting to shoot seismic data far away to minimize interference but with energy still coming through. Seismic interference can also be simultaneous source acquisition operations using for instance random dithering encoding techniques. The represent invention can therefore be used to separate the response from different sources in such simultaneous source operations in a limited frequency band.

[0082] Finally, we note that residual shot noise is signal from a previous shot. The current invention can therefore be used to reduce the time from shot to shot such that the lower most records are removed from the next shot and then

concatenated at the end of the shot to where they belong. Again, a benefit is that the high frequencies typically are not present in the lower parts of the record due to

attenuation in the Earth.

[0083] A particular application of interest is the recording of shear waves and converted waves in seabed seismic or land seismic data for instance. These arrivals propagate slower and shear wave data acquisition is therefore considered expensive as it is necessary to wait until the shear arrivals have arrived before exciting a new shot. The present

invention enables cost-effective acquisition of shear wave data where these data are recorded as residual shot noise on the next shot, separated from the next shot using the present invention and the concatenated to the shot where they belong.

Example

[0084] A synthetic example was created using an acoustic 3D finite-difference synthetic data set mimicking a seismic acquisition geometry over a complex sub-surface model. For simplicity the example is limited to the effect of a single line of shots being recorded on a receiver. The data are shown in Fig. 2 and is free from residual shot noise and therefore serves as our reference data set.

[0085] Fig. 3 shows the data in Fig. 2 after adding residual shot noise from the previous shot. The residual shot noise comprises the data between 5.41s shifted with a random time shift between -0.1s and 0.1s mimicking the effect of a marine seismic survey shooting on position where the vessel speed varies randomly resulting in a +/-0.1s variation in time when the vessel reaches the position where the next shot is fired.

[0086] Fig. 4 shows the ω/t-spectrum of the original data without the residual shot noise shown in Fig. 2.

[0087] Fig. 5 shows the difference between the ω/c-spectra of the residual shot noise (which is known exactly in this example) and the reconstructed residual noise using the invention described herein applied to the data shown in Fig. 3. It is clear that the method has successfully predicted the residual shot noise below a frequency around 14Hz. Above this frequency the method fails to reconstruct the residual noise free data as the problem is under-determined with more

unknowns than equations (the cone constraint is too wide to allow for the unambiguous separation of signal from residual shot noise) . However, in applying the current invention one simply stops removing residual shot noise above this

frequency .

[0088] The example presented in this disclosure was based on using synthetic data. In real data on the other hand the residual shot noise will contain mostly low frequencies as propagation through the Earth leads to attenuation of higher frequencies with time. Residual shot noise is therefore predominantly a problem in the low frequency part of the spectrum where the current invention is most effective.

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

List of cited References

[Abma et al . , 2012] Abma, R., Q. Zhang, A., Arogunmati, and G. Beaudoin, 2012, An overview of BP' s marine Independent

Simultaneous Source field trials: 82nd Annual International Meeting, SEG, Expanded Abstracts.

[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) .

[Akerberg et al . , 2008] Akerberg, P., Hampson, G., Rickett, J., Martin, H., and Cole, J., 2008, Simultaneous source separation by sparse Radon transform: 78th Annual

International Meeting, SEG, Expanded Abstracts, 2801-2805, doi : 10.1190/1.3063927.

[Andersson et al . , 2017] Andersson, F., Robertsson, J. 0. A., van Manen, D. J., Wittsten, J., Eggenberger, K., and Amundsen, L., 2017, Express Letters: Flawless diamond separation in simultaneous source acquisition by seismic apparition:

Geophysical Journal International, 209 (3), 1735-1739. DOI: https : //doi , org/10.1093/gi l/gqxl24

[Beasley et al . , 1998] Beasley, C. J., Chambers, R. E., and Jiang, Z., 1998, A new look at simultaneous sources: 68th Annual International Meeting, SEG, Expanded Abstracts, 133- 136.

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

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

US20140278119A1 (2014) . [Hager, 2016] Hager, E., 2016, Marine Seismic Data: Faster, Better, Cheaper?: GeoExpro, Vol. 13.

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

[Jiang and Abma, 2010] Jiang, Z. and R. Abma, 2010, An

analysis on the simultaneous imaging of simultaneous source data: 80th Annual International Meeting, SEG, Expanded

Abstracts .

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

[Lynn et al . , 1987] Lynn, W., Doyle, M., Larner, K., and

Marschall, R., 1987, Experimental investigation of

interference from other seismic crews: Geophysics, 52, 1501- 1524.

[Moldoveanu et al . , 2008] Moldoveanu, N.,Kapoor, J., and Egan, M., 2008, Full-azimuth imaging using circular geometry

acquisition: The Leading Edge, 27(7), 908-913. doi:

10.1190/1.2954032

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

(2015) .

[Robertsson et al . , 2012] Robertsson, J. 0. A., Halliday, D., van Manen, D. J., Vasconcelos, I., Laws, R., Ozdemir, K., and Gr0naas, H., 2012, Full-wavefield, towed-marine seismic acquisition and applications: 74th Conference and Exhibition, EAGE, Extended Abstracts.

[Robertsson et al . , 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. [Robertsson et al . , 2016] Robertsson, J. 0. A., Amundsen, L., and Pedersen, A. S., 2016, Express Letter: Signal apparition for simultaneous source wavefield separation: Geophys . J.

Int., 206(2), 1301-1305: doi: 10.1093/g i/ggw210.

[Shipilova et al . , 2016] Shipilova, E., Barone, I., Boelle, J. L., Giboli, M., Piazza, J. L., Hugonnet, P., and Dupinet, C., 2016, Simultaneous-source seismic acquisitions: do they allow reservoir characterization? A feasibility study with blended onshore real data: 86th Annual International Meeting, SEG, Expanded Abstracts.

[Stefani et al . , 2007] Stefani, J., Hampson, G., and

Herkenhoff, E. F., 2007, Acquisition using simultaneous sources: 69th Annual International Conference and Exhibition, EAGE, Extended Abstracts, B006.

[Ziolkowski, 1987] Ziolkowski, A. M., 1987, The determination of the far-field signature of an interacting array of marine seismic sources from near-field measurements: Results from the Delft Air Gun experiment: First Break, 5, 15-29.