地球物理学报  2017, Vol. 60 Issue (2): 541-553   PDF    
起伏地形下的高精度反射波走时层析成像方法
张新彦1 , 徐涛2,3 , 白志明2 , 高锐1 , 李秋生1 , 刘有山2 , 张智4     
1. 中国地质科学院地质研究所, 国土资源部深部探测与地球动力学重点实验室, 北京 100037;
2. 中国科学院地质与地球物理研究所, 岩石圈演化国家重点实验室, 北京 100029;
3. 中国科学院青藏高原地球科学卓越创新中心, 北京 100101;
4. 桂林理工大学广西矿治与环境科学实验中心, 广西隐伏金属矿产勘查重点实验室, 广西桂林 541004
摘要: 全球造山带及中国大陆中西部普遍具有强烈起伏的地形条件.复杂地形条件下的地壳结构成像问题像一面旗帜引领了当前矿产资源勘探和地球动力学研究的一个重要方向.深地震测深记录中反射波的有效探测深度可达全地壳乃至上地幔顶部,而初至波通常仅能探测上地壳浅部.为克服和弥补初至波探测深度的不足,本文基于前人对复杂地形条件下初至波成像的已有研究成果,采用数学变换手段将笛卡尔坐标系的不规则模型映射到曲线坐标系的规则模型,并将快速扫描方法与分区多步技术相结合,发展了反射波走时计算和射线追踪的方法.进而利用反射波走时反演,实现起伏地形下高精度的速度结构成像,从而为起伏地形下利用反射波数据高精度重建全地壳速度结构提供了一种全新方案.数值算例从正演计算精度、反演中初始模型依赖性、反演精度、纵横向分辨率以及抗噪性等方面验证了算法的正确性和可靠性.
关键词: 起伏地表      贴体网格      分区多步      射线追踪      走时层析成像     
High-precision reflection traveltime tomography for velocity structure with an irregular surface
ZHANG Xin-Yan1, XU Tao2,3, BAI Zhi-Ming2, GAO Rui1, LI Qiu-Sheng1, LIU You-Shan2, ZHANG Zhi4     
1. Key Laboratory of Earthprobe and Geodynamics, MLR, Institute of Geology, Chinese Academy of Geological Sciences, Beijing 100037, China;
2. State Key Laboratory of Lithospheric Evolution, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
3. CAS Center for Excellence in Tibetan Plateau Earth Sciences, Beijing 100101, China;
4. Guangxi Scientific Experiment Center of Mining, Metallurgy and Environment, Guangxi Key Laboratory of Hidden Metallic Ore Deposits Exploration, Guilin University of Technology, Guangxi Guilin 541004, China
Abstract: Global orogenic belts and the regions of middle and western China have strongly varied topography. Crustal structure imaging in such condition is of great interest in both the mineral resource exploration and the geodynamics study. The probing depth of seismic reflection waves in deep seismic sounding can reach up to the whole crust and even the top of the upper mantle of the earth, which is an effective supplement to the first-arrival refraction waves traveling in shallower depth. Current methods to deal with irregular surface usually take the step shape approximation or use a low velocity layer to overlay the surface. However, these methods may cause accuracy loss or even imaging distortion. This paper aims to introduce a grid-based method for traveltime tomography using the reflection traveltimes, to invert for velocity structure of the whole crust with an irregular surface. To deal with the irregular surface problem, this study is based on an irregular surface flattening scheme. The irregular surface flattening, which involves the transformation between the curvilinear and the Cartesian coordinate systems, can describe rugged terrain with no accuracy loss. Wave-front traveltime is calculated using the fast-sweeping method (FSM) to solve the topography-dependent Eikonal equation, and a multistage technique is then applied to calculate the reflection waves by reinitializing FSM in the incident layer. Ray paths are found by following the steepest traveltime gradient from the receiver to the interface and then to the source. Using the back-projection algorithm, the slowness perturbations are obtained, which have a relationship with the already existing traveltime. The final velocity model is obtained after a few iterations during which the slowness is updated, until the fit of the reflection traveltimes is satisfactory. We take four numerical examples to verify our method. First, the perfect match of traveltimes and ray paths obtained by our method and the shortest path method verifies that our forward modeling of the reflection traveltime tomography is accurate and reliable. Second, we test the dependency of the traveltime inversion on the initial velocity model. All four initial models converge to the true model within fairly good accuracy, which shows that our method is robust and effective with regard to the choice of the initial models. Third, our method is applied to a relatively complex model with strong variation of topography, and the final inverted velocity model gives correct ray paths of the reflections, very low traveltime residual, and high resolution velocity structure. Finally, both the full recovery of checkerboard patterns and the high-precision results inversed with the reflection traveltimes with Gaussian noise added show robustness and effectiveness of the proposed method. In summary, this paper addresses the problem of the irregular surface in seismic traveltime tomography. Based on the body-conforming grid and coordinate transformations, the wave-fronts of the reflection waves beneath the irregular surface are traced with high accuracy using the FWM and the multistage technique. Consequently, the first-arrival traveltime tomography is extended to the reflection traveltime tomography to image the deep structure of the crust in regions with complex surface topography. Numerical experiments show its reliability and effectiveness. This method can thus have a great potential application in mountain regions and basins, where high accuracy and high resolution seismic imaging is crucial to unravel their structure and reveal the tectonic evolution..
Key words: Irregular surface      Body-conforming grid      Multistage technique      Ray tracing      Traveltime tomography     
1 引言

全球造山带及中国大陆中西部普遍具有强烈起伏的地形条件.随着油气能源、矿产资源勘探和地球动力学研究的日渐深入,剧烈起伏的地形变化给盆山耦合区高精度地震探测带来极大挑战(Teng et al., 1987, 2003; Kaila and Krishna, 1992; Kaila and Krishna, 1992; Zeng et al., 1995; Li and Mooney, 1998; Gao et al., 2000, 2005; Artemieva, 2003; Thybo et al., 2003; Carbonell et al., 2004; Wu et al., 2005; 高锐等, 2000, 2004, 2006; Zhao et al., 2006; Zhang and Wang, 2007; Tian et al., 2010; Zhang and Klemperer, 2005, 2010; Zhang et al., 2011; Lei et al., 2006, 2011).传统基于平缓构造的探测方法在处理复杂地表区域采集的地震数据时受到很大限制,强烈起伏的地形常常会造成传统成像方法精度的损失甚至成像的失真.研究复杂地表条件下地震波的走时计算和层析反演问题,对于重建这些地区高精度和高分辨率的地壳地震波速度结构具有重要的意义.

起伏地表问题,通常与浅层勘探密切相关,处理方式多为基于离散介质的波场理论或走时场理论.通常情况下,对介质进行网格离散剖分,包括规则网格和不规则网格.到目前为止,绝大多数可以处理起伏地表的地震波走时计算和层析反演方法都是基于非规则网格的(Fomel, 1997; Kimmel and Sethian, 1998; Sethian, 1999; Sethian and Vladimirsky, 2000; Leliévre et al., 2001; Kao et al., 2008; Qian et al., 2007, 2008).然而,无论是网格剖分阶段还是走时计算阶段,相比于基于规则网格的方法,基于非规则网格的方法需要更多的计算量.传统的基于规则网格的起伏地表处理方法中,典型的是模型扩展方法(Vidale, 1988, 1990; Reshef, 1991; Hole, 1992),该方法的关键是把起伏地表处理成模型内部的分界面,通过在起伏地表上方填充低速介质,把不规则模型扩展为规则模型.其缺点在于起伏地表的阶梯状描述不但会造成精度的损失,而且填充低速的方法可能会产生许多虚假的绕射进而导致成像失真(Ma and Zhang, 2014b).

近几年,在贴体网格模型参数化基础上的地形“平化”策略,广泛应用于处理起伏地表问题(Haines, 1998; Lan and Zhang, 2011, 2012, 2013; 兰海强等, 2011, 2012).借助贴体网格和坐标变换,将笛卡尔坐标系中的物理模型和程函方程变换到曲线坐标系,把物理空间的不规则区域转化为计算空间的规则区域,理论上实现了对起伏地表精度无损的处理(图 1).

图 1 地形“平化”策略:贴体网格下的笛卡尔坐标系和曲线坐标系的转换 Fig. 1 Irregular surface flattening scheme: the transformation between Cartesian coordinates and curvilinear coordinates based on body-conforming grid

与传统连续介质理论上的射线方法(Cervény, 1987; Zelt and Smith, 1992; Zelt et al., 1999; Zelt, 1999; 徐涛等, 2004; Bai and Zhang, 2007; 李飞等,2013; Xu et al., 2006, 2010, 2014)相比,基于离散介质的程函方程数值解方法(Hole, 1992; Aldridge and Oldenburg, 1993; Wang, 2000, 2007Rawlinson and Sambridge, 2003, 2004a, b),采用“矩形网格”划分的模型参数化方法,健壮性好,适应性强,允许对模型进行高精度网格剖分.此外该方法可有效克服传统射线方法难以避免的低速影区问题,获得稳定的初至波走时,因此是非常理想的初至波成像方法.利用贴体网格模型上的初至波成像方法,可以对具有起伏地表的浅层或上地壳速度结构进行准确成像(Lan and Zhang, 2011, 2012; Ma and Zhang, 2014a, b, c, 2015).

初至波走时成像方法通常利用来自浅表沉积盖层的回折波或结晶基底顶部的折射波信息,主要应用于探测上地壳顶部地质体的分布特征(Hole, 1992; 王椿镛等,1995赵烽帆等,2014; Ma and Zhang, 2014a, b, c, 2015),有效探测深度一般局限于上地壳.与此相比,深地震测深记录上通常可观察到来自地壳不同深度界面的清晰反射波.这些反射波的走时数据同时携带了地壳深部介质的速度分布和物性间断面(界面)形态等信息,其探测深度和能力可达到全地壳,乃至上地幔顶部反射界面.因此充分利用地震资料中的反射波信息开展走时层析成像是了解地壳深部介质物理特性和岩石层状态的重要途径,它不仅有效弥补和克服了初至波走时成像探测深度的不足,而且也是对偏移成像方法的天然补充(Farra and Madariage, 1988; Scarascia and Cassinis 1997; Li and Mooney, 1998; Knapp et al., 2004Farquharson et al., 2008; Zhang et al., 2010; Rawlinson and Goleby, 2012).

2004年Rawlinson和Sambridge运用分区多步技术(Multistage scheme+different computation domains),并与快速行进法(Fast Marching Method)相结合,实现了二维层状介质中的多次透射、反射波的追踪和波前数值模拟(Rawlinson, 2004).Bai等将分区多步技术和最短路径方法(Shortest Path Method)相结合,实现了二维和三维的多次透射、反射波的追踪(Bai et al., 2009, 2010; 唐小平和白超英, 2009a, b).

为充分利用深地震测深资料中的反射波信息高精度重建壳幔速度结构,本文将基于前人对初至波成像的已有研究成果,沿用贴体网格模型参数化及曲线坐标系下初至波走时计算和射线追踪的方法,结合分区多步计算技术,完成起伏地表条件下层状介质中反射波的追踪计算,进而发展复杂地表条件下反射波走时反演的方法.

2 贴体网格模型参数化与地形“平化”策略

模型参数化是地震波走时层析成像的基础.对于复杂地表的介质,地形“平化”策略采用贴体网格剖分模型,贴体网格的网格边界与起伏的地表吻合以避免人为地产生边界散射,可以实现对起伏地形的精度无损描述(Thompson, 1985; Hvid, 1994; Lan and Zhang 2013a, b).贴体网格通过由计算空间到物理空间的坐标变换来获得.通过坐标变换将曲线坐标系(q, r)映射到物理空间的笛卡尔坐标系(x, z)(图 1).

生成贴体曲线网格的方法主要有代数法、保角变换法及偏微分方程法等(Thompson, 1985Hvid, 1994; 蒋丽丽,2009;孙建国,2009).其中偏微分方程法因其能方便地控制所生成网格点的疏密,且所生成的网格在边界处有良好的正交性而得到广泛的应用.本文中采用偏微分方程法生成贴体网格.

贴体网格生成之后,笛卡尔坐标系中的网格点与曲线坐标系中的网格点一一对应,即

(1)

由链锁规则,有

(2a)

(2b)

(2c)

(2d)

这里qx表示∂q(x, z)/∂xqzrxrz的意义也类似.这些导数称为度量导数,把公式(2a), (2b)分别代入(3a), (3b),经过化简,可得

(3)

J=xqzr-xrzq是变换的雅克比矩阵.值得注意的是,即使映射关系式(1)有解析的形式,度量导数依然通过数值的方式来计算,以避免当使用守恒形式的动量方程时系数偏导数引起的假源项(Thompson et al., 1985).本文中所有例子的度量导数均是采用二阶差分计算的.

3 反射波射线追踪和走时计算

以下描述基于地形“平化”策略的初至波和反射波走时计算和射线追踪方法.

3.1 与地形相关程函方程计算初至波走时场

传统的程函方程为

(4)

方程(4)为笛卡尔坐标下的各向同性介质情况.方程给出了射线经过慢度为s(x, z)的介质中的点(x, z)时的走时T(x, z).贴体网格模型“平化”策略下,利用方程(2a), (2b)和(4),对应的程函方程变换为(Lan and Zhang, 2013a, b):

(5)

其中,C=,参数A,B和C是与地形相关的.当地表水平时,xr, zq为0,xq, zr为1,此时,雅克比J=1,上述程函方程就简化为笛卡尔坐标系下经典的程函方程.对该程函方程采用Lax-Friedrichs数值哈密尔顿的Gauss-Seidel (G-S)扫描型算法(Kao et al., 2004)来求解,实现了对起伏地形的初至波走时场的模拟.

3.2 初至波走时梯度和射线路径计算

在初至波走时场基础上,需要沿着走时场的负梯度来进行射线追踪.笛卡尔坐标系下的走时梯度∇(T)(Vidale, 1988)为

(6)

其中,ik分别xz和方向走时梯度分量的单位向量.利用链式法则(2)式和度量导数(3)式,可得曲线坐标系下相应的走时梯度分量.

(7)

(8)

其中,xq表示∂x(q, r)/∂q,其余的意义也类似.当地表水平时,xr, zq为0,xq, zr为1,此时,雅克比J=1,该走时梯度分量就简化为笛卡儿坐标系下的形式.

由于初至波射线路径走时最小,所以可以沿着走时梯度最速下降的方向,从检波器逆向追踪到炮点,获得初至波的射线路径(Ma and Zhang, 2014).

3.3 分区多步反射波走时计算和射线追踪原理

分区多步方法是实现离散介质中反射波和多次波计算的高效技术(Rawlinson and Sambridge, 2004; Bai et al., 2009, 2010; 唐小平等, 2009a, b).本文采用快速扫描法处理曲线坐标系下的程函方程数值解(Lan and Zhang, 2013a, b),并结合分区多步方法来实现二维层状介质中反射波的追踪计算.

首先,对于二维模型,我们用一维网格点来描述界面,它与用于模型划分的矩形网格点相互独立,可以与之不重合,反射界面点的走时由该点所在的网格单元的节点走时线性内插得到(图 2a).

图 2 分区多步反射波走时计算和射线追踪示意图 (a)反射界面(黑色圆点)用一系列一维节点描述,并与模型网格点相互独立.背景网格为用于模型参数化的贴体网格.五角星S为炮点,倒三角形R1R2为接收器,O1O2为界面反射点;(b)下行波走时场及射线追踪;(c)上行波走时场及射线追踪. Fig. 2 Schematic diagram for reflectiontraveltime calculation and ray tracing, using Multistage Scheme (a) Interface nodes (black circles) are defined by a serial of 1D nodes, independent of velocity nodes.Grids on background are the body-conforming grids used to parameter the model; The star (S) denotes the energy source and the inverted triangles (R1, R2) indicate the seismic receivers deployedon the irregular topographic surface. The reflected point aremarked by O1 and O2; (b) Traveltime filed of downgoing wave and ray tracing; (c) Traveltime filed of upgoing wave and ray tracing

分区多步计算分为两个步骤:(1)分区:按照速度模型的分布进行物理分区,即按照模型速度界面将速度模型分为相对独立的层状或块状区域;(2)多步计算:在每个分区内进行上行波/下行波走时场的单独计算,然后按照波传统路径将这些独立分区的计算结果通过速度界面有机链接起来实现反射波、多次波等的追踪计算.由于本文只涉及反射波的模拟,故只在单层介质中进行上行波和下行波的计算.

单层介质的反射波分区多步计算步骤为:(1)下行波前的计算(图 2b).由震源开始用快速扫描法扫描波前,当该分区的波前扫描计算结束时,波前停止在该区的反射界面上,同时保存反射界面上各点的走时和该区域的走时场(用于射线追踪计算).(2)以界面离散点构成新的震源,扫描该分区的上行波波前(图 2c),扫描结束后,保存接收点处的走时(记录的反射波走时)和该区域的走时场(用于射线追踪).

快速扫描法中,走时计算和射线追踪通常分开进行.在获得了上行波和下行波的波前走时后,沿着走时梯度最陡下降的方向,由检波器到反射界面再由反射界面到炮点(R→O, O→S),反向追踪获得反射波的射线路径(图 2c, 2b).

4 反射波走时反演 4.1 反射波走时反演

利用深地震测深记录的反射波走时数据反演岩石层速度分布通常分为三个过程:(1)反演岩石层平均速度及其底界面平均深度;(2)固定岩石层底界面形态反演岩石层速度横向变化特征;(3)固定岩石层速度分布,反演其底部界面形态.为简单起见,本论文仅阐述固定岩石层底界面形态反演层速度分布的过程.

在通常的人工源地震探测试验中,观测系统采取地表激发和地表接收的方式.在理论模拟中,我们把震源放在起伏地形的表层,接收器也均匀分布于表层以接收来自各炮点的地震波信息(图 3).

图 3 起伏地表模型中的反射波地震勘探示意图 五角星为炮点,灰色倒三角为接收器,灰色背景为速度场,射线路径和走时场由黑实线和灰色虚线标出. Fig. 3 Schematic diagram of seismic prospecting using reflected waves The star denotes the energy source and the inverted triangles indicate the seismic receivers deployed on the irregular topographic surface; Raypaths and traveltime contours (wavefronts) are shown as black solid line and gray line, respectively

地震波的走时T可表示为

(9)

式中,s(x, z)是慢度,或速度的倒数.因为积分是沿射线路径l[s(x, z)]进行的,而射线路径取决于慢度s(x, z),所以,Ts的关系是非线性的.

若给参考慢度s0(x, z)以微小的慢度扰动δs,即s(x, z)=s0(x, z)+δs(x, z),则根据费马原理,走时扰动值可写为

(10)

在大多数情况下,方程(10)中走时残差和慢度扰动可以写成如下线性关系,即

(11)

其中,L是射线路径矩阵,Δ S是慢度扰动矩阵,Δ T是走时残差矩阵.因此,当获得了走时残差和相应的射线路径之后,慢度扰动可以通过求解方程组(11)计算得到.依此思想,即实现了非线性问题(9)的线性化.

4.2 反投影反演方法

求解此线性问题,本文采用反投影算法(Humphreys and Clayton, 1988; Hole, 1992).同其他的反演算法相比,此方法不进行复杂偏导数的计算和大型矩阵的反演,具有计算简单、占用内存少和计算效率高的优点.这种线性迭代的反演方法使正演与反演计算的成本最小,可密集采样,最终可得到速度结构的高分辨率图像.

其基本思想为.

用矩形网格剖分模型空间,用程函方程有限差分的方法计算初至波走时场,根据走时场追踪射线路径,然后将走时残差平均分配到射线路径上,得到各网格单元内的慢度扰动值,即

(12)

其中,δs为节点慢度修正值,K为穿过该节点的单元射线数,δTk为第k根射线的走时残差,lk为第k根射线的长度.

反演过程中的重采样和平滑可以降低个别网格上的大慢度扰动值而将其扩展到邻近的区域,增强网格间的约束性,减少虚假速度异常的出现.重采样和平滑因子的选择,需要保证模型依然能够很好的参数化和具有良好的分辨率.得到研究区域内各网格单元的慢度扰动值后,将其用于模型更新并作为下次迭代的初始模型.多次迭代直到走时残差达到稳定值且满足要求时停止迭代,此时得到的模型即为反演最终模型.

同时,为了保持速度模型的稳定更新和避免虚假速度异常的出现,反演采取先一维反演再二维反演的步骤.一维反演时保持速度横向不变,且横向光滑长度为整个排列长度,只进行速度在深度上的反演.得到一维反演的速度模型后,将其作为初始模型进行二维反演.

5 数值算例

通过以下数值算例分别验证了正、反演计算的正确性和有效性.

5.1 反射波正演走时计算和射线追踪

为了考察本文的反射波走时计算和射线追踪方法,选择一个起伏地表模型进行实验.模型表面是由四个山丘三个凹陷组成的起伏型地表.模型大小为100 km×50 km,在40 km深度处为一水平反射界面.给该模型以向下增大的速度场,在模型左侧放炮,单边接收.一般地,复杂地表下的反射波走时很难用解析的方法进行刻画,所以我们选择采用最短路径方法(SPM)(Moster, 1991)计算出其数值解并和本方法计算的数值解对比,来考察本方法求解反射波走时和射线路径的正确性.图 4a为对该模型的贴体网格剖分图,网格间隔为1 km.最短路径方法用模型扩展和规则网格剖分,网格间隔为0.2 km,以视作真实解.将两种算法的计算结果进行对比,结果见图 4b.比较两个结果可得,虽然地形“平化”方法和最短路径方法的理论基础和计算方法不同,但两者走时和射线路径计算结果均吻合较好,其中反射波计算走时误差最大为0.2 s左右.与最短路径方法计算结果的对比,验证了地形“平化”策略下反射波正演算法的正确性,这为我们进行可靠的反演计算提供了保障.

图 4 正演走时计算模型试验 (a)四个山丘三个凹陷模型的贴体网格剖分图.为清楚起见,网格密度减小为原来的1/3,模型大小为100 km×50 km, 网格大小为101 km×51 km;(b)起伏地表模型地形“平化”方法和最短路径方法正演走时计算和射线追踪结果对比:(b1)走时计算结果红色十字叉:最短路径方法,黑色菱形:地表“平化”方法;(b2)射线追踪结果红实线:最短路径方法,黑实线:地形“平化”法,背景场:速度场. Fig. 4 Accuracy test of forward calculation (a) The body-conforming grids in the model with four hills and three valleys. For clarify, the grids are displayed with a reducing density factor of 3; (b) The comparison of the calculatedtraveltime and ray paths obtained by the irregular surface flattening method and the shortest path method: (b1) Calculated traveltimes. Red cross denotesthe result of the shortest path method while black diamond denotes the irregular surface flattening result; (b2) Ray paths. Red lines denote the ray paths obtained by the shortest path method while black lines present the ray paths of the irregular flattening method.
5.2 初始模型试验

线性迭代的走时反演方法是一种局部最优化方法,因此初始模型的选择可能会影响最终的反演结果.为了测试初始模型对反演结果的影响性,选择一组起伏地表模型进行实验.模型依然采用5.1节中的山丘凹陷模型,模型大小和网格剖分均与5.1节相同.10个炮点和46个检波器均匀置于起伏地表上.赋予该模型以四种不同速度场作为初始模型反演计算理论模型.图 5a中黑色实线和4个虚线分别为理论模型和4个初始模型的速度梯度.其中初始模型1(红色虚线)整体速度偏低,初始模型2(浅蓝虚线)整体速度偏高,初始模型3和4(深蓝虚线)速度上高下低,但3在下层区域与理论模型偏差较大,4在上层区域偏差较大.这四种模型几乎包含了初始模型与理论模型的所有大小关系,因此能够对初始模型的影响性进行全面分析.

图 5 不同初始模型测试结果图 (a)一维初始模型(彩色虚线)和理论速度(黑色实线);(b)一维反演模型(彩色实线)和理论模型(黑色实线);(c)二维理论模型;(d-g)二维反演模型.其中,(c-g)右侧为前一步的一维模型:理论模型(绿实线)、初始模型(黑虚线)和一维反演模型(黑实线). Fig. 5 Test of different initial models (a) The 1 dimension (1D) initial models (color dotted lines) and theoretical model (black solid line); (b) The 1D inversed model (color solid line) and theoretical model (black solid line); (c) The 2D theoretical model; (d-g) the 2D inversed model. In (c-g), the right panels are 1D models obtained in the former step: the theoretical model (green line), initial models (black dotted lines) and 1D inversed models (black solid lines).

反演过程均采用先一维后二维的反演流程,且在一维反演过程横向光滑长度为整个排列长度.图 5b为一维反演的结果,可以看到,与初始模型相比,反演模型均缩小了与理论模型的速度差.图 5(c-g)右侧图为理论模型(绿线)、初始模型(虚线)和一维反演模型(实线).左侧图为理论模型(c)和各模型二维反演结果(d-g).由图可知,不同的初始模型,经过先一维再二维的反演过程之后,均实现了对理论模型很好的重建,反演误差最大为4%左右.图 6为4个初始模型反演迭代过程中的均方根走时残差收敛情况,各模型走时残差均很快的收敛到0.036 ms左右并保持稳定.该实验验证了本算法的稳健性和较小的初始模型依赖性.

图 6 不同初始模型迭代过程中均方根走时残差图,最终0.036 ms左右保持稳定 Fig. 6 Rms traveltime residuals versus number of iterations for the four numerical examples
5.3 复杂起伏地形模型试验

为了考察本文的反射波走时反演方法,选取一个地表剧烈起伏的模型进行实验.模型大小为100 km×40 km,在20~35 km深度范围为一起伏反射界面.用贴体网格剖分模型,网格间隔为1 km×1 km,大小为101×40(图 7),10个炮点和46个接收点均匀分布于起伏地表上.反演过程中选择重新抽样(网格化)因子为(3, 3),滑动平均滤波器为(4, 2),用于水平和垂直方向的慢度扰动计算,在走时残差收敛且稳定时停止迭代.

图 7 贴体网格剖分图 为表述清楚,进行了稀疏采样,稀疏因子为3.炮点(五角星)和接收器(倒三角)均匀分布于起伏地表上. Fig. 7 Body-conforming grids in the model with strongly varied topography For clarify, the grids are displayed with a reducing density factor of 3. Shot points and receivers are placed on the irregular surface uniformly.

由反演结果(图 8)可以看出,即使在初始模型(图 8b)和理论模型(图 8a)相差较大时,通过迭代之后,反演结果在整个模型上也能实现对理论模型很好的恢复(图 8c).值得注意的是,在紧邻反射界面的下方区域,由于受到界面上方区域在反演时的平滑作用,所以也进行了模型更新,但由于没有真实数据约束,故而界面以下的模型不可靠.

图 8 反射波走时层析成像方法获得的反演结果及模型相对误差图 (a)理论模型; (b)初始模型; (c)最终反演模型; (d)与理论模型相比的初始模型相对误差图; (e)与理论模型相比的最终反演模型相对误差图. Fig. 8 The inversed result by the reflection traveltime tomography (a) The theoretical model; (b) The initial model; (c) The finally inversed model; (d) Velocity′s relativeerror of the initial model; (e) Velocity′s relative error of the final inversed model.

速度误差分布图(图 8d-e)表明,与理论模型相比,初始模型与理论模型速度差较大,最大相对速度差超过10%(图 8d).经过反演迭代之后,最终反演模型上,速度差降低到3%以下(图 8e),尤其是中层区域接近于零误差反演.图 9显示了最终反演模型上的射线分布和射线在各网格单元内的覆盖数,从中可以看出,在深度10 km左右的中层区域射线交叉较多,这与速度误差图上所显示的零误差是相一致的.由迭代走时残差图(图 10)可以看出,均方根走时残差由最初的0.83 s降为一维反演的0.03 ms,并经过二维反演后,降为0.002 ms左右保持稳定,说明反演结果已经收敛.对比图 10图 8(d-e),说明走时残差的收敛和速度误差的收敛是一致的.

图 9 (a)最终反演模型上第三炮反射波射线路径图; (b)最终反演模型上各网格单元内射线覆盖数 Fig. 9 (a) Ray paths of the 3th shot on the final inversed model; (b) Number of rays that intersect each grid cell on the final inversed model
图 10 (a)迭代均方根走时残差;(b)第1~9次迭代走时残差 Fig. 10 (a) Rms traveltime residuals versus number of iterations; (b) Rms traveltime residuals for interation 1~9

为了评价算法的纵横向分辨力,我们进行了检测板测试实验(图 11a).在背景速度场上加以间隔为7.5 km×7.5 km、大小为[0.25×sin (x)×sin (z)] km·s-1的速度异常场合成理论检测板,用背景速度场作为初始模型来反演理论检测板(图 11a).反演结果(图 11b)表明:在模型两侧和反射界面凹陷的区域由于射线覆盖较稀疏和方向交叉较少,检测板恢复稍差;在其他区域纵向和横向上速度异常均恢复较好.验证了该算法的较高纵向和横向分辨能力.

图 11 反射波走时层析成像检测板测试图 (a)理论检测板,速度异常值为[0.25×sin (x)×sin (z)] km·s-1,空间间隔为7.5 km×7.5 km;(b)反演检测板上速度异常场. Fig. 11 Checkerboard test for the reflection traveltime tomography scheme (a) The theoreticalcheckerboard, with velocity anomaly to be [0.25×sin (x)×sin (z)] km·s-1, and space scale to be 7.5 km×7.5 km; (b) The velocity anomaly on the inversed checkerboard.

在真实模型合成的走时数据上加以标准差为0.1 s的高斯白噪声(σ=0.1 s),用相同的初始模型进行反演(图 8b),最终得到的反演模型和速度相对误差分布见图 12b.最终反演模型与真实模型(图 8a)较为一致,与无噪声反演的结果(图 8c)具有可比性.说明该算法具有较高的抗噪能力.

图 12 噪声测试,在图 8a模型的反射波走时上加以标准差σ=0.1 s的高斯白噪声 (a)最终反演速度模型;(b)最终模型与真实模型相比的速度相对误差分布图. Fig. 12 Noise tolerance test, random Gaussian noise with standard deviation of σ=0.1 s was added to the exact traveltime data computed from the model inFig. 8(a). (a) The finally inversed model; (b) The relative velocity difference of the final inversed model, compared with the true model.

值得强调的是,与常规基于射线追踪方法重建地壳速度结构(Zelt et al., 1992, 1995, 1999)相比,我们发展的正反演方法由于采用起伏地形条件下的程函方程进行走时场的正反演计算,避免了射线方法容易出现的低速影区问题等,使得走时成像更为精确可靠.同时我们的方法允许对层间速度模型进行非常精细地剖分,这也是传统射线反演方法难以比拟的.典型如当前测深资料处理解释领域经常使用的Seis81系列(Cervény, 1981, 1987)及Rayinvr技术(Zelt, 1992, 1995)参数化模型均采用在速度层顶界面和底界面取速度和深度节点的参数化方式,难以对某一速度层内局部速度异常进行准确刻画.我们发展的方法则不存在此问题,可以根据射线覆盖及研究需要进行任意精细剖分,这使得本文方法相比之下具有显著优越性.可以预见,该方法在深地震测深资料处理解释领域有广泛应用空间.

全球造山带及中国大陆中西部普遍具有强烈起伏的地形条件,发展起伏地表下的地震波走时正反演计算方法对于盆山耦合区高精度地震探测和结构成像具有重要意义.本文借助于贴体网格和数学变换,在与地形有关的程函方程及其计算的初至波走时场基础上,结合分区多步技术,发展了起伏地形条件下反射波走时计算和射线追踪的方法,实现了复杂地表条件下高精度的速度结构成像,为起伏地形下利用反射地震资料重建地壳精细速度结构提供了一种全新方案.该方案在深地震测深资料处理和解释领域有广阔应用前景和价值.数值算例验证了正演走时计算、射线追踪及反演结果的正确性和可靠性,并通过检测板实验和抗噪实验检验了算法的分辨率和稳定性.

致谢

谨以此文纪念中国科学院地质与地球物理研究所张忠杰研究员(1964-2013).感谢滕吉文院士的悉心指导.感谢兰海强博士和马婷博士在算法上提供的支持和帮助.

参考文献
Aldridge D F, Oldenburg D W. 1993. Two-dimensional tomographic inversion with finite-difference traveltimes. Journal of Seismic Exploration, 2(3): 257-274.
Artemieva I M. 2003. Lithospheric structure, composition, and thermal regime of the East European Craton:implications for the subsidence of the Russian platform. Earth and Planetary Science Letters, 213(3-4): 431-446. DOI:10.1016/S0012-821X(03)00327-3
Bai C Y, Tang X P, Zhao R. 2009. 2-D/3-D multiply transmitted, converted and reflected arrivals in complex layered media with the modified shortest path method. Geophysical Journal International, 179(1): 201-214. DOI:10.1111/gji.2009.179.issue-1
Bai C Y, Huang G J, Zhao R. 2010. 2-D/3-D irregular shortest-path ray tracing for multiple arrivals and its applications. Geophysical Journal International, 183(3): 1596-1612. DOI:10.1111/j.1365-246X.2010.04817.x
Bai Z M, Zhang Z J, Wang Y H. 2007. Crustal structure across the Dabie-Sulu orogenic belt revealed by seismic velocity profiles. Journal of Geophysics and Engineering, 4(4): 436-442. DOI:10.1088/1742-2132/4/4/009
Carbonell R, Simancas F, Juhlin C, et al. 2004. Geophysical evidence of a mantle derived intrusion in SW Iberia. Geophysical Research Letters, 31: L11601.
Červený V. 1987. Ray tracing algorithms in three-dimensional laterally varying layered structures.//Nolet G. Seismic Tomography:With Applications in Global Seismology and Exploration Geophysics. Netherlands:Springer, 99-133.
Farquharson C G, Ash M R, Miller H G. 2008. Geologically constrained gravity inversion for the Voisey's Bay ovoid deposit. The Leading Edge, 27(1): 64-69. DOI:10.1190/1.2831681
Farra V, Madariage R. 1988. Non-linear reflection tomography. Geophysical Journal International, 95(1): 135-147. DOI:10.1111/j.1365-246X.1988.tb00456.x
Fomel S. 2000. A variational formulation of the fast marching eikonal solver. SEP-95:Stanford Exploration Project, 455-476.
Gao R, Huang D D, Lu D Y, et al. 2000. Deep seismic reflection profile across the juncture zone between the Tarim Basin and the West Kunlun Mountains. Chinese Science Bulletin, 45(24): 2281-2286. DOI:10.1007/BF02886369
Gao R, Huang D D, Lu D Y, et al. 2000. Deep seismic reflection profiling across the transitional region between western Kunlun orogenic belt and Tarim basin. Chinese Science Bulletin (in Chinese), 45(17): 1874-1879.
Gao R, Dong S W, He R Z, et al. 2004. Subduction process of the Yangtze continental block from Moho reflection image, South China. Earth Science Frontiers (in Chinese), 11(3): 43-49.
Gao R, Lu Z W, Li Q, et al. 2005. Geophysical survey and geodynamic study of crust and upper mantle in the Qinghai-Tibet Plateau. Episodes, 28(4): 263-273.
Gao R, Wang H Y, Ma Y S, et al. 2006. Tectonic relationships between the Zoigê Basin of the Song-pan Block and the West Qinling Orogen at lithosphere scale:Results of deep seismic reflection profiling. Acta Geoscientica Sinica (in Chinese), 27(5): 411-418.
Haines A J. 1988. Multi-source, multi-receiver synthetic seismograms for laterally heterogeneous media using F-K domain propagators. Geophysical Journal International, 95(2): 237-260. DOI:10.1111/j.1365-246X.1988.tb00465.x
Hole J A. 1992. Nonlinear high-resolution three-dimensional seismic travel time tomography. Journal of Geophysical Research, 97(B5): 6553-6562. DOI:10.1029/92JB00235
Humphreys E, Clayton R W. 1988. Adaptation of back projection tomography to seismic travel time problems. Journal of Geophysical Research, 93(B2): 1073-1085. DOI:10.1029/JB093iB02p01073
Hvid S L. 1994. Three dimensional algebraic grid generation[Ph.D.thesis]. Denmark:Technical University of Denmark.
Jiang L L, Sun J G. 2008. Source terms of elliptic system in grid generation. Global Geology (in Chinese), 27(3): 298-305.
Kaila K L, Krishna V G. 1992. Deep seismic sounding studies in India and major discoveries. Current Science, 62(1-2): 117-154.
Kao C Y, Osher S, Qian J L. 2004. Lax-Friedrichs sweeping scheme for static Hamilton-Jacobi equations. Journal of Computational Physics, 196(1): 367-391. DOI:10.1016/j.jcp.2003.11.007
Kimmel R, Sethian J A. 1998. Computing geodesic paths on manifolds. Proceedings of the National Academy of Sciences of the United States of America, 95(15): 8431-8435. DOI:10.1073/pnas.95.15.8431
Knapp C C, Knapp J H, Connor J A. 2004. Crustal-scale structure of the South Caspian Basin revealed by deep seismic reflection profiling. Marine and Petroleum Geology, 21(8): 1073-1081. DOI:10.1016/j.marpetgeo.2003.04.002
Lan H Q, Zhang Z J. 2011a. Comparative study of the free-surface boundary condition in two-dimensional finite-difference elastic wave field simulation. Journal of Geophysics and Engineering, 8(2): 275-286. DOI:10.1088/1742-2132/8/2/012
Lan H Q, Zhang Z J. 2011b. Three-dimensional wave-field simulation in heterogeneous transversely isotropic medium with irregular free surface. Bulletin of the Seismological Society of America, 101(3): 1354-1370. DOI:10.1785/0120100194
Lan H Q, Liu J, Bai Z M. 2011. Wave-field simulation in VTI media with irregular free surface. Chinese J. Geophys. (in Chinese), 54(8): 2072-2084. DOI:10.3969/j.issn.0001-5733.2011.08.014
Lan H Q, Zhang Z J. 2012. Seismic wavefield modeling in media with fluid-filled fractures and surface topography. Applied Geophysics, 9(3): 301-312. DOI:10.1007/s11770-012-0341-5
Lan H Q, Zhang Z, Xu T, et al. 2012. Influences of anisotropic stretching of boundary conforming grid on traveltime computation by topography-dependent eikonal equation. Chinese J. Geophys., 55(5): 564-579. DOI:10.1002/cjg2.v55.5
Lan H Q, Zhang Z, Xu T, et al. 2012a. Effects due to the anisotropic stretching of the surface-fitting grid on the traveltime computation for irregular surface by the coordinate transforming method. Chinese Journal of Geophysics (in Chinese), 55(10): 3355-3369. DOI:10.6038/j.issn.0001-5733.2012.10.018
Lan H Q, Zhang Z, Xu T, et al. 2012b. A comparative study on the fast marching and fast sweeping methods in the calculation of first-arrival traveltime field. Progress in Geophysics (in Chinese), 27(5): 1863-1870. DOI:10.6038/j.issn.1004-2903.2012.05.005
Lan H Q, Zhang Z J. 2013a. Topography-dependent eikonal equation and its solver for calculating first-arrival traveltimes with an irregular surface. Geophysical Journal International, 193(2): 1010-1026. DOI:10.1093/gji/ggt036
Lan H Q, Zhang Z J. 2013b. A high-order fast-sweeping scheme for calculating first-arrival travel times with an irregular surface. Bulletin of the Seismological Society of America, 103(3): 2070-2082. DOI:10.1785/0120120199
Lei J S, Zhao D P. 2006. Global P-wave tomography:On the effect of various mantle and core phases. Physics of the Earth and Planetary Interiors, 154(1): 44-69. DOI:10.1016/j.pepi.2005.09.001
Lei J S, Zhao D P, Fu X R, et al. 2011. An attempt to detect temporal variations of crustal structure in the source area of the 2006 Wen-An earthquake in North China. Journal of Asian Earth Sciences, 40(4): 958-976. DOI:10.1016/j.jseaes.2010.07.003
Lelièvre P G, Farquharson C G, Hurich C A. 2011. Computing first-arrival seismic traveltimes on unstructured 3-D tetrahedral grids using the Fast Marching Method. Geophysical Journal International, 184(2): 885-896. DOI:10.1111/gji.2011.184.issue-2
Li F, Xu T, Wu Z B, et al. 2013. Segmentally iterative ray tracing in 3-D heterogeneous geological models. Chinese J. Geophys. (in Chinese), 56(10): 3514-3522. DOI:10.6038/cjg20131026
Li S L, Mooney W D. 1998. Crustal structure of China from deep seismic sounding profiles. Tectonophysics, 288(1-4): 105-113. DOI:10.1016/S0040-1951(97)00287-4
Ma T, Zhang Z J. 2014a. Calculating ray paths for first-arrival travel times using a topography-dependent eikonal equation solver. Bulletin of the Seismological Society of America, 104(3): 1501-1517. DOI:10.1785/0120130172
Ma T, Zhang Z J. 2014b. A model expansion criterion for treating surface topography in ray path calculations using the eikonal equation. Journal of Geophysics and Engineering, 11(2): 025007. DOI:10.1088/1742-2132/11/2/025007
Ma T, Zhang Z J, Wang P, et al. 2014. Upper crustal structure under Jingtai-Hezuo profile in Northeastern Tibet from topography-dependent eikonal traveltime tomography. Earthquake Science, 27(2): 137-148. DOI:10.1007/s11589-013-0033-0
Ma T, Zhang Z J. 2015. Topography-dependent eikonal traveltime tomography for upper crustal structure beneath an irregular surface. Pure and Applied Geophysics, 172(6): 1511-1529. DOI:10.1007/s00024-014-0984-7
Moser T J. 1991. Shortest path calculation of seismic rays. Geophysics, 56(1): 59-67. DOI:10.1190/1.1442958
Qian J L, Zhang Y T, Zhao H K. 2007. Fast sweeping methods for eikonal equations on triangular meshes. SIAM Journal on Numerical Analysis, 45(1): 83-107. DOI:10.1137/050627083
Rawlinson N, Sambridge M. 2003. Irregular interface parametrization in 3-D wide-angle seismic traveltime tomography. Geophysical Journal International, 155(1): 79-92. DOI:10.1046/j.1365-246X.2003.01983.x
Rawlinson N, Sambridge M. 2004a. Multiple reflection and transmission phases in complex layered media using a multistage fast marching method. Geophysics, 69(5): 1338-1350. DOI:10.1190/1.1801950
Rawlinson N, Sambridge M. 2004b. Wave front evolution in strongly heterogeneous layered media using the fast marching method. Geophysical Journal International, 156(3): 631-647. DOI:10.1111/gji.2004.156.issue-3
Rawlinson N, Goleby B R. 2012. Seismic imaging of continents and their margins:New research at the confluence of active and passive seismology. Tectonophysics, 572-573: 1-6. DOI:10.1016/j.tecto.2012.07.021
Reshef M. 1991. Depth migration from irregular surfaces with depth extrapolation methods. Geophysics, 56(1): 119-122. DOI:10.1190/1.1442947
Scarascia S, Cassinis R. 1997. Crustal structures in the central-eastern Alpine sector:a revision of the available DSS data. Tectonophysics, 271(1-2): 157-188. DOI:10.1016/S0040-1951(96)00206-5
Sethian J A, Vladimirsky A. 2000. Fast methods for the Eikonal and related Hamilton-Jacobi equations on unstructured meshes. Proceedings of the National Academy of Sciences of the United States of America, 97(11): 5699-5703. DOI:10.1073/pnas.090060097
Sethian J A. 1999. Fast marching methods. SIAM Review, 41(2): 199-235. DOI:10.1137/S0036144598347059
Sun J G. 2007. Methods for numerical modeling of geophysical fields under complex topographical conditions:a critical review. Global Geology (in Chinese), 26(3): 345-362.
Sun J G, Jiang L L. 2009. Orthogonal curvilinear grid generation technique used for numeric simulation of geophysical fields in undulating surface condition. OGP (in Chinese), 44(4): 494-500.
Sun Z Q, Sun J G, Zhang D L. 2009. 2D DC electric field numerical modeling including surface topography using coordinate transformation method. Journal of Jilin University:Earth Science Edition, 39(3): 528-534.
Taillandier C, Noble M, Chauris H, et al. 2009. First-arrival traveltime tomography based on the adjoint-state method. Geophysics, 74(6): WCB1-WCB10. DOI:10.1190/1.3250266
Tang X P, Bai C Y. 2009a. Multiple ray tracing within 3-D layered media with the shortest path method. Chinese J. Geophys. (in Chinese), 52(10): 2635-2643. DOI:10.3969/j.issn.0001-5733.2009.10.024
Tang X P, Bai C Y. 2009b. Multiple ray tracing in 2D layered media with the shortest path method. Progress in Geophysics (in Chinese), 24(6): 2087-2096. DOI:10.3969/j.issn.1004-2903.2009.06.022
Teng J W, Wei S Y, Sun K Z, et al. 1987. The characteristics of the seismic activity in the Qinghai-Xizang (Tibet) Plateau of China. Tectonophysics, 134(1-3): 129-144. DOI:10.1016/0040-1951(87)90253-8
Teng J W, Zeng R S, Yan Y F, et al. 2003. Depth distribution of Moho and tectonic framework in eastern Asian continent and its adjacent ocean areas. Science in China Series D:Earth Sciences, 46(5): 428-446. DOI:10.1360/03yd9038
Thompson J F, Warsi Z U A, Mastin C W. 1985. Numerical Grid Generation:Foundations and Applications. New York:Elsevier North-Holland.
Thybo H, Janik T, Omelchenko V D, et al. 2003. Upper lithospheric seismic velocity structure across the Pripyat Trough and the Ukrainian Shield along the EUROBRIDGE'97 profile. Tectonophysics, 371(1-4): 41-79. DOI:10.1016/S0040-1951(03)00200-2
Tian X B, Teng J W, Zhang H S, et al. 2011. Structure of crust and upper mantle beneath the Ordos Block and the Yinshan Mountains revealed by receiver function analysis. Physics of the Earth and Planetary Interiors, 184(3-4): 186-193. DOI:10.1016/j.pepi.2010.11.007
Vidale J E. 1988. Finite-difference calculation of travel times. Bulletin of the Seismological Society of America, 78(6): 2062-2076.
Vidale J E. 1990. Finite-difference calculation of traveltimes in three dimensions. Geophysics, 55(5): 521-526. DOI:10.1190/1.1442863
Wang C Y, Lin Z Y, Chen X B. 1995. Comprehensive study of geophysics on geoscience transect from Menyuan, Qinghai province, to Ningde, Fujian province, China. Acta Geophysica Sinica (in Chinese), 38(5): 590-598.
Wang C Y, Zhang X K, Ding Z F, et al. 1997. Finite-difference tomography of upper crustal structure in Dabieshan orogenic belt. Acta Geophysica Sinica (in Chinese), 40(1): 495-502.
Wang C Y, Zeng R S, Mooney W D, et al. 2000. A crustal model of the ultrahigh-pressure Dabie Shan orogenic belt, China, derived from deep seismic refraction profiling. Journal of Geophysical Research, 105(B5): 10857-10869. DOI:10.1029/1999JB900415
Wang C Y, Han W B, Wu J P, et al. 2007. Crustal structure beneath the eastern margin of the Tibetan Plateau and its tectonic implications. Journal of Geophysical Research, 112: B07307.
Wu C L, Harris J M, Nihei K T, et al. 2005. Two-dimensional finite-difference seismic modeling of an open fluid-filled fracture:Comparison of thin-layer and linear-slip models. Geophysics, 70(4): T57-T62. DOI:10.1190/1.1988187
Xu T, Xu G M, Gao E G, et al. 2004. Block Modeling and Shooting Ray Tracing in Complex 3-D media. Chinese J. Geophys. (in Chinese), 47(6): 1118-1126.
Xu T, Xu G M, Gao E G, et al. 2006. Block modeling and segmentally iterative ray tracing in complex 3D media. Geophysics, 71(3): T41-T51. DOI:10.1190/1.2192948
Xu T, Zhang Z J, Gao E, et al. 2010. Segmentally iterative ray tracing in complex 2D and 3D heterogeneous block models. Bulletin of the Seismological Society of America, 100(2): 841-850. DOI:10.1785/0120090155
Xu T, Li F, Wu Z B, et al. 2014. A successive three-point perturbation method for fast ray tracing in complex 2D and 3D geological models. Tectonophysics, 627: 72-81. DOI:10.1016/j.tecto.2014.02.012
Xu T, Zhang Z J, Tian X B, et al. 2014. Crustal structure beneath the Middle-Lower Yangtze metallogenic belt and its surrounding areas:Constraints from active source seismic experiment along the Lixin to Yixing profile in East China. Acta Petrologica Sinica (in Chinese), 30(4): 918-930.
Zelt C A, Smith R B. 1992. Seismic traveltime inversion for 2-d crustal velocity structure. Geophysical Journal International, 108(1): 16-34. DOI:10.1111/gji.1992.108.issue-1
Zelt C A. 1999. Modelling strategies and model assessment for wide-angle seismic traveltime data. Geophysical Journal International, 139(1): 183-204. DOI:10.1046/j.1365-246X.1999.00934.x
Zelt C A, Hojka A M, Flueh E R, et al. 1999. 3D simultaneous seismic refraction and reflection tomography of wide-angle data from the central Chilean margin. Geophysical Research Letters, 26(16): 2577-2580. DOI:10.1029/1999GL900545
Zeng R S, Ding Z F, Wu Q J. 1995. A review on the lithospheric structures in the Tibetan Plateau and constraints for dynamics. Pure and Applied Geophysics, 145(3-4): 425-443. DOI:10.1007/BF00879582
Zhang J, Toksoz M N. 1998. Nonlinear refraction traveltime tomography. Geophysics, 63(5): 1726-1737. DOI:10.1190/1.1444468
Zhang J H, Wang W M, Wang S Q, et al. 2010. Optimized Chebyshev Fourier migration:A wide-angle dual-domain method for media with strong velocity contrasts. Geophysics, 75(2): S23-S34. DOI:10.1190/1.3350861
Zhang Z J, Klemperer S L. 2005. West-east variation in crustal thickness in northern Lhasa block, central Tibet, from deep seismic sounding data. Journal of Geophysical Research:Solid Earth (1978-2012), 110: B09403.
Zhang Z J, Klemperer S L. 2010. Crustal structure of the Tethyan Himalaya, southern Tibet:New constraints from old wide-angle seismic data. Geophysical Journal International, 181(3): 1247-1260.
Zhang Z J, Deng Y F, Teng J W, et al. 2011. An overview of the crustal structure of the Tibetan plateau after 35 years of deep seismic soundings. Journal of Asian Earth Sciences, 40(4): 977-989. DOI:10.1016/j.jseaes.2010.03.010
Zhao D P, Lei J S, Inouea T, et al. 2006. Deep structure and origin of the Baikal rift zone. Earth and Planetary Science Letters, 243(3-4): 681-691. DOI:10.1016/j.epsl.2006.01.033
Zhao F F, Ma T, Xu T. 2014. A review of the travel-time calculation methods of seismic first break. Progress in Geophysics (in Chinese), 29(3): 1102-1113. DOI:10.6038/pg20140313
高锐, 黄东定, 卢德源, 等. 2000. 横过西昆仑造山带与塔里木盆地结合带的深地震反射剖面. 科学通报, 45(17): 1874–1879.
高锐, 董树文, 贺日政, 等. 2004. 莫霍面地震反射图像揭露出扬子陆块深俯冲过程. 地学前缘, 11(3): 43–49.
高锐, 王海燕, 马永生, 等. 2006. 松潘地块若尔盖盆地与西秦岭造山带岩石圈尺度的构造关系-深地震反射剖面探测成果. 地球学报, 27(5): 411–418.
蒋丽丽, 孙建国. 2008. 基于Poisson方程的曲网格生成技术. 世界地质, 27(3): 298–305.
兰海强, 刘佳, 白志明. 2011. VTI介质起伏地表地震波场模拟. 地球物理学报, 54(8): 2072–2084. DOI:10.3969/j.issn.0001-5733.2011.08.014
兰海强, 张智, 徐涛, 等. 2012a. 贴体网格各向异性对坐标变换法求解起伏地表下地震初至波走时的影响. 地球物理学报, 55(10): 3355–3369. DOI:10.6038/j.issn.0001-5733.2012.10.018
兰海强, 张智, 徐涛, 等. 2012b. 地震波走时场模拟的快速推进法和快速扫描法比较研究. 地球物理学进展, 27(5): 1863–1870. DOI:10.6038/j.issn.1004-2903.2012.05.005
李飞, 徐涛, 武振波, 等. 2013. 三维非均匀地质模型中的逐段迭代射线追踪. 地球物理学报, 56(10): 3514–3522. DOI:10.6038/cjg20131026
孙建国. 2007. 复杂地表条件下地球物理场数值模拟方法评述. 世界地质, 26(3): 345–362.
孙建国, 蒋丽丽. 2009. 用于起伏地表条件下地球物理场数值模拟的正交曲网格生成技术. 石油地球物理勘探, 44(4): 494–500.
孙章庆, 孙建国, 张东良. 2009. 二维起伏地表条件下坐标变换法直流电场数值模拟. 吉林大学学报:地球科学版, 39(3): 528–534.
唐小平, 白超英. 2009a. 最短路径算法下三维层状介质中多次波追踪. 地球物理学报, 52(10): 2635–2643. DOI:10.3969/j.issn.0001-5733.2009.10.024
唐小平, 白超英. 2009b. 最短路径算法下二维层状介质中多次波追踪. 地球物理学进展, 24(6): 2087–2096. DOI:10.3969/j.issn.1004-2903.2009.06.022
王椿镛, 林中洋, 陈学波. 1995. 青海门源-福建宁德地学断面综合地球物理研究. 地球物理学报, 38(5): 590–598.
王椿镛, 张先康, 丁志峰, 等. 1997. 大别造山带上部地壳结构的有限差分层析成像. 地球物理学报, 40(1): 495–502.
徐涛, 徐果明, 高尔根, 等. 2004. 三维复杂介质的块状建模和试射射线追踪. 地球物理学报, 47(6): 1118–1126.
徐涛, 张忠杰, 田小波, 等. 2014. 长江中下游成矿带及邻区地壳速度结构:来自利辛-宜兴宽角地震资料的约束. 岩石学报, 30(4): 918–930.
赵烽帆, 马婷, 徐涛. 2014. 地震波初至走时的计算方法综述. 地球物理学进展, 29(3): 1102–1113. DOI:10.6038/pg20140313