Optimization of Cell Parameterizations for Tomographic Inverse Problems

—?We develop algorithms for the construction of irregular cell (block) models for parameterization of tomographic inverse problems. The forward problem is defined on a regular basic grid of non-overlapping cells. The basic cells are used as building blocks for construction of non-overlapping irregular cells. The construction algorithms are not computationally intensive and not particularly complex, and, in general, allow for grid optimization where cell size is determined from scalar functions, e.g., measures of model sampling or a priori estimates of model resolution. The link between a particular cell j in the regular basic grid and its host cell k in the irregular grid is provided by a pointer array which implicitly defines the irregular cell model. The complex geometrical aspects of irregular cell models are not needed in the forward or in the inverse problem. The matrix system of tomographic equations is computed once on the regular basic cell model. After grid construction, the basic matrix equation is mapped using the pointer array on a new matrix equation in which the model vector relates directly to cells in the irregular model. Next, the mapped system can be solved on the irregular grid. This approach avoids forward computation on the complex geometry of irregular grids. Generally, grid optimization can aim at reducing the number of model parameters in volumes poorly sampled by the data while elsewhere retaining the power to resolve the smallest scales warranted by the data. Unnecessary overparameterization of the model space can be avoided and grid construction can aim at improving the conditioning of the inverse problem. We present simple theory and optimization algorithms in the context of seismic tomography and apply the methods to Rayleigh-wave group velocity inversion and global travel-time tomography.     

Ray tracing based on Fermat's principle in irregular grids1

A new method to trace rays in irregular grids based on Fermat's principle of minimum time is introduced. Besides the usual transmitted and reflected waves, refracted, diffracted and converted waves can also be simulated. The proposed algorithm is fast and stable, and respects the reciprocity principle between source and receiver better than procedures adopting the shooting method. It is particularly suited to form part of a traveltime inversion procedure. The use of irregular grids allows adaptation of the earth discretization to the available acquisition geometry and ray distribution, to obtain more stable and reliable tomographic images.     

分层坐标变换法起伏自由地表弹性波叠前逆时偏移

传统有限差分方法在处理起伏地表时存在一些困难,而坐标变换法可将起伏地表映射为水平地表以克服此缺点.但同时,地下构造被变换得更加复杂,导致了波传播和成像的不准确.本文提出了一种分层的坐标变换方法,并应用到了弹性波逆时偏移中,此方法既可以克服起伏地表的影响,又可以不破坏地下构造.波场正向延拓、逆时延拓和分离是在辅助坐标系下完成的,而成像是在笛卡尔坐标系下完成的.通过对简单起伏模型和中原起伏模型的试算证明了本文提出方法的准确性.同时,对两种极端起伏地层高程不准确的情况进行测试可以看出:分层坐标变换逆时偏移方法的成像效果远好于传统坐标变换方法.     

On the generation of coastline-following grids for ocean models—trade-off between orthogonality and alignment to coastlines

Regional ocean models usually utilize orthogonal curvilinear grids that are fit to the coastline of the modeled regions. While the orthogonality of the grid is required from the perspective of the numerical algorithms, the alignment to the irregular coastlines improves the characterization of the land-sea distribution and the ocean simulation. In this article, we carry out fractal analysis of two representative coastal regions and discuss the trade-offs between the orthogonality and coastline alignment during the grid generation of these regions. A new grid generation method based on Schwarz-Christoffel conformal mappings is proposed, with automatic coastal boundary retrieval algorithm that generates resolution dependent boundary for grid generation and alleviates the human efforts involved in traditional methods. We show that for the southeastern Pacific region, the coastline is smooth with low fractal dimension and there exists effective trade-off with a coastline boundary that adjusts to the desired grid resolution. On the contrary, there is no effective trade-off for southeast China seas where the coastline is of higher fractal dimension, and a coarser coastline boundary is recommended for better orthogonality with little loss in coastline alignment. Further numerical study of coastal trapped Kelvin waves for the typical regions demonstrate that the new coastline-fitting grids achieve smaller error in numerical dispersion and higher accuracy. Through analysis, we conclude that for grid generation for regional ocean modeling, modelers should bring into consideration of the multi-scale fractal characteristics of the coastline.     

各向异性介质弹性波传播的三维不规则网格有限差分方法

提出一种新的三维空间不规则网格有限差分方法,模拟具有地形构造的非均匀各向异性介质中弹性波传播过程. 该方法通过具有二阶时间精度和四阶空间精度的不规则交错网格差分算子来近似一阶弹性波动方程,与多重网格不同,无需在精细网格和粗糙网格间进行插值,所有网格点上的计算在同一次空间迭代中完成. 针对具有复杂物性参数和复杂几何特征的地层结构,使用精细不规则网格处理粗糙界面、断层和空间界面等复杂几何构造, 理论分析和数值算例表明,该方法不但节省了大量计算机内存和计算时间,而且具有令人满意的稳定性和精度.     

On the statistical distribution of seismic velocities in Earth's deep mantle

The existence of uncoupled shear (S) and compression (P) wave velocity variations in Earth's mantle is a characteristic that might only be explained by the presence of significant chemical and/or phase heterogeneity, with important implications for the dynamics and evolution of Earth's interior. While making a one-to-one comparison between tomographic models for P and S velocity (V_P and V_S) variations for a particular geographic region is ill-posed, their global statistical distributions reveal several robust characteristics indicative of the nature of uncoupled V_P and V_S in the deep mantle. We find that all of the V_P and V_S model distributions at a given depth are Gaussian-like throughout the lowermost mantle. However, a distinct low velocity feature is present in V_S distributions below ≈ 2200 km depth that is not present or is relatively weak in V_P models. The presence of anomalously low V_S material cannot be explained as an artifact, nor can the absence of a similarly strong feature in P models be ascribed to under-resolution. We propose that this feature can be partly explained by laterally variable occurrences of post-perovskite (pPv) lenses in the D″ layer, however, the persistence of significantly slow V_S regions at heights up to ≈ 700 km or more above the core–mantle boundary is likely to be incompatible with a pPv origin and might only be explained by the presence of a laterally discontinuous layer of chemically distinct material and/or some other kind of phase heterogeneity. There also exist significant discrepancies between tomographic models with respect to the width of the distributions as well as differences between the modeled peak values. We propose a scheme for comparison between different seismic models in which the widths of the dominant features in their statistical distributions is exploited.     

空间自相似性与地球物理和地球化学场的分解方法

场的分解是应用勘查地球物理和地球化学数据所必不可少的数理处理工作，这项工作往往是为了分解不同的场来研究和了解不同的地质以及相应的地质过程，比如划分不同岩石类型以及识别和圈定矿化蚀变睡.闰定地质过程所产生的同类型的地质体并不简单地对应于相同的场值和频率，场的分布特征往往还与地质体的几何形态，。埋藏深度以及与其周围地质体的相对差异有关。然而，由于地质过程往往具有多期次和空间相关性等特点，所产生的地质体或者相应的场往往会呈现空间自相似性或统计自相似性。这种自相似性可以帮助分解不同地质过程所产生的同类地质体，文中所介绍的多重分形（multifractal)滤波方法（S－A）根据场的能谱分布的不同空间自相似性形成不规则甚至分形滤波器，并将场分解为不同同的组份，如局部异常和背景场，应用该方法对多处地球化学和地球物理场的处理结果表明，该方法对区分与矿有关的局部异常与背景场是很有效的，对加拿大NOVA SCOTIA省西南部航空放射性铀钍钾（U-Tn-K)数据处理所分解的异常对反映岩体中与铀－钨－锡（U－W－Sn)矿化有关的围岩饰变具有明显的效果。     

Three-dimensional local grid refinement for block-centered finite-difference groundwater models using iteratively coupled shared nodes: a new method of interpolation and analysis of errors

This paper describes work that extends to three dimensions the two-dimensional local-grid refinement method for block-centered finite-difference groundwater models of Mehl and Hill [Development and evaluation of a local grid refinement method for block-centered finite-difference groundwater models using shared nodes. Adv Water Resour 2002;25(5):497–511]. In this approach, the (parent) finite-difference grid is discretized more finely within a (child) sub-region. The grid refinement method sequentially solves each grid and uses specified flux (parent) and specified head (child) boundary conditions to couple the grids. Iteration achieves convergence between heads and fluxes of both grids. Of most concern is how to interpolate heads onto the boundary of the child grid such that the physics of the parent-grid flow is retained in three dimensions. We develop a new two-step, “cage-shell” interpolation method based on the solution of the flow equation on the boundary of the child between nodes shared with the parent grid. Error analysis using a test case indicates that the shared-node local grid refinement method with cage-shell boundary head interpolation is accurate and robust, and the resulting code is used to investigate three-dimensional local grid refinement of stream-aquifer interactions. Results reveal that (1) the parent and child grids interact to shift the true head and flux solution to a different solution where the heads and fluxes of both grids are in equilibrium, (2) the locally refined model provided a solution for both heads and fluxes in the region of the refinement that was more accurate than a model without refinement only if iterations are performed so that both heads and fluxes are in equilibrium, and (3) the accuracy of the coupling is limited by the parent-grid size—a coarse parent grid limits correct representation of the hydraulics in the feedback from the child grid.     

Calculation of Ground Vertical Displacements,Tilts and Their Distribution Characteristics around an Irregular Load

Based on the load model of a uniform isotropic semi-infinite elastic medium,we deduced a calculation of vertical displacement and tilt and proposed a method of calculation of vertical displacements and tilts caused by irregular load on the ground or underground at a certain point with two-dimensional and three-dimensional shapes. We compared the difference between the simplified model and the irregular model. Finally,the vertical displacements near the irregular load and the distribution of horizontal tilt are presented.The results show that,compared with the point simplified model,the irregular load model has certain advantages for describing the near field. The establishment of a twodimensional irregular load model can help with the calculation of the modal vector superposition after load scattering. The three-dimensional irregular load model can redistribute load through different weights given to the scattered points after the load scattering,and then obtain displacement with the vector calculation method. The results of vector superposition calculation from the scattered irregular load both in two-dimensions and three-dimensions are all convergent obviously as grids become denser,and it is shown that the calculation method is correct and feasible.     

Mesh type tradeoffs in 2D hydrodynamic modeling of flooding with a Godunov-based flow solver

The effect of mesh type on the accuracy and computational demands of a two-dimensional Godunov-type flood inundation model is critically examined. Cartesian grids, constrained and unconstrained triangular grids, constrained quadrilateral grids, and mixed meshes are considered, with and without local time stepping (LTS), to determine the approach that maximizes computational efficiency defined as accuracy relative to computational effort. A mixed-mesh numerical scheme is introduced so all grids are processed by the same solver. Analysis focuses on a wide range of dam-break type test cases, where Godunov-type flood models have proven very successful. Results show that different mesh types excel under different circumstances. Cartesian grids are 2–3 times more efficient with relatively simple terrain features such as rectilinear channels that call for a uniform grid resolution, while unstructured grids are about twice as efficient in complex domains with irregular terrain features that call for localized refinements. The superior efficiency of locally refined, unstructured grids in complex terrain is attributable to LTS; the locally refined unstructured grid becomes less efficient using global time stepping. These results point to mesh-type tradeoffs that should be considered in flood modeling applications. A mixed mesh model formulation with LTS is recommended as a general purpose solver because the mesh type can be adapted to maximize computational efficiency.     

UK 1-D regional velocity models by analysis of variance of P-wave travel times from local earthquakes

A method is presented for deriving 1-D velocity depth models from earthquake bulletin data. The models can be used as initial models for more advanced modelling techniques such as tomographic inversion. The method is useful when there is little or no refraction and long-range reflection survey data. The bulletin travel times are subjected to an analysis of variance, where they are separated into source, distance, and receiving station terms. The distance terms describe the variation of travel time with distance, and the associated trend lines allow 1-D velocity models for the crustal layers to be determined. The velocity models provide an average crustal model for the region derived from local data. This does not include superficial layers which are necessarily poorly determined. Earthquake bulletin P-wave data from propagation paths across three different regions of the UK are employed to illustrate the use of the technique.     

A Stable and Efficient Numerical Algorithm for Unconfined Aquifer Analysis

The nonlinearity of equations governing flow in unconfined aquifers poses challenges for numerical models, particularly in field-scale applications. Existing methods are often unstable, do not converge, or require extremely fine grids and small time steps. Standard modeling procedures such as automated model calibration and Monte Carlo uncertainty analysis typically require thousands of model runs. Stable and efficient model performance is essential to these analyses. We propose a new method that offers improvements in stability and efficiency and is relatively tolerant of coarse grids. It applies a strategy similar to that in the MODFLOW code to the solution of Richard's equation with a grid-dependent pressure/saturation relationship. The method imposes a contrast between horizontal and vertical permeability in gridblocks containing the water table, does not require "dry" cells to convert to inactive cells, and allows recharge to flow through relatively dry cells to the water table. We establish the accuracy of the method by comparison to an analytical solution for radial flow to a well in an unconfined aquifer with delayed yield. Using a suite of test problems, we demonstrate the efficiencies gained in speed and accuracy over two-phase simulations, and improved stability when compared to MODFLOW. The advantages for applications to transient unconfined aquifer analysis are clearly demonstrated by our examples. We also demonstrate applicability to mixed vadose zone/saturated zone applications, including transport, and find that the method shows great promise for these types of problem as well.     

3D multi-fractal analysis of porous media using 3D digital images: considerations for heterogeneity evaluation

Pore structure heterogeneity is a critical parameter controlling mechanical, electrical and flow transport behaviour of rock. Multi-fractal analysis method was used for a heterogeneity comparison of three-dimensional rock samples with different lithology. Six real digital samples, containing three sandstones and three carbonates, were used. Based on the mercury injection capillary pressure test on these samples, we found that the carbonate samples are more heterogeneous than sandstones, but primary results of multi-fractal behaviours for all samples were similar. We show that if multi-fractal is used to evaluate and compare heterogeneity of different samples, one needs to follow some considerations such as (1) all samples must have the same size in pixel, (2) samples volume must be bigger than representative volume element, (3) multi-fractal dimensions should be firstly normalized to a determined porosity value and (4) multi-fractal results should be interpreted based on resolution of the imaging tool (effects of fine scale sub-resolution pores are missed). Results revealed that using normalized fractal dimensions, the real samples were divided to less and high heterogeneous groups. Moreover, the study of scale effect also showed that porous structures of these samples are scale invariant in a wide range of scales (from one to eight times bigger).     

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.     

自适应非结构有限元MT二维起伏地形正反演研究

在山区进行MT勘探时,用规则网格有限元方法模拟起伏地形会受到限制.本文采用非结构三角网格可以有效地模拟任意二维地质结构,如起伏地形、倾斜岩层和多尺度构造等.正演引入自适应有限元方法,其在网格剖分过程中能根据单元误差自动细化网格,保证了正演结果的精度.将自适应有限元与Occam算法结合,且引用并行处理技术提高正反演计算速度.通过对比两个理论模型,讨论了地形对MT正演响应的影响;其次进行了不同地电模型带地形反演展示了本文算法的正确性和适用性;最后将该方法应用于实测MT数据处理,证明了自适应非结构有限元方法是复杂地形下处理MT数据的有力工具.     

起伏地表条件下各向异性地震波最短路径射线追踪

在地震波正反演研究中,考虑起伏地表和地震各向异性具有非常重要的理论意义和实际应用价值.本文在前人研究的基础上,将最短路径追踪算法引入到起伏地表各向异性介质模型的地震波走时计算中.模型剖分时,整体模型划分成正方形单元,起伏边界附近以不规则网格逼近,进而采用非规则节点布置实现非规则网格处的最短路径计算.追踪计算中采用Sena群速度近似公式,得到各向异性地震波的走时,实现了复杂地表情况下各向异性介质模型中地震波的射线追踪.理论模型计算结果显示,本文方法能够可靠地应用于复杂各向异性介质模型,具有较高的计算精度.     

面向目标自适应三维大地电磁正演模拟

本文将面向目标的自适应算法应用于三维大地电磁数值模拟.使用基于非结构网格的矢量有限单元法对起伏地表大地电磁正演模拟问题进行求解.使用利用垂向电流密度在物性界面上的连续性对后验误差进行估算的算法指导网格优化.由于全局自适应算法针对观测点优化网格的能力较差,本文通过求解正演问题的对偶问题计算后验误差的加权系数,并对相关加权系数进行改进,从而实现了面向目标的自适应算法.与传统基于结构化网格的电磁正演算法相比,采用非结构网格能够更好地拟合起伏地表和地下不规则异常体.由于使用了面向目标的自适应算法,本文能够使用更少的网格达到较高的计算精度.通过对比本文模拟结果与半空间响应和全局自适应算法计算结果,并通过对比使用改进前和改进后加权系数得到的网格剖分结果验证了本文算法的有效性.     

EGM2008地球重力模型数据在中国大陆地区的精度分析

本文介绍了5′×5 ′的EGM2008地球重力模型及其在全球的精度评价.按照地形变化规律,将中国大陆大致分为7个区域,在10 km网度上,将EGM2008地球重力模型数据与中国地面实测空间重力网格数据进行了对比.由于数据源的问题,中国大陆的模型数据精度普遍低于北美和欧洲.二种数据在地形平坦的东部地区差别较小,向西随着地形复杂程度的增加,二种数据之间的标准差从小于10 mGal增大到50多mGal.畸变点分析表明精度极低的网格点均分布在地形起伏大的地区.总体而言,5′×5′的EGM2008地球重力模型数据在中国大陆将近80％的面积上的精度可达10 mGal之内,可用于小比例尺重力编图和构造研究.在地形起伏较大的中国西部以青藏高原为例进一步比较了EGM2008重力模型和重力测点数据,结果表明在重力点分布稀疏不均匀的地区,平面网格数据难以准确表达重力场信息.由于缺少地面重力数据控制,EGM2008重力模型数据在中国西部精度较低,但模型数据依然在很大程度上提高了空间重力异常信息的丰富程度.将中国区域重力调查成果数据应用于地球模型的构建是一项有意义的工作.