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 d.sub.B and a monitor dataset d.sub.M, calculating a baseline filter .sub.B and a monitor filter .sub.M based on a same common reflectivity r of the subsurface and corresponding remigrated baseline data m.sub.B1 and remigrated monitor data m.sub.M1 so that the base filter .sub.B applied to the remigrated baseline data m.sub.B1 equals the monitor filter .sub.M applied to the remigrated monitor data m.sub.M1, applying the baseline filter .sub.B to raw migrated baseline data m.sub.B0 and applying the monitor filter .sub.M to raw migrated monitor data m.sub.M0 to generate LSM baseline data m.sub.B and LSM monitor data m.sub.M, and generating the 4D image of the subsurface based on the LSM baseline data m.sub.B and the LSM monitor data m.sub.M.

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 d.sub.B and a monitor dataset d.sub.M, with the monitor dataset d.sub.Mbeing taken later in time than the baseline dataset d.sub.B, for the same subsurface; calculating a baseline filter .sub.B and a monitor filter .sub.M based on a same common reflectivity r of the subsurface and corresponding remigrated baseline data m.sub.B1 and remigrated monitor data m.sub.M1 so that the base filter .sub.B applied to the remigrated baseline data m.sub.B1 equals the monitor filter .sub.M applied to the remigrated monitor data m.sub.M1; applying the baseline filter .sub.B to raw migrated baseline data m.sub.B0 and applying the monitor filter .sub.M to raw migrated monitor data m.sub.M0 to generate LSM baseline data m and LSM monitor data m.sub.M; and generating the 4D image of the subsurface based on the LSM baseline data mg and the LSM monitor data m.sub.M.

2. The method of claim 1, wherein the step of calculating comprises: performing raw migration, with a baseline migration operator L.sup.T.sub.B on the baseline dataset d.sub.B and with a monitor migration operator L.sup.T.sub.M on the monitor dataset d.sub.M, to calculate the baseline raw migrated data m.sub.B0and the monitor raw migrated data m.sub.M0.

3. The method of claim 2, further comprising: calculating the remigrated baseline data m.sub.B1 by applying a baseline demigration operator L.sub.B on the reflectivity r, followed by applying the baseline migration operator L.sup.T.sub.B so that the remigrated baseline data m.sub.B1=L.sup.T.sub.BL.sub.Br; and calculating the remigrated monitor data m.sub.M1 by applying a monitor demigration operator L.sub.M on the reflectivity r, followed by applying the monitor migration operator L.sup.T.sub.M so that the remigrated monitor data m.sub.M1=L.sup.T.sub.ML.sub.Mr.

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 m.sub.B1 with the filter transform , from the image domain to the filter domain; and transforming the remigrated monitor data m.sub.M1 with the filter transform , from the image domain to the filter domain.

5. The method of claim 4, further comprising: computing a baseline filter f.sub.B in the filter domain based on a filter domain representation (r) of the reflectivity r, a filter domain representation (m.sub.B1) of the remigrated baseline data m.sub.B1, and a filter domain representation (m.sub.B1) of the remigrated monitor data m.sub.M1; and computing a monitor filter f.sub.M in the filter domain based on the filter domain representation (r) of the reflectivity r, the filter domain representation (m.sub.B1) of the remigrated baseline data m.sub.B1, and the filter domain representation (m.sub.M1) of the remigrated monitor data m.sub.M1.

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 f.sub.B in the filter domain and the monitor filter f.sub.M in the filter domain to produce a same result when applied to the remigrated baseline data m.sub.B1 and to remigrated monitor data m.sub.M1, respectively.

8. The method of claim 7, further comprising: transforming (1) the raw migrated baseline data m.sub.B0 to transformed raw migrated baseline data (m.sub.B0) and (2) the raw migrated monitor data m.sub.M0 to transformed raw migrated monitor data (m.sub.M0) in the filter domain with the filter transform , and the step of applying comprises: applying the baseline filter f.sub.B in the filter domain, to the transformed raw migrated baseline data (m.sub.B0), and applying the monitor filter f.sub.m in the filter domain, to the transformed raw migrated monitor data (m.sub.M0), to generate the LSM baseline data f.sub.B(m.sub.B0) in the filter domain and the LSM monitor data f.sub.M(m.sub.M0) in the filter domain; and transforming back, with an inverse transform .sup.?1, (1) the LSM baseline data f.sub.B(mg) to obtain the LSM baseline data m.sub.B=.sup.?1(f.sub.B(m.sub.0 )) and (2) the LSM monitor data f.sub.M(m.sub.M) to obtain the LSM monitor data m.sub.M=.sup.?1(f.sub.M(m.sub.M0)).

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 d.sub.B and a monitor dataset d.sub.M, with the monitor dataset d.sub.Mbeing taken later in time then the baseline dataset d.sub.B, for the same subsurface; and a processor connected to the interface and configured to, calculate a baseline filter .sub.B and a monitor filter .sub.M based on a same common reflectivity r of the subsurface and corresponding remigrated baseline data m.sub.B1 and remigrated monitor data m.sub.M1 so that the base filter .sub.B applied to the remigrated baseline data m.sub.B1 equals the monitor filter .sub.M applied to the remigrated monitor data m.sub.M1; apply the baseline filter .sub.B to raw migrated baseline data m.sub.B0 and applying the monitor filter .sub.M to raw migrated monitor data m.sub.M0 to generate LSM baseline data m.sub.B and LSM monitor data m.sub.M; and generate the 4D image of the subsurface based on the LSM baseline data m.sub.B and the LSM monitor data m.sub.M.

12. The computing device of claim 11, wherein the processor is further configured to: perform raw migration, with a baseline migration operator L.sup.T.sub.B on the baseline dataset d.sub.B and with a monitor migration operator L.sup.T.sub.M on the monitor dataset d.sub.M, to calculate the monitor raw baseline data m.sub.B0 and the monitor raw migrated data m.sub.M0.

13. The computing device of claim 12, wherein the processor is further configured to: calculate the remigrated baseline data m.sub.B1 by applying a baseline demigration operator L.sub.B on the reflectivity r, followed by applying the baseline migration operator L.sup.T.sub.B so that the remigrated baseline data m.sub.B1=L.sup.T.sub.BL.sub.Br; and calculate the remigrated monitor data m.sub.M1 by applying a monitor demigration operator L.sub.M on the reflectivity r, followed by applying the monitor migration operator L.sup.T.sub.M so that the remigrated monitor data m.sub.M1=L.sup.T.sub.ML.sub.Mr.

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 m.sub.B1 with the filter transform , from the image domain to the filter domain; and transform the remigrated monitor data m.sub.M1 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 f.sub.B in the filter domain based on a filter domain representation (r) of the reflectivity r, a filter domain representation (m.sub.B1) of the remigrated baseline data m.sub.B1, and a filter domain representation (m.sub.M1) of the remigrated monitor data m.sub.M1; and compute a monitor filter f.sub.M in the filter domain based on the filter domain representation (r) of the reflectivity r, the filter domain representation (m.sub.B1) of the remigrated baseline data m.sub.B1, and the filter domain representation (m.sub.M1) of the remigrated monitor data m.sub.M1.

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 f.sub.B in the filter domain and the monitor filter f.sub.M in the filter domain to produce a same result when applied to the remigrated baseline data m.sub.B1 and to remigrated monitor data m.sub.M1, respectively.

18. The computing device of claim 17, wherein the processor is further configured to: transform (1) the raw migrated baseline data m.sub.B0 to transformed raw migrated baseline data (m.sub.B0) and (2) the raw migrated monitor data m.sub.M0 to transformed raw migrated monitor data (m.sub.M0) in the filter domain with the filter transform ; apply the baseline filter f.sub.B in the filter domain, to the transformed raw migrated baseline data (m.sub.B0), and applying the monitor filter f.sub.M in the filter domain, to the transformed raw migrated monitor data (m.sub.M0), to generate the LSM baseline data f.sub.B(m.sub.B) in the filter domain and the LSM monitor data f.sub.M(m.sub.M) in the filter domain; and transform back, with an inverse transform .sup.?1, (1) the LSM baseline data f.sub.B(mg) to obtain the LSM baseline data m.sub.B=.sup.?1(f.sub.B(m.sub.B0)) and (2) the LSM monitor data f.sub.M(m.sub.M) to obtain the LSM monitor data m.sub.M=.sup.?1(f.sub.M(m.sub.M0)).

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 d.sub.B and a monitor dataset d.sub.M, with the monitor dataset d.sub.Mbeing taken later in time then the baseline dataset d.sub.B, for the same subsurface; iteratively updating a least square migration baseline data m.sub.B and a least square migration monitor data m.sub.M based on descent directions of (1) a first baseline residual between baseline migrated seismic data L.sup.T.sub.Bd.sub.B of the baseline dataset d.sub.B and baseline remigrated seismic data L.sup.T.sub.BL.sub.Bm.sub.B of the least square migration baseline data m.sub.B, and (2) a second monitor residual between monitor migrated seismic data L.sup.T.sub.Md.sub.M of the monitor dataset d.sub.Mand monitor remigrated seismic data L.sup.T.sub.ML.sub.Mm.sub.M of the least square migration monitor data m.sub.M; migrating the seismic data d with the updated least square migration baseline data m.sub.B and the updated least square migration monitor data m.sub.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.sub.B and the updated least square migration monitor data m.sub.M, wherein the least square migration baseline data mg and the least square migration monitor data m.sub.M are preconditioned with a baseline filter .sub.B and a monitor filter .sub.M, respectively, which are calculated based on a same common reflectivity r of the subsurface and corresponding remigrated baseline data and remigrated monitor data.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0014] 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:

[0015] 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;

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

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

[0018] 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;

[0019] 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;

[0020] 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

[0021] 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

[0022] 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.

[0023] 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.

[0024] 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.

[0025] 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:

[00003]$\begin{array}{cc}C\left(m\right)={.Math.d-\mathrm{Lm}.Math.}^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.

[0026] The main difficulty in solving equation (2) is in the evaluation of the operator (L.sup.TL).sup.?1, which is the inverse of the Hessian. The second part of equation (2) is the raw migrated image m.sub.0=L.sup.Td.

[0027] The author in [4] proposed to estimate the operator (L.sup.TL).sup.?1 (the inverse of the Hessian) by creating another image realization m.sub.1 by applying the de-migration operator L and the re-migration operator L.sup.T, that is:

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

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

[00005]$\begin{array}{cc}D\left(?\right)=\arg {\min }_{?}.Math.m_0-?*m_1.Math.,& \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.sub.0 and m.sub.1, and * denotes a convolution.

[0029] The filter then approximates the inverse of the Hessian operator (L.sup.TL).sup.?1 i.e.,

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

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

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

[0030] 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.

[0031] 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:

[00008]$\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}$

[0032] For computing the re-migrations m.sub.B1 and m.sub.M1, instead of de-migrating each raw migration m.sub.B0 and m.sub.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.

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

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

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

[00010]$\begin{array}{cc}{?}_B\left(m_{B1}\right)?r,{?}_M\left(m_{M1}\right)?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.

[0035] 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

[00011]$\begin{array}{cc}{.Math.\left(m_0\right)-f\left(m_1\right).Math.}^2+{?}^2{.Math.f.Math.}^2.& \left(13\right)\end{array}$

[0036] A direct solution to this equation is given by:

[00012]$\begin{array}{cc}f=\frac{\left(m_0\right){\left(m_1\right)}^{*}}{{.Math.\left(m_1\right).Math.}^2+{?}^2}.& \left(14\right)\end{array}$

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

[00013]$\begin{array}{cc}?\left(m_0\right)={?}^{-1}\left(\mathrm{??}\left(m_0\right)\right).& \left(15\right)\end{array}$

[0038] Note that equation (15) indicates that the transform effectively transforms the raw migration data m.sub.0 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 .sup.?1. 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.sup.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:

[0039] 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.

[0040] 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.

[0041] 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.

[0042] 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.

[0043] 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.

[0044] 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.

[0045] 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.

[0046] 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.

[0047] 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.

[0048] 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.

[0049] 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.

[0050] 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.

[0051] 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.

[0052] 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.

[0053] 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.

[0054] Returning to equation (12), the two filters .sub.B and .sub.M are required to satisfy the relation:

[00014]$\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.sup.T.sub.BL.sub.B is captured in the monitor filter .sub.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.

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

[00015]$\begin{array}{cc}{m}_B?{?}_B\left(m_{B0}\right),m_M?{?}_M\left(m_{m0}\right)& \left(17\right)\end{array}$

[0056] The two filters .sub.B and .sub.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.sub.B and f.sub.M (in the curvelet space in this instance):

[00016]$\begin{array}{cc}{.Math.?\left(r\right)-{?}_B?\left(m_{1B}\right).Math.}^2+{?}_B^2{.Math.{?}_B.Math.}^2+{.Math.?\left(r\right)-f_M?\left(m_{1M}\right).Math.}^2+{?}_M^2{.Math.f_M.Math.}^2-\left(.Math.?,f_B?\left(m_{1M}\right)-f_M?\left(m_{1M}\right).Math.+c.c\right).& \left(18\right)\end{array}$

[0057] 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:

[00017]$\begin{array}{cc}{f}_B=\frac{2{.Math.?\left(m_{1M}\right).Math.}^2{?\left(m_{1B}\right)}^{*}?\left(r\right)}{2{.Math.?\left(m_{1B}\right).Math.}^2{.Math.?\left(m_{1M}\right).Math.}^2+{?}_B^2{.Math.?\left(m_{1M}\right).Math.}^2+{?}_M^2{.Math.?\left(m_{1B}\right).Math.}^2},f_M=\frac{{.Math.?\left(m_{1M}\right).Math.}^2{?\left(m_{1B}\right)}^{*}?\left(r\right)}{2{.Math.?\left(m_{1B}\right).Math.}^2{.Math.?\left(m_{1M}\right).Math.}^2+{?}_B^2{.Math.?\left(m_{1M}\right).Math.}^2+{?}_M^2{.Math.?\left(m_{1B}\right).Math.}^2}.& \left(19\right)\end{array}$

The above solutions can be contrasted with those obtained through independent evaluations of the filters, f.sub.B.sup.0 and f.sub.M.sup.0, as in equation (14):

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

It is worth noticing a relation between these solutions:

[00019]$\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.sub.1, x.sub.2) is the harmonic mean of x.sub.1 and x.sub.2. This leads to a generic formula for obtained the joint filters based on any generalized mean. In particular, for the Lehmer mean L.sub.p [10]:

[00020]$\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.

[0058] Then, the solution for the 4D single iteration LSM is

[00021]$\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 .sup.?1.

[0059] 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.

[0060] 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.

[0061] 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:

[00022]$\begin{array}{cc}{?}_B^i\left(m_{B1}^i\right)?r,{?}_M^i\left(m_{M1}^i\right)?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:

[00023]$\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}$

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

[00024]$\begin{array}{cc}{m}_B^i?{?}_B^i\left(L_B^{\mathrm{iT}}{d}_B\right),m_M^i?{?}_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:

[00025]$\begin{array}{cc}{?}_B^i\left(m_{B1}^i\right)?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}$

[0063] 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.,

[00026]$\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?{?}_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).

[0064] 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.sub.B associated with a first seismic survey, the baseline survey, and a step 102 of receiving seismic data d.sub.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)).

[0065] 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.

[0066] The seismic data d.sub.B and d.sub.M are then migrated with the migration operator L and LM, respectively, in steps 104 and 106, for obtaining the raw migrated data m.sub.B0 and m.sub.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.sub.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.sup.T.sub.B is applied). This means that steps 110 and 112 are equivalent to operation L.sup.T.sub.BL.sub.Br, which results in calculating the remigration data m.sub.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.sub.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.sup.T.sub.M is applied). This means that steps 114 and 116 are equivalent to operation L.sup.T.sub.ML.sub.Mr, which results in calculating the remigration data m.sub.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.

[0067] 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.

[0068] Next, the method transforms in step 118 the reflectivity r (from step 108) and the remigration data m.sub.B1 and m.sub.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 .sub.B and .sub.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 .sub.B and .sub.M are applied to calculate LSM migration data m.sub.B?F.sub.B(m.sub.B0) and m.sub.M?.sub.M(m.sub.M0). In step 124, the baseline LSM migration data m.sub.B is transformed back to the image domain and in step 126, the monitor LSM migration data m.sub.M is also transformed back to the image domain. Based on the LSM migration data m.sub.B and m.sub.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.

[0069] 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.

[0070] 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.

[0071] 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.sub.B and a monitor dataset d.sub.M, with the monitor dataset d.sub.Mbeing taken later in time then the baseline dataset d.sub.B, for the same subsurface, a step 502 of calculating a baseline filter .sub.B and a monitor filter .sub.M based on a same common reflectivity r of the subsurface and corresponding remigrated baseline data m.sub.B1 and remigrated monitor data m.sub.M1, a step 504 of applying the baseline filter .sub.B to raw migrated baseline data m.sub.B0 and applying the monitor filter .sub.M to raw migrated monitor data m.sub.M0 to generate LSM baseline data m.sub.B and LSM monitor data m.sub.M, and a step 506 of generating the image of the subsurface based on the LSM baseline data m.sub.B and the LSM monitor data m.sub.M.

[0072] The step of calculating may include performing raw migration, with a baseline migration operator L.sup.T.sub.B on the baseline dataset d.sub.B and with a monitor migration operator L.sup.T.sub.M on the monitor dataset d.sub.M, to calculate the monitor raw migrated data m.sub.B0 and the baseline raw migrated data m.sub.M0. This step may further include calculating the remigrated baseline data m.sub.B1 by applying a baseline demigration operator L.sub.B on the reflectivity r, followed by applying the baseline migration operator L.sup.T.sub.B so that the remigrated baseline data m.sub.B1=L.sup.T.sub.BL.sub.Br, and calculating the remigrated monitor data m.sub.M1 by applying a monitor demigration operator L.sub.M on the reflectivity r, followed by applying the monitor migration operator L.sup.T.sub.M so that the remigrated monitor data m.sub.M1=L.sup.T.sub.ML.sub.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.sub.B1 with the filter transform , from the image domain to the filter domain, and transforming the remigrated monitor data m.sub.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.sub.B in the filter domain based on a filter domain representation (r) of the reflectivity r, a filter domain representation (m.sub.B1) of the remigrated baseline data m.sub.B1, and a filter domain representation (m.sub.M1) of the remigrated monitor data m.sub.M1, and computing a monitor filter f.sub.M in the filter domain based on the filter domain representation (r) of the reflectivity r, the filter domain representation (m.sub.B1) of the remigrated baseline data m.sub.B1, and the filter domain representation (m.sub.M1) of the remigrated monitor data m.sub.M1.

[0073] 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.sub.B in the filter domain and the monitor filter f.sub.M in the filter domain to produce a same result when applied to the remigrated baseline data m.sub.B1 and to remigrated monitor data m.sub.M1, respectively, and/or transforming (1) the raw migrated baseline data m.sub.B0 to transformed raw migrated baseline data (m.sub.B0) and (2) the raw migrated monitor data m.sub.M0 to transformed raw migrated monitor data (m.sub.M0) in the filter domain with the filter transform . The step of applying may include applying the baseline filter f.sub.B in the filter domain, to the transformed raw migrated baseline data (m.sub.B0), and applying the monitor filter f.sub.M in the filter domain, to the transformed raw migrated monitor data (m.sub.M0), to generate the LSM baseline data f.sub.B(m.sub.B ) in the filter domain and the LSM monitor data f.sub.M(m.sub.M) in the filter domain, and transforming back, with an inverse transform .sup.?1, (1) the LSM baseline data f.sub.B(m.sub.B ) to obtain the LSM baseline data m.sub.B=.sup.?1(f.sub.B(m.sub.B0)) and (2) the LSM monitor data f.sub.M(m.sub.M) to obtain the LSM monitor data m.sub.M=.sup.?1(f.sub.M(m.sub.M0)). The inverse transform .sup.?1 takes a quantity from the curvelet domain to the image domain.

[0074] 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 .sub.B and .sub.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

[00027]$\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 .sub.B and .sub.M before computing the updates for the reflectivities.

[0075] 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).

[0076] 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.sub.M are updated in the above steps.

[0077] 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.

[0078] 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.

[0079] 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.

[0080] 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.

[0081] 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.

[0082] 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.

[0083] 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.

[0084] 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.

[0085] 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.

[0086] 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.

[0087] 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

[0088] 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.sup.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.

[0089] [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.