Wavefront construction using waverays

WAVERAYS AND WAVEFRONTS

The product of merging the two previously discussed ray tracing methods is sketched in the following pseudo-code:

1 		 Allocate memory for wavefronts and output 

2 		 for each source 

3 		 		 initialize  

4 		 		 while number of points in wavefront $> 4$ 

5 		 		 		 propagate wavefront 

6 		 		 		 check for self crossing 

7 		 		 		 check for rays that go out of bounds 

8 		 		 		 eliminate self crossing and outofbound rays 

9 		 		 		 calculate amplitudes 

10 		 		 		 gridding 

11 		 		 		 interpolate new rays 

12 		 		 		 movwavf

The algorithm defines the data structure ``cube'' at line 1. In it, the ray parameters are stored as wavefronts propagate.

		 struct heptagon {     

struct point {    

float x;      

float z; } x0, x1;   

float angle;       

float ampl;        

char cf; } $*$cube;

where:

• cube[ii].x0 contains the starting position of a ray on a wavefront at a time $\tau$.
• cube[ii].x1 contains the arriving position of a ray on a wavefront at a time step later, $\tau + \Delta \tau$.
• cube[ii].angle gives the arriving angle at point cube[ii].x1.
• cube[ii].ampl is the computed amplitude at point cube[ii].x1.
• cube[ii].cf is a flag that defines if point cube[ii].x1 is in bounds, out of bounds, or belongs to the inner section of a self crossing wavefront.

The index ii runs over the points of a wavefront. Since there is no a priori way of determining how large can a wavefront grow over a velocity model, a predefined limit (nrmax) of cube elements is established at the beginning of the algorithm, which defines the total memory allocated for cube. In other words, nrmax is the maximum number of points that a wavefront can have. If the wavefront grows bigger than nrmax points the algorithm is stopped and an error message is produce indicating that a bigger value for nrmax should be used. This integer depends directly on the maximum allowed separation DSmax between two contiguous points in the wavefront. For the examples shown on Figure 9 through 12, the program ran with a value of nrmax=1300. DSmax was set to 21 meters for those examples.

In the case that the algorithm is run in a ``all arrivals'' mode, which could be done by eliminating line 6 out of the algorithm, the number of crossing points on a wavefront could become considerably large. As the wavefronts crosses and crosses many times over itself, for a velocity model with strong variations (as for example the Marmousi model), the number of crossing points can easily reach the 6 digits figure. This translates directly into a bigger need of computer resources, in use of memory and time.

Subroutine initialize defines the initial wavefront. It assigns an initial amplitude and take-off angle to the points on the initial wavefront.

The propagate wavefront subroutine ray traces using Lomax's waverays. The waverays are traced starting at cube[ii].x0 with a take-off angle cube[ii].angle during one time step at a certain frequency. The time step and the frequency are user predefined.

Subroutine on line 6 checks for self crossed wavefronts and flags the points that belong to the inner crossed section of the wavefront.

Line 7 checks for points that fall out of boundaries, raising a flag. Notice from Figure 8 that the rays cross over the boundaries of the model. This is done in order to obtain arrivals at the receivers that lie on the boundaries.

Subroutine at line 8 eliminates the points on a wavefront that are flagged for laying out of bounds or belong to self crossed wavefronts.

Subroutines on lines 9, 10 and 11 are implemented as previously explained for the NORSAR method on Figures 7, 6 and 4 respectively. On line 11 the number of rays that may be interpolated between any two contiguous rays, is given by the number of times the distance between the two rays is bigger than the maximum allowed distance DSmax.

Subroutine gridding is a very time consuming, due to the irregular distribution of the data in the model. First, the subroutine checks for receivers inside the ray cells of two contiguous wavefronts. If a receiver is found, the ray parameters are interpolated to it.

Subroutine movwavf prepares the structure cube for a new wavefront to be propagated. It takes as the new starting point the previous arriving point (cube[ii].x0 = cube[ii].x1).

Wavefront construction using waverays