Maximum-variance gradiometer technique for removal of spacecraft-generated disturbances from magnetic field data

In situ measurement of the magnetic field using spaceborne instruments requires a magnetically clean platform and/or a very long boom for accommodating magnetometer sensors at a large distance from the spacecraft body. This significantly drives up the costs and the time required to build a spacecraft. Here we present an alternative sensor configuration and a technique allowing for removal of the spacecraft-generated AC disturbances from the magnetic field measurements, thus lessening the need for a magnetic cleanliness programme and allowing for shorter boom length. The final expression of the corrected data takes the form of a linear combination of the measurements from all sensors, allowing for simple onboard software implementation. The proposed technique is applied to the Service Oriented Spacecraft Magnetometer (SOSMAG) on board the Korean geostationary satellite GeoKompsat-2A (GK2A). In contrast to other missions where multi-sensor measurements were used to clean the data on the ground, the SOSMAG instrument performs the cleaning on board and transmits the corrected data in real time, as needed by space weather applications. The successful elimination of the AC disturbances originating from several sources validates the proposed cleaning technique.

1 Introduction

Since very early in space exploration it has been clear that the main limitation in performing accurate magnetic field measurements comes not from the instruments themselves but rather from the strong artificial magnetic fields generated by the spacecraft carrying them. It was recognized that there are three possible approaches to mitigating this problem: one could limit the electromagnetic emissions coming from the spacecraft by going through a rigorous magnetic cleaning procedure. This is a costly and complicated engineering task and introduces limitations on building and operating other onboard instruments; see, for example, Narvaez (2004) for details on the magnetic cleanliness programme for the Cassini magnetic field experiment (Dougherty et al., 2004). Another approach is to accommodate the magnetometer at a large distance from the spacecraft, usually at the end of a long boom, such as the 12 m long Kaguya boom (Kato et al., 2010) or the 13 m long Voyager boom (Behannon et al., 1977). This introduces constrains on the spacecraft operations and still requires a certain degree of magnetic cleanliness of the spacecraft in order to keep the boom at a reasonable length. A third way is to accept the presence of spacecraft-generated disturbances in the measured magnetic field and to remove the artificial contributions afterwards on board or on the ground through special techniques (Mehlem, 1978; Georgescu et al., 2008; Pope et al., 2011). An extreme case where neither magnetic cleanliness nor boom was provided is, for example, the magnetic field experiment on the MASCOT lander (Herčík et al., 2017). In most cases, however, a combination of two or all of the approaches above is employed. For instance, Cluster (Escoubet et al., 1997) and THEMIS (Angelopoulos, 2008) are magnetically clean spacecraft carrying magnetometers on relatively long booms. For normal science investigations, the stray magnetic field from these spacecraft is well below the required accuracy, and no further steps to remove it are usually necessary. Venus Express (Titov et al., 2006), on the other hand, was a magnetically dirty spacecraft with two magnetometers (Zhang et al., 2006) mounted on a short boom for which extensive data cleaning efforts had to be undertaken (Pope et al., 2011). A comprehensive overview of the instrumentation and challenges related to measuring magnetic fields in space is given by Balogh (2010). In this work we focus on the third approach: removal of the contribution of the spacecraft-generated magnetic field from the measured data, without the need of extensive information on potential spacecraft disturbance sources.

One of the first studies on using multi-sensor measurements to clean magnetic field data measured on board spacecraft came from Ness et al. (1971). The proposed method was then successfully applied in a simplified manner to Mariner 10 magnetic field data (Ness et al., 1974) assuming one single dipole disturber source. Neubauer (1975) gave a detailed error analysis of the Ness et al. (1971) method and discussed the optimum placement of collinear sensors. The more recent cleaning procedure used by Pope et al. (2011) for Venus Express, though based on the same principle, is much more sophisticated, allowing removal of disturbances from several different sources. However, additional information about the spacecraft operation and fuzzy logic had to be used to distinguish between the disturbance sources. Such a complex algorithm would be difficult to implement for onboard data cleaning. Our aim is a correction method which reduces to a linear (or at most quadratic) combination of the magnetic field values measured by several sensors without input from other sources and is therefore easy to implement on board.

Similarly to the Ness et al. (1971) and Pope et al. (2011) methods, the disturbance removal method described in the following sections is based on the fact that the magnetic field measured by each sensor is the sum of the ambient magnetic field, and the artificial magnetic field generated by the spacecraft. Because the ambient field is the same for all sensors, it vanishes in the difference between the measurements from any two sensors, similar to the gradiometer working principle. The difference is entirely determined by the artificial magnetic field sources from the spacecraft, preserving their time dependence. Magnetic disturbances generated by time-dependent currents flowing through simple mechanically fixed current loops keep constant direction; therefore in general the disturbance affects only one component of the measured field. If the variation of the disturbing magnetic field is much larger than the variation of the ambient magnetic field during the time interval selected to determine the cleaning parameters, the direction of the strongest disturbance will coincide with the principal component (maximum-variance component) of the measured field, allowing application of the correction only to the affected component. This is the type of magnetic disturbances which can be treated using the method described in the next sections. If the direction of the disturbance changes in time – as is the case for instance for magnetic fields produced by flywheels or other moving mechanisms – then another approach must be used.

The proposed method is applied to the Service Oriented Spacecraft Magnetometer (SOSMAG) instrument (Auster et al., 2016; Magnes et al., 2020), which, together with the Particle Detector experiment (Seon et al., 2020), is part of the Korea Space wEather Monitor (KSEM) (Oh et al., 2018) on board the GeoKompsat-2A (GK2A) geostationary spacecraft. SOSMAG consists of four three-axial magnetic field sensors, two of them mounted on a short boom extended from the spacecraft and the other two placed near strong magnetic disturbance sources within the spacecraft. Once the correction coefficients are determined on the ground, they are uploaded to the spacecraft and are used by the onboard software to correct the magnetic field measurements in flight. This enables accurate magnetic field measurements which are delivered in near-real time to the ground stations without the need of passing through a magnetic cleanliness programme before launch. The quick data delivery is essential in the context of space weather monitoring.

The remainder of the paper is organized as follows: in Sect. 2 we discuss the gradiometer principle on which our method is based. Section 3 outlines the proposed principal-component gradiometer (PiCoG) method to remove spacecraft-generated disturbances from the measured magnetic field data. Section 4 describes how the PiCoG method is applied to clean the GK2A SOSMAG data. The limitations of the proposed method are discussed in Sect. 5. Section 6 summarizes our work.

2 Disturbances from known sources

This section gives the analytical expressions for disturbances when the exact locations of the magnetic field sources and of the sensors are known. While in most cases the direct application of these expressions is not practical, this section outlines the general principle used by gradiometer-based disturbance cleaning methods, namely the possibility to express the spacecraft-generated disturbances in terms of differences between measurements taken at distinct places. The relations derived here constitute the basis of the PiCoG technique detailed in Sect. 3. They are valid for both AC and DC disturbances, though the PiCoG technique only deals with AC disturbances. Because higher multipole moments attenuate strongly with the distance to the source and become negligible even for short booms, we will concentrate only on the dipole and quadrupole contributions.

2.1 Single disturbance source

The magnetic field produced at the position $\mathit{r}=r\widehat{\mathit{r}}$ by a dipole characterized by a slowly varying, time-dependent magnetic moment M(t) is given by

$$\begin{array}{}\text{(1)}& \mathit{b}(\mathit{r},t)={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }{r}^{\mathrm{3}}}}\left(\mathrm{3}\mathbf{X}\right(\widehat{\mathit{r}})-\mathbf{I})\mathit{M}\left(t\right),\end{array}$$

where the elements $X_{kl}={\widehat{r}}_k{\widehat{r}}_l$ of the matrix X are given by the product between the components of the position versor $\widehat{\mathit{r}}$, and I is the 3×3 identity matrix. The subscripts k and l refer to Cartesian components. Knowing the magnetic field at the position rⁱ, one can compute the magnetic field at any position r^j without knowledge about the source magnetic moment M(t):

$$\begin{array}{}\text{(2)}& \mathit{b}\left({\mathit{r}}^j,t\right)={\mathbf{T}}^{\text{dip}}\left({\mathit{r}}^i,{\mathit{r}}^j\right)\mathit{b}\left({\mathit{r}}^i,t\right),\end{array}$$

where the superscripts i and j denote the measurement positions, and the time-independent linear transformation T^dip is

$$\begin{array}{}\text{(3)}& {\mathbf{T}}^{\text{dip}}\left({\mathit{r}}^i,{\mathit{r}}^j\right)={\left({\displaystyle \frac{r^i}{r^j}}\right)}^{\mathrm{3}}\left(\mathrm{3}{\mathbf{X}}^j-\mathbf{I}\right){\left(\mathrm{3}{\mathbf{X}}^i-\mathbf{I}\right)}^{-\mathrm{1}}.\end{array}$$

The inverse ${\left(\mathrm{3}\mathbf{X}-\mathbf{I}\right)}^{-\mathrm{1}}$ always exists and is equal to $\left(\mathrm{3}\mathbf{X}/\mathrm{2}-\mathbf{I}\right)$. To derive this, we used the fact that X is an idempotent matrix, X²=X.

Assuming that the ambient magnetic field is generated by distant sources and thus it is the same at the positions rⁱ and r^j, it is possible to separate the contribution b(rⁱ,t) due to a nearby dipole from the ambient field by computing the difference between the measured magnetic field at the two positions

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\mathit{b}\left({\mathit{r}}^i,t\right)=& {\left({\mathbf{T}}^{\text{dip}}\left({\mathit{r}}^i,{\mathit{r}}^j\right)-\mathbf{I}\right)}^{-\mathrm{1}}\\ & \cdot \left({\mathit{B}}_{\text{measured}}\left({\mathit{r}}^j,t\right)-{\mathit{B}}_{\text{measured}}\left({\mathit{r}}^i,t\right)\right),\end{array}\end{array}$$

where the total measured magnetic field ${\mathit{B}}_{\text{measured}}(\mathit{r},t)=\mathit{B}\left(t\right)+\mathit{b}(\mathit{r},t)$ contains both the position-independent ambient magnetic field B(t) and the position-dependent disturbance magnetic field b(r,t). Sensor-specific disturbances such as sensor noise and sensor offset will be considered later.

Note that the T^dip matrix only depends on the position vectors rⁱ and r^j. It is independent of the dipole M(t) and performs a similar function to the propagator operator in quantum mechanics. Equation (4) shows that, once the T^dip matrix is determined for a pair of sensors, measurements from those two points are sufficient to separate the contribution from a single magnetic field source with arbitrary time variation from the ambient magnetic field. This is the theoretical justification for our method. Also note that for Eq. (4) to be satisfied it is not strictly necessary that the disturbance has a dipole character. It is enough that a (time-independent) linear relation exists between the disturbing magnetic field affecting the sensors at the positions rⁱ and r^j.

Similar relations can be written for a time-dependent quadrupole defined by its moment Q(t):

$$\begin{array}{}\text{(5)}& {\displaystyle }\mathit{b}\left(\mathit{r},t\right)& {\displaystyle }={\displaystyle \frac{{\mathit{\mu }}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }{r}^{\mathrm{4}}}}\left(\mathrm{5}\mathbf{X}\right(\widehat{\mathit{r}})-\mathrm{2}\mathbf{I})\mathbf{Q}\left(t\right)\widehat{\mathit{r}},\text{(6)}& {\displaystyle }\mathit{b}\left({\mathit{r}}^j,t\right)& {\displaystyle }={\mathbf{T}}^{\text{quad}}\left({\mathit{r}}^i,{\mathit{r}}^j\right)\mathit{b}\left({\mathit{r}}^i,t\right),\text{(7)}& {\displaystyle }{\mathbf{T}}^{\text{quad}}\left({\mathit{r}}^i,{\mathit{r}}^j\right)& {\displaystyle }={\left({\displaystyle \frac{r^i}{r^j}}\right)}^{\mathrm{4}}\left(\mathrm{5}{\mathbf{X}}^j-\mathrm{2}\mathbf{I}\right){\mathbf{G}}^{ji}{\left(\mathrm{5}{\mathbf{X}}^i-\mathrm{2}\mathbf{I}\right)}^{-\mathrm{1}},\end{array}$$

where ${\mathbf{G}}^{ji}={\mathbf{QR}}^{ij}{\mathbf{Q}}^{-\mathrm{1}}$. R^ij is the rotation matrix which transforms the versor ${\widehat{\mathit{r}}}^i$ to the versor ${\widehat{\mathit{r}}}^j$, and ${\left(\mathrm{5}\mathbf{X}-\mathrm{2}\mathbf{I}\right)}^{-\mathrm{1}}$ is equal to $\left(\mathrm{5}\mathbf{X}/\mathrm{6}-\mathbf{I}/\mathrm{2}\right)$. While there are instances when T^quad is independent of the quadrupole moment Q(t) (e.g. when the quadrupole source and the two sensors are aligned), in general T^quad depends on it.

For our purposes, however, it is important that T^quad does not depend of time. A common situation in which this happens is when only the magnitude of the quadrupolar disturbance depends on time. Then the time dependence can be separated as an independent scalar factor in the expression of the quadrupole moment, Q(t)=q(t)Q₀, and therefore T^quad becomes time independent. A relation similar to Eq. (4) can then also be written for the quadrupole, and the disturbance can be removed using the same procedure as for a dipole disturbance. In what follows we assume this kind of time variation for quadrupole sources.

2.2 Multiple disturbance sources

The contributions from more than one simultaneously active, arbitrary placed source with arbitrary time dependence cannot be separated from the ambient field in a simple way. However, if multiple sensors are arranged in a suitable configuration and if specific properties of the disturbers, such as known polarization or time dependence, are used, it is possible to remove disturbances generated by multiple sources.

Two magnetometers represent the minimal configuration needed to eliminate stray spacecraft magnetic fields. Many spacecraft carry two magnetometers attached at different positions along one boom. If the boom is long enough such that the distances between the disturbance sources are much smaller compared to the distances to the measurement points and if the disturbances all have either pure dipole or quadrupole character, then their T matrices will be the same and their collective disturbance can be separated from the ambient field in one step using only two sensors as was done, for example, by Ness et al. (1974). Of course, a collection of dipoles will in general produce multipole moments. For the procedure to work, the quadrupole and higher-order contributions must be much weaker than the dipole contribution at both sensors. If, however, both dipole and quadrupole contributions are present at the same time with comparable strengths, then their T matrices will differ due to the different attenuation with the distance. In this case, one must rely on specific properties of the disturbance to eliminate the quadrupole contribution.

In contrast to the minimum two-magnetometer configuration, one can imagine a configuration such that for each disturber there is a sensor placed much closer to it than to all other disturbers, plus an additional sensor far away from all disturbers. Then for each sensor the far disturbers can be assimilated to the ambient field and the problem becomes the single-disturber problem discussed at the beginning of this section. Each contribution can then be separated from the ambient field independently. Such a sensor configuration is ideal and can be attained with a number of sensors placed within the spacecraft plus one sensor placed on a short boom.

If the disturbing magnetic field has a time-dependent magnitude but does not change its direction, i.e. if its variation is linearly polarized, then up to three independent, simultaneously active disturbances with mutually orthogonal variance directions can be separated using two sensors. This is done by projecting Eq. (4) on the direction of each disturbance. The direction of each disturbance can be determined using principal-component analysis as described in Sects. 3 and 4. This kind of linearly polarized disturbances produced by fixed-configuration time-dependent currents is commonly encountered. The PiCoG cleaning method assumes this type of linearly polarized disturbances. If more than three disturbances, or disturbances with their polarization directions not mutually orthogonal, are present, then information from more sensors is necessary. Different sensor pairs will correct different disturbances.

The SOSMAG configuration on board the GK2A spacecraft lies somewhere in between the ideal configuration above and the minimum two-magnetometer configuration. It consists of two high-accuracy magnetometers placed on a relatively short boom and a number of resource-saving magnetometers placed inside the spacecraft. As we will show in Sect. 4, this configuration is well suited to apply the PiCoG cleaning method.

3 The principal-component gradiometer technique

The PiCoG cleaning technique is based on the fact that, while the ambient magnetic field does not change over the spacecraft scale, the magnitude of a spacecraft-generated disturbance in the magnetic field decreases with the distance to the disturbance source. Therefore, the disturbance can be detected – and subsequently removed from the useful signal – by comparing measurements from sensors placed at different distances to the disturbance source as outlined in Sect. 2.

If the precise positions of the disturbers and of the sensors are known, then the transformation matrices T, which allow the separation of disturbances generated by the spacecraft, can be computed directly. The measurements are then cleaned using Eq. (4) and the equivalent equation for quadrupole disturbers. However, this is in general not the case. Here we describe the derivation of the T matrices under certain assumptions but without prior knowledge about the exact positions of the disturbance sources.

The magnetic field measured by the sensor i can be written as the sum of the ambient magnetic field, B(t); the sum of the disturbances ${\mathit{b}}^q(t,{\mathit{r}}^{iq})={\mathit{b}}^{qi}\left(t\right)$ created by N sources placed at relative positions ${\mathit{r}}^{iq}={\mathit{r}}^i-{\mathit{r}}^q$ from the sensor i; and a term containing the sensor-specific disturbance (noise and time-dependent offset), Zⁱ(t):

$$\begin{array}{}\text{(8)}& {\mathit{B}}^{\mathrm{0},i}\left(t\right)=\mathit{B}\left(t\right)+\sum _{q=\mathrm{1}}^N{\mathit{b}}^{qi}\left(t\right)+{\mathit{Z}}^i\left(t\right),\end{array}$$

where the index 0 on the left side indicates the initially measured magnetic field.

We can eliminate the ambient field by subtracting the measurements from two sensors placed at distinct positions:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}\mathrm{\Delta }{\mathit{B}}^{\mathrm{0},ij}\left(t\right)& ={\mathit{B}}^{\mathrm{0},i}\left(t\right)-{\mathit{B}}^{\mathrm{0},j}\left(t\right)\\ & =\sum _{q=\mathrm{1}}^N\mathrm{\Delta }{\mathit{b}}^{qij}\left(t\right)+\mathrm{\Delta }{\mathit{Z}}^{ij}\left(t\right).\end{array}\end{array}$$

If we neglect the sensor-specific disturbances, for single disturbers, the correction to be applied to the measurements consists of a linear combination of the components of the difference ΔB^0,ij(t) between the measured magnetic field at each sensor position:

$$\begin{array}{}\text{(10)}& {\mathit{B}}_{\text{corrected}}^i\left(t\right)={\mathit{B}}^{\mathrm{0},i}\left(t\right)+{\mathbf{A}}^{ij}\mathrm{\Delta }{\mathit{B}}^{\mathrm{0},ij}\left(t\right).\end{array}$$

Using Eq. (4) we find that the matrix A^ij is the equal to $-\left(\mathbf{T}\right({\mathit{r}}^i,{\mathit{r}}^j)-\mathbf{I}{)}^{-\mathrm{1}}$. For each sensor pair i,j a matrix A^ij must be determined. This may of course be computed if we know the exact coordinates of the sensors and of the disturbers and if the disturbers are pure single dipoles or quadrupoles. This is in general not true; therefore we will derive the correction matrix A^ij directly from the measurements.

3.1 First-order correction

We now assume that one of the terms in Eq. (9) is much larger than the others. This is true if one of the disturbance sources is much stronger or much closer to one of the sensors than to the others. In this case the corrected measurements are given by Eq. (10). Note that the small distance between the disturbance source and a sensor does not imply significant contribution from higher-order multipoles. It merely implies that the dominant dipole/quadrupole term produced by the source in question at the sensor location is much larger than the contribution from the other disturbance sources. Also note that, even though the positions of the disturbers do not enter the PiCoG formalism, some rough information about the positions of major disturbers can help in optimizing the accommodation of the sensors by placing them near major disturbers.

Equations (8) and (9) reduced to the single-disturber form are

$$\begin{array}{}\text{(11)}& {\displaystyle }{\mathit{B}}^{\mathrm{0},i}\left(t\right)=\mathit{B}\left(t\right)+{\mathit{b}}^i\left(t\right)+{\mathit{Z}}^i\left(t\right),\text{(12)}& {\displaystyle }\mathrm{\Delta }{\mathit{B}}^{\mathrm{0},ij}\left(t\right)=\mathrm{\Delta }{\mathit{b}}^{ij}\left(t\right)+\mathrm{\Delta }{\mathit{Z}}^{ij}\left(t\right),\end{array}$$

where we drop the disturbance source index, q.

For many spacecraft, including GK2A, many artificial disturbances are produced by simple fixed-geometry currents without phase delays, and thus their magnetic moments are fixed in direction with only their modules changing in time (M(t)=m(t)M₀ and/or Q(t)=q(t)Q₀). Therefore, in the proper coordinate system, only one component of the measured field is affected by one disturbance source. This is a key condition for applying the PiCoG technique.

To find the direction of the disturbance at the sensor positions, we need to assume that the variance due to the disturbance at the sensor positions determines the maximum-variance direction of the measured magnetic field. This holds either when the variance of the disturbance is much larger than the variance of the ambient field or when the variance of the ambient field does not have a preferred direction. In this case, the direction of the disturbance at both sensors can be estimated through variance analysis (Sonnerup and Scheible, 1998; Song and Russell, 1999) of the 3D time series from each sensor. The principal components at each sensor are then the magnetic field components along the maximum-variance directions. The variance is used as a measure of how strong the AC disturbance is in each direction. The maximum-variance direction identifies the strongest component (the only component in the case of linear polarization) of the disturbance for both regular and random disturbances.

Our strategy is first to isolate the disturbance as the maximum-variance component of the differences ΔB^0,ij and then to use it to correct only the maximum-variance component of the measurements B^0,i. Since in general the direction of the disturbance varies from sensor to sensor, different reference systems must be used for different sensors and for the measurements differences.

The components of the magnetic field at the sensor i, corrected using measurements from the sensor j, can be written in the variance principal system (VPS) of the sensor i measurements as

$$\begin{array}{}\text{(13a)}& {\displaystyle }{B}_x^{\mathrm{1},ij}& {\displaystyle }=B_x^{\mathrm{0},i}-{\mathit{\alpha }}^{\mathrm{0},ij}(\mathrm{\Delta }{\mathit{B}}^{\mathrm{0},ij}{)}_x,\text{(13b)}& {\displaystyle }{B}_y^{\mathrm{1},ij}& {\displaystyle }=B_y^{\mathrm{0},i},\text{(13c)}& {\displaystyle }{B}_z^{\mathrm{1},ij}& {\displaystyle }=B_z^{\mathrm{0},i}.\end{array}$$

The superscript “1” in Eq. (13) stands for the first-order correction. Note that, while the left-hand sides and the first term of the right-hand sides of Eq. (13) are represented in the VPS of the measurements at the sensor i, (ΔB^0,ij)_x on the right-hand side of Eq. (13a) is represented in the VPS of the difference ΔB^0,ij. The VPS has the x axis aligned with the maximum variance and the z axis aligned with the minimum variance. Equation (13a) reflects the fact that the correction for the maximum-variance component of the measured magnetic field B^0,i is proportional to the maximum-variance component of the difference ΔB^0,ij. The other two components of the measured magnetic field remain unaffected by the correction.

Since the difference (ΔB^0,ij)_x is proportional to the disturbance to be cleaned, the scaling factor α^0,ij in Eq. (13a) is the ratio between the amplitude of the difference and the amplitude of the disturbance at the position of the sensor i. Assuming that most of the variance of the magnetic field measured by the sensor i is due to the disturbance to be cleaned, a first estimation of the α^0,ij factor is given by the variance of the measurements:

$$\begin{array}{}\text{(14)}& {\mathit{\alpha }}^{\mathrm{0},ij}=\pm \sqrt{{\displaystyle \frac{\mathrm{Var}\left({\left({\mathit{B}}^{\mathrm{0},i}\right)}_x\right)}{\mathrm{Var}\left({\left(\mathrm{\Delta }{\mathit{B}}^{\mathrm{0},ij}\right)}_x\right)}}}\end{array}$$

The ± sign above is due to the fact that, while the orientation of the x axis of the VPS is determined from variance analysis, its sense remains arbitrary. If necessary, the scaling factor value computed using Eq. (14) can be refined, for example, by minimizing the correlation between the corrected magnetic field $B_x^{\mathrm{1},ij}$ and the difference (ΔB^0,ij)_x.

If R^0,i is the rotation matrix from the sensor system to the VPS of the measurements from the sensor i, and R^0,ij is the rotation matrix from the sensor system to the VPS of the difference ΔB^0,ij, then in the sensor system Eq. (13) takes the form

$$\begin{array}{}\text{(15)}& B_k^{\mathrm{1},ij}=B_k^{\mathrm{0},i}-{\mathit{\alpha }}^{\mathrm{0},ij}\left({\left({\mathbf{R}}^{\mathrm{0},i}\right)}^{-\mathrm{1}}{)}_{kx}\right({\mathbf{R}}^{\mathrm{0},ij}\mathrm{\Delta }{\mathit{B}}^{\mathrm{0},ij}{)}_x\end{array}$$

for $k=\mathrm{1},\mathrm{\dots },\mathrm{3}$.

In matrix form the above relation can be written as

$$\begin{array}{}\text{(16)}& {\mathit{B}}^{\mathrm{1},ij}={\mathit{B}}^{\mathrm{0},i}+{\mathbf{A}}^{\mathrm{0},ij}\mathrm{\Delta }{\mathit{B}}^{\mathrm{0},ij},\end{array}$$

where the matrix A with elements

$$\begin{array}{}\text{(17)}& A_{kl}^{\mathrm{0},ij}=-{\mathit{\alpha }}^{\mathrm{0},ij}\left({\left({\mathbf{R}}^{\mathrm{0},i}\right)}^{-\mathrm{1}}{)}_{kx}\right({\mathbf{R}}^{\mathrm{0},ij}{)}_{xl}\end{array}$$

is the correction matrix for the first (strongest) disturber. Note that there is no implicit summation over repeating indices.

3.1.1 The collinear case

While not required, the special case in which the disturbance source is collinear with the two sensors is instructive. In this case, the direction of a linearly polarized disturbance will be the same at both sensors, therefore the same coordinate system will be used for Eq. (13). Substituting $B_x^{\mathrm{0},i}$ in Eq. (13a) using Eq. (11), we obtain

$$\begin{array}{}\text{(18a)}& \begin{array}{rl}{B}_x^{\mathrm{1},ij}=B_x& +\left(a-{\mathit{\alpha }}^{\mathrm{0},ij}(a-\mathrm{1})\right)b_x^j\\ & +Z_x^i-{\mathit{\alpha }}^{\mathrm{0},ij}(Z_x^i-Z_x^j),\end{array}\end{array}$$

$$\begin{array}{}\text{(18b)}& \begin{array}{rl}{B}_x^{\mathrm{1},ji}=B_x& +\left(\mathrm{1}+{\mathit{\alpha }}^{\mathrm{0},ji}(a-\mathrm{1})\right)b_x^j\\ & +Z_x^j-{\mathit{\alpha }}^{\mathrm{0},ji}\left(Z_x^j-Z_x^i\right),\end{array}\end{array}$$

where we made use of the proportionality between the spacecraft-generated disturbances at the sensors i and j: $b_x^i=a{b}_x^j$. For a dipolar disturber at distance rⁱ from the sensor i, and r^j from the sensor j, $a=(r^j/r^i{)}^{\mathrm{3}}$. For a quadrupolar disturber $a=(r^j/r^i{)}^{\mathrm{4}}$.

Since the corrected magnetic field should be independent of the disturbing magnetic field $b_x^j$, it results that the factors multiplying $b_x^j$ in Eq. (18) must be 0; therefore

$$\begin{array}{}\text{(19)}& {\mathit{\alpha }}^{\mathrm{0},ij}={\displaystyle \frac{a}{a-\mathrm{1}}}\text{and}{\mathit{\alpha }}^{\mathrm{0},ji}={\displaystyle \frac{-\mathrm{1}}{a-\mathrm{1}}}.\end{array}$$

This shows that in the collinear case the sum of the α coefficients is equal to 1:

$$\begin{array}{}\text{(20)}& {\mathit{\alpha }}^{\mathrm{0},ij}+{\mathit{\alpha }}^{\mathrm{0},ji}=\mathrm{1}.\end{array}$$

A consequence of the above is that the difference between the corrected measurements at the two sensors is always 0:

$$\begin{array}{}\text{(21)}& \begin{array}{rl}\mathrm{\Delta }{B}_x^{\mathrm{1},ij}& =B_x^{\mathrm{1},ij}-B_x^{\mathrm{1},ji}\\ & =\left(\mathrm{1}-{\mathit{\alpha }}^{\mathrm{0},ij}-{\mathit{\alpha }}^{\mathrm{0},ji}\right)\left(b_x^i-b_x^j+Z_x^i-Z_x^j\right)\equiv \mathrm{0}.\end{array}\end{array}$$

In other words, the corrected field is the same regardless of which sensor is used as the “primary” sensor: $B_x^{\mathrm{1},ij}=B_x^{\mathrm{1},ji}$.

For α obeying Eq. (19) the corrected field given by Eq. (18) is

$$\begin{array}{}\text{(22)}& B_x^{\mathrm{1},ij}=B_x+Z_x^i-{\mathit{\alpha }}^{\mathrm{0},ij}\left(Z_x^i-Z_x^j\right).\end{array}$$

Comparing the above with Eq. (11) shows that, apart from eliminating the spacecraft-generated disturbance $b_x^i$, the procedure introduces an additional disturbance which mixes the two sensor-specific disturbances $Z_x^i$ and $Z_x^j$ scaled by α^0,ij, potentially increasing the noise in the corrected measurements. This effect was also noted by Delva et al. (2002). However, if α^0,ij approaches unity (disturbance source much closer to sensor i), the i sensor-specific noise is replaced by the j sensor-specific noise, which might lead to reduced noise.

3.2 Higher-order corrections

Further corrections can be iteratively applied as long as the stray fields from different disturbers do not have the same direction at the magnetometer location. The iteration relation from order n−1 to order n is

$$\begin{array}{}\text{(23)}& {\mathit{B}}^{n,ij}={\mathit{B}}^{n-\mathrm{1},ij}+{\mathbf{A}}^{n-\mathrm{1},ij}\mathrm{\Delta }{\mathit{B}}^{n-\mathrm{1},ij},\end{array}$$

with ${\mathit{B}}^{\mathrm{0},ij}={\mathit{B}}^{\mathrm{0},i}$ and

$$\begin{array}{}\text{(24)}& A_{kl}^{n-\mathrm{1},ij}=-{\mathit{\alpha }}^{n-\mathrm{1},ij}\left({\left({\mathbf{R}}^{n-\mathrm{1},i}\right)}^{-\mathrm{1}}{)}_{kx}\right({\mathbf{R}}^{n-\mathrm{1},ij}{)}_{xl}.\end{array}$$

The α^n,ij coefficient is estimated from the variance of the field corrected up to order n. The rotation matrices R^n,i and R^n,ij refer to the order n corrected field.

Using Eqs. (23) and (24), we find the corrected magnetic field in the second and third order written as linear combinations of the difference of the measurements taken at the two sensors:

$$\begin{array}{}\text{(25)}& {\displaystyle }\begin{array}{rl}{\mathit{B}}^{\mathrm{2},ij}& ={\mathit{B}}^{\mathrm{0},i}+({\mathbf{A}}^{\mathrm{0},ij}+{\mathbf{A}}^{\mathrm{1},ij}\\ & +{\mathbf{A}}^{\mathrm{1},ij}\left({\mathbf{A}}^{\mathrm{0},ij}+{\mathbf{A}}^{\mathrm{0},ji}\right))\mathrm{\Delta }{\mathit{B}}^{\mathrm{0},ij},\end{array}\text{(26)}& {\displaystyle }\begin{array}{rl}{\mathit{B}}^{\mathrm{3},ij}& ={\mathit{B}}^{\mathrm{0},i}+({\mathbf{A}}^{\mathrm{0},ij}+{\mathbf{A}}^{\mathrm{1},ij}+{\mathbf{A}}^{\mathrm{2},ij}\\ & +{\mathbf{A}}^{\mathrm{1},ij}\left({\mathbf{A}}^{\mathrm{0},ij}+{\mathbf{A}}^{\mathrm{0},ji}\right)\\ & +{\mathbf{A}}^{\mathrm{2},ij}\left({\mathbf{A}}^{\mathrm{0},ij}+{\mathbf{A}}^{\mathrm{0},ji}+{\mathbf{A}}^{\mathrm{1},ij}+{\mathbf{A}}^{\mathrm{1},ji}\right)\\ & +{\mathbf{A}}^{\mathrm{2},ij}\left({\mathbf{A}}^{\mathrm{1},ij}+{\mathbf{A}}^{\mathrm{1},ji}\right)\left({\mathbf{A}}^{\mathrm{0},ij}+{\mathbf{A}}^{\mathrm{0},ji}\right))\mathrm{\Delta }{\mathit{B}}^{\mathrm{0},ij}.\end{array}\end{array}$$

The corrected field B^n,ij determined for the sensor i can now replace the measured field B^0,i in a similar procedure involving the next (third) sensor until the measurements from all sensors are used.

Ideally, the hardware should consist of a “main”, least disturbed sensor and additional sensors close to each major disturbance source as described in Sect. 2.2. Then, only the first-order correction for each sensor pair containing the main sensor is necessary to clean the data. However, other sensor configurations can also be used as described in the next section.

Figure 1 The components in the OB sensor system of the uncorrected measurements taken by the four magnetometers onboard GK2A on 4 March 2019. (a–b) FGMO, FGMI, AMR1, and AMR2.

4 Application to GK2A SOSMAG measurements

The GK2A spacecraft launched on 4 December 2018 in a 128.2^∘ E geostationary orbit is operated by the Korea Aerospace Research Institute (KARI) and provides meteorological and space weather monitoring over the Asia-Pacific region. The magnetic field vector is measured by the SOSMAG instrument (Magnes et al., 2020) at four locations on board the spacecraft. Two high-accuracy 3-axis flux gate magnetometers (FGMs) with a design similar to the THEMIS FGM instruments (Auster et al., 2008) are placed at the end (outboard sensor, FGMO) and 80 cm from the end (inboard sensor, FGMI) of an approximatively 1 m long boom. The other two magnetometers are 3-axis anisotropic magnetic resistance (AMR) (Brown et al., 2012) solid-state magnetometers placed on the body of the spacecraft.

The placement of the sensors can be seen in Fig. 6 of Magnes et al. (2020). Compared with the spacecraft dimensions, $\mathrm{290}\times \mathrm{240}\times \mathrm{460}$ cm, the magnetometer boom is relatively short, leading to strong spacecraft-generated disturbances at both FGM sensors.

As far as magnetic cleanliness is concerned, GK2A is a black box; i.e. no access to spacecraft operation time tables and to satellite-specific housekeeping data is available to aid the cleaning of the magnetic field data. Therefore the cleaning process must be based exclusively on the magnetic field measurements. Our goal is to eliminate the time-dependent spacecraft-generated disturbances from the FGM measurements. The strategy we adopt in order to take maximum advantage of the high accuracy of the FGMs and of the placement of the AMRs close to the disturbance sources is to first use the AMR measurements to clean the data from both FGMs and then use these corrected measurements to clean each other.

When a disturbance is much stronger at one sensor – as is the case for the AMR sensors – the scaling factor α is roughly given by the ratio of the magnitudes of the disturbance at the two sensors. This ratio is about 40 for AMR1, and 5 for AMR2, when paired with any of the FGMs. Since the sensor-specific noise for the flux gate magnetometers is lower by a factor of 20 compared to the AMR sensor noise, according to Eq. (22), the correction using the AMR sensors will introduce roughly the AMR noise divided by α. In particular, for the AMR2 one-fifth of its noise would be introduced in the corrected measurements. Since the same main disturbance is seen by both AMR sensors, no extra information is present in the AMR2 measurements; therefore we decided not to use the AMR2 sensor for removing the stray time-dependent spacecraft magnetic field. The 1∕40 from the AMR1 noise is much more favourable; therefore we will use this sensor to clean both FGM sensors' measurements.

4.1 FGM outboard and FGM inboard cleaning using the AMR1

Figure 1 shows the uncorrected measurements taken by the outboard (OB) FGM, the inboard (IB) FGM, and the two AMR sensors on 4 March 2019. We choose this day because it is representative for the routine operations, all the disturbance sources are active, and the ambient field shows little variance. Both step-like and spike-like disturbances can easily be seen in the picture. Among them, a prominent step-like disturbance between about 15:00 and 16:00 UTC is clearly detected by all four sensors, showing a very large magnitude at the AMR1. Note that the disturbance, which starts shortly after 15:00 UTC, affects the measurements until around 20:00 UTC. Because at 15:00 UTC the spacecraft is close to local midnight, we call this disturbance “midnight disturbance” (MD) to distinguish it from the other step-like disturbances. In 2019 this disturbance appears daily at the beginning and at the end of the year for about 14 weeks in total. We begin the cleaning of the data by first removing this disturbance from the FGM sensors' measurements using the AMR1 data.

For the sake of clarity, in the following we use the index s for the outboard FGM, the index t for the inboard FGM, and the index a for the AMR1 sensor. Equation (16), giving the magnetic field measured by the FGM sensors, corrected in the first order using the AMR1 sensor yields

$$\begin{array}{}\text{(27a)}& {\displaystyle }{\mathit{B}}^{\mathrm{1},sa}& {\displaystyle }={\mathit{B}}^{\mathrm{0},s}+{\mathbf{A}}^{\mathrm{0},sa}\left({\mathit{B}}^{\mathrm{0},s}-{\mathit{B}}^{\mathrm{0},a}\right),\text{(27b)}& {\displaystyle }{\mathit{B}}^{\mathrm{1},ta}& {\displaystyle }={\mathit{B}}^{\mathrm{0},t}+{\mathbf{A}}^{\mathrm{0},ta}\left({\mathit{B}}^{\mathrm{0},t}-{\mathit{B}}^{\mathrm{0},a}\right),\end{array}$$

with the matrices A^0,ja (where $j=s,t$) given by Eq. (17).

We select the time interval [15:10, 16:15] to isolate the targeted disturbance and use it to determine the variance directions of the disturbance and the scaling factors which give us the correction matrices. To lift the indetermination of the sign of the scaling factor α in Eq. (14), we compute the corrected fields Eq. (27) for both signs and keep the sign for which the disturbance is successfully removed. Equation (14) gives a very good estimation for the scaling factor. However, since this estimation uses the measured magnetic field, which includes the ambient magnetic field, it may slightly deviate from the correct value. To improve the precision, one may use the scaling factor determined from Eq. (14) as the initial value for a minimization procedure of the correlation between (ΔB^0,ja)_x and the corrected (B^1,ja)_x. While we found this to improve the determination of α for cleaning the FGMI data using FGMO measurements for days with disturbed ambient magnetic field, for cleaning the FGM data using AMR1 measurements on 4 March 2019 the minimization does not significantly change the value of α.

The angle between the direction of the disturbance at the AMR1 sensor and the direction of the disturbance at the inboard FGM sensor is 31^∘. For the outboard FGM sensor this angle is 25^∘, indicating that the disturbance source is not collinear with either of the sensor pairs. This is not surprising given the placement on the spacecraft body of the AMR sensors. Even so, the sum of the α coefficients differs from unity by less than 0.005.

Figure 2 The difference ΔB^st before the first-order correction was applied, represented on components in its corresponding VPS (the x axis is aligned with the maximum-variation direction of ΔB^st). The mean values were subtracted from all components.

The higher-order corrections should identify and eliminate disturbances roughly ordered by their strength at the AMR1 location. However, attempting the second-order correction only introduces spurious data in the FGMs measurements, increasing their variance. This is because the noise level of the AMR sensors is higher than the noise level of the FGM sensors and the AMR1 noise is added to the corrected measurements according to Eq. (22). Consequently we limit the AMR1 corrections to the first order.

Since the data cleaning on board the spacecraft should not require frequent updates of the correction parameters once uploaded to the spacecraft, it is necessary that the determined A correction matrices remain stable in time. In order to confirm this, we checked the stability of the cleaning parameters by using the same procedure once for every week showing the targeted disturbance in 2019. The standard deviation for the maximum-variance directions was below 1^∘, while the standard deviation for the scale factors was below 10⁻³. These low values are not surprising since for a given source the cleaning parameters depend only on its multipole character and on the geometry of the sources–sensors system. Other factors such as the intensity of the current generating the magnetic disturbance or the temperature do not influence the cleaning parameters. Applying the correction using the determined set of parameters removes the disturbance throughout the entire year of 2019.

Figure 3 The cleaning parameters resulting from the sliding window scan for the first-order correction of the outboard FGM. The top panel shows the number density of the maximum-variance direction (θ_w,φ_w) on a $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid. The bottom panel shows a histogram of the statistical distribution of the α_w coefficients. The grey filled bars are the coefficients corresponding to directions within 2.5^∘ from most probable direction. The red vertical line marks their mean value.

4.2 FGM cleaning using the AMR1-corrected data

We now use the AMR1-corrected FGMO and FGMI measurements given by Eq. (27) as a starting point in the iteration Eq. (23) for cleaning the FGMO data using the FGMI data and vice versa.

Unlike the single-step disturbance we dealt with in Sect. 4.1, the disturbances to be removed now show a repetitive pattern over the entire day, apparent in Fig. 1. Apart from the removed large-magnitude disturbance, one can visually identify at least two other types of disturbances in Fig. 1: step-like disturbances at a timescale of over 1 h and spike-like disturbances at timescales of minutes. To determine the correct cleaning parameters, the length of the analysis interval has to be chosen so as to contain many samples of the targeted disturbance but avoid including other disturbances. This can be accomplished by first eliminating the highest-frequency disturbances using a small enough interval length. Then the interval length is increased to encompass the next-highest-frequency disturbance.

In order to increase the precision of the cleaning and to have an indication as to the stability of the determined parameters, we compute the cleaning parameters using sliding windows covering the entire 24 h interval. For each window w we find the scaling factor α_w, elevation angle θ_w, and azimuth angle φ_w of the maximum-variance direction. After we scan the entire day interval, we determine the most probable direction (θ,φ) of the maximum variance which determines the rotation matrices R in Eq. (24). For this direction we select the corresponding coefficients α_w, and we compute their average value. At the end, the correction matrix A is computed using θ,φ, and α.

Figure 4 The initial AMR1-corrected FGM inboard measurements represented in the inboard VPS are plotted with the black lines. The first-order correction is plotted with red. Mean values were subtracted.

The disturbances can be much better identified in the difference $\mathrm{\Delta }{\mathit{B}}^{\mathrm{0},sta}={\mathit{B}}^{\mathrm{0},sa}-{\mathit{B}}^{\mathrm{0},ta}$ plotted in Fig. 2. The difference was first rotated in the VPS corresponding to a window length of 100 s, smaller than the time interval between the spike-like disturbances. In this coordinate system, different disturbance types tend to sort themselves on components.

The spike-like disturbances appear now in the x and y components with a cadence of 10 min and a magnitude larger than 10 nT during the entire interval. The step-like disturbances with slightly smaller magnitudes than the spikes are present in the y and z components. The duration between upward and downward variations of the step-like disturbances is 80 to 90 min, not as regular as the timing for the spikes. A new type of disturbance, not evident in Fig. 1, is now clearly apparent as a variation at higher frequencies (periods less than 1 min) than the cadence of the steps or the spikes. A closer investigation shows that this disturbance is irregular, with a maximum peak-to-peak amplitude of up to 4 nT in the x component and with its spectral power spread up to the Nyquist frequency.

The much smaller amplitude of the higher-frequency disturbance in the y and z components indicates its linear polarization. This was the disturbance which determined the orientation of the VPS used to plot the differences in Fig. 2. However, the spike-like disturbance has a large contribution to the x component; therefore its maximum-variance direction is not orthogonal to the maximum-variance direction of the high-frequency disturbance. In fact, the angle between the maximum-variance directions of the spike-like disturbance and of the high-frequency disturbance is 25^∘, which grossly violates the orthogonality condition. As a consequence, if the two sources producing the high-frequency and the spike-like disturbances have different scaling factors, the PiCoG method will not be able to remove both disturbances from the x-component using one single pair of sensors. The 75^∘ angle between the directions of the spike-like disturbance and the step-like disturbance is more favourable, but it will still prevent the complete removal of these disturbances simultaneously unless they have the same scaling factors. The closest to orthogonality is the angle between the directions of the high-frequency disturbance and of the step-like disturbance, which is 87^∘. Since the orthogonality condition is not fulfilled, to proceed further we must assume that the disturbances to be removed come from a small volume compared with the distances between the sensors, and therefore their scaling factors are not very different from each other. The results of the cleaning will either confirm or refute our assumption.

For the first-order correction we target the highest-frequency disturbance by choosing the same window length of 100 s used to compute the VPS for the difference plotted in Fig. 2. The statistical distribution for the resultant direction (θ_w,φ_w) of the maximum variance and a histogram of the α_w values is shown in Fig. 3. Both distributions exhibit clear isolated maxima, which is a strong indication that the targeted disturbance does not change its characteristics during the day interval. The angle between the disturbance directions at the two sensors is 15^∘, closer to collinearity than for the AMR1 correction.

Since the disturbances are larger at the inboard sensor, the effect of the correction is better illustrated for it than for the outboard sensor. The first-order correction of the inboard measurements for the first 4 h of the day is plotted in Fig. 4 with red over the initial AMR-corrected inboard measurements represented in the inboard VPS. The targeted high-frequency disturbance is eliminated from the x component. As apparent from the top panel of Fig. 4, between 00:40 and 00:50 UTC the high-frequency disturbance was switched off. One can see that the disturber also introduces a constant offset of about 5 nT, which is removed by the applied correction.

Figure 5 The first- (black) and the second-order correction (red) for the inboard FGM sensor. Mean values were subtracted.

The magnitude of the spike-like disturbance is much reduced in the x component of the corrected magnetic field in Fig. 4, so we conclude that the sources of both high-frequency and step-like disturbances are close to each other and are therefore removed together from the maximum-variance component. This justifies the application of the PiCoG method in this particular case when the directions of the two disturbances are far from orthogonal.

For the second-order correction we target the remaining spike-like disturbance by choosing a window width of 700 s. Figure 5 shows the result of the second-order correction for the inboard sensor. Both the targeted spike-like disturbance and the step-like disturbance are removed from the x component by this correction step, showing that indeed the distances between the sources of all three disturbances are much smaller than the distances between the disturbance sources and the FGM sensors, confirming our previous assumption.

The step-like disturbance and traces of the spike-like disturbance still remain in the y and z components in Fig. 5. To eliminate them, we select a window width of 16 000 s, enough to always include at least one step-like disturbance sample. As seen in Fig. 6, the correction removes the targeted disturbance and strongly reduces the remnants of the other two disturbance types from the x component. A leftover step-like disturbance, with a magnitude of about 1 nT, is still visible in the intermediate-variance component. This is due to the fact that, even with carefully chosen window lengths, the maximum-variance directions are still influenced by all present disturbances and therefore do not perfectly coincide with the polarization direction of the targeted disturbances. This leads to remanent disturbance on the other components. In our case, leftovers from the high-frequency disturbance interfered with the determination of the step-like disturbance polarization direction. The result is the further reduction of the high-frequency disturbance at the cost of not completely removing the step-like disturbance.

Figure 6 The second- (black) and the third-order correction (red) for the inboard FGM sensor. Mean values were subtracted.

We made use of the different characteristic timescales of the three disturbances treated in this section to help decouple them from one another even if their maximum-variance directions were not orthogonal and even if the amplitudes of the spike-like disturbances were not much different from the amplitudes of the step-like disturbances. If the disturbances had had the same timescales, these non-ideal conditions would have prevented the PiCoG cleaning method from working, unless some other specific properties of the disturbances could have been used to help decouple them.

To check the stability of the cleaning parameters, we determine them for every Sunday in 2019 with available data. The procedure produces very similar results apart from three instances when the ambient magnetic field was very disturbed. After eliminating the three outliers, we computed the standard deviations for the principal-component directions and for the scale factors, displayed in Table 1. The table also shows the corresponding maximum change in the corrected magnetic field on 4 March 2019 due to changes in the parameters equal to the standard deviations. The last row displays the maximum change due to the deviations in the parameters for one single order, while the parameters for the other orders are kept constant. Similarly, the last column displays the maximum change related to variations either in one single direction or in one scale factor. The last value in the table is the maximum change in the corrected magnetic field corresponding to all computed deviations, $B_{\text{max}}^{\text{dev}}=\mathrm{0.186}$ nT. This is the expected error due to the variations in the ambient magnetic field. However, the main error source is related to disturbers which do not fit our assumptions, such as collinearity, or to the presence of higher multipoles, as discussed in Sect. 5.

The standard deviations for the first two orders are very small, indicating very stable cleaning parameters for the high-frequency disturbance and for the spike disturbance. The third order, used to clean the step-like disturbance, shows larger deviations, especially for the outboard maximum-variance direction. This is because of the small contribution of the step-like disturbance at the outboard sensor which makes the procedure susceptible to the influence of the ambient magnetic field.

Table 1 Standard deviations for the directions and scale factors of the correction of FGMI data using FGMO measurements, together with the corresponding maximum deviation of the corrected magnetic field.

Figure 7 The final combined correction result for 4 March 2019 in the sensor system. The black lines show the original measurements taken by the outboard FGM; the red lines show the corrected data. The DC offset was restored to the value before the correction.

4.3 Parameters for spacecraft upload

Since the onboard correction is designed as a one-step linear combination of the measurements from different sensors, it cannot follow the iterative procedure described in Sect. 4. Therefore, we have to write the final correction in the form

$$\begin{array}{}\text{(28)}& {\mathit{B}}^{c,s}={\mathbf{M}}^s{\mathit{B}}^{\mathrm{0},s}+{\mathbf{M}}^t{\mathit{B}}^{\mathrm{0},t}+{\mathbf{M}}^a{\mathit{B}}^{\mathrm{0},a},\end{array}$$

where the superscript c stands for the combined correction, and the matrices M^j have constant coefficients given by the A correction matrices determined on the ground from the third-order correction, Eq. (26), applied to the first-order AMR correction of the FGM measurements, Eq. (27):

$$\begin{array}{}\text{(29a)}& {\displaystyle }{\mathbf{M}}^s& {\displaystyle }=\left(\mathbf{I}+{\mathbf{C}}^s\right)\left(\mathbf{I}+{\mathbf{A}}^{\mathrm{0},sa}\right),\text{(29b)}& {\displaystyle }{\mathbf{M}}^t& {\displaystyle }=-{\mathbf{C}}^s\left(\mathbf{I}+{\mathbf{A}}^{\mathrm{0},ta}\right),\text{(29c)}& {\displaystyle }{\mathbf{M}}^a& {\displaystyle }={\mathbf{C}}^s{\mathbf{A}}^{\mathrm{0},ta}-\left(\mathbf{I}+{\mathbf{C}}^s\right){\mathbf{A}}^{\mathrm{0},sa}.\end{array}$$

Here I denotes the identity matrix, and the matrix C^s has the form

$$\begin{array}{}\text{(30)}& \begin{array}{rl}{\mathbf{C}}^s& ={\mathbf{A}}^{\mathrm{0},st}+{\mathbf{A}}^{\mathrm{1},st}+{\mathbf{A}}^{\mathrm{2},st}\\ & +{\mathbf{A}}^{\mathrm{1},st}\left({\mathbf{A}}^{\mathrm{0},st}+{\mathbf{A}}^{\mathrm{0},ts}\right)\\ & +{\mathbf{A}}^{\mathrm{2},st}\left({\mathbf{A}}^{\mathrm{0},st}+{\mathbf{A}}^{\mathrm{0},ts}+{\mathbf{A}}^{\mathrm{1},st}+{\mathbf{A}}^{\mathrm{1},ts}\right)\\ & +{\mathbf{A}}^{\mathrm{2},st}\left({\mathbf{A}}^{\mathrm{1},st}+{\mathbf{A}}^{\mathrm{1},ts}\right)\left({\mathbf{A}}^{\mathrm{0},st}+{\mathbf{A}}^{\mathrm{0},ts}\right).\end{array}\end{array}$$

The AC correction described in Sect. 4 introduces a constant offset in the corrected data. This corresponds to the sources whose disturbances were removed. This DC offset can be determined by subtracting the mean value of the corrected measurements from the mean value of the original measurements:

$$\begin{array}{}\text{(31)}& {\mathit{G}}^s=\langle {\mathit{B}}^{\mathrm{0},s}\rangle -\left({\mathbf{M}}^s\langle {\mathit{B}}^{\mathrm{0},s}\rangle +{\mathbf{M}}^t\langle {\mathit{B}}^{\mathrm{0},t}\rangle +{\mathbf{M}}^a\langle {\mathit{B}}^{\mathrm{0},a}\rangle \right),\end{array}$$

where 〈⋯〉 denotes the average over a time interval longer than the timescale of the corrected AC disturbances.

In practice, there are additional DC offsets affecting the measurements, which are treated in a separate cleaning step. The vector G^s can be used to restore the original DC offset if a pure AC correction is desired.

$$\begin{array}{}\text{(32)}& {\mathit{B}}_{\text{pure AC}}^{c,s}={\mathbf{M}}^s{\mathit{B}}^{\mathrm{0},s}+{\mathbf{M}}^t{\mathit{B}}^{\mathrm{0},t}+{\mathbf{M}}^a{\mathit{B}}^{\mathrm{0},a}+{\mathit{G}}^s\end{array}$$

It follows from Eq. (29) that the sum of the M matrices is equal to the unit matrix:

$$\begin{array}{}\text{(33)}& {\mathbf{M}}^s+{\mathbf{M}}^t+{\mathbf{M}}^a=\mathbf{I}.\end{array}$$

A consequence of Eq. (33) is that an arbitrary vector added to the measurements B^0,s, B^0,t, and B^0,a in the expression of the offset G^s in Eq. (31) vanishes; therefore G^s is independent of the ambient magnetic field. This is to be expected because the magnetic field measurements enter the correction only as differences between distinct sensors; hence the correction – and therefore also the offset due to the correction – is determined only by the spacecraft-generated disturbances. This makes G^s a useful tool for monitoring changes in the DC offsets.

Applying Eq. (32) to the FGMO measurements yields the combined AMR1–FGMI correction to the outboard FGM measurements. We plot the original outboard FGM measurements in sensor system with black lines and the result of the combined correction with red lines in Fig. 7.

The M matrices were uploaded on GK2A 4 months after its launch. Since then the magnetic field measurements are corrected on board and transmitted to the ground stations within minutes from acquisition. The stability of the correction parameters is monitored, and a new set of parameters will be computed and uploaded in case changes in the spacecraft operation require a change in the parameters.

5 Errors and limitations

Even though we were able to eliminate most of the magnetic field disturbances on board the GK2A spacecraft, we need to be aware of the limitations the proposed method is subject to. We have already seen that, due to other disturbances or due to the ambient magnetic field variations, the maximum-variance direction might not coincide with the polarization direction of the disturbance to be removed. This difference will cause non-zero projections of the disturbance on the intermediate- and minimum-variance direction components which are not removed by the current applied correction. They will be reduced, however, by the next correction if the targeted sources lie close to each other. A disturbed ambient magnetic field may also interfere with the determination of the scaling factors. While there are ways to mitigate these effects, they are not within the scope of the present work.

Figure 8 The error due to sensor-specific disturbances and to the quadrupole disturbance introduced by PiCoG dipole disturbance correction. The x axis represents the position $x=r^t/r^s$ of the inboard sensor relative to the outboard sensor. Sensor noise is the same for both sensors. Each line corresponds to a fixed value of the quadrupole disturbance at the outboard position. Red lines, bottom to top: $b_q^s=$(0, 1 and 10)×Z^s. The black lines in between correspond to (0.2, 0.4, 0.6, and 0.8)×Z^s and (2, 4, 6, and 8)×Z^s, respectively.

An important benefit of the PiCoG method is the ability to treat up to three separate disturbance sources using measurements from two sensors. In order to be able to decouple the individual disturber contributions, two conditions must be satisfied: the disturbances must have well-defined polarization directions, and these directions must be orthogonal to each other. This may seem a strong condition to impose. However, apart from moving mechanisms such as reaction wheels, many, if not most, of the magnetic disturbances from a spacecraft come from current loops without phase delays and are therefore linearly polarized. The orthogonality, on the other hand, is not guaranteed. Even in the non-orthogonal case, disturbances coming from sources close to each other compared to the distance to the sensors share the same scaling factor (if both are either dipoles or quadrupoles) and are therefore removed together. A possible way to treat non-orthogonal disturbances coming from positions separated by large distances compared to the distances to the sensors is first transforming the data to a non-orthogonal system with its axes aligned with the maximum-variance directions of the three largest disturbers. This exercise is left for future examination.

Figure 9 (a) Power spectral density of the initial FGMO measurements (black line) compared with the PSD of the cleaned data (red line). (b) The ratio between the initial and cleaned PSDs for 4 March 2019

For each correction order, the disturbance to be removed has to be decoupled from the other disturbances. This is the case if the targeted disturbance amplitude is much larger than the amplitudes of the other disturbances, as assumed in Sect. 3.1. Another situation in which the disturbances can be decoupled is when they have different characteristic timescales. Then one may either use windowing in the time domain, as done in Sect. 4.2, or use band-pass filtering in the frequency domain. If the disturbances cannot be decoupled, then the cleaning procedure only works if the distances between the disturbance sources which cannot be separated are small compared to the distances to the sensors used for cleaning. In this case the cleaning parameters are similar, and the coupled disturbances are cleaned together.

One important class of error sources is additional disturbances which do not follow the determined scaling factor α or that are present at one sensor only. Among these are the sensor-specific noise; temperature effects, which sometimes cause sensor offset oscillations; and multipoles of higher order than the targeted disturbance. These disturbances are introduced into the cleaned magnetic field data either reduced or enhanced, depending on the sensor positions. In particular for GK2A, sensor offset oscillations triggered by large temperature gradients are quite significant, reaching peak-to-peak amplitudes up to 5 nT in the cleaned data (Magnes et al., 2020).

To estimate the error introduced by the sensor-specific noise combined with a quadrupole contribution additional to a dipole disturbance to be removed, let us assume a simple collinear geometry: a disturber placed in the origin of the coordinate system producing a disturbance characterized by both a dipole moment M and a quadrupole moment Q, an inboard sensor placed at the distance r^t characterized by a sensor-specific noise Z^t, and an outboard sensor placed at the distance r^s characterized by a sensor-specific noise Z^s. In these conditions, the correction of the dipole disturbance will introduce an error stemming from the quadrupole disturbance and the sensor-specific disturbances. The magnitude of the error will depend on the relative positions of the two sensors, on the sensor-specific noise, and on the strength of the quadrupole disturbance. After projecting on the principal-component direction, the magnetic field measured by the outboard sensor is (dropping the x-component index)

$$\begin{array}{}\text{(34)}& B^{\mathrm{0},s}=B+b_d^s+b_q^s+Z^s,\end{array}$$

where $b_d^s$ and $b_q^s$ represent the disturbance dipole and quadrupole contributions at the outboard sensor. A similar expression can be written for the inboard sensor.

The corrected field is obtained by applying Eq. (13a):

$$\begin{array}{}\text{(35)}& B^{\mathrm{1},s}=B+{\mathit{ϵ}}^{\mathrm{1},s},\end{array}$$

with the error ϵ^1,s given by

$$\begin{array}{}\text{(36)}& {\mathit{ϵ}}^{\mathrm{1},s}=(\mathrm{1}-\mathit{\alpha })Z^s+\mathit{\alpha }{Z}^t+(\mathrm{1}-\mathit{\alpha })b_q^s+\mathit{\alpha }{b}_q^t.\end{array}$$

Using the notation

$$\begin{array}{}\text{(37)}& x={\displaystyle \frac{r^t}{r^s}}<\mathrm{1}\end{array}$$

and keeping in mind that $b_d^s=x^{\mathrm{3}}{b}_d^t$ and $b_q^s=x^{\mathrm{4}}{b}_q^t$, Eq. (36) becomes

$$\begin{array}{}\text{(38)}& {\mathit{ϵ}}^{\mathrm{1},s}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-x^{\mathrm{3}}}}[\mathrm{1}+x^{\mathrm{3}}{\displaystyle \frac{Z^t}{Z^s}}+\left({\displaystyle \frac{\mathrm{1}}{x}}-\mathrm{1}\right){\displaystyle \frac{b_q^s}{Z^s}}]Z^s.\end{array}$$

Table 2 The initial magnitudes of the disturbances at all sensors and the final magnitudes in the corrected data for 04/03/2019. For the MD and for the spikes the sign shows the direction of the disturbance. AMR1 does not detect a quiet interval; therefore we cannot estimate the high-frequency disturbance magnitude at AMR1. The MD and the step magnitudes are defined as the size of their ramps. The magnitudes of the spikes are equal to the spike heights/depths. The high-frequency-disturbance magnitude is defined as the peak-to-peak amplitude. Samples of the disturbances affecting the z component of the outboard sensor (black lines), together with the corrected measurements (red lines) over 10 min intervals, are illustrated in the second column.

Similar to the findings of Neubauer (1975), the optimum position x results from a trade-off between the error due to the sensors, Z, and the error due to higher-order multipoles, b_q. We plot the error given by Eq. (38) for a number of quadrupole strengths in Fig. 8. The bottom red line corresponds to the pure dipole case. In this case, minimum error, equal to the outboard-sensor-specific noise, is obtained for x=0, i.e. for the inboard sensor placed at the position of the dipole disturber. As soon as a higher multipole is present, the inboard sensor must be moved away from the disturbance source in order to minimize the error. Already for a quadrupole disturbance at the outboard position equal to 1∕10 of the sensor noise, the optimum position of the inboard sensor is almost at the mid-distance between the disturber and the outboard sensor. When the quadrupole disturbance becomes equal to the sensor noise, the optimum distance becomes about 0.6x (middle red line). If the boom is very short, the quadrupole disturbance at the outboard sensor can reach very large values. The topmost red line in Fig. 8 corresponds to a quadrupole contribution 10 times as large as the outboard-sensor noise. In this case the optimum position of the inboard sensor approaches even more the outboard-sensor position (0.8x).

A way to estimate the overall performance of the cleaning is to compare the power spectral densities (PSDs) of the initial measurements with the power spectral densities of the cleaned data as shown in Fig. 9. The spectra in the figure were computed as the average over the entire day of 4 March 2019 using a sliding window of 512 s. Both PSDs contain not only the (remaining) disturbances but also the ambient magnetic field. Their difference shows the absolute total power of the removed disturbances, while their ratio represents the minimum factor with which the power of the disturbances is reduced. The mean of this factor for the 24 h interval shown in Fig. 9 over the frequency range covering periods from 2 s to 1 min is equal to 7.8. For lower frequencies, in the range covering periods from 1 min to 6 h we obtain a factor of 3.9 from the PSDs computed without windowing.

The success of the cleaning procedure can also be estimated for each individual disturbance class. The initial magnitudes of the disturbances targeted for cleaning are shown in Table 2 for each sensor. Values are given for each component in the outboard sensor reference system and for the module. The last column shows the remnants of the disturbances in the corrected data.

For the midnight disturbance we separated the leading ramp occurring around 15:00 UTC from the abrupt trailing ramp about 1 h later. The magnitude is computed as the difference between the median over 1.5 min of the field before and after the ramp. The leading ramp is reduced from about 34 nT in the FGMO measurements to less than 2 nT in the corrected measurements. The trailing ramp is reduced from 40 nT to about 1 nT. For the components, a positive sign denotes an upward ramp, and a negative sign a downward ramp.

The ramps of the step-like disturbances are symmetric; therefore we do not differentiate between the leading and the trailing ramps. The magnitudes are computed in the same way as for the MD. The mean step magnitude is reduced from 2.5 to 1.3 nT. However, note that the x component is more than doubled, from 0.4 to 1 nT. This is a necessary compromise we have to make because the polarization directions of the disturbances are not orthogonal, as discussed in Sect. 4.2.

The magnitude for the spikes was computed as the difference between the value of the peak of the spike and the median over 20 s intervals 5 s before and 5 s after the peak of the spike. For 4 March 2019 we obtain a mean magnitude of 13.8 nT for the initial FGMI measurements, 4.9 nT for the initial FGMO measurements, and 0.3 nT for the corrected measurements. For the components, a positive sign denotes upward spikes, and a negative sign downward spikes.

To estimate the reduction of the high-frequency disturbance, we use as disturbance-free etalon the quiet 10 min interval visible in Fig. 4 between 00:40 and 00:50 UTC. The magnitude of the high-frequency disturbance is computed as the difference between the mean peak-to-peak amplitude ($\mathrm{2}\sqrt{\mathrm{2}{〈B^{\mathrm{2}}〉}_{\text{time}}}$) of the measurements during the reference quiet interval – which is 0.2 nT for the corrected measurements – and the mean peak-to-peak amplitude over the adjacent interval between the next two spikes. The result is below 0.1 nT for all components and for the module. Note that, while AMR1 does not detect a quiet interval, it is still affected by a disturbance in the high-frequency range of about 18 nT peak-to-peak amplitude, possibly coming from another source(s). Despite the large amplitude of this disturbance at AMR1, the increase of the disturbance in the high-frequency range of the FGMO and FGMI measurements after the AMR correction is below 0.1 nT. All other discussed disturbances apart from the MD are lower at AMR1, which combined with the large-scale factor used for correcting the MD assures minimum transfer of these disturbances to the corrected data.

6 Summary and conclusions

We propose a multi-sensor method for removing spacecraft-generated AC disturbances from magnetic field data. The method employs principal-component analysis to decouple multiple disturbance sources and minimize the introduction of artefacts to the components free of the targeted disturbance.

A pair of sensors can resolve up to three independent disturbers. While no prior knowledge on the disturber source is required, linear polarization of the disturbance is assumed, and the polarization direction of different disturbers should ideally be mutually orthogonal. The method is robust enough to provide sensible results even if these assumptions are not strictly met. Of course, specific situations may provide additional opportunities to help separate distinct disturbers. One example is using the different characteristic timescales of the disturbances to determine the window lengths in Sect. 4.

There are, however, situations, such as non-orthogonal disturbances from sources with large spatial separation compared with the distance to the sensors, in which two sensors are not enough to remove the disturbances with the described method. Non-linearly polarized disturbances, as those produced by reaction wheels, need special treatment not covered by this work.

We applied the PiCoG cleaning method to the GK2A SOSMAG sensor configuration by first using the spacecraft-body-mounted AMR sensor measurements to remove large disturbances from the two boom-mounted FGM sensors. Three distinct types of disturbances were then removed using the two FGM sensor measurements: high-frequency disturbance in the range of less than 1 min, spikes occurring every 10 min, and steps occurring at intervals above 1 h.

We proved that on a specific day the method was able to reduce the spectral power of magnetic field disturbances by at least a factor of 7.8 in the period range of 2 s to 1 min and by 3.9 in the period range of 1 min to 6 h. These values are representative for the performance of the method over the entire year of 2019.

The final correction takes the form of a linear combination of the different sensor readings whose coefficients were determined on the ground. These coefficients were uploaded to the GK2A spacecraft, allowing for in-flight removal of spacecraft disturbances and near-real-time delivery of cleaned magnetic field data, essential for space weather applications. In the future we shall apply the PiCoG method for post-processing of data from other spacecraft, for example from BepiColombo (Benkhoff et al., 2010) and Cluster.