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.     

Joint parameter estimation from magnetic resonance and vertical electric soundings using a multi‐objective genetic algorithm

Magnetic resonance sounding (MRS) has increasingly become an important method in hydrogeophysics because it allows for estimations of essential hydraulic properties such as porosity and hydraulic conductivity. A resistivity model is required for magnetic resonance sounding modelling and inversion. Therefore, joint interpretation or inversion is favourable to reduce the ambiguities that arise in separate magnetic resonance sounding and vertical electrical sounding (VES) inversions. A new method is suggested for the joint inversion of magnetic resonance sounding and vertical electrical sounding data. A one‐dimensional blocky model with varying layer thicknesses is used for the subsurface discretization. Instead of conventional derivative‐based inversion schemes that are strongly dependent on initial models, a global multi‐objective optimization scheme (a genetic algorithm [GA] in this case) is preferred to examine a set of possible solutions in a predefined search space. Multi‐objective joint optimization avoids the domination of one objective over the other without applying a weighting scheme. The outcome is a group of non‐dominated optimal solutions referred to as the Pareto‐optimal set. Tests conducted using synthetic data show that the multi‐objective joint optimization approximates the joint model parameters within the experimental error level and illustrates the range of trade‐off solutions, which is useful for understanding the consistency and conflicts between two models and objectives. Overall, the Levenberg‐Marquardt inversion of field data measured during a survey on a North Sea island presents similar solutions. However, the multi‐objective genetic algorithm method presents an efficient method for exploring the search space by producing a set of non‐dominated solutions. Borehole data were used to provide a verification of the inversion outcomes and indicate that the suggested genetic algorithm method is complementary for derivative‐based inversions.     

Practical issues in imaging hydraulic conductivity through hydraulic tomography

Hydraulic tomography has been developed as an alternative to traditional geostatistical methods to delineate heterogeneity patterns in parameters such as hydraulic conductivity (K) and specific storage (S(s)). During hydraulic tomography surveys, a large number of hydraulic head data are collected from a series of cross-hole tests in the subsurface. These head data are then used to interpret the spatial distribution of K and S(s) using inverse modeling. Here, we use the Sequential Successive Linear Estimator (SSLE) of Yeh and Liu (2000) to interpret synthetic pumping test data created through numerical simulations and real data generated in a laboratory sandbox aquifer to obtain the K tomograms. Here, we define "K tomogram" as an image of K distribution of the subsurface (or the inverse results) obtained via hydraulic tomography. We examine the influence of signal-to-noise ratio and biases on results using inverse modeling of synthetic and real cross-hole pumping test data. To accomplish this, we first show that the pumping rate, which affects the signal-to-noise ratio, and the order of data included into the SSLE algorithm both have large impacts on the quality of the K tomograms. We then examine the role of conditioning on the K tomogram and find that conditioning can improve the quality of the K tomogram, but can also impair it, if the data are of poor quality and conditioning data have a larger support volume than the numerical grid used to conduct the inversion. Overall, these results show that the quality of the K tomogram depends on the design of pumping tests, their conduct, the order in which they are included in the inverse code, and the quality as well as the support volume of additional data that are used in its computation.     

Three-dimensional geostatistical inversion of flowmeter and pumping test data

We jointly invert field data of flowmeter and multiple pumping tests in fully screened wells to estimate hydraulic conductivity using a geostatistical method. We use the steady-state drawdowns of pumping tests and the discharge profiles of flowmeter tests as our data in the inference. The discharge profiles need not be converted to absolute hydraulic conductivities. Consequently, we do not need measurements of depth-averaged hydraulic conductivity at well locations. The flowmeter profiles contain information about relative vertical distributions of hydraulic conductivity, while drawdown measurements of pumping tests provide information about horizontal fluctuation of the depth-averaged hydraulic conductivity. We apply the method to data obtained at the Krauthausen test site of the Forschungszentrum Jülich, Germany. The resulting estimate of our joint three-dimensional (3D) geostatistical inversion shows an improved 3D structure in comparison to the inversion of pumping test data only.     

Cost Comparisons of Aquifer Heterogeneity Characterization Methods

The characterization of heterogeneity in hydraulic conductivity (K) is a major challenge for subsurface remediation projects. There are a number of field studies that compare the K estimates obtained using various techniques, but to our knowledge, no field‐based studies exists that compare the performance of estimated K heterogeneity fields or the associated characterization costs. In this paper, we compare the costs of characterizing the three‐dimensional K heterogeneity and its uncertainty estimates of a glaciofluvial aquifer‐aquitard sequence at a 15 m × 15 m × 18 m field site situated on the University of Waterloo campus. We compare geostatistical analysis of high resolution permeameter K data obtained from repacked core samples in five boreholes and hydraulic tomography analysis of four pumping tests consisting of up to 41 monitoring points per test. Aside from the comparison of costs, we also assess the performance of each method by predicting several pumping tests. Our analysis reveals that hydraulic tomography is somewhat more costly than the geostatistical analysis of high resolution permeameter K data due to the higher capital costs associated with the method. However, the equipment may be reused at other sites; hence these costs may be recovered over the life of the equipment. More significantly, hydraulic tomography is able to capture the most important features of the aquifer‐aquitard sequence leading to more accurate predictions of independent pumping tests. This suggests that more robust remediation systems may be designed if site characterization is performed with hydraulic tomography.     

介质密度反演偏导矩阵的精确计算方法

实现反演偏导矩阵的计算是基于导数最优化反演方法的关键,然而目前的地震反演几乎都是基于Zoeppritz方程近似实现的,使计算精度和适应范围受到限制.本文利用Zoeppritz方程建立了反射系数对地层介质密度比偏导方程,导出了Zoeppritz方程矩阵元对介质密度比的导数.通过求解偏导方程获得了反射系数对介质密度比偏导数的精确计算(考虑了速度中含介质密度的问题).利用数值算例分析了反射系数对介质密度比偏导数的变化特点.本文采用直接解法求解偏导矩阵方程组,获得了快的计算速度和高的计算精度,为实现地层介质密度反演(包括大角度反演)提供了偏导矩阵的计算方法.     

大地电磁三维数据空间反演并行算法研究

目前大地电磁三维反演实际应用的主要问题是计算效率低.在对大地电磁三维数据空间反演算法进行深入分析的基础上,本文提出了基于频点和矩阵划分的大粒度并行反演方案和具体实现步骤,并在曙光TC5000A高性能计算平台上实现了基于MPI的大地电磁三维数据空间反演并行算法.该算法实现了包括三维正演、灵敏度矩阵、叉积矩阵以及模型改正量的并行执行,不仅计算效率高,而且每个节点机上灵敏度矩阵的存储空间只需原来微机上的2/N(N是参加并行计算的节点机个数),大大地减少了内存开销.通过两个理论模型合成的数据对实现的三维数据空间反演并行算法进行试算,对比分析了多个节点机下程序的执行效率.测试结果表明,所实现的三维数据空间反演并行算法是可行的、高效的,与单机相比,不仅可以提高运行速度,缩短计算时间,而且还可以扩大计算规模,极大地推动了大地电磁三维反演的实用化.     

Hydraulic Tomography: Continuity and Discontinuity of High‐K and Low‐K Zones

Hydraulic tomography is an emerging field and modeling method that provides a continuous hydraulic conductivity (K) distribution for an investigated region. Characterization approaches that rely on interpolation between one‐dimensional (1D) profiles have limited ability to accurately identify high‐K channels, juxtapositions of lenses with high K contrast, and breaches in layers or channels between such profiles. However, locating these features is especially important for groundwater flow and transport modeling, and for design and operation of in situ remediation in complex hydrogeologic environments. We use transient hydraulic tomography to estimate 3D K in a volume of 15‐m diameter by 20‐m saturated thickness in a highly heterogeneous unconfined alluvial (clay to sand‐and‐gravel) aquifer with a K range of approximately seven orders of magnitude at an active industrial site in Assemini, Sardinia, Italy. A modified Levenberg‐Marquardt algorithm was used for geostatistical inversion to deal with the nonlinear nature of the highly heterogeneous system. The imaging results are validated with pumping tests not used in the tomographic inversion. These tests were conducted from three of five clusters of continuous multichannel tubing (CMTs) installed for observation in the tomographic testing. Locations of high‐K continuity and discontinuity, juxtaposition of very high‐K and very low‐K lenses, and low‐K “plugs” are evident in regions of the investigated volume where they likely would not have been identified with interpolation from 1D profiles at the positions of the pumping well and five CMT clusters. Quality assessment methods identified a suspect high‐K feature between the tested volume and a lateral boundary of the model.     

Accelerating the T‐matrix approach to seismic full‐waveform inversion by domain decomposition

A seismic variant of the distorted Born iterative inversion method, which is commonly used in electromagnetic and acoustic (medical) imaging, has been recently developed on the basis of the T‐matrix approach of multiple scattering theory. The distorted Born iterative method is consistent with the Gauss–Newton method, but its implementation is different, and there are potentially significant computational advantages of using the T‐matrix approach in this context. It has been shown that the computational cost associated with the updating of the background medium Green functions after each iteration can be reduced via the use of various linearisation or quasi‐linearisation techniques. However, these techniques for reducing the computational cost may not work well in the presence of strong contrasts. To deal with this, we have now developed a domain decomposition method, which allows one to decompose the seismic velocity model into an arbitrary number of heterogeneous domains that can be treated separately and in parallel. The new domain decomposition method is based on the concept of a scattering‐path matrix, which is well known in solid‐state physics. If the seismic model consists of different domains that are well separated (e.g., different reservoirs within a sedimentary basin), then the scattering‐path matrix formulation can be used to derive approximations that are sufficiently accurate but far more speedy and much less memory demanding because they ignore the interaction between different domains. However, we show here that one can also use the scattering‐path matrix formulation to calculate the overall T‐matrix for a large model exactly without any approximations at a computational cost that is significantly smaller than the cost associated with an exact formal matrix inversion solution. This is because we have derived exact analytical results for the special case of two interacting domains and combined them with Strassen's formulas for fast recursive matrix inversion. To illustrate the fact that we have accelerated the T‐matrix approach to full‐waveform inversion by domain decomposition, we perform a series of numerical experiments based on synthetic data associated with a complex salt model and a simpler two‐dimensional model that can be naturally decomposed into separate upper and lower domains. If the domain decomposition method is combined with an additional layer of multi‐scale regularisation (based on spatial smoothing of the sensitivity matrix and the data residual vector along the receiver line) beyond standard sequential frequency inversion, then one apparently can also obtain stable inversion results in the absence of ultra‐low frequencies and reduced computation times.     

地震层析成像LSQR算法的并行化

讨论了地震层析成像的LSQR算法(最小二乘QR分解). 在建立偏导数矩阵方程组时,对区内地震在方程中保留震源项,引入正交投影算子进行参数分离,对区外远震采用传统的平滑处理方式,用LSQR法求解联立的方程组. 由于区内地震的正交分解处理和区外远震的平滑处理,使得偏导数矩阵中的非零元素成倍增加,对于大型反演问题,这些非零元素常常达到几十GB到几百GB的数量级,巨量的内存占用成为LSQR算法的瓶颈. 针对这一问题,本文研究了偏导数矩阵中非零元素的分布规律,设计出合理的存储结构,采用分布式存储进行矩阵计算,提出了LSQR算法的并行化方案,并在联想深腾6800超级计算机上实现. 导出了LSQR算法的并行效率估算公式. 对两个地区的实际地震层析成像数据进行了效率测试.     

时间域航空电磁2.5维非线性共轭梯度反演

对于时间域航空电磁法二维和三维反演来说,最大的困难在于有效的算法和大的计算量需求.本文利用非线性共轭梯度法实现了时间域航空电磁法2.5维反演方法,着重解决了迭代反演过程中灵敏度矩阵计算、最佳迭代步长计算、初始模型选取等问题.在正演计算中,我们采用有限元法求解拉式傅氏域中的电磁场偏微分方程,再通过逆拉氏和逆傅氏变换高精度数值算法得到时间域电磁响应.在灵敏度矩阵计算中,采用了基于拉式傅氏双变换的伴随方程法,时间消耗只需计算两次正演,从而节约了大量计算时间.对于最佳步长计算,二次插值向后追踪法能够保证反演迭代的稳定性.设计两个理论模型,检验反演算法的有效性,并讨论了选择不同初始模型对反演结果的影响.模型算例表明:非线性共轭梯度方法应用于时间域航空电磁2.5维反演中稳定可靠,反演结果能够有效地反映地下真实电性结构.当选择的初始模型电阻率值与真实背景电阻率值接近时,能得到较好的反演结果,当初始模型电阻率远大于或远小于真实背景电阻率值时反演效果就会变差.     

基于不等式约束的最小二乘法三维电阻率反演及其算法优化

基于光滑约束的最小二乘法是三维电阻率反演的主要方法,但该方法在某些情况下存在着多解性较强的问题,且普遍耗时较长,严重制约了三维反演方法的推广与发展.为改善上述问题,将表征模型参数变化范围的不等式约束作为先验信息引入最小二乘线性反演方法中,有效地改善了反演结果的精度,降低了反演的多解性问题.为了解决耗时较长的问题,基于预条件共轭梯度(PCG)算法和Cholesky分解法的特点提出了一套优化三维电阻率反演计算效率的计算方案.在该方案中,Cholesky分解法被用来求解敏感度矩阵计算中的多个点源场的正演问题,Cholesky分解法只需对总体系数矩阵进行一次分解,然后对不同的右端向量进行回代即可.将预条件共轭梯度法引入到三维电阻率反演方程的求解中,将雅可比迭代中的对角阵作为预处理矩阵,其具有求逆方便、无需内存空间的特点,有效地加快了收敛速度.对合成数据以及实测数据的反演算例表明,借助不等式约束和反演效率优化方案,最小二乘反演方法可得到较为精确的反演结果,有效地提高了反演计算效率,具有良好的推广前景.     

基于数据空间和稀疏约束的三维重力和重力梯度数据联合反演

随着重力和重力梯度测量技术的日趋成熟,基于重力和重力梯度数据的反演技术得到了广泛关注.针对反演多解性严重、计算效率低和内存消耗大等难点问题,本文开展了三维重力和重力梯度数据的联合反演研究,该方法结合重力和重力梯度两种数据,将L0范数正则化项加入到目标函数中,并在数据空间下采用改进的共轭梯度算法求解反演最优化问题.同时,...     

Frequency‐domain acoustic‐elastic coupled waveform inversion using the Gauss‐Newton conjugate gradient method

We developed a frequency‐domain acoustic‐elastic coupled waveform inversion based on the Gauss‐Newton conjugate gradient method. Despite the use of a high‐performance computer system and a state‐of‐the‐art parallel computation algorithm, it remained computationally prohibitive to calculate the approximate Hessian explicitly for a large‐scale inverse problem. Therefore, we adopted the conjugate gradient least‐squares algorithm, which is frequently used for geophysical inverse problems, to implement the Gauss‐Newton method so that the approximate Hessian is calculated implicitly. Thus, there was no need to store the Hessian matrix. By simultaneously back‐propagating multi‐components consisting of the pressure and displacements, we could efficiently extract information on the subsurface structures. To verify our algorithm, we applied it to synthetic data sets generated from the Marmousi‐2 model and the modified SEG/EAGE salt model. We also extended our algorithm to the ocean‐bottom cable environment and verified it using ocean‐bottom cable data generated from the Marmousi‐2 model. With the assumption of a hard seafloor, we recovered both the P‐wave velocity of complicated subsurface structures as well as the S‐wave velocity. Although the inversion of the S‐wave velocity is not feasible for the high Poisson's ratios used to simulate a soft seafloor, several strategies exist to treat this problem. Our example using multi‐component data showed some promise in mitigating the soft seafloor effect. However, this issue still remains open.     

Pitfalls of 1D inversion of small-loop electromagnetic data for detecting man-made objects

In order to interpret field data from small-loop electromagnetic (EM) instruments with fixed source–receiver separation, 1D inversion method is commonly used due to its efficiency with regard to computation costs. This application of 1D inversion is based on the assumption that small-offset broadband EM signals are insensitive to lateral resistivity variation. However, this assumption can be false when isolated conductive bodies such as man-made objects are embedded in the earth. Thus, we need to clarify the applicability of the 1D inversion method for small-loop EM data. In order to systematically analyze this conventional inversion approach, we developed a 2D EM inversion algorithm and verified this algorithm with a synthetic EM data set. 1D and 2D inversions were applied to synthetic and field EM data sets. The comparison of these inversion results shows that the resistivity distribution of the subsurface constructed by the 1D inversion approach can be distorted when the earth contains man-made objects, because they induce drastic variation of the resistivity distribution. By analyzing the integrated sensitivity of the small-loop EM method, we found that this pitfall of 1D inversion may be caused by the considerable sensitivity of the small-loop EM responses to lateral resistivity variation. However, the application of our 2D inversion algorithm to synthetic and field EM data sets demonstrate that the pitfall of 1D inversion due to man-made objects can be successfully alleviated. Thus, 2D EM inversion is strongly recommended for detecting conductive isolated bodies, such as man-made objects, whereas this approach may not always be essential for interpreting the EM field data.     

Geostatistical seismic Amplitude‐versus‐angle inversion

Seismic inversion plays an important role in reservoir modelling and characterisation due to its potential for assessing the spatial distribution of the sub‐surface petro‐elastic properties. Seismic amplitude‐versus‐angle inversion methodologies allow to retrieve P‐wave and S‐wave velocities and density individually allowing a better characterisation of existing litho‐fluid facies. We present an iterative geostatistical seismic amplitude‐versus‐angle inversion algorithm that inverts pre‐stack seismic data, sorted by angle gather, directly for: density; P‐wave; and S‐wave velocity models. The proposed iterative geostatistical inverse procedure is based on the use of stochastic sequential simulation and co‐simulation algorithms as the perturbation technique of the model parametre space; and the use of a genetic algorithm as a global optimiser to make the simulated elastic models converge from iteration to iteration. All the elastic models simulated during the iterative procedure honour the marginal prior distributions of P‐wave velocity, S‐wave velocity and density estimated from the available well‐log data, and the corresponding joint distributions between density versus P‐wave velocity and P‐wave versus S‐wave velocity. We successfully tested and implemented the proposed inversion procedure on a pre‐stack synthetic dataset, built from a real reservoir, and on a real pre‐stack seismic dataset acquired over a deep‐water gas reservoir. In both cases the results show a good convergence between real and synthetic seismic and reliable high‐resolution elastic sub‐surface Earth models.     

Radial basis function neural network for hydrologic inversion: an appraisal with classical and spatio-temporal geostatistical techniques in the context of site characterization

This paper investigates three techniques for spatial mapping and the consequential hydrologic inversion, using hydraulic conductivity (or transmissivity) and hydraulic head as the geophysical parameters of concern. The data for the study were obtained from the Waste Isolation and Pilot Plant (WIPP) site and surrounding area in the remote Chihuahuan Desert of southeastern New Mexico. The central technique was the Radial Basis Function algorithm for an Artificial Neural Network (RBF-ANN). An appraisal of its performance in light of classical and temporal geostatistical techniques is presented. Our classical geostatistical technique of concern was Ordinary Kriging (OK), while the method of Bayesian Maximum Entropy (BME) constituted an advanced, spatio-temporal mapping technique. A fusion technique for soft or inter-dependent data was developed in this study for use with the neural network. It was observed that the RBF-ANN is capable of hydrologic inversion for transmissivity estimation with features remaining essentially similar to that obtained from kriging. The BME technique, on the other hand, was found to reveal an ability to map localized lows and highs that were otherwise not as apparent in OK or RBF-ANN techniques.     

Development and application of a coupled-process parameter inversion model based on the maximum likelihood estimation method

The coupled flow-mass transport inverse problem is formulated using the maximum likelihood estimation concept. An evolutionary computational algorithm, the genetic algorithm, is applied to search for a global or near-global solution. The resulting inverse model allows for flow and transport parameter estimation, based on inversion of spatial and temporal distributions of head and concentration measurements. Numerical experiments using a subset of the three-dimensional tracer tests conducted at the Columbus, Mississippi site are presented to test the model's ability to identify a wide range of parameters and parametrization schemes. The results indicate that the model can be applied to identify zoned parameters of hydraulic conductivity, geostatistical parameters of the hydraulic conductivity field, angle of hydraulic conductivity anisotropy, solute hydrodynamic dispersivity, and sorption parameters. The identification criterion, or objective function residual, is shown to decrease significantly as the complexity of the hydraulic conductivity parametrization is increased. Predictive modeling using the estimated parameters indicated that the geostatistical hydraulic conductivity distribution scheme produced good agreement between simulated and observed heads and concentrations. The genetic algorithm, while providing apparently robust solutions, is found to be considerably less efficient computationally than a quasi-Newton algorithm.     

Preconditioned non‐linear conjugate gradient method for frequency domain full‐waveform seismic inversion

We present preconditioned non‐linear conjugate gradient algorithms as alternatives to the Gauss‐Newton method for frequency domain full‐waveform seismic inversion. We designed two preconditioning operators. For the first preconditioner, we introduce the inverse of an approximate sparse Hessian matrix. The approximate Hessian matrix, which is highly sparse, is constructed by judiciously truncating the Gauss‐Newton Hessian matrix based on examining the auto‐correlation and cross‐correlation of the Jacobian matrix. As the second preconditioner, we employ the approximation of the inverse of the Gauss‐Newton Hessian matrix. This preconditioner is constructed by terminating the iteration process of the conjugate gradient least‐squares method, which is used for inverting the Hessian matrix before it converges. In our preconditioned non‐linear conjugate gradient algorithms, the step‐length along the search direction, which is a crucial factor for the convergence, is carefully chosen to maximize the reduction of the cost function after each iteration. The numerical simulation results show that by including a very limited number of non‐zero elements in the approximate Hessian, the first preconditioned non‐linear conjugate gradient algorithm is able to yield comparable inversion results to the Gauss‐Newton method while maintaining the efficiency of the un‐preconditioned non‐linear conjugate gradient method. The only extra cost is the computation of the inverse of the approximate sparse Hessian matrix, which is less expensive than the computation of a forward simulation of one source at one frequency of operation. The second preconditioned non‐linear conjugate gradient algorithm also significantly saves the computational expense in comparison with the Gauss‐Newton method while maintaining the Gauss‐Newton reconstruction quality. However, this second preconditioned non‐linear conjugate gradient algorithm is more expensive than the first one.     

Seismic envelope inversion: reduction of local minima and noise resistance

Waveform inversion met severe challenge in retrieving long‐wavelength background structure. We have proposed to use envelope inversion to recover the large‐scale component of the model. Using the large‐scale background recovered by envelope inversion as new starting model, we can get much better result than the conventional full waveform inversion. By comparing the configurations of the misfit functional between the envelope inversion and the conventional waveform inversion, we show that envelope inversion can greatly reduce the local minimum problem. The combination of envelope inversion and waveform inversion can deliver more faithful and accurate final result with almost no extra computation cost compared to the conventional full waveform inversion. We also tested the noise resistance ability of envelope inversion to Gaussian noise and seismic interference noise. The results showed that envelope inversion is insensitive to Gaussian noise and, to a certain extent, insensitive to seismic interference noise. This indicates the robustness of this method and its potential use for noisy data.