Method of redatuming geophysical data

Abstract

A method of redatuming geophysical data, wherein there is provided multi-component geophysical data, and the method includes obtaining at least one focussing function and/or at least one Green's function from the multi-component geophysical data.

Claims

1. A method of prospecting for and producing hydrocarbons, the method comprising: obtaining geophysical data collected using an acquisition setup comprising a plurality of geophysical receivers and a plurality of geophysical sources, wherein at least one geophysical source is not at the same location as at least one geophysical receiver; redatuming the geophysical data, wherein there is provided multi-component geophysical data, and redatuming geophysical data comprises obtaining at least one focussing function and/or at least one Green's function from the multi-component geophysical data, wherein the provided multicomponent geophysical data has not been and is not processed to account for free surface multiples; and using the redatumed geophysical data to determine a location, a direction, or a depth of a well; and producing hydrocarbons from the well based on the determined location, direction, or depth.

2. A method as claimed in claim 1, wherein the provided multi-component geophysical data has not been and is not processed to account for a source wavelet.

3. A method as claimed in claim 1, wherein using the redatumed geophysical data to prospect for hydrocarbons by selecting a location, a direction, or a depth of a well comprises using the obtained at least one focussing function and/or at least one Green's function to generate an image of a geophysical structure and using the image of the geophysical structure to select the location, direction, or depth of the well.

4. A method as claimed in claim 1, wherein the plurality of geophysical receivers are located at or near a first depth and the plurality of geophysical sources are located at or near a second depth.

5. A method as claimed in claim 4, wherein the plurality of geophysical receivers are located proximate a sea bed and the plurality of geophysical sources are located proximate a sea surface.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) More details of the theory behind this invention and some worked examples are provided below with reference to following Figures, which are provided by way of example only.

(2) FIG. 1 shows a schematic view of how the present method (FIG. 1(b)) has fewer requirements than the known Marchenko method (FIG. 1(a));

(3) FIG. 2 shows a schematic view of seismic acquisition explaining the theory behind the present method;

(4) FIG. 3 shows an example of a test model used to test the present method;

(5) FIG. 4 shows comparative results of the test run on the model of FIG. 3 showing the effectiveness of the present method;

(6) FIG. 5 shows more comparative results of the test run on the model of FIG. 3 showing the effectiveness of the present method for irregular and simultaneous sources;

(7) FIG. 6 shows more comparative results of the test run on the model of FIG. 3 showing the effectiveness of the present method with respect to source location repeatability.

DETAILED DESCRIPTION OF THE INVENTION

1. Introduction to Theory

(8) As is mentioned in the introduction of this specification, retrieving the correct focusing and/or Green's function from a virtual source in a subsurface using data measured at the surface is a key component for any imaging and inversion algorithm. Yet, producing full wavefield responses that contain internal multiples with correct amplitudes has long presented a challenge. Recent research on Marchenko-type equations revealed that such virtual responses (the “Green's functions”), as well as special wavefields that focus at any chosen point in the subsurface (the “focussing functions”), can be constructed by iterative substitution (Broggini et al., 2012; Wapenaar et al., 2014a) or direct inversion (van der Neut et al., 2015a) of Marchenko-type equations. Knowledge of these Green's functions and/or the focussing functions enables accurate imaging of deep subsurface targets with no artefacts from overburden-related internal multiples, and opens the way to new techniques for target-oriented imaging (van der Neut et al., 2013) and inversion (Vasconcelos et al, 2014; Broggini and Robertsson, 2014), microseismic source localization (Behura et al., 2013), and primary prediction (Meles et al., 2015).

(9) However, despite first encouraging results reported on field data (Ravasi et al., 2015; van der Neut et al., 2015b), the coupled Marchenko equations require some very strict acquisition requirements, as discussed in the introduction of this specification. For instance, the sources and receivers must be coincident (i.e. at the same depth and location); accurate knowledge/calculation of the source wavelet is required; deconvolution of the source wavelet from the recorded data prior to redatuming is required; and removal of direct wave, source- and receiver-side ghosts and surface related multiples is needed.

(10) While Singh et al. (2015) have shown that free-surface multiples can also be included in the redatuming scheme, yielding an additional integral in the Marchenko equations that requires knowledge of a ghost model, methods to handle the other two strict requirements have not been presented to date.

(11) The present method addresses all three of these strict requirements, and more. As is mentioned above, and as is discuss in more detail later, the inventor has discovered that, provided the availability of seismic data that can multi-component, such as dual-sensor (e.g., pressure and vertical particle velocity) data, these strict requirements can be overcome.

(12) As is discussed in detail below, the inventor's discovery originated from a realisation that the coupled Marchenko equations can be combined with a one-way version of the Rayleigh integral representation to provide a new redatuming scheme that handles band limited seismic data from any acquisition system that presents arbitrarily located sources (preferably above a line of regularly sampled receivers), while also including free-surface effects. Here, the “band limited” data is data that is shaped by a source wavelet. As has been discussed above, using the present method such a wavelet does not need to be removed (and therefore known) prior to redatuming.

(13) The Green's function(s) and the focussing function(s) can then be found by solving this newly-found relationship, for instance by an optimisation procedure such as a gradient-based optimisation procedure, using multi-component recorded seismic data.

(14) As a schematic view of the fewer acquisition requirements needed for the present method, FIG. 1 is provided. In FIG. 1, each star indicates a source and each triangle indicates a receiver. The arrows show exemplary seismic energy paths reflecting off boundaries of a geophysical structure from source to receiver.

(15) FIG. 1(a) shows the ideal scenario for acquisition geometry required by the existing technology. Thus, all sources and receivers are located at the same depth as each other and are regularly and evenly spaced from one another. In reality, this ideal scenario is not physically possible, so data gathered in the real scenario is required to be pre-processed to transfer the real data to the idealised scenario. However, there are inevitable inaccuracies and approximations involved in doing this pre-processing, and the pre-processing is costly. Further, the waveform shown adjacent one of the sources is present in the Figure to emphasise that in the existing technology knowledge about the source wavelet is required for redatuming.

(16) FIG. 1(b) shows a typical acquisition set-up that works for the current method. As can be seen, the requirements have been much relieved in comparison to FIG. 1(a). For instance, the sources no longer need to be regularly spaced or at the same depth as the receiver. Further, the question mark (?) adjacent the source is present in the Figure to emphasise that no knowledge about the source wavelet is required for redatuming. Further still, a free surface (e.g. the sea surface) above the sources is shown to emphasise that the present method can handle free surface multiples and does not require it to be removed prior to redatuming, unlike the existing technology.

2. Background to the Marchenko Equations Theory of the Prior Art

(17) The existing technology (the Marchenko equations) is derived from considering two different states (states A and B), each being idealised/simplified models of a geophysical structure, such as the Earth's crust or sedimentary layer. The model of state A is different to the model of state B, as is described below.

(18) FIG. 2(a) shows state A. State A is the ideal scenario that has been discussed at length above. The model of the geophysical structure here comprises an inhomogeneous lossless medium below a transparent boundary (i.e., the acquisition level Λ.sub.R) that has substantially identical properties to the geophysical structure below depth Λ.sub.R. Of course, these substantially identical properties are somewhat simplified in comparison to the real geophysical structure. A second transparent boundary is located at an arbitrary depth level Λ.sub.F. Below the second transparent boundary, the model of state A is also identical to the physical medium of the geophysical structure below the level Λ.sub.R. Above the acquisition level Λ.sub.R there is no reflection surface or indeed anywhere in the model. By considering a source (the star) at x.sub.R′ just above the level Λ.sub.R, the reflection response R from source x.sub.R′ to receiver (the triangle) x.sub.R is shown in FIG. 2(a). The reflection response R is the data that is recorded at receiver x.sub.R when source x.sub.R is fired in the described state.

(19) A virtual source is located in the model of state A at x.sub.F (which is a focussing location, and is located at the arbitrary depth level Λ.sub.F). The up-radiating and down-radiating Green's functions (g.sup.− and g.sup.+) from source x.sub.F to receiver x.sub.R are defined in state A. Since state A is a simplified model of the actual geophysical substructure, the Green's functions and the reflection response R in state A are as they would be in the simplified geophysical substructure. Making such simplified models of geophysical structures is common in the field and need not be explained in detail here.

(20) FIG. 2(b) shows state B. The model of the geophysical structure here is simpler than that of state A. This model shares the same medium properties of the physical medium (and hence state A) only between depths Λ.sub.R and Λ.sub.F, but it is homogeneous above Λ.sub.R and below Λ.sub.F. Also, like state A, there is no reflection surface above Λ.sub.R, or indeed anywhere in the model. Again, similar to the model of state A, the skilled person is well aware of how to construct such a simplified model of a geophysical structure.

(21) In state B, by considering a source x.sub.F just below the level Λ.sub.F, upgoing f.sup.− and downgoing f.sup.+ focussing functions can be defined. The Green's function representations that undergird the coupled Marchenko equations (Wapenaar et al., 2014a) relate the various quantities in states A and B and can be written in the frequency domain as:

(22) $\begin{array}{cc}{g}^{-}\left(x_R,x_F\right)=\underset{{\Lambda }_R}{\int }R\left(x_R,x_R^{\prime }\right)f^{+}\left(x_R^{\prime },x_F\right){\mathrm{dx}}_R^{\prime }-f^{-}\left(x_R,x_F\right)& \left(1\right)\\ -g^{+*}\left(x_R,x_F\right)=\underset{{\Lambda }_R}{\int }{R}^{*}\left(x_R,x_R^{\prime }\right)f^{-}\left(x_R^{\prime },x_F\right){\mathrm{dx}}_R^{\prime }-f^{+}\left(x_R,x_F\right)& \left(2\right)\end{array}$

(23) Here, * denotes complex conjugation (i.e., time-reversal in the time domain), and the integrals in equations 1 and 2 correspond to multidimensional convolution in time domain. The broadband reflection response of the medium R(x.sub.R,x.sub.R′) must be here interpreted as twice the pressure recording at x.sub.R from a vertical dipole source at x.sub.R′ (Wapenaar and Berkhout, 1989—Appendix B). This gives an expression for the reflection of response:

(24) $\begin{array}{cc}R\left(x_R,x_R^{\prime }\right)=\frac{2}{j\omega \rho \left(x_R^{\prime }\right)}{\partial }_{z^{\prime }}{g}^{-}\left(x_R,x_R^{\prime }\right)& \left(3\right)\end{array}$

(25) where ∂.sub.z represents the spatial derivative along the depth axis, ω is the angular frequency, ρ the density at x.sub.R′, and j the imaginary unit.

(26) To solve equations 1 and 2, it is convenient to introduce a discrete framework of matrix-vector multiplications (van der Neut et al., 2015a):

(27) $\begin{array}{cc}\left[\begin{array}{c}-g^{-}\\ g^{+*}\end{array}\right]=\left(I-\left[\begin{array}{cc}0& R\\ R^{*}& 0\end{array}\right]\right)\left[\begin{array}{c}{f}^{-}\\ f^{+}\end{array}\right]& \left(4\right)\end{array}$

(28) where g.sup.± and f.sup.± are vectors that represent the downgoing and upgoing Green's functions and focussing functions respectively and in which the seismic traces are concatenated in the time-space domain, while R and R* can be seen as operators that apply multidimensional convolution and correlation. Matrices I and 0 are the identity matrix and a matrix filled with zeros.

(29) It is g.sup.+, g.sup.−, f.sup.+ and f.sup.− that are desirable to be known for redatuming, imaging and inversion. However, in equation 4 there are effectively four unknowns (g.sup.+, g.sup.−, f.sup.+ and f.sup.−) in two equations, which is obviously not directly solvable.

(30) To avoid having to solve an underdetermined system of two equations and four unknowns (g.sup.+, g.sup.−, f.sup.+ and f.sup.−), a causality argument can be invoked: by noting that the Green's function contains a direct arrival followed by a scattering coda, a window Θ is designed such that all events after t.sub.d(x.sub.R,x.sub.F) (the time take for the Green's function direct wave to arrive at receiver x.sub.R from source x.sub.F) and before −t.sub.d(x.sub.R,x.sub.F), including the Green's function direct wave itself, are cancelled (i.e., Θg.sup.+=0 and Θg.sup.−=0). Fortunately, focussing functions respond to the action of the window in a completely different manner, i.e., Θf.sup.+=f.sub.d.sup.+ and Θf.sup.−=f.sup.− (Wapenaar et al., 2014a). Here, it is assumed that f.sup.+=f.sub.d.sup.++f.sub.m.sup.+, that is, f.sup.+ is composed of a direct wave f.sub.d.sup.+ followed by a subsequent coda f.sub.m.sup.+. The application of the window Θ to equation 4 leaves us with the coupled Marchenko equations:

(31) $\begin{array}{cc}\left[\begin{array}{c}\Theta {\mathrm{Rf}}_d^{+}\\ 0\end{array}\right]=\left(I-\left[\begin{array}{cc}0& \Theta R\\ \Theta R^{*}& 0\end{array}\right]\right)\left[\begin{array}{c}{f}^{-}\\ f_m^{+}\end{array}\right]& \left(5\right)\end{array}$

(32) which can be solved by direct inversion (van der Neut et al., 2015a) or expanded as a Neumann series by noting that equation 5 is a discretized Fredholm integral equation of the second kind (Slob et al., 2014).

3. Theory Behind the Present Method

(33) With regard to FIG. 2, the inventor constructed a new state, state C shown in FIG. 2(c). State C represents the real, physical acquisition scenario, and (unlike state A) has not been idealised. This state is the same as state A as far as it is a model whose properties are substantially identical (though preferably simplified in comparison to the real geophysical structure) to the actual properties of the geophysical structure, and it has the same two transparent depth levels Λ.sub.R and Λ.sub.F. However, unlike state A, there is also a free surface above the above the depth level Λ.sub.R at depth Λ.sub.S; and the source x.sub.S is located at a depth above the depth level Λ.sub.R (but below depth level Λ.sub.S). Receivers x.sub.R and x.sub.R′ are located at the same depth Λ.sub.R.

(34) With the aim of including the effect of having the source at different depth to the receivers and including the free surface heterogeneity in equations 1 and 2, the inventor used the one-way version of Rayleigh integral representation (Wapenaar and Berkout, 1989) to link the reflection response R in state A to a pair of upgoing and downgoing Green's functions in state C:

(35) $\begin{array}{cc}{G}^{-}\left(x_R^{\prime },x_S\right)=-\underset{{\Lambda }_R}{\int }\frac{2}{j\omega \rho \left(x_R\right)}{\partial }_z{G}^{+}\left(x_R,x_S\right)G^{-}\left(x_R,x_R^{\prime }\right){\mathrm{dx}}_R& \left(6\right)\end{array}$

(36) Here, G.sup.−(x′.sub.R,x.sub.S) is the Green's function response measured by receiver x.sub.R′ originating from source x.sub.S in state C; G.sup.+(x.sub.R,x.sub.S) is the Green's function response measured by receiver x.sub.R originating from source x.sub.S in state C; and G.sup.−(x.sub.R,x′.sub.R) is the Green's function response measured at receiver x.sub.R originating from a source at receiver x.sub.R′ in state A. Similarly to equation 3 above, here ∂.sub.z represents the spatial derivative along the depth axis, ω is the angular frequency, ρ the density at x.sub.R, and j the imaginary unit.

(37) Multiplying each side of equation 6 by s(ω)(jωρ(x′.sub.R)).sup.−1θ.sub.z′ where s(ω) is the source wavelet spectrum, the inventor obtained:

(38) $\begin{array}{cc}\underset{\underset{-v_z^{-}\left(x_R^{\prime },x_S\right)}{︸}}{\frac{s\left(\omega \right)}{j\mathrm{\omega \rho }\left(x_R^{\prime }\right)}{\partial }_{z^{\prime }}{G}^{-}\left(x_R^{\prime },x_S\right)}=\underset{{\Lambda }_R}{\int }\underset{\underset{v_z^{+}\left(x_R,x_S\right)}{︸}}{-\frac{s\left(\omega \right)}{j\mathrm{\omega \rho }\left(x_R\right)}{\partial }_z{G}^{+}\left(x_R,x_S\right)}\underset{\underset{R\left(x_R,x_R^{\prime }\right)}{︸}}{\left(\frac{2}{j\omega \rho \left(x_R^{\prime }\right)}{\partial }_{z^{\prime }}{G}^{-}\left(x_R,x_R^{\prime }\right)\right)}{\mathrm{dx}}_R& \left(7\right)\end{array}$

(39) As can be seen from the annotations below equation 7, by starting with equation 6 and manipulating it using the factor s(ω)(jωρ(x′.sub.R)).sup.−1θ.sub.z′ the inventor has found a relationship between upgoing particle velocity, downgoing particle velocity and the response function R of the Marchenko equations. This is because the expression on the left hand side of equation 7 is equal to the negative of the upgoing particle velocity at receiver x.sub.R′ from a wavefield originating from source x.sub.S (−v.sub.z.sup.−(x′.sub.R,x.sub.S)); the first expression on the right hand side of equation 7 is equal to the downgoing particle velocity at receiver x.sub.R from a wavefield originating from source x.sub.S(v.sub.z.sup.+(x.sub.R,x.sub.S)); and the second expression on the right hand side of equation 7 is equal to the reflection response R of the measured reflected seismic wavefield measured at receiver x.sub.R from a source x.sub.R′ in state A (this expression is equivalent to equation 3).

(40) Thus, equation 7 relates the reflection response R to the vertical particle velocity measurements from a source at x.sub.S, decomposed into their up- and downgoing constituents at the receiver level Λ.sub.R, for the actual physical scenario (state C).

(41) If each side of equation 1 is multiplied by v.sub.z.sup.+(x.sub.R,x.sub.S) and v.sub.z.sup.+*(x.sub.R,x.sub.S) is multiplied to equation 2, after performing integration over receivers and using equation 7, a new set of Green's function representations can be derived such that the focussing functions in state B and the Green's functions in state A are now linked via the upgoing and downgoing vertical particle velocity measurements v.sub.z.sup.− and v.sub.z.sup.+ measurements in state C. These equations are given below.

(42) $\begin{array}{cc}\underset{{\Lambda }_R}{\int }{v}_z^{+}\left(x_R,x_S\right)G^{-}\left(x_R,x_F\right){\mathrm{dx}}_R=-\underset{{\Lambda }_R}{\int }{v}_z^{-}\left(x_R^{\prime },x_S\right)f^{+}\left(x_R^{\prime },x_F\right){\mathrm{dx}}_R^{\prime }-\underset{{\Lambda }_R}{\int }{v}_z^{+}\left(x_R,x_S\right)f^{-}\left(x_R,x_F\right){\mathrm{dx}}_R& \left(8\right)\\ -\underset{{\Lambda }_R}{\int }{v}_z^{+}\left(x_R,x_S\right)G^{+*}\left(x_R,x_F\right){\mathrm{dx}}_R=-\underset{{\Lambda }_R}{\int }{v}_z^{-*}\left(x_R^{\prime },x_S\right)f^{-}\left(x_R^{\prime },x_F\right){\mathrm{dx}}_R^{\prime }-\underset{{\Lambda }_R}{\int }{v}_z^{+*}\left(x_R,x_S\right)f^{+}\left(x_R,x_F\right){\mathrm{dx}}_R& \left(9\right)\end{array}$

(43) In summary, by linking the response function R to vertical upgoing and downgoing particle velocities (i.e. up/down-separated multi-component seismic data) using the Rayleigh relationship of equation (6), the inventor has been able to eliminate the response function R from equations 1 and 2, which means that there is no longer any need to account for the source wavelet. Instead, in equations 8 and 9, the Green's functions and the focussing functions are related only via actual measured multi-component seismic data.

(44) Further, in contrast to equations 1 and 2, equations 8 and 9 no longer require the source(s) and the receivers to be at the same depth, since the integrations are not performed over the source location x.sub.S. This is unlike equations 1 and 2 where the integration is performed over the source location x.sub.R′. Similarly, the sources no longer need to be regularly spaced.

(45) In the existing Marchenko method, the focussing functions are described sometimes as a functions of x.sub.R (the focussing functions outside the integrals of equations 1 and 2) and sometimes as functions of x.sub.R′ (the focussing functions inside the integrals of equations 1 and 2). However, in equations 1 and 2, x.sub.R refers to a receiver and x.sub.R′ refers to a source. Because of this, in order to transition from equations 1 and 2 to equations 4 and 5 (the latter being the equation to be solved when performing the existing Marchenko method) it is required to have sources (x.sub.R′) and receivers (x.sub.R) at the same locations as each other. In contrast, the present method does not require sources and receivers to be at the same location. The focussing functions in equations 8 and 9 are only functions of receiver locations (x.sub.R′ and x.sub.R). There is therefore no constraint on the source sampling or source location.

(46) Equations 8 and 9 can be written as:

(47) 0$\begin{array}{cc}\left[\begin{array}{c}-{\stackrel{\_}{g}}^{-}\\ {\stackrel{\_}{g}}^{+*}\end{array}\right]=\left[\begin{array}{cc}{v}_z^{+}& v_z^{-}\\ v_z^{-*}& v_z^{+*}\end{array}\right]\left[\begin{array}{c}{f}^{-}\\ f^{+}\end{array}\right]& \left(10\right)\end{array}$

(48) Here, v.sub.z.sup.− and v.sub.z.sup.+(v.sub.z.sup.−* and v.sub.z.sup.+*) are operators that apply multidimensional convolution (correlation) with the upgoing and downgoing vertical velocity measurements. Moreover, by noting that the terms on the left hand side of equations 8 and 9 can be interpreted as the upgoing and downgoing Green's functions with free-surface effects included from the focal point x.sub.F, to sources x.sub.S, the inventor has written these terms as g.sup.+=v.sub.z.sup.+g.sup.+ and g.sup.−=v.sub.z.sup.+g.sup.−.

(49) To solve equation 10, a causality argument is invoked, in a similar way as for equation 5. However, the specific causality argument in this case is different.

(50) Since a step of multidimensional convolution (or correlation) with upgoing and downgoing vertical velocity measurements is applied to Green's and focussing functions, the causality argument holds only for the Green's functions, i.e. and not for the focussing functions.

(51) In particular, if a window Θ that mutes energy before −t.sub.d(x.sub.S,x.sub.F) and after t.sub.d(x.sub.S,x.sub.F) (where t.sub.d(x.sub.S,x.sub.F) is the time taken for the Green's function direct wave to arrive at source location x.sub.S from virtual source x.sub.F) is applied to equation 10, focussing functions can be obtained by solving the following system of equations

(52) $\begin{array}{cc}\left[\begin{array}{c}-\stackrel{\_}{\Theta }{v}_z^{-}{f}_d^{+}\\ -\stackrel{\_}{\Theta }{v}_z^{+*}{f}_d^{+}\end{array}\right]=\left[\begin{array}{cc}\stackrel{\_}{\Theta }{v}_z^{+}& \stackrel{\_}{\Theta }{v}_z^{-}\\ \stackrel{\_}{\Theta }{v}_z^{-*}& \stackrel{\_}{\Theta }{v}_z^{+*}\end{array}\right]\left[\begin{array}{c}{f}^{-}\\ f_m^{+}\end{array}\right]& \left(11\right)\end{array}$

(53) where the upgoing focussing function f.sup.+=f.sub.d.sup.++f.sub.m.sup.+. That is, r is composed of a direct wave f.sub.d.sup.+ followed by a subsequent coda f.sub.m.sup.+.

(54) Note that equation 11 is a Fredholm integral of the first kind which cannot be expanded into a Neumann series. In other words, this equation can be solved by direct inversion. It may not be possible to solve equation 11 for the focussing functions using iterative substitution, unlike equation 5.

(55) Thus, as should be appreciated from the above, while equation 5 is only valid when sources and receivers are co-located, equation 11 holds for any (or even a single) source located above Λ.sub.R: when multiple sources are available, solving the ensemble of equations becomes a better posed problem. Of course, state C and equations 6-11 could be reformulated so that the location of the receiver is arbitrary, or so that the location of the source is at the same depth or a lower depth than the receiver.

4. Numerical Example of Current Method, and Comparison with Existing Marchenko Method

(56) In order to demonstrate that the present method of using multi-component seismic data works for a more relaxed acquisition set-up, the following comparison experiment was conducted.

(57) In order to perform this comparison, a model was used as shown in FIG. 3. FIG. 3 shows a constant velocity, variable density subsurface model. A solid black line 1 indicates the source array, and a solid grey line 2 refers to the receiver array. A cross is used to indicate the subsurface focal point.

(58) The model is 2D representation of a simplified geophysical structure comprising different layers of different densities. It includes an Earth's surface-sea boundary at a depth of 200 m. The density of the different layers in the substructure can be seen as different shades, as can be seen from the density key in FIG. 3. Thus, as can be seen several layers are present with differing shapes, thicknesses and densities. The model is constructed so that velocity is constant throughout (v=2400 m/s). A free surface (i.e. a sea-air boundary) is included in the model at 0 m depth.

(59) Dual-sensor (i.e., pressure and vertical particle velocity) data are modelled in the model. The seismic energy is input into the model using a line 1 of 201 (two hundred and one) sources located at 10 m depth. The seismic energy is recorded in the model using a line 2 of 201 (two hundred and one) receivers at Λ.sub.R=50 m depth. A focal point (focussing location) is selected at x.sub.F=1200 m, 1040 m (at the cross on FIG. 3). The modelling was done using a finite-difference vector-acoustic scheme (e.g., Vasconcelos et al., 2014), with a 20 Hz Ricker wavelet.

(60) The dual-sensor data is recorded, and is then processed so that it is decomposed into up/down-separated multi-component data (i.e. the upgoing particle velocity and downgoing particle velocity components).

(61) The upgoing v.sub.z.sup.− and downgoing v.sub.z.sup.+ constituents are found using v.sub.z.sup.±=v.sub.z∓k.sub.z/(ρω)ρ, where k.sub.z is the vertical wavenumber (Wapenaar, 1998), ρ is the density of the geophysical structure, ω is the frequency of the seismic energy recorded and p is the pressure of the seismic energy recorded.

(62) The direct part of the focussing function (f.sub.d.sup.+) is constructed by reversing in time the first arriving wave from the focussing location x.sub.F to the receiver depth level Λ.sub.R, computed via finite-difference in a constant velocity (v.sub.0=2400 m/s), constant density (ρ.sub.0=2400 kg/m.sup.3) subsurface model. Note this is the only information that needs to be known with respect to the model.

(63) To construct upgoing and downgoing focussing functions and subsequently the scattering coda of the Green's function from the measured data, the measured multi-component seismic data recorded by the receivers at Λ.sub.R is used in equation 11. Equation 11 is inverted using a gradient-based optimization scheme, such as LSQR (Pauge and Saunders, 1982). The outputs of this inversion are the focussing functions and Green's functions calculated by the present method.

(64) Using the model of FIG. 3, focussing functions for the existing Marchenko redatuming method can also be calculated. This is done by inverting the original Marchenko equation 5 using the same number of iterations as for the present method's inversion of equation 11. In both cases a preconditioning operator has been added as suggested by van der Neut et al. (2015), ensuring that the retrieved focussing functions contain events only inside the window Θ.

(65) It should be noted that, in order to meet the strict acquisition requirements necessary to perform the existing Marchenko method, in order to calculate the focussing functions for the existing Marchenko method, a new seismic dataset had to be produced using the model of FIG. 3. In this new model, the geophysical structure is the same as for when the data for the present method was produced. However, the sources were placed at the same location as the receivers; the free-surface was replaced by an absorbing boundary; the source wavelet was known and deconvolved from the recorded data; and the direct wave has been removed.

(66) Since it is known that the existing Marchenko outputs correct focussing functions as long as the strict acquisition requirements are met, if the focussing functions output from the present method match the focussing functions of the existing Marchenko method then it is evidence that the present method can successfully construct focussing functions without the need for the strict acquisition requirements. [As has been mentioned above, the focussing functions output from the present method should be, from the theory at least, the very same focussing functions output from the existing Marchenko method for the same focussing location in the same geophysical structure.]

(67) With regard to FIG. 4(a) this shows a comparison of the downgoing focussing functions calculated after 30 iterations using the existing Marchenko method (solid grey lines, marked “Abs”) and the new method (dotted lines, marked “Fs”). The plot shows receiver position along the horizontal axis against time along the vertical axis. For clarity, only a certain number of traces are shown (e.g. only every 10.sup.th receiver trace is shown); this is only done for the sake of clarity—if all receivers were shown the diagram would be very unclear.

(68) As can be seen from the close correlation between the solid and dotted lines, the downgoing focussing function calculated using the present method (using a less-strict acquisition set-up) is very similar—almost identical—to the downgoing focussing function calculated using the existing Marchenko method (using a very-strict acquisition set-up). This is evidence that the present method can be used to obtain focussing functions with fewer acquisition requirements or without the need for pre-processing acquired data into the idealised scenario.

(69) With regard to FIG. 4(b) this shows a comparison of the upgoing focussing functions calculated after 30 iterations using the existing Marchenko method (solid grey lines, marked “Abs”) and the new method (dotted lines, marked “Fs”). The plot shows receiver position along the horizontal axis against time along the vertical axis. For clarity, only a certain number of traces are shown (e.g. only every 10.sup.th receiver trace is shown); this is only done for the sake of clarity—if all receivers were shown the diagram would be very unclear.

(70) As can be seen from the close correlation between the solid and dotted lines, the upgoing focussing function calculated using the present method (using a less-strict acquisition set-up) is very similar—indeed almost identical—to the upgoing focussing function calculated using the existing Marchenko method (using a very-strict acquisition set-up). This is evidence that the present method can be used to obtain focussing functions with fewer acquisition requirements or without the need for pre-processing acquired data into the idealised scenario.

(71) The Green's function itself for the focussing location x.sub.F was also calculated using the present method and the model of FIG. 3. This was done similarly to how the focussing functions were found: by using the recorded multi-component seismic data in equation 11.

(72) In order to show that the Green's function calculated by the present method are accurate, despite the flexible acquisition geometry used to gather the seismic data as well as minimal processing, a Green's function for the source x.sub.F was also created in a separate data run in the model. If the Green's functions output from the present method match the actual Green's functions of the geophysical structure when a source is located at the focussing location then it is evidence that the present method can successfully construct Green's functions without the need for accurate velocity models and strict acquisition/processing requirements prior to redatuming. [As has been mentioned above, the Green's functions output from the present method should be, from the theory at least, the very same Green's functions as would be expected to be measured if a source were actually located at the focussing location (i.e. the “actual” Green's function). This cannot be said for the existing Marchenko method where the Green's functions calculated are not the same as the “actual” Green's functions; instead they are the Green's functions in the ideal scenario of state A.]

(73) With regard to FIG. 4(c) this shows a comparison of the Green's functions using the source at x.sub.F (solid grey lines, marked “True”) and the present method (dotted lines, marked “Estimate”). The plot shows receiver position along the horizontal axis against time along the vertical axis. For clarity, only a certain number of traces are shown (e.g. only every 10.sup.th receiver trace is shown); this is only done for the sake of clarity—if all receivers were shown the diagram would be very unclear.

(74) As can be seen from the close correlation between the solid and dotted lines, the Green's functions calculated using the present method (using a less-strict acquisition set-up) are very similar—indeed almost identical—to the “actual” Green's functions. This is evidence that the present method can be used to obtain Green's functions with fewer acquisition requirements than for the Marchenko method or without the need for pre-processing acquired data into the idealised scenario.

(75) As an aside about this numerical test, it should be noted that because both the data and the initial focussing functions are band-limited by a 20 Hz Ricker wavelet, the retrieved Green's function (and thus also the “actual” modelled Green's function) is shaped by the autocorrelation of a 20 Hz Ricker wavelet. While the wavelet that acts on the data has been considered unknown, any other wavelet could be chosen for the initial focussing function. It could be then deconvolved deterministically from the retrieved Green's function to broaden its bandwidth. Lastly, given the high computational cost in the application of the multidimensional convolution (correlation) operators, it is important that the gradient-based inversion converges quickly to the correct solution. It has been found that the results of the iteration converge quickly to the correct values when finding the focussing and/or Green's functions.

(76) As a further extension of the present method, to show that focussing functions and Green's functions can be found using data gathered from an even more flexible acquisition set-up and without requiring pre-processing to an ideal scenario, two further comparison were carried out.

(77) Firstly, with regard to FIGS. 5(a) and (b), the new method was performed using data gather from the model of FIG. 3, except where the previously evenly-distributed source array was randomly decimated, in this case by a factor of two. This means that there are fewer (half of the previously available) equations available to solve for the focussing functions. Further, this means that the sources are not regularly distributed (as well as not being coincident with, and at the same depth as, the receivers).

(78) With regard to FIG. 5(a) this shows a comparison of the focussing functions calculated after 30 iterations using the new method with no sources removed (solid grey lines, marked “full”) and the new method with half the sources randomly decimated (dotted lines, marked “Missing”). The plot shows receiver position along the horizontal axis against time along the vertical axis. For clarity, only a certain number of traces are shown (e.g. only every 10.sup.th receiver trace is shown); this is only done for the sake of clarity—if all receivers were shown the diagram would be very unclear.

(79) As can be seen from the close correlation between the solid and dotted lines, the focussing functions calculated using the present method having half the sources randomly decimated is very similar—almost identical—to the focussing functions calculated using the present method but using all of the sources. This is evidence that the present method can be used to obtain accurate focussing functions with having fewer sources and by having non-regularly spaced sources and without the need for pre-processing acquired data into the idealised scenario.

(80) With regard to FIG. 5(b) this shows a comparison of the Green's functions calculated from the focusing functions using the source at x.sub.F (solid grey lines, marked “True”) and the present method but with half the sources randomly decimated (dotted lines, marked “Estimate”). The plot shows receiver position along the horizontal axis against time along the vertical axis. For clarity, only a certain number of traces are shown (e.g. only every 10.sup.th receiver trace is shown); this is only done for the sake of clarity—if all receivers were shown the diagram would be very unclear.

(81) As can be seen from the close correlation between the solid and dotted lines, the Green's functions calculated using the present method having half the sources randomly decimated is very similar—almost identical—to the Green's function calculated using the present method but using all of the sources. This is evidence that the present method can be used to obtain accurate Green's functions with having fewer sources and by having non-regularly spaced sources and without the need for pre-processing acquired data into the idealised scenario.

(82) Secondly, with regard to FIGS. 5(c) and (d), the new method was performed using data gather from the model of FIG. 3, except where the previously sequentially-fired sources were fired such that two adjacent sources are shot simultaneously (with a random ignition time of t.sub.i˜U(0,0.25)s.

(83) With regard to FIG. 5(c) this shows a comparison of the focussing functions calculated after 30 iterations using the new method with purely sequential sources (solid grey lines, marked “full”) and the new method with simultaneous sources (dotted lines, marked “Sim”). The plot shows receiver position along the horizontal axis against time along the vertical axis. For clarity, only a certain number of traces are shown (e.g. only every 10.sup.th receiver trace is shown); this is only done for the sake of clarity—if all receivers were shown the diagram would be very unclear.

(84) As can be seen from the close correlation between the solid and dotted lines, the focussing functions calculated using the present method using simultaneous sources is very similar—almost identical—to the focussing functions calculated using the present method but using purely sequential sources. This is evidence that the present method can be used to obtain accurate focussing functions whilst using simultaneous sources, which can cut down on acquisition time.

(85) In this case, solving for the focusing functions is effectively equivalent to solving a linear combination equation 11 that has been solved for sequential sources. In other words, direct inversion of the modified Marchenko equations is also responsible for source deblending.

(86) With regard to FIG. 5(d) this shows a comparison of the Green's functions using the source at x.sub.F (solid grey lines, marked “True”) and the present method but with simultaneous sources (dotted lines, marked “Estimate”). The plot shows receiver position along the horizontal axis against time along the vertical axis. For clarity, only a certain number of traces are shown (e.g. only every 10.sup.th receiver trace is shown); this is only done for the sake of clarity—if all receivers were shown the diagram would be very unclear.

(87) As can be seen from the close correlation between the solid and dotted lines, the Green's function calculated using the present method with simultaneous sources is very similar—almost identical—to the Green's function calculated using a source at x.sub.F. This is evidence that the present method can be used to obtain accurate Green's functions whilst using simultaneous sources and without the need for pre-processing acquired data into the idealised scenario.

(88) As a further extension of the present method, to show that focussing functions (and hence Green's functions) can be found using data gathered from an even more flexible acquisition set-up and without requiring pre-processing to an ideal scenario, a further comparison was carried out.

(89) FIG. 6 demonstrates the repeatability of the present method from one survey to another subsequent survey, even if the two surveys do not have the same source locations (as is typically the case when taking subsequent seismic surveys, such as in 4D seismic monitoring).

(90) With regard to FIG. 6(a), this shows a schematic illustration of how sources of a first survey (the black stars) are in different locations to the sources of a second survey (the white stars).

(91) In the modelling, firstly the sources at the first locations (black stars) were fired, the data was acquired and the focussing functions were calculated. The calculated focussing functions are shown in the solid grey lines of FIG. 6(b).

(92) Secondly, the sources were moved to the second locations (white stars), which are offset by 15 m from the first locations (black stars) and hence are not at the same location as the first locations. The sources were then fired, the data was acquired and the focussing functions were calculated. The calculated focussing functions are shown in the dotted lines of FIG. 6(b).

(93) As can be seen from FIG. 6(b), the focussing functions from the first acquisition are largely identical to the focussing functions from the second acquisition, despite the sources being in fired in different locations. This is evidence that the present method can be used to obtain accurate focussing functions and Green's functions from one survey to the next, whilst having the sources at different locations from one survey to the next, and without the need for pre-processing acquired data into the idealised scenario.

5. Summary

(94) The inventor has found a very flexible redatuming scheme, based on the coupled Marchenko equations, that uses multi-component seismic data (e.g. from dual-sensor data) without the need for pre-processing to an idealised scenario prior to redatuming. The redatuming is robust and can handle any acquisition scenario with sources located at a different depth to the receiver (or indeed the same depth, depth is no longer an issue). Further, the sources or the receivers are no longer required to be located at the same horizontal position, and sources do not need to be regularly spaced. Further, the source wavelet need not be known nor accounted for, and free-surface effects need not be separately accounted for. The application of the present method may have benefits in redatuming and imaging of ocean-bottom and borehole seismic data from sequential, or even simultaneous, sources.

(95) It should be apparent that the foregoing relates only to the preferred embodiments of the present application and the resultant patent. Numerous changes and modification may be made herein by one of ordinary skill in the art without departing from the general spirit and scope of the invention as defined by the following claims and the equivalents thereof.

6. REFERENCES

(96) 1. Behura J., and Snieder, R., 2013. Imaging direct as well as scattered events in microseismic data using inverse scattering theory. SEG Technical Program, Expanded Abstracts. 2. Berryhill, J. R., 1984. Wave-equation datuming before stack. Geophysics, 49(11), 2064-2066. 3. Broggini, F., Snieder, R., and Wapenaar, K., 2012. Focusing the wavefield inside an unknown 1D medium: Beyond seismic interferometry. Geophysics, 77, A25-A28. 4. Broggini, F., and Robertsson, J. O. A., 2014. FD injection utilizing the wavefields generated by Marchenko redatuming: A target-oriented approach. SEG Technical Program Expanded Abstracts. 5. Meles, G. A., Wapenaar, K., Curtis A., 2015, Synthesising Primary reflections by Marchenko Redatuming and convolutional interferometry, EIP Report, 21. 6. van der Neut, J., C. A. Vidal, N. Grobbe, and K. Wapenaar, 2013, Turning onesided illumination into two-sided illumination by target-enclosing interferometric redatuming. 75th EAGE Conference & Exhibition, Extended Abstracts. 7. van der Neut, J., Thorbecke, J., Wapenaar, K., and Slob, E., 2015a. Inversion of the multidimensional Marchenko equation. 77th annual EAGE meeting. 8. van der Neut, J., Wapenaar, K., Thorbecke, J., and Slob, E., 2015b. Practical challenges in adaptive Marchenko imaging, 85th annual SEG meeting. 9. Pauge, C. C. and Saunders, M. A., 1982. LSQR: An algorithm for sparse linear equations and sparse least squares. Transactions on Mathematical Software, 8(1), 43-71. 10. Ravasi, M., Vasconcelos, I., Kritski, A., Curtis, A., da Costa Filho, C. & Meles, G. A., 2015. Marchenko imaging of Volve Field, North Sea. 77th EAGE Conference and Exhibition, Extended Abstracts. 11. Singh, S., Snieder, R., Behura, J., van der Neut, J., Wapenaar, K., and Slob, E., 2015. Marchenko imaging: Imaging with primaries, internal multiples, and free-surface multiples. Geophysics, 80 (5), S165-S174. 12. Vasconcelos, I., M. Ravasi, and J. van der Neut, 2014, An interferometry-based, subsurface-domain objective function for waveform inversion: 76th EAGE Conference & Exhibition, Extended Abstracts. 13. Wapenaar, C. P. A., and Berkhout, A. J., 1989, Elastic wave field extrapolation: Redatuming of single- and multi-component seismic data: Elsevier Science Publ. Co., Inc. 14. Wapenaar, K., 1998, Reciprocity properties of one-way propagators: Geophysics, 63, 1795-1798 15. Wapenaar, K., Thorbecke, J., van der Neut, J., Broggini, F., Slob, E., and Snieder, R., 2014a, Marchenko imaging: Geophysics, Vol. 79 (3), WA39-WA57. 16. Wapenaar, K., Slob, E, 2014b, On the Marchenko equation for multicomponent single-sided reflection data, Geophysics Journal International, 199 pp. 1367-1371. 17. Van der Neut, J., Wapenaar, K., Thorbecke, J., Vasconcelos, I., 2014, Internal multiple suppression by adaptive Marchenko redatuming, SEG Denver 2014 Annual Meeting, pp. 4055-4059. 18. Liu, Y., Van der Neut, J., Arnsten, B., Wapenaar, K., 2015, Data-driven deep local imaging using both surface and borehole seismic data, SEG New Orleans Annual Meeting 2015, pp 5554-5559.