Neighborhood discontinuities in bivariate interpolation of scattered observations

Because of the need for computational efficiency, bivariate interpolation methods applied to scattered observations often involve two stages. Initially the variable is estimated at regular grid nodes using a running subset of data (usually of fixed number). This, however, will produce discontinuities in the interpolated surface. Thus a second stage, curvilinear interpolation technique, is applied to estimated values to smooth out the effect of discontinuities. Such problems can be overcome efficiently in processing large data sets by interpolating over natural neighbor subsets. Interpolation procedures that generate discontinuities in the interpolated surface are inappropriate for geological applications, where dislocations due to structural complications may be present.     

An approximate nearest neighbors search algorithm for low-dimensional grid locations

We propose a new algorithm for the problem of approximate nearest neighbors (ANN) search in a regularly spaced low-dimensional grid for interpolation applications. It associates every sampled point to its nearest interpolation location, and then expands its influence to neighborhood locations in the grid, until the desired number of sampled points is achieved on every grid location. Our approach makes use of knowledge on the regular grid spacing to avoid measuring the distance between sampled points and grid locations. We compared our approach with four different state-of-the-art ANN algorithms in a large set of computational experiments. In general, our approach requires low computational effort, especially for cases with high density of sampled points, while the observed error is not significantly different. At the end, a case study is shown, where the ionosphere dynamics is predicted daily using samples from a mathematical model, which runs in parallel at 56 different longitude coordinates, providing sampled points not well distributed that follow Earth’s magnetic field-lines. Our approach overcomes the comparative algorithms when the ratio between the number of sampled points and grid locations is over 2849:1.     

二维输运方程高精度数值模拟

采用剖开算子法,把二维输运问题剖分为两个子初值问题(对流分步、扩散分步)。在任意三角形网格中,分别对不同性质的算子采用各自适合的算法,即采用特征线法求解对流分步,采用半隐式有限元法求解扩散分步。重点探讨了对流插值问题,给出了一种完全对称三次插值模式,有效地减少了数值阻尼。为了克服高阶插值数值震荡问题,计算中保证了函数及其一阶偏导数连续。算例表明,数值方法模拟结果与精确解吻合较好。该算法在求解输运方程(包括纯对流输运方程)时,既能有效减少数值阻尼,也能保证计算中不出现数值震荡。     

Moving Surface Spline Interpolation Based on Green’s Function

Some commonly used interpolation algorithms are analyzed briefly in this paper. Among all of the methods, biharmonic spline interpolation, which is based on Green’s function and proposed by Sandwell, has become the mainstream method for its high precision, simplicity and flexibility. However, the minimum curvature method has two flaws. First, it suffers from undesirable oscillations between data points, which is solved by interpolation with splines in tension. Second, the computation time is approximately proportional to the cube of the number of data constraints, making the method slow for situations with dense data coverage. Focusing on the second problem, this paper introduces the moving surface spline interpolation method based on Green’s function, and the interpolation error equations are deduced. Because the proposed method only chooses the nearest data points by using the merge sort algorithm for interpolating, the computation time is greatly decreased. The optimal number of the nearest points can be determined by using the interpolation error estimation equation. No matter how many data points there are, this method can be implemented without difficulty. Examples show that the proposed method can obtain high interpolation precision and high computation speed at the same time.     

地球物理勘探中几种二维插值方法的误差分析

为了研究采样和网格化方法对地球物理数据成图精度的影响,为野外数据采集布设提供一定的依据,采用数值模拟确定重力异常场场值,通过不同采样间距和不同插值方法计算重力异常绝对误差均方根值和节点处的绝对误差值,对比不同插值方法的误差,得到了如下认识：1）对于同一插值方法而言,存在小间距绝对误差均方根值小于大间距绝对误差均方根值的关系。2）对不同的插值方法而言:当采样间距小于最小异常地质体尺度时,绝对误差均方根值由小到大的顺序是径向基函数法、改进的谢别德法、克里金插值法、自然邻点法、反距离加权插值法、最近邻点法、最小曲率法,并且线性插值三角网法与自然邻点法具有几乎相同的数值;当采样间距大于最小异常地质体尺度时,绝对误差均方根值由小到大的顺序是径向基函数法、改进的谢别德法、克里金插值法、自然邻点法、最小曲率法、最近邻点法、反距离加权插值法,并且线性插值三角网法和自然邻点法具有几乎相同的数值。3）从绝对误差均方值看,径向基函数方法、改进的谢别德方法和克里金方法数值较小,其中径向基函数值绝对误差均方根值最小。4）从节点处绝对误差值来看,径向基函数方法、克里金方法、改进的谢别德方法相对其他插值方法具有更小的误差,不存在局部误差较小或较大的情况,是相对较好的插值方法,并且径向基函数方法是最好的。     

Quasi-Monte Carlo integration on the grid for sensitivity studies

In this paper we present error and performance analysis of quasi-Monte Carlo algorithms for solving multidimensional integrals (up to 100 dimensions) on the grid using MPI. We take into account the fact that the Grid is a potentially heterogeneous computing environment, where the user does not know the specifics of the target architecture. Therefore parallel algorithms should be able to adapt to this heterogeneity, providing automated load-balancing. Monte Carlo algorithms can be tailored to such environments, provided parallel pseudorandom number generators are available. The use of quasi-Monte Carlo algorithms poses more difficulties. In both cases the efficient implementation of the algorithms depends on the functionality of the corresponding packages for generating pseudorandom or quasirandom numbers. We propose efficient parallel implementation of the Sobol sequence for a grid environment and we demonstrate numerical experiments on a heterogeneous grid. To achieve high parallel efficiency we use a newly developed special grid service called Job Track Service which provides efficient management of available computing resources through reservations.     

Enhancement of backtracking search algorithm for identifying soil parameters

In this paper, an enhanced backtracking search algorithm (so-called MBSA-LS) for parameter identification is proposed with two modifications: (a) modifying the mutation of original backtracking search algorithm (BSA) considering the contribution of current best individual for accelerating convergence speed and (b) novelly incorporating an efficient differential evolution (DE) as local search for improving the quality of population. The proposed MBSA-LS is first validated with better performance than the original BSA and some other typical state-of-the-art optimization algorithms on a benchmark of soil parameter identification in terms of effectiveness, efficiency, and robustness. Then, the efficiency of the MBSA-LS is further illustrated by two representative cases: identifying soil parameters from both laboratory tests and field measurements. All comparisons demonstrate that the proposed MBSA-LS algorithm can give accurate results in a short time. Finally, to conveniently solve the problems of parameter identification, a practical tool ErosOpt for parameter identification is developed by integrating the proposed MBSA-LS and some other efficient algorithms for readers to conduct the parameter identification using optimisation algorithms.     

基于贴体全交错网格的起伏地表正演模拟影响因素

贴体网格有限差分正演模拟算法不仅能够精确模拟任意起伏地形下的波场特征,且计算效率较高,是一种很有应用前景的处理西部复杂地表问题的方法;然而,目前求解波动方程时常用的同位网格和标准交错网格,在处理贴体网格起伏地表正演模拟时存在诸多问题。为此,将全交错网格引入到曲线坐标系下,避免了标准交错网格的插值误差和同位网格中奇偶失联引起的高频振荡现象,提高了模拟精度,减小了算法实现的复杂度。在自由边界条件实施时,采用牵引力镜像法计算速度分量,速度自由边界条件配合紧致交错差分格式更新应力分量,得到了较好的效果。随后,重点研究了贴体全交错网格正演模拟算法的影响因素,考虑了网格正交性、网格间距和网格拼接等的影响,并取得了如下认识：算法对网格的正交性没有过分要求;网格间距的突变会引起虚假反射的产生;不同类型的网格拼接对模拟结果不会造成明显的影响。     

基于统计模型的气象数据无损压缩新方法

对当前广泛使用的气象格点数据结构进行了统计分析，通过分析常用气象要素格点资料相邻格点之间的相关性及计算要素场的符号熵和信息冗余度，认为气象格点数据中存在着明显的信息冗余，具有很高的可压缩性，且相关性越好，可压缩性越高。在此分析基础上建立了气象格点数据的二维线性预测统计模型，剔除冗余信息，并结合Huffman编码，提出了一种气象格点数据无损压缩新方法。该方法可极大提高气象格点数据的压缩率，且能保证在有效精度内数据完全无损。最后对常用气象格点资料进行了压缩对比试验，结果表明，该方案压缩效果明显优于当前国际通用的气象数据压缩编码格式(如GRIB和netCDF码)，从而能够大大提高气象以及地球科学中海量数据采集、存储和传输交换的业务应用效率。     

基于Chebyshev走时逼近的三维多次反射射线计算

在利用有限差分等基于网格的数值分析方法解地震波走时所满足的程函方程时,由于速度模型的网格化离散等原因,会使走时在各网格节点之间不具有计算射线路径所要求的光滑性,即走时在邻近网格节点之间不具有连续的一阶导数。因此,直接利用网格节点走时计算射线路径会使最终的射线路径不光滑。为解决这个问题,已有研究者提出了基于B样条插值的逆向梯度方案（法）。然而,在速度发生突变时,B样条逆向梯度法所计算出的射线路径会具有较大的误差。针对这个问题,首先采用适合于解最小零偏差逼近及最佳平方逼近问题的Chebyshev多项式取代B样条对来自于分区多级计算方案的网格节点走时进行最佳逼近,得到在最小平方意义下的最优走时公式;然后采用与B样条逆向梯度法类似的计算过程得到光滑的射线路径。数值实验表明,利用Chebyshev多项式逼近走时可以得到具有很高精度的多次反射射线路径,在多次波偏移成像研究中具有潜在的价值。     

Parallel viewshed analysis on a PC cluster system using triple-based irregular partition scheme

Using digital elevation models (DEMs), viewshed analysis algorithms determine the visibility of each point on the terrain at a given location in space. As a data-parallel algorithm, real-time viewshed analysis from grid DEM poses a practical challenge to personal computer (PC) users, particularly when dealing with higher resolution and accuracy of large terrain data. Therefore, this paper presents a universal domain decomposition algorithm based on an equal-area strategy for the parallel viewshed analysis on a PC cluster system. The approach uses a scan-line filling method for data partitioning of the irregular bounding polygon of the terrain. The terrain data are divided into sectors of the same area that are connected by the viewpoint and the region vertices, ignoring the null value (or NODATA) points. Furthermore, each sector is assigned to one processor and is organized in the form of triples composed of location and elevation at one point. An index of triples is built to store all of the locations of terminal vertices row-by-row and thus the random access of any point is achieved by using the offsets in each row. Two commonly applied viewshed algorithms, namely, “reference plane” and “Xdraw” algorithms are employed to verify the performance. In addition, two experiments focus on evaluating the efficiency performance and comparing traditional implementation, respectively. Experimental results demonstrate a significant performance improvement compared with the sequential computing method. The memory usage gradually decreases as the number of processors increases. Based on the equal-area decomposition, partitions in terms of sectors can guarantee a suitable load balance. Additional benefits of the proposed solution also include high storage efficiency and program portability.     

基于三次插值的大地电磁自适应有限元二维正演模拟

大地电磁(MT)数值模拟中通常使用有限单元法,通过伽辽金(Galerkin)法将微分方程转化为与其等价的泛函形式,对泛函求取极值并在单元上定义插值基函数,得到节点上电磁场值的线性方程组,最终形成大型复对称稀疏矩阵。要达到较高的有限元计算精度,一般采用密集的网格或高次插值的方法,这样做大大的减慢了正演的速度。结合两者的优点利用三次插值和h-型自适应相结合的有限元法来实现MT的正演算法。首先从一个粗网格出发并利用三次插值,通过后验误差估计方法局部加密网格,在计算量较小的情况获得较高的计算精度。这种方法可以针对目标区域和介质分界面发生突变处进行网格加密,不需要全局加密网格。最后通过对国际标准模型COMMEMI-2D1的模拟,分别比较二次插值与三次插值的自适应网格数量和数值模拟结果,证明了三次插值自适应有限元算法的可行性。     

旅行时线性插值射线追踪提高计算精度和效率的改进方法

旅行时线性插值（LTI）射线追踪算法是基于线性假设的,在向前处理过程中仅用按行或按列扫描的方法来计算节点旅行时没有考虑逆向传播射线,导致其计算精度与网格剖分大小有关,在处理复杂介质时会使得追踪出来的射线路径不一定满足最短旅行时。因此,笔者提出了两点改进措施：在向前处理时需采用全方位循环的方法来计算节点最小旅行时;在网格边界加入次生节点。模型试算结果表明：采用全方位循环的LTI法考虑了来自各个方向的射线,可提高其对复杂模型的适应能力;在节点间距相同的情况下,网格边界插入次生节点的LTI法较传统的LTI法计算精度至少可以提高一个数量级,同时,计算速度也更快;随着节点间距剖分的越精细,计算耗时下降也越明显,计算速度较传统的方法可提高n~10n倍。     

Modeling blocky nonlinear mechanics behavior on plasticity and progressive fractures with the enriched mesh free DDA method

This study explores a new form of Discontinuous Deformation Analysis (DDA) method, which uses mesh free displacement functions instead of linear polynomial ones adopted in the original DDA, hence it can effectively describe complex block displacements and deformations. Moreover, the capability of describing a block’s nonlinear mechanical behaviors, i.e., plasticity and fracture, are developed under classical fracture and increment plasticity mechanics theories. With consideration of computation efficiency and convenience, the Sibson natural neighbor interpolation technique for block plasticity analysis and the enriched Moving Least-Squares (MLS) approximation for block fracture analysis are employed, respectively. Numerical results show the applicability of the proposed mesh free DDA method.     

Application of FFT-based Algorithms for Large-Scale Universal Kriging Problems

Looking at kriging problems with huge numbers of estimation points and measurements, computational power and storage capacities often pose heavy limitations to the maximum manageable problem size. In the past, a list of FFT-based algorithms for matrix operations have been developed. They allow extremely fast convolution, superposition and inversion of covariance matrices under certain conditions. If adequately used in kriging problems, these algorithms lead to drastic speedup and reductions in storage requirements without changing the kriging estimator. However, they require second-order stationary covariance functions, estimation on regular grids, and the measurements must also form a regular grid. In this study, we show how to alleviate these rather heavy and many times unrealistic restrictions. Stationarity can be generalized to intrinsicity and beyond, if decomposing kriging problems into the sum of a stationary problem and a formally decoupled regression task. We use universal kriging, because it covers arbitrary forms of unknown drift and all cases of generalized covariance functions. Even more general, we use an extension to uncertain rather than unknown drift coefficients. The sampling locations may now be irregular, but must form a subset of the estimation grid. Finally, we present asymptotically exact but fast approximations to the estimation variance and point out application to conditional simulation, cokriging and sequential kriging. The drastic gain in computational and storage efficiency is demonstrated in test cases. Especially high-resolution and data-rich fields such as rainfall interpolation from radar measurements or seismic or other geophysical inversion can benefit from these improvements.     

A two-grid search scheme for large-scale 3-D finite element analyses of slope stability

It is well known that the trial process for seeking the safety factor in the shear strength reduction finite element method (SSRFEM) is quite expensive, particularly for large 3-D slope stability analyses. The search algorithm for the safety factor is crucial to the entire solution process for the shear strength reduction finite element method, but few studies have attempted to exploit it. Among search algorithms, the commonly used bracketing and bisection search has not been fully optimised. Consequently, to improve the search scheme for the safety factor associated with the shear strength reduction finite element method, two strategies are suggested. First, a generalised bisection search algorithm is proposed to reduce the possibility of encountering non-convergence from a statistical point of view. To further improve the efficiency, a new two-grid scheme, characterised by a coarse mesh search and followed by a fine mesh search, is developed. Based on the drained or undrained analyses of the 3-D slope examples, the new search algorithm can markedly outperform the commonly used bisection search algorithms based on a single finite element mesh.     

应用自然邻接点插值法的块体非连续变形分析

传统的非连续变形分析法(DDA)采用线性位移模式存在诸多缺陷。为准确计算块体应力场,传统上一般直接增加位移函数的多项式阶次,或进行子块体划分或耦合有限元等改进措施,但应用上仍不够方便有效。建议引进无网格节点位移插值模式,采用自然单元法中的自然邻接点插值(NNI)法,具有插值特性,易于准确实施边界条件或材料连续性条件,且具有无网格特征和良好的计算精度,计算更快效。可在此基础上进一步分析大块体弯曲、裂纹扩展破坏形式等,以解决线性位移模式等的不足。     

Mapping Curvilinear Structures with Local Anisotropy Kriging

Interpolating geo-data with curvilinear structures using geostatistics is often disappointing. Channels, for example, become disconnected sets of lakes when interpolated from point data. In order to improve the interpolation of geological structures (e.g., curvilinear structures), we present a new form of kriging, local anisotropy kriging (LAK). Local anisotropy kriging combines a gradient algorithm from image analysis with kriging in an iterative way. After an initial standard kriging interpolation, the gradient algorithm determines the local anisotropy for each cell in the grid using a search area around the cell. Subsequently, kriging is carried out with the spatially varying anisotropy. The anisotropy calculation and subsequent kriging steps will then succeed until the result is satisfactory in the way of reproducing the curvilinear structures. Depending on the size of the search area more or less detail in the geological structures can be reproduced with LAK. Using test examples we show that LAK interpolates data with curvilinear structures more realistically than standard kriging. In a real world case, using bathymetric data of the Oosterschelde estuary, LAK also proves to be quantitatively superior to standard kriging. Absolute interpolation errors are decreased by 23%. Local anisotropy kriging only uses information from point data, which makes the method very objective, it only presents “what the data can tell.”     

Addressing Conditioning Data in Multiple-Point Statistics Simulation Algorithms Based on a Multiple Grid Approach

Multiple-point statistics (MPS) allows simulations reproducing structures of a conceptual model given by a training image (TI) to be generated within a stochastic framework. In classical implementations, fixed search templates are used to retrieve the patterns from the TI. A multiple grid approach allows the large-scale structures present in the TI to be captured, while keeping the search template small. The technique consists in decomposing the simulation grid into several grid levels: One grid level is composed of each second node of the grid level one rank finer. Then each grid level is successively simulated by using the corresponding rescaled search template from the coarse level to the fine level (the simulation grid itself). For a conditional simulation, a basic method (as in snesim) to honor the hard data consists in assigning the data to the closest nodes of the current grid level before simulating it. In this paper, another method (implemented in impala) that consists in assigning the hard data to the closest nodes of the simulation grid (fine level), and then in spreading them up to the coarse grid by using simulations based on the MPS inferred from the TI is presented in detail. We study the effect of conditioning and show that the first method leads to systematic biases depending on the location of the conditioning data relative to the grid levels, whereas the second method allows for properly dealing with conditional simulations and a multiple grid approach.     

Interpolation by regularized spline with tension: II. Application to terrain modeling and surface geometry analysis

A general approach to the computation of basic topographic parameters independent of the spatial distribution of given elevation data is developed. The approach is based on an interpolation function with regular first and second order derivatives and on application of basic principles of differential geometry. General equations for computation of profile, plan, and tangential curvatures are derived. A new algorithm for construction of slope curves is developed using a combined grid and vector approach. Resulting slope curves better fulfill the condition of orthogonality to contours than standard grid algorithms. Presented methods are applied to topographic analysis of a watershed in central Illinois.