Method of obtaining reflection properties of a subsurface region

Abstract

There is provided a computer-implemented method 70 of obtaining reflection properties of a subsurface region. The method 70 comprises determining 72 an initial dictionary constructed from a starting pulse. The method 70 further comprises inverting 74 seismic data of a subsurface region using the initial dictionary to obtain a model of subsurface reflectivities, and the method further comprises updating 76 the initial dictionary based on the model of subsurface reflectivities according to set criteria to generate a learned dictionary; and repeating the inverting 76 of seismic data using the learned dictionary until a convergence criteria is met.

Description

[0001] Method of Obtaining Reflection Properties of a Subsurface Region

[0002] Field

[0003] The subject disclosure generally relates to geophysics, and in particular to a method of obtaining reflection properties of a subsurface region.

[0004] Background

[0005] One of the goals of seismic interpretation is to relate seismic amplitude variations to rock and fluid properties from variations in acoustic impedance. However, when geological layers are thinner than the seismic wavelength, interference from the top and bottom (base) of a layer distorts the relationship between seismic amplitude and acoustic impedance.

[0006] For example, during CCS operations injected CO2 moves away from the well, typically forming thin layers in a subsurface formation. Efficient management of this process requires being able to accurately map the spatial distribution of the fluids. Seismic methods are commonly employed for this purpose, but a fundamental problem is that the CO2 layer is typically beneath seismic resolution; we cannot see a separation between the top and bottom of the layer. To overcome this limitation, geophysicists typically interpret the interference patterns between the top and bottom reflections - referred to as “tuning analysis”.

[0007] Tuning analysis assumes that the top and bottom reflections have the same waveforms, which is likely to be the case when there is no attenuation of the waves. However, CO2 injection may give rise to significant attenuation. This attenuation has the potential to change the shape of the reflected waveforms, which in turn leads to errors in the tuning analysis and potentially incorrect estimations of the layer thickness. Attempts to apply tuning analysis to CO2 injection projects have led to thickness estimates which do not match observations from the well. It is not possible to correct these estimates without knowledge of the attenuation coefficient, which is not usually available.

[0008] Accordingly, it is an object of this disclosure to at least partially address these issues and / or more accurately estimate subsurface structure. This background serves only to set a scene to allow a person skilled in the art to better appreciate the following description. Therefore, none of the above discussion should necessarily be taken as an acknowledgement that that discussion is part of the state of the art or is common general knowledge. One or more aspects / embodiments of the disclosure may or may not address one or more of the background issues.

[0009] Summary

[0010] The subject disclosure proposes a method of obtaining reflection properties of a subsurface region. The method comprises determining an initial dictionary constructed from a starting pulse. Seismic data of a subsurface region is inverted using this initial dictionary, providing a model of subsurface reflectivities. The reflectivities are then used to update the dictionary according to a set criteria. The seismic data is inverted again with the new dictionary. The process is repeated until a convergence criterion is met.

[0011] The outcome of the method is a set of spatially varying reflection coefficients and wavelet phase factors. The method accounts for the impact of attenuation on wavelet phase. Thus, the positions of the reflectivities more closely match the physical positions of subsurface geological layers in the region as convergence is approached. This may facilitate more accurate layer thickness determination. Additionally, the wavelet phase factors may be transformed into subsurface attenuation measurements.

[0012] According to an aspect of the disclosure there is provided a computer-implemented method of obtaining reflection properties of a subsurface region, the method comprising: determining an initial dictionary constructed from a starting pulse; inverting seismic data of a subsurface region using the initial dictionary to obtain a model of subsurface reflectivities; updating the initial dictionary based on the model of subsurface reflectivities according to set criteria to generate a learned dictionary; and repeating the inverting of seismic data using the learned dictionary until a convergence criteria is met.

[0013] The starting pulse may define a starting pulse shape. The starting pulse may comprise a starting wavelet. The starting wavelet may be a zero-phase wavelet. The starting wavelet may comprise a Ricker Wavelet. The Ricker Wavelet may be a zero-phase Ricker Wavelet. The Ricker Wavelet may be a 35 Hz zero-phase Ricker Wavelet.

[0014] The subsurface reflectivities may comprise a reflectivity series, i.e., a series of reflectivities. The reflectivities may be generated or reflected from the subsurface region.

[0015] The method may further comprise determining layer thickness of the subsurface region based on the model of subsurface reflectivities.

[0016] The method may further comprise mapping a distribution of fluids injected into or preexisting in the subsurface region based on the determined layer thickness.

[0017] The seismic data may be associated with reflections from top and bottom boundary layers of a seismic layer.

[0018] The seismic layer may be formed by CO2 injected into a wellbore. The seismic data may be obtained after CO2 is injected into a wellbore. That is to say, the seismic data may represent seismic data of a subsurface region into which CO2 has been injected.

[0019] Inverting may comprise basis pursuit (BP) and / or sparse spike inversion (SSI).

[0020] Inverting may comprise minimising the following expression to determine the model: where r is a reflectivity series, where s is a seismic trace, where Do is the initial dictionary, and where Acis a regularisation parameter.

[0021] The initial dictionary may be constructed from standard pulse shapes, or from wavelets extracted from specific datasets. For example, the standard pulse may comprise a Ricker Wavelet.

[0022] The initial dictionary comprises a matrix. The learned dictionary comprises a matrix. Different values of the regularisation parameters may be used for different updates. The regularisation parameter may comprise a first parameter for coarse sparse coding and a second parameter for fine sparse coding. The method may comprise coarse sparse coding using a first parameter. The method may further comprise fine sparse coding using a second parameter.

[0023] The first parameter may be approximately ten times larger than the second parameter.

[0024] Optimising may comprise iteratively updating the dictionary to generate the learned dictionary.

[0025] Iteratively updating may comprise modifying non-zero elements of the dictionary to generate the learned dictionary.

[0026] Modifying non-zero elements may comprise filtering elements of the dictionary to generate the learned dictionary. The filtering may comprise introducing a constant phase shift using a Hilbert transform. Elements may be phase shifted through a phase range at a pre-set interval. The interval may comprise -60 to +60 degrees, -50 to +50 degrees or -100 to +100 degrees, or any other suitable range. The interval may comprise 20, 10, 5 or 2 degrees, or any other suitable interval.

[0027] The convergence criterion may comprise a minimum sparsity of the learned dictionary.

[0028] The minimum sparsity may correspond to an optimal phase shift of the recovered reflectivity.

[0029] Updating the dictionary may comprise convolving at least one element of the dictionary with a phase shifted wavelet. Reconstructing the elements may comprise reconstructing the i-th element of the learned dictionary according to the following formula: where DL is the learned dictionary, where Itis the i-th atom or element of the identity matrix of ]RWxWwhere N denotes the length of the seismic data, and where a>e. is a wavelet filtered by a suitable operator that mimics the effect of dispersion and attenuation. For example, filtering can be by introducing a constant phase rotation by an amount of using a Hilbert transform.

[0030] Updating may comprise calculating a sparsity of a trace reconstructed using the learned dictionary. Optimising may comprise iteratively reconstructing and calculating the sparsity until a desired sparsity is calculated. The desired sparsity may correspond to a convergence criterion. Calculating may comprise calculating the sparsity according to following formula: where r is the reflectivity, where s is the seismic data, where DL is the learned dictionary, where the formula is updated for every j-th version of the dictionary, and where Af is the fine regularisation parameter.

[0031] The method may further comprise obtaining an optimal phase from the learned dictionary.

[0032] The method may further comprise estimating attenuation based on the model of subsurface reflectivities. The method may further comprise estimating inter-layer attenuation of the subsurface region based on the optimal phase.

[0033] The method may further comprise controlling injection of CO2 into the wellbore according to determined parameters. The determined parameters may comprise at least one of the recovered reflectivity, the determined layer thickness, the obtained optimal phase, and the estimated inter-layer attenuation. The parameters may be more accurately determined than conventional methods. As such, controlling injection may be more accurate than conventional methods. This may allow for greater levels of CO2 to be injected. Additionally, this may improve the safety and / or reduce the environmental impact of the injection. The method may further comprise detecting seismic data of the subsurface region using a sensor array. The method may further comprise receiving seismic data of the subsurface region using a sensor array.

[0034] According to another aspect there is provided a non-transitory computer-readable medium having computer-program code stored thereon, the code executable by a processor to perform the described method.

[0035] According to another aspect there is provided a computing system comprising: a processor; and a memory having computer-program code stored thereon, the code executable by the processor to perform the described method.

[0036] According to another aspect there is provided a controller adapted to perform the described method. The controller may comprise a processor and memory comprising computer-readable code executable by the processor to perform the described method. The controller may take the form of one or more microcontrollers (MCUs).

[0037] Further elements of the aspects described may include one or more examples, embodiments or features in isolation or in various combinations whether or not specifically stated (including claimed) in that combination or in isolation.

[0038] Brief Description of the Drawings

[0039] These and other aspects of the present disclosure will now be described, by way of example only, with reference to the accompanying Figures, in which:

[0040] Figures 1a and 1b are plots of a 35 Hz zero-phase Ricker Wavelet in the time and frequency domains, respectively;

[0041] Figure 2 is a plot of a dictionary of the Ricker Wavelet of Figures 1a and 1 b;

[0042] Figure 3 is a plot of phase shifted wavelets of the Ricker Wavelet of Figures 1 a and 1 b; Figures 4a1-d3 are plots of the Ricker Wavelet of Figures 1a and 1b, and SSI recovered traces with various regularisation parameters;

[0043] Figures 5a1-d3 are plots of a 30 degree phase shifted Ricker Wavelet and SSI recovered seismic traces with various regularisation parameters; Figures 6a1-d3 are plots of a 50 degree phase shifted Ricker Wavelet and SSI recovered seismic traces with various regularisation parameters;

[0044] Figure 7 is a flowchart of a method of obtaining reflection properties according to an aspect of the disclosure;

[0045] Figures 8a1-d2 are plots of a 30 degree phase shifted Ricker Wavelet, SSI recovered seismic traces, and seismic traces recovered with a method according to an aspect of the disclosure;

[0046] Figures 9a1-d2 are plots of a 50 degree phase shifted Ricker Wavelet, SSI recovered seismic traces, and seismic traces recovered with a method according to an aspect of the disclosure;

[0047] Figures 10a1-d3 are plots of the inversion results of a method according to an aspect of the disclosure and conventional SSI for 0 degree phase shift at various noise levels;

[0048] Figures 11a1-d3 are plots of the inversion results of a method according to an aspect of the disclosure and conventional SSI for +30 degree phase shift at various noise levels; and

[0049] Figures 12a1-d3 are plots of the inversion results of a method according to an aspect of the disclosure and conventional SSI for -30 degree phase shift at various noise levels.

[0050] Detailed Description

[0051] The foregoing summary, as well as the following detailed description of certain embodiments will be better understood when read in conjunction with the accompanying drawings. As will be appreciated, like reference characters are used to refer to like elements throughout the description and drawings. As used herein, an element or feature recited in the singular and preceded by the word "a" or "an" should be understood as not necessarily excluding a plural of the elements or features. Further, references to "one example" or “one embodiment” are not intended to be interpreted as excluding the existence of additional examples or embodiments that also incorporate the recited elements or features of that one example or one embodiment. Moreover, unless explicitly stated to the contrary, examples or embodiments "comprising", "having" or “including” an element or feature or a plurality of elements or features having a particular property might further include additional elements or features not having that particular property. Also, it will be appreciated that the terms “comprises”, “has” and “includes” mean “including but not limited to” and the terms “comprising”, “having” and “including” have equivalent meanings. As used herein, the term “and / or” can include any and all combinations of one or more of the associated listed elements or features.

[0052] It will be understood that when an element or feature is referred to as being “on”, “attached” to, “connected” to, “coupled” with, “contacting”, etc. another element or feature, that element or feature can be directly on, attached to, connected to, coupled with or contacting the other element or feature or intervening elements may also be present. In contrast, when an element or feature is referred to as being, for example, “directly on”, “directly attached” to, “directly connected” to, “directly coupled” with or “directly contacting” another element of feature, there are no intervening elements or features present.

[0053] It will be understood that spatially relative terms, such as “under”, “below”, “lower”, “over”, “above”, “upper”, “front”, “back” and the like, may be used herein for ease of describing the relationship of an element or feature to another element or feature as depicted in the figures. The spatially relative terms can however, encompass different orientations in use or operation in addition to the orientation depicted in the figures.

[0054] Reference herein to “example” means that one or more feature, structure, element, component, characteristic and / or operational step described in connection with the example is included in at least one embodiment and or implementation of the subject matter according to the present disclosure. Thus, the phrases “an example,” “another example,” and similar language throughout the present disclosure may, but do not necessarily, refer to the same example. Further, the subject matter characterizing any one example may, but does not necessarily, include the subject matter characterizing any other example.

[0055] Reference herein to “configured” denotes an actual state of configuration that fundamentally ties the element or feature to the physical characteristics of the element or feature preceding the phrase “configured to”.

[0056] Unless otherwise indicated, the terms “first,” “second,” etc. are used herein merely as labels, and are not intended to impose ordinal, positional, or hierarchical requirements on the items to which these terms refer. Moreover, reference to a “second” item does not require or preclude the existence of lower-numbered item (e.g., a “first” item) and / or a higher-numbered item (e.g., a “third” item).

[0057] As used herein, the terms “approximately” and “about” represent an amount close to the stated amount that still performs the desired function or achieves the desired result. For example, the terms “approximately” and “about” may refer to an amount that is within less than 10% of, within less than 5% of, within less than 1% of, within less than 0.1% of, or within less than 0.01% of the stated amount.

[0058] In geophysics, an acoustic impulse generated by explosive dynamite, a seismic vibrator, airgun or acoustic source propagates through Earth materials. A portion of the energy is reflected back to surface at layer boundaries. The boundaries are the top and bottom boundaries between layers in the substructure. A sensor array at surface may observe the reflected signal, e.g., a seismogram, as a function of time. The basic goal in reflection seismology is to recover the layered structure of the subsurface as a function of depth from observed seismic data.

[0059] The sensor array observes seismic data referred as a seismic trace - a curve representing ground movement. According to the convolutional model, the seismic trace s[ may be modelled as the convolution of a reflectivity series r(t) and a wavelet w(t) plus additional noise n(t) as set out in Equation 1 below. s(t) = w(t) * r(t) + n(t) Equation 1

[0060] The operator * is the one-dimensional linear convolution operator. The convolutional model assumes the subsurface structure may be represented by planar horizontal layers of elastic constant impedance, so that reflections are generated at the interfaces (top / bottom of adjacent layers) between two adjacent layers where the acoustic impedances are different. After time discretisation, Equation 1 can be written in the matrix-vector from as Equation 2 below. NX1 ~ ^NxMrM+l +nNxl Equation 2 where WNxMe ]RWx ]RMEquation 3

[0061] This is known as the wavelet dictionary, i.e., the elements which compose the elements of the wavelet.

[0062] Seismic inversion is a process of transforming the seismic reflection data r(t), captured by the surface sensor array, into the layered structure of the subsurface. Seismic inversion generally requires determination of the reflectivity r(t), given a wavelet w(t) estimated from the captured seismic data s(t).

[0063] Several methods of seismic inversion exist such as sparse spike inversion (SSI) and basis pursuit inversion (BPI). SSI assumes the reflectivity series r(t) is sparse. Generally, a sparse matrix is a matrix in which the great majority of matrix elements are zero. SSI has been shown to be capable of handling random noise suppression, but has difficulty when dealing with phase distorted waveforms. An increase in phase distortion may result in an increase in error in the recovery of SSI parameters, e.g., regularisation parameter, and may distort the waveform.

[0064] In SSI, WNxMis the convolution matrix formed by the time-variant seismic wavelet w(t). This is often chosen to be constructed using a zero-phase wavelet given the advantages of zero-phase, e.g., wavelet symmetry, minimal ambiguity of wavelet shape in correlation, centre and maximum amplitude of the wavelet coinciding in time with reflection interface, and best resolution among wavelets with the same amplitude spectrum.

[0065] Based on the sparse assumption of the reflectivity series r(t), the inverse problem of finding rMxlfrom a noisy measurement of the seismic waveform SNxlmay be formulated as a minimisation problem of a suitably chosen objective function F which depends on the distance between the reconstructed and seismic traces as well as the sparseness of the reconstructed trace.

[0066] F may be reformatted as a basis pursuit, matching pursuit or other algorithm that combines a least square penalty as well as a sparseness penalty. In particular, F may be written as a basis pursuit function as per Equation 4 below. Equation 4

[0067] A regularisation parameter is introduced in Equation 4 to balance signal misfit, reflectivity sparsity and noise. Increasing the regularisation parameter may decrease the sparsity of the solution, while increasing the parameter may increase noise, i.e., noise amplification. The parameter may be determined empirically. To obtain the parameter empirically, the inversion is performed for different values of the regularisation parameter and the optimal parameter is selected based on the fewest number of spikes when convolving with the wavelet that can fit the data within a given tolerance. Hence the main goal of SSI is to find the number of spikes that, when convolved with the seismic wavelet, fit the data within a given tolerance.

[0068] Solving Equation 4 is thus essentially an optimisation problem which can be solved with various techniques such as Linear Program, Least Operator Shrinkage and Selector Operator (LASSO) or Convex Optimisation Technique (CVX). These are described in Chen, S. S., D. L. Donoho, and M. A. Saunders (2001), “Atomic decomposition by basis pursuit’, In: SIAM review 43.1, pp. 129-159; Zhang, R. (2010), “Seismic reflection inversion by basis pursuit’, University of Houston, isbn: 1124504680; Rozenberg, I., I. Cohen, and A. Vassiliou (2014), “Performance measures for sparse spike inversion vs. basis pursuit inversion", In: 2014 IEEE 28th Convention of Electrical & Electronics Engineers in Israel (IEEEI). IEEE, pp. 1-5. isbn: 147995988X; and Pereg, D., I. Cohen, and A. Vassiliou (2017), “Sparse seismic timevariant deconvolution using Q attenuation modef’, In: SEG Technical Program Expanded Abstracts 2017. Society of Exploration Geophysicists, pp. 508-513. isbn: 1949-4645, respectively.

[0069] However, SSI may not be suitable for all use cases. For example, CO2 injection in subsurface structures may result in the formation of thin layers. These thin layers may be interbed between adjacent elastic layers. Such thin interbed layers may have anelastic or dispersive properties. Conventional seismic inversion techniques assume top and bottom reflections from these layers have the same waveforms, which is likely to be the case when there is no attenuation of the reflections. The form of the reflected waves may not change with distance travelled if the layers are elastic. However, if the layers are anelastic and / or dispersive, the reflected wave may change with distance travelled, i.e., wave attenuation. In particular, due to attenuation, reflected seismic waves from interfaces between elastic and dispersive layers may be affected by amplitude and / or phase distortions.

[0070] Theory and observation suggest that CO2 injection gives rise to significant attenuation. This attenuation has the potential to change the shape of the reflected waveforms, which in turn leads to errors in seismic inversion and incorrect estimates of the layer thickness. As mentioned, SSI may not suitably account for such amplitude and / or phase attenuation.

[0071] Performance of SSI may be evaluated using an exemplary 35 Hz zero-phase Ricker Wavelet which is illustrated in Figure 1. A dictionary of this wavelet is constructed and illustrated in Figure 2. Attenuation is simulated on the Ricker Wavelet from +100 to - 100 degrees as illustrated in Figure 3.

[0072] As previously described the regularisation parameter S must be chosen. The performance of the inversion for various values of the parameter are illustrated in Figures 4-6.

[0073] Figures 4a1-d3 are plots of the Ricker Wavelet and SSI recovered traces with various regularisation parameters. In particular, the synthetic data (i.e., Ricker Wavelet) is identified as the Model curve, while two other curves with various regularisation parameters are illustrated, namely regularisation parameter s = 0.05 and regularisation parameter ^ = 0.01. The curves are illustrated for signal to noise ratios (SNR) of: (a1) SNR = 100dB, (b1) SNR = 20dB, (c1) SNR = 10dB and (d1) SNR = 5dB. True reflectivity and inverted reflectivity is illustrated with the regularisation parameter S = 0.05 in (a2), (b2), (c2), and (d2). True reflectivity and inverted reflectivity is illustrated with the regularisation parameters = 0.01 in (a3), (b3), (c3), and (d3).

[0074] Figures 5a1-d3 are plots of a +30 degree phase shifted Ricker Wavelet and SSI recovered seismic traces with the same regularisation parameters as Figures 4a1-d3. In particular, the synthetic data (i.e., Ricker Wavelet phase shifted by +30 degrees) is identified as the Model curve, while two other curves with various regularisation parameters are illustrated, namely regularisation parameter s = 0.05 and regularisation parameters = 0.01. The curves are illustrated for SNR of: (a1) SNR = 100dB, (b1) SNR = 20dB, (c1) SNR = 10dB and (d1) SNR = 5dB. True reflectivity and inverted reflectivity with A = 0.05 in (a2), (b2), (c2) and (d2). True reflectivity and inverted reflectivity with A = 0.01 in (a3), (b3), (c3) and (d3).

[0075] Figures 6a1-d3 are plots of a +50 degree phase shifted Ricker Wavelet and SSI recovered seismic traces with the same regularisation parameters as Figures 4a1-d3. In particular, the synthetic data (i.e., Ricker Wavelet phase shifted by +50 degrees) is identified as the Model curve, while two other curves with various regularisation parameters are illustrated, namely regularisation parameter A = 0.05 and regularisation parameters = 0.01. The curves are illustrated for SNR of: (a1) SNR = 100dB, (b1) SNR = 20dB, (c1) SNR = 10dB, (d1) SNR = 5dB. True reflectivity and inverted reflectivity with S = 0.05 in (a2), (b2), (c2), and (d2). True reflectivity and inverted reflectivity with S = 0.01 in (a3), (b3), (c3), and (d3).

[0076] In general, the inversion quality decreases as noise level increases. Traces recovered with smaller S fit to the input data better than those with larger S. But smaller A would also recover more noisy spikes when SNR decreases. In the zero-phase case (Figures 4a1-d3_, the location of spikes is correct, although the amplitude of inverted reflectivity is underestimated which is a common problem with tv-norm regularization.

[0077] For the phased wavelets (Figures 5 and 6), the location of the spikes is recovered several samples away from their true locations and with an underestimation of reflectivity amplitude. The larger phase the wavelet carries (i.e., the larger the phase attenuation), the more samples the spikes are shifted. This may be due to the peak phase wavelets being shifted several samples from the zero centre. Since the actual waveform is distorted due to dispersion (which causes the attenuation), SSI is unable to recover the actual spike positions. This illustrates the limitations of SSI when dealing with attenuation which may include phase and / or amplitude shifts.

[0078] The subject disclosure proposes a method of at least partially addressing these limitations. The method comprises determining an initial dictionary constructed from a starting pulse, then inverting seismic data of a subsurface region using the initial dictionary to obtain a model of subsurface reflectivities. The initial dictionary is updated based on the model of subsurface reflectivities according to set criteria to generate a learned dictionary. The inverting is repeated using the learned dictionary until a convergence criterion is met. The resulting output of the method is a more accurate determination of reflections (reflectivities) than conventional methods. The reflections may be used to provide a more accurate estimate of interface positions. The interface positions may be used to provide a more accurate estimate of layer thickness. This may be particularly beneficial when attenuation is significant such as for example in thin layers of CO2 injected into a wellbore. Additionally, the described method recovers phase information which provides an indication of inter-layer attenuation. This may be useful in determining additional geophysical properties of the subsurface structure.

[0079] Turning now to Figure 7, a flowchart of the method 70 of obtaining reflection properties is illustrated. The method 70 comprises determining 72 an initial dictionary constructed from a starting pulse. The seismic data is collected or received from a sensor array located at a surface or subsea location of the surface region. The seismic data is associated with reflections from top and bottom boundary layers of a seismic layer. When used in carbon capture and storage (CCS), the seismic layer may be formed by CO2 injected into a wellbore. The wellbore may be an abandoned oil and gas wellbore or a wellbore specifically drilled for CCS.

[0080] The method 70 further comprises inverting 74 seismic data of a subsurface region using the initial dictionary to obtain a model of subsurface reflectivities. The inverting 74 comprises solving the Equation 4 with a coarse regularisation parameter. The resulting formula, Equation 5, is reproduced below. Equation 5 where r is a reflectivity series, where s is a seismic trace, where Do is the initial dictionary, and where Acis a coarse regularisation parameter.

[0081] The coarse regularisation parameter Acmay be determined via trial and error. As such, the method 70 may comprise determining the coarse regularisation parameter of Equation 5. Equation 5 is solved to obtain a reflectivity model m. The set of non-zero elements of the matrix of the dictionary is given by Equation 6 which written below.

[0082] B = {j|mj #= 0} Equation 6

[0083] The set is ranked by its absolute amplitude in a descent sequence. The phase is then determined by modifying the corresponding wavelets in the specific dictionary elements indicated in the set B. This is performed in the following updating 76 step of the method 70.

[0084] As such, the method 70 further comprises updating 76 the initial dictionary based on the model of subsurface reflectivities according to set criteria to generate a learned dictionary. The updating 76 is performed based on the fact that the correct wavelet (i.e., the wavelet having the correct phase) would approximate the reflection with fewer spikes in SSI. The optimal phase may be determined using a grid search over a predefined phase range. The phase range may be from -60 degrees to +60 degrees in 10 degree intervals. This range is only exemplary and other ranges may be used. For examples, the range may be from -50 degrees to +50 degrees and the intervals may be 5 degrees. As one of skill in the art will appreciate other ranges may be used.

[0085] The updating 76 comprises iteratively updating the dictionary by modifying the nonzero elements of the dictionary. Specifically, by phase shifting the non-zero elements, i.e., wavelets, of the dictionary. For every non-zero reflector, the corresponding wavelet in the dictionary is updated with every value of phase in the given range. The ith element or atom of the dictionary is then reconstructed according to Equation 7 written below.

[0086] BLl= It * o)e. Equation 7 where DL is the learned dictionary, where It is the i-th atom or element of the identity matrix of ]RWxW, where N denotes the length of the seismic data, and where o)0is a wavelet filtered by a suitable operator that mimics the effect of dispersion and attenuation. For example, filtering can be by introducing a constant phase rotation by an amount of using a Hilbert transform. The sparsity of the now learned dictionary is calculated according to Equation 8 written below. Equation 8 where r is the reflectivity, where s is the seismic data, where DL is the learned dictionary, where the formula is updated for every j-th version of the dictionary, and where Af is the fine regularisation parameter.

[0087] The fine regularisation parameter Af may be determined via trial and error. As such, the method 70 may comprise determining the fine regularisation parameter of Equation 8. It has been found that the fine regularisation parameter is approximately ten times smaller than the coarse regularisation parameter. The exact parameters may depend on the amplitude of the seismic data and noise levels.

[0088] As shown in Figure 7, the method 70 comprises iteratively repeating the inverting 74 and updating 76 until a convergence criterion is reached. In other words, the wavelet is phase shifted for every phase in the range according to the interval until the sparsity according to Equation 8 is minimised.

[0089] The final output of the method 70 will take the form of Equation 9 written below.

[0090] From this equation, several entities may be determined. In particular, the reflectivity; the optimal phase, as per the minimised sparsity; and the learned dictionary. The learned dictionary includes the wavelet with the optimal phase shift selected on a minimum sparsity basis. The optimal phase is given by Equation 10 written below and the learned dictionary is given by Equation 11 written below. Equation 10

[0091] DLl= It * wg. Equation 11

[0092] The sparsity calculation (Equation 10) of the initial dictionary is performed with a coarse regularisation parameter Ac. Then, around each major reflection, i.e. , each jth version of the updated dictionary (learned dictionary), the sparsity is calculated with a smaller fine regularisation parameter Af.

[0093] The described method may be referred to as phased wavelet inversion (PWI). Performance of the PWI compared to conventional SSI is illustrated in Figures 8-12. In particular, Figures 8 and 9 illustrate phase shifted Ricker Wavelet, SSI recovered seismic traces, and PWI recovered seismic traces. In Figures 8a1-d2, the phase shift is +30 degrees. In Figures 9a1-d2, the phase shift is +50 degrees. The zero phase is not shown as this would be the same as illustrated in the plots of Figure 4.

[0094] The curves in Figures 8a1-d2 are illustrated for SNRs of: (a1) SNR = 100dB, (b1) SNR = 20dB, (c1) SNR = 10dB, and (d1) SNR = 5dB. Inverted reflectivities are illustrated for various regularisation parameters, in particular: (a2) A = 0.01 , (b2) A = 0.01 , (c2) A = 0.03, (d2) A = 0.03. In these plots, the synthetic data (i.e., ground truth model of the Ricker Wavelet) curves are identified as “Model”. SSI sparse coding, i.e., coding with a zero phase wavelet dictionary, curves are identified as “Do”. Sparse coding according to PWI, i.e., with a phase wavelet dictionary, curves are identified as “Dp”.

[0095] The curves in Figures 9a1-d2 are illustrated for SNRs of: (a1) SNR = 100dB, (b1) SNR = 20dB, (c1) SNR = 10dB, and (d1) SNR = 5dB. Inverted reflectivities are illustrated for various regularisation parameters, in particular: (a2) A = 0.01 , (b2) A = 0.01 , (c2) A = 0.05, (d2) A = 0.06. Similar to Figure 8, in these plots, the synthetic data (i.e., ground truth model of the Ricker Wavelet) curves are identified as “Model”. SSI sparse coding, i.e., coding with a zero phase wavelet dictionary, curves are identified as “Do”. Sparse coding according to PWI, i.e., with a phase wavelet dictionary, curves are identified as “D Up ” . Figures 8 and 9 illustrates that the spike locations of inverted reflectivity with zero phase wavelet dictionary (Do), i.e., conventional SSI, are shifted several samples away from ground truth (Model) no matter what noise level. However, in the dictionary according to the PWI (Dp), the location of inverted reflectivity is accurate, although the amplitude is still underestimated. Further, there are no small spikes around major reflections which are attempting to compensate for the difference in the waveform due to the phase change, as in the case with the zero phase wavelet dictionary (Do). This illustrates the significance of accurate wavelet (phase) in seismic reflectivity inversion especially for temporal layer thickness estimation using reflectivity arrival time.

[0096] Turning now to Figures 10-12, plots of the inversion results of PWI and conventional SSI are illustrated for various phase shifts at various noise levels. In particular, Figures 10a1-d3 are plots of PWI and SSI for a 0 degree phase shift, Figures 11a-d3 are plots of PWI and SSI for a +30 degree phase shift, and Figures 12a-d3 are plots of PWI and SSI for a -30 degree phase shift. In these Figures, the black curves denote ground truth model. The dotted curves denote signals recovered by SSI and dashed curves denote signals recovered by PWI, i.e., according to the described method.

[0097] Figures 10a-d3 illustrate the inversion results of zero-phase thin-bed (thickness = 10 ms) with top reflectivity (0 degree phase) being 0.1 and bottom reflectivity (0 degree phase) being -0.1. In (a1) and (a2), the ground truth model (synthetic data) has an SNR = 100 dB. In (a1), the inverted reflectivity by SSI has a regularisation parameter of A = 0.05. In (a2), the inverted reflectivity by PWI has regularisation parameters of Ac= 0.05, and Af = 0.003. In (b1) and (b2), the ground truth model (synthetic data) has an SNR = 20 dB. In (b1), the inverted reflectivity by SSI has a regularisation parameter of A = 0.05. In (b2), the inverted reflectivity by PWI has regularisation parameters of Ac= 0.05, and Af = 0.003. In (c1) and (c2), the ground truth model (synthetic data) has an SNR = 10 dB. In (c1), the inverted reflectivity by SSI has a regularisation parameter of A = 0.05. In (c2), the inverted reflectivity by PWI has regularisation parameters of Ac= 0.05, and Af = 0.003. In (d1) and (d2), the ground truth model (synthetic data) has an SNR = 5 dB. In (d1), the inverted reflectivity by SSI has a regularisation parameter of A = 0.06. In (d2), the inverted reflectivity by PWI has regularisation parameters of Ac= 0.06, and Af = 0.003. The curves in (a3),(b3),(c3) and (d3) illustrate show the phase model (dotted lines) and PWI estimated phases (dashed lines). Figures 11a-d3 illustrate the inversion results of zero-phase thin-bed (thickness = 10 ms) with top reflectivity (zero phase) being 0.1 and bottom reflectivity (+30 degree phase shifted) being -0.1. In (a1) and (a2), the ground truth model (synthetic data) has an SNR = 100 dB. In (a1), the inverted reflectivity by SSI has a regularisation parameter of A = 0.05. In (a2), the inverted reflectivity by PWI has regularisation parameters of Ac= 0.05, and Af = 0.003. In (b1) and (b2), the ground truth model (synthetic data) has an SNR = 20 dB. In (b1), the inverted reflectivity by SSI has a regularisation parameter of A = 0.06. In (b2), the inverted reflectivity by PWI has regularisation parameters of Ac= 0.06, and Af = 0.003. In (c1) and (c2), the ground truth model (synthetic data) has an SNR = 10 dB. In (c1), the inverted reflectivity by SSI has a regularisation parameter of A = 0.06. In (c2), the inverted reflectivity by PWI has regularisation parameters of Ac= 0.06, and Af = 0.003. In (d1) and (d2), the ground truth model (synthetic data) has an SNR = 5 dB. In (d 1 ), the inverted reflectivity by SSI has a regularisation parameter of A = 0.06. In (d2), the inverted reflectivity by PWI has regularisation parameters of Ac= 0.06, and Af = 0.003. The curves in (a3),(b3),(c3) and (d3) illustrate show the phase model (dotted lines) and PWI estimated phases (dashed lines).

[0098] Figures 12a-d3 illustrate the inversion results of zero-phase thin-bed (thickness = 10 ms) with top reflectivity (zero phase) being 0.1 and bottom reflectivity (-30 degree phase shifted) being -0.1. In (a1) and (a2), the ground truth model (synthetic data) has an SNR = 100 dB. In (a1), the inverted reflectivity by SSI has a regularisation parameter of A = 0.06. In (a2), the inverted reflectivity by PWI has regularisation parameters of Ac= 0.06, and Af = 0.003. In (b1) and (b2), the ground truth model (synthetic data) has an SNR = 20 dB. In (b1), the inverted reflectivity by SSI has a regularisation parameter of A = 0.06. In (b2), the inverted reflectivity by PWI has regularisation parameters of Ac= 0.06, and Af = 0.003. In (c1) and (c2), the ground truth model (synthetic data) has an SNR = 10 dB. In (c1), the inverted reflectivity by SSI has a regularisation parameter of A = 0.06. In (c2), the inverted reflectivity by PWI has regularisation parameters of Ac= 0.06, and Af = 0.003. In (d1) and (d2), the ground truth model (synthetic data) has an SNR = 5 dB. In (d 1 ), the inverted reflectivity by SSI has a regularisation parameter of A = 0.06. In (d2), the inverted reflectivity by PWI has regularisation parameters of Ac= 0.06, and Af = 0.003. The curves in (a3),(b3),(c3) and (d3) illustrate show the phase model (dotted lines) and PWI estimated phases (dashed lines).

[0099] As illustrated in Figures 10-12, the inverted reflectivities by PWI are of higher accuracy than those estimated by SSI, especially in terms of locations of spikes. Even in the zero phase situation shown in Figure 10, we still see the improvement of recovered spike locations from PWI. In phased examples (Figures 11 and 12), the wavelet ’peak’ for the bottom reflector arrives several samples later than the reflectivity spike in Figure 11, but it arrives earlier in Figure 12. Both examples show that PWI performs better in estimation locations of reflectivity. This is due to the updating of the dictionary which optimises the phase of the wavelets to correct for waveform distortions. The recovered phases fit with the ground truth with an acceptable error in Figures 11 and 12. While Figures 11 d1-d3 and 12d1-d3 illustrate that higher noise levels decrease the quality of phase estimation. However, the method further include suppressing random and linear noise before processing to reduce noise to acceptable levels.

[0100] These examples illustrate the improved results when using PWI compared to conventional inversion methods.

[0101] It should be understood that the examples provided are merely exemplary of the present disclosure, and that various modifications may be made thereto.

Claims

CLAIMS:

1. A computer-implemented method of obtaining reflection properties of a subsurface region, the method comprising: determining an initial dictionary constructed from a starting pulse; inverting seismic data of a subsurface region using the initial dictionary to obtain a model of subsurface reflectivities; updating the initial dictionary based on the model of subsurface reflectivities according to set criteria to generate a learned dictionary; and repeating the inverting of seismic data using the learned dictionary until a convergence criterion is met.

2. The method of claim 1, further comprising determining layer thickness of the subsurface region based on the model of subsurface reflectivities.

3. The method of claim 2, further comprising mapping a distribution of fluids injected into or pre-existing in the subsurface region based on the determined layer thickness.

4. The method of any preceding claim, wherein the seismic data is associated with reflections from top and bottom boundary layers of a seismic layer.

5. The method of claim 4, wherein the seismic layer is formed by CO2 injected into a wellbore.

6. The method of any preceding claim, wherein inverting comprises basis pursuit (BP) and / or sparse spike inversion (SSI).

7. The method of any preceding claim, wherein inverting comprises minimising the following expression to determine the model:where r is a reflectivity series, where s is a seismic trace, where Do is the initial dictionary, andwhere Acis a regularisation parameter.

8. The method of claim 7, wherein the regularisation parameter comprises a first parameter for coarse sparse coding and a second parameter for fine sparse coding.

9. The method of claim 8, further comprising coarse sparse coding using the first parameter and fine sparse coding using the second parameter.

10. The method of claim 8 or 9, wherein the first parameter is approximately ten times larger than the second parameter.

11. The method of any preceding claim, wherein updating comprises iteratively updating the dictionary to generate the learned dictionary.

12. The method of claim 11 , wherein iteratively updating comprises modifying nonzero elements of the dictionary to generate the learned dictionary.

13. The method of claim 12, wherein modifying non-zero elements comprises phase shifting elements of the dictionary to generate the learned dictionary.

14. The method of any preceding claim, wherein the convergence criterion comprises a minimum sparsity of the learned dictionary.

15. The method of claim 14, wherein the minimum sparsity corresponds to an optimal phase shift of the subsurface reflectivities.

16. The method of any preceding claim, wherein updating comprises reconstructing elements of the learned dictionary by convolving at least one element of the learned dictionary with a phase shifted wavelet.

17. The method of claim 16, wherein reconstructing comprises reconstructing the i- the element of the learned dictionary according to the following formula:where DL is the learned dictionary,where Itis the i-th atom or element of the identity matrix of ]RWxW, where N denotes the length of the seismic data, and where o)0is a filtered wavelet18. The method of claim 17, wherein o)0is filtered by introducing a constant phase rotation by an amount of using a Hilbert transformation.

19. The method of any preceding clam, wherein updating comprises calculating a sparsity of a trace reconstructed using the learned dictionary.

20. The method of claim 19, wherein calculating the sparsity comprises calculating the sparsity according to the following formula:where r is the reflectivity, where s is the seismic data, where DL is the learned dictionary, where the formula is updated for every j-th version of the dictionary, and where Af is the fine regularisation parameter.

21. The method of any preceding claim, further comprising obtaining an optimal phase from the learned dictionary.

22. The method of claim 21, further comprising estimating inter-layer attenuation of the subsurface region based on the optimal phase.

23. The method of any preceding claim, further comprising detecting seismic data of the subsurface region using a sensor array.

24. The method of any preceding claim, further comprising receiving seismic data of the subsurface region using a sensor array.

25. A non-transitory computer-readable medium having computer-program code stored thereon, the code executable by a processor to perform the method of any preceding claim.

26. A computing system comprising: a processor; and a memory having computer-program code stored thereon, the code executable by the processor to perform the method of any of claims 1 to 24.

Images