J. Cent. South Univ. Technol. (2007)01-0100-04
DOI: 10.1007/s11771-007-0020-5
Fast ray tracing method in 3-D structure and its proof of positive definiteness
GAO Er-gen(高尔根)1, Uk HAN2, TENG Ji-wen(滕吉文)3
(1. School of Earth and Space Sciences, University of Science and Technology of China, Hefei 230026, China;
2. Department of Environmental Sciences, Korea Military Academy, Seoul 139799, Korea;
3. Institute of Geophysics, Chinese Academy of Sciences, Beijing 100101, China)
Abstract: Based on Fermat’s principle, two-point ray tracing method was studied in three-dimensional structure. By means of first order Taylor’s incomplete series expansion (i.e. no expansion to the length of the ray), a symmetry block tridiagonal matrix equation set was deduced. Further, the positive definiteness of coefficient matrix was discussed, and the positive definiteness was accurately proved in a mathematical way. It assured that the algorithm was well-posed. Associated with iterative method, the solution to ray tracing can be got through step-by-step linearized iteration of the nonlinear problem. An algorithm of the whole path iterative ray tracing method in three-dimensional velocity structure was obtained. This method shows a clear and simple as well as explicit computation formula, which makes ray tracing computation easily applicable in practice. The correction vector is obtained through finding the solution to the positive definite block tridiagonal equation set, which ensures the method is robust convergence. This study offers a new kind of feasible and efficient ray tracing method for three dimensional seismic migration and tomography. Meanwhile, it also provides the prerequisite guarantee to design a fast algorithm.
Key words: Fermat’s principle; ray tracing; ray path; positivity; seismic migration; tomography
1 Introduction
Ray tracing method takes a very important role in fields such as seismology and seismic exploration data processing. It is one of the effective means to study the propagation of the seismic waves in situations of random velocity distribution and it has been widely used in amplitude versus offset (AVO) analysis, prestack depth migration, seismic velocity analysis based on models, tomography and earthquake locating. It can be applied to obtain, when the underground medium velocity has already known, the trajectory of ray path and travel time through the tomographic zone. The basic principle of the ray method had important applications long ago in the study of the property and behavior of light and it was comparatively late used in seismology. The early application of ray tracing method in seismology and seismic exploration was limited to kinematics. Until KARAL et al[1] applied the results of the electro- magnetics research to the area of elastic waves and the further developments made by ?ERVEN? et al[2], the theory of ray progression method was formed and ray method theories were greatly enriched. In past decades, many geophysicists were engaged in numerical simulated computations and method investigation concerning with ray tracing. A large number of data had been accumulated[1-18], which is significant in promoting the development of seismology and the improvement of the methods and theories in geophysical explorations.
In the aspect of the application of seismic waves such as earthquake locating and the studies of the subsurface velocity structure by means of seismic data tomography, geophysicists have to find solutions to the problem of two point ray tracing, which has made two points ray tracing occupy a very important position in ray tracing. However, it is a nonlinear problem when the two points are fixed such as source point (shot point) and the receiver point (detector point), and then, it is very difficult to be solved. FARRA et al[3-9] conducted the researches relating to the two point ray tracing and proposed their practical computation methods which have made contributions to the improvement of the effect of tomography and the quality of migration. However, no solution could settle all the problems, although, dealing with special medium distributions, each method may have its own advantage. In terms of special medium distribution, it can be convinced that relative ray tracing method is needed before the achievement of the highest efficiency of tracing computations. Therefore, it is necessary to bring forth a suitable computing method for ray tracing with a view to the actual medium distribution. In this study, the focus is to deal with a computation method for two point ray tracing in three-dimension structures—the whole path iterative ray tracing computation method and it tries to make further studies of the positive definiteness and to provide accurate mathematical proof. This method can be suitable to block model-building for velocity distribution and arbitrary interface connection between blocks.
2 Principle of method
For example, taken transmission waves as example, the principle is the same in terms of reflected waves. As shown in Fig.1, the v1, v2, …, and vn are layer velocities; f1(x, y), f2(x, y), …, and fn(x, y) are interface functions; S and R are shot point and receiver point, respectively, and they lie in (x0, y0, z0) and (xn+1, yn+1, zn+1). It is assumed that the lengths of RP1, P1P2, …, PnS are l1, l2, …,ln+1, respectively. The total travel time is given by the following equation:
(1)
where zi=fi(xi, yi), i=1, 2, …, n.
Fig.1 Diagram of whole path iterative ray tracing
The RP1P2…PnS is initial ray path, and RP′1P′2…P′nS is the ray path after one iteration.
According to the Fermat’s principle, we can get:
(2)
where i=1, 2, …, n.
Assumed that the initial ray path point (intersection between ray path and each subsurface) is (xi, yi, fi(xi, yi)), the shortest time ray path corresponds to (xi+Δxi, yi+Δyi, fi(xi+Δxi, yi+Δyi)), and replace xi，yi and zi in Eqn.(2) with xi+Δxi，yi+Δyi and fi(xi+Δxi, yi+Δyi), respectively and expand Eqn.(2) in first order Taylor’s incomplete series in the initial ray path point, i.e. carry out no expansion as to the length of the ray. Let: ai=x-xi-1，bi=yi-yi-1，ci= fi(xi, yi)-fi-1(xi-1, yi-1), θi=?fi(xi, yi)/?x, ψi=?fi(xi, yi)/?y, γi=1/(vili), the following result can be finally deduced from Eqn.(2) :
(3)
where Ai and Bi are both second-order matrixes; the main diagonal element of Ai is
=
, 
=
，the negative diagonal element is 
=
=
; the main diagonal element of Bi is
=
, 
= 
， the negative diagonal element is 
=
,
=
；the block vector of the right-end term is di=((ai+ciθi)γi-(ai+1+ci+1θi)γi+1, (bi+ciψi)γi-(bi+1+ci+1ψi)γi+1)T，the correction value of various ray path points is mi=(Δxi, Δyi)T, i=1, 2, …, n.
After obtaining the solutions to Δxi and Δyi through Eqn.(3), we use (xi+Δxi, yi+Δyi, fi(xi+Δxi, yi+Δyi)) in place of the original(xi, yi, zi). Going on repeat the above process until the sum of a certain norm of Δxi and Δyi on the whole path can satisfy the precision requirement of ray tracing.
The above result of deduction is obtained under the circumstance of transmission waves. In the event of reflected waves, it remains the same. All that needed is to give a serial number, according to progressive increases and from bottom upward following the reflection point, to the interface function and the velocity distribution. If in media there is lateral velocity variables or there exists the complex features of geological structure such as faults, nips, overlaps and phase transformation. Such method can be used as coordinate transformation: conducing mathematical description of geological models by means of the block models, treating each layer as an interface; for the layers with bigger obliquity, if necessary replacing z with x or z with y. Go on with the above deduction without changing the form of the formula.
3 Symmetrical positivity of coefficient matrix and its proof
In the above derivation of equations concerning the whole path iterative ray tracing method, as a result of adopting such methods as an one order Taylor’s incompletion series expansion, the coefficient matrix obtained is a block tridiagonal symmetric positivity matrix. It, therefore, has a well-posed property. In the process of using whole path iterative method to obtain the ray paths, the robust convergence can be still achieved, even though there are more times of iterations, when the initial paths are not in a fine state. Of course, in seismic prospecting, as the receiver points and shooting points are ranged according to a certain order, it will be effective to deal with the undesirable state of the initial path if the previous ray path is used as the next initial ray path. But as far as the research method is concerned, it is of great significance to study the robust computation method. The following is for the proof of the positivity of block tridiagonal symmetric matrix.
To prove that the symmetric coefficient matrix of Eqn.(3) is a positivity matrix, we may as well assume the coefficient matrix of Eqn.(3) to be 
，then, according to the property of matrix, we only prove that there is 
in each given nonzero vector 
, and mk=(αk, βk)T，k=1, 2, …, n.}. Therefore, we provide the following block matrix and expression for the block vectors:
, 
where k=2, 3, …, n.
The following result can be obtained by taking advantage of the rule for the multiplication of block matrix:
(4)
In the same way as in Eqn.(4), expand
again and divide the last term into two terms，that is:
(5)
With a view to the fact that each γi=1/(vili)>0, we can obtain:
(βiψi-βi-1ψi-1+αiθi-αi-1θi-1)2≥
(i=2, 3, …, n) (6)
Again because:
≥0 (i=1, 2, …, n) (7)
If substituting Eqns.(6) and (7) into Eqn.(5), we can obtain the following expression:
≥
Go on using the results of Eqns.(6) and (7), and conducing the computation in above ways, we can obtain:
≥
≥
≥…≥
≥0
To make the above equations sound, it’s necessary that each mk=0, k=1, 2, …, n, which is contradictory to the requirement that 
should be nonzero vector. Therefore, for the nonzero vector, there is:
As a result, the positivity of coefficient matrix is proved.
4 Discussion
In mathematics, the problem of ray tracing between two points is actually the boundary value problem. The numerical value solution to the boundary value problems often leads to the necessity of finding a solution to a nonlinear equation set, as is the case in Eqn.(2). The iterative method is currently applied to the numerical value solution to a nonlinear equation set and the conventional iterative method requires that the iterative initial value should be selected from certain neighbor- hood of the solution. Otherwise, the process of iteration may not converge to the true solution. However, we don’t know the solutions themselves to the nonlinear problems and at the same time it is very difficult to measure the neighborhood, and therefore it is difficult to determine the solution to the neighborhood. If the boundary value problems are transformed into the initial value problems, it is required, in terms of ray tracing between two points, that the angle of incidence of the given ray at the shooting point should be determined so that the ray can reach the receiver point. In three- dimension complex geological structures, the angle of incidence has to be tried many time, moreover, large amount of time has to be spent on calculation. Even so, it may sometimes occur that convergence can not be achieved.
5 Conclusions
1) The whole path iterative ray tracing method in three-dimensional velocity structure presented shows a clear and simple as well as explicit computation formula which makes ray tracing computation easily applicable in practice.
2) The correction vector is obtained through finding the solution to the positive definite block tridiagonal equation set, which on the one hand ensures the method is robust convergence and on the other hand provides mathematical support for the designing of fast ray tracing computation formula.
3) This study offers a new kind of feasible and efficient ray-tracing method for three dimensional seismic migration and tomography.
References
[1] KARAL F C, KELLER J B. Elastic wave propagation in homogeneous and inhomogeneous media[J]. J Accoust Soc Am, 1959, 31(6): 694-705.
[2] ?ERVEN? V, MOLOTKOV I A, P?EN??K I. Ray Method in Seismology[M]. Prague: Charles University Press, 1977.
[3] FARRA V. Bending method revisited: a Hamiltonian approach[J]. Geophys J Int, 1992, 109(1): 138-150.
[4] GAO Er-gen, XU Guo-ming, ZHAO Yi. Segmentally-iterative ray tracing method for any interface[J]. Oil Geophysical Prospecting, 1998, 33(1): 54-60.(in Chinese)
[5] GAO Er-gen, XU Guo-ming. A new kind of step by step iterative ray-tracing method[J]. Acta Geophysica Sinica, 1996, 39(S): 302-308.(in Chinese)
[6] GAO Er-gen, XU Guo-ming, Li Guang-pin, et al. A new total iterative ray-tracing method in random interface[J]. Acta Custica, 2002, 27(3): 282-287.(in Chinese)
[7] SAMBRIDGE M S, KENNETT B L N. Boundary value ray tracing in heterogeneous medium: a simple and versatile algorithm[J]. Geophys J Int, 1990, 101(1): 157-168.
[8] UM J, THURBER C H. A fast algorithm for two-point seismic ray tracing[J]. Bull Seis Soc Am, 1987, 77(3): 972-986.
[9] XU Guo-ming, WEI Shan, GAO Er-gen, et al. Block model-building and ray-tracing in 2-D complicated medium[J]. Oil Geophysical Prospecting, 2001, 36(2): 213-219.(in Chinese)
[10] CHIU S K L, KANASEWICH E R, PHADKE S. Three-dimensional determination of structure and velocity by seismic tomography[J]. Geophysics, 1986, 51(8): 1559-1571.
[11] FARRA V. Ray tracing in complex media[J]. J Appl Geophys, 1993, 30(1): 55-73.
[12] LIU Hong, MENG Fan-lin, LI You-ming. The interface grid method for seeking global minimum travel-time and the correspondent ray path[J]. Acta Geophysica Sinica, 1995, 38(6): 823-832.(in Chinese)
[13] VIDALE J E. Finite different calculation of travel-time in three dimensions[J]. Geophysics, 1990, 55(5): 521-526.
[14] MA Zheng-ming, LI Yan-da. Two-step ray tracing method[J]. Acta Geophysica Sinica, 1991, 34(4): 501-508. (in Chinese)
[15] YANG Chang-chun, LENG Chuan-bo, LI You-ming. Fast and accurate ray tracing in 3-D media[J]. Acta Geophysica Sinica, 1997, 40(3): 414-420. (in Chinese)
[16] YANG Wen-cai, LI You-ming, LI Shi-xiong, et al. Applied seismic Tomography[M]. Beijing: Geological Press, 1993. (in Chinese)
[17] ZHANG Lin-bin, YAO Zhen-xing, JI Chen. Finite difference calculation of seismic first-arrival travel-time[J]. Progress in Geophysics, 1996, 11(4): 47-52. (in Chinese)
[18] ZHOU Zhu-sheng, ZHANG Sai-min, CHEN Ling-jun. Seismic ray tracing calculation based on parabolic travel-time interpolation[J]. J Cent South Univ Technol, 2004, 11(2): 199-205.
(Edited by YANG You-ping)
Foundation item: Project(40674071) supported by the National Natural Science Foundation of China; project(KFAS2002-2003) supported by Korea Foundation for Advanced Studies
Received date: 2006-04-21; Accepted date: 2006-06-18
Corresponding author: GAO Er-gen, Associated professor, PhD; Tel:+86-551-3606871; E-mail: grg@ustc.edu.cn
