In this study, a new finite-difference technique is
designed to reduce the number of grid points needed in
frequency-space domain modeling. The new algorithm
uses optimal nine-point operators for the approximation
of the Laplacian and the mass acceleration terms. The
coefficients can be found by using the steepest descent
method so that the best normalized phase curves can be
obtained.
ABSTRACT
This method reduces the number of grid points per
wavelength to 4 or less, with consequent reductions of
computer memory and CPU time that are factors of tens
less than those involved in the conventional secondorder
approximation formula when a band type solver is
used on a scalar machine

Jo, Churl-Hyun

Shin, Changsoo

Suh, Jung Hee

Geophysics

SNU Open Repository and Archive

GEOPHYSICS, VOL. 61, NO. 2 (MARCH-APRIL 1996); P. 529-537, 11 FIGS.An optimal 9-point, finite-difference, frequency-space,2-D scalar wave extrapolatorChurl-Hyun Jo*, Changsoo Shin*, and Jung Hee Suh‡In this study, anew finite-difference technique isdesigned to reducethe number of grid points needed infrequency-spacedomain modeling. The new algorithmuses optimal nine-point operators for the approximationof the Laplacian and the mass acceleration terms. Thecoefficients can be found by using the steepest descentmethod so that the best normalized phase curves can beobtained.ABSTRACTThis method reduces the number of grid points perwavelength to 4 or less, with consequent reductions ofcomputer memory and CPU time that are factors of tensless than those involved in the conventional second-order approximation formula when a band type solver isused on a scalar machine.MOTIVATION  Numerical modeling of the wave equation in the frequency-space domain was pioneered by Lysmer and Drake (1972) andhas been developed by Marfurt (1984) Shin (1988), Marfurtand Shin (1989) Pratt and Worthington (1990) and Pratt(1990a, b). Finite-element modeling in the frequency domainhas been extensively upgraded by Marfurt and his colleagues(Marfurt 1988, personal communication).Wave equation modeling in the time domain is popularbecause of its easy implementation and accuracy, compared tofrequency domain modeling. The advantage of frequency-domain modeling over time-domain modeling is that multi-experiment seismic data can be simulated economically bydirect multiplication once the triangular factors of theimpedance matrix are calculated. Modeling the effects ofattenuation is more flexible in the frequency domain than inthe time domain (Pratt, 1990b), because in the frequencydomain one can directly input the attenuation coefficient asa function of frequency. Furthermore, for certain geome-tries, only a few frequency components are required toperform wave-equation inversion and tomography (Prattand Worthington, 1990).Because of numerical dispersion and anisotropy, the numer-ical solution of the wave equation can lead to inaccurateresults. Such errors can be overcome by generating a suffi-ciently fine grid compared to the wavelength. In time-domainmodeling, higher-order differencing operators are used toreduce the required spatial sampling (Dablain, 1986). For theconventional second-order central finite-difference operatorused by Pratt and Worthington (1990), more grid points perwavelength are needed in frequency-domain solutions to ob-tain an accuracy that is comparable to that achieved withtime-domain solutions. Moreover, the fact that frequency-domain modeling can be done only implicitly makes theproblem harder because of the huge resultant matrix.In this paper, several techniques are employed to increasethe accuracy of frequency-domain modeling; these being (1)exploiting the extended Laplacian difference operator to re-duce numerical anisotropy, (2) combining the mass accelera-tion term into the point collocation and weighted average term(Marfurt, 1984) and (3) finding averaging coefficients usingthe steepest descent method so that the best normalized phasecurves can be obtained.ACCURACY OF NUMERICAL MODELINGThe scalar wave equation for a homogeneous isotropicmedium in the frequency domain is     ( 1 )where P = the pressure wavefield, = the angular frequency,and v = velocity.When the conventional second-order, central finite-differ-ence, five-point approximation is applied to the Laplacian termin two dimensions, equation (1) becomesManuscript received by the Editor May 5, 1995; revised manuscript received August 9, 1995.*Korea Institute of Geology, Mining, and Materials, P.O. Box 111, Taedok Science Town, Korea 305-350.‡Department of Petroleum and Mineral Engineering, Seoul National University, Seoul, Korea, 151-742.© 1996 Society of Exploration Geophysicists. All rights reserved.529530 Jo et al.          (2)where Pm,n represents the pressure of the wavefield at thelocation             Figure 1shows this grid.Since exploration seismology requires relating events in timeto horizons in depth, one concern is to minimize numericalvelocity or dispersion errors (Alford et al., 1974; Marfurt,1984). Dispersion analysis (Appendix A) allows an estimationof the accuracy of the solution in the frequency domain.Figures 2 and 3 compare dispersion curves corresponding tofrequency-domain and time-domain algorithms. Since disper-sion in the frequency domain is greater than that in the timedomain, more grid points per wavelength should be used in thefrequency domain to obtain an accuracy that is comparable tothat of explicit time-domain extrapolation.A NEW FINITE-DIFFERENCE SCHEMETo get more accurate results with the same, or a smallernumber of, grid points per wavelength, a different implemen-tation is suggested.The first step is to generalize the Laplacian term      (3)where         and the (x', z') coordinate system is rotated 45” relative to the(x, z) system. Finite-difference approximations to these oper-ators can be expressed as           and P +    +  + where A = AX = AZ.This generalization makes a nine-point Laplacian operator(Figure 4). With this formulation, the dispersion curves arecalculated with a = 0.5. The result (Figure 5) shows relativelysmall numerical anisotropy with respect to propagation angleeven though the method is, in general, less accurate than theconventional technique shown in Figure 2. The dispersion ismaximum at a propagation angle of 45” and minimum at 0”.This difference can be explained by the fact that the gridinterval of the 45” rotated coordinate system is larger than thatof the 0” system.The second step is to modify the technique used by Marfurt(1984) in which he considered the mass acceleration term to bea linear combination of the lumped mass matrix and theconsistent mass matrix. Using this technique, the finite-differ-ence approximation of the collocation point P of the massFIG. 1. The pressure fields P at the collocation point (m, n) andits eight nearest neighbors in a 2-D medium.FIG. 2. Dispersion curves for finite-difference solution of thescalar 2-D wave equation in the frequency-space domain.Numerical phase velocity Vph and numerical group velocity Vgrare normalized with respect to the true velocity V0 and plottedversus wavenumber        where G is thenumber of grid points per wavelength.f-x Domain Wave Extrapolator 531FIG. 3. Normalized phase Vph and group Vgr velocity curves for different propagation angles with respect to the grid, for the explicitsecond-order central finite-difference scheme in time-space domain when the stability limit,      where c is thevelocity.FIG. 4. Finite-difference stars for the Laplacian operator. (a) Conventional second-order central difference(five-point) star, (b) 45” rotated star, (c) nine-point star combining (a) and (b).532 Jo et al.acceleration term   in equation (1) can be representedas a linear combination of the points corresponding to theLaplacian operators, so that            In practice the averaging coefficient b can be found by optimi-zation (See Appendix B) in a manner similar to that used inHolberg (1987) to obtain optimal differencing operators intime domain modeling. Dispersion curves based on the opti-mized constant, b = 0.737, are shown in Figure 6. In this figure,dispersion curves are clustered on either side of the linesVp h/ v= 1 and Vgr/v = 1, which means that the numericalvelocity does not diverge significantly from the true mediumvelocity.Finally, by combining these two techniques we can expectmore accurate results. The quantity P in the mass accelerationterm  P, in conjunction with the nine-point Laplacianoperator, is                 (5)where c + 4d + 4e = 1.Substitution of equations (5) and (3) into equation (1) givesthe following difference equation                                                (6)The coefficients are optimized by minimizing the numericaldispersion error of the phase velocity (See Appendix B).FIG. 5. Normalized phase Vph and group Vgr velocity curves forfinite-difference solutions in the frequency domain using thenine-point finite-difference formulation to approximate theLaplacian operator when a = 0.5 in equation (3).FIG. 6. Normalized phase Vph and group Vgr velocity curves forfrequency-domain, finite-difference solutions using an averageterm with b = 0.7370 in equation (4).f-x Domain Wave Extrapolator 533FIG. 7. Normalized phase Vph and group Vgr velocity curves for frequency-domain, finite-difference solutionsusing the nine-point finite-difference formulation with a = 0.5461, c = 0.6248, and d = 0.9381 X 10-l  inequation (6).FIG. 8. A homogeneous half-space model to test the frequency-domain modeling. The velocity of the medium is3000 m/s. The symbol * denotes the shot point; every tenth receiver is shown with the symbol  534 Jo et al.Numerical integration is done to calculate the residual error Our new scheme uses an optimal nine-point, finite-differ-[equation (B-l)], using an interval of integration from 0” to 45” ence approximation of the Laplacian and mass accelerationfor angle  and 1/512 to 1/4 for wavenumber k*. The optimized terms. We admit that the frequency-domain modeling tech-coefficients are a = 0.5461, c = 0.6248, and d = 0.9381 X 10-l. nique cannot compete with other modeling techniques, such asDispersion curves based on these optimized coefficients are the Fourier method (Kosloff and Baysal, 1982) in modelingshown in Figure 7. When Gmin, the number of grid points per conventional synthetic sections. However, for certain problemsshortest wavelength, is 4, the phase velocity error will be (e.g., modeling structures with frequency-dependent physicalbounded within ±0.5%, and when Gmin is 3.2, the phase properties and wave-equation tomography where only a smallvelocity error is bounded within ±1%. For group velocity, this number of frequency components are used), this frequency-value will produce about 3% error. For a comparable degree of domain method will make the problem less expensive. Further-accuracy, the conventional five-point frequency-domain algo- more, the optimal finite-difference formulation, with a morerithm would have Gmin = 13, while the explicit time-domain, efficient sparse matrix solver such as nested dissection, canfinite-difference method would have Gmin = 10. further reduce the computational cost of such modeling.We use a single active column matrix solver (Zienkiewicz,1978) to factorize and solve the impedance matrix. Fastersolutions would result if more efficient matrix solution meth-ods, such as nested dissection (George and Liu, 1981) areused.ACKNOWLEDGMENTSTime-domain seismograms are generated by Fourier synthe-sis; complex frequencies can be used to suppress the wrap-around effect of the Fast Fourier transform (FFT).EXAMPLESTo examine the fidelity of solutions generated with theoptimal nine-point, finite-difference formulation developed inthis study, the conventional frequency-domain solution of Prattand Worthington (1990) and the time-domain solution ofAlford et al. (1974) are compared for the same homogeneousisotropic half-space model shown in Figure 8. The optimalaveraging coefficients used are the same as those used inF i g u r e  7 ;  G m i n= 4 in the frequency-domain calculations,Gmin = 9 in the time domain, and the high-cut frequency of thesource wavelet is 60 Hz. For the time-domain solution, second-order central finite differences are used to approximate thespatial derivatives. Time integration is performed explicitly andrun at the stability      The seismogramgenerated by the conventional five-point frequency-domainscheme (Figure 9a) shows evidence of dispersion because ofinsufficient nodal points, while the seismograms calculated bythe optimal nine-point operator (Figure 9b), as well as those byexplicit time-domain modeling (Figure 9c), show no apparentdispersion. Another example for a syncline model (Figure 10)illustrates the reflection characteristics of the optimal nine-point, finite-difference formulation. Results of time-domainmodel ing are presented for comparison. The resul ts(Figure 11) look similar to each other, except for unstablesignals at late times for near-shot receivers in the frequency-domain modeling (Figure 1 la), which are caused by timealiasing of energy before time zero.The research was inspired by Dr. K. R. Kelly of Amoco. Wewish to thank Prof. B. K. Hyun and C. E. Baag of SeoulNational University and Dr. S. H. Cheong of KIGAM forhelpful discussions during the study. Dr. K. J. Marfurt ofAmoco was kind enough to review the manuscript and to offermany suggestions that greatly improved the study. We are alsograteful to Dr. R. G. Pratt of Imperial College and Prof. E. A.Robinson for their constructive comments on the manuscript.REFERENCESAlford, R. M., Kelly, K. R., and Boore, D. M., 1974, Accuracy offinite-difference modeling of the acoustic wave equation: Geophys-ics, 39, 834-842.Cont, S. D., and de Boor, C., 1980, Elementary numerical analysis, 3rded.: McGraw-Hill Book Co.Dablain, M. A., 1986, The application of high-order differencing to thescalar equation: Geophysics, 51, 54 - 66.George, A., and Liu, J. W., 1981, Computer solution of large sparsepositive definite systems: Prentice-Hall, Inc.Holberg, O., 1987, Computational aspect of the choice of operator andsampling interval for numerical differentiation in large-scale simu-lation of wave phenomena: Geophys. Prosp., 35, 629-655.Kosloff, D. D., and Baysal, E., 1982, Forward modeling by a Fouriermethod: Geophysics, 47, 1402-1412.Lysmer and Drake, 1972, A finite-element method for seismology, inBolt, B. A., Eds., Methods in computational physics, Volume 11:Seismology: Surface waves and earth oscillations: Academic PressInc., 181-216.Marfurt, K. J., 1984, Accuracy of finite-difference and finite-elementmodeling of the scalar and elastic wave equations: Geophysics, 49,533-549.Marfurt, K. J., and Shin, C. S., 1989, The future of iterative modelingin geophysical exploration, in Eisner, E., Ed., Handbook of geo-physical exploration: I Seismic Exploration, Volume 21: Supercom-puters in seismic exploration: Pergamon Press, 203-228.Pratt, R. G., 1990a, Inverse theory applied to multi-source cross-holetomography, Part II: Elastic wave-equation method: Geophys.Prosp., 38, 287-310.CONCLUSION— 1990b, Frequency-domain elastic wave modeling by finitedifferences: A tool for crosshole seismic imaging: Geophysics, 55,626-632.Numerical modeling in the frequency domain has not beenwidely used since the large matrix equations require largecomputer memories. These matrices can be reduced in size byusing a more accurate scheme.Pratt, R. G., and Worthington, M. H., 1990, Inverse theory applied tomulti-source cross-hole tomography, Part I: Acoustic wave-equationmethod: Geophys. Prosp., 38, 287-310.Shin, C. S., 1988, Nonlinear elastic wave inversion by blocky parame-terization: Ph.D. thesis, Univ. of Tulsa.Zienkiewicz, 0. C., 1978, The finite-element methods, 3rd ed.:McGraw-Hill Book Co.f-x Domain Wave Extrapolator 5 3 5FIG. 9. Synthetic seismograms generated by the finite-difference solutions of the 2-D scalar waveequation for thehomogeneoushalf-space model shown in Figure 8 by (a) the conventional frequency-domain method, (b) the frequency-aomam formula using theoptimal nine-point finite differences, and (c) the explicit time-domain scheme (second-order central difference).536 Jo et al.FIG. 10. Syncline model showing receiver and source geometry. V1 = 3000 m/s, V2 = 5000 m/s.FIG. 11. Synthetic seismograms generated from thennit -difference solutions for the syncline model shown in Figure 10 by (a) thefrequency-domain formula using the optimal nine-point finite differences, and (b) the explicit time-domain scheme. The shot pointis located at 750 m.537f-x Domain Wave ExtrapolatorAPPENDIX ADISPERSION ANALYSISSubstituting a plane harmonic wave solution of the form    into equation (6) gives the dispersion equation                                              V (A-1)where      is the grid interval, c + 4d + 4e = 1, and  is the propagation angle from the z-axis. From the definition ofphase and group velocities,       the normalized numerical phase and group velocities can be written asandwhere    (A-2b)                (A-2a)v is the medium velocity,   the number of grid points per wavelength,               and+ sin  cos          + sin  cosAPPENDIX BOPTIMIZATIONSince numerical dispersion in 2-D depends on propagation To obtain the optimal values of a, c, and d, the method ofangle as well as on grid size, the averaging coefficients a, c, and steepest descent had been applied to the minimization of E. Byd in equations (6) and (A-2a) should be determined such that a general iterative rule (Cont and De Boor, 1980), the param-the L2 norm of the residual error is minimized. When the error eters a, c, and d are updated until the minimum of the objectiveof numerical phase velocity is defined as function  (B-2)is reached, where k is the iteration number, p is one of a, c, orthe  norm of the residual to be minimized can be expressedasd, and  c, d) =   a, c, d)  is the gradient vector of E(x).Using equation (B-2), the averaging coefficients can bewhere  =  =  = determined, minimizing the error defined as equation (B-l).

An optimal 9-point, finite-difference, frequency-space, 2-D scalar wave extrapolator

https://s-space.snu.ac.kr/bitstream/10371/6111/1/gphy1996-01.pdf

An optimal 9-point, finite-difference, frequency-space, 2-D scalar wave extrapolator

Abstract

Similar works

Full text

Available Versions

SNU Open Repository and Archive