A Study of Joint Inversion of Direct Current Resistivity, Transient Electromagnetic and Magnetotelluric Sounding Data

Geoelectric methods have been widely used in the field of engineering and mineral exploration. Each method has its own advantage and disad­ vantage in interpreting depth and electrical properties of underlying sub­ surface. Combination of sounding data from various methods has been proved to be an effective way to improve reliability of interpretation as cited from many published literatures. Among them, several authors de­ veloped intriguing schemes to simultaneously invert two sounding data col­ lected from two different surveying methods into a geoelectric section. These studies show that the interpreted sections based on joinMnversion scheme are superior to those from inversion of single data set alone. This research is an extension of their work by inverting coincident loop transient electromagnetic (TEM) data, magnetotelluric (MT) data and Schlumberger sounding (DC) data simultaneously. An appropriate joint inversion scheme has been developed. Synthetic sounding data sets for DC, MT, and TEM based on a six-layer earth model were used to evaluate the performance of our scheme. Our study shows that the estimated model parameters from joint inversion of the triple data sets are more consistent with the hypothetical model than those derived by any inversion result us­ ing DC, TEM, or MT data alone. Furthermore, our result also indicates that our proposed scheme is superior to the joint inversion scheme based on DC and MT, DC and TEM, or TEM and MT.


INTRODUCTION
The most common model used in the geoelectric interpretation is the horizontally layered earth. In a geoelectric survey, the depth of investigation is always limited by the topographic surface and/or available area for spread arrangement. Therefore, the direct current (DC) resis tivity and transient electromagnetic (TEM) method are suitable for investigating the subsur-face at a depth less than 1000 m. The magnetotelluric (MT) method reveals information from deeper parts of the earth than DC or TEM method does. In addition, resolution of the layered structure in geoelectric interpretation as in general is affected by the electrical properties of the geoelectric layers. In order to overcome these shortcomings from each method, joint inver sion schemes had been developed to estimate subsurface resistivity distribution. For example, Vozoff and Jupp (1975) developed an algorithm for joint inverting MT and Schlumberger sounding data. Joint inversion is also suggested by Gomez-Trevino and Edwards (1983) for controlled source EM and Schlumberger soundings, and Raiche et al. (1985) for coincident loop TEM and Schlumberger soundings. All these investigation show improvement of resolu tion of the interpreted layered earth structure.
To extend the work described above, we developed an algorithm of joint inversion of triple data sets of coincident loop TEM data, MT data, and Schlumberger DC sounding data.
These data are inverted simultaneously in order to provide a better estimation of subsurface resistivity distribution. The algorithm has been programmed and tested by synthetic data generated from a horizontal six-layer earth model. Our results demonstrate the advantage of the triple joint inversion. Comparison made a:mong the results from triple joint inversion, single inversion, and double inversion techniques indicates that our proposed triple inversion scheme is better in describing model parameters than the other inversion schemes. The algo rithm developed by this study represents an efficient and reliable inversion approach.

INVERSION METHOD
The application of inversion techniques to geoelectric methods has been described before (e.g., Gomez-Trevino and Edwards, 1983;Inman, 1975;Raiche et al., 1985;Vozoff and Jupp, 197 5). The Jupp-Vozoff algorithm, which is an iterative second-order Marquardt least-squares scheme, is used here. We present only the basic equations as they apply to the geoelectric problem.
The vector 11d is the difference between the measured apparent resistivity, phase, or transient voltage data and the theoretical appa!¥nt resistivity, phase, or transient voltage of the initial model, respectively. The vector /1 x is the difference between the unknown model parameters and the initial model parameters. In our case, the model parameter is the resistivity of each block. The matrix g, referred to as the Jacobian matrix, is the matrix of partial deriva tives of the data values with respect to the model parameters.
Because both the model parameters and the resistivity data values may vary over several orders of magnitude, and also because they must be constrained to positive values, we use logarithmic fitting. Thus, let D =Ind, H =In h, and X = ln x. Using these equations and incorporating a data weighting matrix W into the Jacobian matrix and equation (2), and (4) The inverse procedure is to find the model parameters which minimize AD . Since 4For a joint inversion of different sou11Wng data, as shown in equation (5), the error vector where the elements in each subvector are the differences between the measured physical pa rameters and the theoretical physical parameters for a specific sounding method. The se-quence of superscript shown in equati�n (5) represents the number of survey methods being used. Computed 4 \he model update!J.X by equation (5), we calculate each element of the error vector, !J. D , and Jacobian matrix M 1 for different sounding surveys separately, then substitute the final results again into equation (5). Model parameters can be achieved by several iterations similar to conventional single inversion schemes.
Gloub and Reinsch (1970) decomposed the Jacobian matrix d into its row and column eigenvectors, and the associated singular values, as (6) where U is an M by N data eigenvector matrix, V is an N by N solution eigenvector matrix, and A is an N by N diagonal singular value matrix. -4 The parameter improvement vector !J. X is obtained by substituting d from equation (6) into equation (5). The solution is The problem is that when small singular values are present, the estimate for !J. X is grossly contaminated by numerical noise. To overcome this problem, a damped N by N diagonal matrix I is added to equation (7); thus, The elements of I are (9) where µ is known as the relative singular value threshold and kj = A.j I A, . A.j is the j th singular value and A, is the maximum of the singular values. The estimate is further stabilized by initially including only the largest singular values in the estimate. As the fitting error decreases, the singular values of less importance are also included in the estimate. This is performed by initially giving µ a high value of 0.2 and then, as the fitting error decreases, permitting µ to be decreased from iteration to iteration until it reaches a minimum allowed value.

INITIAL DATA PREPARATION
To process efficiently Schlumberger sounding curves or Wenner field curves, Zohdy (1989) proposed a fast iterative automatic interpretation method. The algorithm is based on inter-preted depths and resistivities obtained from shifted electrode spacings and adjusted apparent resistivities, respectively. The advantage of this algorithm is that it does not require an initial guess of the number of layers, their thicknesses, or their resistivities. In this paper, the itera tive procedure to interpret DC field sounding curves is based on Zohdy' s algorithm. Further more, a modified Zohdy algorithm suggested by the authors can be easily implemented to iterate the MT data, Le., using sounding frequency shifting instead of shifting electrode spac ings to obtain an initial guess of the number of layers, their thicknesses, and their resistivities for MT data inversion.

FORWARD MODELING
A synthetic test on noisy artificial sounding data, in which we add 5 % random Gaussian noise to the DC sounding and 10 % to the TEM and MT soundings, had been undertaken to evaluate our joint inversion algorithm. Figurel shows the six-layer test model. The synthetic data used for inversion are the computed apparent resistivities (for DC and MT data) and transient voltages (for TEM data). Figure2 shows the Schlumberger sounding curve. Figure3 shows the computed transient voltage curve and the corresponding apparent resistivity curve of the TEM response for a coincident l(lop configuration, The method of Knight and Raiche ( 1982) was used to calculate the TEM response. The computed MT apparent resistivity sound ing curve is shown in Figure 4. From these synthytic " observed " sounding data, it is obvious that electrical properties of the first to the third layer are more influential on the DC data, the third to the fifth layer on the TEM data, and the fourth to the sixth layer on the MT data.

S. INVERSION OF ARTIFICIAL DATA
Joint inversion of the DC, TEM, and MT responses by using the  inversion scheme was carried out, a second-order Marquardt method was used to stabilize the        inversion process. The results are shown in Table 1. The results from inverting single DC; TEM; or MT; joint DC-MT; and joint DC-TEM data using the same starting model are also shown in Table 1. Table I indicates that the model parameters derived from the joint inversion of DC-TEM MT data are more consistent with the original model than those from any single or double data inversion result. Figure 5 shows the eigenvalues and the corresponding eigenvectors, eigen values of model and data space. The eigenvectors represent a measure of the inversion perfor mance. The well resolved shallow layer parameters are recognized by the solution eigenvec tors of the eigenvalues 3.545, 2.992, and 1.090. At the eigenvalues of 6.334, 5.213, and 1.720, it shows a good resolution of the medium deep layer parameters by the TEM sounding data. At the eigenvalues 3.323 and 2.021, the well resolved deeper layer parameters by the MT sounding data can be recognized.

CONCLUSIONS
In this paper an algorithm for the joint inversion of DC, TEM, and MT sounding data has been presented to solve some basic problems inherent in each geoelectric data interpretation. The shortcoming of each method may partially be compensated by the others. As pointed out by Raiche et al. (1985), resolution of layer resistivity and thickness may be achieved by using an eigenvalue analysis. Based on our model study, we have shown that the joint inversion of triple data is proved to be more effective in parameter estimation than inversion of single or inversion of double data sets.