Broadband constant-coefficient propagators

The phase error between the real phase shift and the Gazdag background phase shift, due to lateral velocity variations about a reference velocity, can be decomposed into axial and paraxial phase errors. The axial phase error depends only on velocity perturbations and hence can be completely removed by the split‐step Fourier method. The paraxial phase error is a cross function of velocity perturbations and propagation angles. The cross function can be approximated with various differential operators by allowing the coefficients to vary with velocity perturbations and propagation angles. These variable‐coefficient operators require finite‐difference numerical implementation. Broadband constant‐coefficient operators may provide an efficient alternative that approximates the cross function within the split‐step framework and allows implementation using Fourier transforms alone. The resulting migration accuracy depends on the localization of the constant‐coefficient operators. A simple broadband constant‐coefficient operator has been designed and is tested with the SEG/EAEG salt model. Compared with the split‐step Fourier method that applies to either weak‐contrast media or at small propagation angles, this operator improves wavefield extrapolation for large to strong lateral heterogeneities, except within the weak‐contrast region. Incorporating the split‐step Fourier operator into a hybrid implementation can eliminate the poor performance of the broadband constant‐coefficient operator in the weak‐contrast region. This study may indicate a direction of improving the split‐step Fourier method, with little loss of efficiency, while allowing it to remain faster than more precise methods such as the Fourier finite‐difference method.     

Direct interpretation of 2D potential fields for deep structures by means of the quasi-singular points method

The founder of the Russian school of direct interpretation of potential fields (with minimal prior geological‐geophysical information) was V.M. Berezkin, who introduced the operator of total normalized gradient for the 2D interpretation of profile gravity data sets. This operator was successfully applied in searches of hydrocarbon reservoirs. The further development of this approach (the so‐called quasi‐singular points method) has allowed solution also to various structural problems, using mathematical criteria for the transition from extremes of total normalized gradient fields to coordinates of anomalous sources. The main numerical evaluation strategy is based on stabilized downward continuation of field derivatives and specific use of the filtration properties of Fourier series approximation. The characteristic properties of the quasi‐singular points method are: 1) presentation of a more general total normalized gradient function through additional parameters (derivative order m, form of smoothing function Q, number of Fourier coefficients N* with maximal N), optimum values being chosen during a peak‐spectrum analysis of the interpreted function; 2) calculation of the set of total normalized gradient fields for various values of N*/N, representing coordinate systems {x,N*/N} as an ‘axes tree’ of extrema, where each 2D total normalized gradient field is representationally compressed in a 1D line, permitting a) immediate overview of the positions of the axes in all variants of the calculated fields and b) reduction of the retained information, as required in subsequent interpretation; 3) development of two criteria for transition from extrema of total normalized gradient fields to the coordinates of anomaly sources. The quasi‐singular points method is intended for tracing limiting gently‐sloping boundaries, if their micro‐relief features are sources of the interpreted anomaly but sub‐vertical contacts may also be traced. The method has been tested in delineating various geological structures. One of the most challenging, successfully achieved, was tracing of the Moho discontinuity and study of the upper mantle, using only Bouguer anomaly data along interpretation profiles. This is attested in an example of two regional profiles intersecting the European part of Russia. The central part of one of them coincides with the results from a deep seismic profile.     

Application of randomized sampling schemes to curvelet‐based sparsity‐promoting seismic data recovery

Reconstruction of seismic data is routinely used to improve the quality and resolution of seismic data from incomplete acquired seismic recordings. Curvelet‐based Recovery by Sparsity‐promoting Inversion, adapted from the recently‐developed theory of compressive sensing, is one such kind of reconstruction, especially good for recovery of undersampled seismic data. Like traditional Fourier‐based methods, it performs best when used in conjunction with randomized subsampling, which converts aliases from the usual regular periodic subsampling into easy‐to‐eliminate noise. By virtue of its ability to control gap size, along with the random and irregular nature of its sampling pattern, jittered (sub)sampling is one proven method that has been used successfully for the determination of geophone positions along a seismic line. In this paper, we extend jittered sampling to two‐dimensional acquisition design, a more difficult problem, with both underlying Cartesian and hexagonal grids. We also study what we term separable and non‐separable two‐dimensional jittered samplings. We find hexagonal jittered sampling performs better than Cartesian jittered sampling, while fully non‐separable jittered sampling performs better than separable jittered sampling. Two other 2D randomized sampling methods, Poisson Disk sampling and Farthest Point sampling, both known to possess blue‐noise spectra, are also shown to perform well.     

3D magnetotelluric modelling including surface topography

An edge finite‐element method has been applied to compute magnetotelluric (MT) responses to three‐dimensional (3D) earth topography. The finite‐element algorithm uses a single edge shape function at each edge of hexahedral elements, guaranteeing the continuity of the tangential electric field while conserving the continuity of magnetic flux at boundaries. We solve the resulting system of equations using the biconjugate gradient method with a Jacobian preconditioner. The solution gives electric fields parallel to the slope of a surface relief that is often encountered in MT surveys. The algorithm is successfully verified by comparison with other numerical solutions for a 3D‐2 model for comparison of modelling methods for EM induction and a ridge model. We use a 3D trapezoidal‐hill model to investigate 3D topographic effects, which are caused mainly by galvanic effects, not only in the Zxy mode but also in the Zyx mode. If a 3D topography were approximated by a two‐dimensional topography therefore errors occurring in the transverse electric mode would be more serious than those in the transverse magnetic mode.     

求解二阶解耦弹性波方程的低秩分解法和低秩有限差分法

时间域的波场延拓方法在本质上都可以归结为对一个空间-波数域算子的近似.本文基于一阶波数-空间混合域象征,提出一种新的方法求解解耦的二阶位移弹性波方程.该方法采用交错网格,连续使用两次一阶前向和后向拟微分算子,推导得到了解耦的二阶位移弹性波方程的波场延拓算子.由于该混合域象征在伪谱算子的基础上增加了一个依赖于速度模型的补偿项,可以补偿由于采用二阶中心差分计算时间微分项带来的误差,有效地减少模拟结果的数值频散,提高模拟精度.然而,在非均匀介质中,直接计算该二阶的波场延拓算子,每一个时间步上需要做N次快速傅里叶逆变换,其中N是总的网格点数.为了减少计算量,提出了交错网格低秩分解方法;针对常规有限差分数值频散问题,本文将交错网格低秩方法与有限差分法结合,提出了交错网格低秩有限差分法.数值结果表明,交错网格低秩方法和交错网格低秩有限差分法具有较高的精度,对于复杂介质的地震波数值模拟和偏移成像具有重要的价值.     

最优化辛格式广义离散奇异核褶积微分算子地震波场模拟

将波动方程变换至Hamilton体系,构造了一种新的保结构算法,即最优化辛格式广义褶积微分算子(OSGCD). 在时间离散上,首先引入了Lie算子设计二级二阶辛格式,基于最小误差原理得到了优化的辛格式. 在空间离散上,引入广义离散奇异核褶积微分算子计算空间微分,提出了一种有效方法优化GCD并得到了稳定的算子系数. 针对本文发展的新方法,给出了OSGCD稳定性条件. 在数值实验中,将OSGCD与多种方法比较,从精度和计算效率两方面分析了OSGCD的计算优势,计算结果也表明OSGCD长时程以及非均匀介质中地震波模拟亦具有较强能力.     

Tide‐Induced Air Pressure Fluctuations in a Coastal Unsaturated Zone: Effects of Thin Low‐Permeability Pavements

This paper presents an analytical case study to explore one‐dimensional subsurface air pressure variation in a coastal three‐layered unsaturated zone. The upper layer is thin and much less permeable than the middle layer, and water table is located in the very permeable lower layer. An analytical solution was derived to describe the air pressure variation caused by tide‐induced water table fluctuations. We revisited the case study at Hong Kong International Airport conducted by Jiao and Li (2004) who used a two‐dimensional numerical model. The analytical prediction using the parameter values equivalent to the two‐dimensional numerical model agreed very well with the observed air pressure, indicating the validity and applicability of our one‐dimensional model in approximating the actual situation in this coastal zone with adequate accuracy. The analysis revealed that the asphalt pavement played an important role in causing air pressure fluctuations below it. Abnormally high air pressure can be caused beneath the surface pavement when the air permeability decreases due to rainfall infiltration, which may lead to heaving problems during rising tides.     

A nearly analytic exponential time difference method for solving 2D seismic wave equations

In this paper, we propose a nearly analytic exponential time difference (NETD) method for solving the 2D acoustic and elastic wave equations. In this method, we use the nearly analytic discrete operator to approximate the high-order spatial differential operators and transform the seismic wave equations into semi-discrete ordinary differential equations (ODEs). Then, the converted ODE system is solved by the exponential time difference (ETD) method. We investigate the properties of NETD in detail, including the stability condition for 1-D and 2-D cases, the theoretical and relative errors, the numerical dispersion relation for the 2-D acoustic case, and the computational efficiency. In order to further validate the method, we apply it to simulating acoustic/elastic wave propagation in multilayer models which have strong contrasts and complex heterogeneous media, e.g., the SEG model and the Marmousi model. From our theoretical analyses and numerical results, the NETD can suppress numerical dispersion effectively by using the displacement and gradient to approximate the high-order spatial derivatives. In addition, because NETD is based on the structure of the Lie group method which preserves the quantitative properties of differential equations, it can achieve more accurate results than the classical methods.     

Wavefield interpolation in 3D large‐step Fourier wavefield extrapolation

Extrapolating wavefields and imaging at each depth during three‐dimensional recursive wave‐equation migration is a time‐consuming endeavor. For efficiency, most commercial techniques extrapolate wavefields through thick slabs followed by wavefield interpolation within each thick slab. In this article, we develop this strategy by associating more efficient interpolators with a Fourier‐transform‐related wavefield extrapolation method. First, we formulate a three‐dimensional first‐order separation‐of‐variables screen propagator for large‐step wavefield extrapolation, which allows for wide‐angle propagations in highly contrasting media. This propagator significantly improves the performance of the split‐step Fourier method in dealing with significant lateral heterogeneities at the cost of only one more fast Fourier transform in each thick slab. We then extend the two‐dimensional Kirchhoff and Born–Kirchhoff local wavefield interpolators to three‐dimensional cases for each slab. The three‐dimensional Kirchhoff interpolator is based on the traditional Kirchhoff formula and applies to moderate lateral velocity variations, whereas the three‐dimensional Born–Kirchhoff interpolator is derived from the Lippmann–Schwinger integral equation under the Born approximation and is adapted to highly laterally varying media. Numerical examples on the three‐dimensional salt model of the Society of Exploration Geophysicists/European Association of Geoscientists demonstrate that three‐dimensional first‐order separation‐of‐variables screen propagator Born–Kirchhoff depth migration using thick‐slab wavefield extrapolation plus thin‐slab interpolation tolerates a considerable depth‐step size of up to 72 ms, eventually resulting in an efficiency improvement of nearly 80% without obvious loss of imaging accuracy. Although the proposed three‐dimensional interpolators are presented with one‐way Fourier extrapolation methods, they can be extended for applications to general migration methods.     

Out‐of‐plane effects in 2D borehole‐to‐surface resistivity tomography and applications in mineral exploration

In this paper, we discuss the effects of anomalous out‐of‐plane bodies in two‐dimensional (2D) borehole‐to‐surface electrical resistivity tomography with numerical resistivity modelling and synthetic inversion tests. The results of the two groups of synthetic resistivity model tests illustrate that anomalous bodies out of the plane of interest have an effect on two‐dimensional inversion and that the degree of influence of out‐of‐plane body on inverted images varies. The different influences are derived from two cases. One case is different resistivity models with the same electrode array, and the other case is the same resistivity model with different electrode arrays. Qualitative interpretation based on the inversion tests shows that we cannot find a reasonable electrode array to determine the best inverse solution and reveal the subsurface resistivity distribution for all types of geoelectrical models. Because of the three‐dimensional effect arising from neighbouring anomalous bodies, the qualitative interpretation of inverted images from the two‐dimensional inversion of electrical resistivity tomography data without prior information can be misleading. Two‐dimensional inversion with drilling data can decrease the three‐dimensional effect. We employed two‐ and three‐dimensional borehole‐to‐surface electrical resistivity tomography methods with a pole–pole array and a bipole–bipole array for mineral exploration at Abag Banner and Hexigten Banner in Inner Mongolia, China. Different inverse schemes were carried out for different cases. The subsurface resistivity distribution obtained from the two‐dimensional inversion of the field electrical resistivity tomography data with sufficient prior information, such as drilling data and other non‐electrical data, can better describe the actual geological situation. When there is not enough prior information to carry out constrained two‐dimensional inversion, the three‐dimensional electrical resistivity tomography survey is the better choice.     

三维非均匀介质中真振幅地震偏移算子研究

利用三维非均匀介质中的波动方程,进行振幅保真波场偏移算子分解,得到用于真振幅偏移的单程波方程. 经过数学推理,导出裂步Fourier法真振幅偏移和Fourier有限差分法真振幅偏移的算子方程,并给出其具体的实现过程.     

Field experiments constraining the probability distribution of particle travel distances during natural rainstorms on different slope gradients

Rain splash erosion is an important soil transport mechanism on steep hillslopes. The rain splash process is highly stochastic; here we seek to constrain the probability distribution of splash transport distances on natural hillslopes as a function of hillslope gradient and total precipitation depth. Field experiments were conducted under natural precipitation events to observe splash travel on varying slope gradients. The downslope fraction of splash transport on 15°, 25° and 33° gradients were 85%, 96% and 96%, respectively. Maximum splash transport (L_max) was related to the rain splash detachment of soil particles and slope gradient. An empirical relationship of L_max to the precipitation depth and gradient was obtained; it is linearly proportional to hillslope gradient and logarithmically related to precipitation depth. Measured splash distances were calibrated to the fully two‐dimensional (2D) model of splash transport of Furbish et al. (Journal of Geophysical Research 112 : F01001, 2007) that is based on the assumption that radial splash distances are exponentially distributed; calibrated values of mean splash transport distances are an order of magnitude greater than those previously determined in a controlled laboratory setting. We also compared measured data with several one‐dimensional (1D) probability distributions to asses if splash transport distances could be better explained by a heavy‐tailed probability distribution rather than an exponential probability distribution. We find that for hillslopes of 15° and 25°, although a log‐normal probability distribution best describes the data, we find its likelihood is nearly indistinguishable from an exponential distribution based on computing maximum likelihood estimators for all 1D distributions (exponential, log‐normal and Weibull). At 33°, however, we find stronger evidence that measured travel distances are heavy‐tailed. Copyright © 2011 John Wiley & Sons, Ltd.     

一种全局优化的隐式交错网格有限差分算法及其在弹性波数值模拟中的应用

在数值模拟中,隐式有限差分具有较高的精度和稳定性.然而,传统隐式有限差分算法大多由于需要求解大型矩阵方程而存在计算效率偏低的局限性.本文针对一阶速度-应力弹性波方程,构建了一种优化隐式交错网格有限差分格式,然后将改进格式由时间-空间域转换为时间-波数域,利用二范数原理建立目标函数,再利用模拟退火法求取优化系数.通过对均匀模型以及复杂介质模型进行一阶速度-应力弹性波方程数值模拟所得单炮记录、波场快照分析表明:这种优化隐式交错网格差分算法与传统的几种显式和隐式交错网格有限差分算法相比不但降低了计算量,而且能有效的压制网格频散,使弹性波数值模拟的精度得到有效的提高.     

CALCULATION OF THE VERTICAL GRADIENT OF THE GRAVITY FIELD USING THE FOURIER TRANSFORM*

Fourier transform techniques have been used to calculate the theoretical filter (amplitude) response function of Nth order vertical derivative continuation operation. The amplitude response functions of the vertical gradient and its continuation follow from the same. These response functions are subsequently used to calculate the weighting coefficients suitable for two dimensional equispaced data. A shortening operator has been incorporated to limit the extent of the operator. For comparative study, some of the developed coefficient sets and the one presented in this paper are analysed in the frequency domain and their merits and demerits are discussed.     

High‐accuracy practical spline‐based 3D and 2D integral transformations in potential‐field geophysics

Potential, potential field and potential‐field gradient data are supplemental to each other for resolving sources of interest in both exploration and solid Earth studies. We propose flexible high‐accuracy practical techniques to perform 3D and 2D integral transformations from potential field components to potential and from potential‐field gradient components to potential field components in the space domain using cubic B‐splines. The spline techniques are applicable to either uniform or non‐uniform rectangular grids for the 3D case, and applicable to either regular or irregular grids for the 2D case. The spline‐based indefinite integrations can be computed at any point in the computational domain. In our synthetic 3D gravity and magnetic transformation examples, we show that the spline techniques are substantially more accurate than the Fourier transform techniques, and demonstrate that harmonicity is confirmed substantially better for the spline method than the Fourier transform method and that spline‐based integration and differentiation are invertible. The cost of the increase in accuracy is an increase in computing time. Our real data examples of 3D transformations show that the spline‐based results agree substantially better or better with the observed data than do the Fourier‐based results. The spline techniques would therefore be very useful for data quality control through comparisons of the computed and observed components. If certain desired components of the potential field or gradient data are not measured, they can be obtained using the spline‐based transformations as alternatives to the Fourier transform techniques.     

地震信号内插与噪音剔除（一）

信号内插与噪音剔除是相互有着有机联系的两个方面。首先讨论地震信号的内插方法,提出了一种“最佳内插算子”,从而导出了另一种“检噪算子”,它可以有效地把干扰波从记录中识别出来,并将其“剔除”出去。进而讨论空间域数据的内插问题--道内插问题。在一定的条件下,空间域数据内插与时间域内插是完全等效的,因此道内插也可以采用与时间域相同的“最佳内插算子”。地震记录上的干扰波其绝大部分表现为空间域的脉冲式干扰,所以用“检噪算子”在空间域作褶积,便能识别干扰,从而对干扰加以剔除。本文解决了一种迭代的剔除方法,取得了好的效果。     

Operator-oriented CRS interpolation

In common‐reflection‐surface imaging the reflection arrival time field is parameterized by operators that are of higher dimension or order than in conventional methods. Using the common‐reflection‐surface approach locally in the unmigrated prestack data domain opens a potential for trace regularization and interpolation. In most data interpolation methods based on local coherency estimation, a single operator is designed for a target sample and the output amplitude is defined as a weighted average along the operator. This approach may fail in presence of interfering events or strong amplitude and phase variations. In this paper we introduce an alternative scheme in which there is no need for an operator to be defined at the target sample itself. Instead, the amplitude at a target sample is constructed from multiple operators estimated at different positions. In this case one operator may contribute to the construction of several target samples. Vice versa, a target sample might receive contributions from different operators. Operators are determined on a grid which can be sparser than the output grid. This allows to dramatically decrease the computational costs. In addition, the use of multiple operators for a single target sample stabilizes the interpolation results and implicitly allows several contributions in case of interfering events. Due to the considerable computational expense, common‐reflection‐surface interpolation is limited to work in subsets of the prestack data. We present the general workflow of a common‐reflection‐surface‐based regularization/interpolation for 3D data volumes. This workflow has been applied to an OBC common‐receiver volume and binned common‐offset subsets of a 3D marine data set. The impact of a common‐reflection‐surface regularization is demonstrated by means of a subsequent time migration. In comparison to the time migrations of the original and DMO‐interpolated data, the results show particular improvements in view of the continuity of reflections events. This gain is confirmed by an automatic picking of a horizon in the stacked time migrations.     

Full waveform inversion and the truncated Newton method: quantitative imaging of complex subsurface structures

Full waveform inversion is a powerful tool for quantitative seismic imaging from wide‐azimuth seismic data. The method is based on the minimization of the misfit between observed and simulated data. This amounts to the solution of a large‐scale nonlinear minimization problem. The inverse Hessian operator plays a crucial role in this reconstruction process. Accounting accurately for the effect of this operator within the minimization scheme should correct for illumination deficits, restore the amplitude of the subsurface parameters, and help to remove artefacts generated by energetic multiple reflections. Conventional minimization methods (nonlinear conjugate gradient, quasi‐Newton methods) only roughly approximate the effect of this operator. In this study, we are interested in the truncated Newton minimization method. These methods are based on the computation of the model update through a matrix‐free conjugate gradient solution of the Newton linear system. We present a feasible implementation of this method for the full waveform inversion problem, based on a second‐order adjoint state formulation for the computation of Hessian‐vector products. We compare this method with conventional methods within the context of 2D acoustic frequency full waveform inversion for the reconstruction of P‐wave velocity models. Two test cases are investigated. The first is the synthetic BP 2004 model, representative of the Gulf of Mexico geology with high velocity contrasts associated with the presence of salt structures. The second is a 2D real data‐set from the Valhall oil field in North sea. Although, from a computational cost point of view, the truncated Newton method appears to be more expensive than conventional optimization algorithms, the results emphasize its increased robustness. A better reconstruction of the P‐wave velocity model is provided when energetic multiple reflections make it difficult to interpret the seismic data. A better trade‐off between regularization and resolution is obtained when noise contamination of the data requires one to regularize the solution of the inverse problem.     

Optimized staggered-grid finite-difference operators using window functions

The staggered-grid finite-difference (SGFD) method has been widely used in seismic forward modeling. The precision of the forward modeling results directly affects the results of the subsequent seismic inversion and migration. Numerical dispersion is one of the problems in this method. The window function method can reduce dispersion by replacing the finite-difference operators with window operators, obtained by truncating the spatial convolution series of the pseudospectral method. Although the window operators have high precision in the low-wavenumber domain, their precision decreases rapidly in the high-wavenumber domain. We develop a least squares optimization method to enhance the precision of operators obtained by the window function method. We transform the SGFD problem into a least squares problem and find the best solution iteratively. The window operator is chosen as the initial value and the optimized domain is set by the error threshold. The conjugate gradient method is also adopted to increase the stability of the solution. Approximation error analysis and numerical simulation results suggest that the proposed method increases the precision of the window function operators and decreases the numerical dispersion.