Network Distributed Parallel Processing of Pre-Stack Layer-Stripping Reverse-Time Migration

Pre-stack layer-stripping reverse-time migration is successful for struc­ tures with steep dips and strong velocity contrasts. This technique uses reverse-time method to migrate seismic sections through constant or smoothly varying velocity layers, one layer at a time. The final migration result is composited from the individual layers' images. The most important advan·­ tage of the pre-stack layer-stripping reverse-time migration method is that it makes interpretation become a step of migration. However, although su­ per computers have made computation less of a problem, applying pre-stack layer-stripping reverse-time migration to a large seismic model still requires significant computation expense. With the help of modern computer network techniques, we have suc­ cessfully developed a cheap and high efficient strategy to perform pre-stack layer-stripping reverse-time migration. This technique uses several networked workstations to share the works of migration. Migration is distributed to dif­ ferent machines through network and all machines share the same data file. The networked distributed parallel processing technique greatly increases the computational efficiency, reduces the computational cost, and shares compu­ tational resources available. Comparing the performance of migration on 6 clustered VAX workstations and a CONVEX C220 mini super computer, the total computer time on both systems are almost identical. Using the network distributed parallel processing techniqμe, we may develop time consuming data processing algorithm with small computer systems. This technique is particularly important for 3-P data processing.


INTRODUCTION
It has been long recognized that CMP stacking is exact only for horizontal layer me dia.In areas of complex structures, the stacking procedure may discriminate against steeply dipping events, because different dipping events have different apparent stacking velocities.Some of the CMP stacking effects can be corrected by applying dip moveout (DMO) correc tion to common-offset section after it has been NMO corrected (Yilmaz and Claerbout, 1980;Deregowski and Rocca, 1981;Bolondi et a� 1982;Hale, 1984).However, the use of DMO correction is restricted either by the asSuil1ption of small-offset (Yilmaz and Claerbout, 1980) or constant velocity medium (Hale, 1984).The most important drawback of the DMO tech nique is that it still can not perform accurately in areas having mild or strong lateral velocity variations.Under these circumstances, pre-stack depth migration becomes the ultimate tool for imaging subsurface structures correctly.
The pre-stack migration result can be obtained by common-shot gather migration fol lowed by stacking.In each common-shot gather migration, only the wavefield representing the receivers is extrapolated; migration image is formed at the reflector rather than at time zero.The complete migration result is a superposition of these individual migration images.
Traditionally, wavefield extrapolations are achieved by downward continuation of sur face data toward increasing depth into an earth model.In that case, the knowledge of the surface recording alone is insufficient to provide a unique solution to solve the scalar wave equation.In addition to the surface recording; its normal derivative is needed as well.For practical use, the scalar wave equation has to be approximated in some ways.For exam ple, Claerbout applied the parabolic approximations to the wave equations as a wavefield extrapolation operator, which allows migration to be performed using the surface recorded information only.However.using the parabolic wave equation, the wavefi eld can only be propagated accurately within a limited.range of propagation angles.
Reverse-time migration is a depth migration technique.This method permits a different approach to migration that is based on reverse-time wavefi eld calculation instead of depth extrapolation.This migration technique has shown superior imaging capabilities over the conventional dip-limited finite-difference depth migration algorithms (Gazdag and Carrizo, 1986;Kosloff and Baysal, 1983;Loewenthal and Mufti, 1983).However, the high expense of reverse-time migration prohibits its use as a general migration algorithm at this time, particularly if strong velocity contrasts are present in the model.The other drawback of reverse-time migration, conventional migration methods as well, is that the users have to provide a complete velocity model before migration.But the detailed velocity is not always available.
In order to increase the efficiency of reverse-time migration and make interpretation a step of migration, Shih (1990Shih ( , 1991Shih ( , and 1992) ) has successfully developed a layer-stripping reverse-time migration technique.The method uses the reverse-time algorithm to migrate the seismic section layer by layer.The final migrated image is composited from the individual layer migrations.The pre-stack layer-stripping reverse-time migration technique has proven to be successful for structures with steep dips and strong velocity contrasts.In general, pre stack migration deals with a large amount of seismic data.However, although supercomputers make computation less of a problem, applying this algorithm to a large seismic model requires significant computational expense.
In this paper, pre-stack layer-stripping reverse-time migration is presented using a newly developed technique, named network distributed parallel processing pre-stack layer-stripping reverse-time migration.Using this new technique, we are able to distribute migration through computer network to different computer CPU and migrate simultaneously.All machines use the same migration program and share the same input and output data files .This algorithm has successfully ported on several networked VAX 3100 and 4100 VMS system workstations at the Exploration & Development Research Center, Chinese Petroleum Corp. in Miaoli, Taiwan.An outline of the migration algorithm is described in the following sections.

PRE-STACK LAYER-STRIPPING REVERSE-TIME MIGRATION
Assume that a two-dimensional velocity structure can be adequately described by a series of layers, each of which has constant or smoothly varying velocity together with a laterally varying boundary.The layer-stripping migration technique involves numerically solving a boundary value problem for the source-free wave equation over a series of layers.
In each layer the extrapolated wavefield P ( x, z, t) satisfies the scalar wave equation and 8 2 P(x,z,t) 8 2 P(x,z,t) 1 8 2 P(x,z,t) vi(x, z)2 8t2 for (x, z) r/.R, P(x,z,t) = PR(x,z,t), for (x, z) ER. (1) ( 2 ) Here R is the top boWidary of the layer i, Vi ( x, z) is a constant or slowly chang ing velocity field of that layer with horizontal coordinate x and vertical coordinate z, and PR(x, z, t) is the recorded wavefield along the boundary R.
In principle, the observed wavefield PR(x, z, t) and its normal derivative are needed to calculate the extrapolated wavefield P( x, z, t ).From a practical point of view, the extrapolated P ( x, z, t) can be computed in a number of ways by using the observed boundary condition PR(x, z, t) only (Berryhill, 1979;Wiggins, 1984;Yilmaz and Lucas. 1986;Hill and Wuensche!, 1985;Clayton and McMechan, 1981).Here we use the reverse-time method to extrapolate the wavefield.The reverse-time method is implemented by runn ing a finite difference algorithm backward in time and using the time reversed observed wavefield as boundary value at the receivers.The wavefield extrapolated to the redefined boundary is represented by the amplitude of the wavefield along that boWidary at every computational time step.Let Ei denote a wavefield extrapolation operator a 2 a2 i a 2 Ei -( -+ -) ---,... -- Wavefield extrapolated layer by layer in the layer-stripping reverse-time migration scheme can be expressed as -� Ei P(x, z,t)=O,for i=l, ... ,N, where N is the number of layers.Suppose that the finite-difference spatial and temporal grid are defined as follows: Xm = mD.. x,m = 0,1, ... ,M, Zn = nD.. z, n = 0, 1, ... , N, and ti = 16.t,l = 0, 1, ... , L-, (5) where Ax = Az = Ah is spatial grid spacing and At is time step.Based on the 2nd-order central-difference approximation, wave equation (1) can be written as .v(m,n)2At2 P(m, n, I + 1) = 2P(m, n, 1) -P(m, n, I -1) -

Ah2
[P(m + 1, n , l) + P(m -1, n, l) + P(m,n + 1,1) + P(m,n-1,1)-4P(m,n,l)].( 6) Equation ( 6) can be used for migration with the boundary values recorded at the earth's surface.For migration, the surface recording P( x, z = 0, t), for 0 < t < TL, is used as the boundary condition and is inserted at the physical boundary z=O, and extrapolated backward in time.The boundary values at TL and TL+ 1 are also needed to solve Equation (6), where TL is the last recorded time sample.It is assumed that P(x, z, t) = 0 for all t > TL.After time TL, the wave has traveled away from the medium, and the earth model below the last seismic arrival is assumed to be homogeneous and infinite.Equation ( 6) is formulated from the exact two-way wave equation.Solutions of this equation allow energy to be propagated in all directions.In this case, artifi cial interlayer multiples are generated from where the velocity varies rapidly.Since only a constant or smoothly varying velocity function is assigned to each layer migration, the layer-stripping method prevents artificial interlayer multiples.
Practically, layer-stripping reverse-time migration consists of layer definition, reverse time migration and boundary determination three steps.For a velocity model with strong velocity contrasts between layers, the initial velocity model is subdivided into several dis crete layers with approximate boundaries (see Figure 1 ).Each layer has a constant or slowly varying velocity.The boundaries can have any geometry, and are assumed to be interfaces separating layers with large velocity contrasts.Before migration, only the uppermost bound ary is assumed to be known exactly.The bottom boundary of the layer is approximated.After defining the layer.we migrate the seismic data in each layer with the use of the .reverse-time migration method.The migration application involves reverse-time wavefield extrapolation and followed by the imaging fonnation.To extrapolate the wavefield in each layer.we use a 10th-order space.4th-order time explicit finite-diff erence scheme to approxi mate the two way acoustic wave equation (Dablain, 1986).Besides, we apply a sponge layer transmission boundary (Israeli and Orzag, 1981) around the computation region to prevent reflections from the edge of the grid.In order to maintain nondispersive propagation in the grid interior, we need to sample the 10th-order finite-difference scheme at 3 grid points per wavelength.Based on the above criterion, the layer-stripping method allows the grid space to be scaled appropriately for the velocity of each layer.The grid space can be broadened in the high velocity layers, thus reducing the cost of the computation.To honor the original experiment geometry, we keep the horizontal grid equal to the receiver group spacing in migration.
After propagating the surface recordings back to the earth medium, we form the mi gration image •according to an imaging condition.To construct the migration image, we may use a ray tracing method to calculate the travel time from the source to each finite difference grid point, and then sample the reflected wavefi eld at the downgoing wavefront position.This imaging condition bases on Claerbout's (1971) idea that reflectors exist where upcoming and downgoing are time coincident.The final migrated section will be given by R( x, z) = U ( x, z, t( x, z )), where R represents the reflectivity distribution, U denotes the upgoing wavefi eld, and t( x, z) is the one-way travel time from the source to the finite difference grid point ( x, z ).

THE NETWORK DISTRIBUTED PARALLEL PROCESSING MET HOD
Due to the robustness of reverse-time migration, we have chosen the reverse-time wavefield exploration method for layer-stripping migration algorithm.Unfortunately, the band-limiting and stability criteria of the finite-difference scheme limit the computational efficiency for migration.Using a higher-order accurate finite-difference scheme in a forward modeling problem, we can broaden the grid spacing and increase the calculation efficiency (Dablain, 1986).In migration, the trace spacing is kept the same as the recording group spac ing.Therefore, reverse-time migration needs to use a higher-order accurate finite-difference scheme to prevent grid dispersion.Because of the long oper ator length, using a high-order finite-difference scheme makes migration become more time• consuming.
The presented migration method is implemented with an explicit scheme 10th-order space, 4th-order time accurate finite-difference approximation to the scalar wave equation.The pre-stack layer-stripping reverse-time migration program was first developed on an ELXSI 6400 mini super computer, and then paned on a CONVEX C220 mini super com puter.Although super computers have made computation less of a problem, applying this algorithm to a large seismic model still requires signifi cant computer CPU time, memory and disk storage space.For some research institutes or industry corporations, an expensive super computer is not always available, therefore, the usefulness of pre-stack layer-stripping reverse-time migration is strictly limited.
With the help of modem computer network techniques, we have successfully developed a cheap and high efficient strategy to perform pre-stack layer-stripping reverse-time migration.This method uses several networked workstations to share the work of pre-stack migration.Each shot record migration is distributed, through computer network, to different workstations and migrated simultaneously.In addition, all machines share the same migration program and seismic data.
The proposed network distributed parallel processing technique can be treated as a manual computer parallelization method.The basic requirement of the computer system is that the computer file system has to be network clustered together.In other words, the disk storage could be shared from network, if the problem we are solving can be subdivided into independent tasks.For instance, in pre-stack migration, each shot record can be migrated independently.Then we distribute the separated tasks to different computers assigned and run all jobs at the same time.In this case, the total CPU time can be shared by different machines.Since the disk space is shared by all machines, only one single copy is needed for all machines.As long as we track the data positions correctly, only one data file is needed.That algorithm greatly reduces the disk space.
We have successfully apply the network distributed parallel processing pre-stack layer stripping reverse-time migration program to numbers of VAX 3100 and 4100 VMS system workstations.The first step of this technique is to log on to every machine available and stan the same migration job.Once the jobs are submitted, migrations on each machine will stan from the first shot record one at a time.Since all machines are working on the same data set, the fi rst shot record will be accessed more than once.To make sure all the shot records are processed once only, we provide one shot record only to one machine and first come first served.To do this, we have also created a flag file together with the seismic data file.In that flag file, a flag one is assigned to each shot record.If one shot record is in use or migrated by a machine, the flag associated with that shot record will be turned to zero.Once the flag is turned, other machines won't be able to access that shot record and will automatically skip to the next shot record available.When a shot record is migrated the program writes the migration image to disk and search for the next unmigrated shot record.The job is completed until all shot records are migrated.The same as the input data, if we could track the data position correctly, only one fi le for the output data is needed. .The advantages of the networked distributed parallel processing pre-stack layer-stripping reverse-time migration technique are: first, it saves the cost of buying and maintaining super computers; second, it greatly increases the computational efficiency; third, it shares all usable computer resources; and the last, it makes time consuming processes be able to be developed on a small computer system, which is particularly important for 3-D data processing.
The networked distributed parallel processing pre-stack layer-stripping reverse-time mi gration program can be run on one single machine or as many machines as possible.The more machines we use, the higher computational efficiency we gain.The success of the networked distributed parallel processing pre-stack layer-stripping reverse-time migration scheme will be demonstrated by migration of a synthetic example in the next section.The illustration fo cuses on obtaining high computer efficiency, reflector position and shape rather than detailed stratigraphic structural interpretation.

SYNTHETIC DATA EXAMPLE
In this example, 46 2-D synthetic shot records were generated from the Sierra ray-tracing forward modeling programs QUIKSHOT and QUIKDIF.QUIKSHOT performs non-normal incidence raytracing for simulation of field shot records.QUIKDIF adds on the simulation of true diffraction amplitude from faulted or pinched surfaces.The complex velocity model is shown in Figure 2. The model is 18 km wide and 6 km deep, and is composed of 10 irregular layers with velocities from 2,840 m/sec at the surface to 4,320 m/sec at the bottom.The principal structural features in the model are an anticline and a thrust fault.46 synthetic shot records were generated, each 4 seconds in length with a 4 ms sampling interval.The first shot is at 3570 m surface position, the shot interval is 240 m.Receivers are evenly distributed at 60 m intervals.This geometry simulates a 30 m midpoint interval CMP survey.
Figure 3 shows one of the synthetic shot records.A spherical divergence correction and a 40 Hz low-pass filter have been applied to each shot record before migration.Figure 4 shows the CMP stacked section for this model, average fold is 15.Prominent features on the CMP section are: diffraction hyperbolas which originate from the tip of the anticlines and the intersections of two events and velocity pull-up below the thrust fault.( where Ts( m, n) represents the one-way travel time from the shot point to point (m.6.x,n.6.z).
in the next layer, T8( i, j) is the one-way travel time from the shot point to point (i.6.x,j .6.z) on the layer boundary, and Tt1:in is the travel time between boundary point (i.6.x,j.6.z) and point (m.6.x,n.6.z).The same method is repeated for each layer.In migration, when the calculated imaging time lies between two computational time steps, a linear interpolation is used to detennine the amplitude of the wavefield at that specific time.Repeating the same migration procedure, Figure 8 shows the migration images for the second layers velocity model, when velocity 3600 m/sec for the second layer was used .The horizontal grid spacing used for migration in this layer is kept the same as 60 m and the vertical grid size is 30 m.Since the velocity function is closer than the true one and most of the destructive noises are cancelled out, the migration image is better focused.
Using the same migration sequence, a migration composited from 3 layers is displayed in Figure 9.The velocity used for the third layer model is 4,000 m/sec.In Figure 9, the overall image is reasonably migrated correctly.
The above migration was migrated by using five VAX-3100 and one VAX-4000 VMS system workstations.To see the performance of the network distributed parallel processing pre-stack layer-stripping reverse-time migration program on the clustered workstations.We have also migrated the same data set on a CONVEX C220 mini super computer at National ChungCheng University.The result shows that the both systems produce the same migration images and use almost the same total computer time.Since the cost of workstations is much cheaper than a super computer, the use of the networked distributed parallel processing pre-stack layer-stripping reverse-time migration becomes very practical.

Fig. 1 .
Fig. 1.A velocity model with strong velocity contrasts between layers is sub divided into several discrete layers with approximate boundaries.Each layer has a constant or slowly varying velocity.The boundaries can be any geometry, and are assumed to be interfaces separating layers with large velocity contrasts.
hnage of the pre-stack migration consists of the amplitude of the wavefield at the direct wave arri val time from the source to each finite-difference grid.The final migration result is the superposition of these collected amplitude values.If the migration result is considered satisfactory, we then go to the step of boundary detennination.Otherwise, we readjust the velocity and migrate this layer again until a satisfactory image is obtained.Although we have chosen the bottom boundary of each layer when defining the initial model, the initial locations of the bottom boundaries are approximated.The last step of layer stripping reverse-time migration is to detennine the proper position of the bottom boundary for a given layer.This boundary is also the top boundary for the next layer.To find the exact boundary at the bottom of each layer, we scan a zone of seismic data around the initial guess for the interface.Picking the maximum square value of the migrated wavefield for each trace determines the proper interface position.Where the migrated signal is weak or ambiguous, the user can guide the choice of the boundary location.In field data migration, the boundary definition procedure requires a small amount of user smoothing.After detennining the bottom boundary of a layer, we collect a new seismic section along the interface.The collected seismic time section serves as the boundary value for migration in the next layer.Repeating the above procedures, layer by layer, we composite the individual layer images to obtain the final migration image.The above boundary determination step allows the user some control over the migration image.

Fig. 2 .Fig. 3 .Fig. 4 .
Fig. 2. The 10 layers velocity model used to generate synthetic shot records The model is 18 km wide and 6 km deep with velocities from 2,840 m/sec at the surface to 4,320 m/sec at the bottom of that velocity modeL

Figure 5
Figure5shows the velocity model used for pre-stack layer-stripping reverse-time migra tion.This velocity model consists of 3 layers with velocities 3,000 m/sec.3,600 m/sec.and 4,000 m/sec, respectively.The model is chosen to approximate the original velocity model shown in Figure2.Applying the layer-stripping scheme, we first assigned the velocity model as one single layer, 18 km wide and 6 km deep.The constant velocity reverse-time migration is performed using a rectangular grid.The horizontal grid size is 60 m, the receiver interval, and the vertical grid spacing for migration in the first layer is chosen as 30 m.As mentioned above, common-shot ga�er migration requires an imaging condition different from time zero.The imaging condition is the one-way travel time from the shot point to each finite-difference grid points and can be obtained by ray tracing.Because constant velocity layers are used in the migrations presented here, Fennat's principle is expressed as the simple square-root equationv'x2 + z2t= 'Vi

Fig. 5 .
Fig. 5.The 3 layers velocity model used for pre-stack layer-stripping reverse-time migration.The velocities are 3,000 m/sec, 3,600 m/sec, and 4,000 m/sec for the first.second and third layer, respectively.

Figure 6
Figure6shows migration images from the individual shot gathers, when the velocity 3000 mfw; is used.Observation from a single shot point provides a partial image.To obtain a complete image, we composite all shots.However, as seen in the individual migration ' results, migration smiles and grid dispersions from the lateral traveling waves mask the true migration image.The composited image will be strongly distorted if these large amplitude images are also added.To obtain a clear migration result, we have muted images at the large offsets before stacking.The same mute is applied to each partial image and the choice of the muting range is not critical.A number of mute patterns are applied to the individual migration images until the best image is obtained.Figure7shows the composite image for the first layer from 46 shot records.The shape and the position of the reflector is correctly imaged.Since the surface recordings are migrated to their midpoint position in the spatial domain, interfaces on both ends of the velocity model are invisible.Because most of the shot points are in the center of the velocity model, the right and left flanks of the anticline has smaller amplitude.Following the method we have described in the previous section, we determined the location of the lower boundary of the layer.The event representing the redefined boundary is treated as part of the migration image for this layer.For migration in the next layer, seismic signals prior to the one-way travel time from the source to the newly defined boundary are muted in the new boundary condition.A cosine taper is applied to the end of the mute pattern.After being scaled with the transmission coefficient, the extrapolated time sections are used as the boundary values for migration in the next layer.

Fig. 8 .
Fig. 8.The migration images of the 2 layers velocity model.The velocities for the first and second layer are 3000 rn/sec and 3600 rn/sec, respectively.

Fig. 9 .
Fig. 9.The final image of the 3 layers velocity model.The velocities for the first, second and third layer are 3000 rn/sec, 3600 rn/sec and 4000 rn/sec, respectively.