Background noise estimation of the geomagnetic signal

A fast Fourier transform was applied to fit the geomagnetic diurnal variation. Fitting results showed that when the polynomial degree was greater than 160, the residual error was close to 0 nT. White noise is the main component of the residual error when the polynomial degree was greater than 160, so this method was adopted to calculate the background noise of the geomagnetic field. Spectrum analysis further demonstrated that the noise estimation result is reliable.

1 Introduction

The geomagnetic field is critical for solar activity monitoring, space weather detection, and some investigations into crustal motions that lead to geomagnetic field changes (Chapman and Bartels, 1940; Campbell, 1997). Because of the various factors influencing the geomagnetic field, the sources of geomagnetic signal noise are also diverse. According to previous studies, the noise of the geomagnetic signal can be divided into two types (Ren, 2006). Variable noise changes over time and comes from fluctuations within conductive fault zones or from instability within the observation environment, which acts as an antenna to couple with the external geomagnetic field. Background noise is more stable; it originates within a stable observation environment and from instrumental responses, such as thermal noise and other electronic noise. The noise within geomagnetic data plays an important role in evaluating the quality of geomagnetic data and also has an impact on scientific research (Yao et al., 1995).

Historically, variable noise in geomagnetic observations was calculated by first difference. In previous studies, some results such as temporal characteristics, spatial distribution features, and influencing factors of variable noise were achieved (Yan et al., 2013; Wang et al., 2015). According to previous research results, the intensity of variable noise is extremely weak and usually lower than an instrumental resolution of 0.1 nT. This suggests that background noise may be the main component of geomagnetic noise. However, previous research also showed that it is difficult to calculate background noise quantitatively because of the inseparability of the geomagnetic signal and noise. In seismology, signal-to-noise ratio (SNR) estimation on noisy data is mainly obtained through the energy superposition method, spectrum analysis, or power spectrum calculation (Zhang et al., 2009). Because of the nonstationarity of the geomagnetic field, these methods are not suitable for geomagnetic SNR estimation. Zhu et al. (2012, 2013) and Wang et al. (2015) applied principal component analysis (PCA) to suppress noise in airborne electromagnetic data. Nevertheless, PCA needs observation data from at least three groups at the same observatory, which is improbable for most geomagnetic observatories. Jiang et al. (2013) used maximum likelihood estimation to calculate geomagnetic noise through multiple iterations. However, many previous studies have suffered from deficiencies in testing.

In this paper, the diurnal geomagnetic data were fitted through a fast Fourier transform (FFT), and then the residual error between original and fitted data was obtained to estimate the background noise of geomagnetic data. In previous studies, researchers applied an FFT to geomagnetic diurnal variation (Han et al., 2009; Zhao et al., 2014; Koch and Kuvshinov, 2015; Yamazaki and Maute, 2017). However, almost all of them focused on the Sq diurnal variation, so the polynomial degree is no greater than 6. In general, when geomagnetic disturbances are absent, the first four harmonics are sufficient to capture most of the variability in a daily record of the geomagnetic field. We suggest analysis of geomagnetic diurnal variation by the use of FFT with a degree of 250; the residual error may represent changes in background noise. The testing result showed that this approach to estimate background noise in geomagnetic data is effective.

2 Data processing

FFT is the most widely used method of spectrum analysis. Any periodic signal can be decomposed into several components such as first harmonics (T), second harmonics (T∕2), third harmonics (T∕3) and more (Cooley and Tukey, 1965) through the FFT.

A time-series signal can be expressed as a function of sine and cosine as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }f\left(l\right)=\sum _{m=\mathrm{0}}^{M-\mathrm{1}}\left[a\right(m)\cos m{\mathit{\lambda }}_l+b(m\left)\sin m{\mathit{\lambda }}_l\right],\text{(2)}& {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_l={\displaystyle \frac{\mathrm{2}\mathit{\pi }}{M}}l,\end{array}$$

where a(m) and b(m) are coefficients of sine and cosine functions, respectively. λ_l is a function of l, and l indicates the sequence number of the data series. M represents the total number of data point.

Normal daily variation in the geomagnetic field mainly comprises the first six harmonic components (Fig. 1); these components represent signals of period 24, 12, 8, 6, 4.8, and 4 h, and results of higher degree (N) achieve closer fits.

Figure 1 Fitting results of first six harmonic components.

Figure 2 Error results of 10th–250th harmonic components.

The background noise of the geomagnetic vertical component (Z) is more complex because this component is more susceptible to the change in observation environment. Furthermore, in order to reduce the influence of the external geomagnetic field, data from the quietest days of 2013 were chosen and an FFT was applied to fit them for the 10th–250th harmonic components. The residual error (Error in Eq. 3) between original data and fitted data was calculated by as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }\mathrm{Error}={\displaystyle \frac{\sum _{l=\mathrm{1}}^M\left(ff\left(l\right)-f\left(l\right)\right)}{M}}(l=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots },\mathrm{86}\mathrm{400}).\end{array}$$

Here, ff(l) represents the original geomagnetic data, and f(l) indicates the fitted data through FFT. The sampling rate of the original data is 1 Hz, so the total number (M) of samples in 1 day is 86 400. As displayed in Fig. 2, the residual error is less than 1.0 nT when the polynomial degree is greater than 10, and smaller residual errors generally correlated with larger polynomial degrees. When the polynomial degree is greater than 160, the residual error approaches 0 nT. Based on previous analysis results, it is reasonable to assume that the fitted data of the 160th degree could represent the original signal in which background noise is not contained.

Figure 3 shows the original signal and the FFT-fitted data with a polynomial degree of 160 on 29 May 2013 at the LYH (37.40 N, 114.70 E) observatory as an example; a constant is added between them for comparison. The fitted curve is almost identical to the original curve, though the fitted curve is smoother. This implies that the background noise of the geomagnetic signal could be estimated through FFT fitting with a polynomial degree of 160. The background noise of the geomagnetic vertical component is marked as Z_noise and is obtained from Eq. (4). ff(l) in Eq. (4) represents the original geomagnetic data, and f(l) indicates the fitted data through FFT with a polynomial degree of 160. Figure 4 shows the estimated background noise on 29 May 2013 at the LYH observatory as an example. It is randomly distributed between −0.2 and 0.2 nT with a mean value of 0 nT.

$$\begin{array}{}\text{(4)}& {\displaystyle }Z\mathit{\_}\mathrm{noise}=ff\left(l\right)-f\left(l\right)\end{array}$$

White noise is a random signal with a mean value of 0. Based on the characteristics of Z_noise obtained from Fig. 4, we think the main composition of Z_noise may be white noise, which associated with the geomagnetic instrument and the observation environment.

Figure 3 Original signal and fitting data of the Z component.

Figure 4 Estimated noise of the Z component.

3 Testing process

Standard white noise is a random signal with a mean value of 0, and its autocorrelation function is close to 0 when the lag (τ) is not equal to 0. To confirm that the main composition of Z_noise is white noise, the autocorrelation function of Z_noise is calculated from Eq. (5).

$$\begin{array}{}\text{(5)}& {\displaystyle }{R}_{\mathit{\tau }}={\displaystyle \frac{E\left[(X_t-\mathit{\mu })(X_{t+\mathit{\tau }}-\mathit{\mu })\right]}{{\mathit{\sigma }}^{\mathrm{2}}}}\end{array}$$

R_τ is the autocorrelation function of signal, E is the expected value operator, X_t indicates data of Z_noise at time t, μ and σ² represent the mean value and variance of Z_noise, and τ is the lag. Figure 5 shows the autocorrelation function of background noise on 29 May 2013 at the LYH observatory. The autocorrelation function of Z_noise clearly reaches up to 1 when τ=0 and is close to 0 when τ≠0, the same as the autocorrelation of white noise. Therefore, it is demonstrated that the main composition of Z_noise is white noise, and the background noise of the Z component of the geomagnetic field (Z_noise) could be obtained through FFT fitting with a polynomial degree of 160.

Figure 5 Autocorrelation function of background noise in the Z component on 29 May 2013 at LYH.

The SNR of the geomagnetic signal can be calculated as Eq. (6) when background noise is obtained. The result shows that the SNR of geomagnetic data from the LYH observatory is about 47.

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathrm{SNR}}_{\mathrm{0}}=\mathrm{10}\cdot {\log }_{\mathrm{10}}\left({\displaystyle \frac{\sum f^{\mathrm{2}}\left(l\right)}{\sum {\left(ff\left(l\right)-f\left(l\right)\right)}^{\mathrm{2}}}}\right)\end{array}$$

4 Spectrum analysis

To contrast the original geomagnetic data and the noise-free geomagnetic data, their frequency spectra were analyzed. Waveforms of the geomagnetic vertical component (Z) during each interval of 30 min were subjected to FFT spectrum analysis as follows:

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }X\left(k\right)=\sum _{k=\mathrm{0}}^{N-\mathrm{1}}x\left(k\right){\mathit{e}}^{-\mathrm{2}\mathit{\pi }nki/N}k=\mathrm{0},\mathrm{1},\mathrm{\dots },N-\mathrm{1},\text{(8)}& {\displaystyle }& {\displaystyle }\left|F(u,\mathit{\nu })\right|={\left[R^{\mathrm{2}}\left(u,\mathit{\nu }\right)+I^{\mathrm{2}}(u,\mathit{\nu })\right]}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

Here, $i=\sqrt{-\mathrm{1}}$, and R(u, ν) and I(u, ν) indicate the real and imaginary parts of FFT result, respectively. F(u, ν) represents the amplitude of the FFT spectrum.

Figure 6 Spectrum of the original geomagnetic signal and the filtered data (Z component).

The upper panel of Fig. 6 shows amplitude spectrum of original signal on 29 and 30 May 2013 at the LYH observatory. The amplitude spectrum of the original data is clearly more irregular; its background noise can be seen as scattered points in the spectrum, and the intensity of these points is not related to frequency or time. Because of the influence of background noise, the amplitude spectrum of the original geomagnetic signal does not show any significant changes over this interval. The bottom panel of Fig. 6 represents the amplitude spectrum of the geomagnetic data after removing background noise. It is different from the upper panel. Irregularly scattered points representing background noises are removed almost completely, and, as a result, the amplitude spectrum is more regular and informative. From the bottom panel of Fig. 6, some significant characteristics are obtained. First, geomagnetic energy changes with frequency: the lower the frequency of the signal is, the higher the geomagnetic energy. Second, geomagnetic activity in each 30 min block is different. In the first 24 half hours of 29 May 2013 and the first 24 half hours of 30 May 2013 the geomagnetic activity is more intense than any others. In contrast, these characteristics are not immediately obvious in the amplitude spectrum of the original data, implying that FFT filtering is an effective method for geomagnetic background noise estimation.

5 Discussion and conclusion

The background noise of the geomagnetic signal can be obtained through FFT filtering with a polynomial degree of 160 on the vertical component (Z) during the quietest days. The main conclusions are summarized as follows:

1. The residual error between FFT-filtered data and the original signal approaches 0 nT as the polynomial degree is greater than 160, and it has been confirmed that a residual error with a degree of 160 could represent the background noise of the geomagnetic signal.

2. Geomagnetic background noise is a random signal distributed between −0.2 and 0.2 nT.

3. The autocorrelation function of background noise showed that it reaches up to 1.0 when lag (τ) is 0 and is close to 0 in other cases, which confirms that white noise is the main component of geomagnetic background noise.

4. Spectrum analysis further confirms that FFT filtering is an effective method for geomagnetic background noise estimation, and some geomagnetic changes are more remarkable after filtering.

FFT-filtered data with a polynomial degree greater than 160 could represent the original geomagnetic signal with a period of less than 540 s. Any signal with a period of less than 540 s in the original data will be removed completely in the filtered data. To avoid overprocessing, data from the quietest days were chosen, avoiding short-period variations such as pulsations or geomagnetic bays. In addition, because the vertical component (Z) contains more noise information than other components and is not as susceptible to the external geomagnetic field, it was chosen as the analysis object in this paper. Because the main factors influencing background noise are the local observation environment and the instrumental response, the geomagnetic background noise at different observatories differs. For any one particular observatory, background noise is usually nearly invariable due to the stability of the observation environment and the instrument condition.