Simultaneous inversion of 4D data

Abstract

A method of processing time-lapse seismic data comprises providing a forward model configured to generate synthetic seismic data for a given combination of velocity ratio and density ratio, and for each of a plurality of given sampling depths of the subsurface region, iterating over different combinations of velocity ratio and density ratio to determine a given combination of velocity ratio and density ratio that minimises an objective function including a misfit term between (i) synthetic seismic data computed using the forward model and (ii) corresponding seismic data of the second seismic trace, wherein said ratios are defined at said sampling depth. Also disclosed is a method of determining a wavelet associated with seismic data as an analytical wavelet, wherein the analytical wavelet is an inverse Fourier transform of a product of at least two sigmoid functions.

Description

Technical field The present invention relates to a computer-implemented method for processing timelapse seismic data. The invention also relates to a computer program product with instructions for carrying out the same on a processor, and computer system for performing the same. Background Seismic data (e.g., a seismic trace) can be used to build a geophysical model (e.g., a velocity model) of a subsurface. In layman terms, seismic data is a representation of the reflectivity of the subsurface at different sampling or sample depths. The different sampling depths are captured based on an arrival time of the seismic waves that propagate from the acoustic signal source, through the subsurface, and which are reflected back to the acoustic signal receiver. Seismic data can be described as being 2D, 3D, or 4D. 2D seismic data is an image of a vertical “slice” of a subsurface (i.e., a 2D plane). 3D seismic data is an image of a horizontal stack of vertical “slices” of the subsurface. 4D seismic data is also referred to as time lapse seismic data and involves seismic data (2D, or 3D) recorded at different times. The skilled reader will understand that seismic data is representative of an image of the subsurface. The image may represent the subsurface in the spatial (e.g., real or reciprocal), time, or frequency domain. Time lapse seismic data is especially useful for better understanding hydrocarbon reservoir dynamics, such as fluid movement and pressure development, and to update existing geological models based on the observed differences in the time lapse seismic data. For this reason, time lapse seismic data is important for reservoir monitoring. For example, it is commonplace to record 3D seismic before production begins and to record one or more further 3D seismic thereafter, either periodically (e.g., on a yearly basis) or before well planning activities to be used as basis for a drill decision. Figure 1A shows an example of time lapse seismic data 100. The y-axis represents the amplitude of the recorded seismic signal, and the x-axis represents the travel time for the recorded seismic signal. The time lapse seismic data includes baseline seismic data 102 (herein baseline), denoted by the continuous line, and monitor seismic data 104 (herein monitor), denoted by the dotted line. A monitor survey refers to the acquisition survey of the monitor, and a baseline survey refers to the acquisition survey of the baseline. Generally, the baseline refers to seismic data that was acquired before the monitor, although this is not essential. A difference between the monitor and baseline indicates that the subsurface being imaged has changed. In particular, a difference in the measured arrival or travel time might indicate that the seismic waves are travelling at different velocities, whereas differences in the measured amplitude might indicate that the seismic waves are being reflected differently. This information can be used to infer how the subsurface has changed over time. Conventionally, the differences between a monitor and baseline are attributed to two independent sources: (i) differences in travel time (i.e. the velocity of the seismic wave); and (ii) differences in amplitude (i.e. the reflectivity of the subsurface). In a first step, time shifts are calculated to minimise the magnitude of the differences between a time-shifted monitor and the baseline. This is seen in Figure 1A as an offset in the x-direction between the monitor and baseline. In the specific example shown, the baseline appears to lag behind the monitor, and the time shifts assume positive values. Figure 1B shows the monitor following the time shift operation. In a second, subsequent step, amplitude shifts or amplitude differences are calculated to minimise the difference between the amplitude-and-time-shifted monitor and the baseline. The amplitude shifts required can be seen in Figure 1B as an offset in the y-direction between the monitor and baseline. In the specific example shown, the monitor appears to have, in general, a lower amplitude (i.e. lower reflectivity) than the baseline, and the amplitude shifts, therefore, assume a generally positive value. Figure 1C shows the monitor following the amplitude-and-time-shift operations. This approach leads to inaccuracies because differences in travel time and amplitude between the monitor and baseline cannot be regarded as being independent: velocity changes affect both travel time and reflectivity (i.e. amplitude). As a result, any changes in the model of the subsurface based on the time lapse seismic data analysis may be inaccurate. An improved approach for analysing 4D seismic data is desirable. Summary According to aspects of the present invention, there is provided a method of processing time-lapse seismic data, a computer system, a computer program product and a method of determining a wavelet for seismic data. According to a first aspect, there is provided a method of processing time-lapse seismic data. The time-lapse seismic data includes a first seismic trace (e.g., a baseline seismic trace) and a second seismic trace (e.g., a monitor seismic trace) captured of a subsurface region at different times (e.g., a first and second time, respectively). The second time may, in general, be after the first time. The method comprises providing a forward model (e.g. a “travel time and amplitude” forward model), which is configured to generate synthetic seismic data for a given combination of velocity ratio and density ratio; and iterating over different velocity and density ratio vector combinations (i.e., over a plurality of velocity and density ratio vector pairs) in order to determine a velocity and density ratio vector pair (i.e., an optimal pair combination) that minimises an objective function (e.g., as defined in Equation 9) that includes a misfit term between (i) a synthetic seismic trace computed according to a forward model (and which receives, as an input, the given combination of density and velocity ratio vector pair) and (ii) the second seismic trace. Equivalently, the method comprises, for each of a plurality of sampling depths of the subsurface region, iterating over different velocity ratio and density ratio combinations (i.e. over a plurality of velocity and density ratio pairs) that minimises an objective function (e.g., as defined in Equation 9) including a misfit term between (i) synthetic seismic data computed according to or using the forward model (which receives, as an input, the given combination of velocity and density ratio pair) and (ii) corresponding seismic data of the second seismic trace (i.e. the corresponding sampling point in the second seismic trace). Each of the velocity and density ratio pairs is associated with a corresponding sampling interval of the subsurface region. The sampling interval corresponds to the depth interval defined by neighbouring sampling points from the seismic trace. Within the sampling interval, the velocity and density ratio is taken to assume a constant value. The velocity and density ratio can be associated with any one of the sampling depths within that sampling interval. Each of the combinations is an estimate for a “true” or “absolute” velocity and density ratio combination, which, of course, cannot be fully resolved. Each velocity ratio is an estimate of a ratio between (i) a seismic wave velocity at a sampling depth i at a first time (i.e., when the first seismic trace is acquired) and (ii) a seismic wave velocity at a sampling depth i at a second time (i.e., when the second seismic trace is acquired). The seismic wave velocity may be referred to as being notional, as the velocity ratio can be deduced without experimentally obtaining a seismic wave velocity. Each density ratio is an estimate of a ratio between (i) the density of a subsurface at a sampling depth, i, at a first time (i.e., when the first seismic trace is acquired) and (ii) the density of the subsurface at the same sampling depth, i, at a second time (i.e., when the second seismic trace is acquired). Put differently, each velocity ratio is an estimate of a ratio between (i) a seismic wave velocity through a sampling interval (or at a sampling depth) of the subsurface, which has been captured, obtained, or imaged by the second seismic trace (i.e., during a monitor seismic trace survey) and (ii) a seismic wave velocity through the same sampling interval (or at the same sampling depth) of the subsurface, which has been captured, obtained, or imaged by the first seismic trace (i.e., during a baseline seismic trace survey). As such, the velocity ratio describes how the seismic wave velocity (e.g., p-wave velocity) through the subsurface has changed. If the value differs from one, the seismic wave velocity has changed. Each density ratio is an estimate of a ratio between (i) a density of the subsurface at a sampling depth (or a mean density through a sampling interval) obtained, captured, or imaged by the second seismic trace (i.e., during a monitor seismic trace survey) and (ii) a density of the subsurface at the same sampling depth (or a mean density through the same sampling interval) obtained, captured, or imaged by the first seismic trace (i.e., during a baseline seismic trace survey). As such, the density ratio describes how the density of the subsurface through the subsurface has changed. If the value differs from one, the density has changed. Between the iterations referred to above, both the velocity and density ratio are adjusted together (i.e. simultaneously). The method may further include, providing a further forward model (i.e. the travel time only forward model) configured to generate seismic data for a given velocity ratio), and, for each of the given plurality of sampling depths, iterating over different velocity ratios to determine an initialised velocity ratio that minimises an objective function which includes a misfit term between (i) synthetic seismic data computed according to a further forward model (i.e. the travel time only forward model) and (ii) corresponding seismic data of the second seismic trace (i.e. the corresponding sampling point in the second seismic trace). Each of the ratios are defined at the respective sampling depth in the plurality of sampling depths. The method further includes, iterating over different density ratios to determine an initialised density ratio that minimises an objective function which includes a misfit term (i) synthetic seismic data computed according to the forward model based on the corresponding initialised velocity ratio (i.e., the initialised velocity ratio computed for that sampling depth); and (ii) corresponding seismic data of the second seismic trace (i.e. the corresponding sampling point in the second seismic trace), wherein said ratios are defined at said sampling depth. The velocity ratio is not changed between iterations, i.e., it is held constant. In some examples, the initialised velocity ratio and the initialised density ratio are then used as an initial combination of velocity ratio and density ratio for the step of iterating over different combinations of velocity ratio and density ratio referred to above. The time-lapse seismic data may include a plurality of seismic trace pairs, each pair comprising a first and a second seismic trace of the same subsurface region, obtained or captured at different times. Initialised velocity and density ratios can be determined for each of the seismic trace pairs. The method may further comprise, for each of the plurality of sampling depths, iterating over different velocity ratios to determine an exploratory velocity ratio that minimises an objective function which includes a misfit term between (i) synthetic seismic data computed according to the first seismic trace; and (ii) corresponding seismic data of the second seismic trace (e.g., seismic data associated with a corresponding sampling point). Each of the ratios is defined at the corresponding sampling depth. The method then comprises determining a synthetic wavelet as a statistical or analytical wavelet based on the first or second seismic data. The synthetic wavelet comprises a wavelet scalar and a wave function and the forward model and the further forward model includes, as an input term, the determined synthetic wavelet. In some examples, the time lapse data includes a plurality of seismic trace pairs, each pair comprising a first and a second seismic trace of the same subsurface region, obtained at different times, and the method steps of the exploratory stage are carried out on a subset of said plurality of seismic trace pairs for which it is determined that a 4D anomaly is present. The wavelet scalar may be determined by iterating over different wavelet scalar values in order to determine a wavelet scalar value that minimises a summation of objective functions (i.e. a coupled objective function), wherein each objective function corresponds with one of the seismic trace pairs in the subset (i.e., which has been minimised according to the determined exploratory velocity ratios). In some examples, the further forward model (i.e., the travel time only forward model) assumes that, for each of the plurality of sampling depths, an acoustic impedance of a subsurface within the subsurface region remains unchanged between the first time and the second time. The method may include constructing a model of the subsurface region using the determined combinations of velocity ratio and density ratio. The subsurface region may be a combination of a hydrocarbon reservoir and a geological formation. The model may relate to a hydrocarbon reservoir and include determining a reservoir saturation and / or reservoir pressure, according to the model of the subsurface region. The synthetic wavelet may be computed as an analytical wavelet and the analytical wavelet is an inverse Fourier transform of a product of at least two sigmoid functions. In a specific example, the analytical wavelet is an inverse Fourier transform of a product of two sigmoid functions only. The method may then comprise: determining an amplitude spectrum of, or based on, the first or second seismic trace; determining fitting coefficients for each of the sigmoid functions (e.g., of the two sigmoid functions) based on the amplitude spectrum; and determining, using the fitting coefficients, the analytical wavelet as the inverse Fourier transform of the product of the sigmoid functions. The amplitude spectrum may be determined as the mean of the amplitude spectrums for each of the first or second seismic traces in the plurality of seismic trace pairs referred to above. That is, the amplitude spectrum may be representative of all the first or second seismic traces in the time lapse seismic data. According to a second aspect, there is provided a method of determining a wavelet associated with seismic data (e.g., a seismic trace) as an analytical wavelet, wherein the analytical wavelet is an inverse Fourier transform of a product of at least two sigmoid functions. The analytical wavelet may be an inverse Fourier transform of a product of two sigmoid functions only. The method may include determining an amplitude spectrum for a seismic trace in the seismic data; determining fitting coefficients for each of the sigmoid functions according to the amplitude spectrum; and determining the analytical wavelet as the inverse Fourier transform of the product of the sigmoid functions, using the fitting coefficients. In some examples, the amplitude spectrum is a mean of the amplitude spectrums for each seismic trace in the seismic data. That is, the amplitude spectrum may be representative of all the seismic trace associated with the seismic data. The fitting coefficients for each of sigmoid functions (e.g., a logistic function) can be determined based on a “best match” criterion, i.e., the fitting coefficients that minimise the misfit (e.g., the L2-norm) between the amplitude spectrum and the product of the sigmoid functions. The method of the second aspect may be incorporated into the method of the first aspect to, for example, determine its synthetic wavelet. According to third aspect, there is provided a computer system, comprising one or more processors configured to perform the methods according to the first and / or second aspects. According to a fourth aspect, there is provided a computer program product comprising instructions, when executed by a computer processor, cause the computer processor to perform any of the methods according to the first and / or second aspects. Brief Description Figure 1A, 1B, and 1C are graphical representations of seismic data; Figure 2 is a method flow diagram of an exploratory stage of a processing method; Figure 3 is a method flow diagram of an initialisation stage of the processing method; and Figure 4 is a method flow diagram of a simultaneous inversion stage of the processing method. Detailed Description In general terms, the present disclosure relates to a computer-implemented method of processing time lapse seismic data. The method includes three main stages: an exploratory stage, an initialisation stage and a simultaneous inversion stage. In the exploratory stage, a synthetic wavelet comprising a wavelet scalar and wave function is determined. Regularization parameters may also be computed in this stage. In the initialization stage, initial estimates for a velocity and density ratio are determined sequentially. These estimates are used as a starting point in the simultaneous inversion stage. In the simultaneous inversion stage, the optimal velocity and density ratios are determined together (rather than sequentially) by minimising an objective function over different combinations of ratios simultaneously. The method makes use of two different forward models: one which accounts only for travel time changes between the time lapse seismic data and which assumes the acoustic impedance of the subsurface remains unchanged, and one which accounts for both travel time and amplitude changes. The disclosure also proposes a computer program product with instructions for carrying out the same on a processor, and a computer system for performing the same. The disclosure also proposes a method of determining a wavelet as an analytical wavelet according to a product of at least two sigmoid functions. Forward model A forward model refers to an approach in which a synthetic seismic (e.g., a monitor trace) is generated from a synthetic seismic wave (referred to as a wavelet) and assumed geophysical parameters (e.g., acoustic impedance) of the subsurface being imaged. Two forward models are proposed; a “travel time only” forward model, and a “travel time and amplitude” forward model. These models are described in more detail below. For simplicity, the forward models below are described with reference to a pair of seismic traces, comprising a baseline and a monitor trace. It will be understood, however, that the forward models can be applied to a plurality of seismic trace pairs, for example, the plurality of seismic trace pairs may include a plurality of monitor traces to a given baseline trace or a plurality of baseline traces and a plurality of respective seismic traces. In mathematical terms, it can be shown that: (1) (f) — ) * ™(kMfU)) 2 dt where S<m) is the synthetic monitor seismic trace (herein synthetic monitor) with the superscript m indicating that the synthetic seismic trace relates to the monitor, w(t) is the wavelet function, which is a function of the travel time, t, (in Equation 1 w also includes a wavelet scalar term, c), Al(m) is the modelled acoustic impedance of the monitor subsurface with the superscript m again indicating that the acoustic impedance relates to the subsurface from the monitor survey, t is travel time, and * denotes the convolution function. Equation 1 shows that a synthetic monitor, S<m>, can be calculated by convolving a wavelet, w, with a band-limited reflectivity series. The band-limited reflectivity series is related to the rate of change of the modelled acoustic impedance, In general, the acoustic impedance of the subsurface being imaged or captured is not known. Instead of modelling the synthetic monitor from an acoustic impedance, Al(m), in absolute terms, it is more convenient to model the observed travel time and amplitude differences from an existing baseline. By imposing relative velocity and density changes to an experimental baseline, a synthetic monitor including both the travel time and amplitude difference can be modelled. In other words, the observed experimental monitor d(m) can be modelled using a function S(m)(t, d(b); qV, qp), where S(m) is the synthetic monitor, t is travel time, d(b) is the observed experimental baseline, and qV, qp are the model parameters, as defined in Equation 7. Modelling travel times only According to a first model, differences between the monitor and baseline are attributed to travel time changes only. This model is referred to as the “travel time only” model. In particular, it can be shown that, for the travel time only model, Equation 1 becomes: (2) — ~w(r) * ™ln( 2 ' at \ / where is the synthetic monitor (as with Equation 1), with the superscript m indicating that the synthetic seismic trace relates to the monitor survey, w(t) is the wavelet function, which is a function of travel time, t, (as with Equation 1, w also includes a wavelet scalar term, c), is the modelled acoustic impedance of the monitor subsurface but which is derived from the acoustic impedance of the subsurface being captured or imaged by the baseline survey (as denoted by the superscript b), t is travel time, and * denotes the convolution function. According to the first model, the synthetic monitor samples the same depths as the experimental baseline, i.e., the seismic trace acquired by the baseline survey. This means that each depth measurement for the observed experimental baseline corresponds to a depth measurement for the synthetic monitor, i.e., Z(t[b)) = Z(tfm>), where Z is a depth measurement made at travel time, t'b> or t[m), the subscript, I, denotes the sample depth, i, and the superscript b and m denote the baseline and monitor, respectively. The first model also assumes that the acoustic impedance of the subsurface being imaged by the baseline is identical to the acoustic impedance of the subsurface being imaged by the monitor. Put differently, the first model assumes there is no change in acoustic impedance. It follows, from these assumptions, that the acoustic impedance, Al<b)(t^b)), of the subsurface being imaged by the baseline at travel time, tfb) is equal to the acoustic impedance, Al<m>(ti<m>), of the subsurface being imaged by the monitor at travel time, tfm). The subscript i denotes the sample depth or sampling interval, i, which as noted above corresponds to equivalent depths. These constraints are denoted in Equations 3 and 4. It will also be understood that each expression for the acoustic impedance in Equation 3 and 4 is equal to any other one of the expressions because the travel time t[b) corresponds with the depth Z? and the travel time t(m) corresponds with the depth ZF. (4) It will be understood that tj(b) and t(m) are not equal to one another for the first model because it assumes that differences in travel times exist between the synthetic monitor and the experimental baseline. A synthetic monitor computed from an experimental baseline, according to the first model, will be defined in a different time domain. In particular, the time domain of the experimental baseline and monitor is {to(b), ti(b>, t2(b>, ..., yb)}, and the time domain of the synthetic monitor is {to(m>, ..., Um)], where (N+1) denotes the total number of samples (i.e. sample depths). As the sampling rate of the experimental baseline is typically constant, the time domain of the experimental baseline is generally uniform. This means that the time domain of the synthetic monitor will typically be non-uniform, following its computation from the experimental baseline. Put differently, this is because, according to the first model, the seismic wave velocities at a given travel time for the baseline and monitor are different, but the distance that the seismic waves have travelled is the same. In concrete terms, for normal incidence seismic traces, the travel time, t[m), for the synthetic monitor trace to sample depth i, can be expressed as: (5) tF - AtV . * - h 2. X.... iV where At is the constant sampling interval from the experimental baseline survey and qVj is a velocity ratio for sampling interval, j. It will be understood that, if there are N sample depths, there will be (N-1) sampling intervals. Each sampling interval, j, corresponds to a time interval, t / m)], and so sampling interval “1” relates to the time interval [to(m), ti<m>]. For this reason, the sampling interval j is referred to herein simply as sampling interval, i. The velocity ratio is defined as the ratio between the p-wave velocity, Vi(m)A / i(b), where Vj(m) is the p-wave velocity for the sampling interval, i, at the corresponding travel time of the monitor, and Vi(b) is the p-wave velocity for the sampling interval, i, at the corresponding travel time of the baseline, i, i.e. the p-wave velocities for the monitor and baseline at equivalent depths. It will be understood that the velocity ratio can instead be defined as the inverse of the ratio above. Equation 5 assumes that at time (t = 0 seconds), the depth sampled is 0 metres. In general, this will not be the case but the offset can be accounted for by including an offset term, to, to the right-hand-side of Equation 5. The general assumption is then that to = tdml. For normal incidence seismic traces, the travel time for the baseline trace to sample depth, i, which is sampled at regular time intervals, is / At, where At is the sampling rate from the experimental baseline survey. In other words, according to this model, it is possible to transform the experimental baseline to a synthetic monitor by mapping the amplitudes from the baseline time domain to a new time domain, according to the velocity ratios, qV\, or vice versa according to 1 / qV,. That is, each sampling depth i for the baseline is mapped to a corresponding sampling depth i (which samples the same depth) according to the corresponding velocity ratio, qV\„ associated with the corresponding sampling interval, i. As has already been noted, this leads to a warping of the time domain from a uniform one (e.g., of the experimental baseline) to a non-uniform one (e.g., of the synthetic monitor). Referring back to Equation 3, if the values for Al(b> at {ti(b), t2, ..., tN} are known, it follows that the corresponding values for Al<m) at t2(m).....t^m>} are also known, where 1, 2, ..., N denote the sample depths, i. At the same time, it has been demonstrated that the mapping between the baseline travel times to the monitor travel times can be determined according to the velocity ratios. As a result, according to the travel time only model, the acoustic impedance of the subsurface being imaged by the baseline can be mapped from the uniformly spaced time domain to the non-uniformly spaced travel times, according to the velocity ratios {qVi, qVz, ..., qV^-i} for each of the (N-1) sampling intervals. The velocity ratios can be determined by an inversion technique. Possible inversion techniques are described in more detail below. In practice terms, it is desirable that the experimental monitor and the synthetic monitor are defined on the same or common time domain. Any sampling in the time domain could be used, however, the sampling of the experimental baseline and monitor is a convenient choice (herein baseline time domain). In other examples, the same or common time domain has a denser sampling (i.e. a greater number of samples) than the baseline time domain. The acoustic impedance of the baseline, following the mapping operation and the interpolation operation back to the baseline time domain, is denoted as in Equation 2. The tilde denotes that the acoustic impedance has been interpolated from the synthetic monitor time domain to the baseline time domain. Interpolation advantageously reduces warping of the seismic data. Possible interpolation techniques are described in more detail below. As this model assumes that the difference between the experimental baseline and the synthetic monitor are attributed to travel times only, the experimental baseline itself, d(b) can be transformed into an approximate synthetic monitor by mapping the samples from the baseline to the monitor time domain, according to the velocity ratios. Unlike some conventional techniques, no attempt is made to directly recover Al and AJ^ from the experimental baseline. After the mapping operation, the synthetic monitor can then be re-sampled in the baseline time domain, using interpolation techniques. The experimental baseline, following the mapping and interpolation operations, is referred to as d^. The tilde denotes that an interpolation operation has been performed. It is approximately equal to the synthetic monitor from Equation 2. Modelling travel times and amplitude changes According to a second model, differences between the monitor and baseline are attributed to both travel time changes and amplitude changes. This model is referred to as the “travel time and amplitude” model. In particular, it can be shown that, for the travel time and amplitude model, Equation 1 becomes: where S(m) is the synthetic monitor seismic trace (herein synthetic monitor), with the superscript m indicating that the synthetic seismic trace relates to the monitor, w(t) is the wavelet function, which is a function of the travel time, f, c is a wavelet scalar, qAI is a ratio of the acoustic impedances of the subsurface being captured or imaged by the monitor and the baseline for a given sampling interval or depth, t is travel time, is the synthetic monitor according to the first model as noted above, and * denotes the convolution function. It will be understood that the ratio of the acoustic impedances can instead be defined as the inverse of the ratio above. The second term on the right hand side of Equation 6 would then change sign. As with the first forward model, the second forward model assumes that the synthetic monitor samples the same depths as the baseline. Unlike the travel time only model however, the second model assumes that both the acoustic impedances of the subsurface being imaged in the monitor and baseline, and the travel times from the monitor and baseline, are different. Equations 3 and 4, therefore, no longer apply. Instead, the acoustic impedance of the monitor at a particular sampling interval or depth in Equation 1 is given by the product of the acoustic impedance of the baseline and the acoustic impedance ratio between the monitor and baseline, qAI, at that sampling interval or depth. Equation 6 represents a refinement of the model in Equation 2, in which variations in acoustic impedance between the monitor and baseline are taken into account. It is worth noting that, in Equation 2, it is assumed that there are no such variations and so the ratio between the monitor and baseline acoustic impedances is equal to one and the second term on the right-hand side of Equation 6 is zero. As there are variations in acoustic impedance in this model, the estimation of the wavelet scalar, c, becomes important. Different techniques for estimating the wavelet scalar are described in more detail below. Referring to Equation 6, the inventors have noted that, it is possible to arrive at any given synthetic monitor S<m> from different combinations of wavelet scalar and acoustic impedance ratio. For example, the effect of a large wavelet scalar and a low acoustic impedance ratio in the second model is indistinguishable from the effect of a small wavelet scalar and a large acoustic impedance ratio. In other words, the output of the forward model in Equation 6 is nonunique. Approximating the wavelet scalar is of lesser importance for the “travel time only” model because the wavelet scalar does not affect travel times and the model assumes there are no amplitude variations between the monitor and the baseline. The first model can, therefore, be used to determine a realistic wavelet scalar, whereas much more care is required in the second model. The acoustic impedance ratio, qAI, is defined as: where the superscript (m) and (b) denote the monitor and baseline, respectively, and i denotes the sampling interval. The acoustic impedance ratio is also equal to the product of the density ratio, qpi, and the velocity ratio, qV,. The density ratio, qp,, is defined as the ratio between subsurface density of the subsurface at a sampling interval or sampling depth i at a second time (e.g., as obtained, captured or imaged by the monitor survey) and the subsurface density of the subsurface at the same sampling interval or sampling depth i at a first time (e.g., as obtained, captured or imaged by the baseline survey). The density ratio may alternatively be defined as the inverse of the ratio above. Combining Equations 6 and 7, leads to: Equation 8 demonstrates that, in addition to the contribution made by velocity changes to the “travel time only” model, velocity changes also effect the amplitude of the synthetic monitor, which is parameterised by the velocity ratio, qV. A further contribution to the synthetic monitor is made by density changes, which is parameterised by the density ratio, qp. More specifically, an increase in each of the velocity and / or density ratio causes a corresponding increase in the amplitude of the synthetic monitor. As noted in the background, the contribution made by velocity changes to the amplitude of the monitor trace has not been fully considered in the prior art. The travel time and amplitude forward model is an improvement on these earlier techniques. As with the travel time only model, it is preferable, but by no means essential, to reevaluate the synthetic monitor into a different time domain, such as the time domain of the baseline. Correspondingly, the density ratio, qp and velocity ratio, qV, may also be resampled. Possible interpolation techniques are described in more detail below. The tilde in Equations 6 and 8 denotes parameters which have been re-sampled by interpolation. Inversion Inversion refers to an approach in which changes in geophysical properties of a subsurface are calculated from a synthetic seismic wave (referred to as a wavelet) and an experimental monitor (i.e., a monitor trace that has been acquired through a monitor survey). In the context of this disclosure, the geophysical properties of interest are qV, qp and / or qAI. In inversion, an objective function, which describes the misfit between the experimental and synthetic seismic, is minimised, according to a model parameter, m. A regularization or penalty term may also be added to the objective function to promote certain characteristics of the model parameter. The synthetic monitor can be determined according to any one of the forward models or by mapping an experimental baseline, as proposed above. The model parameter m may be qV, qp, or [qV, qp]. For example, if the synthetic monitor in the objective function is computed according to the travel time only forward model, then m may be qV, i.e., [qVi, qV2, ..., qVN-i], whereas if the synthetic monitor in the objective function is computed according to the travel time and amplitude forward model, then m may be qp or [qV, qp], i.e., [qpi, qp2, .... qpN-i], or [[qVi, qV2, ..., qVN-i], [qpi, qp2, ..., qpN-i]], wherein 1, 2, ..., N-1 denote the sampling interval with which the velocity and density ratio correspond with. In mathematical terms, the objective function, l(m), is parameterised by model parameter, m, and may be given by: (9) — ||Sm(m) - dm||$ 4' where dm is the experimental monitor, Sm is the synthetic monitor as obtained by (i) transforming the experimental baseline trace according to a velocity ratio vector, (ii) the travel time only forward model or (iii) the travel time and amplitude forward model described above, A is a regularization parameter, L is a regularization operator, and p indicates the order (e.g., 0, 1, 2) of norm of the regularization term. A||Lm||p is, collectively, referred to as the regularization term, and ||Sm - dm||22 is referred to the L2-norm of the residual vector, Sm - dm . In general terms, the regularization term imposes a penalty term to help ensure that the optimal solution to the objective function is unique and to steer the optimisation problem away from unrealistic solutions. As inversion involves determining a difference between an experimental and synthetic monitor, it will be understood that the time domains with which the traces respectively sample are common or the same. This can be achieved using interpolation techniques, which are described in more detail below. The inverse problem can be posed as a minimisation of the sum of the L2-norm of the residual vector and the regularization term (as shown in Equation 9). Example methods to perform these minimisation steps include: steepest descent, conjugate gradient, Gauss-Newton, Levenberg-Marquardt and the like. These methods are known to the skilled reader, per se. In a specific example, the total variation regularization technique is applied to solve these inverse problems. The regularization operator, L, becomes a forward gradient operator, with p equal to 1. Total variation is known to the skilled reader, per se. The regularization parameter, A, can be determined by, for example, generalised cross validation or L-curve techniques. These techniques are known to the skilled reader. It will be understood that alternative or additional regularization terms can be used, as needed. Each additional regularization term may be subject to an additional constraint, e.g., constraints to lateral variation and other statistical norms, such as the L2 norm. Although the objective function in Equation 9 includes an L2-norm between the experimental and synthetic monitor, the skilled reader will understand that other statistical norms may be used instead, such as the L1-norm. Interpolation As has already been noted, parameters with a tilde denote that they have been resampled by an interpolation operation into a different time domain. Typically, although not always, this is the baseline time domain. Different interpolation techniques can be applied to modify the domain in which a given parameter is sampled. The technique deployed may be parameter-dependent. For example, piecewise linear interpolation can be used, when interpolating the geophysical parameters: qV, qp, and qAI, whereas cubic interpolation can be used to interpolate seismic data (e.g., a seismic trace). This is because subsurface properties, such as velocity and density, may change abruptly, which makes piecewise linear interpolation a better choice. Cubic interpolation, on the other hand, gives a smooth result and generally provides a better estimate for waveform-like data, such as seismics. In cubic interpolation, boundary conditions are used to solve for the polynomial coefficients. A common set of boundary conditions are the so-called “natural” boundary conditions, in which the second derivative at each end of the range is set to be zero. With these boundary conditions, 4N polynomial coefficients can be uniquely determined, wherein N is the number of intervals (i.e. interpolation intervals). The use of boundary conditions to solve for these polynomial coefficients is known to the skilled reader. Wavelet scalar approximation In order to model the observed amplitude differences and acoustic impedance changes accurately, an appropriate wavelet scalar is preferably computed. As previously stated, there is an inherent non-uniqueness with the travel time and forward model, as the effect of a large wavelet scalar and a low acoustic impedance ratio is, in essence, indistinguishable from the effect of a small wavelet scalar and a large acoustic impedance ratio. If the wavelet scalar is incorrectly defined, errors can be introduced into the forward model. On the other hand, the wavelet scalar is of lesser importance for the travel time only forward model, as it has much less of an effect on travel times. This means that an estimate of the velocity ratio, qV, can be obtained more accurately using the travel time only forward model. By consequence, the wavelet scalar can be estimated more accurately by solving an inversion problem that depends on the use of the travel time only forward model (i.e., in which the model parameter, m, is qV). The wavelet scalar can be regarded as a hyper-parameter that is common or shared for the objective functions described above. The wavelet scalar facilitates the solving of an inversion problem in the initialization and simultaneous inversion stages of the workflow described below. In an example, the wavelet scalar can be estimated according to Equations 6 or 8. Alternatively, the wavelet scalar can be estimated by determining a maximum and minimum qAI, according to amplitude differences between the seismic trace pair (after the travel time differences have been accounted for) and attributing the amplitude difference to a corresponding maximum / minimum reflectivity change. As the experimental monitor generally comprises a plurality of seismic traces, the wavelet scalar can be approximated as being a parameter that is common to all of the seismic traces. As a result, an optimised value for the wavelet scalar can be computed by minimising a coupled objective function, given as the sum of the “individual” objective functions. Each “individual” objective function being the objective function used to optimise a model parameter, m, corresponding to a given seismic trace of the experimental monitor. In yet another alternative, an optimised value for the wavelet scalar is calculated using the generalized cross validation approach. Workflow Figures 2, 3, and 4 are method flow diagrams that describe how time lapse data is processed according to different stages of the proposed workflow. The time lapse seismic data includes an experimental baseline seismic trace (herein experimental baseline) and an experimental monitor seismic trace (herein experimental seismic). The experimental baseline represents a first image captured of a subsurface region. The experimental monitor represents a second image captured of the same subsurface region at a different time. In general, the experimental baseline and monitor are different due to changes to the subsurface of the subsurface region. In a specific use case, the changes to the subsurface are attributed to changes to a reservoir of fluid (e.g., pressure and / or saturation changes in a hydrocarbon reservoir), while the subsurface geology remains substantially the same. Figure 2, 3, and 4 respectively describe the (i) exploratory; (ii) initialisation; and (iii) simultaneous inversion stages of the proposed workflow. The three main stages are described in more detail below, with reference to these Figures. The purpose of the exploratory stage is to estimate a wavelet, establish a wavelet scalar and, if relevant, determine the regularization parameters. The exploratory stage is performed on a subset of the seismic data in which 4D anomalies are determined to be present. The exploratory stage then involves solving a simultaneous inversion problem (potentially many times) in order to evaluate the wavelet scalar and, if applicable, regularization parameters.. The purpose of the initialization stage is to determine an initial estimate for the velocity and density ratios, for use as a starting point or solution for the simultaneous inversion stage. In this stage, both subsets of the seismic data which do and do not show 4D anomalies, are used. The initialization stage makes use of the travel time and amplitude forward model, but the velocity and density ratios are determined one after the other. In brief, it involves 1) solving for qV for the whole dataset, 2) accounting for the travel time differences between the experimental baseline and monitor according to 1) and 3) keeping qV fixed and solving for either qp or qAI. The purpose of the simultaneous inversion stage is to determine optimal velocity and density ratios concurrently. This advantageously means that the effect of velocity changes on the amplitude the monitor seismic can be accounted for. Exploratory stage In optional step 202, a subset of the time lapse seismic data is selected for processing. The subset corresponds to a collection of experimental baseline / monitor seismic trace pairs for which the amplitude difference between the experimental baseline and experimental monitor exceeds a preset threshold. This set is referred to as a 4D anomaly region, since it is an indication of a subsurface area in which the subsurface properties have changed between the baseline and the monitor survey. The 4D signal to noise ratio, SNR, in 4D anomaly regions is high, compared to other regions of the time lapse seismic data. The use of 4D anomaly regions means that the parameters determined in the exploratory stage are less vulnerable to noise. That said, the method of Figure 2 can instead be applied to all the time lapse seismic data, or any subset thereof. In step 204, an objective function, which is parameterised by the model parameter -velocity ratio vector, qV, only, is minimised by determining a synthetic monitor trace, computed from the experimental baseline trace, which minimises the objective function. As described above, this may include minimising the sum of (i) a misfit or residual term between the synthetic monitor trace and the experimental monitor trace and (ii) a regularization term. In some examples, the objective function is taken as being minimised, according to a convergence criteria or criterion. As an example, an objective function may be taken as minimised, if, between successive iterations, a fractional difference between the value assumed by the respective objective functions, is less than a predetermined threshold. The synthetic monitor traces determined in step 204 are referred to as exploratory synthetic monitor traces, and they can be computed by transforming the experimental baseline trace, according to different velocity ratio vectors (as described in relation to the travel time only model). As noted above, this transformation constitutes a mapping from the baseline time domain to a different time domain, according to the velocity ratio vector. The velocity ratio vector, which minimises the objective function, is referred to as, qVtt, and the exploratory synthetic monitor trace computed according to this velocity ratio vector is referred to as dtt. The subscript “tt” denotes that the velocity ratio vector and synthetic monitor trace account for differences in travel time only. Optionally, the exploratory synthetic monitor trace may be interpolated into yet another different time domain (e.g., back into time domain of the experimental baseline). The velocity ratio and exploratory synthetic monitor trace, following the interpolation operation, are referred to as qV\t and d^t, where the tilde denotes the interpolation operation. Suitable interpolation techniques have been described in detail above. In step 206, a wavelet is estimated as a statistical wavelet or an analytical wavelet using the experimental baseline and / or the experimental monitor data. The statistical or analytical wavelet can be calculated by calculating the autocorrelation of each seismic trace and transforming the result into the frequency domain using a Fourier Transform, such as a FFT. This gives the power spectrum of each trace. From these power spectra, an average amplitude spectrum, Aavg(f), can be calculated. A statistical wavelet is evaluated by transforming the average amplitude spectrum into the time domain, using an Inverse Fourier Transform, such as a IFFT. The statistical wavelet can be filtered or unfiltered, as necessary. Conventional filtering techniques, known to the skilled reader, can be applied either in the frequency domain or in the time domain. Alternatively, an analytical wavelet is computed by determining, which analytical wavelet, best matches the average amplitude spectrum, Aavg. Example analytical wavelets include Ricker, Ormsby, Klauder, Morlet, or Butterworth. Preferably, although not necessarily, the average amplitude spectrum can be modelled as the product of two logistic (sigmoid) functions, each sigmoid function being given by Equation 10: A (10) x where A, k, and fo are fitting coefficients and f, is the frequency of the wavelet. The analytical wavelet is given by Equation 11 and exhibits zero phase: (11) w(t) J*''^1( / )^2(--^ where F1 is the Inverse Fourier Transform operator, and 01 and 02 are each according to Equation 10, having different sets of fitting coefficients. In particular, with one of the logistic functions being adjusted to reflect the high-end of the amplitude spectrum and the other being adjusted to reflect the low-end of the amplitude spectrum. The fitting coefficients can be determined as those which best fit the average amplitude spectrum, Aavg(f). In step 208, a wavelet scalar is computed, according to any one of the approaches described above. In an example, the wavelet scalar is estimated by determining a maximum and minimum qAI according to amplitude differences between the seismic trace pairs. As noted above, this estimate takes place after travel time only amplitude differences have been accounted for, i.e., after the experimental baseline has been transformed by the velocity ratio vector obtained from step 204. Step 208 may comprise minimising an objective function, which includes the model parameters - velocity ratio vector, qV, and density ratio vector, qp, by determining a synthetic monitor trace according to a combination of velocity and density ratio vectors. That is, for each different combination of velocity and density ratio vector, both the velocity and density ratio vector are adjusted simultaneously. As described above, this may include a) minimising the sum of (i) a misfit or residual term between the synthetic monitor trace and the experimental monitor trace and (ii) a regularization term. The synthetic monitor trace is determined using the travel time and amplitude forward model, with the wavelet, wavelet scalar, and the exploratory synthetic monitor trace, dtt, as the synthetic monitor, db, in Equation 6 or 8. The velocity ratio vector obtained in step 204 can be used as an initial guess for the minimisation problem posed by step 208. The velocity ratio vector and density ratio vector, which minimise the objective function, are referred to herein as the optimal exploratory velocity and density ratio vector, qVe, qpe, where the subscript “e” denotes the exploratory stage. In some examples, the experimental baseline and experimental monitor seismic data includes a plurality of seismic traces. Each seismic trace from the baseline seismic data is then associated with a corresponding seismic trace from the monitor seismic data so as to form a plurality of seismic trace pairs. Steps 202 to 206 or steps 202 to 208 may be repeated for each seismic trace pair, or for a subset of these seismic trace pairs. Inversion algorithms to solve the minimisation problem posed by steps 204 and 208 have been described in detail above. The inversion problem posed in steps 204 and 208 may include one or more regularization parameters. These regularization parameters may be used, if, for example, convergence fails. In steps 204 and 208, the synthetic monitor traces computed are resampled in a different time domain (e.g., back to the time domain of the experimental baseline) prior to inversion taking place. Initialization stage In the exploratory stage, inversion is carried out on a subset of the seismic data in which it has been determined that 4D anomalies are present. Conversely, in the initialization stage, seismic data, both which shows 4D anomalies and which does not, is used to compute an estimate for the optimal values for the model parameters, qV, and qp, where the subscript i denotes the initialization stage. In step 302, an objective function, which is parameterised by the velocity ratio vector, qV, only, is minimised by determining a synthetic monitor trace according to the travel time only model that minimises the objective function. As described above, this may include a) minimising the sum of (i) a misfit or residual between the synthetic monitor trace and the experimental monitor trace and (ii) a regularization term. In some examples the objective function for the velocity ratio is taken as being minimised if, between successive iterations, the difference between the value assumed by the respective objective functions is less than a predetermined threshold. More specifically, synthetic monitor traces are computed using the travel time only forward model, using the wavelet scalar, c, wavelet function, w, and regularization parameters (if applicable) determined in the exploratory stage, for different velocity ratio vectors, qV. The optimal exploratory velocity ratio vector, qVe, can be used as an initial value for the velocity ratio vector during this minimisation process (for the subset of traces used in the exploratory stage). The velocity ratio vector that minimises the objective function is referred to as the optimal initialised velocity ratio vector, qV,. In step 304, an objective function, which is parameterised by the density ratio vector, qp, only, is minimised by determining a synthetic monitor trace, according to the travel time and amplitude forward model that minimises the objective function. As described above, this may include a) minimising the sum of (i) a misfit or residual between the synthetic monitor trace and the experimental monitor trace and (ii) a regularization term. The synthetic monitors are computed according to the travel time and amplitude forward model (e.g. Equation 8). The value for the velocity ratio, qV, in Equation 8 is taken as the initialised optimal velocity ratio computed in step 302 - it is held constant throughout step 304. The density ratio vector which minimises the objective function is referred to as the initialised density ratio vector, qps. In some examples, the objective function for the density ratio is taken as being minimised according to a convergence criterion or criteria, which, as mentioned above, may be that a fractional difference between the value assumed by the respective objective functions is less than a predetermined threshold. It will be understood that the optimal initialised density ratio vector can also be determined by computing the acoustic impedance ratio vector (in substantially the same way as described in step 304) and then dividing the result by the velocity ratio output in step 302. Inversion algorithms to solve the minimisation problem posed by steps 302 and 304 have been described in detail above. Unlike steps 204, 208, and 302, step 304 is a linear inversion problem, which is easier to solve than a non-linear inversion problem. This is essentially because it only involves varying one of the velocity and density vector ratios. In this stage, the objective functions may each include a regularization parameter determined from the exploratory stage. In steps 302 and 304, the synthetic monitor traces computed are resampled in a different time domain (e.g., back to the time domain of the experimental baseline) prior to inversion taking place. Simultaneous inversion In step 402, an objective function, which is parameterised by the velocity ratio vector, qV, and the density ratio vector, qp. is minimised by computing a synthetic monitor trace, according to the travel time and amplitude forward model, that minimises the objective function. As described above, this may include a) minimising the sum of (i) a misfit or residual between the synthetic monitor trace and the experimental monitor trace and (ii) a regularization term. In some examples, the objective function is taken as being minimised according to a convergence criteria (as described above in the exploratory and initialisation stages). Unlike in step 302 and 304, where the different synthetic monitors are computed by changing only one of the velocity or density ratio vectors, in step 402, the synthetic monitor traces are computed by changing the velocity and density ratio in combination, i.e. simultaneously. That is, for each different combination of velocity and density ratio vector, both the velocity and density ratio vector are adjusted, or more specifically, both the velocity ratio and density ratio for each sampling depth or interval are adjusted between successive iterations. In this way, the velocity and density ratio vectors are optimised together, and the optimal values are determined simultaneously. The optimal velocity and density ratios are denoted qV and qp, respectively. The synthetic monitor traces are computed using the travel time and amplitude forward model, i.e., using Equation 8, with the wavelet scalar, c, wavelet function, w, and regularization parameters, if applicable, determined from the exploratory stage, for different velocity and density ratio vectors. The optimal initialised velocity and density ratio vectors output in the initialisation stage can be used as a starting point for the minimisation problem. The initialised values for the velocity and density ratio help to ensure that the minimisation procedure starts near a local minimum. Ultimately, this reduces the number of iterations required during the inversion process (compared with selecting initial values at random) and steers inversion away from unrealistic solutions. As with the previous stages, the synthetic monitor traces computed are resampled in a different time domain (e.g., back to the time domain of the experimental baseline) prior to inversion taking place. Inversion algorithms to solve the minimisation problem posed by steps 402 have been described in detail above. In optional step 404, a velocity and / or a density model of the subsurface formation being captured by the monitor survey is constructed based on the optimal velocity and / or optimal density ratio vectors. With the proposed approach, more accurate combinations of velocity and density ratio can be determined. This, in turn, improves the accuracy of the model being constructed, which is critical for making subsequent operational decisions. The approach is particularly useful for monitoring hydrocarbon reservoirs because it provides a better basis for assessing the contribution made to the 4D anomaly by changes in reservoir saturation and pressure. The method is not limited to hydrocarbon wells, however. More generally, it finds use in any processes which involve the injection and / or extraction of fluid into or from a subsurface formation, for example, in geological carbon sequestration. The ordering of the steps in Figures 2 to 4 is not intended to impose a strict order in which the method steps are performed. The skilled reader will understand that the method steps can be performed in other working orders. The predetermined thresholds in the exploratory, initialisation, and simultaneous inversion stages described in Figure 2, 3, and 4 may be different or the same. The skilled reader will understand that each of these thresholds can be set based on a known trade-off between computational cost and accuracy. Although the invention has been described in terms of preferred embodiments as set forth above, it should be understood that these embodiments are illustrative only and that the claims are not limited to those embodiments only. Features from each of the different stages may be combined or omitted as appropriate to form other working examples.

Claims

1. A computer-implemented method of processing time-lapse seismic data, wherein the time-lapse seismic data includes a first seismic trace for a subsurface region obtained at a first time and a second seismic trace for the same subsurface region obtained at a second time that is different from the first time, the method comprising:providing a forward model configured to generate synthetic seismic data for a given combination of velocity ratio and density ratio;for each of a plurality of given sampling depths of the subsurface region, iterating over different combinations of velocity ratio and density ratio to determine a given combination of velocity ratio and density ratio that minimises an objective function including a misfit term between (i) synthetic seismic data computed using the forward model and (ii) corresponding seismic data of the second seismic trace, wherein said ratios are defined at said sampling depth, wherein a velocity ratio is an estimate of a ratio between (i) a seismic wave velocity at the given sampling depth at a second time and (ii) a seismic wave velocity at the given sampling depth at a first time, andwherein a density ratio is an estimate of a ratio between (i) a density of a subsurface at the given sampling depth at a second time and (ii) a density of a subsurface at the given sampling depth at a first time.

2. The method of claim 1, wherein, for each successive iteration, both the velocity ratio and density ratio are adjusted.

3. The method of claim 1 or 2, further comprising:providing a further forward model configured to generate synthetic seismic data for a given velocity ratio;for each of the plurality of given sampling depths:iterating over different velocity ratios to determine an initialised velocity ratio that minimises an objective function including a misfit term between (i) synthetic seismic data computed using the further forward model and (ii) corresponding seismic data of the second seismic trace, wherein said ratios are defined at said sampling depth;iterating over different density ratios to determine an initialised density ratio that minimises an objective function including a misfit term (i) synthetic seismic data computed using the forward model based on the corresponding initialised velocity ratio; and (ii) corresponding seismic data of the second seismic trace, wherein said ratios are defined at said sampling depth; andusing the initialised velocity ratio and the initialised density ratio as an initial combination of velocity ratio and density ratio for the step of iterating over different combinations of velocity ratio and density ratio.

4. The method according to claim 3, in which time-lapse seismic data includes a plurality of seismic trace pairs, each pair comprising a first and a second seismic trace of the same subsurface region, obtained at different times, and the method steps of claim 3 are performed for each of the seismic trace pairs.

5. The method according to claim 3 or 4, further comprising:determining a synthetic wavelet as a statistical or an analytical wavelet based on the first or second seismic trace,wherein, the synthetic wavelet comprises a wave function and the forward model and the further forward model include, as an input, the determined synthetic wavelet.

6. The method according to claim 5, comprising:for each of the plurality of given sampling depths,iterating over different velocity ratios to determine an exploratory velocity ratio that minimises an objective function parameterised by a misfit between (i) synthetic seismic data computed according to the first seismic trace; and (ii) seismic data associated with a corresponding sampling point of the second seismic trace, wherein said ratios are defined at said sampling depth.

7. The method according to claim 5 or 6, wherein the time lapse data includes a plurality of seismic trace pairs, each pair comprising a first and a second seismic trace of the same subsurface region, obtained at different times, and the method steps of claim 5 are carried out on a subset of said plurality of seismic trace pairs for which it is determined that a 4D anomaly is present.

8. The method according to claim 7, in which the wavelet scalar is determined by:iterating over different wavelet scalar values in order to determine a wavelet scalar value that minimises a summation of objective functions, wherein each objective function corresponds with one of the seismic trace pairs in the subset..

9. The method according to any of claims 3 to 8, when dependent on claim 3, in which the further forward model assumes that, for each of the plurality of given sampling depths, an acoustic impedance of the subsurface within the subsurface region remains unchanged between the first time and the second time.

10. The method according to any one of the preceding claims, further comprising: constructing a model of the subsurface region using the determined combinations of velocity ratio and density ratio.

11. The method according to claim 10, in which the model relates to a hydrocarbon reservoir, and the method further comprises:determining a reservoir saturation and / or reservoir pressure, using the model of the subsurface region.

12. The method according to any one of claims 5 to 11, wherein the synthetic wavelet is computed as an analytical wavelet and the analytical wavelet is an inverse Fourier transform of a product of at least two sigmoid functions.

13. The method of claim 12, in which the analytical wavelet is an inverse Fourier transform of a product of two sigmoid functions only.

14. The method of claim 12 or 13, comprising:determining an amplitude spectrum based on the first or second seismic trace;determining fitting coefficients for each of the sigmoid functions based on the amplitude spectrum; anddetermining, using the fitting coefficients, the analytical wavelet as the inverse Fourier transform of the product of the sigmoid functions.

15. The method according to claim 14, in which the amplitude spectrum is a mean of the amplitude spectrums for each of the first or second seismic traces in the plurality of seismic trace pairs.

16. A computer system comprising one or more processors configured to perform the method steps of any one of claims 1 to 15.

17. A computer program product comprising instructions, which, when executed by a computer processor, cause the computer processor to perform the method steps of any one of claims 1 to 15.

18. A method of determining a wavelet associated with seismic data as an analytical wavelet, wherein the analytical wavelet is an inverse Fourier transform of a product of at least two sigmoid functions.

19. A method of claim 18, in which the analytical wavelet is an inverse Fourier transform of a product of two sigmoid functions only.

20. A method of claim 18 or 19, comprising:determining an amplitude spectrum of a seismic trace in the seismic data;determining fitting coefficients for each of the sigmoid functions based on the amplitude spectrum; anddetermining the analytical wavelet as the inverse Fourier transform of the product of the sigmoid functions, using the fitting coefficients.

21. A method according to claim 20, in which the amplitude spectrum is a mean of the amplitude spectrums for each seismic trace in the seismic data.

Images