Simultaneous Structure-Coupled Joint Inversion of Gravity and Magnetic Data Based on a Damped Least-Squares Technique

The structure-coupled joint inversion method of gravity and magnetic data is a powerful tool for developing improved physical property models with high resolution and compatible features; however, the conventional procedure is inefficient due to the truncated singular values decomposition (SVD) process at each iteration. To improve the algorithm, a technique using damped least squares is adopted to calculate the structural term of model updates, instead of the truncated SVD. This produces structural coupled density and magnetization images with high efficiency. A so-called coupling factor is introduced to regulate the tuning of the desired final structural similarity level. Synthetic examples show that the joint inversion results are internally consistent and achieve higher resolution than separated. The acceptable runtime performance of the damped least squares technique used in joint inversion indicates that it is more suitable for practical use than the truncated SVD method.


Introduction
The joint inversion technique (i.e., imaging) is essential in interpreting collocated gravity and magnetic data.
However, handling the combination of geophysical data types, each of which is sensitive to different material properties, poses a major challenge.In the early stage, Menichetti and Guillen [1] proposed a generalized inverse method to recover 2.5-D models with fixed densities and susceptibilities.Serpa and Cook [2] inverted the densities, susceptibilities and the vertices of 2.5-D models with a priori information added.In recent years, some advanced numerical solutions were developed.Zeyen and Pous [3] were the first to study 3-D joint inversion of magnetic and gravity data, presenting a quasi-Newton algorithm for a designed model composed of vertical rectangular prisms with fixed lateral boundaries parallel to the local coordinate system.Gallardo-Delgado et al. [4] proposed a versatile algorithm for joint 3-D inversion of gravity and magnetic data, with the depths to the tops and bottoms as unknowns to be determined.Pilkington [5] jointly processed the gravity and magnetic data in terms of a model consisting of an interface separating two layers with constant density contrast and magnetization.Gallardo et al. [6] described and applied an automated refinement technique for 3-D multilayer models conditioned by gravity and magnetic data and by meaningful geometrical and physical constraints.The objective function they had used was a quadratic program in a stable iterative scheme, which produced models that correlated with the surface geology and revealed the subsurface features.Besides the layered model, a model including a large set of combined prisms is also preferred to approximately simulate the spatial variations of the physical properties.Fregoso and Gallardo [7] presented a cross-gradients joint framework applicable to gravity and magnetic data.The kernel technique and structural similarity criterion were first applied to electrical and seismic data [8], and were then soon developed for other types of data [9] [10].Zhdanov et al. [11] developed Gramian constraints for gravity and magnetic data to illustrate the effectiveness of the approach.Shamsipour et al. [12] presented a novel stochastic joint inversion method for gravity and magnetic data based on cokriging.To sum up, various methodologies of joint inversion for simultaneously handling these two kinds of data have been investigated, which is an indication of its importance and broad potential application.
In this study we present a joint inversion framework for gravity and magnetic data using a structural similarity criterion, which is proposed as an improved version of the Fregoso and Gallardo [7] approach when taking its practical application into consideration.First, the structural similarity criterion is reviewed in detail.Second, a cost function with cross-gradients component constraints is given, and iterative formula is presented based on the general nonlinear least-squares technique.The key issue in solving this problem concerns the structural vector updates at each iteration.Instead of conventional truncated singular value decomposition (SVD) adopted in previous work, we propose a damped least-squares algorithm to carry out controlled structure resemblance and efficient computation.We then test the method with several synthetic examples and compare the results obtained using mono-inversion of each dataset, and joint inversion of the integrated dataset, to demonstrate the reliability and efficiency of the proposed algorithm.

The Structural Similarity Criterion
The principle of the structural similarity criterion is that collocated multiple physical property models that reflect the same geological structure are considered implicitly to have common boundaries.Therefore the structural similarity criterion is required to quantitatively evaluate the structural resemblance of one or two physical property models.Haber and Oldenburg [13] introduced a structural operator for conducting constraint inversion processes.Gallardo and Meju [8] extended the formula to the cross-product of two given model gradients, termed the "cross gradient": , , , , , , where 1 m and 2 m are the two participating physical property values at an arbitrary subsurface location in 3-D space, and ∇ is the gradient operator.The model is a naturally continuous variable of x-, y-and z-geometrical measures; however, it is discretized with an aggregated prism set, which is suitable for carrying out numerical computations.The cross-gradients vector, τ , indicates the structural similarity at the current point.The vector can obviously be written as three orthometric components: ; .
which express the structural similarity of the two models on their orthometric vertical projection surfaces.When the model gradients have the same or reverse direction, that is, their gradients are parallel, these three components tend to have very small values; otherwise they would be non-zero quantities.So, for the whole model volume beneath the study area, a built-up array τ τ τ is used to measure the structural similarity throughout the whole space.When the structures are consistent, the array τ should approach a null space 0. This property is what we have used as a structural constraint in the cost function in order to obtain the desired results with a good structural resemblance feature.

Cost Function of Joint Inversion
As described by Gallardo and Meju [8], joint inversion promises to reduce the set of solutions by combining various data in a single scheme, and driving the updated model to explain all the requirements included in the cost function, such as data fitting, model smoothness, the distance of the model from prior information, and the structural similarity.Thus the cost function is expressed as: ) where subscripts 1 and 2 represent the gravity and magnetic methods; subscript 0 represents reference vector information; d is the observational data vector; p is the transformed model vector; ĝ is the transformed forward-modeling operator; the transform function between p and m is defined by the user to satisfy some known physical information; Z is the depth-weighting matrix [14]; L is the Laplacian matrix; C is the covariance matrix with respect to subscripts d (data), L (smoothness) and p (smallness), which are terms in the cost function; and τ is the assembled cross-gradients three-component array.The cost function satisfies the equality constraint, which is mainly where it differs from unconstrained optimization frameworks.The structural constraints provide the correlation between model parameters that formulate the complete algorithm.When the two models are coupled in structure, the operator tends to reach null space 0.
Following the robust statistical optimization algorithm framework developed by Fregoso and Gallardo [7], the nonlinear generalized inversion problem that minimizes this cost function is solved so that the data and discrete model can be submitted into a general framework and can be processed simultaneously: ( ) where p p is the combined model vector; N and n are a matrix and an array which involve the matrices described above; B is the cross-gradients Jacobian matrix with respect to 1 m and 2 m ; and k is the iteration number.The intermediate variables are expressed by: ( ) ( ) where J is the combined and transformed forward modeling Jacobian matrix.The inverse of N may be calculated by any conventional algorithm, and the inverse of  BN B , which does not have to be computed directly, is presented in detail below.The model update in the kth iteration is then obtained by: The final models 1 m and 2 m are extracted from p by an inverse transform function when the iteration procedure satisfies the given threshold concerned with data misfit and structural similarity.

Damped Least-Squares Technique
The inverse of matrix is usually difficult to compute because the production of the several involved matrices makes it singular.As described by Frogoso and Gallardo [7], the truncated SVD is an effective way of calculating its pseudo inverse.For the kth iteration, the decomposition formula is: ( ) where U and V are orthonormal matrices, and Λ is the corresponding singular value matrix.Then the reconstitution is expressed as: ( ) where c is the truncated index of the non-null singular sequence 1 2 c λ λ λ ≥ ≥ ≥ .However, this technique is not suitable for large-scale computations.Especially in the present case, the vectors are ready to compute in a combined manner, which enlarges the dimensions of both the model and the data.We suggest the use of an alternative algorithm to handle this problem.A damping factor is introduced to help improve the ill-posed of the matrix.The equation becomes: ( ) where I is the identity matrix, and β is the denominator of the damping factor.The second term added in brackets has the effect of making the total matrix well-posed.Although some structural perturbation information is lost in such a scheme, the computation stability is enhanced.This linear system is easy to compute with an iterative conjugate gradient solver, which accelerates the computing rates compared to truncated SVD.Then model updates of the equation at the kth iteration is solved using: which is the equation adopted in this paper.In this scheme, β is termed a coupling factor, and has an impact on the final achieved structural resemblance such that large values of β indicate greater model similarity.In other words, the structural coupling level may be readily controlled by tuning β by trial and error.A synthetic test below outlines the principle of choosing an optimal value of β.

Synthetic Examples
We first describe a synthetic example that illustrates the basic features of the proposed joint technique, using a typical synthetic model similar to that of Frogoso and Gallardo [7].The subsurface was divided into a mesh of 20 × 20 × 10 rectangular elements with a 50 m space in the x-, y-and z-directions.The geological target was located at a roof depth of 100 m beneath the center of the study area.The density contrast was 1.0 g/cm 3 and magnetization amplitude was 1.0 A/m. Figure 1 shows the theoretically modeled gravity and total magnetic intensity anomaly contaminated by 2.5% Gaussian noise, and the geomagnetic field was set as declination D = −5˚ and inclination I = 50˚.The separate inversion experiment was carried out by omitting the second term on the right-hand side of Equation (10), and the joint inversion was completed with 4 5.0 10 β = × ; the selection of this value is ex- plained below.For gravity and magnetic inversions, the depth-weighting technique should be included to counteract anomaly concentrations near the surface.The reference models for both density and magnetization model are set to 0. The maximum number of iterations was 6.
Figure 2 shows the inversion results for both the separate and joint inversion cases.For reasons of symmetry, only the vertical profile on northing = 500 m is presented.All the results reflect the true geological target in a fuzzy manner.Considered separately, better recovery was obtained from the magnetization model than from the density model due to the greater sensitivity of the magnetic method; however, the density values were lower, and these results were distorted with a broad tail at depth.In the joint inversion, both the density model and magnetization model improved, with more focused features than when separate.A significant characteristic was that the resultant density and magnetization structure were consistent.Not only did the resultant compatible images enhance the resolution, they also reduced the inherent non-uniqueness of the inverse problems.It is clear that the joint inversion technique incorporating the structural similarity constraint produced improved images with higher resolution.
To illustrate the use of the coupling factor β, we assigned different values to compute the cross-gradients L-2 norm at each iteration as shown in   structures are consistent with each other, while a large value accounts for the incompatible feature.Clearly, a large β value achieved better convergence than a small value; however, if too big, the stability of the computation may be compromised.The choice often entails a trial-and-error process.A relatively large number may be chosen provided it achieves satisfactory convergence of the cross-gradients norm and that the numerical computation is stable.
We then tested and compared the efficiency of the proposed damped least-squares technique and the conventional truncated SVD.Two different mesh discretization scales were adopted with various model dimensions to demonstrate runtime performance on the same computer platform.Table 1 shows that, for the smaller model scale, the damped method (final cross-gradients norm 8.3 × 10 −8 ) was slightly faster than SVD (final cross-gradients norm 3.9 × 10 −8 ), but for the larger model scale, the runtime of the truncated SVD (final cross-gradients norm 6.0 × 10 −8 ) dramatically increased at each iteration except the initial computation, compared to the damped method (final cross-gradients norm 2.4 × 10 −8 ).It was seen that the truncated SVD is not suitable for large-scale problems.The damped least-squares method obtains the same results but with better runtime performance.

Conclusions
Efficient, simultaneous 3-D structure-coupled joint inversion of collocated gravity and magnetic data using a damped least-squares technique is presented.It is consistent with the previous view that the cross-gradients constraints have a significant capacity to obtain more compatible models with higher structural resemblance than separated inversion.
The damped least-squares technique is introduced to perform efficient computation of the structural perturbation term in the linear system at each iteration.A coupling factor is proposed to regulate the final structural similarity level achieved.A relatively large value of this factor, chosen by trial and error, is suggested to achieve desired cross-gradients norm convergence on condition that the numerical computation is stable.
The synthetic experiments show that the joint inverted density and magnetization results are consistent, which reveals the geological setting.Furthermore, the runtime performance is more acceptable than that of truncated SVD, indicating its suitability for practical use.

Figure 3 .Northing 5 Figure 1 .Figure 2 .
Figure 1.Schematic of synthetic gravity data (a) and magnetic data (b).The black boxes indicate projected view of geological targets with different physical property feature.

Table 1 .
Runtime comparison of damped least squares and truncated SVD.Program) of China (No. 2013AA063901-4 and 2013AA063905-4), R&D of Key Instruments and Technologies for Deep Resources Prospecting (The National R&D Projects for Key Scientific Instruments) (No. ZDYZ2012-1-02-04), Constrained multi-parameter inversion of geophysical technology and software systems (No. 2014AA06A613), and The Fundamental Research Funds for the Central Universities.