SYSTEM AND METHOD FOR LEAST-SQUARES MIGRATION OF TIME-LAPSE SEISMIC DATA

Abstract

A least-square migration, LSM, based method for generating a 4D image of a subsurface, the method including receiving seismic data d related to the subsurface, the seismic data d including a baseline dataset dB and a monitor dataset dM, calculating a baseline filter B and a monitor filter M based on a same common reflectivity r of the subsurface and corresponding remigrated baseline data mB1 and remigrated monitor data mM1 so that the base filter B applied to the remigrated baseline data mB1 equals the monitor filter M applied to the remigrated monitor data mM1, applying the baseline filter B to raw migrated baseline data mB0 and applying the monitor filter M to raw migrated monitor data mM0 to generate LSM baseline data mB and LSM monitor data mM, and generating the 4D image of the subsurface based on the LSM baseline data mB and the LSM monitor data mM.

Description

BACKGROUND OF THE INVENTION

Technical Field

Embodiments of the subject matter disclosed herein generally relate to a system and method for migrating seismic data associated with a subsurface, and more particularly, to a time-lapse least-square migration (LSM) method that applies different non-stationary filters, for the baseline and monitor surveys, while using a common reflectivity. The method may be a single iteration LSM or a plural iteration LSM.

Discussion of the Background

Seismic data is acquired to generate a profile (image) of the geophysical structure (subsurface) under the surface. While this profile does not provide an accurate location for a desired resource (e.g., oil and gas) reservoir, it suggests, to those trained in the field, the presence or absence of the resource reservoirs. Thus, providing a high-resolution image of the subsurface is an ongoing process for the exploration of natural resources.

The seismic data used for generating the image of the subsurface may be acquired in various ways, for example, with land equipment, marine equipment, ocean bottom equipment, autonomous underwater equipment, or aerial equipment. The type of sensors and sources used in any of these scenarios are known in the art and thus, not repeated herein. After the seismic data is acquired, it needs to be processed to eventually generate the image of the subsurface. There are many known algorithms for processing the seismic data. Migration is one such algorithm that is used for seismic data analysis. Migration allows for the correct positioning of subsurface structures and reflectivity. The migration algorithms are ray-based, for example, Kirchhoff migration as described in Schneider, 1978, (Schneider, W., 1978, Integral formulation for migration in two and three dimensions, Geophysics, 43, 49-76) and wave equation based, for example, Reverse-time migration (RTM) Hemon, 1978 (Hemon, C., 1978, Equations d'onde et modeles. Geophysical Prospecting, 26, 790-821).

Kirchhoff migration and Reverse-time migration are inversion methods where the subsurface reflectivity is extracted from the recorded seismic data. They are influenced by the acquisition setup and the geological complexity of the subsurface. The images may suffer from lack of resolution, bad event continuity and zones with dimmed amplitude, due to band-limited sources and uneven illumination patterns of seismic waves.

A particular application for migration is time-lapse (4D) processing where images of the same underground region are obtained from two or more different datasets, e.g., baseline and monitors datasets, which are recorded at different times, usually for the purpose of reservoir monitoring (i.e., to monitor the changes in the reservoir by determining the changes between the monitors and the baseline). The level of the 4D signal (signal that corresponds to the differences between the baseline and monitor datasets) is usually small and can easily blend in with noise or problems related to the acquisition.

These limitations can be overcome by using Least-Squares Migration (LSM), which is a more robust inversion method based on minimizing the misfit between the observed data d and the synthetic data generated by the inverted reflectivity [1]. This method has been proposed to seek an inverted image, which generates simulated data that best matches the recorded seismic data d. This method effectively improves the resolution of the seismic image, compensates for inadequacies for the recorded seismic data due to acquisition, compensates for image illumination, and attenuates imaging artefacts.

LSM is traditionally solved by iteration, where each iteration typically includes: a migration, a de-migration (modelling) and a comparison to the real seismic data being imaged. More specifically, a misfit C(m) between the observed data d and the synthetic data Lm generated by the inverted reflectivity is given by:

$\begin{array}{cc}C\left(m\right)=\frac{1}{2}{d-\mathrm{Lm}}^2,& \left(1\right)\end{array}$

where d is the data, m is the reflectivity, and L the Born-modeling or de-migration operator. Minimization of this cost function results in the reflectivity m:

$\begin{array}{cc}m={\left(L^{T}L\right)}^{-1}{L}^{T}d,& \left(2\right)\end{array}$

where L^T is the migration operator and L^TL is the Hessian matrix.

In real scenarios, directly inverting the Hessian matrix L^TL is computationally unfeasible and some approximation must be performed. Iterative schemes [2-4] can be used to invert the Hessian, but they suffer from high computational cost for common 3D surveys. Alternative solutions have been proposed to make the process more feasible, for example the work described in [5-7].

However, all these algorithms have their own limitations. For example, the algorithm described in [7] uses a single iteration approach, approximating the inverse Hessian with non-stationary filters. Although it is cost-efficient, this method may leave enough imprint of the acquisition. In the context of 4D imaging, if there is not repeatability between the baseline and monitor acquisition geometries, this can destroy the 4D signal. The other known algorithms have their own limitations that are known in the art and thus, not discussed herein.

Thus, there is a need for a new process that handles the recorded seismic data (especially the 4D signal) in a more efficient way, reduces a time spent by the computing device for generating the image of the subsurface, and overcomes the limitations noted above with the existing methods. Accordingly, it would be desirable to provide systems and methods with such capabilities.

SUMMARY OF THE INVENTION

According to an embodiment, there is a least-square migration, LSM, based method for generating a 4D image of a subsurface and the method includes receiving seismic data d related to the subsurface, where the seismic data d includes a baseline dataset d_B and a monitor dataset d_M, with the monitor dataset d_M being taken later in time than the baseline dataset d_B, for the same subsurface, calculating a baseline filter _B and a monitor filter _M based on a same common reflectivity r of the subsurface and corresponding remigrated baseline data m_B1 and remigrated monitor data m_M1 so that the base filter _B applied to the remigrated baseline data m_B1 equals the monitor filter _M applied to the remigrated monitor data m_M1, applying the baseline filter _B to raw migrated baseline data m_B0 and applying the monitor filter _M to raw migrated monitor data m_M0 to generate LSM baseline data m_b and LSM monitor data m_M, and generating the 4D image of the subsurface based on the LSM baseline data m_B and the LSM monitor data m_M.

According to another embodiment, there is a computing device for generating a 4D image of a subsurface based on a least-square migration, LSM, based method, and the computing device includes an interface for receiving seismic data d related to the subsurface, wherein the seismic data d includes a baseline dataset d_B and a monitor dataset d_M, with the monitor dataset d_M being taken later in time then the baseline dataset d_B, for the same subsurface and a processor connected to the interface. The processor is configured to calculate a baseline filter _B and a monitor filter _M based on a same common reflectivity r of the subsurface and corresponding remigrated baseline data m_B1 and remigrated monitor data m_M1 so that the base filter _B applied to the remigrated baseline data m_B1 equals the monitor filter _M applied to the remigrated monitor data m_M1, apply the baseline filter _B to raw migrated baseline data m_B0 and applying the monitor filter _M to raw migrated monitor data m_M0 to generate LSM baseline data m_B and LSM monitor data m_M, and generate the 4D image of the subsurface based on the LSM baseline data mg and the LSM monitor data m_M.

According to yet another embodiment, there is a least-square migration, LSM, based method for generating an image of a subsurface, and the method includes receiving seismic data d related to the subsurface, wherein the seismic data d includes a baseline dataset d_B and a monitor dataset d_M, with the monitor dataset d_Mbeing taken later in time then the baseline dataset d_B, for the same subsurface, iteratively updating a least square migration baseline data m_B and a least square migration monitor data m_M based on descent directions of (1) a first baseline residual between baseline migrated seismic data L^T_Bd_B of the baseline dataset d_B and baseline remigrated seismic data L^T_BL_Bm_B of the least square migration baseline data m_B , and (2) a second monitor residual between monitor migrated seismic data L^T_Md_Mof the monitor dataset d_Mand monitor remigrated seismic data L^T_ML_Mm_M of the least square migration monitor data m_M, migrating the seismic data d with the updated least square migration baseline data m_B and the updated least square migration monitor data m_M, and generating the 4D image of the subsurface based on the migrated seismic data d with the updated least square migration baseline data m_B and the updated least square migration monitor data m_M. The least square migration baseline data m_B and the least square migration monitor data m_M are preconditioned with a baseline filter _B and a monitor filter _M, respectively, which are calculated based on a same common reflectivity r of the subsurface and corresponding remigrated baseline data and remigrated monitor data.

BRIEF DESCRIPTION OF THE DRAWINGS

For a more complete understanding of the present invention, reference is now made to the following descriptions taken in conjunction with the accompanying drawings, in which:

FIG. 1 is a flow chart of a method for generating an image of a subsurface based on single iteration LSM using distinct filters for a baseline dataset and a monitor dataset;

FIG. 2 schematically illustrates a geological formation located into a subsurface and how a source probes the formation with seismic waves;

FIGS. 3A and 3B illustrate seismic data corresponding to two different vintages for a same subsurface;

FIG. 4A illustrates the image differences between a traditional imaging method and a traditional LSM method, FIG. 4B illustrates the image difference between the traditional imaging method and the constrained LSM method, and FIG. 4C shows the image differences between the traditional LSM method and the constrained LSM method;

FIG. 5 is another flow chart of a method for generating an image of the subsurface based on single iteration LSM using filters for both a baseline dataset and a monitor dataset;

FIG. 6 is a flow chart of a method for generating an image of the subsurface based on plural iterations LSM using filters for preconditioning descent directions used to update a velocity model; and

FIG. 7 is a schematic diagram of a data computing system for generating an image of the subsurface based on LSM using filters.

DETAILED DESCRIPTION OF THE INVENTION

The following description of the embodiments refers to the accompanying drawings. The same reference numbers in different drawings identify the same or similar elements. The following detailed description does not limit the invention. Instead, the scope of the invention is defined by the appended claims. The following embodiments are discussed, for simplicity, with regard to a single iteration time-lapse LSM method for processing the 4D data for generating an image of the subsurface. However, the embodiments to be discussed next are not limited to a single iteration LSM or to generating an image of the subsurface, but may be applied to plural iteration methods and/or for generating other data from the seismic data.

The methods discussed herein may be applied to the field of subsurface exploration, for example, hydrocarbon exploration and development, geothermal exploration and development, and carbon capture and sequestration, or other natural resource exploration and exploitation. They could also be employed for surveying and monitoring for windfarm applications, both onshore and offshore, and also for medical imaging applications.

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

To circumvent the problems associated with the existing LSM methods, single iteration methods based on approximating the inverse Hessian with non-stationary filters have been proposed [7, 8]. The LSM method, as solved by [7], is posed in a least squares sense by minimizing a cost function C(m), which is given by the following equation:

$\begin{array}{cc}C\left(m\right)={d-\mathrm{Lm}}^2,& \left(3\right)\end{array}$

where d is the recorded seismic data, L is the demigration ‘modelling’ operator and m is the reflectivity, and Lm is the simulated data. In this problem, d and L are known and m is calculated. The least squares solution to this problem is given by equation (2) discussed above.

The main difficulty in solving equation (2) is in the evaluation of the operator (L^TL)⁻¹, which is the inverse of the Hessian. The second part of equation (2) is the raw migrated image m₀=L^Td.

The author in [4] proposed to estimate the operator (L^TL)⁻¹ (the inverse of the Hessian) by creating another image realization m₁ by applying the de-migration operator L and the re-migration operator L^T, that is:

$\begin{array}{cc}{m}_1=\left(L^{T}L\right)L^{T}d.& \left(4\right)\end{array}$

In other words, this method calculates the re-migration m₁=L^TLm₀ by de-migrating m₀ with operator L and then migrating the result of the de-migration step with operator L^T, thus relating to the raw migration m₀ through the Hessian matrix L^TL. A matching filter between the two image realizations m₀ and m₁ is calculated using the following equation:

$\begin{array}{cc}D\left(ℱ\right)=\arg {\min }_{ℱ}m_0-ℱ*m_1,& \left(5\right)\end{array}$

where D() is a cost function to be minimized, is the matching filter (it has many components) derived by matching the two realizations m₀ and m₁, and ‘*’ denotes a convolution.

The filter then approximates the inverse of the Hessian operator (L^TL)⁻¹ i.e.,

$\begin{array}{cc}ℱ\left(m_1\right)\approx m_0,\mathrm{or}ℱ\approx {\left(L^{T}L\right)}^{-1}.& \left(6\right)\end{array}$

Filter F is then applied to the raw reflectivity m₀ , which results in the following approximation of equation (2), i.e., the LSM is approximated by:

$\begin{array}{cc}m={\left(L^{T}L\right)}^{-1}{L}^{T}d\approx ℱ*L^{T}d=ℱ\left(m_0\right).& \left(8\right)\end{array}$

As discussed above, this approach may leave enough imprint of the acquisition and, in the context of 4D imaging, if there is not repeatability between the baseline and monitor acquisition geometries, this can destroy the 4D signal.

According to an embodiment, a method that achieves repeatability between the baseline and monitor acquisitions, which is based on the non-stationary filter F is now discussed. Subscripts B and M are introduced to label the baseline and monitor related data, filters, and operators, respectively. The LSM approach discussed above, when applied for each survey dataset, satisfies the relationships:

$\begin{array}{cc}{m}_B={\left(L_B^T{L}_B\right)}^{-1}{m}_{B0},m_M={\left(L_M^T{L}_M\right)}^{-1}{m}_{M0},\mathrm{where}& \left(9\right)\end{array}$ $\begin{array}{cc}{m}_{B0}=L_B^T{d}_B,m_{M0}=L_M^T{d}_M.& \left(10\right)\end{array}$

For computing the re-migrations m_B1 and m_M1, instead of de-migrating each raw migration m_B0 and m_M0, according to this embodiment, a common reflectivity r is used. This common reflectivity r acts as a constraint on the LSM method [11]. The common reflectivity r can be estimated/selected based on the experience of the operator running the simulation for best describing the subsurface of interest. In one application, the common reflectivity r is calculated based on the data from a previous seismic survey. In yet another application, the common reflectivity r is calculated, based on other methods, either from the baseline data or the monitor dataset. Any known method for calculating the reflectivity may be used.

In this case, the re-migrations are given by:

$\begin{array}{cc}{m}_{B1}=L_B^T{L}_{B}r,m_{M1}=L_M^T{L}_{M}r.& \left(11\right)\end{array}$

Similar to the independent LSM method discussed above with regard to equations (1-8), filters _B and _M, for each of the baseline and monitor datasets are calculated. The filters _B and _M are constructed so that each of them satisfies a corresponding relationship:

$\begin{array}{cc}{ℱ}_B\left(m_{B1}\right)\approx r,{ℱ}_M\left(m_{M1}\right)\approx r,& \left(12\right)\end{array}$

such that each filter approximates the respective inverse Hessian matrix. These filters can be computed in any domain, by applying a linear transform on (1) the reflectivity and (2) the re-migrations before the filter computation.

For example, in one application, a curvelet transformation , which is a linear transform from the image domain is used. Other transformations may be used, as discussed later. This transform is named herein a “filter transform” because the filter is calculated in this domain, which is different from the image domain and the data domain. The corresponding domain may be called “filter domain.” Following the approach in [8], the image is decomposed in the curvelet domain (based on [9]), where the seismic events are decomposed in wavenumber scales and dips, allowing for better illumination compensation for different components of the data. In this domain (i.e., a curvelet domain), a filter f is sought that minimizes the cost-function

$\begin{array}{cc}{\left(m_0\right)-f\left(m_1\right)}^2+{\alpha }^2{f}^2.& \left(13\right)\end{array}$

A direct solution to this equation is given by:

$\begin{array}{cc}f=\frac{\left(m_0\right){\left(m_1\right)}^{*}}{{❘\left(m_1\right)❘}^2+{\alpha }^2}.& \left(14\right)\end{array}$

Then, the solution for the single iteration LSM method is given by:

$\begin{array}{cc}ℱ\left(m_0\right)={𝒞}^{-1}\left(\mathrm{𝒻𝒞}\left(m_0\right)\right).& \left(15\right)\end{array}$

Note that equation (15) indicates that the transform effectively transforms the raw migration data m₀ to the curvelet domain, performs the filtering f in that domain, and then transforms the results back into the image domain with the inverse curvelet transform ⁻¹. While this embodiment transformed the data (raw migration data) from the image domain to the curvelet domain, it is possible to use other domains for applying the filter. In this regard, note that the term “image domain” is used in this context to refer to the transformed representation of subsurface structures and geological features that result from the application of migration algorithms L^T to the seismic data d. These algorithms aim to create accurate and interpretable images of the subsurface. Some examples of image domains used in seismic imaging, that are appropriate for this embodiment include:

Time-Migrated Image Domain: This is a common type of image domain where the seismic data are processed in the time domain. It involves positioning the seismic reflections at their correct time positions in the subsurface, creating an image that represents the depth at which geological features are located. Time-migrated images are widely used for structural interpretation and subsurface mapping.

Depth-Migrated Image Domain: Depth migration transforms seismic data from the time domain to the depth domain. It accurately positions seismic reflections at their correct spatial depths in the subsurface. Depth-migrated images are valuable for understanding the three-dimensional geometry of subsurface structures.

Common-Offset Gather Image Domain: In some cases, pre-stack gathers (common-offset gathers) are transformed into an image domain. This can involve stacking seismic traces at each offset to enhance signal-to-noise ratio and improve resolution.

Common-Azimuth (Angle) Gather Image Domain: In areas with complex subsurface geometries, images can be transformed into common-azimuth (angle) gather domains. This type of imaging helps reveal the directional variation of subsurface features and is important for understanding anisotropy effects.

Amplitude-Preserved Image Domain: Some migration algorithms aim to preserve the amplitude information of seismic reflections, resulting in images that accurately represent the subsurface reflectivity and potential reservoir properties.

Anisotropic Image Domain: Anisotropic migration takes into account the directional dependence of seismic wave propagation in anisotropic media. The resulting images account for the varying velocities and complexities introduced by anisotropy.

Amplitude-Versus-Offset (AVO) Image Domain: AVO inversion and imaging techniques transform seismic data to highlight the variations in seismic amplitudes with offset, providing insights into subsurface fluid content, lithology, and other reservoir properties.

Multi-Dimensional Image Domain: Seismic images can be presented in multi-dimensional forms, incorporating attributes beyond just amplitude and time or depth. These attributes might include attributes derived from seismic inversion or other advanced processing techniques.

Each of these image domains serves a specific purpose and can provide unique insights into the subsurface. The choice of image domain depends on the geological setting, the goals of the study, and the specific challenges posed by the seismic data and subsurface conditions.

The curvelet transform noted with regard to equations (13-15) is a mathematical technique used in seismic imaging, and image processing in general, where data is often characterized by curves and edges at various scales and orientations. It is an extension of the wavelet transform, designed to capture highly directional and anisotropic features in data. In the context of seismic imaging, the curvelet transform is employed to analyze seismic data in a way that efficiently represents features of subsurface structures. Seismic data often contains information about subsurface layers, faults, and other geological features that have different orientations and scales. The curvelet transform excels at capturing these features by using a multiscale and multidirectional decomposition approach.

While other transforms may be used, for example, Fourier transform, ridgelet, contourlet, seislet, linear Radon, parabolic Radon, hyperbolic Radon, shifted-hyperbolic Radon, wavelet, complex wavelet, etc., some of which are described in Yilmaz, Öz (2001), “Seismic data analysis,” Society of Exploration Geophysicists, ISBN 1-56080-094-1, the curvelet transform has the following advantages.

Multiscale Decomposition: Similar to the wavelet transform, the curvelet transform decomposes data into different scales (frequencies). However, it goes beyond this by also considering different directions. This allows the transform to capture features that are not aligned with the standard horizontal and vertical orientations.

Directional Sensitivity: The curvelet transform is particularly effective at capturing features with directional information, which is desired in seismic imaging. Geological structures often exhibit directional patterns, such as faults, fractures, and layer boundaries. The curvelet transform can represent these features more accurately than traditional transforms.

Sparse Representation: The curvelet transform produces a sparse representation of data, meaning that only a few coefficients are significant while the rest are close to zero. This can greatly reduce the amount of data needed to represent complex features accurately, leading to more efficient storage and processing.

Anisotropic Features: The curvelet transform is well-suited for detecting and representing anisotropic features, which are features that exhibit different properties in different directions. Seismic data often contains such anisotropic features due to the varied orientations of subsurface structures.

Returning to equation (12), the two filters _B and _M are required to satisfy the relation:

$\begin{array}{cc}{ℱ}_B\left(m_{B1}\right)={ℱ}_M\left(m_{M1}\right),& \left(16\right)\end{array}$

which ensures consistency between both filters. With this constraint, the null space of the baseline Hessian L^T_BL_B is captured in the monitor filter _M and vice-versa. The null space (also known as the kernel) of a matrix or a linear transformation is the set of all vectors in the domain that are mapped to the zero vector in the codomain. In other words, it is the set of solutions to the homogeneous equation Ax=0, where A is a given matrix and x is a vector of variables.

Thus, the resulting 4D LSM migrated data is described by:

$\begin{array}{cc}{m}_B\approx {ℱ}_B\left(m_{B0}\right),m_M\approx {ℱ}_M\left(m_{m0}\right)& \left(17\right)\end{array}$

The two filters _B and _M are joint filters because they satisfy relationship (16), which relates the two filters to each other. A method for calculating the joint filters is now illustrated. Because of equation (16), the method seeks to minimize the following cost function for two filters f_B and f_M (in the curvelet space in this instance):

$\begin{array}{cc}{𝒞\left(r\right)-{𝒻}_B𝒞\left(m_{1B}\right)}^2+{\alpha }_B^2{{𝒻}_B}^2+{𝒞\left(r\right)-f_M𝒞\left(m_{1M}\right)}^2+{\alpha }_M^2{f_M}^2-\left(\langle \lambda ,f_B𝒞\left(m_{1M}\right)-f_M𝒞\left(m_{1M}\right)\rangle +c.c\right).& \left(18\right)\end{array}$

The first two terms in equation (18) are the independent minimization of the baseline and monitor filters, respectively, corresponding to equation (15). The third term in equation (18) is coupling the constraint of equation (16) imposed by the Lagrange multiplier λ. The solution to equation (18) is given by:

$\begin{array}{cc}{f}_B=\frac{2{❘𝒞\left(m_{1M}\right)❘}^2{𝒞\left(m_{1B}\right)}^{*}𝒞\left(r\right)}{2{❘𝒞\left(m_{1B}\right)❘}^2{❘𝒞\left(m_{1M}\right)❘}^2+{\alpha }_B^2{❘𝒞\left(m_{1M}\right)❘}^2+{\alpha }_M^2{❘𝒞\left(m_{1B}\right)❘}^2},f_M=\frac{{❘𝒞\left(m_{1M}\right)❘}^2{𝒞\left(m_{1B}\right)}^{*}𝒞\left(r\right)}{2{❘𝒞\left(m_{1B}\right)❘}^2{❘𝒞\left(m_{1M}\right)❘}^2+{\alpha }_B^2{❘𝒞\left(m_{1M}\right)❘}^2+{\alpha }_M^2{❘𝒞\left(m_{1B}\right)❘}^2}.& \left(19\right)\end{array}$

The above solutions can be contrasted with those obtained through independent evaluations of the filters, f_B⁰ and f_M⁰, as in equation (14):

$\begin{array}{cc}{f}_B^0=\frac{𝒞\left(r\right){𝒞\left(m_{1B}\right)}^{*}}{{❘𝒞\left(m_{1B}\right)❘}^2+{\alpha }_B^2},f_B^0=\frac{𝒞\left(r\right){𝒞\left(m_{1M}\right)}^{*}}{{❘𝒞\left(m_{1M}\right)❘}^2+{\alpha }_M^2}.& \left(20\right)\end{array}$

It is worth noticing a relation between these solutions:

$\begin{array}{cc}{f}_B=\frac{\mathrm{HM}\left(f_B^0𝒞\left(m_{1B}\right),f_M^0𝒞\left(m_{1M}\right)\right)}{𝒞\left(m_{1B}\right)},f_M=\frac{\mathrm{HM}\left(f_B^0𝒞\left(m_{1B}\right),f_M^0𝒞\left(m_{1M}\right)\right)}{𝒞\left(m_{1M}\right)},& \left(21\right)\end{array}$

where HM(x₁, x₂) is the harmonic mean of x₁ and x₂. This leads to a generic formula for obtained the joint filters based on any generalized mean. In particular, for the Lehmer mean L_p [10]:

$\begin{array}{cc}{f}_B=\frac{L_p\left(f_B^0𝒞\left(m_{1B}\right),f_M^0𝒞\left(m_{1M}\right)\right)}{𝒞\left(m_{1B}\right)},f_M=\frac{L_p\left(f_B^0𝒞\left(m_{1B}\right),f_M^0𝒞\left(m_{1M}\right)\right)}{𝒞\left(m_{1M}\right)}.& \left(22\right)\end{array}$

Here p is a free parameter. Special cases are p=0 for the harmonic mean, p=0,5 for the geometric mean and p=1 for the arithmetic mean.

Then, the solution for the 4D single iteration LSM is

$\begin{array}{cc}{ℱ}_B\left(m_{0B}\right)={𝒞}^{-1}\left(f_B𝒞\left(m_{0B}\right)\right),{ℱ}_M\left(m_{0M}\right)={𝒞}^{-1}\left(f_M𝒞\left(m_{0M}\right)\right),& \left(23\right)\end{array}$

which is expressed in the image domain due to the inverse curvelet transform ⁻¹.

Having the LSM data, it is then possible to generate an image of the subsurface, to visualize changes in the reservoir associated with the baseline and monitor datasets, i.e., the 4D data. In one embodiment, the output of the migration process is the subsurface image that represents the geological structures, faults, and other features. This image is typically displayed in a 2D or 3D format, where the vertical axis represents depth or time, and the horizontal axes correspond to spatial coordinates. This image is used by those skilled in the art of image interpretation, for example, geophysicists and/or geologists, to identify potential subsurface reservoirs, geological structures, and other features that might indicate the presence of hydrocarbons, minerals, or other valuable resources.

The method discussed herein may be applied in any pre-stack domain common to both baseline and monitor datasets. In seismic analysis, the pre-stack domain refers to the stage of data processing and analysis where seismic data is manipulated and interpreted before it is stacked. This domain retains the individual seismic traces acquired at various source-receiver offsets and angles, allowing for more detailed and accurate subsurface imaging and interpretation. Working in the pre-stack domain is particularly useful in areas with complex geological structures, anisotropic media, or when performing advanced analysis techniques.

For a pre-stack application, let i label an instance of such dataset in a given domain. The joint filters are computed separately for each i, but all using the same reflectivity r:

$\begin{array}{cc}{ℱ}_B^i\left(m_{B1}^i\right)\approx r,{ℱ}_M^i\left(m_{M1}^i\right)\approx r,{ℱ}_B^i\left(m_{B1}^i\right)={ℱ}_M^i\left(m_{M1}^i\right),& \left(24\right)\end{array}$

where the re-migrations are given by:

$\begin{array}{cc}{m}_{B1}^i=L_B^{\mathrm{iT}}{L}_B^{i}r,m_{M1}^i=L_M^{\mathrm{iT}}{L}_M^{i}r.& \left(25\right)\end{array}$

The 4D pre-stack LSM are then given by:

$\begin{array}{cc}{m}_B^i\approx {ℱ}_B^i\left(L_B^{\mathrm{iT}}{d}_B\right),m_M^i\approx {ℱ}_M^i\left(L_M^{\mathrm{iT}}{d}_M\right).& \left(26\right)\end{array}$

This approach can also be applied for 3D pre-stack LSM. For example, for a single vintage, all pre-stack remigrations are performed with the same reflectivity r. If there are n instances of the pre-stack domain, it is possible to compute joint filters satisfying:

$\begin{array}{cc}{ℱ}_B^i\left(m_{B1}^i\right)\approx r,{ℱ}_B^i\left(m_{B1}^i\right)={ℱ}_B^j\left(m_{B1}^j\right),i,j=1\mathrm{to}n.& \left(27\right)\end{array}$

It is worth noting that in some 4D projects, more than two surveys, i.e., vintages, can be considered. The method discussed above can be extended and it applies naturally to such situations, by generalizing the conditions from equations (16) and (17), i.e.,

$\begin{array}{cc}{ℱ}_i\left(m_{i1}\right)={ℱ}_j\left(m_{j1}\right),i,j=1\mathrm{to}n,& \left(28\right)\end{array}$ $\begin{array}{cc}{m}_i\approx {ℱ}_i\left(m_{i0}\right)& \left(29\right)\end{array}$

where i and j represent different seismic surveys i.e., different 4D vintages, and n is the number of surveys under consideration. Expressions for the filters in the filter domain can be found through the generalized Lehmer mean formulation of equation (22).

The method of single iteration, time-lapse LSM discussed above can be summarized as now discussed with regard to FIGS. 1 and 2. Note that this method does not use multiple iterations, as other methods do, but uses a single iteration for migrating the data using the filters discussed above. The method starts in step 100 with receiving seismic data d_B associated with a first seismic survey, the baseline survey, and a step 102 of receiving seismic data d_Massociated with a second seismic survey, the monitor survey, as illustrated in FIG. 1. The two seismic surveys acquire seismic data over a same subsurface 210, as illustrated in FIG. 2, but at different times, for monitoring a reservoir or formation 202. FIG. 2 shows a land seismic survey in which a source S emits seismic energy into the subsurface. A wavefield 212 is shown propagating toward the reservoir 202 and then being reflected, and the reflected wavefield 214 is recorded by a corresponding seismic sensor R1. The source S may be any type of seismic source, examples include airgun, pinger, sparker, boomer, marine vibrator, land vibrator, dynamite, etc. The source may be a single source (e.g., single airgun) or an array of sources (e.g., airgun array). The seismic sensor may be one of hydrophones, geophones, particle motion sensors, particle velocity sensors, accelerometers, near-field hydrophones, near-field accelerometers or other sensor configured to detect seismic energy. The data may be from towed streamer, ocean-bottom sensor, OBS, (node or cable) acquisition, land acquisition, transition zone campaign, or borehole (e.g., vertical seismic profile (VSP), distributed acoustic sensing (DAS)).

FIG. 2 further shows another wavefield 216 that is reflected from an interface 220 between layers having different impedances, and the reflected wavefield 218 is recorded by another sensor R2, or by the first sensor R1. The seismic acquisition system 200 is shown in FIG. 2 being a land system, located on the surface 201 of the earth. However, the system 200 may be a marine system or any other system used in the field. The system 200 may include plural sources S and plural sensors Ri, with i taking a value up to tens of thousands.

The seismic data d_B and d_M are then migrated with the migration operator L′ and L′M, respectively, in steps 104 and 106, for obtaining the raw migrated data m_B0 and m_M0, respectively. In step 108, the reflectivity r of the subsurface is obtained, as previously discussed, either from a preexisting survey, or from being estimated by the operator of the method, based on prior knowledge, and the details of the subsurface. In step 110, the reflectivity r is demigrated (i.e., operation L_Br is performed) for the baseline dataset, i.e., for the baseline acquisition geometry, based on equation (11), and in step 112 the results of step 110 are re-migrated (i.e., operator L^T_B is applied). This means that steps 110 and 112 are equivalent to operation L^T_BL_Br, which results in calculating the remigration data m_B1 for the baseline dataset. Similar steps are performed for the monitor dataset, i.e., in step 114, the same reflectivity r as in step 110 is demigrated (i.e., operation L_Mr is performed) for the monitor dataset (for the baseline acquisition geometry), based on equation (11), and in step 116 the results of step 114 are re-migrated (i.e., operator L^T_M is applied). This means that steps 114 and 116 are equivalent to operation L^T_ML_Mr, which results in calculating the remigration data m_M1 for the monitor dataset. Note that steps 110 and 112 and steps 114 and 116 are performed independent of each other. This means that these steps may be performed in parallel or in series.

As various domains are used for processing the data for the method illustrated in FIG. 1, it is noted that steps 100 and 102 are performed in the data domain (e.g., space-time domain), steps 104, 106 and 108 are performed in the image domain, steps 110 and 114 are performed in the data domain, and steps 112 and 116 are performed in the image domain.

Next, the method transforms in step 118 the reflectivity r (from step 108) and the remigration data m_B1 and m_M1 from steps 112 and 116 into the filter domain (for example, a curvelet domain using a curvelet transform, as described by equation (18)). Then, the joint filters _B and _M are calculated in step 120, in the filter domain, based on equations (19) and (20). In step 122, the raw migrations from steps 104 and 106 are transformed into the filter domain and the calculated filters _B and _M are applied to calculate LSM migration data m_B≈F_B(m_B0) and m_M≈_M(m_M0). In step 124, the baseline LSM migration data m_B is transformed back to the image domain and in step 126, the monitor LSM migration data m_M is also transformed back to the image domain. Based on the LSM migration data m_B and m_M, the method generates in step 128 an image of the subsurface, which might identify the reservoir 202, so that the user of the image can decide, if necessary, where to drill a new well, where to focus the production effort, etc.

The term “image” has a broader meaning than a two-dimensional grid of pixel values, being rather a 3-dimensional map of one or more attribute values. The most frequently mapped attribute is seismic wave propagation velocity. The attribute values characterize (and thus allow identification of) different material volumes inside the subsurface and therefore enable estimation of the location and shape of interfaces between different materials. While this profile/image does not provide an accurate location for natural resources such as (but not limited to) oil and/or gas, it enables those trained in the field to determine the presence or absence thereof.

The method discussed above with regard to FIG. 2 has been tested on the synthetic data illustrated in FIGS. 3A and 3B. FIG. 3A shows the seismic data of a first vintage for a given subsurface and FIG. 3B shows the seismic data for a second vintage for the same subsurface. FIG. 4A illustrates the image differences between the traditional LSM method and the image of the subsurface obtained with a traditional process, and FIG. 4B illustrates the image differences between the constrained LSM of FIG. 2 and the traditional image. FIG. 4C shows the true differences between the images shown in FIGS. 4A and 4B. One skilled in the art would recognize that the image generated by the constraint LSM of FIG. 2 is superior to the traditional LSM method.

An LSM based method with constraint for generating an image of a subsurface is schematically illustrated in FIG. 5 and includes a step 500 of receiving seismic data d related to the subsurface, wherein the seismic data d includes a baseline dataset d_B and a monitor dataset d_M, with the monitor dataset d_Mbeing taken later in time then the baseline dataset d_B, for the same subsurface, a step 502 of calculating a baseline filter _B and a monitor filter _M based on a same common reflectivity r of the subsurface and corresponding remigrated baseline data m_B1 and remigrated monitor data m_M1, a step 504 of applying the baseline filter _B to raw migrated baseline data m_B0 and applying the monitor filter _M to raw migrated monitor data m_M0 to generate LSM baseline data m_B and LSM monitor data m_M, and a step 506 of generating the image of the subsurface based on the LSM baseline data m_B and the LSM monitor data m_M.

The step of calculating may include performing raw migration, with a baseline migration operator L^T_B on the baseline dataset d_B and with a monitor migration operator L^T_M on the monitor dataset d_M, to calculate the monitor raw migrated data m_B0 and the baseline raw migrated data m_M0. This step may further include calculating the remigrated baseline data m_B1 by applying a baseline demigration operator L_B on the reflectivity r, followed by applying the baseline migration operator L^T_B so that the remigrated baseline data m_B1=L^T_BL_Br, and calculating the remigrated monitor data m_M1 by applying a monitor demigration operator L_M on the reflectivity r, followed by applying the monitor migration operator L^T_M so that the remigrated monitor data m_M1=L^T_ML_Mr. In one application, the step may further include transforming the reflectivity r with a filter transform , from an image domain to a filter domain, transforming the remigrated baseline data m_B1 with the filter transform , from the image domain to the filter domain, and transforming the remigrated monitor data m_M1 with the filter transform , from the image domain to the filter domain. In addition, the step of calculating may also include computing a baseline filter f_B in the filter domain based on a filter domain representation (r) of the reflectivity r, a filter domain representation (m_B1) of the remigrated baseline data m_B1, and a filter domain representation (m_M1) of the remigrated monitor data m_M1, and computing a monitor filter f_M in the filter domain based on the filter domain representation (r) of the reflectivity r, the filter domain representation (m_B1) of the remigrated baseline data m_B1, and the filter domain representation (m_M1) of the remigrated monitor data m_M1.

In one application, the filter transform is a curvelet transform and the filter domain is a curvelet domain. The method may further include subjecting the baseline filter f_B in the filter domain and the monitor filter f_M in the filter domain to produce a same result when applied to the remigrated baseline data m_B1 and to remigrated monitor data m_M1, respectively, and/or transforming (1) the raw migrated baseline data m_B0 to transformed raw migrated baseline data (m_B0) and (2) the raw migrated monitor data m_M0 to transformed raw migrated monitor data (m_M0) in the filter domain with the filter transform . The step of applying may include applying the baseline filter f_B in the filter domain, to the transformed raw migrated baseline data (m_B0), and applying the monitor filter f_M in the filter domain, to the transformed raw migrated monitor data (m_M0), to generate the LSM baseline data f_B(m_B ) in the filter domain and the LSM monitor data f_M(m_M) in the filter domain, and transforming back, with an inverse transform ⁻¹, (1) the LSM baseline data f_B(m_B ) to obtain the LSM baseline data m_B=⁻¹(f_B(m_B0)) and (2) the LSM monitor data f_M(m_M) to obtain the LSM monitor data m_M=⁻¹(f_M(m_M0)). The inverse transform ⁻¹ takes a quantity from the curvelet domain to the image domain.

The above method is performed by doing a single iteration LSM approach. However, it is possible to use multiple iterations for the LSM method. If multiple iterations are used, then the steps of determining the filters _B and _M may be adapted to be used as a preconditioning for iterative LSM. For example, at the k-th iteration of the method, the descent directions

$\begin{array}{cc}{g}_B^{\left(k\right)}=L_B^T\left(d_B-L_B{m}_B^{\left(k\right)}\right),g_M^{\left(k\right)}=L_M^T\left(d_M-L_M{m}_M^{\left(k\right)}\right)& \left(30\right)\end{array}$

can be preconditioned by the joint filters _B and _M before computing the updates for the reflectivities.

The iterative LSM aims to improve the accuracy of seismic images by iteratively updating the migrated image to minimize the difference between the observed seismic data and the modeled data. According to this approach, and as illustrated in FIG. 6, the method receives in step 600 the seismic data for data preprocessing. Similar to steps 100 and 102, the method may receive the data from baseline and monitor surveys. Preprocessing the raw seismic data may include shot gathers and receiver gathers, removing noise, correcting for static corrections, and performing data conditioning. In step 602, the method builds the velocity model. A preliminary velocity model of the subsurface is either inputted by the user, based on the knowledge of the surveyed subsurface, e.g., using well logs, seismic tomography, or other available geological and geophysical information, or is randomly initiated. In step 604 an initial migration is performed using the initial velocity model to create an initial subsurface image. Next, in step 606, forward modeling is performed, i.e., using the initial migrated image and the initial velocity model, and the method computes synthetic seismic data (modeled data) through the process of forward modeling. This involves simulating how seismic waves would propagate through the subsurface based on the given velocity model. In step 408, the method computes the residuals, i.e., calculates the residuals by subtracting the observed seismic data from the modeled data. Residuals represent the differences between the actual recorded data and the data predicted by the initial migration. In this step, the method may also calculate the gradient of the objective function (misfit) with respect to the velocity model parameters. The gradient points in the direction of the steepest increase in the objective function. In the context of optimization, taking the negative of the gradient gives the direction of steepest decrease (descent direction).

In step 610, the velocity model is iteratively updated, in the direction of the descent direction (see equations (30)), using, for example, inversion techniques, to minimize the misfit between the observed data and the modeled data. By iteratively following the descent direction provided by the negative gradient, the optimization algorithm seeks to find the subsurface velocity model that best fits the observed seismic data. Convergence occurs when the algorithm reaches a point where further adjustments in the velocity model do not lead to significant reductions in the misfit. This step involves solving an optimization problem to update the velocity model parameters. Steps 604 to 610 are then reiterated by performing migration using the updated velocity model, generating new modeled data, computing residuals, and updating the velocity model until convergence is achieved or a predefined stopping criterion is met in step 612. Once convergence is achieved or the iterative process is completed in step 612, the method performs a final migration step 614 using the refined velocity model, to generate in step 616 an improved subsurface image. In an optional step 618, the method may apply additional post-processing steps to the final migrated image, such as amplitude scaling, noise suppression, and edge enhancement, to improve the interpretability of the subsurface features. In one application, instead of updating the velocity model, the least square migration baseline data mg and the updated least square migration monitor data m_M are updated in the above steps.

The above-discussed procedures and methods may be implemented in a computing device as illustrated in FIG. 7. Hardware, firmware, software or a combination thereof may be used to perform the various steps and operations described herein. The computing device 700 is suitable for performing the activities described in the above embodiments and may include a server 701. Such a server 701 may include a central processor (CPU) 702 coupled to a random access memory (RAM) 704 and to a read-only memory (ROM) 706. ROM 706 may also be other types of storage media to store programs, such as programmable ROM (PROM), erasable PROM (EPROM), etc. Processor 702 may communicate with other internal and external components through input/output (I/O) circuitry 708 and bussing 710 to provide control signals and the like. Processor 702 carries out a variety of functions as are known in the art, as dictated by software and/or firmware instructions.

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

Server 701 may be coupled to other devices, such as seismic sources, seismic sensors, or any other data imaging systems. The server may be part of a larger network configuration as in a global area network (GAN) such as the Internet 728, which allows ultimate connection to various landline and/or mobile computing devices.

As described above, the apparatus 700 may be embodied by a computing device. However, in some embodiments, the apparatus may be embodied as a chip or chip set. In other words, the apparatus may comprise one or more physical packages (e.g., chips) including materials, components and/or wires on a structural assembly (e.g., a baseboard). The structural assembly may provide physical strength, conservation of size, and/or limitation of electrical interaction for component circuitry included thereon. The apparatus may therefore, in some cases, be configured to implement an embodiment of the present invention on a single chip or as a single “system on a chip.” As such, in some cases, a chip or chipset may constitute means for performing one or more operations for providing the functionalities described herein.

The processor 702 may be embodied in a number of different ways. For example, the processor may be embodied as one or more of various hardware processing means such as a coprocessor, a microprocessor, a controller, a digital signal processor (DSP), a processing element with or without an accompanying DSP, or various other processing circuitry including integrated circuits such as, for example, an ASIC (application specific integrated circuit), an FPGA (field programmable gate array), a microcontroller unit (MCU), a hardware accelerator, a special-purpose computer chip, or the like. As such, in some embodiments, the processor may include one or more processing cores configured to perform independently. A multi-core processor may enable multiprocessing within a single physical package. Additionally or alternatively, the processor may include one or more processors configured in tandem via the bus to enable independent execution of instructions, pipelining and/or multithreading.

In an example embodiment, the processor 702 may be configured to execute instructions stored in the memory device 704 or otherwise accessible to the processor. Alternatively or additionally, the processor may be configured to execute hard coded functionality. As such, whether configured by hardware or software methods, or by a combination thereof, the processor may represent an entity (e.g., physically embodied in circuitry) capable of performing operations according to an embodiment of the present invention while configured accordingly. Thus, for example, when the processor is embodied as an ASIC, FPGA or the like, the processor may be specifically configured hardware for conducting the operations described herein. Alternatively, as another example, when the processor is embodied as an executor of software instructions, the instructions may specifically configure the processor to perform the algorithms and/or operations described herein when the instructions are executed. However, in some cases, the processor may be a processor of a specific device (e.g., a pass-through display or a mobile terminal) configured to employ an embodiment of the present invention by further configuration of the processor by instructions for performing the algorithms and/or operations described herein. The processor may include, among other things, a clock, an arithmetic logic unit (ALU) and logic gates configured to support operation of the processor.

The term “about” is used in this application to mean a variation of up to 20% of the parameter characterized by this term. It will be understood that, although the terms first, second, etc. may be used herein to describe various elements, these elements should not be limited by these terms. These terms are only used to distinguish one element from another. For example, a first object or step could be termed a second object or step, and, similarly, a second object or step could be termed a first object or step, without departing from the scope of the present disclosure. The first object or step, and the second object or step, are both, objects or steps, respectively, but they are not to be considered the same object or step.

The terminology used in the description herein is for the purpose of describing particular embodiments and is not intended to be limiting. As used in this description and the appended claims, the singular forms “a,” “an” and “the” are intended to include the plural forms as well, unless the context clearly indicates otherwise. It will also be understood that the term “and/or” as used herein refers to and encompasses any possible combinations of one or more of the associated listed items. It will be further understood that the terms “includes,” “including,” “comprises” and/or “comprising,” when used in this specification, specify the presence of stated features, integers, steps, operations, elements, and/or components, but do not preclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and/or groups thereof. Further, as used herein, the term “if” may be construed to mean “when” or “upon” or “in response to determining” or “in response to detecting,” depending on the context.

The disclosed embodiments provide an LSM method for addressing time-lapse seismic data by using filters selected for each survey. It should be understood that this description is not intended to limit the invention. On the contrary, the embodiments are intended to cover alternatives, modifications and equivalents, which are included in the spirit and scope of the invention as defined by the appended claims. Further, in the detailed description of the embodiments, numerous specific details are set forth in order to provide a comprehensive understanding of the claimed invention. However, one skilled in the art would understand that various embodiments may be practiced without such specific details.

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

This written description uses examples of the subject matter disclosed to enable any person skilled in the art to practice the same, including making and using any devices or systems and performing any incorporated methods. The patentable scope of the subject matter is defined by the claims, and may include other examples that occur to those skilled in the art. Such other examples are intended to be within the scope of the claims.

REFERENCES

The entire content of all the publications listed herein is incorporated by reference in this patent application.[1] Tarantola, A., 1987, Inverse problem theory: Methods for data fitting and model parameter estimation: Elsevier Science Publishing Company.[2] Schuster, G. T., 1993, Least-squares crosswell migration: 63^rd Annual International Meeting, SEG, Expanded Abstracts, 110-113, doi: 10.1190/segam2012-1425.1.[3] Nemeth, T., C. Wu, and G. T. Schuster, 1999, Least-squares migration of incomplete reflection data: Geophysics, 64, 208-221, doi: 10.1190/1.1444517.[4] Tang, Y., 2008, Wave-equation Hessian by phase encoding: 78th Annual International Meeting, SEG, Expanded Abstracts, 2201-2205.[5] Claerbout, J., and D. Nichols, 1994, Spectral preconditioning, in SEP-82: Stanford Exploration Project, 183-186.[6] Rickett, J., 2003, Illumination-based normalization for wave-equation depth migration: Geophysics, 68, 1371-1379.[7] Guitton, 2004, Amplitude and kinematic corrections of migrated images for non-unitary imaging operators: geophysics, 69, 1017-1024.[8] Ping Wang, Shouting Huang, and Ming Wang, (2017), “Improved subsalt images with least-squares reverse time migration,” Interpretation 5: SN25-SN32.

[9] Candès, E., L. Demanet, D. Donoho, and L. Ying, 2006, Fast Discrete Curvelet Transforms: Multiscale Modeling & Simulation, 5, 861-899.

[10] P. S. Bullen. Handbook of means and their inequalities. Springer, 1987.[11] Zhao Wang, Rongxin Huang, Yi Xuan, Qing Xu, Scott Morton, and Brad Kuntz, (2017), “Improving OBS-streamer 4D imaging by least-squares migration,” SEG Technical Program Expanded Abstracts: 5819-5823.

Claims

1. A least-square migration, LSM, based method for generating a 4D image of a subsurface, the method comprising: receiving seismic data d related to the subsurface, wherein the seismic data d includes a baseline dataset dB and a monitor dataset dM, with the monitor dataset dMbeing taken later in time than the baseline dataset dB, for the same subsurface;calculating a baseline filter B and a monitor filter M based on a same common reflectivity r of the subsurface and corresponding remigrated baseline data mB1 and remigrated monitor data mM1 so that the base filter B applied to the remigrated baseline data mB1 equals the monitor filter M applied to the remigrated monitor data mM1;applying the baseline filter B to raw migrated baseline data mB0 and applying the monitor filter M to raw migrated monitor data mM0 to generate LSM baseline data m and LSM monitor data mM; andgenerating the 4D image of the subsurface based on the LSM baseline data mg and the LSM monitor data mM.

2. The method of claim 1, wherein the step of calculating comprises: performing raw migration, with a baseline migration operator LTB on the baseline dataset dB and with a monitor migration operator LTM on the monitor dataset dM, to calculate the baseline raw migrated data mB0and the monitor raw migrated data mM0.

3. The method of claim 2, further comprising: calculating the remigrated baseline data mB1 by applying a baseline demigration operator LB on the reflectivity r, followed by applying the baseline migration operator LTB so that the remigrated baseline data mB1=LTBLBr; andcalculating the remigrated monitor data mM1 by applying a monitor demigration operator LM on the reflectivity r, followed by applying the monitor migration operator LTM so that the remigrated monitor data mM1=LTMLMr.

4. The method of claim 3, further comprising: transforming the reflectivity r with a filter transform , from an image domain to a filter domain;transforming the remigrated baseline data mB1 with the filter transform , from the image domain to the filter domain; andtransforming the remigrated monitor data mM1 with the filter transform , from the image domain to the filter domain.

5. The method of claim 4, further comprising: computing a baseline filter fB in the filter domain based on a filter domain representation (r) of the reflectivity r, a filter domain representation (mB1) of the remigrated baseline data mB1, and a filter domain representation (mB1) of the remigrated monitor data mM1; andcomputing a monitor filter fM in the filter domain based on the filter domain representation (r) of the reflectivity r, the filter domain representation (mB1) of the remigrated baseline data mB1, and the filter domain representation (mM1) of the remigrated monitor data mM1.

6. The method of claim 4, wherein the filter transform is a curvelet transform and the filter domain is a curvelet domain.

7. The method of claim 5, further comprising: subjecting the baseline filter fB in the filter domain and the monitor filter fM in the filter domain to produce a same result when applied to the remigrated baseline data mB1 and to remigrated monitor data mM1, respectively.

8. The method of claim 7, further comprising: transforming (1) the raw migrated baseline data mB0 to transformed raw migrated baseline data (mB0) and (2) the raw migrated monitor data mM0 to transformed raw migrated monitor data (mM0) in the filter domain with the filter transform ,and the step of applying comprises:applying the baseline filter fB in the filter domain, to the transformed raw migrated baseline data (mB0), and applying the monitor filter fm in the filter domain, to the transformed raw migrated monitor data (mM0), to generate the LSM baseline data fB(mB0) in the filter domain and the LSM monitor data fM(mM0) in the filter domain; andtransforming back, with an inverse transform −1, (1) the LSM baseline data fB(mg) to obtain the LSM baseline data mB=−1(fB(m0 )) and (2) the LSM monitor data fM(mM) to obtain the LSM monitor data mM=−1(fM(mM0)).

9. The method of claim 1, wherein the LSM is performed in a pre-stack domain.

10. The method of claim 1, wherein the LSM method performs a single iteration.

11. A computing device for generating a 4D image of a subsurface based on a least-square migration, LSM, based method, the computing device comprising: an interface for receiving seismic data d related to the subsurface, wherein the seismic data d includes a baseline dataset dB and a monitor dataset dM, with the monitor dataset dMbeing taken later in time then the baseline dataset dB, for the same subsurface; anda processor connected to the interface and configured to,calculate a baseline filter B and a monitor filter M based on a same common reflectivity r of the subsurface and corresponding remigrated baseline data mB1 and remigrated monitor data mM1 so that the base filter B applied to the remigrated baseline data mB1 equals the monitor filter M applied to the remigrated monitor data mM1;apply the baseline filter B to raw migrated baseline data mB0 and applying the monitor filter M to raw migrated monitor data mM0 to generate LSM baseline data mB and LSM monitor data mM; andgenerate the 4D image of the subsurface based on the LSM baseline data mB and the LSM monitor data mM.

12. The computing device of claim 11, wherein the processor is further configured to: perform raw migration, with a baseline migration operator LTB on the baseline dataset dB and with a monitor migration operator LTM on the monitor dataset dM, to calculate the monitor raw baseline data mB0 and the monitor raw migrated data mM0.

13. The computing device of claim 12, wherein the processor is further configured to: calculate the remigrated baseline data mB1 by applying a baseline demigration operator LB on the reflectivity r, followed by applying the baseline migration operator LTB so that the remigrated baseline data mB1=LTBLBr; andcalculate the remigrated monitor data mM1 by applying a monitor demigration operator LM on the reflectivity r, followed by applying the monitor migration operator LTM so that the remigrated monitor data mM1=LTMLMr.

14. The computing device of claim 13, wherein the processor is further configured to: transform the reflectivity r with a filter transform , from an image domain to a filter domain;transform the remigrated baseline data mB1 with the filter transform , from the image domain to the filter domain; andtransform the remigrated monitor data mM1 with the filter transform , from the image domain to the filter domain.

15. The computing device of claim 14, wherein the processor is further configured to: compute a baseline filter fB in the filter domain based on a filter domain representation (r) of the reflectivity r, a filter domain representation (mB1) of the remigrated baseline data mB1, and a filter domain representation (mM1) of the remigrated monitor data mM1; andcompute a monitor filter fM in the filter domain based on the filter domain representation (r) of the reflectivity r, the filter domain representation (mB1) of the remigrated baseline data mB1, and the filter domain representation (mM1) of the remigrated monitor data mM1.

16. The computing device of claim 14, wherein the filter transform is a curvelet transform and the filter domain is a curvelet domain.

17. The computing device of claim 15, wherein the processor is further configured to: subject the baseline filter fB in the filter domain and the monitor filter fM in the filter domain to produce a same result when applied to the remigrated baseline data mB1 and to remigrated monitor data mM1, respectively.

18. The computing device of claim 17, wherein the processor is further configured to: transform (1) the raw migrated baseline data mB0 to transformed raw migrated baseline data (mB0) and (2) the raw migrated monitor data mM0 to transformed raw migrated monitor data (mM0) in the filter domain with the filter transform ;apply the baseline filter fB in the filter domain, to the transformed raw migrated baseline data (mB0), and applying the monitor filter fM in the filter domain, to the transformed raw migrated monitor data (mM0), to generate the LSM baseline data fB(mB) in the filter domain and the LSM monitor data fM(mM) in the filter domain; andtransform back, with an inverse transform −1, (1) the LSM baseline data fB(mg) to obtain the LSM baseline data mB=−1(fB(mB0)) and (2) the LSM monitor data fM(mM) to obtain the LSM monitor data mM=−1(fM(mM0)).

19. The computing device of claim 11, wherein the LSM is performed in a pre-stack domain.

20. A least-square migration, LSM, based method for generating an image of a subsurface, the method comprising: receiving seismic data d related to the subsurface, wherein the seismic data d includes a baseline dataset dB and a monitor dataset dM, with the monitor dataset dMbeing taken later in time then the baseline dataset dB, for the same subsurface;iteratively updating a least square migration baseline data mB and a least square migration monitor data mM based on descent directions of (1) a first baseline residual between baseline migrated seismic data LTBdB of the baseline dataset dB and baseline remigrated seismic data LTBLBmB of the least square migration baseline data mB, and (2) a second monitor residual between monitor migrated seismic data LTMdM of the monitor dataset dMand monitor remigrated seismic data LTMLMmM of the least square migration monitor data mM;migrating the seismic data d with the updated least square migration baseline data mB and the updated least square migration monitor data mM; andgenerating the 4D image of the subsurface based on the migrated seismic data d with the updated least square migration baseline data mB and the updated least square migration monitor data mM,wherein the least square migration baseline data mg and the least square migration monitor data mM are preconditioned with a baseline filter B and a monitor filter M, respectively, which are calculated based on a same common reflectivity r of the subsurface and corresponding remigrated baseline data and remigrated monitor data.