Ray tracing in Matlab PDF

Preview

Summary

matlab...

Description

Ray tracing in Matlab Ruiqing He

ABSTRACT Ray tracing is an important procedure that is needed often in geophysics. Matlab is an excellent scientific and engineering software tool. This report introduces a ray tracer using Matlab. The ray tracer is so flexible that it can trace transmissions, reflections and multiples easily in both 2D and 3D. Several ray tracing experiments were conducted and the computed traveltimes were compared with the theoretical travel times and those from eikonal ray tracer. Results showed this ray tracer to be more capable than the eikonal ray tracer. This ray tracer has been successfully used in both tomography and migration algorithms applied to both synthetic and real data.

INTRODUCTION Raytracing is so vital to geophysics that it has enjoyed a long history of research. There are many ray tracing methods that can be categorized as initial value ray tracing (like the analytical method and shooting method) and boundary value ray tracing (like the bending method and perturbation method). For application purposes, we often need both travel times and ray paths, sometimes, reflections and multiples as well. To simplify the ray tracing problem, suppose we are given a digitized velocity model, how can we find all of the first arrivals and their corresponding travel paths to each node? Moser (1991), Fischer (1993) and others used a coarse grid to model velocity while using a finer grid to model raypaths. Realizing a finer velocity model is more important than a finer raypath model in ray tracing, I use the opposite approach as Moser and Fischer did, i.e. the velocity grid is equal or finer than the raypath grid (in this report I use equal grid). Without losing much accuracy, I assume that any ray path consists of one or several line segments which have the nodes as end points. A longer line segment is usually not necessary, so we can limit the line segment length to several grid interval (also called the sub-ray tracing length). This means any ray traveling from the source to a geophone must rest on nodes within the sub-ray tracing length. With the above assumption, the method I am going to present is based on the general Fermat’s principle, i.e. rays travel along minimum travel time paths.

1

RAYTRACING PROCESS The proposed ray tracing process iteratively computes sub rays, searches for other possible raypaths and then executes a decision. It has only four simple steps: Step 1: starting from the original source, we first trace the straight line paths to the nodes within the sub ray tracing length. All the visited nodes record the direct travel time, and remember the source coordinates. Step 2: then we can decide among these visited nodes, the one with the minimum travel time must has the first arrival it records, and the ray comes from the source it remembers. Step 3: use the decided node as a secondary source and repeat the sub-ray tracing. Any node that is visited twice is compared with its two travel times. Only the minimum one is kept, and the remembered source is always the one with the minimum travel time, no matter if it is the original source or a secondary source. Step 4: repeat from step 2 until all interested nodes (or geophones) are decided. This ray tracing process is so flexible that it makes reflection and multiple ray tracing very easy and straightforward. Suppose we are required to trace the reflections from one subsurface layer back to the ground, first we can apply flags to identify the nodes beneath the reflector so that they are not to be traced. Do one ray tracing until all the nodes on the reflector are decided. Stop here and keep only the travel times on the reflector nodes while refreshing all other nodes’ travel times, (like an exploding reflector, but they are triggered at their first arrival time). Do one more ray tracing directly from step 2 in the above loop. So one reflection rayt racing appears as two ray tracing procedures. A similar approach can be taken to trace multiples, so that one first-order multiple ray tracing appears as 4 transmission ray tracing procedures, and one second-order multiple ray tracing appears as 6 transmission ray tracing procedures, and so on. A similar strategy can be applied in 3D ray tracing.

EXAMPLES I now present some ray tracing examples. The first example is a linear gradient model. The grid size is 100x100 at 1m interval. The top layer velocity is 1000 m/s, and the gradient gives a 10 m/s increase at every 1m depth. Theoretically we can calculate the first arrivals on the ground. Figure 1 illustrates the errors of this ray tracer and that of the eikonal ray tracer. From it we can see, the proposed ray tracer has almost the same accuracy as the eikonal ray tracer which used a finer grid spacing. Figure 2 displays the travel time field of the first experiment. On our AMD 2.2 GHZ computer, it takes less than 2 seconds. Figure 3 shows the ray paths to the four model boundaries. The ray paths may cross each other and not enough curve. That is because of numerical modeling errors. Figure 4 shows a ray tracing of reflections in the same model and assuming there is a dipping reflector in the middle. Figure 5 shows a ray tracing of a multiple reflection. Only one multiple is shown although all 2

multiples are obtained simultaneously. Figure 6 shows a simple example of 3D ray tracing. Figure 7 and Figure 9 shows the ray tracing result in a complex salt dome model. The grid size is 150x645, and it takes about 20 seconds. Figure ?? shows the CPU time used for tracing all grid points in different size models.

Figure 1: Comparison of the traveltimes computed by the ray tracer and the eikonal finite difference method.

CONCLUSION In this report I proposed a flexible ray tracing method using Matlab that is versatile, accurate, and affortablly fast. In order to gain the same speed as the eikonal ray tracer , we may need to implement it in C or Fortran. The obvious advantage of this ray tracer is its capability to find ray paths and trace reflections and multiples. This ray tracer has already been successfully used in our tomography and migration algorithms for both synthetic and real data.

ACKNOWLEDGMENT I thank the members of the 2002 University of Utah Tomography and Modeling/Migration (UTAM) Consortium for the financial support of this research.

3

Figure 2: Travel time field in a linear gradient model.

Figure 3: Ray paths in a linear gradient model.

4

Figure 4: Reflection ray tracing in a dipping layer model with a velocity gradient.

Figure 5: Same as the previous figure except multiple ray tracing.

5

Figure 6: 3D ray tracing example.

Figure 7: Raytracing of traveltime field in a salt dome model.

6

Figure 8: Ray paths in the previous model.

Figure 9: CPU time vs model size.

7

REFERENCES Moser, T.j., 1991, Shortest path calculation of seismic rays, Geophysics, Vol 56, Issue 1, 59-67. Fischer, R., and Lees, J., 1993, Shortest path ray tracing with sparse graphs, Geophysics, Vol. 58, NO. 7, 987-996. R. T. Langan, 1985, Tracing of rays through heterogeneous media: An accurate and efficient procedure, Geophysics, Vol. 50, No. 9, 1456-1465. Sun, H., 1999, 2-D and 3-D Shooting Ray Tracing: UTAM 1999 Annual Report, 401-412.

8...