1 Introduction

The magnetotelluric (MT) method is geophysical induction method, using measurements of the naturally time varying geomagnetic and telluric field components at the Earth’s surface, to determine the Earth’s subsurface electrical resistivity (Chave & Jones, 2012). However, both the electric and magnetic fields, which are recorded in an MT measurement, are susceptible to distortion from anthropogenic noise sources, when in close proximity to human infrastructure (e.g., Szarka, 1988). The dominant ‘source’ of anthropogenic distortion is often the power grid utility frequency and associated higher harmonics (50 Hz in Europe, Asia, Africa, Oceania; 60 Hz in North and South America; and both in locals like Japan). Integer multiples of the utility frequency (i.e., 100, 150, 200 Hz and higher for a 50 Hz utility frequency) typically have a disproportionately large effect on the measured signal. The closer such a noise source is to the MT observation location, the stronger the effect of this noise is in the time series. The presence of utility frequency noise propagates through the naturally occurring source field as harmonic oscillations distorting the entire high frequency range in MT studies.

Power line generated noise effects are most noticeably seen in high frequency magnetotelluric response (50/60–10 kHz) preventing accurate recovery of near surface structure. Multiple approaches of varying complexity have been applied to minimize and remove utility frequency noise in MT observations (e.g., Borah et al., 2015; Chen, 2008; Fontes et al., 1988; Trad & Travassos, 2000), the simplest of which is application of a notch filter in frequency domain (Fontes et al., 1988), however, this approach distorts and/or attenuates neighboring frequencies (Butler & Don Russel, 1993). Among the many different notch filters available the finite impulse response notch filter (Borah et al., 2015) applied in the frequency domain has been commonly used in MT data analysis. A more complex approach to filtering the power line utility frequency response is a delay-line filter (Chen, 2008) in the time domain, while more complex numerically, a more accurate response is normally recovered. Wavelet filters are based on combined time–frequency analysis in the ‘wavelet’ domain (Trad & Travassos, 2000), and can be used to isolate the utility frequency of the power grid, but as in the case of other methods, do not deal with time changes of utility frequency.

Filtering approaches such as notch filter, work well, when the utility frequency is stable exactly around one particular frequency, e.g., 50 Hz, otherwise wider frequency band must be removed and more of the natural signal is lost. Power grids are however, inherently unstable with variations in the utility frequency resulting from differences in the demand and supply of electricity. Power grids in high-income countries (e.g., Europe, North America) are usually well balanced with small variations in the utility frequency (e.g., ± 0.5 Hz), and simple removal approaches remove only narrow frequency band (e.g. 49–51 Hz). While, less stable power grids (e.g., Africa, Asia) or renewable energy sources, may produce large utility frequency variations (e.g., ± 2 Hz or more), and simple filtering approaches may result in heavily distorted responses or removal of wider frequency range of natural signal (e.g. 47–53 Hz). Thus, we apply the method of Larsen et al. (2014), who used time varying sinusoidal subtraction for removal of power line utility frequency and its higher harmonics in the time domain and we utilize it for MT signal. Presented approach show possible way of cleaning frequency band around utility frequency in case of time varying noise, which is useful for very noisy datasets, when choose of target frequencies might be difficult. Five test localities within different countries (Table 1) are used to explore the efficacy of sinusoidal subtraction and for comparison with the standard notch filter and remote reference method.

Table 1 Summary of used MT stations
Full size table

2 Utility Frequency of the Power Line Network

Power grids rely on the ability to instantaneously balance the demand (draw by customers) to the supply (contribution from power suppliers), when demand surges beyond supply capacity, the deficit is accommodated via a decrease in the rotational energy of the generators resulting in a decrease in the frequency of supplied power (i.e., a fluctuation and reduction in the utility frequency).

Regulations governing the operation of electricity distribution networks vary by country or by regions linked to interconnected grids. For example, the European standard (EN 50160:2010) requires the utility frequency (50 Hz) to be provided at a stability level of ± 0.5 Hz over the distribution network. Minor variations are permitted in the utility frequency of the transmission system during standard operations (Table 2). However, extreme events can cause larger deviations in the utility frequency, these can range from weather events maintenance and failures of power generation infrastructure; and space weather and geo-magnetic storms disrupting the electrical grid. These extreme events are rare, as such the utility frequency is stable near to 50 Hz the majority of the time. Though, it may naturally drift within a narrow frequency band centered at 50 Hz (Fig. 1).

Table 2 Power line utility frequency standards according to statute EN 50160:2010, specifying frequency stability requirements for European transmission systems
Full size table
Fig. 1
figure 1

Five minutes record of the utility power line frequency within the European grid (10.3.2020 15:08–15:13 UTC) illustrating the utility frequency variations are predominantly restricted to the narrow band 49.95–50.05 Hz. Data accessed from www.mainsfrequency.com on 10.3.2020

Full size image

Fluctuations of the utility frequency under a strict rules typical for European power grids (49.5–50.5 Hz) affects narrower frequency band, than in developing countries, where looser regulations are applied to grid networks. For example, in Vietnam 49–50.5 Hz (Duong et al., 2018), in Bangladesh 48.9–51.2 Hz, with spikes to 48.7 Hz or to 51.5 Hz (Chawdhury & Chattopadhyay, 2018); in Zimbabwe 48.5–51.25 Hz (Zimbabwe grid code, 2017), and Kenya 45–52 Hz (Kenya electricity grid code, 2008). Major fluctuations of grid frequency can be caused by renewable energy sources, which may not be as stable as more established power generations methods. Wind generators in South Africa only need to maintain a utility frequency band of 47–52 Hz (Chown et al., 2018). On the Indonesian island of Sulawesi varying wind speeds result in a variable generation loads supplied to the grid and hence significant short-lived decreases (~ 48.2 Hz) of the utility frequency can occur (Arief et al., 2020). In close proximity to renewable power sources (e.g., wind farms) or in regions where the power grid stability is low, a wider frequency range must be usually removed from MT signal by simpler filtration techniques.

3 Magnetotelluric Method and Cultural Noise

The Earth’s naturally occurring time varying electromagnetic field serve as a source field in the magnetotelluric method. Origin of these variations can be largely attributed to global lighting activity for higher frequencies (i.e., > 1 Hz), while lower frequency (i.e., < 1 Hz) fluctuations result from the interaction of solar winds with the Earth’s magnetosphere (Chave & Jones, 2012). The resulting electromagnetic field propagate down to the Earth’s surface, with partial transmission into the Earth. The transmitted magnetic and electric fields within the Earth produce secondary fields, which can be used for inferring the subsurface electrical structure (Chave & Jones, 2012), because their amplitude and orientation depends directly on resistivity. Surface measurements of the orthogonal horizontal components of the electric (E), and magnetic (H) fields are used to determine the complex impedance tensor (Z) (Cagniard, 1953; Tikhonov, 1950)

$$\mathbf{E}=\mathbf{Z}\mathbf{H}$$
(1)

and the horizontal and vertical components of the magnetic (H) field are used to determine the magnetic field transfer function (K) (Parkinson, 1962)

$${\mathrm{H}}_{\mathrm{z}}=-\mathbf{K}\mathbf{H}$$
(2)

Time series observations of the different components E and H observed at the surface are translated to the frequency domain via Fourier transform and the stacking of calculated power spectra, which are then in turn used to compute the Z and K transfer functions. The impedance tensor can be presented in terms of apparent resistivity (ρz) and phase (ϕ) (Jones et al., 1989). The broad frequency range (10,000–0.0001 Hz) of natural variations make MT a valuable tool for investigating the solid earth over multiple scales of interest from first tens of meters to hundreds of kilometers. Thus, MT is commonly used for geothermal resource evaluation (e.g., Lautze et al., 2020), ore deposit exploration (e.g., Lahti et al., 2019), near surface and crustal investigations of magmatic systems (e.g., Hill et al., 2020; Hogg et al., 2018), hydrogeology (Xu et al., 2020), basin research (e.g., Campanyà et al., 2019), active tectonic processes (e.g., Wannamaker et al., 2017) or investigations of the asthenospheric mantle (e.g., Robertson et al., 2016).

Naturally occurring electromagnetic fields are weak in comparison to those generated by anthropogenic activity, as such accurate MT measurements close to industrialized and even residential areas are challenging. Commonly human induced electromagnetic fields are generated by imbalance in the three-phase transmission network, transmission power grid networks, heavy industrial operations (e.g., mining camps), electric railways (especially DC), and/or electric fences (Szarka, 1988). There has been variable success in recovering the natural signal from the distorted time series using filtering approaches and stacking strategies during processing (e.g., Chave & Thompson, 2004; Gamble et al., 1979; Kappler, 2012; Krings, 2007; Kütter, 2015; Weckmann et al., 2005), when the distortion is comparable in amplitude to the naturally occurring signal. However, there is still inability to effectively separate the small natural part from a large amplitude anthropogenic component of the observed MT response even with complex multi-site acquisition and processing strategies (Egbert, 2002; Oettinger et al., 2001; Wannamaker et al., 2009).

The 50/60 Hz utility frequency of power movement within power grids is the most pervasive cultural noise source encountered in electromagnetic geophysical studies (Chave & Jones, 2012). AC power transmission networks are typically three-phase systems, where individual conducting lines carry voltage separated by one-third of a cycle. In an ideally balanced three phase system, the sum of the phased currents is zero, and the noise produced in respect to electromagnetic geophysical measurements is small, resulting from induction current within the conductors. However, ideally balanced loads are not often the actuality of transmission systems, where the current sum is near to, but not identically zero. These unbalanced non-zero current sum transmission systems produce ‘noise’ for electromagnetic geophysical surveying via current flow on the neutral line and current seepage from the transmission system into the Earth (Chave & Jones, 2012). The amplitude of the utility frequency current is dependent upon the time varying loading of the transmission grid and requirement to meet varying power requirements throughout the day, effectively causing a modulation and broadening of spectral peaks. The dynamic nature of transmission networks may also produce slight variations in the utility frequency (Szarka, 1988) making it difficult to isolate in geophysical measurements. It is not unusual for large transmission networks to produce in excess of tens higher harmonics of the utility frequency (Pellerin et al., 2004) affecting frequencies up to 1 kHz and above, reaching the dead band (range between 1 and 5 kHz, where is MT signal naturally attenuated) for shallow looking approaches like Audio-magnetotellurics (AMT) (Hoover et al., 1978).

Typically, even with the high sampling rates of modern MT instruments (512 Hz and greater) MT time series obtained in developed regions (where power grids are present) are dominated by anthropogenic sinusoidal 50/60 Hz signal to the extent of masking all of the target naturally occurring signal and even other anthropogenically induced noise components. On Fig. 2 can be seen observed time series and stacked power spectra by robust technique in Mapros (Friedrichs, 2004) of a representative MT sounding from Mýtina maar, Czechia. It shows both the dominance of the utility frequency over the natural sourced signal and how the dominance of the utility frequency propagates through the signal spectrum with large peak amplitudes occurring at the higher harmonics in the stacked power spectra.

Fig. 2
figure 2

A Example of observed time series in length of 0.125 s (recorded with a 4096 Hz sampling rate) affected by the power grid utility frequency from Mýtina maar, Czechia. The naturally occurring time varying natural electromagnetic field response is dominated by the large amplitude 50 Hz sinusoidal response of the power grid utility frequency. B Power spectra response (stacked in 4096 sample windows) of the example time series showing the increased energy associated with the 50 Hz utility frequency and higher harmonics

Full size image

4 Power line Frequency Treatment

Following the approach of Larsen et al. (2014) who used sinusoid subtraction (Butler & Don Russel, 1993) to suppress power line noise in time-series from magnetic resonance soundings, we apply this approach to MT soundings. The power line harmonic components occur at number of higher harmonic frequencies mf0, where m is positive integer and f0 utility frequency, while the power line interference can be modelled as a sum of sinusoidal signals with frequencies given by mf0, each having an independent amplitude and phase. When the utility frequency, amplitude and phase have been determined, the powerline noise can be constructed and subtracted from the recorded data.

The utility frequency constantly fluctuates around its base frequency, so the sinusoidal subtraction has to be applied in individual time windows. These time windows must contain sufficient data for processing MT time series, while, being short enough to isolate changes in the utility frequency (for discussion see chapter 6.1). The power line signal for a time window k can be described as

$${S}_{p}\left(k\right)={\sum }_{m}{A}_{m}^{p}\mathit{cos}\left(2\pi m\frac{{f}_{0}}{{f}_{s}}k+{\Phi }_{m}^{p}\right),$$
(3)

where f0 is the power line utility frequency, fs the sampling frequency, and \({A}_{m}^{p}\) & \({\Phi }_{m}^{p}\) are the amplitude and phase of the mth harmonic component (first m is power line utility frequency). Summation over m extends to all excited harmonics in the time window of investigation.

In solving Eq. (3), model parameters Am, Φm and f0 must be determined first through non-linear optimization, and can be simplified by separating the cosine

$${A}_{m}\mathit{cos}\left(2\pi m\frac{{f}_{0}}{{f}_{s}}k+{\Phi }_{m}\right)={\alpha }_{m}\mathit{cos}\left(2\pi m\frac{{f}_{0}}{{f}_{s}}k\right)+{\beta }_{m}\mathit{sin}\left(2\pi m\frac{{f}_{0}}{{f}_{s}}k\right),$$
(4)

where αm and βm are coefficients used for estimation of amplitude and phase of mth harmonic given by:

$${A}_{m}=\sqrt{{\alpha }_{m}^{2}+{\beta }_{m}^{2}},$$
(5)
$$\mathit{tan}\left({\Phi }_{m}\right)=-\frac{{\beta }_{m}}{{\alpha }_{m}}.$$
(6)

Coefficients αm and βm are chosen such that a sinusoidal frequency mf0 fits the harmonic oscillations in MT signal (within time window k) as closely as possible in a least-squares sense (Butler & Don Russel, 1993).

The only remaining non-linear parameter is f0 after simplification in Eq. (4). As f0 is close to the utility power line frequency, the simplest method for determining the model parameters is to treat f0 as a parameter and solve the linear equation system. Employing a standard least-squares approach with a large number of initial f0 values surrounding the utility frequency, the optimum (true) f0 can be determined. Once the (true) f0 is acquired, the modelled power line signal is subtracted from the associated individual time window k. The measured electric and magnetic components of the natural electromagnetic field inherently have different responses to power line noise (e.g., due to orientation, and/or sensor types), so the (true) f0 must be determined individually and each component unique response subtracted from the recorded data within each time window.

5 Removing Power line Interference from MT data

With the high sampling rates typical in MT data collection (> 512 Hz) the resulting time series have millions of samples (individual data). To improve the computational efficiency, when working with large data sets such as MT, the optimization search for f0 may be separated into two steps. The first step involves a roughly sampled sweep using a 0.025 Hz increment between 49.5 and 50.5 Hz (in case of stable European power grid, otherwise signal power spectra need to be inspected for adequate frequency range) to identify the approximate minima f0a. Then a finer sweep increment of 0.001 Hz is used over the band f0a ± 0.02 Hz to reach the absolute minima of f0. The result of the two step sweep and identification of f0 for the example MT station from Czechia are shown in Fig. 3 identifying both the approximate f0a and absolute minima f0.

Fig. 3
figure 3

Identification of the optimum utility frequency f0 through the multi stage minimization residual power sweep of the MT time series window of interest from site in Czechia. A First stage coarse sweep located the ‘approximate’ f0a in range 49.95 ± 0.02 Hz. B Second stage fine sweep centered around the ‘approximate’ f0a subsequently identify the absolute f0

Full size image

The efficiency of this approach for noise removal is dependent on the length of the considered time window, resulting in an inherent trade-off between computation time and accuracy. The longer the time window used the less computational resources required, but the larger the compromise (i.e., variation of the estimate over the entirety of the time window) is in optimization of the f0 minima introduced by temporal variations of the utility frequency.

Representative examples of the efficacy of the sinusoidal subtraction removal approach and standard notch filter from two test locations in Czechia and Mongolia are shown in Fig. 4. Same as all other following data presented in this study, original and filtered time series data using both notch and sinusoidal removal of f0 were processed by Mapros (Friedrichs, 2004) using robust stacking algorithm (Egbert & Booker, 1986). Different noise distortion levels are observed at the two locations depending on proximity to infrastructure and the regional population density. A significant improvement in the final transfer functions are clearly visible in Czechia (Fig. 4A), while in Mongolia (Fig. 4B) an improved response is recovered but not to the same extent probably due to close proximity to residential houses.

Fig. 4
figure 4

Representative examples of noise removal MT response curves by both sinusoidal subtraction and notch filter: A Mýtina maar (Czechia) and B Tsenkher (Mongolia). Improved results are recovered by both methods for both locations, though, site from Mongolia has a lesser improvement due to higher level of cultural noise. The ideal time window length (in seconds) for sinusoidal subtraction is noted for both sites

Full size image

Sinusoidal subtraction was also tested on MT site from Norway, which was recorded with remote reference station. On Fig. 5 are compared different possible processing schemes: robust stacking, sinusoidal subtraction filtering, notch filtering and processing by remote reference method. The best improvement was achieved by sinusoidal subtraction and then by notch filter, while remote reference station recovered only slightly better sounding curves.

Fig. 5
figure 5

Noise removal MT response curves by sinusoidal subtraction, notch filter and remote reference (RR) station from Norway site. Improved results are achieved especially by sinusoidal subtraction and in lesser extent by notch filter. Remote reference method improved sounding curves the least

Full size image

For testing the ability to deal with time changes of utility frequency, the sinusoidal subtraction was tested on data from Egypt and Ethiopia. Example soundings from both countries contained only weak 50 Hz peaks in the power spectra, without any higher harmonics and data processed by robust methods were free of scattered resistivity/phase points (Egypt station Fig. 6A; Ethiopia not shown) probably due to data collection in sparsely populated areas. However, sinusoidal subtraction found frequencies f0 within 1 s time windows in wide range of 48–52 Hz for Egypt and 49.5–50.5 Hz for Ethiopia during entire recordings, suggesting that a notch filter needs to remove wider frequency range in areas with more unstable grid networks.

Fig. 6
figure 6

Results of synthetic test, where was added sinusoidal noise with wide span 48–52 Hz together with nine odd higher harmonic components into original clean data from Egypt presented by sounding curves and sample stacked spectra of Ex component (Y axis has same range for A–F). A original data, B data with added noise, C–E data filtered by sinusoidal subtraction with different window lengths, F data filtered by notch filter. A, B, F were processed by 2 s long windows, while C–E by indicated window lengths. Note: Y axis of stacked spectra is in same scale for all four plots

Full size image

For further validation of sinusoidal subtraction for power grids with a wide utility frequency range (i.e., peaks), we conducted a synthetic test. In almost ‘clean’ data from Egypt with a sampling frequency of 4096 Hz (Fig. 6A), we added sinusoidal noise representing the wide utility frequency range with span ± 2 Hz centered around 50 Hz (and nine odd higher harmonic components with same span ± 2 Hz centered around 150, 250⋯950 Hz). Each harmonic component had a different amplitude and phase, while the base frequency between 48–52 Hz, amplitude and phase changed every 2 s (8192 sample long windows). Whereas original data are almost free of noise as shown in power spectra (Fig. 6A), the added noise significantly affect both spectra and processed data (Fig. 6B). Sinusoidal subtraction applied over a 1 s and 2 s time windows successfully removed the distortion introduced by the sinusoidal noise and recovered original sounding curves (Fig. 6C, D). Sinusoidal subtraction applied over a 4 s window was not able to deal with faster time changes of utility frequency and failed (Fig. 6E). The standard notch filter designed for removing frequency range of 48–52 Hz removes only a narrow band of the distorted utility frequency in given range (Fig. 6F). For complete removal of observed peaks a wider frequency range have to be employed by notch filter due to width of frequency peaks close to limit values 48 and 52 Hz.

6 Discussion

6.1 Effect of Chosen Window Length on Sinusoidal Subtraction

The application of the stacked sinusoidal noise removal approach shows improved final transfer function responses; however, the trade-off or ‘cost’ of these improved responses is the numeric computational cost, which can require large computational times (ca. 1 h for 30 min recording by 4096 Hz sampling rate). The approach allows for a very narrow time varying frequency band representing the power grid utility frequency to be identified and removed, maintaining significantly more of the source signal than the application of the standard static notch filters. The ‘width’ of the sinusoidal subtraction filter applied here was 10–4 Hz, though this can be varied based on local conditions.

Limitations of the presented approach are largely related to the window length selection, over which f0 is swept and identified. The selected window length is then fixed through subsequent processing and power spectra stacking. While selection of the fixed window lengths for the filtering is a necessity, this can be undesirable for power spectra stacking (Borah et al., 2015). However, there is not an elegant workaround for the fixed window length requirement of the filtering – for further processing where having the fixed window length is detrimental, the only solution is to perform the utility frequency filtering multiple times using different window lengths.

The power grid utility frequency can be considered ‘constant’ or ‘static’, when the time window is small (i.e., 20 ms). However, such short time window is not suitable for MT time series analysis. The requirement for longer time windows in the processes of moving from the time to frequency domain of the MT response, results in the optimum f0 utility frequency dependency on the length of the time window (Fig. 7) within the first coarse frequency sweep (49.5–50.5 Hz). Shorter time windows are both more precise and accurate in identification of the power grid utility frequency f0 (Fig. 7A-B). When the length of the time window is increased (> 0.5 s) the utility frequency experiences smaller drift, as such, the absolute accuracy of the f0 identification decreases and represents the average peak energy frequency for the selected time window (Fig. 7C–F).

Fig. 7
figure 7

Determination of the peak energy frequency f0 using 15 min of Ex observations from Czechia site, illustrating the dependence of time window length selection in the f0 identification sweeps. Time window length was varied from 0.125 to 4.0 s. Shorter time windows A–B allow more precise identification of utility frequency, while identification in long windows C–F more likely represents the average peak energy frequency

Full size image

The choice of window length for the sinusoidal filtering can be seen as a trade-off between the accuracy of utility frequency identification & removal, and longer window lengths, which produce better results in MT time series processing (Borah et al., 2015). Shorter window lengths identify the utility frequency for removal more accurately, but do not provide enough data required for estimation of the MT transfer functions. Improved results are obtained in all cases for different time window lengths between 0.25 and 4.0 s (Fig. 8A-E, original data from site in Czechia are in Fig. 4A). In general, longer time windows produce more stable estimates of the MT response for time windows up to ~ 4 s (with a sampling rate of 4096 Hz), while effectively removing the higher harmonic effects of the utility frequency successfully. For longer time windows (> 4 s) the drift in the utility frequency outweighs the benefit of the longer time window for estimation of the MT response. For the 8 s time window the sounding curves become scattered due to variations of the utility frequency within the time window. Analyzing the results from Fig. 8, the optimal length of time window for the Czechia data is between 2 and 4 s.

Fig. 8
figure 8

Post-sinusoidal sum utility frequency f0 noise removal MT curves for the different time window lengths used to identify f0 at Czechia site shown in Fig. 4A. The best sounding curves are produced by sinusoidal subtraction within window lengths 2–4 s (D–E)

Full size image

The stability of the power grid utility frequency is important in identifying the optimum time window length to sweep for f0. More stable utility frequency permit to use longer time windows, which offer better performance during MT transfer function estimates. However, test proved, that even short windows as 1 s are sufficiently long for MT processing and lead to improved MT sounding curves. Determination of the optimal time window to be used is numerically expensive given the trial and error approach used in its determination.

6.2 Comparison Between Sinusoidal Subtraction and Notch Filter

Differences between sinusoidal subtraction and notch filter are generally small, however, sinusoidal subtraction removes from original signal just one exact frequency (e.g. 49.596 Hz) compared to frequency band by notch filter (e.g. 49–51 Hz). For deeper comparison of both approaches, stacked power spectra from the site in Czechia were carefully tested (Fig. 9). Almost no changes can be observed within time window of 0.25 s between both methods (Fig. 9A). However, as the time window is increased to 2 s the benefit of sinusoidal subtraction becomes clear; the notch filter produces energy drops around the notch location and sinusoidal subtraction preserved more of the original signal (Fig. 9B, see magnified parts of spectra). Inaccurate removal of 50 Hz peak in case of Ey channel might be caused by non-ideal time window length of 2 s for this channel or by some additional noise source present in time series. With application of 8 s long window sinusoidal subtraction is failing due to already too large variations of utility frequency and notch filter becomes again dominant (Fig. 9C).

Fig. 9
figure 9

Stacked power spectra for original (top line), notch filter (middle line) and sinusoidal subtraction (bottom line) for three different time window lengths of Czechia site. Power spectra have same scale, but are vertically shifted for better visibility. In case of 2 s long window sinusoidal subtraction gives generally better results as it can be seen in magnified parts of 50 Hz peak and preserves more of the original signal. On contrary with application of 8 s window, the drift in the utility frequency outweighs benefits of subtraction filtering and notch filter gives better result

Full size image

The benefit of sinusoidal subtraction method can be seen with the synthetic noise added example with a wide range of utility frequency variations (Fig. 6). Here sinusoidal subtraction recover original response in case of applied time windows shorter than time changes of utility frequency (Fig. 6C, D vs F) better than standard notch filter. In case of longer time window of 4 s applied over synthetic utility frequency changing every 2 s, sinusoidal subtraction failed (Fig. 6E). However, synthetic signal was designed in a simple way, that utility frequency could change in given span (48–52 Hz) without any limitations. In real world example, utility frequency will probably drift over time (as in Fig. 1 and just faster) from e.g., 48 Hz to 51 Hz and not change instantly in one or two seconds. Thus, sinusoidal subtraction will be probably able to deal with such changes. It should be noted that real world examples will likely not contain such a well-defined sinusoidal noise component and the results in such cases would be somewhat reduced. Synthetic example proved, that sweeping multiple window lengths to optimize the process of identifying f0 is important especially in locations with wide utility frequency peaks.

7 Conclusion

The sinusoidal subtraction method for removal of the time varying power line utility frequency signal and higher harmonics in MT data significantly improves observational data in urban areas. The adaptive approach allows for removal of a single frequency compared to standard notch filtering, resulting in the useable response more closely honoring the observed natural signal. Analysis of power spectra from test soundings within Europe show that sinusoidal subtraction preserves more of the original signal than a standard notch filter. Presented approach should be beneficial especially in areas, where large variations of utility frequency occur and standard notch filters need to remove wider frequency band from original data. This can cause issues in case of very noisy datasets, when choose of target frequencies might be difficult or in shallow MT studies, where structures of interest lies in depths corresponding to 50 Hz frequency or its higher harmonic components. A trade-off between optimum time series processing window length and identification precision of the utility frequency must be managed when using sinusoidal subtraction. This is currently achieved by trialing multiple time window lengths to optimize the trade-off due to the instability of the utility frequency, which can be numerically expensive.