The problem of model recovering in the presence of impulse noise on the data is considered for the magnetotelluric (MT) inverse problem. The application of total variation regularization along with L1-norm penalized data fitting (TVL1) is the usual approach for the impulse noise treatment in image recovery. This combination works poorly when a high level of impulse noise is present on the data. A nonconvex operator named smoothly clipped absolute deviation (TVSCAD) was recently applied to the image recovery problem. This operator is solved using a sequence of TVL1 equivalent problems, providing a significant improvement over TVL1. In practice, TVSCAD requires the selection of several parameters, a task that can be very difficult to attain. A more simple approach to the presence of impulse noise in data is presented here. A nonconvex function is also considered in the data fitness operator, along with the total variation regularization operator. The nonconvex operator is solved by following a half-quadratic procedure of minimization. Results are presented for synthetic and also for field data, assessing the proposed algorithm’s capacity in model recovering under the influence of impulse noise on data for the MT problem.
Electromagnetic sounding methods have been applied for a long time to investigate the interior of the Earth from depths of a few meters to hundreds of kilometers. The magnetotelluric (MT) sounding method is one of the most popular today, and it is aimed at estimating the electrical conductivity distribution at depth associated with subsurface geological structures. It has been successfully applied in shallow exploration problems and to infer the Earth’s deep electrical conductivity distribution, temperature regimes, and geological structure. The method uses measurements of natural electric and magnetic fields recorded at the surface of the Earth. If it is assumed that the underground structure is one-dimensional (1D), the inverse problem consists of recovering the vertical conductivity distribution from the ratio of surface electric and magnetic field measurements taken at different frequencies. One of the main characteristics of the MT inverse problem is its nonlinearity. Several methods have been developed to cope with this problem : linearization, asymptotic behavior, and exact methods. Constable et al.  consider smooth models using a procedure called Occam’s inversion. The method consists of Tikhonov’s regularization applied to the linearized equations of the 1D MT problem. A fundamental characteristic of Occam’s models is that the resulting subsurface conductivity distribution is a continuous function. This is very convenient from a mathematical point of view . Due to the roughness penalizer inclusion of Constable et al. , the developed model tends to smear discontinuities, an undesirable feature in most situations. In an effort to generate models with sharp boundaries, several techniques have been applied to this problem: piecewise continuous formulations as in ref. [4,5], l1 norm minimization [6,7,8], total variation penalizers in ref.  and , among others. An important requirement for the regularizer is to allow the recovery of edges and smooth the homogeneous parts. As is well known, total variation  is now the standard approach to meet this requirement. TV consists on applying a convex regularizer, ensuring existence and uniqueness of a solution. In the case of regular functions, the total variation of u is the norm of the gradient, i.e., .
There are two basic methods for MT transfer function estimation from field data: the single-site method and the remote reference method. In the first method, one or two biased estimates at each frequency are obtained for each measuring direction . In the second method, only one unbiased estimate can be obtained at one site simultaneously using a magnetic field measurement at a remote site [13,14]. The latter can be used effectively only for uncorrelated noise between local and remote sites. The remote reference method is a standard procedure in MT data acquisition because is the only way to remove bias. The single-site method is still used when a second instrument is not available or is faulty. The MT linear relationship of the impedance tensor is given by
Results are presented for the MT problem, but can also be applied to other types of electromagnetic soundings.
Tikhonov’s regularization method is the standard technique applied to obtain models of the subsurface conductivity distribution from electric or electromagnetic measurements. The model proposed for the conductivity distribution (m) is obtained by minimizing the functional
In an effort to generate models with sharp boundaries, several techniques have been applied to this problem: piecewise continuous formulations as in ref. [4,5], norm minimization [6,7,8], total variation penalizers in ref. [9,10], among others. An important requirement for the regularizer is to allow the recovery of edges and smooth the homogeneous parts. As is well known, Total Variation  is now the standard approach to meet this requirement. Considering in (2), a more realistic model is recovered, but the resulting model can be severely affected if the Gaussian’s assumption of data noise is not accomplished. The algorithm solving (2) with and norm for fitting will be called TVL2. The model developed by solving TVL2 is observed in Figure 1a for the synthetic data presented in Figure 1b, obtained by solving the forward equations for a four-layered model.
For the case of impulse noise, a preferred model is TVL1, replacing the norm in the data fitness term with , and considering TV as regularizer :
Nikolova  noticed that the solutions of the TVL1 model substantially deviate from both the data-acquisition and the prior models. Gu et al.  proposed to combine TV regularization with nonconvex smoothly clipped absolute deviation penalty for data fitting called TVSCAD. A difference of convex functions is adopted to solve the non-convex TVSCAD.
Cai et al.  proposed a two-phase approach to restore images corrupted by blur and impulse noise called CTVL1. The first phase is to detect possibly corrupted pixels, then in the second phase, they restore images using only the noncorrupted pixels. The adaptive correction procedure obtains an initial estimator of the image by TVL1 and then apply a corrective step by blurring the initially estimated image . The connections among TVSCAD and CTVL1 are observed by Gu et al. . They notice that the main difference is in the function considered in the correction step of CTVL1. They claim to obtain comparable results to TVL1. An important issue of TVSCAD is the selection of SCAD function parameters. As observed by Gu et al., they are problem-dependent, but the best choices are very difficult to determine.
Recently, Zhang et al.  proposed a nonconvex approach to image restoration contaminated with impulse noise. They consider a nonconvex data fitting term and TV for regularization operator. The nonconvex functions used are an exponential type and the Geman function . They proposed a proximal linearized minimization algorithm for the nonconvex part and the alternate direction method of multipliers (ADMM) for the total variation part. They sow convergence for the case of a concave function .
In this paper, we also consider applying a nonconvex function to the impulse noise treatment. The nonconvex part is solved using a half-quadratic minimization algorithm proposed by Geman and Yang . The total variation regularized is also considered with the ADMM method. The algorithm is easy to implement and the convergence properties were established by Charbonnier et al. for a nonconvex function satisfying the edge-preserving properties . The application is to the MT nonlinear inverse problem. As mentioned earlier, other algorithms (TVSCAD, CTVL1) present some difficulties in parameter selection, and others like  require a concave function to converge.
A smooth, nonconvex regularization algorithm was proposed by Charbonnier et al. , considering the following operator:
Charbonnier et al. apply an iterative algorithm developed by Geman and Yang  called “half-quadratic regularization” for the minimization of a modified version of the original functional (2) with P = PEPR. This algorithm is basically applying a convex conjugate to replace the original nonconvex function with a new one defined over an extended domain. Charbonnier et al. showed that for a potential function satisfying several conditions, the sequence generated by this algorithm is convergent.
In this paper, we propose to apply the operator to the fitness term in (2), solving
The algorithm is aimed to minimize the functional
The solution with respect to v is given by a shrinkage operator,
The Geman–McClure function is applied here, but several other functions can also be considered, as observed in ref. .
Parameter selection is not as complicated as in other methods; and can be selected by some of the regularization parameter selection methods available . The penalty parameter γ is obtained by trial and error, and as noticed by Getreuer , the algorithm is not so sensitive to the fair value of γ. Parameter β is selected as a scaling term according to the model’s dynamic range.
The algorithm is presented in Table 1, minimization is stopped when the model does not change in some min value.
|update with (10) and (11)|
It should be noticed that, as the model evolves from a homogeneous single layer to a multilayered system, the residuals will be large in the first steps. Because of that the potential function should be adapted, to consider those large residuals at the first steps. To implement that parameter β is initialized to a value so that a quadratic error fitness function is approximated for a few (ten) steps. Then, the desired value of β is applied to obtain a non-convex fitness function. The algorithm parameters are observed in Table 2.
TVL1 and TVL2 were stopped when was less than 0.04 and 0.01, respectively. All algorithms were asked to stop when min was 0.01, but TVL1 would not converge for the parameters shown in Table 2. TVL1 is implemented using ADMM for both, TV and L1 minimization, it uses other parameter for the L1 part, its value was also 10.
For TVL2, both regularization parameters were the same because the standard deviation is used as a weight for the measurements; a standard deviation of 0.05 for magnitude and 0.5 degrees for phase were considered in the example. The standard deviation affects the selection of μ, which is the reason for a smaller μ in that algorithm. The developed model by TVNC2 and the corresponding fitness are shown in Figure 4(a and b). It can be observed that the recovery of the model improves significantly over the previous cases. Notice that the resistivity of the third layer is underestimated. This is not so much a drawback of the algorithm but an intrinsic feature of the MT method, which detects much better good conductors than good resistors. This is because in good resistors, the induced currents are smaller than in good conductors. Consequently, the former does not alter significantly the electric and magnetic fields on the surface.
Figure 5 presents a model developed by TVNC2 for field data of site PC5008 in the COPROD2 data set . The data are the “uncorrected” original data provided by Jones . The field data are observed in Figure 6(a and b). A large deviation can be observed in the first two lectures.
It should be noticed that TVL2 does not converge for this uncorrected data. A correction was made by Jones , obtaining the data observed in Figure 7(b and c). The model developed by TVL2 is observed in Figure 7a. Notice the similarity with the one obtained by the TVNC2 algorithm when using the uncorrected data.
The MT sounding data are usually contaminated with large errors. The robust treatment of MT data has been generally applied in the transfer function components determination. In this paper, the case of large errors on the resulting apparent resistivity and phase data is considered.
A comparison is made among the usual least-squares inversion, a robust technique based on the L1 norm applied on the fitness term, and a nonconvex operator proposed in the image restoration some time ago for image noise filtering.
The application of a nonconvex operator for the fitness term and TV for regularization for impulse noise treatment on MT data is proposed. The inverse problem is solved by using a half-quadratic procedure for the nonconvex operator, along with a split Bregman procedure, ensuring convergence. The proposed algorithm parameters are easy to select.
Results are presented for data obtained from a synthetic model, observing the capability of the proposed method in model recovering when large deviations are inserted on the data. Application to field COPROD2 data is presented too, assessing the capability of the algorithm to recover a model of raw, uncorrected data.
The authors would like to thank P. J. Savage of PanCanadian Petroleum Limited of Calgary, Alberta, for making the raw (uncorrected) data available to the Geological Survey of Canada.
 Whittall KP, Oldenburg DW. Inversion of magnetotelluric data for a one-dimensional conductivity. Tulsa: Society of Exploration Geophysicists; 1992 Jan 1. Search in Google Scholar
 Constable SC, Parker RL, Constable CG. Occam’s inversion: A practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics. 1987 Mar;52(3):289–300. Search in Google Scholar
 Gómez-Treviño E. Nonlinear integral equations for electromagnetic inverse problems. Geophysics. 1987 Sep;52(9):1297–302. Search in Google Scholar
 Hidalgo H, Marroquin JL, Gómez-Treviño E. Piecewise smooth models for electromagnetic inverse problems. IEEE Trans Geosci Remote Sens. 1998 Mar;36(2):556–61. Search in Google Scholar
 Varentsov IM. A general approach to the magnetotelluric data inversion in a piecewise-continuous medium. Izv. Phys Solid Earth C/C Fiz Zemli-Rossiiskaia Acad. Nauk. 2002 Nov 1;38(11):913–34. Search in Google Scholar
 Esparza FJ, Gómez-Treviño E. Electromagnetic sounding in the resistive limit and the Backus-Gilbert method for estimating averages. Geoexploration. 1987 Dec 1;24(6):441–54. Search in Google Scholar
 Esparza FJ, Gómez-Treviño E. 1-D inversion of resistivity and induced polarization data for the least number of layers. Geophysics. 1997 Nov 1;62(6):1724–9. Search in Google Scholar
 Hidalgo H, Gómez-Treviño E, Pérez-Flores MA. Linear programs for the reconstruction of 2d images from geophysical electromagnetic measurements. Subsurface Sens Technol Appl. 2004 Apr 1;5(2):79–96. Search in Google Scholar
 Dobson DC. Recovery of blocky images in electrical impedance tomography. In Inverse problems in medical imaging and nondestructive testing. Vienna: Springer; 1997. p. 43–64. Search in Google Scholar
 Vogel CR, Oman ME. Fast, robust total variation-based reconstruction of noisy, blurred images. IEEE Trans Image Process. 1998 Jun;7(6):813–24. Search in Google Scholar
 Rudin LI, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms. Phys D Nonlinear Phenom. 1992 Nov 1;60(1–4):259–68. Search in Google Scholar
 Sims WE, Bostick Jr, FX, Smith HW. The estimation of magnetotelluric impedance tensor elements from measured data. Geophysics. 1971 Oct;36(5):938–42. Search in Google Scholar
 Gamble TD, Goubau WM, Clarke J. Magnetotellurics with a remote magnetic reference. Geophysics. 1979 Jan;44(1):53–68. Search in Google Scholar
 Clarke J, Gamble TD, Goubau WM, Koch RH, Miracky R. Remote-reference magnetotellurics: equipment and procedures. Geophys Prospecting. 1983 Jan 1;31(1):149–70. Search in Google Scholar
 Trad DO, Travassos JM. Wavelet filtering of magnetotelluric data. Geophysics. 2000 Mar;65(2):482–91. Search in Google Scholar
 Bai M, Zhang X, Shao Q. Adaptive correction procedure for TVL1 image deblurring under impulse noise. Inverse Probl. 2016 Jun 16;32(8):085004. Search in Google Scholar
 Nikolova M. A variational approach to remove outliers and impulse noise. J Math Imaging Vis. 2004 Jan 1;20(1–2):99–120. Search in Google Scholar
 Gu G, Jiang S, Yang J. A TVSCAD approach for image deblurring with impulsive noise. Inverse Probl. 2017 Nov 10;33(12):125008. Search in Google Scholar
 Nikolova M. Model distortions in Bayesian MAP reconstruction. Inverse Probl & Imaging. 2007;1(2):399. Search in Google Scholar
 Cai JF, Chan RH, Nikolova M. Two-phase approach for deblurring images corrupted by impulse plus Gaussian noise. Inverse Probl Imaging. 2008;2(2):187. Search in Google Scholar
 Zhang X, Bai M, Ng MK. Nonconvex-TV based image restoration with impulse noise removal. SIAM J Imaging Sci. 2017;10(3):1627–67. Search in Google Scholar
 Geman D, Yang C. Nonlinear image recovery with half-quadratic regularization. IEEE Trans Image Process. 1995 Jul;4(7):932–46. Search in Google Scholar
 Charbonnier P, Blanc-Féraud L, Aubert G, Barlaud M. Deterministic edge-preserving regularization in computed imaging. IEEE Trans Image Process. 1997 Feb;6(2):298–311. Search in Google Scholar
 Hidalgo H, Gómez-Treviño E, Marroquin JL, Esparza FJ. Piecewise continuous models for resistivity soundings. IEEE Trans Geosci Remote Sens. 2001 Dec;39(12):2725–8. Search in Google Scholar
 Geman S, McClure D. Statistical methods for tomographic image reconstruction. Bull Int Stat Inst. 1987;52:5–21. In Proceedings of the 46th Session of the International Statistical Institute. Search in Google Scholar
 Vogel CR. Computational methods for inverse problems. Philadelphia: Society for Industrial and Applied Mathematics; 2002 Jan 1. Search in Google Scholar
 Getreuer P. Rudin–Osher–Fatemi total variation denoising using split Bregman. Image Process Line. 2012 May 19;2:74–95. Search in Google Scholar
 Jones AG. The COPROD2 dataset: Tectonic setting, recorded MT data, and comparison of models. J Geomag Geoelec. 1993 Sep 20;45(9):933–55. Search in Google Scholar