A FINITE DIFFERENCE ALGORITHM FOR TWO-DIMENSIONAL INVERSION OF GEOELECTROMAGNETIC DATA

Abstract

Data inversion is an interplay of physico-mathematical operators devised to extract meaningful information about a system from an observed data set and to appraise the quality of an inverse solution. In geoelectromagnetic methods, where sources are natural electromagnetic (EM) fields, the earth is parameterized in terms of electrical resistivity which is of special significance as itcarries information about the lithology, porefluid, temperatureand chemical variations. The presentstudyis an effort to improve the data inversion capabilities of EM data. Forthis purpose an efficient 2-D inversion algorithm, EM2INV, for geoelectromagnetic data has been developed. The EM field is a non-linear function of subsurface resistivity distribution. As a result the inverse problem is quasi-linearized and solved iteratively. For each inversion iteration, a new forward problem, yielding the response of current resistivity model, has to be solved. Therefore, the forward algorithm is a prerequisite for an inversion algorithm. For generation of EM response, a boundary value problem is solved analytically or numerically. However, for the problems involving complex geometries one has to seek numerical solutions only. Due to its simple mathematics and easy implementation, finite difference method has been chosen over other numerical techniques for solving the EM boundary value problem. The research work was initiated with the implementation of finite difference formulation of the forward EMproblem (Brewitt-Taylor&Weaver, 1976). Since the use of Dirichlet boundary conditions results in a large study domain, special finite domain, integral and asymptotic boundary conditions are implemented. In the present work, an integrated formulation of these boundary conditions has been developed. The quasi-linearization of non-linear problem results in a matrix equation which is solved using Bi-Conjugate Gradient Method (BCGM), a semi-iterative matrix solver that dispenses with the necessity of explicit computation of Jacobian matrix. To fix the in number of unknown block resistivities for all frequencies and throughout the inversion process, a superblock notion has been developed. The initial guess is made on the basis of observed anomaly and other a priori information. The inversion algorithm EM2INV is the culmination of research that started with the development of a primitive algorithm. The algorithm has been written in FORTRAN 77 and implemented on an IBM compatible EISA based PC-486 machine with 32 MB RAM and 383 MB hard disk, using the SVR 4.0 version of Unix operating system and F78 FORTRAN compiler. For a typical model, having 31 x 15 nodes, the algorithm takes about 3 minutes for 10 inversion iterations. EM2INV comprises 6120 lines, 42 subroutines and 3 function subprograms. The main program has two basic modules - Forward and Inverse. Its special efficiency features which result in cost effectiveness are - (i) Optimal grid generation based on grid design thumb rules, (ii) Finitedomain boundary conditions, (iii) Interpolation matrix that permits generation of response at observation points, (iv) Gaussian elimination, the forward matrix solver, which enables reuse of already decomposed coefficient matrix, (v) Use of logarithm of resistivity to ensure positive values of estimated parameters, (vi) Superblock notion that reduces the number of blocks with unknown resistivities and, in turn, the size of Jacobian matrix and (vii) BCGM matrix solver for inverse problem. Besides beingefficient, EM2INV is versatileon account ofitsfeatures like - (i) Inversion with field/synthetic data, (ii) Error free/erroneous synthetic data, (iii) Inversion of MT/GDS data and (iv) Inversion of profiling/sounding data. The algorithm has been rigorously tested by setting up exercises of diverse nature and practical significance. Forestablishing the validity offorward computations, the published results of various models have been reproduced after carrying out the no contrast and mesh convergence studies. Similarly, for checking validity of inversion computations, the synthetic anomalies have been inverted and compared with those of the true model. The stability of the algorithm has been established by inverting the synthetic response corrupted with Gaussian noise. EM2INV has been employed in two different sets of experiments designed to study the nature of forward and inverse problems. The forward experiments aim at studying the impact of parameters like depth of burial, resistivity contrast, separation between two bodies on model responses. Although a preliminary quantitative discriminant analysis was attempted to design thumb rules for estimating size and iv resistivity of target yet it did not succeed. However, qualitative inferences have been drawn. The inversion experiments performed are aimed at gauging - (i) Relative performances of response functions, (ii) Inversion quality fo two modes of polarization, (iii) Efficacy of single and multifrequency inversions and (iv) Minimum number of frequencies and observation points needed for successful data inversion. The inversion of MT data provides better estimates of vertical position of the target, whereas the inversion of GDS data deciphers the horizontal variations better. It has been observed that the conductive and resistive bodies are better resolved by inversion of E- and Bpolarization data respectively. The results of multifrequency inversion imply that increase in number of frequencies does not necessarily enhance the inversion quality es ecially when the spread of observation points is sufficiently large to sense the target. The study of minimum number of observation points highlights the importance of single point inversion which furnishes useful information about the inhomogeneity. After the design exercises, EM2INV has been exhaustively tested by inverting synthetic data, field data, as well as data derived from models based on field studies. Few geologically significant models are picked up from the literature for generation of synthetic data. For these models, initially 1-D inversion is carried out at each point of the profile which is then stacked to get the starting model for 2-D inversion. The comparison of inverted model with the 1-D stacked model leads us to conclude that 2-D inversion substantially improves the quality of the inverted model. Next, a study has been carried out on models derived from GDS or MT field studies. The reliability of the estimates of resistivity is evident in the goodness of fit of the computed and observed responses. Lastly, the algorithm EM2INV is tested on Trans Himalayan conductor and COPROD2, GDS and MTfield data respectively. The inverted models are in broad agreement with the published results. This supported the confidence in the utility of the algorithm. The results of various experiments and those of inversion of synthetic/field geoelectromagnetic data in terms of resistivity model have established the veracity of the algorithm and also amply displayed the capabilities of the inversion algorithm. Also discussed, is the possible scope of future work in various directions for its upgradation and extension to 3-D environment.