Incorporation of Topography into Two-dimensional de Resistivity and CSAMT Joint Inversion

One approach to increasing the reliability of geoelectric data interpre­ tation involves the application of multiple techniques and processing the data by using a joint inversion scheme. However, most of the current joint inversion techniques are based on one-dimensional models with a flat ground surface. Topographic effects are addressed with a correction factor. Therefore, good results can only be obtained for simple geologic structures. To avoid the influence of a possible incomplete topographic correction, the topography should be incorporated in the inversion model before being undertaken. Here, a 2-D joint inversion algorithm for direct current (de) resistivity and controlled source audio-frequency magnetotellurics (CSAMT) data in which finite source effects are not considered is developed. This algorithm is able to invert sounding data collected on an irregular earth surface with­ out presuming a flat earth surface. De resistivity and CSAMT modeling are based on the finite-element method; and, the iterative joint inversion scheme is derived from the second-order Marquardt damped least-squares method. The algorithm has been tested on both synthetic and field de resis­ tivity and CSAMT data with explicit incorporation of topography in the joint inversion. Both synthetic and field studies indicate that the technique is computationally efficient in providing an improved geologic interpreta­ tion for complex subsurface structures. (


INTRODUCTION
Electrical resistivity methods are used widely to study environmental and engineering problems.Often, one-dimensional (1-D) horizontal-layered resistivity models can not fit the observed data satisfactorily, because of strong lateral inhomogeneity.In order to obtain detailed, reliable models whose response fits close to the data, interpretations of two-dimensional (2-D) de resistivity data by the inversion techniques are investigated.Pelton et al. (1978) developed an inexpensive computer algorithm for estimating the location and resistivity of a single 2-D prism with a rectangular cross-section.Oristaglio and Worthington (1980) and Tripp et al. (1984) also developed 2-D inversion algorithms for the line-source electromagnetic method and the de resistivity sounding method, respectively.Since the resistivity surveys sometimes were carried out in mountainous areas, where topographic effects may produce spurious anoma lies in the resistivity section, Fox et al. (1980) showed a terrain-correction procedure to reduce topographic effect in resistivity data.Holcombe and Jiracek (1984) further suggested that the corrected data can be interpreted as if they were the response of a flat-earth model.Tong and Yang (1990) developed an algorithm based on the Tripp et al. (1984) approach with a refine ment that the topography is incorporated in the inversion model.
The most common models used in geoelectric interpretation assume a horizontally lay ered earth, with each layer having a constant resistivity.In geoelectric surveys, the depth of penetration is limited by the surface topography and/or available area for the array configuration.Different geoelectrical methods have their own optimum ranges of penetration depth.Therefore, in the sampling of information at different depths by different data sets, the need for joint inversion is quite obvious.Joint survey and inversion schemes have to be developed to esti mate the subsurface resistivity distribution.For example, Vozoff and Jupp (1975) developed an algorithmfor jointly inverting magnetotelluric (MT) and Schlumberger sounding data, Gomez-Trevino and Edwards (1983) for controlled source EM and Schlumberger soundings, and Raiche et al.(1985) for coincident loop time domain electromagnetics (TEM) and Schlumberger soundings.Yang and Tong (1988) developed an algorithm for inverting the de, TEM, and MT responses simultaneously.Although these joint inversion methods are based on 1-D models, they did improve the resolution of the interpretive layered earth models.
In order to extend and connect the works described above, we have developed an algo rithm of 2-D joint inversion of de and control source audio-frequency magnetotelluric (CSAMT) sounding data to get a better estimate of the subsurface resistivity distribution.Joint de and CSAMT inversion was done by Sasaki (1989), in which both synthetic and field data did not involve the effects of surface topography.In addition, we used TE mode in stead of TM mode as shown in Sasaki's paper.Tong and Yang (1990) made a comparison between a conven tional topography-corrected resistivity pseudosection of the results of a fault model for dipole dipole and resistivity pseudosection obtained by resistivity inversion with topography, which reveals substantial differences.In addition, variation of the receiver relative to the transmitter will produce elevation errors in electric or magnetic fields along a traverse line, relative to the fields that would be observed over a flat surface (Tripp et al. 1978).Therefore, we used a finite-element method for both the de and CSAMT modeling, explicitly incorporating an ir regular earth surface instead of assuming a flat earth in the inversion scheme.The proposed inversion technique was applied to both synthetic and field data.Our results, as described below, demonstrate that the algorithm can perform efficiently, and produce reliable estimates of subsurface resistivity.

2-D INTERPRETATION
Two-dimensional (2-D) interpretation of the apparent resistivity pseudosections is highly desirable for a full understanding of data collected with 2-D coverage.The data sometimes were measured over a rough terrain surface.The joint inversion technique used here is quite similar to Tong and Yang's scheme (1990).The errors caused by assuming a flat ground surface are removed by incorporating topography in the inversion.Here is a brief description of this technique.
The first step is to input the de apparent resistivity data and the CSAMT apparent resistiv ity and phase data at nodes on the mesh of an element.The initial model is based on the inversion model of 1-D de and CSAMT data for each location.After this the forward model responses are calculated.The model responses are then compared with both de and CSAMT field data, and the model resistivities are automatically adjusted by using an iterative second order Marquardt damped least-squares scheme (Tong and Yang 1990).The forward calcula tion is repeated.After each iteration the r.m.s.value of the differences between field data and model response is calculated.The process is repeated and the model is adjusted until an ac ceptable fit is achieved.

FORWARD MODELLING
The forward program is a modified version of the finite element code for modeling 2-D structures incorporating topography, presented previously by Tong and Yang (1990) except CSAMT responses were also calculated.For de cases, this technique involves solving a three dimensional potential distribution due to point sources located in a two-dimensional half space, often referred to as a 2.5-D modeling.For the CSAMT case, the data were previously made of near field corrections, so they could be treated as purely plane-wave EM data.The modeled subsurface is divided into discrete elements, in which the element sizes can vary, but their shapes are always rectangular.The element can be coarsened successively towards the outer edges in order to simulate infinite extent.Surrounding edge elements are assigned with values of the nearest pseudosection data points.The resistivity distribution is defined for each element at the nodes of the mesh, and therefore can be varied arbitrarily between nodes.The mesh was designed as a trade-off between lateral coverage, sufficient depth penetration and time-memory requirements for the forward calculation.Topographic effects are included in the model.The geometry of the subsurface determined by the resistivity values is estimated in the inversion.Details of the modeling approach can be referred in Tong and Yang (1990).

INVERSION METHOD
The application of the inversion technique to geoelectrical methods has been described before (e.g., Inman 1975;Vozoff and Jupp 1975).The Jupp-Vozoff algorithm, which is an iterative second-order Marquardt least-squares scheme, is used here .. We present only the ba sic equations as they apply to the de resistivity problem.
Assume N unknown model parameters x are related to the M resistivity data d by a nonlin ear functional h defined as where c denotes the electrode separation.This nonlinear equation is linearized by expansion in a Taylor's series at each electrode separation.Elimination of the higher order terms leads to a linear set of M equations in N unknowns, where and ' ' (2) The vector L1d is the difference between the measured apparent resistivity or phase data and the theoretical apparent resistivity or phase of the initial model.The vector ,1x 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 f!, sometimes known as the Jacobian matrix, is the matrix of partial derivatives of the apparent resistivities with re spect 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 = ln d, H = ln h, and X = ln x.Using these equations and adding a data weighting matrix W into the Jacobian matrix and equation ( 2), and Detailed examples of partial derivatives of de and MT responses can be cited from Tong and Yang (1990) and Sasaki (1989).The inverse procedure is to find the model parameters which minimize !ID .Since equation ( 4) was deduced by linearizing a nonlinear system, the solution requires several iterations.At each step, equation ( 4) is solved for LiX and yields a new set of model parameters.This procedure is repeated until an accepted minimum-square residual is reached.Gloub and Reinsch (1970) decomposed the Jacobian matrix A into its row and column eigenvectors, and the associated singular values, as (5) where U is an M by M data eigenvector matrix, l:'. is an N by N solution eigenvector matrix, and J is an N by N diagonal singular value matrix.
The parameter improvement vector LiX is obtained by substituting A from equation ( 5) into equation ( 4).The solution is (6 The problem is that when small singular values are present, the estimate for LiX is grossly contaminated by numerical noise.To overcome this problem, a damped N by N diagonal ma trix 'Lis added to equation ( 6); thus,

The elements of T are
LiX =VT A-I U T LID.
---t j = ky11 f(ky11 + µ 2 11 ).for j = 1,2, ...... N where µ is known as the relative singular value threshold and K j = A j I Ai.A j is the j th singular value and Ai is the maximum of the singular values.For T/ = 1, equation ( 8) repre sents the conventional Marquardt's technique, and for larger values of' T/' the singular value truncation method with relative threshold µ.We adopted a value T/ = 2 for this study.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 by initially giving µ a high value of 0.2 and then as the fitting errors decrease, per mitting µ to be decreased from iteration to iteration until it reaches a minimum allowed value.More details about this were cited by Tong and Yang (1990).
In the above description, there is no implicit restriction on the choice of inversion model parameters.In our application, the resistivity parameter is the resistivity of each block.The geometry of each block is specified and the topography is included explicitly.We estimated the resistivity of each block in the inversion process.The result was then the resistivity distri bution under the specified irregular surface.

NUMERICAL TESTING AND VERIFICATION
Numerical tests on both noise-free and noise-contaminated 2-D synthetic de resistivity data and CSAMT data calculated along an irregular surface were performed to evaluate the capabilities and limitations of our inversion algorithm.The synthetic data which were used are the computed apparent resistivities for a Schlumberger array and scalar CSAMT TE mode on a traverse perpendicular to the strike of a 2-D earth model.
A 90-degree dipping fault model incorporating a step-like ground surface (Fig. 1), was used to demonstrate our joint inversion approach.The vertkal steep surface coincides with the fault trace.Based on this fault model, the noise-free electric responses for both de resistiv ity and CSAMT data were computed and are displayed in Fig. 2. Figure 2a shows the modeled apparent resistivity pseudosection of de-resistivity, and Figs.2b, c show the modeled apparent resistivity and the impedance phase pseudosection for the CSAMT data respectively.The fault affecting the apparent resistivity distribution near the center of the section can be recognized from both figures.The block resistivities of the inverted model, whose responses fit the syn thetic noise-free data of the model to an acceptable degree, is shown in Fig. 3.This inversion routine requires three to ten iterations to converge to the final model.Although the resistivity of each block is slightly different from the true resistivity, one can easily identify the subsur face structure from the obvious resistivity contrasts.
Furthermore, in order to test the effect of noise-contamination on the inversion, 10 per cent random Gaussian noise was added to the synthetic data obtained from the previous fault model shown in Fig. 2. Figure 4a shows the modeled apparent resistivity pseudosection of     , c show the modeled apparent resistivity and also the impedance phase pseudosection for the CSAMT data.Again, the faults affecting the apparent resistivity distributions near the center of the section are still recognizable from both sections.The sub surface structure delineated by this inversion is obvious.The final model inverted from these data is shown in Fig. 5.It is very close to the synthetic fault model.The numerical modeling studies show that our joint inversion algorithm is able to resolve complex 2-D resistivity structures.

INVERSION OF FIELD DATA
In order to demonstrate the applicability of our inversion algorithm to interpret actual 2-D resistivity data measured along a rough topographic surface, we applied the method to the field data acquired at a site located at a gravel-covered area off the southern part of National Central University (NCU) campus, NW Taiwan (Fig. 6).The topographic feature of this site is quite similar to the synthetic fault model.The lithology of the subsurface at a depth of less than 150 m can be inferred from seven water wells located on campus.The top layer, with a thick ness ranging from 3 to 5 m, is laterite .The second layer consists of gravel, with a thickness ranging from 5 to 15 m.The bottom layer consists of alternating sand and clay layers, with thicknesses being greater than 130 m.There is a deep dried gas well located in the Pin-Jing area at about 1 km.from the test site.The lithologic units obtained from the well can be grouped into three formations.The Toukou-shan Formation (mainly laterite, gravel, sand and clay) with a thickness of 671 mis at top.At a depth from 671 to 1342 mis the Choulan Formation.The presence of shale in this formation is reflected by its lower resistivity than the overburden shown in E-log.Below the depth of 1342 mis Chinshui Formation.The dominant presence of shale in this formation makes its resistivity even lower than that of the Choulan Formation.Eight dipole-dipole sounding locations were selected with a maximum spread of 400 m.The survey line is marked by a straight line with a star as shown in Fig. 6.The inter-spacing of each sounding location is 50 m.The topography has a 10 m elevation change at the x-coordi nate of 25 m.In order to fit the 2-D assumption, the electrodes of each station along this traverse were spread perpendicular to the suspected fault strike.Figure 7a shows the observed  8b) reflects the presence of thick clay layers of the Toukou shan formation, and the low resistivity zone below the depth of 700 m can be attributed to the effect of Choulan Formation.Since there are no obvious lateral discontinuities of resistivity layers within a depth of 100 m, the view that the step-like ground surface was the end product of faulting is not compelling.

CONCLUSIONS
We have presented an algorithm for automatic 2-D joint inversion of de resistivity data and CSAMT data measured on an irregular surface.The main features of this scheme is to incorporate topography in the inversion model, so the measured apparent resistivity data can be inverted directly to determine the subsurface structure.The results of both numerical and field studies show that this inversion algorithm can provide a rapid and efficient analysis of the resistivity data.Our algorithm can deal with any four-electrode configuration for the de data and can be extended to study 2-D joint inversion of other geoelectric data.Some limitations of this 2-D joint inversion technique are as follows: (1) The geometry of the subsurface must be specified in advance.(2) The total number of inversion parameters allowed may be limited.(3) The process of dividing the subsurface into blocks must be done carefully, because improper block sizes will lead to poorer resolution of the resistivity and increase the number of iterations.(4) The surface elevation and the distance between elec trodes must be measured as accurately as possible; otherwise, the erroneous elevation and distance will lead to biased interpretations.The traditional topographic correction procedures have the same problem.(5) It is not possible to evaluate 3-D topography and subsurface struc tures from the 2-D joint inversion results.However, the use of this algorithm process is practical, rapid and reliable for delineating subsurface structures when subjected to the possible limita tions mentioned above.This approach can be used to analyze 2-D resistivity from joint geoelectric surveys carried out on either flat or irregular topographic surfaces.

Fig
Fig. 1.A vertical fault model incorporating a step-like ground surface.

Yang
et al.

Fig. 2 .
Fig. 2. Synthetic noise-free pseudosections for (a) de apparent resistivity (b) CSAMT apparent resistivity and (c) impedance phase, Computed data are based on the model shown in Fig. 1.