Use of VFSA for resolution, sensitivity and uncertainty analysis in 1D DC resistivity and IP inversion

We present results from the resolution and sensitivity analysis of 1D DC resistivity and IP sounding data using a non-linear inversion. The inversion scheme uses a theoretically correct Metropolis–Gibbs' sampling technique and an approximate method using numerous models sampled by a global optimization algorithm called very fast simulated annealing (VFSA). VFSA has recently been found to be computationally efficient in several geophysical parameter estimation problems. Unlike conventional simulated annealing (SA), in VFSA the perturbations are generated from the model parameters according to a Cauchy-like distribution whose shape changes with each iteration. This results in an algorithm that converges much faster than a standard SA. In the course of finding the optimal solution, VFSA samples several models from the search space. All these models can be used to obtain estimates of uncertainty in the derived solution. This method makes no assumptions about the shape of an a posteriori probability density function in the model space. Here, we carry out a VFSA-based sensitivity analysis with several synthetic and field sounding data sets for resistivity and IP. The resolution capability of the VFSA algorithm as seen from the sensitivity analysis is satisfactory. The interpretation of VES and IP sounding data by VFSA, incorporating resolution, sensitivity and uncertainty of layer parameters, would generally be more useful than the conventional best-fit techniques.     

Appraisal of equivalence and suppression problems in 1D EM and DC measurements using global optimization and joint inversion*

Global optimization with very fast simulated annealing (VFSA) in association with joint inversion is performed for 1D earth structures. The inherent problems of equivalence and suppression in electromagnetic (EM) and direct current (DC) resistivity methods are studied. Synthetic phase data from multifrequency sounding using a horizontal coplanar coil system and synthetic apparent resistivity data from Schlumberger DC resistivity measurements are inverted individually and jointly over different types of layered earth structures. Noisy data are also inverted. The study reveals that global optimization of individual data sets cannot solve inherent equivalence or suppression problems. Joint inversion of EM and DC measurements can overcome the problem of equivalence very well. However, a suppression problem cannot be solved even after combination of data sets. This study reveals that the K-type earth structure is easiest to resolve while the A-type is the most difficult. We also conclude that the equivalence associated with a thin resistive layer can be resolved better than that for a thin conducting layer.     

电阻率成像反演中的模拟退火方法

电阻率成像反演是一种典型的多参数、非线性反演问题。模拟退火方法将组合优化问题与统计力学中的热平衡问题相类比,开辟了求解组合优化问题的新途径。本文通过数字模拟及对野外资料的处理,证明快速模拟退火方法在二维电阻率成像反演中具有较好的应用前景。     

Background velocity estimation using non-linear optimization for reflection tomography and migration misfit

We show that it is possible to estimate the background velocity for prestack depth migration in 2D laterally varying media using a non-linear optimization technique called very fast simulated annealing (VFSA). We use cubic splines in the velocity model parametrization and make use of either successive pairs of shot gathers or several constant-offset sections as input data for the inversion. A Kirchhoff summation scheme based on first-arrival traveltimes is used to migrate/model the input data during the velocity analysis. We evaluate and compare two different measures of error. The first is defined in the recorded data or (x,t) domain and is based on a reflection-tomography criterion. The second is defined in the migrated data or (x,z) domain and is based on a migration-misfit criterion. Depth relaxation is used to improve the convergence and quality of the velocity analysis while simultaneously reducing the computational cost. Further, we show that by coarse sampling in the offset domain the method is still robust. Our non-linear optimization approach to migration velocity analysis is evaluated for both synthetic and real seismic data. For the velocity-analysis method based on the reflection-tomography criterion, traveltimes do not have to be picked. Similarly, the migration-misfit criterion does not require that depth images be manually compared. Interpreter intervention is required only to restrict the search space used in the velocity-analysis problem. Extension of the proposed schemes to 3D models is straightforward but practical only for the fastest available computers.     

重力与地震资料的模拟退火约束联合反演

联合反演是综合地球物理研究的重要定量解释手段.本文在总结和分析重力与地震资料联合反演的研究现状基础上,利用改进的全局寻优的快速模拟退火算法,实现了重力和地震资料的约束同步联合反演.针对性地设计了密度和速度界面不完全一致的模型,理论模型的试验说明了方法的效果和适用性.结合最近完成的广东徐闻地区实际资料的处理和解释,表明该方法可准确确定复杂构造物性界面的密度和速度结构,在该地区的油气勘探中发挥了作用.在先验信息约束下,该联合反演方法要明显优于单独的重力反演.     

Solutions of the inherent problem of the equivalence in direct current resistivity and electromagnetic methods through global optimization and joint inversion by successive refinement of model space

The problem of equivalence in direct current (DC) resistivity and electromagnetic methods for a thin resistive and conducting layer is well‐known. Attempts have been made in the past to resolve this problem through joint inversion. However, equivalence still remains an unresolved problem. In the present study, an effort is made to reduce non‐uniqueness due to equivalence using global optimization and joint inversion by successive refinement of the model space. A number of solutions derived for DC resistivity data using very fast simulated annealing global inversion that fits the observations equally well, follow the equivalence principle and show a definite trend. For a thin conductive layer, the quotient between resistivity and thickness is constant, while for a resistive one, the product between these magnitudes is constant. Three approaches to obtain very fast simulated annealing solutions are tested. In the first one, layer resistivities and thicknesses are optimized in a linear domain. In the second, layer resistivities are optimized in the logarithmic domain and thicknesses in the linear domain. Lastly, both layer resistivities and thicknesses are optimized in the logarithmic domain. Only model data from the mean models, corresponding to very fast simulated annealing solutions obtained for approach three, always fit the observations. The mean model defined by multiple very fast simulated annealing solutions shows extremely large uncertainty (almost 100%) in the final solution after inversion of individual DC resistivity or electromagnetic (EM) data sets. Uncertainty associated with the intermediate resistive and conducting layers after global optimization and joint inversion is still large. In order to reduce the large uncertainty associated with the intermediate layer, global optimization is performed over several iterations by reducing and redefining the search limits of model parameters according to the uncertainty in the solution. The new minimum and maximum limits are obtained from the uncertainty in the previous iteration. Though the misfit error reduces in the solution after successive refinement of the model space in individual inversion, it is observed that the mean model drifts away from the actual model. However, successive refinement of the model space using global optimization and joint inversion reduces uncertainty to a very low level in 4–5 iterations. This approach works very well in resolving the problem of equivalence for resistive as well as for conducting layers. The efficacy of the approach has been demonstrated using DC resistivity and EM data, however, it can be applied to any geophysical data to solve the inherent ambiguities in the interpretations.     

Global Optimisation of Time Domain Electromagnetic Data Using Very Fast Simulated Annealing

--Efficacy of the global optimisation technique is demonstrated in the inversion of time domain electromagnetic data. Transient EM responses observed using a coincident loops system along several profiles and different time channels over a plate like conducting body are inverted, using very fast simulated annealing (VFSA) as an optimisation tool. Three time channels for each profile are considered in the inversion. Study reveals that only one run of global inversion considering a single multi-channel profile is not enough to resolve all the nine model parameters of a plate-like conducting body. However, the global inversion of a single multi-channel profile with several runs yields a mean model that is quite close to the true model. Considering many profiles and time channels together in the global inversion can yield reliable estimates of all the parameters. Computationally, this is not an efficient procedure. Analysis of the results shows, however, that two distant profiles are enough to yield all the model parameters reliably, even after one run of the global inversion. Both noise-free and noisy synthetic data are used in the inversion. Finally, field data are also inverted to study the performance of the global inversion.     

电阻率和速度随机分布的MT与地震联合反演

在已有研究成果的基础上,为了适应物性参数剧烈变化的复杂模型并满足联合反演的要求,开发了速度和电阻率随机分布共网格单元模型的建模技术.基于这种统一的物性随机分布的网格介质模型,利用有限元方法和改进的射线追踪法分别正演计算大地电磁场和地震走时,结合改进的模拟退火算法,研究实现了电阻率和速度随机分布条件下的大地电磁与地震资料的同步联合反演.对物性界面不完全一致和物性变化剧烈的带地形复杂模型的试验,表明了该方法在精细反演复杂电阻率和速度结构方面的效果,克服了以往研究局限于简单模型的不足.对地震资料品质差的地区开展的实际资料联合反演,表明了方法的适用性,先验信息约束下的联合反演提高了反演精度.     

基于非结构网格的电阻率三维带地形反演

地表起伏地形在野外矿产资源勘察中不可避免,其对直流电阻率法勘探影响巨大.近年来,电阻率三维正演取得诸多进展,特别是应用非结构网格我们能够进行任意复杂地形和几何模型的电阻率三维数值模拟,但面向实际应用的起伏地形下电阻率三维反演依然困难.本文基于非结构化四面体网格,并考虑到应用GPS/GNSS时,区域地球物理调查中可非规则布设测网的实际特点,实现了任意地形(平坦或起伏)条件下、任意布设的偶极-偶极视电阻率数据的不完全Gauss-Newton三维反演.合成数据的反演结果表明了方法的有效性,可应用于复杂野外环境下的三维电法勘探.     

地震与大地电磁测深数据的自适应正则化同步联合反演

在传统的联合反演研究中,地球物理学者往往更多地关注数据拟合,很少涉及正则化理论.本文在电阻率和速度随机分布的大地电磁测深(MT)与地震联合反演研究的基础之上,将正则化思想引入到同步联合反演中,加入先验信息进行模型约束,选取最小模型为稳定泛函,并首次采用自适应正则化算法来确定联合反演的正则化因子.根据以往研究成果,采用非...     

The use of reference models from a priori data to guide 2D inversion of electrical resistivity tomography data

Clay-rich till plains cover much of the UK. Such sites are attractive locations for landfills, since clay aquitards lower the risk of landfill leachate entering groundwater. However, such tills often contain sand and gravel bodies that can act as leachate flow routes. Such bodies may not be detected by conventional site investigation techniques such as drilling boreholes and trial pitting. A method of guided inversion, where a priori data are used to construct structural reference models for use in inverting electrical resistivity tomography data, was proposed as a tool to improve the detection of sand and gravel bodies within clay-rich till sequences.
Following a successful 2D guided inversion synthetic modelling study, a field study was undertaken. Wenner 2D electrical resistivity tomography lines, resistivity cone penetrometry bores and electromagnetic induction ground resistivity data were collected over a site on the East Yorkshire coast, England, where sand and gravel lenses were known to exist from cliff exposures. A number of equally valid geoelectrical models were constructed using the electromagnetic and resistivity cone data. These were used as structural reference models in the inversion of the resistivity tomography data. Blind inversion using an homogenous reference model was also carried out for comparison.
It was shown for the first time that the best solution model produced by 2D inversion of one data set with a range of structural reference models could be determined by using the l ₂ model misfit between the solution models and associated reference models (reference misfit) as a proxy for the l ₂ misfit between the solution models and the synthetic model or 'best-guess' geoelectrical model (true misfit). The 2D methodology developed here is applicable in clay-rich till plains containing sand and gravel bodies throughout the UK.     

Inversion of DC resistivity data using neural networks

The inversion of geoelectrical resistivity data is a difficult task due to its non-linear nature. In this work, the neural network (NN) approach is studied to solve both 1D and 2D resistivity inverse problems. The efficiency of a widespread, supervised training network, the back-propagation technique and its applicability to the resistivity problem, is investigated. Several NN paradigms have been tried on a basis of trial-and-error for two types of data set. In the 1D problem, the batch back-propagation paradigm was efficient while another paradigm, called resilient propagation, was used in the 2D problem. The network was trained with synthetic examples and tested on another set of synthetic data as well as on the field data. The neural network gave a result highly correlated with that of conventional serial algorithms. It proved to be a fast, accurate and objective method for depth and resistivity estimation of both 1D and 2D DC resistivity data. The main advantage of using NN for resistivity inversion is that once the network has been trained it can perform the inversion of any vertical electrical sounding data set very rapidly.     

Two‐dimensional joint inversions of cross‐hole resistivity data and resolution analysis of combined arrays

In this study, a new two‐dimensional inversion algorithm was developed for the inversion of cross‐hole direct current resistivity measurements. In the last decades, various array optimisation methods were suggested for resistivity tomography. However, researchers have still collected data by using classical electrode arrays in most cross‐hole applications. Therefore, we investigated the accuracy of both the individual and the joint inversion of the classical cross‐hole arrays by using both synthetic and field data with the developed algorithm. We showed that the joint inversion of bipole–bipole, pole–bipole, bipole–pole, and pole–tripole electrode arrays gives inverse solutions that are closer to the real model than the individual inversions of the electrode array datasets for the synthetic data inversion. The model resolution matrix of the suggested arrays was used to analyse the inversion results. This model resolution analysis also showed the advantage of the joint inversion of bipole–bipole, pole–bipole, bipole–pole, and pole–tripole arrays. We also used sensitivity sections from each of the arrays and their superpositions to explain why joint inversion gives better resolution than the any individual inversion result.     

1D joint multi‐offset inversion of time‐domain marine controlled source electromagnetic data

The accurate estimation of sub‐seafloor resistivity features from marine controlled source electromagnetic data using inverse modelling is hindered due to the limitations of the inversion routines. The most commonly used one‐dimensional inversion techniques for resolving subsurface resistivity structures are gradient‐based methods, namely Occam and Marquardt. The first approach relies on the smoothness of the model and is recommended when there are no sharp resistivity boundaries. The Marquardt routine is relevant for many electromagnetic applications with sharp resistivity contrasts but subject to the appropriate choice of a starting model. In this paper, we explore the ability of different 1D inversion schemes to derive sub‐seafloor resistivity structures from time domain marine controlled source electromagnetic data measured along an 8‐km‐long profile in the German North Sea. Seismic reflection data reveal a dipping shallow amplitude anomaly that was the target of the controleld source electromagnetic survey. We tested four inversion schemes to find suitable starting models for the final Marquardt inversion. In this respect, as a first scenario, Occam inversion results are considered a starting model for the subsequent Marquardt inversion (Occam–Marquardt). As a second scenario, we employ a global method called Differential Evolution Adaptive Metropolis and sequentially incorporate it with Marquardt inversion. The third approach corresponds to Marquardt inversion introducing lateral constraints. Finally, we include the lateral constraints in Differential Evolution Adaptive Metropolis optimization, and the results are sequentially utilized by Marquardt inversion. Occam–Marquardt may provide accurate estimation of the subsurface features, but it is dependent on the appropriate conversion of different multi‐layered Occam model to an acceptable starting model for Marquardt inversion, which is not straightforward. Employing parameter spaces, the Differential Evolution Adaptive Metropolis approach can be pertinent to determine Marquardt a priori information; nevertheless, the uncertainties in Differential Evolution Adaptive Metropolis optimization will introduce some inaccuracies in Marquardt inversion results. Laterally constrained Marquardt may be promising to resolve sub‐seafloor features, but it is not stable if there are significant lateral changes of the sub‐seafloor structure due to the dependence of the method to the starting model. Including the lateral constraints in Differential Evolution Adaptive Metropolis approach allows for faster convergence of the routine with consistent results, furnishing more accurate estimation of a priori models for the subsequent Marquardt inversion.     

套管井井壁附近地层横波速度径向分布反演

套管外地层受异常地应力、油气开采的影响,在径向上表现出非均质性;采用声波测井可以对该非均质性进行探测,利用偶极子横波测井数据可以对横波速度的径向分布进行反演.本文建立了套管井外地层横波速度径向分层参考模型,采用修正的微扰法计算了该模型的偶极弯曲波频散曲线,建立了横波速度径向分布反演目标函数,采用高斯牛顿法和快速模拟退火法对目标函数进行了求解,得到了套管井外地层横波速度的径向分布.分析了偶极弯曲波频段、套管横波速度对反演结果的影响,对比了高斯牛顿法和快速模拟退火法对反演过程的影响.分析对比结果表明,采用偶极弯曲波激发强度较高的频段与采用全频段的反演结果相近;套管的横波速度准确度越高,反演结果越准确;高斯牛顿法和快速模拟退火法计算精度相同,都可以得到高精度的横波速度径向分布;快速模拟退火法的计算效率略低于高斯牛顿法,但其收敛性对初始值依赖更小,实际处理中应选择快速模拟退火法.     

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.     

航空电磁拟三维模型空间约束反演

为了克服时间域航空电磁数据单点反演结果中常见的电阻率或层厚度横向突变造成数据难以解释的问题,通过引入双向约束实现航空电磁拟三维空间约束反演.除考虑沿测线方向相邻测点之间的横向约束外,同时还考虑了垂直测线方向测点在空间上的相互约束.为此,首先设计拟三维模型中固定层厚和可变层厚两种空间约束反演方案,然后通过在目标函数中引入沿测线和垂直测线方向上的模型参数约束矩阵,并使用L-BFGS算法使目标函数最小化,获得最优拟三维模型空间反演解.基于理论模型和实测数据反演,对单点反演与两种空间约束反演方案的有效性进行比较,证明本文空间约束反演算法对于噪声的压制效果好,反演的界面连续光滑,同时内存需求和反演时间少,是一种快速有效的反演策略.     

基于交叉梯度约束的时间域激发极化法二维同步反演

时间域激发极化法(Time-domain induced polarization method,简称为TDIP)已有的反演算法采用的是分步反演的思路,即先由视电阻率资料反演电阻率,固定电阻率再由视极化率资料反演极化率,这样就存在极化率结果严重依赖于电阻率反演结果的问题.为了有效解决这一问题,本文实现了TDIP二维数据空间分步反演算法,提出了基于交叉梯度约束的TDIP二维同步反演策略,实现了交叉梯度约束的电阻率和极化率二维同步反演算法.分别用电阻率和极化率结构一致和不一致的二维模型合成数据进行了分步和同步反演试算,对不同模型试算结果进行了对比分析.结果表明:对于电阻率和极化率结构一致和不一致模型,同步反演结果比分步反演结果能更好地确定异常体的空间分布范围,反演得到的电阻率和极化率值更接近真值.理论模型算例表明本文提出的同步反演算法有效解决了分步反演的问题,优于分步反演算法,具有更好的实用性.     

Two-dimensional joint inversion of radiomagnetotelluric and direct current resistivity data

An algorithm for the two-dimensional (2D) joint inversion of radiomagnetotelluric and direct current resistivity data was developed. This algorithm can be used for the 2D inversion of apparent resistivity data sets collected by multi-electrode direct current resistivity systems for various classical electrode arrays (Wenner, Schlumberger, dipole-diplole, pole-dipole) and radiomagnetotelluric measurements jointly. We use a finite difference technique to solve the Helmoltz and Poisson equations for radiomagnetotelluric and direct current resistivity methods respectively. A regularized inversion with a smoothness constrained stabilizer was employed to invert both data sets. The radiomagnetotelluric method is not particularly sensitive when attempting to resolve near-surface resistivity blocks because it uses a limited range of frequencies. On the other hand, the direct current resistivity method can resolve these near-surface blocks with relatively greater accuracy. Initially, individual and joint inversions of synthetic radiomagnetotelluric and direct current resistivity data were compared and we demonstrated that the joint inversion result based on this synthetic data simulates the real model more accurately than the inversion results of each individual method. The developed 2D joint inversion algorithm was also applied on a field data set observed across an active fault located close to the city of Kerpen in Germany. The location and depth of this fault were successfully determined by the 2D joint inversion of the radiomagnetotelluric and direct current resistivity data. This inversion result from the field data further validated the synthetic data inversion results.