Finite Element Reverse Time Imaging for Earthquake Sources and Scatterers

In this study, we applied the reversibility of elastic wave to a discussion of the reverse time propagation of wave equation. Based on the backward propga­ tion concept, we evalu�ted the feasibility of reconstructing the image of elastic wave sources and scatterers. The finite element. method was used to calculate elastic wave propagation in the medium. The same computational code was also used to reconstruct the. reverse time image. It was found that after our careful process­ ing of the absorption boundary for artificial reflectio,n, we could use the time-trace records at the surface to invert the locations of seismic sources or underground scat­ terer accurately. Our results show that point source, 'finite fault and multi-sources are well reconstructed. Although direct waves are much stronger than scattering wavefields, the imaging of the underground anomaly is still well reconstructed from the reverse time process. The image reconstruction technique based on the finite element method is flexible and easily to use for many applications.


INTRODUCTION 17
Earthqua ke sources and earth structure inversion are active research top ics in seismology.They may be determined by the extrapolation of observed wavefield backward in time to obtain an image of sources or niedia scatterers using numerical methods.The principle upon which the imaging of sources and scatterers is based is the reversibility �f the wave equation for both forward and backward propagation.Forward mo deling . of seislllic data .is an initia) value problem.A source is specified in a given model arid the respcmse is the wave field for any desired observation point.Since the wave equation is reversible, one can propagate the wavefield backward from the observation points in ; order to recover the wavefield at previous stages in both time and space.Reverse .time imaging, a kind of backward propagation, solves the wavefield of a given space model based on the wave equation using time traces for the observation points.

TAO
Vol.2, No.l The numerical code used for forward modeling can also be used for reverse time imaging (McMechan, 1982).N umerical modeling of a wave propagation has been well discussed.Boore (1972) used the finite difference method (FDM) to calculate wave propagation, Alterman and Aboudi (1970) simulated a dislocation source based on an equiv alent body force.Smith (1975) used the finite element method (FEM) to calcu late elastic wave propagation.Cerveny et al. (1982) used the Gaussian beam method (GBM) to trace the ray path.Benz and Smith (1987), and Huang (1989) used the FEM to simulate kinematic fault rupture processes.Kosloff et al. (1984) used the pseudo-spectrum method (PSM) to model wave propaga tion.
Reverse imaging for media reflectors is well developed in explosion seis mology.The wavefi.eldreverse time image technique for geophysical applications was mostly developed by McMechan and his co-workers (Chang and McMechan, 1987;Zhu and McMechan, 1988).However, most of their attention was focussed on petroleum geophysical applications and assumed explosive sources.Their targets were reflections from underground structures.They used the FD M to compute the wave equation because it saved computation time and was easy to calculate.In this study, our work is the determination of source image and media scatterers for earthquake seismology.A method for a source parameter inversion based on reverse time propagation was first proposed by McMechan (1982) and by McMechan et al. (1985) for real earthquake data, and extended to elastic cases by Hu and McMechan (1988).In order to consider the source complexity, structure heterogeneity, model flexibility and required numerical accuracy, we choose the FEM to formulate the elastic wave equation.In this study a finite element code for image reconstruction is developed first.Second, some interesting numerical cases are proposed and discussed.Finally, we dis cussed the possibility of using the reverse time imaging technique on real data (such as the SMART-1 array data) for source rupture inversion.

THEORY
Since forward modeling and reverse time imaging use the same concepts of wave propagation and share the same finite element codes, the forward modeling for the elastic wave propagation will first be presented, which then followed by its inverse.
a.The forward problem: Synthetic seismograms for 2-D in-plane models (P-SV and Rayleigh waves) can be deduced from the following equations of motion (Aki and Richards, 1980, p781), (I) where t denotes time, (u,w) is the displacement vector of horizontal (x or dis tance) and vertical (z or depth) components.f:x and fz are the body forces on the x and z components.a(x, z) is the compression velocity ((,\ + 2µ)/ p)1!2 and j3 ( x, z) is the shear velocity (µ / p ) 1 I 2 , where ,\ and µ are the Lame con stants, and p(x, z) is density.Equation ( 1) can be expressed as a finite element representation (Smith, 1975) as follows, . .

[M][
where [M] is the mass matrix, [CJ is the damping matrix, [K] is the stiffness matrix and [F] is a vector of point forces acting on the nodes.[A], [V] and [U] are respectively the acceleration, velocity and displacement vectors of the nodes.The wavefield can be extrapolated in its time domain with either explicit or implicit integration methods (Bathe and Wilson, 1976).Equation ( 2) is im plemented, for example in this study, as the central difference approximation (Wylie, 1975) and its finite element code is well checked with analytic method and other finite element code (Huang, 1989).To keep numerical dispersion down to an acceptable level, 10 grid points per wavelength at least are needed at the dominant source frequency for the S wave velocity and a damping matrix is used with a damping coefficient d defined by Archuleta and. Frazier (1978) from the 0.1 to 0.2 range.In order to suppress numerical reflections from the ar tificial boundaries of the finite element model, the absorbing boundary method proposed by Smith (1974) is used.This method involves the superposition of independently calculated wave equation solutions with Neumann or Dirichlet conditions.In theory, the summed wavefield completely eliminates artificial reflections.Using the absorbing boundary method (Smith, 1974), wave prop agation allows the computation on a small finite element model without being distorted by numerical reflections from artificial boundaries.The numerical re flections in Fig. l (b) are efficiently eliminated but do not affect the real signals (Fig. 1 (a)).As an important difference from the FDM and PSM, the FEM does not set any artificial boundary conditions on the free surface, the traction free boundary condition is automatically satisfied in the FEM formulation.Fig. 1 Comparisoii. of seismic profiles between (a) with, and (b) without absorption boundary condition (Smith, 1974) for a half-space and a dislocation point source. b.
The inverse problem: The reconstruction of an image of the source (inverse problem,) , utilizes the reverse of the forward modeling process.This forms a boundary value prob lem in which the boundary conditions are defined at the observation point {free surface or downhole receivers) .The complete boundary conditions need all of the explicitly varying values in time, that is, all of the information from seis mograms.The procedure for image reconstruction using a the finite element code is described as follows.The finite element code for reverse time imaging is the same as that for forward modeling but the input source locations are substituted by the receiver locations and the input source time function is re placed by the reversed seismic : time traces.First, reverse the recorded seismic traces in time (both horizontal (x) and vertical (z) component) as input data for the imaging code.Second, extrapolate the input data from the observation points at free surface back to the medium by pushing the horizontal and vertical component of finite element meshes respectively.Finally, the source image is identified in the model space from the wavefield concentration pattern of the snapshots (or time slices) .The snapshot display method is the key technique for reverse time imaging.The time slice display is necessary for the evaluation of image concentration.Fo llowing the reverse time imaging calculation, the wave equation carries all the energy backward along the same paths over which it was previously propagated.More examples and discussion about reverse time imaging can be found in a paper by McMechan (1982) for acoustic waves.The reconstruction of scatterer images differs only slightly from that of source im ages.Since, the direct P, S and conversion phases (where arrival times are predictable) do not contain any information about scatterers, we ignore those phases from seismograms first to isolate the scatterer wavefields in seismograms.
Then, the processed time traces are treated as input data for source imaging to get the scatterer image.

NUMERICAL EXAMPLES
a.
Point source imaging: The first example of this section, the forward modeling of the fre� surface displacements due to a dislocation point source with a 45° dip-slip angle was formulated by the split node technique (Melosh and Raefsky, 1981) is shown in Fig. 2. The grid size was 121 points horizontally and 61 points vertically, and the grid interval was 10 meters for the model in Fig. 2. A source time function of Gaussian pulse is used.:r" he two seismic profiles of Fig. 2 give the observations on a free surface which indicates the propagated wavefronts at different locations at different times.These seismograms include information about the relative locations of the source to observations.Figure 3(a) shows the outward propagating seismic wavefronts in a vertical plane at three widely separated time steps with time increasing from the lower panel to the upper one.
For the purposes of source image reconstruction, the seismic profiles of Fig. 2 are processed using the reverse time inverse procedure.The results are shown in Fig. 3(b) which shows the x-z wavefield three times during backward propagation.The snapshots have the same propagation times as for the forward modeling.The energy of the imaging wavefields converges to focus in time and space as a reverse of the forward � odeling (comparing sn � pshots in Fig. 3), but these wavefields used from inversion .are simply the upper half of those from the forward modeling, the lower half of the forward spherical wavefront is missing because only the energy propagated upward is recorded in the seismograms on a free surface.Therefore, the downward waves can rtot be recovered simply from the obserVa.tionsat a free surface.Additionally, diffraction tails can be seen in the reverse wavefields.This is because the, seismic profile data is of a finite aperture for the wave equation and the tails are produced by truncation of the spatial data.When the backward propagation continues through focused imaging to a negative time (pre-origin), the wavefields diverge again as shown in Fig. 4. The continuing image procedure will propagate the energy downgoing from the focal point.Therefore, using reverse time imaging, some criterion is required to determine an optimal source location and origin time.
Having demonstrated the characteristics of the imaging concept with a simple example, we shall now illustrate three numerical examples with more seismological meaning in the latter part of this section.These include earth-  b.Finite fault imaging: A first example of source imaging is the case of a line source of a 45° dip angle thrust fault.The velocity parameters are the same as for the point source model (Fig. 2). Figure 5 shows four representative time snapshots of forward and reverse propagation.The forward propagation snapshots (Fig. 5( a)) show that the major wavefront energy observed at a free surface is on the side of the hanging wall.The S wavefront is clearer than the P wavefront.The reverse propagation snapshots (Fig. 5(b)) show that the major reverse wavefronts came from the right hand side (hanging wall) and a backward propagating wavefronts approach to the true source location (from the right upper corner).The fault dip angle and its length are well reconstructed.Although the inverse results include only the upper half wavefield of source at the original time, the direction of faulting (or soµrce rupture direction) is well reproduced.The dominatant energy of reverse time snapshots are S wavefronts.The diffraction tails are weaker than that for point source example.It is also found that the finite seismic profile span is enough to identify the orientation and spatial extent of the rupture fault.The source with a finite spatial extent is usually observed in near distance strong motion records (Spudich and Frazer, 1984).This example illustrates that the method is suitable for finite fault image reconstruction.

c. Multiple source imaging:
T.Q,is �xample is used to illustrate the resolution of the imaging technique for.ip.ultiple point sources.Figu.re 6 contains representative forward and in verse wavefields for two dislocation point sources (left upper and right lower locations in the model).The model velocity parameters are the same as for the point source model (Fig. 2).The four snapshots from bottom to top in Fig. 6(a) show the wavefields extrapolated over time.The energy of ho.th the point source wavefronts interferes and superimposes on the free surface.Based on reverse time imaging, the interference seismograms are used to drive finite element meshes . .'t� reconstruct the source images.The inversion results of this back ward propag • ation show that the image sources are well separated and correctly focused in tirri.e and spac� as seen in Fig. 6(b)� This can be interpreted as the superposition property of wave propagation.Although the wavefields observed at the free surface"are mixed in the records, they can be extrapolated backward along individual paths to foci.It is found that the resolution of the source image was affected by the frequency content, i.e. the higher the frequency, the greater the resolution.Because of the limitation of a .finiterecording aperture,  -----------------.:i: 0.2 x :i:: .... Q. w 0.4 c 0. 6 '---'--J----'--,    niques for the geophysical engineers.Fig. 7 The model for scatterers imaging problems.The square anomaly has a density of 0.1 g/cm3, P wave velocity of 2.03 km/sec and S wave velocity of 0.01 km/sec.The point source at a depth of 0.02 km simulated the surface control source.---'-----' '---'--'----'-'---'-----' : Fig. 9 Snapshots of the wavefield during reverse time propagation.The input data are the muted seismograms of scattering waves as shown in Fig. 8.Many interesting waves are identified.The symbols are defined as P for direct P wave, S for direct S wave, PS for P conyerts to a S wave, PP for P converts to a P wave, SP for S converts to a P wave, SS for S converts to a S wave, R for Rayleigh wave.
late the control source.Two component seismic records calculated from the forward modeling are used as input data for reverse time imaging.No anomaly is assumed, thus a homogeneous half space model is used for the inverse procedure.From our computation, the wavefields are gradually extrap olated into the model as shown in Fig. 9.The eight snapshots in Fig. 9 from lower panel to the upper one and from left to right show the discrete time steps of the wavefield backwardly propagated from pre-focus to post-divergence.As discussed in the point source example, there is no absolute time for the focus, hence, the wavefields will diverge after focusing ..The image of the most mas sive wavefield is recognized as the location for the scatterer or anomaly (at the time slice t = 0.405 sec ).One forward modeling snapshot is shown in Fig. 10 ..Although the direct waves (phase P, S, PS, R as seen in Fig. 10) are much stronger than scattering wavefields.(phase PP, SP and SS), the imaging of the anomaly is still well reconstructed by the reverse time process (Fig. 9).The computational accuracy of free surface response is important for forward modeling and for reverse time imaging.The traction boundary condition of the free surface is automatically satisfied in the finite element formulation for wave equations but is not always so suitable for the FDM or PSM (Kosloff et al., 1990).The accuracy of the free surface amplitude response computed by the FEM is better than for the FDM.The discretization of complex models for wave propagation using the FEM is more suitably accomplished than using the FDM.Therefore, the advantages of FEM for reverse time inversion primarily lie in the flexibility of its model description and the accuracy of energy recovery for complex cases.The numerical cases for the point source, finite fault and multi-source are well reconstructed using the FEM.This example using multiple sources is important, as when we treat the multiple sources as secondary sources in media, for discussing the resolution of scatterers wavefield inversion.Currently, the seismic refraction or reflection methods> are used to sense underground holes or basement topography, but the two methods are only use ful over certain suitable distance ranges and are subject to other limitations, as discussed by Sharma (1976).The seismic data processes using the ray tracing and refraction methods which are usually used for geophysical interpretations only applies the first arrival times.However, seismogram information is still not well utilized by these two methods, thus, the seismic record waveform carries not only a travel time but information about underground structures.Herein, the reverse time imaging technique proposes a concept to apply the seismic waveform.This method is based on the single shot data and processes seismic data directly without regard to phases identification.The information, it em ploys, is the waveform (multiple phases and amplitudes).Otherwise, the image reconstruction technique based on the finite element method is easily used in other applications.For example, underwater moving source detection, high res olution underground anomaly detection, nondestructive testing and biomedical applications.
One of the important purposes of this study is to evaluate the feasibility of the reverse time image technique for earthquake source inversion.The reso lution of the images produced in the examples above are, in some ways, better than might be expected for real data under the less-than-ideal conditions of the profile span angle, larger offset and complex rupture processes.It has also been found that, in elastic cases, the S waveform is very important for reverse time imaging.Most of the focus energy in a snapshot is in the S wavefield.If the ob servations do not include or only include a part (for example, only the vertical component) of the S wave energy, reverse time imaging is not complete and im-age resolution is poor.Therefore, 3-component seismograms including complete S wave waveforms are important for this method.The other limitation of the reverse time imaging method is the large computation time and storage mem ory needed to compute wave propagation by numerical methods� This weak point also restrlcts• the high frequency signal inversion and the resolution of the image targets.In any' case, the image from the reverse time reconstruction is only a picture, which can not support to give us the quantitative resolution of the image.The statistical properties of the image require a systematic in version analysis.Based on the numerical examples and discussion above, we suggest that before applying this method to process SMART-1 seismic array data• for source imaging and ruptures, some difficulties need to be overcome.First, the basin structure of the Ilan plain needs detailed investigation.Second, the aperture . of this array needs to be extended across the basin.Finally, for the complex 3-D structure of the Ilan plain (Chiang, 1976), a 3-D reverse time imaging is needed.Even so, with development of computer engineering tech niques and further installation of a near fault high dense array for source study, the imaging technique will surely be widely applied in seismology in the near future.
In summary: (1) The image reconstruction technique based on the finite element method is flexible and easily used for many application fields.
(2) The earthquake source is well imaged using the reverse time image tech nique.
(3) Some difficulties need to be overcome to apply the method to the SMART-1 array data.
Fig. 2 The forward modeling for a dislocation point source in a half-space.(a) horizontal, (b) vertical displacement profiles, and (c) model parameters and source location.The observations are along the line of the free surface. Fig.3 Snapshots of (a) forward and (b) reverse propagation for a point source model (Fig.2).The elastic wavefields are expressed in a vector format.The same as other snapshots display in this study.Stratigraphic snapshots of the energy wavefield during backward propagation.The wavefield at the negative time represents propagation over the origin time.TAO Vol.2, No.1 quake sources (line source, multiple point source) and scatterers (underground anomalies) imaging; The input data illustrated for the reverse time imaging are all generated from forward modeling by the same finite element code.
Fig. 5 Snapshots of (a) forward and (b) reverse propagation in line source model.The model parameters are shown in Fig. 2(c).The source is a thrust fault and located at the center of the model with a dip angle 45° to right.

d.
Scatterer imaging:•The follow�ng example demonstrates the application of reverse time tech-
Figure 8 shows the horizontal component profile.The strongest signals can be recog nized as direct waves which propagate from the source to receivers at the free surface (Fig. 8 (a)) .The direct waves are muted with a linear taper as shown in Fig. 8 (b).The data, horizontal (Fig. 8 (b)) and processed vertical seismic profile (not shown) , are synchronously driving the finite element meshes, respectively.
AND CONCLUSIONSMost observations of natural earthquakes are made on the free surface.