Improved Variational Mode Decomposition for Magnetotelluric Data Denoising Combined with Transformer Neural Network

1. Introduction

Magnetotelluric (MT) is a natural source electromagnetic (EM) geophysics method, which is widely used in subsurface resistivity prospecting. Generally, timeseries of orthogonal EM fields are observed; and frequency dependent transfer functions, including impedance tensors and induction arrows, are calculated using linear regression algorithms [1]. As a prerequisite and basic guarantee, the data quality of transfer functions plays an import role for a final reliable resistivity model. And approaches in both time and frequency domain were developed in the last decades [2,3,4,5,6,7,8,9,10,11].

Technically, a general MT data processing includes following steps: (1) timeseries are subdivided into overlapping data segments where cascade decimations are also involved. (2) Each data segment is pre-whited and tapered with a Hanning or Hamming window. (3) Fourier transforms are employed to yield Fourier coefficients, which are then calibrated. (4) Spectral density matrix is calculated and pre-selecting corresponding to coherence and other criteria is employed. (5) Regression procedures generally regarded as robust estimation techniques are further conducted to get optimal transfer functions. Corresponding to the basic steps, various statistical methods and timeseries denoising algorithms were proposed. In practice, robust estimation techniques, basing on either dependent variable (electric field data) [2,4,7,12] or predictor variable (magnetic field data) [13], and remote reference methods [5,6] were employed in combination [2,14,15]. While synchronously observed array data are available, multivariate linear regression methods, including multivariate robust remote reference estimator [16], frequency-domain independent component analysis [17,18] and principal component analysis [3,19,20], can be practical to improve processing quality.

When EM noise persists throughout the entire observation period, large proportion of dataset would be severely contaminated, leading requirements of breakdown point cannot be satisfied for the robust estimation algorithms. Approaches of timeseries denoising can contribute to improve data quality to some extent. For instance, wavelet transform and/or Hilbert-Huang transform were employed for data segments pre-selecting [21,22]. Methods using wiener filtering with reference data [23] and compressive sensing [24] were used to replace or construct noised data segment. These approaches may be invalid while noisy data segments dominated.

Furthermore, mode decomposition techniques have become effective tools for MT data denoising because of its solid mathematical foundation and excellent capability in processing non-stationary signals. Chen et al. [25] utilized empirical mode decomposition (EMD) and Li et al. [26] further employed variational mode decomposition (VMD) to enhance robustness. Unlike EMD, VMD formulates signal decomposition as a constrained variational optimization problem and simultaneously estimates the center frequencies and bandwidths of all modes, resulting in improved decomposition stability and reduced mode mixing [27]. Since the denoising performance of VMD strongly depends on the selection of the mode number and penalty factor, adaptive optimization methods such as particle swarm optimization (PSO), sparrow search algorithm (SSA), and beluga whale optimization (BWO) have been introduced to improve parameter selection [28,29,30]. Wang et al. [31] combined VMD with wavelet thresholding (WT) and mathematical morphology filtering (MMF) to further suppress residual high-frequency noise while recovering useful low-frequency components discarded during decomposition. Nevertheless, VMD-based methods still rely heavily on parameter optimization, and the separation of strong interference from weak low-frequency MT signals remains challenging under complex noise conditions.

In addition, artificial neural networks show advantages in timeseries analysis and processing, including noise identification and removal, signal reconstruction and prediction, suggesting significant potential in geoscience [32]. Generally, a supervised learning requires manual identification and labeling of features (e.g., data quality, signal correlation, etc.) [11]. Typical neural networks, such as convolutional neural network (CNN), long-short term memory (LSTM), U-net, gated recurrent unit (GRU), etc., were utilized in an either independent or combination way to handle noise contaminated MT data and suppress EM noise with typical morphological features [9,10,33]. For instance, Li et al. [10] proposed a two-steps method using CNNs for signal-to-noise separation and LSTM for model training with denoised data from CNNs. Using the LSTM network, Tian et al. [34] utilized synchronously observed data from a reference station deployed in an underground laboratory as a model training set in a data-driven manner to reconstruct noisy segments in local stations. Multiple channels can be used in a combination way for neural networks training and noise depressing [35]. Furthermore, sparse coding and/or dictionary learning algorithms are employed in EM noise suppression studies. Feng et al. [36] developed a K-SVD dictionary learning algorithm to improve marine MT signal processing. A combination of deep convolutional learning classification and dictionary learning denoising shows advantages to achieve fully adaptive sparse decomposition of MT data [37]. Similar methods comprehensively using mathematical morphology filtering, convolutional networks, sparse coding, and dictionary learning were developed in previous studies [38,39,40,41]. Generally, multiple experiments are required to appropriately define specified coefficients, which are essential to achieve optimal results for these methods, remaining challenging in practice. More recently, Transformer-based [42] models have shown strong potential in geophysical and vibration signal processing due to their ability to capture long-range dependencies and global temporal correlations. As a typical limitation of the above methods, misfits may occur when continuous noise or long persistent noise data exist in timeseries, or noisy and clean datasets exhibit similar morphological features, leading ambiguous in longer period data. Besides, training datasets encompassing principal signal and noise morphological features control the efficacy and reliability of the noise depressing procedure.

In this study, we propose a joint MT denoising method that combines Transformer-based noise identification with a PSO optimized VMD algorithm. The Transformer model is first used to identify and locate strong interference noise in the original MT time series, and the abnormal segments are removed and repaired by linear interpolation. Meanwhile, PSO-VMD is applied to the original untreated signal, where PSO adaptively optimizes the mode number $K$ and penalty factor $\alpha$ to obtain the optimal decomposition result. Based on the optimized modal components, the noise-dominated component is removed, and the remaining effective high-frequency components are further filtered. Finally, the denoised MT signal is reconstructed by adding the filtered high-frequency components to the linearly interpolated signal. By combining the strong temporal noise recognition ability of Transformer with the adaptive decomposition capability of PSO-VMD, the proposed method can collaboratively suppress strong interference noise and random noise while preserving effective low-frequency information, thereby improving the stability of MT response functions.

2. Theory and Method

2.1. Variational Mode Decomposition

The VMD is an adaptive decomposition method for non-stationary signals [27]. It decomposes an input signal $f\left(t\right)$ (i.e., electric and magnetic components for MT timeseries) into $K$ band-limited intrinsic mode functions $u_k\left(t\right)$ with different center frequencies ${\omega }_k$, which can be written as:

$$f\left(t\right)={\displaystyle \sum _{k=1}^K}{u}_k\left(t\right),$$

(1)

The corresponding constrained variational problem is:

$$\begin{array}{c}min\\ \left\{u_k\},\{{\omega }_k\right\}\end{array}{\displaystyle \sum _{k=1}^K}{‖{\partial }_t\left[\left(\delta \left(t\right)+{\displaystyle \frac{j}{\pi t}}\right)u_k\left(t\right)\right]e^{-j{\omega }_{k}t}‖}_2^2,s.t.{\displaystyle \sum _{k=1}^K}{u}_k\left(t\right)=f\left(t\right),$$

(2)

By introducing a quadratic penalty factor and a Lagrange multiplier, the constrained problem can be solved iteratively using an algorithm named as alternate direction method of multipliers (ADMM) [27,43,44]. Finally, the original signal is adaptively decomposed into several modal components with different frequency bands, which provides a basis for subsequent noise identification and signal reconstruction.

The decomposition performance of VMD is strongly affected by the mode number $K$ and the penalty factor $\alpha$. If these parameters are selected improperly, under-decomposition, over-decomposition, or poor mode separation may occur. Therefore, PSO is introduced to adaptively search for the optimal VMD parameter combination.

2.2. PSO-Based Optimization of VMD Parameters

2.2.1. Basic Principle of PSO

The PSO is a population-based intelligent optimization algorithm inspired by the foraging behavior of bird flocks [45]. In PSO, each candidate solution is regarded as a particle in the search space. Each particle has two attributes: position and velocity. During iteration, particles update their search directions according to their own historical best position and the global best position of the swarm[45].

Assume that the position and velocity of the $i$-th particle at the $t$-th iteration are $x_i^t$ and $v_i^t$, respectively. The individual best position is denoted as $p_i$, and the global best position is denoted as $g$. The velocity and position are updated as follows:

$$v_i^{t+1}=v_i^t+c_1{r}_1\left(p_i-x_i^t\right)+c_2{r}_2\left(g-x_i^t\right),$$

(3)

$$x_i^{t+1}=x_i^t+\beta v_i^{t+1},$$

(4)

where $c_1$ and $c_2$ are learning factors, $r_1$ and $r_2$ are random numbers uniformly distributed in [0,1], and $\beta$ is a step-size factor. Through iterative updating, the particle swarm gradually approaches the optimal solution.

2.2.2. PSO-VMD Parameter Optimization

In the proposed PSO-VMD method, each particle represents a candidate parameter combination:

$$x_i=\left[K_i,{\alpha }_i\right],$$

(5)

where $K_i$ is the candidate mode number and ${\alpha }_i$ is the candidate penalty factor. For each particle, VMD is performed using the corresponding parameter combination, and the decomposition quality is evaluated by a fitness function.

First, the reconstructed signal is obtained by summing all decomposed intrinsic mode functions (IMFs), and the reconstruction error is calculated as:

$$E_{rec}= {\displaystyle \frac{1}{L}}{\sum }_{t=1}^L{\left(f\left(t\right) - {\displaystyle \sum _{k=1}^K}{u}_k\left(t\right)\right)}^2,$$

(6)

where $L$ is the signal length, $f\left(t\right)$ is the input signal, and $u_k\left(t\right)$ is the $k$-th IMF obtained by VMD.

To evaluate whether the noise component is concentrated in a specific IMF, the noise concentration index is defined as:

$$C= {\rho }_{k_n}{R}_n\left(1-{\bar{\rho }}_{\mathrm{o}}\right),$$

(7)

where ${\rho }_k=\mid \mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}\left(u_k\left(t\right),n\left(t\right)\right)\mid$, $t\in {\mathsf{\Omega }}_n$, is the absolute correlation coefficient between the $k$-th IMF and the detected noise pattern $n\left(t\right)$ within the detected noise region ${\mathsf{\Omega }}_n$. The dominant noise mode is determined by $k_n=\mathrm{a}\mathrm{r}\mathrm{g}{\mathrm{m}\mathrm{a}\mathrm{x}}_k{\rho }_k$, and ${\rho }_{k_n}$ is its corresponding maximum correlation coefficient. The term:

$$R_n={\displaystyle \frac{\sum t\in {\mathsf{\Omega }}_n u{k}_{n}2\left(t\right)}{{\displaystyle \sum _{k=1}^K}\sum t\in {\mathsf{\Omega }}_n u_k^2\left(t\right)}},$$

(8)

denotes the energy proportion of the dominant noise mode within the detected noise region, and ${\bar{\rho }}_{\mathrm{o}}$ is defined as:

$${\bar{\rho }}_{\mathrm{o}}={\displaystyle \frac{1}{K-1}}\sum k\ne 1 {\rho }_k,$$

(9)

The modal energy entropy is calculated as:

$$H=-{\displaystyle \sum _{k=1}^K}{p}_k\ln p_k,$$

(10)

$$p_k={\displaystyle \frac{E_k}{{\displaystyle \sum _{j=1}^K}{E}_j}},$$

(11)

$$E_k={\displaystyle \sum _{t=1}^L}{u}_k^2\left(t\right),$$

(12)

where $p_k$ denotes the normalized energy proportion of the $k$-th IMF. A smaller energy entropy indicates a more concentrated modal energy distribution.

Accordingly, the fitness function is formulated as:

$$F=w_1{E}_{rec}+w_2\left(1-C\right)+w_{3}H,$$

(13)

where the weighting coefficients are set to $w_1=0.3,w_2=0.5,w_3=0.2.$ Finally, the optimal VMD parameters are obtained by minimizing the fitness function:

$$\left(K^{*},{\alpha }^{*}\right)=arg\underset{K,\alpha }{\min }F,$$

(14)

In this study the particle population size is set to 18, and the maximum number of iterations is 10. The learning factors are set to $c_1=c_2=1.5$, and the position update step-size factor is set to 0.5. The search ranges of the mode number and penalty factor are $K\in \left[\mathrm{4,9}\right]$ and $\alpha \in [2000,8000]$, respectively. The corresponding velocity ranges are $v_K\in [-\mathrm{2,2}]$ and $v_{\alpha }\in [-2000,2000]$. For VMD, the parameters are set to $\tau =0$, $DC=0$, $init=1$, and $tol={10}^{-7}$. These settings are consistent with the uploaded implementation. The detailed procedure is summarized in Algorithm 1.

Algorithm 1: PSO Optimization for VMD Parameters

$$\begin{gathered} \mathrm{Input}\text{:} \mathrm{MT} \mathrm{time} \mathrm{series} f\left(\mathrm{t}\right), \mathrm{search} \mathrm{ranges} \mathrm{of} \mathrm{K} \mathrm{and} \mathsf{\alpha }, \mathrm{population} \mathrm{size} \mathrm{N}, \mathrm{maximum} \mathrm{iteration} \mathrm{number} \mathrm{G}, \mathrm{learning} \mathrm{factors} c_1 \mathrm{and} c_2 \\ \mathrm{Output}\text{:} \mathrm{Optimal} \mathrm{VMD} \mathrm{parameters} K^{*},{ \alpha }^{*}. \\ 1. \mathrm{Initialize} \mathrm{particle} \mathrm{positions} x_i^0=[K_i^0,{\alpha }_i^0] \mathrm{and} \mathrm{velocities} v_i^0, i=\mathrm{1,2},\dots ,N; \\ 2. \mathrm{For} \mathrm{each} \mathrm{particle}, \mathrm{perform} \mathrm{VMD} \mathrm{decomposition} \mathrm{using} x_i^0=[K_i^0,{\alpha }_i^0]; \\ \mathrm{br}-\mathrm{to}-\mathrm{break} F_i=w_1{E}_{rec}+w_2\left(1-C\right)+w_{3}H;4. \mathrm{Set} \mathrm{the} \mathrm{personal} \mathrm{best} \mathrm{position} p_i=x_i^0 \\ 5. \mathrm{while} \mathrm{g}<G \\ 7. \mathrm{for} i=1:N \\ \mathrm{br}-\mathrm{to}-\mathrm{break} v_i^{\mathrm{g}}=v_i^{\mathrm{g}-1}+c_1{r}_1\left(p_i-x_i^{\mathrm{g}-1}\right)+c_2{r}_2\left(xbest -x_i^{\mathrm{g}-1}\right);\mathrm{br}-\mathrm{to}-\mathrm{break} x_i^{\mathrm{g}}=x_i^{\mathrm{g}-1}+\beta v_i^{\mathrm{g}};11. \mathrm{Apply} \mathrm{boundary} \mathrm{constraints} \mathrm{to} K_i^{\mathrm{g}} \mathrm{and} {\alpha }_i^{\mathrm{g}}; \\ 12. \mathrm{Perform} \mathrm{VMD} \mathrm{decomposition} \mathrm{using} \mathrm{the} \mathrm{updated} x_i^{\mathrm{g}}=[K_i^g,{\alpha }_i^g]; \\ 13. \mathrm{Calculate} \mathrm{the} \mathrm{current} \mathrm{fitness} \mathrm{value} F_i^g; \\ 15. \mathrm{for} i=1:N \\ 16. \mathrm{if} F_i^g<F\left(p_i\right) \\ 17. p_i=x_i^g \\ 19. \mathrm{if} F_i^g<F\left(xbest \right) \\ 20. xbest =x_i^g \\ 24. \mathrm{Return} xbest =[K^{*},{\alpha }^{*}]. \end{gathered}$$

, step-size factor β.

3. Calculate the fitness value:

, and set the global best position as the particle with the minimum fitness value;

 and the stopping criterion is not satisfied do

6.    g = g + 1

 do

8.        Update particle velocity:

9.        Apply velocity constraints;

10.       Update particle position:

14.   end for

 do

then

18.      end if

then

21.      end if

22.    end for

23. end while

Finally, the optimal particle $g$ is selected as the optimal VMD parameter combination $\left[K^{*},{\alpha }^{*}\right]$. Compared with empirical parameter selection, the PSO-based strategy can automatically search for suitable VMD parameters, thereby improving the adaptivity and stability of VMD decomposition for non-stationary MT signals.

2.3. Intelligent Filtering

To further suppress transient spike interference while preserving the underlying MT signal, an intelligent filtering strategy is applied to the reconstructed signal obtained from VMD. Instead of directly removing abnormal samples, the proposed method first detects transient spike regions using multiple complementary criteria and then reconstructs these regions according to the local statistical characteristics of neighboring normal samples. Specifically, the reconstructed signal is first obtained by:

$$\mathrm{y}\left(t\right)={\displaystyle \sum _{k=2}^K}{u}_k\left(t\right),$$

(15)

where $u_k\left(t\right)$ denotes the $k$-th IMF. The first IMF, which mainly contains low-frequency noise interference, is excluded from the reconstruction. To identify abnormal samples, amplitude, local statistical deviation, and signal gradient are jointly considered. The local trend is first estimated using a median filter:

$$m\left(t\right)=MedFilt\left\{\mathrm{y}\left(t\right)\}\right.,$$

(16)

and three detection masks are constructed according to amplitude, local deviation, and gradient, respectively. The final spike mask is obtained by:

$$M\left(t\right)=M_a\left(t\right)\bigvee M_s\left(t\right)\bigvee M_g\left(t\right),$$

(17)

where $M_a\left(t\right)$, $M_s\left(t\right)$, and $M_g\left(t\right)$ denote the amplitude-, statistical-, and gradient-based detection results, respectively. The detected regions are subsequently expanded using morphological dilation to ensure complete coverage of the transient disturbance.

For each detected abnormal interval, the local trend is estimated using a median filter, while the local fluctuation is synthesized using a first-order autoregressive (AR(1)) process estimated from the surrounding normal signal:

$$z\left(t\right)= \varphi z\left(t-1\right)+ \epsilon \left(t\right),$$

(18)

where $\varphi$ is the autoregressive coefficient and $\epsilon \left(t\right)$ is a Gaussian innovation term. The reconstructed signal within the abnormal region is then generated as:

$$\widehat{y}\left(t\right)= m\left(t\right)+z\left(t\right),$$

(19)

where $m\left(t\right)$ denotes the local trend estimated from neighboring normal samples. Finally, a weighted transition together with a Savitzky–Golay filter [46,47] is applied near the boundaries of each reconstructed interval to ensure smooth continuity. The final filtered signal is expressed as:

$${\mathrm{y}}_f\left(t\right)=\left\{\begin{array}{c}\widehat{y}\left(t\right), t\in {\mathsf{\Omega }}_n\\ y\left(t\right), t\notin {\mathsf{\Omega }}_n\end{array}\right.,$$

(20)

where ${\mathsf{\Omega }}_n$ denotes the detected abnormal region.

In this study, the median filter window, local statistical window, and Savitzky–Golay smoothing window were set to 5001, 1001, and 101 samples, respectively. The amplitude threshold was set to 5000, while the gradient threshold was determined as the 99.5th percentile of the gradient distribution. The detected spike regions were expanded by 50 morphological dilation iterations before reconstruction, and 1500 neighboring samples on each side of each detected interval were used for AR(1) parameter estimation.

2.4. Transformer-Based Noise Identification

An Anomaly-Transformer model [48] is used to identify strong interference noise in MT time series. Given an input sequence $X\in {\mathbb{R}}^{N\times d}$, token embedding and positional embedding are first applied to obtain high-dimensional temporal features. The embedded sequence is then fed into stacked Transformer encoder layers, which consist of anomaly-attention modules, feed-forward networks, residual connections, and layer normalization. The structure of the Anomaly-Transformer used for time series noise identification is shown in Figure 1.

To enable the Transformer model to learn the intrinsic temporal characteristics of MT signals, a dedicated noise identification dataset was constructed before model training. The original MT observations were first normalized using Z-score normalization to eliminate amplitude differences among different channels. Subsequently, a sliding-window strategy was adopted to segment the continuous time series into fixed-length samples, where each sample contains local temporal information and long-range correlation characteristics.

The constructed dataset mainly consists of MT signals with weak interference, which are regarded as normal patterns for unsupervised anomaly detection. The model learns the temporal association characteristics of normal MT signals during training without requiring manual labels. During testing, typical strong interference signals, including abrupt disturbances and pulse-like noise, are introduced and labeled to evaluate the capability of the trained model in identifying abnormal segments.

After training, the Transformer model generates anomaly scores based on the discrepancy between prior association and series association. The segments with high anomaly scores are considered strong interference noise and are further processed in the subsequent denoising procedure.

The core component of the model is the anomaly-attention mechanism, which simultaneously models the global temporal dependency and the local neighborhood dependency. Specifically, the hidden feature from the previous layer is linearly projected to generate the query, key, value, and learnable scale parameter: $Q=X{W}_Q$, $K=X{W}_K$, $\sigma =X{W}_{\sigma }$. Based on the learnable scale parameter, a Gaussian kernel is employed to construct the prior association, which explicitly characterizes the local correlation around each timestamp:

$$P_{i,j}^l={\displaystyle \frac{1}{\sqrt{2\mathsf{\pi }}{\mathsf{\sigma }}_i}}\exp \left(-{\displaystyle \frac{{\left|j-i\right|}^2}{2{\mathsf{\sigma }}_i^2}}\right),$$

(21)

Meanwhile, the series association is obtained using the standard self-attention mechanism to capture long-range temporal dependency,

$$S=Softmax\left({\displaystyle \frac{Q{K}^T}{\sqrt{d_{model}}}}\right),$$

(22)

For normal MT signals, long-range temporal dependency dominates the sequence, resulting in a noticeable discrepancy between the prior association and the series association. In contrast, strong interference noise is generally characterized by impulsive and locally continuous patterns, causing the learned series association to concentrate on local neighborhoods and become more similar to the prior association. Therefore, the discrepancy between the two associations can effectively distinguish abnormal interference from normal MT signals. Following the original Anomaly Transformer, the association discrepancy is quantified using the symmetric Kullback-Leibler (KL) divergence[49],

$$\mathrm{AssDis}\left(P,S\right)={\displaystyle \frac{1}{L}}{\sum }_{l=1}^L\left(\mathrm{KL}\left(P^l\parallel S^l\right)+\mathrm{KL}\left(S^l\parallel P^l\right)\right),$$

(23)

To improve the robustness of noise identification, the association discrepancy is further combined with the sequence reconstruction error to calculate the anomaly score:

$$Score\left(X\right)=Softmax\left(-AssDis\right)\odot |X-\widehat{X}{|}_2^2,$$

(24)

where $\widehat{X}$ denotes the reconstructed sequence. Time points with anomaly scores exceeding a predefined threshold are identified as strong interference noise. The detected noise intervals are subsequently provided as prior information for the following denoising stage, enabling the proposed framework to suppress interference while preserving the intrinsic characteristics of MT signals.

2.5. Workflow of the Proposed Method

The proposed framework consists of three main stages: Transformer-based strong interference noise recognition and preprocessing, PSO-optimized adaptive VMD decomposition and filtering, and final signal reconstruction. The workflow is shown in Figure 2.

First, the original MT time series is input into the trained Transformer-based anomaly detection model. The model identifies strong interference segments according to anomaly scores, and the detected abnormal segments are removed and repaired by linear interpolation to suppress impulsive noise while maintaining signal continuity.

Second, the original signal is decomposed by VMD, where PSO is introduced to adaptively optimize the mode number $K$ and penalty factor $\alpha$. The optimized VMD decomposition results are analyzed according to spectral characteristics, and noise-dominated components are removed while effective components are retained for further filtering.

Finally, the filtered VMD reconstruction components are combined with the repaired Transformer outputs to reconstruct the final denoised MT signal. The proposed framework integrates Transformer-based strong interference recognition and PSO-VMD adaptive decomposition, enabling simultaneous suppression of impulsive interference and random high-frequency noise.

3. Implementation

3.1. Synthetic Experiment

We employed example MT data from the open-source EMTF package [50] to evaluate the proposed method. Five orthogonal electric and magnetic channels corresponding to a uniform 100 Ωm half-space model, with a sampling rate of 1 Hz and a total duration of 10,000 s were used. The horizontal electric field component Ex was selected as the test signal. To simulate strong interference conditions commonly encountered in practical MT observations, triangular and tilted-top square noise clusters were artificially introduced into the clean Ex time series. The four noise segments were located at approximately 550–1664 s, 3150–4299 s, 5750–6816 s, and 8350–9450 s, with amplitudes of approximately +500, −450, +830, and −530 mV/km, respectively. These artificial disturbances were designed to represent strong transient interference with different polarities and amplitudes in field MT data.

Figure 3 presents the Transformer-based identification results of strong interference noises and the subsequent signal reconstruction process. As shown in Figure 3, the Transformer model identifies interference segments with different waveform characteristics, polarity, and amplitude, demonstrating its capability to capture abnormal temporal patterns and distinguish strong interference from normal MT signal fluctuations. Based on the identified noise locations, the abnormal segments were removed and the missing intervals were subsequently reconstructed using linear interpolation, as shown in Figure 3(b). After noise removal, the reconstructed Ex signal eliminates the large amplitude deviations introduced by strong interference while maintaining the continuity and general variation trend of the original MT response. The results indicate that the proposed Transformer-based method can not only effectively locate strong interference noises in MT data but also provide reliable preprocessing for subsequent signal analysis and inversion.

Before VMD decomposition, the key parameters of VMD, including the number of modes $K$ and the penalty factor $\alpha$, were optimized using the particle swarm optimization (PSO) algorithm. The particle population was set to 18 with a maximum of 10 iterations. The search ranges of $K$ and $\alpha$ were defined as 4–9 and 2000–8000, respectively. The fitness function was constructed based on the decomposition performance of VMD, where a lower fitness value indicates better modal separation.

Figure 4 illustrates the convergence process of the PSO optimization. During the initialization stage, particles with different combinations of $K$ and $\alpha$ were randomly generated within the predefined search space. Through iterative particle updating, the optimization process rapidly converged within the first few iterations, and the fitness value gradually reached a stable minimum. Finally, the optimal parameters were determined as $K=8$ and $\alpha =2000$. These optimized parameters were subsequently adopted for VMD decomposition of the noisy Ex signal.

After determining the optimal VMD parameters through PSO, the repaired MT Ex signal was decomposed using the PSO-VMD method. As shown in Figure 5, the signal was adaptively decomposed into eight intrinsic mode functions (IMFs), where each component represents different frequency characteristics of the original signal. The first mode (IMF1, Figure 5(a)) mainly contains the low-frequency trend component. It exhibits significant amplitude variations corresponding to the locations of the artificially introduced strong interference, indicating that the abnormal disturbances are mainly concentrated in this low-frequency component. Therefore, IMF1 was regarded as the interference-dominated trend component and excluded during signal reconstruction. The intermediate modes (IMF2–IMF4, Figure 5b-5d) mainly contain medium-frequency oscillatory information and residual transient characteristics. These components preserve the dynamic variations of the MT signal while gradually reducing the influence of strong interference. The higher-order modes (IMF5–IMF8, Figure 5e-5h) mainly represent high-frequency fluctuations and fine-scale signal characteristics. With increasing mode order, the oscillation frequency gradually increases and the amplitude decreases, demonstrating that PSO-VMD can effectively separate the non-stationary MT signal into components with different temporal-frequency characteristics.

After identifying the interference-dominated IMF1, the remaining modal components (IMF2–IMF8) were reconstructed to preserve the useful signal information distributed over different frequency bands (Figure 6). However, the reconstructed signal still contains residual high-frequency fluctuations caused by weak interference and random noise. Therefore, an additional filtering operation was applied to the reconstructed high-frequency signal to suppress these residual noise components before the final signal reconstruction. Figure 6(b) compares the clean signal, the reconstructed high-frequency signal before filtering, and the filtered high-frequency signal, demonstrating that the filtering operation effectively suppresses residual random noise while maintaining the useful high-frequency characteristics required for subsequent local signal reconstruction.

After obtaining the filtered high-frequency components from PSO-VMD, they were selectively superimposed onto the signal repaired by Transformer-based interference identification and linear interpolation, as illustrated in Figure 7. Specifically, the filtered high-frequency components were added only within the strong interference intervals identified by the Transformer model, where the original signal had been removed and reconstructed by linear interpolation. Outside these repaired intervals, the original MT signal was retained without any modification. Consequently, the proposed strategy reconstructs only the contaminated segments while preserving the original waveform characteristics of the uncontaminated regions.

The final denoising results are presented in Figure 7(b) and 7(c). Compared with the noisy MT signal, the proposed method effectively removes strong interference while further suppressing residual random fluctuations. As shown in the enlarged comparison in Figure 7(c), the reconstructed signal agrees well with the clean signal in both waveform and amplitude within the repaired interval, indicating that the proposed local reconstruction strategy effectively restores the contaminated segments without introducing unnecessary distortion into the surrounding signal.

These results demonstrate that the proposed hybrid denoising framework effectively combines Transformer-based temporal interference localization with PSO-VMD-based adaptive random noise suppression. The Transformer model accurately identifies and repairs strong interference intervals, whereas PSO-VMD extracts and filters the useful high-frequency components, which are subsequently used to reconstruct only the repaired intervals. By restricting the reconstruction process to the detected interference regions, the proposed method effectively suppresses noise while preserving the original characteristics of uncontaminated MT signals.

To further evaluate the effectiveness of the proposed denoising method in MT response estimation, the same noise contamination and denoising procedures were applied to all five MT components. The impedance tensors were then estimated using the EMTF package [50] with the clean, noisy, and denoised five-component datasets. Figure 8 shows the apparent resistivity and impedance phase curves obtained from these three datasets. The apparent resistivity and impedance phase curves exhibit obvious distortion for the noisy data, where the resistivity responses of both xy and yx modes deviate significantly from the theoretical value of 100 Ωm over the whole frequency range, and the impedance phases show apparent deviations from the expected 45° response, especially in the longer period band (Figure 8), indicating a near-field effect. In contrast, the results obtained from the denoised data show strong agreement with those from the clean data in both apparent resistivity and phase responses. These results demonstrate that the proposed Transformer and PSO-VMD combined denoising method can effectively suppress strong interference and random noise while preserving the useful MT response characteristics.

3.2. Field Data Experiment

To further verify the applicability of the proposed method in real complex observation environments, field MT time-series data observed in southern Tibet using LEMI-417 were used for denoising experiments. Five orthogonal components data, including three magnetic channels and two electric channels, were recorded with a 1 Hz sampling rate. The east-west orientated electric field (Ey) was contaminated by random noises probably caused by unstable non-polarized electrode.

Figure 9 presents the Transformer-based identification results of strong interference noise in the observed Ey time series. As shown in Figure 9(a), the field MT data contain various types of non-stationary interference, including long-duration step-like disturbances, impulsive spikes, and local abnormal fluctuations. Compared with the synthetic experiments, the interference characteristics are more complicated because multiple noise types coexist in the observed data. The Transformer model accurately identifies the abnormal segments with different durations and amplitudes, demonstrating its capability to capture complex temporal patterns in real MT observations. An enlarged view is shown in Figure 9(b), where the detected interference segments are highlighted. It can be observed that both long-duration step-like disturbances and short impulsive anomalies are effectively identified, providing reliable noise localization for the subsequent signal reconstruction process.

Based on the identified noise locations, the abnormal segments were removed and the missing intervals were reconstructed using linear interpolation, as shown in Figure 10. Figure 10(a) compares the original and reconstructed Ey time series, while Figure 10(b) provides an enlarged view of the reconstructed interval. After removing the detected interference, the reconstructed signal effectively eliminates the large-amplitude distortions caused by strong cultural noise while preserving the continuity and long-term variation trend of the original MT response. These preprocessing results provide a reliable input for the subsequent PSO-VMD decomposition and random noise suppression.

Before VMD decomposition, the number of modes K and the penalty factor α were adaptively optimized using the PSO algorithm to improve the decomposition performance for the field MT signal. Figure 11 illustrates the convergence process of the PSO optimization. During the initialization stage, particles with different combinations of $K$ and $\alpha$ were randomly generated within the predefined search space. Through iterative particle updating, the optimization process gradually converged, and the fitness value decreased to 15000 after the second iteration. Finally, the optimal parameters were determined as $K=9$ and $\alpha =2000$, achieving the minimum fitness value of 13000. These optimized parameters were subsequently adopted for the subsequent VMD decomposition and random noise suppression of the field Ey signal.

The PSO-optimized VMD method was further applied to the observed Ey signal for adaptive modal decomposition. As shown in Figure 12, the field signal was decomposed into nine intrinsic mode functions with distinct frequency characteristics. IMF1 mainly contains the low-frequency trend together with the dominant step-like interference observed in the field data, indicating that most of the large-scale non-stationary disturbances are concentrated in this mode. Therefore, IMF1 was regarded as the interference-dominated trend component and removed during the subsequent signal reconstruction. The remaining modes from IMF2 to IMF9 mainly reflect oscillatory components in different frequency bands. With the increase of mode order, the modal frequency gradually increases and the amplitude decreases, indicating that PSO-VMD can adaptively separate the complex non-stationary field signal into different frequency-scale components.

After removing the interference-dominated IMF1, the remaining modes (IMF2 – IMF9) were reconstructed to obtain the high-frequency components containing useful signal information together with residual random noise. As shown in Figure 13(a), the reconstructed high-frequency signal still exhibits noticeable impulsive fluctuations in several interference-affected intervals, indicating that part of the residual noise remains after modal reconstruction. A filtering operation was therefore applied to the reconstructed high-frequency components to suppress these fluctuations while retaining the useful high-frequency information. Figure 13(b) presents an enlarged comparison before and after filtering. This selective reconstruction strategy suppresses residual random noise within the repaired segments while avoiding unnecessary alteration of the uncontaminated portions of the field MT signal.

After obtaining the filtered high-frequency components from the reconstructed IMF2 – IMF9 modes (Figure 13), these components were selectively superimposed onto the linearly interpolated signal reconstructed after Transformer-based interference identification and removal (Figure 10) to generate the final denoised field MT signal (Figure 14). Specifically, the filtered high-frequency components were only added to the intervals where the strong interference had been removed and repaired by linear interpolation, whereas the original observed MT signal outside these repaired intervals was preserved without any modification. Consequently, the proposed reconstruction strategy restores only the interference-contaminated segments while maintaining the original waveform characteristics of the uncontaminated portions of the field MT signal.

As illustrated in Figure 14(a), the proposed method effectively suppresses large-amplitude step-like interference and local abnormal fluctuations in the field MT signal while maintaining the overall temporal variation trend. Figure 14(b) presents an enlarged comparison within a representative interference interval. It can be observed that the reconstructed signal effectively restores the waveform continuity after removing the strong interference and suppresses residual random fluctuations without introducing noticeable distortion into the surrounding uncontaminated signal. These results demonstrate that the proposed local reconstruction strategy successfully repairs the interference-contaminated intervals while preserving the original characteristics of the remaining field MT signal.

We further employed the RRMC algorithm [51] to estimate impedance tensors using both original observed and denoised data. Figure 15 shows the apparent resistivity and impedance phase curves of the original observed and denoised data. The results show that the resistivity curve of the yx mode from the observed data is severely distorted in the period of > 4000 s, and the phase curve approaches to be ${-180}^{\circ }$ in the longer period. After the denoising processing, the distorted sounding curve was corrected, demonstrating that the proposed method can effectively improve the quality and stability of field MT transfer functions. Additionally, the xy modes derived from both the original and denoised data show phases out-of-quadrant at periods greater than 9000s, which are probably due to electrical anisotropy, galvanic distortion, and/or local 3D conductive bodies [52], but not due to the influence of typical background EM noise.

4. Conclusions

To address the contamination of MT time series by strong cultural interference and random noise in complex observation environments, this study proposed a joint denoising method combining Transformer-based noise identification and PSO-VMD adaptive decomposition. The method was verified using both synthetic and field MT data. The main conclusions are as follows.

• The Transformer model can accurately identify strong interference noise segments in MT time series, including step-like anomalies, impulsive peaks, and complex non-stationary noise. By removing abnormal segments and repairing them through linear interpolation, the low-frequency trend of the original signal can be effectively restored, providing a stable basis for subsequent signal reconstruction;
• PSO can adaptively optimize the VMD parameters, thereby improving the modal decomposition ability and stability for complex non-stationary signals. The PSO-optimized VMD decomposition can better distinguish low-frequency noise-dominated trends from useful high-frequency information and reduce mode mixing caused by empirical parameter selection;
• After removing the noise-dominated mode, the remaining high-frequency components are further located and filtered to suppress random high-frequency noise and residual impulsive interference. The final reconstruction combines the effective high-frequency components with the interpolated low-frequency signal, achieving coordinated suppression of strong interference and random noise;
• Synthetic data experiments show that the proposed method can effectively recover the true variation characteristics of contaminated MT time series, improve the signal-to-noise ratio, and enhance the continuity and stability of apparent resistivity and impedance phase curves. Field data experiments further demonstrate that the method has good applicability and stability in complex cultural-noise environments and can improve the reliability of low-frequency effective information and MT response estimation.

Overall, the proposed Transformer and PSO-VMD joint denoising method combines the advantages of deep learning in abnormal noise identification and VMD in adaptive decomposition of non-stationary signals, providing an effective technical approach for MT data processing under complex noise conditions.

However, some limitations remain. The recognition performance of the Transformer model depends on the quality of training samples and the coverage of noise types. When field data contain more complex or unknown noise patterns, the generalization ability of the model still needs further improvement. In addition, the PSO-VMD optimization process requires a certain computational cost, which may affect processing efficiency for long-duration or large-scale MT datasets.