Recursive impedance inversion of ground-penetrating radar data in stochastic media

The travel time and amplitude of ground-penetrating radar (GPR) waves are closely related to medium parameters such as water content, porosity, and dielectric permittivity. However, conventional estimation methods, which are mostly based on wave velocity, are not suitable for real complex media because of limited resolution. Impedance inversion uses the reflection coefficient of radar waves to directly calculate GPR impedance and other parameters of subsurface media. We construct a 3D multiscale stochastic medium model and use the mixed Gaussian and exponential autocorrelation function to describe the distribution of parameters in real subsurface media. We introduce an elliptical Gaussian function to describe local random anomalies. The tapering function is also introduced to reduce calculation errors caused by the numerical simulation of discrete grids. We derive the impedance inversion workflow and test the calculation precision in complex media. Finally, we use impedance inversion to process GPR field data in a polluted site in Mongolia. The inversion results were constrained using borehole data and validated by resistivity data.     

基于Pade逼近的Cole-Cole频散介质GPR有限元正演

针对Cole-Cole频散介质中的复介电常数是jω的分数次幂函数,传统的时域有限元法难以离散及计算时间域分数阶导数,本文采用Pade逼近算法将含有时间分数阶导数的Cole-Cole频散介质电磁波方程推导为一组整数阶辅助微分方程,提出了一种适用于Cole-Cole频散介质的GPR有限元正演模拟算法.在复数伸展坐标系下,通过在频率域Cole-Cole频散介质电磁波方程中引入2个中间变量,并将其变换到时间域,从而以变分形式将PML边界条件加载到Cole-Cole频散介质GPR有限元方程组中,并给出了详细的求解公式.在此基础上,编制了基于Pade逼近的Cole-Cole频散介质GPR有限元正演程序,利用该程序对均匀模型进行计算,并与解析解进行对比,验证了本文构建的GPR有限元正演算法的正确性和有效性.设计了一个复杂Cole-Cole频散介质GPR模型,利用本文构建的GPR有限元正演算法进行模拟并与非频散介质模型的模拟结果进行对比,分析了电磁波在Cole-Cole频散介质中传播衰减增强、子波延伸,分辨率降低等传播规律,有助于实测雷达资料更可靠、更准确的解释.模拟结果表明,基于Pade逼近的GPR有限元正演算法可用于复杂Cole-Cole频散介质结构模拟,且具有较高的计算精度.     

Simultaneous elastic parameter inversion in 2-D/3-D TTI medium combined later arrival times

Traditional traveltime inversion for anisotropic medium is, in general, based on a “weak” assumption in the anisotropic property, which simplifies both the forward part (ray tracing is performed once only) and the inversion part (a linear inversion solver is possible). But for some real applications, a general (both “weak” and “strong”) anisotropic medium should be considered. In such cases, one has to develop a ray tracing algorithm to handle with the general (including “strong”) anisotropic medium and also to design a non-linear inversion solver for later tomography. Meanwhile, it is constructive to investigate how much the tomographic resolution can be improved by introducing the later arrivals. For this motivation, we incorporated our newly developed ray tracing algorithm (multistage irregular shortest-path method) for general anisotropic media with a non-linear inversion solver (a damped minimum norm, constrained least squares problem with a conjugate gradient approach) to formulate a non-linear inversion solver for anisotropic medium. This anisotropic traveltime inversion procedure is able to combine the later (reflected) arrival times. Both 2-D/3-D synthetic inversion experiments and comparison tests show that (1) the proposed anisotropic traveltime inversion scheme is able to recover the high contrast anomalies and (2) it is possible to improve the tomographic resolution by introducing the later (reflected) arrivals, but not as expected in the isotropic medium, because the different velocity (qP, qSV and qSH) sensitivities (or derivatives) respective to the different elastic parameters are not the same but are also dependent on the inclination angle.     

偏离大圆路径传播对四川西部面波相速度成像的影响

面波成像是研究地壳上地幔横波速度结构的一种重要方法．通常,面波相速度或群速度成像假设面波沿大圆路径传播．但是,在地下介质速度结构变化较大时,面波会偏离大圆路径传播,从而导致基于大圆路径假设下的面波成像结果存在一定的误差. 我们采用基于射线追踪的面波成像方法,研究了面波的偏离大圆路径传播对四川西部地区面波相速度成像结果的影响．使用快速行进法(fast marching method)进行面波传播路径的射线追踪,采用子空间反演法(subspace inversion)进行迭代反演,对理论模型合成数据和川西台阵的短周期背景噪声相速度频散数据进行成像分析,并与使用大圆路径传播的成像结果进行对比．对理论模型的测试结果表明,当速度结构变化较大时,基于偏离大圆路径传播的面波成像能够更好地恢复模型异常．对川西台阵的真实数据反演结果显示:在短周期为6 s时,基于偏离大圆路径传播的反演方法较基于大圆路径传播的反演方法所获得的相速度异常的幅度更大些,在四川盆地区域两者的差异接近0.2 km/s;在周期为10 s时,两种反演方法的差异显著减小,基本都在0.1 km/s以内．这主要是因为6 s周期的面波相速度对复杂的上地壳浅层结构更为敏感,从而使得面波的偏离大圆路径传播效应对反演结果的影响更为显著．本文结果表明,当某一周期不同路径的面波相速度测量值变化较大,例如相对于平均相速度的异常超过10％时,则需考虑采用基于偏离大圆路径传播的面波成像方法,否则速度异常较大区域的反演结果可能会造成较大的偏差．      

Ray tracing in elliptical anisotropic media using the linear traveltime interpolation (LTI) method applied to traveltime seismic tomography

The linear traveltime interpolation (LTI) method is a suitable ray‐tracing technique for modelling first‐arrival times in isotropic media. LTI is extended to elliptical anisotropic media and applied to a tomographic inversion procedure. A theoretical formulation is first derived and then LTI implementation is discussed in terms of source–receiver arrays and cell size. The method is then combined with the tomographic inversion procedure adopted. The matching of the ray tracing with inversion in elliptical anisotropic media posed a double non‐linear problem. Thus two assumptions were made: the velocity in each cell is uniform and the main directions of anisotropy are known. To take into account the geometrical characteristics of the area under investigation (depth and velocity of the weathering, and thickness of the inner media), cells of varying size were considered. No hypothesis was made on anisotropy weakness.
The algorithm was first tested on synthetic models and then applied to a field survey. On comparing the results of the synthetic models and the field survey with those obtained with a linear raypath approximation, it was found that there were fewer data misfits.     

2-D Frequency-domain Waveform Inversion of Coupled Acoustic-Elastic Media with an Irregular Interface

In order to correctly interpret marine exploration data, which contain many elastic signals such as S waves, surface waves and converted waves, we have developed both a frequency-domain modeling algorithm for acoustic-elastic coupled media with an irregular interface, and the corresponding waveform inversion algorithm. By applying the continuity condition between acoustic (fluid) and elastic (solid) media, wave propagation can be properly simulated throughout the coupled domain. The arbitrary interface is represented by tessellating square and triangular finite elements. Although the resulting complex impedance matrix generated by finite element methods for the acoustic-elastic coupled wave equation is asymmetric, we can exploit the usual back-propagation algorithm used in the frequency domain through modern sparse matrix technology. By running numerical experiments on a synthetic model, we demonstrate that our inversion algorithm can successfully recover P- and S-wave velocity and density models from marine exploration data (pressure data only).     

跨孔雷达全波形反演成像方法的研究

跨孔雷达全波形反演是一种使用全波形信息反演两钻孔之间地下信息的层析成像技术.常规的层析成像反演大部分采用射线追踪方法,其中基于初至时的射线追踪方法可以反演出速度剖面（介电常数剖面）,基于最大振幅的层析成像可以反演出衰减剖面（电导率剖面）.常规射线追踪方法有许多不足,究其原因是该方法仅使用了小部分的信号信息.为了进一步提高成像分辨率,本文全面推导了全波形跨孔雷达层析成像反演方法,该方法利用雷达波全幅度相位信息能够反演出地下高分辨率的介电常数和电导率图像.本文通过基于局域网的分布式并行算法,有效地解决了巨量数据正演计算问题.文中首先建立了基于单轴各向异性介质完全匹配层的时间域有限差分二维正演算法,进而通过应用包括时间维度在内的全波场信息与残场逆向传播的全波场信息乘积来计算梯度方向,通过求取以步长为自变量的目标函数的极值确定步长公式,并提出以第一次介电常数反演作为同步反演的初始模型,能够有效提高收敛速度.本文对多组模型进行成像实验,取得了较好的反演效果.     

Adaptation of prestack migration to multi-offset ground-penetrating radar (GPR) data

In adapting the prestack migration technique used in seismic imaging to the inversion of ground‐penetrating radar (GPR) from time‐ to depth‐sections, we show that the theoretical integral formulation of the inversion can be applied to electromagnetic problems, albeit with three assumptions. The first two assumptions concern the electromagnetic characteristics of the medium, primarily that the medium must be perfectly resistive and non‐dispersive, and the third concerns the antennae radiation pattern, which is taken to be 2D. The application of this adaptation of the inversion method is confirmed by migrating actual GPR measurements acquired on the test site of the Laboratoire Central des Ponts et Chaussées. The results show good agreement with the geometry of the structures in the medium and confirm that the possible departure from the assumption of a purely resistive medium has no visible effect on the information concerning the geometry of scattering and reflecting structures. The field experiments also show that prestack migration processing is sufficiently robust with regard to the assumption of a non‐dispersive medium. The assumption of a 2D antennae radiation pattern, however, produces artefacts that could be significant for laterally heterogeneous media. Nevertheless, where the medium is not highly laterally heterogeneous, the migration gives a clear image of the scattering potential due to the geometry of structural contrasts in the medium; the scatterers are well focused from diffraction hyperbolae and well localized. Spatial geometry has limited dimensional accuracy and positions are located with a maximum error equal to the minimum wavelength of the signal bandpass. Objects smaller than one wavelength can nevertheless be detected and well focused if their dielectric contrasts are sufficiently high, as in the case of iron or water in gneiss gravels. Furthermore, the suitability of multi‐offset protocols to estimate the electromagnetic propagating velocity and to decrease the non‐coherent noise level of measurements is confirmed. Our velocity estimation is based on the semblance calculation of multi‐offset migrated images, and we confirmed the relevance of this quantification method using numerical data. The signal‐to‐noise ratio is improved by summing multi‐offset results after the addition of random noise on measurements. Thus the adaptation of prestack migration to multi‐offset radar measurements significantly improves the resolution of the scattering potential of the medium. Limitations associated with the methods used here suggest that 3D algorithms should be applied to strongly laterally heterogeneous media and further studies concerning the waveform inversion are necessary to obtain information about the electric nature of the medium.     

Interpretation of shallow refraction seismic data by reflection/refraction tomography

A new algorithm for tomographic inversion of traveltimes of reflected and refracted seismic waves is developed. The inversion gives interface configurations and velocity distributions in layers. The important features of the algorithm are: (a) the inclusion of shot time delays in the list of unknown parameters; (b) the regularization is applied in such a way that the most probable model is characterized by the similarity of neighbouring interfaces. As the problem under consideration is non-linear, several iterations are necessary in order to obtain the final model. In the case of a very inexact initial model, a 'layer-by-layer' inversion strategy is recommended as a first inversion step. The inversion program is supplied with a user interface, thanks to which one can: (a) pick interactively and identify seismic traveltimes; (b) build and edit depth/velocity models; and (c) display calculated traveltime curves and compare them with picked traveltimes as well as with the original seismic sections. The efficiency of the inversion software developed is illustrated by a numerical example and a field example in which shallow seismic data are considered. Application to wide-aperture reflection/refraction profiling (WARRP) data is also possible.     

Estimating Crustal Heterogeneity from Double-difference Tomography

Seismic velocity parameters in limited, but heterogeneous volumes can be inferred using a double-difference tomographic algorithm, but to obtain meaningful results accuracy must be maintained at every step of the computation. MONTEILLER et al. (2005) have devised a double-difference tomographic algorithm that takes full advantage of the accuracy of cross-spectral time-delays of large correlated event sets. This algorithm performs an accurate computation of theoretical travel-time delays in heterogeneous media and applies a suitable inversion scheme based on optimization theory. When applied to Kilauea Volcano, in Hawaii, the double-difference tomography approach shows significant and coherent changes to the velocity model in the well-resolved volumes beneath the Kilauea caldera and the upper east rift. In this paper, we first compare the results obtained using MONTEILLER et al.'s algorithm with those obtained using the classic travel-time tomographic approach. Then, we evaluated the effect of using data series of different accuracies, such as handpicked arrival-time differences (``picking differences'), on the results produced by double-difference tomographic algorithms. We show that picking differences have a non-Gaussian probability density function (pdf). Using a hyperbolic secant pdf instead of a Gaussian pdf allows improvement of the double-difference tomographic result when using picking difference data. We completed our study by investigating the use of spatially discontinuous time-delay data.     

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 .     

Simulation and analysis of GPR signal based on stochastic media model with an ellipsoidal autocorrelation function

Target detection using ground penetrating radar (GPR) is based on the contrast between the electrical parameters of the target and the background medium, such as dielectric permittivity, conductivity and permeability. The application mainly concentrates on the detection of the medium interface and the target shape. In any theoretical study, a simulation model is built with a homogeneous medium. However, real detection encounters heterogeneous media which might produce scattering and diffraction at electrical interfaces and distort the radar pulse shape and affect the detection resolution. In this paper, we build multi-scale random media model with an ellipsoidal autocorrelation function and use FDTD method to simulate the GPR signal response. We then estimate and analyze the arrival time, layer thickness, permittivity and the physics relation in different scale random models according to the S transform method and the transmission wave method. The results demonstrate that we can use GPR to obtain geophysical information of multi-scale heterogeneous media, and provide a foundation for real media detection and complex media inversion.     

Feasibility of joint 1D velocity model and event location inversion by the Neighbourhood algorithm

Using a set of synthetic P‐ and S‐wave onsets, computed in a 1D medium model from sources that mimic a distribution of microseismic events induced by hydrofrac treatment to a monitoring geophone array(s), we test the possibility to invert back jointly the model and events location. We use the Neighbourhood algorithm for data inversion to account for non‐linear effects of velocity model and grid search for event location. The velocity model used is composed of homogeneous layers, derived from sonic logging. Results for the case of one and two monitoring wells are compared. These results show that the velocity model can be obtained in the case of two monitoring wells, if they have optimal relative position. The use of one monitoring well fails due to the trade‐off between the velocity model and event locations.     

Ray path of head waves with irregular interfaces

Head waves are usually considered to be the refracted waves propagating along flat interfaces with an underlying higher velocity. However, the path that the rays travel along in media with irregular interfaces is not clear. Here we study the problem by simulation using a new approach of the spectral-element method with some overlapped elements (SEMO) that can accurately evaluate waves traveling along an irregular interface. Consequently, the head waves are separated from interface waves by a time window. Thus, their energy and arrival time changes can be analyzed independently. These analyses demonstrate that, contrary to the case for head waves propagating along a flat interface, there are two mechanisms for head waves traveling along an irregular interface: a refraction mechanism and transmission mechanism. That is, the head waves may be refracted waves propagating along the interface or transmitted waves induced by the waves propagating in the higher-velocity media. Such knowledge will be helpful in constructing a more accurate inversion method, such as head wave travel-time tomography, and in obtaining a more accurate model of subsurface structure which is very important for understanding the formation mechanism of some special areas, such as the Tibetan Plateau.     

三维地震与地面微地震联合校正方法

由于地面微地震监测台站布设在地表,会受到地表起伏、低降速带厚度和速度变化的影响,降低了微地震事件的识别准确度和定位精度,限制了地面微地震监测技术在复杂地表地区的应用.因此,将三维地震勘探技术的思路引入到地面微地震监测中,提出了三维地震与地面微地震联合校正方法,将油气勘探和开发技术更加紧密地结合在一起.根据三维地震数据和低降速带测量数据,通过约束层析反演方法建立精确的近地表速度模型,将地面微地震台站从起伏地表校正到高速层中的平滑基准面上,有效消除复杂近地表的影响.其次,根据射孔数据和声波测井速度信息,通过非线性反演方法建立最优速度模型,由于已经消除复杂近地表的影响,在进行速度模型优化时不需要考虑近地表的影响,因而建立的速度模型更加准确.最后,在精确速度模型的基础上,通过互相关方法求取剩余静校正量,进一步消除了复杂近地表和速度模型近似误差的影响.三维地震与地面微地震联合校正方法采用逐步校正的思路,能够有效消除复杂近地表的影响,提高微地震数据的品质和速度模型的精确度,保证了微地震事件的定位精度,具有良好的应用前景.     

3D waveform inversion of downhole microseismic data for transversely isotropic media

3D anisotropic waveform inversion could provide high-resolution velocity models and improved event locations for microseismic surveys. Here we extend our previously developed 2D inversion methodology for microseismic borehole data to 3D transversely isotropic media with a vertical symmetry axis. This extension allows us to invert multicomponent data recorded in multiple boreholes and properly account for vertical and lateral heterogeneity. Synthetic examples illustrate the performance of the algorithm for layer-cake and ‘hydraulically fractured’ (i.e. containing anomalies that simulate hydraulic fractures) models. In both cases, waveform inversion is able to reconstruct the areas which are sufficiently illuminated for the employed source-receiver geometry. In addition, we evaluate the sensitivity of the algorithm to errors in the source locations and to band-limited noise in the input displacements. We also present initial inversion results for a microseismic data set acquired during hydraulic fracturing in a shale reservoir.     

微测井约束分步层析静校正技术在古峰庄地区的应用

针对地表剧烈起伏,速度纵、横向变化大的复杂区,层析静校正较以往的折射波静校正方法有明显的优势,但是受初至时间拾取精度、炮检距的选择、近地表模型约束等问题的影响,层析反演的精度还不能满足低幅度构造预测的需求,为此提出微测井约束分步层析的静校正方法,即将浅层速度模型与最终近地表模型分步进行层析反演,并应用微测井信息约束浅层速度模型层析反演,有效地提高了近地表速度模型反演的精度,通过在古峰庄地区的应用较好地解决了复杂地表条件下的静校正问题,低幅度构造预测精度得到提高.     

Soil moisture estimation using tomographic ground penetrating radar in a MCMC–Bayesian framework

In this study, we focus on a hydrogeological inverse problem specifically targeting monitoring soil moisture variations using tomographic ground penetrating radar (GPR) travel time data. Technical challenges exist in the inversion of GPR tomographic data for handling non-uniqueness, nonlinearity and high-dimensionality of unknowns. We have developed a new method for estimating soil moisture fields from crosshole GPR data. It uses a pilot-point method to provide a low-dimensional representation of the relative dielectric permittivity field of the soil, which is the primary object of inference: the field can be converted to soil moisture using a petrophysical model. We integrate a multi-chain Markov chain Monte Carlo (MCMC)–Bayesian inversion framework with the pilot point concept, a curved-ray GPR travel time model, and a sequential Gaussian simulation algorithm, for estimating the dielectric permittivity at pilot point locations distributed within the tomogram, as well as the corresponding geostatistical parameters (i.e., spatial correlation range). We infer the dielectric permittivity as a probability density function, thus capturing the uncertainty in the inference. The multi-chain MCMC enables addressing high-dimensional inverse problems as required in the inversion setup. The method is scalable in terms of number of chains and processors, and is useful for computationally demanding Bayesian model calibration in scientific and engineering problems. The proposed inversion approach can successfully approximate the posterior density distributions of the pilot points, and capture the true values. The computational efficiency, accuracy, and convergence behaviors of the inversion approach were also systematically evaluated, by comparing the inversion results obtained with different levels of noises in the observations, increased observational data, as well as increased number of pilot points.     

各向异性TI介质qP反射波走时层析成像

地震走时层析成像是反演地层各向异性参数分布的有效方法,但是关于地震各向异性介质走时层析成像的研究并不多,其技术远远没有达到成熟的阶段.在野外数据采集时,地表反射波观测方式相对井间和垂直地震剖面观测方式的成本更低,利用qP反射波走时反演各向异性参数具有更加广泛的实用价值.本文实现的TI介质地震走时层析成像方法结合了TI介质反射波射线追踪算法、走时扰动方程和非线性共轭梯度算法,它可以对任意强度的TI介质模型进行反演,文中尝试利用qP反射波走时重建TI介质模型的参数图像.利用qP反射波对层状介质模型和块状异常体模型进行走时反演,由于qP波相速度对弹性模量参数和Thomsen参数的偏微分不同,所以可以分别反演弹性模量参数和Thomsen参数.数值模拟结果表明:利用qP反射波可以反演出TI介质模型的弹性模量参数与Thomsen参数,不同模型的走时迭代反演达到了较好的收敛效果,与各向同性介质走时反演结果相比较,各向异性介质走时反演结果具有较好的识别能力.     

Three-dimensional geostatistical inversion of synthetic tomographic pumping and heat-tracer tests in a nested-cell setup

A main purpose of groundwater inverse modeling lies in estimating the hydraulic conductivity field of an aquifer. Traditionally, hydraulic head measurements, possibly obtained in tomographic setups, are used as data. Because the groundwater flow equation is diffusive, many pumping and observation wells would be necessary to obtain a high resolution of hydraulic conductivity, which is typically not possible. We suggest performing heat tracer tests using the same already installed pumping wells and thermometers in observation planes to amend the hydraulic head data set by the arrival times of the heat signals. For each tomographic combinations of wells, we recommend installing an outer pair of pumping wells, generating artificial ambient flow, and an inner well pair in which the tests are performed. We jointly invert heads and thermal arrival times in 3-D by the quasi-linear geostatistical approach using an efficiently parallelized code running on a mid-range cluster. In the present study, we evaluate the value of heat tracer versus head data in a synthetic test case, where the estimated fields can be compared to the synthetic truth. Because the sensitivity patterns of the thermal arrival times differ from those of head measurements, the resolved variance in the estimated field is 6 to 10 times higher in the joint inversion in comparison to inverting head data only. Also, in contrast to head measurements, reversing the flow field and repeating the heat-tracer test improves the estimate in terms of reducing the estimation variance of the estimate. Based on the synthetic test case, we recommend performing the tests in four principal directions, requiring in total eight pumping wells and four intersecting observation planes for heads and temperature in each direction.