网格法射线追踪研究及应用
|
摘要
射线追踪方法作为一种快速有效的波场近似计算方法，对于地震波理论研究以及地震波反演及偏移成像等过程具有重要意义。其理论基础是，在高频近似条件下，地震波场的主要能量沿射线轨迹传播。传统的射线追踪方法包括初值问题的试射法和边值问题的弯曲法。试射法根据由源出发的一束射线到达接收点的情况对射线出射角及密度进行调整，最后由最靠近接收点的两条射线走时内插求出接收点处走时；弯曲法则是从源与接收点之间的一条假想初始路径开始，根据最小走时准则对路径进行扰动，从而求出接收点处的走时及射线路径。
近年来，随着Kirchhoff积分叠前深度偏移在解决复杂构造成像中获得一系列成功，作为其算法基础之一的射线追踪方法也得到了很大的促进和发展，出现了大量不同于传统方法的新型算法。这些方法的主要特点在于不再局限于地震波的射线路径描述,而是直接从Huygens原理或者Fermat原理出发，采用等价波前描述地震波场的特征。在导师的指导下，我在研究传统的射线追踪方法以及新发展出来的这些新的射线追踪方法后，根据相关文献里的某种思路，实现了一种网格法射线追踪，并在不同的模型上对其进行了实践。
传统的射线追踪都是从射线方程出发，对介质进行了离散化之后运用射线方程求出走时和射线，但是传统的射线追踪方法存在着很多不足之处，主要问题在于:
① 难于处理介质中较强的速度变化；
② 难于求出多值走时中的全局最小走时；
③ 计算效率较低；
④ 阴影区内射线覆盖密度不足。
然而按照目前的观点，仅考虑地震波所有走时中最小走时无疑具有很大局限性，即使从射线偏移的角度来看，要获得较好的成像效果，只考虑地震波走时中的最小走时也是远远不够的。因而最近几年，关于射线追踪方法的研究主要集中在多值走时计算方面,研究进展主要体现在:
传统的试射法及弯曲法的基础上的改进,如各类波前重建方法，除多值走时外,还较好地解决了计算效率及阴影区覆盖不足的问题；
对最小走时算法的改进,使之可适应多值走时计算,如慢度匹配法，可认为是最短路径方法的推广；
③ 传统方法与最小走时算法的结合,如ＨＷＴ方法，则是通过波前传播计算射线路径。
本文阐述了两种比较常用的射线追踪方法——有限差分法和有序波前重建法。它们都解决了传统射线追踪方法的不足。有限差分法提出的计算旅行时的程函方程有限差分法，对于速度差较大的模型也能得到正确结果，但对于速度结构较复杂的模型，仍像Vidale算法一样找不到最小旅行时，因此仍需寻求一种更有效的搜索方案使有限差分法能真正实现波前面追踪。有序波前重建方法有如下特点：
（1）用实际波前扩展代替扩展方阵方式进行射线追踪计算，较好地把握了射线走向，提高了计算精度；
(2)网格节点最小走时计算采用搜索次级源的计算策略，充分考虑了全波场信息。走时计算的同时记录射线，实现了基于两步法思想的同步计算，降低了计算量；
(3)算法原理清晰，易于实现，波前走时计算采用双曲线近似计算，近似程度高，实用性强。
可是我们必须看到有序波前重建法的不足之处：由于离散方法本身使得走时曲线不可能足够光滑，故也不适于叠前速度分析和叠前速度偏移。为了
解决这个难题，Popovici和Sethian又引进一种三维快速推进（FM）技术，它所解的程函方程是沿着扩展波前而不是扩展方阵。这种FM技术保持了原有的有限差分法的效率。不幸的是，这个进步的有限差分方法的精确性还是难以满足近地表的实际应用，比如在初至层析成像中。
鉴于上面所说的几种射线追踪方法的特点，我们引入一种非常直观的网格法射线追踪技术（GRT）。它结合了有限差分法和波前重建法的优点，并且更加结合了近年来一种非常先进的FM技术把FM法和波前构造方法的优点结合在一起。它非常灵活，为精确的走时和射线计算在近地表应用提供了有效的工具。网格法射线追踪具有以下特点：
（1）网格法结合了有限差分法，在计算走时和射线路径时都很精确而且效率也很高，同时通过结合FM方法，避免了有限差分法无法建立把射线转换成初至层析成像所需要的模型的缺陷。
（2） 网格法运用了快速波前追踪算法，在射线计算方面效率很高。同时避免了影子区问题，可以在任何介质中建立弯曲曲线的模型。
（3）网格法是以网格射线追踪和局部波前构造为基础的，所以它可以在不需要大量的进行球状波前构造和走时的插值就能精确的对波动向量进行描述。
通过网格法在不同模型中的实践，检验出它在不同模型中的精确性和高效率性，其中重点我们把网格法与比较先进的FM法进行了比较，发现无论在二维还是在三维的模型中，网格法在精确性和效率上都优于FM方法。
As a fast and effective approximative calculate method of wave field,raytracing has a very important meaning to the research of the seism wave theory and the seism wave’s invertion and excursion imaging process. The theory foundation of the raytracing is that under the conditions of high-frequency approximate the major energy of the seism wave transmits along the radial trace. The traditional raytracing methods comprise the try-shooting method of the initial value problem and the bending method of the border value problem. According to the circs of a bunch of rays set up from the source to the receive point, the try-shooting method gets the traveltime at the receive point by the interpolation between two traveltime of rays close to the receive point. The bending method begins with a imaginary radial road from the source to the receive point, then disturb the ray road base on the minimum traveltime rules and at last get the traveltime and the raytrace.
In recent years, the Kirchhoff’s intergral fore-stack depth migrations get much useful in the complex formation imaging. As one of its algorithm basises, the raytracing metods develop very fast. Large number of new algorithms different from the traditional raytracing methods have appeared. The main features of these methods are that they aren’t limited by the descriptions of the travetime and raytracing of the seism wave. They set out directly from the principle of Huygens and the Fermat theory and then descript the features of the seism wave fields by adopting the equivalent fore-waves. Under my supervisor’s guiding, I have studied traditional raytracing methods and their new developments and then get a grid raytracing method according to the ideas in correlate documents. At last I practise the method on difference models.
All traditional raytracing methods set out from the ray-equations and after divergence to the medium they get the traveltime and raytace from the ray-equations. But they have many deficiencies and main problem is that:
Hard to handle the stronger change of velocity of the medium.
Hard to work out the global minimum traveltime of hypervalue traveltimes.
Calculation efficiency is lower.
The shortage of the rays overlay density in the shadow areas.
However, it is very limited to consider the minimum traveltime in the seism wave traveltimes only. Even if considering the ray migration, it is far not enough. As a result, in the recent years, it concentrates primarily the calculations of the multivalue traveltimes. The development of the research is that:
The improvement on the foundation between traditional try-shooting method and bending method, eg. Various wave front
reconstruction methods. In addition to multivalue traveltimes, they have resolved the calculation efficiency and the shadow areas problems.
The improvement of the minimum traveltime algorithms, for example, slowness matching method. It is a extend of the shortest path methods.
The combinative method between the traditional methods and the shortest traveltime method, such as HWT method, it calculate ray path by the wave front transmitting.
The text expands two raytracing methods in common use, that is the finite difference method and the ordered wave front reconstruction method. They all have resolved the defects of the traditional raytracing methods. The finite difference method brings up eikonal equation finite difference method for the traveltime calculation. It can get exact result even in greater speed model. But for complicated velocity model in flat-out construction, it can’t find the minimum traveltime similar to the Vidale algorithm. Therefore it still needs to find a kind of efficient search scheme to let the finite difference method realize the wave front tracing really. The ordered wave front reconstruction method has features as follows:
a. Put up the raytracing calculation by the actual wave front expand instead of the expand square matrix. Better to hold the raytrace trend. Increase the calculation precision.
b. The calculations of the minim
近年来，随着Kirchhoff积分叠前深度偏移在解决复杂构造成像中获得一系列成功，作为其算法基础之一的射线追踪方法也得到了很大的促进和发展，出现了大量不同于传统方法的新型算法。这些方法的主要特点在于不再局限于地震波的射线路径描述,而是直接从Huygens原理或者Fermat原理出发，采用等价波前描述地震波场的特征。在导师的指导下，我在研究传统的射线追踪方法以及新发展出来的这些新的射线追踪方法后，根据相关文献里的某种思路，实现了一种网格法射线追踪，并在不同的模型上对其进行了实践。
传统的射线追踪都是从射线方程出发，对介质进行了离散化之后运用射线方程求出走时和射线，但是传统的射线追踪方法存在着很多不足之处，主要问题在于:
① 难于处理介质中较强的速度变化；
② 难于求出多值走时中的全局最小走时；
③ 计算效率较低；
④ 阴影区内射线覆盖密度不足。
然而按照目前的观点，仅考虑地震波所有走时中最小走时无疑具有很大局限性，即使从射线偏移的角度来看，要获得较好的成像效果，只考虑地震波走时中的最小走时也是远远不够的。因而最近几年，关于射线追踪方法的研究主要集中在多值走时计算方面,研究进展主要体现在:
传统的试射法及弯曲法的基础上的改进,如各类波前重建方法，除多值走时外,还较好地解决了计算效率及阴影区覆盖不足的问题；
对最小走时算法的改进,使之可适应多值走时计算,如慢度匹配法，可认为是最短路径方法的推广；
③ 传统方法与最小走时算法的结合,如ＨＷＴ方法，则是通过波前传播计算射线路径。
本文阐述了两种比较常用的射线追踪方法——有限差分法和有序波前重建法。它们都解决了传统射线追踪方法的不足。有限差分法提出的计算旅行时的程函方程有限差分法，对于速度差较大的模型也能得到正确结果，但对于速度结构较复杂的模型，仍像Vidale算法一样找不到最小旅行时，因此仍需寻求一种更有效的搜索方案使有限差分法能真正实现波前面追踪。有序波前重建方法有如下特点：
（1）用实际波前扩展代替扩展方阵方式进行射线追踪计算，较好地把握了射线走向，提高了计算精度；
(2)网格节点最小走时计算采用搜索次级源的计算策略，充分考虑了全波场信息。走时计算的同时记录射线，实现了基于两步法思想的同步计算，降低了计算量；
(3)算法原理清晰，易于实现，波前走时计算采用双曲线近似计算，近似程度高，实用性强。
可是我们必须看到有序波前重建法的不足之处：由于离散方法本身使得走时曲线不可能足够光滑，故也不适于叠前速度分析和叠前速度偏移。为了
解决这个难题，Popovici和Sethian又引进一种三维快速推进（FM）技术，它所解的程函方程是沿着扩展波前而不是扩展方阵。这种FM技术保持了原有的有限差分法的效率。不幸的是，这个进步的有限差分方法的精确性还是难以满足近地表的实际应用，比如在初至层析成像中。
鉴于上面所说的几种射线追踪方法的特点，我们引入一种非常直观的网格法射线追踪技术（GRT）。它结合了有限差分法和波前重建法的优点，并且更加结合了近年来一种非常先进的FM技术把FM法和波前构造方法的优点结合在一起。它非常灵活，为精确的走时和射线计算在近地表应用提供了有效的工具。网格法射线追踪具有以下特点：
（1）网格法结合了有限差分法，在计算走时和射线路径时都很精确而且效率也很高，同时通过结合FM方法，避免了有限差分法无法建立把射线转换成初至层析成像所需要的模型的缺陷。
（2） 网格法运用了快速波前追踪算法，在射线计算方面效率很高。同时避免了影子区问题，可以在任何介质中建立弯曲曲线的模型。
（3）网格法是以网格射线追踪和局部波前构造为基础的，所以它可以在不需要大量的进行球状波前构造和走时的插值就能精确的对波动向量进行描述。
通过网格法在不同模型中的实践，检验出它在不同模型中的精确性和高效率性，其中重点我们把网格法与比较先进的FM法进行了比较，发现无论在二维还是在三维的模型中，网格法在精确性和效率上都优于FM方法。
As a fast and effective approximative calculate method of wave field,raytracing has a very important meaning to the research of the seism wave theory and the seism wave’s invertion and excursion imaging process. The theory foundation of the raytracing is that under the conditions of high-frequency approximate the major energy of the seism wave transmits along the radial trace. The traditional raytracing methods comprise the try-shooting method of the initial value problem and the bending method of the border value problem. According to the circs of a bunch of rays set up from the source to the receive point, the try-shooting method gets the traveltime at the receive point by the interpolation between two traveltime of rays close to the receive point. The bending method begins with a imaginary radial road from the source to the receive point, then disturb the ray road base on the minimum traveltime rules and at last get the traveltime and the raytrace.
In recent years, the Kirchhoff’s intergral fore-stack depth migrations get much useful in the complex formation imaging. As one of its algorithm basises, the raytracing metods develop very fast. Large number of new algorithms different from the traditional raytracing methods have appeared. The main features of these methods are that they aren’t limited by the descriptions of the travetime and raytracing of the seism wave. They set out directly from the principle of Huygens and the Fermat theory and then descript the features of the seism wave fields by adopting the equivalent fore-waves. Under my supervisor’s guiding, I have studied traditional raytracing methods and their new developments and then get a grid raytracing method according to the ideas in correlate documents. At last I practise the method on difference models.
All traditional raytracing methods set out from the ray-equations and after divergence to the medium they get the traveltime and raytace from the ray-equations. But they have many deficiencies and main problem is that:
Hard to handle the stronger change of velocity of the medium.
Hard to work out the global minimum traveltime of hypervalue traveltimes.
Calculation efficiency is lower.
The shortage of the rays overlay density in the shadow areas.
However, it is very limited to consider the minimum traveltime in the seism wave traveltimes only. Even if considering the ray migration, it is far not enough. As a result, in the recent years, it concentrates primarily the calculations of the multivalue traveltimes. The development of the research is that:
The improvement on the foundation between traditional try-shooting method and bending method, eg. Various wave front
reconstruction methods. In addition to multivalue traveltimes, they have resolved the calculation efficiency and the shadow areas problems.
The improvement of the minimum traveltime algorithms, for example, slowness matching method. It is a extend of the shortest path methods.
The combinative method between the traditional methods and the shortest traveltime method, such as HWT method, it calculate ray path by the wave front transmitting.
The text expands two raytracing methods in common use, that is the finite difference method and the ordered wave front reconstruction method. They all have resolved the defects of the traditional raytracing methods. The finite difference method brings up eikonal equation finite difference method for the traveltime calculation. It can get exact result even in greater speed model. But for complicated velocity model in flat-out construction, it can’t find the minimum traveltime similar to the Vidale algorithm. Therefore it still needs to find a kind of efficient search scheme to let the finite difference method realize the wave front tracing really. The ordered wave front reconstruction method has features as follows:
a. Put up the raytracing calculation by the actual wave front expand instead of the expand square matrix. Better to hold the raytrace trend. Increase the calculation precision.
b. The calculations of the minim
引文
[1] 杨文采. 地球物理反演和地震层析成象 北京：地质出版社，1989.
[2] 杨宝俊等. 勘探地震学导论. 吉林科学技术出版社，1992.
[3] 杨文采等. 应用层析成像. 北京：地质出版社，1993.
[4] G.Nolet（荷兰）.地震层析成象及应用. 王椿镛、吴宁远等译. 北京：学术书刊出版社，1989.
[5] 何樵登、熊维纲. 应用地球物理教程——地震勘探. 北京：地质出版社，1991
[6] 赫尔曼. 层析成象和反演问题的基本方法. 北京：石油工业出版社，1997
[7] 常旭、刘伊克. 地震正反演与成像[Ｍ]. 北京：华文出版社，2001.
[8] 赵文敏. 有限差分初步. 台湾：台湾徐氏基金，1977
[9] 赵成斌等．地震波走时和射线的有限差分计算．华北地震科学，1995， 13（4），41～46
[10] 王新民、董小刚. 工程中德数值方法. 长春：吉林人民出版社，2001
[11] 董敏煜. 地震勘探信号分析. 东营：石油大学出版社，1989
[12] 王华忠等. 有限差分法地震波走时计算. 同济大学学报，1997，25（3），318～321
[13] 张霖斌等. 地震反射波波前法射线追踪. 石油地球物理勘探，1995，30（1），50～56
[14] 赵连锋等. 有序波前重建法的射线追踪. 地球物理学报，2003，46（3），415～420
[15] 张钊等. 射线追踪方法的发展现状. 地球物理学进展，2000，15（1），36～45
[16] 刘洪等. 计算最小走时和射线路径的界曲网全局方法. 地球物理学报，1995，38，821～831.
[17] 陈裕明等. 最小时间两点射线追踪. 石油地球物理勘探，2002，37（5），455～459
[18] 高尔根等. 一种任意界面的逐段迭代射线追踪方法. 石油地球物理勘探，1998，33（1），54～60.
[19] 朱介寿等. 任意介质中的射线追踪方法. 地震学中的计算方法,地震出版社，1998，603～630.
[20] 黄联捷等. 用于图像重建的波前法射线追踪. 地球物理学报，1992，35（2），223～
233.
[21] 马争鸣等. 二步法射线追踪. 地球物理学报，1991，34（4），501～508.
[22] 张建中等. 最短路径射线追踪方法及其改进. 地球物理学进展，2003，18（1），146～150.
[23] 谢海兵等. 一种新程函方程差分法层析成象. 石油地球物理勘探，1997，32（5），623～628.
[24] 杨长春等. 射线追踪的微变网格法. 地球物理学报，1996，39（1），97～102.
[25] 鲁传恒、顾贤明. 井间地震射线层析成象. CT理论与应用研究，1994，3（3），18～25.
[26] 赵改善等. 基于旅行时线性插值的地震射线追踪算法. 石油物探，1998，37（2），14～24.
[27] 张霖斌等. 有限差分法射线追踪. 石油地球物理勘探，1993，29（6），673～677.
[28] 刘殿秘、韩立国. 一种改进的网络射线追踪方法. 吉林大学学报，2004，34，77～80.
[29] Cao S and Greenhalgh S A. Finite-difference solution of the eikonal equation using an efficient, firstarrival wavefront tracking scheme. Geophysics, 1994, 59, 632～643.
[30] Jackson M J and Tweeton D R. Migratom-geophysical tomography using wavefront migration andfuzzy constraints, Bureau of Mines Report of Investigation9497, United States Department of the Interior.
[31] Vidale, J. E., Finite-difference calculation of traveltimes, Bull. Seism. Soc. AM. 78, 2062-2076.
[32] Qin F, Luo Y, Olsen K B, Schuster G T. Finite-difference solution of the eikonal equation along expanding wavefronts. Geophysics, 1992, 57， 478～487.
[33] Popovici A M，Sethian J A. Three dimensional traveltime computation using the fast marching method. 67th Annual Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1997, 1778～1781.
[34] Vidale J E. Finite-difference calculation of traveltimes in three dimensions. Geophysics, 1990, 55: 521～526.
[35] Zhu T. and Brown L. Two-dimensional velocity inversion and synthetic seismogram computation. Geophysics, 1987, 52, 37～50.
[36] Mose r T. Shortest path calculation of seismic rays. Geophysics，1991，
56，59～67.
[37] William A and Schneider Jr. Robust a md efficient upwind finite一difference traveltime calculations in three dimension. Geophysics，1995，60，1108～1117.
[38] Trier J V. Up-wind finite difference calculation of seismic traveltimes. Presented at 60th SEG Annual.
[39] Klimes L and Kvasnicha M. 3-D network ray tracing[ J]，Ceophys.J.Int.，1994，116 ，726～738.
[40] U m J，Thurber C. H. A fast algorithm for two-point seism is ray tracing.Bull. Seism.Soc.am，1987，77，972～986.
[41] Podvin P and Lecomte I . Finite-difference computation of traveltimes in
very contrasted velocity models : A massively parallel approach and its
associated tools. Ceophys.J.Int.,1991，105，271～284
[2] 杨宝俊等. 勘探地震学导论. 吉林科学技术出版社，1992.
[3] 杨文采等. 应用层析成像. 北京：地质出版社，1993.
[4] G.Nolet（荷兰）.地震层析成象及应用. 王椿镛、吴宁远等译. 北京：学术书刊出版社，1989.
[5] 何樵登、熊维纲. 应用地球物理教程——地震勘探. 北京：地质出版社，1991
[6] 赫尔曼. 层析成象和反演问题的基本方法. 北京：石油工业出版社，1997
[7] 常旭、刘伊克. 地震正反演与成像[Ｍ]. 北京：华文出版社，2001.
[8] 赵文敏. 有限差分初步. 台湾：台湾徐氏基金，1977
[9] 赵成斌等．地震波走时和射线的有限差分计算．华北地震科学，1995， 13（4），41～46
[10] 王新民、董小刚. 工程中德数值方法. 长春：吉林人民出版社，2001
[11] 董敏煜. 地震勘探信号分析. 东营：石油大学出版社，1989
[12] 王华忠等. 有限差分法地震波走时计算. 同济大学学报，1997，25（3），318～321
[13] 张霖斌等. 地震反射波波前法射线追踪. 石油地球物理勘探，1995，30（1），50～56
[14] 赵连锋等. 有序波前重建法的射线追踪. 地球物理学报，2003，46（3），415～420
[15] 张钊等. 射线追踪方法的发展现状. 地球物理学进展，2000，15（1），36～45
[16] 刘洪等. 计算最小走时和射线路径的界曲网全局方法. 地球物理学报，1995，38，821～831.
[17] 陈裕明等. 最小时间两点射线追踪. 石油地球物理勘探，2002，37（5），455～459
[18] 高尔根等. 一种任意界面的逐段迭代射线追踪方法. 石油地球物理勘探，1998，33（1），54～60.
[19] 朱介寿等. 任意介质中的射线追踪方法. 地震学中的计算方法,地震出版社，1998，603～630.
[20] 黄联捷等. 用于图像重建的波前法射线追踪. 地球物理学报，1992，35（2），223～
233.
[21] 马争鸣等. 二步法射线追踪. 地球物理学报，1991，34（4），501～508.
[22] 张建中等. 最短路径射线追踪方法及其改进. 地球物理学进展，2003，18（1），146～150.
[23] 谢海兵等. 一种新程函方程差分法层析成象. 石油地球物理勘探，1997，32（5），623～628.
[24] 杨长春等. 射线追踪的微变网格法. 地球物理学报，1996，39（1），97～102.
[25] 鲁传恒、顾贤明. 井间地震射线层析成象. CT理论与应用研究，1994，3（3），18～25.
[26] 赵改善等. 基于旅行时线性插值的地震射线追踪算法. 石油物探，1998，37（2），14～24.
[27] 张霖斌等. 有限差分法射线追踪. 石油地球物理勘探，1993，29（6），673～677.
[28] 刘殿秘、韩立国. 一种改进的网络射线追踪方法. 吉林大学学报，2004，34，77～80.
[29] Cao S and Greenhalgh S A. Finite-difference solution of the eikonal equation using an efficient, firstarrival wavefront tracking scheme. Geophysics, 1994, 59, 632～643.
[30] Jackson M J and Tweeton D R. Migratom-geophysical tomography using wavefront migration andfuzzy constraints, Bureau of Mines Report of Investigation9497, United States Department of the Interior.
[31] Vidale, J. E., Finite-difference calculation of traveltimes, Bull. Seism. Soc. AM. 78, 2062-2076.
[32] Qin F, Luo Y, Olsen K B, Schuster G T. Finite-difference solution of the eikonal equation along expanding wavefronts. Geophysics, 1992, 57， 478～487.
[33] Popovici A M，Sethian J A. Three dimensional traveltime computation using the fast marching method. 67th Annual Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1997, 1778～1781.
[34] Vidale J E. Finite-difference calculation of traveltimes in three dimensions. Geophysics, 1990, 55: 521～526.
[35] Zhu T. and Brown L. Two-dimensional velocity inversion and synthetic seismogram computation. Geophysics, 1987, 52, 37～50.
[36] Mose r T. Shortest path calculation of seismic rays. Geophysics，1991，
56，59～67.
[37] William A and Schneider Jr. Robust a md efficient upwind finite一difference traveltime calculations in three dimension. Geophysics，1995，60，1108～1117.
[38] Trier J V. Up-wind finite difference calculation of seismic traveltimes. Presented at 60th SEG Annual.
[39] Klimes L and Kvasnicha M. 3-D network ray tracing[ J]，Ceophys.J.Int.，1994，116 ，726～738.
[40] U m J，Thurber C. H. A fast algorithm for two-point seism is ray tracing.Bull. Seism.Soc.am，1987，77，972～986.
[41] Podvin P and Lecomte I . Finite-difference computation of traveltimes in
very contrasted velocity models : A massively parallel approach and its
associated tools. Ceophys.J.Int.,1991，105，271～284