Regularized magnetotelluric inversion based on a minimum support gradient stabilizing functional

Abstract

Regularization is used to solve the ill-posed problem of magnetotelluric inversion usually by adding a stabilizing functional to the objective functional that allows us to obtain a stable solution. Among a number of possible stabilizing functionals, smoothing constraints are most commonly used, which produce spatially smooth inversion results. However, in some cases, the focused imaging of a sharp electrical boundary is necessary. Although past works have proposed functionals that may be suitable for the imaging of a sharp boundary, such as minimum support and minimum gradient support (MGS) functionals, they involve some difficulties and limitations in practice. In this paper, we propose a minimum support gradient (MSG) stabilizing functional as another possible choice of focusing stabilizer. In this approach, we calculate the gradient of the model stabilizing functional of the minimum support, which affects both the stability and the sharp boundary focus of the inversion. We then apply the discrete weighted matrix form of each stabilizing functional to build a unified form of the objective functional, allowing us to perform a regularized inversion with variety of stabilizing functionals in the same framework. By comparing the one-dimensional and two-dimensional synthetic inversion results obtained using the MSG stabilizing functional and those obtained using other stabilizing functionals, we demonstrate that the MSG results are not only capable of clearly imaging a sharp geoelectrical interface but also quite stable and robust. Overall good performance in terms of both data fitting and model recovery suggests that this stabilizing functional is effective and useful in practical applications.

Introduction

The MT sounding is widely used in geophysical research and plays an important role in the exploration of oil and gas resources and the basic research of deep geological structures. Inversion of MT data has a direct impact on geological interpretation and thus is a key point in the MT exploration. Currently, the most widely used inversion methods include the followings: Occam’s inversion, which yields a model with the smallest roughness for a specified misfit, providing a stable and rapidly convergent solution (Constable et al. 1987; de Groot-Hedlin and Constable 1990); rapid relaxation inversion (RRI), in which horizontal derivative terms are approximated from the fields of previous iterations (Smith and Booker 1991); reduced basis Occam’s inversion (REBOCC), which transforms the linearized inverse problem from the model space to the data space (Siripunvaraporn and Egbert 2000); and nonlinear conjugate gradients (NLCG), which avoid excessive evaluations of the full Jacobian matrix and the complete solution of a linearized inverse problem at each step of iteration (Rodi and Mackie 2001). These inversion methods usually cast the mathematical constraint of the model with the minimum model, the flattest model, or the smoothest model into the objective functional to relieve the nonuniqueness of the inverse problem. As a result, these algorithms provide stable and smooth inversion results. However, it is difficult for these methods to clearly image the shape and sharp boundary of geological structures, although such an imaging is required in many practical situations (Marcuello-Pascual et al. 1992; Smith et al. 1999; Zhdanov 2004, 2009, 2010).

Aiming for a clear imaging of electrical interfaces, Last and Kubik (1983) first considered the use of a MS functional that provides models with the minimum domain volume with anomalous parameter distribution. Rudin et al. (1992) introduced a method for reconstruction of noisy and blurred images based on total variation (TV). Farquharson and Oldenburg (1998) and Farquharson (2008) modified the typical minimum-structure inversion with an iteratively reweighted least-squares procedure and used l ₁-type measure to generate blocky, piecewise-constant earth models. Portniaguine and Zhdanov (1999) demonstrated that the use of a MGS functional minimizes the area with a finite gradient of the model parameters and helps to generate focused images. The MGS has been successfully applied to the inversion of gravitational, magnetic, and electromagnetic data. de Groot-Hedlin and Constable (2004) developed a linearized method that uses smooth boundaries to define sharp resistivity contrasts. This method penalizes variations in the boundary depth rather than in the contrast between the resistivity of adjacent blocks and yields more geologically realistic results. They concluded that this method is also suitable for situations in which the value of the target anomaly is known to have a strong conductivity contrast with the background medium. However, it needs to build a layered model to invert both resistivity and thickness of each layer for every horizontal mesh, which is very complicated in the model construction and inversion calculation. Meanwhile, Zhang et al. (2009, 2010) added diagonal gradient support to the objective functional, which consists of both smoothing and focusing constraints, and obtained better inversion results in imaging an inclined geoelectrical interface. Recently, Grayver and Kuvshinov (2016) applied arbitrary norms for residual and regularization terms to produce classes of equivalent solutions of smooth or compact models, in which a L _1.5 norm was shown to be most efficient for a stable inversion with sharp boundary.

Here our research intends on focusing a sharp boundary, following general Tikhonov regularization inversion within the framework of L₂ norm, which is justified by general Gaussian distribution of EM data (Fournier and Febrer 1976; Weaver et al. 2000). In this paper, we propose a new stabilizing functional called the MSG functional, which adds stable constraints to help focusing a sharp boundary, and apply it to a unified form of the objective functional in the regularized inversion. We perform the inversion of synthetic data to compare the results with those obtained by the conventional stabilizing functionals. The implications of the proposed functional are also discussed.

Theory and algorithm

A general way to solve the nonuniqueness of geophysical inverse problems is first described in the regularization theory by Tikhonov and Arsenin (1977), which uses both data misfit and model stabilizing functionals to construct the parametric objective functional. The objective functional of a regularized inversion can generally be written as

$$P({\mathbf{m,d}})^{\alpha } = f({\mathbf{m}},{\mathbf{d}}) + \alpha s({\mathbf{m}}),$$

(1)

where P(m,d)^α is the parametric objective functional; m is a vector of model parameters; f(m, d) is the data objective functional and can be given by $f({\mathbf{m,d}}) = \left\| {{\mathbf{W}}_{\text{d}} A({\mathbf{m}}) - {\mathbf{W}}_{\text{d}} {\mathbf{d}}} \right\|^{2}$, here ${\mathbf{W}}_{\text{d}}$ is the data weighting matrix, $\left\| {\mathbf{g}} \right\|^{2} = {\mathbf{g}}^{\text{T}} g$ means the two norm, and ^T denotes the transpose; A is an operator of the theoretical data from model parameters; d is a vector containing the observation data; s(m) is the model stabilizing functional; and α is the regularization parameter. Logarithmic scaling is applied to both the observational data (MT impedance elements) and a model parameter (logarithm of the electrical resistivity in the present case).

A stabilizing functional is used to select an appropriate class of models from the space of all possible models. The following stabilizing functionals have been applied successfully in many inversion studies to yield stable, smooth solutions (Constable et al. 1987; de Groot-Hedlin and Constable 1990; Smith and Booker 1991; Uchida 1993; Rodi and Mackie 2001). The simplest is the minimum model stabilizing functional (MM), which is based on the least-squares criterion and uses the minimum norm of the difference from the a priori model m _apr. The MM stabilizing functional is given as

$$s_{\text{MM}} ({\mathbf{m}}) = \sum\limits_{i = 1}^{N} {\left( {m^{(i)} - m_{\text{apr}}^{(i)} } \right)^{2} } ,$$

(2)

where m ⁱ is the ith model parameter and N is the total number of model parameters. If we consider the spatial gradient $\nabla {\mathbf{m}}$ of the model parameters to enhance the stability of the inversion, we obtain the flattest model (FM) stabilizing functional. We calculate the gradient in each direction separately for 2D or 3D case, but the weight on the gradient in each direction is the same so as to compare the results from different stabilizers using spatial gradients in the same condition. Here we describe generally 1D form of each functional. In this formulation, the FM stabilizer is given as:

$$s_{\text{FM}} ({\mathbf{m}}) = \sum\limits_{i = 1}^{N} {\left( {\nabla \left( {m^{(i)} - m_{\text{apr}}^{(i)} } \right)} \right)^{2} .}$$

(3)

The norm of the spatial Laplacian $\nabla^{2} {\mathbf{m}}$ of the model parameters is used as another stabilizing functional to produce the smoothest model (SM):

$$s_{\text{SM}} ({\mathbf{m}}) = \sum\limits_{i = 1}^{N} {\left( {\nabla^{2} \left( {m^{(i)} - m_{\text{apr}}^{(i)} } \right)} \right).}$$

(4)

The modified TV stabilizing functional requires the distribution of model parameters of bounded variation (Acar and Vogel 1994) which is given as:

$$s_{\text{TV}} ({\mathbf{m}},\beta ) = \sum\limits_{i = 1}^{N} {\left( {\sqrt {\nabla \left( {m^{(i)} - m_{\text{apr}}^{(i)} } \right)^{2} + \beta^{2} } } \right)} ,$$

(5)

where β is a small value called the focusing parameter. It tends to decrease the bounds of variation of model parameters and therefore is effective in resolving sharp boundary to some extent (Zhdanov 2002). Farquharson (2008) generated sharp interfaces through adding first-order finite-difference diagonal weighting operators based on L ₁ measure. Small β tends toward L ₁ norm and large β tends like a scaled sum-of-squares measure, so different β will produce different inversion results.

The MS stabilizing functional is proposed to build a model that minimizes the volume with anomalous model parameters (Last and Kubik 1983), focusing on the volumes in which the model parameters differ from the given a priori model. In this case, the MS stabilizing functional is given as

$$s_{\text{MS}} ({\mathbf{m}},\beta ) = \sum\limits_{i = 1}^{N} {\left( {\frac{{m^{(i)} - m_{\text{apr}}^{(i)} }}{{\sqrt {(m^{(i)} - m_{\text{apr}}^{(i)} )^{2} + \beta^{2} } }}} \right)^{2} } ,$$

(6)

where β is the focusing parameter. If the difference in model parameters is much smaller than β, the functional becomes significantly smaller than N, and it approaches N when the difference is much greater than β.

Analogously, Portniaguine and Zhdanov (1999) developed the MGS functional to minimize the volume where the model parameter varies steeply and discontinuously. The MGS stabilizing functional is written as

$$s_{\text{MGS}} ({\mathbf{m}},\beta ) = \sum\limits_{i = 1}^{N} {\left( {\frac{{\nabla (m^{(i)} - m_{\text{apr}}^{(i)} )}}{{\sqrt {(\nabla (m^{(i)} - m_{\text{apr}}^{(i)} ))^{2} + \beta^{2} } }}} \right)^{2} } .$$

(7)

When the gradient of a model parameter is much smaller than β, the functional takes a value much smaller than N, whereas it approaches N when the gradient is much larger than β. In the latter case, the solution of the inversion images the sharp boundaries but strongly rely on the selection of the focusing parameter β. Zhang et al. (2012) proposed the estimation of the proper value of β by referring to the cutoff length, but it is difficult to estimate a priori the actual range of the cutoff in general. In this case, tuning by trial and error is necessary to obtain an optimum solution.

This study aims to construct a new stabilizer that is less dependent on the tuning of the focusing parameter. Thus, considering the advantages and disadvantages of each of these stabilizing functionals, we propose a new stabilizing functional, the minimum support gradient, which is defined as

$$s_{\text{MSG}} ({\mathbf{m}},\beta ) = \sum\limits_{i = 1}^{N} {\left( {\nabla \frac{{(m^{(i)} - m_{\text{apr}}^{(i)} )}}{{\sqrt {(m^{(i)} - m_{\text{apr}}^{(i)} )^{2} + \beta^{2} } }}} \right)^{2} } .$$

(8)

The spatial gradient is calculated after the MS stabilizer is obtained using Eq. (6), and thus, we call this new stabilizing functional the MSG. It uses MS value instead of $\left( {m^{(i)} - m_{\text{apr}}^{(i)} } \right)$ in Eq. (3), which leads to a stable constraint of focusing sharp boundaries for inversion. It works on the class of models so as the domain with anomalous parameter distribution occupies the minimum volume, while it is expected to avoid the instability caused by the tuning of focusing parameter at the same time. And the MSG stabilizing functional is affected little by a priori model. Zhdanov (2002) proved that the MS and MGS satisfy the Tikhonov criterion for a regularization stabilizer. Since the MSG is a spatial gradient of the MS, it also satisfies the regularization criterion.

As for the numerical evaluation of spatial gradient $\nabla$, we simply approximate it by a difference between values in neighboring cells of a numerical model. Therefore, this is not exactly a finite-difference approximation of the gradient operator, but we confirmed that little difference is produced. Thus, the stabilizing functionals for the MS and MSG can formally be expressed as the following pseudo-quadratic functional of model parameters:

$$s({\mathbf{m}}){ = }\sum\limits_{i = 1}^{N} {\left( {{\mathbf{W}}_{{\mathbf{e}}} \left( {m^{(i)} - m_{\text{apr}}^{(i)} } \right)} \right)^{2} } ,$$

(9)

where W _e is a N × N matrix consisting of m, m _apr, and β. We then use the matrix notation of each stabilizing functional to build the unified form of the objective functional for the discrete model parameters (Zhdanov 2002). Using this notation, the discrete form of Eq. (1) is

$$P({\mathbf{m}},{\mathbf{d}})^{\alpha } = \sum\limits_{j = 1}^{\text{ND}} {\left( {A(m)^{(j)} - d^{(j)} } \right)^{2} } + \alpha \sum\limits_{i = 1}^{N} {\left( {{\mathbf{W}}_{{\mathbf{e}}} \left( {m^{(i)} - m_{\text{apr}}^{(i)} } \right)} \right)^{2} } .$$

(10)

where d ^j is the jth observed data and ND is the total number of observed data.

Here, we examine the characteristics of the stabilizing functionals suggested above. The stabilizing matrix for MS in this case could be written as

$${\mathbf{W}}_{{\mathbf{e}}}^{{{\mathbf{MS}}}} = \left( {\begin{array}{*{20}c} {\frac{1}{{\sqrt {\left( {m_{k}^{\left( 1\right)} - m_{\text{apr}}^{(1)} } \right)^{2} + \beta^{2} } }}} & {} & {} & {} & 0 \\ {} & \ldots & {} & {} & {} \\ 0 & {} & {\frac{1}{{\sqrt {\left( {m_{k}^{\left( i \right)} - m_{\text{apr}}^{\left( i \right)} } \right)^{2} + \beta^{2} } }}} & {} & 0 \\ {} & {} & {} & \ldots & {} \\ 0 & {} & {} & {} & {\frac{1}{{\sqrt {\left( {m_{k}^{\left( N \right)} - m_{\text{apr}}^{(N)} } \right)^{2} + \beta^{2} } }}} \\ \end{array} } \right) ,$$

(11)

where $m_{k}^{(N)}$ is the Nth model parameter of the kth iteration and β is the focusing parameter. The stabilizing matrix of the MSG for 1D case can then be derived as

$${\mathbf{W}}_{{\mathbf{e}}}^{{{\mathbf{MSG}}}} = \left( {\begin{array}{*{20}c} 0 & 0 & {} & {} & {} & 0 \\ {\frac{ - 1}{{\sqrt {\left( {m_{k}^{\left( 1 \right)} - m_{\text{apr}}^{\left( 1 \right)} } \right)^{2} + \beta^{2} } }}} & {\frac{1}{{\sqrt {\left( {m_{k}^{\left( 2 \right)} - m_{\text{apr}}^{\left( 2 \right)} } \right)^{2} + \beta^{2} } }}} & {} & {} & {} & 0 \\ {} & \ldots & \ldots & {} & {} & {} \\ 0 & 0 & {\frac{ - 1}{{\sqrt {\left( {m_{k}^{\left( i \right)} - m_{\text{apr}}^{\left( i \right)} } \right)^{2} + \beta^{2} } }}} & {\frac{1}{{\sqrt {\left( {m_{k}^{{\left( {i + 1} \right)}} - m_{\text{apr}}^{{\left( {i + 1} \right)}} } \right)^{2} + \beta^{2} } }}} & {} & 0 \\ {} & {} & {} & \ldots & \ldots & {} \\ 0 & 0 & {} & {} & {\frac{ - 1}{{\sqrt {\left( {m_{k}^{{\left( {N - 1} \right)}} - m_{\text{apr}}^{{\left( {N - 1} \right)}} } \right)^{2} + \beta^{2} } }}} & {\frac{1}{{\sqrt {\left( {m_{k}^{\left( N \right)} - m_{\text{apr}}^{\left( N \right)} } \right)^{2} + \beta^{2} } }}} \\ \end{array} } \right) .$$

(12)

In 1D case, the gradient is taken just in the vertical direction. In 2D or 3D situation, however, we need to separate the gradient calculation in each direction.

We examine how different values of β affect the stabilizing performance by using a simple 1D model which is given by

$$m(x) = \left\{ {\begin{array}{cl} {0,} \hfill &\quad {0 \le x \le 1} \hfill \\ {{\text{e}}^{{( - 32*(x - 1.5)^{2} )}} ,} \hfill &\quad {1 < x \le 1.5} \hfill \\ {1,} \hfill &\quad {1.5 < x \le 2.5} \hfill \\ {{\text{e}}^{{( - 32*(x - 2.5)^{2} )}} ,} \hfill &\quad {2.5 < x \le 3} \hfill \\ {0,} \hfill &\quad {3 < x \le 4} \hfill \\ \end{array} } \right.,$$

(13)

where x is the logarithm of the depth and has a value that changes from 0 to 4 in increments of 0.04, and m _apr is assumed to be 0 over the whole range of x. The behaviors of the different stabilizing functionals after normalization are shown in Fig. 1. We see that stabilizing functionals, such as the TV, MGS, and MSG functionals, tend to have a narrow peak around the boundary of the model parameter, indicating effective focusing. This property is expected to help us obtain a sharp image of the electrical interface. If the value of the focusing parameter is changed, features of the TV and MGS are significantly altered; the effect on the sharp imaging becomes weak when the focusing parameter is increased. Conversely, its influence on the amplitude of MSG is not very large in comparison with those of the TV or MGS. The MSG stabilizer shows a similar effect on focusing the boundary even if the value of focusing parameter is varied over a wide range (Fig. 1). Different values of the focusing parameter will have different MS values, so that the position of MS boundary should change accordingly. The gradient of the MSG is calculated after the calculation of the MS, so the boundary positions of MSG are expected to change with changes in focusing parameter. Nevertheless, it still focuses the right boundary position (x-axis from 1 to 1.5 and 2.5 to 3).

Fig. 1

Comparison of mathematical model distributions of different stabilizing functionals with square of focusing parameters β ² of a 1, b 0.1, and c 0.01

We can thus choose different optimization methods to minimize the objective functional so that inversion can be performed by using different stabilizing functionals in the same framework. In the next section, we perform synthetic tests using 1D-layered and 2D wedge models to compare the performance of inversion with different stabilizing functionals. Here, we applied Occam’s inversion (Constable et al. 1987) with a self-adaptive conjugate gradient weighted optimization scheme for 1D case. We used the ratio of data misfit and the value of each stabilizing functional as the initial value of regularization parameter, we kept the same regularization parameter if the data misfit decreases during the inversion, while we set $\alpha_{n + 1} = 0.9*\alpha_{n}$ if the data misfit did not decrease (Zhdanov 2002). We applied the original 2D Occam’s inversion with the Cholesky decomposition to obtain inverse matrix (de Groot-Hedlin and Constable 1990) for 2D case. The regularization parameters were determined exactly as suggested in the original paper (Constable et al. 1987).

Synthetic model test

One-dimensional model

Here we test 1D models with sharp and smooth boundaries to observe the effect of different stabilizing functionals, i.e., smoothing functionals (MM, FM, and SM) and focusing functionals (TV, MS, MGS, and MSG). We use 80 frequencies from 1000 to 0.001 Hz for the forward calculation to generate synthetic data. Then, 1% Gaussian noise was added to the synthetic data and the observation error was assumed also to be 1%. If we applied higher noise and error, the result of each stabilizing functional easily reached the target misfit so that all results showed little difference. Then, we calculated apparent resistivity, and a logarithmic scaling was applied.

The synthetic 1D-layered model is shown by a gray thick line in Fig. 2a (model A) and b (model B). The models with opposite contrast are given by the function

$$m(x) = \left\{ {\begin{array}{*{20}l} {2,} \hfill &\quad {x < - 1} \hfill \\ {3,} \hfill &\quad { - 1 \le x \le 0} \hfill \\ {1,} \hfill &\quad {0 < x \le 1} \hfill \\ {2,} \hfill & \quad {x > l} \hfill \\ \end{array} } \right.,$$

(14)

and

$$m(x) = \left\{ {\begin{array}{*{20}l} {2,} \hfill &\quad {x < - 1} \hfill \\ {1,} \hfill &\quad { - 1 \le x \le 0} \hfill \\ {3,} \hfill &\quad {0 < x \le 1} \hfill \\ {2,} \hfill & \quad {x > l} \hfill \\ \end{array} } \right.,$$

(15)

where x is the logarithm of the depth, and it consists of one high-resistivity anomalous layer and one low-resistivity anomalous layer embedded in a uniform half-space with a resistivity of 100 Ω m. The resistivity and thickness of the shallower anomaly are 1000 Ω m and 0.9 km ranging from 0.1 to 1 km in depth, and those of the deeper anomaly are 10 Ω m and 9 km ranging from 1 to 10 km. The common boundary of the two anomalous layers is 1 km (x = 0) in depth. The synthetic sounding curve of each model is shown in Fig. 2c, d. The model is discretized into 40 layers, and the maximum depth is 100 km.

Fig. 2

1D-layered model inversion results of a model A, b model B, c synthetic sounding curve of the model A, and d synthetic sounding curve of model B. x in horizontal axis is the logarithm of the depth, and the unit of depth is km. Thick gray line indicates a profile of the synthetic model, and thin lines with different colors correspond to different stabilizing functionals as shown in the legend on the top

The initial model is a 100-Ω m uniform half-space, and the a priori model is the same as the initial model. The dependence on the starting model is mainly from the different optimization methods to minimize the objective functional. Here we use Occam’s inversion method to solve the inverse problem whose solution is well known to be affected little by the choice of the initial model. The choice of a priori model, m _apr, is independent from the choice of the initial model. It was confirmed that Occam’s inversion is stable so far as m _apr is not too far from the true background value (100 Ω m in this case).

We calculated the vertical gradient of model parameters to determine the position and sharpness of the interface (Fig. 3c, d). To evaluate how the model responses fit the synthetic data, we define the root-mean-square (RMS) data misfit as

$${\text{RMS}}_{\text{d}} = \sqrt {\frac{{\sum\nolimits_{j = 1}^{\text{ND}} {\left[ {(d_{j}^{\text{obs}} - d_{j}^{\text{cal}} )/{\text{err}}_{j} } \right]^{2} } }}{\text{ND}}} ,$$

(16)

where ND is the number of data and err _j is the observation error. For Occam’s inversion, the target misfit is set to 1. Here we set the error floor equal to the assumed noise level, so it is possible to reduce RMS misfit to 1 by regularized or nonregularized inversion. However, in some cases, the solution suffered from overfitting. To evaluate how the synthetic model parameters are recovered by inversion, we also calculate the RMS model recovery (Zhang et al. 2012) defined as

$${\text{RMS}}_{\text{m}} = \sqrt {\frac{{\sum\nolimits_{k = 1}^{N} {\left( {m_{k}^{\text{inv}} - m_{k}^{\text{model}} } \right)^{2} } }}{N}} .$$

(17)

Fig. 3

1D-layered model inversion results of different stabilizing functionals with β ² = 0.1: a inversion results of the model A, b inversion results of the model B, c, d the gradient of the each results, e, f normalized gradient of each results, and g, h convergence of RMS data misfit of each model

The obtained RMS_d and RMS_m values are given in Tables 1 and 3. For a detailed discussion, we also separately calculated the RMS model recovery of each electrical boundary (for the range of x from − 1.5 to − 0.5, from − 0.5 to 0.5, and from 0.5 to 1.5, respectively) for the two models, as given in Tables 2 and 4.

Table 1 RMS data misfit and RMS model recovery of different stabilizing functionals for the 1D-layered model A (Fig. 2a)

Table 2 RMS model recovery around each boundary of the 1D-layered model A (Fig. 2a)

Table 3 RMS data misfit and RMS model recovery of different stabilizing functionals for the 1D-layered model B (Fig. 2b)

Table 4 RMS model recovery around each boundary of the 1D-layered model B (Fig. 2b)

Both for models A and B, all inversion results reached the target misfit even with different values of β ² in both two cases. When β ² was set to 0.1 (Fig. 3a, b), TV, MS, MGS, and MSG had a similar RMS model recovery for model A, significantly smaller than smooth results. All of them had a better recovery for the boundary around x = 0 (Table 2). In model B, the MGS had the smallest RMS model recovery of 0.2464. However, the RMS model recovery of the MS and MSG was also as small as 0.3628 and 0.3307, respectively, while TV had a poor RMS model recovery as 0.3926. The MGS results showed an accurate boundary at x = 0 (Fig. 3c, d), and the MSG results showed a more accurate boundary location at x = − 1 and x = 1.

Next we varied the focusing parameter value to 0.001 and 0.0001 and observed the influence. With decreasing β ² to 0.001, the focusing stabilizers had a better effect on imaging the sharp interface than those results with smooth model constraints (Fig. 4a, b). With the MSG and MGS stabilizers, the position of the boundary was imaged more accurately than that with the MS and TV at x = 1. At x = 0, MGS showed a more sharp boundary. However, the MS and MGS results appear to be unstable and produce some false features around x = − 0.5 and x = 0.8 (black circles in Fig. 4c). For model B, the MSG result had the smallest RMS model recovery of 0.2396, and it had a better performance for imaging every boundary (Fig. 4d and Table 4). However, TV and MGS produced false features around x = − 0.5 and x = 0.5 (Fig. 4d), respectively. When β ² was decreased further to 0.0001 (Fig. 5a, b), for model A, the MS results produce more false features around x = 0.5, while MGS had large model recovery because of its false features occurred around x = − 0.5 to x = 1.2. The MSG result still remains stable, accurately recovering both the resistivity values and interface depths, especially x = 0 (Fig. 5c, e). For model B, TV and MGS produced lots of false features around x = − 0.8 and x = 0.2 (Fig. 5d). From the comparison of models A and B, the boundary changing from resistive to conductive (x = 0 in Figs. 3a, 4a, 5a and x = − 1 in Figs. 3b, 4b, 5b) is imaged more sharply than the boundary of opposite contrast.

Fig. 4

The same as Fig. 3 with β ² = 0.001

Fig. 5

The same as Fig. 3 with β ² = 0.0001

Although the TV, MS, and MGS yield the sharp boundaries in some cases (boundaries of x = 1 in Figs. 3a, 4a and x = 0 in Figs. 3b, 4b), they also produce numerous false structures (at depths of x = − 0.5 in Figs. 4a, 5a and x = 0.5 in Figs. 4b, 5b), suggesting the inversion is unstable due to the small value of the focusing parameter. For a given stabilizer, we also calculated the difference between each pair of inversion results produced by different β ² values (β ₁ and β ₂) as

$${\text{Diff}}_{\text{m}} = \sum\limits_{k = 1}^{N} {\left( {m_{k}^{{\beta_{1} }} - m_{k}^{{\beta_{2} }} } \right)^{2} } .$$

(18)

These differences for different stabilizers are summarized in Tables 5 and 6. The MSG results had the smallest differences when selecting different values of β ², which means the MSG inversion is more robust against variation in β ² than those of the TV, MS, or MGS. It is a great advantage of using the MSG stabilizer that fine-tuning of β ² is not necessary to obtain a stable solution.

Table 5 Differences [calculated using Eq. (18)] between the inversion results of model A (Fig. 2a) obtained using different values of β ²

Table 6 Differences (calculated using Eq. (18)) between the inversion results of model B (Fig. 2b) obtained using different values of β ²

The third 1D synthetic model shows spatially smooth variation in the electrical conductivity and is represented by a thick gray line in Fig. 6. The model is given by the function

$$\text{m} (x) = \left\{ {\begin{array}{cl} {1,} \hfill & {x < - 1} \hfill \\ {\sin ((x + 1) * \pi * 10/7) + 1,} \hfill & { - 1 \le x \le - 0.3} \hfill \\ {1,} \hfill & { - 0.3 < x \le 0.4} \hfill \\ {\sin ((x + 0.2) * \pi * 5/3) + 1,} \hfill & {0.4 < x \le 1} \hfill \\ {1,} \hfill & {x > \text{l} } \hfill \\ \end{array} } \right. .$$

(19)

Thus, the anomaly changes at the logarithm depths ranging from 0.1 to 0.5 km and from 2.5 to 10 km in a background resistivity of 10 Ω m. The initial and a priori models are both in a 10-Ω m half-space.

Fig. 6

1D smooth boundary model inversion results with different stabilizing functionals. Inversion results with β ² of a 0.1, b 0.001, and c 0.0001

The inversion results from each stabilizing functional with different β ² values are shown in Fig. 6. The data misfit and model recovery are given in Table 7. The TV, MS, MGS, and MSG results have discontinuities where the synthetic model parameters change smoothly. The FM and SM results show excellent performance for conductive anomaly, as their RMS data misfit values both reached 1 and their RMS model recovery values became small, reaching 0.1707 and 0.1744, respectively. Thus, the results imply that the FM or SM is a better choice for cases with smooth boundaries, although results using all focusing stabilizers also show an acceptable performance both in terms of RMS misfit and model recovery.

Table 7 RMS data misfit and RMS model recovery of different stabilizing functionals for 1D smooth boundary model

Two-dimensional wedge model

For further discussion, we test all five stabilizing functionals using 2D synthetic models with a resistive and conductive wedge (Fig. 7a, b). The 2D wedge model configuration is the same as that used by de Groot-Hedlin and Constable (2004) and Zhang et al. (2009). The upper boundary of the wedge has a slope of 5.7°, and the lower boundary has a slope of 16.6°. In total, 24 synthetic observation sites are distributed regularly at intervals of 0.5 km. For the modeling, 92 mesh grids are defined in the horizontal direction with a 0.125-km spacing, and 100 layers are set in the vertical direction, which is discretized with the same increment. We used the finite-element code in Occam’s inversion for the modeling. The mesh spacing for the inversion is set twice as large as that in the forward modeling, using 46 meshes in the horizontal direction and 50 layers in the vertical direction. We used 20 frequencies from 4 to 0.0063 Hz for the resistive wedge model and 24 frequencies from 400 to 0.01 Hz for the conductive wedge model to perform the modeling that generates synthetic data (de Groot-Hedlin and Constable 2004). Both transverse electric (TE) and transverse magnetic (TM) modes were used in the calculation. The synthetic sounding curves of apparent resistivity and phase from site 10 (red triangle in Fig. 7a) of the resistive wedge case are shown in Fig. 8 as an example.

Fig. 7

Synthetic 2D models. a Resistive wedge model, b conductive wedge model. Black triangle indicates the site location for which the synthetic TE and TM responses are shown in Fig. 8

Fig. 8

Synthetic 2D apparent resistivity and phase curves of site 10 for the resistive model: a TE apparent resistivity, b TE phase, c TM apparent resistivity, and d TM phase. Observation errors corresponding to 1% of impedance amplitude are also shown

In the 2D resistive wedge model, the background resistivity is 1 Ω m and that of the wedge-shaped body is 100 Ω m. We allow the resistivity values to change from 0.01 to 500 Ω m during the inversion. The starting model is a 1-Ω m half-space, and the a priori model is also a 1-Ω m uniform half-space. Again, the different initial model affects little, and the regularization works well when a priori model is not too far from the true background resistivity (1 Ω m in this case). We add 1% random noise in the synthetic data and set the error floor as 1% in the inversion. Then, it is not easy to get the data misfit small. We also tried cases with larger error floor. In those cases, it was quite easy to get data misfit down, but resulting model was much less sharp than the present results. It depends on the noise and error given to the synthetic data whether we can easily get the data misfit down to an RMS value of 1.

We calculated the RMS data misfit and the model recovery, the definitions of which are the same as in the 1D case (Table 8). The inversion using the smooth stabilizing functionals, such as the MM (Fig. 9a), FM (Fig. 9b), and SM (Fig. 9c), reached the target misfit, and their RMS model recovery values are reduced only a little to 0.4584, 0.3583, and 0.3456, respectively. This means that none of these stabilizing functionals helps yield the right shape and value of the wedge. Here the result with TV is not shown because the inversion did not converge, and most probably because the TV stabilizer constrains L₁ norm of the spatial gradient, which is not consistent with the framework of Occam’s inversion with L ₂ norm.

Table 8 RMS data misfit and RMS model recovery of different stabilizing functionals for 2D resistive model

Fig. 9

Synthetic inversion results for the 2D resistive model (Fig. 6a) using different stabilizing functionals a MM, b FM, and c SM. Color scale is the same as Fig. 6a

Results using focusing (MS, MGS, and MSG) functionals and variations for different focusing parameters are shown in Figs. 10, 11, 12. They are generally better than the results using smooth constraints (Table 8). We then compared the results using different focusing parameters. When β ² is 0.0001, the MS result (Fig. 10a) has a model recovery of 0.4302 and shows some false structures below the wedge. As a whole, the shape of the wedge body is not restored well. The MGS (Fig. 10b) and MSG results (Fig. 10c) both have smaller RMS model recovery values (0.3242 and 0.2499, respectively) than those of the results with smooth constraints. However, the bottom boundary of the wedge for the MGS result is not as accurate as that for the MSG result.

Fig. 10

Synthetic inversion results for the 2D resistive model (Fig. 6a) using different stabilizing functionals with β ² = 0.0001: a MS, b MGS, c MSG results, and d convergence of RMS data misfit. Color scale is the same as Fig. 6a

Fig. 11

Synthetic inversion results for the 2D resistive model (Fig. 6a) using different stabilizing functionals with β ² = 0.001: a MS, b MGS, c MSG results, and d convergence of RMS data misfit. Color scale is the same as Fig. 6a

Fig. 12

Synthetic inversion results for the 2D resistive model (Fig. 6a) using different stabilizing functionals with β ² = 0.01: a MS, b MGS, c MSG results, and d convergence of RMS data misfit. Color scale is the same as Fig. 6a

If β ² is increased to 0.001, each inversion result varies in different degrees. The MS (Fig. 11a) and MGS results (Fig. 11b) exhibit considerable changes. Both of them perform better than previous results in terms of the model recovery (0.2983, 0.3164, respectively). The results of the MSG stabilizer (Fig. 11c) are more stable and acceptable than those of other stabilizing functionals, as they reached the target misfit and achieved the smallest model recovery (0.2634). Further increasing β ² to 0.1, the MS and MGS results changed drastically, but the MSG results still achieve stability and remain the best among results using different functionals (model recovery of 0.3005). Additionally, the MS (Fig. 12a) and MGS results (Fig. 12b) produced some false structures, in both the shape and resistivity of the wedge. Conversely, the MSG results (Fig. 12c) successfully reproduced the shape and the resistivity of the wedge. As a whole, the MSG stabilizer shows the best performance as can be seen in a general view of the RMS model recovery with different focusing parameters (Fig. 13).

Fig. 13

The model recovery curve with different focusing parameters for the 2D resistive model (Fig. 6a)

Thus, the MSG stabilizing functional was shown to be able to not only describe the sharp boundary but also estimate the model parameter more accurately than other functionals. In addition, it is very stable in the regularized inversion and is influenced by the focusing parameter value much less than the MS and MGS stabilizers.

We also performed synthetic tests for a 2D conductive wedge model (Fig. 7b). Figure 14 shows a result of inversion using the MSG regularization. Again we can conclude that the MSG stabilizer provides a good performance among all stabilizers tested. The synthetic inversion results for 1D case in “One-dimensional model” section showed that MT is insensitive to the sharpness in the vertical direction. The 2D results are sharper than 1D results, indicating that the sharpness of inclined interface of 2D cases is mostly due to the horizontal gradient.

Fig. 14

A synthetic (a) inversion result with the MSG stabilizer for the 2D conductive wedge model (Fig. 6b). Color scale is the same as that in Fig. 6b

Conclusion

In this paper, to consider the stability on sharp boundaries, we proposed a new stabilizing functional called the MSG stabilizing functional for imaging sharp interfaces of resistivity structure by magnetotelluric inversion. We compared its performance with other stabilizers under the same inversion framework through the unified objective functional form of a discrete weighted matrix. By comparing 1D and 2D synthetic inversed results using different stabilizers, we confirmed that the MSG stabilizer was not only able to resolve a geoelectrical interface but also stable for a wide range of the focusing parameter value. The MSG stabilizer also showed good performance regarding both data misfit and model recovery from synthetic data inversion with sharp boundaries. The overall results of the present study suggest that this stabilizing functional would be effective in practical applications.