The strategy to calculate impedance using RMHHT estimator for continuous and discontinuous broadband magnetotelluric time-series data

: MT method is widely used in exploration surveys worldwide but hardly applicable in urban areas because of the large amount of artificial electromagnetic noise. The MT time-series data at night is much quieter than in the daytime because most of the electric equipments are shut down for several hours during nighttime. Therefore, we focused on calculating the MT impedance using the quiet time-series data. In this research, the data observed by the Phoenix System was used. The data contain continuous and discontinuous time-series data. As an alternative way, we introduced a robust impedance estimator based on the Hilbert-Huang transform (RMHHT). We indicated that this technique needs 4-hour data to get a reliable impedance up to 1,000-second in the numerical simulation. This short measurement time makes it possible to carry out MT surveys in urban areas with strong noises. This paper demonstrated a strategy to process the broadband MT time-series data properly by the RMHHT estimator. Finally, a successful case study, using the nighttime data to get reliable impedance tensor in the areas contaminated by heavy noises, was demonstrated.


INTRODUCTION
The 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. MT method is widely used in exploration surveys around the world. Still, along with urban constructions, artificial disturbances to the electromagnetic observations are becoming more serious. MT exploration is hardly applicable in urban areas with a large amount of artificial electromagnetic noise.
The field data consists of natural sources and local noise. Szarka(1988) and Junge(1996) summarized active and passive noise sources observed in MT measurements. When the data segment contains heavy noise, it is difficult to differentiate between the MT signals coming from the natural sources and the noise signals produced by the local environment (e.g., electric transmission lines, electric fences, trains). If most of the time-series data contain strong noise, any impedance estimator will fail to get reliable results. But in one day, there might be several hours' time window at nighttime that most of the electric equipment is shutdown. The time-series data in the period is much quieter than in the daytime. In this paper, we consider using the quiet period data only to calculate the impedance tensor.
The first step in MT data processing is to estimate the frequency domain's impedance tensor from the measured time-series data. In conventional MT data processing routines, the original time-series is divided into short segments and then transforms the segments into the frequency domain using the Fourier Transform. Here we should know one segment is one dataset when doing the impedance estimation in the frequency domain. MT time-series data occasionally contain transient noise, which is not suitable for the basic requirements of conventional methods based on the Fourier transform. In the Fourier transform, a part of the time-series data contaminated by transient noise will influence the whole spectrum.
The purpose of impedance estimation is to get the impedance tensor that depends on frequency.
We also can transform the time-series into the frequency domain data by the time-frequency transform technique like HHT. Huang et al. (1998) introduced empirical mode decomposition (EMD) in the framework of the Hilbert-Huang Transform (HHT). This novel time-series analysis tool is data-adaptive and suitable for nonlinear and nonstationary data. The time-series data is decomposed into intrinsic mode functions (IMFs) that can be represented as a mono-component, containing a single frequency. Huang et al. (1998) argued that the IMF allowed for a meaningful computation of its instantaneous parameter (IP) with the Hilbert Transform. Chen et al. (2012) used the IPs for impedance estimation. Neukirch and Garcia (2013) firstly proved that IPs could be used directly for impedance estimation in theory. Neukirch et al. (2014) also showed that when the signal source is nonstationary (NS), the method based on the HHT performs better than the traditional method based on windowed FFT. The estimator based on Hilbert-Huang Transform can estimate MT response functions even in NS signal and transient noise.
In this paper, the data sets observed by the MTU-A Geophysical Instruments was used. The time-series data are stored in three files. Two of them store the high and middle-frequency bands (2,400 and 150 Hz) discontinuously. The other file stores the low-frequency data (15 Hz) continuously. As an alternative, we introduced a robust impedance estimator based on the Hilbert-Huang transform (RMHHT) and demonstrate a proper strategy based on RMHHT to calculate the full-scale broadband MT time-series data, including the continuous and discontinuous data. The property of the HHT, which could express the spectrum variation in the time-frequency domain well, makes it possible to directly connect the discontinuous time-series and process the combined discontinuous time-series in the same way as the continuous timeseries. The biggest problem of the estimator based on HHT is time-consuming. By this strategy, the computation time can be shortened much.
In the synthetic data test, we have shown that there is continuous 4-hour noise-free time-series data. We could get a reliable impedance until 1,000-second. This strategy is suitable to get a reliable result by selecting the data in the time domain. That makes it possible to use the nighttime data only to get a reliable impedance.
Finally, we compared the performance of RMHHT with the conventional method Bounded Influence Remote Reference Processing (Chave and Thomson, 2004) and the SSMT 2000 program and showed the successful case study using the nighttime data to get a reliable result from the areas contaminated by heavy noisy.

METHOD
In this section, three different impedance tensor estimators, BIRRP, SSMT-2000, and RMHHT, are introduced to compare them using the field data in subsection 4.3.

BIRRP
Bounded Influence Remote Reference Processing (Chave and Thomson, 2004) is a typical conventional robust estimator to calculate the impedance tensor based on windowed FFT.
The first estimator (least-squares estimator, Sims et al., 1971) divides the original time-series into the short segments and then transforms the segments into the frequency domain using Fourier Transform. Finally, the impedance tensor at a specific frequency can be calculated as follow: The objective is to minimize the squares residual sum by the least-squares theory, and the solution of impedance tensor at a specific frequency is given by: where E and H are the horizontal electric and magnetic field components at a specific frequency, and Z is the MT impedance tensor. The superscript † denotes the complex conjugate transpose.
The M-estimator (Egbert, 1986) gives a weight w based Huber weight function to the outlier depending on the residual between the output (electric field) of the least-square estimator and the observation data. The solution of M-estimator at a specific frequency is given by: 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 bounded influence (BI) estimator combines the robust M-estimator with leverage weight ( ) 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: 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 crossspectral 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.

SSMT-2000
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. It produces Fourier coefficients in an intermediate step to calculate the impedance tensor using the robust routine.
Before impedance estimation, coherency and resistivity variance are used to filter out the noisy data. Resistivity variance gives a larger weight to data points with smaller error bars. Ordinary coherency gives a larger 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. It gives error estimates that are systematically smaller than other methods, especially at high frequencies.

RMHHT
The purpose of impedance estimation is to get the impedance tensor depends on frequency. We also can transform the time-series into the frequency domain data by the time-frequency transform technique like HHT.
The impedance tensor derived from Eq. 1 is defined purely in the frequency domain and cannot be used for the instantaneous spectra obtained by the HHT. Neukirch (2013) showed that the instantaneous parameters (IPs) obtained by the HHT could be used to calculate the impedance directly. The impedance is calculated in the time-frequency domain as follows: and where ℎ ( ), ℎ ( ), ℎ ( ), ℎ ( ) denote the instantaneous values of the analytic signal corresponding to the electromagnetic field, which has the same instantaneous frequency ω at the same time point. denotes the n th time point. The symbol † denotes the complex conjugate transpose. N is the total number of sampling points. The procedure of RMHHT method is as follows: (1) Transforming the time-series data into the frequency domain data by HHT (2) The least-squares regression with the Huber weight function is repeated to estimate MT impedances using the IPs obtained by HHT.
The detail of the RMHHT is as follow: Step1: Time-frequency transformation Huang et al. (1998) proposed the Hilbert-Huang transform (HHT). It is the combination of the empirical mode decomposition (EMD) and the Hilbert spectral analysis.
EMD is a technique to adaptively decompose a signal into a finite number of intrinsic mode functions (IMFs), with a trend using a process called the sifting algorithm. IMFs are time-varying mono-component functions, which contain a single frequency. To obtain meaningful IPs, the IMFs have to satisfy the following two conditions (Huang et al., 1998): (1) in the whole dataset, the number of extrema and the number of zero-crossings must either equal or differ at most by one (2) at any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero Based on these two conditions, a signal x(t) is decomposed into IMFs through an iterative process termed sifting (Huang et al., 1998) that is described as follows: (1) Identify the local extrema(maxima and minima) of x(t).
(2) Interpolate all local maxima and local minima by cubic splines to obtain upper and lower envelopes, respectively.
(3) Subtract the mean m(t) of the upper and lower envelopes from x(t) to obtain the signal h(t) = x(t) − m(t).
(4) If h(t) fulfils the two conditions of IMF, the first IMF c 1 (t) = h(t). Otherwise, steps 1-3 are applied on h(t) (instead of x(t)) and are repeated i-1 times until h i = h (i−1) (t) − m i (t) fulfills the conditions of IMF. c 1 (t) = h i (t) is considered as the first IMF.
(5) Subtract the first IMF c 1 (t) from x(t) and treat the obtained signal as a new time-series x(t) = x(t) − c 1 (t).
(6) Steps 1 to 5 are repeated to determine the subsequent IMFs from the remaining signal x(t).
The decomposition process is stopped after n iterations when the residue becomes a monotonic function, which can not be decomposed further. The final signal can then be expressed as a sum of all IMFs and the residual.
The second step of HHT is the Hilbert transform, the complex conjugate y(t) of any real-valued function x(t) can be determined as: Then the time complex number can be created as: where a(t)=√ (t) 2 + (t) 2 , θ(t)=arc tan( ). Here, a(t) is the instantaneous amplitude, and θ is the instantaneous phase function. The instantaneous frequency is simply defined as: Finally, we apply the Hilbert transform to each IMFs, and we can get the IPs in the timefrequency domain. Huang et al. (1998) presented the HHT to univariate data, but MT data consists of at least four data channels depending on each other. Rehman and Mandic (2009)  scale for all channels. Neukirch (2014) showed that the instantaneous frequencies (IFs) obtained from IMFs with a similar scale can be considered the same. We will extract the instantaneous parameter in the same frequency band based on the IFs later.
There are different ways to derive instantaneous frequencies (IF) and instantaneous amplitudes (IA) from IMFs. Chen et al. (2012) discussed applying the direct quadrature method to MT data and showed that this method is more accurate than Hilbert transform. A detailed description of the method is given by Huang et al. (2009). The direct quadrature method is implemented by an empirical AM-FM decomposition based on a spline fitted normalization scheme to produce a unique and smooth empiric envelope(AM) and a unity-valued carrier(FM).  The blue line denotes the original data. The red line denotes the data's absolute values. At first, the empiric envelope(en(t)) is obtained from the maxima of the data's absolute values by spline fitting (black line in Fig. 1), the AM part is defined first. It is exactly A(t) = e(t), and then the envelope is used to normalize the data by F(t) = x(t)/en(t), where F(t) is the normalized data. This procedure is repeated until all of the F(t) extrema are unity value (grey line in fig. 1). F(t) is then designated as the empirical FM part of the data. The mono-component signal can be expressed as Its quadrature part is simply defined as x q (t) = A(t) sinφ(t); and The phase angle is computed by θ(t) = arctan (F(t)/ (1 -F 2 (t))). The instantaneous frequency can be calculated for each IMF by eq. (10). For the complex time-series We will use the instantaneous parameters to calculate the impedance tensor.
Step 2: Extract the IPs In the paper, we depend on the EMD's (Neukirch et al. 2014) spectrum estimation and extraction. Before the regression, the independent instantaneous parameters calculated by HHT were gathered, and the same frequency band's IPs were extracted to calculate the impedance tensor. Neukirch et al. (2014) describe how to extract the IPs for the different frequency bands in detail.
We divide the frequency into different frequency bands. As the impedance tensor changes smoothly with frequency, we can group the IPs based on the IFs in the same frequency band to do further processing at different frequency bands.
Step 3: Selection scheme by hat matrix Before the estimation of impedance, the leverage point should be removed. The hat matrix is N by N matrix and defined as follows: Chave (2003) suggested that the hat matrix's diagonal element, which is more than several times 2/N, is problematic. As the noisy data is energetic, the hat matrix's diagonal element's statistical analysis is useful to detect the noisy data. N denotes the number of data. Before the regression, the data where the hat matrix's diagonal element is more than two times 2/N should be removed by the selection strategy to remove the leverage point.
Step 4: Estimation of the impedance tensor 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=1,2,…, 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 Mestimator. Both methods give a weight to reject the outlier in the electric field. They are the same in theory.
When x equals the remote reference N by two matrices [ , ], the interstation transfer function is calculated at first as follow: Then impedance tensor is calculated using the relationship between the local transfer function and interstation transfer function as follow: This strategy 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 sampling rate is 1,000 Hz, and each segment is 4 seconds. Both of the time segments are added with 10% Gaussian noises. The combination of the two segments is shown in Fig. 2

. The
HHT spectrum is shown in Fig. 3. The color denotes the value of 10·log10(amp/max), where "amp" denotes the amplitude of the IPs, "max" denotes the maximum amplitude of the IPs. Fig.3 showed that the HHT could express the spectrum variation in the time-frequency domain very well.
In this section, we showed the capability of the HHT to express the spectrum variation in the time-frequency domain. Therefore, we could combine the two-segment time-series directly, which are in a different property, and get a reliable spectrum by HHT. Because of this kind of property, we could combine the discontinuous time-series data directly and process the combined discontinuous time-series data in the same way as the continuous time-series data when we compute the impedance by RMHHT. Fig. 2 The time-series x 1 (t)=sin(2π·50t)+sin(2π·200t) and x 2 (t)=sin(2π·25t)+sin (2π·100t) +sin(2π·250t). Both of the segments are added with 10% Gaussian noises. Fig. 3 The HHT spectrum of time-series data, which is shown in Fig. 1. The color denotes the value of 10·log10(amp/max), where "amp" denotes the amplitude of the IPs, "max" denotes the maximum of the amplitude of the IPs.

The Test with Different Length of Synthetic Time-series data
In this subsection, we tried to compare the performance of RMHHT by different lengths of synthetic MT time-series data. The way to create synthetic time-series is similar to the method proposed by Chen (2012). The 1-day magnetic time-series data from Memambetsu (MMB) station was used as the synthetic magnetic time-series data. MMB is one of the Magnetic Observatory stations where geomagnetic and geoelectric observations are performed in Japan.
The sampling rate is 1 second, and its unit is nT. At first, the time-series magnetic field h(t) is transferred to the frequency domain using the Fourier transform. Then, the impedance Z(ω), calculated from the simple 1-D model, is multiplied with H(ω) to determine E(ω). Subsequently, the electric field spectra are transformed back into the time domain to obtain the electric timeseries e(t).
We tested the effectiveness of data-length using 1-hour, 2-hour, 3-hour, 4-hour, 5-hour, 6hour, 12-hour, and 24-hour synthetic MT time-series data. The 1-hour time-series data can only get the maximum reliable period is 150 seconds. The maximum reliable period is determined by the error bar and the impedance result. 4-hour time-series data is relatively good enough to get the maximum reliable period of about 1,000-second. From this result, we confirmed that if there were 2-hour noise-free time-series data, we could get reliable results at least with a period of 300-second. If there are approximately 4-hour noise-free time-series data, we could get reliable data with a period of 1,000-second by RMHHT.   The observation period is shown in Table 1. In this research, we used the data observed on June 26, 2016. The field data consists of natural sources and local noise. The natural MT source coming from the magnetosphere or global lighting is far enough from the local site following the plane wave's assumption. All other parts of the measured electric and magnetic fields are considered as noise.

The Characteristic of the Field MT Data Observed by Phoneix Equipment
We can subdivide the noise and signals into three parts: the correlated noise ( CN ), where the noise is coherent between E and H field, uncorrelated noise ( UN ), and MT signal ( MT ). With this subdivision, we can rewrite the fields E and H as follows: Szarka (1988) and Junge (1996) summarized the active and passive noise sources observed in MT measurements. Most of the artificial noise is caused by earthing currents originated from the electrical equipment. Usually, the remote reference magnetic field has a high correlation with the local magnetic field. Therefore, the correlation between the local magnetic field and remote reference magnetic field is a good indicator of whether the noise is strong or not. show the daytime data, and the right figures show the night data. Each data segment is for ten minutes. The unit of the magnetic field is nT. Fig. 7 shows the Hx component time-series measured simultaneously 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 nighttime at site L6274 and L6276. The variation of the time-series at other sites is almost keeping the same scale. Fig. 8 shows the Hx component, whose sampling rate is 15 Hz, measured simultaneously in different sites. The left figures show the data in the daytime, and the right figures show the data at nighttime. Each segment is ten minutes. The unit of the magnetic field is nT. The variation at site L6274 and L6276 is about 100 times the variation at other sites in the daytime and similar at nighttime. That means the noise to signal ratio is 100 in the daytime; this will bias the impedance tensor a lot. Fig. 9 shows the correlation variation between the local and remote magnetic field (Hx) at L6274 and Y0625 from 3:00:00 to 22:00:00 in UCT time on June 26. The correlation is based on the Pearson correlation using the time-series segment, and each segment length is 10 minutes.
The blue denotes a negative correlation, and the red indicates a positive correlation. Due to local noise, the correlation becomes low or negative in the daytime, and the correlation becomes higher at nighttime. It means that the data in the local night is much quiet. Moreover, we also introduce another parameter, polarization directions, to estimate the background noise in a different frequency band. The polarization directions of the electric field ( ) and magnetic field ( ) (Fowler et al., 1967) at a specific frequency are defined as: The polarization directions describe the time-harmonic variation of two orthogonal field components with constant phasing. 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).
In this research, we analyzed the polarization directions at different frequencies for all sites.
As the signal strength at the dead band (1-10s) is very low, it is easy to influence the local noise.
The typical polarization directions for the dead band is shown in Fig.10. This figure shows the result of polarization directions calculated by the site L7158 at 6-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.
This data set is contaminated by typical industrial noise from a manufacturing plant from the analysis above. The correlation of time-series between the local site and remote site, at site L6274 and L6276, is much higher at nighttime. The magnetic field's polarization has a preferred direction in the dead-band at all sites in the daytime. That means the data is contaminated by strong noise in the dead band at all sites in the daytime. Garcia et al. (2002) also found that the high-frequency band's MT signal was stronger at nighttime. The signal to noise ratio is higher at nighttime than in the daytime at all frequency bands in this data set. Therefore, we consider using the data at nighttime only to calculate the impedance tensor. Fig. 10 The polarization direction for a period of 6-second at site L7158 from 3:00:00 to 22:00:00 in UCT time. The upper figure shows the electric field's polarization directions, and the lower is the polarization direction of the magnetic field.

The Strategy to Calculate the Wideband Time-series Data Based on RMHHT
The processing of broadband MT time-series data by RMHHT is demonstrated in this section.
The data set observed using the MTU-A geophysical instruments was used. The time-series data are stored in three files. The high-frequency band (2,400 Hz) and the middle-frequency band (150 Hz) were sampled intermittently. The length of each segment is 8 seconds and 1 second separately. The low-frequency band (15 Hz) was sampled continuously. Table 2 shows the strategy for the discontinuous time-series. For the windowed FFT, we should transform each segment into the frequency domain and group the target frequency data in the frequency domain. In contrast, the HHT can express the spectrum variation well in the timefrequency domain, as shown in subsection 3.1. We can combine the intermittent time-series data directly and process the discontinuous data in the same way as the continuous time-series data.
Considering the accuracy of the complex coefficient data, we set the effective frequency to 4 times of sampling rate. Using the high-frequency (2,400 Hz) data, calculate the period between 4/2,400 to 4/150 seconds. One continuous segment (1 second Table 3 shows the strategy for continuous time-series. As the computation of HHT is timeconsuming. We design a strategy that can decrease the computation time much. The 15-Hz data was used to calculate the period data between 4/15 to 5 seconds. There are 54,000 data points in 1-hour, and 100,000 points (about 2 hours) are used to calculate the impedance, and we downsampled the continuous 15-Hz data to 1-Hz to compute the period larger than 5 seconds. As shown in the simulation, 4-hour noise-free time-series data can get a relatively reliable result until 1,000-second. We can calculate the full scall period from high frequency to 1,000 seconds by four hours by combining all the results.
This strategy allows us to select data in the time domain. There may be several hours at nighttime that most of the electric equipment is shutdown. The time-series data in these periods are much quiet. This strategy is suitable to use only the quiet period time-series data to calculate the impedance. The BIRRP code calculated the red curves. We only used the BIRRP code to process continuous time-series. The time-bandwidth 3 for the Slepian sequence taper is applied to each section. The result was also calculated by the nighttime data from 16:00:00 to 20:00:00 in UTC on June 26.
The SSMT-2000 calculated the blue curves. The results were calculated from one-day data.
In the SSMT-2000 result, the XY component phase was close to 0º, and the apparent resistivity act as a straight line in the log scale in the long period. That is the phenomenon of artificial noise (Zonge and Hughes,1987 Similar results were obtained at the L6178 and L7158. Fig.14 and Fig.15 show the result calculated using the data at site L6178 and L7158, respectively. The RMHHT result performs well at the full-scale. In the dead band, as we used the relatively quiet data at the Y0625 as a remote site, it performs as well as BIRRP, but in the long period of around 1,000-second, the error bar calculated by RMHHT becomes lagger. Suppose there is another quiet nighttime data; we can combine them directly. In this way, we can get more data at the impedance estimation step in the frequency domain, and the result will become more stable. Fig. 14 The magnetotelluric sounding curves calculated using the data at site L6178. SSMT-2000 was used to calculate the blue curves using one-day data, RMHHT was used to calculate the black curves using the nighttime data from 16:00:00 to 20:00:00 in UTC on June 26 at the remote site Y0625. BIRRP calculates the red curves using the nighttime data from 16:00:00 to 20:00:00 in UTC on June 26. Fig. 15 The magnetotelluric sounding curves calculated using the data at site L7158. SSMT-2000 was used to calculate the blue curves using one-day data, RMHHT was used to calculate the black curves using the nighttime data from 16:00:00 to 20:00:00 in UTC on June 26 at the remote site Y0625. BIRRP was used to calculate the red curves using the nighttime data from 16:00:00 to 20:00:00 in UTC on June 26. At nighttime, most of the electrical equipment is shutdown. Usually, the data is much quieter than in the daytime. We tried to use the quiet period time-series data only to calculate the impedance tensor. We introduced a robust impedance estimator based on the HHT (RMHHT) and demonstrated a strategy to process the broadband MT time-series data properly. It has been shown that we could get a reliable result until 1,000-second from continuous 4-hour noise-free time-series data in the simulation. This short measurement time makes it possible to do MT exploration in strong electromagnetic noise areas. This paper also demonstrated that a time window of 4-hour length was almost noise-free during the local nighttime at site L7158, L6178, and Y0625. We used the data from this quiet period alone to estimate impedances and get a reliable result.
Moreover, if there are several nighttime time-series data, we can combine them directly. In this way, we can increase the number of data sets at the impedance estimation step. The result will become more stable.
This paper showed that RMHHT could get a similar result with the conventional result based on windowed FFT from the quiet data set. We have not shown the method's superiority based on the HHT than the conventional method based on FFT. Any method needs the normal data at the last step of impedance estimation. When transforming the time-series data into the frequency domain, a part of the time-series data, which is contaminated by the transient noise, will influence the whole spectrum calculated by the Fourier transform. In the worst case, every segment contains transient noises. The method based on the FFT is impossible to get a reasonable impedance.
In contrast, HHT can express the spectrum variation in the time-frequency domain well. The transient noises will affect the signal contaminated time-span, and the other parts are not affected.
We can remove the IPs in the contaminated period in the frequency domain before the impedance estimation and get a reliable result. The method based HHT can reject the noisy data maximum in the time-frequency domain. We think it is promising to pay more attention to identifying the noise data in the time-frequency domain based on the HHT.

Availability of data and materials
The magnetic time-series data used to construct the synthetic MT data was downloaded from the INTERMAGNET (International Real-time Magnetic Observatory Network). The data used in this study belong to the Institute of Geophysical and Geochemical Exploration, China Geological Survey under the regional geophysical survey technology research and information system platform construction, DD20190030-05. Alan Chave provided the BIRRP code. Maik Neukirch provided the EMT code.

Competing interests
There are no conflicts of interest associated with this publication.

Funding
The data used in this study belong to the Institute of Geophysical and Geochemical Exploration, China Geological Survey under the regional geophysical survey technology research and information system platform construction, DD20190030-05.