1. Patent Search
  2. Details

SYSTEM AND METHOD FOR SEISMIC INVERSION

20180210101 ยท 2018-07-26

    Inventors

    Cpc classification

    International classification

    Abstract

    A method for elastic wave modeling up to tilted orthorhombic symmetry in anisotropy involves applying a curved grid adapting to a free-surface topography and transforming this the curved grid to a rectangular grid. Calculations using numerical discretization methods such as finite-difference are performed. The boundary condition for any free-surface is .sub.ij.Math.n.sub.j=0, which requires vanishing the global stresses' normal components to the free-surface. This resulting model more accurately simulates surface effects and better inversion for subsurface earth properties.

    Claims

    1. A method for simulating seismic waves, comprising: reading a plurality of parameters of a grid and a plurality of seismic acquisition data into a model; reading a plurality of parameter of a medium of a earth formation into the model; applying a free-surface boundary condition to the model, wherein the free-surface boundary conditions are applied in a curved grid; converting the free-surface boundary condition in the curved grid to a rectangular computational grid; constructing a Hooke's law equation (4) in the rectangular computational grid; constructing momentum conservation equations (15)-(17) in the rectangular computational grid; solving the Hooke's law equation and the momentum conservation equations in the rectangular computation grid by applying the free-surface boundary conditions; and obtaining a resulting data, wherein the output data are seismograms or snapshots of seismic pressure, or stress components, or particle velocities of the seismic wave.

    2. The method of claim 1, wherein the curved grid has a planar bottom and a curvilinear top that coincides with the free-surface topology.

    3. The method of claim 1, wherein the free-surface boundary condition is vanishing normal traction at any point on a surface of the curved grid.

    4. The method of claim 1, wherein the medium of the earth formation is an elastic medium of anisotropy of tilted orthorhombic symmetry.

    5. The method of claim 1, wherein a numerical discretization method is employed in solving the Hooke's law equation and the momentum conservation equations in the rectangular computation grid.

    6. The method of claim 5, where in the numerical discretization method is a finite-difference (FD) method, a finite element (FE) method, a spectral element (SE) method.

    7. The method of claim 1, further comprising inputting the resulting data into a seismic imaging method, a full waveform inversion (FWI) method, an elastic reverse time migration method.

    Description

    DESCRIPTIONS OF DRAWINGS

    [0009] FIG. 1 is a flow chart illustrating a method of the current disclosure.

    [0010] FIG. 2 illustrates in 2-D the 3-D transform from the curved grid (x, y, z) to the rectangular computational grid (, , ).

    [0011] FIGS. 3A-3F are seismic snapshots of the Vx (three images in the top row) and Vz (three images in the bottom row) tangential and normal component particle velocities for a Ricker explosion source released at the middle of a tilted, plane free-surface, covering a homogeneous, isotropic medium.

    [0012] FIGS. 4A-4F are pressure snapshots along a vertical plane of a Ricker point source released at the focus of a parabola.

    [0013] FIG. 5A and FIG. 5B are pressure snapshots using acoustic parameters in elastic tilted orthorhombic modeling, in which Ricker point source is released at depth of a homogeneous, acoustic tilted orthorhombic medium. For the purpose of comparison, FIG. 5C presents simulation results in a previously published article when a pure acoustic program code is employed.

    [0014] FIGS. 6A-6G compare Vx horizontal particle velocity snapshots along a vertical xz-plane of the SEAM isotropic topography model for with FIGS. 6H-6N, which are Vx horizontal particle velocity snapshots for a homogeneous model covered by the same topography. The same vertical Klauder wavelet is released at the surface of the model in each case.

    [0015] FIGS. 7A, 7B, and 7C are pressure seismograms generated using the heterogeneous SEAM isotropic topography model and respectively show the cross section of the middle, the left, and the right cross section of the same topography shown in FIGS. 6A-6N.

    [0016] FIGS. 8A-8D show representative corresponding seismograms from the simulation of FIGS. 6A-6N, with FIGS. 8A and 8C from the full (heterogeneous) SEAM topography model, and FIGS. 8B and 8D from the corresponding homogeneous model.

    DETAILED DESCRIPTION OF THE EMBODIMENT

    [0017] The present invention may be described and implemented in the general context of a system and computer methods to be executed by a computer. Such computer-executable instructions may include programs, routines, objects, components, data structures, and computer software technologies that can be used to perform particular tasks and process abstract data types. Software implementations of the present invention may be coded in different languages for application in a variety of computing platforms and environments. It will be appreciated that the scope and underlying principles of the present invention are not limited to any particular computer software technology.

    [0018] Moreover, those skilled in the art will appreciate that the present invention may be practiced using any one or combination of hardware and software configurations, including but not limited to a system having single and/or multiple computer processors, hand-held devices, programmable consumer electronics, mini-computers, mainframe computers, and the like. The invention may also be practiced in distributed computing environments where tasks are performed by servers or other processing devices that are linked through a one or more data communications network. In a distributed computing environment, program modules may be located in both local and remote computer storage media including memory storage devices.

    [0019] Also, an article of manufacture for use with a computer processor, such as a CD, pre-recorded disk or other equivalent devices, may include a computer program storage medium and program means recorded thereon for directing the computer processor to facilitate the implementation and practice of the present invention. Such devices and articles of manufacture also fall within the scope of the present disclosure.

    [0020] The Figures (FIG.) and the following description relate to the embodiments of the present disclosure by way of illustration only. It should be noted that from the following discussion, alternative embodiments of the structures and methods disclosed herein will be readily recognized as viable alternatives that may be employed without departing from the principles of the claimed inventions.

    [0021] Referring to the drawings, embodiments of the present disclosure will be described. Various embodiments can be implemented in numerous ways, including for example as a system (including a computer processing system), a method (including a computer implemented method), an apparatus, a non-transitory computer readable medium, a computer program product, a graphical user interface, a web portal, or a data structure tangibly fixed in a non-transitory computer readable memory. Several embodiments of the present invention are discussed below. The appended drawings illustrate only typical embodiments of the present disclosure and therefore are not to be considered limiting of its scope and breadth.

    [0022] Reference will now be made in detail to several embodiments of the present disclosure(s), examples of which are illustrated in the accompanying figures. It is noted that wherever practicable similar or like reference numbers may be used in the figures and may indicate similar or like functionality. The figures depict embodiments of the present disclosure for purposes of illustration only. One skilled in the art will readily recognize from the following description that alternative embodiments of the structures and methods illustrated herein may be employed without departing from the principles of the disclosure described herein.

    [0023] In a seismic acquisition, a shot is deployed at the source location, generating a wavefield that propagates from the source to the subsurface structure. The reflection from the subsurface structure propagates back to the surface and is detected by the receivers. The receivers convert the seismic signals to voltage signals. The voltage signals are then transmitted to a computer in a recording station (not shown) to be processed and converted into seismic data. The seismic data can be stored and transmitted for further processing.

    [0024] FIG. 1 shows a flow chart of the process of elastic seismic wave simulations with free-surface topography of the current disclosure. First, the parameters of the grid and the seismic acquisition data for seismic inversion are read and analyzed. Such parameters/data may include grid size (nx, ny, nz), directions (dx, dy, dz), cube dimensions (including the start and end values), shot location, run time, absorption layer thicknesses, absorption attenuation, receiver increments in x and y-directions, receiver depth, snapshot start, end and increment times, wanted output wave fields, wanted snapshot planes, etc.

    [0025] Next, the data regarding the rock physics and other input data are obtained. Such data can be one or more among P-velocity (Vp), S-velocity (Vs), formation density (), Thomsen epsilon parameters (.sub.1, .sub.2), Thomsen delta parameters (.sub.1, .sub.2, .sub.3), gamma ratios (.sub.1, .sub.2), porosity (), the incident angle (), and the opening angle between the symmetry axes in the rotated xy-plane ().

    [0026] After loading all necessary material parameters, a decision is made whether a free-surface or an absorbing surface should be simulated at the top of the model. If a free-surface is confirmed, then a choice is made whether a surface topography is chosen at the free-surface. If the free-surface is chosen as the surface topography, its data is read and modeled as the top boundary of the medium; if not, a free plane surface is applied. The standard output of any simulation is seismic pressure seismograms, with optional output of seismic pressure snapshots, or seismograms or snapshots of stress components. Alternatives not listed in the flow chart are snapshots and seismograms of particle velocities.

    [0027] FIG. 2 is a 2D illustration of the 3D transform from the curved grid (x, y, z) (presented in its xz-plane) to the rectangular grid (, , ) (presented in its -plane). The numerical discretizations are performed in the rectangular grid. The top of the curved grid coincides with a free-surface topography being modeled. The bottom of the curved grid is planar. The undulations in the vertical direction reduce with depth towards the bottom. The wave equations and boundary conditions are originally imposed in this physically conformed curved grid.

    [0028] When modeling a free-surface topography, it is assumed that elastic wave equations are given in a curved 3D grid. The mapping function from the curved (x, y, z) to the rectangular (, , ) system is assumed to have positive downward direction, and it can be written as (Puente et al., 2014)

    [00001] z ( , , ) = max .Math. ( z 0 ( , ) + m ) - z 0 ( , ) , ( 1 )

    where .sub.max is the maximum value (the depth) of the rectangular grid; z.sub.0 (,) is the elevation (topography) function, which is positive upwardsthe opposite of the vertical coordinates z and . m is the maximum depth of the physical model (the curved grid z(,,)). The chain rule is apply to express the spatial derivatives in the curved (x,y,z) grid by the ones in the rectangular (,,) grid, which gives

    [00002] x .fwdarw. .Math. x + .Math. x = + A ( , , ) .Math. ; .Math. .Math. A ( , , ) = x = max - z 0 ( , ) + m .Math. h x , .Math. y .fwdarw. .Math. y + .Math. y = + B ( , , ) .Math. ; .Math. .Math. B ( , , ) = y = max - z 0 ( , ) + m .Math. h y , .Math. z .fwdarw. .Math. z = C ( , ) .Math. ; .Math. C ( , ) = z = max z 0 ( , ) + m , ( 2 )

    where

    [00003] h x = z 0 ( , )

    is the topograpny slope in the x (or ) direction and

    [00004] h y = z 0 ( , )

    is the topography slope in the y (or ) direction. Details of derivation of the coordinate partial derivatives A(,,),B(,,) and C(,) can be found in Hestholm and Ruud (1998).

    [0029] Let and with components .sub.ij and .sub.i, i,j=1 . . . 3, be stresses and particle velocities in a vertical orthorhombic medium. Let and , with components .sub.ij and .sub.i, be the stresses and particle velocities directed along the coordinates of the symmetry axes of a tilted orthorhombic medium. The rotation of velocities from vertical to tilted orthorhombic system can then be written =R, and the rotation of stresses from vertical to tilted orthorhombic system can be written =M, where M is the Bond matrix, consisting of products and sums of the components R.sub.ij of rotation matrix R. The rotation matrix is given by

    [00005] R = ( cos .Math. .Math. .Math. .Math. cos .Math. .Math. .Math. .Math. cos .Math. .Math. - sin .Math. .Math. .Math. .Math. sin .Math. .Math. - cos .Math. .Math. .Math. .Math. cos .Math. .Math. .Math. .Math. sin .Math. .Math. - sin .Math. .Math. .Math. .Math. cos .Math. .Math. cos .Math. .Math. .Math. .Math. sin .Math. .Math. sin .Math. .Math. .Math. .Math. cos .Math. .Math. .Math. .Math. cos .Math. .Math. + cos .Math. .Math. .Math. .Math. sin .Math. .Math. - sin .Math. .Math. .Math. .Math. cos .Math. .Math. .Math. .Math. sin .Math. .Math. + cos .Math. .Math. .Math. .Math. cos .Math. .Math. sin .Math. .Math. .Math. .Math. sin .Math. .Math. - sin .Math. .Math. .Math. .Math. cos .Math. .Math. sin .Math. .Math. .Math. .Math. sin .Math. .Math. cos .Math. .Math. ) , ( 3 )

    where , , and , respectively, are the angles of azimuth, tilt, and opening angle between the symmetry axes in the rotated xy-plane.

    [0030] In an embodiment of the current disclosure, when simulating a tilted orthorhombic system, the stress-strain relationship (Hooke's law) is given by


    =C,(4)

    with C being the full 66 stiffness matrix in the rotated system, which is given in terms of the Bond matrix M as


    C=MCM.sup.T.(5)

    [0031] and are stress and strain in the tilted orthorhombic system, and C is the stiffness matrix in the vertical orthorhombic system. C, being a full matrix, increases the amount of computations and/or storage from those of a vertical orthorhombic system. Furthermore, when the curved grid is introduced and transformed to the rectangular grid according to transform (2), the stress-strain relation described by equation (4) will have time derivatives of strain components .sub.i, i=1 . . . 6, given by

    [00006] .Math. 1 t = v x + A .Math. v x , ( 6 ) .Math. 2 t = v y + B .Math. v y , ( 7 ) .Math. 3 t = C .Math. v z , ( 8 ) .Math. 4 t = C .Math. v y + v z + B .Math. v z , ( 9 ) .Math. 5 t = C .Math. v x + v z + A .Math. v z , ( 10 ) .Math. 6 t = v x + B .Math. v x + v y + A .Math. v y . ( 11 )

    [0032] Similarly, Newton's 2.sup.nd law (the momentum conservation equations) in a curved coordinate system,

    [00007] v x t = 1 .Math. ( xx x + xy y + xz z ) , ( 12 ) v y t = 1 .Math. ( xy x + yy y + yz z ) , ( 13 ) v z t = 1 .Math. ( xz x + yz y + zz z ) , ( 14 )

    when transformed to the rectangular (,,)-system through coordinate transform (2) becomes

    [00008] v x t = 1 .Math. ( xx + A .Math. xx + xy + B .Math. xy + C .Math. xz ) , ( 15 ) v y t = 1 .Math. ( xy + A .Math. xy + yy + B .Math. yy + C .Math. yy ) , ( 16 ) v z t = 1 .Math. ( xz + A .Math. xz + yz + B .Math. yz + C .Math. zz ) . ( 17 )

    [0033] At the free-surface, the rectangular grid vertical coordinate =0, so the equations (18) can be derived from equations (2)

    [00009] A ( , , ) = max z 0 ( , ) + m .Math. h x = C ( , ) .Math. h x , .Math. B ( , , ) = max z 0 ( , ) + m .Math. h y = C ( , ) .Math. h y , .Math. C ( , ) = max z 0 ( , ) + m . ( 18 )

    [0034] In the embodiment of the current disclosure, the boundary condition for any free-surface is vanishing normal traction, .sub.ij.Math.n.sub.j=0, where .sub.ij are the stress tensor components, and n.sub.j are the components of the inward directed (unit) normal vector n=(h.sub.x, h.sub.y, 1) to the surface, with h.sub.x and h.sub.y the topography slopes as before. Component-wise this becomes, after differentiating with respect to time,

    [00010] h x .Math. xx t + h y .Math. xy t + xz t = 0 , .Math. h x .Math. xy t + h y .Math. yy t + yz t = 0 , .Math. h x .Math. xz t + h y .Math. yz t + zz t = 0. ( 19 )

    [0035] Hooke's law, equation (4), is differentiated with respect to time, and the time differentiated stresses from it are inserted into boundary conditions (19). Moving all vertical derivatives of the particle velocities to the left hand sides of the equations, leads to the following boundary conditions in terms of the velocities .sub.i in the tilted orthorhombic system,

    [00011] v x [ h x 2 .Math. Cc 11 + h y 2 .Math. Cc 66 + h x .Math. h y .Math. C ( c 16 + c 61 ) + h x .Math. C ( c 15 + c 51 ) + h y .Math. C ( c 56 + c 65 ) + Cc 55 ] + v y [ h x 2 .Math. Cc 16 + h y 2 .Math. Cc 62 + h x .Math. h y .Math. C ( c 12 + c 66 ) + h x .Math. C ( c 14 + c 56 ) + h y .Math. C ( c 64 + c 52 ) + Cc 54 ] + v z [ h x 2 .Math. Cc 15 + h y 2 .Math. Cc 64 + h x .Math. h y .Math. C ( c 14 + c 65 ) + h x .Math. C ( c 13 + c 55 ) + h y .Math. C ( c 63 + c 54 ) + Cc 53 ] = - h x ( c 11 .Math. v x + c 12 .Math. v y + c 14 .Math. v z + c 15 .Math. v z + c 16 ( v x + v y ) ) - h y ( c 61 .Math. v x + c 62 .Math. v y + c 64 .Math. v z + c 65 .Math. v z + c 66 ( v x + v y ) ) - c 51 .Math. v x - c 52 .Math. v y - c 54 .Math. v z - c 55 .Math. v z - c 56 ( v x + v y ) , ( 20 ) v x [ h x 2 .Math. Cc 61 + h y 2 .Math. Cc 26 + h x .Math. h y .Math. C ( c 21 + c 66 ) + h x .Math. C ( c 41 + c 65 ) + h y .Math. C ( c 25 + c 46 ) + Cc 45 ] + v y [ h x 2 .Math. Cc 66 + h y 2 .Math. Cc 22 + h x .Math. h y .Math. C ( c 26 + c 62 ) + h x .Math. C ( c 46 + c 64 ) + h y .Math. C ( c 24 + c 42 ) + Cc 44 ] + v z [ h x 2 .Math. Cc 65 + h y 2 .Math. Cc 24 + h x .Math. h y .Math. C ( c 25 + c 64 ) + h x .Math. C ( c 45 + c 63 ) + h y .Math. C ( c 23 + c 44 ) + Cc 43 ] = - h x ( c 61 .Math. v x + c 62 .Math. v y + c 64 .Math. v z + c 65 .Math. v z + c 66 ( v x + v y ) ) - h y ( c 21 .Math. v x + c 22 .Math. v y + c 24 .Math. v z + c 25 .Math. v z + c 26 ( v x + v y ) ) - c 41 .Math. v x - c 42 .Math. v y - c 44 .Math. v z - c 45 .Math. v z - c 46 ( v x + v y ) , .Math. .Math. and ( 21 ) v x [ h x 2 .Math. Cc 51 + h y 2 .Math. Cc 46 + h x .Math. h y .Math. C ( c 41 + c 56 ) + h x .Math. C ( c 31 + c 55 ) + h y .Math. C ( c 36 + c 45 ) + Cc 35 ] + v y [ h x 2 .Math. Cc 56 + h y 2 .Math. Cc 42 + h x .Math. h y .Math. C ( c 52 + c 46 ) + h x .Math. C ( c 54 + c 36 ) + h y .Math. C ( c 32 + c 44 ) + Cc 34 ] + v z [ h x 2 .Math. Cc 55 + h y 2 .Math. Cc 44 + h x .Math. h y .Math. C ( c 45 + c 54 ) + h x .Math. C ( c 35 + c 53 ) + h y .Math. C ( c 34 + c 43 ) + Cc 33 ] = - h x ( c 51 .Math. v x + c 52 .Math. v y + c 54 .Math. v z + c 55 .Math. v z + c 56 ( v x + v y ) ) - h y ( c 41 .Math. v x + c 42 .Math. v y + c 44 .Math. v z + c 45 .Math. v z + c 46 ( v x + v y ) ) - c 31 .Math. v x - c 32 .Math. v y - c 34 .Math. v z - c 35 .Math. v z - c 36 ( v x + v y ) , ( 22 )

    [0036] where the A's and B's have been replaced by Ch.sub.x and Ch.sub.y respectively, from relations (18) valid at the free-surface. Equations (20)-(22) are boundary conditions in terms of particle velocities for free-surface topography on top of tilted orthorhombic or simpler anisotropic elastic media.

    [0037] For a plane, horizontal free-surface, the curved grid reduces to a rectangular grid, so we can set C1, and the slopes h.sub.x and h.sub.y vanish. Then the boundary conditions reduce to

    [00012] v x .Math. c 55 + v y .Math. c 54 + v z .Math. c 53 = - c 51 .Math. v x - c 52 .Math. v y - c 54 .Math. v z - c 55 .Math. v z - c 56 ( v x + v y ) , ( 23 ) v x .Math. c 45 + v y .Math. c 44 + v z .Math. c 43 = - c 41 .Math. v x - c 42 .Math. v y - c 44 .Math. v z - c 45 .Math. v z - c 46 ( v x + v y ) , .Math. .Math. and ( 24 ) v x .Math. c 35 + v y .Math. c 34 + v z .Math. c 33 = - c 31 .Math. v x - c 32 .Math. v y - c 34 .Math. v z - c 35 .Math. v z - c 36 ( v x + v y ) . ( 25 )

    [0038] A regular standard grid or a staggered grid may also be used. However, regular standard grids are preferable to staggered grids, since often regular staggered grids do not naturally adapt to correct physical locations when describing anisotropic wave equations and topography boundary conditions. Standard grid discretization tends to be less stable in application. Therefore, Fully Staggered Grids (FSG; Puente et al., 2014) is preferable to ensure accuracy and stability in complex media in spite of increased memory and computational requirement associated with FSG.

    [0039] The time derivatives of equations (4) were solved for the stresses and equations (15)-(17) for the velocities. These nine field variables become 36 (a factor 4) due to the use of FSG. The same factor theoretically applies to the increase of computational cost, however due to loop parallelization (using open MP) this factor is closer to 3 in practice. The wave equations (4) and (15)-(17) are solved as first-order partial differential equations (PDEs) for propagation in an elastic tilted orthorhombic medium. They are discretized by eight-order staggered FDs (Fornberg, 1988) in space and second-order staggered FDs in time. Nine full grids of material parameters (the nine stiffness parameters) are required for elastic vertical orthorhombic modeling, and the 21 stiffnesses required for tilted orthorhombic modeling can be calculated on the fly from these to save memory. Then in addition the three rotation angles for dip, azimuth and need to be stored in full 3D grids.

    [0040] In an embodiment of the current disclosure, the method involves storing all 21 stiffnesses but reducing the number of other computations. On the basis of definitions given in Tsvankin (1987), relationships can be calculated for transition from material input parameters 's, 's and 's mentioned in FIG. 1, to the material stiffness matrix Cwhich again determines the (an)isotropic speeds Vp and Vs and density .

    [0041] Boundary conditions (20)-(22) (alternatively in their form of a free plane-surface, equations (23)-(25)), are implemented at the free-surface. Full (8th) FD-order can be maintained in the implementation of these free-surface boundary conditions by assuming the local particle velocities to be constant above the free-surface in a layer of nodal depth of the FD-operator's half order (4 in examples in this disclosure; Jastram and Behle, 1993). Experience shows no qualitative degeneration of results by assuming such constant fields above the free-surface, when compared to the alternative, which is gradual reduction of FD-order from interior medium (8.sup.th order FDs) towards minimum 2.sup.nd order FDs at the free-surface. In this process, one would go via 6.sup.th and 4.sup.th order central FDs. For the stresses, mirror conditions around the free-surface are imposed on the normal traction, implying that its surface value is zero.

    [0042] For absorbing boundary conditions along all boundaries (except the top, in free-surface cases), the sponge/exponential damping method (Cerjan et al., 1985) was used, for which the damping coefficient for optimal absorption has to be found empirically for the relevant wave propagation scheme.

    [0043] FIGS. 3A-3F show a same homogeneous elastic isotropic medium covered by a tilted, plane free-surface. The tilt in this case is negative 30. Vp=3500 m/s. Vs=1750 m/s. The density is 2084.5 kg/m.sup.3. A Ricker explosion source of peak frequency 15 Hz is released at the middle of the tilted free-surface. Dimensions in x- and y-directions are 4000 m and model depth (in rectangular grid) is 6000 m. The grid lengths are 10 m uniformly. Twenty (20) absorbing layers are added to each grid edge except at the free-surface. The top three panels (FIGS. 3A-3C) are snapshots in the xz-plane whiles the bottom three panels (FIGS. 3D-3F) are snapshots in the yz-plane. Components shown are tangential (top) and normal (bottom) particle velocities to the free-surface. The P-, S-, and head waves connecting them from the surface, as well as the Rg-waves at the surface are observed and are labeled in the two snapshots on the right. Note that the Rg-waves are slightly slower than the S-wave front.

    [0044] FIGS. 4A-4F are pressure snapshots along a vertical plane of a Ricker point source released at the focus of a parabolic free-surface, and the resulting plane wave reflected back from it into the medium. The model has homogeneous P- and S-velocities of 4500 and 2600 m/s, and the grid distance is 20 m in the horizontal directions and 40 m in the vertical direction. The result of the simulation is consistent with a known effect in the particular case of a free-surface parabola covering a homogeneous, isotropic medium. That is, the reflection of a point source released at the focus of the parabola from the free-surface back into the medium is expected to be a plane wave, which is approximately true for the reflected pressure (P-) wave.

    [0045] In FIGS. 4A-4F, the parabola is a free-surface topography over a homogeneous, isotropic medium. The pressure (P-) wave reflected from the free-surface back into the medium is approximately a plane wave, as expected. The deviation from the circle around the P-wave in the three panels (FIGS. 4A-4C) in the top row are numerical grid artifacts from the curved grid. The artifacts in three panels (FIGS. 4D-4F) in the bottom row along the parabola free-surface after the plane P-wave reflection back from it are the physical phenomena of the slower S/Rg (Rayleigh) wave train propagating along the surface following the horizontal plane P-wave. There is a wave proceeding the plane P-wave reflection, appear like a hook on either side of the plane P-wave are due to grid artifacts from the compression of the grid nodes at either side (the steepest part) of the parabola.

    [0046] FIGS. 5A and 5B show a homogeneous acoustic tilted orthorhombic model covered by a plane horizontal free-surface, using acoustic parameters in the elastic code. The model has Vp=4500 m/s, Vs=2250 m/s, density=2200 kg/m.sup.3, .sub.1=0.2, .sub.2=0.12, .sub.1=.sub.2=0.06, .sub.3=0, .sub.1=.sub.2=0, =35, =40 and =25. The dimensions are 4000 m horizontally in both directions and 6000 in depth. Grid distances are 10 m uniformly, and a Ricker point source of peak frequency 15 Hz is excited at 2500 m depth. For comparison, the bottom panel FIG. 5C shows P-waves from explosion source in acoustic simulation from Chu (2012). As one can readily appreciate, the P-wave form resulting from a simulation according to this disclosure is consistent with those in Chu, further validating the applicability of this modeling method.

    [0047] FIGS. 6A-6G are Vx particle velocity component snapshots along a vertical xz-plane through the SEAM (SEG Advance Modeling Program) heterogeneous, isotropic topography model while FIGS. 6H-6N from the same model but with homogeneous interior parameters Vp=3500 m/s, Vs=1750 m/s and density 2200 kg/m.sup.3. A vertical Klauder wavelet of central frequency 15 Hz is released at the same location of the free-surface of both models. The two cases are compared for the same free-surface topography. Common wave features can be confirmed. However, snapshots for the heterogeneous medium (panels on the left, FIGS. 6A-6G) exhibit more wavefield complexities, e.g., scattering from interior inhomogeneities.

    [0048] FIGS. 7A, 7B, and 7C respectively shows pressure seismograms from middle, left and right cross sections from the above SEAM topography (heterogeneous) example of FIGS. 6A-6N. P-wave, S-wave, as well as scattering are labeled.

    [0049] FIGS. 8A-8D are some representative corresponding seismic sections from the full heterogeneous isotropic SEAM topography model (FIGS. 8A and 8C) and the homogeneous model (FIGS. 8B and 8D) of the simulations of FIGS. 6A-6N. Both exhibit similar waveforms. However, the heterogeneous isotropic SEAM topography model provides far more details in the heterogeneities add far more wavefield complexity.

    [0050] While in the foregoing specification this invention has been described in relation to certain preferred embodiments thereof, and many details have been set forth for purpose of illustration, it will be apparent to those skilled in the art that the invention is susceptible to alteration and that certain other details described herein can vary considerably without departing from the basic principles of the invention. In addition, it should be appreciated that structural features or method steps shown or described in any one embodiment herein can be used in other embodiments as well.

    REFERENCES

    [0051] Cerjan, C., Kosloff, D., Kosloff R., and Reshef, M. (1985) Short note on a nonreflecting boundary condition for discrete acoustic-wave and elastic-wave equations, Geophysics, 50, pp. 705-708. [0052] Fornberg, B. (1988) Generation of Finite Difference Formulas on Arbitrarily Spaced Grids, Mathematics of Computation, 51, no. 184, pp. 699-706. [0053] Hestholm, S., and Ruud, B. (1998) 3-D finite-difference elastic wave modeling including surface topography, Geophysics, 63 no. 2, pp. 613-622. [0054] S. Hestholm (2002) Composite Memory Variable Velocity-Stress Viscoelastic Modelling, Geophys. J. Int., 148, pp. 153-162. [0055] S. Hestholm and B. Ruud (2002) 3-D Free-Boundary Conditions for Coordinate-Transform Finite-Difference Seismic Modelling, Geophys. Prosp., 50, pp. 463-474. [0056] Jastram, C., and Behle, A. (1993) Elastic Modeling by Finite-Difference and the Rapid Expansion Method (REM), Expanded abstract, ST3.6, Soc. Of Exploration Geophysics, pp. 1573-1576. [0057] de la Puente, J., Ferrer, M., Hanzich, M., Castillo, J. E., and Cela, J. M. (2014) Mimetic seismic wave modeling including topography on deformed staggered grids, Geophysics, 79, no. 3, pp. T125-T141. [0058] Tessmer, E., and Kosloff, D. (1994) 3-D elastic modeling with surface topography by a Chebychev spectral method, Geophysics, 59, no. 3, pp. 464-473. [0059] Tsvankin, I. (1997) Anisotropic parameters and P-wave velocity for orthorhombic media, Geophysics, 62, no. 4, pp. 1292-1309.