地球物理学报  2019, Vol. 62 Issue (5): 1704-1715   PDF    
模型扩展和地形平化起伏地形处理方案在地震走时层析成像中的对比研究
郭高山1,2, 兰海强1,2, 陈凌1,2,3     
1. 中国科学院地质与地球物理研究所, 岩石圈演化国家重点实验室, 北京 100029;
2. 中国科学院大学, 北京 100049;
3. 中国科学院青藏高原地球科学卓越创新中心, 北京 100101
摘要:中国大陆中西部普遍具有强烈的地形起伏,起伏地形会对地震资料的处理分析产生严重干扰.精细处理起伏地形成为高精度地震成像的必然要求.传统方法通过填充低速介质将不规则模型扩展为规则模型来处理起伏地形.近年来,借助坐标变换将物理空间不规则模型转换为计算空间规则模型的地形平化方法,为解决起伏地形问题提供了新思路.本文基于经典的模型扩展和新发展的地形平化方法分别处理起伏地形,从走时正演、射线追踪和反演成像三个方面,全面细致地评判了两种地形处理方法在起伏地形层析成像中的适用性和有效性.结果表明,模型扩展中阶梯状近似和填充介质速度参与计算,会造成起伏地形走时计算精度损失,出现虚假射线路径和错误出射角,导致反演分辨率降低,成像结果模糊甚至失真;地形平化中采用贴体网格参数化,能够保证离散模型完全匹配起伏地形,并且保持起伏地形在物理空间和计算空间中均为自由表面.在此基础上发展的层析成像技术具有高度的保真性,有效地处理了地形起伏效应,为起伏地形区域精细速度成像提供了有力的技术保障.
关键词: 模型扩展      地形平化      起伏地形      层析成像     
A comparative study of both model expansion and topography flatten in topography handling in seismic traveltime tomography
GUO GaoShan1,2, LAN HaiQiang1,2, CHEN Ling1,2,3     
1. State Key Laboratory of Lithospheric Evolution, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
2. University of the Chinese Academy of Sciences, Beijing 100049, China;
3. Chinese Academy of Sciences Center for Excellence in Tibetan Plateau Earth Sciences, Beijing 100101, China
Abstract: Irregular surfaces are common situations in seismic exploration, especially at regions of middle and western China, making seismic processing and interpretation rather difficult in both oil/gas reservoir surveys and deep seismic soundings. The fine treatment of irregular topography becomes an inevitable requirement in high-precision seismic imaging. Traditionally, the problem of irregular topography is usually settled by model expansion scheme, in which the topography is regarded as an inner discontinuity after adding a low-velocity layer above the surface to expand an irregular physical space to a regular computational space. In recent years, a new method called topography flatten scheme, which involves the transformation between Cartesian and curvilinear coordinate systems, has been proposed to provide an alternative way to deal with rugged surfaces. In this paper, the both treatment schemes for irregular surfaces are briefly summarized. Then, a series of comparative tests are performed to comprehensively evaluate the accuracy, effectiveness and applicability of both methods from the first-arrival traveltimes, raypaths and imaging results of different synthetic models. These experiments indicate that the transition from a free surface to an inner discontinuity and the use of low-velocity media in model expansion scheme will cause many undesirable problems, such as loss of accuracy in calculating traveltimes and physically unrealistic seismic rays or spurious take-off angles, which can reduce seismic resolution and distort or smear structural images. For topography flatten scheme, the body-conforming grid can completely describe rugged surfaces and the topography-dependent eikonal equation (TDEE) is capable of exactly charactering the travel-time fields under complex surfaces, which establish a firm foundation for seismic imaging free from the effect of irregular surfaces. Therefore, we believe that topography flatten scheme has an exciting potential application in reconstructing high-quality seismic structures under complex topographical conditions.
Keywords: Model expansion scheme    Topography flatten scheme    Irregular topography    Seismic tomography    
0 引言

起伏地形是地球浅层油气矿产资源勘探(Reshef,1991孙建国,2007Farquharson et al., 2008Lelièvre et al., 2012Jiang and Zhang, 2016)及地球深部结构探测(Carbonell et al., 2004Zhang and Klemperer, 2010)等过程中普遍遇到的问题.起伏地形观测的地震资料存在信噪比低(Neuberg and Pointer, 2000)、干扰类型复杂(Bean et al., 2008)、静校正不准(Bevc,1997)等问题,给地震数据的处理和解释造成极大困难(Chávez-García et al., 1996Gray and Marfurt, 1995),一定程度上制约了地球内部精细结构探测的发展.走时层析成像是地震数据处理中最为简单且使用最为广泛的方法之一,利用它重建地球速度结构是地球内部结构探测和动力学研究的重要内容(Zhao et al., 1992, 1994Van Der Hilst et al., 1997Zhang and Thurber, 2005Rawlinson et al., 2006; Rawlinson and Fishwick, 2012Roecker et al., 2004, 2006Zhao et al., 2009Lv and Chen, 2017; Guo et al., 2019).发展能够有效处理起伏地形的层析成像方法,对于高精度的近地表速度重建和地球内部精细结构探测均具有重要意义.

根据走时计算的不同,走时层析成像方法可以分为射线类(Sambridge,1990Zhao et al., 1992Thurber and Eberhart-Phillips, 1999Xu et al., 2014俞贵平等,2017)和程函求解类(Hole,1992Rawlinson et al., 2006Sun et al., 2011孙章庆等,2012Lan and Zhang, 2013a, b).程函方程求解类(走时层析成像)方法在计算精度和效率方面具有显著优势,有效避免了射线类方法容易出现低速影区和焦散面的问题,自提出以来发展迅速,逐渐成为地球内部速度结构探测的主要方法之一(Hole et al., 2001Roecker et al., 2004, 2006郭飚等,2009Rawlinson et al., 2006; Rawlinson and Fishwick, 2012Li et al., 2009Huang and Bellefleur, 2011张风雪等,2014毛慧慧等,2016).经典程函方程求解类成像方法采用模型扩展方法直接处理起伏地形下程函方程的求解问题(Hole,1992Taillandier et al., 2009Huang and Bellefleur, 2012李德勇等,2016),其核心思想是将起伏地形当作介质内部不连续面,通过填充低速介质将物理空间地形起伏的不规则模型扩展为计算空间的规则模型来处理.然而,该处理过程中采用的起伏地形的阶梯状表示及虚拟填充介质可能导致计算精度的损失(在后续内容中将进一步探讨和证实).地形平化是近年逐渐发展起来处理起伏地形的新方法,借助贴体网格和坐标变换,将物理空间的不规则模型映射为计算空间的规则模型,通过数值求解曲线坐标系中与地形相关的程函方程(详细说明见附录)精确获取起伏地形的走时场和射线路径,在此基础上实现起伏地形下的地震层析成像(Ma and Zhang, 2015Zhang et al., 2018).然而,两种起伏地形的处理方法到底哪种效果更优?截至目前,尚无定量的结论.针对这一问题,本文对上述两种处理起伏地形的方法,从层析成像中的三个环节,即走时计算、射线追踪和层析成像,进行全面系统的测试分析、对比评判,研究他们对起伏地形处理的有效性和适用性,为起伏地形地震资料高精度速度结构成像的方法选择提供借鉴依据和实践基础.

1 方法介绍

起伏地形模型参数化作为走时计算和反演的基础,是获得起伏地形条件下高精度成像结果的关键步骤(Zelt and Smith, 1992).目前,适用于起伏地形的建模方式以非结构网格为主(Sethian and Vladimirsky, 2000Zhang and Thurber, 2005Qian et al., 2007Kao et al., 2008Bai et al., 2011Lelièvre et al., 2012Zhang et al., 2017),然而从计算量和算法复杂度来看,基于结构网格的方法更加简洁便利,计算效率更高.发展基于结构网格的起伏地形地震数据处理方法对日益增长的地震探测数据的高效处理显得尤为重要.如前文所述,基于结构网格的起伏地形处理方法主要有两种,包括经典的模型扩展和新发展的地形平化.如图 1a所示,模型扩展将地形最高点作为模型扩展边界,将起伏地形视为介质内部界面,在起伏地形以上的空间填充低速介质,将不规则的模型扩展为规则模型,进而用规则矩形网格来离散扩充所获得的规则模型.这种方法的误差来源于规则网格离散导致的对起伏地形的阶梯状近似;此外,该方法需要对起伏地形之上的空间填充低速介质,这会产生虚假的射线路径和错误的出射角等(将在后文中证实)问题.不同于模型扩展,地形平化使用形状各异、大小不等的四边形结构化贴体网格离散模型(如图 1b所示),通过坐标变换将笛卡尔坐标下的不规则四边形网格映射为曲线坐标系下的规则矩形网格,此时笛卡尔坐标系下的经典程函方程变换为曲线坐标系的与地形相关的程函方程(详细介绍见附录).由于在物理空间和计算空间中都保持起伏地形为自由表面,理论上可以完备、精确地描述复杂的地质模型同时兼有规则结构化网格求解高效的特点(Lan and Zhang, 20112013bb).

图 1 起伏地形模型参数化方案 (a)模型扩展策略(Hole,1992);(b)地形平化策略(Lan and Zhang, 2013ab).浅灰色部分(V1)表示物理模型,深灰色部分(V0)表示低速填充介质. Fig. 1 Schematic diagrams illustrating two different irregular surface treatments (a) Model expansion scheme (Hole, 1992); (b) Topography flatten scheme (Lan and Zhang, 2013ab). The light and dark gray denote the physical model (V1) and the low-velocity layer (V0), respectively.

基于贴体网格参数化,兰海强等分别发展了求解各向同性与各向异性程函方程的快速扫描法(Lan et al., 2014; Lan and Chen, 2018),结合拟牛顿L-BFGS反演算法,开展了起伏地形下地震层析成像研究.在本文的对比研究中,我们采用兰海强等人发展的走时计算和反演迭代方法,保持正反演所用的方法和其他参数相同,只是处理起伏地形的方式不同,通过具体实例进行详细的分析测试、对比评判,给读者定量直观的认识,为地震学成像中起伏地形的处理提供借鉴依据.

2 数值试验

下面我们通过数值实例,分别从走时计算、射线追踪、反演成像三个方面来定量评判两种起伏地形处理模式的精度及有效性.

2.1 起伏地形条件下初至波走时计算和射线追踪

我们设计了两组数值模型来评判基于模型扩展和地形平化的起伏地形初至波走时计算和射线追踪效果.第一组为斜坡地形均匀介质模型,用以考察两种方法处理起伏地形的精度;第二组为不同起伏地形的连续介质模型,用以测试两种方法处理不同起伏地形的有效性.

2.1.1 精度测试

第一组模型的起伏地形为斜坡地形,为方便与真实走时和射线路径对比,模型取速度为4 km·s-1的均匀介质,模型扩展过程采取广泛使用的空气层填充介质(0.34 km·s-1)(Hole,1992Taillandier et al., 2009Huang and Bellefleur, 2012李德勇等,2016Ryberg and Haberland, 2018).分别使用相同网格间距(0.05 km)的贴体网格和矩形网格离散模型.使用高阶Lax-Friedrichs快速扫描法(Lan and Zhang, 2013b)求解程函方程模拟初至波的走时场,沿走时梯度最速下降方向,从检波器逆向追踪到炮点,获得初至波的射线路径.从斜坡地形走时可见,模型扩展获得的地形走时呈现为锯齿状,计算走时与精确走时有明显偏差且一般较精确解大(图 2a);地形平化获得的地形走时连续光滑,与精确解吻合得很好(图 2b).在走时场正演的基础上计算射线路径,为清楚起见,将起伏地形的射线路径放大显示.图 3a表示模型扩展的射线路径,与解析解对比可见,填充低速介质后斜坡地形射线出射角畸变,模拟射线严重偏离真实路径.图 3b表示地形平化的射线路径,对比结果显示模拟射线与真实射线结果一致.可以初步断定:模型扩展方式计算斜坡模型的走时误差较大,射线路径与真实路径存在明显偏差;而地形平化可以精确计算走时和射线路径.

图 2 斜坡地形均匀介质模型地形走时数值解与解析解对比 (a)模型扩展;(b)地形平化.黑色点代表走时数值解,灰色点代表走时解析解. Fig. 2 Comparison of numerical and analytical traveltimes at the topography of homogeneous model with a slope surface (a) The results obtained by model expansion scheme; (b) The results obtained by topography flatten scheme. The black and gray dotted lines denote numerical and analytical traveltimes, respectively.
图 3 斜坡地形均匀介质模型计算射线与真实射线对比 (a)模型扩展;(b)地形平化.黑色和白色实线分别表示计算射线和真实射线.图(a)中深灰色部分表示物理模型,浅灰色部分表示填充介质. Fig. 3 Comparison of numerical and analytical raypaths for the homogeneous model with a slope surface (a) The results obtained by model expansion scheme; (b) The results obtained by topography flatten scheme. The black and white solid lines denote numerical and analytical raypaths, respectively. The dark gray media indicate physical model and light gray media show low-velocity layer above the slope surface.
2.1.2 有效性测试

进一步测试两种起伏地形处理方法对山峰、凹陷、台阶及山峰和凹陷组合的起伏地形模型的有效性,模型均取速度随深度均匀增加(2~5 km·s-1)的连续介质.图 47表示不同模型的计算走时场和射线路径,同样放大显示起伏地形的结果.如图所示,模型扩展中低速介质的填充导致起伏地形处出现虚假射线和错误出射角;地形平化有效消除了不同起伏地形对模型内部计算产生的影响,射线路径在连续介质中表现为弧形且垂直于走时等值线,在地形起伏较大处射线路径会有轻微的曲折,这与网格疏密及剖分的有效性有关.分析可知,模型扩展将起伏地形处理为模型内部速度间断面,使用矩形网格剖分模型,对起伏地形进行阶梯状近似,为满足有限差分的计算要求,走时计算和射线追踪需要使用填充介质的速度值,从而引发上述问题;地形平化的优势在于贴体网格能够精确描述地形的不同起伏形态,差分格式中的节点全部为实际模型离散点,通过求解与地形相关的程函方程可以精确模拟起伏地形模型的走时场和射线路径.上述模型实例证明:针对不同类型的地形,模型扩展难以消除起伏地形对走时和射线路径计算的影响,而地形平化均能够有效地消除起伏地形效应.

图 4 山丘地形连续介质模型的地震波走时场和射线路径 (a)模型扩展;(b)地形平化.其中白色等值线表示走时场,黑色实线表示射线路径. (a)图中灰色部分表示填充的低速介质. Fig. 4 Seismic traveltime contours and raypaths for the constant-gradient velocity model with a hill surface (a) The results obtained by model expansion scheme; (b) The results obtained by topography flatten scheme. White isolines and black solid lines show traveltime contours and raypaths, respectively. The dark gray media in figure (a) denote a low-velocity layer above the hill surface.
图 5 凹陷地形连续介质模型地震波走时场和射线路径.图例同图 4 Fig. 5 Seismic traveltime contours and raypaths for the constant-gradient velocity model with a valley surface. Same legend as Fig. 4
图 6 台阶地形连续介质模型地震波走时场和射线路径.图例同图 4 Fig. 6 Seismic traveltime contours and raypaths for the constant-gradient velocity model with a step surface. Same legend as Fig. 4
图 7 复杂地形连续介质模型地震波走时场和射线路径.图例同图 4 Fig. 7 Seismic traveltime contours and raypaths for the constant-gradient velocity model with a wavy surface. Same legend as Fig. 4
2.2 起伏地形条件下初至波层析成像

在实际数据处理工作中,为方便起见,可能忽略地形起伏的影响,而把起伏地形采集的数据用水平地形模型进行处理,这种水平地形假设在地震层析成像中的合理性尚待商榷.因此,我们首先评价水平地形假设导致的走时误差,进而考虑对反演结果的影响程度,探讨层析成像过程中处理地形的必要性.然后,以复杂地形模型为例评判不同地形处理方法在地震层析成像过程中的处理效果.

2.2.1 水平地形假设误差评价

以平缓起伏地形模型为例,评估水平地形假设对走时计算的影响程度.设置模型如图 8所示,地形起伏幅度为-0.15~0.2 km,模型深度为4.0 km,最大起伏程度仅为模型深度的5%,接近水平地形模型.该模型的背景速度为2~5 km·s-1的连续介质,存在大小3 km×1 km的高、低速度异常体,速度异常幅度±0.5 km·s-1.首先以如图 8a所示的3个炮点和位于起伏地形的4个检波器记录的走时为例,使用地形平化方法处理起伏地形,利用迎风格式快速扫描法(Lan and Chen, 2018; Lan et al., 2018)计算走时,分析水平地形假设造成的走时误差(表 1).可以看到,不考虑地形因素会造成明显的走时误差,该误差会对成像过程中异常体的恢复产生影响(这点在下文得到验证).

图 8 水平地形假设对异常体模型的成像影响测试 (a)真实模型;(b)异常体扰动;(c)基于水平地形假设的反演结果;(d)基于地形平化处理的反演结果.其中五角星表示震源,三角形表示检波器,用于标注表 1中的观测系统.①、②两条灰色虚线表示所抽取速度剖面的位置,将用于下图展示地形起伏对成像结果的具体影响. Fig. 8 The effect of the flat-surface approximation on anomalous-body inversion (a) The true model; (b) The anomalous-body perturbation; (c) The imaging result based on flat-surface approximation; (d) The imaging result based on topography flatten scheme. The stars mark the sources, and triangles indicate the receivers in Table 1. The two grey lines (① and ②) are used to mark the positions of velocity profiles, which will be employed to show the influence of topography on inversion in the next figure.
表 1 水平地形假设导致的走时百分比误差(%) Table 1 The relative percent error of the traveltimes induced by flat-surface approximation

接下来,测试水平地形假设对成像结果的影响,观测系统为均匀分布在地形的39个炮点和199个检波器.分别以起伏地形和水平地形背景速度模型为初始模型,用L-BFGS反演算法更新速度模型,异常体恢复结果如图 8c8d所示.为定量对比,抽取地形附近和通过异常体中心的速度曲线(图 9).对比可知,基于水平地形假设,忽略地形影响的成像结果中异常体的形态和数值恢复产生较大误差,地形附近出现明显的虚假异常;考虑地形起伏因素的成像结果中,异常体边界和数值得到有效恢复.可以推断,即使在地形起伏较平缓的情况下,水平地形假设也会产生明显的虚假异常,降低速度反演的分辨率.因此,地震层析中考虑地形起伏效应是高精度成像的必然需求.

图 9 不同深度的速度曲线对比 (a)和(b)分别对应图 8b中的速度剖面①和②. Fig. 9 Comparison of the velocity logs extracted from different depths (a) and (b) correspond to the positions of grey lines ① and ② in Fig. 8b, respectively.
2.2.2 地形处理方法在地震层析成像中的效果评价

通过起伏地形条件下的简单异常体、梯形速度界面和标准检测板的反演,对比分析地形平化和模型扩展在起伏地形层析成像过程中的有效性.沿用2.1中设计的复杂地形模型,观测系统仍为39个炮点和199个检波器均匀布置于起伏地形.首先,仍以图 10所示的简单异常体反演为例,以连续介质背景速度模型为初始模型,图 10c10d分别为采用模型扩展和地形平化的反演结果,图 11为穿过异常体中心不同方向的速度曲线.可观察到在起伏程度剧烈的情况下,低速介质填充导致反演结果中地形附近人为高速异常,异常体周围产生的明显拖尾等干扰影响异常体边界的恢复;地形平化仍能获取可靠的信息,对近地表速度恢复准确,异常体边界刻画清晰.然后,设计如图 12a所示的梯形起伏地层模型,速度随深度均匀增加(2~5.5 km·s-1),在梯形起伏间断面发生0.5 km·s-1的速度跳跃.初始模型(图 12b)与真实模型速度分布相同,但速度间断面水平,以此观察地形处理方法对速度界面反演的影响.反演结果如图 12c12d所示,从中可见模型扩展的反演界面形态难以得到有效恢复,而地形平化反演的梯形速度界面形态清晰,两者的数值误差在穿过梯形地层的速度曲线中体现得更为明显(图 13).最后,如图 14所示,以检测板形式的速度扰动恢复为例,速度异常仍为大小1.33 km×1.33 km,振幅0.3 km·s-1的高速、低速相间分布的模式,用以观测地形处理方法对分辨率的影响.反演的初始模型为连续介质背景速度模型,图 14c14d分别为基于模型扩展和地形平化起伏地形处理方案的反演结果,图 15为不同方向的速度曲线.从图中可以观察到,模型扩展处理方案的成像结果中低速填充介质对检测板异常幅值和形状的恢复产生明显干扰,导致浅部成像结果模糊,分辨率降低,深部射线稀疏且交叉较少使得深部检测板恢复较差;地形平化中检测板的浅部异常形态及幅值均恢复得较好.从上述数值反演实例可知,模型扩展过程中的低速填充介质对近地表速度的恢复和形态的刻画有较大影响;地形平化能够完全消除层析成像中的地形效应,对于速度异常具有较高的保真性.

图 10 复杂地形异常体模型反演结果 (a)真实模型;(b)异常体扰动;(c)模型扩展反演结果;(d)地形平化反演结果. ①、②、③灰色虚线表示所抽取速度剖面的位置,将用于下图中展示两种方法的具体差别. Fig. 10 Inversion results for the anomalous-body model with a wavy surface (a) The true model; (b) The anomalous-body perturbation; (c) The imaging result obtained by model expansion scheme; (d) The imaging result obtained by topography flatten scheme. The three grey dotted lines (①, ②, ③) are used to mark the positions of velocity profiles, which will be employed to show the specific difference between the both methods in the next figure.
图 11 通过异常体中心不同方向的速度曲线对比 (a)、(b)和(c)分别对应图 10b中的速度剖面①、②和③. Fig. 11 Comparison of the velocity logs extracted from different directions through the center of anomalous bodies (a), (b) and (c) correspond to the positions of grey dotted lines ①, ② and ③ in Fig. 10b, respectively.
图 12 复杂地形梯形地层模型反演结果 (a)真实模型;(b)初始模型;(c)模型扩展反演结果;(d)地形平化反演结果.图例同图 10. Fig. 12 Inversion results for the trapezoidal-interface model with a wavy surface (a) The theoretical model; (b) The initial model; (c) The imaging result obtained by model expansion scheme; (d) The imaging result obtained by topography flatten scheme. Same legend as Fig. 10.
图 13 通过梯形地层的速度曲线对比.该剖面对应图 12a中的速度剖面.图例同图 11 Fig. 13 Comparison of the velocity logs extracted through the trapezoidal interface. The velocity profile corresponds to the position of grey dotted line in Fig. 12a. Same legend as Fig. 11
图 14 复杂地形检测板模型反演结果 (a)真实模型;(b)检测板扰动;(c)模型扩展的反演结果;(d)地形平化的反演结果.图例同图 10. Fig. 14 Inversion results for the checkboard model with a wavy surface (a) The true model; (b) The checkboard perturbation; (c) The imaging result obtained by model expansion scheme; (d) The imaging result obtained by topography flatten scheme. Same legend as Fig. 10.
图 15 通过检测板不同方向的速度曲线对比 (a)、(b)和(c)分别对应图 14b中的速度剖面①、②和③.图示同图 11. Fig. 15 Comparison of the velocity logs extracted from different directions through the checkerboard (a), (b) and (c) correspond to the positions of grey dotted lines ①, ② and ③ in Fig. 14b, respectively. Same legend as Fig. 11.
3 结论与讨论

起伏地形条件下走时层析成像方法的研究为起伏地形区域精细速度结构反演提供理论基础和技术保障.本文从走时正演、射线路径追踪和层析成像三个方面,对模型扩展和地形平化两种(地震层析成像中)起伏地形处理方案进行了详细的分析评判,得到以下两点认识:

(1) 模型扩展将起伏地形处理成模型内部不连续性界面,起伏地形的阶梯状表示和填充介质速度参与计算是走时误差的主要来源,地形走时精度的损失导致虚假射线路径和错误出射角,使得近地表成像中出现虚假异常,导致反演分辨率降低,结果失真.

(2) 地形平化采用的贴体网格离散保证离散模型能够完全匹配起伏地形,将物理空间中的起伏地形映射为计算空间的水平地形,在此基础上发展的层析成像技术具有高度的保真性,有效地处理了地形起伏效应,为起伏地形区域精细速度结构成像提供了新的处理思路.

目前本研究仅限于利用初至波(单震相)进行近地表速度结构成像,对深部结构的成像约束稍显不足.下一步拟将该方法拓展至后续震相(多震相),并将起伏地形处理的思路和方法拓展用于深部起伏界面的处理,发展地球内部结构的高精度成像方法,为深部探测研究提供便利.

附录  经典程函方程和与地形相关的程函方程

笛卡尔坐标系各向同性介质中的经典程函方程为

(A1)

上述方程描述了射线经过慢度为s(x, z)的介质中的点(x, z)时的走时T(x, z),模型扩展方案通过求解该方程可以获得规则填充模型的走时场.

曲线坐标系各向同性介质的程函方程为

(A2)

上述方程可由笛卡尔坐标系的经典程函方程推导变换得到(详见Lan and Zhang, 2013a).方程中这些参数表征了起伏地形的局部特征.该方程可以描述任意地形起伏模型中地震波传播的走时规律,被称之为与地形相关的程函方程(Topography-dependent eikonal equation,TDEE).这里xr表示意义类似,J=xqzrxrzq.当地形水平时,xrzq退化为0,xqzr退化为1,此时,J为1,上述程函方程简化为经典的程函方程,可见经典程函方程为与地形相关程函方程在水平地形下的特殊形式.对于该方程的求解,具体可详见Lan and Zhang(2013a, b),Lan and Chen(2018).

References
Bai C Y, Li X L, Tang X P. 2011. Seismic wavefront evolution of multiply reflected, transmitted, and converted phases in 2D/3D triangular cell model. Journal of Seismology, 15(4): 637-652. DOI:10.1007/s10950-011-9242-y
Bean C, Lokmer I, O'Brien G. 2008. Influence of near-surface volcanic structure on long-period seismic signals and on moment tensor inversions:Simulated examples from Mount Etna. Journal of Geophysical Research:Solid Earth, 113(B8): B08308. DOI:10.1029/2007JB005468
Bevc D. 1997. Flooding the topography:Wave-equation datuming of land data with rugged acquisition topography. Geophysics, 62(5): 1558-1569. DOI:10.1190/1.1444258
Carbonell R, Simancas F, Juhlin C, et al. 2004. Geophysical evidence of a mantle derived intrusion in SW Iberia. Geophysical Research Letters, 31(11): L11601. DOI:10.1029/2004GL019684
Chávez-García F J, Sánchez L R, Hatzfeld D. 1996. Topographic site effects and HVSR.A comparison between observations and theory.. Bulletin of the Seismological Society of America, 86(5): 1559-1573.
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.
Gray S H, Marfurt K J. 1995. Migration from topography:Improving the near-surface image. Canadian Journal of Exploration Geophysics, 31(1-2): 18-24.
Guo B, Liu Q Y, Chen J H, et al. 2009. Teleseismic P-wave tomography of the crust and upper mantle in Longmenshan area, West Sichuan. Chinese Journal of Geophysics (in Chinese), 52(2): 346-355.
Guo B, Chen J H, Liu Q Y, et al. 2019. Crustal structure beneath the Qilian Orogen Zone from multiscale seismic tomography. Earth and Planetary Physics, 3(2): 1-11.
Hole J A. 1992. Nonlinear high-resolution three-dimensional seismic travel time tomography. Journal of Geophysical Research:Solid Earth, 97(B5): 6553-6562. DOI:10.1029/92JB00235
Hole J A, Catchings R D, St Clair K C, et al. 2001. Steep-dip seismic imaging of the shallow san andreas fault near parkfield. Science, 294(5546): 1513-1515. DOI:10.1126/science.1065100
Huang J W, Bellefleur G. 2011. Joint transmission and reflection traveltime tomography in imaging permafrost and thermokarst lakes in Northwest Territories, Canada. Leading Edge, 30(11): 1304-1312. DOI:10.1190/1.3663404
Huang J W, Bellefleur G. 2012. Joint transmission and reflection traveltime tomography using the fast sweeping method and the adjoint-state technique. Geophysical Journal International, 188(2): 570-582. DOI:10.1111/gji.2012.188.issue-2
Jiang W B, Zhang J. 2016. First-arrival traveltime tomography with modified total-variation regularization. Geophysical Prospecting, 65(5): 1138-1154.
Kao C Y, Osher S, Qian J L. 2008. Legendre-transform-based fast sweeping methods for static Hamilton-Jacobi equations on triangulated meshes. Journal of Computational Physics, 227(24): 10209-10225. DOI:10.1016/j.jcp.2008.08.016
Lan H Q, Zhang Z J. 2011. 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, 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 traveltimes with an irregular surface. Bulletin of the Seismological Society of America, 103(3): 2070-2082. DOI:10.1785/0120120199
Lan H Q, Chen J Y, Zhang Z J. 2014. A fast sweeping scheme for calculating p wave first-arrival travel times in transversely isotropic media with an irregular surface. Pure and Applied Geophysics, 171(9): 2199-2208. DOI:10.1007/s00024-014-0836-5
Lan H Q, Chen L. 2018. An upwind fast sweeping scheme for calculating seismic wave first-arrival travel times for models with an irregular free surface. Geophysical Prospecting, 66(2): 327-341. DOI:10.1111/gpr.2018.66.issue-2
Lan H Q, Chen L, Badal J. 2018. A hybrid method for calculating seismic wave first-arrival traveltimes in two-dimensional models with an irregular surface. Journal of Applied Geophysics, 155: 70-77. DOI:10.1016/j.jappgeo.2018.05.011
Lelièvre P G, Farquharson C G, Hurich C A. 2012. Joint inversion of seismic traveltimes and gravity data on unstructured grids with application to mineral exploration. Geophysics, 77(1): K1-K15.
Li Y D, Dong L G, Liu Y Z. 2016. First arrival traveltime tomography based on an improved scattering-integral algorithm. Chinese Journal of Geophysics (in Chinese), 59(10): 3820-3828. DOI:10.6038/cjg20161026
Li Z W, Roecker S, Li Z H, et al. 2009. Tomographic image of the crust and upper mantle beneath the western Tien Shan from the MANAS broadband deployment:Possible evidence for lithospheric delamination. Tectonophysics, 477(1-2): 49-57. DOI:10.1016/j.tecto.2009.05.007
Lv Y, Chen L. 2017. Upper crustal P-wave velocity structure beneath two volcanic areas in northern Iran. Science China:Earth Sciences, 60(4): 786-795. DOI:10.1007/s11430-016-9005-7
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
Mao H H, Lei J S, Teng J W. 2016. Teleseismic P-wave tomography of the upper mantle along the north-south profile under the northern Ordos basin. Chinese Journal of Geophysics (in Chinese), 59(6): 2056-2065. DOI:10.6038/cjg20160612
Neuberg J, Pointer T. 2010. Effects of volcano topography on seismic broad-band waveforms. Geophysical Journal International, 143(1): 239-248.
Qian J L, Zhang Y T, Zhao H K. 2007. Fast sweeping methods for eikonal equations on triangular meshes. Society for Industrial and Applied Mathematics, 45(1): 83-107.
Rawlinson N, Reading A M, Kennett B L N. 2006. Lithospheric structure of Tasmania from a novel form of teleseismic tomography. Journal of Geophysical Research:Solid Earth, 111(B2): B02301. DOI:10.1029/2005JB003803
Rawlinson N, Fishwick S. 2012. Seismic structure of the southeast Australian lithosphere from surface and body wave tomography. Tectonophysics, 572-573: 111-122. DOI:10.1016/j.tecto.2011.11.016
Reshef M. 1991. Depth migration from irregular surfaces with depth extrapolation methods. Geophysics, 56(1): 119-122. DOI:10.1190/1.1442947
Roecker S, Thurber C, Mcphee D. 2004. Joint inversion of gravity and arrival time data from Parkfield:New constraints on structure and hypocenter locations near the SAFOD drill site. Geophysical Research Letters, 31(12): L12S04. DOI:10.1029/2003GL019396
Roecker S, Thurber C, Roberts K, et al. 2006. Refining the image of the San Andreas Fault near Parkfield, California using a finite difference travel time computation technique. Tectonophysics, 426(1-2): 189-205. DOI:10.1016/j.tecto.2006.02.026
Ryberg T, Haberland C. 2018. Bayesian inversion of refraction seismic traveltime data. Geophysical Journal International, 212(3): 1645-1656. DOI:10.1093/gji/ggx500
Sambridge M S. 1990. Non-linear arrival time inversion:constraining velocity anomalies by seeking smooth models in 3-D. Geophysical Journal International, 102(3): 653-677. DOI:10.1111/gji.1990.102.issue-3
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
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, Sun Z Q, Han F X. 2011. A finite difference scheme for solving the eikonal equation including surface topography. Geophysics, 76(4): T53-T63. DOI:10.1190/1.3580634
Sun Z Q, Sun J G, Han F X. 2012. Traveltime computation using the upwind finite difference method with nonuniform grid spacing in a 3D undulating surface condition. Chinese Journal of Geophysics (in Chinese), 55(7): 2441-2449. DOI:10.6038/j.issn.0001-5733.2012.07.028
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
Thurber C, Eberhart-Phillips D. 1999. Local earthquake tomography with flexible gridding. Computers & Geosciences, 25(7): 809-818.
Van Der Hilst R D, Widiyantoro S, Engdahl E R. 1997. Evidence for deep mantle circulation from global tomography. Nature, 386(6625): 578-584. DOI:10.1038/386578a0
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
Yu G P, Xu T, Zhang M H, et al. 2017. Nonlinear travel-time inversion for 3-D complex crustal velocity structure. Chinese Journal of Geophysics (in Chinese), 60(4): 1398-1410. DOI:10.6038/cjg20170414
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
Zhang F X, Wu Q J, Li Y H. 2014. A traveltime tomography study by teleseismic S wave data in the Northeast China area. Chinese Journal of Geophysics (in Chinese), 57(1): 88-101. DOI:10.6038/cjg20140109
Zhang H J, Thurber C. 2005. Adaptive mesh seismic tomography based on tetrahedral and Voronoi diagrams:Application to Parkfield, California. Journal of Geophysical Research:Solid Earth, 110(B4): B04303. DOI:10.1029/2004JB003186
Zhang M H, Xu T, Bai Z M, et al. 2017. Ray tracing of turning wave in elliptically anisotropic media with an irregular surface. Earthquake Science, 30(5-6): 219-228. DOI:10.1007/s11589-017-0192-5
Zhang X Y, Bai Z M, Xu T, et al. 2018. Joint tomographic inversion of first-arrival and reflection traveltimes for recovering 2-D seismic velocity structure with an irregular free surface. Earth and Planetary Physics, 2(3): 1-11.
Zhang Z J, Klemperer S. 2010. Crustal structure of the Tethyan Himalaya, southern Tibet:new constraints from old wide-angle seismic data. Geophysical Journal International, 181(3): 1247-1260.
Zhao D P, Hasegawa A, Horiuchi S. 1992. Tomographic imaging of P and S wave velocity structure beneath northeastern Japan. Journal of Geophysical Research:Solid Earth, 97(B13): 19909-19928. DOI:10.1029/92JB00603
Zhao D P, Hasegawa A, Kanamori H. 1994. Deep structure of Japan subduction zone as derived from local, regional, and teleseismic events. Journal of Geophysical Research:Solid Earth, 99(B11): 22313-22329. DOI:10.1029/94JB01149
Zhao L, Allen R M, Zheng T Y, et al. 2009. Reactivation of an Archean craton:Constraints from P-and S-wave tomography in North China. Geophysical Research Letters, 36(17): L17306. DOI:10.1029/2009GL039781
郭飚, 刘启元, 陈九辉, 等. 2009. 川西龙门山及邻区地壳上地幔远震P波层析成像. 地球物理学报, 52(2): 346-355.
李勇德, 董良国, 刘玉柱. 2016. 基于改进的散射积分算法的初至波走时层析. 地球物理学报, 59(10): 3820-3828. DOI:10.6038/cjg20161026
毛慧慧, 雷建设, 滕吉文. 2016. 鄂尔多斯盆地北缘南北向剖面上地幔远震P波层析成像. 地球物理学报, 59(6): 2056-2065. DOI:10.6038/cjg20160612
孙建国. 2007. 复杂地表条件下地球物理场数值模拟方法评述. 世界地质, 26(3): 345-362. DOI:10.3969/j.issn.1004-5589.2007.03.015
孙章庆, 孙建国, 韩复兴. 2012. 三维起伏地表条件下地震波走时计算的不等距迎风差分法. 地球物理学报, 55(7): 2441-2449. DOI:10.6038/j.issn.0001-5733.2012.07.028
俞贵平, 徐涛, 张明辉, 等. 2017. 三维复杂地壳结构非线性走时反演. 地球物理学报, 60(4): 1398-1410. DOI:10.6038/cjg20170414
张风雪, 吴庆举, 李永华. 2014. 中国东北地区远震S波走时层析成像研究. 地球物理学报, 57(1): 88-101. DOI:10.6038/cjg20140109