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.
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)
1D-layered model inversion results of different stabilizing functionals with β 2 = 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 RMSd and RMSm 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.
Both for models A and B, all inversion results reached the target misfit even with different values of β 2 in both two cases. When β 2 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 β 2 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 β 2 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.
The same as Fig. 3 with β 2 = 0.001
The same as Fig. 3 with β 2 = 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 β 2 values (β 1 and β 2) 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 β 2, which means the MSG inversion is more robust against variation in β 2 than those of the TV, MS, or MGS. It is a great advantage of using the MSG stabilizer that fine-tuning of β 2 is not necessary to obtain a stable solution.
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.
1D smooth boundary model inversion results with different stabilizing functionals. Inversion results with β 2 of a 0.1, b 0.001, and c 0.0001
The inversion results from each stabilizing functional with different β 2 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.
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.
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
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 L1 norm of the spatial gradient, which is not consistent with the framework of Occam’s inversion with L 2 norm.
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 β 2 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.
If β 2 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 β 2 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).
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.
