用改进的光滑约束最小二乘正交分解法实现电阻率三维反演

对三维电阻率反演问题进行了深入研究,提供了一种利用地表观测数据实现三维反演的实用算法.该方法应用有限差分求正演解,并通过对粗糙度矩阵元素进行适当改进,使之适用于各种情况下粗糙度矩阵的求取,进而建立在模型的总粗糙度极小条件下的反演方程.对反演方程采用收敛速度快且稳定的最小二乘正交分解（LSQR）法进行迭代求解,在迭代求解过程中只需利用偏导数矩阵和其转置矩阵乘以一个向量的结果,回避了直接求偏导数矩阵的繁琐计算,节省了内存,加快了反演的计算速度.不同的计算实例表明上述方法是求解大规模三维电阻率反演问题的有效方法.     

An efficient non-linear least-squares 1D inversion scheme for resistivity and IP sounding data

Non-linear least-squares inversion operates iteratively by updating the model parameters in each step by a correction vector which is the solution of a set of normal equations. Inversion of geoelectrical data is an ill-posed problem. This and the ensuing suboptimality restrict the initial model to being in the near vicinity of the true model. The problem may be reduced by introducing damping into the system of equations. It is shown that an appropriate choice of the damping parameter obtained adaptively and the use of a conjugate-gradient algorithm to solve the normal equations make the 1D inversion scheme efficient and robust. The scheme uses an optimal damping parameter that is dependent on the noise in the data, in each iterative step. The changes in the damping and relative residual error with iteration number are illustrated. A comparison of its efficacy over the conventional Marquardt and simulated annealing methods, tested on Inman's model, is made. Inversion of induced polarization (IP) sounding is obtained by inverting twice (true and modified) DC apparent resistivity data. The inversion of IP data presented here is generic and can be applied to any of the IP observables, such as chargeability, frequency effect, phase, etc., as long as these observables are explicitly related to the DC apparent resistivity. The scheme is used successfully in inverting noise-free and noisy synthetic data and field data taken from the published literature.     

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

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

Assessment of the reliability of 2D inversion of apparent resistivity data

The reliability of inversion of apparent resistivity pseudosection data to determine accurately the true resistivity distribution over 2D structures has been investigated, using a common inversion scheme based on a smoothness‐constrained non‐linear least‐squares optimization, for the Wenner array. This involved calculation of synthetic apparent resistivity pseudosection data, which were then inverted and the model estimated from the inversion was compared with the original 2D model. The models examined include (i) horizontal layering, (ii) a vertical fault, (iii) a low‐resistivity fill within a high‐resistivity basement, and (iv) an upfaulted basement block beneath a conductive overburden. Over vertical structures, the resistivity models obtained from inversion are usually much sharper than the measured data. However, the inverted resistivities can be smaller than the lowest, or greater than the highest, true model resistivity. The substantial reduction generally recorded in the data misfit during the least‐squares inversion of 2D apparent resistivity data is not always accompanied by any noticeable reduction in the model misfit. Conversely, the model misfit may, for all practical purposes, remain invariant for successive iterations. It can also increase with the iteration number, especially where the resistivity contrast at the bedrock interface exceeds a factor of about 10; in such instances, the optimum model estimated from inversion is attained at a very low iteration number. The largest model misfit is encountered in the zone adjacent to a contact where there is a large change in the resistivity contrast. It is concluded that smooth inversion can provide only an approximate guide to the true geometry and true formation resistivity.     

三维大地电磁自适应正则化有限内存拟牛顿反演

有别于传统基于梯度信息的反演方法在正则化约束中用总梯度逼近海塞逆矩阵的技术,本文将正则化约束问题的数据拟合项和模型光滑项分开考虑,只利用数据拟合函数的梯度信息对数据拟合项的海塞矩阵进行逼近,通过求解类高斯牛顿下降方向方程得到不依赖前几次迭代正则化因子的更精确下降方向,在求解当前迭代下降方向的过程中,通过保证右端项中两个向量的二范数在同一数量级的原则,实现了正则化因子的自动更新.对理论模型的试算表明这种自适应正则化反演方案可以在拟牛顿反演框架下基本达到OCCAM的算法稳定性,反演结果对初始模型依赖性较小,同时又无需在一次迭代中多次搜索最佳正则化因子.本文还基于此算法讨论了大地电磁各参数对于反演结果的影响,由于本文的反演结果能得到充分的正则化约束,因而在此框架下讨论阻抗和倾子在反演中的作用相对更为客观.

     

ABP法在高密度电阻率法反演中的应用

非线性反演方法作为地球物理反演的一个重要分支,在地球物理反演中发挥着特有的作用.近年来学者对非线性联合反演研究较多,但目前仍未有实质性的研究进展;本文尝试利用BP(Back Propagation)神经网络优化方法与蚁群算法联合演算,实现高密度电阻率法的电阻率二维非线性反演.通过两组模型的结果比较,BP与ABP法的反演较传统反演法优势较为突出,而且ABP(Ant colony optimization-Back Propagation)方法明显优于BP神经网络反演法,它可以克服BP神经网络反演方法的不足、减少迭代次数、节约计算时间,获得更好的反演结果.     

利用拟牛顿法和信赖域法联合反演震中分布与一维速度结构

本文将拟牛顿法与信赖域法相结合的反演方法用于地震震中和一维速度结构的联合反演. 通过数值试验比较发现：用拟牛顿法构造Hessian矩阵,虽然仅利用目标函数的一阶导数信息, 但由于保证了正定性和拟牛顿条件, 比解析求得Hessian矩阵更加高效稳定, 通过与信赖域方法的结合,使得该联合反演方法与经典的联合反演方法相比不仅具有较高的计算效率和较高的反演精度,同时也具有更强的总体收敛性.利用本文方法对台湾集集地震重新定位结果（纬度23.86°N, 经度120.80°E, 深度10km, 发震时刻1999年9月21日17:47:16.0）与前人主流结果一致.此外,该方法对不同初始速度模型以及不同的地震到时数据集均可给出较一致的结果,进一步表明该方法的优越性.本文还利用1998～2001年间华北遥测地震台网40个台站记录到的64次地震的971条P直达波和首波到时资料,采用本文提出的方法吩北京及其邻区的一维地壳P波速度结构进行联合反演, 结果表明该地区浅层地壳速度较低,中部地壳速度较高, 下部地壳速度逐渐递增.     

Robust inversion of vertical electrical sounding data using a multiple reweighted least-squares method

The root cause of the instability problem of the least-squares (LS) solution of the resistivity inverse problem is the ill-conditioning of the sensitivity matrix. To circumvent this problem a new LS approach has been investigated in this paper. At each iteration, the sensitivity matrix is weighted in multiple ways generating a set of systems of linear equations. By solving each system, several candidate models are obtained. As a consequence, the space of models is explored in a more extensive and effective way resulting in a more robust and stable LS approach to solving the resistivity inverse problem. This new approach is called the multiple reweighted LS method (MRLS). The problems encountered when using the L ₁- or L ₂-norm are discussed and the advantages of working with the MRLS method are highlighted. A five-layer earth model which generates an ill-conditioned matrix due to equivalence is used to generate a synthetic data set for the Schlumberger configuration. The data are randomly corrupted by noise and then inverted by using L ₂, L ₁ and the MRLS algorithm. The stabilized solutions, even though blurred, could only be obtained by using a heavy ridge regression parameter in L ₂- and L ₁-norms. On the other hand, the MRLS solution is stable without regression factors and is superior and clearer. For a better appraisal the same initial model was used in all cases. The MRLS algorithm is also demonstrated for a field data set: a stable solution is obtained.     

某煤矿典型CSAMT法视电阻率曲线的一维模拟

了解煤矿开采区的地层地质情况,指明断层、破碎带、陷落柱等构造及其含水性,并分析其对开采过程中可能造成的安全威胁是煤矿安全生产的重要保证．可控源音频大地电磁测深(CSAMT)法是浅层地球物理勘探的一种有力的手段,可以用来探明从地表到地下1500米范围的地质情况．在开采区完成三维地震勘探,对构造有了较系统的控制后,利用电法CSAMT手段重点来分析构造及地层的含水及富水特性,是一种科学的工作流程．本文采用一维正演方法对该矿区典型的CSAMT视电阻率曲线进行模拟并与实际曲线进行了对比．通过大量的模型拟合,得到了较好的拟合效果,这对反演解释提供了有用信息。     

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.     

时间域电磁法视电阻率的一种数值计算方法

主动源时间域电磁测深方法的主要优势在于其稳定的电磁响应信号和优越的勘探深度,特别是其对低阻薄层的敏感反映.但传统的视电阻率计算方法是采用极端情况进行定义的,因而对极端条件的合适度与电磁响应结果的数据处理和解释如何保持一致是一个很难确定的前提.本文建议采用一种直接利用观测结果进行2阶数值微分的方法来求解对应地点断面的视电阻率.这种方法计算所得视电阻率能很好的与地点断面的真实电阻率对应.这就为时间域电磁测深结果的直接解释或为反演初始模型参数选择提供了良好的基础.     

Rapid inversion of data from 2D resistivity surveys with electrode displacements

Resistivity monitoring surveys are used to detect temporal changes in the subsurface using repeated measurements over the same site. The positions of the electrodes are typically measured at the start of the survey program and possibly at occasional later times. In areas with unstable ground, such as landslide‐prone slopes, the positions of the electrodes can be displaced by ground movements. If this occurs at times when the positions of the electrodes are not directly measured, they have to be estimated. This can be done by interpolation or, as in recent developments, from the resistivity data using new inverse methods. The smoothness‐constrained least squares optimisation method can be modified to include the electrode positions as additional unknown parameters. The Jacobian matrices with the sensitivity of the apparent resistivity measurements to changes in the electrode positions are then required by the optimisation method. In this paper, a fast adjoint‐equation method is used to calculate the Jacobian matrices required by the least squares method to reduce the calculation time. In areas with large near‐surface resistivity contrasts, the inversion routine sometimes cannot accurately distinguish between electrode displacements and subsurface resistivity variations. To overcome this problem, the model for the initial time‐lapse dataset (with accurately known electrode positions) is used as the starting model for the inversion of the later‐time dataset. This greatly improves the accuracy of the estimated electrode positions compared to the use of a homogeneous half‐space starting model. In areas where the movement of the electrodes is expected to occur in a fixed direction, the method of transformations can be used to include this information as an additional constraint in the optimisation routine.     

电阻率二维神经网络反演

由于非线性特性地球物理反演一直以来都是一个比较困难的问题. 近十年来,非线性反演方法如人工神经网络、遗传算法在地球物理数据解释中得到越来越多的应用,但目前基本仍限于一维反演问题. 对于二维反问题,反演参数较多,神经网络反演运用较少. 本文利用BP神经网络优化方法,实现了电阻率二维非线性反演. 与传统线性化的迭代反演比较,神经网络反演能够克服传统方法的不足、获得更好的反演结果.     

多通道瞬变电磁法频率域比值算法及视电阻率计算研究

视电阻率在电磁法勘探中起重要作用,可以作为初步解释结果,多通道瞬变电磁法通过大地脉冲响应的峰值时刻定义视电阻率.鉴于传统方法计算峰值时刻视电阻率需要恢复完整大地脉冲响应,本文通过研究大地脉冲响应及其频谱特征,提出一种新的频率域比值法.利用参考谱及其与一般谱的伸缩因子计算一般谱的峰值时刻,进而得到视电阻率.该方法只需大地脉冲响应单个频点幅度谱,从而不用估计整个大地脉冲响应.实际计算时,可以用大地脉冲响应多个频点幅度谱得到平均视电阻率.模拟数据与实测数据结果表明,在一维介质中,频率域比值法可以得到合理的视电阻率.

     

One-dimensional magnetotelluric inversion using an adaptation of Zohdy's resistivity method

An iterative refinement method for determining a layered resistivity model from a Schlumberger or Wenner sounding curve is adapted to determine a layered resistivity model by using apparent resistivity and phase derived from the magnetotelluric impedance. Magnetotelluric observations presented as a function of period are first converted to an approximate resistivity–depth profile using Schmucker's transformation and this is used to construct an initial guess (starting) model. A two-stage procedure is then invoked. Keeping resistivities constant, layer boundaries are first adjusted to give a minimum misfit between measured data and responses and this is followed by resistivity adjustments with fixed layer boundaries to reduce the misfit further. The method is illustrated by application to some synthetic data both exact and with added noise, to a real field data set and to some magnetotelluric profile data obtained in a survey over the Carnmenellis granites in south Cornwall. The method is validated by recovering conductivity models from the exact and noisy 1D synthetic data. For complicated three-dimensional data at a single site and along a profile of stations, the method is shown to produce acceptable solutions which may be used as starting models in further two- or three-dimensional studies.     

地震地电阻率数据处理方法

本文总结了消除地电阻率现测资料中干扰的数据处理方法和地震地电阻率异常的判定方法,分析了它们的优缺点及适用条件,可为地电阻率观测资料的系统分析处理提供参考,有助于更好地利用地电阻率观测来为地震监测预报服务.     

Piecewise 1D laterally constrained inversion of resistivity data

In a sedimentary environment, layered models are often capable of representing the actual geology more accurately than smooth minimum structure models. Furthermore, interval thicknesses and resistivities are often the parameters to which non‐geophysicist experts can relate and base decisions on when using them in waste site remediation, groundwater modelling and physical planning. We present a laterally constrained inversion scheme for continuous resistivity data based on a layered earth model (1D). All 1D data sets and models are inverted as one system, producing layered sections with lateral smooth transitions. The models are regularized through laterally equal constraints that tie interface depths and resistivities of adjacent layers. Prior information, e.g. originating from electric logs, migrates through the lateral constraints to the adjacent models, making resolution of equivalences possible to some extent. Information from areas with well‐resolved parameters will migrate through the constraints in a similar way to help resolve the poorly constrained parameters. The estimated model is complemented by a full sensitivity analysis of the model parameters, supporting quantitative evaluation of the inversion result. Examples from synthetic 2D models show that the model recognition of a sublayered 2D wedge model is improved using the laterally constrained inversion approach when compared with a section of combined 1D models and when compared with a 2D minimum structure inversion. Case histories with data from two different continuous DC systems support the conclusions drawn from the synthetic example.     

1D inversion of DC resistivity data using a quality-based truncated SVD

Many DC resistivity inversion schemes use a combination of standard iterative least-squares and truncated singular value decomposition (SVD) to optimize the solution to the inverse problem. However, until quite recently, the truncation was done arbitrarily or by a trial-and-error procedure, due to the lack of workable guidance criteria for discarding small singular values. In this paper we present an inversion scheme which adopts a truncation criterion based on the optimization of the total model variance. This consists of two terms: (i) the term associated with the variance of statistically significant principal components, i.e. the standard model estimate variance, and (ii) the term associated with statistically insignificant principal components of the solution, i.e. the variance of the bias term. As an initial model for the start of iterations, we use a multilayered homogeneous half-space whose layer thicknesses increase logarithmically with depth to take into account the decrease of the resolution of the DC resistivity technique with depth. The present inversion scheme has been tested on synthetic and field data. The results of the tests show that the procedure works well and the convergence process is stable even in the most complicated cases. The fact that the truncation level in the SVD is determined intrinsically in the course of inversion proves to be a major advantage over other inversion schemes where it is set by the user.     

基于全区视电阻率的半航空瞬变电磁法对低阻薄层的探测能力探讨

针对使用感应电动势定义全区视电阻率会出现"二值性"的问题,本文采用与电阻率一一对应的垂直磁场分量结合二分搜索算法计算全区视电阻率,这使得全区视电阻率具有唯一性和可靠性.将全区视电阻率的概念引入到半航空瞬变电磁法中,在此基础上,建立以中间层为低阻薄层的三层地电模型,分别绘制了不同低阻薄层埋深与层厚比、围岩与低阻薄层电阻率比的全区视电阻率异常曲线图,并分析其全区视电阻率异常曲线形态及变化规律,通过计算其相对误差,分析并归纳了半航空瞬变电磁法对低阻薄层的探测能力.结果 表明全区视电阻率对地下低阻薄层分辨率较好,能较好反映地电断面信息,这对半航空瞬变电磁法有效探测低阻薄层具有实际的指导意义.     

基于交叉梯度约束的电阻率法和背景噪声法三维联合反演研究

电阻率法和背景噪声法是通过获得地下介质的电阻率和速度参数的分布来探究地球内部物质分布的非均匀性特征,联合反演可以有效地发挥两种方法的优势,减小反演多解性.本文采用有限内存拟牛顿（L-BFGS）算法实现了电阻率法和背景噪声法的单方法三维反演,然后基于电阻率和速度模型结构耦合约束,采用交替迭代方式实现了电阻率法和背景噪声法的三维联合反演.通过设计规则体组合模型和不规则体组合模型进行理论模型合成数据三维反演试算,结果表明：联合反演可以获得同时满足数据拟合差和结构更为相似的速度-电阻率模型;联合反演可以压制背景噪声单方法反演出现的假异常,获得更精细的速度结构;联合反演获得的电阻率模型对倾斜异常体、高阻覆层下方异常体、圈闭的高/低阻体等边界信息有明显的提升,有效克服电阻率法单方法反演的局限,提高深部电阻率的分辨率.