A reduction method for magnetic disturbance caused by DC railway system

The stray current of direct current (DC) railway systems causes magnetic disturbance in geomagnetic measurements, which may complicate the identification of useful information. The magnetic disturbance exhibits broadband characteristics in the frequency domain. In this paper, we propose a noise reduction method based on the adaptive Kalman filter to extract useful signals from the geomagnetic data with a high level of noise. The covariance matrixes of both the process noise (Q) and measurement noise (R) can be adaptively estimated to improve the performance of the adaptive Kalman filter. The proposed method is adopted to process the geomagnetic data collected at the Beijing Geomagnetic Observatory (BJI), which is affected by the DC railway system. The magnetic disturbance is largely reduced, and the signal-to-noise ratios of the horizontal and vertical components of the geomagnetic field are improved by more than 14 dB and 27 dB, respectively. The K-indices are calculated to evaluate the performance of the adaptive Kalman filter method. To assess the influence of the adaptive Kalman filter on natural signals, the geomagnetic data that contain rapid variations are processed. The denoising results show that the adaptive Kalman filter can effectively reduce the magnetic disturbance caused by DC railway system without large impact on the natural geomagnetic rapid variations.


Introduction
Geomagnetic data have many potential scientific uses, which can be used, for example, to study the geomagnetic variations and space weather conditions (e.g., Finlay et al. 2020;Love and Finn 2017), monitor volcanic and seismic activities (e.g., Masci and Di Persio 2012;Takla et al. 2014), study ocean tidal signals (e.g., Maus and Kuvshinov 2004;Schnepf et al. 2018) and interactions between the lithosphere and atmosphere(e.g., Kherani et al. 2012). In recent years, with economic growth and metropolitan developments, direct current (DC) railway systems have become popular in many cities due to their faster speed and greater travel capacity for relieving traffic pressure (Chen et al. 2017b;Ibrahem 2017), which will cause some negative impacts on geomagnetic measurements (Chen et al. 2017a). In general, the DC traction loop consists of a third rail (or overhead copper wire), a running rail and a power station. When the train operates, the DC power station supplies traction current, and the third rail carries the traction current from the power station to the train, the current then returns to the power station through the running rail. Since the running rail cannot be completely insulated over the ground, a part of the traction current, which is called the stray current, will flow into the earth through the running rail and return to the power station by an earth diversion (Wang et al. 2018). Hence the magnetic disturbance is caused by the imbalance of current in the third rail and the running rail (Lowes 2009).
Several studies have addressed magnetic disturbances caused by DC railway systems. Dupouy (1950) was probably the first researcher to study the magnetic disturbance of the stray current. Since then, the magnetic disturbance caused by DC railway systems has been Open Access *Correspondence: yuguo@ouc.edu.cn 1 College of Marine Geoscience, Ocean University of China, Qingdao 266100, China Full list of author information is available at the end of the article Ding et al. Earth, Planets and Space (2021) 73:107 simulated and analyzed (Lowes 2009;Pirjola 2011;Pirjola et al. 2007;Tokumoto and Tsunomura 1984). Linington (1974) reported the influences of disturbance caused by DC railway systems on archaeological magnetic surveys. Sun and Wang (1984) surveyed the magnetic disturbance caused by electric railways in Beijing. Lowes (1987) observed the vertical component of magnetic fields close to the railway and found that the individual train movements can be resolved by analyzing the magnetic data. Chen et al. (2017a) observed magnetic disturbance caused by the mass rapid transit system in Taiwan and analyzed the disturbance in the frequency domain.
Most previous studies have focused on the simulation of magnetic disturbance of DC railway systems. To the best of our knowledge, no specific study has focused on the characteristic analysis and reduction method of magnetic disturbance caused by DC railway systems. In this paper, we analyze the characteristics of the magnetic disturbance of DC railway systems in both the frequency domain and time domain. The adaptive Kalman filter is adopted to reduce the magnetic disturbance, and its performance is evaluated by calculating the K-indices and processing the geomagnetic data that contain natural rapid variations. The rest of this paper is organized as follows. In "Data and analysis" section, we investigate the influence of magnetic disturbance caused by DC railway system on geomagnetic measurements. The magnetic data are analyzed in both the frequency domain and time domain. In "Reduction of magnetic disturbance using the adaptive Kalman filter" section, we briefly describe the adaptive Kalman filter method and apply it to the geomagnetic data collected on both the geomagnetically quiet and geomagnetically disturbed days. To assess the influence of the adaptive Kalman filter on the natural signals, the geomagnetic data that contain rapid variations were also processed. Finally, in "Conclusions" section we draw the conclusions of this paper.

Data and analysis
The magnetic disturbance caused by the DC railway system was observed near the Beijing and Qingdao metro lines. In this section, we will analyze the characteristics of magnetic disturbance using geomagnetic data collected in Beijing and Qingdao.
The geomagnetic data recorded at the Beijing Geomagnetic Observatory (BJI), which is located 3 km away from one of the Beijing metro lines, are affected by the DC railway system. If a geomagnetic observatory is far from metro lines, the geomagnetic data will not be affected by the DC railway system, such as the geomagnetic data collected at the Beijing Ming Tombs Geomagnetic Observatory (BMT). Figure 1 shows the geomagnetic data collected at BJI and BMT on May 3, 2019, which are located 33.4 km away from each other. The data sampling rate is 1/60 Hz. From Fig. 1, one can see that the geomagnetic data at the two observatories have the same variation trend. It is obvious that the geomagnetic data collected at BJI are affected by the DC railway system, which operates between 04:30 am and 00:30 am the next day (local time) (Zhang et al. 2018), while the magnetic disturbance is very weak from 00:30 am to 04:30 am, when the metro is not in operation. The vertical component of the geomagnetic field is much noisier than the horizontal component, and the vertical and horizontal components of magnetic disturbance are about 6 nT and 3 nT, respectively. The signal-to-noise ratios (SNR) of the vertical and horizontal components of the geomagnetic field are − 29.57 dB and − 14.51 dB, respectively. The SNR is defined as follows: where P s and P n represent the power of the noise-free geomagnetic signal and noise, respectively. Here, the geomagnetic data collected at BMT, which is not affected by the DC railway system, are regarded as noise-free signals.
To analyze the characteristics of magnetic disturbance in the frequency domain, we use the Welch's method (MATLAB, Version R2019b) to estimate the power spectra density (PSD). Figure 2 displays the PSD of the geomagnetic data from both BJI and BMT. From Fig. 2, one can see that the disturbance caused by the DC railway system shows broadband characteristics in the frequency domain, which means that it will affect broadband geomagnetic signals. It is obvious that the frequency ranges of the horizontal and vertical components of the magnetic disturbance are 0.0006 Hz ~ 0.08 Hz and 0.0001 Hz ~ 0.08 Hz, respectively. The vertical component of the magnetic disturbance has higher energy than the horizontal one. Compared with the geomagnetic data collected at BMT, the PSD of the vertical component of the geomagnetic data collected at BJI is increased by 30.99 dB at a frequency of 0.08 Hz, while the PSD of the horizontal component is increased by 18.86 dB.
The magnetic disturbance caused by the DC railway system has also been observed near one of Qingdao metro lines. Figure 3 shows the geomagnetic data collected in Qingdao for a total of 3.5 h. A magnetometer was installed on our campus and is 680 m away from the metro railway. The period of this experiment covers three states of metro operation. The trains run on two parallel (1) SNR = 10log 10 P s P n tracks in opposite directions during period (a), while the train runs only on one of tracks in one direction during period (b). In period (c), the power is turned off for maintenance. From Fig. 3, one can see that the magnetic disturbance varies over time and is obviously related to the metro operation states. When the trains operate in both directions, the magnetic field varies between 51,700 nT and 52,270 nT, and the maximum magnetic disturbance is about 200 nT, while the magnetic field varies between 52,000 nT and 52,100 nT, and the maximum magnetic disturbance is about 100 nT when the train operates only in one direction. The magnetic disturbance almost disappears when the railway is powered off. The SNR of the data collected in the two-track and one-track operation states are − 37.46 dB and − 26.12 dB, respectively. This indicates that the magnetic disturbance is larger in the two-track operation state, in which a greater number of trains are running.

Reduction of magnetic disturbance using the adaptive Kalman filter
The magnetic disturbance caused by the DC railway system in the geomagnetic data will complicate the identification of useful information. We employ the adaptive Kalman filter to reduce the magnetic disturbance. To evaluate the performance of the proposed method, we processed not only the geomagnetic data of both the geomagnetically quiet and geomagnetically disturbed days, but also the data containing storm sudden commencement (SSC), sudden impulse (SI), and magnetic effect of a solar flare (Sfe).

Adaptive Kalman filter
The Kalman filter is a kind of linear optimal estimation method that is regarded as one of the most famous Bayesian filters (Li et al. 2015). It was first proposed by Kalman (1960), and has been utilized in diverse areas of science and engineering, such as noise reduction in magnetotelluric, seismic and gravity data (e.g., Leśniak et al. 2009;Deng et al. 2016;Wang et al. 2019), target tracking and navigation (e.g., Yang et al. 2010), and digital image processing (e.g., Piovoso and Laplante 2003). The status equation and observation equation of the Kalman filter at time k are expressed as follows: (2) X(k) = �X(k − 1) + W (k) process noise and observation noise vector, respectively, assuming that W and V are zero mean Gaussian white noise and independent of each other.
The process of the Kalman filter consists of two steps: Predict and Update (Akhlaghi et al. 2017;Li et al. 2015). They will be briefly described in the following.
In the prediction step, an estimate of the state at the current time step is predicted using it at the previous step (Akhlaghi et al. 2017;Li et al. 2015): in which X(k|k − 1) and P(k|k − 1) are the predicted state estimate and its covariance matrix, respectively.Q(k) is the covariance matrix of process noise. X(a|b) indicates the estimate of X at time a given observations up to and including at time b ≤ a.
In the update step, the current prediction is combined with current observation information to refine the state estimate (Akhlaghi et al. 2017;Li et al. 2015):  where d(k) , S(k),R(k) , K (k) , X(k|k) and P(k|k) are measurement innovation, innovation covariance, covariance matrix of observation noise, Kalman gain, posteriori state estimate and posteriori covariance matrix, respectively.
The covariance matrix of process noise Q and observation noise R need to be determined in the filter process, which have a significant impact on the performance of the filter (Akhlaghi et al. 2017;Ding et al. 2007). Improper choice of Q and R may influence the performance of the filter method. Many estimation methods have been proposed to determine Q and R in the filter process (e.g., Akhlaghi et al. 2017;Basso et al. 2017;Ding et al. 2007;Liebich et al. 2017;Saho and Masugi 2015). With the use of adaptive filtering approach, we estimate Q and R , and it can be expressed as (Akhlaghi et al. 2017)  where ε(k) = Y (k) − HX(k|k) is the difference between the observed value and estimated one. α is a forgetting factor between 0 and 1, and a larger α puts more weights on previous estimates and longer time delays to catch up with changes (Akhlaghi et al. 2017). α is set to 0.7 in the examples presented in this paper. (11)

Data of both geomagnetically quiet and geomagnetically disturbed days
The adaptive Kalman filter is adopted to treat the geomagnetic data collected at the Beijing Geomagnetic Observatory (BJI) on May 5 (geomagnetically quiet day) and May 14 (geomagnetically disturbed day), 2019. The processed horizontal and vertical components of the geomagnetic field are shown in Fig. 4. For comparison, the raw data is also plotted. From Fig. 4, one can see that the magnetic disturbance caused by DC railway system has been effectively reduced and the SNR has been improved. The SNR of the horizontal components of geomagnetic data on May 5 and May 14 are improved by 19.58 dB and 14.8 dB, respectively, and the SNR of the vertical  Note that the adaptive Kalman filter may not completely predict rapidly changing geomagnetic data such as those marked by blue rectangle in Fig. 4.
The geomagnetic K-index is a quasi-logarithmic local index of geomagnetic activities, which quantifies disturbances in the horizontal component of geomagnetic data with an integer in the range 0-9 with 1 being calm and 5 or more indicating a geomagnetic storm (Uwamahoro and Habarulema 2014). The K-indices are derived from the maximum fluctuations of the horizontal component of geomagnetic data in a 3-h interval. To evaluate the performance of the adaptive Kalman filter, we adopted the Finnish Meteorological Institute (FMI) method (V1.0) to calculate the K-indices, as shown in Fig. 5. For comparison, the K-indices for the geomagnetic data at BMT were also plotted. From Fig. 5a, one can see that on the geomagnetically quiet day (May 5, 2019), the K-indices for the raw data at BJI (black bars) except for the last two 3-h intervals are larger than those for the data at BMT (blue bars), while after the application of the adaptive Kalman filter, the K-indices (red bars) are reduced and become equal to those for the data at BMT. Actually, the maximum fluctuations in the last two 3-h intervals are also reduced, as shown in Fig. 6a. This may indicate that the adaptive Kalman filter is useful to improve the reliability of the calculated K-indices.
From Fig. 6b, one can see that on the geomagnetically disturbed day (May 14, 2019), the maximum fluctuations for the raw data at BJI are larger than that for the data at BMT, while after the application of the adaptive Kalman filter, the maximum fluctuations are reduced and become close to those for the data at BMT. Note that the K-indices for both the raw data and the processed data have no much difference (Fig. 5b). This is because the magnitudes of the natural magnetic variations are larger than those of the magnetic disturbances caused by DC railway system.

Data with naturally rapid magnetic variation
The above example demonstrates that the adaptive Kalman filter method can effectively reduce the magnetic disturbance caused by DC railway system. In the following example, we will assess the influence of the adaptive Kalman filter on natural rapid variations, including storm sudden commencement (SSC), sudden impulse (SI) and magnetic effect of a solar flare (Sfe) (Curto 2020).
We randomly selected 6 days of geomagnetic data collected at BJI from January 2016 to May 2019, which contain rapid variations (SSC, SI and Sfe), and processed them using the adaptive Kalman filter. The processed Fig. 8 Difference between the raw data and processed data by the adaptive Kalman filter at BJI (a). Difference between BMT and BJI before and after processing for horizontal component (b) and vertical component (c) horizontal component is shown in Fig. 7. For comparison, the raw data is also presented. From Fig. 7, one can see that the magnetic disturbances caused by DC railway system are reduced and the natural rapid variations are retained. This means that the adaptive Kalman filter will not affect the recognition of rapid variations.

Influence of the adaptive Kalman filter on natural signals
The adaptive Kalman filter has been adopted to reduce the noise in the geomagnetic data collected at BJI, as shown in Fig. 1 (black line). After the application of the adaptive Kalman filter, the SNR of the vertical and horizontal components of the geomagnetic data are improved by 27.5 dB and 21.35 dB, respectively. To assess the influence of the adaptive Kalman filter on natural signals, we calculate the difference between the raw data and the filtered data at BJI (from 01:00 to 04:00 am, when the metro was powered off ), as shown in Fig. 8a. The differences between them lie between − 0.6 nT and 1.4 nT. We also calculated the difference between BMT and BJI before and after processing, as shown in Fig. 8b, c. The high-frequency magnetic disturbance has been reduced. However, the natural signals are also affected by the filter. The impact of the filter on the horizontal component is larger than that on the vertical component.

Conclusions
The geomagnetic data collected at Beijing and Qingdao shows that DC railway system can cause magnetic disturbance, which will complicate the identification of useful information. The vertical component of the magnetic disturbance is larger than the horizontal component, and the magnetic disturbance is related to the metro operation state. The magnetic disturbance collected at "twotrack operation state" is larger than that collected at "one-track operation state", and the magnetic disturbance almost disappears when the railway is powered off.
To improve the signal-to-noise ratio, we conducted studies of noise reduction method based on adaptive Kalman filter. The results of practical data processing show that the adaptive Kalman filter is effective in reducing the magnetic disturbance caused by DC railway system. After adaptive Kalman filtering, the magnetic disturbance due to DC railway system can be largely reduced without affecting the recognition of the natural magnetic variations, including storm sudden commencement (SSC), sudden impulse (SI) and magnetic effect of a solar flare (Sfe). The adaptive Kalman filter has some impact on the natural signals, especially the geomagnetic data collected on geomagnetically disturbed days.