Method for seismic acquisition and processing

Abstract

A simultaneous sources seismic acquisition method is described that introduces notch diversity to improve separating the unknown contributions of one or more sources from a commonly acquired set of wavefield signals while still allowing for optimal reconstruction properties in certain diamond-shaped regions. In particular, notch diversity is obtained by heteroscale encoding.

Claims

1. A wavefield acquisition and processing method, comprising: obtaining wavefield recordings based on activation of at least two sources along one or more activation lines or grids while varying the source activation times of the sources such that at least one shifted replica of a related version of the wavefield of at least one of the sources appears in the frequency- wavenumber domain; separating a contribution of at least one of the at least two sources to the obtained wavefield recordings as generated by the at least two sources individually in absence of other sources; generating subsurface representations of structures or Earth media properties using the separated contribution of the at least one of the at least two sources; and outputting the generated subsurface representations, wherein the varying of the source activation times of the at least two sources comprises time shifts that vary from one source activation location to another source activation location and that are representable as a product of at least two separate source activation time functions corresponding to variations in the activation times on at least two different length scales along the one or more activation lines or grids.

2. The method of claim 1, wherein the varying of the source activation times on at least two different length scales are chosen before or during the acquisition.

3. The methods of claim 1, wherein the at least one shifted replica of the related version of the wavefield of the at least one of the sources is a wavefield with limited, cone-shaped support.

4. The method of claim 1, wherein the varying of the source activation times on the at least two different length scales comprises variations on a first scale including a periodic part, and variations on a second scale that introduce further notch diversity with avoidance of non-diverse notches associated with the periodic part.

5. The method of claim 4, wherein the variations on the first scale give rise to substantially non-overlapping diamond-shaped regions in the frequency-wavenumber domain and the variations on the second scale are chosen to be slow as a function of space such that the corresponding separate source activation time functions have relatively small numerical support in a spatial transform domain.

6. The method of claim 1, wherein the variations in the source activation times on the at least two different length scales enables computation of a stable inverse matrix, which is required in the separation of the contribution.

7. The method of claim 1, wherein the variations in the source activation times on the at least two different length scales enables changing, away from zero, a determinant of a matrix whose inverse is required in the separation of the contribution.

8. The method of claim 1, wherein the step of separating the contribution is carried out in the frequency-space domain.

9. The method of claim 5, wherein the variations on the second scale are smooth enough so that convolutions they cause in the frequency-wavenumber domain do not substantially change the non-overlapping diamond-shaped regions in the frequency-wavenumber associated with varying the source activation times only on a single first scale.

10. The method of claim 1, wherein the obtained wavefield recordings for multilple receivers are processed jointly in a multi-dimensional or higher-dimensional space where aliasing ambiguity is reduced.

11. The method of claim 1, wherein the separating step further comprises seismic data processing steps including at least one of imaging, migration, and demultiple for which dealiasing allows for better data quality.

12. The method of claim 1, wherein the obtained wavefield recordings comprise multiple components.

13. The method of claim 12, wherein one or more of the multiple components have a contribution separated independently, or two or more of the multiple components have a contribution separated jointly.

14. The method of claim 12, wherein one or more of the multiple components are combined before one or more products of the combination have a contribution separated.

15. The method of claim 1 applied to land seismic data, marine seismic data, seabed seismic data, permanent monitoring seismic data, time-lapse seismic data, transition zone seismic data or borehole seismic data with near surface or downhole placed receivers, and sources such as VSP, 3D VSP, or distributed acoustic sensing seismic data.

16. The method of claim 1, wherein one of the at least two different length scales corresponds to two, three, four, five, six, seven, or eight times a distance between source activation locations along one or more directions of the activation lines or grids, or a representative average of distances between source activation locations.

17. The method of claim 16, wherein a second or higher length scale of the at least two different length scales corresponds to N times a distance between source activation locations along the one or more directions of the activation lines or grids, or a representative average of distances between source activation locations, and wherein N is an integer larger than a total number of source activation points along one or more directions of the source activation lines or grids.

18. The method of claim 1, wherein the variations in the activation times on the at least two different length scales along the one or more activation lines or grids result in the source activation times being non-periodic.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

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

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

(3) FIG. 2 illustrates one embodiment of time shifts at various scales as a function of shot position;

(4) FIG. 3 illustrates heteroscale encoding sequence derived from the example in FIG. 2;

(5) FIG. 4 shows original data in the f-x domain;

(6) FIG. 5 shows reconstruction using periodic single-scale encoding;

(7) FIG. 6 illustrates reconstruction using heteroscale encoding;

(8) FIG. 7 shows the general practice of marine seismic surveying;

(9) FIG. 8 shows the general practice of land seismic surveying;

(10) FIG. 9 summarizes key steps for one embodiment of the methods disclosed herein; and

(11) FIG. 10 illustrates how the methods herein may be computer-implemented.

DETAILED DESCRIPTION

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

(13) The slowest observable (apparent) velocity of a signal along a line of recordings in any kind of wave experimentation is identical to the slowest physical propagation velocity in the medium where the recordings are made. As a result, after a spatial and temporal Fourier transform, large parts of the frequency-wavenumber (ωk) spectrum inside the Nyquist frequency and wavenumber tend to be empty. In particular, for marine reflection seismic data (Robertsson et al., 2015), the slowest observable velocity of arrivals corresponds to the propagation velocity in water (around 1500 m/s).

(14) FIG. 1, (A) part (A) illustrates how all signal energy when represented in or transformed into the frequency-wavenumber (ωk) domain sits inside a “signal cone” centered at k=0 and bounded by the propagation velocity of the recording medium.

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

(16) Furthermore, the terms “cone” and “cone-shaped” used herein are used to indicate the shape of the “bounded” or “effective numerical” support of the data of interest (e.g., the data that would be recorded firing the sources individually [i.e. non-simultaneously]) in the frequency-wavenumber domain. In many cases, it will still be possible to apply the methods described herein if the actual support is approximately conic or approximately cone-shaped. For example, at certain frequencies or across certain frequency ranges the support could be locally wider or less wide than strictly defined by a cone. Such variations are contemplated and within the scope of the appended claims. That is, the terms “cone” and “cone-shaped” should be understood to include approximately conic and approximately cone-shaped. In addition, in some cases we use the terms “bounded support” or “limited support” and “effective numerical support” to refer to data with “conic support” or “cone-shaped support” even though in the strict mathematical sense a “cone” is not bounded (as it extends to infinite temporal frequency). In such cases, the “boundedness” should be understood to refer to the support of the data along the wavenumber axis/axes, whereas “conic” refers to the overall shape of the support in the frequency-wavenumber domain.

(17) Note that the term “cone-shaped support” or similar refers to the shape of the support of e.g. the data of interest (in the frequency-wavenumber domain), if it were regularly sampled along a linear trajectory in 2D or Cartesian grid in 3D. That is, it refers only to the existence of such a support and not to the actual observed support of the data of interest in the simultaneous source input data or of the separated data of interest sampled as desired. The support of both of these depends on the chosen regularly or irregularly sampled straight or curved input (activation) and output (separation) lines or grids. Such variations are within the scope of the appended claims.

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

(19) In a wavefield experiment it may be that a source is excited sequentially for multiple source locations along a line while recording the reflected wavefield on at least one receiver. The source may be characterized by its temporal signature. In the conventional way of acquiring signals representing a wavefield the source may be excited using the same signature from source location to source location, denoted by integer n. Next, consider the alternative way of acquiring such a line of data using a periodic sequence of source signatures: every second source may have a constant signature and every other second source may have a signature which can for example be a scaled or filtered function of the first source signature. Let this scaling or convolution filter be denoted by a(t), with frequency-domain transform A(ω). Analyzed in the frequency domain, using for example a receiver gather (one receiver station measuring the response from a sequence of sources) recorded in this way, can be constructed from the following modulating function m(n) applied to a conventionally sampled and recorded set of wavefield signals:
m(n)=½[1+(−1).sup.n]+½A[1−(−1).sup.n],
which can also be written as
m(n)=½[1+e.sup.iπn]+½A[1−e.sup.iπn].  (0.1)

(20) By applying the function m in Eq. 0.1 as a modulating function to data ƒ(n) before taking a discrete Fourier transform in space (over n), F(k)=(ƒ(n)), the following result can be obtained:

(21) $\begin{array}{cc}ℱ\left(f\left(n\right)m\left(n\right)\right)=\frac{1+A}{2}F\left(k\right)+\frac{1-A}{2}F\left(k-k_N\right)=A_{+}F\left(k\right)+A_{-}F\left(k-k_N\right),& \left(0.2\right)\end{array}$
which follows from a standard Fourier transform result (wavenumber shift) (Bracewell, 1999).

(22) Eq. 0.2 shows that the recorded data ƒ will be replicated in two places in the spectral domain as illustrated in FIG. 1(B) by means of related versions of ƒ where the relationship between the original ƒ and the two related versions is due to the multiplicative operators A.sub.+ and A.sub.− quantified in Tab. I for different choices of A(ω).

(23) TABLE-US-00001 TABLE I Mapping of signal to cone centered at k = 0 (A.sub.+) and cone centered at k = k.sub.N (A.sub.−) for different choices of A(ω) for signal separation or signal apparition in Eq. (0.2). A(ω) A.sub.− = (1 − A)/2 A.sub.+ = (1 + A)/2 1 0 1 −1 1 0 0 ½ ½ ½ ¼ ¾ e.sup.iωT (1 − e.sup.iωT)/2 (1 + e.sup.iωT)/2 1 + e.sup.iωT −e.sup.iωT/2 1 + e.sup.iωT/2

(24) Part of the data will remain at the signal cone centered around k=0 (denoted by A.sub.+ in FIG. 1, part (B) (b)) and part of the data will be scaled and replicated to a signal cone centered around k.sub.N (denoted by A.sub.−). It can be observed that by only knowing one of these parts of the data it is possible to predict the other.

(25) This process may be referred to as “wavefield apparition” or “signal apparition” in the meaning of “the act of becoming visible”. In the spectral domain, the wavefield caused by the periodic source sequence is nearly “ghostly apparent” and isolated.

(26) 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 ƒ=g+h are recorded at a receiver location.

(27) In the frequency-wavenumber domain, where the recorded data are denoted by F=G+H, the H-part is partitioned into two components H.sub.+ and H.sub.− with H=H.sub.++H.sub.− where the H.sub.−-component is nearly “ghostly apparent” and partially isolated around the Nyquist-wavenumber [FIG. 1(B)], whereas G and H.sub.+ are overlapping wavefields around k=0. Furthermore, H.sub.− completely determined by H and the multiplicative operator .sub.−. This means that H can be reconstructed from H.sub.− for all values of H.sub.− such that A.sub.−(ω) is non-zero and such that H.sub.−(ω,k) and H.sub.−(ω,k-k.sub.N) do not overlap. For the partial reconstruction of, the corresponding partial reconstruction of G can be obtained through the relationship G=F−H yielding the separated partially recovered wavefields g and h in the time-space domain.

(28) Although the above description has focused on acquisition along essentially straight lines, the methodology applies equally well to curved trajectories such as coil-shaped trajectories, circles, or other smoothly varying trajectories or sequences of source activations.

(29) The concept may be extended to the simultaneous acquisition of more than two source lines by choosing different modulation functions for each source and it can be applied to higher dimensional source sampling in space (van Manen et al., 2016c).

(30) 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 A(ω) and then repeating the pattern again until the full source line has been acquired, will generate additional signal cones centered around ±k.sub.N/2.

(31) In the following, such additional scaled or related signal cones, which may occur at various “shifted” locations along the wavenumber axis, or axes (in case of higher dimensional spatial encoding), are also referred to individually as a “shifted replica of a related version of the wavefield of one of the sources”. The support of these shifted replica, may or may not be theoretically bounded and may or may not be strictly shaped as a cone as explained in more detail elsewhere herein. For completeness, the introduction of such shifted replicates of related versions at various locations along the wavenumber axis/axes is also referred to herein using the terms “injecting” and “partially injecting”. Furthermore, as should be clear from the text and Tab. I, the words “scaled” and “related” are used herein to describe the action of a multiplicative operator which can have, for example, a temporal frequency dependence.

(32) FIG. 1, (B)Part(B) also illustrates a possible limitation of signal apparition. The H.sub.+ and H.sub.− parts are separated within the respective diamond-shaped regions in FIG. 1, (B)part(B). In the triangle-shaped parts they interfere and may no longer be separately predicted without further assumptions and van Manen et al. (2016c) describe methods how to address this. In the example shown in FIG. 1, part(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.

(33) As is evident from Tab. I, the special case A=1 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.sub.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 A(ω) for wave experimentation seems to let every second source be excited a time shift T later compared to neighbouring recordings, that is, select A=e.sup.iωT.

(34) The above description assumes a simple modulating sequence m(n), and thus generating the wavefield h. In practice it is difficult to obtain perfectly periodic time shifts from a measurement setup. It is for example common practice for seismic vessels to shoot or trigger their sources at predetermined (essentially equidistant) positions, and due to practical variations (vessel velocity etc.) it will be difficult to realize shots at both predetermined locations and times.

(35) Deviations, or also called perturbations, from a nominal encoding setup in the simultaneous source acquisition, for example such as described above, can be dealt with using a variety of methods, for example such as the methods discussed in van Manen et al. [2016a].

(36) Following such deterministic and periodic shot modulation, recorded data are partially injected into a small number of “new” signal cones offset along the wavenumber axis thereby populating the frequency-wavenumber space much better compared to conventional single-source acquisition and enabling exact simultaneous source separation below a certain frequency.

(37) A critical observation and insight in the signal apparition approach is that partially injecting energy along the wavenumber axis is sufficient as long as the source variations are known as the injected energy fully predicts the energy that was left behind in the “conventional” signal cone.

(38) [Andersson et al., 2017b] demonstrate, by mathematical proof, that the signal apparition simultaneous source approach overcomes a fundamental limitation in non-periodic (e.g., random dithered) simultaneous source acquisition and is able to exactly recover a region in frequency-wavenumber space that is twice as large compared to the case of random, dithered simultaneous source acquisition. We refer herein to such a region as a “flawless diamond” or to this process and its result as “flawless reconstruction” on a “diamond-shaped” area.

(39) A known drawback, for the exemplified periodic time shift modulation functions, is that the periodicity introduces notches in the frequency spectrum which prevent direct exact recovery as outlined above, and this poses challenges to source separation and reconstruction. Thus it is an objective of the current invention to provide methods that do not suffer this drawback. To this end, a methodology for heteroscale encoding and decoding (e.g., acquisition, respectively, processing) is introduced which introduces notch diversity while allowing for flawless reconstruction in diamond-shaped regions that are only slightly smaller compared to those enabled by encoding using periodic time shifts.

(40) Let us start by recapitulating some one-dimensional properties of the Fourier transform. We will use the notation
.sub.{circumflex over (ƒ)}=∫.sub.−∞.sup.∞ƒ(x)e.sup.−2πixξdx
for the Fourier transform in one variable, and consequently .sub.{circumflex over (ƒ)}(ω,ξ) for the Fourier transform of two dimensional function ƒ(t,x) with a time (t) and spatial (x) dependence, and .sub.{circumflex over (ƒ)}.sup.t(ω,x) to denote the Fourier transform of ƒ solely with respect to the time variable, and correspondingly we will use .sub.{circumflex over (ƒ)}.sup.x(ω,x) for the Fourier transform solely with respect to the spatial variable.

(41) We will make use of the following version of the Poisson sum formula

(42) $\stackrel{\infty }{\underset{k=-\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 }^{\hat{f}}}{\overset{\infty }{.Math.}}\left(\xi +\frac{k}{{\Delta }_x}\right)e^{-2\pi i{x}_0\frac{k}{{\Delta }_x}}$
given some spatial sampling parameter Δ.sub.x and fixed spatial shift x.sub.0, [Andersson et al., 2016].

(43) Suppose that ƒ.sub.1(t,x), . . . , ƒ.sub.M(t,x) are M represent seismic data from M different sources and consider the case where sources are sampled simultaneously using the pattern

(44) $\begin{array}{cc}d\left(t,{\Delta }_{x}k\right)=\underset{n=1}{\overset{M}{.Math.}}{a}_{k,n}{f}_n\left(t-{\tau }_{k,n},{\Delta }_{x}k\right),& \left(1\right)\end{array}$
where the time shifts τ.sub.k,n are of the form
τ.sub.k,n=h(k)τ.sub.k,n.sup.p
where τ.sub.k,n.sup.p as well as the amplitude factors a.sub.k,n vary periodically with period M, i.e.,
τ.sub.k+M,n.sup.p=τ.sub.k,n.sup.p, and a.sub.k+M,n=a.sub.k,n.

(45) Let us define

(46) $D\left(\omega ,\xi \right)=\underset{k=-\infty }{\overset{\infty }{.Math.}}{\int }_{-\infty }^{\infty }d\left(t,{\Delta }_{x}k\right)e^{-2\pi i\left(t\omega +{\Delta }_{x}k\xi \right)}\mathrm{dt}.$

(47) We may rewrite this as

(48) $\begin{array}{cc}D\left(\omega ,\xi \right)=\underset{m=1}{\overset{M}{.Math.}}\underset{k=-\infty }{\overset{\infty }{.Math.}}{\int }_{-\infty }^{\infty }d\left(t,{\Delta }_x\left(\mathrm{Mk}+m\right)\right)e^{-2\pi i\left(t\omega +{\Delta }_x\left(\mathrm{Mk}+m\right)\xi \right)}=\underset{m=1}{\overset{M}{.Math.}}\underset{n=1}{\overset{M}{.Math.}}{a}_{m,n}\underset{k=-\infty }{\overset{\infty }{.Math.}}{\int }_{-\infty }^{\infty }\left(f_n\left(t-{\tau }_{m,n}^{p}h\left(\mathrm{kM}+m\right),{\Delta }_x\left(\mathrm{kM}+m\right)\right)e^{-2\pi i\left(t\omega +{\Delta }_x\left(\mathrm{Mk}+m\right)\xi \right)}\right)\mathrm{dt}& \left(2\right)\end{array}$

(49) Let us look more carefully at the inner sum above. To begin with, we may evaluate the Fourier integral to obtain terms of the form

(50) ${\hat{f}}_n^t\left(\omega ,{\Delta }_x\left(\mathrm{kM}+m\right)\right)e^{-2\pi i{\tau }_{m,n}^{p}h\left(\mathrm{kM}+m\right)}{e}^{-2\pi i\left(t\omega +{\Delta }_x\left(\mathrm{Mk}+m\right)\xi \right)}$

(51) By introducing

(52) ${\hat{g}}_{m+n}^t\left(\omega ,\mathrm{kM}+m\right)={\hat{f}}_n^t\left(\omega ,{\Delta }_x\left(\mathrm{kM}+m\right)\right)e^{-2\pi i{\tau }_{m,n}^p\left(\mathrm{kM}+m\right)}$
we can express the inner sum of using the Poisson sum formula as

(53) $\underset{k=-\infty }{\overset{\infty }{.Math.}}{\hat{g}}_{m+n}^t\left(\omega ,\mathrm{kM}+m\right)e^{-2\pi i\left(t\omega +{\Delta }_x\left(\mathrm{Mk}+m\right)\xi \right)}=\underset{k=-\infty }{\overset{\infty }{.Math.}}{\hat{g}}_{m+n}\left(\omega ,\xi +\frac{k}{M{\Delta }_x}\right)e^{-2\pi i\frac{\mathrm{mk}}{M}}=\underset{k=-\infty }{\overset{\infty }{.Math.}}\left({\hspace{0.17em}}_{f_n}^{^}\stackrel{x}{*}{e}^{-2\pi i{\tau }_{m,n}h(.Math.)}\right)\left(\omega ,\xi +\frac{k}{M{\Delta }_x}\right)e^{-2\pi i\frac{\mathrm{mk}}{M}}$

(54) We may thus conclude that

(55) $\begin{array}{cc}D\left(\omega ,\xi \right)=\underset{k=-\infty }{\overset{\infty }{.Math.}}\times \underset{n=1}{\overset{M}{.Math.}}\underset{m=1}{\overset{M}{.Math.}}\left({\hspace{0.17em}}_{f_n}^{^}\stackrel{x}{*}{e}^{-2\pi i{\tau }_{m,n}h(.Math.)}\right)\left(\omega ,\xi +\frac{k}{M{\Delta }_x}\right)e^{-2\pi i\frac{\mathrm{mk}}{M}}& \left(3\right)\end{array}$

(56) Let us for the moment assume that h(x)=1. Since the functions ƒ.sub.r1ƒ.sub.M describe seismic measurements, their Fourier transforms have support a conic domain of the form
={(ω,ξ):ω.sup.2≥c.sup.2ξ.sup.2},
where c denotes the velocity of sound in the media, e.g., water.

(57) For values of (ω,ξ) in the diamond shaped region

(58) 0$𝒟=𝒞\backslash \left(\left(𝒞-\left(0,\frac{1}{M{\Delta }_x}\right)\right).Math.\left(𝒞+\left(0,\frac{1}{M{\Delta }_x}\right)\right)\right),$
it thus holds that the supports of

(59) ${\hspace{0.17em}}_{f_n}^{^}\left(.Math.,.Math.+\frac{k}{M{\Delta }_x}\right)$
are non-intersecting for different values of k (and n fixed).

(60) This, in combination with the fact that D(ω,ξ) is periodic in ξ with a period of 1/(MΔ.sub.x), implies that if we define

(61) $D_m\left(\omega ,\xi \right)=D\left(\omega ,\xi +\frac{m}{M{\Delta }_x}\right),$
for each (ω.sub.ξ)−, there is a linear relationship between D.sub.m(ω,ξ) and .sub.{circumflex over (ƒ)}n(ω,ξ).

(62) Let us now turn our attention back to the function h. For general choices of h, the argument about non-overlapping parts of the diamond shaped regions break down. However, if h is chosen to be of slow variation (i.e. that .sub.ĥ has small numeric support), then the convolutions

(63) ${\hspace{0.17em}}_{g_n}^{^}\stackrel{x}{*}{e}^{-2\pi i{\tau }_{m,n}h(.Math.)}$
will effectively have support in domains that are only slightly larger than .

(64) It then approximately holds that

(65) $D_m\left(\omega ,\xi \right)\approx D_m^a\left(\omega ,\xi \right)=\underset{n=1}{\overset{M}{.Math.}}\underset{k=1}{\overset{M}{.Math.}}\left({\hspace{0.17em}}_{f_n}^{^}\stackrel{x}{*}{e}^{-2\pi i{\tau }_{k,n}h(.Math.)}\right)\left(\omega ,\xi \right)e^{2\pi i\frac{\mathrm{mk}}{M}}.$

(66) If we now apply a Fourier transformation with respect to the spatial variable we obtain the relationship

(67) $D_m^a\left(\omega ,x\right)=\underset{n=1}{\overset{M}{.Math.}}\underset{k=1}{\overset{M}{.Math.}}{\hat{f}}_n^t\left(\omega ,x\right)e^{-2\pi i{\tau }_{k,n}h\left(x\right)}{e}^{2\pi i\frac{\mathrm{mk}}{M}}.$

(68) Upon inspection, we can see that the above expression can be simplified if expressed in terms of matrices. Hence, let define the column vectors
D.sub.ω,x=(D.sub.1.sup.a(ω,x), . . . ,D.sub.M.sup.a(ω,x)).sup.T,
and
F.sub.ω,x=({circumflex over (ƒ)}.sub.1.sup.t(ω,x), . . . ,{circumflex over (ƒ)}.sub.M.sup.t(ω,x)).sup.T,
it holds that
D.sub.ω,x=BA.sub.ω,xF.sub.ω,x,
where A.sub.ω,x and B are M×M matrices defined by

(69) $A_{\omega ,x}\left(m,n\right)=a_{m,n}{e}^{-2\pi i\omega {\tau }_{m,n}^{p}h\left(x\right)},\mathrm{and}$$B\left(m,n\right)=e^{2\pi i\frac{\mathrm{mn}}{M}}.$
The above system is invertible as long as A.sub.ω,x is invertible, since the matrix B is a discrete Fourier matrix.

(70) Let us now (without practical limitation) assume that the time shift τ.sub.m,n.sup.p are rational numbers. They can then clearly be represented as τ.sub.m,n.sup.p=T.sub.m,nΔ.sub.T for integer values of T.sub.m,n and some Δ.sub.T. Let ω.sub.x=ωh(x)Δ.sub.T. Note that the elements in A.sub.ω,x are now periodic with respect to ω with a period defined by h(x)Δ.sub.T. A condition about invertibility can then be posed using polynomial as follows: Let P.sub.T(z) be the polynomial matrix with monomial elements

(71) ${\left(P_T\right)}_{m,n}\left(z\right)=a_{m,n}{z}^{T_{m,n}},$
and define
p.sub.T(z)=det(P.sub.T(z)).  (4)

(72) It clearly holds that

(73) $\det \left(A_{\omega ,x}\right)=p_T\left(e^{2\pi i\omega h\left(x\right){\Delta }_T}\right).$
The individual sources can now be reconstructed from the relation

(74) $F_{\omega ,x}=\frac{1}{p_T\left(e^{2\pi i\omega h\left(x\right){\Delta }_T}\right)}\mathrm{adj}\left(P_T\left(e^{2\pi i\omega h\left(x\right){\Delta }_T}\right)\right)B^{-1}{D}_{\omega ,x},$
where adj denotes the matrix adjugate.

(75) The problem is solvable as long as p.sub.T is non-zero. A special case of practical interest is when there is no amplitude variation, i.e., when a.sub.m,n=1. In this case p.sub.T does always have at least one root on the unit circle, namely at x=1. This happens for ω=0, and the polynomial matrix P.sub.T(1) then becomes a matrix with only ones, which clearly have zero determinant. Note that this also implies that there is a notch at ω=1/(h(x)Δ.sub.t).

(76) In case it is not possible for practical reasons to choose Δ.sub.x small enough to avoid non-zero notches, the above procedure shows that it is possible to introduce notch diversity by introducing a factor h that varies slowly. The sampling thus contains two types of time shifts, one that is varying rapidly and one that is varying slowly. The presence of variation in sampling delays at different scales motivates the term heteroscale encoding.

(77) A person skilled in the art will realize that it would be possible to obtain such heteroscale encoding sequences by a variety of means, e.g., by superposition of modulation functions, by varying acquisition parameters other than the source activation time (e.g., varying the relative depth of the simultaneous sources etc.). In this context, it should be mentioned that a related beneficial concept of notch diversity has been introduced previously [Soubaras and Lafet, 2013], where the depth of individual sources or receivers is varied intentionally to in order to directly control and vary the frequency of the notch introduced by sea-surface reflection. Note that in this case the notch is directly the result of the destructive interference of the upgoing and downgoing, sea-surface reflected waves and not due to variations in source activation times (or variations in the relative tow depths of the simultaneous sources, as suggested herein).

(78) A schematic outline of the reconstruction procedure is given below: a. Conduct a two dimension Fourier transform of the M-source simultaneous data to compute D with N.sub.t samples in time and N.sub.s samples in the, e.g., inline source coordinate. b. Reshape the two dimensional data D into a data cube of size N.sub.t×N.sub.s/M×M. c. Apply FFT in the third dimension, i.e., apply B.sup.−1. d. Apply FFT in the second dimensions, i.e., compute a spatial inverse Fourier transform e. Apply the adjugate of the polynomial matrix P.sub.T.

(79) f. Apply a regularized version of 1/p.sub.T, e.g.

(80) 0$\frac{\stackrel{\_}{p_T}}{\max \left({.Math.p_T.Math.}^2,\mu \right)}$ for some choice of regularization parameter μ.

(81) FIGS. 1-3 provide an example where the notch effects from different encoding patterns are illustrated. FIG. 1 shows the original data in the f-x domain, whereas FIG. 2 shows a reconstruction using periodic single-scale encoding with a single time-shift modulation function applied. FIG. 3 shows reconstruction using heteroscale encoding where the notch frequency varies with source location and which thus introduces notch diversity.

(82) The dealiasing methods introduced by van Manen et al. (2016b) and Andersson et al. (2017a) can also be applied in conjunction with the current invention and the scope of the current invention is not limited to sufficiently sampled data. Thus, for example, it is possible to introduce directionality penalties to reduce aliasing ambiguity in the separation at higher frequencies by means of estimating local directionality structure using partial reconstructions comprising only the flawlessly reconstructed diamond-shaped regions. Alternatively, for example, the flawlessly reconstructed diamond-shaped regions can be used to form conjugated phase factors from an analytic part of these regions and combined with an analytic part of the obtained wavefield recordings to derive essentially non-aliased functions that can be filtered and recombined using non-conjugated phase factors to yield essentially dealiased and separated wavefield recordings.

(83) The current invention is not limited to data acquired without deviations or perturbations in the nominal (that is, intended) heteroscale modulation functions and the methods introduced by van Manen et al. [2016a] can be applied in conjunction with the current invention when such deviations or perturbations are present on top of the nominal heteroscale modulation functions to yield improved results. For example, provided the actual heteroscale modulation functions are known (that is, including the deviations/perturbations), the Fourier transform of the actual modulation functions can be used to set up a so-called cyclic convolution matrix which models the actual scaling and replication of the signal along the wavenumber axis/axes. This can be used to set up and and invert a model of the actual obtained wavefield recordings.

(84) Provided the perturbations/deviations in the nominal heteroscale modulation functions are sufficiently small, one or more of the advantages of the current invention will at least be partially retained and may be improved by using the methods introduced by van Manen et al. [2016a], including as explained herein.

(85) The methods described herein have mainly been illustrated using so-called common receiver gathers, i.e., all seismograms recorded at a single receiver. Note however, that the method can be applied straightforwardly over one or more receiver coordinates, to individual or multiple receiver-side wavenumbers. Processing in such multi-dimensional or higher-dimensional spaces can be utilized to reduce data ambiguity due to sampling limitations of the seismic signals.

(86) We note that further advantages may derive from applying the current invention to two- or three-dimensional shot activation grids, where beyond, e.g., the (inline) x-location of the simultaneous sources, the shot grids also extend in (crossline) y-location and/or the vertical (z or depth) direction (van Manen et al., 2016c). Furthermore, the methods described herein could be applied to different two-dimensional shot activation grids, such as shot activation grids in the x-z plane or y-z plane. The vertical wavenumber is limited by the dispersion relation and hence the encoding and decoding can be applied similarly to 2D or 3D shot activation grids which involve the z (depth) dimension, including by making typical assumptions in the dispersion relation. In these and other cases, the 2D modulation functions may have one or more short scales complemented by one or more longer scales which may or may not align with the axes of the shot activation grids or be fully parallel with each other.

(87) As should be clear to one possessing ordinary skill in the art, the methods described herein apply to different types of wavefield signals recorded (simultaneously or non-simultaneously) using different types of sensors, including but not limited to; pressure and/or one or more components of the particle motion vector (where the motion can be: displacement, velocity, or acceleration) associated with compressional waves propagating in acoustic media and/or shear waves in elastic media. When multiple types of wavefield signals are recorded simultaneously and are or can be assumed (or processed) to be substantially co-located, we speak of so-called “multi-component” measurements and we may refer to the measurements corresponding to each of the different types as a “component”. Examples of multi-component measurements are the pressure and vertical component of particle velocity recorded by an ocean bottom cable or node-based seabed seismic sensor, the crossline and vertical component of particle acceleration recorded in a multi-sensor towed-marine seismic streamer, or the three component acceleration recorded by a microelectromechanical system (MEMS) sensor deployed e.g. in a land seismic survey.

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

(89) Further, it should be understood that the various features, aspects and functionality described in one or more of the individual embodiments are not limited in their applicability to the particular embodiment with which they are described, but instead can be applied, alone or in various combinations, to one or more of the other embodiments of the invention.

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

(91) In FIG. 9, the key steps for one embodiment of the methods disclosed herein are summarized. In a first step, 901, wavefield recordings are obtained based on the activation of at least two sources along one or more activation lines or grids varying the source activation times of the sources such that at least one shifted replica of a related version of the wavefield of at least one of the sources appears in the frequency-wavenumber domain and where the varying the source activation times of the sources comprise time shifts that vary from source activation location to source activation location and that are representable as a product of at least two separate source activation time functions corresponding to variations in the activation times on at least two different length scales along the one or more activation lines or grids. This is done in accordance with the general practice of marine or land seismic acquisition and/or the methods disclosed herein. In a second step, 902, A contribution of at least one of the at least two sources to the obtained wavefield recordings as generated by the at least two sources individually in the absence of the other sources is separated. In a third step, 903, 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 fourth step, 904, the generated subsurface representations are output.

(92) 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:

(93) (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;

(94) (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

(95) (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.

(96) 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.

(97) 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. 10. In FIG. 10, 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.

(98) 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.

(99) 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

(100) [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). [Andersson et al., 2016] Andersson, F., Eggenberger, K., van Manen, D. J., Robertsson, J., & Amundsen, L. (2016). Seismic apparition dealiasing using directionality regularization. In SEG Technical Program Expanded Abstracts 2016 (pp. 56-60). Society of Exploration Geophysicists. [Andersson et al., 2017a] Andersson, F., van Manen, D. J., Robertsson, J. O. A., Wittsten, J., and Eggenberger, K., 2017, Method for dealiasing data: GB patent application No. 1700520.8 filed on 12 Jan. 2017. [Andersson et al., 2017b] Andersson, F., Robertsson, J., van Manen, D. J., Wittsten, J., Eggenberger, K., and Amundsen, L., 2017, Flawless diamond reconstruction for simultaneous source separation: 79th Conference and Exhibition, EAGE, Extended Abstracts. [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. [Barnes, 1992] A. E. Barnes, GEOPHYSICS, 57(5), 749-751 (1992). [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). [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). Laws, R. M., Halliday, D. F., Özbek, A., and Hopperstad, F., 2016, Exploiting the control of phase in marine vibrators: The First Break, 34(11), pp: 65-74. [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., Özdemir, K., and Grønaas, H., 2012, Full-wavefield, towed-marine seismic acquisition and applications: 74th Conference and Exhibition, EAGE, Extended Abstracts. [Robertsson et al., 2015] Robertsson, J. O. 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.sup.nd edition (ed. G. Schubert), Elsevier-Pergamon, Oxford. [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. [Robertsson et al., 2016b], Robertsson, J. O. A., Eggenberger, K., van Manen, D. J., and Andersson, F., 2016, Simultaneous source acquisition and separation method, GB Patent application 1619037.3 filed on Nov. 10, 2016. [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. [Soubaras and Lafet, 2013] Soubaras, R., and Lafet, Y.,2013, Variable-depth streamer acquisition: Broadband data for imaging and inversion, Geophysics, 78(2), WA27-WA39. [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. [Stockwell, 1996] R. G. Stockwell, L. Mansinha, and R. P. Lowe. Signal Processing, IEEE Transactions on 44(4), 998-1001 (1996). [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. [van Manen et al., 2016a] van Manen, D. J., Andersson, F., Robertsson, J. O. A., and Eggenberger, K., 2016, Source separation method: GB Patent Application No. 1603742.6 filed on 4 Mar. 2016. [van Manen et al., 2016b] D. J. van Manen, F. Andersson, J. O. A. Robertsson, K. Eggenberger, 2016, De-aliased source separation method: GB patent application No. 1605161.7 filed on 28 Mar. 2016. [van Manen et al., 2016c] van Manen, D. J., Andersson, F., Robertsson, J. O. A., and Eggenberger, K., 2016, Simultaneous source acquisition and separation on general related sampling grids: GB patent application No. 1608297.6 filed on 12 May 2016