基于三维模型中的射线追踪方法测定鲁甸M_s6.5地震震源深度

2014年8月3日云南鲁甸发生M_S6.5地震,造成重大人员伤亡和财产损失。此次地震发震构造复杂,引起了地震学界的广泛关注。本文基于三维建模方法,建立鲁甸地区三维层状非均匀速度模型,采用逐段迭代射线追踪方法对鲁甸M_S6.5地震进行三维射线追踪走时计算;利用震中距150km范围内的近震Pg波震相走时数据,通过射线追踪走时拟合获得鲁甸地震的震源深度约为12km,与前人研究成果基本一致,表明了本文采用的三维建模和射线追踪方法的有效性。     

三维地质模型中地震波共轭梯度非线性走时反演

地震体波走时层析成像是探测地球内部速度结构的重要方法之一。基于三维块状建模以及三角形拼接的界面描述方式,结合快速高效的逐段迭代射线追踪方法,获得三维复杂地质模型中的地震射线路径与走时信息,采用共轭梯度非线性反演算法,进行地震波走时反演。实验结果表明共轭梯度反演算法在三维层状模型中具有较高的有效性。     

基于图形结构的三维射线追踪方法

在地震层析成像研究中,为了克服最小走时射线路径追踪方法存在的问题,对该方法计算过程中的关键步骤进行了改进．在节点走时的计算中引入Ｂｒｅｓｅｎｈａｍ画线算法;在最小走时节点查寻中,结合使用快速排序算法与插入排序算法,替代以往方法中多采用的堆排序算法;所采用的节点设置方式,可以引入速度界面,还可以实现反射波射线追踪．模型计算证明,改进的最小走时射线路径方法具有精度高,速度快的特点,所提出的三维空间反射波射线追踪算法简便易行。     

基于LTI和网格界面剖分的三维地震射线追踪算法

将二维线性走时插值射线追踪算法(LTI)推广应用至三维模型,并结合网格界面剖分方式,提出了一种三维射线追踪算法．该算法既可获得高精度的全局最小走时和射线路径,又具有快速稳定的特点．三维模型计算结果表明,在模型参数包括网格密度完全相同情况下,本文算法较传统的三维最短路径算法在计算效率、走时和射线的计算精度上均有明显改进．     

地震波走时和射线的有限差分计算

以往都是采用射线追踪的方法计算地震波的走时和射线，但是当速度模型复杂时这种方法存在一些问题。本文提出另一种计算地震波走时和射线的方法。该方法从程函方程出发，利用互换原理和Ｆｅｒｍａｔ原理计算出各种波的到时和射线。解决了射线追踪方法存在的问题。为地震波走时和射线的计算以及地震波走时反演开辟了一条新途径。     

三维复杂地形近地表速度估算及地震层析静校正

在地表一致性模型的基础上提出一种可适用于宽线剖面、弯曲测线、传统的二维和目前广泛使用的三维地震观测.在地形及近地表低降速带地质结构复杂的探区，低降速带厚度及速度估算的精度是静校正处理的关键.本研究根据三维地震观测的初至走时数据，利用最小平方与QR分解相结合的算法，在三维空间重建近地表低降速带速度模型，根据重建速度模型实现了静校正长波长分量与短波长分量的同步计算.分析了复杂的近地表低降速带模型初至波的性质，在观测值的自动拾取以及理论值的计算中充分考虑了可能成为初至波的直达波、折射波和反射波的利用，提高了低降速带速度模型反演的精度.在初至走时观测数据的拾取中，本研究采用分形算法克服了初至波波形差异以及折射波相位反转导致的拾取误差，实现了三维初至拾取的大规模全自动化运算.在射线路径与初至波理论走时的计算中，本研究采用一种计算量与模型复杂程度无关的三维射线追踪方法，该方法以最小走时射线路径保证了与观测数据有同等意义的初至波的射线追踪及理论走时的计算.野外实际资料的处理结果表明了方法的有效性.     

运用改进的地震定位交切法重定位2001—2010年云南地震

传统地震定位方法利用震源轨迹确定震源位置,但基于均匀或横向均匀介质模型必然导致定位误差。为此对传统方法进行改进,发展适用于三维复杂地壳速度模型的地震定位交切法。利用最小走时树射线追踪技术,以离散方式准确计算三维复杂地壳速度模型中的震源轨迹,将震源定位于震源轨迹交汇的密集点。将该方法应用于云南地区地震重定位,得到较高定位精度。     

基于三维层状介质模型地震层析成像正演网格研究

地震层析成像技术是研究地层速度结构的一个重要手段,该技术是利用震源到接收点的地震波走时或波形的观测数据,结合建立的数学物理模型来重建地层速度结构。而层析成像正演模拟的精度将直接影响反演出的速度与地层真实速度的拟合度,因此找出一种精度较高的层析成像正演模拟算法是非常重要的。本文针对三维层状介质,通过网格化模型,采用最小走时射线追踪方法展开层析成像正演研究,分别采用三角形网格、矩形网格和六边形网格进行模型参数化,并对比分析各自正演结果的精度,总结出针对三维层状介质模型的最佳网格划分方式。      

青藏高原东部及其边缘地区的地壳上地幔三维速度结构

利用在青藏高原东部及其邻近地区记录到的１万余条近震到时资料,反演该地区的地壳上地幔三维速度结构。采用网格点模型描述三维速度结构,模型维数为２２２２６,网格点间距水平向为１００ｋｍ,垂直向为２０ｋｍ,网格点之间的速度值通过线性插值给出。采用改进了的快速三维射线追踪方法,确定三维非均匀介质中的地震射线路径和理论走时。反演结果显示,青藏高原南部的上地壳中（３０ｋｍ左右的深度）存在一低速区,这和面波反演的结果一致,羌塘块体下地壳有明显的低速异常带,青藏公路沿线的垂直速度剖面显示出岩石层受挤压增厚的构造特征。     

利用三维弯曲界面的反演方法重建唐山震区莫霍界面三维构造形态

研究利用反射波的走时反演三维弯曲界面和介质层速度的计算方法．各层界面利用分段非完全三次多项式描述．正问题采用一种快速的三维射线追踪方法．反演过程采用变阻尼最小二乘法．数值模拟结果表明，解很快收敛到真模型．处理了通过唐山震区的实测资料，重建了该震区莫霍界面三维构造形态，并揭示了该区域构造与地震活动的关系．      

Seismic traveltime inversion of 3D velocity model with triangulated interfaces

Seismic traveltime tomographic inversion has played an important role in detecting the internal structure of the solid earth. We use a set of blocks to approximate geologically complex media that cannot be well described by layered models or cells. The geological body is described as an aggregate of arbitrarily shaped blocks, which are separated by triangulated interfaces. We can describe the media as homogenous or heterogeneous in each block. We define the velocities at the given rectangle grid points for each block, and the heterogeneous velocities in each block can be calculated by a linear interpolation algorithm. The parameters of the velocity grid positions are independent of the model parameterization, which is advantageous in the joint inversion of the velocities and the node depths of an interface. We implement a segmentally iterative ray tracer to calculate traveltimes in the 3D heterogeneous block models. The damped least squares method is employed in seismic traveltime inversion, which includes the partial derivatives of traveltime with respect to the depths of nodes in the triangulated interfaces and velocities defined in rectangular grids. The numerical tests indicate that the node depths of a triangulated interface and homogeneous velocity distributions can be well inverted in a stratified model.     

Non‐uniqueness with refraction inversion – a syncline model study

Non‐uniqueness occurs with the 1D parametrization of refraction traveltime graphs in the vertical dimension and with the 2D lateral resolution of individual layers in the horizontal dimension. The most common source of non‐uniqueness is the inversion algorithm used to generate the starting model. This study applies 1D, 1.5D and 2D inversion algorithms to traveltime data for a syncline (2D) model, in order to generate starting models for wave path eikonal traveltime tomography. The 1D tau‐p algorithm produced a tomogram with an anticline rather than a syncline and an artefact with a high seismic velocity. The 2D generalized reciprocal method generated tomograms that accurately reproduced the syncline, together with narrow regions at the thalweg with seismic velocities that are less than and greater than the true seismic velocities as well as the true values. It is concluded that 2D inversion algorithms, which explicitly identify forward and reverse traveltime data, are required to generate useful starting models in the near‐surface where irregular refractors are common. The most likely tomogram can be selected as either the simplest model or with a priori information, such as head wave amplitudes. The determination of vertical velocity functions within individual layers is also subject to non‐uniqueness. Depths computed with vertical velocity gradients, which are the default with many tomography programs, are generally 50% greater than those computed with constant velocities for the same traveltime data. The average vertical velocity provides a more accurate measure of depth estimates, where it can be derived. Non‐uniqueness is a fundamental reality with the inversion of all near‐surface seismic refraction data. Unless specific measures are taken to explicitly address non‐uniqueness, then the production of a single refraction tomogram, which fits the traveltime data to sufficient accuracy, does not necessarily demonstrate that the result is either ‘correct’ or the most probable.     

DETERMINATION OF LATERAL INHOMOGENEITIES IN REFLECTION SEISMICS BY INVERSION OF TRAVELTIME RESIDUALS*

Lateral inhomogeneities generate fluctuations in the traveltime of seismic waves. By evaluation of these traveltime fluctuations from different source and receiver positions, lateral inhomogeneities can be located using a pseudo inverse matrix method (Aki, Christoffersson and Husebye 1977). The formulation of the problem is possible for transmitted waves as well as for reflected and refracted waves. In reflection seismics this method is of importance, if no reflections from the inhomogeneities themselves, but only reflections from lower boundaries can be observed. The basic assumptions for the mathematical formulation are (1) the average velocities and depths of the reflecting horizons are known already from standard processing methods, and (2) the traveltime residuals are due to lateral velocity changes between different reflectors or between reflectors and the surface. The area of the earth to be considered is divided into layers and the layers into rectangular blocks. The parallel displacement of a ray after passing a disturbed block is neglected, only the traveltime residual is taken into account. In this paper the method and its application to data obtained with two-dimensional models are described.     

基于初至波走时层析成像的Tikhonov正则化与梯度优化算法

初至波走时层析成像是利用地震初至波走时和其传播的射线路径来反演地下介质速度的技术.该问题本质上是一个不适定问题,需要使用正则化方法并辅之以适当的最优化技巧.本文从数值优化的角度介绍了初至波走时层析成像的反演原理,建立了Tikhonov正则化层析成像反演模型并提出求解极小化问题的加权修正步长的梯度下降算法.该方法可以从速度模型的可行域中迭代找到一个最优解.数值试验表明,该方法是可行和有应用前景的.     

三维非均匀地质模型中的逐段迭代射线追踪

地震射线追踪是地震定位、层析成像、偏移等领域的重要正演环节.随着这些领域研究的深入,针对传统的网格结构和层状结构在描述复杂地质模型遇到的很大困难,我们采用大小不等、形状各异的地质块组成的集合体来描述三维复杂地质模型,并用三角形面片来描述地质块之间的物性间断面,理论上可以描述任意复杂的地质模型.为适应任意非均匀速度分布的地质模型,基于费马原理,本文发展了与之相适应的逐段迭代射线追踪方法.该方法属于弯曲法范畴,对路径点采用一阶显式增量修正,相对于传统的迭代法,高效省时.数值试验表明,联合逐段迭代法和伪弯曲法的射线追踪扰动修正方案在三维复杂非均匀块状模型中有适用性和高效性.     

Secondorder interpolation of traveltimes

To carry out a 3D prestack migration of the Kirchhoff type is still a task of enormous computational effort. Its efficiency can be significantly enhanced by employing a fast traveltime interpolation algorithm. High accuracy can be achieved if secondorder spatial derivatives of traveltimes are included in order to account for the curvature of the wavefront. We suggest a hyperbolic traveltime interpolation scheme that permits the determination of the hyperbolic coefficients directly from traveltimes sampled on a coarse grid, thus reducing the requirements in data storage. This approach is closely related to the paraxial ray approximation and corresponds to an extension of the wellknown     method to arbitrary heterogeneous and complex media in 3D. Application to various velocity models, including a 3D version of the Marmousi model, confirms the superiority of our method over the popular trilinear interpolation. This is especially true for regions with strong curvature of the local wavefront. In contrast to trilinear interpolation, our method also provides the possibility of interpolating source positions, and it is 56 times faster than the calculation of traveltime tables using a fast finitedifference eikonal solver.     

Eigenrays in 3D heterogeneous anisotropic media,Part I: Kinematics

We present a new ray bending approach, referred to as the Eigenray method, for solving two‐point boundary‐value kinematic and dynamic ray tracing problems in 3D smooth heterogeneous general anisotropic elastic media. The proposed Eigenray method is aimed to provide reliable stationary ray path solutions and their dynamic characteristics, in cases where conventional initial‐value ray shooting methods, followed by numerical convergence techniques, become challenging. The kinematic ray bending solution corresponds to the vanishing first traveltime variation, leading to a stationary path between two fixed endpoints (Fermat's principle), and is governed by the nonlinear second‐order Euler–Lagrange equation. The solution is based on a finite‐element approach, applying the weak formulation that reduces the Euler–Lagrange second‐order ordinary differential equation to the first‐order weighted‐residual nonlinear algebraic equation set. For the kinematic finite‐element problem, the degrees of freedom are discretized nodal locations and directions along the ray trajectory, where the values between the nodes are accurately and naturally defined with the Hermite polynomial interpolation. The target function to be minimized includes two essential penalty (constraint) terms, related to the distribution of the nodes along the path and to the normalization of the ray direction. We distinguish between two target functions triggered by the two possible types of stationary rays: a minimum traveltime and a saddle‐point solution (due to caustics). The minimization process involves the computation of the global (all‐node) traveltime gradient vector and the traveltime Hessian matrix. The traveltime Hessian is used for the minimization process, analysing the type of the stationary ray, and for computing the geometric spreading of the entire resolved stationary ray path. The latter, however, is not a replacement for the dynamic ray tracing solution, since it does not deliver the geometric spreading for intermediate points along the ray, nor the analysis of caustics. Finally, we demonstrate the efficiency and accuracy of the proposed method along three canonical examples.