Method and apparatus for anisotropic PP / PS tomography without well information

Abstract

A seismic exploration method 1100 includes obtaining 1110 PP and PS seismic datasets representing P-waves and S-waves emerging from the subsurface formation. PP kinematic invariants and PS kinematic i

Description

[0001] Embodiments of the subject matter disclosed herein generally relate to a system and method for exploring a subsurface formation by building an anisotropic PP / PS velocity model using seismic data without using well information. DISCUSSION OF THE BACKGROUND

[0002] Seismic tomography is a technique used to create images of a subsurface formation by analyzing seismic data (such as, seismic waves’ travel times, amplitudes, etc.) acquired by recording waves emerging from therefrom in response to injecting seismic waves into the subsurface formation. Essentially, the seismic tomography yields a 3D model of at least one parameter, such as velocity, by tracking how traveling through different materials changes the seismic waves.

[0003] The seismic waves traveling in the probed underground structure are typically P-waves (where P stands for primary or pressure waves). The P-waves are compressional waves that are longitudinal (meaning that the direction of particle movement is mainly along the direction of wave propagation), the wave causing a trail of compressions and dilatations in the medium they move through. At an interface between layers of different materials in which the seismic waves move with different velocities, the P-waves may be reflected back in the medium they come from (PP reflection), transmitted or refracted in the medium beyond the interface, and may also be converted to S-waves (as PS reflection or transmission). The S-waves (where S stands for secondary or shear) are transversal waves, meaning that the direction of particle movement of an S wave is mainly perpendicular to the direction of wave propagation. The S-waves travel slower than the P-waves and propagate only in solids.

[0004] When performing seismic imaging from reflected waves, the aim is to reposition and quantify the model heterogeneities generating the reflected waves. In seismic data processing, "migration" refers to a technique that geometrically repositions seismic events (i.e., P-wave or S-wave reflections) identified in seismic data set to their actual location, creating a more accurate image of underground structures by correcting distortions caused by dipping reflectors and wave propagation paths. Migration essentially “moves” the seismic events to where they actually originated rather than where they were recorded at the surface. Stacking is a method to enhance the signal carrying information about the subsurface formation by adding traces of several shots, offset (i.e., the distance between a seismic source and a receiver) or aperture angle at reflectors to enhance signal to noise ratio. A background velocity model is necessary to reposition the events though the migration process. The process of analyzing the travel time of the reflected events and deriving a background velocity model is what we call here tomography. When only PP reflected events are analyzed, only a P-wave velocity model is necessary. When PS reflected events are considered, a P-wave and a S-wave velocity model are necessary. Such P-wave and S-wave velocity models may also account for anisotropy, that is, the change of the velocity with the direction of propagation.

[0005] A tomography is a method that allows to recover the parameters characterizing the background velocity model from the kinematic parameters observed and picked on seismic data. It is an optimization process. The kinematic parameters can be directly picked in the data domain (e.g., travel times, slopes, curvatures), or, more frequently, in the migrated image domain (e.g., dips, curvatures, residual move out (RMO), etc.) because it helps with the interpretation and denoising. The kinematic parameters picked in the migrated image domain depend on the velocity model used for the migration. Therefore, the kinematic parameters picked in the migrated image domain are not kinematic invariants, and the corresponding tomographic optimization process is frequently done based on a local linearized approximation while the relationship between the model parameters and kinematic data is non-linear.

[0006] A non-linear optimization can be performed by converting the kinematic parameters characterizing the picked migrated data into kinematic invariants. This provides an efficient and flexible approach for velocity model building (i.e. building the background velocity model). An efficient class of non-linear tomographic approaches is the family of non-linear slope tomography methods. These methods are based on a concept of kinematic invariants defined as locally coherent events characterized by their position and slopes in the un-migrated pre-stack domain. Kinematic invariants offer a versatile tool for velocity model building. They can be derived by a process called kinematic demigration from dip and residual move-out picks made either in pre-stack depth migrated (preSDM) or pre-stack time migrated (preSTM) domains, or even directly from picks in the seismic data domain (the un-migrated domain). Note that unlike time migration that assumes that background velocity varies vertically but only smoothly laterally, depth migration allows for velocity to vary both vertically and horizontally.

[0007] The kinematic invariants being in the un-migrated domain, the picking only needs to be done once. The classical iterative velocity update made of several iterations of residual move-out (RMO) picking, pre-stack migration and velocity update can be replaced by a more efficient sequential approach involving a single preSDM and a single picking step (typically dip and RMO picking followed by kinematic demigration) followed by a non-linear tomographic inversion. In seismic data processing, RMO refers to the remaining depth (or time when considering time migration) shift to apply at each reflected event in the migrated traces gathered at the same location for various offsets, shot position, aperture angle, etc. to ensure its flatness along offsets, shot or receiver position, aperture angle, etc. Essentially, the RMO optimizes image stack quality in seismic data processing by applying an adjustment needed to properly align seismic reflections from different offsets before stacking.

[0008] Estimating an accurate anisotropic background velocity model is a challenge in seismic imaging (i.e., generating an image of the subsurface formation via an accurate model thereof). The simplest type of anisotropy is the case of transverse isotropy. In this configuration the P-wave and SV wave velocities (SV stand for vertically polarized S-wave, i.e. those having a vertical component in their polarization vector) depends in addition to the vertical P and S-wave velocity on the two anisotropic parameters e and 6 [7], The RMO at short and large offsets of the PP events does not allow to decouple vertical P-wave velocity, e and 6, but provide two coupled relationships between these three parameters. Adding RMO of PS event does not solve the issue but provides an additional coupled relationship involving S-wave velocity. Estimating an anisotropic background velocity model can’t then be done accurately when only PP or even PP and PS RMOs are available. But it can be done using the additional information of the vertical positioning of specific horizons identified on well data. Thus, there is a need for new methods and devices able to generate an accurate anisotropic velocity model of an explored subsurface formation thereby allowing locating sought-after resources (such as oil and gas) when no well data is available. SUMMARY OF THE INVENTION According to an embodiment, there is a seismic exploration method for locating targeted resources in a subsurface formation. This method includes (1) obtaining PP and PS seismic datasets representing P-waves and S-waves emerging from the subsurface formation in response to injecting P-waves into the subsurface formation, (2) extracting PP kinematic invariants and PS kinematic invariants from PP and PS seismic dataset, respectively, and (3) obtaining a complete anisotropic model of the subsurface formation depending on P-wave velocity, 1^, S-wave velocity, 1^, and Vertically or Tilted Transverse Isotropy (VTI or TTI) Thomsen anisotropic parameters e and 6, so that seismic PP and PS migrated images predicted using the complete anisotropic model fulfills predetermined criteria. The predetermined criteria include the RMO flattening of PP and PS kinematic invariants and the registration of PP and PS migrated events (registration is the process that ensures that the same reflected event observed in the PP and PS is migrated at the same position in the image). The complete anisotropic model can be obtained through the optimization of a cost function involving terms describing the RMO flattening of PP and PS kinematic invariant and the quality of the registration of PP and PS migrated events. For example, the joint inversion can be done through a local optimization (gradient based) through iterative loops. One loop may consist of performing a first PP tomography using PP kinematic invariants, a PS tomography using PS kinematic invariants, and then a PP and PS events registration to update Vp, and the Thomsen anisotropy parameters e, 6. The anisotropic model is not constrained using well information acquired from one or more wells dig into the subsurface formation. The output criterion may be convergence of the optimization of the cost function.

[0010] In some embodiments, the anisotropic model is updated by iteratively performing a sequence of steps of a loop. For example, in each loop, the variables Vp and e of an initial anisotropic model are updated by performing the first PP tomography using the PP kinematic invariants to obtain a modified anisotropic model; variable of the modified anisotropic model by performing the PS tomography using the PS kinematic invariants to obtain a complete anisotropic model; variable Vp of the complete anisotropic model is then updated using the Vp / Vs derived from the PP and PS image registration to obtain an updated anisotropic model; and, last, variable 6 and e of the updated anisotropic model are updated by performing a second PP tomography while keeping Vp,Vs fixed to obtain a next-iteration anisotropic model. The next-iteration anisotropic model is usable as the initial model for repeating the first PP tomography, the PS tomography and, if the complete anisotropic model is not output, the image registration and the second PP tomography (i.e., another loop). When a global output criterion is fulfilled, the complete anisotropic model is output as the model of the subsurface formation.

[0011] A seismic processing apparatus for locating targeted resources in a subsurface formation has a communication interface for obtaining PP and PS seismic datasets representing P-waves and S-waves emerging from the subsurface formation in response to injecting P-waves into the subsurface formation, and a processor connected to the communication interface configured: (1) to extract PP kinematic invariants and PS kinematic invariants from PP and PS seismic dataset, respectively, and (2) to obtain an anisotropic model of the subsurface formation depending on P-wave velocity, 1^, S-wave velocity, 1^, and Thomsen anisotropy parameters e and 6, so that seismic data predicted using the anisotropic model fulfills an output criterion based on a comparison of the predicted seismic data with the PP and PS seismic datasets. The processor is configured to perform a PP and PS tomography constrained by a registration of PP and PS migrated events for obtaining the anisotropic model. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0013] FIG. 1 is a schematic diagram of seismic data acquisition using ocean bottom nodes;

[0014] FIG. 2 in a graphical representation of the process of the registration between PP and PS migrated images.

[0015] FIG. 3 schematically illustrates a non-linear depth tomography technique with the kinematic invariants obtained by kinematic de-migration of locally coherent events;

[0016] FIGs. 4A and 4B illustrate observed elementary kinematic information used for picking RMO and dip;

[0017] FIG. 5 graphically illustrates several concepts used in non-linear slope tomography technique;

[0018] FIG. 6 illustrates the positioning in the migrated image of a kinematic invariant when kinematic migration is performed for two different velocity models;

[0019] FIG. 7 schematically illustrates an iterative anisotropic velocity model optimization scheme according to an embodiment;

[0020] FIG. 8 schematically illustrates a specific embodiment of the workflow of an anisotropic PP / PS tomography method without using well information according to an embodiment;

[0021] FIGs. 9A and 9B illustrate a TTI velocity model of a subsurface formation used to test the anisotropic PP / PS tomography method without using well information;

[0022] FIGS. 10A, 10B, 10C, and 10D graphically illustrate the results of applying the anisotropic PP / PS tomography method without using well information;

[0023] FIG. 11 is a flowchart of a method according to an embodiment; and

[0024] FIG. 12 is a schematic diagram of an apparatus configured to perform anisotropic PP / PS tomography methods that do not use well information according to various embodiments. DETAILED DESCRIPTION OF THE INVENTION

[0025] 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 seismic data acquired using ocean bottom nodes (OBN). However, the techniques to be discussed next are not limited to such data acquisition setup being pertinent to seismic data acquired on land or using ocean bottom cables, or as soon as PP and PS reflected events can be observed.

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

[0027] The methods described herein could 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. The methods could also be employed for surveying and monitoring for windfarm applications, both onshore and offshore, and also for medical imaging applications.

[0028] Before focusing on various embodiments, FIG. 1 illustrates seismic data acquisition using ocean bottom nodes. Water-bottom and land seismic acquisition systems use multi-component detectors 100 (only one labeled) for recording seismic data at locations on the top surface 118 (i.e., water bottom) of the explored subsurface formation 126. These multi-component detectors are typically able to detect both P and S waves, for example, by including a hydrophone and a three-dimensional (3D) geophone or accelerometer. One or more sources 120 (e.g., towed by a vessel 121 above the detectors 100 in FIG. 1 along sail lines 110) generate seismic P-waves waves 122 that penetrate the top surface 118 to then propagate into the subsurface formation 126. Inside subsurface formation 126, these injected P-waves waves are reflected, refracted and / or generate S-waves at various interfaces or reflectors R. Some of these waves (such as 123) emerge from the subsurface formation and are recorded by one of the multi-component detectors. The seismic waves propagating through water are P waves (S-waves do not propagate through water) and the waves emerging from the subsurface formation are both P and S waves (yielding PP datasets and PS datasets). The detectors record traces (i.e., amplitude versus time series) based on which reflector’s location and other characteristics (e.g., incidence and azimuth angles) may be calculated. For example, a traditional OBN dataset is acquired with an inline distance ID (along the sail lines) between the nodes of about 100 m, and a cross-line distance CD between the nodes of about 400 m. The axes X and Y in FIG. 1 indicate the inline and the crossline directions. The recorded seismic data is subsequently processed in multiple steps (e.g., pre-processing to separate the PP dataset and the PS dataset from the recorded seismic data, velocity model building, and imaging) to enable identifying and locating sought-after resources such as oil and gas, water, minerals, etc. in the subsurface formation.

[0029] The seismic data processing techniques according to embodiments described in this section obtain an anisotropic velocity model of the subsurface by processing multi-component seismic data acquired, for example, as in FIG. 1 or equivalent data acquisition geometries. In contrast to isotropic velocity models, the anisotropic velocity models (such as, a tilted transverse isotropy (TTI) or even vertical transverse isotropy (VTI) velocity model) assume that wave propagation velocities are not direction independent, and thereby knowledge about anisotropy parameters needs to complement the velocity values. A VTI velocity model is a 4-parameter model e, 6) yielding information inside the subsurface formation probed by seismic waves. Here, Vp is the vertical P-wave velocity, Vs is the vertical S-wave velocity, e and 6 are Thomsen anisotropy parameters defined in [7], The Thomsen anisotropy parameter e represents the difference between horizontal and vertical P-wave velocities, and the Thomsen anisotropy parameter 6 quantifies the difference between vertical P-wave velocity and the apparent P-wave velocity deduced from the NMO curvature of the reflected event, vpNM0. The TTI assumes that rock layers in the subsurface formation have a preferred wave propagation direction (i.e., a direction along which the wave propagation direction is largest) that is not necessarily vertical. This assumption allows for tilted bedding planes within the subsurface formation and, thereby, provides a more accurate representation of complex structures.

[0030] The significance and relationship between the TTI velocity model and observables obtained from seismic data, for a single reflector at the bottom of a homogeneous VTI anisotropic layer, are illustrated next. Note that vp and Vy are vertical velocities in case of VTI or velocities perpendicular to geological bed in case of TTI. For a homogeneous velocity and a flat horizontal reflector, the PP and PS NMO velocities (vpNM0 and vpsNM0) and the horizontal P-wave velocity (yph) have the following relationships with the anisotropic parameters: Vpnmo = V1 + 2<5 (1) ^ph — VpNMOy / 1 T 2?) (2) VrsNm = + (3) with t] = y0 = and a = t] Yo (1 + 25). Here, vpNM0 (extracted from the PP data) is a move-out corrected vp, vsNM0 (extracted from the PS data) is a move-out corrected 14, and vph (extracted from long offset PP data) is the horizontal P-wave velocity. Estimating four model parameters e, 3) from the three data parameters VpNM0, Vph, VpsNM0 is not possible, as well as estimating (Vp,e, 3) from the kinematic parameter describing the PP reflected events, VpNM0, VPh- One way to solve the challenge is to complement data with well information. Well information is based on measurements (e.g., composition, porosity, etc.) of rock samples obtained at specific depths inside a well or on measurements (e.g., natural gamma ray, electrical, acoustic, stimulated radioactive responses, electromagnetic, nuclear magnetic resonance, pressure and other properties of the rocks and their contained fluids) made by lowering an instrument into a well dig into the explored subsurface formation. However, digging wells is expensive and time consuming and may still miss significant features of the subsurface formation. Therefore, the following methods do not use well information but use information obtained by registration between the PP and PS migrated events, that is, by collocating PP and PS reflectors (i.e., △Zpp-Ps=0 as illustrated in FIG. 2). This registration information can be converted into a volumetric volume of y0 as shown in [8],

[0031] Four decoupled quantities may be extracted from the data: VpNM0, Vph, VpsNM0 and PP-PS registration information (as for example of y0). The VTI velocity model has also four model parameters (yp,Vs, e, d). The relationship between the four decoupled quantities extracted from the data and the four model parameters is not linear but the problem of obtaining the model parameters from the quantities extracted from the data is sufficiently well constrained to be solved by a joint nonlinear optimization scheme.

[0032] FIG. 3 illustrates a non-linear slope tomography technique with the kinematic invariants obtained by de-migration of locally coherent events as described in [3] (a complete list of references invoked in this description is provided at the end of this section). In this method, dRMO and dip are picked 320 from the acquired and then migrated seismic data 310 (i.e., preSTM or preSDM); that is, each locally coherent event is defined by its structural dip and dRMO.

[0033] The kinematic invariants then calculated at 340 after kinematic de-migration 330 include the shot and receiver positions (S and R), the two-way time (T) and the two-way time slope (slope) computed in a constant-offset un-migrated seismic cube as described later. The kinematic invariants are velocity model independent and are similar to the observed data entering a travel-time slope tomography (as described in [4]). As a result, both the time consuming preSTM or preSDM 310 and picking stages 320 have been removed from the iterative velocity updating loop 300 (i.e. steps 350-380) and a non-linear depth tomography has emerged yielding a velocity model at 390. In this method, a cost function used as a criterion to exit the loop at 370 involves an assessment of the local flattening of the local derivative of RMO (i.e., dRMO), migrated facets (see dRMO in Fig. 5), which means the RMO has been minimized.

[0034] The non-linear slope tomography method needs an accurate picking in preSTM or preSDM domains. Such picking of locally coherent events in the preSDM or preSTM domain may involve an automated dense picking of both the structural dip and the dRMO. For example, FIGs 4A and 4B illustrate observed elementary kinematic information used for dRMO and dip picking. Specifically, FIG. 4A illustrates picking RMO (z(h)) from a common image gather (i.e., a gather of migrated traces obtained at the same position for different offset, shot, receiver or aperture angles), and FIG. 4B illustrates picking a dip corresponding to the picked RMO from the image sub-stack.

[0035] Kinematic de-migration and migration used in non-linear slope tomography are illustrated in FIG. 5. A preamble of kinematic de-migration is (time or depth) migration of seismic data (at 310) using an initial velocity model. The migrated data is illustrated on the right side of FIG. 5. In a slope tomography approach, a common offset h subset of a (time or depth) migrated facet is characterized by its central position (X, Z) and a structural dip together with the RMO slope (dRMO) for the offset h at the central position. Note that alternatively to common offset data can be a common shot or receiver position, or aperture angle gather as described in

[13] and

[14] , The dRMO slope may be derived from a picked high-order RMO curve.

[0036] After kinematic de-migration (by specular offset ray tracing upwards from the image point in the depth case), the data in the unmigrated domain is illustrated on the left side of FIG. 5. In this unmigrated domain, a locally coherent event is associated with a set of kinematic invariants including: (a) a pair of source and receiver locations S and R, respectively, (both locations being at the water surface or the water bottom or any other location), (b) a two-way time (T) and a slope (slopem) in the common offset (h) gather and (c) a slope (sloped in a common midpoint (CMP) gather.

[0037] As explained in [5], using a first order Taylor expansion of the preSDM imaging equations in the neighborhood of the picked event, the slopeh invariant can be expressed as a linear combination of the measured dRMO and of ray parameters obtained from the specular de-migration, for example: dRMO = (4)

[0038] The computed kinematic invariants do not depend on the velocity model as long as the same velocity model is used for the migration and kinematic de-migration. Note that the kinematic invariant can be also directly picked in the unmigrated domain (and then no kinematic demigration step is necessary).

[0039] Returning now to FIG. 3, a complete velocity model is obtained by performing at least once steps 350-360-370-380. At 350, the local kinematic invariants may be kinematically re-migrated 350 with a current velocity model (if other than the initial velocity model) for then predicting 360 the dRMO corresponding to the current velocity model. The dRMO prediction 360 involves a non-linear remodeling of the wave path. The current velocity model is then compared 370 with the initial velocity model, without having to re-run a preSTM or preSDM 310 and dip and RMO picking for the new updated velocity model. When the comparison indicates 370 that the current velocity model has minimized dRMO (e.g., using a cost function), the current velocity model is output as the result of this tomographic technique. Otherwise, the invariants re-migrated at 350 may be used to perform 380 a linearized velocity update using the Frechet derivatives computed along the modeled ray-paths. The velocity model updated at 380 becomes the current velocity model for repeating step 350 and following steps as needed. If the initial dip and RMO picking are dense and accurate enough, a complete velocity model (i.e., the velocity model output at 370), as accurate as it can be in the illuminated area, is obtained directly from a single picking step.

[0040] FIG. 6 illustrates the velocity model update when using invariants picked from the data. The source 610 and the receiver 620 are fixed from one iteration to the next, but initial facet location 630 is moved to a next iteration location 640. Note that improving quality of the velocity model means that dRMO (i.e., derivative of depth z relative to offset h) approaches zero (vanishes).

[0041] In the conventional anisotropic tomography, well information is used to guide the model updates. This information can be typically the depth of some specific reflectors that are identified in the data or migrated image and matched to the depth observed in the well. The following embodiments perform anisotropic tomography without using well information. In the conventional methods, Vp e, 6 are updated matching the short and long offset dRMO of PP reflected event and matching at the wells the depth positioning of a set of selected PP events. is update afterwards considering dRMO of PS reflected events or co-depthing reflectors in PP and PS images (i.e., matching location of PP migrated events and location of PS migrated events. In the methods according to various embodiments described below, without using any well information, the co-depthing of PP migrated events and location of PS migrated events is directly used to estimate (Vp,Vs, e, 3) in a joint PP-PS non-linear slope tomographic process.

[0042] FIG. 7 illustrates an iterative velocity model optimization scheme according to one embodiment. The values variables (Vp,Vs, e, 3) of the (initial or later) velocity model 710 are non-linearly updated 720 by minimizing dRMOpp, dRMOps, and AZpp-ps based on input data including PP kinematic invariants 700, PS kinematic invariants 701, obtained by dip and dRMO picking and kinematic demigration and registration of PP and PS events 702. In other words, the iterative optimization scheme starts with initial velocity model 710 (yp,Vs, e, 3), which is iteratively updated based on the PP kinematic invariants 700, PS kinematic invariants 701, and a set of parameters AZpp.ps describing the registration of a set of PP and PS migrated events 702. The minimization of dRMOpp, dRMOps and AZpp.ps is assessed using a predefined joint cost function. When no significant variation of the cost function is observed, the complete velocity model 730 is output. The registration parameters may be introduced in the process by converting the set of AZpp.psobserved between short offset PP and PS migrate image into a volumetric y0 = — (as shown in [8]) that can be used directly as a volumetric constrain during Vs the iterative e, 6) model update from PP and PS kinematic invariants.

[0043] FIG. 8 schematically illustrates a specific embodiment of the workflow of an anisotropic PP / PS tomography method without using well information according to an embodiment. Starting with an initial velocity model 810 (Vpo> e> Sq), a first PP tomography (“PP tomol”) updates Vp and e with Vs, y, and 6 fixed (which means updating VpNM0 = + 26 and t] = |^). This first PP tomography uses the dip and dRMOpp computed from the PP kinematic invariants to update Vp (to Vpl) and e (to d) while keeping unchanged Vs, y, and 6 (i.e., fixed) at the values corresponding to the initial velocity model 810. The first PP tomography outputs a modified velocity model 820 (ypl,Vs0, d, 6o). Further starting with the modified velocity model 820 (Vpl,Vs0, d, <50), a PS tomography (“PS tomo”) updates (to 7S1) while keeping Vp, e and 6 unchanged at the values corresponding to the modified velocity model 820. The PS tomography uses the dip and dRMOps computed from the PS kinematic invariants. The PS tomography obtains a complete velocity model 830 (Vpl,Vsl, d, <50). If the complete model’s quality meets an output criterion (e.g., a cost function converges, being adequately improved after few iterations of the loop in Fig. 8), the complete velocity model 830. The output criterion (e.g., a cost function optimization) ensures convergence of the iterative process (e.g., there are no more significant variations in the various data components of the cost function (dRMOpp, dRMOps, and AZpp.ps).

[0044] When the output criterion is not met, a further improvement of the velocity model is pursued by updating the complete model 830 based on data yielded by PP and PS event registration. For example, a volumetric y0 model may be extracted from PP to PS event registration by co-depthing PP reflectors to PS reflectors. That is, the model parameter Vp is updated to Vp2 based on y0 while maintaining fixed all the other model parameters, thus obtaining an updated velocity model 840 e', <50).

[0045] A second PP tomography (“PP tomo2”) technique updating 8 (to and e (to e") (according to a fixed value of 77) is then applied to the updated velocity model 840 to obtain a next velocity model 850. The next velocity model 850 (Vp2, Vsi, e", 6-^) becomes the initial model for a next-iteration of the techniques illustrated in FIG. 8.

[0046] The method illustrated in Fig. 8 has been tested using a horizontal homogeneous layer with PP and PS invariants obtained by kinematic demigration for the exact velocity model. The exact y0 = 2 is also used in the iterative inversion process. This first test yielded the results illustrated in Table 1. Attribute Initial Recovered Exact 3500 2996 3000 Fs (ms-1) 1750 1498 1500 6 (%) 30.0 15.3 15.0 5(%) 15.0 7 5 Table 1

[0047] For a second test, the same experiment is done except that now in the exact model 6 = -3.0%. This second test yielded the results illustrated in Table 2. Attribute Initial Recovered Exact 3500 3001 3000 (ms-1) 1750 1500 1500 30.0 14.9 15,0 5 (%) -7.0 —3.0 —3.0 Table 2

[0048] Further, a more complex test aimed to reconstruct the VTI velocity model for a 6 km deep subsurface formation characterized by a one dimensional velocity variation with depth (within the range of 1.75 to 6.41 km / s) and horizontal reflectors at every 40 m. FIG. 9A graphically illustrates the true velocity vp (continuous line) and 14 (dashed line) as functions of depth z. FIG. 9B graphically illustrates the anisotropy parameters (e and 6 using continuous line and dashed line, respectively) as functions of depth z.

[0049] FIGS. 10A, 10B, 10C and 10D are graphical representation of this complex test’s results. In these graphs of Vp,Vs, e, and 6 as functions of depth, the start model is illustrated using large-dashed lines, the actual model is illustrated using continuous lines and the method results are illustrated using small-dashed lines. The methods result accurately retrieved the exact model.

[0050] FIG. 11 is a flowchart of a seismic exploration method 1100 according to an embodiment. The method includes obtaining 1110 PP and PS seismic datasets representing P-waves and S-waves emerging from the subsurface formation in response to injecting P-waves into the subsurface formation. As previously mentioned, these PP and PS seismic datasets result from pre-processing seismic data acquired by multi-component receivers placed on the top surface of the explored subsurface formation (e.g., water bottom or land surface). The method 1100 obtains the anisotropic model without using well information acquired from one or more wells dig into the subsurface formation.

[0051] The method 1100 further includes extracting 1120 PP kinematic invariants and PS kinematic invariants from the PP dataset and the PS seismic dataset, respectively. The kinematic invariant may be automatically extracted as explained in the descriptions of FIGS. 4 and 5.

[0052] The method 1100 then includes obtaining 1130 a complete anisotropic model () of the subsurface formation depending on P-wave velocity, 1^, S-wave velocity, 1^, and Thomsen anisotropy parameters e and 6, so that seismic data predicted using the complete anisotropic model fulfills an output criterion based on an optimum value of the cost function representing a joint minimization of dRMOpp, dRMOps and AZpp.ps. When no significant variation of the cost function is observed from one iteration to the other the iterative updating loop is stopped. The obtaining of the complete anisotropic model of the subsurface formation includes performing a first PP tomography, a PS tomography, and a PP and PS data registration before updating the Thomsen anisotropy parameter 6 by performing a second PP tomography as described relative to FIG. 8, or for generally to perform jointly or sequentially PP tomography, PS tomography and PP-PS event registration.

[0053] In one embodiment, obtaining of the complete anisotropic model of the subsurface formation includes: (1) updating Vp and e of an initial anisotropic model by performing the first PP tomography using the PP kinematic invariants to obtain a modified anisotropic model; (2) updating of the modified anisotropic model by performing the PS tomography using the PS kinematic invariants to obtain a complete anisotropic model; (3) updating Vp of the complete anisotropic model by performing the PP and PS data registration matching reflector locations as reconstructed using P-waves and S-waves to obtain an updated anisotropic model; and (4) updating 6 (to 8') and e (to e") (according to a fixed value of tj) of the updated anisotropic model by performing a second PP tomography to obtain a next-iteration anisotropic model usable as the initial model for repeating performing the first PP tomography and the PS tomography. Note that the complete anisotropic model is output when a convergence criterion is met.

[0054] Thus, the method 1100 may also include outputting the complete anisotropic model when an output criterion is fulfilled. Otherwise, the method continues with performing the PP and PS data registration, the second PP tomography, the first PP tomography using the next-iteration anisotropic model as the initial anisotropic model, and the PS tomography to obtain a new version of the complete anisotropic model.

[0055] In view of the tomographic technique described relative to FIG. 3, the first PP tomography may be performed starting from a current anisotropic model being the initial anisotropic model, until P-wave RMO predicted based on the current anisotropic model is minimized: (i) performing a kinematic remigration of the PP kinematic invariants using the current anisotropic model, (ii) calculating the P-wave RMO, and (iii) if the P-wave RMO has been minimized, outputting the current anisotropic model as the modified anisotropic model, otherwise updating the current anisotropic model by varying the vp and e of the current anisotropic model.

[0056] The PS tomography may be performed starting from a current anisotropic model being the modified anisotropic model, until S-wave residual moveout, RMO, predicted based on the current anisotropic model is minimized, by: (i) performing a kinematic remigration of the PS kinematic invariants using the current anisotropic model, (ii) calculating the S-wave RMO, and (iii) if the S-wave RMO has been minimized, outputting the current anisotropic model as the complete anisotropic model, otherwise further updating the current anisotropic model by varying the 14 of the current anisotropic model.

[0057] The second PP tomography may be performed starting from a current anisotropic model being the updated anisotropic model, until P-wave residual moveout, RMO, predicted based on the current anisotropic model is minimized, by: (i) performing a kinematic remigration of the PP kinematic invariants using the current anisotropic model, (ii) calculating the P-wave RMO, and (iii) if the P-wave RMO has been minimized, outputting the current anisotropic model as the next-iteration anisotropic model, otherwise update the current anisotropic model by varying the vp of the current anisotropic model.

[0058] The method 1100 may be performed by a seismic exploration apparatus 1200 (i.e., an apparatus configured to perform anisotropic PP / PS tomography methods that do not use well information according to various embodiments) as illustrated in FIG. 12. Hardware, firmware, software ora combination thereof may be used to perform the various steps and operations described herein. The computing device 1200 is suitable for performing the activities described in the above embodiments and may include a server 1201. Such a server 1201 may include a central processor (CPU) 1202 coupled to a random access memory (RAM) 1004 and to a read-only memory (ROM) 1206. ROM 1206 may also be other types of storage media to store programs, such as programmable ROM (PROM), erasable PROM (EPROM), etc. Processor 1202 may communicate with other internal and external components through input / output (I / O) circuitry 1208 and bussing 1210 to provide control signals and the like. Processor 1202 carries out a variety of functions as are known in the art, as dictated by software and / or firmware instructions.

[0059] Server 1201 may also include one or more data storage devices, including hard drives 1212, CD-ROM drives 1214 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 1216, a USB storage device 1018 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 1214, disk drive 1212, etc. Server 1201 may be coupled to a display 1220, 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 1222 is provided, including one or more user interface mechanisms such as a mouse, keyboard, microphone, touchpad, touch screen, voice-recognition system, etc.

[0060] Server 1201 may be coupled to other devices, e.g., 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 1228, which allows ultimate connection to various landline and / or mobile computing devices.

[0061] As described above, the apparatus 1200 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.

[0062] The processor 1202 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.

[0063] In an example embodiment, the processor 1202 may be configured to execute instructions stored in the memory device 1204 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.

[0064] The term “about” is used in this application to mean a variation of up to 20% of the parameter characterized by this term.

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

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

[0067] The disclosed embodiments provide a system for hyperspectral imaging using a supercontinuum light source. 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.

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

[0069] 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 and reflects the current level of knowledge of a person skilled in seismic exploration. [1] Alkhalifah, T. and Tsvankin, I., 1995, Velocity analysis for transversely isotropic media, Geophysics Vol. 60 No. 5, pp. 1550-1566. [2] Audebert, F., Granger, P-Y., Herrenschmidt, A., 1999. CCP-Scan technique: true common conversion point sorting and converted wave velocity analysis solved by PP and PS Pre-Stack Depth Migration. SEG 1999 Expanded Abstracts. [3] Audebert, F., Granger, P-Y., Gerea, C., Herrenschmidt, A., 2001. Can joint PP and PS velocity analysis manage to corner 5, the anisotropic depthing parameter? EAGE 64th Conference and Exhibition, Florence, Italy, 2002. [4] Billette, F., Lambare, G. Velocity macro-model estimation by stereotomography Geophysical Journal International (1998) 135, pp. 671 -690. [5] Chauris, H., Noble, M., Lambare, G., and P. Podvin, 2002. Migration velocity analysis from locally coherent events in 2-D laterally heterogeneous media, Part I: Theoretical aspects. Geophysics, Vol. 67 No. 4, pp. 1213-1224. [6] Guillaume, P., Audebert, F., Berthet, P., David, B., Herrenschmidt, A. and Zhang, X. 2001.3D Finite-offset tomographic inversion of CRP-scan data, with or without anisotropy, 71st SEG International Exhibition &Annual Meeting, pp. 31-34. [7] Lambare, G., 2008. Stereotomograpy, Geophysics Vol. 73 No. 5, pp. VE25-VE34. [8] Perrone, F., Guillaume, P., and Seher, T., 2015. “Dynamic-image Warping and Volumetric VpA / s Constraint for Nonlinear PP / OS Tomography, 77th EAGE Conference &Exhibition , Tu N116 05. [9] Thomsen, L., 1986, Weak elastic anisotropy, Geophysics Vol. 51 No. 10, pp. 1954-1966.

[10] Thomsen L., 1999, Converted-wave reflection seismology over inhomogeneous, anisotropic media, Geophysics, 64 No. 3, pp. 678-690.

[11] Tsvankin I. and Thomsen L., 1994. Nonhyperbolic reflection moveout in anisotropic media, Geophysics, Vol. 59 No. 8, pp. 1290-1304.

[12] Woodward, M.J., Nichols, D., Zdraveva, 0., Whitfield, P., and Johns, T., 2008, A decade of tomography, Geophysics, 73 No. 5, pp. VE5-VE11.

[13] Montel, J.P., and G. Lambare, 2019. Kinematics of common-image gathers — Part 1: Theory, Geophysics, Vol. 84 No. 5, pp. S437-S447.

[14] Montel, J.P., and G. Lambare, 2019. Kinematics of common-image gathers — Part 2: Tomographic ray tracing and applications, Geophysics, Vol 84 No. 5, pp. S501-S510.

Claims

25WHAT IS CLAIMED IS:

1. A seismic exploration method (1100) for locating targeted resources in a subsurface formation, the method comprising:obtaining (1110) PP and PS seismic datasets representing P-waves and S-waves emerging from the subsurface formation in response to injecting P-waves into the subsurface formation;extracting (1120) PP kinematic invariants and PS kinematic invariants from PP and PS seismic dataset, respectively; andobtaining (1130) a complete anisotropic model of the subsurface formation depending on P-wave velocity, Vp, S-wave velocity, Vs, and Thomsen anisotropy parameters e and 8, so that model seismic data predicted using the complete anisotropic model fulfills an output criterion, wherein the output criterion is based on an optimum value of the cost function representing a joint minimization of dRMOpp, dRMOps and AZpp.ps,wherein the obtaining of the complete anisotropic model of the subsurface formation includes performing a PP tomography, a PS tomography using the PP kinematic invariants and PS kinematic invariants, respectively, and a PP-PS data registration.

2. The method of claim 1, wherein the obtaining of the anisotropic model is not constrained using well information acquired from one or more wells dig into the subsurface formation.08 12 253. The method of claim 1, wherein the obtaining comprises:updating Vp and e of an initial anisotropic model by performing a first PP tomography using the PP kinematic invariants to obtain a modified anisotropic model with same IZ? and 6 values as the initial anisotropic model;updating I / s of the modified anisotropic model by performing the PS tomography using the PS kinematic invariants to obtain a complete anisotropic model with same Vp, e and 6 values as the modified anisotropic model;updating Vp of the complete anisotropic model by performing the PP and PS data registration matching reflector locations as reconstructed using P-waves and S-waves to obtain an updated anisotropic model with same Vs, e and 6 values as the modified anisotropic model; andupdating 6 of the updated anisotropic model by performing a second PP tomography to obtain a next-iteration anisotropic model usable as the initial model for repeating the first PP tomography followed by the PS tomography, the nextiteration anisotropic model having same Vp and values as the updated anisotropic model.

4. The method of claim 3, wherein the output criterion is based on convergence of a cost function depending on minimizing P-wave and S-wave residual moveouts and differences between the reflector locations as reconstructed using the P-waves and using the S-waves,the method further comprising:when the output criterion is fulfilled08 12 25outputting the complete anisotropic model as the complete anisotropic model of the subsurface formation,otherwise, performing the PP and PS data registration, the second PP tomography, the first PP tomography using the next-iteration anisotropic model as the initial anisotropic model, and the PS tomography.

5. The method of claim 4, wherein the first PP tomography comprises: starting from a current anisotropic model being the initial anisotropic model, until P-wave residual moveout, RMO, predicted based on the current anisotropic model is minimizedperforming a kinematic remigration of the PP kinematic invariants using the current anisotropic model,calculating the P-wave RMO, andif the P-wave RMO has been minimized, outputting the current anisotropic model as the modified anisotropic model, otherwise updating the current anisotropic model by varying vp and e thereof.

6. The method of claim 4, wherein the PS tomography comprises:starting from a current anisotropic model being the modified anisotropic model, until S-wave residual moveout, RMO, predicted based on the current anisotropic model is minimized,performing a kinematic remigration of the PS kinematic invariants using the current anisotropic model,08 12 25calculating the S-wave RMO, andif the S-wave RMO has been minimized, outputting the current anisotropic model as the complete anisotropic model, otherwise updating the current anisotropic model by varying thereof.

7. The method of claim 4, wherein the second PP tomography comprises: starting from a current anisotropic model being the updated anisotropic model, until P-wave residual moveout, RMO, predicted based on the current anisotropic model is minimized,performing a kinematic remigration of the PP kinematic invariants using the current anisotropic model,calculating the P-wave RMO, andif the P-wave RMO has been minimized, outputting the current anisotropic model as the next-iteration anisotropic model, otherwise update the current anisotropic model by varying vp thereof.

8. A seismic data processing apparatus (1200) for locating targeted resources in a subsurface formation, the apparatus comprising:a communication interface (1220) for obtaining PP and PS seismic datasets representing P-waves and S-waves emerging from the subsurface formation in response to injecting P-waves into the subsurface formation; anda processor (1210) connected to the communication interface and configured:to extract PP kinematic invariants and PS kinematic invariants from PP08 12 25and PS seismic dataset, respectively, andto obtain an anisotropic model of the subsurface formation depending on P-wave velocity, Vp, S-wave velocity, Vs, and Thomsen anisotropy parameters e and 8, so that model seismic data predicted using the anisotropic model fulfills an output criterion, wherein the output criterion is based on an optimum value of the cost function representing a joint minimization of dRMOpp, dRMOps and AZpp.ps, wherein the processor is configured to perform a PP tomography, a PS tomography using the PP kinematic invariants and PS kinematic invariants, respectively, and a PP and PS data registration for obtaining the anisotropic model.

9. The seismic data processing apparatus of claim 8, wherein the processor obtains the anisotropic model by:updating Vp and e of an initial anisotropic model via the first PP tomography using the PP kinematic invariants to obtain a modified anisotropic model with same 14 and 6 values as the initial anisotropic model;updating I / s of the modified anisotropic model via the PS tomography using the PS kinematic invariants to obtain a complete anisotropic model with same Vp, e and 6 values as the modified anisotropic model;updating Vp of the complete anisotropic model via the PP and PS data registration matching PP and PS reflector locations as reconstructed using P-waves08 12 25and S-waves to obtain an updated anisotropic model with same Vs, e and 8 values as the modified anisotropic model; andupdating 6 and e of the updated anisotropic model via a second PP tomography to obtain a next-iteration anisotropic model usable as the initial model for repeating the first PP tomography followed by the PS tomography, the nextiteration anisotropic model having same Vp and values as the updated anisotropic model.

10. The seismic data processing apparatus of claim 9, whereinthe output criterion is based on convergence of a cost function depending on minimizing P-wave and S-wave residual moveouts and differences between the reflector locations as reconstructed using the P-waves and using the S-waves, andthe processor is further configured, when the output criterion is fulfilled, to output the complete anisotropic model, otherwise, to repeat the PP and PS data registration, the second PPtomography, the first PP tomography using the next-iteration anisotropic model as the initial anisotropic model, and the PS tomography.

11. The seismic data processing apparatus of claim 9, wherein the processor is configured to perform the first PP tomography by:starting from a current anisotropic model being the initial anisotropic model, until P-wave residual moveout, RMO, predicted based on the current anisotropic model is minimized to:08 12 25perform a kinematic remigration of the PP kinematic invariants using the current anisotropic model,calculate the P-wave RMO, andif the P-wave RMO has been minimized, output the current anisotropic model as the modified anisotropic model, otherwise update the current anisotropic model by varying vp and e thereof.

12. The seismic data processing apparatus of claim 9, wherein the processor is configured to perform the PS tomography by:starting from a current anisotropic model being the modified anisotropic model, until S-wave residual moveout, RMO, predicted based on the current anisotropic model is minimized to:perform a kinematic remigration of the PS kinematic invariants using the current anisotropic model,calculate the S-wave RMO, andif the S-wave RMO has been minimized, output the current anisotropic model as the complete anisotropic model, otherwise update the current anisotropic model by varying v; thereof.

13. The seismic data processing apparatus of claim 9, wherein the processor is configured to perform the second PP tomography by:08 12 25starting from a current anisotropic model being the updated anisotropic model, until P-wave residual moveout, RMO, predicted based on the current anisotropic model is minimized to:perform a kinematic remigration of the PP kinematic invariants using the current anisotropic model,calculate the P-wave RMO, andif the P-wave RMO has been minimized, output the current anisotropic model as the next-iteration anisotropic model, otherwise update the current anisotropic model by varying vp thereof.

14. The seismic data processing apparatus of claim 8, wherein the obtaining of the anisotropic model is not constrained using well information acquired from one or more wells dig into the subsurface formation.

15. A computer-readable recording medium (1030) non-transitorily storing executable codes, which, when executed by a processor, make the processor perform a seismic explorations method (1100) for locating targeted resources in a subsurface formation, the method comprising:obtaining (1110) PP and PS seismic datasets representing P-waves and S-waves emerging from the subsurface formation in response to injecting P-waves into the subsurface formation;extracting (1120) PP kinematic invariants and PS kinematic invariants from PP and PS seismic dataset, respectively; and08 12 25obtaining (1130) a complete anisotropic model of the subsurface formation depending on P-wave velocity, Vp, S-wave velocity, Vs, and Thomsen anisotropy parameters e and 8, so that model seismic data predicted using the complete anisotropic model fulfills an output criterion, wherein the output criterion is based on an optimum value of the cost function representing a joint minimization of dRMOpp, dRMOps and AZpp.ps,wherein the obtaining of the complete anisotropic model of the subsurface formation includes performing a PP tomography, a PS tomography using the PP kinematic invariants and PS kinematic invariants, respectively, and a PP-PS data registration.

16. The computer-readable recording medium of claim 15, wherein the obtaining of the anisotropic model is not constrained using well information acquired from one or more wells dig into the subsurface formation.

17. The computer-readable recording medium of claim 15, wherein the obtaining comprises:updating Vp and e of an initial anisotropic model by performing a first PP tomography using the PP kinematic invariants to obtain a modified anisotropic model with same 14 and 6 values as the initial anisotropic model;updating I / s of the modified anisotropic model by performing the PS tomography using the PS kinematic invariants to obtain a complete anisotropic model with same Vp, e and 6 values as the modified anisotropic model;08 12 25updating Vp of the complete anisotropic model by performing the PP and PS data registration matching reflector locations as reconstructed using P-waves and S-waves to obtain an updated anisotropic model with same Vs, e and 6 values as the modified anisotropic model; andupdating 6 of the updated anisotropic model by performing a second PP tomography to obtain a next-iteration anisotropic model usable as the initial model for repeating the first PP tomography followed by the PS tomography, the nextiteration anisotropic model having same Vp and values as the updated anisotropic model.

18. The computer-readable recording medium of claim 17, wherein the output criterion is based on convergence of a cost function depending on minimizing P-wave and S-wave residual moveouts and differences between the reflector locations as reconstructed using the P-waves and using the S-waves,the method further comprising:when the output criterion is fulfilledoutputting the complete anisotropic model as the complete anisotropic model of the subsurface formation,otherwise, performing the PP and PS data registration, the second PP tomography, the first PP tomography using the next-iteration anisotropic model as the initial anisotropic model, and the PS tomography.08 12 2519. The computer-readable recording medium of claim 18, whereinthe first PP tomography comprises:starting from a current anisotropic model being the initial anisotropic model, until P-wave residual moveout, RMO, predicted based on the current anisotropic model is minimizedperforming a kinematic remigration of the PP kinematic invariants using the current anisotropic model,calculating the P-wave RMO, andif the P-wave RMO has been minimized, outputting the current anisotropic model as the modified anisotropic model, otherwise updating the current anisotropic model by varying vp and e thereof;the PS tomography comprises:starting from a current anisotropic model being the modified anisotropic model, until S-wave residual moveout, RMO, predicted based on the current anisotropic model is minimized,performing a kinematic remigration of the PS kinematic invariants using the current anisotropic model,calculating the S-wave RMO, andif the S-wave RMO has been minimized, outputting the current anisotropic model as the complete anisotropic model, otherwise updating the current anisotropic model by varying v; thereof; and the second PP tomography comprises:08 12 25starting from a current anisotropic model being the updated anisotropic model, until P-wave residual moveout, RMO, predicted based on the current anisotropic model is minimized,performing a kinematic remigration of the PP kinematic invariants using the current anisotropic model,calculating the P-wave RMO, andif the P-wave RMO has been minimized, outputting the current anisotropic model as the next-iteration anisotropic model, otherwise update the current anisotropic model by varying vp thereof.

20. The computer-readable recording medium of claim 15, wherein the obtaining of the anisotropic model is not constrained using well information acquired from one or more wells dig into the subsurface formation.

Images