GB2553890A - Improvement to seismic processing based on predictive deconvolution - Google Patents

Abstract

A seismic processing method 200 for determining a subsurface model comprises receiving 210 seismic data recorded in a space-time domain. A gapped deconvolution prediction operator- is derived 220 from the received seismic data in a model domain which is different from said space-time domain. The deconvolution prediction operator is converted 222 to the space-time domain using a linear conversion operator. The subsurface model is determined 223 by convolving said gapped deconvolution prediction operator with the received seismic data. The linear conversion operator may be a reverse linear Radon operator, reverse parabolic Radon operator, reverse hyperbolic Radon operator, reverse migration operator, reverse curvelet transform, reverse AVP operator, reverse tilted hyperbolic Radon operator. The subsurface model may comprise primaries, multiples or both.

Description

(71) Applicant(s):

CGG Services

Avenue Carnot, Massy 91300,

France (including Overseas Departments and Territori es) (72) Inventor(s):

Gordon Poole (51) INT CL:

G01V1/36 (2006.01) G01V 1/38 (2006.01) (56) Documents Cited:

WO 2015/118409 A2

SEG Technical Programme Expanded Abstracts, 2007, Broadhead et al, Predictive deconvolution by frequency domain Wiener filtering (58) Field of Search:

INT CL G01V

Other: WPI, EPODOC, Patents Fulltext, Geophysics (74) Agent and/or Address for Service:

Mewburn Ellis LLP

City Tower, 40 Basinghall Street, LONDON, Greater London, EC2V 5DE, United Kingdom

Title of the Invention: Improvement to seismic processing based on predictive deconvolution Abstract Title: Convolving seismic data with a gapped deconvolution prediction operator derived domain other than the space-time domain

A seismic processing method 200 for determining a subsurface model comprises receiving 210 seismic data recorded in a space-time domain. A gapped deconvolution prediction operator'is derived 220 from the received seismic data in a model domain which is different from said space-time domain. The deconvolution prediction operator is converted 222 to the space-time domain using a linear conversion operator. The subsurface model is determined 223 by convolving said gapped deconvolution prediction operator with the received seismic data. The linear conversion operator may be a reverse linear Radon operator, reverse parabolic Radon operator, reverse hyperbolic Radon operator, reverse migration operator, reverse curvelet transform, reverse AVP operator, reverse tilted hyperbolic Radon operator. The subsurface model may comprise primaries, multiples or both.

200

FIG«8

First Gs-rfer Solace Side Makipfe >

PrisnssV”

T

Seco-^ad

Order

Receiver

Side Mu&ipfe n

i

I ,4 / >

f

4-.

V ii i ; !· 1 ί / u >

FIG.2a

FIG.2c

Ο

F F

X (Offset)

t (time) t (time)

x (Offset)

t (time)

200 fN fN

τ~Ί τ~Ί

IMPROVEMENT TO SEISMIC PROCESSING BASED ON PREDICTIVE

DECONVOLUTION fiBOSS^EFERENCEJO^ELAISD^APP|JCATIONS [0001] This application claims priority and benefit from U.S. Provisional Patent Application No. 62/358,200, filed July 5, 2016, for Improved predictive deconvolution, the entire content of which is incorporated in its entirety herein by reference.

BACKGROUND izchnxcalfield [0002] Embodiments of the subject matter disclosed herein generally relate to processing seismic data using a subsurface model. More specific, primaries and multiple reflections are modelled and separated as part of a process for generating an image of the subsurface.

DISCUSSION OF THE BACKGROUND [0003] Hydrocarbon exploration and development uses waves (e.g., seismic waves or electromagnetic waves) to explore the structure of underground formations on land and at sea (i.e., formations under the seafloor). As schematically illustrated in Figure 1, waves emitted by a source 1 at a known location penetrate an explored formation 2 and are reflected at interfaces 20, 21, 22 that separate the formation's layers having different layer impedances. Sensors 3 detect the reflected waves. The detected waves include primary reflections such as wave 4 which travels directly from formation interface 22 to sensor 3, and multiple reflections such as wave 5, which are reflected at least one additional time inside the formation before being detected. Note that, as used herein, the term formation refers to any geophysical structure into which source energy is used to perform seismic surveying, e.g., land or water based, such that a formation will include a water layer when the context is marine seismic surveying.

[0004] There are various types of multiples, e.g., surface-related multiples and interbed multiples. In the case of surface-related multiples, when the energy reflected from the subsurface reaches the water surface it will be reflected back downwards into the wafer column and subsurface. This produces a second set of reflected energy containing spurious events. Interbed multiples are similar, but in this case the downward reflecting surface is a rock interface in the subsurface.

[0005] Moreover, multiples can be characterized as belonging to different orders, e.g., first order, second order, third order, etc., based on the number of additional reflections involved, For example, a primary has a single reflection between a source S and a receiver R as shown in Figure 2a. By way of contrast a first order multiple (shown in Figure 2b) can have two additional reflections relative to the primary, whereas a second order multiple (shown in Figure 2c) can have four additional reflections.

[0006] In order to understand the structure of the explored underground formation, it is preferable that primaries and multiples be separated as part of the processing of the recorded seismic data, and frequently it is preferable to remove the multiples. Accordingly, there are also numerous processing techniques which have been developed to attenuate or suppress multiples in the recorded seismic data.

[0007] Many surface related multiple attenuation algorithms model multiples by convolving recorded data with a primary estimate. Probably the best known of these approaches is Surface-Related Multiple Estimation or SRME, where the primary estimate is initialized by the data itself (see for instance Estimation of multiple scattering by iterative inversion, Part 1: Theoretical considerations, Berkhout and Verschuur, 1997, Geophysics, 62, 158661595). [0008] The multiple model from SRME is typically adaptively subtracted from the data, a process which corrects for inaccuracies in the multiple model relating to the source wavelet, cross-talk between multiples and missing or under sampled data. When the multiple generator has not been sufficiently recorded SRME can break down.

[0009] In such cases, the primary estimate may be provided by a reflectivity image modelling as disclosed for instance by Pica, A., Poulain, G., David, B.,

Magesan, M., Baldock, S., Weisser, T., Hugonnet, P., and Herrmann, P (2005, 3D Surface-Related Multiple Modeling, Principles and Results, 75^th annual international meeting, SEG, Expanded abstracts, p. 2080-2083) or Green's function multiple modelling (see Wang, P. Jin, H., Xu, S. And Zhang, Y., 2011, Model-based water-layer demultiple, 81^st annual international meeting, SEG, expanded abstracts, p, 3551-3555).

[0010] Gapped deconvolution approaches such as tau-p deconvolution or 2D deconvolution offer an alternative where information about the multiple generator is derived from the periodicity of multiples in the data itself (see for instance Biersteker J., 2001, MAGIC: Shell's surface multiple attenuation technique, 70^th annual international meeting, SEG, expanded abstracts 13011304).

[0011] Backus (1959, Water reverberations - their nature and elimination, Geophysics, Vol. XXIV, N°2, pages 233-261) describes a second order correction for predictive deconvolution which improves the amplitude consistency of the multiple prediction. The approach has been adapted for the cases of SRME (see Hugonnet, P, 2002, Partial Surface Related Multiple Elimination, 72^nd annual international meeting, SEG, expanded abstracts 2102-2105) and Green's function multiple modelling (See Cooper, 3., Poole, G., Wombell, R. and Wang, P., 2015, Recursive Model-based Water-layer Demultiple, 77^th EAGE conference).

[0012] The following focuses on the shallow water demultiple based on predictive deconvolution:

[0013] 2D deconvolution employs correlations to predict the multiple generators (e.g. the water bottom or any intermediate layer on which multiple reflections may happen). Demultiple process normally consist of two steps, a first step of designing the deconvolution prediction operator (i.e. finding the prediction operator which, when convolved with the recorded data, will give only the multiples), and a second step of applying the deconvolution prediction operator in order to estimate the multiples. The estimated multiples can then be subtracted from recorded data to deliver the primaries.

[0014] A 2D deconvolution equation may be given by :

d - (S + R)g - Cg (Equation 1) where:

d is a vector representative of the estimated input data (excepted the primary);

g is a vector representative of the spatial deconvolution prediction operator;

S is a matrix corresponding to the spatial convolution with the data on the source side;

R is a matrix corresponding to the spatial convolution with the data on the receiver side, and

C is a spatial convolution matrix resulting in the addition of S and R.

[0015] The spatial convolutions hence consist in a number of convolutions between data and prediction operator within a defined spatial aperture followed by a summation. The summation of multiple convolutions may be performed within a Multiple Contribution Gather (MCG), i.e. a gather of convolved trace pairs that, when summed, predicts surface multiples for a target trace (as defined by Drafgoset et al. in « A perspective on 3D surfacerelated multiple elimination » GEOPHYSICS, VOL. 75, NO. 5 SEPTEMBEROCTOBER 2010; P. 75A245-75A261, 10.1190/1.3475413. The process may include frequency dependent operations to compensate for phase changes and/or Fresnel zone compensation. The Fresnel zone compensation may relate to low frequencies stacking in more strongly than higher frequencies. The compensation may be based on one or more spatial dimension on source and/or receiver sides. For example this may be based on 2πΙί, where f is frequency in Hertz, or in the case of 2D, the square-root thereof, In addition, an obliquity correction may optionally be applied based on the cosine of the angle of reflection.

[0016] This linear operation may also be described by the following pseudocode:

- Loop through all traces;

- Define a spatial aperture relating to a MCG;

- Loop through all traces in the MCG;

® Convolve the trace with the associated trace in the prediction operator;

® Accumulate the convolved data; and - Correct the accumulated data for phase and/or amplitude errors, [0017] The deconvolution prediction operator must have an appropriate gap followed by an active part with an appropriate operator length to prevent predictions of the primaries. The gap may be dependent on the data (which in the shallow section may approximate the kinematics of the deconvolution operator), water depth or water bottom time.

[0018] The least squares solution may be found, for example, by solving the following equation:

(^Td c'^!Cg (Equation 2) [0019] Equation 2 may be solved using conjugate gradients, Cholesky factorization, Lower-Upper (LU) decomposition, singular value decomposition (SVD), or other methods.

[0020] While 2D deconvolution is better suited in shallow water environments for the purpose of multiple attenuation, some errors are encountered in specific situations regarding the amplitudes of some predicted multiples, as will be apparent in the following description. Accordingly, it remains desirable to improve the multiple prediction using spatial deconvolution prediction operators.

SUMMARY [0021] According to one embodiment, there is a seismic processing method for determining a subsurface model comprising:

receiving seismic data recorded in a space-time domain;

defining a gapped deconvolution prediction operator from the received seismic data; and determining the subsurface model by convolving said gapped deconvolution prediction operator with the received seismic data, wherein said gapped deconvolution prediction operator is defined by deriving said gapped deconvolution prediction operator in a model domain which is different from said space-time domain, and by converting the deconvolution prediction operator derived in the model domain to the space-time domain using a linear conversion operator.

[0022] According to another embodiment, there is a data processing apparatus for determining a subsurface model, the apparatus comprising:

an interface configured to receive seismic data recorded in a spacetime domain; and a data processing unit configured to define a gapped deconvolution prediction operator from the received seismic data; and determine the subsurface model by convolving said gapped deconvolution prediction operator with the received seismic data, wherein said gapped deconvolution prediction operator is defined by deriving said gapped deconvolution prediction operator in a model domain which is different from said space-time domain, and by converting the deconvolution prediction operator derived in the model domain to the space-time domain using a linear conversion operator.

BRIEF DESCRIPTION OF THE DRAWINGS [0023] The accompanying drawings, which are incorporated in and constitute a part of the specification, illustrate one or more embodiments and, together with the description, explain these embodiments. In the drawings:

Figure 1 is a schematic illustration of exploring underground formations using waves including a primary and a multiple;

Figures 2a-2c illustrate a primary, a first order source side peg-leg multiple, and a second order receiver side peg-leg multiple, respectively;

Figure 3 shows actual amplitude terms of multiples taken directly from seismic data;

Figures 4a and 4b show predicted amplitude terms of multiples for the same seismic data as in Figure 3 which were predicted using a conventional prediction deconvolution approach, respectively at the source and at the receiver side;

Figure 4c shows the resulting summation of 4a and 4b;

Figure 5 illustrates an example unconstrained spatial prediction operator determined according to prior art;

Figure 6 illustrates a first exemplary spatial prediction operator determined according to the invention;

Figure 7 illustrates a second exemplary spatial prediction operator determined according to the invention;

Figure 8 is a flow chart of a seismic processing method for determining a subsurface model according to an embodiment;

Figure 9 is a schematic diagram of a computing device configured to implement one or more of the methods discussed herein.

DETAILED DESCRIPTION [0024] The following description of the exemplary embodiments refers to the accompanying drawings. The same reference numbers in different drawings identify the same or similar elements. The following detailed description does not limit the invention. Instead, the scope of the invention is defined by the appended claims. The following embodiments are discussed using terminology of seismic data processing. However, the described methods may be used for other wave data processing (e.g., electro-magnetic wave data).

[0025] Reference throughout the specification to one embodiment or an embodiment means that a particular feature, structure or characteristic described in connection with an embodiment is included in at least one embodiment of the subject matter disclosed. Thus, the appearance of phrases in one embodiment or in an embodiment in various places throughout the specification is not necessarily referring to the same embodiment, Further, the particular features, structures or characteristics may be combined in any suitable manner in one or more embodiments.

[0026] One problem associated with using the predictive deconvolution technique to attenuate multiples in received seismic data is that it tends to overestimate the amplitude of peg-leg multiples of order two or higher, as is now explained.

[0027] The combined source and receiver 2D gapped deconvolution equations may be given by:

^dik = Pik *^djk ⁺yL^dij*Sjk (Equation 3) j j where input data, d, is predicted by a multi-channel convolution between the data and surface consistent deconvolution prediction operator, g, with indices i, j, and k relating to spatial coordinates for the spatial summation plus primary, p. Index i relates to the shot position, index k relates to the receiver position, and index j relates to the spatial position the convolution between data and the operator takes place.

[0028] The surface consistent prediction operator may be shared for source and receiver sides with an appropriate gap and active operator length to prevent prediction of the primaries, p, and/or shaping of the wavelet. The gap may be constant, or may vary as a function of space and/or time. The linear equations may be solved using least squares inversion to find the prediction operator g, following which a multiple estimate is made by convolving the prediction operator with the recorded data.

[0029] To explore the accuracy of the multiple prediction, one may consider a subset of arrivals in the input data consisting of primaries followed by peg-leg multiples, d^^li, as illustrated in Figure 3. This dataset may be described in terms of a primary, p, source only side multiples, CRP , receiver only side multiples, P/^R' , and mixed side multiples, Operators S and R relate to convolution operations (provided by multiplications in the frequency domain) by the multiple generator on source and receiver sides respectively, a and b relate to the multiple order on source and receiver sides respectively.

co co <7 = 1

7/ co co λ

(Equation 4) [0030] Including source and receiver side convolutions, the 2D deconvolution will produce:

/(a) (Equation 5) ⁺²Σ^^αΣ^’ b---l a=l b----1 j

While source oniy and receiver only side multiples have been correctly 5 predicted, the mixed side multiples have been double counted. Figures 4a and 4b show the multiples predicted by the source and receiver side convolutions. Figure 4c shows the combined multiple estimate, highlighting the double prediction of the mixed side multiples.

[0032] To correct for the over prediction it is proposed here to iteratively 10 modify the data used for the convolutions using the prediction operator from the previous iteration. In the following equation, g^ⁿ~^ relates to the operator from the previous iteration and

is a new operator founded by least squares inversion:

(Equation 6) where:

and ,(/) , 1 v¹ in-ί} , .

3k ~^akk fk (Equation 6a) [0033] The factor of a half prevents over-prediction of the mixed side multiples. Note that the data on the left hand side of Equation 6 is still the original input data with no such modification. Hence, we are predicting original data from modified data.

[0034] In one embodiment, an iterative deconvolution prediction operator estimation scheme may consist of:

a) Receiving input data;

b) Deriving a standard mulfi-dimensional prediction operator;

c) Modifying input data using the prediction operator as indicated in

Equation 6a;

d) Calculating an improved 2D prediction operator using the modified input data (Equation 6);

e) Optionally going back to step c) to further improve the prediction operator.

[0035] Using the same notations as for Equation 1 given in the background section, the previous scheme may involve the following pseudo-code:,

1. Initialise operator go ~ 0

2. Calculate an updated prediction operator g_n by solving linear equation:

where the data d_s and dr used respectively for source and receiver convolutions (in operators S and R) are:

3. Go back to step 2 using the updated operator [0036] A multiple model computed using the aforementioned scheme may be more accurate than the standard multi-dimensional deconvolution approach, and as such may require less adaptive subtraction.

[0037] One may now consider multiples that relate to a primary corresponding to an event represented by the prediction operator. Considering the case where the prediction operator predicts only water bottom multiples, this would relate to arrivals that have only travelled in the wafer layer. This sub-set of the data may be described as a primary followed by receiver side multiples:

) [0038] Applying the modified prediction scheme (Equation 6) for peg-leg multiples, one can find:

h=l ό=2 I b---2 (Equation 7) [0039] Equation 7 shows that while multiple orders two and onwards have been predicted accurately, the first order multiple has been double counted. To solve this problem, it is proposed to scale the data in the interval relating to the multiple generator in Equation 7 by a half.

[0040] Thus, the final algorithm may be stated as follows:

(«)

(Equation 8) where h is the input data muted to the time interval of the prediction operator. [0041] A multiple model computed using the aforementioned scheme may be more accurate than the standard multi-dimensional deconvolution approach, and as such may require less adaptive subtraction.

[0042] The algorithm may be applied in two or more spatial dimensions where data sampling allows, A receiver-side 3D application may be used for towed streamer data where receiver side multiples are predicted using shot domain operators. The operators may be derived shot by shot, on a sliding spatial window of shots (for example 3 or 5 shots in the window), on overlapping spatial windows of shots along a sail line or on all shots within a sail line. Using more shots in the operator design may result in an operator that is more robust to the prediction of primary events, however it may also lack detail and as such may not model the multiples accurately enough. In this case, and by making equivalence between source and receiver side multiples, one may use equation:

d,

IK

Pik .-/09

--4, (Equation 9) [0043] This 3D scheme may be implemented for the case of equations 1, 3, 6, or 8.

[0044] De-aliasing of data in-between the streamers or extrapolation of data away from the streamers (e.g, reconstructing near offsets or extrapolating streamers outside the cable spread) may be implemented with data interpolation, for example using NMO copy and spatial weighting, hyperbolic Radon interpolation, tilted hyperbolic radon interpolation, fx interpolation, or another interpolation.

[0045] The algorithm may be used to output the full multiple model, the second order correction term and the term where the primary relates to the same interval as the multiple generator which may be used for joint or cascaded adaptive subtraction.

[0046] Another drawback of the deconvolution approach is directly linked to the deconvolution operator design. As stated above, the standard 2D gapped deconvolution process is modelled by Equation 4:

j j where input data, d, is predicted by a multi-channel convolution between the data and surface consistent prediction operator, g, with indices i, j, and k relating to spatial coordinates for the spatial summation plus primary, p.

[0047] The deconvolution prediction operator, g, is traditionally found by solving Equation 3, e.g. by inversion.

[0048] Prediction operator g is thus a function of time and space. Figure 5 gives an exemplary prediction operator g where the x-axis relates to space (the position relating to space at 0 relating to the apex of the operator), and the y-axis is time.

[0049] Depending on the length of the deconvolution operator (length of the gap and/or length of the active part), primaries can be either predicted or damaged. It is thus necessary to find solutions to constrain the deconvolution operator and make it less susceptible to the prediction of primaries.

[0050] The prediction operator may be constrained with sparseness weights to enforce a gap and operator length both of which may vary with aperture offset and/or with space (e.g. aperture apex position). In addition, the sparseness weights may be further modified using iteratively re-weighted least squares inversion. In one embodiment, a least squares solution is found to Equation 1, following which the envelope of the deconvolution operator is calculated and used to define sparseness weights for future sparseness iterations. The use of sparseness weights may reduce the likelihood of primary damage and wavelet truncation.

[0051] Equation 2 modified for operator domain sparseness weights may take the form:

A

W¹ C¹ d-^W¹ C^TCW g (Equation 10)

A where g = — w [0052] In addition to operator domain sparseness weight, data domain sparseness weights may also be used. Data domain sparseness weights may penalize missing or erroneous input data to avoid such data being detrimental to the quality of the output operators. Sparse inversion is explained in detail for instance in Trad et al., Latest views of the sparse Radon transform, 2003, Geophysics, 68, 386-399.

[0053] In another embodiment, it is proposed to derive the spatial deconvolution prediction operator directly in a model domain which is different to the space-time domain. This can be done by adding a reverse model transform to the deconvolution equation.

[0054] Linear Equations 1, 3, 6, 8, 9, or other prediction equations discussed herein, may be combined with a second linear transform L which converts the prediction operator derived from the new domain to the space-time domain, according to the following general equation.

£xt — Lg dom (Equation 11) where g_dom is the vector representing the deconvolution prediction operator derived in the new model domain.

g_xt is the vector representing the deconvolution prediction operator in the classical x-t domain; and

L is the linear operator for the conversion between the two domains. [0055] The reverse model transform operator L may relate to any operator transforming the model from a domain other than a space-time domain to the space-time domain. Linear operator L can be any one of:

® Reverse linear Radon operator ® Reverse parabolic Radon operator ® Reverse hyperbolic Radon operator ® Reverse migration operator (i.e. demigration) o Kirchhoff, beam, one way wave equation, RTM, etc o Constant or variable velocity field ® Reverse curvelet transform ® Reverse wavelet transform ® Reverse FK transform ® Reverse AVO operator (e.g. reforming a gather from R0 and G parameters) [0056] In one exemplary embodiment, the reverse model transform is a reverse hyperbolic Radon transform which models or represents data in the space-time (x-t) domain by a series of hyperbolic or hyperbolical events, for example defined in the (tO-V) domain (also known as tau-v domain). In this case, Equation 11 is now:

8xt — Lg-w where g,_v is the vector representing the deconvolution prediction operator g in the tau-v domain (also known as the hyperbolic Radon domain) g_xt is the vector representing the deconvolution prediction operator g in the classical x-f domain; and

Operator L is a matrix with (rx v) columns and (x x t) rows, each (x,t) position being a function of all points in the tau-v domain, as follows:

L(x, v, τ, v) = FFTfJreq - v to τ-ν domain)e~l™^J^FFT(τ - v to freq-v domain) where Δ/ ----- Jr² +^-, v being the RMS velocity of the sub-surface (m/s). V N

As both Equations 1 and 11 are linear operators, they can be applied in the following sequence as part of a least squares inversion using conjugate gradients:

1. Transform prediction operator from tau-v domain to x-t domain (Equation 11).

2. Convolve prediction operator with the data to estimate multiples (Equation 1).

[0057] The (tau-v) domain may be considered as a matrix of values, columns relating to v and rows relating to tau, each representing a different hyperbola in the (x-t) domain. The range of v may be from 1400 m/s to 6000 m/s, and may relate to the rms-velocity of sound in the subsurface. Alternatively, where an estimate of the velocity of the subsurface is known, the velocity may relate to a range of scaled versions of the velocity field estimate for example in the range 70% to 105% with 5% increment. As the velocity of sound for the multiple generator may be known (for example, dose to water velocity), the range of values in the (tau-v) domain may be reduced accordingly. In addition, the range of v as a function of tau may be variable. The reverse hyperbolic transform may consist of a mapping of energy from each (tau-v) position to each trace in the operator. The travel-time in each corresponding to t in the above equation (i.e. sqrt both sides of the above equation). As the travel-time value t will not fall exactly on a sample or time interval of the recorded trace, if may be necessary to interpolate the energy, for example using a sine function or Fourier interpolation.

[0058] By defining the prediction operator in the hyperbolic Radon domain, the prediction operator may only consist of events relating to hyperbolic moveout, which allows the operator to be well constrained, thus reducing the scope of predicting/damaging primaries. Figure 6 is an example of 2D prediction operator constrained by using the hyperbolic Radon transform.

[0059] The amplitude of the prediction operator may vary with space, and as such one embodiment may define a modified model domain called the AVO (amplitude variation with offset) hyperbolic Radon domain. Such a domain may now consist of two matrices. The first matrix may define the amplitude of the hyperbola at zero offset (i.e. for the central trace in the previous figure, also known as the intercept R0), and the second matrix may define how the amplitude of the hyperbolic event changes with space (known as the gradient, grad). In this embodiment, the reverse hyperbolic Radon transform is now modified so that the amplitude of the energy mapping to the (x-t) domain is defined by the following Equation:

Amp(f, x) = (R0 + x * grad)

- 2πϊίΔΐ where At is the hyperbolic shift as introduced earlier. The above equation may be reverse Fourier transformed to form data in the (x-t) domain. Instead of defining the variation of amplitude with space by an intercept amplitude and gradient, we may alternatively define an amplitude at one or more offsets, for example R0 (amplitude at zero offset), R250 (amplitude at 250m offset), and R500 (amplitude at 500m). The amplitude of the prediction operator in between the defined offsets may be calculated by a linear combination or linear interpolation between the derived values.

[0060] Figure 7 is an example of 2D prediction operator constrained by using the AVG hyperbolic Radon transform.

[0061] As stated above in reality a plurality of different transforms may be used to constrain the operators. In one embodiment, the operator may be defined as a parabolic Radon representation in the migrated domain. In another embodiment, the operator may be defined as a tau-px-qh representation in the migrated domain, where px is the slowness along a common offset section of the migrated operators, and qh is the parabolic radon representation in the offset direction.

[0062] Sparseness weights may also be used when the operator is defined in a model domain.

[0063] A new transform called the tilted hyperbolic Radon domain can also be used. The concept of a tilted hyperbola is consistent with travel times relating to a dipping plane in a constant velocity medium. For a horizontal recording datum, the linear operator for the tilted hyperbolic transform may be given by:

+

v + s x + s v (Equation 12)

X V^J where:

s is a model space parameter relating to the slowness in x-direction (s/m) s is a model space parameter relating to the slowness in y-direction (s/m) v

J x is the offset in the x-direction (m) of a given input trace y is the offset in the y-direction (m) of a given input trace v is a model space parameter that relates to the velocity of sound in the subsurface, e.g, RMS velocity (m/s) τ is a model space parameter relating to the zero offset arrival time (s) [0064] This might be implemented in 2D or 3D, The tilted hyperbolic domain may be found by solving the following equation:

d - ^Lthy(0)^mthy (Equation 13) where:

d is input data in the space-time domain, relating to a number of recordings (traces) that are a function of time, each having their own position in x and y. This data may relate any recording in a space-time domain, e.g. primaries only, multiples only, or a mixture of primaries and multiples;

m_th is a tilted hyperbolic domain for zero receiver depth datum, a function of τ, velocity, x-slowness, and y-slowness.;

^Lthy(O) is the linear operator reverse transforming from the tilted hyperbolic Radon domain to the space time domain, based on the travel-times given in Equation 12, with, for example, linear or sparse inversion (11, 10, Cauchy, or another norm, for example solved with conjugate gradients iteratively reweighted least squares).

[0065] The linear operator may be considered of a mapping of energy from the tilted hyperbolic Radon domain to a time on each trace in the space-time domain. For example, the mapping of energy may include the use of sine functions, linear interpolation, point to point mapping or phase shifts in a Fourier domain. Schemes used for such a mapping have been used in the case of hyperbolic Radon, as described above, or in Sacchi, M., and T. Ulrych, 1995,

High-resolution velocity gathers and offset space reconstruction: Geophysics, 60, 1169-1177.

[0066] The equation may be modified to account for data recorded on an arbitrary datum (e.g. non-horizontai) using at least one receiver depth (z) and vertical slowness (s ):

^k z ⁷ d = Lthy(z)m_thy Equation 14 where linear operator , L_thy(_Z), uses traveltimes:

i 2 X + N*

i.., ( > = h ( _λ + s x + s y + s z = . τ I--------------_s----------h s x + s y + s z thy(z) hy(x,y) x y-> z T _γ2 x z [0067] The vertical slowness, s , may be calculated using the relation:

dt , A ( , hy(*,y) ₊ dx +

/ where v is the speed of sound in medium where the receiver or source is positioned, for example the speed of sound in water, dtp^v/x (x ) and where derivatives--—— and ————may be calculated analytically.

dx dy

It should be noted that each receiver might have its own unique depth and the receivers need not all be at the same depth.

[0068] In one embodiment, the tilted hyperbolic model may be used for recreating data at new positions. In order to reconstruct the wavefield at new positions, we may reverse transform the tilted hyperbolic model by applying a linear operator derived using the output coordinates. This may be combined with any other embodiment.

[0069] The inversion equation 14 may be modified for the purposes of up/down separation or deghosting. In this embodiment linear operators for up-going and down-going arrivals are combined into a single inversion problem, for example:

^{d L}thy(z)t ^+RLthy(z)^] ^mthy (Equation 15) where:

R is the surface reflectivity (usually -1),

I.._thy(2)t is the reverse ti.ted hyperbolic transform for up-going wavefield, ^Lthy(z>

wavefield, and is the reverse tilted hyperbolic transform for down-going

Thy(z)t hy(x,y) x y

t.

t, ! X + S X + S V + s z ny(x,y) x v-' z ^Lthy(z)4r

3] The model satisfies a single wavefield as it reflects at the free-surface that is simultaneously consistent with the up-going (primary) and down-going (ghost) wavefields in the input data, [0071] Once the model has been found we may:

- Reverse transform the model to estimate the receiver ghost by applying operator RL,_hy(z)4..

- Reverse transform the model to estimate the up-going energy by applying operator L

The reverse transformed data may be subtracted from the input data. In one embodiment, the tilted hyperbolic model may be used for modeling source designature. We may modify the previous equation to optionally model the source farfleid signature, F:

d = F L . + RL , v thy(z): thy(z)v thy

In the case a source is made up of an array of N sources, each represented by a notional source, F_n, the equation becomes:

Σ n=l thy(z) + RL thv(z) thy [007S] The offsets used in the linear operator may also vary as a function of n, for example based on the position of an individual source element relative to the centre of the source.

[0076] Using this equation, the inversion may find a tilted hyperbolic model of the data relating to a point isotropic source, which is consistent with the recorded data relating to a source array below a free surface.

[0077] Once the model has been transformed, it may be:

o Reverse transformed using a Dirac source; so as to estimate data relating to a point source without a free surface, o Model the difference between the Dirac source and the source with the full array and free surface, This energy may be subtracted from the input data to correct for array and free surface effects.

o Model the difference between the source with bubble and free surface reflection and the source without bubble and free surface reflection. A model of the ghost and bubble energy may be subtracted from the input data.

[0078] In one embodiment, the tilted hyperbolic model may be used for demultiple. This may be achieved, for example, by muting primary energy in the tilted hyperbolic Radon domain, reverse transforming multiples, and subtracting the multiples from the input data.

[0079] Alternatively, we may output a wavefield relating to a new receiver datum to simulate the timing of peg-leg mutliples relating to a sea bed multiple generator. These multiples may then be subtracted (for example using adaptive subtraction) from the input data.

[0080] Whether the tilted hyperbolic transform for horizontal or variable datum has been used, the resulting model makes a continuous representation of the data and may be used for a variety of purposes such as:

® Data interpolation (e.g. cable interpolation prior to demuifipie or imaging);

® Data regularization;

« Data extrapolation (e.g. near offset extrapolation prior to demultiple;

® Redatuming;

® Wavefield separation:

o Deghosting o Up/Down separation ® Source designature;

® Noise attenuation.

[0081] The tilted hyperbolic transform may not adequately model all energy. As such, it may be necessary to subtract the modelled data from the received data to calculate residual data. The residual data may be further processed. For example, in the case of deghosting the residual may be processed using regular receiver deghosting algorithms (e.g. tau-p based), or in the case of redatuming the residual energy may be redatumed using normal moveout correction. Once processed, the processed residual data may be combined with data reverse transformed from the tilted hyperbolic Radon transform. [0082] Any known linear transform may be tilted by adding additional slowness terms in x- and / or y-directions; for example:

- Tilted parabolic Radon

- Tilted sinusoidal Radon

- Tilted 3D parabolic Radon; e.g. Q-R-S (Hugonnet et al., 2009, 3D high resolution parabolic Radon filtering, EAGE conference proceedings)

- Tilted conical Radon transform

- Tilted cylindrical Radon transform [0083] Figure 8 gives a flowchart of a method 200 implementing the principle of using a model domain other than the time-space domain. The method is particularly suitable for determining a subsurface model, such as a multiple model, in shallow water environments. The method applies equally well to fresh, brackish and salt water environments.

[0084] Method 200 begins with step 210, wherein seismic data are received and recorded in the space-time domain. Seismic data may be received by one or more receivers, such as sensor 3 shown in figure 1, after one or more sources, such as source 1 in figure 1, have transmitted a plurality of source wavelets.

[0085] In step 220, a gapped deconvolution prediction operator is defined from the received data. This is achieved by first, deriving a gapped prediction operator in a model domain which is different from said space-time domain (step 221), then converting the prediction operator derived in the model domain to the space-time domain using a linear conversion operator (step 222). The subsurface model can then be classically calculated, by convolving the space-time domain gapped prediction operator with the received seismic data in the space-time domain (step 223), [0086] As already disclosed above, a seismic processing method for determining a multiple mode! comprises:

receiving seismic data d containing primaries p and multiples m; deriving a gapped deconvolution prediction operator g from the received data (either directly in the space-time domain, or in a model domain as disclosed above); and determining the multiple model by convolving said gapped prediction operator with the received seismic data.

[0087] Once the multiple model has been determined, the estimated multiples m are generally subtracted to recorded data to deliver the primaries p, A straight subtraction or adaptive subtraction may be used as previously discussed.

[0088] This can be mathematically expressed by rearranging for instance Equation 3 as :

ik §ij jk

d..

' ij >jk (Equation 16) [0089] Improvement in the way primaries p, or more generally a set of data can be predicted by using the prediction operator g and the multiple model m.

[0090] In one embodiment, any set of data is predicted from a multiple by solving linear Equation 16, i.e. running an inversion on linear Equation 16.

[0091] Running an inversion means to finding set of data which, when convolved with the known prediction operator g, produces the known multiples m, i.e. the adapted multiples from the adaptive subtraction step.

[0092] Same reasoning can be made starting from Equation 6 or Equation 8 in case over prediction of either second order or higher order multiples and/or first order multiples has been corrected.

[0093] In one embodiment, finding a set of data r consisting of primary reflections from received data d containing primaries and multiples is achieved by running an inversion on the following equation:

d = Gr where G is an operator representing the convolution of the prediction operator with the data r on the source and the receiver sides.

For instance, G is a multi-order prediction operator of the type:

where ί_ς is the multiple order on the source side and i is the multiple order on the receiver side. Index zero on source and receiver sides corresponds to the primary. The non-zero indices encode all orders of multiple on source and receiver sides. It should be understood that while the above equation encodes all peg-leg multiples with consistent amplitude, any multiple relating to a convolution between the deconvolution operator and its corresponding primary will not be correct in terms of amplitude.

In this case the prediction operator may be used to encode all orders of multiple from the primary, r. For practical purposes we may only need to encode a finite number of multiple orders, for example 4. Note that instead of encoding each multiple individually, it may be faster to compound the calculations of different multiples.

[0094] An alternative approach is to receive multiples, m, a multi-dimensionai deconvolution operator, p, and run an inversion on the following equation to find data, d:

m --- Pd where P = P_g + P_R

This formulation will not be fully amplitude consistent for the reasons discussed earlier.

This approach could be implemented in the following workflow:

1) Receive input data containing primaries and multiples;

2) Estimate a multiple prediction operator using a multi-dimensional predictive deconvolution approach ;

3) Estimate a multiple model using the input data and the multiple prediction operator;

4) Optionally process the multiple model, for example using adaptive subtraction; and

5) Estimate primary reflection data using the multiple model and the prediction operator, where the step of estimating includes solving an inversion.

In particular this flow may be of interest in the case near offsets have not been recorded, for example:

1) Receive input data;

2) Reconstruct near offsets using a kinematic approach, e.g, NMO copy;

3) Invert for the prediction operator;

4) Convolve the operator with the data, including the near offset NMO copy traces;

5) Adapt the multiple model to the multiples recorded in the data;

6) Reconstruct near offset data (including primaries and multiples) from the adapted multiple model.

[0095] As an alternative in which over prediction of second or higher order multiples are corrected, finding a set of data d consisting of primary reflections and multiples from received data m containing multiples is achieved by running an inversion on equation:

where P

	P„			pc Ί
P - Pc	I--*-	+ P„	I -	s
s	2	R		2

[0096] As another alternative in which over prediction of first order multiple is achieved, data d is scaled by factor 1/2 for the time window relating to the duration of the prediction operator.

[0097] Another embodiment aims to find a set of data r consisting of primaries from a dataset m only including multiples. This is achieved by running an inversion on the following equation:

m = (G - l)r where G is one of the above-defined multi-order prediction operators, and subtracting I from G removes the primary from the output to leave just the multiples, [0098] Another embodiment can be used to find a dataset r consisting of primaries from a dataset that contains only a particular order of peg-leg multiple, mo. This is achieved by running an inversion on the following equation:

nig = Gr where

G = Gg + G_r in case the particular order of peg-leg multiple is the first order; or

G_SG_R + G_rG_r

GpGp in case the particular order of peg-leg multiple is the second order;

Etc,., [0099] The above-discussed procedures and methods may be implemented in a computing device as illustrated in Figure 9. Hardware, firmware, software or a combination thereof may be used to perform the various steps and operations described herein.

[00100]Exemplary computing device 1100 suitable for performing the activities described in the exemplary embodiments may include a server 1101. Such a server 1101 may include a central processor (CPU) 1102 coupled to a random access memory (RAM) 1104 and to a read-only memory (ROM) 1106. ROM 1106 may also be other types of storage media to store programs, such as programmable ROM (PROM), erasable PROM (EPROM), etc. Processor 1102 may communicate with other internal and external components through input/output (I/O) circuitry 1108 and bussing 1110 to provide control signals and the like. Processor 1102 carries out a variety of functions as are known in the art, as dictated by software and/or firmware instructions.

[00101]Server 1101 may also include one or more data storage devices, including hard drives 1112, CD-ROM drives 1114 and other hardware capable of reading and/or storing information, such as DVD, etc. In one embodiment, software for carrying out the above-discussed steps may be stored and distributed on a CD- ROM or DVD 1116, a USB storage device 1118 or other form of media capable of portably storing information. These storage media may be inserted into, and read by, devices such as CD-ROM drive 1114, disk drive 1112, etc. Server 1101 may be coupled to a display 1120, which may be any type of known display or presentation screen, such as LCD, plasma display, cathode ray tube (CRT), etc. A user input interface 1122 is provided, including one or more user interface mechanisms such as a mouse, keyboard, microphone, touchpad, touch screen, voice-recognition system, etc.

[00102]Server 1101 may be coupled to other devices, such as sources, detectors, etc. The server may be part of a larger network configuration as in a global area network (GAN) such as the Internet 1128, which allows ultimate connection to various landline and/or mobile computing devices.

[00103]Computing device or computing apparatus 1100 can be configured to implement any of the above-discussed procedures and methods, including combinations thereof. For instance, interface 1008 is configured to receive seismic data recorded in a space-time domain; and data processing unit 1002 may be configured to define a gapped deconvolution prediction operator from the received seismic data; and to determine the subsurface model, especially the multiple model, by convolving said gapped deconvolution prediction operator with the received seismic data, wherein said gapped deconvolution prediction operator is defined by deriving said gapped deconvolution prediction operator in a model domain which is different from said space-time domain, and by converting the deconvolution prediction operator derived in the model domain to the space-time domain using a linear conversion operator.

[00104]The disclosed exemplary embodiments provide a computing device, software instructions and a method for seismic data processing. If should be understood that this description is not intended to limit the invention. On the contrary, the exemplary embodiments are intended to cover alternatives, modifications and equivalents, which are included in the scope of the invention as defined by the appended claims. Further, in the detailed description of the exemplary embodiments, numerous specific details are set forth in order to provide a comprehensive understanding of the claimed invention. However, one skilled in the art would understand that various embodiments may be practiced without such specific details.

[00105]Although the features and elements of the present embodiments are described in the embodiments in particular combinations, each feature or element can be used alone without the other features and elements of the embodiments or in various combinations with or without other features and elements disclosed herein.

[00106]This written description uses examples of the subject matter disclosed to enable any person skilled in the art to practice the same, including making and using any devices or systems and performing any incorporated methods.

The patentable scope of the subject matter is defined by the claims, and may include other examples that occur to those skilled in the art. Such other examples are intended to be within the scope of the claims.

Claims (15)

1. A seismic processing method (200) for determining a subsurface model comprising:

receiving (210) seismic data recorded in a space-time domain; defining (220) a gapped deconvolution prediction operator from the received seismic data; and determining (223) a subsurface model by convolving said gapped deconvolution prediction operator with the received seismic data, wherein said gapped deconvolution prediction operator is defined by deriving (221) said gapped deconvolution prediction operator in a model domain which is different from said space-time domain, and by converting (222) the deconvolution prediction operator derived in the model domain to the spacetime domain using a linear conversion operator.

2. The is a reverse

3. The is a reverse

4. The is a reverse

5. The is a reverse

6. The is a reverse

7. The is a reverse seismic method of claim 1, wherein said linear conversion operator linear Radon operator.

seismic method of claim 1, wherein said linear conversion operator parabolic Radon operator.

seismic method of claim 1, wherein said linear conversion operator hyperbolic Radon operator.

seismic method of claim 1, wherein said linear conversion operator migration operator.

seismic method of claim 1, wherein said linear conversion operator curvelet transform.

seismic method of claim 1, wherein said linear conversion amplitude variation with offset, AVO, operator.

operator

8. The seismic method of claim 1, wherein said linear conversion operator is a reverse amplitude variation with offset, AVO, hyperbolic Radon operator.

9. The seismic method of claim 1, wherein said linear conversion operator is a reverse tilted hyperbolic Radon operator.

10. The seismic method of claim 1, wherein said subsurface model is a multiple model.

11. The seismic method of claim 1, wherein said subsurface model is a primary model.

5

12. The seismic method of claim 1, wherein said subsurface model contains modelled primaries and multiples.

13. The seismic method of claim 1, wherein the subsurface model is subtracted from the received seismic data.

14. The seismic method of claim 1, wherein the seismic data is modified so

10 as to improve the prediction of the multiples.

15. A data processing apparatus (1000) for determining a subsurface model, the apparatus comprising:

an interface (1008) configured to receive seismic data recorded in a space-time domain; and

15 a data processing unit (1002) configured to define a gapped deconvolution prediction operator from the received seismic data; and to determine the subsurface model by convolving said gapped deconvolution prediction operator with the received seismic data, wherein said gapped deconvolution prediction operator is defined by deriving

20 said gapped deconvolution prediction operator in a model domain which is different from said space-time domain, and by converting the deconvolution prediction operator derived in the model domain to the space-time domain using a linear conversion operator.

Intellectual

Property

Office

Application No: GB1710708.7 Examiner: Stephen Jennings