Method for seismic data acquisition and processing

Abstract

Methods for separating the unknown contributions of two or more sources from a commonly acquired set of wavefield signals based on varying parameters at the firing time, location and/or depth of the individual sources in a lateral 2D plane.

Claims

1. A method for separating sources, comprising: obtaining, by at least one receiver, wavefield recordings of an underlying continuous wavefield that has conic support in a temporal-spatial frequency-wavenumber domain, based on a synchronous activation of at least two separate sources moved from shotpoint to shotpoint along a first spatial direction, and repeating the activation along two or more activation lines, wherein the two or more activation lines are distributed along a second spatial direction and each activation line comprises separate lines corresponding to the at least two separate sources; varying at least one parameter between the at least two sources from one activation to a following activation along at least one of the first and second spatial directions, the at least one parameter being at least one of source signal amplitude, source signal spectrum, source activation time, source location at activation time, and source depth, such that the varying causes support of the frequency-wavenumber representation of the obtained wavefield to be a union of at least two shifted replicas of the support of the frequency-wavenumber representation of the continuous wavefield in at least one of the first and second spatial directions, wherein the frequency-wavenumber representation of information contained in the continuous wavefield is distributed, by means of the varying, to the obtained wavefield so that a periodicity of the obtained wavefield is different than a periodicity as determined by the sampling parameters when the varying is omitted; using knowledge of the varying, at least partially resolving ambiguity present when the varying is omitted, to generate separated contributions of at least one of the at least two sources; generating subsurface representations of structures or Earth media properties using the separated contributions of the at least one of the at least two sources; and outputting the generated subsurface representations.

2. The method of claim 1, further comprising reconstructing the separated contributions of the at least one of the at least two sources up to a temporal frequency that is higher than a corresponding temporal frequency when the varying is omitted.

3. The method of claim 1, further comprising reconstructing the separated contributions of the at least one of the at least two sources up to a temporal frequency that is higher than a corresponding temporal frequency when considering only data along a single activation line.

4. The method of claim 1, wherein the varying of the at least one parameter is periodic.

5. The method of claim 4, wherein the varying of the at least one parameter comprises time shifting by time shifts that are smaller than 100 ms with respect to a firing time of individual sources.

6. The method of claim 1, wherein the varying of the at least one parameter comprises superposing multiplying at least two varying functions, wherein at least one of the at least two varying functions is periodic.

7. The method of claim 1, wherein the varying of the at least one parameter is one of random, naturally random by acquisition, pseudo-random, and quasi-random.

8. The method of claim 7, wherein the varying of the at least one parameter is a time dither greater than 100 ms.

9. The method of claim 1, further comprising at least one of dealiasing and reconstructing to recover the continuous wavefield beyond the temporal frequency below which the ambiguity has been resolved.

10. The method of claim 9, wherein the at least one of dealiasing and reconstructing comprises forming an analytic part of recorded wavefield information, extracting a non-aliased representation of a part of the recorded wavefield information, forming a phase factor from a conjugate part of an analytic part of the non-aliased representation, combining the analytic part of the recorded wavefield information with the phase factor to derive an essentially non-aliased function, applying a filtering operation to the non-aliased function, and recombining the filtered non-aliased function with the non-conjugated phase factor to reconstruct a representation of essentially dealiased recorded wavefield information.

11. The method of claim 9, wherein the at least one of the dealiasing and reconstructing further comprises iteratively forming analytic parts from parts of the at least one of the dealiasing and reconstruction, and using the analytic parts in solvers of linear systems to increase a range of the parts of the at least one of the dealiasing and reconstruction.

12. The method of claim 9, wherein a quaternion analytic part is used instead of the analytic part.

13. The method of claim 1, wherein the obtaining step comprises obtaining the wavefield recordings, which are marine seismic data or seabed seismic data, wherein the sources are towed by a same vessel.

14. The method of claim 13, wherein airguns or marine vibroseis devices are used as seismic sources.

15. The method of claim 1, wherein the obtaining step comprises obtaining the wavefield recordings, which are marine seismic data or seabed seismic data, and wherein the sources are towed by at least two different vessels.

16. The method of claim 1, wherein at least one of the at least two sources is towed at a greater depth than other sources such that each source exhibits different notch frequencies to provide a broadband source effect for all of the at least two sources for wavelengths corresponding to and below a certain frequency.

17. The method of claim 1, wherein the wavefield is obtained at locations corresponding to a Bravais grid.

18. The method of claim 1, wherein a number of the activation lines is eight or more, and a number of the shot points along an activation line is eight or more.

19. The method of claim 1, wherein the varying of the at least one parameter between the at least two sources from one activation to a next is source activation time, and a difference in source activation time between at least two of the sources is smaller than 50 ms to enable constructive interference of an emitted signal of the at least two sources at frequencies below 6 Hz.

20. The method of claim 1, wherein the at least two sources are located within 750 m from each other to enable constructive interference of an emitted signal of the at least two sources below 4 Hz.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) In the following description reference is made to the attached figures, in which:

(2) FIG. 1, parts (A) and (B) illustrate how in a conventional marine seismic survey along a single shot line 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;

(3) FIG. 2 shows an illustration of a frequency plane for a receiver gather with sources shot on a 2D (x,y)-grid. The example illustrates a well-sampled case without any aliasing or interference of contributions from different sources during simultaneous source acquisition. The “*” and the “+” mark two features of interest within the signal cone;

(4) FIG. 3 shows an illustration of a frequency plane for a receiver gather with sources shot on a 2D (x,y)-grid. The example illustrates an under-sampled case where a single source was shot on a four times as sparse grid as that in FIG. 2, generating overlapping regions of replicated and shifted signal cones;

(5) FIG. 4 shows an illustration of a frequency plane for a receiver gather with sources shot on a 2D (x,y)-grid. The example illustrates a case of under-sampling in the y-direction but where apparition encoding has been carried out in the x-direction;

(6) FIG. 5 shows a zoom in of the region of interest in FIG. 3, where we would like to remove overlap between cones;

(7) FIG. 6 shows a zoom in of the region of interest in FIG. 4, where we would like to remove overlap between cones;

(8) FIG. 7 shows the fk transform of the reference data in the data example along a line corresponding to a single shotline in a common receiver gather;

(9) FIG. 8 shows the fk transform of 8 simultaneous sources in the data example along a line corresponding to a single shotline in a common receiver gather;

(10) FIG. 9 shows the fk transform of the difference between the reference data in FIG. 7 and the data in FIG. 8 decoded by the method of signal apparition along a single line;

(11) FIG. 10 shows the fk transform of the difference between the reference data in FIG. 7 and the data in FIG. 8 decoded by means of the method described in the present invention using source lines covering a 2D lateral plane;

(12) FIG. 11 shows a linlog plot of the cumulative error in FIG. 10 summed over all inline and crossline wavenumbers.

(13) FIG. 12 Shows the general practice of marine seismic surveying;

(14) FIG. 13 Shows the general practice of land seismic surveying;

(15) FIG. 14 Summarizes key steps for one embodiment of the methods disclosed herein;

(16) FIG. 15 Illustrates how the methods herein may be computer-implemented;

DETAILED DESCRIPTION

(17) The following examples may be better understood using a theoretical overview as presented below.

(18) It is herein proposed to use acquisition patterns in a plane instead of just a line to reduce ambiguity in the overlapping region thereby causing the area of exact or nearly exact separation to be larger.

(19) By a plane it is meant a portion of the (infinite) horizontal plane covering an area corresponding to the length of the shotlines and the length in the crossline direction occupied by adjacent shotlines. Note that the number of sampling locations in the inline and crossline directions will limit the resolution properties of the Fourier transform. In practice the number of points should be greater than or equal to 8 or more. Specifically, for a Bravais lattice (x.sub.1,x.sub.2)=(u.sub.1,1k.sub.1+u.sub.1,2k.sub.2,u.sub.2,1k.sub.1+u.sub.2,2k.sub.2) thus including rectangular sampling, staggered grid and hexagonal sampling, K.sub.1.sup.min≤k.sub.1<K.sub.1.sup.max, K.sub.2.sup.min≤k.sub.2<K.sub.2.sup.max, it should hold that K.sub.1.sup.max−K.sub.1.sup.min≥8 and K.sub.2.sup.max−K.sub.2.sup.min≥8.

(20) Method for separating wavefields in the overlapping region includes utilizing the analytic part or the quaternion analytic part of the wavefields. In the following we describe the procedures in one spatial coordinate although it generalizes to two coordinates or more. We will use the notation
{circumflex over (ƒ)}(ξ)=∫.sub.−∞.sup.∞ƒ(x)e.sup.−2πixξdx
for the Fourier transform.

(21) Let denote the cone ={(ω,ξ):|ω|>|ξ|}, and let denote the non-aliased (diamond shaped) set =\({(ω,ξ):|ω|>|ξ−½|}∪{(ω,ξ):|ω|>|ξ+½|}).

(22) Suppose that
d(t,j)=ƒ.sub.1(t,j)+ƒ.sub.2(t−Δ.sub.t(−1).sup.j,j)  (1)
is a known discrete sampling in x=j recorded at a first sampling interval. Note that due to this type of apparition sampling, the data d will always have aliasing effects present if the data is band unlimited in the temporal frequency direction.

(23) If ƒ.sub.1 and ƒ.sub.2 represent seismic data recorded at a certain depth, it will hold that supp()⊂ and supp()⊂. We will assume that the source locations of ƒ.sub.1 and ƒ.sub.2 are relatively close to each other. Let

(24) $D_1\left(\omega ,\xi \right)={\int }_{-\infty }^{\infty }\underset{j=-\infty }{\overset{\infty }{.Math.}}d\left(t,j\right)e^{-2\pi i\left(j\xi +t\omega \right)}\mathrm{dt}\mathrm{and}$$D_2\left(\omega ,\xi \right)={\int }_{-\infty }^{\infty }\underset{j}{.Math.}d\left(t+{{\Delta }_t\left(-1\right)}^j,j\right)e^{-2\pi i\left(j\zeta +t\omega \right)}dt.$

(25) It is shown in Andersson et al. (2016) that

(26) $\begin{array}{cc}{D}_1\left(\omega ,\zeta \right)={\Sigma }_{k=-\infty }^{\infty }\left(\omega ,\xi +k\right)+{\Sigma }_{k=-\infty }^{\infty }\left(\omega ,\xi +\frac{k}{2}\right)\frac{1}{2}\left(e^{-2\pi i{\Delta }_t\omega }+{\left(-1\right)}^k{e}^{2\pi i{\Delta }_t\omega }\right).& \left(2\right)\end{array}$

(27) For each pair of values (ω,ξ)∈, most of the terms over k in (2) vanish (and similarly for D.sub.2), which implies that (ω,ξ) and (ω,ξ) can be recovered through

(28) $\begin{array}{cc}\left(\omega ,\xi \right)=\frac{D_1\left(\omega ,\xi \right)-\cos \left(2\pi {\Delta }_t\omega \right)D_2\left(\omega ,\xi \right)}{{\sin }^2\left(2\pi {\Delta }_t\omega \right)},\left(\omega ,\xi \right)=\frac{D_2\left(\omega ,\xi \right)-\cos \left(2\pi {\Delta }_t\omega \right)D_1\left(\omega ,\xi \right)}{{\sin }^2\left(2\pi {\Delta }_t\omega \right)},& \left(3\right)\end{array}$
given that sin(2πΔ.sub.tω)≠0.

(29) By including an amplitude variation in (1), the last condition can be removed. For values of (ω,ξ).Math.\ it is not possible to determine the values of (ω,ξ) and (ω,ξ) without imposing further conditions on the data.

(30) Given a real valued function ƒ with zero average, let
ƒ.sup.a(t)=2∫.sub.0.sup.∞ƒ.sup.∞.sub.−∞ƒ(t′)e.sup.2πi(t−′)ωdt′dω.

(31) The quantity is often referred to as the analytic part of ƒ, a description that is natural when considering Fourier series expansions in the same fashion and comparing these to power series expansions of holomorphic functions. It is readily verified that Re(ƒ.sup.a)=ƒ.

(32) As an illustrative example, consider the case where
ƒ(t)=cos(2πt)
for which it holds that
ƒ.sup.a(t)=e.sup.2πit.

(33) Now, whereas |ƒ(t)| is oscillating, |ƒ.sup.a(t)|=1, i.e., it has constant amplitude. In terms of aliasing, it can often be the case that a sampled version of ≡ƒ.sup.a| exhibits no aliasing even if ƒ and |ƒ| do so.

(34) Let us now turn our focus back to the problem of recovering (ω,ξ) and (ω,ξ) for (ω,ξ).Math.. We note that due to linearity, it holds that
d.sup.a(t,j)=ƒ.sub.1.sup.a(t,j)+ƒ.sub.2.sup.a(t−Δ.sub.t(−1).sup.j,j).

(35) A natural condition to impose is that the local directionality is preserved through the frequency range. The simplest such case is when ƒ.sub.1 and ƒ.sub.2 are plane waves (with the same direction), i.e, when
ƒ.sub.1(t,x)=h.sub.1(t+bx), and ƒ.sub.2(t,x)=h.sub.2(t+bx).

(36) Without loss of generality, we assume that b>0. We note that

(37) $\begin{array}{c}{\hat{f}}_1\left(\omega ,\xi \right)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{h}_1\left(t+bx\right)e^{-2\pi i\left(t\omega +x\zeta \right)}d\mathrm{tdx}\\ =\left\{s=t+bx\right\}\\ ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{h}_1\left(s\right)e^{-2\pi i(s\omega +x\left(\xi -b\omega \right)}\mathrm{dsdx}\\ =\left(\omega \right)\delta \left(\xi -b\omega \right).\end{array}$
A similar formula holds for .

(38) Let us now assume that ω<½. Inspecting (2) we see that if, e.g., −½<ξ<0 then all but three terms disappear and therefore the blended data satisfies

(39) $\begin{array}{c}\left(\omega ,\xi \right)=(\left(\omega \right)+\cos \left(2\pi {\Delta }_t\omega \right)\left(\omega \right))\delta \left(\xi -b\omega \right)-\\ i\sin \left(2\pi {\Delta }_t\omega \right)\left(\omega \right)\delta \left(\xi -\left(b\omega -1\text{/}2\right)\right)\\ =\left(\omega \right)\delta \left(\xi -b\omega \right)+\left(\omega \right)\delta \left(\xi -\left(b\omega -1\text{/}2\right)\right).\end{array}$

(40) Let w.sub.h, and w.sub.l be two filters (acting on the time variable) such that w.sub.h has unit L.sup.2 norm, and such that w.sub.h has a central frequency of ω.sub.0 and ω.sub.l has a central frequency of ω.sub.0/2. For the sake of transparency let w.sub.h=w.sub.l.sup.2. Suppose that we have knowledge of (ω,ξ) and (ω,ξ) for ω<ω.sub.0, and that the bandwidth of w.sub.l is smaller than ω.sub.0/2.

(41) Let g.sub.1=ƒ.sub.1.sup.a*w.sub.l and g.sub.2=ƒ.sub.2.sup.a*w.sub.l. Note that, e.g.,
g.sub.1(t,x)=h.sub.1.sup.a*w.sub.l(t+bx),

(42) so that g.sub.1 is a plane wave with the same direction as ƒ.sub.1. Moreover, |g.sub.1| will typically be mildly oscillating even when ƒ.sub.1 and |ƒ.sub.1| are oscillating rapidly.

(43) Let

(44) $p_1=\frac{g_1}{.Math.g_1.Math.}$
be the phase function associated with g.sub.1, and define p.sub.2 in a likewise manner as phase function for g.sub.2. If w.sub.l is narrowbanded, g.sub.1 will essentially only contain oscillations of the form

(45) 0$\left(\frac{{\omega }_0}{2}\right)e^{2\pi i\frac{{\omega }_0}{2}\left(t+bx\right)}$
i.e., |g.sub.1(t,x)| is more or less constant.

(46) Under the narrowband assumption on w.sub.l (and the relation w.sub.h=w.sub.l.sup.2), we consider
d.sub.h(t,x)=(d.sup.a*w.sub.h)(t,x)≈h.sub.1(ω.sub.0)e.sup.2πiω.sup.0.sup.(t+bx)+h.sub.2(ω.sub.0)e.sup.2πi(ω.sup.0.sup.t+(bω.sup.0.sup.−1)x).

(47) By multiplication by p.sub.1p.sub.2, which may be considered as an example of a conjugated phase factor, we get
d.sub.h(t,x)=p.sub.1(t,x)p.sub.2(t,x)≈h.sub.1(ω.sub.0)+h.sub.2(ω.sub.0)e.sup.−πix.  (4)

(48) In the Fourier domain, this amounts to two delta-functions; one centered at the origin and one centered at (0,−½). Here, we may identify the contribution that comes from ƒ.sub.2 by inspecting the coefficient in front of the delta-function centered at (0,−½). By the aid of the low-frequency reconstructions g.sub.1 and g.sub.2, it is thus possible to move the energy that originates from the two sources so that one part is moved to the center, and one part is moved to the Nyquist wavenumber. Note that it is critical to use the analytic part of the data to obtain this result. If the contributions from the two parts can be isolated from each other, it allows for a recovery of the two parts in the same fashion as in (3). Moreover, as the data in the isolated central centers is comparatively oversampled, a reconstruction can be obtained at a higher resolution than the original data. Afterwards, the effect of the phase factor can easily be reversed, and hence a reconstruction at a finer sampling, i.e. at a smaller sampling interval than the original can be obtained.

(49) A similar argument will hold in the case when the filters w.sub.l and w.sub.h have broader bandwidth. By making the approximation that
p.sub.1(t,x)≈e.sup.πiω.sup.0.sup.(t+bx)
we get that
d.sub.h(t,x)=p.sub.1(t,x)p.sub.2(t,x)≈h.sub.1(ω)e.sup.2πi(ω−ω.sup.0.sup.)(t+bx)+h.sub.2(ω)e.sup.−πix+e.sup.2πi(ω−ω.sup.0.sup.)(t+bx),
and since w.sub.h is a bandpass filter, (ω,ξ) will contain information around the same two energy centers, where the center around (0,−½) will contain only information about ƒ.sub.2. By suppressing such localized energy, we can therefore extract only the contribution from ƒ.sub.2 and likewise for ƒ.sub.1.

(50) The above procedure can now be repeated with ω.sub.0 replaced by ω.sub.1=βω.sub.0 for some β>1. In this fashion we can gradually recover (dealias) more of the data by stepping up in frequency. We can also treat more general cases. As a first improvement, we replace the plane wave assumption by a locally plane wave assumption, i.e., let φ.sub.α be a partition of unity (Σ.sub.αφ.sub.α.sup.2=1), and assume that
ƒ.sub.1(t,x)φ.sub.α.sup.2(t,x)≈h.sub.1,α(t+b.sub.ax)φ.sub.α.sup.2(t,x).

(51) In this case the phase functions will also be locally plane waves, and since they are applied multiplicatively on the space-time side, the effect of (4) will still be that energy will be injected in the frequency domain towards the two centers at the origin and the Nyquist wavenumber.

(52) Now, in places where the locally plane wave assumption does not hold, the above procedure will not work. This is because as the phase function contains contributions from several directions at the same location, the effect of the multiplication in will no longer correspond to injecting the energy of D.sub.1 (and D.sub.2) towards the two centers around the origin and the Nyquist number. However, some of this directional ambiguity can still be resolved.

(53) In fact, upon inspection of (2), it is clear that the energy contributions to a region with center (ω.sub.0,ξ.sub.0) must originate from regions with centers at (ω.sub.0,ξ.sub.0+k/2) for some k. Hence, the directional ambiguity is locally limited to contributions from certain directions. We will now construct a setup with filters that will make use of this fact.

(54) Let us consider the problem where we would like to recover the information from a certain region (ω.sub.0,ξ.sub.0). Due to the assumption that ƒ.sub.1 and ƒ.sub.2 correspond to measurements that take place close to each other, we assume that ƒ.sub.1 and ƒ.sub.2 have similar local directionality structure. From (2) we know that energy centered at (ω.sub.0,ξ.sub.0) will be visible in measurements D.sub.1 at locations (ω.sub.0,ξ.sub.0+k/2). We therefor construct a (space-time) filter w.sub.h, that satisfies

(55) $\begin{array}{cc}\left(\omega ,\xi \right)={\Sigma }_k\hat{\Psi }\left(\omega -{\omega }_0,\xi -{\xi }_0-\frac{k}{2}\right).& \left(5\right)\end{array}$

(56) We now want to follow a similar construction for the filter w.sub.l. Assuming that there is locally only a contribution from a direction associated with one of the terms over k above, we want the action of multiplying with the square of the local phase to correspond to a filtration using the part of w.sub.h that corresponds to that particular k.

(57) This is accomplished by letting

(58) $\begin{array}{cc}\left(\omega ,\xi \right)={\Sigma }_k\hat{\Psi }\left(\omega -\frac{{\omega }_0}{2},\xi -\frac{{\xi }_0}{2}-\frac{k}{4}\right),& \left(6\right)\end{array}$
where {circumflex over (Ψ)}={circumflex over (ψ)}*{circumflex over (ψ)}.

(59) Under the assumption that ƒ.sub.1*w.sub.h and ƒ.sub.2*w.sub.j has a local plane wave structure, we may now follow the above procedure to recover these parts of ƒ.sub.1 and ƒ.sub.2 (by suppressing localized energy as described above). We may then completely recover ƒ.sub.1 and ƒ.sub.2 up to the temporal frequency ω.sub.0 by combining several such reconstructions, and hence we may proceed by making gradual reconstructions in the ω variable as before.

(60) Let be the quaternion algebra (Hamilton, 1844). An element g∈ can be represented as q=q.sub.0+iq.sub.1+jq.sub.2+kq.sub.3, where the q.sub.j are real numbers and i.sup.2=j.sup.2=k.sup.2=ijk=1. We also recall Euler's formula, valid for i,j,k:
e.sup.iθ=cos θ+i sin θ,e.sup.jθ=cos θ+j sin θ,e.sup.kθ=cos θ+k sin θ.

(61) Note that although i,j,k commute with the reals, quaternions do not commute in general. For example, we generally have e.sup.iθe.sup.jϕ≠e.sup.jϕe.sup.iθ which can easily be seen by using Euler's formula. Also recall that the conjugate of q=q.sub.0+iq.sub.1+jq.sub.2+kq.sub.3 is the element q*=q.sub.0−iq.sub.1−jq.sub.2−kq.sub.3. The norm of q is defined as ∥q∥=(qq*).sup.1/2=(q.sub.0.sup.2+q.sub.1.sup.2+q.sub.2.sup.2+q.sub.3.sup.2).sup.1/2.

(62) Given a real valued function ƒ=ƒ(t,x), we define the quaternion Fourier transform (QFT) of ƒ by
ƒ(ω,ξ)=∫.sup.∞.sub.−∫.sup.∞.sub.−∞e.sup.−2πitωƒ(t,x)e.sup.−2πjxξdtdx.
Its inverse is given by
ƒ(t,x)=.sup.−1(ƒ)(t,x)=∫.sup.∞.sub.−∫.sup.∞.sub.−∞e.sup.−2πitωƒ(ω,ξ)e.sup.2πjxξdωdξ.

(63) In a similar fashion, it is possible to extend the Fourier transform to other hypercomplex representations, e.g., octanions (van der Blij, 1961), sedenions (Smith, 1995) or other Cayley or Clifford algebras. A similar argument applies to other well-known transform domains (e.g., Radon, Gabor, curvelet, etc.).

(64) Let

(65) $\chi \left(\omega ,\xi \right)=\{\begin{array}{cc}4& \mathrm{if}\omega >0\mathrm{and}\xi >0,\\ 0& \mathrm{otherwise}.\end{array}$

(66) Using χ we define ƒ.sup.q:.sup.2.fwdarw. as
ƒ.sup.q=.sup.−1χƒ.

(67) We call ƒ.sup.q the quaternion part of ƒ. This quantity can be seen as a generalization of the concept of analytic part. For the analytic part, half of the spectrum is redundant. For the case of quaternions, three quarters of the data is redundant.

(68) In a similar fashion, it is possible to extend the analytic part to other hypercomplex representations, e.g., octanions (van der Blij, 1961), sedenions (Smith, 1995) or other Cayley or Clifford algebras.

(69) The following results will prove to be important: Let ƒ(t,x)=cos u, where u=2π(at +bx)+c with a>0. If b>0 then
ƒ.sup.q(t,x)=cos u+i sin u+j sin u−k cos u,
and if b<0 then
ƒ.sup.q(t,x)=cos u+i sin u−j sin u+k cos u.

(70) The result is straightforward to derive using the quaternion counterpart of Euler's formula. Note that whereas |ƒ(t,x)| is oscillating, ∥ƒ.sup.q∥=√{square root over (2)}, i.e., it has constant amplitude. In terms of aliasing, it can often be the case that a sampled version of ∥ƒ.sup.q∥ exhibits no aliasing even if ƒ and |ƒ| do so.

(71) Assume that ƒ(t,x)=cos(u), and that g(t,x)=cos(ν), where u=2π(a.sub.1t+b.sub.1x)+e.sub.1 and ν=2π(a.sub.2t+b.sub.2x)+e.sub.2 with a.sub.1,a.sub.2≥0. It then holds that

(72) $\stackrel{\_}{g^q}{f}^q\stackrel{\_}{g^q}=\{\begin{array}{cc}2\left(\cos \left(u-2v\right)\left(1+k\right)+\sin \left(u-2v\right)\left(i+j\right)\right)& \mathrm{if}{b}_1,b_2\ge 0,\\ 2\left(\cos \left(u\right)\left(-1-k\right)+\sin \left(u\right)\left(i+j\right)\right)& \mathrm{if}{b}_1\ge 0\mathrm{and}{b}_2<0,\end{array}$
with similar expressions if b.sub.1<0.

(73) Let us describe how to recover (ω,ξ) and (ω,ξ) us the quaternion part. We will follow the same procedure as before, and hence it suffices to consider the case where ƒ.sub.1=h.sub.1(t+bx), ƒ.sub.2=h.sub.2(t+bx), with b>0, and ω<½. Let w.sub.h and w.sub.l be two (real-valued) narrowband filters with central frequencies of ω.sub.0 and ω.sub.0/2, respectively, as before. From (2), it now follows that
d*w.sub.h(t,x)≈c.sub.1 cos(2πω.sub.0(t+bx)+e.sub.1)+c.sub.2 cos(2πω.sub.0(t+{tilde over (b)}x)+e.sub.2),

(74) for some coefficients c.sub.1,c.sub.2, and phases e.sub.1,e.sub.2, and with {tilde over (b)}<0.

(75) Since ƒ.sub.1 and ƒ.sub.2 are known for ω=ω.sub.0/2, we let
g.sub.1(t,x)=ƒ.sub.1*w.sub.l(t,x)≈c.sub.3 cos(πω.sub.0(t+bx)+e.sub.3).
We compute the quaternion part g.sub.1.sup.q of g.sub.1, and construct the phase function associated with it as

(76) $p_1=\frac{g_1^q}{.Math.g_1^q.Math.},$
and define p.sub.2 in a likewise manner as the phase function for g.sub.2.

(77) Let d.sup.q be the quaternion part of d*w.sub.h. It then holds that after a left and right multiplication of conjugate of the phase factors p.sub.1 and p.sub.2

(78) $\stackrel{\_}{p_1}{d}^q\stackrel{\_}{p_1}\approx c_{1}2\left(\cos \left(e_1-2e_3\right)\left(1+k\right)+\sin \left(e_1-2e_3\right)\left(i+j\right)\right)+c_{2}2\left(\cos \left(2\pi {\omega }_0\left(t+\tilde{b}x\right)+e_1\right)\left(-1-k\right)+\sin \left(2\pi {\omega }_0\left(t+\tilde{b}x\right)+e_1\right)\left(i+j\right)\right).$

(79) This result is remarkable, since the unaliased part of d is moved to the center, while the aliased part remains intact. Hence, it allows for a distinct separation between the two contributing parts.

(80) It is herein proposed to increase the maximum frequency for fully recovering the contributions in recorded a wavefield from individual sources in a simultaneous source survey. More precisely, this is achieved by considering sampling the same wavefield along a plurality of parallel sail-lines to unambiguously delineate regions of overlap in a Fourier-wavenumber plane of at least two spatial dimensions.

(81) Let us begin by introducing notation and recapitulating the theory for regular seismic apparition. We will use the notation
{circumflex over (ƒ)}(ξ)=∫.sup.∞.sub.−∞ƒ(x)e.sup.−2πixξdx
for the Fourier transform in one variable, and consequently {circumflex over (ƒ)}(ω,ξ.sub.1,ξ.sub.2) for the Fourier transform of two dimensional function ƒ(t,x.sub.1,x.sub.2) with a time (t) and spatial (x.sub.1,x.sub.2) dependence. When referring to Fourier transforms with respect to only part of the variables we will use the {circumflex over ( )} notation and use the variable names to implicitly specify with regards to which dimensions the transform is applied.

(82) Suppose that ƒ(t,x.sub.1,x.sub.2) is a function with frequency support in two cones of the form

(83) $\frac{{\omega }^2}{c^2}\ge \left({\xi }_1^2+{\xi }_2\right).$

(84) The constraint comes from assuming that the function ƒ represent the recording of a wavefield at time t at a fixed receiver coordinate, source coordinates (x.sub.1,x.sub.2), and fixed depth, where the recorded wave field at this depth is a solution to the homogeneous wave equation with a velocity c. The wavefields are generated at discrete values of (x.sub.1,x.sub.2) which we assume to be equally spaced with sampling distances Δ.sub.x.sub.1 and Δ.sub.x.sub.2 meters. For the sake of transparency we will assume the two spacing coordinates to be equal with a distance set to 1 in the presentation that follows, i.e. Δ.sub.x.sub.1=Δ.sub.x.sub.2=1. The derived expressions hold in a more general setting, and it is fairly straightforward to incorporaty kind of Bravais sampling lattice (x.sub.1,x.sub.2)=(u.sub.1,1k.sub.1+u.sub.1,2k.sub.2,u.sub.2,1k.sub.1+u.sub.2,2k.sub.2), thus including rectangular sampling, staggered grid and hexagonal sampling.

(85) We now consider the case where several sources are used simultaneously. In this presentation we assume that the M sources are aligned sequentially, although the theory can easily be adjusted to cover more general sampling formations. Moreover, in this presentation, we will assume that the different sources are fired with a small internal time delays, that vary with period M. This variation can be of more general form, for instance varying with a delay pattern that includes a heteroscale variation (several superimposed or products of periodic or non-periodic functions), random or quasi-random time delay patterns or variations in source amplitude for instance. Let us denote the periodic delay pattern by Δ.sub.m.sub.1.sub.,m.sub.2. One formulation for the representation described above is

(86) $d\left(t,M{k}_1+m_1,k_2\right)=\underset{m_2=1}{\overset{M}{.Math.}}f\left(t-{\Delta }_{m_1,m_2},\left(M{k}_1+m_1\right){\Delta }_x,\left(M{k}_2+m_2\right){\Delta }_x\right)=\underset{m_2=1}{\overset{M}{.Math.}}{ℱ}^{*}\left(\hat{f}\left(\omega ,{\xi }_1,{\xi }_2\right)e^{-2\pi i{\Delta }_{m_1,m_2}\omega }\right)\left(t,\left(M{k}_1+m_1\right){\Delta }_x,\left(M{k}_2+m_2\right){\Delta }_x\right),$

(87) The Poisson sum formula
Σ.sub.k=−∞.sup.∞ƒ(k)=Σ.sub.k=−∞.sup.∞ƒ(k)
can be modified (by using the Fourier modulation rule) as

(88) $\underset{k=-\infty }{\overset{\infty }{.Math.}}f\left(k\right)e^{-2\pi i\xi k}=\underset{k=-\infty }{\overset{\infty }{.Math.}}ℱ\left(f(.Math.)e^{-2\pi i\xi .Math.}\right)\left(k\right)=\underset{k=-\infty }{\overset{\infty }{.Math.}}\hat{f}\left(\xi +k\right).$

(89) By the standard properties of the Fourier transform it is straightforward to show that

(90) 0$\underset{k=-\infty }{\overset{\infty }{.Math.}}f\left(k{\Delta }_x+x_0\right)e^{-2\pi i\left(k{\Delta }_x+x_0\right)\xi }=\frac{1}{{\Delta }_x}\underset{k=-\infty }{\overset{\infty }{.Math.}}\hat{f}\left(\xi +\frac{k}{{\Delta }_x}\right)e^{2\pi i{x}_0\frac{k}{{\Delta }_x}}$
hold for the spatial sampling parameter Δ.sub.x and some fixed spatial shift x.sub.0.

(91) Let us now consider Fourier representation of the data. We can assume that the sampling in time is fine enough to use the continuous Fourier transform for the temporal part. For the spatial parts, we will expand the data by means of Fourier series representation, thus obtained a representation that is periodic with respect to the continuous frequency-wavenumber parameters ξ.sub.1 and ξ.sub.2.

(92) The transform thus takes the shape

(93) $D\left(\omega ,{\xi }_1,{\xi }_2\right)={\int }_{-\infty }^{\infty }\underset{k_1=-\infty }{\overset{\infty }{.Math.}}\underset{k_2=-\infty }{\overset{\infty }{.Math.}}d\left(t,k_1,k_2\right)e^{-2\pi i\left(\left(k_1{\xi }_1+k_2{\xi }_{2}M\right)+t\omega \right)}\mathrm{dt}={\int }_{-\infty }^{\infty }\underset{m_1=1}{\overset{M}{.Math.}}\underset{k_1=-\infty }{\overset{\infty }{.Math.}}\underset{k_2=-\infty }{\overset{\infty }{.Math.}}d\left(t,{\mathrm{Mk}}_1+m_1,k_2\right)e^{-2\pi i\left(\left({\mathrm{Mk}}_1+m_1\right){\xi }_1+{\mathrm{Mk}}_2{\xi }_{2}t\omega \right)}\mathrm{dt}={\int }_{-\infty }^{\infty }\underset{m_1=1}{\overset{M}{.Math.}}\underset{m_2=1}{\overset{M}{.Math.}}\underset{k_1=-\infty }{\overset{\infty }{.Math.}}\underset{k_2=-\infty }{\overset{\infty }{.Math.}}f\left(t,{\mathrm{Mk}}_1+m_1,{\mathrm{Mk}}_2+m_2\right)\times e^{-2\pi i\left(\left({\mathrm{Mk}}_1+m_1\right){\xi }_1+\left({\mathrm{Mk}}_2\right){\xi }_2+t\omega +{\Delta }_{m_1,m_2}\omega \right)}\mathrm{dt}$

(94) By taking a one-dimensional Fourier transform with respect to the time variable it holds that

(95) $D\left(\omega ,{\xi }_1,{\xi }_2\right)=\underset{m_1=1}{\overset{M}{.Math.}}\underset{m_2=1}{\overset{M}{.Math.}}\underset{k_1=-\infty }{\overset{\infty }{.Math.}}\underset{k_2=-\infty }{\overset{\infty }{.Math.}}\hat{f}\left(\omega ,M{k}_1+m_1,M{k}_2+m_2\right)\times e^{-2\pi i\left(\left(M{k}_1+m_1\right){\xi }_1+\left(M{k}_2\right){\xi }_2+{\Delta }_{m_1{m}_2}\omega \right)}$

(96) Let us now focus on the part
Σ.sub.k.sub.1.sub.=−∞.sup.∞{circumflex over (ƒ)}(ω,Mk.sub.1+m.sub.1,Mk.sub.2+m.sub.2)e.sup.−2πi(Mk.sup.1.sup.+m.sup.1.sup.)ξ.sup.1

(97) By using the version of the Poisson sum formula above, we can rewrite the expression above as

(98) $\underset{k_1=-\infty }{\overset{\infty }{.Math.}}\hat{f}\left(\omega ,M{k}_1+m_1,M{k}_2+m_2\right)e^{-2\pi i\left(M{k}_1+m_1\right){\xi }_1}=\underset{k_1=-\infty }{\overset{\infty }{.Math.}}\hat{f}\left(\omega ,{\xi }_1+\frac{k_1}{M},{\mathrm{Mk}}_2+m_2\right)e^{2\pi i{m}_1{k}_1/M}$

(99) Using a similar argument for the sum over k.sub.2 we obtain

(100) $\hat{h}\left(\omega ,{\xi }_1,{\xi }_2\right)=\underset{k_1=-\infty }{\overset{\infty }{.Math.}}\underset{k_2=-\infty }{\overset{\infty }{.Math.}}\hat{f}\left(\omega ,M{k}_1+m_1,{\mathrm{Mk}}_2+m_2\right)e^{-2\pi i\left(\left(M{k}_1+m_1\right){\xi }_1+\left(M{k}_2\right){\xi }_2\right)}=\underset{k_2=-\infty }{\overset{\infty }{.Math.}}\underset{k_1=-\infty }{\overset{\infty }{.Math.}}\hat{f}\left(\omega ,{\xi }_1+\frac{k_1}{M},{\xi }_2+\frac{k_2}{M}\right)e^{2\pi i\left(\left(m_1{k}_1+m_2{k}_2\right)/M+m_2\xi \right)}$

(101) It is now important to note that the periodicity of ĥ with respect to ξ.sub.2 is (1/M) and not 1 as it is for the periodicity with respect to ξ.sub.1. To see this, note that for any integer k.sub.2′ it holds that

(102) $\hat{h}\left(\omega ,{\xi }_1,{\xi }_2+k_2^{\prime }/M\right)=\underset{k_2=-\infty }{\overset{\infty }{.Math.}}\underset{k_1=-\infty }{\overset{\infty }{.Math.}}\hat{f}\left(\omega ,{\xi }_1+\frac{k_1}{M},{\xi }_2+\frac{k_2+k_2^{\prime }}{M}\right)e^{2\pi i\left(\left(m_1{k}_1+m_2{k}_2\right)/M+m_2\left(\xi +k_2^{\prime }\right)\right)}=\underset{k_2=-\infty }{\overset{\infty }{.Math.}}\underset{k_1=-\infty }{\overset{\infty }{.Math.}}\hat{f}\left(\omega ,{\xi }_1+\frac{k_1}{M},{\xi }_2+\frac{k_2+k_2^{\prime }}{M}\right)e^{2\pi i\left(\left(m_1{k}_1+m_2\left(k_2+k_2^{\prime }\right)\right)/M+m_2\xi \right)}=\hat{h}\left(\omega ,{\xi }_1,{\xi }_2\right).$

(103) It now follows that

(104) $D\left(\omega ,{\xi }_1,{\xi }_2\right)=\underset{m_1=1}{\overset{M}{.Math.}}\underset{m_2=1}{\overset{M}{.Math.}}\hat{h}\left(\omega ,{\xi }_1,{\xi }_2\right)e^{-2\pi i{\Delta }_{m_1{m}_2}\omega }$
will have period 1 in the ξ.sub.1-direction and period 1/M in the ξ.sub.2-direction. It is now possible to choose the time delays Δ.sub.m.sub.1.sub.,m.sub.2 such that the periodicity of D (ω,ξ.sub.1,ξ.sub.2) in the ξ.sub.1 direction is larger than 1/M (where we in the definition of periodicity use a requirement of independence of ƒ, discarding degenerate cases such as ƒ=0), and typically such that the periodicity becomes 1 (using the normalized setup with Δ.sub.x=1 as mentioned earlier on).

(105) Let us now consider what happens for fixed values of ω. Since the support of {circumflex over (ƒ)} lies inside a cone whose width is determined by the spatial sampling parameters and the propagation medium velocity, it follows that the footprint of {circumflex over (ƒ)}(ω,ξ.sub.1,ξ.sub.2) will be a disk in the (ξ.sub.1,ξ.sub.2)-plane. FIG. 2 illustrates such a disk with two particular features marked “*” and “+”. From the derived expression above it follows that D(ω,ξ.sub.1,ξ.sub.2) consists of the union of several, possibly overlapping, such disks. Whereas the information that is recoverable in a regular sampling scheme would only be parts of disks that do not overlap (see FIG. 3 where the “*” and “+” features are reproduced identically in all overlapping regions), it is possible to recover information even when the disks are overlapping in the case illustrated in this application, because of the fact that there are several disks with different information content overlapping (see FIG. 4 where the “*” and “+” features are now not reproduced identically in the horizontal direction but reproduced identically in the vertical direction). This implies that it is possible to design a set of linear equation to recover information well beyond what would be feasible in a standard sampling setup. It also allows for trading sampling information in the different sampling directions.

(106) FIG. 5 and FIG. 6, show the part of the data that is used to solve for delineating the overlapping parts. In the case of regular sampling this will not be possible since the overlapping parts are all identical (FIG. 5). However, in the case of applying the sampling strategy in the inline direction as described in the present invention we see that the overlapping parts are different such that it is possible to determine where the marked features in the overlapping domains belong, provided that we know how the data were sampled in the inline direction (e.g., what time delay pattern that was used for the simultaneous source acquisition). As described above, the sampling in the inline direction can for instance be carried out using time-delays that vary periodically or randomly for instance in the inline direction. The sampling does not need to be identical in the crossline direction (the case illustrated in FIG. 4 and FIG. 6). However, it will be important to know what time delays that were used along the parallel sail-lines.

(107) One such scenario would be a vessel travelling along a sail-line, where it is possible to collect information along that line fairly densely. To collect information not only along single lines, the vessel will have to travel along several sail lines. To obtain a dense sampling, these sail lines would therefore have to be laid out close to each other, and this results in that the total measurement scenario takes a long time.

(108) Instead, the vessel could carry sources spread out in the direction perpendicular to the sail line. By using the properties described above it is possible to obtain a sampling that provides full resolution up to a certain frequency ω.sub.0 where the overlapping of disk does not prohibit perfect recovery, at an (more) equally spaced sampling lattice than what is obtained by acquiring data along individual sail lines.

(109) An important consideration for the acquisition of data that is rich in low frequencies is related to how we choose to encode the simultaneous sources relative other in for instance choosing the magnitude of the time shifts used for encoding and where the sources are located relative to each other in space.

(110) In the following for simplicity we will discuss the case of two simultaneous sources. However, the same arguments generalizes in a straight forward way to a larger number of sources, for instance with time delays smaller than the time delay considered between the two sources in the following and/or spaced at distances within the two sources considered below.

(111) The combined output of two sources at a certain frequency is A=e.sup.iωt+e.sup.iω(t+Δt), where Δ.sub.t is the time shift between the activation time of the two sources and ω is the (angular) frequency of the emitted amplitude spectrum.

(112) The interference manifests itself in an amplitude effect and a phase effect. Both are unraveled during the decoding. However, the level of the interference of the two sources is relevant when considering the signal-to-noise ratio in the acquired data.

(113) For simplicity we consider an average time delay in activation time between the two sources of 50 ms. The simple equation for the combined output shows that frequencies below 6.7 Hz benefit from constructive interference between the two sources.

(114) The spatial separation of the two sources is also important to achieve constructive interference. The phase shift between two sources that are separated by a small distance relative to the distance of propagation can be computed in a straight forward manner using Pythagoras theorem. If we simply consider a case where the distance of propagation is 5000 m in a medium with propagation velocity 1500 m/s (water), the phase shift between two sources spaced 750 m apart in a direction perpendicular to the plane of propagation is 37 ms. In such a case we will obtain constructive interference below 10 Hz. However, in the broadside direction the limit for constructive interference is lower as the phase shift becomes greater. Below 4 Hz we can expect constructive interference for all azimuths and emergence angles smaller than 15 degrees relative to the vertical which will enable the acquisition for instance of refracted arrivals for tomography and full waveform inversion propagating in the deeper part of the subsurface of interest.

EXAMPLE

(115) A synthetic example was created using an acoustic 3D finite-difference synthetic data set mimicking a seabed seismic acquisition geometry over a complex sub surface model.

(116) In the example a 25 m by 25 m shot grid was generated over a single seabed node receiver. The source had a spectrum up to 30 Hz. FIG. 7 shows the 2D Fourier transform (fk) over a line of shots in a single source acquisition mode. The data map inside a cone bounded by the water velocity.

(117) FIG. 8 shows the 2D Fourier transform (fk) over an 8 simultaneous source acquisition mode where 8 parallel shot lines were acquired simultaneously encoded using the method of signal apparition.

(118) FIG. 9 shows the difference between the reference data in FIG. 7 and the decoded data corresponding to the same source line shown in FIG. 8 for the case where the decoding was carried out in a 2D mode along a single swath line of acquisition (using the 8 parallel simultaneously acquired shot lines only) as described in the prior art (e.g., Robertsson et al., 2016). The lozenge-shaped region where the decoding is exact (or nearly exact) is empty as expected and the error is very large outside this region. FIG. 9 confirms equations (0.3) and (0.4) resulting in a lozenge-shaped region with ƒ.sub.mid=3.75 Hz and ƒ.sub.top=7.5 Hz for this case with 8 sources.

(119) FIG. 10 shows a similar difference plot as in FIG. 9 but this time using the method in the present invention where data from a large number of parallel swaths positioned in the 2D horizontal plane (3D acquisition mode) were used to decode the data from which one of the swaths is shown in FIG. 8. The region of highly accurate decoding is no longer lozenge-shaped but covers a full triangle resulting in exact (or nearly exact) decoding below a certain frequency which in this case is approximately 10.5 Hz.

(120) FIG. 11 shows a logarithmic plot of the cumulative error in FIG. 10 summed over all inline and crossline wavenumbers. The error is very low below 10.5 Hz.

(121) In FIG. 14, the key steps for one embodiment of the methods disclosed herein are summarized. In a first step, 1401, at least two different sources are encoded relative to each other for shotpoints acquired along single simultaneous acquisition lines (together comprising the first activation line) along a first spatial direction to enable the separation of the sources along a single activation line. In a second step, 1402, the at least two different sources are encoded relative to each other along adjacent activation lines distributed along a second spatial direction, by encoding the at least two different sources relative to each other in such a way that the periodicity of the obtained wavefield is different than the periodicity as determined by the sampling parameters in the case that the varying is omitted, using the encoding methods disclosed herein. In a third step, 1403, wavefield recordings are obtained along the activation lines comprising the data generated by the at least two different sources. In a fourth step, 1404, the separated contributions from the at least two different sources acquired along the set of activation lines are generated (decoded) using the methods disclosed herein and in particular using knowledge of the varying along the second spatial direction to at least partially resolve ambiguity which would be present in the case that the varying is omitted. In a fifth step, 1405, subsurface representations of structures or Earth media properties are generated using the separated contribution of at least one of the at least two sources. In a sixth step, 1406, the generated subsurface representations are output.

(122) In the description of the embodiments of the invention it is clear that acquisition lines or source lines do not necessarily need to be straight lines but can correspond to curved trajectories in the acquisition plane. It is also clear that in such a case, the term spatial directions is used to denote the path describing these trajectories.

(123) The methods described herein may be understood as a series of logical steps and (or grouped with) corresponding numerical calculations acting on suitable digital representations of the acquired seismic recordings, and hence can be implemented as computer programs or software comprising sequences of machine-readable instructions and compiled code, which, when executed on the computer produce the intended output in a suitable digital representation. More specifically, a computer program can comprise machine-readable instructions to perform the following tasks:

(124) (1) Reading all or part of a suitable digital representation of the obtained wave field quantities into memory from a (local) storage medium (e.g., disk/tape), or from a (remote) network location;

(125) (2) Repeatedly operating on the all or part of the digital representation of the obtained wave field quantities read into memory using a central processing unit (CPU), a (general purpose) graphical processing unit (GPU), or other suitable processor. As already mentioned, such operations may be of a logical nature or of an arithmetic (i.e., computational) nature. Typically the results of many intermediate operations are temporarily held in memory or, in case of memory intensive computations, stored on disk and used for subsequent operations; and

(126) (3) Outputting all or part of a suitable digital representation of the results produced when there no further instructions to execute by transferring the results from memory to a (local) storage medium (e.g., disk/tape) or a (remote) network location.

(127) Computer programs may run with or without user interaction, which takes place using input and output devices such as keyboards or a mouse and display. Users can influence the program execution based on intermediate results shown on the display or by entering suitable values for parameters that are required for the program execution. For example, in one embodiment, the user could be prompted to enter information about e.g., the average inline shot point interval or source spacing. Alternatively, such information could be extracted or computed from metadata that are routinely stored with the seismic data, including for example data stored in the so-called headers of each seismic trace.

(128) Next, a hardware description of a computer or computers used to perform the functionality of the above-described exemplary embodiments is described with reference to FIG. 15. In FIG. 15, the computer includes a CPU 1000 (an example of “processing circuitry”) that performs the processes described above. The process data and instructions may be stored in memory 1002. These processes and instructions may also be stored on a storage medium disk such as a hard drive (HDD) or portable storage medium or may be stored remotely. Further, the claimed advancements are not limited by the form of the computer-readable media on which the instructions of the inventive process are stored. For example, the instructions may be stored on CDs, DVDs, in FLASH memory, RAM, ROM, PROM, EPROM, EEPROM, hard disk or any other information processing device with which computer communicates, such as a server or another computer.

(129) Further, the claimed advancements may be provided as a utility application, background daemon, or component of an operating system, or combination thereof, executing in conjunction with CPU 1000 and an operating system such as Microsoft Windows 10, UNIX, Solaris, LINUX, Apple MAC-OS and other systems known to those skilled in the art.

(130) The hardware elements in order to achieve the computer can be realized by various circuitry elements, known to those skilled in the art. For example, CPU 1000 can be a Xenon or Core processor from Intel of America or an Opteron processor from AMD of America, or may be other processor types that would be recognized by one of ordinary skill in the art (for example so-called GPUs or GPGPUs). Alternatively, the CPU 1000 can be implemented on an FPGA, ASIC, PLD or using discrete logic circuits, as one of ordinary skill in the art would recognize. Further, CPU 1000 may be implemented as multiple processors cooperatively working in parallel to perform the instructions of the inventive processes described above.

LIST OF CITED REFERENCES

(131) [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, F., Eggenberger, K., van Manen, D. J., Robertsson, J. O. A., and Amundsen, L., 2016, Seismic apparition dealiasing using directionality regularization: 2016 SEG annual meeting, Dallas. [Andersson et al., 2017] Andersson, F., Robertsson, J. O. 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. [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. [Hamilton, 1844] W. R. Hamilton, “Ii. on quaternions; or on a new system of imaginaries in algebra.” The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 25.163: 10-13, (1844). [Ikelle, 2010] L. T. Ikelle, Coding and Decoding: Seismic Data: The Concept of Multishooting. (Elsevier, 2010), Vol. 39. [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. O. A. Robertsson, Geophysics. 80, V133 (2015). [Robertsson et al., 2012] Robertsson, J. O. A., Halliday, D., van Manen, D. J., Vasconcelos, I., Laws, R., Ozdemir, K., and Gronaas, H., 2012, Full-wavefield, towed-marine seismic acquisition and applications: 74th Conference and Exhibition, EAGE, Extended Abstracts. [Robertsson et al., 2016] Robertsson, J. O. A., Amundsen, L., and Pedersen, Å. S., 2016, Express Letter: Signal apparition for simultaneous source wavefield separation: Geophys. J. Int., 206(2), 1301-1305: doi: 10.1093/gji/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. [Smith, 1995]. J. D. H. Smith, “A left loop on the 15-sphere.” Journal of Algebra 176.1:128-138 (1995). [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. [van der Blij, F. 1961] F. van der Blij, “History of the octaves.” Simon Stevin 34 1961: 106-125. [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.