KIRCHHOFF ALGORITHMS

Kirchhoff migration (Schneider, 1978) is applied to the data as a convolution with a space and time variant operator. The impulse response of Kirchhoff migration, although it may be very long, is continuos in space. By exploiting the spatial continuity of the impulse response, an efficient implementation of Kirchhoff migration on a massively parallel computer needs to use only fast nearest neighbor communications. The algorithm is very simple. The seismic traces are laid out local to each processor (Figure ); the input traces are passed around among processors while the output traces are kept fixed (it may be as well the contrary). At each step of this process, each input trace will be in the same processor of an output trace, and thus its contributions to the overlapping output trace can be added to the result by simple in-processor computations.

The following code fragments express in Fortran 90 a Kirchhoff prestack time-migration algorithm. In this example the data are assumed to be sorted in common-offset sections and the traveltimes are computed analytically assuming an RMS slowness function. Furthermore, for the sake of simplicity, I did not include the constant phase shift that should be applied to the results of Kirchoff migration.

First I show how the desired data layout is specified to the Fortran compiler.

C
C arrays declarations
C
      REAL DATAL(N1,N2)
      REAL DATAR(N1,N2)
      REAL MIGRES(NTAU,N2)

C
C arrays layout 
C
C :SERIAL = local axis
C :NEWS = parallel axis
C

CMF$ LAYOUT DATAL(:SERIAL,:NEWS)
CMF$ LAYOUT DATAR(:SERIAL,:NEWS)
CMF$ LAYOUT MIGRES(:SERIAL,:NEWS)

Then, I show the main loop of the algorithm.

C
C initialize buffers for right-side half of operator
C
      DATAR=DATAL

C
C loop over operator lag
C
      DO IOPER=1,LOPER

C
C compute constants on the front-end
C	
	DELMIDPOINT=(IOPER)*D2
	DELMH2=(DELMIDPOINT-HOFFSET)**2
	DELPH2=(DELMIDPOINT+HOFFSET)**2

C
C shift input data to the left and to the right
C
	DATAL=CSHIFT(DATAL,DIM=2,SHIFT=1)
	DATAR=CSHIFT(DATAR,DIM=2,SHIFT=-1)

C
C loop over output times
C
	DO ITAU=3,NTAU

C
C compute traveltimes on the front-end
C
	   TIME=SLOWTAU(ITAU) *
     1          (SQRT(DEPTH2(ITAU)+(DELMH2)) +
     2	        SQRT(DEPTH2(ITAU)+(DELPH2)))
	   TIME=((TIME-O1)*INVD1)+1
	   ITIME=INT(TIME)

C
C compute amplitudes (your favorite amplitude function can go here)
C
    	   AMP=1./TIME

IF (ITIME .LE. (N1-3)) THEN

C
C compute interpolation samples and weights on the front-end,
C using the interpolator table TABINT
C
		ITAB=INT((TIME-ITIME)*LTABINT)+1

INDEX1=ITIME
		INDEX2=ITIME+1
		INDEX3=ITIME-1
		INDEX4=ITIME+2
		INDEX5=ITIME-2
		INDEX6=ITIME+3

INTER1=AMP*TABINT(0,ITAB)
		INTER2=AMP*TABINT(-1,ITAB)
		iNTER3=AMP*TABINT(1,ITAB)
		INTER4=AMP*TABINT(-2,ITAB)
		INTER5=AMP*TABINT(2,ITAB)
		INTER6=AMP*TABINT(-3,ITAB)

C
C compute contributions to the output and sum into output array
C

MIGRES(ITAU,:)=MIGRES(ITAU,:)
     1		     +	INTER1*DATAL(INDEX1,:)
     2		     +	INTER1*DATAR(INDEX1,:)
     3		     +	INTER2*DATAL(INDEX2,:)
     4		     +	INTER2*DATAR(INDEX2,:)
     5		     +	INTER3*DATAL(INDEX3,:)
     6		     +	INTER3*DATAR(INDEX3,:)
     7		     +	INTER4*DATAL(INDEX4,:)
     8		     +	INTER4*DATAR(INDEX4,:)
     9		     +	INTER5*DATAL(INDEX5,:)
     1		     +	INTER5*DATAR(INDEX5,:)
     2		     +	INTER6*DATAL(INDEX6,:)
     3		     +	INTER6*DATAR(INDEX6,:)

END IF

END DO

END DO

In this simple case each processor does exactly the same operation because the traveltimes are independent of the midpoint location. The same would be true if dip moveout were implemented instead of migration. On the contrary, in a depth migration algorithm the impulse response changes with the midpoint location, and thus each processor needs to perform a slightly different operation. More precisely, the indices of input time samples used in the summation differ from processor to processor. This requires the use of the indirect addressing capabilities of the processing hardware, and it usually leads to an increase in run time. Nevertheless, Kirchhoff depth migration can be efficiently run on a massively parallel computer (Van Trier, 1991). If the data are not evenly sampled, as in 3-D surveys, implementing a Kirchhoff migration algorithm would become even more complicated. However, also in this case, the basic idea of exploiting the spatial continuity of the impulse response can be used for deriving efficient algorithms.