GB2605999A

GB2605999A - Method of, and apparatus for, geophysical investigation using seismic signal decomposition

Info

Publication number: GB2605999A
Application number: GB2105706.2A
Authority: GB
Inventors: Wang Yanghua
Original assignee: Individual
Current assignee: Individual
Priority date: 2021-04-21
Filing date: 2021-04-21
Publication date: 2022-10-26
Anticipated expiration: 2041-04-21
Also published as: GB2605999B; GB202105706D0

Abstract

A method of subsurface exploration comprises determining one or more physical properties of a portion of the volume of the Earth from a seismic measurement representative of at least one physical parameter. The method comprises: providing an observed multichannel seismic dataset derived from a seismic measurement of a portion of the volume of the Earth; estimating local dipping information of the portion of the volume of the Earth from the observed multichannel seismic dataset; defining the observed multichannel seismic dataset as a multichannel residual dataset; iteratively performing seismic signal decomposition on each channel of the multichannel residual dataset simultaneously to generate an output seismic dataset comprising a series of decomposed wavelets; and providing the output seismic dataset to determine one or more physical properties of the portion of the volume of the Earth for subsurface exploration.

Description

Method of, and apparatus for, Geophysical Investigation Using Seismic Signal Decomposition The present invention relates to an improved method of geophysical investigation.

More particularly, the present invention relates to an improved method of geophysical investigation using seismic signal decomposition.

There is much interest in surveying the subsurface of the Earth to detect natural mineral resources or other sites of geological interest using seismic surveys. Seismic to surveys are the key means for the petroleum industry to explore the subsurface of the Earth for hydrocarbon reserves.

In general, seismic survey data is acquired and analysed with regard to identifying locations suitable for direct investigation of the sub-surface by drilling. Seismic surveying also has applications within the mining industry and within other industrial sectors that have an interest in details of the subsurface of the Earth.

In a seismic survey, seismic sources are used to generate vibrational energy which is directed into the subsurface of the Earth. Reflected, refracted and other acoustic signals returned from subsurface features are then detected and analysed to map the subsurface of the Earth.

A schematic illustration of an experimental set up 10 for an undersea seismic survey is shown in Figure 1. This example is non-limiting and an equivalent experiment can be carried out on land. The present invention is applicable to subsurface exploration in any suitable environment, for example land or marine measurements of a portion of the subsurface of the Earth. The present invention may be applicable to identification of numerous subsurface resources in any form and is intended to include oil exploration and gas prospecting.

The skilled person would be readily aware of the suitable environments in which data could be gathered for analysis and exploration purposes as set out in the present disclosure.

In this example, the experimental set up 10 comprises a source 12 located on a ship 14. However, this need not be the case and the source may be located on land, or within the sub-surface, or on any other suitable vessel or vehicle.

The source 12 generates acoustic and/or elastic waves having sufficient vibrational energy to penetrate the subsurface of the Earth and generate sufficient return signals to aid useful detection. The source 12 may comprise, for example, an explosive device, or alternatively an air gun or other mechanical device capable of creating sufficient vibrational disturbance. Commonly, for many seismic survey experiments a single source is used which is shot from multiple locations. Naturally occurring sources may also be employed.

A plurality of detectors 16 is provided, where the detectors 16 may comprise any suitable vibrational detection apparatus. Commonly, geophones (which detect particle motion) and hydrophones (which detect pressure variations) are used. It is typical for a large number of detectors 16 to be laid out in lines for two-dimensional (2D) data acquisition although the detectors 16 may be arranged in sets of lines or in a grid for three-dimensional (3D) data acquisition. Detectors 16 may also be placed within the subsurface, for example down boreholes. The detectors 16 are connected to trace acquisition apparatus such as a computer or other electronic storage device. In this example, the acquisition apparatus is located on a further ship 18. However, this need not be the case and other arrangements are possible.

In use, elastic waves 20 generated by the source 12 propagate into the subsurface 22 of the Earth. The subsurface 22, in general, comprises one or more layers or strata 24, 26, 28 formed from rock or other materials. The elastic waves 20 are transmitted and refracted through the layers and/or reflected off the interfaces between them and/or scattered from other heterogeneities in the sub-surface and a plurality of return signals 30 is detected by the detectors 16.

In general, the returning signals 30 comprise elastic waves having different polarisations. Primary or compression waves (known as P-waves) are approximately longitudinally polarised and comprise alternating rarefactions and compressions in the medium in which the wave is travelling. In other words, in an isotropic environment, the oscillations of a P-wave are parallel to the direction of propagation.

Shear or secondary waves (known as S-waves) may also be generated. S-waves have an approximately transverse polarisation. In other words, in an isotropic environment, the polarisation is perpendicular to the direction of propagation. Whilst S-wave analysis is possible and falls within the scope of the present invention, the following description will focus on the analysis of P-waves.

to A seismic survey is typically composed of a large number of individual source excitation events. The Earth's response to these events is recorded at each receiver location, as a seismic trace for each source-receiver pair. For a 2D survey, the tens of thousands of individual traces may be taken. For the 3D case, this number may run into the millions.

A seismic trace comprises a sequence of measurements in time made by one or more of the multiplicity of detectors 16, of the returning reflected, refracted and/or scattered elastic waves 30 originating from the source 12. In general, a partial reflection of the elastic wave 20 occurs at a boundary or interface between two dissimilar materials, or when the elastic properties of a material changes. Traces are usually sampled in time at discrete intervals of the order of milliseconds.

Seismic surveys at the surface or seabed can be used to extract rock properties and construct reflectivity images of the subsurface. Such surveys can, with the correct interpretation, provide an accurate picture of the subsurface structure of the portion of the Earth being surveyed. This may include subsurface features associated with mineral resources such as hydrocarbons (for example, oil and natural gas). Features of interest in prospecting include: faults, folds, anticlines, unconformities, salt domes, reefs. An example of such a feature or anomaly 32 is shown in Figure 1.

Multiple techniques for the analysis of seismic data are known. One technique is seismic matching pursuit. Seismic matching pursuit is a technique to decompose a seismic trace into a stack of wavelets. The decomposition is performed iteratively. At each iteration, an optimal wavelet is chosen so that it has the highest correlation coefficient with the current residual trace. Because the time-frequency spectrum generated by those decomposed wavelets has relative high resolution, matching pursuit method has gained much popularity in seismic field, including seismic wavefield reconstruction, and hydrocarbon reservoir detection.

However, matching pursuit suffers from high computational burden caused by extracting the optimal wavelet. Various approaches to speed up the process are known, for example by using the analytical forms of Morlet wavelets, by extracting multiple basic wavelets through parallel computation, or through the introduction of to quantum swarm evolutionary matching pursuit to boost the efficiency of the conventional matching pursuit.

However, a technical problem with such approaches is that the classic matching pursuit method suffers from the non-uniqueness in seismic trace decomposition and IS meanwhile has high sensitivity to random noise.

The present invention, in embodiments, addresses this technical problem.

According to a first aspect of the present invention there is provided a method of subsurface exploration, the method comprising determining one or more physical properties of a portion of the volume of the Earth from a seismic measurement representative of at least one physical parameter, the method comprising: a) providing an observed multichannel seismic dataset derived from a seismic measurement of a portion of the volume of the Earth; b) estimating local dipping information of the portion of the volume of the Earth from the observed multichannel seismic dataset; c) defining the observed multichannel seismic dataset as a multichannel residual dataset; d) iteratively performing seismic signal decomposition on each channel of the multichannel residual dataset simultaneously to generate an output seismic dataset comprising a series of decomposed wavelets; and e) providing the output seismic dataset to determine one or more physical properties of a portion of the volume of the Earth for subsurface exploration.

In one embodiment, step d) comprises: f) extracting a basic wavelet from the multichannel residual dataset; g) expanding a basic-wavelet matrix; h) calculating an amplitude of each wavelet for each of the multiple channels simultaneously to generate a synthetic multichannel seismic dataset using a dipping constraint derived from said estimated local dipping information; i) forming an updated multichannel residual dataset by subtracting the synthetic multichannel seismic dataset from the observed multichannel seismic dataset; j) repeating steps f) to i) using the updated multichannel residual dataset until a criterion of a threshold parameter is met; and, if said criterion is met; k) generating an output seismic dataset comprising a series of decomposed wavelets.

to In one embodiment, step h) further comprises utilising an objective function having said dipping constraint.

In one embodiment, the dipping constraint is operable to constrain solutions of the objective function towards extracted wavelets which are optimal for seismic signals of IS multiple channels along a dipping orientation derived from said estimated dipping information.

In one embodiment, the dipping constraint is a soft dipping constraint.

In one embodiment, the objective function is defined in terms of a Frobenius Norm.

In one embodiment, the objective function is defined as J=11s,,,r -s],kAkx,t+ii 11(),,T.D,S4 +F("2"Dr(Lt_, where M is the number of data channels, s is the multichannel seismic dataset, and A", is an amplitude matrix in which each column vector is composed of the amplitudes of a single channel.

In one embodiment, step h) further comprises minimizing the objective function.

In one embodiment, minimising the objective function comprises calculate all wavelet amplitudes in a least-squares approach.

In one embodiment, the step of minimising estimates the current wavelet amplitude and refines the amplitudes of wavelets extracted in all previous iterations.

In one embodiment, the dipping constraint comprises one or more local dipping operators applied to the reconstruction of multichannel synthetic signals in the form of -th4"' In one embodiment, the degree of constraint of the dipping constraint is variable.

In one embodiment, the soft dipping constraint comprises first order differential operators D and D, along the time axis and the space axis respectively, defined as: -1 1 -1 1 D= DL= -1 1 where Q is 1.1>«P4 -1 a tridiagonal matrix with the sine elements, Q = cliag{sin01 sin/9z * Q-is a tridiagonal matrix with cosine elements, Q = diag{rosel cosq_ *..

and 0, are the dipping angles along the time taken from a central channel of the multiple channels.

In one embodiment, step f) further comprises extracting the basic wavelet from the defined multichannel residual dataset utilising structure-adaptive multichannel matching pursuit methods.

In one embodiment, the structure-adaptive multichannel matching pursuit comprises the steps of: I) defining the Rh basic wavelet by one or more of the parameters of: mean frequency, wavelet breadth, phase, and time delay; m) setting an initial value of the or each parameters; and n) generating a structure-adaptive value of the multichannel residual dataset comprising the sum along a local dip of a coherent event.

In one embodiment, step m) is set according to the parameters of an average channel from the multichannel residual dataset and step n) generates a structure-adaptive value.

In one embodiment, step g) comprises: o) adding an extra column vector, formed by the basic wavelet extracted at the current iteration of the robust multichannel matching pursuit for seismic signal decomposition.

In one embodiment, in step o) the extra column vector is of the form s",,=s",l1g," where N is the number of samples of each data channel, e,x,, is the basic-wavelet matrix obtained in the previous iteration, g, is the vector of the current basic wavelet, and u is the augment operator which combines the basic-wavelet matrix s and the vector g, . In one embodiment, step i) comprises subtracting the synthetic multichannel dataset from the observed multichannel seismic dataset, =s,,, -§,;t;" , where 1C" is the multichannel residual dataset after the leh iteration.

In one embodiment, the threshold parameter comprises the relative energy of the multichannel residual dataset and the respective criterion is that the relative energy is lower than a pre-set percentage threshold.

In one embodiment, step b) comprises defining the local wavefield as local plane wave and utilising plane wave destruction methods to estimate the local dipping information.

In one embodiment, the method is performed on data in a plurality of windows taken along a length direction.

In one embodiment, step e) further comprises utilising the output seismic dataset for subsurface exploration. In one embodiment, step e) further comprises utilising the output seismic dataset to determine one or more physical properties of a portion of the volume of the Earth as part of subsurface exploration.

According to a second aspect of the present invention, there is provided a computer system comprising a processing device configured to perform the method of the first aspect.

According to a third aspect of the present invention, there is provided a computer readable medium comprising instructions configured when executed to perform the method of the first aspect.

According to a fourth aspect of the present invention, there is provided a computer system comprising: a processing device, a storage device and a computer readable medium of the third aspect.

In embodiments, the present invention enables the basic wavelet in a local region to be obtained along a dipping plane defined by local structure dipping information. This local dipping information is casted as a constraint to a multichannel least-square inversion problem. The constraint may, in embodiments, be a soft constraint.

In embodiments, the present invention provides a method of robust multichannel matching pursuit for seismic signal decomposition. In embodiments, seismic signal decomposition is performed as a multichannel matching pursuit. Any extracted basic wavelet is the optimal basic wavelet suitable to all channels within a group of multiple channels. In embodiments, multichannel matching pursuit is performed as a multichannel inverse problem. The wavelet amplitudes are calculated simultaneously for all wavelets and for all channels within a group of multiple channels.

In embodiments, multichannel seismic signal decomposition is performed as a multichannel inverse problem with dipping constraint. The extracted wavelets are optimal for seismic signals of multiple channels along the dipping orientation in a soft-constrained manner.

In embodiments, the present invention is directed to seismic exploration and, more particularly, to a method of robust and multichannel matching pursuit for seismic signal decomposition.

The same procedure described in embodiments of the present invention may be performed along sliding windows following the space direction, and the output of each sliding window is the decomposition result of a central channel.

In embodiments, the first step of estimating local dipping information from the original multichannel seismic dataset. In embodiments, the present invention exploits a unique feature of multichannel seismic signals, and this unique feature is the spatial coherence of seismic reflection signals. The local wavefield may be treated as local plane wave, and the plane-wave destruction method may be used to obtain the local to dipping information.

In embodiments, the second step of treating the original multichannel seismic dataset as a multichannel residual dataset. At the beginning of iterative implementation, there is no extracted wavelet, the multichannel residual signals are the original multichannel seismic signals.

In embodiments, the third step of extracting a basic wavelet from the current multichannel residual dataset, implemented as structure-adaptive multichannel matching pursuit. The kth basic wavelet gy, is defined by four parameters, n =fon". o'"), 0("). trcul: the mean angular frequency, the parameter controlling the wavelet breadth, the phase, and the time delay. The initial value of this set of parameters is given according to an average channel from the multichannel residual dataset. A structure-adaptive average of the multichannel residual dataset is the sum along the local dip of a coherent event. These four parameters are refined simultaneously over multiple seismic channels.

In embodiments, the fourth step of expanding the basic-wavelet matrix by augmenting an extra column vector, formed by the basic wavelet extracted at the current iteration: lkik-111--igx where N is the number of samples of each seismic channel, (1!"4__il is the basic-wavelet matrix obtained in the previous iteration, g, is the vector of the current basic wavelet, and u is the augment operator which combines the basic-wavelet matrix and the vector In embodiments, the fifth step of calculating all wavelet amplitudes of all multiple channels simultaneously, implemented as a multichannel inverse problem with a soft dipping constraint. The objective function is defined in terms of the Frobenius norm as -@A ±g[11Q-Dis kit ±PL,DAQ-11:], where M is the number of seismic channels, s is the multichannel seismic dataset, and k,", is the amplitude matrix, in which each column vector is composed of the amplitudes of a single channel. The local dipping operators are applied to the reconstruction of synthetic multichannel signals sT!,-x) = * A. The degree of constraint is controlled softly by the trade-off parameter p. In embodiments, the fifth step of calculating all wavelet amplitudes of all multiple channels simultaneously; In the soft dipping constraint, Dr and D, are the first-order differential operators along the time (t) axis and the space (x) axis respectively, D.= -1 1 - (N-1;XN - IM 1) Q is a tridiagonal matrix with the sine elements, Q = diagIsir.61 sin 02, Q 0, is a tridiagonal matrix with cosine elements, g = diag {cost), cost92 cosi& } , and 0, are the dipping angles alongsor the time taken from the central channel of the multiple channels.

In embodiments, the fifth step of calculating all wavelet amplitudes of all multiple channels simultaneously. In the definition of the objective function above, lirt lt, is the Frobenius norm of matrix R,

D

where r are the elements of R. Because it is a sum of all squared elements of a matrix, thus, minimising the objective function means to calculate all wavelet amplitudes in a least-squares creation. The instant invention not only estimates the current wavelet amplitude but also refine the amplitudes of wavelets extracted in all previous iterations.

In embodiments, the sixth step of forming a new multichannel residual dataset by subtracting the synthetic multichannel signals from the original multichannel seismic dataset, 1C"-S,"-§(2,14, , where R4 is the multichannel residual dataset after the kth iteration.

In embodiments, the seventh step of performing steps three to six iteratively until the relative energy of the multichannel residual dataset is lower than a pre-set threshold in pre centage. The scaled per centage threshold reduces the dependency on the original amplitude of the signal, but it depends upon seismic data noise.

In embodiments, the eighth step of outputting a series of decomposed wavelets of the central channel. The same procedure described in the instant invention is performed along sliding windows following x direction, and the output of each sliding window is the decomposition result of the central channel.

The present invention, in embodiments, mitigates the effect of seismic data noise and any errors in estimating dip information, leads to robust and accurate multichannel seismic signal decomposition.

Embodiments of the present invention will now be described with reference to the accompanying figures, in which: Figure 1 shows a schematic diagram of a seismic measurement process to generate seismic trace data; Figure 2 is a flowchart outlining the eight implementation steps of the instant invention; Figure 3a shows a set of seismic channels; Figure 3b is the corresponding amplitude spectrum (at 15 Hz) generated by single channel matching pursuit; Figure 3c is the corresponding amplitude spectrum (at 15 Hz) generated by the to method of an embodiment of the present invention utilising multi-channel matching pursuit with a dipping constraint; Figure 4a is a group of seismic channels, the average channel without considering dipping information and the average channel which take account the dipping information; Figure 4b is the matching pursuit decomposition of the average channel without considering dipping information; Figure 4c is the matching pursuit of the average channel which take account of the dipping information; Figure 5 compares two convergence rates -the moralised residual energy versus iteration number of standard MC-MP (dotted line) and the instant invention, MC-MP with structure adaptative implementation (triangle-dotted line); and Figure 6 compares two residuals, in which Figure 6a is the residual channels after the standard MC-MP (without the structure adaptative implementation) after 100 iterations, and Figure 6b is the residual channels after using the structure-adaptative MC-MP of an embodiment after 20 iterations.

The present invention relates to a method of robust multichannel matching pursuit for seismic signal decomposition. In embodiments, seismic signal decomposition is performed as a multichannel matching pursuit. Any extracted basic wavelet is the optimal basic wavelet suitable to all channels within a group of multiple channels.

In embodiments, multichannel matching pursuit is performed as a multichannel inverse problem. The wavelet amplitudes are calculated simultaneously for all wavelets and for all channels with a group of multiple channels.

Embodiments of the present invention are implemented utilising the following generalised steps: 1) Obtaining an observed multichannel seismic dataset derived from a seismic measurement of the portion of said portion of the volume of the Earth; 2) Estimating local dipping information from the original multichannel seismic dataset; 3) Treating the original multichannel seismic dataset as a multichannel residual dataset, and iterating the following steps for seismic signal decomposition; 4) Basic-wavelet extraction: extracting a basic wavelet from the current multichannel residual dataset, implemented as structure-adaptive multichannel matching pursuit; 5) Basic-wavelet matrix augment: augmenting the basic-wavelet matrix by adding an extra column vector, formed by the basic wavelet extracted at the current iteration; 6) Amplitude calculation: calculating all wavelet amplitudes of all multiple channels simultaneously, implemented as a multichannel inverse problem with a soft dipping constraint; 7) Multichannel residual dataset: forming a new multichannel residual dataset by subtracting the synthetic multichannel signals from the original multichannel seismic dataset 8) Iteration: performing steps 4-7 iteratively until the relative energy of the multichannel residual dataset is lower than a pre-set threshold in percentage terms; and 9) Output: outputting a series of decomposed wavelets. The wavelets may be of the central channel or a different channel.

to It is noted that a multichannel data set is a data set comprising a plurality of data channels.

The same procedure described in the following embodiments may be performed along sliding windows following the space direction, and the output of each sliding window may be the decomposition result of the central channel.

A method according to the present invention will now be described with reference to Figure 2. Figure 2 shows a flow diagram of an embodiment of the present invention. In Figure 2, steps 200 to 208 are illustrated as S200 to S208.

Step 200: Obtain observed seismic data set Initially, it is necessary to obtain a set of experimentally gathered seismic data in order to initiate subsurface exploration. This may be gathered by an experimental arrangement such as the set up shown and described with reference to Figure 1.

The gathered seismic data may be optionally pre-processed in various ways including by propagating numerically to regions of the surface or subsurface where experimental data have not been acquired directly. The skilled person would readily be able to design and undertake such pre-processing as might be necessary or desirable. With or without such pre-processing, the resultant seismic dataset representing experimentally-gathered data is known as an "observed seismic data set".

As shown in Figure 1, a large number of receivers or detectors 16 are positioned at well known positions on the surface of the portion of the Earth to be explored. The detectors 16 may be arranged in a two dimensional (such as a line) or a three dimensional (such as a grid or plurality of lines) arrangement. The physical location of the detectors 16 is known from, for example, location tracking devices such as GPS devices. Additionally, the location of the source 12 is also well known by similar location tracking means.

The observed seismic data set may comprise multiple source 12 emissions known in the art as "shots". The data comprises pressure as a function of receiver position (on to the x-axis) with respect to time (on the t-axis). This is because, in general, a detector such as a hydrophone measures the scalar pressure at its location. However, other arrangements may be used.

The seismic trace data comprises a plurality of observed data points. Each measured discrete data point has a minimum of seven associated location values -three spatial dimensions (x, y and z) for receiver (or detector) position (r), three spatial dimensions (x, y, z) for source location (s), and one temporal dimension measuring the time of observation relative to the time of source initiation, together with pressure magnitude data. The seven coordinates for each discrete data point define its location in space and time.

The seismic trace data also comprises one or more measurement parameters which denote the physical property being measured. In this embodiment, a single measurement parameter, pressure is measured. The observed data set is defined as dobs(r,s,t) and, in this embodiment, is in the time domain. For clarity, the following discussion considers a single source-receiver pair and so r, s are not needed.

However, it is possible to measure other parameters using appropriate technology; for example, to measure the particle velocities or particle displacements in three spatial dimensions in addition to the pressure. The present invention is applicable to the measurement of these additional variables.

The actual gathering of the seismic data set is described here for clarity. However, this is not to be taken as limiting and the gathering of the data may or may not form part of the present invention. The present invention simply requires a real-world observed seismic data set upon which analysis can be performed to facilitate subsurface exploration of a portion of the Earth.

The method now proceeds to step 201.

Step 201: Estimate local dipping information At step 201, local dipping information is estimated from the original multichannel to seismic dataset.

Dipping information relates to the angle to the horizontal of a reflector such as a geophysical layer which results in a systematic change in arrival time due to the dipping layer.

In embodiments, a specific and unique feature of multichannel seismic signals is exploited, namely the spatial coherence of multichannel seismic signals. The local wavefield may be treated as local plane wave, and the plane-wave destruction (PWD) method may be used to obtain the local dipping information.

For example, the local dipping information may be obtained using the plane-wave destruction method where the discretised plane-wave destructor is expressed in a vector-matrix form as: f(cr,)=e7LIC(o-,) Lie =0, where a, is the apparent dipping along the x coordinate, c(0-,) is a two-dimensional operator, u is the matrix containing seismic data samples, [c(0-Jou] is the Hadamard product of two matrices, and e is the all-ones vector. The dipping is estimated by solving this destructor with a least-squares inversion method.

In the alternative, different plane wave destruction methods may be used. For example, linear plane wave destruction may be used which comprises application of explicit finite differences and obtains the slope field by least-squares estimation. A non-linear PWD could also be used which applies a max-flat fractional delay filter to approximate a linear phase operator (or phase-shift operator) and to provide a polynomial equation for the local slope.

When the local dipping information has been estimated, the method proceeds to step 202.

Step 202: Initialisation to At step 202, the original multichannel seismic dataset provided in step 200 is treated as a multichannel residual dataset. The original dataset is then utilized in the following steps for seismic signal decomposition. By this is meant that at the beginning of the iterative process of steps 203 to 207, there is no extracted wavelet and the multichannel residual signals are the original multichannel seismic signals are the original multichannel seismic signals generated or provided in step 200.

The method proceeds to step 203.

Step 203: Basic wavelet extraction At step 203, a basic wavelet is extracted from the current multichannel residual dataset (which at the first run through of the process comprises the original multichannel seismic dataset as described in step 202 above).

This is implemented as structure-adaptive multichannel matching pursuit. The kth basic wavelet g, is defined by four parameters: ={: Q.) which comprise: i) the mean angular frequency, ii) the parameter controlling the wavelet breadth, iii) the phase, and iv) the time delay.

The initial value of this set of parameters is given according to an average channel from the multichannel residual dataset. A structure-adaptive average of the multichannel residual dataset is the sum along the local dip of a coherent event.

These four parameters are refined simultaneously over multiple seismic channels.

The method then proceeds to step 204.

Step 204: Basic wavelet matrix augment At step 204, the basic-wavelet matrix is augmented by adding an extra column vector, formed by the basic wavelet extracted at the current iteration: eArxk-IPAT4k 111-4

A

where N is the number of samples of each seismic channel, 0,0_" is the basic-wavelet matrix obtained in the previous iteration, g7 is the vector of the current basic wavelet, and u is the augment operator which combines the basic-wavelet matrix ILkk_ip and the vector g.

The method proceeds to step 205.

Step 205: Amplitude calculation At step 205, in embodiments, the wavelet amplitudes of all multiple channels are calculated simultaneously. In embodiments, this is implemented as a multichannel inverse problem with a soft dipping constraint.

The objective function is then defined in terms of the Frobenius norm as: -(1),.kAk.112, +11[11Q1iikk%42, where M is the number of seismic channels, s],x, is the multichannel seismic dataset, and is the amplitude matrix, in which each column vector is composed of the amplitudes of a single channel.

In general terms, the Frobenius norm defines an operation in which the columns of the matrix are stacked on top of each other to create a vector of size m x n, and then the vector 2-Norm of the result is taken.

The local dipping operators are applied to the reconstruction of multichannel synthetic signals S, = 4, . The degree of constraint is controlled softly by the trade-off parameter p. The inventor of the present application has recognised, for the first time, that numerous benefits can be derived from, in embodiments, defining the objective function in terms of the Frobenius norm. This has the following advantages: 1. Amplitude determination: The first term of the objective function (equation 4) relates to data fitting. Conventional matching pursuit methods determine the amplitude of each wavelet in sequence. However, in embodiments of the present invention, the amplitudes of all wavelets are estimated simultaneously. Such a simultaneous estimation is achieved by an inversion method which maximises data fit (equivalently minimises the misfit).

2. Multichannel determination: Even for a simultaneous estimation of the amplitudes of all wavelets, conventional methods use a single channel, which is a data vector, rather than a data matrix.

In data fitting terms, the observed dataset is a matrix), the synthetic dataset is generated by a matrix (43,,A"." ). In embodiments, the present invention utilises multiple channels simultaneously, for which the data fitting is defined in terms of the Frobenius norm.

3. Construction of the dipping information: The second term in the objective function comprises a unique approach to constructing the dipping information. In embodiments, the present invention constructs the dipping information in a matrix form, such that it is consistent with the matrices in the data-fitting term.

4. Soft constraint: In this embodiment, the dipping information in successfully constructed in a consistent form as the data matrix. Thus a soft constraint parameter can be assigned to the system. A soft constraint is one which has one or more variable value that is penalised in the objective function if the conditions on the one or more variable is not satisfied.

In embodiments, all wavelet amplitudes of all multiple channels are calculated simultaneously. In the soft dipping constraint, D, and D, are the first order derivative operators along the time axis and the space axis respectively, 1 -1 D= D = -1 1 (1*1-1)xx -1,1xUd I) Q is a tridiagonal matrix with the sine elements, Q sr, = g{s±91 5±102 * * * sjna, } , Q cos is a tridiagonal matrix with cosine elements, Q spri = cliag loose, c.os02 cose,_,} , and 0, are the dipping angles along the time taken from the central channel of the multiple channels.

In embodiments, the calculation of all wavelet amplitudes of all multiple channels is performed simultaneously. In the definition of the objective function above, IR lk is the Frobenius norm of matrix R,

I IR

where r are the elements of R. Because it is a sum of all squared elements of a matrix, thus, minimising the objective function means to calculate all wavelet amplitudes in a least-squares creation. In embodiments, the method not only estimates the current wavelet amplitude but also refine the amplitudes of wavelets extracted in all previous iterations.

The method proceeds to step 206.

Step 206: Multichannel residual dataset At step 206, a new multichannel residual dataset is formed by subtracting the synthetic multichannel signals from the original multichannel seismic dataset, R1 = -S1, where it,^;:<A, is the multichannel residual dataset after the kth iteration.

The method proceeds to step 207.

Step 207: Iteration At step 207, steps 203 to 206 are repeated iteratively until a threshold is met. In embodiments, the threshold metric utilised is that the relative energy of the multichannel residual dataset is lower than a pre-set threshold in percentage.

In embodiments, the scaled per centage threshold reduces the dependency on the original amplitude of the signal. However, it depends upon seismic data noise.

If the threshold has not been met, the method proceeds back to step 203. If the threshold criterion/criteria has/have been met, the method proceeds to step 208.

Step 208: Utilise data Once a threshold is reached in step 208, the method of an embodiment is operable to output a series of decomposed wavelets. These may, in embodiments, be of the central channel.

The same procedure may be performed along sliding windows following the space direction, and the output of each sliding window is the decomposition result of the central channel.

The output data set can then be provided to determine one or more physical properties of a portion of the volume of the Earth for subsurface exploration. The provided data set can be used for subsurface exploration.

In other words, the output from the method is a stack of seismic wavelets which are decomposed from the original seismic data set. These decomposed wavelets may be classified into various groups based on their characteristics. Each of characteristics or selective combination of characteristics reveal the physical properties of the subsurface of the Earth and can be used for subsurface exploration and identification of subsurface features.

For example, with reference to Figure 1, when seismic waves propagate through the anomaly 32 having anomalous physical properties, the decomposed wavelets will have distinguishing characteristics specific to the anomaly.

Based on the distinguished characteristics, the decomposed wavelets thus indicate the existence of the geological anomaly 32. This may be, for example, a petroleum reservoir or a specific type of physical property discontinuation of geophysical interest.

As a further example, the decomposition enables identification of features within a seismic trace. For example, it may be possible to identify a strong reflector (such as a coal seam) and remove this reflector from the data to reveal more subtle substructure having weaker reflections. This may enable identification of natural resources.

As a yet further example, the use of time-frequency spectrum data generated by the present invention may enable identification of signatures representative of subsurface structures such as gas reserves.

In other words, once the threshold criteria in step 207 is met, it is determined that the extracted data is deemed to be sufficiently accurate to be used for subsurface exploration. This may involve the identification of subsurface features such as cavities or channels which may contain natural resources such as hydrocarbons. Examples of such hydrocarbons are oil and natural gas.

Once these features have been identified in the output data, then measures can be taken to investigate these resources. For example, survey vessels or vehicles may be dispatched to drill pilot holes to determine whether natural resources are indeed present in these areas.

The above examples are not intended to be limiting and the skilled person would be readily aware of the suitable environments in which data could be gathered for analysis purposes as set out in the present disclosure.

In summary, in embodiments, the present invention aims to mitigate the effect of seismic data noise and any errors in estimating dip information. This leads to high to accuracy and robustness of multichannel seismic signal decomposition.

The present invention provides, in embodiments, a robust multichannel matching pursuit for seismic signal decomposition. As described, the seismic signal decomposition is performed as a multichannel matching pursuit. Any extracted basic wavelet is the optimal basic wavelet suitable to all channels within a group of multiple channels.

In embodiments, the multichannel matching pursuit is performed as a multichannel inverse problem. The wavelet amplitudes are calculated simultaneously for all wavelets and for all channels with a group of multiple channels.

In embodiments, the multichannel seismic signal decomposition is performed as a multichannel inverse problem with the soft constraint according to the determined dipping information. The extracted wavelets are optimal for seismic signals of multiple channels along the dipping orientation in a soft-constrained manner.

The effectiveness of the above-described embodiments is illustrated with reference to Figures 3 toe.

Figure 3a a set of seismic channels. Figure 3b shows the corresponding amplitude spectrum (at 15 Hz) generated by single channel matching pursuit. Finally, Figure 3c shows the corresponding amplitude spectrum (at 15 Hz) generated using the method of an embodiment of the present invention utilising the multi-channel matching pursuit method with a dipping constraint.

The benefits of the present invention are apparent in that the noisy spectrum data shown in Figure 3a is vastly reduced in the data of Figure 3b. This indicates the instability of known methods and the improvements of the present invention.

Figure 4 further illustrates the advantages of the present invention.

Figure 4a shows a group of seismic channels (left), the average channel without considering the dipping information (centre), and the average channel which takes to into account the dipping information (right).

Figure 4b shows the matching pursuit decomposition of the average channel without considering the dipping information.

Finally, Figure 4c shows the matching pursuit of the average channel which take account the dipping information. This clearly demonstrates the robustness of the structure-adaptive implementation.

Figure 5 illustrates a comparison of two convergence rates. The normalised residual energy versus iteration number is plotted for known MC-MP methods (circles) is plotted and for the methods of the present invention utilising structure adaptative implementation (triangles). This clearly demonstrates that the robust multichannel matching pursuit leads to significantly faster convergence than known methods. As a result, the approach of the present invention enables identification of subsurface features more quickly and utilising fewer computational resources.

Figure 6 shows a comparison of two residuals. Figure 6a shows the residual channels after a known MC-MP process after 100 iterations. Figure 6b shows the residual channels using the method of an embodiment of the present invention utilising the structure-adaptative approach after 20 iterations. This comparison clearly demonstrates that the robust multichannel matching pursuit leads to much higher accuracy than the standard method, because the implementation mitigates the effect of seismic data noise and any errors in estimating dip information.

In aspects, the embodiments described herein relate to a method of subsurface exploration. However, the embodiments described herein are equally applicable as an instruction set for a computer for carrying out said method or as a suitably programmed computer.

The methods described herein are, in use, executed on a suitable computer system or device running one or more computer programs formed in software and/or hardware and operable to execute the above method. A suitable computer system will generally comprise hardware and an operating system.

The term 'computer program' is taken to mean any of (but not necessarily limited to) an application program, middleware, an operating system, firmware or device drivers or any other medium supporting executable program code.

The term 'hardware' may be taken to mean any one or more of the collection of physical elements that constitutes a computer system/device such as, but not limited to, a processor, memory device, communication ports, input/output devices. The term 'firmware' may be taken to mean any persistent memory and the program code/data stored within it, such as but not limited to, an embedded system. The term 'operating system' may taken to mean the one or more pieces, often a collection, of software that manages computer hardware and provides common services for computer programs.

The methods described herein may be embodied in one or more pieces of software and/or hardware. The software is preferably held or otherwise encoded upon a memory device such as, but not limited to, any one or more of, a hard disk drive, RAM, ROM, solid state memory or other suitable memory device or component configured to software. The methods may be realised by executing/running the software. Additionally or alternatively, the methods may be hardware encoded.

The method encoded in software or hardware is preferably executed using one or more processors. The memory and/or hardware and/or processors are preferably comprised as, at least part of, one or more servers and/or other suitable computing systems.

Embodiments of the present invention have been described with particular reference to the examples illustrated. While specific examples are shown in the drawings and are herein described in detail, it should be understood, however, that the drawings and detailed description are not intended to limit the invention to the particular form disclosed. It will be appreciated that variations and modifications may be made to the examples described within the scope of the present invention.

Claims

CLAIMS1. A method of subsurface exploration, the method comprising determining one or more physical properties of a portion of the volume of the Earth from a seismic measurement representative of at least one physical parameter, the method comprising: a) providing an observed multichannel seismic dataset derived from a seismic measurement of a portion of the volume of the Earth; b) estimating local dipping information of the portion of the volume of the Earth to from the observed multichannel seismic dataset; c) defining the observed multichannel seismic dataset as a multichannel residual dataset; d) iteratively performing seismic signal decomposition on each channel of the multichannel residual dataset simultaneously to generate an output seismic dataset comprising a series of decomposed wavelets; and e) providing the output seismic dataset to determine one or more physical properties of a portion of the volume of the Earth for subsurface exploration.
2. The method according to claim 1, wherein step d) comprises: f) extracting a basic wavelet from the multichannel residual dataset; g) expanding a basic-wavelet matrix; h) calculating an amplitude of each wavelet for each of the multiple channels simultaneously to generate a synthetic multichannel seismic dataset using a dipping constraint derived from said estimated local dipping information; i) forming an updated multichannel residual dataset by subtracting the synthetic multichannel seismic dataset from the observed multichannel seismic dataset; j) repeating steps f) to i) using the updated multichannel residual dataset until a criterion of a threshold parameter is met; and, if said criterion is met; and k) generating an output seismic dataset comprising a series of decomposed wavelets.
3. The method according to claim 2, wherein step h) further comprises utilising an objective function having said dipping constraint.
4. The method according to claim 3, wherein the dipping constraint is operable to constrain solutions of the objective function towards extracted wavelets which are optimal for seismic signals of multiple channels along a dipping orientation derived from said estimated dipping information.
5. The method according to claim 4, wherein the dipping constraint is a soft dipping constraint.
6. The method according to claim 3, 4 or 5, wherein the objective function is defined in terms of a Frobenius Norm.
7. The method of claim 6, wherein the objective function is defined as 111Q-DK41-KAA9-11:.] where Mis the number of data channels, s"" is the multichannel seismic dataset, and A,c," is an amplitude matrix in which each column vector is composed of the amplitudes of a single channel.
8. The method according to any one of claims 3 to 7, wherein step h) further comprises minimizing the objective function.
9. The method according to claim 8, wherein minimising the objective function comprises calculate all wavelet amplitudes in a least-squares approach.
10. The method according to claim 9, wherein the step of minimising estimates the current wavelet amplitude and refines the amplitudes of wavelets extracted in all previous iterations.
11 The method according to claim 10, wherein the dipping constraint comprises one or more local dipping operators applied to the reconstruction of multichannel synthetic signals in the form of s-,A), =4,
12. The method according to any one of claims 5 to 11, wherein the degree of constraint of the dipping constraint is variable.
13. The method according to any one of claims 5 to 12, wherein the soft dipping constraint comprises first order differential operators D and D along the time axis and the space axis respectively, defined as: D= -1 1 where Q nir is a tridiagonal matrix with the sine elements, Q rar = ciiag{sinO, sinO, * * * sine"_,}, Q is a tridiagonal matrix with cosine elements, Q diagloos(9, cost92 * ** cnset",} , and 0, are the dipping angles along the time taken from a central channel of the multiple channels.to
14. The method according to any one of claims 2 to 13, wherein step f) further comprises extracting the basic wavelet from the defined multichannel residual dataset utilising structure-adaptive multichannel matching pursuit methods.
15. The method according to claim 14, wherein the structure-adaptive multichannel matching pursuit comprises the steps of: I) defining the Rh basic wavelet by one or more of the parameters of: mean frequency, wavelet breadth, phase, and time delay; m) setting an initial value of the or each parameters; and n) generating a structure-adaptive value of the multichannel residual dataset comprising the sum along a local dip of a coherent event.
16. The method according to claim 15, wherein step m) is set according to the parameters of an average channel from the multichannel residual dataset and step n) generates a structure-adaptive value.
17. The method according to any one of claims 2 to 16, wherein step g) comprises: o) adding an extra column vector, formed by the basic wavelet extracted at the current iteration of the robust multichannel matching pursuit for seismic signal decomposition.
18. The method according to claim 17, wherein in step o) the extra column vector is of the form (D,"x, = ow"11g, ,where Nis the number of samples of each data channel, <DAT_ " is the basic-wavelet matrix obtained in the previous iteration, g" is the vector of the current basic wavelet, and 11 is the augment operator which combines the basic-wavelet matrix ck"" and the vector g, .
19. The method according to any one of claims 2 to 18, wherein step i) comprises subtracting the synthetic multichannel dataset from the observed multichannel seismic dataset, R,, =s,, , where Itit, is the multichannel residual dataset after the km iteration.
20. The method according to any one of claims 2 to 19, wherein the threshold parameter comprises the relative energy of the multichannel residual dataset and the respective criterion is that the relative energy is lower than a pre-set percentage threshold.
21 The method according to any one of the preceding claims, wherein step b) comprises defining the local wavefield as local plane wave and utilising plane wave destruction methods to estimate the local dipping information.
22. The method of any one of the preceding claims, wherein the method is performed on data in a plurality of windows taken along a length direction.
23. The method of any one of the preceding claims, wherein step e) further comprises utilising the output seismic dataset for subsurface exploration.
24. A computer system comprising a processing device configured to perform the method of any one of the preceding claims.
25. A computer readable medium comprising instructions configured when executed to perform the method of any one of claims 1 to 23.
26. A computer system comprising: a processing device, a storage device and a computer readable medium according to claim 25.