Abstract
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 L1norm 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 halfquadratic 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.
1 Introduction
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 onedimensional (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 [1]: linearization, asymptotic behavior, and exact methods. Constable et al. [2] 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 [3]. Due to the roughness penalizer inclusion of Constable et al. [2], 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], l
_{1} norm minimization [6,7,8], total variation penalizers in ref. [9] and [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 [11] 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
There are two basic methods for MT transfer function estimation from field data: the singlesite method and the remote reference method. In the first method, one or two biased estimates at each frequency are obtained for each measuring direction [12]. 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 singlesite 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
where E is the local electric field, Z is the MT transfer function tensor, H is the magnetic field, and N is the noise. The impedance Z is either estimated by the single source or remote reference method. The noise N accounts for higher dimensional sources and short duration nonstationary and instrumental noise. N vanishes for the strict zero wavenumber or plane wave model. Once Z is estimated, the residuals can be computed by substituting their values in equation (1). Residuals are given by the difference between the observed E field and the predicted E field, obtained by substituting the estimated Z in relation (1). The use of classical spectral analysis together with leastsquares regression is warranted if data follow a stationary and Gaussian model. In this case, it is also common to assume that residuals follow a multivariate normal probability distribution. If this is satisfied, the maximum likelihood estimates are obtained. In practice, most data display gross departures, or outliers, from the simple regression model [15]. The main causes are geomagnetic phenomena, thunderstorms, anthropogenic contribution, and instrumental problems. Outliers usually appear as a fraction of the useful observations, with distinct characteristics from the rest of the sample. The phenomenon is independent of the structure of the bulk of geomagnetic field observations. In this paper, we do not consider the robust estimation of the MT transfer functions, the treatment of MT data contaminated with large variations is the main purpose here.
Results are presented for the MT problem, but can also be applied to other types of electromagnetic soundings.
2 MT inversion
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
with
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],
In Figure 2a, a model is obtained by the same algorithm solving (2) but inserting a large variation on the first amplitude data. Notice the effect on the resulting model.
For the case of impulse noise, a preferred model is TVL1, replacing the
where
Nikolova [19] noticed that the solutions of the TVL1 model substantially deviate from both the dataacquisition and the prior models. Gu et al. [18] 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 nonconvex TVSCAD.
Cai et al. [20] proposed a twophase 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 [20]. The connections among TVSCAD and CTVL1 are observed by Gu et al. [18]. 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 problemdependent, but the best choices are very difficult to determine.
Recently, Zhang et al. [21] 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
In this paper, we also consider applying a nonconvex function to the impulse noise treatment. The nonconvex part is solved using a halfquadratic minimization algorithm proposed by Geman and Yang [22]. 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 edgepreserving properties [23]. 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 [21] require a concave function to converge.
2.1 Nonconvex function applied to noise treatment
A smooth, nonconvex regularization algorithm was proposed by Charbonnier et al. [23], considering the following operator:
where
with β > 0. This function does not satisfy assumptions required by Zhang et al. for convergence [21] because it is not concave as
Charbonnier et al. apply an iterative algorithm developed by Geman and Yang [22] called “halfquadratic regularization” for the minimization of a modified version of the original functional (2) with P = P _{EPR}. 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
with ϕ(⋅) applied to the magnitude
The algorithm is aimed to minimize the functional
with two auxiliary variables
with
The solution with respect to v is given by a shrinkage operator,
with
With ϕ′ the derivative of ϕ.
The Geman–McClure function is applied here, but several other functions can also be considered, as observed in ref. [23].
Parameter selection is not as complicated as in other methods;
The algorithm is presented in Table 1, minimization is stopped when the model does not change in some min value.
Algorithm 
Input:

Initialization:

repeat 
evaluate




update


until

Output

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 nonconvex fitness function. The algorithm parameters are observed in Table 2.
Algorithm 


γ  RMS 

TVL2  10  10  10  8.4 
TVL1  100  10  10  8.22 
TVNC2  100  10  10  8.63 
TVL1 and TVL2 were stopped when
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.
2.2 Application to field data
Figure 5 presents a model developed by TVNC2 for field data of site PC5008 in the COPROD2 data set [28]. The data are the “uncorrected” original data provided by Jones [28]. 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 [28], 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.
3 Conclusion
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 leastsquares 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 halfquadratic 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.
Acknowledgments
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.

Author contributions: HHS developed and implemented the algorithms, and EGT assisted with the analysis and interpretation of results.
References
[1] Whittall KP, Oldenburg DW. Inversion of magnetotelluric data for a onedimensional conductivity. Tulsa: Society of Exploration Geophysicists; 1992 Jan 1.10.1190/1.9781560802419Search in Google Scholar
[2] 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.10.1190/1.1442303Search in Google Scholar
[3] GómezTreviño E. Nonlinear integral equations for electromagnetic inverse problems. Geophysics. 1987 Sep;52(9):1297–302.10.1190/1.1442390Search in Google Scholar
[4] Hidalgo H, Marroquin JL, GómezTreviño E. Piecewise smooth models for electromagnetic inverse problems. IEEE Trans Geosci Remote Sens. 1998 Mar;36(2):556–61.10.1109/36.662738Search in Google Scholar
[5] Varentsov IM. A general approach to the magnetotelluric data inversion in a piecewisecontinuous medium. Izv. Phys Solid Earth C/C Fiz ZemliRossiiskaia Acad. Nauk. 2002 Nov 1;38(11):913–34.Search in Google Scholar
[6] Esparza FJ, GómezTreviño E. Electromagnetic sounding in the resistive limit and the BackusGilbert method for estimating averages. Geoexploration. 1987 Dec 1;24(6):441–54.10.1016/00167142(87)900135Search in Google Scholar
[7] Esparza FJ, GómezTreviño E. 1D inversion of resistivity and induced polarization data for the least number of layers. Geophysics. 1997 Nov 1;62(6):1724–9.10.1190/1.1444272Search in Google Scholar
[8] Hidalgo H, GómezTreviño E, PérezFlores MA. Linear programs for the reconstruction of 2d images from geophysical electromagnetic measurements. Subsurface Sens Technol Appl. 2004 Apr 1;5(2):79–96.10.1023/B:SSTA.0000034544.25382.8fSearch in Google Scholar
[9] 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.10.1007/9783709165218_5Search in Google Scholar
[10] Vogel CR, Oman ME. Fast, robust total variationbased reconstruction of noisy, blurred images. IEEE Trans Image Process. 1998 Jun;7(6):813–24.10.1109/83.679423Search in Google Scholar
[11] 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.10.1016/01672789(92)90242FSearch in Google Scholar
[12] Sims WE, Bostick Jr, FX, Smith HW. The estimation of magnetotelluric impedance tensor elements from measured data. Geophysics. 1971 Oct;36(5):938–42.10.1190/1.1440225Search in Google Scholar
[13] Gamble TD, Goubau WM, Clarke J. Magnetotellurics with a remote magnetic reference. Geophysics. 1979 Jan;44(1):53–68.10.1190/1.1440923Search in Google Scholar
[14] Clarke J, Gamble TD, Goubau WM, Koch RH, Miracky R. Remotereference magnetotellurics: equipment and procedures. Geophys Prospecting. 1983 Jan 1;31(1):149–70.10.1111/j.13652478.1983.tb01047.xSearch in Google Scholar
[15] Trad DO, Travassos JM. Wavelet filtering of magnetotelluric data. Geophysics. 2000 Mar;65(2):482–91.10.1190/1.1444742Search in Google Scholar
[16] Bai M, Zhang X, Shao Q. Adaptive correction procedure for TVL1 image deblurring under impulse noise. Inverse Probl. 2016 Jun 16;32(8):085004.10.1088/02665611/32/8/085004Search in Google Scholar
[17] Nikolova M. A variational approach to remove outliers and impulse noise. J Math Imaging Vis. 2004 Jan 1;20(1–2):99–120.10.1023/B:JMIV.0000011920.58935.9cSearch in Google Scholar
[18] Gu G, Jiang S, Yang J. A TVSCAD approach for image deblurring with impulsive noise. Inverse Probl. 2017 Nov 10;33(12):125008.10.1088/13616420/aa9383Search in Google Scholar
[19] Nikolova M. Model distortions in Bayesian MAP reconstruction. Inverse Probl & Imaging. 2007;1(2):399.10.3934/ipi.2007.1.399Search in Google Scholar
[20] Cai JF, Chan RH, Nikolova M. Twophase approach for deblurring images corrupted by impulse plus Gaussian noise. Inverse Probl Imaging. 2008;2(2):187.10.3934/ipi.2008.2.187Search in Google Scholar
[21] Zhang X, Bai M, Ng MK. NonconvexTV based image restoration with impulse noise removal. SIAM J Imaging Sci. 2017;10(3):1627–67.10.1137/16M1076034Search in Google Scholar
[22] Geman D, Yang C. Nonlinear image recovery with halfquadratic regularization. IEEE Trans Image Process. 1995 Jul;4(7):932–46.10.1109/83.392335Search in Google Scholar PubMed
[23] Charbonnier P, BlancFéraud L, Aubert G, Barlaud M. Deterministic edgepreserving regularization in computed imaging. IEEE Trans Image Process. 1997 Feb;6(2):298–311.10.1109/83.551699Search in Google Scholar PubMed
[24] Hidalgo H, GómezTreviño E, Marroquin JL, Esparza FJ. Piecewise continuous models for resistivity soundings. IEEE Trans Geosci Remote Sens. 2001 Dec;39(12):2725–8.10.1109/36.975007Search in Google Scholar
[25] 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
[26] Vogel CR. Computational methods for inverse problems. Philadelphia: Society for Industrial and Applied Mathematics; 2002 Jan 1.10.1137/1.9780898717570Search in Google Scholar
[27] Getreuer P. Rudin–Osher–Fatemi total variation denoising using split Bregman. Image Process Line. 2012 May 19;2:74–95.10.5201/ipol.2012.gtvdSearch in Google Scholar
[28] Jones AG. The COPROD2 dataset: Tectonic setting, recorded MT data, and comparison of models. J Geomag Geoelec. 1993 Sep 20;45(9):933–55.10.5636/jgg.45.933Search in Google Scholar
© 2021 Hugo HidalgoSilva and Enrique GómezTreviño, published by De Gruyter
This work is licensed under the Creative Commons Attribution 4.0 International License.