Seismic Imaging Systems and Methods Employing a 3D Reverse Time Migration with Tilted Transverse Isotropy

Abstract

Seismic survey data is converted into a subsurface data volume that users can view to understand and analyze the subsurface structures in the survey region. At least some of the disclosed systems and methods employ three-dimensional reverse time migration using wave equations adapted for use in formations having tilted transverse isotropy. Relative to existing methods, the disclosed systems and methods rely on fewer approximations and suffer fewer limitations on the circumstances in which they can be employed. Moreover, because the disclosed wave equations are derived from Hooke's law (and consequently they operate on fields that are directly tied to physical quantities), they demonstrate an increased stability relative to existing methods. Survey data analysts employing the disclosed systems and methods should obtain better images of the subsurface and be better able to identify reservoirs and deposits for commercial exploitation.

Description

BACKGROUND

Seismology is used for exploration, archaeological studies, and engineering projects that require geological information. Exploration seismology provides data that, when used in conjunction with other available geophysical, borehole, and geological data, can provide information about the structure and distribution of rock types and their contents. Such information greatly aids searches for water, geothermal reservoirs, and mineral deposits such as hydrocarbon reservoirs and ore veins. Most oil companies rely on exploration seismology to select sites in which to drill exploratory oil wells.

Traditional seismology employs artificially generated seismic waves to map subsurface structures. The seismic waves propagate from a source down into the earth and reflect from boundaries between subsurface structures. Surface receivers detect and record reflected seismic waves for later analysis. Though some large-scale structures can often be perceived from a direct examination of the recorded signals, generally the recorded signals are processed to remove distortion and reveal finer detail in the subsurface image. Various existing processing methods do not sufficiently remove distortion, and further they can be unduly demanding in terms of computational resources. Improved systems and methods are disclosed herein.

BRIEF DESCRIPTION OF THE DRAWINGS

A better understanding of the various disclosed embodiments can be obtained when the following detailed description is considered in conjunction with the attached drawings, in which:

FIG. 1 shows an illustrative seismic survey environment;

FIG. 2 shows an illustrative seismic source and receiver geometry;

FIG. 3 shows an illustrative seismic survey recording system;

FIG. 4 shows illustrative seismic traces;

FIG. 5 shows an illustrative data volume in three dimensions;

FIG. 6 shows an illustrative imaging system; and

FIG. 7 shows a flowchart of an illustrative seismic imaging method.

While the invention is susceptible to various modifications and alternative forms, specific embodiments thereof are shown by way of example in the drawings and will herein be described in detail. It should be understood, however, that the drawings and detailed description are not intended to limit the disclosed embodiments to the particular forms shown, but on the contrary, the intention is to cover all modifications, equivalents and alternatives falling within the scope of the appended claims.

DETAILED DESCRIPTION

Disclosed herein are various systems and methods that convert seismic survey data into a subsurface data volume that users can view to understand and analyze the subsurface structures in the survey region. At least some of the disclosed systems and methods employ three-dimensional reverse time migration using wave equations adapted for use in formations having tilted transverse isotropy. Relative to existing methods, the disclosed systems and methods rely on fewer approximations and suffer fewer limitations on the circumstances in which they can be employed. Moreover, because the disclosed wave equations are derived from Hooke's law (and consequently they operate on fields that are directly tied to physical quantities), they demonstrate an increased stability relative to existing methods. Survey data analysts employing the disclosed systems and methods should obtain better images of the subsurface and be better able to identify reservoirs and deposits for commercial exploitation.

The disclosed systems and methods are best understood when described in an illustrative usage context. Accordingly, FIG. 1 shows an illustrative seismic survey environment, in which surveyors position an array of seismic receivers 102 in a spaced-apart arrangement on the earth's surface 104 to detect seismic waves. (The array typically extends for a couple of kilometers in each direction, and is moved periodically from place to place until the whole survey region has been covered.) The receivers 102 communicate wirelessly or via cable to a data acquisition unit 106 that receives, processes, and stores the seismic signal data collected by the receivers. The surveyors trigger a seismic energy source 108 (e.g., a vibrator truck) at multiple positions to generate seismic energy waves that propagate through the earth 110. Such waves reflect from acoustic impedance discontinuities to reach the receivers 102. Illustrative discontinuities include faults, boundaries between formation beds, and boundaries between formation fluids. The discontinuities appear as bright spots in the subsurface structure representation that is derived from the seismic signal data.

FIG. 1 further shows an illustrative subsurface model. In this model, the earth has three relatively flat formation layers and two dipping formation layers of varying composition and hence varying speeds of sound. Within each formation, the speed of sound can be isotropic (i.e., the same in every direction) or anisotropic. Due to the manner in which rocks are formed, nearly all anisotropic formations are transversely isotropic. In other words the speed of sound in anisotropic formations is the same in every “horizontal” direction, but is different for sound traveling in the “vertical” direction. Note, however, that geological activity can change formation orientations, turning a vertical transversely isotropic (“VTI”) formation into a tilted transversely isotropic (“TTI”) formation. In FIG. 1 shows the third flat layer as being VTI, while the first dipping formation layer is M.

FIG. 2 shows an illustrative geometry for source positions 202 and receiver positions 204, as they might be seen from an overhead view. Viable seismic survey geometries are infinite in variety, and can be gathered on a sector-by-sector basis, a rolling geometry basis, a moving-array basis, and various combinations thereof. The main message here is that the number of receiver signals acquired in response to each firing of the seismic source is fairly large, and when the number of source firings is taken into account, the resulting number of traces can easily reach into the millions.

FIG. 3 shows an illustrative seismic survey recording system having the receivers 102 coupled to a bus 302 to communicate digital signals to data recording circuitry 306. Position information sensors 304 (and optionally sensors for other parameters) are also coupled to the data recording circuitry 306 to enable the data recording circuitry to store additional information useful for interpreting the recorded data. Illustratively, such additional information can include source waveform characteristics, digitization settings, detected faults in the system, etc.

Recording circuitry 306 acquires the high speed data stream(s) from receivers 102 onto a nonvolatile storage medium such as a storage array of optical or magnetic disks. The data is stored in the form of (possibly compressed) seismic traces, each trace being the signal detected and sampled by a given receiver in response to a given shot. (The associated shot and receiver positions are also stored,) Illustrative seismic signals are shown in FIG. 4. The signals indicate some measure of seismic wave energy as a function of time (e.g., displacement, velocity, acceleration, pressure), and they are digitized at high resolution (e.g., 16 to $2 bits) at a programmable sampling rate (e.g., 400 to 1000 Hz) for a fixed duration after each shot (e.g., 30 seconds). Such signals can be grouped in different ways, and when so grouped, they are called a “gather”. For example, a “common midpoint gather” is the group of traces that have a midpoint within a defined region. A “shot gather” is the group of traces recorded for a single firing of the seismic source.

A general purpose data processing system 308 receives the acquired seismic survey data from the data recording circuitry 306. In some cases the general purpose data processing system 308 is physically coupled to the data recording circuitry and provides a way to configure the recording circuitry and perform preliminary processing in the field. More typically, however, the general purpose data processing system is located at a central computing facility with adequate computing resources for intensive processing. The survey data can be transported to the central facility on physical media or communicated via a computer network. Processing system 308 includes a user interface having a graphical display and a keyboard or other method of accepting user input, enabling users to view and analyze the subsurface structure images derived from the seismic survey data.

The recorded seismic survey data is of little use when maintained in the format of FIG. 4. Although it is possible to plot the various recorded waveforms side by side in a plot that reveals large scale subsurface structures, such structures are distorted and finer structures cannot even be seen. Hence the data is processed to create a data volume, i.e., a three dimensional array of data values such as that shown in FIG. 5. The data volume represents some seismic attribute throughout the survey region. The three-dimensional array comprises uniformly-sized cells, each cell having a data value representing the seismic attribute for that cell. Various seismic attributes may be represented, and in some embodiments, each cell has multiple data values to represent multiple seismic attributes. Examples of suitable seismic attributes include reflectivity, acoustic impedance, acoustic velocity, and density. The volumetric data format readily lends itself to computational analysis and visual rendering, and for this reason, the data volume may be termed a “three-dimensional image” of the survey region.

FIG. 6 shows an illustrative computer system 600 for performing seismic data processing including three-dimensional reverse time migration with tilted transverse isotropy. A personal workstation 602 is coupled via a local area network (LAN) 604 to one or more multi-processor computers 606, which are in turn coupled via the LAN to one or more shared storage units 608. Personal workstation 602 serves as a user interface to the processing system, enabling a user to load survey data into the system, to retrieve and view image data from the system, and to configure and monitor the operation of the processing system. Personal workstation 602 may take the form of a desktop computer with a graphical display that graphically shows survey data and 3D images of the survey region, and with a keyboard that enables the user to move files and execute processing software.

LAN 604 provides high-speed communication between multi-processor computers 606 and with personal workstation 602. The LAN 604 may take the form of an Ethernet network.

Multi-processor computer(s) 606 provide parallel processing capability to enable suitably prompt conversion of seismic trace signals into a survey region image. Each computer 606 includes multiple processors 612, distributed memory 614, an internal bus 616, and a LAN interface 620. Each processor 612 operates on an allocated portion of the input data to produce a partial image of the seismic survey region. Associated with each processor 612 is a distributed memory module 614 that stores conversion software and a working data set for the processor's use. Internal bus 616 provides inter-processor communication and communication to the LAN networks via interface 620. Communication between processors in different computers 606 can be provided by LAN 604.

Shared storage units 608 may be large, stand-alone information storage units that employ magnetic disk media for nonvolatile data storage. To improve data access speed and reliability, the shared storage units 608 may be configured as a redundant disk array. Shared storage units 608 initially store a velocity data volume and shot gathers from a seismic survey. The illumination matrix values and/or reflectivity image volumes can be stored on shared storage units 608 for later processing. In response to a request from the workstation 602, the image volume data can be retrieved by computers 606 and supplied to workstation for conversion to a graphical image to be displayed to a user.

In this context, we turn to a discussion of the methods for converting seismic survey data into three dimensional data volumes. The following discussion necessarily omits a significant amount of detail that would be known to one of ordinary skill in the art, but which can be found in the references cited throughout. Part of the discussion assumes that a velocity distribution of the subsurface is known beforehand, as well as the orientation and degree of anisotropy. Many ways exist for determining the velocity distribution from the seismic survey data, including the traditional approach described Jon F. Claerbout, Fundamentals of Geophysical Data Processing, p. 246-56, which is hereby incorporated herein by reference, and an iterative tomographic migration technique such as that described in U.S. application Ser. No. 12/566,885, “Seismic imaging Systems and Methods Employing Tomographic Migration-Velocity Analysis Using Common Angle Image Gathers” by Xia Fan et al, which is hereby incorporated herein by reference.

Given a volumetric velocity distribution, the reverse time migration technique consists of three stages: forward migration, reverse migration, and correlation. In this context, the term “migration” refers to the extrapolation of wave fields to determine, at each point in a volume, a wave intensity as a function of time. In the forward migration stage, a computer extrapolates a wave field from the shot into the volume, using a wave equation to determine how the seismic waves from the shot propagate forward in time. In the reverse migration stage, the computer uses the wave equation to extrapolate the wave field measured by the receivers into the volume, this time using the wave equation to determine how the seismic waves propagate in reverse, i.e., backward in time. In the third stage, the computer combines the two wave fields, in effect obtaining a zero-lag time correlation of the wave intensity function at each point in the volume. A high correlation indicates the presence of a reflector, and a low correlation indicates its absence.

The reverse time migration is repeated for each shot in the survey. The results are combined to build up a three-dimensional map of the subsurface. The quality of the map is strongly dependent on the precise manner in which migration is performed. In many cases, the wave equations that are employed to perform the wave field extrapolations employ assumptions and approximations in an effort to reduce computational requirements. Such assumptions and approximations can obscure or blur the features of the subsurface structure and/or cause those features to appear shifted from their true positions. Ideally, the wave equations used will be computationally tractable without suffering from such problems. Accordingly, there follows a detailed development and description of the wave equations and their usage in the seismic imaging systems and methods of the present disclosure. These equations are believed to alleviate such problems while remaining inexpensive in terms of required computational resources.

E. Duveneck, P. Milcik, P. M. Bakker, “Acoustic VII wave equations and their application for anisotropic reverse-time migration”, SEG Las Vegas 2008, p. 2186-2190, discloses three-dimensional wave equations derived based on Hooke's law, with an assumption of vertical transverse isotropy (VTI) in the subsurface formations:

$\begin{array}{cc}\frac{{\partial }^2{\sigma }_H}{\partial t^2}=V_p^2\left\{\left(1+2\varepsilon \right)\left(\frac{{\partial }^2{\sigma }_H}{\partial x^2}+\frac{{\partial }^2{\sigma }_H}{\partial y^2}\right)+\sqrt{1+2\delta }\frac{{\partial }^2{\sigma }_V}{\partial z^2}\right\}& \left(1\right)\\ \frac{{\partial }^2{\sigma }_V}{\partial t^2}=V_p^2\left\{\left(1+2\delta \right)\left(\frac{{\partial }^2{\sigma }_H}{\partial x^2}+\frac{{\partial }^2{\sigma }_H}{\partial y^2}\right)+\frac{{\partial }^2{\sigma }_V}{\partial z^2}\right\}& \left(2\right)\end{array}$

where σ_H is the horizontal stress field, σ_V is the vertical stress field, V_p is the p-wave velocity, ε and δ are the Thomsen anisotropy parameters (see L. Thomsen, “Weak elastic anisotropy”, Geophysics 51, p. 1954-1966, 1986), and x, y, z, and t are the normal space and time coordinates for the data volume. To account for tilted transverse isotropy (TTI), assume that a VTI formation is aligned with the subject formation by two rotations—a first counterclockwise rotation about the y axis by angle θ and a second counterclockwise rotation about the x axis by angle φ—respectively represented by the following rotation matrices A and B:

$\begin{array}{cc}A=\left[\begin{array}{ccc}\cos \theta & 0& \sin \theta \\ 0& 1& 0\\ -\sin \theta & 0& \cos \theta \end{array}\right]& \left(3\right)\\ B=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \phi & \sin \phi \\ 0& -\sin \phi & \cos \phi \end{array}\right]& \left(4\right)\end{array}$

The resulting 3D wave equations are:

$\begin{array}{cc}\frac{{\partial }^2{\sigma }_H}{\partial t^2}=V_p^2\left\{\begin{array}{c}\left(1+2\varepsilon \right)\left[\begin{array}{c}\begin{array}{c}{\cos }^2\theta \frac{{\partial }^2{\sigma }_H}{\partial x^2}+\left({\sin }^2\theta {\sin }^2\phi +{\cos }^2\phi \right)\frac{{\partial }^2{\sigma }_H}{\partial y^2}+\\ \left({\sin }^2{\mathrm{\theta cos}}^2\phi +{\sin }^2\phi \right)\frac{{\partial }^2{\sigma }_H}{\partial z^2}-{\cos }^2\theta \sin 2\phi \frac{{\partial }^2{\sigma }_H}{\partial y\partial z}-\end{array}\\ \sin 2\theta \cos \phi \frac{{\partial }^2{\sigma }_H}{\partial x\partial z}-\sin 2\theta \sin \phi \frac{{\partial }^2{\sigma }_H}{\partial x\partial y}\end{array}\right]+\\ \sqrt{1+2\delta }\left[\begin{array}{c}\begin{array}{c}{\sin }^2\theta \frac{{\partial }^2{\sigma }_V}{\partial x^2}+{\cos }^2\theta {\sin }^2\phi \frac{{\partial }^2{\sigma }_V}{\partial y^2}+\\ {\cos }^2{\mathrm{\theta cos}}^2\phi \frac{{\partial }^2{\sigma }_V}{\partial z^2}+{\cos }^2\theta \sin 2\phi \frac{{\partial }^2{\sigma }_V}{\partial y\partial z}+\end{array}\\ \sin 2\mathrm{\theta cos\phi }\frac{{\partial }^2{\sigma }_V}{\partial x\partial z}+\sin 2\theta \sin \phi \frac{{\partial }^2{\sigma }_V}{\partial x\partial y}\end{array}\right]\end{array}\right\}& \left(5\right)\\ \frac{{\partial }^2{\sigma }_V}{\partial t^2}=V_p^2\left\{\begin{array}{c}\sqrt{1+2\delta }\left[\begin{array}{c}\begin{array}{c}{\cos }^2\theta \frac{{\partial }^2{\sigma }_H}{\partial x^2}+\left({\sin }^2\theta {\sin }^2\phi +{\cos }^2\phi \right)\frac{{\partial }^2{\sigma }_H}{\partial y^2}+\\ \left({\sin }^2\theta {\cos }^2\phi +{\sin }^2\phi \right)\frac{{\partial }^2{\sigma }_H}{\partial z^2}-{\cos }^2\mathrm{\theta sin}2\phi \frac{{\partial }^2{\sigma }_H}{\partial y\partial z}-\end{array}\\ \sin 2\mathrm{\theta cos}\phi \frac{{\partial }^2{\sigma }_H}{\partial x\partial z}-\sin 2\mathrm{\theta sin}\phi \frac{{\partial }^2{\sigma }_H}{\partial x\partial y}\end{array}\right]+\\ \left[\begin{array}{c}\begin{array}{c}{\sin }^2\theta \frac{{\partial }^2{\sigma }_V}{\partial x^2}+{\cos }^2{\mathrm{\theta sin}}^2\phi \frac{{\partial }^2{\sigma }_V}{\partial y^2}+\\ {\cos }^2{\mathrm{\theta cos}}^2\phi \frac{{\partial }^2{\sigma }_V}{\partial z^2}+{\cos }^2\mathrm{\theta sin\phi }\frac{{\partial }^2{\sigma }_V}{\partial y\partial z}+\end{array}\\ \sin 2\mathrm{\theta cos\phi }\frac{{\partial }^2{\sigma }_V}{\partial x\partial z}+\sin 2\mathrm{\theta sin\phi }\frac{{\partial }^2{\sigma }_V}{\partial x\partial y}\end{array}\right]\end{array}\right\}& \left(6\right)\end{array}$

Various techniques can be used to apply these wave equations for migration in a discretized volume. Three are discussed here: finite difference (FD), pseudospectral (PS), and hybrid. In the finite difference technique, the spatial second derivatives are calculated using high-order difference equations. See, e.g., M. A. Dablain, “The application of high-order differencing to the scalar wave equation”, Geophysics vol. 51, No. 1, January 1986 p. 54-66. Thus, for example:

$\begin{array}{cc}\frac{{\partial }^2\sigma \left(x\right)}{\partial x^2}=\Delta x^{-2}*\left\{-2.92722222\sigma \left(x\right)+1.66666667\left[\sigma \left(x-1\right)+\sigma \left(x+1\right)\right]-0.23809524\left[\sigma \left(x-2\right)+\sigma \left(x+2\right)\right]+0.03968254[\sigma \left(x-3\right)+\sigma \left(x+3\right)]-0.00496031[\sigma \left(x-4\right)+\sigma \left(x+4\right)]+0.00031746\left[\sigma \left(x-5\right)+\sigma \left(x+5\right)\right]\right\}& \left(7\right)\end{array}$

in a tenth order approximation. For the cross-derivatives, a second order finite difference has been found to be adequate, e.g.:

$\begin{array}{cc}\frac{{\partial }^2\sigma \left(x,y\right)}{\partial x\partial y}=\left[\sigma \left(x-1,y-1\right)-\sigma \left(x-1,y+1\right)-\sigma \left(x+1,y-1\right)+\sigma \left(x+1,y+1\right)\right]/\left(4\Delta x\Delta y\right)& \left(8\right)\end{array}$

The spatial derivatives are repeatedly calculated for each point in the volume to step forward or backward in time in accordance with equations (5) and (6).

In the pseudospectral technique, the spatial second derivatives are calculated using the spatial Fourier Transform to find the transform coefficients for each wavenumber, multiplying each transform coefficient by a corresponding value as explained below, and taking the inverse transform. For each partial derivative along a given axis, the transform coefficients c_n are multiplied by ikn, where i=√{square root over (−1)}, k=2π/N, n=1, 2, . . . , N, and N is the number of coefficients of the spatial transform along that axis. The pseudospectral technique employs a spatial transform (and inverse transform) along each of the axes. It yields a simultaneous determination of the spatial second derivative at each point in the volume.

To determine the spatial second derivatives, the hybrid technique uses a combination of the pseudospectral technique and the finite difference technique. The pseudospectral technique is applied to determine partial derivatives along the horizontal (x- and y-axes). The finite-differencing technique is used to determine the partial derivatives along the vertical (z-axis). These operations are applied sequentially to determine the spatial second derivative at each point in the volume.

In each case, the calculated spatial derivatives can then be combined with the previously estimated distributions of velocity and anisotropy to determine the time progression of the horizontal and vertical stress fields in accordance with equations (5) and (6). Systems employing these equations do not suffer any dip limitations. There are, however, some stability considerations.

Equations (5) and (6) need to be recalculated for each cell at each time step in the shot interval. The larger the cells and the time step, the fewer the number of recalculations will be. However, if the cell dimensions or time step becomes too large, the numerical errors in the calculation will dominate, making the results meaningless. The cell dimensions in the horizontal direction (Δx and Δy) should be limited to less than:

$\begin{array}{cc}\Delta x\le \frac{V_{pmin}}{S*f_{\mathrm{Nyquist}}}& \left(9\right)\end{array}$

where V_max is the smallest wave velocity, f is the Nyquist frequency for the seismic signals, and S is a empirically determined value which is less than 10 and Which generally decreases for higher order differencing schemes. See, e.g., M. A. Dablain, “The application of high-order differencing to the scalar wave equation”, Geophysics, vol. 51, no. 1 (January 1986). For 10^th order FD applied to equations (5) and (6), S=4 is adequate to avoid grid dispersion. The optimal time step has been found to be:

$\begin{array}{cc}\Delta t\le \frac{2\Delta d_{min}}{\sqrt{W}V_{p\mathrm{ma}x}}{\left\{\begin{array}{c}\left(1+2{\varepsilon }_{m\mathrm{ax}}\right)\left[2-\left({\cos }^2\theta \sin 2\phi +\sin 2\theta \left(\sin \phi +\cos \phi \right)\right)\right]+\\ \sqrt{1+2{\delta }_{m\mathrm{ax}}}\left[1+\left({\cos }^2\theta \sin 2\phi +\sin 2\theta \left(\sin \phi +\cos \phi \right)\right)\right]\end{array}\right\}}^{-1/2}& \left(10\right)\end{array}$

where d_min is the smallest dimension of a cell in the discretized volume, V_max is the highest wave velocity, W is the sum of the absolute values of the finite difference equation coefficients (for 10^th order FD, W=6.826), ε_max and δ_max are the maximum anisotropy parameters, and the combination of θ and φ is that combination taken from the discretized volume that yields the smallest time step, (For determining stability for the pseudospectral technique, the √{square root over (W)} should be replaced by π.)

Another consideration is the effect of the boundaries on the migration process. Recall that acoustic impedance mismatches cause reflections, so if the boundaries of the region are not treated carefully, extraneous reflections will adversely affect the image. In at least some embodiments, the disclosed systems and methods employ an “absorbing boundary condition” technique such as that taught by Cerjan, et al., “A nonreflecting boundary condition for discrete acoustic and elastic wave equations”, Geophysics v.50, n.4, p. 705-708 April 1985. In this technique, a 20-cell thick layer is provided along each boundary. Stress field values in this layer are multiplied at each time step by an attenuation coefficient that varies with distance from the boundary to reduce reflections from the boundary.

In other embodiments, the disclosed systems and methods employ a “perfectly matched layer” technique such as that taught by Komatitsch and Tromp, “A perfectly matched layer absorbing boundary condition for the second-order seismic wave equation”, Geophys. J. Int. v.154 p. 146-153 (2003). In this technique, the wave equation used in the boundary layers is modified to include a damping factor for migration perpendicular to the boundary. To simplify the calculation, Komatitseh and Tromp split the equation into components. Adapting that approach to the VTI equations (1) and (2) above, the inventors have found that for damping along the x-axis:

σ_H=u₁+u₂+u₃, and  (11)

σ_V=u₄+u₅+u₆  (12)

where u₁ through u₆ are determined with the following equations:

$\begin{array}{cc}{\left(\partial t+D\left(x\right)\right)}^2u_1=V^2\left(1+2\varepsilon \right)\frac{{\partial }^2{\sigma }_H}{\partial x^2}& \left(13\right)\\ {\left(\partial t+D\left(x\right)\right)}^3u_2=-V^2D\left(x\right)\left(1+2\varepsilon \right)\frac{\partial {\sigma }_H}{\partial x}& \left(14\right)\\ \partial t^2u_3=V^2\left(\left(1+2\varepsilon \right)\frac{{\partial }^2{\sigma }_H}{\partial y^2}+\sqrt{1+2\delta }\frac{{\partial }^2{\sigma }_V}{\partial z^2}\right)& \left(15\right)\\ {\left(\partial t+D\left(x\right)\right)}^2u_4=V^2\sqrt{1+2\delta }\frac{{\partial }^2{\sigma }_H}{\partial x^2}& \left(13\right)\\ {\left(\partial t+D\left(x\right)\right)}^3u_5=-V^2D\left(x\right)\sqrt{1+2\delta }\frac{\partial {\sigma }_H}{\partial x}& \left(14\right)\\ \partial t^2u_6=V^2\left(\sqrt{1+2\delta }\frac{{\partial }^2{\sigma }_H}{\partial y^2}+\frac{{\partial }^2{\sigma }_V}{\partial z^2}\right).& \left(15\right)\end{array}$

D(x) is the damping factor that increases as you progress into the boundary layer:

$\begin{array}{cc}D\left(x\right)=-\log \left(R\right)*\frac{4{\mathrm{Vx}}^4}{2{\delta }^5}& \left(16\right)\end{array}$

where R is the reflection coefficient (typically 10⁻³), V is wave velocity, and δ is the thickness of the boundary layer.

Similar equations provide attenuation for boundaries perpendicular to the y and z axes. For overlaps in the boundary layers (e.g., at the coders and edges of the discretized volume), the horizontal and vertical stress values are calculated with each set of equations applicable to the boundaries and then averaged together.

Yet another consideration during migration is the effect of anisotropy on the source field. The disclosed systems and methods can minimize source field artifacts by assuming a small isotropic region around the shot. A region of about 20 cells in diameter has been found to be adequate for most purposes.

Given the foregoing discussion, we turn to FIG. 7 which outlines an illustrative seismic imaging method. In block 702, the seismic imaging system (e.g., a computer running seismic imaging software) retrieves configuration information and seismic trace data from a seismic survey database. The configuration information includes the dimensions of the subsurface volume to be imaged, the granularity of the cells (or grid points) in the volume, the size of the time steps, the locations of the shots and receivers, etc. The seismic trace data itself is generally stored in a structured form in a disk array to enable high bandwidth accessibility by a parallel processing computer system. In block 704 the seismic imaging system retrieves the velocity distribution model, which can include the amount and orientation of formation anisotropy.

In block 706, the system employs the three dimensional tilted transverse isotropy (3D TTI) equations to perform forward migration of a seismic shot, thereby determining for each grid point a horizontal and vertical stress as a function of time, σ_Hs(t) and σ_Vs(t) respectively. The system can use the finite difference technique, the pseudospectral technique, a hybrid technique, or any other suitable migration technique to do this determination.

In block 708, the system employs the 3D TTI equations to perform reverse time migration of the seismic traces from the receivers into the volume, thereby determining for each grid point a second horizontal and vertical stress as a function of time, σ_Hr(t) and σ_Vr(t) respectively. In block 710, the horizontal and vertical stress functions at each grid point are combined to determine an intensity at each grid point, e.g.:

I _H=∫σ_Hs(t)σ_Hr(t)dt  (17)

I _V=∫σ_Vs(t)σ_Vr(t)dt  (18)

I=√{square root over (I_H²+I_V²)}  (19)

And these intensities are added to the intensity values for the cumulative image in block 712.

In block 714, the system checks to determine if more shots exist in the survey data, and if so, the operations of blocks 706-712 are repeated with the next shot. Otherwise, in block 716 the system applies optional filters and makes the image available for a user to view. In some embodiments, the system applies a Laplacian filter as described in Yoon and Zhou, “Depth imaging with Multiples”, Geophysics v.66, n. 1, p. 246-255 (2001), to eliminate low-frequency artifacts and other artifacts of the migration process:

$\begin{array}{cc}{I}_{\mathrm{out}}\left(x,y,z\right)=\frac{{\partial }^2I_{in}}{\partial x^2}+\frac{{\partial }^2I_{in}}{\partial y^2}+\frac{{\partial }^2I_{in}}{\partial z^2}& \left(20\right)\end{array}$

In block 718 the system enables the user to interact with the image to manipulate it and gain a thorough understanding of the subsurface structure and hopefully identify reservoirs that can be developed. In at least some embodiments, the system user is able to apply different filters, rotate and slice the data volume, change colors and opacities, and manipulate the image in other ways as well.

It is contemplated that the operations described in FIG. 7 would be implemented in the form of software, which can be stored in computer memory, in long-term storage media, and/or in portable information storage media. It should be noted that illustrative method of FIG. 7 is provided as an explanatory aid. In practice, the various operations shown in FIG. 7 may be performed in different orders and are not necessarily sequential. For example, seismic data processing can benefit substantially from parallelism. In some processing method embodiments, data from different shots or survey regions can generally be processed independently. In other embodiments, the operations may be “pipelined” so that individual calculation operations are determined in the sequence shown, even though the operations are all occurring concurrently. Additional operations may be added to the illustrative method and/or several of the operations shown may be omitted.

Numerous variations and modifications will become apparent to those skilled in the art once the above disclosure is fully appreciated. It is intended that the following claims be interpreted to embrace all such variations and modifications.

Claims

1. A seismic imaging method that comprises: Obtaining at least one shot gather for a survey region;performing reverse time migration of the shot gather to determine horizontal and vertical stresses throughout a survey region while accounting for formations having tilted transverse isotropy; andconstructing an image based at least in part on said horizontal and vertical stresses.

2. The method of claim 1, further comprising displaying said image.

3. The method of claim 1, further comprising writing said image to an information storage medium.

4. The method of claim 1, further comprising: carrying out forward migration of a shot corresponding to said shot gather to determine horizontal and vertical stresses throughout a survey region while accounting for formations having tilted transverse isotropy,wherein said constructing includes combining said horizontal and vertical stresses from said forward and reverse time migrations.

5. The method of claim 4, wherein said constructing further comprises: repeating said performing, carrying out and combining operations for other shot gathers; andaccumulating results from each shot gather to form said image.

6. The method of claim 4, wherein said constructing further comprises applying a filter to correlations of the horizontal and vertical stresses from the forward and reverse time migrations.

7. The method of claim 4, wherein said forward migration employs an isotropic region around the shot to reduce migration artifacts.

8. The method of claim 1, wherein said migration employs an absorbing boundary layer to reduce extraneous reflections.

9. The method of claim 1, wherein said migration employs a matched boundary layer to reduce extraneous reflections.

10. The method of claim 1, wherein said migration is performed in three dimensions.

11. A seismic imaging system that comprises: at least one storage device that stores at least one shot gather from a seismic survey of a given survey region;a memory that stores seismic imaging software; andat least one processor coupled to the memory to execute the seismic imaging software, wherein the software configures the at least one processor to: retrieve said at least one shot gather from the storage device;perform reverse time migration of the shot gather to determine horizontal and vertical stresses throughout a survey region while accounting for formations having tilted transverse isotropy; andconstruct an image based at least in part on said horizontal and vertical stresses.

12. The system of claim 11, wherein the software further configures the at least one processor to display said image.

13. The system of claim 11, wherein the software further configures the at least one processor to write said image to an information storage medium.

14. The system of claim 11, wherein the software further configures the at least one processor to: implement forward migration of a shot corresponding to said shot gather to determine horizontal and vertical stresses throughout a survey region while accounting for formations having tilted transverse isotropy, andwherein as part of said constructing, the software configures the processor to combine said horizontal and vertical stresses from said forward and reverse time migrations.

15. The system of claim 14, wherein as part of said constructing the software further configures the processor to: repeat said performing, implementing, and combining operations for other shot gathers; andaccumulate results from each shot gather to form said image.

16. The system of claim 14, wherein as part of said constructing the software further configures the processor to apply a filter to correlations of the horizontal and vertical stresses from the forward and reverse time migrations.

17. The system of claim 14, wherein as part of said forward migration the software employs an isotropic region around the shot to reduce migration artifacts.

18. The system of claim 11, wherein as part of said migration, the software employs an absorbing boundary layer to reduce extraneous reflections.

19. The system of claim 11, wherein as part of said migration, the software employs a matched boundary layer to reduce extraneous reflections.

20. The system of claim 11, wherein the processor performs said migration in three dimensions.