WO2002065159A1

WO2002065159A1 - Method of seismic imaging using direct travel time computing

Info

Publication number: WO2002065159A1
Application number: PCT/KR2001/000215
Authority: WO
Inventors: Hae-Ryong Lim; Chang-Soo Shin
Original assignee: Hae-Ryong Lim; Chang-Soo Shin
Priority date: 2001-02-14
Filing date: 2001-02-14
Publication date: 2002-08-22

Abstract

The present invention provides a method of fast computing two-dimensional or three-dimensional travel-times for an arbitrary velocity model in the fields of geophysical data processing. The present invention simply egts rid of the major drawbacks of the shooting ray tracing technique by computing the travel-times with velocities and lengths from a source location to grid points, for a velocity model that a wave is transmitted from the source location to the grid points in straight ray in rectangular coordinates, and transforming the rectangular coordinates into polar coordinates. In the transformed velocity model, all ray paths are fixed with normal incident rays to a radial direction. Therefore, the travel-times and the amplitudes can be easily computed in a manner that the length or unit distance between grids is divided by a specific velocity between cells or grids. The directly transmitted travel-time can be substituted into a maximum energy travel-times in seismic imaging, depth migration. This method is executed on a pc based computer, a main frame or a massively parallel computer because it can compute the travel-time used for the depth migration with a short computing time and small data storage caèacity having no effect on the resolution of image.

Description

METHOD OF SEISMIC IMAGING USING DIRECT TRAVEL TIME COMPUTING

Technical Field

The present invention relates to the exploration of geophysics using seismic signals. More particularly, the present invention relates to a method of computing the direct travel-time of a seismic wave in rectangular coordinates or polar coordinates for an arbitrary velocity model to use the travel-time in a two- dimensional or three-dimensional depth migration.

Background Art

In general, a seismic wave is imaged in a manner that prestack depth migration is performed on the basis of the first inaccurate velocity model and the first image produced by the prestack depth migration is analyzed to correct the velocity model. This procedure is repeated until a satisfactory underground image is obtained. This process requires iterative computations of the travel-time for the corrected velocity models. Accordingly, a lot of computation time and computer memories are needed for calculation of the travel-time, resulting in large amount of computing cost.

There are currently three general methods of computing travel-time tables: ray tracing method, finite difference eikonal equation solver and graph theory method. These methods mainly calculate the minimum travel-times or first arrival times. Prior examples are disclosed in the following United States patents - U.S. PAT. No. 5265068, U.S. PAT. No. 5394325, U.S. PAT. No. 5724310, U.S. PAT. No. 5991695, U.S. PAT. No. 6018499.

Quotated from recent studies (Geoltrain and Brae, 1993, "Can we image complex structures with first-arrival travel-time", Geophysics, 58, p564-575; Nichols, 1996, "maximum energy travel-times calculated in the seismic frequency band", Geophysics, 61, p253-263), maximum energy arrival time gives better images in depth migration rather than first arrival times. However, it is very difficult to calculate maximum energy arrival times and a long period of time is required for the calculation. Furthermore, it is expensive to calculate the maximum energy arrival time because it needs a memory for storing a vast quantity of data and a CPU for processing the data.

Disclosure of Invention

Accordingly, it is an object of the present invention to provide a method of calculating travel-times using intermediate travel-times between the first arrival times and the maximum energy arrival times. The present invention can produce improved images, calculate travel times within a short period of time, and reduce storage space of data, compared with the conventional method of computing the first arrival times. In addition, the invention generates images similar to those produces by the method of computing the maximum energy arrival times while easily calculating the direct travel-times. The shooting ray-tracing method is mainly associated with a directly transmitted wave, therefore calculating the travel-times close to a maximum energy arrival event. However, it suffers from a caustic problem and an interpolation problem in the shadow zone, which requires a fine division of the ray take off angles. This takes a lot of time and cost. Furthermore, even with the fine division of the ray take off angles, the travel-times may not be calculated around the shadow zone.

The present invention simply gets rid of the major drawbacks of the shooting ray tracing technique by computing the travel-times with velocities and lengths from a source location to grid points, for a velocity model that a wave is transmitted from the source location to the grid points in straight ray in rectangular coordinates, and transforming the rectangular coordinates into polar coordinates. In the transformed velocity model, all ray paths are fixed with normal incident rays to a radial direction. Therefore, the travel-times and the amplitudes can be easily computed in a manner that the length or unit distance between grids is divided by a specific velocity between cells or grids. The calculated amplitudes is used as weighted values in the computation of the depth migration so as to be employed for weighted Kirchhoff prestack depth migration that considers even the amplitudes the measured seismic wave.

In case where the travel-times of the present invention are used for the depth migration of seismic wave exploration, the time required for computing the travel-times is shortened and data storage capacity is reduced while the resolution of image is not deteriorated, compared to the method of using the maximum energy arrival times. In addition, the resolution of image becomes close to the actual model, the time required for computing the arrival times is reduced, and data storage capacity becomes small, compared with the method of using the first arrival times.

Brief Description of the Drawings

Further objects and advantages of the invention can be more fully understood from the following detailed description taken in conjunction with the accompanying drawing in which: FIG. 1 is a flow chart outlining the steps of the method according to the present invention;

FIG. 2A depicts a method of computing direct travel-times in rectangular coordinates for a cell-based velocity model in the rectangular coordinates according to the present invention;

FIG. 2B depicts a method of computing the direct travel-times in the rectangular coordinates for a grid-based velocity model in the rectangular coordinates according to the present invention; FIG. 3 depicts a method of computing the direct travel-times when transforming the velocity model in the rectangular coordinates into the polar or the spherical coordinates according to the present invention;

FIG. 4A shows travel-time contours computed in the rectangular coordinate according to the present invention, which are superimposed on the velocity model for a vertical slice in the x direction through the SEG/EAEG Salt Model, where a point source is located at the surface;

FIG. 4B shows travel-time contours computed in the rectangular coordinate according to the present invention, which are superimposed on the velocity model for a vertical slice in the y direction through the SEG/EAEG Salt Model, where a point source is located at the surface;

FIG. 4C shows travel-time contours computed in the rectangular coordinate direct computation according to the present invention, which are superimposed on the velocity model for a horizontal slice at the respective depth through the SEG/EAEG Salt Model, where a point source is located at the surface; FIG. 5A shows travel-time contours computed in the spherical coordinates using the velocity model transform of the present invention, which are superimposed on the velocity model for a vertical slice in the x direction through the SEG/EAEG Salt Model;

FIG. 5B shows travel-time contours computed in the spherical coordinates using the velocity model transform of the present invention, which are superimposed on the velocity model for a vertical slice in the y direction through the SEG/EAEG Salt Model; FIG. 5C shows travel-time contours computed in the spherical coordinates using the velocity model transform of the present invention, which are superimposed on the velocity model for a horizontal slice at the respective depth through the SEG/EAEG Salt Model; FIG. 6A shows a vertical slice of the velocity model in the y direction through the SEG/EAEG Salt Model;

FIG. 6B shows a vertical slice of the final depth migrated image when imaging a shot-line data of SEG/EAEG Salt Model C3-NA data set using the travel-time tables -generated by the present invention; FIG. 7 depicts the velocity model transformed into polar coordinate and ray tracing;

FIG. 8A shows an original four-layered model;

FIG. 8B shows the transformed velocity model, where Δr=70m and Δθ=10 degree; FIG. 8C shows the transformed velocity model, where Δr=7m and Δθ=l degree;

FIG. 8D shows the transformed velocity model, where Δr=0.7m and Δθ=0.1 degree;

FIG. 9 shows the direct arrival times computed by the present invention (black solid line), first arrival times by Vidale's method (black dashed line) and the maximum energy arrival times picked from the wave fields by FEM modeling

(white solid line) superimposed on a snapshot image obtained by FEM modeling for the well-known two dimensional Marmousi Model;

FIG. 10A shows the migrated image for Marmousi Model data using the first arrival times obtained by Vidale's method; and

FIG. 10B shows the migrated image for Marmousi Model data using the direct arrival times according to the present invention. Best Mode for Carrying Out the Invention

The present invention will now be described in detail in connection with preferred embodiments with reference to the accompanying drawings. For reference, like reference characters designate corresponding parts throughout several views.

FIG. 1 shows a flow chart outlining the steps of the method of the invention. If a velocity model in the rectangular coordinates is given (step 101), a travel-time from a given source (step 102) to an underground grid is computed. In general, the velocity model is a 3-D cube array in the rectangular coordinates, describing the distribution of velocities on a grid (grid-based velocity model) or within a cell (cell-based velocity model). In the grid-based model, it is assumed that the velocity between grids can be interpolated using the velocities of the neighboring grids. In the cell-based model, it is assumed that the velocity within a cell is constant.

The present invention computes a travel-time from a source to a grid using a straight ray, transforming the velocity model. There are two ways of the invention in this step. One is directly computing the travel-time in the rectangular coordinates (step 104) and the other is transforming the velocity model into the spherical coordinates using an interpolation (step 105), computing the travel-time at each grid in the spherical coordinates (step 106) and then transforming the travel-time in the spherical coordinates into the rectangular coordinates (step 107).

FIG. 2A depicts the method of directly computing the travel-time in the rectangular coordinates for a cell-based velocity model. A travel-time from a source to a grid j, T_j, is represented as T_J = ∑ _l /V_l (1)

where Lj is the raypath length in the 1^th cell and V_} is the velocity value of the 1^th cell.

FIG. 2B depicts the method of directly computing the travel-time in the rectangular coordinates for a grid-based velocity model. A travel-time from a source to a grid, j, T_j, is represented as

T_j = ^1, 1/7, + VV_M)/2 (2)

where Li is the unit raypath length when the straight line connecting the source and grid j is divided into a predetermined number of sections, and Vj and

Vi₊i are velocities at both end points of the unit raypath, which can be determined by the interpolation of the velocities of the adjacent grids.

In the case of using the velocity model transform into the spherical coordinates, first the velocity model in the rectangular coordinated is transformed to the spherical coordinates having an origin at the source point (step 105). Next, the travel-time at each grid is computed in the spherical coordinates.

FIG. 3 depicts the method of computing the travel-time in the spherical coordinates (step 106). Although two-dimensional polar grids (r,θ) are presented in this figure, the polar coordinates are transformed into the spherical coordinates (r,θ,φ) in case of the three-dimensional velocity model so that there is no significant difference in computing the travel-times between the two-dimensional and three-dimensional models. A travel-time from a source to a grid j, T_j, is given as

E₇ = ∑Δr(l/ , + l/ ;₊₁)/2 (3)

where V; and V_i+ι are velocities at grid points located on the straight line connecting the source and the grid j, and Δr and Δθ are lengths between grids.

After computing the travel-times at all grid points in the spherical coordinates (r, θ,φ), then they are transformed back to the rectangular grids using an interpolation (step 107). After generating a travel-time table for a given source location, then this table is stored as a file in a hard disk or transferred being stored in a memory to the depth migration step (step 108). The steps 101 to 108 are iteratively performed for the next source location (step 109). Upon completion of computing the travel- times for all source locations, computations for new velocity model are iteratively carried out (step 110). The conventional travel-time computing methods store the computed travel-times as files on hard disks, and then read them to perform the depth migration because they need large computer memory capacities and long computing time. However, the present invention requires a small quantity of memory capacity and computing time shorter than the time required for writing the computed three-dimensional travel-times into files and reading them. Accordingly, the travel-time computing process can be made into a sub-program and called by a depth migration program. In other words, the travel-time tables can be directly transferred, not being in the form of file on a hard disk but being stored in a memory, to the depth migration step.

EMBODIMENTS

To verify the suitability of the invention, an experiment on the 3-D SEG/EAEG Salt Model (Aminzadeh, 1995, "3-D Modeling Project: Third Report", The Leading Edge 14, i 25- 128) was carried out. The overall model size is 13.5x13.5x4.2km on a 20m grid. FIGS. 4A-5C show the computed travel-times superimposed on the velocity model. A point source is activated at the surface (x=7220, y=6740) and the 201x201x210 3-D travel times are obtained with 20m spacmg.

FIGS. 4A-4C are obtained by the direct computation in the rectangular coordinates (step 104) and FIGS. 5A-5C are obtained by using the spherical coordinate transform (steps 105 to 107). They show nearly the same value without small discrepancies from an interpolation error in the spherical coordinate. But this small error gives no effects to the depth migration process.

The prestack depth migration was performed for one shot-line data of SEG/EAEG Salt Model C3-NA data set using the travel-time tables generated by the method of the present invention. FIG. 6A shows a vertical slice of the velocity model in the y direction through the SEG/EAEG Salt Model, and FIG. 6B shows a vertical slice of the final depth migrated image when imaging a shot-line data of SEG/EAEG Salt Model C3-NA data set using the travel time tables generated by the present invention. We note that the final migrated image has a good agreement with the original model. The theoretical background of the present invention will be described in detail.

According to the present invention, the travel-time table can be computed very fast with a small memory capacity and very low cost not only in the rectangular coordinates but also in the spherical coordinates using a straight ray- path. All of these two methods in the two coordinates have the same theoretical basis as the transformation of the velocity model. In other words, the direct computation of the travel-times in the rectangular coordinates also based on the same idea of the velocity model transform. We will discuss more detail about the velocity model transform. For the simplicity, we explain by taking the 2-D velocity model transform into the polar coordinates.

Suppose that we transform the velocity model in Cartesian coordinates into a polar grid having an origin at a source point. In this transformation, we invoke a bilinear or bi-cubic interpolation to assign velocities at polar grids. We then assume that the slowness is constant along a finite polar arc and only varies linearly in a radial direction (r) within a unit cone-like element constructed of adjacent four polar grid points as shown in FIG. 7. If we shoot ray in the radial direction, raypath becomes a straight line along the center axis of the cone. In this way we can trivially fix a take-off point and a terminal point of ray, therefore computing travel times from a source to a point in a simple manner. Furthermore, we can also calculate the amplitude by exploiting a transmission coefficient relationship between a layered cake model and a geometrical spreading. Our method gives no travel-time error in the ray tracing for a transformed velocity model but a misfit between the original model and the transformed model exists. FIGS. 8A-8D illustrate how fine polar cell division produces the model discretization error. FIG. 8A shows an original four-layered model, FIG. 8B shows the transformed velocity model, where Δr=70m and Δθ=10 degree, FIG. 8C shows the transformed velocity model, where Δr=7m and Δθ=l degree, and FIG. 8D shows the transformed velocity model, where Δr=0.7m and Δθ=0.1 degree. As we limit Δr→-0 and Δθ-»0, the transformed velocity model approaches to the original model on a rectangular grid.

To demonstrate the theoretical basis of the invention, we migrate the Marmousi data using the travel-times produced by the method of the invention. First, we compared the travel-times computed by the method of the invention with first arrival times by Vidale's algorithm and the maximum energy arrival times picked from the finite element wave equation modeling. FIG. 9 shows a snapshot at 1.1 second of the seismic wave field generated by the time domain finite element modeling. In FIG. 9, the travel-time computed by the method of the invention, Vidale's travel-time and maximum energy arrival time are superimposed on a snapshot. From FIG. 9, we note that the direct arrival times (black solid line) obtained by the method of the invention lies in the intermediate zone between the maximum energy travel-time denoted by a white solid line and the first arrival times denoted by a black dotted line. Therefore, if the direct arrival time obtained by the present invention is used, instead of the maximum energy travel-time, for imaging seismic wave, an image better than that is obtained by using the first arrival time.

FIG. 10A shows the final migrated image using the first arrival times by Vidale's algorithm, whereas FIG. 10B shows the final migrated image using the travel-times obtained by the present invention. From FIGS. 10A and 10B, we note that the travel-time of the invention provides better-focused and well-defined images. Especially, the present invention provides the most improved images in the anticline structures that become exploration targets due to their high possibility of containing petroleum.

While the present invention has been described with reference to the particular illustrative embodiments, it is not to be restricted by the embodiments but only by the appended claims. It is to be appreciated that those skilled in the art can change or modify the embodiments without departing from the scope and spirit of the present invention.

Claims

What Is Claimed Is:

1. A method of seismic imaging using direct travel time computing, comprising: a step of deterrnining a velocity model in rectangular coordinates (step

101); a step of deterrnining a source (step 102); a travel-time computing step of computing a direct travel-time when a wave is transmitted from the source to a grid in straight ray (step 103); a travel-time outputting step of outputting the computed travel-time (step

108); a travel-time output loop step of computing and outputting the travel-time for the next source (step 109); and a velocity model correction step of correcting the velocity model based on the computed travel-times for all sources and repeating the steps 101 to 109 (step 110).

2. The method according to claim 1, wherein the velocity model is a velocity model represented in rectangular coordinates constructed of cells having constant velocities, and the travel-time is obtained by dividing the raypath length in each cell by the velocity of each cell and adding up the divided results.

3. The method according to claim 1, wherein the velocity model is a velocity model represented in rectangular coordinates constructed of grids, and the travel-time is obtained by dividing the straight line connecting the source and the grid into sections, multiplying half the sum of reciprocals of the velocities at both end points of each section by the length of each section, and adding up the multiplied results for the divided sections.

4. The method according to claim 3, wherein the velocities at both end points of each section are obtained by the interpolation of the velocities of adjacent grids.

5. The method according to claim 1, wherein the travel-time is computed using polar coordinates or spherical coordinates in the travel-time computing step (step 103).

6. The method according to claim 5, wherein the step of computing the travel-time using the polar coordinates or spherical coordinates comprising: a step of transforming the velocity model in rectangular coordinates to the grids in the polar coordinates or spherical coordinates having an origin at a source point (step 105); a step of computing the travel-time by adding up time values each of which is calculated by multiplying the unit radial raypath length by half the sum of reciprocals of the velocity at each grid (step 106); and a step of transforming the travel-time computed by using the polar coordinates or spherical coordinates to the rectangular grids using an interpolation (step 107).

7. The method according to claim 1, wherein, in the outputting step (step 108), the travel-time is not outputted in the form of file stored in a hard disk or external storage device but used in a manner that depth migration program calls sub-program configured of the travel-time stored in a memory during the travel- time computing step.