Methods for joint inversion and analysis of EM and other geophysical data

5 Methods for joint inversion and analysis of EM and other geophysical data 5.1 Introduction The ultimate goal of the acquisition and analysis of geophysical data is typically to evaluate the composition, lithology, or petrology of the rocks, their fluid saturation, porosity, permeability, etc. These “hyperparameters” (or second-level parameters) are in some way related to the first-level parametersdthe physical properties of the rocks (electrical conductivity, density, temperature, etc.) which, in turn, are directly or indirectly determined from the measured geophysical data. The interface between the model spaces of the first and second levels is ensured by petrophysical data and geostatistics while between the first-level model space and the data, by the solution of the direct problems for the respective methods (Fig. 5.1). Depending on the objectives of the study, various algorithms are used for constructing the single-method models based on the existing data and prior geological and geophysical information. In doing so, the researcher faces the well-known difficulties: - insufficient resolution of the measured data compared to the target parameters; - relatively low accuracy and low degree of detail of the constructed models associated with insufficient amount and/or quality of the measured data and prior information; - theoretical nonuniqueness and instability of solution of the inverse problems of geophysics; - insufficient efficiency of the existing inversion algorithms for separate data types (which are particularly relevant for the 3-D/4-D problem statements); - the absence of mechanisms for quantitative integration of the entire set of the available prior information and expert estimates in the geophysical inversion.

Computational Geo-Electromagnetics. https://doi.org/10.1016/B978-0-12-819631-1.00005-5 Copyright © 2020 Elsevier Inc. All rights reserved.

133

134

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

Figure 5.1 Data and model spaces during joint inversion of geophysical data. Modified after Bosch, 1999.

These difficulties are, to some degree, coped with in the inversion of an individual data type. However, the final results frequently dissatisfy the geophysicists and geologists; they try to reduce the uncertainties and, accordingly, to improve the reliability of the constructed models by involving different data types (e.g., Bedrosian, 2007; Spichak et al., 2006, 2007; Fullea, 2017). To this effect, attempts are made to build a self-consistent and robust model of a geophysical/geological system by integrating independent data sets of different physical type, with different scales, resolution, quality, etc. Naturally, this approach has its own challenges: - In the general case, it is not guaranteed that the data included in the joint processing have been acquired from the same segment of the study area and/or measured with the same degree of detail. - Different data types may have different resolution with respect to the sought parameters of the medium. - Rock properties or hyperparameters determined from the different geophysical data are not guaranteed to be consistent with each other, etc.

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

If the described risk factors are neglected, the joint inversion of the different data types may even worsen the sought result rather than improve it. Moreover, in the real (rather than model) situation, the researcher may not even be aware of it because there are no objective criteria to judge how much better the model constructed in this way is than the models constructed without combining different data types. Therefore, critical analysis of the methods for joint geophysical inversion is not only of a pure academic interest but also of practical importance. Despite availability of first monographs on this item (Dell’Aversana, 2014; Moorkamp et al., 2016) it is still actual to provide a comparative analysis of different approaches highlighting their advantages and limitations. In this chapter different methods for joint inversion and posterior analysis of the results of singlemethod inversions are considered following a review paper by Spichak (2019).

5.2 Simultaneous inversion The idea of simultaneous inversion is that in the hope to narrow the margins of uncertainty and to increase the accuracy, the models of the first- or second-level parameters are constructed simultaneously (in parallel) over the entire set of the initial data. Each model is built based on its own data, and the relationship between them is implemented through prior information which is artificially introduced into the inversion process. Depending on the mathematical apparatus applied, deterministic and stochastic approaches are distinguished.

5.2.1 Deterministic techniques The simplest way for simultaneous inversion of all the available data types consists in forming the residual functional with its subsequent minimization within the framework of a certain computational method. The idea of this joint inversion is that the researcher may vary, at one’s own discretion, the influence of some group of the data on the final results (typically, firstlevel models) by adjusting the respective weighting coefficients (the regularization parameters) in the residual functional (Lelievre et al., 2012). In publications by Abubakar et al. (2012) and Gao et al. (2012) this approach is implemented based on the modified Newtone Gauss iterative method (Habashy and Abubakar, 2004). These authors consider the so-called petrophysical inversion. Assuming that the second-level parameters (including porosity and fluid

135

136

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

saturation of rocks) are linked with the first-level parameters (electrical conductivity, seismic velocities, and density) by the known empirical relations, they initially include them into the objective residual functional to be minimized. The joint inversion of the electromagnetic and seismic data immediately yields the second-level models: the porosity and fluid saturation ones. The authors show that the accuracy of the two constructed models is higher than if they were simultaneously built based on the single (their “own”) data type. At the same time, the obtained results are indecisive as to the effect of the inversion in terms of the second-level parameters compared to the variant where the inversion is conducted in two stepsdinitially for the first-level parameters and then for the parameters of the second level. In the latter case, there is no need to immediately postulate the empirical dependences of hyperparameters on the first-level parameters which could distort the results of the inversion. In some cases, there are grounds to expect that there is a structural similarity between the models of the first-level or second-level parameters (in particular, a common geometrical model can be hypothesized a priori (Golizdra, 1978)). In this case, the “curvature” or gradients of the physical properties characterize the geometrical peculiarities of the model and are used as prior constraints between two different data sets which are formalized and forcedly taken into account in the inversion. For example, Haber and Oldenburg (1997) defined the structure operator of the model in the following way:

0;
V2 m
< s1



SðmÞ ¼ P5
V2 m
; s1 <
V2 m
< s2 (5.1)
2
1; s2 <
V m
where m is the vector of model parameters and P5 is a polynomial of fifth degree which is selected in such a way that S be a continuous twice differentiable function. As follows from definition (5.1), the structure operator is always positive and normalized to the interval [0, 1]. This guarantees that, first, both the positive and negative variations of the model are taken into account identically and, second, operator S is invariant to the scale of the models. Fig. 5.2 demonstrates two examples of application of this structure operator to the model in the case when s1 ¼ s2 ¼ 105.

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

137

Figure 5.2 Results of application of the structural operator: (A)dmodels, (B)dresulting images (s1 ¼ s2 ¼ 105) (Haber, Oldenburg, 1997).

Correspondingly, the requirement of structural similarity between two models is reduced to solving the minimization problem for the difference between the structures of two models, or penalty function 41 ¼

N X   i  2 S m1  S mi2 0min;

(5.2)

i¼1

where mi1 and mi2 are the values of model parameters in the ith grid cell and N is the total number of cells, provided that the data constraints are satisfied. In this approach, named “structural inversion,” one of the difficulties is that the subjective choice of the thresholds s1 and s2 which quantify the “large difference” may considerably affect the solution. For example, if s1 is too small while s2 is too large, the entire model space will be deemed having a structure. Another shortcoming of this approach was that during the joint inversion, the trend of the changes in the structure was not taken into account. Gallardo and Meju (2003) have overcome these difficulties by introducing a nonnormalized cross-gradient function which, unlike operator (5.1), defines the geometrical similarity between two models as the distribution of the changes in their properties rather than in the quantities themselves: 42 ¼

N X   Vmi1  Vmi2 0min

(5.3)

i¼1

In other words, these authors propose the inversion with structural constraints that are determined by the gradients of the existing models. Note that the zero value of the vector

138

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

product (5.3) can mean both the complete collinearity of the gradients and their absence at least in one model. Fig. 5.3 illustrates the use of this concept of geometrical similarity between two models of the physical properties (p (a) and q (b)) (Gallardo et al., 2005). The vectors on the contour maps (a) and (b) show the gradients of the properties in the corresponding zones of the models. Vectors (1) and (2) on the both maps have significant amplitudes but point in different directions (almost perpendicular to each other) which means the practical absence of structural similarity. Vectors (3) also have large amplitudes; however, they are directed exactly opposite to each other which means the presence of structural similarity. The contour map (c) shows the values of cross-gradient function calculated for models p and q. It can be seen that the largest positive or negative cross-gradient values are confined to the regions with smallest structural similarity. Fig. 5.4 shows the example of joint inversion of magnetotelluric and seismic data with the use of the described approach (Gallardo and Meju, 2007). For the electromagnetic and seismic models (Fig. 5.4A,B, respectively) which have a similar structure with the characteristic structural elements (the sedimentary cover, the Earth’s crust, the upper mantle, an isolated reservoir, an inclined fault, and a horst), the inversions of the corresponding synthetic data contaminated with noise were conducted. The results of the one-dimensional (1-D) inversions are shown in Fig. 5.4C,D, and the results of the respective joint inversion by the crossgradient procedure are illustrated in Fig. 5.4E,F. By comparing

Figure 5.3 Illustration of the concept of geometrical similarity between two images using schematic p (A) and q (B) images, both of arbitrary units from 0 to 10 (Gallardo et al., 2005).

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

139

Figure 5.4 The models obtained from separate and joint inversion of seismic and MT data (Gallardo, Meju, 2007): (A) MT-resistivity model recovered by separate inversion using the formulation of Smith and Booker (1991) implemented in this study. (B) Seismic model obtained by separate inversion of the seismic data using the scheme implemented by Gallardo (2004). (C) MT-resistivity and (D) seismic velocity models obtained by joint inversion using cross-gradient constraints. The thick lines outline the test structures that are the exploration targets (units AeF).

the constructed sections one can see that the joint inversion provides a better vertical resolutions and better resolves the structures beneath the isolated reservoir. Despite numerous publications that use this inversion scheme (Haber and Oldenburg, 1997; Pinheiro et al., 1997; Kaipio et al., 1999; Gallardo and Meju, 2003, 2007; Gallardo et al., 2005; Saunders et al., 2005; Linde et al., 2006; Fregoso and Gallardo, 2009; Hu et al., 2009; Doetsch et al., 2010; Infante et al., 2010; Moorkamp et al., 2011; Gallardo and Meju, 2011; Hamdan et al., 2012; Lochbuhler et al., 2013, etc.), it should be noted that numerical implementation of the suggested approach is associated with a number of difficulties among which the instability of the scheme and the necessity to use a proper regularization are perhaps the main ones. That is why Galardo et al. (2005) recommended to use the simplest possible joint structures in the inversion, which, naturally, reduces the practical value of the method. Since the regularization uses second-order tensor (covariance of the model), the numerical implementation of the algorithm is based on quadratic programming. The latter supports parallel varying of the corresponding properties of two models (the structural similarity) through the use of the Lagrange multipliers and

140

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

does not limit the values of the physical properties: they can be any real numbers. These values form the additional information that constraints the models and can reduce the uncertainties in their reconstruction. The main advantage of the structural approach to the joint inversion is that there is no need in this case to a priori determine the functional relationship between different physical properties (although the known coefficient of correlation can be included into the covariance matrix of the model): they are obtained a posteriori. However, this does not necessarily mean that a useful model has been constructed and that the obtained petrophysical relationships are close to the actual ones. Moreover, in each problem it is necessary to carefully evaluate the validity of the minimization of the cross-gradient function using the physical models of rocks and the information about the geology of the region. But the most critical disadvantage of this approach is the necessity to a priori postulate the structural similarity of the considered models. It is also worth mentioning that neither the Lagrange multipliers nor quadratic programming are suitable for solving the three-dimensional (3-D) inversion problems because in this case the computer memory could be insufficient.

5.2.2 Stochastic techniques Stochastic methods, also frequently referred to as Monte Carlo methods, are based on generating a large number of realizations of a stochastic process and making the decisions on their acceptability. The general idea is to simulate the processes known from nature or physics which have the same probabilistic characteristics (Sambridge and Mosegaard, 2002). For doing this, a frequent practice is to use the Gibbs sampling (Geman and Geman, 1984), genetic algorithms (Holland, 1975; Haupt and Haupt, 2004; Moorkamp et al., 2006, 2007; 2010), and the simulated annealing algorithms (Kirkpatrick et al., 1983; Cerny, 1985; Aarst and Korst, 1989; Bertsimas and Tsitsiklis, 1993; Harris and MacGregor, 2007; Mota and Monteiro-Santos, 2010). The Bayesian statistical approach, first proposed for joint inversion of geophysical data by Goltzman and Kalinina (1973) is perhaps a most flexible instrument which allows the data and the prior geological and geophysical information to be reconciled with constructing the adequate models of the medium (see in this connection Chapter 2). Subsequently, this approach was further developed in Tarantola (1987); Backus (1988); Press (1989); Mosegaard and Tarantola (1995); Bosch (1999); Spichak et al. (1999); Bosch et al. (2001); Munoz et al. (2010); Dell’Aversana et al.

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

(2011); Tondi et al. (2012); JafarGandomi and Binley (2013); Jardani et al. (2013); MacCalman et al. (2014); Mellors et al. (2014); Ren et al. (2017); etc. Within this approach, Bosch (1999) proposed to carry out a joint inversion using sets of geophysical data, petrophysical and geostatistical information, and geological data describing the structure of the model in statistical terms as the input. As the output, the inversion yields the single-method first-level models which are explicitly linked with the geophysical data by the physical theory as well as a lithological (or some other) second-order model. Denote the set of all the first- and second-level parameters by m ¼ {m1, m2,., mN} and the corresponding set of all the model parameters by M (m ˛ M). Let mp be the vector of the first-level model parameters (mp ˛ Mp, p ¼ 1,.,k), and ms be the vector of the second-level model parameters (ms ˛ Ms where s ¼ k þ 1,.,N). Correspondingly, M ¼ Mp.Ms. Different types of the information about the properties of the medium and their spatial structure can be expressed in terms of the probability density function (PDF) defined on the sets of model parameters. The purpose of Bayesian inversion is to calculate the posterior probability density function   ppost ms jmp Þ ¼ cqs=p ðms jmp ppr ðmp ÞLðms Þ; (5.4) where each term corresponds to its own information content: - Prior geological information on the lithological structure of the medium is described by function ppr (mp) defined on the set of first-level model parameters. - qs/p (ms j mp) is the conditional probability that reflects the information on the second-order model parameters, their spatial interrelation, cross-correlations, and dependence on the firstlevel properties. - L (ms) is the joint likelihood function which is a probabilistic measure of the difference between the observations and the data calculated for the joint model. - c is the normalization constant calculated from the probabilistic formulation of the premise postulating the existence of the solution of the inverse problem on the space of all parameters (M). Thus, the resultant probability density function ppost (ms j mp) is composed of two factors corresponding to two different sources of information: the joint likelihood function L (ms)dthe product of the independent likelihood functions associated with each geophysical data setdand the joint prior PDF. In turn, the latter,

141

142

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

in accordance with the partition of the space of model parameters into the primary (physical) and secondary (lithological) subspaces, consists of the marginal PDF ppr (mp) in the primary space of model parameters and the conditional PDF in the secondary space of model parameters. The conditional PDF qs/p (ms j mp) is especially useful for accommodating the petrophysical (empirical or theoretical) laws that link the properties of the rocks with each other. The physical properties of rocks naturally depend on lithology. They reflect macroscopic manifestations of the nature of rocks (their composition, texture, structure, genesis), which is described by the lithological properties. This dependence has been extensively explored for various geological media and can be empirically examined for a particular region. The marginal probability density ppr (mp) is useful for describing the properties that are better constrained by prior information. It is typically represented in lithological terms (geological maps, probable lithotypes, geometric relationships between lithotypes, stratigraphic directions in the sedimentary formations, etc.). Within each lithology, statistical relationships between the properties (e.g., mean values and variograms or, more generally, the marginal and conditional PDFs) can be described more accurately. Under favorable conditions it can also be guessed that these statistical relationships are uniform. This is a vitally important assumption because most models for estimating and modeling the properties of the medium rely on the hypothesis of statistical uniformity, or spatial stationarity. If it is anticipated that within each lithotype the properties are uniform, then the conditional PDF of the secondary parameters will be trivial and the inversion will be reduced to estimating the primary model parameters. In particular, if only one property of the medium and only one type of geophysical data are considered, then the joint PDF will be reduced to the simple form of the posterior PDF for the single-method geophysical inversions (see the examples of the Bayesian inversion of electromagnetic data in Roussignol et al., 1993; Grandis, 1994; Spichak et al., 1999; and references in Chapter 2 of this book). As noted above, depending on problem statement, the firstlevel models can in principle also include other properties that characterize the texture of the rocks or their state (porosity, permeability, fluid saturation, etc.). The important primary properties for describing the medium depend on the scale and on the content of the problem under study: for example, the upper crust

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

is described by the lithotypes; the characteristic attributes of sedimentary basins are lithotypes, porosity, fluid content, and stratigraphic trend. Generally speaking, the selection of the first-level model parameters should meet two criteria: - these parameters should be useful for efficient identifying and characterizing the lithology or the state of the medium; and - the second-order parameters should statistically depend on them. Due to the complicated character of all functions contained in Formula (5.4), the posterior marginal PDF for model parameters does not have a closed form and can only be calculated numerically with the use of the well-known statistical methods. After its calculation, posterior values of the first- and second-level parameters are readily found based on the set of the joint models selected during the Bayesian inversion. Figs. 5.5 and 5.6 illustrate the steps of the joint inversion of magnetic and gravity data for constructing the lithological sections (Bosch et al., 2001). Fig. 5.5 shows the geological map of the study region (top) and the preliminary lithological cross sections along the NWeSE and NEeSW directions (bottom). Fig. 5.6 illustrates the steps of solving the inverse problem on finding the posterior lithology of the study area and the joint distribution of magnetic permeability and density of masses fitting the observed magnetic and gravity data, respectively, and the prior geological and geophysical information. Thus, one may conclude that stochastic approaches prove more suitable for joint interpretation of various types of geophysical data than deterministic approaches, as they provide the formalism for involving the prior geological and geophysical information as well as expert estimates. Moreover, among the output products, the researcher not only obtains the distributions of the sought parameters but also the posterior estimate of their uncertainty at each node of the spatial grid. On the other hand, the efficiency of stochastic approaches directly depends on how adequately the characteristics of the random processes in the discussed algorithms reflect the real situation. Among the challenges associated with this approach we note the necessity to specify prior PDFs for all parameters, the assumption of their Gaussian distribution (though not rigorousd see Section 2.2.1), and the difficulties in its numerical implementation which requires huge memory and computing time of the order of a few weeks even in the case of using the multiprocessor systems.

143

144

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

Figure 5.5 Plan view and two cross sections of a preliminary 3-D model of the study area (Bosch et al., 2001).

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

145

Figure 5.6 Prior and posterior joint models for the Eastern section indicated in Fig. 5.5 (Bosch et al., 2001). Lithotype color codes are same as in Fig. 5.5. The secondary properties are the base 10 logarithm of the mass density (kg/m3) and the base 10 logarithm of the magnetic susceptibility (SI units). Calculated data are shown with a blue line and observed data with a red line; the yellow band indicates 1 standard deviation uncertainty in the observed data.

5.3 Cooperative inversion An alternative approach to the simultaneous inversion methods considered above consists in the so-called “cooperative” inversion when the results of inversion of some data are used as the starting models for inverting other data. This direction in the joint geophysical inversion has been pioneered by Lines et al. (1988). In the beginning of the 2000s, this approach was reincarnated owing to the studies by Dell’Aversana and his followers (Dell’Aversana, 2001, 2006; 2014; Dell’Aversana et al., 2002, 2011; Zhu and Harris, 2011; Paasche et al., 2012). Let us illustrate the ideas of this approach by the example of the joint interpretation of seismic, magnetotelluric, and gravity data (Dell’Aversana, 2001) (Fig. 5.7). The algorithm of iterative operations is following.

146

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

Figure 5.7 Cooperative inversion of seismic and MT data (DellAversana, 2001, Figs. 1e5, 7): (A)dexample of common receiver gather (Global Offset data); (B)dtomographic Vp model by Global Offset data inversion; (C)dempirical relationship between resitivity and velocity, obtained by well log analysis; (D)dresistivity model derived from seismic tomography and borehole information; (E)dfinal magnetotelluric model, which honors the MT data and, at the same time, derives from a starting model based on seismic and borehole data; (FeA) parametric model (in colors) imported into the seismic section (time domain), the colors represent velocities derived by resistivity back-transformation.

(1) Based on the Global Offset data (Fig. 5.7A), the seismic velocity Vp section is built (Fig. 5.7B). (2) The Vp section is transformed into the electrical resistivity (R) section with resistivity logging data taken into account.

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

To this end, the resistivity logging data are analyzed for the presence of the probable correlation between the two physical parameters. The graph of the dependence between velocity VP and resistivity R is constructed (Fig. 5.7C); this dependence is approximated by the analytical formula Vp ¼ aðlnðlnðRÞÞÞ þ b

(5.5)

Based on this formula, seismic velocities of the model in Fig. 5.7 are converted into electrical resistivities (Fig. 5.7D). (3) Modeling of the MT field in the constructed starting resistivity model for separating the TE and TM modes and correcting the static shift. (4) Two-dimensional (2-D) MT inversion is carried out. The inversion starts from the resistivity model selected as a result of numerical modeling at step (3). Here, the resistivity logging data are used as the constraints. The final model is shown in Fig. 5.7F. (5) The quality of the geoelectrical model (Fig. 5.7E) is controlled by gravity modeling. To this end, the resistivity section is converted into the density section by the empirical formula which is derived from the borehole data and geological hypotheses. Although the solution of the gravity problem is fundamentally nonunique, by doing so we may constrain the variations in geometry and density and, thus, to reduce the number of the solutions. (6) The constructed resistivity section (Fig. 5.7E) is transformed into a new seismic velocity section (by Formula 5.5) and then the depth section is converted into the time section (Fig. 5.7F). If necessary, based on the new seismic section, a new resistivity model is constructed, etc. The iterative cycle (1)e(6) continues until the residuals of the inversions get stabilized for all data types. This scheme results in the construction of the “generalized multiparameter geophysical model” which, in the opinion of the author, is easier to analyze compared to the initial model (compare Fig. 5.7A and F). This approach has the advantages that its numerical implementation is relatively simple and that the computer memory requirements are much weaker than in the simultaneous inversion of all data. At the same time, during the iteration process it is necessary to use prior relations between the petrophysical properties which have been established beforehand based on the correlated logging data or results of laboratory experiments on rock samples. A weak point of this approach is that this empirical

147

148

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

relationship established from the logging or laboratory data is not necessarily valid for all the regions of the 2-D and even less the 3-D section. Besides, these data are not always available.

5.4 Classification methods A fundamentally different approach to the joint analysis of the geophysical data consists in the following. Initially, the singlemethod cross sections (the first-level models) are constructed; then, the sections characterized by the uniformity of physical properties are identified using the classification methods; and, finally, the parameters of the second-level models are forecasted with the allowance for the expert estimates and prior geological and geophysical information. This approach is based on the fundamental assumption of geophysics overall that geological objects are characterized by their own physical properties and can be distinguished by their individual measurements. Thus, the structures can be discriminated provided that their distinctions in the parameter space are greater than the variations of the physical parameters within the structures. In this section, we briefly consider the main approaches to the joint posterior analysis of the results of single-method inversions leaving aside the inversion methods for separate types of geophysical data. There are several popular classification methods (e.g., Reimann et al., 2008): cluster analysis, Gaussian classification, K-means clustering, and discriminant analysis. Let us recall their brief characteristics. Cluster analysis. This is a method of grouping based on the closeness of samples in the space of properties. Here, grouping is exclusively based on the properties and leaves out of account the information about the belonging of the cases to a particular lithological group (e.g., Hartigan, 1975; Kaufman and Rousseeu, 2005). This distinguishes this method from the other methods in which a training set is used for constructing the statistical model for classification. Gauss classification. It is assumed that within each lithological group, the parameter vector has n-dimensional Gaussian distribution of the conditional probability density function. The centroid and the covariance matrix are determined from the training set which should exist for each lithotype (e.g., Rasmussen and Williams, 2006). K-means clustering. This is an iterative method for finding a set of centroids that provide the best fitting of the distribution of the training observations (e.g., Kanungo et al., 2002). The

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

149

method attributes each training observation to one of K clusters in such a way that each cluster is represented by the centroid of the corresponding observations and each observation is closer to the centroid of its cluster than to the centroids of the other clusters. Discriminant analysis. This method seeks to find the functions of the properties that optimally (in the mean-square sense) discriminate the groups. The purpose of the method is to select a combination of linear coefficients that maximizes the variance of the group centroids and at the same time minimizes the variance of the cases within the groups (e.g., Huberty, 1994). The discriminant functions are used as an auxiliary coordinate system for data representation - the discriminant space. Among the methods based on the most popular cluster analysis, in turn, two main approaches can be noted: the probabilistic and neural network ones.

5.4.1 Probabilistic clustering According to Bedrosian (2007), for applying probabilistic clustering, the parameters of the models obtained in the singlemethod inversions should be initially interpolated on the common coordinate grid (Fig. 5.8). Here, it is assumed that each

Figure 5.8 Probabilistic approach to structure classification using independent geophysical models (Bedrosian, 2007). Independent models are interpolated (1) onto a common grid. The correlation between the models is examined (2) and significant classes (localized regions of high correlation) are identified (3). Classes are then mapped back onto the depth section (4) and interpreted as distinct lithologies (Bedrosian, 2007).

150

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

first-level model is specified by its discrete values at the nodes of its coordinate grid. Thus, after the interpolation on the common grid, its each node is characterized by the set of the values of two or more physical parameters. At the second step (“correlation”) in the common parameter space, a joint probability density function is determined; this function is used at the third step (“classification”) as the bases for identifying local domains of the increased probability density (identification of classes and their boundaries in the common parameter space). Finally, at the last step (“structure mapping”), these domains are mapped back into the coordinate space in which each class determines the geological structure (lithology). Fig. 5.9 illustrates the application of this approach based on the example of Maercklin et al. (2005). Bedrosian et al. (2007) used the approach based on the Gaussian clustering. The classification in the cited paper was conducted under the assumption that lithotypes are spatially connected domains characterized by uniform physical properties which are normally distributed about a mean. The optimal number of classes (clusters) was determined by analyzing the global

Figure 5.9 Geophysical subsurface structure around the Arava Fault (AF) shown by arrow (Maercklin et al., 2005): (A)dslices through the 3-D velocity Vp model (left) and 2-D electrical resistivity (log r) models (right); (B)d histogram of Vp vs log r (good correlation is shown in black, and ellipses outline clusters); (C)ddata pairs (pixels) from the outlined regions in the histogram (B) remapped into the subsurface (cluster W is shown in dark blue, W in light blue and E in red).

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

misfit as a function of the suitable number of classes. In papers by Munoz et al. (2010) and Jousset et al. (2011), classes were associated with geological objects in accordance with geological survey, borehole data, and independent geophysical studies. The correlation of independent models has shown that no simple empirical relations exist between different physical properties; instead, there are a number of localized correlations each of which expresses separate lithology.

5.4.2 Neural network classification In contrast to the clustering methods considered above, the neural network classification methodsdmaximal correlation similarity technique (Spichak et al., 2006, 2007) and method of SelfOrganizing Maps (SOM) (Kohonen, 2001)dare based on training of the artificial neural network (ANN) by correspondence between the considered physical parameters.

5.4.2.1 Maximal correlation similarity technique Similar to the probabilistic methods discussed above, application of this approach requires preliminary interpolation of the firstlevel models onto a common grid. Then, the domains in the coordinate space are sought in which the considered physical properties have the maximal coefficient of correlation. These clusters can be considered as the regions of spatial stationarity and be used as the basis for the subsequent differentiation of the model into the lithological types. The algorithm of the method is following. Let there be given a sample E ¼ fxg consisting of n pairs x ¼ ðx1 ; x2 Þ of the values of, generally speaking, random quantities. The sample correlation coefficient of these quantities is    P P P n x1 x2  x1 x2 x˛E x˛E x˛E r ¼ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (5.6)   2ffisffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ffi P 2 P 2 P P n x1  n x2  x1 x2 x˛E

x˛E

x˛E

x˛E

It is required to find a subset A4E of m  n points so that the “partial” correlation coefficient    P P P m x1 x2  x1 x2 x˛A x˛A x˛A ffi s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (5.7) rA ¼ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   2 2 P P P 2 P 2 m x1  x1 x2 m x2  x˛A

x˛A

x˛A

x˛A

151

152

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

be maximal. This problem is solved in the following way. Let q ðxÞ be a characteristic function of a set A, i.e.,

0; if x;A; qðxÞ ¼ 1; if x˛A: Then, the formula for the partial correlation coefficient (5.7) can be rewritten in the form:      P P P qðxÞ x1 x2 qðxÞ  x1 qðxÞ x2 qðxÞ x˛E x˛E x˛E x˛E ffi s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rA ¼ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi     2     2ffi P P 2 P P 2 P P qðxÞ x1 qðxÞ  x1 qðxÞ qðxÞ x2 qðxÞ  x2 qðxÞ 

x˛E

P

x˛E

x˛E

x˛E

x˛E

(5.8)

x˛E

Let us now consider the situation when q ¼ q ðx; aÞ in Formula (5.8) is a smooth function of x and of a certain set of parameters a, which takes on the values in the interval [0, 1]. Then, quantity rA is not necessarily exactly equal to the partial correlation coefficient of a certain subset A4E, because for some x values it is possible that 0 < q ðx; aÞ < 1. But, on the other hand, quantity rA becomes now a smooth function rðaÞ of parameter a and we can find the a value at which rðaÞ is maximum using standard algorithms for finding the maximum of a function. After this, the subset can be defined by the following condition: x˛A when and only when q ðx; aÞ > 0:5. Moreover, there is every reason to expect that the corresponding partial correlation coefficient will be close to maximal. Since this mathematical problem clearly has a nonunique solution, the required subset can be found by imposing a natural constraint of one-to-one correspondence between the values of two functions. Fig. 5.10 illustrates the application of this approach for determining the potentially seismically active regions based on the joint analysis of the resistivity model and density of the hypocenters of the previous earthquakes. As can be seen from Fig. 5.10, the regions of maximal correlation (marked by red) elementary cells are highly resistive domains of the brittle crust and the deep faults. Compared to the probabilistic techniques used by Bedrosian (2007), Bedrosian et al. (2007), and Maercklin et al. (2005), the

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

153

Figure 5.10 3-D resistivity model of the northern Tien Shan area (Spichak et al., 2006). Elementary cells (in red) indicate locations of maximal correlation between the bulk resistivity and EQ hypocenters density.

discussed approach to the analysis of geophysical data has the following advantages: - The values of two physical properties at the same nodes of the spatial grid are determined by neural network interpolation rather than by standard kriging procedure, which increases the accuracy of the determination. - The classification is conducted directly in the coordinate space which obviates the need in two steps of the scheme of the socalled probabilistic clustering (the transitions from the coordinate space to the parameter space and back) (Fig. 5.8). In turn, this makes it possible to avoid solving the unnecessary problems arising in this method (in particular, to assume Gaussian distribution of the parameters, which is not always justified). - Unlike in other clustering methods, it is not required to specify the number of clusters a priori: they are obtained automatically. On the other hand, the fact that at most two physical properties could be analyzed simultaneously is the drawback of this technique.

5.4.2.2 Self-organizing map (SOM) technique Another neural network method is based on the use of the socalled self-organizing maps (Kohonen, 2001) or unsupervised artificial neural networks. The idea of the method is to introduce the prior information about the number of clusters and to train the Kohonen ANN to identify spatial segments with the same type of characteristics in the study area. Let us illustrate the substance of this approach by the example from Bauer et al. (2012) (Fig. 5.11).

154

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

Figure 5.11 Concept for the SOM analysis of seismic and magnetotelluric models (Bauer et al., 2012). Three model types and parameters with two hypothetical rock types are considered (A). The discrimination of the rock types is possible in the 3-D parameter space, but separation of clusters will be more complicated with increasing number of rock types and if more than three parameters would be analyzed (B). Data vectors are generated from the models to form the input of the self-organizing map. Rock types are finally identified at the Kohonen layer, which represents a topological mapping of the rock-type properties from the 3-D parameter space to the 2-D Kohonen layer (C).

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

The paradigm of the Kohonen neural network does not include the comparison of the output layer neurons with reference values (which is the case in the well known backpropagation technique of supervised training, see Section 3.2). During training, this neural network is successively inputted by training examples which are the values of physical properties of the rocks at the same nodes of the coordinate grid. It can be said that the system works on the competitive basis. Neurons of the second (cluster) layer representing the hyperparameters compete with each other by the winner-gets-it-all principle, and the neuron element whose weight vector is closest to the input vector (in this case this is a three-component vector) wins this competition. During this process, each input vector is classified to a certain cluster element. After the training process has ended, the Kohonen network can classify the input examples into the groups of similar elements. The entire set of neurons in the output layer exactly models the distribution structure of training samples in the multidimensional parameter space. Thus, the Kohonen selforganizing maps transform the multivariate space of the physical properties into the space of clusters (e.g., Spichak et al., 2015). The approach considered above was applied by Spichak et al. (2008) for building a cluster petrophysical cross section based on three sections of the physical properties of the medium (seismic velocity, effective density, and electrical resistivity) constructed from the geophysical data measured in East Siberia along a segment of regional profile 1-SB (the details are presented in Section 11.2).

5.4.3 Hybrid approaches In a number of studies, the construction of a cluster section is preceded by building of the general geometrical model of the medium. For doing this, Hellman et al. (2017) use the method of cross-gradients considered earlier in Section 5.2.1. In several works (Nikitin et al., 2003; Kaplan et al., 2006; Cheremisina et al., 2004; Galuev and Kaplan, 2009), the preliminary geometrical model is constructed based on locating the regions of sharp changes in the properties of the medium in the single-method models. For better focusing of these properties and for passing to dimensionless units, the authors use the so-called “differentially normalized parameter” (DNP) which is a depth increment of logarithmic seismic stiffness (density/magnetization/ resistivity): K ¼ 1=2 d=dz ðln PÞ dz;

(5.9)

155

156

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

where P, depending on the data under consideration, is V$2 (V is seismic velocity and is 2 density), for gravimetry it is density 2, and for electrical prospecting it is electrical resistivity r. As DNP (5.9) reflects the degree of relative inhomogeneity of the medium, tracking its extreme values along the section can determine the geometry of the interfaces (blocks) characterized by the uniformity of the physical properties. Fig. 5.12 shows the “geometerized” model of the Earth’s crust constructed in the described way for the segment of the reference profile 1-SB (against the single-method models) (Kaplan et al., 2006). At the next step, the statistical characteristics of each property within the identified blocks are estimated (by the K-means technique) and the physicalegeological model with a set of physical properties averaged on the common geometry of block boundaries is constructed. At the final step, each set of the properties in the block is set in correspondence with a first-level model parameter (most frequently, lithotype) (Fig. 5.13). This approach has the advantage that the boundaries of the domains with uniform physical properties are drawn based on the contour maps of DNP extrema which naturally outline zones of smooth variations in these properties. In the cases when the boundaries coincide for the different first-level models, it is possible to construct the so-called “geometerized” model and subsequently determine the physical properties of the corresponding uniform blocks. The disadvantage of this approach is, as is often the case, the flip side of this advantage: if the boundaries of the uniform regions determined in the described way do not coincide (which is likely to be the most frequent situation in practice), the reliability of all the subsequent constructions is called into question. Just as in the cross-gradient method considered above, if this hypothesis is at odds with reality, the revealed “joint structure” and, naturally, its details may have little in common with the real situation. Overall, it should be noted that clustering methods discussed in this section are exclusively based on the statistical correlation of the physical properties in the common parameter space and do not depend on the theoretical or empirical relationships between physical properties of rocks. Besides, the considered methods can also be useful for developing the multivariate statistical description of the relationship between rock properties and lithological groups as well as for selecting the set of the properties and geophysical data suitable for determining the lithology

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

157

Figure 5.12 Models of effective geophysical parameters (in conditional units) derived from seismic velocity (A), density (B), magnetization (C) and electrical resistivity (D) (Kaplan et al., 2006).

within a specific exploration scenario. From the methodological standpoint, this approach highlights the structural conformity of the models and provides a natural framework for regularization of joint inversion of the data which can be implemented on this basis (in particular, using some deterministic or probabilistic technique discussed earlier in Section 5.2).

158

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

Figure 5.13 Lithology clusters derived from joint analysis of geophysical cross sections shown in the Fig. 5.12 (Kaplan et al., 2006): 1dterrigenous-carbonate deposits (2 ¼ 2630 kg/m3, Vp ¼ 5800 m/s); 2dcarbonateterrigenous deposits (2 ¼ 2640e2700 kg/m3> Vp ¼ 5600e6100 m/s); 3dgneisses, shales, igneous rocks of acidic composition (2 ¼ 2670e2700 kg/m3, Vp ¼ 6500 m/s); 4dquartz-biotite, gneisses and shales, igneous rocks of medium composition (grano-diorite, diorite, monzonites) (2 ¼ 2710e2730 kg/m3, Vp ¼ 6700e7100 m/s); 5dbiotitehornblended and amphibole shales, crystalline gneisses, and igneous rocks of medium composition (2 ¼ 2800e2820 kg/m3, Vp ¼ 6500e7100 m/s); 6dgabbro, gabbro-diabases, and amphibolites (2 ¼ 2850e2890 kg/m3, Vp ¼ 7300e7500 m/s); 7dgabbro, gabbro-diabases, and amphibolites (possibly increased fracturing) (2 ¼ 2850e2890 kg/m3, Vp ¼ 6800e7100 m/s); 8dgabbro, basalts, hornblended shales (2 ¼ 2930e2980 kg/m3, Vp ¼ 7600e7800 m/s); 9ddunites, peridotite, pyroxenite (2 ¼ 3000 kg/m3, Vp ¼ 7900e8300 m/s); 10declogites (2 > 3070 kg/m3, Vp > 8300 m/s).

5.5 Conclusions The composition of geophysical data used for joint inversion should depend substantially on the initial statement of the problem in terms of the second-level hyperparameters. In turn, the latter should be determined through the dialogue between the geologist and the geophysicist who may develop the formalized criteria for searching for the solution based on their experience. Therefore, preliminary formulation of the search criteria based on the necessary and sufficient conditions (in terms of macroparameters characterizing a certain object or process) appears to be the most efficient approach. As it was demonstrated, joint geophysical inversion per se is neither a necessary nor a sufficient condition for obtaining good results (in particular, for constructing the models of the medium close to the real situation). The difficulties of using the discussed approaches are associated with the difference in the resolution and spatial scale of the geophysical data involved in the joint inversion. Another challenge in joint inversion of the data that are affected by different physical processes consists in the fact that it is practically impossible to a priori separate the cases when

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

the data are associated with correlated structures in the Earth and the cases when the things are different. When all the data are sensitive to the same physical properties, the process of joint inversion is conceptually transparent and it can be expected that the resulting model will be at least not worse than the individual models. The situation changes if the data are combined in such a way that they correspond to different physical properties (e.g., electromagnetic and seismic data). On one hand, we may hope to retrieve more information about the Earth’s interior, to reduce the number of t he acceptable models, and to suppress the noise effects. On the other hand, this approach may prove totally misleading if the data sets are independent, i.e., affected by the different structures. For example, using the cross-gradient method we may obtain an object that does not really exist, due to that mere fact that we have to postulate the existence of a common structure for all the involved data (even in regions where there is no solid correlation between physical parameters) even at the stage of statement of the problem of joint inversion. Therefore, posterior evaluation and critical analysis of the results of the joint inversion are even more important than in the case of the inversion of separate data types. From the methodological standpoint, it is more correct (at least, safe against the challenge noted above) to carry out a joint posterior analysis of the results of independent single-method inversions which would allow finding the regions of maximal correlation between the different studied parametersdthe potential indicators of certain phenomena or the clusters of petrophysical properties of the medium characterizing a certain lithology. It is important that this analysis uses the transformations of geophysical data that are most sensitive to the studied phenomena or structures. In other words, the joint posterior analysis should be preceded by the analysis of the sensitivity within each method used. In turn, the cluster model constructed by joint posterior analysis of the first-level models can be used for the subsequent, a more substantiated joint geophysical inversion. Based on the considered methods of Bayesian statistical inversion (Chapter 2) and neural network analysis (Chapter 3), a general scheme could be suggested for the joint analysis of geological and geophysical information, data inversion, and forecasting the target parameters (Fig. 5.14).

159

160

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

Figure 5.14 General workflow for geophysical data analysis and inversion.

References Aarst, E., Korst, J., 1989. Simulated Annealing and Boltzman Machines. Wiley Publ., Chichester. Abubakar, A., Gao, G., Habashy, T.M., Liu, J., 2012. Joint inversion approaches for geophysical electromagnetic and elastic full-waveform data. Inverse Probl. 28 https://doi.org/10.1088/0266-5611/28/5/055016. Backus, G.E., 1988. Bayesian inference in geomagnetism. Geophys. J. 92, 125e142. Bauer, K., Munoz, G., Moeck, I., 2012. Pattern recognition and lithological interpretation of collocated seismic and magnetotelluric models using selforganizing maps. Geophys. J. Int. 189, 984e998. Bedrosian, P.A., 2007. MTþ, integrating magnetotellurics to determine earth structure, physical state and processes. Surv. Geophys. 28, 121e167. Bedrosian, P.A., Maercklin, N., Weckmann, U., Bartov, Y., Ryberg, T., Ritter, O., 2007. Lithology-derived structure classification from the joint interpretation of magnetotelluric and seismic models. Geophys. J. Int. 170, 737e748. Bertsimas, D., Tsitsiklis, J., 1993. Simulated annealing. Stat. Sci. 8 (1), 10e15. Bosch, N., 1999. Lithologic tomography: from plural geophysical data to lithology estimation. J. Geoph. Res. 104 (B1), 749e766. Bosch, N., Guillen, A., Ledru, P., 2001. Lithologic tomography: an application to geophysical data from Cadomian belt of northern Brittany, France. Tectonophysics 331, 197e227. Cerny, V., 1985. Thermodynamical approach to the traveling salesman problem: an efficient simulation algorithm. J. Optim. Theory Appl. 45, 41e51. Cheremisina, E.N., Galuev, V.I., Kaplan, S.A., Malinina, S.S., 2004. Technique for identifying the reference deep boundaries marking the changes in the physical properties of rocks for the problems of integration of geophysical data in the regional geophysical studies. Geoinformatika 1, 50e53 (in Russian). Dell’Aversana, P., 2001. Integration of Seismic, MT and Gravity Data in a Thrust Belt Interpretation//First Break 6, pp. 335e341. Dell’Aversana, P., 2006. Joint inversion of seismic, gravity and magnetotelluric data combined with depth seismic imaging. In: Extended Abstr. 18th IAGA WG 1.2 Workshop on EM Induction in the Earth (El Vendrell, Spain).

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

Dell’Aversana, P., 2014. Integrated geophysical models combining rock physics with seismic, electromagnetic and gravity models. In: EAGE Publ. DB Houten. (The Netherlands). Dell’Aversana, P., Bernasconi, G., Miotti, F., Rovetta, D., 2011. Joint inversion of rock properties from sonic, resistivity and density well-log measurements. Geophys. Prospect. 59, 1144e1154. Dell’Aversana, P., Morandi, S., 2002. Depth model building by constrained magnetotelluric inversion. Ann. Geophys. 45 (2), 247e257. Doetsch, J., Linde, N., Coscia, I., Greenhalgh, S.A., Green, A.G., 2010. Zonation for 3D aquifer characterization based on joint inversions of multimethod crosshole geophysical data. Geophysics 75 (6), G53eG64. Fregoso, E., Gallardo, L.A., 2009. Cross-gradients joint 3D inversion with applications to gravity and magnetic data. Geophysics 74 (4), L31eL42. Fullea, J., 2017. On joint modelling of electrical conductivity and other geophysical and petrological observables to infer the structure of the lithosphere and underlying upper mantle. Surv. Geophys. https://doi.org/ 10.1007/s10712-017-9432-4. Gallardo, L.A., 2004. Joint Two-dimensional Inversion of Geoelectromagnetic and Seismic Refraction Data with Cross-gradients Constraint. PhD dissertation. Lancaster University, UK. Gallardo, L.A., Meju, M.A., 2003. Characterization of heterogeneous near-surface materials by joint 2D inversion of DC resistivity and seismic data. Geophys. Res. Lett. 30 (13) https://doi.org/10.1029/2003GL017370. Gallardo, L.A., Meju, M.A., 2007. Joint two-dimensional cross-gradient imaging of magnetotelluric and seismic traveltime data for structural and lithological classification. Geophys. J. Int. 169, 261e1272. Gallardo, L.A., Meju, M.A., 2011. Structure coupled multiphysics imaging in geophysical sciences. Rev. Geophys. https://doi.org/10.1029/2010RG000330. RG1003. Gallardo, L.A., Meju, M.A., Perez-Flores, M.A., 2005. A quadratic programming approach for joint image reconstruction: mathematical and geophysical examples. Inverse Probl. 21, 435e452. Galuev, V.I., Kaplan, S.A., 2009. Complex interpretation of the data of studies on a segment of the reference geological-geophysical profile 2-DV. Razvedka Okhrana Nedr (Prospect Prot. Miner. Resour.) 4, 49e56 (in Russian). Gao, G., Abubakar, A., Habashy, T.M., 2012. Joint petrophysical inversion of electromagnetic and full-waveform seismic data. Geophysics 77 (3). WA3eWA18. Geman, S., Geman, D., 1984. Stochastic relaxation, Gibbs distribution and the Bayesian restoration of images. IEEE Trans. Pattern Anal. Mach. Intell. 721e741. Golizdra, G.Y., 1978. On integrating the gravimetric and seismic methods. Izv. Phys. Solid Earth 6, 26e38. Golzman, F.M., Kalinina, T.B., 1973. Integration of geophysical observations. Izv. Phys. Solid Earth 8, 31e42. tique Bayesienne par la simulation Grandis, H., 1994. Imagerie electromagne , Universite Paris VII. d’une chaine de Markov. Doctoral d’Universite Habashy, T.M., Abubakar, A., 2004. A general framework for constraint minimization for the inversion of electromagnetic measurements. In: Progress in Electromagnetic Research Symp, 46, pp. 265e312. https:// doi.org/10.2528/PIER03100702. Haber, E., Oldenburg, D., 1997. Joint inversion: a structural approach. Inverse Probl. 13, 63e77.

161

162

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

Hamdan, H., Economou, N., Kritikakis, G., Andronikidis, N., Manoutsoglou, E., Vafidis, A., Pangratis, P., Apostolidou, G., 2012. 2D and 3D imaging of the metamorphic carbonates at Omalos plateau/polje, Crete, Greece by employing independent and joint inversion on resistivity and seismic data. Int. J. Speleol. 41 (2), 199e209. Harris, P., MacGregor, L., 2007. Enhancing the resolution of CSEM inversion using seismic constraints. In: Expanded Abstr. SEG San Antonio Annual Meeting. Hartigan, J., 1975. Clastering Algorithms. John Wiley & Sons, Inc., p. 351 Haupt, R.L., Haupt, S.E., 2004. Practical Genetic Algorithms, second ed. John Wiley & Sons, Inc., Hoboken, New Jersey. Hellman, K., Ronczka, M., Gunther, T., Wennermark, M., Rucker, C., Dahlin, T., 2017. Structurally coupled inversion of ERT and refraction seismic data combined with cluster-based model integration. J. Appl. Geophys. 163, 169e181. Holland, J.H., 1975. Adaptation in Natural and Artificial Systems. University of Michigan Press, Ann Arbor, MI. Hu, W., Abubakar, A., Habashy, T.M., 2009. Joint electromagnetic and seismic inversion using structural constraints. Geophysics 74 (6), R99eR109. Huberty, C.J., 1994. Applied Discriminant Analysis. John Wiley and Sons, Inc., New York, 466pp. Infante, V., Gallardo, L.A., Montalvo-Arrieta, J.C., de León, I.N., 2010. Lithological classification assisted by the joint inversion of electrical and seismic data at a control site in northeast Mexico. J. Appl. Geophys. 70, 93e102. JafarGandomi, A., Binley, A., 2013. A Bayesian trans-dimensional approach for the fusion of multiple geophysical datasets. J. Appl. Geophys. 96, 38e54. Jardani, A., Revil, A., Dupont, J.P., 2013. Stochastic joint inversion of hydrogeophysical data for salt tracer test monitoring and hydraulic conductivity imaging. Adv. Water Resour. 52, 62e77. Jousset, P., Haberland, C., Bauer, K., Árnason, K., 2011. Hengill geothermal volcanic complex (Iceland) characterized by integrated geophysical observations. Geothermics 40, 1e24. Kaipio, J.P., Kolehmainen, V., Vauhkonen, M., Somersalo, E., 1999. Inverse problems with structural prior information. Inverse Probl. 15, 713e729. Kanungo, T., Mount, D.M., Netanyahu, N.S., Piatko, C.D., Silverman, R., Wu, A.Y., 2002. An efficient k-means clustering algorithm: analysis and implementation. IEEE Trans. Pattern Anal. Mach. Intell. 24 (7), 881e892. Kaplan, S.A., Galuev, V.I., Pimanova, N.N., Malinina, S.S., 2006. Integrated processing and interpretation of survey data on reference geophysical profiles. Geoinformatika 3, 38e46 (in Russian). Kaufman, L., Rousseeu, P.J., 2005. Finding Groups in Data. John Wiley & Sons, Inc., 368pp. Kirkpatrick, S., Gelatt Jr., C.D., Vecchi, M.P., 1983. Optimization by simulated annealing. Science 220, 671e680. Kohonen, T., 2001. Self-organizing Maps. Springer-Verlag, Berlin, 534pp. Lelievre, P.G., Farquharson, C.G., Hunch, C.A., 2012. Joint inversion of seismic traveltimes and gravity data on unstructured grids with application to mineral exploration. Geophysics 77 (1), K1eK15. Linde, N., Tryggvason, A., Binley, A., Pedersen, L.B., Revil, A., 2006. A structural approach to joint three-dimensional inversion of geophysical data. In: Extended Abstr. 18th IAGA WG 1.2 Workshop on EM Induction in the Earth (El Vendrell, Spain).

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

Lines, L.R., Schultz, A.K., Treitel, S., 1988. Cooperative inversion of geophysical data. Geophysics 53 (1), 8e20. Lochbuhler, T., Doetsch, J., Brauchler, R., Linde, N., 2013. Structure-coupled joint inversion of geophysical and hydrological data. Geophysics 78 (3), ID1eID14. MacCalman, L., O’Callaghan, S.T., Reid, A., Shen, D., Carter, S., Krieger, L., 2014. Distributed Bayesian geophysical inversions. In: Beardsmore, G., Bonilla, E.V., Ramos, F.T. (Eds.), Expanded Abstr. Thirty-Ninth Workshop on Geothermal Reservoir Engineering. Stanford, California. Maercklin, N., Bedrosian, P.A., Haberland, C., Ritter, O., Ryberg, T., Weber, M., Weckmann, W., 2005. Characterizing a large shear-zone with seismic and magnetotelluric methods: the case of the Dead Sea Transform. Geoph. Res. Lett. 32 https://doi.org/10.1029/2005GL022724. Mellors, R.J., Tompson, A., Dyer, K., Yang, X., Chen, M., Wagoner, J., TrainorGuiton, W., Ramirez, A., 2014. Stochastic joint inversion modeling algorithm of geothermal prospects. In: Expanded Abstr. Thirty-ninth Workshop on Geothermal Reservoir Engineering. Stanford, California. Moorkamp, M., Heincke, B., Jegen, M., Roberts, A.W., Hobbs, R.W., 2011. A framework for 3-D joint inversion of MT, gravity and seismic refraction data. Geophys. J. Int. 184, 477e493. Moorkamp, M., Jones, A.G., Eaton, D.W., 2007. Joint inversion of teleseismic receiver functions and magnetotelluric data using a genetic algorithm: are seismic velocities and electrical conductivities compatible? Geophys. Res. Lett. 34, L16311 https://doi.org/10.1029/2007GL030519. Moorkamp, M., Jones, A.G., Fishwick, S., 2010. Joint inversion of receiver functions, surface wave dispersion, and magnetotelluric data. J. Geophys. Res. 115 https://doi.org/10.1029/2009JB006369. Moorkamp, M., Jones, A.G., Rao, C.K., 2006. Joint inversion of MT and seismic receiver function data using a genetic algorithm. In: Extended Abstr. 18th IAGA WG 1.2 Workshop on EM Induction in the Earth (El Vendrell, Spain). Moorkamp, M., Lelièvre, P.G., Linde, N., Khan, A. (Eds.), 2016. Integrated Imaging of the Earth: Theory and Applications. AGU and John Wiley and Sons, Inc. Mosegaard, K., Tarantola, A., 1995. Monte Carlo sampling of solutions to inverse problems. J. Geophys. Res. 100 (B7), 12431e12447. Mota, R., Monteiro-Santos, F.A., 2010. 2D sections of porosity and water saturation from integrated resistivity and seismic surveys. Near Surf. Geophys. 8, 575e584. Munoz, G., Bauer, K., Moeck, I., Schulze, A., Ritter, O., 2010. Exploring the Gross Schonebeck (Germany) geothermal site using a statistical joint interpretation of magnetotelluric and seismic tomography models. Geothermics 39, 35e45. Nikitin, A.A., Kaplan, S.A., Galuev, V.I., Malinina, S.S., 2003. Determination of physicogeometrical properties of the Earth’s crust from combined geophysical data. Geoinformatika 2, 29e38 (in Russian). Paasche, H., Tronicke, J., Dietrich, P., 2012. Zonal cooperative inversion of partially co-located data sets constrained by structural a priori information. Near Surf. Geophys. 10, 103e116. Pinheiro, P.A.T., Loh, W.W., Dickin, F.J., 1997. Smoothness-constrained inversion for two-dimensional electrical resistance tomography. Meas. Sci. Technol. 8, 293e302. Press, S.J., 1989. Bayesian Statistics: Principle, Models and Applications. John Wiley & Sons.

163

164

Chapter 5 Methods for joint inversion and analysis of EM and other geophysical data

Rasmussen, C.E., Williams, C.K.I., 2006. Gaussian Processes for Machine Learning. MIT Press. Reimann, C., Filzmoser, P., Garrett, R., Dutter, R., 2008. Statistical Data Analysis Explained. John Wiley and Sons Ltd., London, 343pp. Ren, H., Ray, J., Hou, Z., Huang, M., Bao, J., Swiler, L., 2017. Bayesian inversion of seismic and electromagnetic data for marine gas reservoir characterization using multi-chain Markov chain Monte Carlo sampling. J. Appl. Geophys. https://doi.org/10.1016/j.jappgeo.2017.10.004. Roussignol, M., Jouanne, V., Menvielle, M., Tarits, P., 1993. Bayesian electromagnetic imaging. In: Hardle, W., Siman, L. (Eds.), Computer Intensive Methods. Physical Verlag, pp. 85e97. Sambridge, M., Mosegaard, K., 2002. Monte Carlo methods in geophysical inverse problems. Rev. Geophys. 40 (3) https://doi.org/10.1029/ 2000RG000089. Saunders, J.H., Herwanger, J.V., Pain, C.C., Worthington, M.N., de Oliveira, C.R.E., 2005. Constrained resistivity inversion using seismic data. Geophys. J. Int. 160, 785e796. Smith, J.T., Booker, J.R., 1991. Rapid inversion of 2-dimensional and 3-dimensional magnetotelluric data. J. Geophys. Res.-Solid Earth Planets 96, 3905e3922. Spichak, V.V., 2019. Modern techniques for joint analysis and inversion of geophysical data. Russ. Geol. Geophys. 60 (12), 23e44. Spichak, V.V., Bezruk, I.A., Goidina, A.G., 2015. Constructing the threedimensional cluster petrophysical models of geological medium based on the combination of geophysical data measured on reference profiles. Razved. Okhr. Nedr (Prospect Prot. Mineral Resour. 4, 41e45 (in Russian). Spichak, V.V., Bezruk, I.A., Popova, I.V., 2008. Constructing the deep cluster petrophysical sections from geophysical data and forecasting the oil and gas bearing capacity of the regions. Geofizika 5, 43e45 (in Russian). Spichak, V.V., Borisova, V.P., Fainberg, E.B., Khalezov, A.A., Goidina, A.G., 2007. Electromagnetic 3D tomography of the Elbrus volcanic center according to magnetotelluric and satellite data. J. Volcanol. Seismol. 1 (1), 53e66. Spichak, V.V., Menvielle, M., Roussignol, M., 1999. Three-dimensional inversion of EM data using Bayesian statistics. In: Spies, B., Oristaglio, M. (Eds.), 3D Electromagnetics. SEG Publ. GD7, Tulsa. USA, pp. 406e417. Spichak, V.V., Rybin, A., Batalev, V., Sizov, Y., Zakharova, O., Goidina, A., 2006. Application of ANN techniques to combined analysis of magnetotelluric and other geophysical data in the northern Tien Shan crustal area. In: Extended Abstr. 18th IAGA WG 1.2 Workshop on EM Induction in the Earth (El Vendrell, Spain). Tarantola, A., 1987. Inverse Problem Theory: Method for Data Fitting and Model Parameter Estimation. Elsevier, New York, p. 613. Tondi, R., Cavazzoni, C., Danecek, P., Morelli, A., 2012. Parallel, ‘large’, dense matrix problems: application to 3D sequential integrated inversion of seismological and gravity data. Comput. Geosci. 48, 143e156. Zhu, T., Harris, J.M., 2011. Iterative joint inversion of P-wave and S-wave crosswell traveltime data. In: Expanded Abstr. SEG San Antonio Annual Meeting, pp. 479e483. USA.