The application of phase difference to analysis the magnetotelluric data


 Magnetotelluric (MT) method is an electromagnetic geophysical method for inferring the earth's subsurface electrical conductivity from measurements of natural geomagnetic and geoelectric field variation at the earth's surface. The ﬁrst step in MT data processing is to estimate the impedance tensor in the frequency domain from the measured time-series data. The initial MT response function estimator is based on the least-square theory; it can be severely disturbed by the cultural noise. In the presence of a small amount of intermittent contaminated data, it can be improved by remote reference technique, robust procedure or combination of them. In the presence of a large amount of contaminated data, it can still succeed with assistance from data analysis to remove the most contaminated data before the impedance tensor estimation. The phase difference is an important parameter to analyze the data in the frequency domain. In this paper, we investigate three parameters(the predicted coherence, remote coherence and polarization direction) correspond the phase difference to analyze the MT data. We demonstrated that the high predicted coherence could indicate a high signal-to-noise ratio(SNR) or strong coherence noise. The polarization direction was useful to visualize the background noise. The remote coherence was a useful parameter to indicate the quality of the data. In this paper, we will introduce a robust M-estimator at first. At last, we showed the effectiveness of the application of remote linear coherence to the selection strategy based on the M-estimator. By this selection strategy, the result can be improved dramatically in the presence of a large amount of intermittent noise.


INTRODUCTION
Magnetotelluric (MT) method is an electromagnetic geophysical method for inferring the earth's subsurface electrical conductivity from measurements of natural geomagnetic and geoelectric field variation at the earth's surface. There is a linear relationship between the geoelectric field and the geomagnetic field in the frequency domain. The first step of MT data processing is to estimate the frequency-dependent impedance tensor from the measured time-series data. The impedance tensor at a specific frequency can be calculated in the frequency domain as follow: where E and H are the horizontal electric and magnetic field components at a specific frequency, and Z means the MT impedance. The x denotes the north-south direction; y denotes the east-west direction. All the MT data interpretation is based on the result of impedance. Therefore, it is important to get reliable impedance. The initial MT impedance estimator is based on the least-square theory (Sims et al., 1971); it can be severely disturbed by the cultural noise. In the presence of a small amount of contaminated data, it can be improved by remote reference (Gamble et al. 1979) technique, robust procedure (Egbert & Booker 1986;Chave et al., 1989;Smirnov, 2003;Chave et al., 2003;Chave and Thomson, 2004) or combination of them.
The remote reference processing can improve the estimator's performance using the cross-spectral instead of the auto-spectral when performing the regression based on the least-squares estimator (Gamble et al., 1979). Robust algorithms aim to detect and reject outliers by a data-adaptive weighting scheme. Both methods need a reasonable proportion of contaminated data to yield reliable results. In the presence of a large amount of contaminated data, the estimator can still succeed with the data analysis to remove most of the contaminated data. It is practical to use a selection strategy before impedance tensor estimation (Weckmann et al., 2005).
The impedance is calculated in the frequency domain. The phase difference is an important parameter to analyze the data in the frequency domain. In this paper, we proposed a new parameter, linear coherence, based on the phase difference. The linear coherency measures the consistency of the phase difference between the time series in the same frequency. In physics, two waves are perfectly correlated in the time domain if their frequency is identical, and their phase difference is 0°. It is used to define the remote coherence between the local and remote magnetic field and the predicted coherence between the predicted electric field and observed electric field.
In this paper, section 2 will introduce the three impedance estimator Bounded Influence Remote Reference Processing (BIRRP, Chave and Thomson, 2004), the SSMT 2000 program and a robust M-estimator. Section 3 introduce the parameters corresponding to the phase difference. Section 4 will test the M-estimator with quiet field MT data, at first, and then demonstrate the performance of the three parameters, the predict coherence, remote coherence and the polarization direction by the field data which contaminated by intermittent noise.
The field data consists of natural sources and local noise. The magnetic field observed at the different site come from the same source, magnetosphere and ionosphere. Usually, the remote magnetic field has a high correlation with the local magnetic field. That also means the phase difference between local and remote magnetic field should be close to 0° at a different frequency. The linear coherence should be close to 1. Szarka(1988) and Junge(1996) summarized active and passive sources of noise that can be observed in MT measurements. The noise signals produced by the local environment (e.g., electric transmission lines, electric fences, trains) don't correlate. Therefore, the remote linear coherence between the local and remote magnetic field is a good indicator of the data's quality.
The predict coherence strongly depends on the initial model of the least square estimator. The high coherence can indicate a high signal-to-noise ratio(SNR) or strong coherence noise. Therefore, it is not suitable to apply to every situation. Fowler(1967) proposed the polarization direction. Weckmann (2005) showed its effectiveness to visualize the background noise. In this paper, we showed the relationship of the phase difference of the two orthogonal electric field and magnetic field with the polarization direction.
Between the three parameters, the remote linear coherence is suitable for the selection strategy. At last, we will show the effectiveness of remote linear coherence to the select strategy based on the M-estimator in section 5. By this selection strategy, the result can be improved dramatically.

THE DATA PROCESSING ALGORITHM
This section will introduce three independent MT impedance estimators, the Bounded Influence Remote Reference Processing ( BIRRP; Chave and Thomson, 2004), the SSMT 2000 and a standard M-estimator, we will use them to compare the performance with each other. At last, we will show the effectiveness of the select strategy based on the M-estimator in section 5.

BIRRP
BIRRP is one kind of the typical robust estimator based on windowed FFT. It divides the original time series into the short segments at first. Each segment is whitened using an autoregressive filter, and the Fourier transforms of the sections are computed after tapering with a prolate window. The segment length is variable, and the selected frequencies within each segment are grouped to calculate the response function in a different frequency.
The first impedance estimator is based on the least-square theory (Sims et al., 1971). The objective is to minimize the sum of squares residuals, and the solution of impedance at a specific frequency is given by: (2) where E and H are the horizontal electric and magnetic field components at a specific frequency, and Z is the MT impedance. The superscript † denotes the complex conjugate transpose. The M-estimator gives a weight w to the outlier based Huber weight function depending on the residual between the output (electric field) of the least-square estimator and the observation data. And the solution of M-estimator at a specific frequency is given by: (3) M-estimator can reduce the influence of unusual data (outliers) in the electric field but are not sensitive to exceptional input (magnetic field) data, which are termed leverage points. The hat matrix is widely used to detect unusual input data. The hat matrix is N by N matrix and defined as follows: (4) Chave (2003) suggested that the hat matrix's diagonal element, which is more than several times 2/N, is problematic. N denotes the number of data. The bounded influence (BI) estimator combines the robust M-estimator with leverage weighting ( ) based on the statistics of the hat matrix diagonal element (Chave et al., 2004). And the solution of BI-estimator at a specific frequency is given by: (5) The detail of the M-estimator and BI estimator is described by Chave (2012).
The remote reference processing can improve the estimator's performance using the cross-spectral instead of the auto-spectral when performing the regression based on the least-squares estimator (Gamble et al., 1979). And the solution of remote reference based on BI-estimator at a specific frequency is given by: where denotes the remote reference magnetic field data. BIRRP is an abbreviation for the Bounded Influence Remote Reference Processing method.
The error bar is calculated by jackknife at a specified frequency. The jackknife is a resampling method used to estimate the bias of a large population. It involves a leave-one-out strategy to estimate the error bar for the data set.

SSMT2000
SSMT-2000 is one of the standard Phoenix software set. This program takes the raw time series files, calibration files, and site parameter files as inputs. In an intermediate step, it produces Fourier coefficients, which are then used to calculate the impedance tensor by the robust routine that is similar to BIRRP we introduce before.
Before impedance estimation, coherency and resistivity variance, are used to filter out the noisy data. Resistivity variance gives more weight to data points with smaller error bars. Ordinary coherency gives more weight to data points with good coherency between E and H channels. This process reduces the size of the error bars and smooths the curves in plots of apparent resistivity.

M-estimator 2.3.1 Spectrum estimation
The cascade decimation proposed by Wight and Bostick (1980) was widely used to calculate the spectrum for a different frequency. This scheme has a significant advantage: it computes the spectrum fast, produces a spectrum representing constant percentage bandwidth, and is evenly spaced on a log-frequency scale. It divides the original time series into segments with the specific window length (for example,32,128,256, or 512). The selected frequencies within each segment are grouped to calculate the impedance tensor in the frequency domain. Data are then passed through a cascaded sequence of low pass filters and downsampled by parameter 2. The sampling frequency is reduced by two whenever the level raises one stage. The same number of data points within each segment is extracted to calculate the impedance tensor at a different frequency.
In this study, a way similar to BIRRP is adopted to estimate the spectrum. Instead of downsampling the original time series. The segment length is variable with each level. The 5 th and 8 th number of each segment's frequency coefficient is grouped to calculate each level's response function. The section length is then reduced by 2, and the process is repeated to calculate the impedance tensor at a different frequency.
For the real MT time series, , , ℎ , ℎ , ℎ . The unit of the electric field is transformed into mv/km; the magnetic field is transformed into nT. For each level, the procedure to produces the accurate spectrum for each section as follow: (1) Using the first difference filter to remove the long-period trends and means, in the case of , yielding the new series: (2) To reduce the bias of spectral estimation, a multi-tapering method is adopted. Multitapering process is an extension of single taper approaches. The multi-tapering method can be carried out by the Matlab intrinsic function "dpss.m". In this study, we specify the half bandwidth product to be 2.5 and multi-tapering to be 4. It produces four orthogonal discrete prolate spheroidals (Slepian) sequences windows and a column weight vector (lambda), the lambda equal in the length to the number of Slepian sequences. Tapper the 4 Slepian sequences to each segment and get four time-series sets .
(3) Fourier transform is carried out to each , and calculate the weighted mean by multiplying the lambda to the four Fourier coefficient sets, yielding the Fourier coefficients [ ], where i=1, 2,..., N/2 is the frequency number. The other components ( , , , ) perform the same procedure. The 5 th and 8 th number of the frequency coefficient within each segment is grouped to do further processing at each level.

Robust estimation for the response function
The regression step is similar to the method that Neukirch et al. (2014) and Campanya et al. (2014) described. The impedance tensor is calculated as follows : where the n denotes the number of the data, denote the ℎ real part of , denote the ℎ imaginary part of ; it is the same to , , denotes the real part of , denotes the imaginary part of ; it is the same as . The symbol j denotes the imaginary number unit. The total equations for all data can be transformed into matrix form as follows: The bivariate complex linear regression problem converts to the real multivariate linear regression problem. The Matlab intrinsic function "robustfit.m" can easily solve the regression problem, and Huber weight is adopted. We also can solve the by Eq. 3, following the M-estimator. They are the same in theory; both of them give a wight to the outlier in the electric field. We have compared the two methods. The result is almost the same.
When x equals the remote reference N by two matrices [ , ]. The interstation transfer function is calculated at first as follow: (11) Then impedance tensor is calculated using the relationship between the local transfer function and interstation transfer function as follow: (12) It can be categorized into one kind of two-stage processing strategy. In this way, separating the calculation into the regression between the remote reference with all channels can avoid a direct effect of coherent noise between local channels. If the remote magnetic field does not contain correlated noise with the local site, it will work well.
Finally, the robust regression was bootstrapped to compute the data-dependent distribution of impedance values and estimate the intrinsic data errors. Empirically, 1,024 iterations were repeated, considering a sufficient trade-off between accuracy and computation time to evaluate results' uncertainty.

THE APPLICATION OF PHASE DIFFERENCE TO MT DATA ANALYSIS
The coherency is a quantitative measure of the consistency of the phase difference(PD) between the two waves. If two sources are coherence, their phases must be either the same or have a constant difference (Marple and Marino, 2004). We divide the two wave into N individual segments. It is defined as the ratio of the mean cross power to the root mean auto powers. For A and B, it is defined as: where the brackets represent the averages of N individual auto-power and cross-power density spectra.
In the polar coordinate, there is: Where and is the i th spectrum in a specific frequency.
denotes the phase of and , respectively. The superscript † denotes the complex conjugate transpose.
The phase difference between two segments at a specific frequency is defined as follow: the denote the angle of the phase difference between two segments.
Suppose A and B are not coherent with each other. The different cross power will counteract with each other. Therefore, the value of coherency will decrease.
There is a particular situation that there are constant PD equal zeros. At this condition, the time series also have a linear relationship in the time domain. This paper will discuss the application of phase difference in this situation to the MT data analyze.

The linear coherence
In this section, a new parameter, linear coherence, is proposed to express the phase difference in another way. The linear coherence is defined as follow: The Lcoh is the linear coherence and denotes the cosine of the phase difference. As shown in Fig.1 Their frequencies should be identical, and their phase difference equals 0°. For a specific frequency, the higher linear relationship in time domain also means a smaller phase difference. The bigger the phase difference is, the small the linear coherence is. Fig.1 The variation of the cosine functions with a different phase difference.

The phase difference between the orthogonal electric and magnetic field
By Faraday's law, a time-varying magnetic field will induce an orthogonal electric field depending on the media's resistivity. In a simple 1D model, E=ZH. the apparent resistivity is defined as: The phase is defined as: φ=arctan( ImZ ReZ ) .
where T denotes the period. φ denotes the phase. Z is a complex number. The phase does not equal to 0° or 180°. The phase difference between the orthogonal electric field and the magnetic field determines the impedance phase. And the amplitude ratio between the orthogonal electric field and the magnetic field determines the apparent resistivity. For the north-south direction, the phase difference between the orthogonal electric field and the magnetic field is calculated as follow: the linear coherence between the orthogonal electric field and the magnetic field(EHLcoh) is defined as: As the impedance correlates with the conductivity distribution of the underground. It will not change with time usually. The EHLcoh also should be stable.

The remote coherence
The natural MT source comes from magnetosphere and ionosphere, which is far enough from the local site that following the assumption of the plane wave. All other parts of the measured electric and magnetic fields are considered as noise. The field data consists of natural sources and local noise. Usually, the remote reference magnetic field has a high correlation with the local magnetic field. Szarka(1988) and Junge(1996) summarized active and passive sources of noise that can be observed in MT measurements. The noise signals produced by the local environment (e.g., electric transmission lines, electric fences, trains) doesn't correlate with each other. Therefore, the coherence between the local and remote magnetic field is a good indicator of the data's quality. For the north-south direction, the linear coherence between the remote and local magnetic field( Lcoh) is defined as follow: Where H x_i and H xr_i is the i th local and remote magnetic field spectrum at a specific frequency. The Lcoh should have high linear coherence. The linear coherence should be close to 1.

The phase difference between the predicted electric field and observed electric field
We also use the phase difference defined the predict coherence between the measured electric field(E) and electric field predicted(E p ) from E p =ZH, where Z is obtained by LS estimator. The predicted coherence(PLcoh) is defined as follow: As the predicted electric field and the observed electric field should be similar. The PLcoh should have high linear coherence which is close to 1.

The polarization directions
The polarization directions describe the time-harmonic variation of two orthogonal field components with a constant phase difference. The polarization directions of the electric field ( _ ) and magnetic field ( _ ) (Fowler et al., 1967) at a specific frequency are defined as: We can rewrite the function as: Where and are _ and _ or _ and _ , respectively. The polarization direction is related to the phase difference and the amplitude ratio between the two orthogonal fields. As a variety of sources generates a natural magnetic signal. These sources generate magnetic fields that vary in their incidence directions, and thus there is no preferred polarization direction for the magnetic field. However, according to a given conductivity distribution of the subsurface, there might be a preferred polarization direction of the induced electric field (Weckmann et al., 2005).

DATA ANALYSIS
In this section, we will test the M-estimator's correctness with quiet field data and investigate the performance of each parameter by the field MT data, which is contaminated by intermittent noise. For the figure correspond to linear coherence, the upper shows the north-south component, the below shows the east-west component. Blue denotes a negative value, and red represents a positive value. The black curve represents the absolute value of the log of the amplitude ratio between the corresponding component.

Field Data Set 1: quiet data
In this section, we tried to test the performance of the M-estimator by the quiet field data observed in the desert. The data is observed by the MTU-A Geophysical Instruments. Fig. 2 shows the MT sounding curve calculate by SSMT 2000 and M-estimator. SSMT 2000 calculates the blue curves. M-estimator calculates the black curves. Both of them coincide with each other. It shows the correctness of the M-estimator and also confirms the correctness of the spectrum estimation. All the parameter estimation will base on the spectrum estimation introduced in 2.2.1.

Field Data Set 2: KAP03
The Kaapvaal 2003 (KAP03) long period 5-component magnetotelluric time-series data was used in this paper. Data were recorded with a 5-second sampling rate for around a month in 2003 at each site using GSC LIMS systems as a part of the SAMTEX project. The 26 long-period sites distributed in a NE-SW profile are shown in Fig. 3. Data at sites in the middle of the profile (KAP127 to KAP145) were heavily contaminated by DC signals due to a power line for trains running between Kimberley and Johannesburg.  The sampling rate is 5 seconds. The unit of the electric field is mV/km, and of the magnetic field is nT. The blue line shows the Dst (disturbance storm time) index. The Dst index is a negative index of geomagnetic activity used to estimate the averaged change of the horizontal component of the earth's magnetic field based on field measurements. It is expressed in nanoteslas. When the Dst index is less than -50 nT, it is categorized as a geomagnetic storm. When the Dst index is less than -100 nT, it is categorized as a strong geomagnetic storm. This dataset includes one geomagnetic storm from Oct. 29 to Oct. 31. During the geomagnetic storm, the data followed the plane wave's assumption, and the variation was much stronger than the data during non-storm periods. That means the single-noise ratio (SNR) is better than the non-storm day data. Fig. 5 shows the polarization direction for the electric field and magnetic field at 84-second from Oct. 26 to Oct. 31. There is preferred polarization direction in the magnetic field between the non-storm day(Oct. 26 to Oct. 29), and become scattered between the geomagnetic storm day(Oct. 29 to Oct. 31). On the contrary, The polarization direction for the electric field is scattered between the non-storm day and has a preferred direction during the geomagnetic storm day. The result shows that the SNR is low between the non-storm day and is high between the storm day. Fig. 6 shows the MT sounding curves calculated using the data on both storm and non-storm days at site 133. BI-estimator was used to calculate the result. The black curves were calculated using the non-storm data from 00:00:00 on Oct. 26 to 00:00:00 on Oct. 29. The red curves were calculated using data from 00:00:00 on Oct. 29 to 00:00:00 on Oct. 31. The upper figures show the apparent resistivity, and the lower figures show the impedance phase. The XY phase calculated by non-storm day data is close to 0º, and the apparent resistivity increases as a line in the log scale. The impedance curves are improved by using the geomagnetic storm data, and the result became more reasonable than the curves using non-storm data. This result also coincides with the result of the polarization analysis that the SNR is low during the non-storm day and is high during the geomagnetic storm day. Fig. 4 Time variations of magnetotelluric data whose sampling rate is 5 seconds at site 133. The red vertical lines show the data gaps, and the black lines show the 5-component MT data. The blue line shows the Dst index of magnetic data. The unit of the electric field is mV/km, and it of the magnetic field is nT. Fig. 5 The polarization direction for the electric field and magnetic field at 84-second. Fig.6 The MT sounding curves calculated using the data on both storm and non-storm days at site 133. BI-estimator was used to calculate the black curve. The upper figures show the apparent resistivity, and the lower figures show the impedance phase. Fig. 7 shows the variation of RLcoh versus the amplitude ratio at 84-second. Site 151 is set as the remote site. The black curve denotes the amplitude ratio between the local and remote site. The RLcoh becomes high during the geomagnetic storm period, which coincides with the SNR is high between the geomagnetic storm period. Fig. 8 shows the variation of PLcoh at 84-second. The PLcoh is low and scatted between the non-storm day. The PLcoh become high during the geomagnetic storm period.
Here the high SNR data during the geomagnetic storm data correspond to the high predicted coherence. Fig. 9 shows the variation of EHLcoh versus the amplitude ratio at 84-second. The black curve denotes the amplitude ratio between the orthogonal electric field and magnetic field. This figure shows the different statistical property of the phase difference and amplitude ratio. So the result calculated by the data between the non-storm day and storm day is systematically different. Fig. 7 The variation of RLcoh versus the amplitude ratio at 84-second. The black curve denotes the absolute value of the log of the amplitude ratio. Fig. 8 The variation of PLcoh versus the amplitude ratio at 84-second. Fig. 9 The variation of EHLcoh versus the amplitude ratio at 84-second. The black curve denotes the absolute value of the log of the amplitude ratio.

Field Data Set 3: L6-7
In this section, the broadband-frequency MT time-series data observed by the Phoenix MTU-A Geophysical Instruments is used. This data belongs to the Institute of Geophysical and Geochemical Exploration, China Geological Survey. In this data set, the high-frequency band (2,400 Hz) sampled 1 second at intervals of 2 minutes from the beginning of the minute, the middle-frequency band (150 Hz) sampled 8 seconds at intervals of 2 minutes from the beginning of the minute, and the low-frequency band (15 Hz) sampled continuously. Fig.10 shows the location map of sites in the study area. The local survey lines are L6-1, L6-2 and L7-1, and the Y0625 is set as the remote reference site. The lower left inset map shows the survey area in China. The upper left insert map shows the detail of survey line L7-1; The lower right insert map shows the detail of survey line L6-1; The upper right insert map shows the detail of survey line L6-2.   11 shows the Hx component time series measured at the same time from 3:00:00 to 22:00:00 on June 26 in different sites. The sampling rate is 15 Hz. The variation in the daytime is much stronger than that at midnight at site L6274 and L6276. The variation of the time series at other sites is almost keeping the same scale. Fig.12 shows the detail of the Hx component in Fig.11. The left figures show the data in the daytime. The right figures show the data at midnight. Each segment is ten munites. The unit of the magnetic field is nT. The variation at site L6274 and L6276 is about 100 times than other sites in the daytime and similar scale at midnight. That means the noise to signal ratio is 100 at the daytime. It will bias the impedance tensor a lot. As the influence of local noise, the correlation of the magnetic field becomes low in the daytime, and the correlation becomes higher at midnight. Fig.11 Time series of Hx component whose sampling rate is 15 Hz at the same time from 3:00:00 to 22:00:00 in UCT time. The unit of the magnetic field is nT. Fig.12 Time-series data of Hx components whose sampling rate is 15 Hz at the same time. The left figures show the daytime data, and the right figures show the night data. Each data segment is for ten munites. The unit of the magnetic field is nT.
As the signal strength at the dead band (1-10s) is very low, it is easy to be influenced by the local noise. The typical polarization directions for the dead band is shown in Fig.13. This figure shows the result of polarization directions calculated by the site L7158 at 8-second. The polarization directions of the magnetic field have a preferred direction from 3:00:00 to 15:00:00. That means the data in that period is contaminated by strong noise. A similar situation occurs at all sites, including the remote site. Fig.14 shows the magnetotelluric sounding curves calculated using the data at site L7158. SSMT-2000 was used to calculate the blue curves from one-day data, BIRRP was used to calculate the black curves using daytime data from 16:00:00 to 20:00:00 in UTC on June 26. M-estimator was used to calculate the red curves from midnight data from 16:00:00 to 20:00:00 in UTC on June 26. The data calculated by the one-day data is biased between the dead band. And the result calculated by the midnight data changes reasonable.
All the results from Fig. 11 to Fig.14 indicate that the data at midnight is quiet. Fig.15 shows the variation of RLcoh at 84-second versus the amplitude ratio at site L7158. The property of the amplitude ratio is quite different between the daytime and midnight time. The RLcoh is low or negative in the day time. It becomes high at midnight. It shows that RLcoh is a good indicator of the data's quality. Fig.16 shows the variation of PLcoh at 84-second. The PLcoh is high in the daytime. PLcoh becomes low or negative at midnight. In this situation, the noise data in the day time correspond to the high predicted coherence. Fig.17 shows the variation of EHLcoh at 8-second versus the amplitude ratio at site L7158. This figure shows the different statistical property of the phase difference and amplitude ratio between the daytime and midnight. So the result calculated by the data between the day time and midnight time is systematically different, as shown in Fig.14.   Fig.13 The polarization direction for the electric field and magnetic field at 8-second. Fig.14 The magnetotelluric sounding curves calculated using the data at site L7158. SSMT-2000 was used to calculate the blue curves from one-day data, BIRRP was used to calculate the black curves using daytime data from 16:00:00 to 20:00:00 in UTC on June 26. M-estimator was used to calculate the red curves from midnight data from 16:00:00 to 20:00:00 in UTC on June 26. Fig.15 The variation of RLcoh at 84-second versus the amplitude ratio at site L7158. The content is the same as Fig.  8. Fig.16 The variation of PLcoh at 84-second versus the amplitude ratio at site L7158. The content is the same as Fig.  9. Fig.17 The variation of EHLcoh at 8-second versus the amplitude ratio at site L7158. The content is the same as Fig.  10.

SELECTION STRATEGY BY REMOTE COHERENCE
From the performance of the three parameters (the predicted coherence, remote coherence and polarization direction) we found that the high predicted coherence can indicate a high signal-to-noise ratio(SNR) or strong coherence noise. Therefore, it is not suitable to apply to every situation. The polarization direction is effective to visualize the contaminated data, but it is not suitable to use to the selection strategy. The coherence between the local and remote magnetic field is a good indicator of the data's quality. When the SNR is high, the remote linear coherence will close to 1. It is suitable to use as the selection scheme.
In this section, we will show the performance of applying the RLcoh to the selection strategy based on the M-estimator introduced in section 2.3. we will use the data from L6-7. We set the threshold of LRcoh to be 0.9. The low remote coherence data will be removed and not participate in further processing. Fig.18 shows the magnetotelluric sounding curves calculated using the data at site L7158. All the result is calculated by one-day data. SSMT-2000 calculated the blue curves; M-estimator calculated the red curves with the remote reference data. M-estimator calculated the black curves with selection scheme by remote linear coherence. The result gets by the selection scheme is improved dramatically.
The similar result can get at site L6174, L6178, L6180, and L7154. Fig.19 shows The magnetotelluric sounding curves calculated using the data at site L6174. Fig.20 shows the magnetotelluric sounding curves calculated using the data at site L6176. Fig.21 shows the magnetotelluric sounding curves calculated using the data at site L6178. Fig.22 shows the magnetotelluric sounding curves calculated using the data at site L6180. Fig.23 shows the magnetotelluric sounding curves calculated using the data at site L7154. Fig.18 The magnetotelluric sounding curves calculated using the data at site L7158. All the result is calculated by one-day data. SSMT-2000 was used to calculate the blue curves; M-estimator was used to calculate the red curves with the remote reference data. M-estimator was used to calculate the black curves with selection scheme by remote linear coherence. Fig.19 The magnetotelluric sounding curves calculated using the data at site L6174. The content is the same as Fig. 19. Fig.20 The magnetotelluric sounding curves calculated using the data at site L6176. The content is the same as Fig. 19. Fig.21 The magnetotelluric sounding curves calculated using the data at site L6178. The content is the same as Fig. 19. Fig.22 The magnetotelluric sounding curves calculated using the data at site L6180. The content is the same as Fig. 19. Fig.23 The magnetotelluric sounding curves calculated using the data at site L7154. The content is the same as Fig. 19.

DISCUSSION
The initial magnetotelluric(MT) response function estimator is based on the least-square theory (Sims et al., 1971); it can be severely disturbed by the cultural noise. In the presence of small parts of intermittent contaminated data, it can be improved by remote reference technique (Gamble et al. 1979), robust procedure (Egbert & Booker 1986;Chave et al., 1989;Smirnov, 2003;Chave and Thomson, 2004) or combination of them.
The remote reference processing can improve the estimator's performance using the cross-spectral instead of the auto-spectral when performing the regression based on the least-squares estimator (Gamble et al., 1979). Robust algorithms aim to detect and reject outliers by a data-adaptive weighting scheme. Both methods need a reasonable proportion of contaminated data to yield reliable results. In the presence of a large amount of contaminated data, the estimator can still succeed with assistance from data analysis to remove the most contaminated data before the impedance tensor estimation. It is practical to use a selection strategy before impedance tensor estimation (Weckmann et al., 2005).
To recover the digital data information, we should sample at least two points in one period, followed by the sampling theory. But to get the accurate complex coefficient from the time series. We suggest that it is better to sample 4 points in one period, and one continuous segment should contain at least four times the period. For instance, the 32-second continuous time series with the 1-second sampling rate. We can get a reliable complex coefficient between 4 to 8 second.
To get a reliable complex coefficient of 1000-second. We need about 4000 seconds time series segment at least. To keep the independence of each data and get more sample data. The overlay rate set to be 50%. If we assume, 20 data sets are enough to get a reliable result from the regression. We need about 40,000 second length time series. For a difficult period, the length of the time series we need is different. For a shorter period, the length is shorter. The length for a different period is shown in table 1. Table 1 The enough length of the time series for a different period Period 10 s 10 2 s 10 3 s The length of time series 400 s 4,000 s 40,000 s We set there are about 20 samples at the first level. When we do the decimation at each level, we can get double samples than the last level. We will get many samples than we need for the high-frequency data. If any data is problematic, we can remove it directly rather than give a small weight to them. Especially, the data between the dead band (1-10s), as the signal strength is very low; it is easy to be influenced by the local noise. And its enough length is about 400 second(about 6 minutes). Usually, the observation period is one day. When the intermittent noise contaminates the dead band data, we can apply the selection strategy to get a reliable result.
Moreover, the continuous noise causes more severe problems. It isn't easy to get a reliable result from current impedance estimator directly. We have to modify the time series before we transfer it into the frequency domain, just like the method noise identification and modifying (Kappler, 2012;Wang et al., 2017;Li et al., 2019).

CONCLUSIONS
In this paper, a robust iteratively reweighted least squares regression with the Huber weight function(M-estimator) is introduced at first. We proposed a new parameter, linear coherence, based on the phase difference. The linear coherency measures the consistency of the phase difference between the two waves. It denotes the linear relationship between the two waves in the time domain. For a specific frequency, the higher the linear coherence is, the higher the linear relationship in the time domain. It is used to define the remote coherence between the local and remote magnetic field and the predicted coherence between the predicted electric field and observed electric field. We also showed the relationship between the phase difference with the polarization direction. We demonstrated the performance of the three parameters, the predict coherence, remote coherence and the polarization direction by two field data.
Usually, the remote reference magnetic field has a high correlation with the local magnetic field, and artificial noise doesn't correlate with each other. Therefore, the remote coherence between the local magnetic field and remote reference magnetic field is a good indicator of whether the noise is strong or not.
The predict coherence strongly depends on the initial model of the least square estimator. The high coherence can indicate a high signal-to-noise ratio(SNR) or strong coherence noise. So it is not suitable to apply to every situation.
The polarization direction is effective to visualize the strong noise. But it is not suitable to apply to the selection strategy.
At last, We showed the effectiveness of using the remote coherence to the select strategy based on the M-estimator. By this selection strategy, the result can be improved dramatically.

Competing interests
We know of no conflicts of interest associated with this publication. We declare that this manuscript is original, has not been published before and is not currently being considered for publication elsewhere.