From the perspective of whether the data follow the linear relationship in Eq. 1, the field data ( E and H ) can be subdivided into three parts as follows: E = E M T , E H L N , E L L N , (2) H = H M T , H H L N , H L L N , (3) where the superscript HLN denotes the data dominated by noise with high linearity, the superscript LLN denotes the data dominated by noise with low linearity, and the superscript MT denotes the high-quality data with high linearity. Noise with low linearity can be identified from the similarity between the observed electric field and that predicted by the linear relationship. It is similar to the single-input/single-output linear model to evaluate the linearity, as shown in Figure 1. The difference between the data with high linearity and low linearity is that most of the predicted and observed values of the output are similar.

Assuming the data are highly linear related, the observed electric field ( E ) should be similar to the predicted electric field ( E p ), where E p = ZH and Z are obtained by the least-squares estimator. We can identify noisy data with low linearity by comparing the measured electric field ( E ) and the predicted electric field ( E p ). The complex number has two properties: the amplitude and phase. In this study, we use the linear coherence defined by the phase difference between the predicted and observed electric fields to confirm the phase similarity and use the amplitude ratio between the predicted and observed electric fields to confirm the amplitude similarity.

The linear coherence (Lcoh) between two spectra A i and B i is defined by the cosine of the phase difference (PD) as follows: cos θ i = R e e j φ A i − φ B i = R e A i B ¯ i A i ‖ B i , (4) where A i and B i denote the spectrum calculated from the i ^th segment, B ¯ i represents a conjugate of B i , and θ i denotes the angle of the phase difference (PD) between A i and B i . According to Euler’s formula, Lcoh equals the real part of e j ( φ A i − φ B i . Re denotes the real part of the complex number. The value of Lcoh lies in the range of (−1,1). When the PD is close to 0 ° , the Lcoh is high and close to 1. In this study, the predicted linear coherence (P L c o h ) between the measured electric field ( E ) and the predicted electric field ( E p ) is calculated as follows: P L c o h = R e Y p i Y ¯ i Y p i ‖ Y i , (5) where Y p i and Y i are the predicted and measured electric fields corresponding to the i ^th segment, and Y is associated with either E x or E y . The predicted electric field and the observed electric field should be similar when the linearity is high, and PLcoh should be close to 1.

The high predicted linear coherence can ensure the phase similarity. We also use the predicted amplitude ratio (PAR) to ensure the amplitude similarity, and it is defined as follows: P A R = { Y p i Y i , Y p i < Y i , o r Y i Y p i , Y i < Y p i . (6)

The PAR ranges from 0 to 1; the higher PAR is, the higher the similarity between the two spectra in terms of the amplitude.

There is a problem that the energy of the signal and noise changes with time, and the linearity may change. Suppose we perform regression with all the available data; the linearity may be low in the presence of a large amount of noise and provide misleading information. It is similar to the single-input/single-output linear model, as shown in Figure 2. There are six datasets, and only dataset 1 has high linearity. When we perform regression with all the available data, the linearity is low, and we cannot extract data with high linearity. To solve this problem, we subdivide the data into small groups and calculate the PAR and PLcoh values separately. We rename the predicted linear coherence and amplitude ratio as P L c o h s z and P A R s z , where the subscript SZ means we calculate the predicted electric field ( E p = ZH ) by the impedance ( Z ) obtained by the subdivided data. In this research, the field data are subdivided into small groups with 20 samples when evaluating the linearity. It is necessary to divide the data into groups when evaluating linearity (see Supplementary Figures S2, S3).

Fowler et al. (1967) proposed the polarization direction, and Weckmann et al. (2005) showed the effectiveness of MPD in detecting coherent noise. The MPD ( α H i ) at a specific frequency is defined as follows: α H i = tan − 1 2 R e H x i H ¯ y i H x i 2 − H y i 2 = tan − 1 2 H y i H x i · cos θ i 1 − H y i H x i 2 , (7) where i (=1,2, … , N) is the number index of the event, H x i and H y i are the spectra of the magnetic field calculated from the i ^th segment, and θ i denotes the PD between H x i and H y i . The polarization direction is related to the PD and amplitude ratio (AR) between the two orthogonal fields. Various sources generate natural magnetic signals that vary in incident directions and energy, and the PD and AR between the two orthogonal fields vary with time; thus, the magnetic field has no preferred polarization direction (Weckmann et al., 2005). In contrast, the local EM noise source usually has a constant location; the incident direction and energy have similar properties that change with time. Suppose there is a preferred polarization direction for the magnetic field; we can consider that the coherent noise contaminates the data. On the other hand, when incoherent noise contaminates the field data, the magnetic field has no preferred polarization direction. Therefore, the polarization direction for the magnetic field can only detect coherent noise.

To quantify the dispersion degree of MPD, the dispersion degree of the polarization directions ( D D p o l ) is proposed as follows: D D p o l = N i n N , (8) where N i n denotes the number of samples falling in the range of ( m i + 30 ° , m i -30 ° ). m i is the median for each α H i with its surrounding 2k samples (k is set to 20 in this study), and it is calculated as follows: m i = m e d i a n α H i − k , α H i − k − 1 , … , α H i , … , α H i + k − 1 , α H i + k . (9)

When the polarization direction has a preferred direction, m i approximately equal to the preferred direction. m i is calculated by the surrounding 2k samples, and there are two sides; we hope half of the data is beyond a specific range, which means the threshold is set as 0.5; therefore, the expected value of D D p o l should be smaller than 0.5, and 1/3 is chosen in this research, which means the range is 60 ° (180 ° × 1/3), and the specific range is set as ( m i + 30 ° , m i -30 ° ). If the polarization directions vary randomly from −90 ° to 90 ° , D D p o l should be close to 1/3, and D D p o l increases when there is a preferred direction. D D p o l can be used to automatically detect coherent noise with a strong polarization direction.

The first case study used the data observed at TVN48 from the USArray project. The time-series data can be downloaded from the Incorporated Research Institutions for Seismology (IRIS) website. The data sampling period is 1 s, and the used times-series data are observed from July 19 to 24 July 2015. First, we examine the variation in the parameters at different frequencies. Figure 4 shows the parameter variation in the period of 13.1 s. Figures 4A, B show the variation in P L c o h s z and P A R s z associated with E x -component, respectively. The events in which both P L c o h s z and P A R s z are larger than 0.8 are shown in red, and the other is shown in blue. Red denotes an event with high linearity, and blue denotes an event with low linearity. Figure 4C shows that the MPD is scattered for all the events. Figure 4D shows the variation in the hat matrix’s diagonal element. The hat matrix is an N by N matrix (N denotes the number of events) defined as follows (Chave and Thomson, 2004; 2003): H h a t = H H † H − 1 H † , (10) where H represents N by two matrices of the horizontal magnetic field ( H x , H y ) at a specific frequency. The superscript † denotes the complex conjugate transpose. The expected value of the hat matrix’s diagonal element is 2/N. The variation in the diagonal elements of the hat matrix has the same trend as the magnetic field amplitude (Chen et al., 2022; Li et al., 2022; Li et al., 2023; Li et al., 2018; Zhang et al., 2022; Zhang et al., 2021). Therefore, we can use the hat matrix to visualize the energy variation in the magnetic field. Figure 4D shows that the red events have high energy. This is caused by the geomagnetic storm (see Supplementary Figure S1). Since the natural EM signal is relatively low in the dead band (0.1–10 s), local noise can easily influence it during non-storm periods. When there is a geomagnetic storm, the natural EM signal strength increases, and high SNR events appear. In conclusion, the blue events are dominated by incoherent noise in the period of 13.1 s, and P L c o h s z and P A R s z can identify incoherent noise. Figure 5 shows the parameter variation in the period of 105.4 s. Most of the events have high linearity, and the MPD is scattered for all the events. This indicates that only a small part of the events are contaminated by incoherent noise. After analyzing the variation in the parameters in different periods, we found that most of the events are dominated by incoherent noise between 5 and 20 s.

Then, we compare the MT sounding curves calculated by the different methods, as shown in Figure 6. All of those responses are estimated by M-estimator (Egbert and Booker, 1986; Neukirch and García, 2014; Maronna et al., 2019). The result using the data preselection strategy with P L c o h s z and P A R s z coincides with the remote reference result, and they are regarded as the true model. It shows that a reliable result can be obtained even if we do not use the remote site data by the preselection method. On the other hand, the apparent resistivity of the robust results is downbiased between 6 and 20 s. According to the analysis of Figure 4, more than half of the events are contaminated by incoherent noise. The underestimation of the apparent resistivity was probably attributed to the auto-power spectra of the noise in the denominator of the response function, which is a well-known limitation of the single-site data processing (Sims et al., 1971; Simpson and Bahr, 2005).

This subsection compares the performance of the related parameters used to evaluate linearity, e.g., multiple coherence (Travassos and Beamish, 1988; Egbert and Livelybrooks, 1996; Bendat and Piersol, 2011) and bivariate coherence (Ritter et al., 1998; Weckmann et al., 2005). All of those parameters can be indicators of the data quality under the assumption that the dataset follows a linear relationship.

Multiple coherence is defined as the ratio of the ideal output spectrum due to the measured inputs in the absence of noise to the total output spectrum, which includes the noise (Bendat and Piersol, 2011). In equation form, the multiple coherence associated with E x is calculated as follows: r m 2 = 1 − E e r r x i E ¯ e r r x i E x i E ¯ x i , (11) where E e r r x i = E p x i − E x i , E p x i   a n d   E x i denote the predicted and observed electric field calculated by the i ^th segment. The bar denotes the conjugate of a complex number. Because the error between the predicted and observed electric field may be larger than the observed electric field. The right part of Eq. 11 may be a negative value. We take the square root of the absolute value of 1 − E e r r x i E ¯ e r r x i E x i E ¯ x i as the multiple coherence ( r m ), and regard the event, which r m larger than 0.8 and smaller than 1, are in high linearity. The red events in Figure 7C show the event in high linearity based on the multiple coherence.

Bivariate coherence is defined as a function of the amplitude ratio and phase difference between the predicted and the observed electric field. In equation form, bivariate coherence associated with E x is calculated as follows (Weckmann et al., 2005): r b 2 = Z x x * H x i E ¯ x i + Z x y * H y i E ¯ x i E x i E ¯ x i = E p x i E ¯ x i E x i E ¯ x i = E p x i E x i · cos θ i , (12) where E x i , H x i and H y i represent the EM field corresponding to the i ^th segment. θ i is the phase difference between E p x i   a n d   E x i . Because the predicted electric field may be larger than the observed electric field, and bivariate coherence may be larger than 1. We regard the event, which r b larger than 0.8 and smaller than 1, are in high linearity. The red events in Figure 7D show the data in high linearity based on the bivariate coherence.

The comparison of the parameters used to evaluate the linearity is shown in Figure 7. First, the events are divided into groups that contain N samples, e.g., 20 samples. We calculate the predicted electric field for each group, separately. If the data quality is high, all of the parameters should be close to one under the assumption that the data follow a linear relationship. All of the parameters can identify the high-quality data corresponding to the magnetic storm. However, the parameters ( P L c o h s z and P A R s z ) proposed in the research and multiple coherence perform better than bivariate coherence. Some parts of the high-quality events, which r b larger but close to 1 are regarded as in low linearity based on the bivariate coherence. The MT sounding curves calculated by the different preselection strategies are shown in Figure 8. All of the preselection strategies can reduce the influence of noise compared with the robust result in Figure 6. While the apparent resistivity in Figure 8B around the period of 7.9 s is downbiased compared with other results in Figures 8A, C. It may be caused by some parts of the high-quality events being removed when using the bivariate coherence to prescreen data.

The second case study uses data observed in northeastern China on 26 June 2020. Phoenix Geophysics Instruments were used to collect the MT time-series data. These data are provided by the Institute of Geophysical and Geochemical Exploration, China Geological Survey. Time-series data from 3:00 to 22:00 UTC were used in this case study. The sampling rate is 15 Hz. The observation area is in the GMT+8 time zone, and the local midnight time is approximately 16:00.

First, we examine the variation in the parameters at different frequencies. Figure 9 shows the variation in the period of 6.7 s. Red denotes events with high linearity, and blue denotes events with low linearity. The previous 2,500 events in the daytime have a preferred polarization direction of approximately −30 ° , as shown in Figure 9C, and the polarization direction becomes scattered at nighttime (the events are approximately 2,500 to 4,000). This indicates that the daytime event is dominated by coherent noise, and most of the events have high linearity; in contrast, the event at nighttime is relatively quiet. Figure 10 shows the variation in the period of 33.4 s. The previous 300 events in the daytime have low linearity, and the polarization direction is scattered. In contrast, the events during the nighttime between 300 and 500 have a high linearity, and the polarization direction is scattered. This indicates that the daytime event is dominated by incoherent noise and that the nighttime events are relatively quiet. After analyzing the parameter variation at different frequencies, we find that most of the events are dominated by coherent noise between 2 and 20 s and dominated by incoherent noise between 20 and 100 s. The field data are contaminated by coherent and incoherent noise simultaneously.

Then, we compare the MT sounding curves calculated by the different methods. First, we compare the result calculated by the SSMT-2000 and the results with and without the preselection strategy based on linearity, as shown in Figure 11. SSMT-2000 is one of the standard Phoenix software sets. After comparing all the results, we think all the impedance results are biased between 2 and 20 s, and there is a rapid rise and fall in the apparent resistivities. The SSMT-2000 result coincides with the preselection strategy result between 20 and 100 s and changes smoothly. The single-site robust result is improved between 20 and 100 s after using the preselection strategy with P L c o h s z and P A R s z . According to the data quality analysis in different periods, most of the events are contaminated by incoherent noise between 20 and 100 s, and P L c o h s z and P A R s z are effective in removing incoherent noise. Most of the events are contaminated by coherent noise, which is highly linear, between 2 and 20 s. We also try to use the remote reference method to improve the result but fail (see Supplementary Figure S5), and it needs a different strategy to suppress the noise.

Next, we try to use the information on MPD to constrain the coherent noise. Figure 12 shows the variation in the MPD and the corresponding variation in D D p o l in the periods of 6.7 and 33.4 s. The expected value of D D p o l is 1/3, and D D p o l increases when the polarization direction has a preferred direction. It shows that D D p o l is effective in differentiating the events with and without a preferred polarization direction.

We also compare the criteria proposed by Platz and Weckmann (2019) which is based on the statistical information on magnetic polarization direction (SMPD). They subdivided the polarization direction into 180 bins with a bin width of 1 ° . In general, the polarization direction is randomly distributed in each bin. Therefore, the expected value ( N u m b e r   o f   e v e n t s 180 ) is the same for each bin. They removed all the events whose polarization directions fall into bins that are much higher than the expected value. Figure 13 shows the statistical analyses of the polarization direction; the threshold k (k = 1.5σ) is used to detect abnormal values, where σ is the standard deviation (Chave and Thomson, 2003). The corresponding abnormal events are drawn in red. According to the criteria based on SMPD, the quiet event may be removed at nighttime (the event from 2,500 to 4,000) and many of the events in the daytime remain.

At last, we compare the MT sounding curves calculated by the different preselection strategies based on the information on MPD, as shown in Figure 14. The robust estimator combining P L c o h s z , P A R s z for the preselection strategy is used to calculate the blue curve, which prescreens the data only based on the linearity. The robust estimator combining P L c o h s z , P A R s z and SMPD for the preselection strategy is used to calculate the red curve. The criteria proposed by Platz and Weckmann (2019) do not improve the result, the rapid rise and rapid fall between 2 and 20 s remain. The robust estimator combining P L c o h s z , P A R s z and D D p o l for the preselection strategy is used to calculate the black curve, and the threshold for D D p o l is set to 0.5. The rapid rise and rapid fall between 2 and 20 s are removed. Comparing the two criteria between D D p o l and SMPD, most of the events in the daytime are removed based on the D D p o l , while many events in the daytime remain based on the SMPD, and those events in the daytime may also correspond to the coherent noise and dominate the regression, making the preselection strategy fail. This shows the superiority of D D p o l for detecting coherent noise with a strong magnetic polarization direction.