快速射线追踪及其在地震层析中的应用

在地震层析研究中,常常涉及地震波在非均匀介质中的传播路径问题,故快速射线追踪是不可缺少的工具。本文讨论了一种以三点法为基础的逐步逼近的射线追踪方法,其优点是计算速度快、精度高,适合地震层析研究的计算要求。     

非线性插值地震射线追踪方法

针对走时线性插值法(LTI)的计算精度问题,利用非线性插值替代传统的线性插值。该方法利用边界上多个节点的走时,使得插值误差减小,从而提高整体边界网格离散点走时精度,同时保留LTI所具有快速、全局性的优点。数值模拟结果表明,非线性插值算法比LTI算法精度高,可以追踪包括直达波、折射波、透射波等多种射线路径,适合速度突变模型。     

水平层状介质中的快速两点间射线追踪方法

提出了一种快速准确的水平层状介质中两点间的射线追踪方法. 该方法对任意复杂的水平层状速度结构、 任意大小的入射角都能稳定、 快速、 准确地收敛. 将本方法分别与伪弯曲法和Kim方法进行了比较. 结果表明： 本方法比伪弯曲法的计算效率和计算精度更高, 但是只适用于一维速度结构; 在入射角不大的情况下, Kim方法与本方法效果相当. 在大入射角的情况下, Kim方法失效, 但本方法能快速稳定地收敛到真解. 本文提出的射线追踪方法为很多实际地震问题提供了强有力的工具, 但适用范围仅限于水平层状介质.     

基于射线追踪的微地震多波场正演模拟

对于微地震正演模拟,本文以射线追踪的原理为基础,对两点间的射线追踪问题进行了研究,应用二分算法、改进二分算法和微变网格算法对水平层状匀速模型、弯曲层状匀速模型和复杂地质模型进行射线追踪,使得计算效率和适用范围都得到了很大的改善.文中对每种算法误差范围和计算效率进行了对比验证,对于不同的地质模型,选用合适的算法才能在计算速度和精度上得到双重保证,最后正演模拟了多波三分量记录.在模型建立上引入了超薄层概念,并在前人模拟的直达波、透射波、反射波基础上拟了折射波,使正演模拟的多波场信息更丰富.文中的应用实例及模型结果表明:与二分法相比,改进二分法能够对弯曲界面进行射线追踪,并能保证结果的精度.弯曲层状模型中,改进二分法与微变网格法相比计算速度有显著提高,能够应用到资料的反演中.     

复杂介质中地震波前及射线追踪综述

本文较为系统综述了国内外在不均匀介质中各种主要和实用的射线追踪方法,例如:基于射线理论的打靶法、弯曲法(伪弯曲法)、高斯射线束算法等;基于网格单元扩展的有限差分解程函方程法(FD)、最短路径算法(SPM);以及结合射线和网格单元扩展的波前构造法等.同时对目前出现的多次波射线追踪技术、以及多值波前追踪技术(如:相空间算法、水平集算法)也进行了分析讨论.同时对基于网格单元扩展算法的优缺点进行了评述,其基本结论是:基于单元模型的SPM要优于FD算法,而基于网格的SPM算法则次之.就传统的射线追踪算法(如:打靶法和弯曲法)而言,其未来的发展方向是实现完全非线性的相应算法,而基于网格单元的算法则主要是扩展功能(如:后续波、多值波前的追踪).射线追踪方法技术未来需要解决的问题主要有:块状模型中多次波的追踪;多值波前及多值射线追踪;走时与振幅的同时追踪计算;以及其它领域新方法的引入.     

射线追踪方法的发展现状

射线追踪方法作为一种快速有效的波场近似计算方法,不仅对于地震波理论研究具有重要意义,而且也直接应用于地震波反演及偏移成像等过程,本文在收集、整理国内外有关研究资料的基础上,介绍了近十年中这一领域的研究现状及最新发展趋势,并针对其中具有代表性的几类方法的基本思路、方法特点及实现步骤等进行评述。     

对旅行时进行抛物型插值的地震射线追踪方法

提出一种对旅行时进行抛物线插值的地震射线追踪方法（简称PTI方法）,它比基于旅行时线性插值方法（简称LTI方法）计算结果更准确.PTI和LTI方法都是基于2D网格单元模型,用于计算地震波的旅行时和射线路径.首先介绍了相关方法的一些基本概念.旅行时和射线路径都是在网格边界上进行计算的,因此,射线路径在同一恒速网格内是直线.其计算过程有两步.第一步,计算旅行时,第二步追踪射线路径.然后给出了LTI算法的基本公式.因为在炮点网格内可能存在折射波,文章也相应导出了其公式.最后详细推导了PTI 方法的公式.通过模型试算对比说明,用PTI方法较LTI算法更精确、更有效,PTI方法是一种很有发展前途的地震射线追踪算法.     

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

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

最短路径射线追踪方法及其改进

综述了用网络最短路径算法求解地震射线追踪问题的原理、方法技术以及存在问题和改进措施。特别介绍了作者在最短路径算法基础上，提出的动态网络最短路径地震射线追踪方法。该方法先采集从炮点到整个模型所有节点上的初至旅行时，其中，在一个单元内，对相邻每对已计算出最小旅行时的节点进行线性插值，并利用Fermat原理计算未知节点的最小旅行时；然后，利用同样的方法，从接收点开始，反向追踪炮点到接收点的射线路径、该方法能适于各种复杂的非均匀介质，极大地提高了射线追踪的精度。     

Energy Group optimization for forward and inverse problems in nuclear engineering: application to downwell-logging problems

Simulating radiation transport of neutral particles (neutrons and γ‐ray photons) within subsurface formations has been an area of research in the nuclear well‐logging community since the 1960s, with many researchers exploiting existing computational tools already available within the nuclear reactor community. Deterministic codes became a popular tool, with the radiation transport equation being solved using a discretization of phase‐space of the problem (energy, angle, space and time). The energy discretization in such codes is based on the multigroup approximation, or equivalently the discrete finite‐difference energy approximation. One of the uncertainties, therefore, of simulating radiation transport problems, has become the multigroup energy structure. The nuclear reactor community has tackled the problem by optimizing existing nuclear cross‐sectional libraries using a variety of group‐collapsing codes, whilst the nuclear well‐logging community has relied, until now, on libraries used in the nuclear reactor community. However, although the utilization of such libraries has been extremely useful in the past, it has also become clear that a larger number of energy groups were available than was necessary for the well‐logging problems. It was obvious, therefore, that a multigroup energy structure specific to the needs of the nuclear well‐logging community needed to be established. This would have the benefit of reducing computational time (the ultimate aim of this work) for both the stochastic and deterministic calculations since computational time increases with the number of energy groups. We, therefore, present in this study two methodologies that enable the optimization of any multigroup neutron–γ energy structure. Although we test our theoretical approaches on nuclear well‐logging synthetic data, the methodologies can be applied to other radiation transport problems that use the multigroup energy approximation. The first approach considers the effect of collapsing the neutron groups by solving the forward transport problem directly using the deterministic code EVENT, and obtaining neutron and γ‐ray fluxes deterministically for the different group‐collapsing options. The best collapsing option is chosen as the one which minimizes the effect on the γ‐ray spectrum. During this methodology, parallel processing is implemented to reduce computational times. The second approach uses the uncollapsed output from neural network simulations in order to estimate the new, collapsed fluxes for the different collapsing cases. Subsequently, an inversion technique is used which calculates the properties of the subsurface, based on the collapsed fluxes. The best collapsing option is chosen as the one that predicts the subsurface properties with a minimal error. The fundamental difference between the two methodologies relates to their effect on the generated γ‐rays. The first methodology takes the generation of γ‐rays fully into account by solving the transport equation directly. The second methodology assumes that the reduction of the neutron groups has no effect on the γ‐ray fluxes. It does, however, utilize an inversion scheme to predict the subsurface properties reliably, and it looks at the effect of collapsing the neutron groups on these predictions. Although the second procedure is favoured because of (a) the speed with which a solution can be obtained and (b) the application of an inversion scheme, its results need to be validated against a physically more stringent methodology. A comparison of the two methodologies is therefore given.     

Conserved quantities and the inverse problem of reflection seismology

We show that the time-dependent wave equation in both one and two spatial dimensions possesses quantities which are globally conserved. We show how these conserved quantities can be used to determine the characteristic impedance, the rock density and the elastic constant of the rock. We also demonstrate that the conserved quantities possess the capability of determining and/or bracketing the unknown component of the direct pressure response, which is required to begin downward continuation algorithms. Further, we demonstrate that the conserved quantities are always available irrespective of the source structure in time. Numerical instability, arising if the filtering due to the source structure is too harsh, can then be controlled to a degree by demanding that the conserved quantities be indeed conserved.     

Effects of bathymetry on tsunami propagation: Application of ray tracing to tsunamis

Ray tracing of seismic surface waves is applied to tsunami propagation to examine bathymetric effect along its propagation path. Computations are made for trans-Pacific tsunamis and for near-field tsunamis in the Japan Sea. For tsunamis across the Pacific Ocean, the comparison to a uniform ocean shows that focusing and defocusing, due to bathymetry, are significant for some combinations of source and receiver. For example, the refraction of rays is predominant at the East Pacific Rise for the tsunami from Chile. The tsunamis in the Japan Sea are strongly affected by the shallow Yamato Rise. The predicted arrival time and amplitude distribution generally agree with the observations from an actual tsunami. Since the computation can be made very quickly, the method is useful for preliminary analysis of tsunami propagation, such as in an operational warning system or in the determination of computational area for finite-difference computation.     

Trefftz method for the solution of three-dimensional forward and inverse problems of geoelectrics

The Trefftz method is applied to the numerical solution of the three-dimensional (3D) forward problem for an electromagnetic field harmonically oscillating with time in a 3D environment, which ensures the possibility to solve a 3D inverse problem. Here, the known idea of simultaneous joint interpretation of the tangential components of electric and magnetic fields measured on the Earth's surface, which was suggested by A.N. Tikhonov for the development of the electromagnetic sounding method, is employed in the context of the numerical-analytical representation of the forward problem solution in accordance with the Trefftz method. Application of this method allows enables one to generalize the basic plane-stratified vertically 1D model of the medium by allowing for lateral variations in electric conductivity.     

A least-squares penalty method algorithm for inverse problems of steady-state aquifer models

Based on the generalized Gauss–Newton method, a new algorithm to minimize the objective function of the penalty method in (Bentley LR. Adv Wat Res 1993;14:137–48) for inverse problems of steady-state aquifer models is proposed. Through detailed analysis of the “built-in” but irregular weighting effects of the coefficient matrix on the residuals on the discrete governing equations, a so-called scaling matrix is introduced to improve the great irregular weighting effects of these residuals adaptively in every Gauss–Newton iteration. Numerical results demonstrate that if the scaling matrix equals the identity matrix (i.e., the irregular weighting effects of the coefficient matrix are not balanced), our algorithm does not perform well, e.g., the computation cost is higher than that of the traditional method, and what is worse is the calculations fail to converge for some initial values of the unknown parameters. This poor situation takes a favourable turn dramatically if the scaling matrix is slightly improved and a simple preconditioning technique is adopted: For naturally chosen simple diagonal forms of the scaling matrix and the preconditioner, the method performs well and gives accurate results with low computational cost just like the traditional methods, and improvements are obtained on: (1) widening the range of the initial values of the unknown parameters within which the minimizing iterations can converge, (2) reducing the computational cost in every Gauss–Newton iteration, (3) improving the irregular weighting effects of the coefficient matrix of the discrete governing equations. Consequently, the example inverse problem in Bentley (loc. cit.) is solved with the same accuracy, less computational effort and without the regularization term containing prior information on the unknown parameters. Moreover, numerical example shows that this method can solve the inverse problem of the quasilinear Boussinesq equation almost as fast as the linear one.In every Gauss–Newton iteration of our algorithm, one needs to solve a linear least-squares system about the corrections of both the parameters and the groundwater heads on all the discrete nodes only once. In comparison, every Gauss–Newton iteration of the traditional method has to solve the discrete governing equations as many times as one plus the number of unknown parameters or head observation wells (Yeh WW-G. Wat Resour Res 1986;22:95–108).All these facts demonstrate the potential of the algorithm to solve inverse problems of more complicated non-linear aquifer models naturally and quickly on the basis of finding suitable forms of the scaling matrix and the preconditioner.     

Nonlinear ray perturbation theory with its applications to ray tracing and inversion in anisotropic media

A comprehensive approach, based on the general nonlinear ray perturbation theory (Druzhinin, 1991), is proposed for both a fast and accurate uniform asymptotic solution of forward and inverse kinematic problems in anisotropic media. It has been developed to modify the standard ray linearization procedures when they become inconsistent, by providing a predictable truncation error of ray perturbation series. The theoretical background consists in a set of recurrent expressions for the perturbations of all orders for calculating approximately the body wave phase and group velocities, polarization, travel times, ray trajectories, paraxial rays and also the slowness vectors or reflected/transmitted waves in terms of elastic tensor perturbations. We assume that any elastic medium can be used as an unperturbed medium. A total 2-D numerical testing of these expressions has been established within the transverse isotropy to verify the accuracy and convergence of perturbation series when the elastic constants are perturbed. Seismological applications to determine crack-induced anisotropy parameters on VSP travel times for the different wave types in homogeneous and horizontally layered, transversally isotropic and orthorhombic structures are also presented. A number of numerical tests shows that this method is in general stable with respect to the choice of the reference model and the errors in the input data. A proof of uniqueness is provided by an interactive analysis of the sensitivity functions, which are also used for choosing optimum source/receiver locations. Finally, software has been developed for a desktop computer and applied to interpreting specific real VSP observations as well as explaining the results of physical modelling for a 3-D crack model with the estimation of crack parameters.     

Modern trends in the solution of forward and inverse 3D electromagnetic induction problems

The main results obtained during the last 5–8 yr in the solution of forward and inverse problems of 3D induction studies are summarized. The up-to-date status of 3D modelling is presented and prospective improvements in the formulation and numerical solution of forward problems are discussed. Approximate techniques and practical aspects of 3D modelling are specially considered. The general scheme of 3D interpretation of electromagnetic geophysical observations is outlined and realistic formalized approaches to solving 3D inverse problems, namely direct inversion and formalized model fitting, are studied.     

一种最短路径射线追踪的快速算法

为提高最短路径射线追踪的精度,需要增加模型的剖分网格和离散节点,并增加子波传播方向,或者采用其他方法改善计算结果,这些处理会带来大量的额外计算.本文的快速算法改进了波前点的管理和子波传播的计算这两项耗时的工作,较大幅度地提高了传统算法的效率.在波前点的管理上,采用按时间步划分区间的方法,实现了波前点的桶排序管理,其效率高于传统方法中常用的堆排序算法. 在子波传播的计算上,利用斯奈尔定律,同时参考来自邻近节点的波的走时,来限定当前子波传播的有效区域,排除大量不需要计算的子波传播方向. 模型实算表明,本文快速算法的计算速度是传统方法的几倍至十多倍.     

Minimax approach to inverse problems of geophysics

A new approach is suggested for solving the inverse problems that arise in the different fields of applied geophysics (gravity, magnetic, and electrical prospecting, geothermy) and require assessing the spatial region occupied by the anomaly-generating masses in the presence of different types of a priori information. The interpretation which provides the maximum guaranteed proximity of the model field sources to the real perturbing object is treated as the best interpretation. In some fields of science (game theory, economics, operations research), the decision-making principle that lies in minimizing the probable losses which cannot be prevented if the situation develops by the worst-case scenario is referred to as minimax. The minimax criterion of choice is interesting as, instead of being confined to the indirect (and sometimes doubtful) signs of the “optimal” solution, it relies on the actual properties of the information in the results of a particular interpretation. In the hierarchy of the approaches to the solution of the inverse problems of geophysics ordered by the volume and quality of the retrieved information about the sources of the field, the minimax approach should take special place.