Next Article in Journal
Rodingites in the Darbut Ophiolitic Mélange, West Junggar: New Insights into Rodingitization and Tectonic Evolution
Next Article in Special Issue
Editorial for the Special Issue “Electromagnetic Exploration: Theory, Methods and Applications”
Previous Article in Journal
Selenium Dissolution from Decopperized Anode Slimes in ClO/OH Media
Previous Article in Special Issue
Magnetotelluric Noise Attenuation Using a Deep Residual Shrinkage Network
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Magnetotelluric Regularized Inversion Based on the Multiplier Method

1
School of Geosciences and Info-Physics, Central South University, Changsha 410083, China
2
Key Laboratory of Metallogenic Prediction of Nonferrous Metals, Ministry of Education, Changsha 410083, China
3
School of Civil Engineering, Central South University, Changsha 410083, China
*
Author to whom correspondence should be addressed.
Minerals 2022, 12(10), 1230; https://doi.org/10.3390/min12101230
Submission received: 4 July 2022 / Revised: 21 September 2022 / Accepted: 24 September 2022 / Published: 28 September 2022
(This article belongs to the Special Issue Electromagnetic Exploration: Theory, Methods and Applications)

Abstract

:
Magnetotellurics (MT) is an important geophysical method for resource exploration and mineral evaluation. As a direct and effective form of data interpretation, MT inversion is usually considered to be a penalty-function constraint-based optimization strategy. However, conventional MT inversion involves a large number of calculations in penalty terms and causes difficulties in selecting exact regularization factors. For this reason, we propose a multiplier-based MT inversion scheme, which is implemented by introducing the incremental Lagrangian function. In this case, it can avoid the exact solution of the primal-dual subproblem in the penalty function and further reduce the sensitivity of the regularization factors, thus achieving the goal of improving the convergence efficiency and accelerating the optimization calculation of the inverse algorithm. In this study, two models were used to verify the performance of the multiplier method in the regularized MT inversion. The first experiment, with an undulating two-layer model of metal ore, verified that the multiplier method could effectively avoid the MT inversion falling into local minimal. The second experiment, with a wedge model, showed that the multiplier method has strong robustness, due to which it can expand the selection range and reduce the difficulty of the regularization factors. We tested the feasibility of the multiplier method in field data. We compared the results of the multiplier method with those of conventional inversion methods in order to verify the accuracy of the multiplier method.

1. Introduction

Magnetotellurics (MT) is a geophysical exploration method that uses natural planar electromagnetic waves as field sources and obtains geoelectric tectonic information by observing mutually orthogonal electric and magnetic field components at the surface [1,2]. MT inversion is the reduction in the electrical distribution in the subsurface through MT responses, such as apparent resistivity, impedance phase, surface impedance, and tipper, and it is a very important part of MT data processing and interpretation [3].
In the last half-century, MT inversion has developed rapidly in terms of global optimization and local optimization. The global optimization methods include simulated annealing algorithm [4], neural network inversion [5,6], genetic algorithm [7], particle swarm optimization [8], Bayesian method [9,10], and deep learning [11,12]. The local optimization methods include Occam inversion [13,14], simplified-based Occam inversion [15], rapid relaxation inversion [16,17], nonlinear conjugate gradient inversion [18,19,20], and Gauss–Newton inversion [21,22,23,24,25]. The inversion method of global optimization prevents the inversion process from falling into the local minimal through global parameter search. However, it increases the time and storage cost of inversion optimization, especially for high-dimensional problems, and there is still significant scope for improvement in this regard [26]. The inversion method of local optimization offers better speed, but the inversion results easily fall into local minimal, and it is difficult to capture a satisfactory global optimal solution [27].
To resolve the problem that local optimization inversion methods tend to fall into local minimal, researchers have used different methods. Newman and Alumbaugh argue that trying multiple initial models can prevent the problem of inversions falling into local minimal [28]. Adaptive weighted averaging in the Niblett–Bostick transform preserves part of the nonlinearity of the original problem, and this technique weakens the dependence of the inverse model on the initial model [29]. To obtain a reliable initial model, an imprinting algorithm based on the weighting method to construct the initial model was proposed [30]. Numerous studies have shown that regularization methods not only deal with the discomfort problem but also positively affect the problem of falling into local extrema for local inversion [31,32,33,34,35]. In the current MT regularization inversion, the main methods commonly used for the selection of the regularization factor are the L-curve method [36], the generalized cross-validation method [37], the adaptive algorithm [38], the MT adaptive regularization [39], and others. Nevertheless, no common method has been found for the selection of the regularization factor.
At present, the multiplier method is a widely used algorithm in solving optimization problems with multiple constraints [40,41,42,43,44]. Groot-Hedlin and Constable [45] developed a linearization algorithm to invert the noisy MT data of the subsurface conductivity structure represented by the sharp resistivity contrast defined by the smooth boundary and proposed a new balance mechanism of model roughness and data-fitting difference based on the Lagrange multiplier method. This mechanism provides a better solution for regularization factor selection, but it does not consider the effect of the multiplier method on the regularization factor. The multiplier method is less commonly used in geophysical inversions [46]. In this paper, the modified total variation is combined with the multiplier method for the first time and applied to the MT inversion. The multiplier method, also known as the Augmented Lagrange Method (ALM), starts from the essence of the problem and transforms the constraints of the constrained problem into one of the objective functions by introducing a strategy of Lagrange multipliers, completing the conversion from a constrained optimization problem to an unconstrained optimization problem. The improvement of the multiplier method on the data-objective function prevents the MT inversion from falling into local extremes and guarantees the stability of the inversion, while the approximate estimation of the solution reduces the dependence of the optimization process on the regularization factor.
In this paper, the multiplier method is applied to the study of MT inversion for the first time. The multiplier method is used to improve the convergence speed of the inversion by preventing the inversion from falling into local extremes to improve its efficiency. In addition, due to the unique method of solving the approximation estimation of the multiplier method, for the first time, the multiplier method is combined with the MT regularization translation to reduce the dependence of the inversion on the regularization term, improve the threshold range of the regularization factor, and increase the regularization range of the selection of factors. In this study, an undulating two-layer model was used to verify that the multiplier method has better convergence performance in MT inversion. In addition, a wedge model was used to verify that the multiplier method has stronger robustness in MT regularized inversion.

2. Methodology

2.1. Conventional Inversion

The MT inversion problem can be expressed as follows.
d = F ( m ) + e
d denotes the observed data vector, F(m) is the forward response operator, m denotes the model parameter vector, and e is the residual vector. Since the inversion problem has multiple solutions, simply speaking, an exceptionally wide range of solutions can meet the fitting requirements of the observed data. Therefore, to obtain the required inversion results, we need to rely on other a priori information to constrain the inversion. The penalty-function method is a classical algorithm for solving constrained problems. The fundamental idea is that the constraint conditions are transformed into some kind of penalty function added to the data-objective function according to their characteristics, thus transforming the optimization problem with constraints into a series of unconstrained optimization problems to be solved. Thus, the optimization problem can be written as
min Φ d ( m ) s . t .   R ( m )
where Φd(m), R(m) and m represent the data-objective function, regularization term, and inversion parameter, respectively. In mathematical form, Equation (2) can also represent the constraint of the data-objective function on the regularization term as
min R ( m ) s . t .   Φ d ( m )  
The conventional MT inversion problem is an optimization problem with constraints, which can be expressed as [47].
Φ ( m ) = Φ d ( m ) + λ R ( m )
of which
Φ d ( m ) = W d ( F ( m ) d ) 2
where Wd in Equation (5) is the data-weighting matrix, which is typically based on the inverse of the measured standard deviation, F(m) is the forward response operator, m denotes the model parameter vector, and d denotes the observed data vector.
The MT inversion problem with constraints is converted into an unconstrained optimization problem by the external penalty function method, which can be expressed as
min R ( m ) s . t .   W d ( F ( m ) d ) 2
solved using the external penalty Equation (2). The process can be expressed by Algorithm 1 as
Algorithm 1 External Penalty Method
Given β > 0, tolerance ɛ0 > 0, starting points m0 and λ0;
for k = 0, 1, 2, …
Find an approximate minimizer mk of Φ(m), starting at mk, and terminating when m Φ ( m ) ε k ;
if final the convergence test satisfied
stop with approximate solution mk;
end if
Choose new penalty parameter λk+1λk;
Choose a new starting point mk+1;
end for
The external-penalty-function method has the feature of a simple solution framework, which is solved by directly calling the general procedure of the unconstrained optimization algorithm, and the procedure is simple to implement. However, the disadvantages are more obvious: (1), for complex geophysical problems, mk is not a feasible point of the problem itself; (2), it is very difficult to obtain an appropriate penalty parameter λk. When it is too small, the penalty term does not exert its proper effect. Conversely, when the penalty parameter is too large, the number of conditions for solving the penalty function also increases, and the equation becomes more difficult to solve; and (3), the most critical point, since the penalty function is generally non-integrable, it is difficult to directly use the derivative of the optimization algorithm, resulting in the slow convergence of the algorithm.

2.2. The Multiplier Method

The multiplier method, also known as the Augmented Lagrange Method (ALM), is based on the Lagrange function of the original problem to start with an appropriate penalty function, thus transforming the original problem into the solution to a series of unconstrained optimization subproblems. The Lagrangian function of the original problem, however, can be obtained from the Lagrange Multiplier Method (LMM), a method for finding the extrema of multivariate functions whose variables are constrained by one or more conditions, which introduces a new scalar unknown, the Lagrange multiplier: the gradient of the constrained equation of coefficients of each vector in the linear combination. The Lagrangian function for the optimization problem of Equation (3) can be written as
L ( m , μ ) = R ( m ) + μ T ( F ( m ) d )
where μ is the Lagrangian multiplier vector. When L(m,μ) = 0, we can obtain the optimal value by calculation. For the most part, the gradient method cannot be used to optimize the Lagrangian function directly because the optimal value of the original problem corresponds to the saddle point of the Lagrangian, and the gradient method cannot guarantee that the saddle point is found.
The multiplier method, also known as the augmented Lagrange method, is shown in Algorithm 2. Moreover, if the gradient of the primary problem ▽Rmk) is a higher-order function, the primary problem can only be solved by an iterative method, such as the Newton method. However, the iterative method requires that the primary problem function satisfies the second-order derivative, and the second-order derivative-solution process is relatively expensive. Furthermore, due to the penalty parameter μ→+∞ in the external penalty function method in ALM, the pathological nature of the augmented objective function constructed in the external penalty function method is enhanced. This pathological nature of the augmented objective function is the main drawback of the external penalty function method, and this drawback is effectively overcome in the multiplier method by the introduction of the Lagrangian function and the addition of an appropriate penalty function. To overcome the drawbacks of the Lagrangian method, the multiplier method is proposed, and it is used to calculate Equation (7), which can be written as [48]:
Algorithm 2 Augmented Lagrangian Method
Given β > 0, tolerance ɛ0 > 0, starting points m0 and μ0;
for k = 0, 1, 2, …
Find an approximate minimizer mk of L A ( m k , μ k ) , starting at mk, and terminating when m L A ( m k , μ k ) ε k ;
if the final convergence test satisfied
stop with approximate solution mk;
end if
Update Lagrange multipliers using (12) to obtain μk+1;
Choose new penalty parameter μk+1 ≥ μk;
Set starting point for the next iteration to mk+1 = mk;
Select tolerance εk+1;
end for
L A ( m , μ ) = R ( m ) + μ T ( F ( m ) d ) + λ F ( m ) d 2 2
where λ > 0 is the penalty parameter. For the constraint problem, this can be written as
min m max μ , λ R ( m ) μ , F ( m ) d λ , F ( m ) d + μ T ( Ax b ) + λ F ( m ) d 2 2
where < > denotes the inner product. Formula (9) can be abbreviated as
min m max μ , λ R ( m ) μ , F ( m ) d + λ F ( m ) d 2 2
Hence, the updated formula for the multiplicative method can be written as:
m k + 1 = arg min   m R ( m ) μ k , F ( m ) d + λ F ( m ) d 2 2
μ k + 1 = μ k λ ( F ( m ) d )
By integrating Equation (8) into Equations (11) and (12), it is not difficult to find that the value range of λ is very wide at this time. Through mathematical analysis of Equations (11) and (12), we find that the product of the weight factor λ and the square in Equation (8) becomes the product of the weight factor λ and the first-order data. This method reduces the weight of the weight factor λ, which can reduce the difficulty of selecting the weight factor in the Lagrangian method.
To explain Equation (12) more thoroughly, we use Figure 1 to represent the process of data correction, in which “a”, “b”, “c”, and “d” denote the extreme values of the problem, and assuming that “d” is the global minimum value, the others are the local extreme values. In addition, the solid black line indicates the original state of the optimization problem, i.e., the real scenario; the short blue line indicates the case in which the conventional optimization method falls into the local extrema; and the red dashed line indicates that through the multiplicative method proposed in this paper, we modify the overall optimization equation so that we do not fall into the local extrema when performing optimization, thus achieving the dual effects of accelerating convergence and the avoidance of falling into local extrema.

3. Magnetotelluric Inversion Based on Multiplier Method

3.1. Data-Objective Function

In this part, we propose a new iterative algorithm for solving the MT inversion problem of Equation (4) based on the multiplicative method. We use the multiplicative method to modify the data-objective function of the optimization process so that the optimization problem is modified under the constrained condition, the constrained optimization problem is transformed into an unconstrained problem, and multi-scale inversion is achieved. At this point, the MT inversion objective function can be written as
min m max μ   W d ( F ( m ) d * ) 2 2 μ , W d ( F ( m ) d * ) + λ R ( m )
μ is the Lagrange multiplier vector and d* is the observed data. Equation (4) is the multiplier method used to iteratively solve this minimal problem. In each iteration of the algorithm, the objective function is first minimized for m, with μ held constant, after which the multiplier μ is updated to maximize the objective. The major advantage of such a solution procedure, as mentioned in Hestenes and Powell’s paper, is that it features a very simple update iteration formula [48].
m k + 1 = arg min m   W d ( F ( m ) d * ) 2 2 μ k , W d ( F ( m ) d * ) + λ R ( m )
μ k + 1 = μ k + λ W d ( F ( m k + 1 ) d * )
where k is the number of iterations and the initial Lagrange multiplier µ0 = 0. Taking the partial derivatives of the model parameters m for (14) results in the following equation:
J ( m ) T [ 2 β W d ( F ( m ) d * ) μ k ] + λ m R ( m )
in which
J ( m ) = F ( m ) m
J is denoted as a Jacobi matrix, and the superscript T denotes the transpose. In Equation (16), the rightmost term in the parentheses, μk, can be considered as the Lagrange multiplier vector, which partially satisfies the optimality conditions of the original problem equation. It follows that a reasonable choice for the multiplier vector for the next iteration should be Equation (15). Furthermore, if we use the variable projection method to reconstruct the extended Lagrangian function in an equivalent dyadic form, in which only the multiplier vector is involved, then the gradient of the dyadic function is 2βWd (F(m) − d) and, hence, the optimal multiplier satisfies F(m) = d [49]. Therefore, the update direction in Equation (15) is considered as the gradient direction (the step length of the gradient is λ) for solving the pairwise problem so that it is locally maximized concerning μ. With this in consideration, we can apply a quasi-Newton approach to this pairwise problem to improve convergence, but at the cost of increased computational complexity.

3.2. Model Updates

To obtain the model update, this paper presents a solution that uses a simplified approach. For Equations (18) and (19), assuming that the variable dk varies with µ/λk and ignoring the constant term associated with µ, Equation (14) is simplified [46], which gives the data-objective function update equation as
m k = arg min m   W d ( F ( m k 1 ) d k 1 ) 2 + λ R ( m k 1 )
d k = d k 1 ( F ( m k 1 ) d * )
where d0 = d*, m0 is the initial inversion model which can be chosen according to the geological background. This paper adopts the 100 Ω∙m uniform model. It is easy to see that only one step in the above iterative method is needed to approximate the minimization equation before updating the multipliers. Proof of the convergence of such algorithms for inexact solutions of the primitive subproblem was provided by Miele [50]. In addition, this method has also been described by Tapia [51], who calls it the diagonal variable method of the conventional method.
Considering that the modified total-variation regularization has the dual effect of ensuring the inversion stability and accurately reproducing the sharp interface [23], the modified total-variation regularization term is adopted in this paper.

4. Numerical Examples

In this part, the performance of the multiplier method in MT inversion is verified by the following numerical algorithms. Aiming to verify the performance of the multiplier method in avoiding local minimal, this study compared the convergence of the multiplier method and the conventional inversion method using an undulating two-layer model. The comparison experiments were carried out in terms of inversion-result plots, convergence-curve plots, inversion-related index values, and the gradient of the inversion root mean square (RMS). The hardware platform used for this inversion study was a personal computer with a basic configuration of Intel(R) Core (TM) i7-7700 CPU @ 3.60 GHz (Interl, Santa Clara, CA, USA).
In this inversion study, for the forward simulation part, we used a non-regular quadrilateral to dissect the grid of the target body and then adopted the finite-cell method for the simulation calculation [52]. To reduce the influence of other factors on this study, we followed the principle of a single variable, and the MT inversion method was the joint method of TE + TM. To make full use of the existing data, we chose the apparent resistivity and phase-synchronization-inversion approach, in which we do not present the inversion results of phase in the numerical model, in order to ensure the compactness of the article, and we present the results of both the apparent resistivity and the phase uniformly in the field data.

4.1. Undulating Two-Layer Model

To verify the performance of the multiplier method in MT inversion, we first constructed a two-layer model in which the first (top) layer has a resistivity of 100 Ω∙m while the second layer has a resistivity of 1000 Ω∙m. Since the electrical properties of metal deposits are generally of low resistance, we designed the model as in Figure 2a. As shown in Figure 2a, we insert two low-resistivity (10 Ω∙m) ore bodies into the top layer. The two anomalies have different shapes. To make the topography realistic, we set the ground surface in an irregular shape. The center of the left quadrilateral low-resistance body was buried at a depth of 1 km with a lateral length of 1.5 km, and the center of the right elliptical low-resistance body was buried at a depth of 1.2 km with a lateral length of 4 km. The total length of the study area was 20 km, with 41 observation points, and the observation interval was 0.5km; the frequency range was 0.01~1000 (Hz), with 10 frequency points: 0.01, 0.0359, 0.1292, 0.4642, 1.6681, 5.9948, 21.5443, 77.4264, 278.2559, and 1000 (Hz). In this study, the finite-element method was used for the forward simulation, and different forward and inverse grids were used during the study. The forward grid was 208 × 72 (the air region was set to 10 layers), and the inverse grid was 85 × 45.
In this paper, we used the control-variables method to investigate the performance of the multiplier method in MT inversion. The initial model was a homogeneous half-space with a resistivity of 100 Ω∙m in the first layer; the conventional total-variance regularization constraint was used, and the optimal results of the two inversions were selected empirically. The modified total-variation regularization was used in this experiment for the following reasons: (1) the modified total-variation regularization is a better constraint on the inversion effect of the model with interfaces, and the performance of the multiplier method in the MT inversion application can be judged more intuitively; and (2) the total-variance regularization has a better stabilization constraint on the inversion. The inversion results of the conventional method and the multiplier method are shown in Figure 2b,c, respectively.
The black dashed lines in Figure 2b,c indicate the locations of anomalies in the orthorectified model. Firstly, comparing the inversion recovery integrity, the inversion results of both methods completely recovered the two low-resistance bodies and the high-resistance laminae, and no key part was missing from the inversion results; secondly, in terms of the inversion accuracy of the two methods, the multiplier method was more accurate than the conventional inversion method. Compared with the conventional method, the multiplier method can clearly inverse the stratigraphic partition interface and better suppress the high-resistance artifacts next to the low-resistance boundary, all these improved support for the interpretation of the field data.
Figure 3 shows the convergence curves of the inversions of both the conventional MT inversion and the multiplier-method MT inversion. The black dashed line indicates the convergence curve of the conventional inversion method, and the red dashed line indicates the convergence curve of the multiplier-method inversion. There are two convergence curves for the multiplier method inversion; the first appears at the 12th inversion iteration, and the second appears at the 22nd inversion iteration. These two iterations are described in detail below. The comparison graph of the inversion convergence curves demonstrates that the multiplier method can accelerate the inversion convergence speed and improve the inversion efficiency when the inversion falls into slow convergence, while the conventional inversion method makes the results satisfy the conditions set in advance by the smooth convergence in the longer convergence and repeated iterations, which increases the cost of the inversion. The model error, RMS error of the simulated data fitting, and inversion time of the two methods are compared in Table 1. The model error is calculated in Equation (20), where minv is the model parameter calculated by the inversion method, mtrue is the true model parameter, and M is the number of model parameters. Equation (20) represents the error between the model calculated by inversion and the true model. The RMS error is calculated in Equation (21), where dfor is the calculated forward response, dobs is the observed data, and N is the number of data. Table 1 shows the comparison of the inversion-evaluation data of the two methods; the multiplier method is better than the conventional inversion method in terms of both accuracy and speed.
E r r o r = 1 M M ( | m i n v m t r u e | | m t r u e | × 100 )
R M S = 1 N N ( d f o r d o b s d o b s ) 2
This experiment was implemented in two steps. Firstly, when the gradient of the inverse convergence curve was less than 5% for three consecutive times, the first multiplicative method was implemented to induce the data-objective function to adjust to the global extremes under the condition of mathematical constraints; the threshold condition for the second implementation of the multiplicative method was that the gradient of the convergence curve was less than 2% for three consecutive times; the implementation effect is shown in Figure 4. CMTI is the conventional magnetotelluric inversion and MMTI is the multipliers magnetotelluric inversion. The solid black line is the gradient of the RMS data fitting error curve of the conventional MT inversion. The solid red line is the gradient of the RMS curve based on multiplicative MT. The red dashed line is the threshold line of 5%. The black dashed line is the threshold line of 2%. The black solid lines “ο” in Figure 4 indicate the number of iterations needed to achieve the conditions required to implement the multiplicative method in this experiment, which are: 9, 10, 11, 19, 20, and 21. The first implementation of the multiplication method is at the 11th iteration, and the second inversion iteration is at the 21st iteration. The solid black lines in the figure form the first point after the triangle with the red dashed line and the black dashed-dotted line, respectively.
After this experiment, the results showed that the multiplier method could effectively improve the inversion efficiency in the application of MT inversion, mainly in terms of accelerating the inversion convergence and improving the inversion accuracy.

4.2. Wedge Model

In this study, the multiplier method was applied to the MT regularized inversion in response to the difficulty of regularization-factor selection in the regularized inversion. Given the advantages of the multiplier method in the regularization-factor selection, the experimental design in this example aims to investigate the performance of the multiplier method in the MT regularized inversion and to verify the effectiveness of the method through a wedge model.
The experimental model uses a wedge model, as shown in Figure 5, in which the slope of the upper boundary of the wedge is 5.7° and the slope of the lower boundary is 16.6°. The details of the model are shown in Figure 5. The red part of the resistivity is a 1000-ohm-meter high-resistance block, the blue area is a 10-ohm-meter low-resistance wedge, the green polygon indicates the background with a resistivity of 100 Ω∙m, and the color scale on the right-hand side in Figure 5 of the model indicates the range of the resistivity. A total of 22 observations were set up for this model, each at a spacing of 0.5 km. The grid profile in the forward simulation was 132 × 60, with 10 frequency points and frequencies of 0.01, 0.0359, 0.1292, 0.4642, 1.6681, 5.9948, 21.5443, 77.4264, 278.2559, 1000 (Hz), and 10 grid layers in the air layer. Different size grids were used for the forward and inverse, and the inverse network was 66 × 30.
To determine the dependence of the two methods on the regularization factors, this comparison experiment was conducted by selecting four different sets of regularization factors, all within the range required to ensure that the regularization factors were valid for the inversion, i.e., the inversion results were within the corresponding error range, in which the four sets of regularization factors were selected from different scales of variation, namely: 0.05, 0.1, 0.3, and 0.5. The inversion results are shown in Figure 6; Figure 6a–d show the results of the conventional inversion, and Figure 6e–h show the results of the multiplier-method inversion. The results suggest that the proposed method can reconstruct the underground structure more accurately than the conventional method. From the perspective of the overall structure recovery of the model and the reconstruction of the resistivity value of each abnormal body, the proposed method has a good effect. It shows that the accuracy of the proposed method is not only guaranteed but also that the details can be reconstructed better. To demonstrate its advantages in the selection of regularization factors, the new method is discussed in detail below.
In this study, the inversion results are compared in detail. The results of the conventional methods have two characteristics. Firstly, the interface of the low-resistance body is not well recovered, and the inversion results have obvious downward trend, which makes the data analysis difficult. Secondly, in the black boxes of Figure 6a,c,d, the background resistivity of the conventional inversion results is higher than that of the real model. Unlike conventional inversion methods, the multiplier method is modified for the data-objective function. Among them, there are two types of corrections: a predefined fixed number of times and a number of times intelligently controlled by setting a threshold value for the data-fit difference. The optimization was carried out through multiple steps and directions to prevent the inversion process from falling into local extremes and enhances the inversion’s detailed portrayal of the target region. Figure 7 shows the convergence curves of the inversion for different methods and different regularization factors. Compared with the conventional method, the convergence curve of the multiplier method has obvious characteristics. The convergence trend of the multiplier-method curve is very similar to that of the multi-scale method. The multiplier method performs optimization calculations in one scale and enters the next scale domain when the convergence condition is reached. This method can accelerate the convergence of optimization calculation and thus improve the inversion efficiency. This indicates that the multiplier method outperforms the conventional inversion method in MT inversion.
Figure 8 shows the single-point data comparison through the dashed line in Figure 5 of the wedge model, i.e., the data of the observation point at the horizontal 5.5 km. The dashed line indicates the inversion result of the conventional inversion method, and the solid line indicates the inversion result of the multiplier method. The dashed line in Figure 8 has two characteristics: (1) the recovery deviation of the high-resistance body is large, and the result shows linear variation, which, in this case, can be interpreted as a manifestation of inversion multi-solvability, where a single layer of the high-resistance body, recovered into a linearly varying high-resistance body, causes the misjudgment of subsequent interpretations of the data and increases the difficulty of interpreting the data information; (2), the boundary differentiation between different anomalies is insufficient to locate the anomalies and judge their scope in a timely and effective manner. However, from the comparison in Figure 8, it is easy to find that the solid line results in the figure define the target boundary well and prevent the problems discussed above.
The most critical aim in designing this experiment was to demonstrate whether there was a difference between the effect of the sensitivities of the conventional method and the multiplier method on the regularization factor. Therefore, the inversion parameters of the two methods were compared, and the results are presented in Table 2, where the inversion model error, the difference in the fit of the data, and the inversion time are presented, respectively. The model error represents the ratio of the difference between the inversion results and the forward model, as detailed in Equation (20). The data-fit difference represents the ratio of the difference between the last positive result and the observed data in the inversion to the observed data, as detailed in Equation (21). The inversion time represents the time spent on the inversion process for the same number of inversion iterations for different regularization factors. As can be seen in Table 2, the difference in time taken by different regularization factors is not significant; however, there is a large difference in their model error and data fit. For this reason, in this paper, the model error and data fit difference are presented in a more intuitive bar chart and a line graph, with the blue block and blue line in Figure 9 indicating the model error and data-fit difference of the inversion results of the conventional method, respectively, and the red block and red line indicating the model error and data-fit error of the inversion results of the multiplier method. The figure shows that the multiplier method has better results than the conventional method, and in the inversion results of the conventional method, the influence of different regularization factors on the inversion results fluctuates greatly. For example, when λ = 0.1 is the best regularization factor for this experiment, when λ is greater or less than 0.1, the model error and data-fit difference are greater than when λ = 0.1. Combining the inversion results in Figure 6, it was easy to find that these two methods have different performances under different regularization factors. In the conventional method, there were large differences in the inversion results, model errors, and data fit between different regularization factors; however, the multiplier method obtained better inversion results, model errors, and data-fit differences for different regularization factors. This experimental procedure demonstrates that the multiplier method has strong robustness and can obtain better inversion results for different regularization factors, which can reduce the difficulty of regularization factor selection.
In summary, the multiplier method, with its excellent robustness, can obtain more reliable inversion results among different regularization factors. Therefore, this approach can reduce the workload involved in selecting the appropriate regularization factor and greatly improve the efficiency of accurate inversion. Furthermore, the inversion iteration time is not significantly different from that of the conventional method, which reflects the advantages of the multiplier method in reducing the difficulty involved in selecting the regularization factor in many respects.

5. Field-Data Example

The numerical examples demonstrate that the proposed algorithm is feasible and effective. Therefore, we verified its efficiency through field data. We chose the data in an open-source program for verification and compared them with the inversion results of the conventional algorithm. The data were obtained from the observation of the Lachlan folded belt in central Australia, with 53 stations in total, and the measured line trend was East–West, with a total length of 180 km. The number of field collection frequencies was 27. The frequencies were: 0.0092, 0.0134, 0.0183, 0.0269, 0.037, 0.054, 0.073, 0.107, 0.146, 0.215, 0.293, 0.43, 0.59, 0.86, 1.17, 1.72, 2.34, 3.4, 4.7, 6.9, 9.4, 13.7, 18.8, 27.5, 40, 66, and 97 (Hz). The remote reference point was 160 km away from the midpoint of the survey line.
For these field data, we compared the MT inversion results of the conventional method and the multiplier method. Specifically, the inversion results of the field data, the final RMS data fitting errors, and the iteration number of the inversion were compared. We selected an initial model with a resistivity of 100 Ω∙m and carried out 30 iterations of optimization. The multiplier method was used to correct the data-objective function. In this study, we used a fixed number of iterations to modify the data-objective function. Through experiments, we found that the use of five iterations was the most appropriate.
Figure 10 shows the inversion results of the two methods. The two results were very close to each other, and we consider our inversion results to be credible. In Table 3 and Figure 11, it is shown that the RMS data fitting error of both methods is relatively small. The RMS data fitting error of the conventional inversion method is 0.0760, while the RMS data fitting error of the multiplier method is 0.0408. Therefore, we suggest that the inversion result of the multiplier method is closer to the real situation. Finally, there are subtle differences between our proposed method and the conventional method. For example, our proposed method presents a separate small anomaly in the results, but the location of this block can be found in the conventional inversion results, as shown by the black rectangular box in Figure 10. Furthermore, combined with the RMS data fitting error values in Table 3, the proposed method has a smaller data-fit difference, and we can conclude that the result of the proposed method better reflects the real situation.

6. Conclusions

(1) In this paper, the modified total variation was combined with the multiplier method for the first time and successfully applied to the MT inversion. To address the problem of the difficult regularization selection and significant workload in inversion, the effect of the multiplier method on the selection of regularization factors was comprehensively discussed in this paper. In general, the multiplicative method can optimize the selection of regularization factors for inversion.
(2) By comparing with the open-source program, it was proven that the application of the method proposed in this paper to the inversion of field data is feasible. Through the research, it was found that the method proposed in this paper can complete not only the inversion of the measured data but also obtain better results than conventional inversion.
(3) As shown in this paper, the multiplier method can accomplish the task of preventing the inversion from falling into local extremes and reducing the difficulty of selecting the regularization factor in MT inversion. Therefore, similar problems can be solved by applying the multiplier method for other complex geophysical inversion problems.
(4) This paper deals with single-parameter inversion, and the multiplier method can perform this task relatively satisfactorily. However, when faced with two-parameter or multi-parameter inversion, such as when the method is applied to the two-parameter inversion of the radio magnetotelluric method, the effectiveness of the method will require further study. Therefore, two-parameter or multi-parameter multiplier inversion will be a future research direction.

Author Contributions

Conceptualization, D.F.; data curation, X.S. and S.D.; formal analysis, X.S.; funding acquisition, D.F. and X.W.; investigation, S.L.; methodology, X.S. and X.W.; project administration, D.F.; resources, D.F. and X.W.; software, X.S.; supervision, X.W.; validation, X.S. and Y.L.; visualization, C.C.; writing—original draft, D.F.; writing—review and editing, S.L. and C.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded in part by the National Natural Science Foundation of China (grant numbers 42074161 and 42104143) and in part by the Natural Science Foundation of Hunan Province, China (grant number 2021JJ0806).

Data Availability Statement

Data associated with this research are available and can be obtained by contacting the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Cagniard, L. Basic theory of the magneto-telluric method of geophysical prospecting. Geophysics 1953, 18, 605–635. [Google Scholar] [CrossRef]
  2. Tikhonov, A. On determining electrical characteristics of the deep layers of the Earth’s crust. Magn. Methods 1986, 73, 295–297. [Google Scholar]
  3. Jupp, D.L.; Vozoff, K. Two-dimensional magnetotelluric inversion. Geophys. J. Int. 1977, 50, 333–352. [Google Scholar] [CrossRef]
  4. Dittmer, J.K.; Szymanski, J. The stochastic inversion of magnetics and resistivity data using the simulated annealing algorithm. Geophys. Prospect. 1995, 43, 397–416. [Google Scholar] [CrossRef]
  5. Spichak, V.; Popova, I. Artificial neural network inversion of magnetotelluric data in terms of three-dimensional earth macroparameters. Geophys. J. Int. 2000, 142, 15–26. [Google Scholar] [CrossRef]
  6. Wang, H.; Liu, M.; Xi, Z.; Peng, X.; He, H. Magnetotelluric inversion based on BP neural network optimized by genetic algorithm. Chin. J. Geophys. 2018, 61, 1563–1575. [Google Scholar] [CrossRef]
  7. Pérez-Flores, M.A.; Schultz, A. Application of 2-D inversion with genetic algorithms to magnetotelluric data from geothermal areas. Earth Planets Space 2002, 54, 607–616. [Google Scholar] [CrossRef]
  8. Shaw, R.; Srivastava, S. Particle swarm optimization: A new tool to invert geophysical data. Geophysics 2007, 72, F75–F83. [Google Scholar] [CrossRef]
  9. Zhou, S.; Huang, Q. Two-dimensional sharp boundary magnetotelluric inversion using Bayesian theory. Chin. J. Geophys. 2018, 61, 3420–3434. [Google Scholar]
  10. Xiang, E.; Guo, R.; Dosso, S.E.; Liu, J.; Dong, H.; Ren, Z. Efficient hierarchical trans-dimensional Bayesian inversion of magnetotelluric data. Geophys. J. Int. 2018, 213, 1751–1767. [Google Scholar] [CrossRef]
  11. Conway, D.; Alexander, B.; King, M.; Heinson, G.; Kee, Y. Inverting magnetotelluric responses in a three-dimensional earth using fast forward approximations based on artificial neural networks. Comput. Geosci. 2019, 127, 44–52. [Google Scholar] [CrossRef]
  12. Liu, Z.; Chen, H.; Ren, Z.; Tang, J.; Xu, Z.; Chen, Y.; Liu, X. Deep learning audio magnetotellurics inversion using residual-based deep convolution neural network. J. Appl. Geophys. 2021, 188, 104309. [Google Scholar] [CrossRef]
  13. Constable, S.C.; Parker, R.L.; Constable, C.G. Occam’s inversion: A practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics 1987, 52, 289–300. [Google Scholar] [CrossRef]
  14. Degroot-Hedlin, C.; Constable, S. Occam’s inversion and the North American Central Plains electrical anomaly. J. Geomagn. Geoelectr. 1993, 45, 985–999. [Google Scholar] [CrossRef]
  15. Siripunvaraporn, W.; Egbert, G. An efficient data-subspace inversion method for 2-D magnetotelluric dataREBOCC Inversion for 2-D MT Data. Geophysics 2000, 65, 791–803. [Google Scholar] [CrossRef]
  16. Smith, J.T.; Booker, J.R. Rapid inversion of two-and three-dimensional magnetotelluric data. J. Geophys. Res. Solid Earth 1991, 96, 3905–3922. [Google Scholar] [CrossRef]
  17. Tan, H.; Yu, Q.; John, B.; Wei, W. Three-dimensional rapid relaxation inversion for the magnetotelluric method. Chin. J. Geophys. 2003, 46, 1218–1226. [Google Scholar] [CrossRef]
  18. Mackie, R.L.; Madden, T.R. Three-dimensional magnetotelluric inversion using conjugate gradients. Geophys. J. Int. 1993, 115, 215–229. [Google Scholar] [CrossRef]
  19. Rodi, W.; Mackie, R.L. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion. Geophysics 2001, 66, 174–187. [Google Scholar] [CrossRef]
  20. Hu, Z.Z.; Hu, X.Y.; He, Z.X. Pseudo-three-dimensional magnetotelluric inversion using nonlinear conjugate gradients. Chin. J. Geophys. 2006, 49, 1111–1120. [Google Scholar] [CrossRef]
  21. Sasaki, M. An introspective account of L2 writing acquisition. In Reflections on Multiliterate Lives; Multilingual Matters: Bristol, UK, 2001; pp. 291–304, Saracino 2004a. [Google Scholar]
  22. Avdeev, D.; Avdeeva, A. 3D magnetotelluric inversion using a limited-memory quasi-Newton optimization. Geophysics 2009, 74, F45–F57. [Google Scholar] [CrossRef]
  23. Feng, D.; Su, X.; Wang, X.; Wang, X. A modified total variation regularization approach based on the Gauss-Newton algorithm and split Bregman iteration for magnetotelluric inversion. J. Appl. Geophys. 2020, 178, 104073. [Google Scholar] [CrossRef]
  24. Xie, J.; Cai, H.; Hu, X.; Long, Z.; Xu, S.; Fu, C.; Wang, Z.; Di, Q. 3-D Magnetotelluric Inversion and Application Using the Edge-Based Finite Element with Hexahedral Mesh. IEEE Trans. Geosci. Remote Sens. 2021, 60, 1–11. [Google Scholar] [CrossRef]
  25. Qin, C.; Liu, X.; Wang, X.; Sun, W.; Zhao, N. Three-dimensional inversion of magnetotelluric based on adaptive finite element method. Chin. J. Geophys. 2022, 65, 2311–2325. [Google Scholar]
  26. Mueller, J.L.; Siltanen, S. Linear and Nonlinear Inverse Problems with Practical Applications; SIAM: Philadelphia, PA, USA, 2012. [Google Scholar]
  27. Zhdanov, M.S.; Fang, S. Three-dimensional quasi-linear electromagnetic inversion. Radio Sci. 1996, 31, 741–754. [Google Scholar] [CrossRef]
  28. Newman, G.A.; Alumbaugh, D.L. Three-dimensional magnetotelluric inversion using non-linear conjugate gradients. Geophys. J. Int. 2000, 140, 410–424. [Google Scholar] [CrossRef]
  29. Esparza, J.R.F.; Treviño, E.G. 2-D Niblett-Bostick magnetotelluric inversion. Geol. Acta: Int. Earth Sci. J. 2010, 8, 15–30. [Google Scholar]
  30. Ye, T.; Chen, X.-B.; Yan, L.-J. Refined techniques for data processing and two-dimensional inversion in magnetotelluric (Ⅲ); using the Impressing Method to construct starting model of 2D magnetotelluric inversion. Chin. J. Geophys. 2013, 56, 3596–3606. [Google Scholar] [CrossRef]
  31. Tikhonov, A.N.; Glasko, V.B. Use of the regularization method in non-linear problems. USSR Comput. Math. Math. Phys. 1965, 5, 93–107. [Google Scholar] [CrossRef]
  32. Last, B.; Kubik, K. Compact gravity inversion. Geophysics 1983, 48, 713–721. [Google Scholar] [CrossRef]
  33. Portniaguine, O.; Zhdanov, M.S. Focusing geophysical inversion images. Geophysics 1999, 64, 874–887. [Google Scholar] [CrossRef]
  34. Grayver, A.V.; Kuvshinov, A.V. Exploring equivalence domain in nonlinear inverse problems using Covariance Matrix Adaption Evolution Strategy (CMAES) and random sampling. Geophys. J. Int. 2016, 205, 971–987. [Google Scholar] [CrossRef]
  35. Su, Y.; Yin, C.; Liu, Y.; Ren, X.; Zhang, B.; Qiu, C.; Xiong, B. 2D magnetotelluric sparse regularization inversion based on curvelet transform. Chin. J. Geophys. 2021, 64, 314–326. [Google Scholar]
  36. Hansen, P.C.; O’Leary, D.P. The use of the L-curve in the regularization of discrete ill-posed problems. SIAM J. Sci. Comput. 1993, 14, 1487–1503. [Google Scholar] [CrossRef]
  37. Haber, E.; Oldenburg, D. A GCV based method for nonlinear ill-posed problems. Comput. Geosci. 2000, 4, 41–63. [Google Scholar] [CrossRef]
  38. Zhdanov, M.S. Geophysical Inverse Theory and Regularization Problems; Elsevier: Amsterdam, The Netherlands, 2002; Volume 36. [Google Scholar]
  39. Chen, X.-B.; Zhao, G.-Z.; Tang, J.; Zhan, Y.; Wang, J.-J. An adaptive regularized inversion algorithm for magnetotelluric data. Chin. J. Geophys. 2005, 48, 937–946. [Google Scholar] [CrossRef]
  40. Jin, B.; Takeuchi, T. Lagrange optimality system for a class of nonsmooth convex optimization. Optimization 2016, 65, 1151–1166. [Google Scholar] [CrossRef]
  41. Bahreininejad, A. Improving the performance of water cycle algorithm using augmented Lagrangian method. Adv. Eng. Softw. 2019, 132, 55–64. [Google Scholar] [CrossRef]
  42. Alikhani Koupaei, J.; Firouznia, M. A chaos-based constrained optimization algorithm. J. Ambient. Intell. Humaniz. Comput. 2021, 12, 9953–9976. [Google Scholar] [CrossRef]
  43. Rockafellar, R.T. Augmented Lagrangians and hidden convexity in sufficient conditions for local optimality. Math. Program. 2022, 1–36. [Google Scholar] [CrossRef]
  44. Wu, W.; Yang, Y.; Zheng, H.; Zhang, L.; Zhang, N. Numerical manifold computational homogenization for hydro-dynamic analysis of discontinuous heterogeneous porous media. Comput. Methods Appl. Mech. Eng. 2022, 388, 114254. [Google Scholar] [CrossRef]
  45. de Groot-Hedlin, C.; Constable, S. Inversion of magnetotelluric data for 2D structure with sharp resistivity contrasts. Geophysics 2004, 69, 78–86. [Google Scholar] [CrossRef]
  46. Gholami, A.; Aghamiry, H.S.; Operto, S. Extended-space full-waveform inversion in the time domain with the augmented Lagrangian method. Geophysics 2022, 87, R63–R77. [Google Scholar] [CrossRef]
  47. Willoughby, R.A. Solutions of ill-posed problems (AN Tikhonov and VY Arsenin). SIAM Rev. 1979, 21, 266. [Google Scholar] [CrossRef]
  48. Hestenes, M.R. Multiplier and gradient methods. J. Optim. Theory Appl. 1969, 4, 303–320. [Google Scholar] [CrossRef]
  49. Luenberger, D.G.; Ye, Y. Linear and Nonlinear Programming; Springer: New York, NY, USA, 1984; Volume 2. [Google Scholar]
  50. Miele, A.; Moseley, P.; Levy, A.; Coggins, G. On the method of multipliers for mathematical programming problems. J. Optim. Theory Appl. 1972, 10, 1–33. [Google Scholar] [CrossRef]
  51. Tapia, R.A. Diagonalized multiplier methods and quasi-Newton methods for constrained optimization. J. Optim. Theory Appl. 1977, 22, 135–194. [Google Scholar] [CrossRef]
  52. Feng, D.; Liu, J.; Wang, X. MT finite element method forward modeling and inversion using irregular quadrilateral mesh of steep topography model. J. Cent. South Univ. 2018, 49, 626–632. [Google Scholar] [CrossRef]
Figure 1. Multiplicative multi-scale correction process for data. a, b, and c are those local minimals; d is the global minimal.
Figure 1. Multiplicative multi-scale correction process for data. a, b, and c are those local minimals; d is the global minimal.
Minerals 12 01230 g001
Figure 2. Inversion model and results: (a) the model, (b) the inversion result of the conventional inversion method, and (c) the inversion result of the multiplier method.
Figure 2. Inversion model and results: (a) the model, (b) the inversion result of the conventional inversion method, and (c) the inversion result of the multiplier method.
Minerals 12 01230 g002
Figure 3. Comparison of RMS data fitting error curves.
Figure 3. Comparison of RMS data fitting error curves.
Minerals 12 01230 g003
Figure 4. Multiplier-method RMS data fitting error.
Figure 4. Multiplier-method RMS data fitting error.
Minerals 12 01230 g004
Figure 5. Diagram of the wedge model.
Figure 5. Diagram of the wedge model.
Minerals 12 01230 g005
Figure 6. Inversion results for different regularization factors.
Figure 6. Inversion results for different regularization factors.
Minerals 12 01230 g006
Figure 7. Comparison of inversion-convergence curves.
Figure 7. Comparison of inversion-convergence curves.
Minerals 12 01230 g007
Figure 8. Comparison of single-point recovery curves.
Figure 8. Comparison of single-point recovery curves.
Minerals 12 01230 g008
Figure 9. Comparison of RMS data fitting and model errors for different regularization factors.
Figure 9. Comparison of RMS data fitting and model errors for different regularization factors.
Minerals 12 01230 g009
Figure 10. Comparison of inversion for field data (the top is the result of MMTI, the bottom is the result of the CMTI).
Figure 10. Comparison of inversion for field data (the top is the result of MMTI, the bottom is the result of the CMTI).
Minerals 12 01230 g010
Figure 11. Comparison of RMS data fitting error between MMTI and CMTI.
Figure 11. Comparison of RMS data fitting error between MMTI and CMTI.
Minerals 12 01230 g011
Table 1. Model error, RMS data fitting error, and inversion time for the heave model.
Table 1. Model error, RMS data fitting error, and inversion time for the heave model.
MethodError (%)RMSTime (s)
Conventional9.570.0065770.4180
Multipliers7.410.0057740.3230
Table 2. Model errors, RMS data fitting error, and inversion times for different regularization factors of the wedge model.
Table 2. Model errors, RMS data fitting error, and inversion times for different regularization factors of the wedge model.
MethodRegularization ParameterError (%)RMSTime (s)
Conventional0.0517.180.0072206.4840
0.110.130.0040206.6550
0.320.560.0227208.2030
0.521.020.0292210.4520
Multipliers0.055.670.0029204.4870
0.14.410.0001201.8690
0.38.980.0034205.6320
0.58.980.0076208.5530
Table 3. RMS data fitting error and time of inversion for the field data.
Table 3. RMS data fitting error and time of inversion for the field data.
MethodRMSTime (s)
Conventional0.07601652.2350
Multipliers0.04081456.3620
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Feng, D.; Su, X.; Wang, X.; Ding, S.; Cao, C.; Liu, S.; Lei, Y. Magnetotelluric Regularized Inversion Based on the Multiplier Method. Minerals 2022, 12, 1230. https://doi.org/10.3390/min12101230

AMA Style

Feng D, Su X, Wang X, Ding S, Cao C, Liu S, Lei Y. Magnetotelluric Regularized Inversion Based on the Multiplier Method. Minerals. 2022; 12(10):1230. https://doi.org/10.3390/min12101230

Chicago/Turabian Style

Feng, Deshan, Xuan Su, Xun Wang, Siyuan Ding, Cen Cao, Shuo Liu, and Yi Lei. 2022. "Magnetotelluric Regularized Inversion Based on the Multiplier Method" Minerals 12, no. 10: 1230. https://doi.org/10.3390/min12101230

APA Style

Feng, D., Su, X., Wang, X., Ding, S., Cao, C., Liu, S., & Lei, Y. (2022). Magnetotelluric Regularized Inversion Based on the Multiplier Method. Minerals, 12(10), 1230. https://doi.org/10.3390/min12101230

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop