Lithological Characterization and Its Application Based on Three-Dimensional Structure-Coupled Joint Inversion of Gravity and Magnetic Data

Incorporating structural-coupling constraint, known as the cross-gradients criterion, helps to improve the focussing trend in cross-plot of multiple physical properties. Based on this feature, a post-processing technique is studied to characterize the lithological types of subsurface geological materials after joint inversion. A simple domain transform, which converts two kinds of participant physical properties into an artificial complex array, is adopted to extract anomalies manually from homogenous host rock. A synthetic example shows that structure-coupled joint inverted results tend to concentrate on the feature trends in the cross-plot, and the main geological targets are recovered well by a radius-azimuth plot. In a field data example, the lithological characterization reveals that the main rock types interpreted in the study area agree with the geological information, thus demonstrating the feasibility of this technique.

measurements.Based on certain established relationships between lithological types and their statistical property ranges, cheap but effective lithological mapping in a 3-D space can be conducted using the geophysical inversion results [1].Kanasewich and Agarwal [2] presented a method of computing the ratio of density to magnetic susceptibility in the wavenumber domain, which was then applied to real data and demonstrated the practical use of this method.Based on the correlation between gravity and magnetic anomaly gradients, Price and Dransfield [3] conducted advanced studies on apparent lithology mapping.Lane and Guillen [4] investigated the lithological classification technique assisted by inverted physical property models such as density and magnetic susceptibility.Williams and Dipple [5] proposed a method of mapping alteration minerals with added borehole data.Kowalczyk et al. [6] inverted gravity and magnetic data and found correlations between lithotypes and physical properties, which were then used to identify the geological units in a 3-D subsurface space.
However, direct transform from ordinary inverted physical properties to lithological types may be incorrect, because the incompatible resultant physical property models pose difficulties in identifying lithotypes [7].To address this problem, the joint inversion framework proposed here for yielding more compatible, high-resolution models simplifies the recognition of lithological types from inverted property models [8] [9].In the present study, we first briefly review the cross-gradients joint inversion scheme of gravity and magnetic data, which lays a solid foundation for the subsequent lithological characterization.Then, we introduce a domain transform method which, unlike the commonly used topological analysis technique, is based on multi-physical cross-plot.The differences among various lithotypes are amplified by this operation.A synthetic example is provided to illustrate the use of this technique.Finally, it is applied to real data, and its practical use and effectiveness are demonstrated.

Structural-Coupled Joint Inversion
As described by Gallardo and Meju [8] [9], the cross-gradients vector is defined as the cross-product of two spatial gradient vectors of two physical property models (e.g., density contrast and magnetization in our case): where m represents the sorted discrete model parameters; x, y are the horizontal coordinates; and z is the vertical coordinate.The subscript number of the vector marks the type of model concerned in the process.For 3-D models, this function has three components, each of which quantifies the structural similarities in planes normal to the component direction.The resemblance is obtained once all of the components approach zero.Combined with the cross-gradients constraint, a reconstructed model is obtained with the minimum value of the corresponding operator L-2 norm.To obtain stable iterations, the cross gradients are incorporated as part of the objective function: Here, ϕ represents the objective function of the data-fitting term and extra model regularization term [10].The cross-gradients criterion requires the minimization problem to satisfy the condition = 0 τ , which means that any spatial changes occurring in both models must lie in the same or in opposite direction [8] [9] [11].This implies that if a common boundary exists, their gradients must be sensed in a parallel orientation regardless of the amplitude of the model parameter changes.
Following the robust statistical optimization algorithm framework developed by Fregoso and Gallardo [11], the inversion problem of minimizing the objective function subject to the given conditions is solved.The data and the discrete models are submitted into a general framework and are processed simultaneously.To construct the normal equations, the physical relationships and the smoothing mechanism should be involved.In addition, the structural constraints provide the correlation between the model parameters that formulate the total algorithm.When the models are coupled in structure, the cross-gradients operator tends towards null space.Since the subsurface is discretized in sequence into a set of prisms, the cross gradients may be approximately computed in a forward difference form.Using linearized Taylor expansion, the corresponding partial derivatives can be readily obtained using the disassembling form of the cross-gradients vectors.

Lithological Characterization Technique
After cross-gradients joint inversion is carried out, the resultant density and magnetization distribution are con-verted into lithology units using the specified statistical information.The inverted models should be examined carefully before converting to eliminate unexpected misleading features due to error propagation.A synthetic example is employed here to explain the whole process.
As Figure 1 shows, a model consists of four geological targets (seen as rectangular boxes) embedded in the half-space.Density contrast and magnetization properties are assigned to the model (listed in Table 1).The anomaly is designed with different property values so that it can be recognized by petrophysical means.The ambient field has an inclination I = 45˚ and a declination D = 45˚.The observed response is modeled using gravity and magnetic forward programs with 2.5% Gaussian noise.
The recovered slices of density contrast and magnetization distribution horizontally located at depth = 200 m are shown in Figure 2. It is obvious that the separately imaged anomaly (top) reflects the basic feature of the rectangular boxes, but the resolution is low.Note that the interval between two adjacent boxes is blurred due to the smoothness hypothesis.A significant improvement in the jointly imaged density model is to distinguish anomalies with the same attributes.The boundary between two boxes with opposite attributes is also sharpened, thanks to the cooperation between the interactive structural constraints.The magnetization image has a similar appearance to the density image, and the image amplitude is enhanced when compared to the separate results, showing that the joint images have an enhanced focussing feature.
The cross-plot of the two properties is a practical tool for demonstrating congruent relationships in the spatial scale (shown in Figure 2(c) and Figure 2(f)).Because of the smoothness constraints of the inversion scheme, the results tend to produce broad tails around the anomalies, which are displayed in the form of fuzzy cluster clouds (shown in Figure 2(c)).Fortunately, the divergence phenomenon is suppressed such that it is compressed close to the feature trends.This characteristic shows the potential for extracting additional information from the results and thus produces images that are both more detailed and more reliable.The joint inversion results show that lithological characterization can readily be conducted for any relationship that is more obvious in the cross-plot.In addition, a domain transform is used to amplify the difference between geological targets, and also the contrast with the homogenous background.That is: , , , , where i is a complex unit; c is an artificial complex number composed of real part 1 m and imaginary part 2 m ; r is the modulus of c; and θ is the azimuth angle of c on the complex plane.To sum up, the convergence of the relationship between two properties makes it possible to classify different types of subsurface medium based on the distinguishing feature line.The domain transform of the physical properties amplifies the differences between the geological targets, and make it easier to differentiate between anomalous bodies and host rock.

Application to Field Data
The actual collocated gravity and magnetic survey data was obtained from a mineral exploration venture in South China that contained 5625 gravity and 4362 magnetic acquisition points approximately located within a study area measuring 7 × 8 km 2 (Figure 4).The data was appropriately gridded to carry out 2-D low pass filtering in order to eliminate the negative influence of high-frequency noise.The study area is located in the effusive rock basin formed during faulting and magmatic activities in the late Yanshan.This district is mainly dominated by andesitic volcanic rock, volcanoclastics, tuff, syenite, etc.The orebody deposit occurs mainly within the tectonically crushed zone and fractured andesite.
Figure 5 shows the cross-plot of joint inverted density and magnetization distribution.By incorporating structural constraints, the results are consistent, as further shown by the distinct trends on the cross-plot.This reveals the potential for conducting advanced analysis using the proposed lithological characteristics technique.The radius-azimuth plot, which is converted from the physical properties plane of Figure 5, is shown in Figure 6.The red rectangular boxes were marked manually in this presentation.In this case three types of rocks were generally

Conclusion
A simple, practical and effective lithological characteristic technique is presented in this paper based on structure-coupled joint inverted physical property models.The good convergence of the structure-coupled models, which is shown as agglomerated feature lines rather than fuzzy cloud clusters, demonstrates the possibility of classification using this method.Synthetic and real data examples show the potential of the technique for correctly identifying different lithological types in subsurface zones and establishing a solid basis for further exploration decisions.

Figure 1 .
Figure 1.Schematic of four geological units embedded in a uniform host medium.

Figure 2 .
Figure 2. Images of density and magnetization distribution at horizontal slicing plane for depth = 200 m and cross-plots of property in spatial scale: (a) separate density image; (b) separate magnetization image; (c) density-magnetization cross-plot for the separate case; (d) joint density image; (e) joint magnetization image; (f) density-magnetization cross-plot for the joint case.

Figure 3 Figure 3 .
shows the radius- azimuth angle from which the anomaly can be identified by the vertical trends at different angles.The four geological targets are recognizable from the bar-graph-like distribution.The low values homogeneously located close to the zero line are regarded as the host rock contaminated by noise.The spatial targets are recovered which is shown in All anomalous types have been divided successfully and their locations match the originals reasonably well when compared with Figure1.

Figure 4 .
Figure 4. 3-D sketch map of disctreted subsurface materials and its observed gravity and magnetic field data.

Figure 5 .
Figure 5. Cross plot of density and magnetization after joint inversion.

Figure 6 .
Figure 6.Anomaly filtering by radius-azimuth.The boxes indicate identified rock types.extracted from the whole subsurface, leading to the recovered lithology model presented in Figure 7. Rock 1 reflects the physical trends near azimuth angle 0˚, which shows slight density and magnetization anomalies.Rock 2 delineates the high density and magnetization volume in the northeast of the study area.Rock 3 reveals negative density contrast and medium magnetization anomalies.Preliminary geological inferences and rock sample analyses suggested that the volcanic rocks stratum rises in the northeastern region, showing both high density contrast and magnetism, and the crater area is located beneath sedimentary overburden in the central northwestern region, showing low density contrast and high magnetism.All these lithological characteristic interpretations match the geological inferences.

Figure 7 .
Figure 7. Lithological characterization analysis extracted from joint imaging results.

Table 1 .
Rock physics of the subsurface medium.