2.5D modelling, inversion and angle migration in anisotropic elastic media

2.5D modelling approximates 3D wave propagation in the dip‐direction of a 2D geological model. Attention is restricted to raypaths for waves propagating in a plane. In this way, fast inversion or migration can be performed. For velocity analysis, this reduction of the problem is particularly useful. We review 2.5D modelling for Born volume scattering and Born–Helmholtz surface scattering. The amplitudes are corrected for 3D wave propagation, taking into account both in‐plane and out‐of‐plane geometrical spreading. We also derive some new inversion/migration results. An AVA‐compensated migration routine is presented that is simplified compared with earlier results. This formula can be used to create common‐image gathers for use in velocity analysis by studying the residual moveout. We also give a migration formula for the energy‐flux‐normalized plane‐wave reflection coefficient that models large contrast in the medium parameters not treated by the Born and the Born–Helmholtz equation results. All results are derived using the generalized Radon transform (GRT) directly in the natural coordinate system characterized by scattering angle and migration dip. Consequently, no Jacobians are needed in their calculation. Inversion and migration in an orthorhombic medium or a transversely isotropic (TI) medium with tilted symmetry axis are the lowest symmetries for practical purposes (symmetry axis is in the plane). We give an analysis, using derived methods, of the parameters for these two types of media used in velocity analysis, inversion and migration. The kinematics of the two media involve the same parameters, hence there is no distinction when carrying out velocity analysis. The in‐plane scattering coefficient, used in the inversion and migration, also depends on the same parameters for both media. The out‐of‐plane geometrical spreading, necessary for amplitude‐preserving computations, for the TI medium is dependent on the same parameters that govern in‐plane kinematics. For orthorhombic media, information on additional parameters is required that is not needed for in‐plane kinematics and the scattering coefficients. Resolution analysis of the scattering coefficient suggests that direct inversion by GRT yields unreliable parameter estimates. A more practical approach to inversion is amplitude‐preserving migration followed by AVA analysis. SYMBOLS AND NOTATION A list of symbols and notation is given in Appendix D .     

三维波动方程时空域混合网格有限差分数值模拟方法

常规高阶和时空域高阶有限差分方法广泛应用于三维标量波动方程的数值模拟,这两种差分方法仅利用笛卡尔坐标系中的坐标轴网格点构建三维Laplace差分算子,相应的差分离散波动方程本质上仅具有2阶差分精度,模拟精度低.本文将三维笛卡尔坐标系中非坐标轴网格点分为两类:坐标平面内的非坐标轴网格点和坐标平面外的非坐标轴网格点,系统推...     

Pure quasi-P wave equation and numerical solution in 3D TTI media

Based on the pure quasi-P wave equation in transverse isotropic media with a vertical symmetry axis (VTI media), a quasi-P wave equation is obtained in transverse isotropic media with a tilted symmetry axis (TTI media). This is achieved using projection transformation, which rotates the direction vector in the coordinate system of observation toward the direction vector for the coordinate system in which the z-component is parallel to the symmetry axis of the TTI media. The equation has a simple form, is easily calculated, is not influenced by the pseudo-shear wave, and can be calculated reliably when δ is greater than ε. The finite difference method is used to solve the equation. In addition, a perfectly matched layer (PML) absorbing boundary condition is obtained for the equation. Theoretical analysis and numerical simulation results with forward modeling prove that the equation can accurately simulate a quasi-P wave in TTI medium.     

2.5D direct‐current resistivity forward modelling and inversion by finite‐element–infinite‐element coupled method

To reduce the numerical errors arising from the improper enforcement of the artificial boundary conditions on the distant surface that encloses the underground part of the subsurface, we present a finite‐element–infinite‐element coupled method to significantly reduce the computation time and memory cost in the 2.5D direct‐current resistivity inversion. We first present the boundary value problem of the secondary potential. Then, a new type of infinite element is analysed and applied to replace the conventionally used mixed boundary condition on the distant boundary. In the internal domain, a standard finite‐element method is used to derive the final system of linear equations. With a novel shape function for infinite elements at the subsurface boundary, the final system matrix is sparse, symmetric, and independent of source electrodes. Through lower upper decomposition, the multi‐pole potentials can be swiftly obtained by simple back‐substitutions. We embed the newly developed forward solution to the inversion procedure. To compute the sensitivity matrix, we adopt the efficient adjoint equation approach to further reduce the computation cost. Finally, several synthetic examples are tested to show the efficiency of inversion.     

Improved expressions for the ellipsoidal harmonic series representing the Earth gravitational potential and its first and second derivatives on and outside the reference ellipsoid of revolution

The Legendre functions of the second kind, renormalized by Jekeli, are considered in the external space on a set of ellipsoids of revolution which are confocal with respect to the normal ellipsoid. Among these ellipsoids a reference one is chosen which bounds the Earth. New expressions for the first and second order derivatives of the Legendre functions are derived. They depend on two very quickly convergent Gauss hypergeometric series which are obtained by transforming the slowly convergent initial hypergeometric series. The derived expressions are applied for constructing the ellipsoidal harmonic series for the Earth disturbing gravitational potential and its derivatives of the first and second orders. Since outside the chosen reference ellipsoid there are no Earth masses (as compared to the normal ellipsoid) then it is more appropriate for constructing the boundary-value equation and solving it on the basis of surface gravity data reduced to this ellipsoid.     

2D acoustic‐elastic coupled waveform inversion in the Laplace domain

Although waveform inversion has been intensively studied in an effort to properly delineate the Earth's structures since the early 1980s, most of the time‐ and frequency‐domain waveform inversion algorithms still have critical limitations in their applications to field data. This may be attributed to the highly non‐linear objective function and the unreliable low‐frequency components. To overcome the weaknesses of conventional waveform inversion algorithms, the acoustic Laplace‐domain waveform inversion has been proposed. The Laplace‐domain waveform inversion has been known to provide a long‐wavelength velocity model even for field data, which may be because it employs the zero‐frequency component of the damped wavefield and a well‐behaved logarithmic objective function. However, its applications have been confined to 2D acoustic media. We extend the Laplace‐domain waveform inversion algorithm to a 2D acoustic‐elastic coupled medium, which is encountered in marine exploration environments. In 2D acoustic‐elastic coupled media, the Laplace‐domain pressures behave differently from those of 2D acoustic media, although the overall features are similar to each other. The main differences are that the pressure wavefields for acoustic‐elastic coupled media show negative values even for simple geological structures unlike in acoustic media, when the Laplace damping constant is small and the water depth is shallow. The negative values may result from more complicated wave propagation in elastic media and at fluid‐solid interfaces. Our Laplace‐domain waveform inversion algorithm is also based on the finite‐element method and logarithmic wavefields. To compute gradient direction, we apply the back‐propagation technique. Under the assumption that density is fixed, P‐ and S‐wave velocity models are inverted from the pressure data. We applied our inversion algorithm to the SEG/EAGE salt model and the numerical results showed that the Laplace‐domain waveform inversion successfully recovers the long‐wavelength structures of the P‐ and S‐wave velocity models from the noise‐free data. The models inverted by the Laplace‐domain waveform inversion were able to be successfully used as initial models in the subsequent frequency‐domain waveform inversion, which is performed to describe the short‐wavelength structures of the true models.     

Multistep inversion workflow for 3D long‐offset damped elastic waves in the Fourier domain

We present a new workflow for imaging damped three‐dimensional elastic wavefields in the Fourier domain. The workflow employs a multiscale imaging approach, in which offset lengths are laddered, where frequency content and damping of the data are changed cyclically. Thus, the inversion process is launched using short‐offset and low‐frequency data to recover the long spatial wavelength of the image at a shallow depth. Increasing frequency and offset length leads to the recovery of the fine‐scale features of the model at greater depths. For the fixed offset, we employ (in the imaging process) a few discrete frequencies with a set of Laplace damping parameters. The forward problem is solved with a finite‐difference frequency‐domain method based on a massively parallel iterative solver. The inversion code is based upon the solution of a least squares optimisation problem and is solved using a nonlinear gradient method. It is fully parallelised for distributed memory computational platforms. Our full‐waveform inversion workflow is applied to the 3D Marmousi‐2 and SEG/EAGE Salt models with long‐offset data. The maximum inverted frequencies are 6 Hz for the Marmousi model and 2 Hz for the SEG/EAGE Salt model. The detailed structures are imaged successfully up to the depth approximately equal to one‐third of the maximum offset length at a resolution consistent with the inverted frequencies.     

三维热导方程求解的波数域方法（MWN3D）及模型实验研究

给出了三维热导方程波数域的基本解和算法（MWN3D）,并通过数值模型的方式验证、讨论了这一方法的可靠性和适应能力. 模型实验表明这一方法具有计算速度快,处理问题灵活的特点. 模拟实验还显示,只有在横向较为均匀的三维构造条件下才能用一维模型替代三维模型. 事实上,由于在地壳和上地幔中物质的横向变化不可忽略,通过本方法构建一个真三维的热结构模型是必要的.     

有限元算法在声波方程数值模拟中的频散分析

针对有限元算法在地震波数值模拟中的数值频散问题,利用集中质量矩阵双线性插值有限元算法,推导了二维声波方程的频散函数.在此基础上采用定量分析方法,对比分析了网格纵横长度比变化时的入射方向、空间采样间隔、地震波频率以及地层速度对数值频散的影响.数值算例和模型正演结果表明：当采用集中质量矩阵双线性插值有限元算法时,为了有效地压制数值频散,在所使用震源子波的峰值频率对应的波长内,采样点数目应不少于20个;减小网格长度的纵横比可以有效地抑制入射角(波传播方向与z轴的夹角)较小的地震波的数值频散;地震波频率越高,传播速度越慢,频散越严重,尤其是当相速度与其所对应的频率比值小于2倍空间采样间隔时,不仅会出现严重的数值频散,还会出现假频现象.     

Analytical solutions of the advection‐dispersion equation for transient groundwater flow. A numerical validation

An analytical transport‐model was developed to simulate the propagation of a contaminant in one‐ and two‐dimensional transient flow in groundwater. It is proved that the distribution of concentration at a given time and for a given discharge is identical to that obtained for a different discharge if the volumetric flux of water is the same in the two cases. The results of simulations have been compared with results obtained using the MT3DMS numerical model. There is good agreement when the calculated concentrations are flux‐weighted concentrations. On the other hand, there is a notable divergence when the resident mode is considered. Resident mode concentrations express the mass per unit volume whereas flux mode concentrations express the ratio of mass flux to fluid flux. The solutions presented in this paper can thus be a useful alternative to code MT3DMS when the objective is to simulate concentrations in transient flow according to a resident mode. Copyright © 2008 John Wiley & Sons, Ltd.     

Lebedev scheme for the numerical simulation of wave propagation in 3D anisotropic elasticity‡

This paper presents a Lebedev finite difference scheme on staggered grids for the numerical simulation of wave propagation in an arbitrary 3D anisotropic elastic media. The main concept of the scheme is the definition of all the components of each tensor (vector) appearing in the elastic wave equation at the corresponding grid points, i.e., all of the stresses are stored in one set of nodes while all of the velocity components are stored in another. Meanwhile, the derivatives with respect to the spatial directions are approximated to the second order on two‐point stencils. The second‐order scheme is presented for the sake of simplicity and it is easy to expand to a higher order. Another approach, widely‐known as the rotated staggered grid scheme, is based on the same concept; therefore, this paper contains a detailed comparative analysis of the two schemes. It is shown that the dispersion condition of the Lebedev scheme is less restrictive than that of the rotated staggered grid scheme, while the stability criteria lead to approximately equal time stepping for the two approaches. The main advantage of the proposed scheme is its reduced computational memory requirements. Due to a less restrictive dispersion condition and the way the media parameters are stored, the Lebedev scheme requires only one‐third to two‐thirds of the computer memory required by the rotated staggered grid scheme. At the same time, the number of floating point operations performed by the Lebedev scheme is higher than that for the rotated staggered grid scheme.     

3D Modeling and Born approximation inversion for the borehole surface electromagnetic method

We present a 3D approach to numerical modeling of the borehole-surface electromagnetic （BSEM） method. The 3D electromagnetic response created by a vertical line current source in a layered medium is modeled using the 3D integral equation method. The modeling results are consistent with analytical solutions. 3D Born approximation inversion of BSEM data is also conducted for reservoir delineation. The inversion method is verified by a synthetic reservoir model.     

井地电磁法三维正演模拟和Born近似反演

本文提出一种数值模拟井地电磁法的方法。用体积分方程法对层状介质中的垂直长导线源三维电磁响应做了三维模拟。模拟的结果与解析解对比误差很小,说明算法是正确的。开发了井地电磁法Born近似反演程序,理论模型合成的数据反演结果非常好。     

强剩磁强退磁条件下的二维井中磁测反演

强剩磁、强退磁改变了总磁化强度的大小和方向,给磁测资料解释带来困难.为此,本文利用二维井中磁测数据反演磁化强度矢量的二维分布.首先利用井中磁测的磁异常模量反演磁化强度大小的分布.然后,在已知磁化强度大小分布的前提下,拟合磁场分量,反演磁化强度方向的分布.其中,磁化强度大小和方向均用共轭梯度法求解,并通过预优矩阵改善磁化强度大小的反演效果.理论模拟说明,该方法能准确获得磁化强度矢量分布.磁化强度矢量反演结果包括感磁、剩磁及退磁的影响,这为研究强剩磁、高磁化率矿床提供了一种有效方法.     

伪二维弹性波联合反演近地表的速度和衰减

利用弹性波的初至波和面波,应用交叉梯度算子,联合反演了近地表的二维纵横波速度和衰减参数,并提出了采用一维弹性波正演模拟,应用二维Tikhonov正则化,同时反演出二维速度模型和衰减模型的方法。理论模型测试和实际数据应用结果均表明本文算法极大地提高了计算效率,同时能够反演出可靠的速度模型和衰减模型。      

An Adaptive Wavenumber Sampling Strategy for 2.5D Seismic-Wave Modeling in the Frequency Domain

Simulation of seismic waves from a 3D point-source in a 2D medium may be performed in the frequency-wavenumber domain (called 2.5D modelling). It involves computing the Fourier-transformed Green's function for a number of frequency (ω) and strike direction wavenumber (k_y) values and doubly inverse transforming to convert to the traveltime and distance space. Such modeling produces a wavefield with 3D features but the computation becomes pseudo 2D (i.e., in the xz-plane) rather than 3D (in the xyz-frame). The common sampling strategy for the wavenumber is inefficient for 2.5D wave modeling because it employs a large number of wavenumbers (k_y). This leads to a high cost of computer time in the linear-equation-solving processing, which detracts from the advantages of 2.5D modeling. In this paper, we use two analytic frequency-wavenumber-domain solutions for seismic waves in a homogeneous medium and an inhomogeneous media (two semi-infinite media in contact) to investigate the properties of the solutions and an efficient sampling strategy for choosing the wavenumbers. We have carried out analytic and numerical experiments with these solutions, and present adaptive Gauss–Legendre abscissae for the wavenumber sampling in terms of a modeling situation. We show that the effective range and the number of sampling points of the wavenumber define the adaptive sampling strategy, and they can be estimated in terms of the wavelength and the maximum source-receiver offset. We apply this sampling strategy to the finite-element method and demonstrate that the range and number of sampling points may be adapted for obtaining significant computational efficiency and satisfactory accuracy for every frequency component. Such 2.5D wave modeling can be readily applied for frequency-domain full-waveform inversion for seismic surface measurements and crosshole seismic waveform tomography.     

Analysis of the nonlinear storage–discharge relation for hillslopes through 2D numerical modelling

Storage–discharge curves are widely used in several hydrological applications concerning flow and solute transport in small catchments. This article analyzes the relation Q(S) (where Q is the discharge and S is the saturated storage in the hillslope), as a function of some simple structural parameters. The relation Q(S) is evaluated through two‐dimensional numerical simulations and makes use of dimensionless quantities. The method lies in between simple analytical approaches, like those based on the Boussinesq formulation, and more complex distributed models. After the numerical solution of the dimensionless Richards equation, simple analytical relations for Q(S) are determined in dimensionless form, as a function of a few relevant physical parameters. It was found that the storage–discharge curve can be well approximated by a power law function Q/(LK_s) = a(S/(L²(? ? θ_r)))^b, where L is the length of the hillslope, K_s the saturated conductivity, ? ? θ_r the effective porosity, and a, b two coefficients which mainly depend on the slope. The results confirm the validity of the widely used power law assumption for Q(S). Similar relations can be obtained by performing a standard recession curve analysis. Although simplified, the results obtained in the present work may serve as a preliminary tool for assessing the storage–discharge relation in hillslopes. Copyright © 2012 John Wiley & Sons, Ltd.     

A numerical procedure for obtaining the mean solution and variance of the two-dimensional stochastic convection equation

Under certain conditions the concentration of a substance moving in a stochastic flow field is described by the stochastic convection equation. A numerical method yielding the mean solution and variance of the two-dimensional problem is described here. First, the differential operator is replaced by a discrete linear operator based on finite differences. The resulting system of stochastic equations is then replaced by a system of equations whose solution is the mean concentration. The variance of the concentration can then be calculated. In addition, and example is given for which an approximate analytical solution and its variance is known. The numerical method is applied to the example and results compared to the approximate analytical solution and variance.     

3D description and inversion of reflection moveout of PS‐waves in anisotropic media

Common‐midpoint moveout of converted waves is generally asymmetric with respect to zero offset and cannot be described by the traveltime series t²(x²) conventionally used for pure modes. Here, we present concise parametric expressions for both common‐midpoint (CMP) and common‐conversion‐point (CCP) gathers of PS‐waves for arbitrary anisotropic, horizontally layered media above a plane dipping reflector. This analytic representation can be used to model 3D (multi‐azimuth) CMP gathers without time‐consuming two‐point ray tracing and to compute attributes of PS moveout such as the slope of the traveltime surface at zero offset and the coordinates of the moveout minimum. In addition to providing an efficient tool for forward modelling, our formalism helps to carry out joint inversion of P and PS data for transverse isotropy with a vertical symmetry axis (VTI media). If the medium above the reflector is laterally homogeneous, P‐wave reflection moveout cannot constrain the depth scale of the model needed for depth migration. Extending our previous results for a single VTI layer, we show that the interval vertical velocities of the P‐ and S‐waves (V_P0 and V_S0) and the Thomsen parameters ε and δ can be found from surface data alone by combining P‐wave moveout with the traveltimes of the converted PS(PSV)‐wave. If the data are acquired only on the dip line (i.e. in 2D), stable parameter estimation requires including the moveout of P‐ and PS‐waves from both a horizontal and a dipping interface. At the first stage of the velocity‐analysis procedure, we build an initial anisotropic model by applying a layer‐stripping algorithm to CMP moveout of P‐ and PS‐waves. To overcome the distorting influence of conversion‐point dispersal on CMP gathers, the interval VTI parameters are refined by collecting the PS data into CCP gathers and repeating the inversion. For 3D surveys with a sufficiently wide range of source–receiver azimuths, it is possible to estimate all four relevant parameters (V_P0, V_S0, ε and δ) using reflections from a single mildly dipping interface. In this case, the P‐wave NMO ellipse determined by 3D (azimuthal) velocity analysis is combined with azimuthally dependent traveltimes of the PS‐wave. On the whole, the joint inversion of P and PS data yields a VTI model suitable for depth migration of P‐waves, as well as processing (e.g. transformation to zero offset) of converted waves.