地球物理学进展  2014, Vol. 29 Issue (1): 255-260   PDF    
引用本文
秦宁, 李振春, 桑运云, 张凯. 初至波走时层析速度建模方法研究[J]. 地球物理学进展, 2014, 29(1): 255-260, doi: 10.6038/pg20140136  
QIN Ning, LI Zhen-chun, SANG Yun-yun, ZHANG Kai. The research of travel time tomographic velocity modeling method on first break. Progress in Geophysics, 2014, 29(1): 255-260, doi: 10.6038/pg20140136   
初至波走时层析速度建模方法研究
秦宁1, 李振春2, 桑运云3, 张凯2    
1. 中石化胜利油田物探研究院和中石化胜利油田博士后科研工作站, 东营 257002;
2. 中国石油大学(华东)地球科学与技术学院, 青岛 266580;
3. 中石油东方地球物理公司研究院华北分院, 任丘 062550
收稿日期: 2013-9-12 修回日期: 2013-12-11 .
基金项目：国家自然科学基金（41204086）、中国石油大学自主创新科研计划资助项目（13CX02041A）、教育部博士点基金项目（20110133120001）、国家863课题（2011AA060301）和国家科技重大专项（2011ZX05006-002）共同资助.
作者简介：秦宁，女，1985年生，博士，中石化胜利油田在站博士后，主要从事叠前深度偏移、层析建模和波形反演等方面的研究.（E-mail：geoqin@163.com）
摘要：层析反演的核心是敏感核函数的求取和旅行时反演方程的建立.在射线理论下， 射线追踪技术是计算敏感核函数和旅行时的主要方法， 因此初至层析的实现需要有一种可以准确追踪初至波路径的射线追踪方法.本文推导了抛物旅行时插值射线追踪公式， 提出了基于抛物旅行时插值的最短路径射线追踪方法.理论模型和实际资料处理结果表明， 利用该方法实现初至层析速度建模可行.
关键词： 射线追踪     抛物旅行时插值     初至波     走时层析     速度建模    
The research of travel time tomographic velocity modeling method on first break
QIN Ning1, LI Zhen-chun2, SANG Yun-yun3, ZHANG Kai2     
1. Geophysical Research Institute of Sinopec Shengli Oilfield, Postdoctoral Scientific Research Workstation of Sinopec Shengli Oilfield, Dongying 257002, China;
2. School of Geosciences, China University of Petroleum, Qingdao 266580, China;
3. Huabei Sub-institute of Research Institute of BGP, CNPC, Renqiu 062550, China
Abstract: The kernel of tomographic inversion is the solution of sensitive kernel function and the establishing set of travel time equations. Based on the ray theory, ray tracing is the main method to calculate the sensitive kernel function and travel time. Therefore the way to trace first break travel path exactly is needed for first break tomography. In this paper the formula of parabolic travel time interpolation is derived based on the previous research. This paper puts forward the dynamic shortest path ray tracing algorithm using parabolic travel time interpolation. The theoretical model and the real data prove the feasibility of the method, which is used to realize the velocity modeling based on the first break tomography.
Key words: ray tracing     parabolic travel time interpolation     first break     travel time tomography     velocity modeling    

0 引 言

沙漠地区表面厚度从50 m到200 m变化不等的低降速带，地表岩层纵横向岩性多变，速度差异也较大，地面接收到的地震波动力学和运动学特征随测线方向激发条件和接收条件的改变而改变.一方面，表层对地震波吸收严重，导致沙漠地区的地震资料信噪比低(程玖兵，2005，2009；陈湛文和张关泉，2004；董春晖，2008，董春晖和张剑锋，2010)；另一方面，低速层的厚度和速度变化较大，静校正问题非常严重(冯波和王华忠，2011；卢宝坤等，2012；黄中玉等，2009).沙漠地区地震勘探，静校正量的准确求取显得尤为重要，表层速度问题亟待解决(李博等，2009；刘守伟等，2008；刘璐等，2011).

初至走时层析可以很好的适应介质的横向或纵向变速(Nakanishi and Yamaguchi，1986；秦宁等，2012)，过程中利用直达波、折射波、回折波和透射波等不同类型的初至信息，是近地表速度建模的常用方法(秦宁等，2012；桑运云等，2013；Stork and Clayton，1991；孙福利等，2011)，而射线追踪是走时层析的核心技术(王昌龙等，2007；王华忠等，2010；王维红等，2011；张浩和张剑锋，2012).1986年Nakanishi和Yamaguchi(1986)首次把Dijstra网络最短路径算法应用于地震最小走时射线追踪(桑运云等，2013)；1991年Stork(1991)的最短路径射线追踪方法，把具有最小旅行时的网格节点的连线作为射线路径(桑运云等，2013).张建中等(2004年)克服了最短路径射线追踪方法中路径只能从一个节点到达另一个节点的限制，提出了动态网络最短路径射线追踪，并应用旅行时线性插值法(LTI)和Fermat原理进行插值，大大提高了射线追踪的精度.

本文提出的基于抛物旅行时插值的最短路径射线追踪(PTISPR)方法(以下简称PTISPR法)克服了最短路径方法中射线“之”字路径的缺陷，抛物插值方法提高了射线追踪的精度.PTISPR法既可以得到精确的初至走时又能准确计算层析核函数，同时节点设置方式使得此方法拥有变网格的效果，可以在大网格下实现层析速度建模.理论模型和实际资料试算证明了方法的可行性和实用性.

1 基于抛物插值的最短路径射线追踪方法

PTISPR法分两步进行：首先从震源开始正向波前扩展计算整个速度场的节点走时，然后从接收点反向追踪射线路径.其核心思想是对某一计算节点首先利用Dijstra算法计算旅行时，然后在计算过程中搜索满足抛物插值的节点并进行插值，将两部分得到的最小旅行时作为该计算节点的最终时间.射线追踪过程利用抛物旅行时插值从接收点网格开始一直搜索到炮点网格为止.

1.1 局部最小旅行时

波前扩展的过程是由已知旅行时的点(子震源)求取与其相连的节点的旅行时，计算常常在同一单元内，得到的最小旅行时是局部最小旅行时.在某一单元内，若已知节点i的走时，求与其相连接的节点j的旅行时公式为

其中dij为节点i，j的欧式距离，s是单元慢度.这种计算仅仅是点对点的计算，这也是最短路径射线追踪误差的来源之一.为克服这种缺陷，我们提出了抛物旅行时插值射线追踪.如图 1所示，在慢度为s的矩形单元内，若A,B,C为同一边界上的已知走时的三个节点，其坐标分别为(xa,za)，(xb,zb)，(xc,zc)，旅行时分别为TA，TB，TC.则他们就构成了节点D的插值段，现在求通过该边界上R点(xr,zr)到达D点(xd,zd)的最小旅行时TD和射线路径.
图 1 射线穿过网格边界及在一个网格单元内插值段的几何关系 Fig. 1 The diagram of rays passing through grids boundary and the geometry of rays interpolation in one grid

如图 1所示，当za=zr=zb=zc时：利用抛物插值方法得到

其中，AB=xa－xb，AC=xa－xc，BC=xb－xc，这样求解的D点的时间是整个网格内的局部最小走时.

1.2 全局最小旅行时

全局最小旅行时是地震波从震源点传播到任一节点时的最小旅行时.给定震源节点的旅行时，从震源开始，用Dijstra方法逐点计算从震源到达任一网格节点的最小旅行时.Dijstra算法的思想是，把全部节点N划分为四个子集.

P：已经是全局最小旅行时的节点，也就是做过子震源的节点子集；Q:已获得至少一次局部最小旅行时的节点，但没做过子震源的节点子集；R:在N中去掉P、Q后剩余节点子集；FS(i)：与震源(或者子震源)i直接相连接的节点子集，每一步在波前扩展时只计算这种节点的走时.

具体实现时：

1)初始化：P=Φ，t(s)=0，s为震源，Φ代表空集；Q=N，t(i)=∞，i∈N .

2)选择：在Q中选择最小旅行时的节点i，i∈Q .

3)更替：计算从i点传到j点的走时，同时由j所在网格单元内的走时插值段计算j节点的走时，把计算的j节点的局部最小走时与原值比较，若比原值小则用此值代替原值，否则原值保持不变.

4)判断迭代：如果P=N或Q=Φ迭代停止，否则转向2).

1.3 射线路径的确定

最短路径射线追踪在扩展波前的过程中，依次记下每一个节点的父节点(即上一级震源节点)，射线路径的确定是从接收点依次向前寻找其父节点，当追踪到震源点时结束，射线路径就是这些节点的连线.由于角度离散和空间离散的限制，射线路径往往呈“之”字形.为克服这种缺陷，从接收点开始通过互易原理，利用节点的全局走时和PTI方法，依次搜寻上级震源的位置，直到追踪到炮点网格，如图 2所示.

图 2 射线追踪路径示意图 Fig. 2 Paths of ray tracing

2 初至波走时层析理论

地震波的旅行时是介质慢度函数沿射线路径的线积分(Zhang，2013)后表示为

其中射线追踪得到的A是和慢度场 s(x,z)相关的射线路径构成的矩阵， ΔS是求解的慢度更新列向量，ΔT为实际资料或理论模型初至时间和正演计算的初至旅行时的残差列向量.

慢度更新量ΔS通过求解上面方程得到，对慢度场 s(x,z)进行更新，然后在得到的新的慢度场上，进行射线追踪建立新的方程，求解新的慢度修正量，迭代处理，直到满足初至层析反演收敛分析的条件，进而得到我们所需要的层析建模速度场，图3是初至走时层析速度建模流程图.

图3 初至走时层析速度建模流程图 Fig. 3 The flow of travel time tomographic velocity modeling on first break
3 模型试算和实际资料处理
3.1 模型试算

图4a为六层起伏模型，速度网格数101*301，横向采样间隔15 m，深度采样间隔为5 m，由浅到深各层速度分别为：500 m/s，800 m/s，1200m/s，1500m/s，2000m/s，3000m/s，低降速带深度范围500 m内，为起伏横向变速地层，高速层顶界面为水平地层，高速层内有两起伏界面和一水平界面.地面地震观测：炮点和检波点都位于起伏地表上，第一炮横向位置为0.0 m，炮间隔为250 m，共7炮，每炮101道接收，第一道横向位置0.0 m，道间隔15.0 m.图5a真实速度对应的等时线和射线路径图可知，初至波路径穿越了低降速带，因此可以通过初至走时层析进行反演.选用的初始速度场的高速层与真实速度场近似，低速层为梯度变化模型，利用加入正则化的LSQR分解方法求解层析方程组(董春晖和张剑锋, 2010)，经过16次迭代后平均走时残差从114.71 ms变为-0.49 ms，得到的层析反演结果如图4b所示.层析反演基本反映了低降速带的变化趋势，真实速度(图4c所示)和层析速度(图4d所示)对应的等时线和射线路径图基本一致.图4e分别对应速度场深度245 m处真实速度、层析速度和初始速度随深度变化曲线，层析速度更接近真实速度，低降速带反演结果在误差范围内.模型试算表明基于PTISPR法的初至走时层析可以应用于复杂近地表速度建模.

图4 (a)真实速度场;(b)层析反演速度场;(c)真实速度对应的射线路径和等时线图;(d)层析建模速度对应的射线路径和等时线图(e)245 m深度处真实速度、层析速度和初始速度曲线对比. Fig. 4 (a) True velocity;(b) tomographic inversion velocity;(c) Ray paths and isochrones corresponding to true velocity;(d) Ray paths and isochrones corresponding to tomographic velocity;(e) The compares of true velocity, tomographic velocity and initial velocity in depth of 245 m.

图5 (a)层析初始速度场和(b)建模速度场;(c)以采样点数表示的各道静校正量;(d)炮记录静校正前和(e)炮记录静校正校正后对比图. Fig. 5 (a) Initial velocity and (b) updated velocity;(c) The static value of each trace in the shot gathers;(d) and (e)The comparison between initial shot gathers and static correction shot gathers.
3.2 实际资料处理

实际资料是某一工区的三维地震资料，炮排列和检波排列东西向距离达到了1475 m，排列沿南北走向，每条排列50炮，炮点南北向起始大地坐标位置(4921787-4942587)，检波点排列南北向起始大地坐标位置(4914612-4946562)，速度场起始位置大地坐标4914600，沿南北向网格间隔100 m，深度采样间隔5 m，网格化后采样点数321*121.根据微测井和小折射资料，构建沿深度递增的梯度初始模型，梯度变化量25 m/s，图5(a)左图为初始速度场，速度范围(300~3300 m/s)，图5(b)为层析建模速度场，速度从南到北高程变化剧烈，呈现南低北高的趋势，北侧低降速带范围埋深更浅.这里选择500 m/s速度作为低降速带和高速层的分解速度，基准面埋深150 m，位于起伏界面以下的一水平界面，以采样点数表示的各道静校正量如图5(c)所示，图中可以看到，该线从南到北静校正量逐渐变大.选取的炮记录5(d)(e)所示，静校正前同相轴偏离双曲规律，抖动比较剧烈，初至同相轴不连续，经过静校正后，抖动消除，同相轴更加符合双曲规律，初至同相轴连续性更好，因此本文提出的方法适用于实际资料处理.

4 结 论

本文把基于抛物旅行时插值的最短路径射线追踪方法成功的应用于初至波走时层析速度建模方法，准确的计算初至走时和层析核函数，保证了层析建模的精度.理论模型试算证明了该方法可以有效的反演近地表速度模型，基于该方法的实际资料静校正处理效果显著.

参考文献
[1] Chen Z W, Zhang G Q. 2004. A tomography method for modification of migration velocity model using γ gathers[J]. Progress in Geophysics, 19(1): 154-160.
[2] Cheng J B, Wang H Z, Ma Z T. 2005. Double square root equation migration methods of narrow azimuth seismic data[J]. Chinese J. Geophys.   (in Chinese), 48(2): 399-405.
[3] Cheng J B, Wang N, Ma Z T. 2009. Table-driven 3-D angle-domain imaging approach for Kirchhoff prestack time migration[J]. Chinese J. Geophys.   (in Chinese), 52(3): 792-800.
[4] Dong C H. 2008. Prestack-Depth-Migration velocity analysis method based on the CFP method[J]. Progress in Geophysics, 23(6): 1866-1871.
[5] Dong C H, Zhang J F. 2010. Prestack time migration incorporated with redatuming[J]. Chinese J. Geophys.   (in Chinese), 53(10): 2435-2441.
[6] Feng B, Wang H Z. 2011. 3D offset plane-wave finite-difference pre-stack time migration[J]. Chinese J. Geophys.   (in Chinese), 54(11): 2916-2925.
[7] Huang Z Y, Qu S L, Wang Y J, et al. 2009. Kirchhoff prestack time migration of PS-wave data for the layered anisotropic medium[J]. Chinese J. Geophys.   (in Chinese), 52(12): 3109-3115.
[8] Li B, Liu G F, Liu H. 2009. A method of using GPU to accelerate seismic pre-stack time migration[J]. Chinese J. Geophys.   (in Chinese), 52(1): 245-252.
[9] Liu S W, Wang H Z, Cheng J B, et al. 2008. Space-time-shift imaging condition and migration velocity analysis[J]. Chinese J. Geophys.   (in Chinese), 51(6): 1883-1891.
[10] Liu L, Liang G H, Fu C, et al. 2011. Bend-ray Kirchhoff pre-stack time migration based on Chebyshev polynomial[J]. Chinese J. Geophys.   (in Chinese), 54(10): 2665-2672.
[11] Lu B K, Zhang J F, Pu B L. 2012. Pre-stack time migration in angle-domain: a one-way propagator approach in inhomogenenous media[J]. Chinese J. Geophys.   (in Chinese), 55(11): 3795- 3804.
[12] Nakanishi I, Yamaguchi K. 1986. A numerical experiment on nonlinear image reconstruction from first-arrival times for two dimensional island arc structure [J].   Journal of Physical Earth, 34:195-201.
[13] Qin N, Li Z C, Yang X D, et al. 2012. Image domain travel-time tomography velocity inversion based on automatic picking[J].   Oil Geophysical Prospecting (in Chinese), 47(3): 392-398.
[14] Sang Y Y, Li Z C, Zhang K. 2013. Shortest path raytracing based on parabolic traveltime interpolation[J].   Oil Geophysical Prospecting (in Chinese), 48(3): 403-409.
[15] Stork C, Clayton R W. 1991. Linear aspects of tomographic velocity analysis[J].   Geophysics, 56(4): 483-495.
[16] Sun F L, Wang Z L, Hao T Y, et al. 2011. Seismic imaging of complicated deep structures in southern South China Sea[J]. Chinese J. Geophys.   (in Chinese), 54(12): 3210-3216.
[17] Wang C L, Zhang S L, Zhao J X, et al. 2007. Synthetic source record interactive residual migration velocity analysis based on controlled illumination[J]. Chinese J. Geophys.   (in Chinese), 50(3): 860-867.
[18] Wang H Z, Cai J X, Kong X N, et al. 2010. An implementation of Kirchhoff integral prestack migration for large-scale data[J]. Chinese J. Geophys.   (in Chinese), 53(7): 1699-1709.
[19] Wang W H, Lin C H, Chen Z D, et al. 2011. Seismic data imaging techniques and their application for igneous rock formation in Gulong fault depression[J]. Chinese J. Geophys.   (in Chinese), 54(2): 310-319.
[20] Zhang H, Zhang J F. 2012. 3D prestack time migration including surface topography[J]. Chinese J. Geophys.   (in Chinese), 55(4): 1335-1344.
[21] Zhang J Z, Chen S J, Xu C W. 2004. A method of shortest path raytracing with dynamic networks[J]. Chinese J. Geophys.   (in Chinese), 47(5): 889-904.
[22] Zhang K, Yin Z, Li Z C, et al. 2013,Wave equation tomographic velocity inversion method based on the Born/Rytov approximation[J].   Applied Geophysics, 10(3): 314-322.
[23] 程玖兵, 王华忠, 马在田. 2005. 窄方位地震数据双平方根方程偏移方法探讨[J].   地球物理学报, 48(2): 399-405.
[24] 程玖兵, 王楠, 马在田. 2009. 表驱三维角度域Kirchhoff叠前时间偏移成像方法[J].   地球物理学报, 52(3): 792-800.
[25] 陈湛文, 张关泉. 2004. 用γ道集信息修正偏移速度模型的层析成像法[J].   地球物理学进展, 19(1): 154-160.
[26] 董春晖. 2008. 共聚焦点叠前深度偏移速度估计方法研究[J].   地球物理学进展, 23(6): 1866-1871.
[27] 董春晖, 张剑锋. 2010. 结合基准面重建的叠前时间偏移方法[J].   地球物理学报, 53(10): 2435-2441.
[28] 冯波, 王华忠. 2011. 三维偏移距平面波有限差分叠前时间偏移[J].   地球物理学报, 54(11): 2916-2925.
[29] 黄中玉, 曲寿利, 王于静,等. 2009. 层状各向异性介质转换波克希霍夫叠前时间偏移[J].   地球物理学报, 52(12): 3109-3115.
[30] 李博, 刘国峰, 刘洪. 2009. 地震叠前时间偏移的一种图形处理器提速实现方法[J].   地球物理学报, 52(1): 245-252.
[31] 刘守伟, 王华忠, 程玖兵,等. 2008. 时空移动成像条件及偏移速度分析[J].   地球物理学报, 51(6): 1883-1891.
[32] 刘璐, 梁光河, 符超,等. 2011. 基于Chebyshev多项式的弯曲射线Kirchhoff叠前时间偏移[J].   地球物理学报, 54(10): 2665-2672.
[33] 卢宝坤, 张剑锋, 蒲泊伶. 2012. 角度域叠前时间偏移-非均匀介质中的单程波方法[J].   地球物理学报, 55(11): 3795-3804.
[34] 秦宁, 李振春, 杨晓东,等. 2012. 自动拾取的成像空间域走时层析速度反演[J].   石油地球物理勘探, 47(3): 392-398.
[35] 桑运云, 李振春, 张凯. 2013. 基于抛物旅行时插值的最短路径射线追踪[J].   石油地球物理勘探, 48(3): 403-409.
[36] 孙福利, 王真理, 郝天珧,等. 2011. 南海南部深部结构的复杂构造地震成像[J].   地球物理学报, 54(12): 3210-3216.
[37] 王昌龙, 张叔伦, 赵景霞,等. 2007. 基于控制照明的合成震源记录交互剩余偏移速度分析[J].   地球物理学报, 50( 3): 860- 867.
[38] 王华忠, 蔡杰雄, 孔祥宁,等. 2010. 适于大规模数据的三维Kirchhoff积分法体偏移实现方案[J].   地球物理学报, 53(7): 1699-1709.
[39] 王维红, 林春华, 陈志德,等. 2011. 古龙断陷深层火山岩地震资料成像方法及应用研究[J].   地球物理学报, 54(2): 310-319.
[40] 张浩, 张剑锋. 2012. 起伏地表采集数据的三维直接叠前时间偏移方法[J].   地球物理学报, 55(4): 1335-1344.
[41] 张建中, 陈世军, 许初伟. 2004. 动态网络最短路径射线追踪[J].   地球物理学报, 47(5): 889-904.