Spatial gradients of geomagnetic temporal variations causing the instability of inter-station transfer functions

Abstract

Spatial gradients in the primary geomagnetic fields directly contribute to both the amplitudes and phases of inter-station transfer functions (IS-TFs). This suggests that, for the analysis of subsurface resistivity structures, IS-TFs should be carefully treated by checking the establishment of the plane-wave assumption. Geomagnetic time-series data include various and complicated characteristics and accordingly, time–frequency domain analysis is suitable for the discussion of spatial gradients of time-varying geomagnetic fields. However, such evaluations are complicated by the huge amount of information contained in the spectrograms from several stations. Therefore, we propose a Multi-Channel Nonnegative Matrix Factorization (MC-NMF) method that can decompose raw spectrograms into several components, allowing the spatial gradient of each geomagnetic temporal variation to be identified. We confirm that such components actually affect the estimation of IS-TFs using data acquired at the Kakioka and Memambetsu magnetic observatories in Japan. We derive the year-to-year changes in IS-TFs from each set of paired stations among Kakioka, Kanoya, and Memambetsu observatories. Although the IS-TFs should exhibit opposite polarities (a negative correlation) when the input and output observatories are swapped; surprisingly, some of them have “identical” polarities. The application of MC-NMF shows that the analyzed geomagnetic data include several components that have various spatial gradients. Although IS-TFs sometimes fail to give the expected implication regarding the spatial gradients of geomagnetic temporal variations, MC-NMF can verify whether the IS-TFs exhibit any spatial gradients. Thus, the use of IS-TFs with MC-NMF can yield better implications regarding subsurface resistivity information.

Introduction

The electromagnetic (EM) responses of Earth, such as magnetotellurics (MT), are often used to evaluate the subsurface resistivity structure. The inter-station transfer function (IS-TF) between horizontal geomagnetic data at two sites is one type of EM response, and has been used to analyze the subsurface structure (Egbert and Booker 1989; Soyer and Brasse 2001; Arora and Rao 2002; Campanyà et al. 2019). As shown by Soyer and Brasse (2001) and Campanyà et al. (2019), the advantage of IS-TFs is that they improve the consistency of the inverted model because the data from different sites are directly related. IS-TFs are also used as indicators of the temporal changes of resistivity structures triggered, for example, by large earthquakes (Honkura and Koyama 1978; Honkura 1979; Beamish 1982; Hattori 2004) or crustal uplifts (Honkura and Taira 1983). However, IS-TFs have not been used as often as the MT responses to analyze subsurface structures.

Frequently, the changes in MT responses and IS-TFs are not caused by the subsurface environment, but by the source field (electric current in the ionosphere or magnetosphere) of natural geomagnetic fluctuations (Egbert et al. 2000; Brändlein et al. 2012; Romano et al. 2014; Vargas and Ritter 2016; Murphy and Egbert 2018; Sato 2020). In particular, IS-TFs may be expected to be strongly affected by the source field because the spatial gradients of the primary geomagnetic fields directly affect both the amplitudes and phases of IS-TFs. This suggests that IS-TFs should be carefully treated by checking the establishment of the plane-wave assumption or spatial homogeneity of geomagnetic temporal variations when using them for the analysis of subsurface resistivity structures. When using both MT responses and IS-TFs that are not biased by the spatial gradients of geomagnetic variations, we can interpret the subsurface resistivity structure in detail, as reported by Soyer and Brasse (2001) and Campanyà et al. (2019).

Geomagnetic time-series data include various complicated information regarding time-varying geomagnetic fields. Time–frequency analysis (i.e., at the stage of spectrograms instead of IS-TFs) would be suitable for evaluating the spatial gradients of geomagnetic variations. However, the evaluation remains difficult because of the huge amount of information contained in the raw spectrograms from several stations. Therefore, we have developed a Multi-Channel Nonnegative Matrix Factorization (MC-NMF) method. This is a new method that decomposes the raw spectrograms of horizontal geomagnetic data into several matrices/components, each with their own spatial gradients.

Egbert (1997) developed a method of estimating IS-TFs using principle component analysis and robust statistics. Later, Egbert (2002) discussed the spatial gradients of geomagnetic variations by applying his method (Egbert 1997) to array data. Raw geomagnetic spectrograms include more information than the IS-TFs derived from the raw data. MC-NMF can extract components with more information than the available data channel, although the maximum number of components extracted by principle component analysis is limited to the number of available channels. Thus, this study uses MC-NMF to evaluate the spatial gradients of geomagnetic temporal variations.

This paper first introduces the general method of calculating IS-TFs and the details of MC-NMF. We then verify that the IS-TFs between two magnetic observatories in Japan (Kakioka and Memambetsu) are shifting temporally. Evaluations using MC-NMF indicate that the presence of large spatial gradients causes such shifting. We show that the year-to-year changes in IS-TFs derived from array data cannot yield the expected implications regarding spatial characteristics such as the gradients of geomagnetic temporal variations, and we elucidate the causes using MC-NMF. Hereafter, we use the term “geomagnetic event” to distinguish each geomagnetic temporal variation; events are not related to real geomagnetic phenomena (e.g., pc1). In particular, we define an “anomalous (geomagnetic) event” as having a characteristic of spatial gradient different from that of the other events.

Inter-station transfer functions

The x and y directions of geomagnetic field data recorded at one site (e.g., site 1) are defined as geographically north and east, respectively, as in the general coordinate systems used in geomagnetism. In the time–frequency domain after transformation using a Short-time Fourier Transform (STFT), which is a sequence of Fourier Transforms of windowed/tapered time-series data, these data have the following relationship with the magnetic field data at another site (e.g., site 2):

$$\left( {\begin{array}{*{20}c} {H_{{{\text{site}}1,x}} \left( {f,t} \right)} \\ {H_{{{\text{site}}1,y}} \left( {f,t} \right)} \\ \end{array} } \right) = \left( {\begin{array}{*{20}c} {T_{xx} \left( f \right)} & {T_{xy} \left( f \right)} \\ {T_{yx} \left( f \right)} & {T_{yy} \left( f \right)} \\ \end{array} } \right)\left( {\begin{array}{*{20}c} {H_{{{\text{site}}2,x}} \left( {f,t} \right)} \\ {H_{{{\text{site}}2,y}} \left( {f,t} \right)} \\ \end{array} } \right),$$

(1)

where $H_{{{\text{site}}1,x}}$ and $H_{{{\text{site}}1,y}}$ are the magnetic fields in the x and y directions at site 1, respectively, and are defined as the output data. $H_{{{\text{site}}2,x}}$ and $H_{{{\text{site}}2,y}}$ are the magnetic fields at site 2 and are defined as the input data. Thus, $T$ in Eq. 1 is the IS-TF of site 2 to site 1. Generally, to derive $T_{xx}$, the least-squares method is applied (Vozoff 1972; Schmucker 1984; Neska 2006):

$$T_{xx} = \frac{{\left\langle {H_{{{\text{site}}1,x}} H_{{{\text{site}}2,x}}^{*} } \right\rangle \left\langle {H_{{{\text{site}}2,y}} H_{{{\text{site}}2,y}}^{*} } \right\rangle - \left\langle {H_{{{\text{site}}1,x}} H_{{{\text{site}}2,y}}^{*} } \right\rangle \left\langle {H_{{{\text{site}}2,y}} H_{{{\text{site}}2,x}}^{*} } \right\rangle }}{{\left\langle {H_{{{\text{site}}2,x}} H_{{{\text{site}}2,x}}^{*} } \right\rangle \left\langle {H_{{{\text{site}}2,y}} H_{{{\text{site}}2,y}}^{*} } \right\rangle - \left\langle {H_{{{\text{site}}2,x}} H_{{{\text{site}}2,y}}^{*} } \right\rangle \left\langle {H_{{{\text{site}}2,y}} H_{{{\text{site}}2,x}}^{*} } \right\rangle }},$$

(2)

where $H^{*}$ denotes the complex conjugate of $H$ and $\left\langle {H_{{{\text{site}}1,x}} H_{{{\text{site}}2,x}}^{*} } \right\rangle$ denotes the averaged value of the cross-spectra. This may be biased by noise at site 2. However, the data available from several geomagnetic observatories provide high-quality (i.e., approximately noise-free) geomagnetic data, and we can estimate the IS-TF accurately from the relevant data. When applying a remote reference method (Gamble et al. 1979), the conjugate spectra in Eq. 2 are replaced by reference data. Although a remote reference method can produce high-accuracy IS-TFs, the determination of which site affects the IS-TFs and the causal relationship may be ambiguous. Thus, we use the standard least-squares method and focus on only two sites, thus preventing any bias from the spatial gradients of geomagnetic temporal variations included in the reference data.

Multi-Channel Nonnegative Matrix Factorization

MC-NMF is a natural expansion of Complex Nonnegative Matrix Factorization (C-NMF: Kameoka et al. 2009; King and Atlas 2011; Kitamura 2019), which decomposes an observed complex spectrogram $\varvec{X}$ ($F \times T$ matrix), obtained using an STFT, with $F$ frequencies and $T$ time windows into “Basis vectors” $\varvec{B}$ ($F \times K$ matrix), “Activations” $\varvec{U}$ ($K \times T$ matrix), and phase terms $\varvec{\varphi }$ ($K$ matrices with $F \times T$ elements). A singular value decomposition provides an orthogonal basis set. However, the magnitude spectra are not always orthogonal. Therefore, C-NMF is used in this study because this method is not limited to orthogonal basis sets, and thus offers more degrees of freedom. Using C-NMF, $\varvec{X}$ is represented as

$$X\left( {f,t} \right) \cong \mathop \sum \limits_{k} B\left( {f,k} \right)U\left( {k,t} \right){\text{e}}^{{j\varphi \left( {k,f,t} \right)}} ,$$

(3)

where $X\left( {f,t} \right)$, $B\left( {f,k} \right)$, $U\left( {k,t} \right)$, and $\varphi \left( {k,f,t} \right)$ are elements of the above matrices, respectively, all elements of $\varvec{B}$ and $\varvec{U}$ are real and nonnegative, and $j$ is an imaginary unit. These “real” and “nonnegative” constraints suggest that the Basis vectors $\varvec{B}$ denote the patterns of magnitude spectra included in the observed spectrogram $\varvec{X}$, and the Activations $\varvec{U}$ denote the temporal changes in the spectral amplitudes.

To analyze several spectrograms, in MC-NMF, we expand $\varvec{X}$, $\varvec{B}$, and $\varvec{\varphi }$ to $\varvec{X}_{m,d}$, $\varvec{B}_{m,d}$, and $\varvec{\varphi }_{m,d}$, where $m \left( {m = {\text{site}}1, \ldots ,{\text{site}}M} \right)$ and $d \left( {d = x,y} \right)$ denote the sites and directions, respectively. The MC-NMF model is illustrated in Fig. 1. Based on the maximum a posteriori estimation, one of Bayesian estimation, the posterior probabilities of $\varvec{B}_{m,d}$, $\varvec{U}$, and $\varvec{\varphi }_{m,d}$ must be maximized to estimate them from the observed spectrograms $\varvec{X}_{m,d}$. Although such posterior probabilities cannot be derived directly, we can write them as

$$p\left( {\varvec{B}_{m,d} ,\varvec{U},\varvec{ \varphi }_{m,d} |\varvec{X}_{m,d} } \right) \propto p\left( {\varvec{X}_{m,d} |\varvec{B}_{m,d} ,\varvec{U},\varvec{ \varphi }_{m,d} } \right)p\left( {\varvec{B}_{m,d} } \right)p\left( \varvec{U} \right)p\left( {\varvec{\varphi }_{m,d} } \right).$$

(4)

Fig. 1

Model of MC-NMF. $F$, $T$, and $K$ denote number of frequencies, time windows, and Basis vectors, respectively. $\varvec{X}$ denotes observed spectrograms, $\varvec{B}$ denotes Basis vectors, and $\varvec{U}$ denotes Activations

$p\left( \varvec{U} \right)$ is generally assumed to be a super-Gaussian (Kameoka et al. 2009; King and Atlas 2011), and written as

$$p\left( \varvec{U} \right) = \mathop \prod \limits_{k,t} \frac{1}{{2\varGamma \left( {1 + \frac{1}{q}} \right)\xi }}\exp \left( { - \frac{{\left| {U\left( {k,t} \right)} \right|^{q} }}{{\xi^{q} }}} \right),$$

(5)

where $\varGamma$ is Gamma function, $0 < q < 2$, and $\xi$ is the scale parameter. $p\left( {\varvec{X} |\varvec{B},\varvec{U},\varvec{ \varphi }} \right)$ denotes the reconstructed errors between the left- and right-hand sides in Eq. 3 and is assumed to follow a Gaussian distribution, as in preceding studies (Kameoka et al. 2009; King and Atlas 2011; Kitamura 2019). Kameoka et al. (2009) and King and Atlas (2011) demonstrated the stability of C-NMF under the assumption that $p\left( \varvec{B} \right)$ and $p\left( \varvec{\varphi } \right)$ follow uniform distributions; we set $p\left( {\varvec{B}_{m,d} } \right)$ and $p\left( {\varvec{\varphi }_{m,d} } \right)$ accordingly.

Consider the analysis of geomagnetic spectrograms. As reported by Vörös et al. (1998), the empirical probability density function of geomagnetic fluctuations has a long tail, which is not easily modeled by a Gaussian distribution. Thus, an assumption of the specific distributions regarding Activations does not conflict with the physical phenomena, because a super-Gaussian distribution has a long tail. In addition, we may use uniform distributions for the Basis vectors and phases because these distributions are suitable for cases with no prior information (Gelman et al. 2013), as noted in preceding studies (Kameoka et al. 2009; King and Atlas 2011; Kitamura 2019). As shown by Makishima et al. (2019), MC-NMF can estimate the Basis vectors just as well with no prior information as with explicitly available prior information. Thus, our assumption of no prior information and uniform distributions is reasonable.

Focusing on the exponential term of $p\left( {\varvec{B}_{m,d} ,\varvec{U},\varvec{ \varphi }_{m,d} |\varvec{X}_{m,d} } \right)$, the objective function to be minimized can be written as

$$J\left( {B,U,\varphi } \right) = \mathop \sum \limits_{m,d,f,t} \left| {X_{m,d} \left( {f,t} \right) - \mathop \sum \limits_{k} B_{m,d} \left( {f,k} \right)U\left( {k,t} \right){\text{e}}^{{j\varphi_{m,d} \left( {k,f,t} \right)}} } \right|^{2} + 2\lambda \mathop \sum \limits_{k,t} \left| {U\left( {k,t} \right)} \right|^{q} ,$$

(6)

where $\lambda$ is a weighting coefficient, which includes information on $\xi$ and $\varGamma$ in Eq. 5. The nonnegative constraint indicates that the nonlinear objective function $J\left( {B,U,\varphi } \right)$ is minimized using the majorize minimize (MM) algorithm (Hunter and Lange 2004), because the algorithm is guaranteed to converge (Kameoka et al. 2009; Kitamura 2019), as shown in Appendix 1. $\varvec{B}_{m,d}$ and $\varvec{U}$ have a scale ambiguity and should be constrained according to $\mathop \sum \limits_{t} \left| {U\left( {k,t} \right)} \right|^{2} = 1$. We define a separated component:

$$Y_{m,d} \left( {k,f,t} \right) = B_{m,d} \left( {f,k} \right)U\left( {k,t} \right){\text{e}}^{{j\varphi_{m,d} \left( {k,f,t} \right)}} ,$$

(7)

which is an element of $\varvec{Y}_{m,d} \left( k \right)$ at a time–frequency slot ($f$, $t$).

Modeling $\varvec{U}$ using a super-Gaussian distribution, $\varvec{Y}_{m,d} \left( k \right)$ follows a super-Gaussian and approaches a sparse matrix (see Eq. 23) because $p\left( {\varvec{B}_{m,d} } \right)$ and $p\left( {\varvec{\varphi }_{m,d} } \right)$ are assumed to follow uniform distributions. Such sparse components modeled by a super-Gaussian distribution can be considered independent of each other. For example, Hyvärinen and Oja (2000) proposed an independent component analysis that extracts sparse and independent components included in observed data and models them as super-Gaussian distributions. When some independent components are mixed, the distributions approach Gaussian distributions according to the central limit theorem. They proved that the independent components can be modeled by super-Gaussian distributions because of the inverse process. Thus, each element $Y_{m,d} \left( {k,f,t} \right)$ of the separated components $\varvec{Y}_{m,d} \left( k \right)$ at a time–frequency slot ($f$, $t$) can be treated independently from the others, because they are sparse and modeled by super-Gaussian distributions. In addition, this independence and sparseness ensure that each $Y_{m,d} \left( {k,f,t} \right)$ does not overlap with other elements at a time–frequency slot $\left( {f, t} \right)$. Instead of optimizing the phase term, we can define the phase term as $e^{{j\varphi_{m,d} \left( {k,f,t} \right)}} = \frac{{X_{m,d} \left( {f,t} \right)}}{{\left| {X_{m,d} \left( {f,t} \right)} \right|}} \left( {k = 0, \ldots ,K - 1} \right)$, as in the experiment reported by Kameoka et al. (2009) and Kitamura (2019). The initial models of the Basis vectors $\varvec{B}_{m,d}$ and Activations $\varvec{U}$ are set based on nonnegative independent component analysis (Kitamura and Ono 2016).

The number of Basis vectors $K$ is based on many criteria, such as Bayes’ information criterion (Owen and Perry 2009) and Bayesian nonparametric modeling (Hoffman et al. 2010). We use the criterion based on the convergence of the root mean square error (RMSE) between the observed data and models, given by

$${\text{RMSE}}\left( \% \right) = 100\sqrt {\frac{{\mathop \sum \nolimits_{m,d,f,t} \left| {X_{m,d} \left( {f,t} \right) - \mathop \sum \nolimits_{k} B_{m,d} \left( {f,k} \right)U\left( {k,t} \right){\text{e}}^{{j\varphi_{m,d} \left( {k,f,t} \right)}} } \right|^{2} }}{{\mathop \sum \nolimits_{m,d,f,t} \left| {X_{m,d} \left( {f,t} \right)} \right|^{2} }}} ,$$

(8)

which is the modified version used by Sawada et al. (2013) and Kameoka et al. (2018). Their criterion regarding $K$ states that the cost-function like RMSE has almost converged and the raw spectrograms are represented clearly by the sum of the separated components $\varvec{Y}_{m,d} \left( k \right)$. If $K$ is greater than the exact number for representing the observed spectrograms, the one component originally represented by one Basis vector is separated into two or more components. In such cases, we must identify the physical meanings of the Basis vectors and Activations (Sawada et al. 2013; Kameoka et al. 2018). In this study, we determine the number of Basis vectors such that the RMSE between the observed data and models does not increase when $K$ increases. The physical meaning of each Basis vectors is then determined as follows.

To identify spatial characteristics such as the gradient of the geomagnetic event reflected by each Basis vector and Activation, we focus on the difference between $\frac{{B_{{m_{1} ,d}} \left( {f,k} \right)}}{{\mathop \sum \nolimits_{l} B_{{m_{1} ,d}} \left( {f,l} \right)}}$ and $\frac{{B_{{m_{2} ,d}} \left( {f,k} \right)}}{{\mathop \sum \nolimits_{l} B_{{m_{2} ,d}} \left( {f,l} \right)}} \left( {m_{1} \ne m_{2} } \right)$. When the difference is large, the geomagnetic event reflected by the $k$th Basis vector and Activation will have a characteristic of spatial gradient different from that of other events, and can be considered as an anomalous event. However, the equality $\frac{{B_{{{\text{site}}1,x}} \left( {f,k} \right)}}{{\mathop \sum \nolimits_{l} B_{{{\text{site}}1,x}} \left( {f,l} \right)}} = \ldots = \frac{{B_{{{\text{site}}M,y}} \left( {f,k} \right)}}{{\mathop \sum \nolimits_{l} B_{{{\text{site}}M,y}} \left( {f,l} \right)}}$ is established (or approximately established) if all geomagnetic events have the same spatial gradients. These concepts are explained in more detail in Appendix 2. Hereafter, we define the Basis vector Rate (${\text{BR}}$) as

$${\text{BR}}_{m,d} \left( {f,k} \right) = \frac{{B_{m,d} \left( {f,k} \right)}}{{\mathop \sum \nolimits_{l} B_{m,d} \left( {f,l} \right)}}.$$

(9)

Relationship between IS-TFs and Basis vectors

We now show that IS-TFs shift following anomalous geomagnetic events, and that we can evaluate such events using MC-NMF. We analyze geomagnetic data acquired at three magnetic observatories in Japan: Kakioka (KAK), Kanoya (KNY) and Memambetsu (MMB). The locations of those observatories are shown in Fig. 2. The sampling rate of geomagnetic data at these three observatories is 1/60 Hz (i.e., one sample per minute). The data used in this study were recorded in January and February of 2011 (59 days). These data can be downloaded from the Kakioka magnetic observatory website (http://www.kakioka-jma.go.jp/, accessed December 26, 2019).

Fig. 2

Location of each (geo)magnetic observatory. (1) Beijing Ming Tombs (BMT). (2) Kakioka (KAK). (3) Kanoya (KNY). (4) Memambetsu (MMB). The original of this map without any annotation can be downloaded from Craft MAP (http://www.craftmap.box-i.net/world.php; Accessed 26 Dec 2019)

A first-order differential filter (a type of high-pass filter) was applied to these geomagnetic data for pre-processing and detrending. The detrended data are hereafter referred to as the raw time-series data. The data were transformed into the time–frequency domain using an STFT. The Fourier Transform length was fixed to 512 samples (i.e., 512 min) and a frequency range from 9/30,720 to 108/30,720 Hz was analyzed (i.e., about 3000–300 s in a period). The spectrograms from the three observatories in the x or y direction are shown in Fig. 3.

Fig. 3

Amplitudes of spectrograms of geomagnetic data with x (upper) and y (lower) direction, acquired at KAK (left), at KNY (middle), and at MMB (right). The vertical axis denotes the frequencies between 9/30,720 Hz and 108/30,720 Hz, and the standard is 9/30,720 Hz (i.e., the value 0 in the vertical axis corresponds to 9/30,720 Hz). The horizontal axis denotes time window number, with a total of 165

We applied MC-NMF to the spectrograms and modeled them using 10 Basis vectors (Fig. 4a), which are the patterns of magnitude spectra included in each observed spectrogram, and 10 Activations (Fig. 4b), which are common factors in each observed spectrogram and denote the temporal changes in the spectral amplitudes. The objective function (Eq. 6) was optimized through 3000 iterations, and the RMSE (Eq. 8) was found to be 3% and convergent. When $K$ was set to 9 or 11, the RMSE was approximately 3%. This seems to uphold the condition that the model in MC-NMF (i.e., the sum of $\varvec{Y}_{m,d} \left( k \right)$) can represent the raw spectrograms exactly. To evaluate the spatial gradient of each geomagnetic event reflected by a Basis vector, we compare ${\text{BR}}_{{m_{1} ,d}} \left( {f,k} \right)$ and ${\text{BR}}_{{m_{2} ,d}} \left( {f,k} \right)$ in Eq. 9, ($m_{1} ,m_{2} =$ KAK, KNY, or MMB, and $m_{1} \ne m_{2}$, $d = x, y$), and summarize the results in Fig. 4c. Based on Fig. 4c, the difference between ${\text{BR}}_{{{\text{KNY}},x}} \left( {f,0} \right)$ and ${\text{BR}}_{{{\text{MMB}},x}} \left( {f,0} \right)$ appears to be large at many frequencies, especially between 30/30,720 and 95/30,720 Hz. To demonstrate how we evaluate the spatial gradient of each geomagnetic event included in the geomagnetic data, we reconstruct spectrograms by multiplying Basis vector 0, Activation 0, and their phases (i.e., $\varvec{Y}_{m,d} \left( 0 \right)$ defined in Eq. 7) from 30/30,720 to 95/30,720 Hz, and then applying an inverse STFT (Fig. 5a). In Fig. 5b, we show the raw time-series data filtered by a band-pass filter between 30/30,720 and 95/30,720 Hz. Note that the time series in Fig. 5 are cumulatively aggregated (i.e., the inverse process of a first-order differential filter) so that their dimensions can be unified and expressed in nanotesla (nT).

Fig. 4

Results obtained by applying MC-NMF to the geomagnetic spectrograms from three observatories. a Basis vectors included in the data with x (upper) and y (lower) directions at KAK (left), at KNY (middle), and at MMB (right). The vertical axis denotes frequencies and the horizontal axis denotes IDs of Basis vectors. b Activations. The vertical axis denotes IDs of Activations (same as Basis vectors) and the horizontal axis denotes time window number. c${\text{BR}}$ with x (upper) and y (lower) directions at KAK (left), at KNY (middle), and at MMB (right). The vertical axis denotes frequencies and the horizontal axis denotes IDs of Basis vectors

Fig. 5

Comparison of geomagnetic time-series data at three observatories. a Time-series data transformed from separated component $Y_{m,d} \left( {0,f,t} \right)$ within the range of 30/30,720–95/30,720 Hz. b Raw time-series data filtered by a band-pass filter with the pass band between 30/30,720 Hz and 95/30,720 Hz

Ordinarily (i.e., except for the annotated peaks of the top three in Fig. 5b), the amplitudes of geomagnetic signals at MMB are larger than those at KAK and KNY. However, focusing on the annotated peaks of the x-direction data around 70,000 s (the top three peaks in Fig. 5b), which correspond to the top three peaks in Fig. 5a, the amplitudes at KAK and KNY are almost the same as and greater than those at MMB, respectively. The amplitudes of the top three peaks in Fig. 5b are − 7 nT (KAK), − 9 nT (KNY), and − 7 nT (MMB). The geomagnetic event corresponding to Basis vector 0 can be considered to have a characteristic of spatial gradient different from the others. From this analysis, we can identify anomalous geomagnetic events in the time–frequency domain using MC-NMF.

We check the relationship between the temporal changes in IS-TFs and the anomalous events evaluated by MC-NMF. The IS-TFs between the geomagnetic data at KAK (output) and at MMB (input) during January/February 2011 are derived using Eq. 2. The IS-TFs of MMB to KAK during January/February from 2000 to 2011 (i.e., 22 months) are calculated to provide standard values. Here, we focus on the Basis vectors reflecting anomalous events, especially those for which (i) the ${\text{BR}}$ is greater than 10% and (ii):

$$\frac{{\left| {{\text{BR}}_{{m_{1} ,d}} \left( {f,k} \right) - \left\{ {{\text{BR}}_{{m_{2} ,d}} \left( {f,k} \right) + {\text{BR}}_{{m_{3} ,d}} \left( {f,k} \right)} \right\}/2} \right|}}{{\left\{ {{\text{BR}}_{{m_{2} ,d}} \left( {f,k} \right) + {\text{BR}}_{{m_{3} ,d}} \left( {f,k} \right)} \right\}/2}} > \varTheta ,$$

(10)

where $m_{1} ,m_{2} ,m_{3} =$ KAK, KNY, or MMB ($m_{1} \ne m_{2}$, $m_{2} \ne m_{3}$, and $m_{3} \ne m_{1}$) and $\varTheta$ is a threshold. Equation 10 represents the standardization distance of ${\text{BR}}$ at site $m_{1}$ from those at other sites. This is a quantitative definition of anomalous events, allowing us to check the relationship between such events and “quantitative” values of IS-TFs. To enable quantitative evaluation, we set the transfer function difference (${\text{TFD}}$) as

$${\text{TFD}}_{d} \left( f \right) = \frac{{\left| {{\text{TF}}_{22M,d} \left( f \right) - {\text{TF}}_{2M,d} \left( f \right)} \right|}}{{\left| {E_{22M,d} \left( f \right) + E_{2M,d} \left( f \right)} \right|}},$$

(11)

where ${\text{TF}}_{22M,d}$ and ${\text{TF}}_{2M,d}$ denote the standard complex IS-TF and 2-month complex IS-TF, respectively. $E_{22M,d}$ and $E_{2M,d}$ denote the estimated errors at the 95% confidence level based on the error propagation under the assumption of a Gaussian distribution. If ${\text{TFD}}$ is greater than 1, the difference between the 2-month IS-TF and the standard one is greater than the estimated error.

Almost all Basis vectors satisfy Eq. 10 for $\varTheta$ less than 0.01, but do not satisfy Eq. 10 for $\varTheta$ greater than 0.20. Thus, we detect the anomalous events based on Eq. 10 by changing the threshold in the range from 0.01 to 0.20 in intervals of 0.01, and derive the IS-TFs by removing the time windows that have the maximum amplitude at each frequency in the Activations $U\left( {k,t} \right)$. These Activations correspond to the Basis vectors satisfying the two conditions stated above: (i) the ${\text{BR}}$ is greater than 10% and (ii) Eq. 10 holds.

For a quantitative evaluation, we separate the frequency into three ranges—low: 9/30,720–41/30,720 Hz, middle: 42/30,720–74/30,720 Hz, and high: 75/30,720–108/30,720 Hz. In each range, we calculated the averaged ${\text{TFD}}$ values for the 2-month $T_{xx}$ and $T_{yy}$ derived from the raw data as follows—low: 1.50 ($T_{xx}$) and 0.84 ($T_{yy}$), middle: 1.87 ($T_{xx}$) and 0.88 ($T_{yy}$), and high: 1.42 ($T_{xx}$) and 1.01 ($T_{yy}$). The averaged ${\text{TFD}}$ values modified by removing the time windows as above are less than those values as far as substituting 0.01–0.20 into $\varTheta$ in Eq. 10. Based on the trial-and-error approach, the best averaged ${\text{TFD}}$ values, obtained with $\varTheta = 0.04$, are as follows—low: 0.75 ($T_{xx}$) and 0.76 ($T_{yy}$), middle: 0.83 ($T_{xx}$) and 0.62 ($T_{yy}$), and high: 1.36 ($T_{xx}$) and 0.71 ($T_{yy}$). The number of removed windows is dependent on the frequency and varies from 1.2 to 4.2% of the total. These IS-TFs are shown in Fig. 6, and the IS-TFs with removed time windows including anomalous events (Fig. 6c) are similar to the standard IS-TFs (Fig. 6a). The averaged ${\text{TFD}}$ of the off-diagonal components $T_{xy}$ and $T_{yx}$ (not shown) derived from the raw data are as follows—low: 0.78 ($T_{xy}$) and 0.86 ($T_{yx}$), middle: 1.00 ($T_{xy}$) and 1.30 ($T_{yx}$), and high: 0.87 ($T_{xy}$) and 1.32 ($T_{yx}$). The averaged values derived from the modified data are as follows—low: 0.68 ($T_{xy}$) and 0.71 ($T_{yx}$), middle: 0.64 ($T_{xy}$) and 0.70 ($T_{yx}$), and high: 0.89 ($T_{xy}$) and 0.76 ($T_{yx}$). The off-diagonal components have smaller amplitudes than the diagonal components, and so we focus on the diagonal elements in this paper. Based on the ${\text{TFD}}$, we check that some of the IS-TFs are biased following anomalous geomagnetic events, and evaluate such events using MC-NMF (i.e., based on Eq. 9). Note that the ${\text{TFD}}$ includes both amplitude and phase information, which suggests that such anomalous events affect the phases although we focus on the amplitudes here.

Fig. 6

Amplitudes of IS-TFs of MMB (input) to KAK (output). a Standard IS-TFs based on 22 months of data. b 2-month IS-TFs. c Two-month IS-TFs derived by removing time windows including anomalous geomagnetic events. Error bars correspond to 95% confidence level

Ordinarily, the geomagnetic amplitude at MMB is greater than those at KAK or KNY, as shown in Fig. 5b. Focusing on the time-series in the x-direction reconstructed from $\varvec{Y}_{m,x} \left( 0 \right)$ (the top three in Fig. 5a), the largest geomagnetic peak occurs at around 70,000 s at KNY, and the peaks at KAK and MMB are almost the same. Therefore, Basis vector 0 can be considered to reflect the large-power anomalous geomagnetic event at KNY. Moreover, based on the locations of the three observatories (Fig. 2), this event can be regarded as more powerful in the south than the other events. This interpretation does not conflict with the 2-month IS-TFs (Fig. 6b), whose $T_{xx}$ components indicate larger values than those of the standard IS-TFs (Fig. 6a).

Year-to-year changes in IS-TFs derived from array data

To evaluate the spatial gradient (geographically north and east) of each geomagnetic event and delineate the effect on the IS-TFs in more detail, we now analyze the horizontal geomagnetic data acquired at four (geo)magnetic observatories (Beijing Ming Tombs (BMT), KAK, KNY, and MMB), as shown in Fig. 2. The data were recorded during October and November from 2000 to 2009 at a sampling rate of 1/60 Hz. These data can be downloaded from INTERMAGNET (http://www.intermagnet.org/index-eng.php, accessed December 26, 2019) except for the data recorded at KNY in 2000, which are available from the Kakioka magnetic observatory website (http://www.kakioka-jma.go.jp/, accessed December 26, 2019). We selected these 2 months because there was no lack of data from all observatories (BMT, KAK, KNY, and MMB).

Based on Eq. 2, we calculate the IS-TFs during each span of 2 months from 2000 to 2009 under the following six cases: (a) $T_{xx}$ and $T_{yy}$ of KNY (input) to KAK (output) defined as $T_{xx}$ and $T_{yy}$ for KAK/KNY, (b) $T_{xx}$ and $T_{yy}$ for KNY/KAK, (c) $T_{xx}$ and $T_{yy}$ for KAK/MMB, (d) $T_{xx}$ and $T_{yy}$ for MMB/KAK, (e) $T_{xx}$ and $T_{yy}$ for KNY/MMB, and (f) $T_{xx}$ and $T_{yy}$ for MMB/KNY. Generally, MT responses and IS-TFs at low frequencies are biased from the source field more strongly than those at high frequencies (Egbert et al. 2000). Thus, we focus on 9/30,720 Hz, which is the lowest frequency considered in this study. The year-to-year changes in the IS-TF amplitudes at a frequency of 9/30,720 Hz are shown in Fig. 7a–f, which correspond to cases (a)–(f) mentioned above. Note that the vertical axes in Fig. 7b, d, f are reversed because these figures are derived by swapping the input and output data of cases (a), (c), and (e), respectively. The confidence level of the estimated errors is 95%. Comparing Fig. 7a, b, the IS-TFs for the KAK and KNY data exhibit a reasonable opposite polarity (i.e., negative correlation), although it appears to be identical because of the reversed vertical axis of Fig. 7b. However, several exceptions appear in the IS-TFs between KAK and MMB (Fig. 7c, d) and between KNY and MMB (Fig. 7e, f). For example, $T_{xx}$ for KAK/MMB in 2004 shifts above the estimated errors in 2002, 2005, and 2007, although $T_{xx}$ for MMB/KAK remains within the estimated errors. Moreover, $T_{xx}$ for KAK/MMB and $T_{xx}$ for MMB/KAK are shifting with an identical polarity through 2003–2004. In Table 1, we summarize all exceptions for which (1) the IS-TFs shift temporally above or up to their estimated errors, although the values derived by swapping the input and output data do not, and (2) the year-to-year shift in the IS-TF has the same polarity as that given by swapping the input and output data.

Fig. 7

IS-TFs at 9/30,720 Hz derived from the data during October/November from 2000 to 2009. a IS-TFs for KAK/KNY. The vertical axis denotes the absolute amplitude of IS-TFs and the horizontal axis denotes years. b KNY/KAK. c KAK/MMB. d MMB/KAK. e KNY/MMB. f MMB/KNY. The red points denote $T_{xx}$ and the light-blue points denote $T_{yy}$. Note that the vertical axes of b, d, and f are reversed. Error bars correspond to 95% confidence level

Table 1 IS-TFs that shift with the same polarity as those derived by swapping the input and output data

Also, we annotate them in Fig. 7.

To evaluate the anomalous geomagnetic events that have characteristics of spatial gradients different from the other events, we apply MC-NMF to the horizontal geomagnetic data from the four observatories (BMT, KAK, KNY, and MMB). In the application of MC-NMF, the detailed STFT condition is described in the previous section, and the analyzed frequency range is from 9/30,720 Hz to 108/30,720 Hz. We obtain the eight geomagnetic spectrograms during October/November of each year, and construct a model with 10 Basis vectors and 10 Activations. The objective function in Eq. 6 is optimized through 3000 iterations, and the RMSE (Eq. 8) is around 3%. We derive ${\text{BR}}_{m,d} \left( {f_{1} ,k} \right)$ ($m =$ BMT, KAK, KNY, or MMB, $d = x,y$) in Eq. 9, where $f_{1} =$ 9/30,720 Hz, from the result in each year. The results are summarized in Fig. 8(a). The Sum of Squared Errors (${\text{SSE}}$) is defined as

$${\text{SSE}}_{d} \left( {m_{1} ,m_{2} } \right) = \mathop \sum \limits_{k} \left| {{\text{BR}}_{{m_{1} ,d}} \left( {f_{1} ,k} \right) - {\text{BR}}_{{m_{2} ,d}} \left( {f_{1} ,k} \right)} \right|^{2} ,$$

(12)

where $m_{1} \ne m_{2}$ and $m_{1} ,m_{2} =$ KAK, KNY, or MMB. The ${\text{SSE}}$ values between each pair of sites are summarized in Fig. 8b. A large ${\text{SSE}}$ indicates geomagnetic data including anomalous events. Figure 8(b) shows that ${\text{SSE}}_{x} \left( {{\text{KNY}},{\text{MMB}}} \right)$ in 2004 is greater than the ${\text{SSE}}_{x} \left( {{\text{KNY}},{\text{MMB}}} \right)$ in the other years, and that ${\text{SSE}}_{x} \left( {{\text{KAK}},{\text{MMB}}} \right)$ in 2004 is the second largest across all years. These large ${\text{SSE}}_{x}$ in 2004 are triggered by Basis vectors 1 and 2, of which the ${\text{BR}}$ values contribute more than 85% of ${\text{SSE}}_{x}$. Moreover, any other Basis vector contributes less than 10%. Thus, we focus on these two Basis vectors, and summarize the corresponding ${\text{BR}}$ values in Table 2.

Fig. 8

Results obtained applying MC-NMF to geomagnetic spectrograms from four observatories. a${\text{BR}}$ at 9/30,720 Hz (left: BMT, left mid: KAK, right mid: KNY, and right: MMB). The uppers show ${\text{BR}}$ of the x component and the lowers show ${\text{BR}}$ of the y component. The vertical axis denotes years during 2000–2009 (the standard is 2000), and the horizontal axis denotes IDs of Basis vectors. b The ${\text{SSE}}\left( {{\text{KAK}},{\text{KNY}}} \right)$ (red), ${\text{SSE}}\left( {{\text{KAK}},{\text{MMB}}} \right)$ (light blue), and ${\text{SSE}}\left( {{\text{KNY}},{\text{MMB}}} \right)$ (orange). The upper denotes the ${\text{SSE}}$ of the x component and the lower denotes the ${\text{SSE}}$ of the y component. The vertical axis denotes ${\text{SSE}}$ and the horizontal axis denotes years

Table 2 ${\text{BR}}_{m,x} \left( 1 \right)$ and ${\text{BR}}_{m,x} \left( 2 \right)$ of BMT, KAK, KNY, and MMB in 2004

Basis vector 1 reflects a large-power anomalous event in the south (i.e., KAK and KNY as shown in Fig. 2) and Basis vector 2 reflects a large-power event in the north (i.e., BMT and MMB). In addition, ${\text{SSE}}_{y} \left( {{\text{KNY}},{\text{MMB}}} \right)$ in 2005 has a large value, and several anomalous events are included as well as in the data for the x direction in 2004.

Discussion

We now discuss the implications of the IS-TFs summarized in Table 1. The effects of anomalous geomagnetic events that have characteristics of spatial gradients different from the others on the IS-TFs are then elucidated.

We focus on the $T_{xx}$ result between KAK and MMB in 2004 (Fig. 7c, d). Based on $T_{xx}$ for MMB/KAK, the spatial gradient of geomagnetic temporal variations in 2004 seems to be the same as in 2002, 2005, and 2007. However, based on $T_{xx}$ for KAK/MMB, the power of the geomagnetic variations at KAK in 2004 appears to be larger than in the years. Therefore, $T_{xx}$ between KAK and MMB in 2004 does not give a precise reflection of the spatial gradients of time-varying geomagnetic fields. The other IS-TFs in Table 1 are similar to the above.

Based on the IS-TFs between sites 1 and 2 derived using Eq. 2, the reason why the IS-TFs do not exactly reflect the spatial gradients of geomagnetic variations can be explained as follows. For simplicity, we assume that the cross-spectra between the x and y components are smaller than those in the same direction, and define $T_{xx}$ for site 2 with respect to site 1 as

$$T_{xx} \left( {1,2} \right) \cong \frac{{\left\langle {H_{{{\text{site}}1,x}} H_{{{\text{site}}2,x}}^{*} } \right\rangle \left\langle {H_{{{\text{site}}2,y}} H_{{{\text{site}}2,y}}^{*} } \right\rangle }}{{\left\langle {H_{{{\text{site}}2,x}} H_{{{\text{site}}2,x}}^{*} } \right\rangle \left\langle {H_{{{\text{site}}2,y}} H_{{{\text{site}}2,y}}^{*} } \right\rangle }} = \frac{{\left\langle {H_{{{\text{site}}1,x}} H_{{{\text{site}}2,x}}^{*} } \right\rangle }}{{\left\langle {H_{{{\text{site}}2,x}} H_{{{\text{site}}2,x}}^{*} } \right\rangle }}.$$

(13)

We also assume that the geomagnetic data $H_{{{\text{site}}1,x}}$ and $H_{{{\text{site}}2,x}}$ include $I$ geomagnetic events:

$$\begin{aligned} H_{{{\text{site}}1,x}} \left( {f,t} \right) = C_{{{\text{site}}1,x}} \left( {f,0} \right)A\left( {0,f,t} \right) + \ldots + C_{{{\text{site}}1,x}} \left( {f,I - 1} \right)A\left( {I - 1,f,t} \right), \hfill \\ H_{{{\text{site}}2,x}} \left( {f,t} \right) = C_{{{\text{site}}2,x}} \left( {f,0} \right)A\left( {0,f,t} \right) + \ldots + C_{{{\text{site}}2,x}} \left( {f,I - 1} \right)A\left( {I - 1,f,t} \right), \hfill \\ \end{aligned}$$

(14)

where $C_{{{\text{site}}1,x}} \left( {f,0} \right), \ldots , C_{{{\text{site}}2,x}} \left( {f,I - 1} \right)$ are the coefficients of geomagnetic events $A\left( {0,t} \right), \ldots ,A\left( {I - 1,t} \right)$ in the geomagnetic data; the equations in Appendix 2 are also considered. The values of $C$ are dependent on (1) the positional relationship between the observed stations and the source field and (2) the resistivity structure of the subsurface at each site. Therefore, if all geomagnetic events have the same spatial gradients, then $C_{{{\text{site}}1,x}} \left( {f,0} \right) = \ldots = C_{{{\text{site}}1,x}} \left( {f,I - 1} \right)$. Here, each geomagnetic event is assumed to be non-overlapping (see Appendix 2). This enables us to consider

$$\begin{aligned} \mathop \sum \limits_{t} A\left( {i_{1} ,f,t} \right)A^{*} \left( {i_{2} ,f,t} \right) = {\text{AP}}\left( {i_{1} ,f} \right) \left( {i_{1} = i_{2} } \right), \hfill \\ \mathop \sum \limits_{t} A\left( {i_{1} ,f,t} \right)A^{*} \left( {i_{2} ,f,t} \right) = 0 \left( {i_{1} \ne i_{2} } \right), \hfill \\ \end{aligned}$$

(15)

where ${\text{AP}}\left( {i,f} \right)$ denotes the power of the $i$th geomagnetic event. As a result, the IS-TF $T_{xx} \left( {1,2} \right)$ from site 2 to site 1 in Eq. 13 can be represented as

$$\begin{aligned} T_{xx} \left( {1,2} \right) & \cong \frac{{\left\langle {H_{{{\text{site}}1,x}} H_{{{\text{site}}2,x}}^{*} } \right\rangle }}{{\left\langle {H_{{{\text{site}}2,x}} H_{{{\text{site}}2,x}}^{*} } \right\rangle }} \\ & = \frac{{\mathop \sum \nolimits_{i} C_{{{\text{site}}1,x}} \left( i \right)A\left( {i,t} \right)\mathop \sum \nolimits_{i} C_{{{\text{site}}2,x}}^{*} \left( i \right)A\left( {i,t} \right)}}{{\mathop \sum \nolimits_{i} C_{{{\text{site}}2,x}} \left( i \right)A\left( {i,t} \right)\mathop \sum \nolimits_{i} C_{{{\text{site}}2,x}}^{*} \left( i \right)A\left( {i,t} \right)}} = \frac{{\mathop \sum \nolimits_{i} C_{{{\text{site}}1,x}} \left( i \right)C_{{{\text{site}}2,x}}^{*} \left( i \right){\text{AP}}\left( i \right)}}{{\mathop \sum \nolimits_{i} C_{{{\text{site}}2,x}} \left( i \right)C_{{{\text{site}}2,x}}^{*} \left( i \right){\text{AP}}\left( i \right)}}, \\ \end{aligned}$$

(16)

where $f$ was omitted to focus on only one frequency. It can be established that $\frac{{C_{{{\text{site}}1,x}} \left( i \right)}}{{C_{{{\text{site}}2,x}} \left( i \right)}} = {\text{constant }}\left( {i = 0, \ldots ,I - 1} \right)$ if all $I$ geomagnetic events have the same spatial gradients at sites 1 and 2. We test the stability of the IS-TF under the simple cases that (i) the power of each geomagnetic event is the same (i.e., ${\text{AP}}\left( 0 \right) = \cdots = {\text{AP}}\left( {I - 1} \right)$), (ii) the number of geomagnetic events is two, (iii) the subsurface structures of sites 1 and 2 are the same, and (iv) the geomagnetic events affect only the amplitude of the IS-TF (i.e., all $C$ in Eq. 16 are real-positive numbers). Using Eq. 16, we calculate $T_{xx} \left( {1,2} \right)$ and $T_{xx} \left( {2,1} \right)$, which are derived by swapping the input and output data of $T_{xx} \left( {1,2} \right)$, and changing the values of $C_{{{\text{site}}1,x}} \left( 0 \right)$, $C_{{{\text{site}}1,x}} \left( 1 \right)$, $C_{{{\text{site}}2,x}} \left( 0 \right)$, and $C_{{{\text{site}}2,x}} \left( 1 \right)$, as summarized in Table 3. We consider four cases: case A: all geomagnetic events are homogeneous and have the same spatial gradients between two sites; case B: all geomagnetic events have different spatial gradients at site 1, but the same spatial gradients at site 2; case C: all geomagnetic events have different spatial gradients between the two sites (event 0 causes variations near site 2 and event 1 causes variations near site 1); and case D: all geomagnetic events produce different spatial gradients between the two sites.

Table 3 Summation of $C_{{{\text{site}}1,x}} \left( 0 \right)$, $C_{{{\text{site}}1,x}} \left( 1 \right)$, $C_{{{\text{site}}2,x}} \left( 0 \right)$, and $C_{{{\text{site}}2,x}} \left( 1 \right)$, and the resultant values of $T_{xx} \left( {1,2} \right)$ and $T_{xx} \left( {2,1} \right)$

We use the IS-TFs derived from case A (i.e., $T_{xx} \left( {1,2} \right) = T_{xx} \left( {2,1} \right) = 1.00$) as the standard values. $T_{xx} \left( {1,2} \right)$ has a value of 1.00 in cases B and D, although these two situations are different: focusing on site 2, in case B, all events have the same gradients as in case A, but in case D, all events have different spatial gradients. As a result, we cannot distinguish $T_{xx} \left( {1,2} \right)$ in case B and in case D. Both $T_{xx} \left( {1,2} \right)$ and $T_{xx} \left( {2,1} \right)$ in case C are less than 1.00. In this case, $T_{xx} \left( {1,2} \right)$ implies that the geomagnetic temporal variations at site 2 are larger than those at site 1; in contrast, $T_{xx} \left( {2,1} \right)$ implies the opposite. Although $T_{xx} \left( {1,2} \right)$ and $T_{xx} \left( {2,1} \right)$ should have opposite polarities, both of these values in case C decrease, and have the same polarity when case A is used as the standard. However, using case D as the standard, $T_{xx} \left( {1,2} \right)$ and $T_{xx} \left( {2,1} \right)$ in case C have opposite polarities. The results of cases B–D indicate that IS-TFs do not reflect the exact situation regarding the spatial gradients of geomagnetic temporal variations if the analyzed data include several anomalous events. The reason is simply because we generally estimate the IS-TFs without considering the mixture model in statistics, and the mathematical background for such a bias in IS-TFs is explained in Appendix 3. In addition, we use a long time-series to derive the IS-TFs based on stacking, and long-span data are likely to include several anomalous events.

The MC-NMF results show that $H_{x}$ in 2004 includes several anomalous events, as summarized in Table 2. We reconstructed the geomagnetic data at KAK and MMB using only (a) Basis vector 1 as $Y_{m,d} \left( {1,f_{1} ,t} \right) = B_{m,d} \left( {f_{1} ,1} \right)U\left( {1,t} \right)e^{{j\varphi_{m,d} \left( {1,f_{1} ,t} \right)}}$, (b) Basis vector 2 as $Y_{m,d} \left( {2,f_{1} ,t} \right)$, and (c) Basis vectors 1 and 2 as $Y_{m,d} \left( {1,f_{1} ,t} \right) + Y_{m,d} \left( {2,f_{1} ,t} \right)$. Then, $T_{xx}$ for KAK/MMB and $T_{xx}$ for MMB/KAK were derived using the reconstructed geomagnetic data of cases (a)–(c), as shown in Fig. 9. $T_{xx}$ for KAK/MMB and $T_{xx}$ for MMB/KAK in cases (a) and (c) have identical polarities, whereas both $T_{xx}$ in case (b) are somewhat different from those derived from the raw data. The anomalous geomagnetic events corresponding to Basis vectors 1 and 2 affect the IS-TFs. These anomalous events are considered to trigger the inconsistent implications between $T_{xx}$ for KAK/MMB and for MMB/KAK in 2004, as explained above using the simulation in Table 3.

Fig. 9

Year-to-year $T_{yy}$ changes between KAK and MMB at a frequency of 9/30,720 Hz. a$T_{xx}$ for KAK/MMB. The circles were derived from the raw geomagnetic data and correspond to the light-blue dots in Fig. 7c. The diamonds, squares, and triangles in 2004 are derived from the reconstructed data using only Basis vector 1, Basis vector 2, and Basis vectors 1 and 2, respectively. b$T_{xx}$ for MMB/KAK. The circles were derived from raw geomagnetic data, and correspond to the light-blue dots in Fig. 7d. Note that the directions of the vertical axes in a and b are reversed

The other cases of $T_{xx}$ in 2004 and $T_{yy}$ in 2005 between KNY and MMB (see Fig. 7e, f) can also be explained by the arguments above. From these results, the data analyzed using MC-NMF can be considered to include anomalous geomagnetic events when the IS-TFs do not exactly reflect the spatial gradients of geomagnetic temporal variations.

However, the inverse has not been established. For example, we focus on the reasonable opposite polarity of $T_{yy}$ between KNY/MMB and MMB/KNY in 2004 under large ${\text{SSE}}$. The large ${\text{SSE}}\left( {{\text{KNY}},{\text{MMB}}} \right)$ is caused by Basis vectors 2 and 6, whose ${\text{BR}}$ values are summarized in Table 4. Note that the other Basis vectors do not contribute to the large ${\text{SSE}}\left( {{\text{KNY}},{\text{MMB}}} \right)$.

Table 4 ${\text{BR}}_{m,y} \left( 2 \right)$ and ${\text{BR}}_{m,y} \left( 6 \right)$ of BMT, KAK, KNY, and MMB in 2004

Basis vector 2 is more conspicuous in the data at BMT and KNY than at KAK and MMB, whereas Basis vector 6 is more dominant in the data at KAK and MMB than at BMT and KNY. Considering the location of each observatory (Fig. 2), Basis vector 2 reflects a large-power anomalous geomagnetic event in the west, whereas Basis vector 6 reflects a large-power event in the east.

Here, we verify the effect of Basis vectors 2 and 6 on $T_{yy}$, as shown in Fig. 9. The geomagnetic data at KNY and MMB are reconstructed using only (a) Basis vector 2 as $Y_{m,d} \left( {2,f_{1} ,t} \right)$, (b) Basis vector 6 as $Y_{m,d} \left( {6,f_{1} ,t} \right)$, and (c) Basis vectors 2 and 6 as $Y_{m,d} \left( {2,f_{1} ,t} \right) + Y_{m,d} \left( {6,f_{1} ,t} \right)$. Then, we calculate $T_{yy}$ for KNY/MMB and $T_{yy}$ for MMB/KNY using the reconstructed geomagnetic data in cases (a)–(c), and show the results in Fig. 10. Figure 10a, b shows that the IS-TFs derived from the data reconstructed using only Basis vector 2 or 6 (the diamonds and squares in Fig. 10) are largely different from the raw data in each year. Therefore, Basis vectors 2 and 6 are actually unstable components for the IS-TFs. In addition, the $T_{yy}$ for KNY/MMB derived from Basis vectors 2 and 6 (triangle in Fig. 10a) does not shift from 2003 and 2005, although $T_{yy}$ for MMB/KNY (triangle in Fig. 10b) shifts above the estimated errors from 2003 and 2005. This corresponds to the case whereby geomagnetic data include several anomalous events, as summarized in Table 3. However, $T_{yy}$ for KNY/MMB and MMB/KNY in 2004, derived from the raw data, exhibit a shift with an opposite polarity and seem to reflect the spatial gradients of geomagnetic temporal variations exactly. Basis vectors 2 and 6 have opposite spatial characteristics, as mentioned above, and can be considered to cancel each other out. Moreover, the sum of ${\text{BR}}$ values of Basis vectors 2 and 6 is smaller than 20% of the total, and is smaller than that of Basis vectors 1 and 2 in the data for the x direction in 2004. As a result, their effect on the IS-TFs becomes small, as shown in Fig. 10.

Fig. 10

Year-to-year changes in $T_{yy}$ between KNY and MMB at a frequency of 9/30,720 Hz. a$T_{yy}$ for KNY/MMB. The circles are derived from raw geomagnetic data, and correspond to the light-blue dots in Fig. 7e. The diamonds/squares/triangles in 2004 are derived from the reconstructed data using only Basis vector 2/Basis vector 6/Basis vectors 2 and 6, respectively. b$T_{yy}$ for MMB/KNY. The circles are derived from raw geomagnetic data, and correspond to the light-blue dots in Fig. 7f. Note that the vertical axis of b is reversed

$T_{yy}$ in 2004 can be considered as implying that the geomagnetic variations have a large power in the east (or at MMB), although they include anomalous events. If Basis vectors 2 and 6 had larger powers (e.g., the same as Basis vectors 1 and 2 of $H_{x}$ in 2004), the implication could possibly be different. Note that the other IS-TFs (e.g., $T_{yy}$ in 2008 and 2009) having a reasonable opposite polarity with those derived by swapping the input and output data under the large ${\text{SSE}}$ can be explained as follows. The anomalous events cancel each other out as well as the case of $T_{yy}$ in 2004, and moreover, the powers of such events, causing large ${\text{SSE}}$, are smaller than those in 2004.

Although the analyzed datasets are limited, we confirm specifically the instability of the IS-TFs caused by the anomalous geomagnetic events. Such anomalous events can be detected by MC-NMF. When all geomagnetic events have the same spatial gradient (i.e., no anomalous events exist), MC-NMF will show that the number of spatial gradients is one, and the IS-TFs provide an exact implication regarding the spatial gradient of geomagnetic variations. In discussing the spatial gradients of time-varying geomagnetic fields, the IS-TFs with MC-NMF are required. We can identify the geomagnetic conditions related to the source field (i.e., the establishment or not of the plane-wave assumption) through combination with MC-NMF. Consequently, better inferences regarding the subsurface resistivity structure can be obtained through IS-TFs. Although we focus on the amplitudes here, analyzing the phases of IS-TFs using MC-NMF would enhance our understanding of the spatial gradients of the geomagnetic variations.

The components extracted by MC-NMF may be seen in the raw spectrograms. One could argue that it would be better and easier to compare the raw spectrograms directly when evaluating the presence of anomalous geomagnetic events. However, the spectrograms contain huge amounts of information; for example, there are 132,000 elements in those from BMT, KAK, KNY, and MMB. Thus, they are not suitable for the comparison. Instead, MC-NMF can extract anomalous events from such huge matrices in a quantitative manner.

Summary

We developed a method for extracting the components included in several geomagnetic spectrograms (MC-NMF). Using the proposed method, we can evaluate anomalous geomagnetic events that have characteristics of spatial gradients different from others.

The IS-TFs between KAK and MMB changed temporally from the standard ones, and analysis using MC-NMF showed that the causes are anomalous geomagnetic events. The shift in IS-TFs could be modified by removing the time windows including such events. This ensures that the IS-TFs shift in accordance with the spatial gradients of geomagnetic events.

MC-NMF was applied to geomagnetic data acquired at BMT, KAK, KNY, and MMB, and used to evaluate the anomalous geomagnetic events. We also derived the year-to-year changes in IS-TFs between KAK and KNY, between KAK and MMB, and between KNY and MMB. Some IS-TFs exhibited the same polarity as those derived by swapping the output and input data, although the polarities should be opposite. The spatial gradients of time-varying geomagnetic fields cannot be evaluated from the IS-TFs. However, using numerical examples we proved that this is because of anomalous geomagnetic events, although the proof assumes a specific condition, and the results of MC-NMF showed that the analyzed data actually include such events.

The IS-TFs can fail to yield the exact implications regarding the spatial gradients of geomagnetic temporal variations when the analyzed data include anomalous geomagnetic events. However, MC-NMF can evaluate such anomalous events. Using IS-TFs alongside MC-NMF allows information on geomagnetic conditions, such as their spatial gradients, to be obtained. This advantage will be useful for checking the establishment of the plane-wave assumption, and as a result, will yield better implications related to subsurface resistivity information. We also will discuss the effect of spatial gradients of geomagnetic temporal variations on the MT responses or geomagnetic depth-sounding responses.

Appendices

Appendices

Appendix 1

Here, we derive the optimized equation of Basis vectors and Activations by minimizing the objective functions in MC-NMF. Hunter and Lange (2004) suggested the MM algorithm, which minimizes the objective function $J\left( \alpha \right)$ as follows: Consider $J^{ + } \left( {\alpha ,\beta } \right)$ such that $\mathop {\hbox{min} }\limits_{\beta } J^{ + } \left( {\alpha ,\beta } \right) = J\left( \alpha \right)$. Then, $J\left( \alpha \right)$ is not increasing under the update $\beta \leftarrow {\text{argmin}}_{\beta } J^{ + } \left( {\alpha ,\beta } \right)$ and so $\alpha \leftarrow {\text{argmin}}_{\alpha } J^{ + } \left( {\alpha ,\beta } \right)$.

Proof

Let $\beta_{\theta + 1} = {\text{argmin}}_{\beta } J^{ + } \left( {x,\beta } \right)$ and $\alpha_{\theta + 1} = {\text{argmin}}_{\alpha } J^{ + } \left( {\alpha ,\beta } \right)$. The relationship $J\left( {\alpha_{\theta } } \right) = J^{ + } \left( {\alpha_{\theta } ,\beta_{\theta + 1} } \right)$ at the update from $\theta$ to $\theta + 1$ follows the definition of $\mathop {\hbox{min} }\limits_{\beta } J^{ + } \left( {\alpha ,\beta } \right) = J\left( \alpha \right)$. In addition, $\alpha_{\theta + 1} = {\text{argmin}}_{\alpha } J^{ + } \left( {\alpha ,\beta_{\theta + 1} } \right)$ implies that $J^{ + } \left( {\alpha_{\theta } ,\beta_{\theta + 1} } \right) \ge J^{ + } \left( {\alpha_{\theta + 1} ,\beta_{\theta + 1} } \right)$. Therefore, we can obtain $J\left( {\alpha_{\theta } } \right) = J^{ + } \left( {\alpha_{\theta } ,\beta_{\theta + 1} } \right) \ge J^{ + } \left( {\alpha_{\theta + 1} ,\beta_{\theta + 1} } \right) = J\left( {\alpha_{\theta + 1} } \right)$ from $\theta$ to $\theta + 1$. $\square$

The objective function in Eq. 6 is $J\left( {B,U,\varphi } \right) = \mathop \sum \limits_{m,d,f,t} \left| {X_{m,d} \left( {f,t} \right) - \mathop \sum \limits_{k} B_{m,d} \left( {f,k} \right)U\left( {k,t} \right){\text{e}}^{{j\varphi_{m,d} \left( {k,f,t} \right)}} } \right|^{2} + 2\lambda \mathop \sum \limits_{k,t} \left| {U\left( {k,t} \right)} \right|^{q}$, and we require $J^{ + } \left( {B,U,\varphi ,y} \right)$ such that $\mathop {\hbox{min} }\limits_{y} J^{ + } \left( {B,U,\varphi ,y} \right) = J\left( {B,U,\varphi } \right)$. From Jensen’s inequality, we obtain the relationship $\mathop \sum \limits_{k} \gamma_{k} G\left( {x_{k} } \right) \ge G\left( {\mathop \sum \limits_{k} \gamma_{k} x_{k} } \right)$ under the condition that $G$ is a convex function and $\mathop \sum \limits_{k} \gamma_{k} = 1 \left( {\gamma_{k} > 0} \right)$. The first term in Eq. 6 can be transformed into

$$\begin{aligned} & \left| {X_{m,d} \left( {f,t} \right) - \mathop \sum \limits_{k} B_{m,d} \left( {f,k} \right)U\left( {k,t} \right){\text{e}}^{{j\varphi_{m,d} \left( {k,f,t} \right)}} } \right|^{2} \\ & \quad = \left| {\mathop \sum \limits_{k} \frac{1}{K}\frac{1}{{\frac{1}{K}}}\left\{ {X_{m,d}^{ + } \left( {k,f,t} \right) - B_{m,d} \left( {f,k} \right)U\left( {k,t} \right){\text{e}}^{{j\varphi_{m,d} \left( {k,f,t} \right)}} } \right\}} \right|^{2} \\ & \quad \quad \le \mathop \sum \limits_{k} \frac{1}{K}\left| {\frac{{X_{m,d}^{ + } \left( {k,f,t} \right) - B_{m,d} \left( {f,k} \right)U\left( {k,t} \right){\text{e}}^{{j\varphi_{m,d} \left( {k,f,t} \right)}} }}{{\frac{1}{K}}}} \right|^{2} \\ & \quad = \mathop \sum \limits_{k} K\left| {X_{m,d}^{ + } \left( {k,f,t} \right) - B_{m,d} \left( {f,k} \right)U\left( {k,t} \right){\text{e}}^{{j\varphi_{m,d} \left( {k,f,t} \right)}} } \right|^{2} , \\ \end{aligned}$$

(17)

where $\mathop \sum \limits_{k} X_{m,d}^{ + } \left( {k,f,t} \right)$ = $X_{m,d} \left( {f,t} \right)$. Because $2\left| {U\left( {k,t} \right)} \right|^{q}$ is tangent to $q\left| {U^{ + } \left( {k,t} \right)} \right|^{q - 2} \left| {U\left( {k,t} \right)} \right|^{2} + 2\left| {U^{ + } \left( {k,t} \right)} \right|^{q} - q\left| {U^{ + } \left( {k,t} \right)} \right|^{q}$ at $U\left( {k,t} \right) = \pm U^{ + } \left( {k,t} \right)$, the second term in Eq. 6 can be transformed into

$$2\left| {U\left( {k,t} \right)} \right|^{q} \le q\left| {U^{ + } \left( {k,t} \right)} \right|^{q - 2} \left| {U\left( {k,t} \right)} \right|^{2} + 2\left| {U^{ + } \left( {k,t} \right)} \right|^{q} - q\left| {U^{ + } \left( {k,t} \right)} \right|^{q} ,$$

(18)

where $0 < q < 2$ and $U^{ + } \left( {k,t} \right)$ is arbitrary. Therefore, we obtain $J^{ + } \left( {B,U,\varphi ,X^{ + } ,U^{ + } } \right)$ as

$$\begin{aligned} & J^{ + } \left( {B,U,\varphi ,X^{ + } ,U^{ + } } \right) \\ & \quad = \mathop \sum \limits_{k,d,f,t} K\left| {X_{m,d}^{ + } \left( {k,f,t} \right) - B_{m,d} \left( {f,k} \right)U\left( {k,t} \right){\text{e}}^{{j\varphi_{m,d} \left( {k,f,t} \right)}} } \right|^{2} \\ & \quad \quad + \lambda \mathop \sum \limits_{k,t} \left\{ {q\left| {U^{ + } \left( {k,t} \right)} \right|^{q - 2} \left| {U\left( {k,t} \right)} \right|^{2} + 2\left| {U^{ + } \left( {k,t} \right)} \right|^{q} - q\left| {U^{ + } \left( {k,t} \right)} \right|^{q} } \right\}. \\ \end{aligned}$$

(19)

The equality $J^{ + } \left( {B,U,\varphi ,X^{ + } ,U^{ + } } \right) = J\left( {B,U,\varphi } \right)$ can be established when the equalities in Eq. 17 and Eq. 18 are established as

$$\begin{aligned} & X_{m,d}^{ + } \left( {k,f,t} \right) \\ & \quad \leftarrow B_{m,d} \left( {f,k} \right)U\left( {k,t} \right){\text{e}}^{{j\varphi_{m,d} \left( {k,f,t} \right)}} + \frac{1}{K}\left\{ {\left| {X_{m,d} \left( {f,t} \right)} \right| - \mathop \sum \limits_{k} B_{m,d} \left( {f,k} \right)U\left( {k,t} \right)} \right\}, \\ \end{aligned}$$

(20)

$$U^{ + } \left( {k,t} \right) \leftarrow U\left( {k,t} \right).$$

(21)

$J\left( {B,U,\varphi } \right)$ can be minimized after updating $X_{m,d}^{ + } \left( {k,f,t} \right)$ and $U^{ + } \left( {k,t} \right)$ as shown in Eqs. 20 and 21. We need only calculate $\frac{{\partial J^{ + } \left( {B,U,\varphi ,X^{ + } ,U^{ + } } \right)}}{{\partial B_{m,d} \left( {f,k} \right)}}$ and $\frac{{\partial J^{ + } \left( {B,U,\varphi ,X^{ + } ,U^{ + } } \right)}}{{\partial U\left( {k,t} \right)}}$, and under the “nonnegative” constrain, the updates of $B_{m,d} \left( {f,k} \right)$ and $U\left( {k,t} \right)$ are as follows:

$$B_{m,d} \left( {f,k} \right) \leftarrow \frac{{\mathop \sum \nolimits_{t} U\left( {k,t} \right)\left| {X_{m,d}^{ + } \left( {k,f,t} \right)} \right|}}{{\mathop \sum \nolimits_{t} U\left( {k,t} \right)^{2} }},$$

(22)

$$U\left( {k,t} \right) \leftarrow \frac{{\mathop \sum \nolimits_{m,d,f} B_{m,d} \left( {f,k} \right)\left| {X_{m,d}^{ + } \left( {k,f,t} \right)} \right|}}{{\mathop \sum \nolimits_{m,d,f} B_{m,d} \left( {f,k} \right)^{2} + \frac{\lambda }{K}q\left| {U\left( {k,t} \right)} \right|^{q - 2} }}.$$

(23)

The term of $J^{ + }$ regarding $\varphi$ can be transformed as

$$\begin{aligned} & \left| {X_{m,d}^{ + } \left( {k,f,t} \right) - B_{m,d} \left( {f,k} \right)U\left( {k,t} \right){\text{e}}^{{j\varphi_{m,d} \left( {k,f,t} \right)}} } \right|^{2} \\ & \quad = \left| {X_{m,d}^{ + } \left( {k,f,t} \right)} \right|^{2} + \left| {B_{m,d} \left( {f,k} \right)U\left( {k,t} \right)} \right|^{2} \\ & \quad \quad - 2{\text{Real}}\left[ {\left\{ {X_{m,d}^{ + } \left( {k,f,t} \right)} \right\}^{*} B_{m,d} \left( {f,k} \right)U\left( {k,t} \right){\text{e}}^{{j\varphi_{m,d} \left( {k,f,t} \right)}} } \right]. \\ \end{aligned}$$

(24)

Therefore, the phase update is

$${\text{e}}^{{j\varphi_{m,d} \left( {k,f,t} \right)}} \leftarrow \frac{{X_{m,d}^{ + } \left( {k,f,t} \right)}}{{\left| {X_{m,d}^{ + } \left( {k,f,t} \right)} \right|}}.$$

(25)

This optimization rule corresponds to the C-NMF given by Kameoka et al. (2009). In this study, we substitute 1.2 into $q$ because Kameoka et al. (2009) have demonstrated that $q = 1.2$ ensures stability and $\mathop \sum \limits_{m,d,f,t} \frac{{\left| {X_{m,d} \left( {f,t} \right)} \right|^{2} }}{{10^{4.5} }}$ into $\lambda$, which is ten times of the value used by Kameoka et al. (2009), for additional weighting of the super-Gaussian term.

Appendix 2

Assuming that the geomagnetic fields $H_{{{\text{site}}1,x}} \left( {f,t} \right)$ and $H_{{{\text{site}}1,y}} \left( {f,t} \right)$ include $I$ geomagnetic events $A\left( {i,f,t} \right)$$\left( {i = 0, \ldots I - 1} \right)$, $H_{{{\text{site}}1,x}} \left( {f,t} \right)$ and $H_{{{\text{site}}1,y}} \left( {f,t} \right)$ can be represented by the coefficients $C_{{{\text{site}}1,x}} \left( {f,i} \right)$, $C_{{{\text{site}}1,y}} \left( {f,i} \right)$, and $A\left( {i,f,t} \right)$. The geomagnetic fields at other sites can be summarized as

$$\left( {\begin{array}{*{20}c} {H_{{{\text{site}}1,x}} \left( {f,t} \right)} \\ \vdots \\ {H_{{{\text{site}}K,y}} \left( {f,t} \right)} \\ \end{array} } \right) = \left( {\begin{array}{*{20}c} {\begin{array}{*{20}c} {C_{{{\text{site}}1,x}} \left( {f,0} \right)} \\ \vdots \\ {C_{{{\text{site}}K,y}} \left( {f,0} \right)} \\ \end{array} } & \ldots & {\begin{array}{*{20}c} {C_{{{\text{site}}1,x}} \left( {f,I - 1} \right)} \\ \vdots \\ {C_{{{\text{site}}K,y}} \left( {f,I - 1} \right)} \\ \end{array} } \\ \end{array} } \right)\left( {\begin{array}{*{20}c} {A\left( {0,f,t} \right)} \\ \vdots \\ {A\left( {I - 1,f,t} \right)} \\ \end{array} } \right).$$

(26)

We define the geomagnetic events in this study as follows: (i) one $A\left( {i,f,t} \right)$ of frequency $f$ is non-overlapping with the others at $t$ as well as an element of $\varvec{Y}_{{{m},{d}}}$ at a time–frequency slot ($f$, $t$) introduced above, and (ii) each $\varvec{A}\left( i \right)$, whose elements are $A\left( {i,f,t} \right)$, has a constant magnitude spectrum pattern. These geomagnetic events are not “real”, but are defined to allow us to derive the relation with events reflected by the Basis vectors of MC-NMF and to assign physical meanings such as spatial gradients to the Basis vectors. Definition (i) allows us to consider that at most one event has a power at a time window $t$. Thus, we obtain

$$\begin{aligned} & \left( {\begin{array}{*{20}c} {\left| {H_{{{\text{site}}1,x}} \left( {f,t} \right)} \right|} \\ \vdots \\ {\left| {H_{{{\text{site}}M,y}} \left( {f,t} \right)} \right|} \\ \end{array} } \right) \\ & \quad = \left( {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\left| {C_{{{\text{site}}1,x}} \left( {f,0} \right)} \right|} \\ \vdots \\ {\left| {C_{{{\text{site}}M,y}} \left( {f,0} \right)} \right|} \\ \end{array} } & \ldots & {\begin{array}{*{20}c} {\left| {C_{{{\text{site}}1,x}} \left( {f,I - 1} \right)} \right|} \\ \vdots \\ {\left| {C_{{{\text{site}}M,y}} \left( {f,I - 1} \right)} \right|} \\ \end{array} } \\ \end{array} } \right)\left( {\begin{array}{*{20}c} {\left| {A\left( {0,f,t} \right)} \right|} \\ \vdots \\ {\left| {A\left( {I - 1,f,t} \right)} \right|} \\ \end{array} } \right). \\ \end{aligned}$$

(27)

Definition (ii) ensures that the rank of each $\varvec{A}\left( i \right)$ is one, and so $A\left( {i,f,t} \right)$ can be represented as

$$\left| {A\left( {i,f,t} \right)} \right| = A_{s} \left( {f,i} \right)A_{a} \left( {i,t} \right),$$

(28)

where $A_{s} \left( {f,i} \right)$ denotes a constant magnitude spectrum pattern and $A_{a} \left( {i,t} \right)$ denotes the temporal changes of amplitude; both are nonnegative. Definition (ii) also allows us to consider that geomagnetic temporal variations with the same magnitude spectrum pattern are triggered by the same geomagnetic event. Substituting Eq. 28 into Eq. 27, we obtain

$$\left( {\begin{array}{*{20}c} {\left| {H_{{{\text{site}}1,x}} \left( {f,t} \right)} \right|} \\ \vdots \\ {\left| {H_{{{\text{site}}M,y}} \left( {f,t} \right)} \right|} \\ \end{array} } \right) = \left( {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\left| {C_{{{\text{site}}1,x}} \left( {f,0} \right)} \right|} \\ \vdots \\ {\left| {C_{{{\text{site}}M,y}} \left( {f,0} \right)} \right|} \\ \end{array} } & \ldots & {\begin{array}{*{20}c} {\left| {C_{{{\text{site}}1,x}} \left( {f,I - 1} \right)} \right|} \\ \vdots \\ {\left| {C_{{{\text{site}}M,y}} \left( {f,I - 1} \right)} \right|} \\ \end{array} } \\ \end{array} } \right)\left( {\begin{array}{*{20}c} {A_{s} \left( {f,0} \right)A_{a} \left( {0,t} \right)} \\ \vdots \\ {A_{s} \left( {f,I - 1} \right)A_{a} \left( {I - 1,t} \right)} \\ \end{array} } \right).$$

(29)

Substituting $A_{s} \left( {f,i} \right)$ and the power of $A_{a} \left( {i,t} \right)$ (i.e., $\mathop \sum \limits_{t} \left| {A_{a} \left( {i,t} \right)} \right|^{2}$) into $\left| {C_{{{\text{site}}1,x}} \left( {f,i} \right)} \right|$, we can represent Eq. 29 as the MC-NMF model:

$$\left( {\begin{array}{*{20}c} {\left| {H_{{{\text{site}}1,x}} \left( {f,t} \right)} \right|} \\ \vdots \\ {\left| {H_{{{\text{site}}M,y}} \left( {f,t} \right)} \right|} \\ \end{array} } \right) = \left( {\begin{array}{*{20}c} {\begin{array}{*{20}c} {B_{{{\text{site}}1,x}} \left( {f,0} \right)} \\ \vdots \\ {B_{{{\text{site}}M,y}} \left( {f,0} \right)} \\ \end{array} } & \ldots & {\begin{array}{*{20}c} {B_{{{\text{site}}1,x}} \left( {f,I - 1} \right)} \\ \vdots \\ {B_{{{\text{site}}M,y}} \left( {f,I - 1} \right)} \\ \end{array} } \\ \end{array} } \right)\left( {\begin{array}{*{20}c} {U\left( {0,t} \right)} \\ \vdots \\ {U\left( {I - 1,t} \right)} \\ \end{array} } \right),$$

(30)

where $\mathop \sum \limits_{t} \left| {U\left( {i,t} \right)} \right|^{2} = 1$ and $B_{m,d} \left( {f,i} \right)U\left( {i,t} \right) = \left| {C_{m,d} \left( {f,i} \right)} \right|A_{s} \left( {f,i} \right)A_{a} \left( {i,t} \right)$ ($m = {\text{site}}1, \ldots ,{\text{site}}M$ and $d = x$ or $y$). If all $A\left( {i,f,t} \right)$ have the same spatial gradients (this includes the case where all events are spatially homogeneous), then the coefficients of each event in the geomagnetic field at one site are equal (i.e., $\left| {C_{m,d} \left( {f,0} \right)} \right| = \ldots = \left| {C_{m,d} \left( {f,I - 1} \right)} \right| = \left| {C_{m,d} \left( f \right)} \right|$). In other words, the coefficient $C$ depends on only the site, direction, and frequency, and does not depend on the events. As a result, Eq. 30 can be represented as

$$\begin{aligned} & \left( {\begin{array}{*{20}c} {\left| {H_{{{\text{site}}1,x}} \left( {f,t} \right)} \right|} \\ \vdots \\ {\left| {H_{{{\text{site}}M,y}} \left( {f,t} \right)} \right|} \\ \end{array} } \right) \\ & \quad = \left( {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\left| {C_{{{\text{site}}1,x}} \left( f \right)} \right|} \\ \vdots \\ {\left| {C_{{{\text{site}}M,y}} \left( f \right)} \right|} \\ \end{array} } & \ldots & {\begin{array}{*{20}c} {\left| {C_{{{\text{site}}1,x}} \left( f \right)} \right|} \\ \vdots \\ {\left| {C_{{{\text{site}}M,y}} \left( f \right)} \right|} \\ \end{array} } \\ \end{array} } \right)\left( {\begin{array}{*{20}c} {A_{s} \left( {f,0} \right)A_{a} \left( {0,t} \right)} \\ \vdots \\ {A_{s} \left( {f,I - 1} \right)A_{a} \left( {I - 1,t} \right)} \\ \end{array} } \right) \\ & \quad = \left( {\begin{array}{*{20}c} {\begin{array}{*{20}c} {B_{{{\text{site}}1,x}} \left( {f,0} \right)} \\ \vdots \\ {B_{{{\text{site}}M,y}} \left( {f,0} \right)} \\ \end{array} } & \ldots & {\begin{array}{*{20}c} {B_{{{\text{site}}1,x}} \left( {f,I - 1} \right)} \\ \vdots \\ {B_{{{\text{site}}M,y}} \left( {f,I - 1} \right)} \\ \end{array} } \\ \end{array} } \right)\left( {\begin{array}{*{20}c} {U\left( {0,t} \right)} \\ \vdots \\ {U\left( {I - 1,t} \right)} \\ \end{array} } \right). \\ \end{aligned}$$

(31)

Here, the relationship $B_{m,d} \left( {f,i} \right)U\left( {i,t} \right) = \left| {C_{m,d} \left( {f,i} \right)} \right|A_{s} \left( {f,i} \right)A_{a} \left( {i,t} \right)$ gives

$$\begin{aligned} & \left| {B_{m,d} \left( {f,i} \right)} \right|^{2} = \mathop \sum \limits_{t} \left| {B_{m,d} \left( {f,i} \right)U\left( {i,t} \right)} \right|^{2} = \left| {C_{m,d} \left( f \right)} \right|^{2} \left| {A_{s} \left( {f,i} \right)} \right|^{2} \mathop \sum \limits_{t} \left| {A_{a} \left( {i,t} \right)} \right|^{2} , \\ & {\text{or}} \\ & B_{m,d} \left( {f,i} \right) = \left| {C_{m,d} \left( f \right)} \right|A_{s} \left( {f,i} \right)\sqrt {\mathop \sum \limits_{t} \left| {A_{a} \left( {i,t} \right)} \right|^{2} } . \\ \end{aligned}$$

(32)

Dividing the Basis vector $B_{{{\text{site}}1,d}} \left( {f,i} \right)$ by the summation of Basis vectors included in the same spectrogram, we obtain the relationship:

$$\begin{aligned} & \frac{{B_{{{\text{site}}1,d}} \left( {f,i} \right)}}{{\mathop \sum \nolimits_{l} B_{{{\text{site}}1,d}} \left( {f,l} \right)}} = \frac{{\left| {C_{{{\text{site}}1,d}} \left( f \right)} \right|A_{s} \left( {f,i} \right)\sqrt {\mathop \sum \nolimits_{t} \left| {A_{a} \left( {i,t} \right)} \right|^{2} } }}{{\mathop \sum \nolimits_{l} \left| {C_{{{\text{site}}1,d}} \left( f \right)} \right|A_{s} \left( {f,l} \right)\sqrt {\mathop \sum \nolimits_{t} \left| {A_{a} \left( {l,t} \right)} \right|^{2} } }} \\ & \quad = \frac{{A_{s} \left( {f,i} \right)\sqrt {\mathop \sum \nolimits_{t} \left| {A_{a} \left( {i,t} \right)} \right|^{2} } }}{{\mathop \sum \nolimits_{l} A_{s} \left( {f,l} \right)\sqrt {\mathop \sum \nolimits_{t} \left| {A_{a} \left( {l,t} \right)} \right|^{2} } }} = \frac{{\left| {C_{{{\text{site}}M,d}} \left( f \right)} \right|A_{s} \left( {f,i} \right)\sqrt {\mathop \sum \nolimits_{t} \left| {A_{a} \left( {i,t} \right)} \right|^{2} } }}{{\left| {C_{{{\text{site}}M,d}} \left( f \right)} \right|\mathop \sum \nolimits_{l} A_{s} \left( {f,l} \right)\sqrt {\mathop \sum \nolimits_{t} \left| {A_{a} \left( {l,t} \right)} \right|^{2} } }} = \frac{{B_{{{\text{site}}M,d}} \left( {f,i} \right)}}{{\mathop \sum \nolimits_{l} B_{{{\text{site}}M,d}} \left( {f,l} \right)}}. \\ \end{aligned}$$

(33)

Therefore, if all geomagnetic events have the same spatial gradients, the equality $\frac{{B_{{{\text{site}}1,x}} \left( {f,i} \right)}}{{\mathop \sum \nolimits_{l} B_{{{\text{site}}1,x}} \left( {f,l} \right)}} = \ldots = \frac{{B_{{{\text{site}}M,y}} \left( {f,i} \right)}}{{\mathop \sum \nolimits_{l} B_{{{\text{site}}M,y}} \left( {f,l} \right)}}$ is established. Based on this contraposition, when this equation does not hold, the geomagnetic data must include anomalous geomagnetic events. If geomagnetic event $\varvec{A}\left( i \right)$ has a spatial gradient different from that of other events, the difference between $\frac{{B_{{m_{1} ,d}} \left( {f,i} \right)}}{{\mathop \sum \nolimits_{l} B_{{m_{1} ,d}} \left( {f,l} \right)}}$ and $\frac{{B_{{m_{2} ,d}} \left( {f,i} \right)}}{{\mathop \sum \nolimits_{l} B_{{m_{2} ,d}} \left( {f,l} \right)}} \left( {m_{1} \ne m_{2} } \right)$ will become large. This is because the denominator is the summation of Basis vectors, and is less affected by one event, and the numerator is the Basis vector corresponding to $\varvec{A}\left( i \right)$.

Appendix 3

Under the same condition as the examples in Table 3, Eq. 16 can be represented as follows:

$$\begin{aligned} T_{xx} \left( {1,2} \right) & \cong \frac{{C_{{{\text{site}}1,x}} \left( 0 \right)C_{{{\text{site}}2,x}}^{*} \left( 0 \right) + C_{{{\text{site}}1,x}} \left( 1 \right)C_{{{\text{site}}2,x}}^{*} \left( 1 \right)}}{{C_{{{\text{site}}2,x}} \left( 0 \right)C_{{{\text{site}}2,x}}^{*} \left( 0 \right) + C_{{{\text{site}}2,x}} \left( 1 \right)C_{{{\text{site}}2,x}}^{*} \left( 1 \right)}} \\ & = \frac{{C_{{{\text{site}}1,x}} \left( 0 \right)C_{{{\text{site}}2,x}}^{*} \left( 0 \right)}}{{C_{{{\text{site}}2,x}} \left( 0 \right)C_{{{\text{site}}2,x}}^{*} \left( 0 \right)}} \cdot \frac{{1 + \frac{{C_{{{\text{site}}1,x}} \left( 1 \right)C_{{{\text{site}}2,x}}^{*} \left( 1 \right)}}{{C_{{{\text{site}}1,x}} \left( 0 \right)C_{{{\text{site}}2,x}}^{*} \left( 0 \right)}}}}{{1 + \frac{{C_{{{\text{site}}2,x}} \left( 1 \right)C_{{{\text{site}}2,x}}^{*} \left( 1 \right)}}{{C_{{{\text{site}}2,x}} \left( 0 \right)C_{{{\text{site}}2,x}}^{*} \left( 0 \right)}}}} \\ & = \frac{{C_{{{\text{site}}1,x}} \left( 0 \right)}}{{C_{{{\text{site}}2,x}} \left( 0 \right)}} \cdot \frac{{1 + \frac{{C_{{{\text{site}}1,x}} \left( 1 \right)C_{{{\text{site}}2,x}}^{*} \left( 1 \right)}}{{C_{{{\text{site}}1,x}} \left( 0 \right)C_{{{\text{site}}2,x}}^{*} \left( 0 \right)}}}}{{1 + \frac{{C_{{{\text{site}}2,x}} \left( 1 \right)C_{{{\text{site}}2,x}}^{*} \left( 1 \right)}}{{C_{{{\text{site}}2,x}} \left( 0 \right)C_{{{\text{site}}2,x}}^{*} \left( 0 \right)}}}} \\ & = \frac{{C_{{{\text{site}}1,x}} \left( 0 \right)}}{{C_{{{\text{site}}2,x}} \left( 0 \right)}} \cdot \frac{{1 + \frac{{C_{{{\text{site}}1,x}} \left( 1 \right)C_{{{\text{site}}2,x}}^{*} \left( 1 \right)}}{{C_{{{\text{site}}1,x}} \left( 0 \right)C_{{{\text{site}}2,x}}^{*} \left( 0 \right)}} \cdot \frac{{C_{{{\text{site}}2,x}} \left( 0 \right)C_{{{\text{site}}2,x}} \left( 1 \right)}}{{C_{{{\text{site}}2,x}} \left( 1 \right)C_{{{\text{site}}2,x}} \left( 0 \right)}}}}{{1 + \frac{{C_{{{\text{site}}2,x}} \left( 1 \right)C_{{{\text{site}}2,x}}^{*} \left( 1 \right)}}{{C_{{{\text{site}}2,x}} \left( 0 \right)C_{{{\text{site}}2,x}}^{*} \left( 0 \right)}}}} \\ & = \varOmega \left( 0 \right) \cdot \frac{{1 + \frac{\varOmega \left( 1 \right)}{\varOmega \left( 0 \right)}\left| {\varPi_{{{\text{site}}2,x}} } \right|^{2} }}{{1 + \left| {\varPi_{{{\text{site}}2,x}} } \right|^{2} }}, \\ \end{aligned}$$

(34)

where $\varOmega \left( 0 \right) = \frac{{C_{{{\text{site}}1,x}} \left( 0 \right)}}{{C_{{{\text{site}}2,x}} \left( 0 \right)}}$, $\varOmega \left( 1 \right) = \frac{{C_{{{\text{site}}1,x}} \left( 1 \right)}}{{C_{{{\text{site}}2,x}} \left( 1 \right)}}$, and $\varPi_{{{\text{site}}2,x}} = \frac{{C_{{{\text{site}}2,x}} \left( 1 \right)}}{{C_{{{\text{site}}2,x}} \left( 0 \right)}}$. Similarly, $T_{xx} \left( {2,1} \right)$ can be written as

$$\begin{aligned} T_{xx} \left( {2,1} \right) & \cong \frac{{C_{{{\text{site}}2,x}} \left( 0 \right)C_{{{\text{site}}1,x}}^{*} \left( 0 \right) + C_{{{\text{site}}2,x}} \left( 1 \right)C_{{{\text{site}}1,x}}^{*} \left( 1 \right)}}{{C_{{{\text{site}}1,x}} \left( 0 \right)C_{{{\text{site}}1,x}}^{*} \left( 0 \right) + C_{{{\text{site}}1,x}} \left( 1 \right)C_{{{\text{site}}1,x}}^{*} \left( 1 \right)}} \\ & = \frac{{C_{{{\text{site}}2,x}} \left( 0 \right)C_{{{\text{site}}1,x}}^{*} \left( 0 \right)}}{{C_{{{\text{site}}1,x}} \left( 0 \right)C_{{{\text{site}}1,x}}^{*} \left( 0 \right)}} \cdot \frac{{1 + \frac{{C_{{{\text{site}}2,x}} \left( 1 \right)C_{{{\text{site}}1,x}}^{*} \left( 1 \right)}}{{C_{{{\text{site}}2,x}} \left( 0 \right)C_{{{\text{site}}1,x}}^{*} \left( 0 \right)}}}}{{1 + \frac{{C_{{{\text{site}}1,x}} \left( 1 \right)C_{{{\text{site}}1,x}}^{*} \left( 1 \right)}}{{C_{{{\text{site}}1,x}} \left( 0 \right)C_{{{\text{site}}1,x}}^{*} \left( 0 \right)}}}} \\ & = \frac{{C_{{{\text{site}}2,x}} \left( 0 \right)}}{{C_{{{\text{site}}1,x}} \left( 0 \right)}} \cdot \frac{{1 + \frac{{C_{{{\text{site}}2,x}} \left( 1 \right)C_{{{\text{site}}1,x}}^{*} \left( 1 \right)}}{{C_{{{\text{site}}2,x}} \left( 0 \right)C_{{{\text{site}}1,x}}^{*} \left( 0 \right)}} \cdot \frac{{C_{{{\text{site}}2,x}}^{*} \left( 0 \right)C_{{{\text{site}}2,x}}^{*} \left( 1 \right)}}{{C_{{{\text{site}}2,x}}^{*} \left( 1 \right)C_{{{\text{site}}2,x}}^{*} \left( 0 \right)}}}}{{1 + \frac{{C_{{{\text{site}}1,x}} \left( 1 \right)C_{{{\text{site}}1,x}}^{*} \left( 1 \right)}}{{C_{{{\text{site}}1,x}} \left( 0 \right)C_{{{\text{site}}1,x}}^{*} \left( 0 \right)}} \cdot \frac{{C_{{{\text{site}}2,x}} \left( 0 \right)C_{{{\text{site}}2,x}}^{*} \left( 0 \right)}}{{C_{{{\text{site}}2,x}} \left( 1 \right)C_{{{\text{site}}2,x}}^{*} \left( 1 \right)}} \cdot \frac{{C_{{{\text{site}}2,x}} \left( 1 \right)C_{{{\text{site}}2,x}}^{*} \left( 1 \right)}}{{C_{{{\text{site}}2,x}} \left( 0 \right)C_{{{\text{site}}2,x}}^{*} \left( 0 \right)}}}}. \\ & = \frac{1}{\varOmega \left( 0 \right)} \cdot \frac{{1 + \left( {\frac{\varOmega \left( 1 \right)}{\varOmega \left( 0 \right)}} \right)^{*} \left| {\varPi_{{{\text{site}}2,x}} } \right|^{2} }}{{1 + \left| {\frac{\varOmega \left( 1 \right)}{\varOmega \left( 0 \right)}} \right|^{2} \left| {\varPi_{{{\text{site}}2,x}} } \right|^{2} }}. \\ \end{aligned}$$

(35)

$T_{xx} \left( {1,2} \right)$ and $T_{xx} \left( {2,1} \right)$ can be represented by two complex values ($\varOmega \left( 0 \right)$ and $\varOmega \left( 1 \right)$) and one real value ($\left| {\varPi_{{{\text{site}}2,x}} } \right|$). Note that, given the other two values, one of these three values is not uniquely determined. For example, given $\varOmega \left( 0 \right) = 1$ and $\varOmega \left( 1 \right) = 1$, $\varPi_{{{\text{site}}2,x}}$ can take a value of 1, 2, or 3. As in the simulation described in Table 3, we consider the simple case in which all $C$ in Eqs. 34 and 35 are real-positive numbers. Differentiating $T_{xx} \left( {1,2} \right)$ and $T_{xx} \left( {2,1} \right)$ with respect to $\varOmega \left( 0 \right)$ gives

$$\frac{{\partial T_{xx} \left( {1,2} \right)}}{\partial \varOmega \left( 0 \right)} \cong \frac{\partial }{\partial \varOmega \left( 0 \right)}\left\{ {\varOmega \left( 0 \right) \cdot \frac{{1 + \frac{\varOmega \left( 1 \right)}{\varOmega \left( 0 \right)}\left| {\varPi_{{{\text{site}}2,x}} } \right|^{2} }}{{1 + \left| {\varPi_{{{\text{site}}2,x}} } \right|^{2} }}} \right\} = \frac{1}{{1 + \left| {\varPi_{{{\text{site}}2,x}} } \right|^{2} }},$$

(36)

and

$$\begin{aligned} \frac{{\partial T_{xx} \left( {2,1} \right)}}{\partial \varOmega \left( 0 \right)} & \cong \frac{\partial }{\partial \varOmega \left( 0 \right)}\left\{ {\frac{{\varOmega \left( 0 \right) + \varOmega \left( 1 \right)\left| {\varPi_{{{\text{site}}2,x}} } \right|^{2} }}{{\varOmega \left( 0 \right)^{2} + \varOmega \left( 1 \right)^{2} \left| {\varPi_{{{\text{site}}2,x}} } \right|^{2} }}} \right\} \\ & = \frac{{\left\{ {\varOmega \left( 0 \right)^{2} + \varOmega \left( 1 \right)^{2} \left| {\varPi_{{{\text{site}}2,x}} } \right|^{2} } \right\} - 2\varOmega \left( 0 \right)\left\{ {\varOmega \left( 0 \right) + \varOmega \left( 1 \right)\left| {\varPi_{{{\text{site}}2,x}} } \right|^{2} } \right\}}}{{\left\{ {\varOmega \left( 0 \right)^{2} + \varOmega \left( 1 \right)^{2} \left| {\varPi_{{{\text{site}}2,x}} } \right|^{2} } \right\}^{2} }} \\ & = \frac{{\left\{ { - \varOmega \left( 0 \right)^{2} - 2\varOmega \left( 0 \right)\varOmega \left( 1 \right)\left| {\varPi_{{{\text{site}}2,x}} } \right|^{2} + \varOmega \left( 1 \right)^{2} \left| {\varPi_{{{\text{site}}2,x}} } \right|^{2} } \right\}}}{{\left\{ {\varOmega \left( 0 \right)^{2} + \varOmega \left( 1 \right)^{2} \left| {\varPi_{{{\text{site}}2,x}} } \right|^{2} } \right\}^{2} }}. \\ \end{aligned}$$

(37)

Equation 36 shows that $T_{xx} \left( {1,2} \right)$ is a monotonically increasing function of $\varOmega \left( 0 \right)$. However, Eq. 37 shows that $T_{xx} \left( {2,1} \right)$ is not a monotonically decreasing function of $\varOmega \left( 0 \right)$. When $\varOmega \left( 0 \right)$ is less than $- \varOmega \left( 1 \right)\left| {\varPi_{{{\text{site}}2,x}} } \right|^{2} + \varOmega \left( 1 \right)\left| {\varPi_{{{\text{site}}2,x}} } \right|\sqrt {1 + \left| {\varPi_{{{\text{site}}2,x}} } \right|^{2} }$, $\frac{{\partial T_{xx} \left( {2,1} \right)}}{\partial \varOmega \left( 0 \right)}$ is greater than 0. Although the numerical examples in Table 3 are not calculated under the condition that only one of $\varOmega \left( 0 \right)$, $\varOmega \left( 1 \right)$, and $\varPi_{{{\text{site}}2,x}}$ is varying, we can prove that $T_{xx} \left( {1,2} \right)$ and $T_{xx} \left( {2,1} \right)$ do not always shift with the opposite polarity.