Robust Regression Estimators in Magnetotelluric Data Processing: Performance & Evaluation

1. Introduction

As a conventional electromagnetic (EM) geophysical method, magnetotelluric (MT) is widely used to investigate subsurface geo-electrical structures, including crustal and/or lithospheric architecture, deep geodynamic processes, resource exploration targets, and orogenic deformation systems. In particular, subsurface electrical structures provide important constraints on metallogenic systems because conductive anomalies are commonly associated with magma reservoirs, fluid migration pathways, partial melts, and fault-controlled weak zones that play critical roles in mineral deposits-forming processes. Mounts of classical MT studies associated with mineral deposits mechanism were published in the last decades [1,2,3,4,5,6], one of these studies, Heinson et al. [3,4] show an impressive resistivity model in the Olympic Dam district, Australia, demonstrating that MT-derived electrical structures can provide critical information for understanding of mineralization. However, it is generally challenging to observe ideal EM fields data in mining areas, that are usually dominated by complex and strong industrial background noises. Consequently, the reliability of geo-electrical model and the following geological interpretations are limited by strongly distorted MT transfer functions.

MT transfer functions are derived from orthogonal naturally time-varying EM fields. Partial transfer functions and power spectral density matrices are computed from discrete spectra obtained using short-time Fourier transforms within individual time windows. Under the assumption of noise-free conditions, the frequency-dependent impedance tensor $\mathit{Z}\left(\omega \right)$ is defined as [7–9],

$\begin{array}{c}\left(\begin{array}{c}{E}_x\left(\omega \right)\\ E_y\left(\omega \right)\end{array}\right)=\left(\begin{array}{cc}{Z}_{xx}\left(\omega \right)& Z_{xy}\left(\omega \right)\\ Z_{yx}\left(\omega \right)& Z_{yy}\left(\omega \right)\end{array}\right)\left(\begin{array}{c}{H}_x\left(\omega \right)\\ H_y\left(\omega \right)\end{array}\right)\left(1\right)\end{array}$

[2] where $H_x$ and $H_y$ denote the magnetic field components, and $E_x\mathrm{ }$ and $E_y$ denote the electric field components. In practice, EM signals are generally contaminated with background noises and a linear regression system is typically expressed as,

$\begin{array}{c}r\left(\omega \right)=E\left(\omega \right)-Z\left(\omega \right)H\left(\omega \right)\left(2\right)\end{array}$

where $\mathit{r}\left(\omega \right)$ is residual of the dependent variable (i.e., $\mathit{E}\left(\omega \right)$ ). The ordinary least-squares (LS) method minimizes the Euclidean norm of the residuals and is therefore sensitive to outliers when the Gaussian noise assumption is violated. Since the early development of MT, various statistical approaches have been proposed for robust impedance estimation [10–17] .

On the other hand, outliers associated with predictor variables, i.e., magnetic field data, commonly referred to as leverage points [18,19], cannot be effectively handled by conventional robust approaches. The remote reference technique was proposed to suppress local magnetic noises that are uncorrelated with the noise at reference stations [20,21,22]. In cases where a remote reference station is available, both robust and remote reference methods are employed in combination [12,15,23,24]. However, if both local and reference channels are contaminated by correlated noises, the results from the remote reference method will be distorted [25,26]. Different from the classical remote reference approach, in which the cross-power spectral density between the local and reference magnetic fields is introduced into the spectral density matrix, other algorithms have been proposed to suppress leverage points associated with noise in the local magnetic field. The two-stage bounded influence estimator is one of the effective methods used to remove correlated noise in local EM fields using clean remote stations [18,27]. Since the remote reference technique can be considered a two-input–multiple-output linear system, Usui et al. [17] proposed a robust remote reference estimator based on robust multivariate linear regression. Ogawa et al. [28] separated signals from noise-affected components by applying a frequency-domain independent component analysis to local EM data and reference magnetic data. A similar approach based on frequency-domain independent component analysis was proposed by Sato et al. [29]. It is worth noting that high-quality reference data are generally required by both conventional and derivative remote reference methods. Furthermore, when array MT stations and synchronous data are available, a robust multivariate errors-in-variables estimator that utilizes principal component analysis is proven to be a practical approach for extracting signals corresponding to the two polarizations of plane-wave MT fields [30,31,32]. Smirnov and Egbert [33] further extended the multiple-station technique [31] to be applicable to large incomplete arrays, performing practical deployment strategies such as the “rolling array”.

Pre-screening of data segments in the frequency domain is generally employed prior to robust estimation and/or remote reference methods to exclude severely noisy data. Different pre-selection strategies have been suggested in previous studies, and here we list some of the proposed criteria. To deal with low signal-to-noise ratios in the dead band, Egbert and Livelybrooks [14] selected only data segments with high electric-to-magnetic field multiple coherence. Signals meeting a required magnitude that is significantly higher than the self-noise (sensitivity) of magnetic sensors were selected as effective data [34]. In auroral zones, limited vertical magnetic field variations were considered as data contaminated by non-uniform sources [35]. In a more general case, two partial coherences calculated from orthogonal electric and magnetic components were suggested as criteria for removing highly noisy segments [16]. Moreover, physical properties, including spectral power densities, polarization directions, and coherences, were also evaluated to improve MT transfer functions [36], and a further automated pre-selection approach based on Mahalanobis distance and magnetic field constraints was developed by Platz and Weckmann [37]. A multi-criteria data sorting strategy was suggested in the study by Wang et al. [38].

In addition, various approaches have been employed in the time domain to enhance data quality, including time-series denoising using mathematical transforms [39,40], decomposition methods [41,42], digital signal filtering [43], autoregressive models [16], and signal reconstruction algorithms based on compressive sensing [44], as well as many other algorithms not listed in this paper. It should be acknowledged that artificial neural networks offer distinct advantages in time-series data cleaning [45,46]. Conventional artificial neural networks have been applied to MT data processing and perform well for typical synthetic noise [47,48,49,50,51,52]; however, misfits are inevitable when the data contain complex, persistent, or long-duration noise. Furthermore, supervised learning–based approaches typically require time-consuming manual identification and labeling [45], and the robustness of the trained models depends on relatively large datasets, which can be challenging to obtain. Unsupervised learning–based methods typically require careful selection of parameters to achieve optimal performance [53,54,55], which remains challenging in practical applications.

Although various approaches in both the time and frequency domains are generally employed to obtain final transfer functions, robust techniques combined with the remote reference method remain the most practical and widely adopted approaches for addressing the frequency-dependent linear regression relationships between orthogonal EM fields. In practice, the most widely used approach is M-estimation [56], which generalizes the maximum likelihood estimation framework. The M-estimation employs the median absolute deviation as a univariate scale estimator, which is independent of the overall data distribution [57]. An alternative approach is the S-estimator [58], which is based on robust scale estimators. Usui et al. [17] applied the S-estimator to their multiple-output system and illustrated the differences between the M- and S-estimators using a simple linear regression problem. To combine the advantages of both M- and S-estimation, the robustness and high breakdown point of S-estimation and the efficiency of M-estimation, MM-estimation was subsequently proposed by Yohai [59]. Theoretically, the performance of different estimators is expected to vary across typical statistical scenarios.

It is difficult to evaluate the contribution of each step (time- and/or frequency-dependent algorithms) to noise suppression from time-series data to transfer functions because of complex and varying EM noise. Although a combination of all possible approaches in both the time and frequency domains seems to be an optimal strategy, it is still meaningful to analyze the performance of a specified method. In this study, we developed a robust estimator based on the MM-estimation algorithm and applied the M-, S-, and MM-estimators to both synthetic and filed MT datasets to evaluate their performance. Our study suggests that a comprehensive comparative analysis of results obtained from different estimators combined with remote reference techniques should be conducted in practical applications, enhancing studies on mineralization mechanism.

2. Theory and Methods

2.1. MT Data Processing

Here, we briefly describe the basic procedures of MT data processing, while details are introduced in many previous studies [16,17,27,60]. Time-series data of observed orthogonal EM fields are generally checked and/or corrected using time-domain algorithms, including filtering, artificial neural networks, decomposition methods, etc. Pre-cleaned time-series data are then subdivided into overlapping data segments (time windows), and subsequent cascade decimation with low-pass filtering is generally performed. In general, analysis and processing proceed for each data segment as follows: (1) Pre-whitening is performed to remove long-period trends and mean values. (2) The segment is tapered using a Hanning or Hamming window to reduce spectral leakage. (3) Fast Fourier transforms (FFTs) are performed to obtain Fourier coefficients for all channels, and calibration correction is also applied. (4) A spectral density matrix consisting of auto- and cross-spectral values is calculated. (5) Frequency-domain pre-selection is carried out using coherence sorting and other methods.

As a result of the above steps, spectral density matrices from all segments at specified evaluation frequencies are prepared for the subsequent regression procedures.

2.2. Robust Estimation

M-estimation [56], developed within the maximum-likelihood framework, provides a robust generalization of the classical LS method. For the linear regression problem defined in Equations (1) and (2), M-estimation introduces a general loss function, ρ, and minimizes the following objective function,

$\begin{array}{c}{\widehat{\mathsf{\theta }}}_M=\underset{\widehat{\mathsf{\theta }}}{\min }{\displaystyle \sum _{i=1}^N}\rho \left({\mathit{r}}_i\right)\left(3\right)\end{array}$

where “M” denotes “maximum likelihood–type,” and N denotes the total number of observations, defined as the product of the number of observation equations (for either $E_x\left(\omega \right)$ or $E_y\left(\omega \right)$ ) and the number of time windows. And a scale estimate based on the median absolute deviation is typically used, which is largely independent of the overall dataset [57].

In practice, the residuals are standardized as shown in Equation (4) [61],

$\begin{array}{c}{\displaystyle \sum _{i=1}^N}\psi \left({\displaystyle \frac{{\mathit{r}}_i}{\widehat{\sigma }}}\right){\mathit{H}}_i=0\left(4\right)\end{array}$

where $\psi$ is the derivative of $\rho$ , and the Huber function [56] is commonly employed. $\widehat{\sigma }$ is a scale parameter,

$\begin{array}{c}\widehat{\sigma }={\displaystyle \frac{MAD}{0.6745}}={\displaystyle \frac{median\left|{\mathit{r}}_i-median\left({\mathit{r}}_i\right)\right|}{0.6745}}\left(5\right)\end{array}$

To solve the nonlinear optimization problem in Equation (4), the iteratively reweighted least squares (IRLS) algorithm is typically used. The core idea is to transform the objective function of robust estimation into a weighted least-squares problem, where the weights depend on residuals from the previous iteration. Specifically, the weight matrix at the $k$ -th iteration is defined as $W^{\left(k\right)}=diag\left\{\omega \left(r_1^{\left(k\right)}\right),\dots ,\omega \left(r_n^{\left(k\right)}\right)\right\}$ , where $\omega$ (⋅) is a weight function of standardized residuals,

$\begin{array}{c}{\displaystyle \sum _{i=1}^N}{\omega }_i\left({\displaystyle \frac{{\mathit{r}}_i}{\widehat{\sigma }}}\right){\mathit{H}}_i=0\left(6\right)\end{array}$

Consequently, the updated estimator in Equation (6) can be explicitly formulated as follows,

$\begin{array}{c}{\widehat{\mathit{Z}}}^{\left(k+1\right)}={\left({\mathit{H}}^T{W}^{\left(k\right)}\mathit{H}\right)}^{-1}{\mathit{H}}^T{W}^{\left(k\right)}E\left(7\right)\end{array}$

This iterative procedure continues until convergence, effectively attenuating the influence of outliers by assigning them lower weights according to their relative residuals.

In this study, we employed Tukey’s bisquare weight function (Equation (8)) [62] as the weighting scheme in the estimation procedure,

$\begin{array}{c}{\omega }_i\left({\mathit{u}}_i\right)=\left\{\begin{array}{c}{\left[1-{\left({\displaystyle \frac{{\mathit{u}}_i}{c}}\right)}^2\right]}^2,\left|{\mathit{u}}_i\right|\le c\\ 0,\left|{\mathit{u}}_i\right|>c\end{array}\right.\left(8\right)\end{array}$

where ${\mathit{u}}_i$ denotes the standardized residual, i.e., ${\mathit{u}}_i={\displaystyle \frac{{\mathit{r}}_i}{\widehat{\sigma }}}$ . Under the assumption of a normal distribution, the estimator achieves approximately 95% asymptotic efficiency when $\mathrm{c}=4.685$ .

To overcome the limitations of the median-based scale in M-estimation, S-estimation employs the residual standard deviation [57] and is formulated to minimize the scale of the residuals as shown in Equation (9) [61]:

$\begin{array}{c}{\widehat{\mathsf{\theta }}}_s=arg\underset{\widehat{\mathsf{\theta }}}{\min }{\widehat{\sigma }}_s\left\{{\mathit{r}}_1\left(\mathsf{\theta }\right),\dots ,{\mathit{r}}_n\left(\mathsf{\theta }\right)\right\}\left(9\right)\end{array}$

where ${\widehat{\sigma }}_s$ is obtained by solving Equation (10),

$\begin{array}{c}{\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N}\rho \left({\displaystyle \frac{{\mathit{r}}_i}{{\widehat{\sigma }}_s}}\right)=K\left(10\right)\end{array}$

where $K$ represents the expected value of the chosen loss function $\rho$ under the standard normal distribution. For high-breakdown-point S-estimation, a commonly used value is $K\approx 0.199$ . The initial estimate of ${\widehat{\sigma }}_s$ is obtained from Equation (5), and the subsequent iterative update formula is given as follows,

$\begin{array}{c}{{\widehat{\sigma }}_s}^{\left(k+1\right)}={{\widehat{\sigma }}_s}^{\left(k\right)}\sqrt{{\displaystyle \frac{1}{NK}}{\displaystyle \sum _{i=1}^N}\rho \left({\displaystyle \frac{{\mathit{r}}_i}{{{\widehat{\sigma }}_s}^{\left(k\right)}}}\right)}\left(11\right)\end{array}$

S-estimation [58] incorporates a robust scale estimator and attenuates the influence of outliers by minimizing a scale-based criterion of the residuals. Due to the non-convex nature of the S-estimation objective function, a multi-start strategy is employed, where each run is initialized from a random subset. The final solution is selected as the one yielding the minimum scale.

One key distinction between S-estimation and M-estimation is that S-estimation minimizes its objective function by iteratively updating the scale parameter ${\widehat{\mathsf{\sigma }}}_{\mathrm{s}}$ , until no further decrease is observed. Accordingly, Equation (4) can be reformulated into the form given in Equation (12),

$\begin{array}{c}{\displaystyle \sum _{i=1}^n}\psi \left({\displaystyle \frac{{\mathit{r}}_i}{{\widehat{\sigma }}_s}}\right){\mathit{H}}_i=0\left(12\right)\end{array}$

M-estimation exhibits good robustness against outliers but remains sensitive to high-leverage points. In contrast, S-estimation achieves extremely high robustness by maximizing the breakdown point through constrained scale estimation; however, this robustness comes at the cost of increased computational complexity. Therefore, S-estimation is typically used as the initial stage of MM-estimation, while the subsequent refinement is carried out using M-estimation. MM-estimation [59] combines the advantages of both approaches, achieving the high breakdown robustness of S-estimation while retaining the high statistical efficiency of M-estimation, thereby overcoming the limitations associated with using either method alone.

Specifically, MM-estimation begins by computing a high-breakdown scale parameter ${\widehat{\mathsf{\sigma }}}_{\mathrm{s}}$ (Equation (10)) using S-estimation. Based on this initial scale estimate, an M-estimation is subsequently constructed using the weight function in Equation (8), thereby ensuring high statistical efficiency. The final MM-estimator is then obtained by solving the following estimating equation,

$\begin{array}{c}{\displaystyle \sum _{i=1}^n}\psi \left({\displaystyle \frac{{\mathit{r}}_i}{{\widehat{\sigma }}_{MM}}}\right){\mathit{H}}_i=0\left(13\right)\end{array}$

In practice, the three approaches are typically implemented using an iterative strategy, in which specific scale definitions and estimating functions are applied for M- and S-estimation, while MM-estimation is achieved by combining the M-estimator with the scale obtained from S-estimation.

Figure 1 illustrates the algorithmic workflows of the three robust estimation methods. The detailed procedures are described as,

• Initialization: The impedance tensor ${\mathbf{Z}}_0$ is initialized using an ordinary least-squares ( O LS) method.
• Residual computation: For the current estimate ${\mathbf{Z}}_{\mathit{i}}$ , the residuals ${\mathit{r}}_i$ are calculated from the observed magnetic (H) and electric (E) field components.
• Scale estimation and residual updating:
  M-estimation: Compute the scale parameter $\widehat{\sigma }$ using the median absolute deviation (MAD), and update ${\mathit{r}}_i$ using the M-estimator with $\widehat{\sigma }$ .
  S-estimation: Initialize the scale parameter ${\widehat{\sigma }}_s$ using MAD, then update ${\widehat{\sigma }}_s$ according to Equation (10). The residuals ${\mathit{r}}_i$ are updated using the S-estimator with the updated ${\widehat{\sigma }}_s$ .
  MM-estimation: Update ${\mathit{r}}_i$ using the M-estimator together with the scale parameter ${\widehat{\sigma }}_s$ obtained from S-estimation.
• Impedance tensor update: The impedance tensor ${\mathbf{Z}}_{\mathit{i}}$ is updated based on the weighted residuals ${\mathit{r}}_i$.
• Iteration: Reiterate steps 2–4 until the number of iterations reaches a defined maximum value.

3. Implementation

3.1. Synthetic Experiments

3.1.1. Synthetic Data

To better demonstrate that the synthetic experiments are fully driven by artificially imposed noise, we prepared a clean time-series dataset by computing the electrical field response of a theoretical model to observed magnetic field data. Raw MT time-series data were sampled at 0.5 Hz with a total duration of approximately 10 days using a LEMI-423 instrument. The measured magnetic flux density $\mathit{B}\left(\mathrm{t}\right)$ was first converted into magnetic field intensity $\mathit{H}\left(\mathrm{t}\right)$ using

$\begin{array}{c}B\left(t\right)={\mu }_{0}H\left(t\right)\left(14\right)\end{array}$

where ${\mu }_0$ is the magnetic permeability of free space. In the frequency domain, the horizontal electric field components were computed using the impedance tensor,

$\begin{array}{c}{E}_x\left(\omega \right)=Z_{xy}\left(\omega \right)H_y\left(\omega \right)\left(15\right)\end{array}$

$\begin{array}{c}{E}_y\left(\omega \right)=Z_{yx}\left(\omega \right)H_x\left(\omega \right)\left(16\right)\end{array}$

We defined a homogeneous half-space model, and the impedance is defined as,

$\begin{array}{c}Z\left(\omega \right)=\sqrt{{\displaystyle \frac{i\omega {\mu }_0}{\sigma }}}\left(17\right)\end{array}$

where $\sigma$ is electrical conductivity.

To mitigate spectral bias caused by long time series, the data were processed using a segmented windowing approach. Each segment was transformed into the frequency domain via FFT, multiplied by the impedance, and then transformed back using inverse FFT. The final time series were reconstructed using weighted overlap-add, where the weights are derived from the window functions.

The theoretical magnetic fields were then reconstructed from the synthetic electric field. A stabilized inverse formulation was adopted,

$\begin{array}{c}H\left(\omega \right)={\displaystyle \frac{\mathit{E}\left(\omega \right){\mathit{Z}}^{*}\left(\omega \right)}{{\left|\mathit{Z}\left(\omega \right)\right|}^2+\lambda }}\left(18\right)\end{array}$

where λ is a small regularization parameter to ensure numerical stability.

The vertical magnetic field was obtained using tipper transfer functions,

$\begin{array}{c}{H}_z\left(\omega \right)=T_x\left(\omega \right)H_x\left(\omega \right)+T_y\left(\omega \right)H_y\left(\omega \right)\left(19\right)\end{array}$

and both horizontal and vertical magnetic fields were subsequently converted back to magnetic flux density.

3.1.2. Robust Estimations

To evaluate the robustness of the proposed estimation methods, different types of noise, including Gaussian, square, peak, and sawtooth noise, were added to clean synthetic MT time-series data. The amplitude–frequency characteristics of the noise were controlled using the signal-to-noise ratio (SNR),

$\begin{array}{c}SNR=20{log}_{10}\left({\displaystyle \frac{{\sigma }_{signal}}{{\sigma }_{noise}}}\right)\left(20\right)\end{array}$

To simulate realistic non-stationary EM contaminations, noise was superimposed on randomly selected subsets of samples. Table 1 shows the different parameter settings for the four types of noise added to the time series, which were subsequently normalized using a unified SNR control scheme. Specifically, 30%, 50%, and 70% of sample points were randomly selected using a fixed random seed to ensure reproducibility, and noise was added to all channels. This strategy preserved the temporal correlation of noise across different channels and more realistically simulated EM interference in practical environments. Two SNR levels were set to 30 and 5, respectively.

Figure 2 presents comparisons between clean and noise-contaminated time series with an SNR of 30 dB. Regular preprocessing procedures mentioned above were conducted to both synthetic and observed datasets in this study. The mean value of each channel was removed to suppress DC components and reduce low-frequency spectral leakage. Seven levels of decimation were employed. At each level, the sampling frequency was reduced by a factor of 8. The evaluation frequencies were logarithmically distributed with seven frequencies per octave to ensure adequate spectral resolution across the investigated bandwidth. The decimated time series were segmented into overlapping windows of 512 samples with 25% overlap. A Hamming window was applied prior to Fourier transform to suppress spectral leakage and reduce sidelobe effects. Fourier transforms were performed on each windowed segment to obtain complex spectra for subsequent impedance tensor estimation. Frequency-domain quality control was primarily based on correlation analysis between the electric and magnetic field components.

Figure 3 presents the MT sounding results for three synthetic noise contamination cases under an SNR of 30 dB, where Gaussian, square-wave, peak, and sawtooth noise were superimposed onto the original clean time series. The contaminated proportions of the time series were set to 30%, 50%, and 70%, respectively. Overall, as the noise contamination ratio increases from 30% to 70%, the apparent resistivity and phase curves obtained using the three robust estimation methods remain relatively stable and consistent. Only the S-estimation exhibits a slight apparent resistivity deviation in the high-frequency band, while no significant distortion is observed in the phase results. These observations indicate that, under an SNR of 30 dB, all three robust estimation methods can effectively suppress the influence of mixed noise.

In comparison, the S-estimation exhibits relatively more pronounced scatter at several frequency points, followed by the M-estimation, whereas the MM-estimation shows the fewest outliers. It is noteworthy that, although the noise contamination ratio increases continuously from 30% to 70%, the phase curves do not exhibit obvious collapse or large-scale distortion. This indicates that all three robust estimation algorithms can effectively reduce the weights of anomalous spectral windows and non-stationary noise segments under low-noise conditions, thereby suppressing the influence of transient interference and intermittent noise on impedance tensor estimation. Among them, the MM-estimation demonstrates comparatively superior robustness.

Figure 4 presents the final weight distributions of the impedance tensors at ~181 s using different robust estimation methods. The weights were computed based on the linear regression models in Equations (1) and (2), and were iteratively updated during the optimization process of each robust estimator. For the $Z_{xy}$ and $Z_{yx}$ components, a common characteristic across all three cases is that, as the mixed-noise contamination ratio increases, all three methods are still able to maintain most high-weight points clustered within relatively stable magnitude and phase regions. However, for components with phases around 45$°$ and $-$135$°$, the MM-estimation assigns the largest number of high-weight points. Combined with the previously presented apparent resistivity and phase results (Figure 3), these observations demonstrate that the MM-estimation exhibits superior robustness and noise resistance under high mixed-noise contamination conditions.

To compare the robustness of the three methods more clearly, the noise intensity was further increased by setting the SNR to 5 dB, while keeping all other conditions unchanged. Figure 5 presents a comparison between clean and noise-contaminated time series with an SNR of 5 dB.

Figure 6 presents the MT sounding results under the condition of SNR = 5 dB for different levels of mixed-noise contamination. Panels (a-c) correspond to 30%, 50%, and 70% of the time series being contaminated by a combination of Gaussian, square-wave, peak, and sawtooth noise, respectively. Compared with the results obtained at SNR = 30 dB, the apparent resistivity curves under the SNR of 5 dB exhibit more fluctuations and scattered points for all three robust estimation methods, particularly in the high-frequency band. This indicates that, as the background noise level increases, the influence of mixed noise on impedance estimation becomes increasingly significant. Nevertheless, despite the intensified noise contamination, all three methods are still able to preserve the overall structural characteristics of the apparent resistivity and phase curves, indicating that the robust estimation procedures retain the ability to suppress anomalous data under strong-noise conditions. Specifically, as the noise contamination ratio increases, the apparent resistivity curves obtained using the MM-estimation remain smoother than those produced by the M- and S-estimation methods. It is noteworthy that when the contamination ratio increases from 30% to 70%, the short-period apparent resistivity responses are affected much more strongly than the mid- and long-period responses. This is likely because short-period responses are more sensitive to high-frequency transient noise, while the peak and square-wave components of the mixed noise introduce stronger interference in local spectral estimation. Under such low-SNR conditions, the MM-estimation demonstrates relatively better overall stability and stronger resistance to outliers, indicating that this method possesses strong robustness under severe mixed-noise contamination.

Figure 7 presents the weight distributions of the impedance tensors at ~181 s under a noisy condition of SNR = 5 dB. Compared with the results at SNR = 30 dB (Figure 4), the scatter distributions in Figure 7 exhibit stronger dispersion, with a substantial increase in peripheral low-weight points, indicating enhanced instability of impedance estimation under stronger mixed-noise interference. Nevertheless, all three robust estimation methods still maintain most high-weight points concentrated within relatively stable regions, demonstrating their ability to identify anomalous spectral windows and suppress noise contamination. Among them, the MM-estimation exhibits the most pronounced clustering behavior. Its high-weight points are mainly concentrated near the stable phase direction, whereas the peripheral anomalous scatter points are generally assigned relatively low weights, resulting in a clearer weight-separation pattern. Combined with the results shown in Figure 4, it can be concluded that the MM-estimation achieves a better balance between anomalous spectral-window rejection and valid spectral-window preservation. Because the MM-estimation combines the high breakdown point of the S-estimation with the high statistical efficiency of the M-estimation, it can more effectively distinguish valid subsurface responses from noise-dominated anomalous spectral windows under strong mixed-noise conditions, thereby maintaining more stable impedance estimation results. These observations further demonstrate that the MM-estimation exhibits superior robustness and noise suppression capability under high mixed-noise contamination conditions.

3.2. Application to Observed Data

To further evaluate performance of the robust estimation methods, a remote reference technique was introduced in the processing of field data. The experimental data are from a northwest – southeast oriented broadband MT profile deployed in the northeastern china [63]. The profile starts from Songliao block and goes into the Xing’an block, regions mineral deposits and oil and gas fields are developed. The field data were acquired in April 2016 in Inner Mongolia, China, using the ADU instrument. A station located approximately 15 km away from the local station was used as a remote reference station, and a total of 16 hours of data were recorded synchronously. Our synthetic experiments focused on long-period MT data with a low sampling rate of 0.5 Hz. Therefore, field data recorded at higher sampling rates of 128 Hz and 4096 Hz were used to more comprehensively evaluate the performance of the three algorithms, suggesting that the proposed methods are applicable to both broadband and long-period MT data, and can effectively handle multi-scale electromagnetic signals.

Figure 8 compares the results from different methods. Overall, the four methods exhibit consistent trends over most of the investigated period range; however, noticeable differences can still be observed in the mid- to long-period band. Among the three robust estimation approaches, M-estimation exhibits the most pronounced outliers and local fluctuations, particularly in the apparent resistivity curves around the intermediate-period range. S-estimation also shows a relatively large number of scattered points and oscillatory variations, although its deviations are generally smaller than those of M-estimation. In contrast, the MM-estimation results are highly consistent with those obtained using the remote reference method, with the corresponding apparent resistivity and phase curves almost completely overlapping over most of the period band. In particular, the MM-estimation results exhibit smoother transitions and fewer abnormal deviations in the mid- to low-frequency range, indicating a more stable recovery of the MT transfer functions under complex field-noise conditions.

The observed differences among the methods are closely related to their respective robust weighting strategies and statistical properties. These robust estimations are primarily based on iterative residual weighting, where the influence of outliers is gradually reduced through a weighting function [64]. However, due to the relatively mild weight adjustment strategy of M-estimation, its ability to suppress outliers remains limited. S-estimation also exhibits a relatively large number of deviated points, which is likely because, although it achieves a high breakdown point by minimizing a robust scale, its statistical efficiency remains relatively low [58]. When a substantial portion of the frequency band is contaminated by noise, the scale estimate can become distorted. In real field measurements, where noise is typically more complex than in synthetic datasets, this limitation becomes more pronounced, leading to the widespread oscillations and scattered points observed in the curves. By contrast, MM-estimation combines the advantages of both S-estimation and M-estimation. The initial S-estimation provides a high-breakdown, noise-resistant scale that prevents early contamination, while the subsequent M-estimation step refines the solution with higher efficiency and improved smoothness. It is noteworthy that the MM-estimation results are nearly identical to those obtained using the remote reference method, indicating that MM-estimation can effectively recover subsurface electrical responses consistent with remote reference processing without requiring external reference station information. The remote reference method essentially reduces local uncorrelated noise by introducing remote reference signals, whereas MM-estimation actively suppresses anomalous spectral windows through robust statistical weighting functions. These results demonstrate that MM-estimation provides the highest robustness in processing real MT data and represents an effective choice for field records affected by complex background noise.

4. Discussions

4.1. Theoretical Analysis

Since natural EM signals observed in practical MT surveys are often non-stationary, time-varying, truncated, and contaminated by noise, the assumption of an ideal normal noise distribution is generally not satisfied for the linear regression model of MT transfer functions (Equations (1) and (2)). The primary distinction among the estimators lies in the definition of their objective functions: M-estimation minimizes a residual-based objective function, whereas S-estimation minimizes a scale-based objective function. These differences lead to distinct stability behaviors under varying SNR conditions. Under the condition of SNR = 30 dB, the results obtained using S-estimation exhibit relatively stronger dispersion (Figure 3 and Figure 4). When the SNR decreases to 5 dB, all three methods exhibit larger fluctuations and deviations in the mid- and long-period ranges (Figure 6 and Figure 7). This indicates that the robust estimation methods respond differently to changes in noise intensity, reflecting a clear trade-off between robustness and statistical efficiency. S-estimation employs a high-breakdown-point scale estimation strategy that strongly down-weights spectral windows with large residuals, thereby improving resistance to severe noise contamination. Since the noise level does not completely destroy the overall stability of the data, this aggressive weighting mechanism may simultaneously suppress some valid spectral windows, thereby reducing the utilization efficiency of useful data under relatively weak noise conditions. In contrast, the weight adjustment in M-estimation is relatively moderate; however, its ability to suppress outliers remains limited, such that some severely contaminated spectral windows may still retain non-negligible weights.

The MM-estimation results are generally consistent with those obtained using the remote reference method for the field data (Figure 8), indicating that the single-station processing approach based on robust statistical weighting can achieve noise suppression performance comparable to that of remote reference processing under the investigated conditions, without requiring additional reference stations. It is noteworthy that, since MM-estimation relies on the residuals of the dependent variables (electric field data), methods such as the remote reference technique or bounded-influence estimators [18,27] are still required to suppress outliers originating from the predictor variables (magnetic field data). We further suggest that a combination of MM-estimation and the remote reference method could provide advantages in handling both electric-field and magnetic-field noise at local stations. The superior performance of MM-estimation can be attributed to its two-stage structure, which combines the high breakdown point of S-estimation with the high statistical efficiency of M-estimation [59]. The S-step effectively suppresses high-amplitude transients, such as peak noise, by down-weighting outliers through the scale parameter, thereby providing a more reliable initial solution. The subsequent M-step refines the estimate using a smooth redescending weight function, improving convergence and reducing sensitivity to high-leverage points.

In addition, Tukey’s bisquare weight function [62] was employed in this study, although the use of alternative weight functions also merits consideration. Weight functions with a ψ-shape can effectively balance robustness and efficiency. During the iterative procedure, the bisquare scheme progressively suppresses the influence of outliers, enabling the estimator to focus on the dominant structure of the clean data [65]. Notably, the optimal selection of constants for the weight functions is largely based on classical experience and remains an open problem, requiring adaptive strategies for field datasets with varying noise distributions.

4.2. Limitations and Suggestions

Two major limitations of this study should be acknowledged: (1) the synthetic noise scenarios may not fully capture the complexity of industrial EM noise, and natural noise sources—such as long-period geomagnetic variations, instrument drift, and site-specific ground coupling effects—can also produce outliers; (2) the synthetic experiments were conducted by superimposing noise onto randomly selected segments of the time-series. This configuration simplifies the noise distributions, ensuring that sufficiently clean windows remain available for robust estimation. As a result, the evaluated performance may represent an optimistic scenario compared with real-field conditions. In practical MT surveys, cultural EM interference, instrument grounding issues, and environmental transients are generally persistent, widespread, or intermittently distributed throughout the entire recording. Under such circumstances, the breakdown point of the estimator, the availability of uncontaminated windows, and the convergence behavior of the iterative weighting scheme may all be adversely affected.

In this study, the primary purpose of the synthetic experiments is to validate the proposed method. It would be valuable to incorporate a more comprehensive and quantitative assessment of the estimation approaches by introducing noise disturbances with diverse frequency-domain signatures and a broader range of noise types. Evaluating robust estimators under these more realistic and complex noise conditions will help clarify their performance limits and inform the development of noise-specific weighting strategies. In addition, it is recognized that further research should focus on integrating adaptive weight selection procedures, multivariate robust estimators, and hybrid approaches that combine frequency-domain robust regression with time-domain machine-learning-based noise suppression.

5. Conclusions

We propose an iterative MM-estimation algorithm for the frequency-domain linear regression procedure in MT data processing. Together with the traditional M-estimation and S-estimation methods, we systematically evaluated the performance of all three robust regression estimators using both noise-contaminated synthetic data and field datasets. The results demonstrate that all estimators can produce reliable impedance tensors when the data are contaminated by a small amount of noise. For cases with a high proportion of noisy segments, MM-estimation demonstrates superior performance in suppressing outliers. We suggest that MM-estimation is less sensitive to industrial EM noise, which is increasingly becoming a key factor limiting the applicability of the MT method. We propose that MM-estimation is a strong candidate for routine MT processing workflows, particularly when data quality is highly variable.