A comparison of different radiomagnetotelluric data inversion methods for buried waste sites

This study deals with two-dimensional (2D) inversions of synthetic and observed radiomagnetotelluric (RMT) data on typical buried conductive waste sites in Europe, and with the practical aspects of different inversion algorithms. In the inversion calculations, we used smoothing and L2-norm stabilizers and compared the results. The resolution of the geometry of the highly conductive waste site, in particular, was investigated. In the inversion with the L2-norm stabilizer, we used the least-squares solution with singular value decomposition (LSSVD) and conjugate gradient (CG), whereas only the conjugate gradient solver was used in the 2D-inversion with the smoothing stabilizer. The inversion results of the synthetic data showed a better resolution of the geometry of the highly conductive waste site when a L2-norm stabilizer was applied in the inversion; in particular, a better detection of the bottom of the waste deposit was achieved. Additional model studies were carried out using synthetic RMT data in order to investigate the 2D inversion of RMT data observed on a 3D structure; these studies showed that the use of TM mode data yields a better resolution of the structure than joint inversion of TE and TM modes.2D inversions of RMT data on a waste site near Cologne showed that the inversion of the TM mode could resolve well the geometry, especially the bottom of the waste site, if information about the background conductivity structure was available. In this case study, inversion with the L2-norm stabilizer produced a sharper image of the waste site than inversion with the smoothing stabilizer, as indicated also by the inversion study that used synthetic data.     

大地电磁法三维共轭梯度反演研究

Based on the analysis of the conjugate gradient algorithm, we implement a threedimensional （3D） conjugate gradient inversion algorithm with magnetotelluric impedance data. During the inversion process, the 3D conjugate gradient inversion algorithm doesn＇ t need to compute and store the Jacobian matrix but directly updates the model from the computation of the Jacobian matrix. Requiring only one forward and four pseudo-forward modeling applications per frequency to produce the model update at each iteration, this algorithm efficiently reduces the computation of the inversion. From a trial inversion with synthetic magnetotelluric data, the validity and stability of the 3D conjugate gradient inversion algorithm is verified.     

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.     

大地电磁二维各向异性反演及其在青藏高原北部的应用

地球介质的电性特征在不同深度上都可能表现出各向异性,识别电各向异性有助于深入理解地球内部物质状态及变形环境,促进对动力学模型的理解.本研究在已实现的任意电各向异性大地电磁有限差分正演算法的基础之上,通过构建常规的各向异性反演目标函数,完成对目标函数及其梯度的计算.借助成熟的非线性共轭梯度反演技术,实现了对大地电磁各向异性介质的反演过程.在此基础上,构建了二维方位各向异性理论模型,来验证反演过程的稳定性.理论模型的反演结果直观地显示了大地电磁各向异性反演对模型参数中垂直电导率恢复上的局限性;在各向异性结构分析上需要将电导率张量作为整体考虑.对青藏高原北部东昆仑—柴达木盆地西段的实测大地电磁数据进行反演,并对比早期二维各向同性反演结果发现,各向同性结构中位于祁漫塔格山脉下方上地幔顶部存在低阻异常带,在相同的位置,各向异性结构中表现出明显的方位各向异性,其各向异性低阻主轴电阻率指示的剪切带走向与地表造山带走向不一致,这暗示该处的应力环境复杂,除了受印度—欧亚板块碰撞控制外,还可能受到邻近阿尔金走滑断裂带左行剪切运动的影响.

     

基于交叉梯度耦合的大地电磁与地震走时资料三维联合反演

为了有效解决目前大地电磁和地震走时资料单方法反演结果一致性不好的问题,同时克服基于岩石不同物性参数间关系耦合约束联合反演的局限性,本文研究了基于交叉梯度耦合约束的大地电磁与地震走时资料的三维联合反演算法.以较为成熟的天然地震走时资料三维正反演和大地电磁三维正反演算法为基础,实现了具有共同的反演网格,以交叉梯度结构耦合约束,并能同时获得电阻率和速度模型的三维联合反演算法.分别利用单棱柱体模型和双棱柱体模型合成数据进行了联合反演试算.结果表明:无论是单棱柱体模型还是双棱柱体模型,联合反演结果比单独反演对异常体的空间形态都有更好的恢复,其中单棱柱体模型反演的异常体电阻率更接近于真实电阻率,双棱柱体模型的联合反演结果不仅消除了围岩的部分电阻率假异常,而且增强了对异常体深部速度结构特征的恢复程度.联合反演还能同时改善电阻率和速度反向变化异常体的单独反演结果,进一步证明交叉梯度耦合不依赖于岩石物性关系,而强调地下结构的相似性,具有更普遍的适用性.     

倾子资料三维共轭梯度反演研究

在对倾子响应和共轭梯度算法深入分析的基础上,我们实现了倾子资料三维共轭梯度反演算法.基于倾子资料的三维共轭梯度反演研究,探讨了利用倾子资料进行三维反演定量解释的方法.通过对理论模型合成数据进行反演试算,验证了所实现的倾子资料三维共轭梯度反演算法的有效性和稳定性.该反演算法可用于对大地电磁测深和地磁测深(地震地磁台站进行的三分量地磁观测)所整理出的倾子资料进行三维定量反演,获得地下三维模型的电阻率结构.     

Hybrid genetic algorithms in view of the evolution theories with application for the electrical sounding method

The classical genetic algorithm is a stochastic process which operates by natural selection. Although the algorithm may localize a point around the global minimum of the misfit function, it is not efficient at finding the precise solution. This paper suggests some hybrid genetic algorithms, derived from evolution theories, to overcome this problem. Firstly, sexual selection has been incorporated in the classical genetic algorithm to obtain a full representation of the Darwinist evolution concept. The simulation of sexual selection is performed by assigning a higher probability of surviving to some parameters that satisfy some algebraic relationships. This method is called the ‘marked constraints’ algorithm since it permits us to insert geological and geophysical constraints into the problem. The algorithm implementation is realized by progressively shrinking the parameter search space through successive generations. In this way, the genetic algorithm gains some degree of determinism. Secondly, since the evolution theory of Lamarck postulates that the acquired traits are passed on to the next generation, a hybrid use of the damped least‐squares method and the genetic algorithm is called Lamarckian inversion. Lamarckian inversion involves some improvement procedures that simulate the reduction of the misfit with the help of a derivative‐based method between two generations. Finally, although there is no correspondence in nature, Lamarckian and Darwinist evolution concepts are combined to strengthen the deterministic part of the solution algorithm. This is called the Lamarckian‐marked‐constraint algorithm. The merits and behaviours of the suggested algorithms are discussed using two examples. The first is a hypothetical example affected by a multiminima problem. The second examines the equivalence problem using vertical electrical sounding data.     

大地电磁全信息资料三维共轭梯度反演研究(英文)

在对张量阻抗数据、倾子数据和共轭梯度算法深入分析的基础上,我们实现了大地电磁全信息资料三维共轭梯度反演算法。基于全信息资料的三维共轭梯度反演研究,探讨了同时利用五个电磁场分量整理得到的大地电磁资料进行三维反演定量解释的方法以及全信息数据在三维反演中的作用。理论模型合成数据的反演结果表明,在三维反演中使用张量阻抗和倾子数据结合的全信息数据的反演结果优于只使用张量阻抗数据(或只使用倾子数据)的反演结果,提高了反演结果的分辨率和可信度。合成数据的反演算例也验证了所实现的大地电磁全信息资料三维共轭梯度反演算法的正确性和稳定性。     

Linearized seismic waveform inversion using a multiple re-weighted least-squares method with QR preconditioning

Linearized inversion methods such as Gauss‐Newton and multiple re‐weighted least‐squares are iterative processes in which an update in the current model is computed as a function of data misfit and the gradient of data with respect to model parameters. The main advantage of those methods is their ability to refine the model parameters although they have a high computational cost for seismic inversion. In the Gauss‐Newton method a system of equations, corresponding to the sensitivity matrix, is solved in the least‐squares sense at each iteration, while in the multiple re‐weighted least‐squares method many systems are solved using the same sensitivity matrix. The sensitivity matrix arising from these methods is usually not sparse, thus limiting the use of standard preconditioners in the solution of the linearized systems. For reduction of the computational cost of the linearized inversion methods, we propose the use of preconditioners based on a partial orthogonalization of the columns of the sensitivity matrix. The new approach collapses a band of co‐diagonals of the normal equations matrix into the main diagonal, being equivalent to computing the least‐squares solution starting from a partial solution of the linear system. The preconditioning is driven by a bandwidth L which can be interpreted as the distance for which the correlation between model parameters is relevant. To illustrate the benefit of the proposed approach to the reduction of the computational cost of the inversion we apply the multiple re‐weighted least‐squares method to the 2D acoustic seismic waveform inversion problem. We verify the reduction in the number of iterations in the conjugate'gradient algorithm as the bandwidth of the preconditioners increases. This effect reduces the total computational cost of inversion as well.     

Determination of the regularization level of truncated singular-value decomposition inversion: The case of 1D inversion of MT data

The total mean‐square error (MSE) of the estimated model, defined as the sum of the standard model variance and the bias variance, is used to define the truncation level of the singular‐value decomposition to give a reasonable balance between model resolution and model variance. This balance is determined largely by the data and no further assumptions are necessary except that the bias terms are estimated sufficiently well. This principle has been tested on the 1D magnetotelluric inverse problem with special emphasis on high‐frequency radio magnetotelluric (RMT) data. Simulations clearly demonstrate that the method provides a good balance between resolution and variance. Starting from a homogeneous half‐space, the best solution is sought for a fixed set of singular values. The model variance is estimated from the sum of the inverse eigenvalues squared, up to a certain threshold, and the bias variance is estimated from the model projections on the remaining eigenvectors. By varying the threshold, the minimum of the MSE is found for an increasing number of fixed singular values until the number of active singular values becomes greater than or equal to the estimated number. As a side‐effect, the depth of penetration of a given set of measurements can be estimated very efficiently by simply noting at which depth the final model deviates little from the starting homogeneous half‐space model. A suite of synthetic data is inverted and an example of inversion of one site is shown to illustrate how the truncation is carried out as the non‐linear inversion process proceeds. A field example with a profile across a plume of contaminated groundwater in the Netherlands shows good agreement with the electrical resistivity obtained in a nearby borehole.     

一种并行的大地电磁场非线性共轭梯度三维反演方法

本文改进并验证了大地电磁测深数据的三维反演算法和并行计算程序,程序对计算机物理内存和CPU速度及数量要求较低,使普通家用机进行三维反演计算成为可能.本文在Newman和Alumbaugh(2000)提出的三维非线性共轭梯度算法和Rodi和Mackie(2001)给出的大地电磁场二维NLCG反演预处理方法的基础上实现了大地电磁场NLCG三维反演算法,改进了的预处理方法,将反演计算对初始模型的依赖性降到最低,并且通过理论模型验证了程序的正确性,并根据日本KAYABE地区实测数据的反演结果验证了算法的实用性.     

二维大地电磁尖锐边界反演研究

尖锐边界反演是大地电磁及其他地球物理反演中的一个较难解决的问题,目前有不少专家正对这方面进行研究.我们在分析了Zhdanov（2004）提出的基于最小支撑泛函聚焦反演方法的基础上,在目标函数中引入对角梯度支撑以改善倾斜电性分界面的反演效果.同时,为了提高计算效率,本文还使用互易定义计算灵敏度矩阵及使用正则化共轭梯度法优化目标函数.通过模型试验并与未引入对角梯度支撑的反演结果进行了比较,表明本文的方法反演效果更好而且算法稳定和实用.     

基于有限差分正演的带地形三维大地电磁反演方法

本研究实现了一套基于有限差分（FD）方法的大地电磁测深数据带地形三维反演算法及代码.其中,在大地电磁场正演数值模拟方面,开发了起伏地形条件下基于交错网格剖分、有限差分方法的大地电磁测深三维正演代码;在满足平面波场假设的前提下,使用长方体网格剖分模拟三维起伏地形,实现了带地形三维正演计算;并设计理论模型进行试算,经试算结果与前人的有限元法计算结果对比,验证了所研发的带地形三维正演计算的正确性与可靠性.在反演方面,本研究基于非线性共轭梯度方法编写了大地电磁测深带地形三维反演代码,试验了不同的共轭梯度搜索因子β,避免了目标函数对海森矩阵（参数二次导数矩阵）的显式计算和存储,初步实现了大地电磁资料的带地形三维反演.最后,对一系列理论模型进行正演计算,利用其生成的合成数据模拟实测数据进行反演,并与现有的不带地形大地电磁测深三维反演结果比较,检验了所研发的带地形三维反演计算的可靠性与稳定性.     

On regularized time varying gravity field models based on grace data and their comparison with hydrological models

Determination of spherical harmonic coefficients of the Earth’s gravity field is often an ill-posed problem and leads to solving an ill-conditioned system of equations. Inversion of such a system is critical, as small errors of data will yield large variations in the result. Regularization is a method to solve such an unstable system of equations. In this study, direct methods of Tikhonov, truncated and damped singular value decomposition and iterative methods of ν, algebraic reconstruction technique, range restricted generalized minimum residual and conjugate gradient are used to solve the normal equations constructed based on range rate data of the gravity field and climate experiment (GRACE) for specific periods. Numerical studies show that the Tikhonov regularization and damped singular value decomposition methods for which the regularization parameter is estimated using quasioptimal criterion deliver the smoothest solutions. Each regularized solution is compared to the global land data assimilation system (GLDAS) hydrological model. The Tikhonov regularization with L-curve delivers a solution with high correlation with this model and a relatively small standard deviation over oceans. Among iterative methods, conjugate gradient is the most suited one for the same reasons and it has the shortest computation time.     

大地电磁相位张量和倾子的三维主轴各向异性反演

地壳和上地幔普遍存在的电性各向异性是大地电磁测深方法在实践应用中的关键和难点问题,而由近地表局部地电异常引起的电流型畸变可进一步加剧各向异性反演的多解性和不确定性.为压制电流型畸变影响,本文实现了一种利用大地电磁相位张量和倾子数据获取地下介质主轴电性各向异性结构的三维反演算法.该算法正演问题的求解采用交错采样有限差分法,反演问题的求解采用非线性共轭梯度法.理论模型测试表明,相位张量和倾子的三维各向异性反演不仅恢复了真实的各向异性异常分布特征,且压制了为拟合电流型畸变出现的假异常.由于相位张量和倾子数据不含幅值信息,先验模型选取对反演模型电阻率值的恢复效果影响较大,但初始模型的改变基本不影响反演结果.为进一步验证该算法的有效性,将其应用于USArray项目在美国西北部采集的109个长周期大地电磁实测数据.反演结果表明,相对于阻抗数据,相位张量和倾子数据的各向异性反演结果对浅层结构的刻画更为准确,与地表地质特征更为吻合.

     

融入宽范围物性约束策略的大地电磁和重力的Gramian约束联合反演

宽范围物性约束方式可以促进物性参数在一定范围内实现耦合, 具有一定的容错性, 实现简单.但目前该方式只在全局优化算法联合反演中得以实现, 如何将其推广至梯度优化算法联合反演至今没有明确的策略.此外, Gramian约束虽然对先验信息依赖低, 但是对明确的物性关联信息利用率低.基于此, 本文提出了适用于梯度优化算法的宽范围物性约束的策略, 即"岩石物性关联+范围约束+耦合项", 并实现了基于宽范围物性约束策略的大地电磁测深(MT)和重力共轭梯度Gramian约束联合反演新算法.模型试验表明, 本文提出的算法可以提高先验物性关联信息的利用率, 且有效限定物性耦合的范围.我们将本文方法应用于安徽省茶亭铜金矿床实际资料处理之中, 联合反演结果清晰地揭示了铜金矿床的空间分布.     

龙门山断裂带南段地壳电性特征: 来自速度结构约束下大地电磁反演的证据

龙门山断裂带区域内大震频发、构造运动活跃, 其地史演化与隆升动力学机制一直未有定论.2013年4月20日雅安大地震后, 于小金至雅安段布设了一条垂直于龙门山构造带长约200 km的大地电磁剖面.为了克服龙门山地形起伏影响, 同时增强大地电磁反演的垂向分辨率, 本文采用宽角折/反射地震走时反演获得的速度模型作为结构约束, 通过交叉梯度项引入大地电磁非线性共轭梯度二维反演, 实现了考虑地形的速度结构约束大地电磁二维算法.设计地堑-地垒模型进行合成数据反演试算, 结果显示基于速度结构约束的带地形二维大地电磁反演算法能够减少由地形影响引起的假异常, 同时对异常体边界轮廓的刻画更加清晰.实测数据未加地震资料约束的大地电磁二维反演结果, 具有明显"低-高-低"的宏观电性特征, 与前人结果基本吻合.本文将新算法应用于实测数据反演的结果表明: 松潘—甘孜地块下方包含两部分低阻异常带, 并有相互连通趋势.西侧低阻体埋深更深, 分布于15~45 km范围内, 顶部有通道延伸至地表; 东侧低阻体相对变浅, 向东逐渐延伸至高阻异常体之上, 具有向扬子地体逆冲推覆的趋势.龙门山构造带薄皮盖层下方大规模高阻异常体厚约50 km, 以Moho面为底界, 且呈现北西倾向的特征, 与人工地震资料推测的龙门山三条主断裂的深部延伸的产状一致.     

Three-Dimensional Magnetotelluric Inversion: An Introductory Guide for Developers and Users

In the last few decades, the demand for three-dimensional (3-D) inversions for magnetotelluric data has significantly driven the progress of 3-D codes. There are currently a lot of new 3-D inversion and forward modeling codes. Some, such as the WSINV3DMT code of the author, are available to the academic community. The goal of this paper is to summarize all the important issues involving 3-D inversions. It aims to show how inversion works and how to use it properly. In this paper, I start by describing several good reasons for doing 3-D inversion instead of 2-D inversion. The main algorithms for 3-D inversion are reviewed along with some comparisons of their advantages and disadvantages. These algorithms are the classical Occam’s inversion, the data space Occam’s inversion, the Gauss–Newton method, the Gauss–Newton with the conjugate gradient method, the non-linear conjugate gradient method, and the quasi-Newton method. Other variants are based on these main algorithms. Forward modeling, sensitivity calculations, model covariance and its parallel implementation are all necessary components of inversions and are reviewed here. Rules of thumb for performing 3-D inversion are proposed for the benefit of the 3-D inversion novice. Problems regarding 3-D inversions are discussed along with suggested topics for future research for the developers of the next decades.     

Assumptions and goals for least squares migration

Least squares migration uses the assumption that, if we have an operator that can create data from a reflectivity function, the optimal image will predict the actual recorded data with minimum square error. For this assumption to be true, it is also required that: (a) the prediction operator must be error-free, (b) model elements not seen by the operator should be constrained by other means and (c) data weakly predicted by the operator should make limited contribution to the solution. Under these conditions, least squares migration has the advantage over simple migration of being able to remove interference between different model components. Least squares migration does that by de-convolving or inverting the so-called Hessian operator. The Hessian is the cascade of forward modelling and migration; for each image point, it computes the effects of interference from other image points (point-spread function) given the actual recording geometry and the subsurface velocity model. Because the Hessian contains illumination information (along its diagonal), and information about the model cross-correlation produced by non-orthogonality of basis functions, its inversion produces illumination compensation and increases resolution. In addition, sampling deficiencies in the recording geometry map to the Hessian (both diagonal and non-diagonal elements), so least squares migration has the potential to remove sampling artefacts as well. These (illumination compensation, resolution and mitigating recording deficiencies) are the three main goals of least squares migration, although the first one can be achieved by cheaper techniques. To invert the Hessian, least squares migration relies on the residual errors during iterations. Iterative algorithms, like conjugate gradient and others, use the residuals to calculate the direction and amplitudes (gradient and step size) of the necessary corrections to the reflectivity function or model. Failure of conditions (a), (b) or (c) leads the inversion to calculate incorrect model updates, which translate to noise in the final image. In this paper, we will discuss these conditions for Kirchhoff migration and reverse time migration.     

(T2,D)二维核磁共振测井混合反演方法与参数影响分析

一维核磁共振(1D NMR)测井技术在流体识别中具有一定的局限性.二维核磁共振(2D NMR)测井能同时测量到多孔介质中横向弛豫时间(T₂)和扩散系数(D)等信息,利用这两个参数区分流体性质,较一维核磁共振测井技术具有明显的优越性.针对梯度场下的2D NMR测井弛豫机理和数学模型,提出了适用于求解大型稀疏矩阵方程的反演方法-基于非负最小二乘法(LSQR)和截断奇异值分解(TSVD)法的混合算法.为验证方法的有效性,先根据多回波观测模式合成回波串数据,然后再用混合反演算法进行反演,反演得到横向弛豫时间(T₂)和扩散系数(D),并构建T₂-D二维谱图.结果对比表明,该混合反演算法得到的T₂-D二维谱与流体模型一致性好,计算精度均比单一反演方法有较大改善,表明该混合反演方法可用于油气储集层2D NMR测井的反演和流体识别.此外,分别对油水同层和气水同层模型进行了正演模拟和反演实验, 系统考察了不同磁场梯度、不同回波间隔组合对反演效果的影响,为2D NMR参数设计提供依据.