Probabilistic approach to Gramian inversion of multiphysics data

We consider a probabilistic approach to the joint inversion of multiphysics data based on Gramian constraints. The multiphysics geophysical survey represents the most effective technique for geophysical exploration because different physical data reflect distinct physical properties of the various components of the geological system. By joint inversion of the multiphysics data, one can produce enhanced subsurface images of the physical properties distribution, which improves our ability to explore natural resources. One powerful method of joint inversion is based on Gramian constraints. This technique enforces the relationships between different model parameters during the inversion process. We demonstrate that the Gramian can be interpreted as a determinant of the covariance matrix between different physical models representing the subsurface geology in the framework of the probabilistic approach to inversion theory. This interpretation opens the way to use all the power of the modern probability theory and statistics in developing novel methods for joint inversion of the multiphysics data. We apply the developed joint inversion methodology to inversion of gravity gradiometry and magnetic data in the Nordkapp Basin, Barents Sea to image salt diapirs.


Introduction
Different geophysical methods provide information about various physical properties of rock formations and mineralization. In many cases, this information is mutually complementary. At the same time, inversion of the data for a particular survey is subject to considerable uncertainty and ambiguity as to causative body geometry and intrinsic physical property contrast. One productive approach to reducing uncertainty is to jointly invert several types of data. In the framework of joint inversion, the goal is to fit different types of predicted data to the corresponding observed data while keeping specific relationships between the different physical properties of the rocks. This can be achieved by using the known petrophysical relationships between different physical properties of the rocks within the framework of the inversion process (Afnimar et al., 2002;Hoversten et al., 2006;Moorkamp et al., 2011;Gao et al., 2012;Zhdanov, 2015;Moorkamp et al., 2016;Giraud et al., 2017;Giraud et al., 2019).
The alternative approach is based on using the Gramian constraints introduced in (Zhdanov et al., 2012a;Zhdanov et al., 2012b). The advantage of Gramian constraints is that they enforce the functional relationships between different physical parameters without a priori knowledge of the specific form of these relationships. In the original publications (Zhdanov et al., 2012a;Zhdanov et al., 2012b;Zhdanov, 2015), the Gramian constraints were introduced in the framework of the deterministic approach to the solution of the inverse problem, which considers the data and model parameters characterized by specific functions or vectors with certain (maybe unknown) values. However, there is a probabilistic approach to inverse problems where the observed data and model parameters are treated as realizations of some random variables. This approach was introduced in the pioneering papers of (Foster, 1961), (Franklin, 1970), (Jackson, 1972), (Tarantola and Valette, 1982), and (Tarantola, 1987;Tarantola, 2005).
It can be demonstrated that both these approaches result in similar numerical solutions of the inverse problem (Zhdanov, 2002). At the same time, deterministic or probabilistic interpretation of the observed data and model parameters emphasize different aspects of the inversion algorithms. This also helps understand better the properties of the inversion parameters.
This paper introduces a novel approach to the joint inversion where the Gramian constraints are represented in the probabilistic form as the determinant of the covariance matrix between the different model parameters.

Gramian constraints
Gramian constraints enforce the correlation between different parameters and their transforms or attributes (Zhdanov et al., 2012a;Zhdanov, 2015). By imposing the requirement of the minimum of the Gramian mathematical operator in regularized inversion, we obtain multi-modal inverse solutions with enhanced correlations between the different model parameters or their attributes.
Let us consider forward geophysical problems for multiple geophysical data sets. These problems can be described by the following operator relationships: (1) where, in a general case, A (i) is a non-linear operator, d (i) (i = 1, 2, 3, …, n) are different observed data sets (which may have different physical natures and/or parameters), and m (i) (i = 1, 2, 3, …, n) are the unknown sets of model parameters.
Note that, diverse model parameters may have different physical dimensions (e.g., density is measured in g/cm 3 , resistivity is measured in Ohm-m, etc.).
It is convenient to introduce the dimensionless weighted model parameters,m (i) , defined as follows: where W (i) m is the corresponding linear operator of the model weighting (Zhdanov, 2002;Zhdanov, 2015).
Note that, in the following sections of the paper we will omit "wave" above the symbols of the model parameters to simplify the notations. However, we always assume that we work with the dimensionless model parameters.
These parameters can be described by integrable functions of a radius-vector r = (x, y, z) defined within some volume V of a 3D space. The set of these functions forms a complex, in a general case, Hilbert space of the model parameters, M, with a L 2 norm, defined by the corresponding inner product: where asterisk "*" denotes the complex conjugate value.
Similarly, different data, as a rule, have distinct physical dimensions as well. Therefore, it is convenient to consider dimensionless weighted data,d (i) , defined as follows: where W (i) d is the corresponding linear operator of data weighting. However, to simplify the notations, we will drop the "wave" above the symbols of the data, assuming nevertheless that the data are dimensionless.
The selection of the model and data weights was discussed in many publications on inversion theory. For example, in the framework of the probabilistic approach (Tarantola, 1987), the weights are determined as inverse data covariance or model covariance matrices. In the framework of the deterministic approach (Zhdanov, 2002(Zhdanov, , 2015, the weights are determined as inverse integrated sensitivity matrices.
We also assume that the data belong to some complex Hilbert space of the data, D, with the L 2 norm, defined by the corresponding inner product: where S is an observation surface. Joint inversion of multiphysics data can be reduced to minimization of the following parametric functional, where α ∈ [0, ∞) is the regularization parameter. Misfit functionals, φ (m (i) ), are defined as follows, where A (i) (m (i) ) are the predicted data, The term S G (m (1) , m (2) , ….m (n) ) is the constraint stabilizing functional. This functional can be introduced as a determinant of the Gram matrix of a system of model parameters, m (1) , m (2) , ….m (n) , which is called a Gramian (Zhdanov et al., 2012a;Zhdanov, 2015): where (-,-) stands for the inner product in the corresponding Hilbert space of model parameters, defined by Eq. 3.
Note that, in case of discrete model parameters, expression (3) for inner product takes the following form: where m (i) l are the discrete values of the corresponding parameter; L is the total number of every discretized model parameter within the modeling domain V, assuming that we use the same discretization for all parameters.
The main property of the Gramian (Eq. 8) is that it is a nonnegative functional, (e.g., Everitt, 1958;Barth, 1999): The equality holds in (Eq. 10) if the system of functions (m (1) , m (2) , …., m (n) ) is linearly dependent. It was demonstrated in (Zhdanov et al., 2012a) and (Zhdanov, 2015) that one could apply the Gramian constraints to different transforms of the model parameters. For example, we may consider differential operations, like gradient, applied to the spatial distribution of the model parameters. This will result in the shared Earth model characterized by similar behavior of the gradients of different physical properties. In this case, Gramian constraint, similarly to the cross-gradient method (Gallardo and Meju, 2003, 2004, 2007Gallardo, 2007), enforces the structural similarity between different physical inverse models . We can also use various functions of the model parameters, e.g., logarithms or trigonometric functions, to enforce the correlation between the transformed parameters.

Covariance representation of the probabilistic Gramian stabilizer
In the framework of the probabilistic approach to the solution of the inverse problem, we consider the observed data and the model parameters as realizations of some random variables. The joint inversion requires the enforcement of some relations between different model parameters, which can be achieved by adding the term containing the covariance matrix between different model parameters, which serves as an analog of Gramian in a deterministic approach (Zhdanov et al., 2012a).
The probabilistic Gramian stabilizer, S G σ , can be introduced as the determinant of the covariance matrix between different model parameters in a probabilistic approach to the inversion theory: It is well known that the determinant of the covariance matrix is always non-negative, and it is equal to zero if and only if the random variables, m (1) , m (2) , …, m (n) , are linearly dependent (e.g., Ross, 2010). This property of the probabilistic Gramian is similar to that of the deterministic Gramian defined by Eq. 8. The key difference is in the way how we treat the model parameters. In the framework of the deterministic approach, they are described by specific (though unknown) functions. In the framework of the probabilistic approach, the model parameters are the realizations of some unknown random variables.
We should note also that the direct analogy between Eqs. 8 and 11 holds when the random model parameters have zero mean values. Indeed, in the case of discrete and real model parameters, the statistical estimate of the covariance is as follows: (12) where ⟨⋅⟩ indicates the mean. If the mean values are zero, ⟨m (i,j) ⟩ = 0, then according to Eq. 9, we have In Supplementary Appendix S1 we discuss the Gramian space of random variables. One of the important properties of this space is that the distance between two variables is determined by the corresponding Gramian.
Let us consider the Gramian space Γ (n) of random variables, representing model parameters with the core elements formed by model parameters m (1) , m (2) , …, m (n−1) . Then according to definition of norm in this space, equations S1 and S8 in the Supplementary Appendix we can write probabilistic Gramian as the norm of parameter m (n) : We will use this property in formulating the principles of joint inversion using the probabilistic Gramian.

Regularization of the inverse problem using probabilistic Gramian
Following the general principles of regularization theory, we can now introduce a probabilistic parametric functional, P α σ as a linear combination of combined misfit functional, n ∑ i=1 φ (m (i) ), and probabilistic Gramian: where α ∈ [0, ∞) is regularization parameter. We should note that the parametric functional is similar to Eq. 6 if the joint stabilizing functional, S G , is selected equal to the probabilistic Gramian, S G σ .
There exist several standard gradient type methods which can be used to solve the minimization problem-steepest descent, Newton, and conjugate gradient methods, for example,. All of them require calculation of the steepest ascent direction (gradient) of the corresponding functional. In the following section we will consider this problem in details.

Steepest ascent direction of the probabilistic parametric functional
The first variation of the parametric functional, P α σ , with the probabilistic Gramian stabilizer can be calculated as follows: where Γ (1) D is the probabilistic Gramian space of data with inner product defined by covariance according to formula S12 in the Supplementary Appendix.
The variation of the probabilistic Gramian stabilizer, according to Eq. 14, can be expressed as follows: where we took into account that, according to S11 in the Supplementary Appendix.
Using the properties of the inner product operation in the Gramian spaces of random variables, Γ (i) , we can calculate the variation of the norm square of m (i) as follows: where G m σij is the corresponding minor of Gram matrix G σ (m (1) , m (2) , ….m (n−1) , m (n) ) formed by eliminating column i and row j; and elements l (i) G σ are the directions of the steepest ascent for the probabilistic Gramian stabilizing functionals, We introduce now the probabilistic Gramian space of the model parameters, Γ M , with the inner product defined by covariance according to formula S12 in the Supplementary Appendix. Therefore, Eq. 19 can be written as follows: Substituting the last formula in Eq. 17, we arrive at the following expression for the first variation of the probabilistic Gramian: Let us examine again Eq. 16 for the variation of the parametric functional.
Considering that operators A (i) are differentiable, we can write: m is a linear operator of the Fréchet derivative of A (i) . Therefore, the inner product in the first term in Eq. 16 can be modified as follows: where F (i)⋆ m are the adjoint Fréchet derivative operators. Substituting Eqs 23, 22 in the first and second terms of Eq. 16, we obtain the following important result: where l α(i) σ (m (1) , m (2) , ….m (n) ) are the directions of the steepest ascent of the probabilistic functional P α σ : and r (i) are the residuals between the predicted and observed data:

Steepest descent method of joint inversion
The expression for the steepest ascent direction, introduced above, can be used in constructing the computational schemes for the different gradient type methods of solving the minimization Eq. 15.
We begin with the most simple steepest descent method. Let us select where k α is some positive real number, and l α(i) σ (m (1) , m (2) , ….m (n) ) are the directions of the steepest ascent of the functional P α σ defined by Eq. 25. Substituting Eq. 27 into Eq. 24, we have: so, the vector, describes the "direction" of increasing (ascent) of the functional P α σ , in other words, the direction of "climbing" on the hill.  To simplify the notations, we introduce vector m formed by different model parameters: and vector F m formed by the corresponding Fréchet derivative operators: We can construct an iteration process for the regularized steepest descent method as follows, where the coefficient k α n is found by using the minimization of the parametric functional with respect to k α n : In particular, applying the linear line search, we find that the minimum of the probabilistic parametric functional is reached if k α n is determined by the following formula: The iterative process Eq. 32 is terminated at n = N when the combined misfit functional reaches the given level ɛ 0 : Frontiers in Earth Science 05 frontiersin.org

Conjugate gradient method of joint inversion
The conjugate gradient method is based on the same ideas as the steepest descent, and the iteration process is very similar to the last one: where δm n = −k α nl α σ (m n ) .
However, the "directions" of ascentl α σ (m n ) are selected in a different way. At the first step we use the "direction" of the steepest ascent: At the next step the "direction" of ascent is a linear combination of the steepest ascent at this step and the "direction" of ascentl Determination of the length of iteration step, a coefficientk α n , can be based on the linear or parabolic line search: Solution of this minimization problem gives the following best estimation for the length of the step using a linear line search: One can use a parabolic line search also (Fletcher, 1995) to improve the convergence rate of the RCG method. The CG method requires that the vectorsl α σ (m n ) introduced above will be mutually conjugate. This requirement is fulfilled if the coefficients β n are determined by the formula (Zhdanov, 2002(Zhdanov, , 2015 Using Eqs 34-37, we can obtain m iteratively. The iterative process Eq. 34 is terminated when the combined misfit functional reaches the given level ɛ 0 :

Example: Nordkapp Basin
The developed algorithm of joint probabilistic Gramian inversion based on the conjugate gradient method was implemented in the computer code and carefully tested on synthetic models. We have also applied this novel method to joint 3D inversion of the gravity gradiometry, and total magnetic intensity (TMI) data gathered over the Nordkapp Basin, the principal salt-producing basin in the western Barents Sea. Finally, we compare results obtained from standalone, deterministic Gramain and the developed probabilistic Gramian inversion approaches.

Geologic setting of the Nordkapp Basin
The evolution of the Barents Sea and Nordkapp Basin is the consequence of a succession of tectonic events and climatic fluctuations that have affected the Barents Sea from the Late Devonian to the present day (Dengo and Røssland, 1992). There is a complicated distribution of structural highs, domes, platforms, and basins in the western Barents Sea, including the Nordkapp Basin.
Seismic exploration is the most prevalent method for locating the hydrocarbon (HC) deposits. Due to weak primaries, strong multiples, and diffraction noise, however, seismic imaging of salt diapirs is extremely difficult. The salt structures are surrounded by a "shadow zone" where continuous seismic reflectors are difficult to understand (Hokstad et al., 2011). As the top of the diapir is exposed to meteoric water under terrestrial conditions, the halite portion of the diapir will disintegrate, leaving behind a disordered, cumulated cap rock layer. These fragments of changing density and velocity will generate random reflections that exacerbate the difficulty of recovering the geomorphology of the diapirs using only the seismic method. Proper imaging of the diapirs from top to bottom is critical, as large salt bodies might contain small hydrocarbon volumes, but small salt bodies with an overhang can contain large hydrocarbon volumes.
To overcome this difficulty, a multiphysics approach is required to properly illuminate the complete geologic picture. Gravity gradiometry, TMI, and EM data have been acquired and inverted to fill in the gaps left by the incomplete seismic picture (Gernigon et al., 2011;Stadtler et al., 2014;Paoletti et al., 2020;Tao et al., 2021;Zhdanov, 2021, 2022). We examine the Uranus diapir in this example, which ultimately failed to present the HC deposit after drilliing. While the salt volume was large, the overhang critical to trap hydrocarbons was not present.
The G zz component of the gravity gradiometry data used in the inversions is shown in Figure 1. The black square indicates the inversion area, and the profile AA' traversing the Uranus diapir is shown by the black line in Figure 1.

Inversion parameters
The 3D voxel-based inversions were all carried out with a parallelized GPU code utilizing a moving sensitivity domain to minimize computation (Cuma et al., 2012;Cuma and Zhdanov, 2014). Voxel size was 200 m laterally with a logarithmic vertical discretization of 40-800 m spaced over 32 layers. Details of the forward modeling methodologies can be found in . All components of the gravity tensor and TMI data were inverted towards density and magnitude of magnetization vector.
The inversions were all terminated when the respective data misfits reached the error level of about 5%. The plot of data misfit versus iteration shown in Figure 2 demonstrates the probabilistic Gramian approach converged to a solution in roughly half the number of iterations as the deterministic Gramian approach. The standalone inversions converged to a solution in less than ten iterations; however, these models exhibited less structural correlation than the models derived from joint inversion. Figures 3, 4, 5 show the comparison of observed data and data predicted from the different inversion methodologies. The Frontiers in Earth Science 08 frontiersin.org data misfits of approximately 5% are comparable for all inversion methodologies.

Standalone inversion results
Vertical profiles extracted from the 3D voxel models of density and magnetization obtained from standalone inversions of the gravity gradiometry and TMI data, respectively, are shown in Figure 6. A typical density of the base tertiary rocks in the area of investigation is within 2.30-2.38 g/cm 3 . The Uranus salt diapir is characterized usually by negative density anomalies, which is clearly seen in Figure 6. The negative magnetization within the salt diapir indicates that it is opposite to the direction of the inducing magnetic field, which corresponds to a diamagnetic property of the salt structures. At the same time, the magnetization is positive outside the diapir which is typical for paramagnetic minerals present in Cretaceous sea-bottom layers of the host formations (Paoletti et al., 2020;Tao et al., 2021). Thus, the volume distribution of the density and magnetization, produced by inversion is indicative about the salt diapir structure in the Nordkapp basin. However, while the anomalies in the respective models are roughly coindicent, the gross morphology of the bodies varies widely.

Deterministic Gramian inversion results
The deterministic Gramian inversion improved the structural correlation of the inverse density and magnetization models. Compared to the models obtained from standalone inversions, we see a stronger agreement in gross morphology across the respective models. The vertical profiles extracted from the 3D voxel models of density and magnetization obtained from deterministic Gramian inversion of the gravity gradiometry and TMI data, respectively, are shown in Figure 7).

Probabilistic Gramian inversion results
Finally, the same vertical sections extracted from the 3D voxel models of density and magnetization obtained from probabilistic Gramian inversion of the gravity gradiometry and TMI data, respectively, are shown in Figure 8. These models show the strongest degree of correlation across the density and magnetization models, giving a more definitive picture of the Uranus diapir. We can see, indeed, that the structural overhangs necessary for a hydrocarbon trap are absent despite a large salt volume. The poor candidacy for hydrocarbon exploration was confirmed by drilling.
The probabilistic Gramian approach proved to be more numerically stable and efficient, and produced models with enhanced structural correlation.  Petrophysically guided inversion can be a useful tool in refining inversion results (Sun and Li, 2015); however, the a priori petrophysical relations may be unknown. One advantage of the Gramian approach, deterministic or probabilistic, is the oppotunity to see these petrophysical relations in a parametric cross plot of the inverted models. We present cross plots of density versus magnetization for each of the inversion methodologies in Figure 9.
The "cloud" of correlations shown by the blue dots in Figure 9 makes it difficult to determine any a priori relations in the standalone inversions. The calculated correlation coefficient for the standalone inversions in the area surrounding the Uranus diapir is 0.7. The cross plots for the deterministic and probabilistic Gramian inversions are shown by the green and red dots in Figure 9, respectively. Compared to the standalone inversions, a clear trend is now discernible in the scattered correlations. The calculated correlation coefficient for the deterministic Gramian is 0.73, which is a marginal improvement over the standalone inversions, and 0.93 for the probabilistic Gramian, which is a significant improvement over the other inversion methodologies.

Conclusion
Joint inversion of multiphysics data based on the Gramian constraints represents a powerful tool for integrated interpretation. The Gramian is defined as the determinant of the Gram matrix formed by models of different physical properties of the geologic formations, e.g., density, electrical conductivity, seismic velocity, magnetization, etc. This definition appears naturally in the framework of the deterministic approach to the inverse problem. At the same time, a significant body of research is based on the probabilistic approach to the inversion theory. In this case, the observed data, and the physical models are treated as the realizations of some random variables. The advantage of the probabilistic approach over the deterministic one is that the former provides some valuable statistics about the probability and accuracy of the physical models produced by the inversion.
We have developed and presented joint inversion methods based on the probabilistic approach, which provides insight into the properties of the inversion and helps improve the reliability of the inversion results. This was illustrated in the example from the Nordkapp Basin, where inverted models had a higher structural correlation and the numerical scheme was more efficient compared to the deterministic approach.
We have demonstrated how the joint Gramian-based inversion could be reformulated using the probabilistic approach. The Gramian is an analog of the determinant of the covariance matrix between the different physical properties representing the geologic formations. This helps understand better the role of the Gramian in enforcing the relationships between different physical models. It presents an alternative numerical implementation of the Gramiantype constraints by using the statistical estimates of the components of the covariance matrix.
The incorporation of the covariance in the Gramian constraint is an example of a Gaussian process, and was used to determine the relationship between the model parameters and data with no a priori knowledge. This opens the door to more modern probabilistic approaches, such as empirical Bayes, where model parameters can be estimated directly from the data in areas with little a priori knowledge. Such a prediction will, of course, depend on the amount and veracity of a priori knowledge or training models. Gaussian processes, Bayesian statistics, and empirical Bayes paired with rapidly advancing joint inversion strategies present a promising future for multi-modal geophysical inversion.