Robust Adaptive Spacecraft Array Derivative Analysis

Multispacecraft missions such as Cluster, Themis, Swarm, and MMS contribute to the exploration of geospace with their capability to produce gradient and curl estimates from sets of spatially distributed in situ measurements. This paper combines all existing estimators of the reciprocal vector family for spatial derivatives and their errors. The resulting framework proves to be robust and adaptive in the sense that it works reliably for arrays with arbitrary numbers of spacecraft and possibly degenerate geometries. The analysis procedure is illustrated using synthetic data as well as magnetic measurements from the Cluster and Swarm missions. An implementation of the core algorithm in Python is shown to be compact and computationally efficient so that it can be easily integrated in the various free and open source packages for the Space Physics and Heliophysics community.


Introduction
Exploration of the Earth's space environment and the heliosphere is driven by in situ measurements of fields and particles along the trajectories of dedicated spacecraft. Spatiotemporal ambiguities hamper the analysis and the proper interpretation of single-spacecraft data. In the past two decades, this issue has been addressed by a growing number of multisatellite missions such as Cluster, Themis, Swarm, and MMS. In the context of the four-spacecraft ESA mission Cluster, considerable attention was paid to the development of analysis tools for multispacecraft data (e.g., Paschmann & Daly, 1998. Of particular relevance is the multipoint analysis of spatial gradients and related derivatives such as the divergence and, very importantly, the curl of a vector field which gave rise to the term curlometer (Dunlop et al., 1988). The latter allows for an estimation of the electric current through Ampère's law (see Vogt, 2014, andDunlop et al., 2018, for recent reviews).
The quality of multipoint gradient estimates crucially depends on geometrical properties of the spacecraft array. Robert et al. (1998a) reviewed geometric quality parameters proposed for the four-spacecraft Cluster constellation and concluded that the most complete geometric characterization is through the eigenvalues of the so-called volumetric tensor R vol formed from dyadic products of the spacecraft position vectors. The array shape is captured by two parameters E (elongation) and P (planarity), defined in terms of eigenvalue ratios of R vol . The impact of the tetrahedral shape parameters E and P on the errors of curl and gradient estimates was studied by Robert et al. (1998b) and Chanteur (2000). De Keyser et al. (2007) and also Hamrin et al. (2008) applied least squares inversion with error control to more general satellite configurations over a finite period of time, with localization enforced through the choice of covariance matrices or weighting functions. Shen et al. (2012) discussed stable gradient and curl estimators in terms of the effective array geometry resulting from the nonzero eigenvalues of the volumetric tensor. This paper is concerned with a family of gradient and curl estimators that are based on a local least squares principle and formulated in terms of so-called reciprocal vectors q , = 1, 2, … , S, where S denotes the

Earth and Space Science
10.1029/2019EA000953 number of spacecraft in the configuration. Originally formulated for the four-satellite case (Chanteur, 1998), reciprocal vectors were generalized to arrays of more than four spacecraft in noncoplanar geometries (Vogt et al., 2008a), to three-satellite constellations (Vogt et al., 2009), and also to virtual planar arrays with four or more spacecraft (Vogt et al., 2013), a case that is relevant for Swarm measurements in the auroral zone. Using measurements b , = 1, 2, … , S of a vectorial observable b, the reciprocal vector concept allows for compact and suggestive estimator representations of spatial derivatives: The gradient matrix is estimated as ∑ q b t , the curl through ∑ q × b , and the divergence as ∑ q · b . For two-dimensional configurations, these estimators reproduce only the planar contributions to the spatial derivatives.
In section 2, we derive the Robust Adaptive Spacecraft Array Derivative Analysis (rasada) algorithm, that is, a compact recipe that integrates the various members of the reciprocal vector family of gradient and curl estimators and allows for efficient control of errors and degenerate array geometries. As demonstrated in a supplemented jupyter notebook, a fast and efficient implementation of the rasada algorithm in Python requires only a few lines of code, thus facilitating its adoption in free and open software toolboxes dedicated to the analysis of spacecraft data from geospace and heliospheric missions such as SpacePy (Morley et al., 2011), pysat , or HelioPy (Stansby et al., 2019), see also the recent review by Burrell et al. (2018). In section 3, the algorithm is validated and demonstrated through a series of applications to synthetic data and also to Cluster and Swarm measurements. All examples are transparently and reproducibly documented in the supplemented jupyter notebook, in support of the Geoscience Paper of the Future initiative (Gil et al., 2016).

Derivation of the Algorithm
Without losing generality and to ease the notation, the spacecraft positions are assumed to be given in a mesocentric reference frame: ∑ r = 0. The generalized reciprocal vectors q are solutions of R pos q = r , = 1, 2, … , S, where R pos = ∑ r r t denotes the so-called position tensor (Vogt et al., 2008a;Vogt, 2014), that is related to the volumetric tensor through R vol = 1 S R pos . Note that elsewhere, the position tensor is written without a subscript, but in this paper we reserve the plain symbol R for the position matrix, that is, the (S × 3) matrix R with the mesocentric spacecraft position vectors as its rows: R = (r ) , where = 1, 2, … , S and = 1, 2, 3. This implies that R pos = R t R.
The various flavors of reciprocal vectors allow for intuitive and compact representations of local linear least squares estimators for spatial derivative operators such as the gradient ∇b, the divergence ∇ · b, and the curl ∇ × b of a vector field b (Chanteur, 1998;Vogt et al., 2008a) using measurements b of b at the positions r : Here the gradient matrix G = ∇b is defined through the differential db = G t dr, to be consistent with the gradient vector g = ∇h of a scalar field h: dh = g · dr = g t dr. In Cartesian coordinates x , = 1, 2, 3, the matrix elements of G = ∇b are G k = b k , k = 1, 2, 3, and the gradient ∇b k of a Cartesian component b k is a column vector.

Reciprocal Matrix
Suppose the singular value decomposition (SVD) of the position matrix R is given by where U is a (S × 3) matrix and V a (3 × 3) matrix satisfying U t U = E and V t V = E with the (3 × 3) identity matrix E. The matrix V further satisfies VV t = E and is thus orthogonal. The elements of the diagonal matrix S are the singular values of R. Then R t = VSU t , and the position tensor can be written in the form The column vectors of the orthogonal matrix V are the eigenvectors of R pos . The corresponding eigenvalues R (1) pos ≥ R (2) pos ≥ R (3) pos ≥ 0 (note that R pos is positive semidefinite) are the squares of the singular values of R: 10.1029/2019EA000953 The set of linear equations R pos q = r , = 1, 2, … , S, is condensed into a single matrix equation R pos Q t = R t . After transposition and since R pos is symmetric, QR pos = R. The rows of the matrix Q are formed by the generalized reciprocal vectors q , = 1, 2, … , S. Note that elsewhere, the symbol Q is used to denote the tensor ∑ q q t , also called the (generalized) reciprocal tensor. For consistency and to avoid confusion, we write Q rec = ∑ q q t and reserve the plain symbol Q for the matrix of reciprocal vectors defined above, then Q rec = Q t Q. For brevity, we call Q the reciprocal matrix. Note that Q is the pseudoinverse of R.

Right-multiplying QVS
This is the defining equation for the reciprocal matrix Q, based on the SVD of the position matrix R, and corresponding to the set of equations R pos q = r , = 1, 2, … , S.

Nondegenerate Spacecraft Arrays in Three Dimensions
In the case of a regular three-dimensional array where not all S ≥ 4 spacecraft are in one plane, the position tensor R pos is invertible, and the generalized reciprocal vectors are q = R −1 pos r (Vogt et al., 2008a). In matrix representation, In order to obtain the matrix form of the estimator G ≃ ∑ q b t for the gradient matrix G = ∇b of a vectorial observable b, the measurements b , = 1, 2, … , S, are assumed to form the rows of a data matrix B.
The curl c = ∇ × b is estimated either through the skew-symmetric part 1 The gradient vector g = ∇h of a scalar observable h is estimated as g ≃ Q t h. Here the measurements h , = 1, 2, … , S, are assumed to form a data column vector h.

Array Geometry and Error Control
In the two-dimensional case where at least three spacecraft reside in one plane but do not degenerate to a linear configuration, R (1) pos ≥ R (2) pos > 0 and R (3) pos = 0. The third eigenvector is normal to the spacecraft plane and constitutes a singular dimension and without further information the normal gradient component cannot be estimated. Although a full inversion of the set of equations R pos q = r , = 1, 2, … , S, is not possible, we may compute minimum-norm solutions and still obtain gradient estimates in the spacecraft plane as Vogt et al., 2009Vogt et al., , 2013.
In matrix representation, we first compute a stable pseudoinverse of S through and then Q = US − V t (Vogt et al., 2009(Vogt et al., , 2013 with row vectors q , = 1, 2, … , S, as before. Finally, for one-dimensional arrays, R (1) pos > 0 and R (2) pos = R (3) pos = 0. In analogy with the two-dimensional case, we write In practice, the eigenvalues R (2) pos and/or R (3) pos may be nonzero but still very small, producing unstable estimates of spatial derivatives in the affected directions, that is, results that are very susceptible to instrumental and other random errors. Meaningful estimates can only be constructed in the subspace spanned by the eigenvectors corresponding to the larger eigenvalues, leading to an effective error control strategy closely associated with the geometry of the spacecraft array. Earlier studies following a similar approach include the reports of De Keyser et al. (2007), De Keyser (2008), Hamrin et al. (2008), and Shen et al. (2012).
Since this paper aims at a robust yet simple implementation of error control, we adopt a single tolerance parameter to stabilize the inversion of the position matrix. A direction along an eigenvectorv is taken into account only if the corresponding singular value ratio R ( ) ∕R (1) relative to the largest value R (1) is larger than the prescribed tolerance parameter . As a canonical value, we suggest = 0.1, corresponding to an 10.1029/2019EA000953 elongation of E = 0.9, and also a planarity of P ∼ 0.9 if the first two singular values are of comparable size. In the two-dimensional case where the third direction is excluded a priori by construction (Vogt et al., 2013), the coefficient number of the two-dimensional position tensor corresponds to the inverse square tolerance parameter: R (1) pos ∕R (2) pos = −2 = 100. In the case of planar or even linear arrays, measurements do not offer information on spatial derivatives along missing dimensions. Additional information in the form of geometrical or physical constraints can be utilized to construct the missing components (see Vogt et al., 2009Vogt et al., , 2013. Integration of constraint equations in the rasada algorithm should be straightforward but is beyond the scope of the present study and thus left for future work.

Error Estimation
Error analysis of gradient and curl estimators from the reciprocal vector family has been discussed, for example, by Chanteur (1998), Vogt and Paschmann (1998), and Vogt et al., (2008a, 2009. It turns out that positional errors are typically small compared to instrumental uncertainties. Another important contribution to the statistical error are random deviations from the supposedly linearly varying magnetic field profile, for example, small-scale magnetic fluctuations that do not show a correlation between the different sensors. Assuming the effective (instrumental plus random deviations from linearity) statistical errors are isotropic and mutually uncorrelated, they can be captured by a scalar quantity b, and the error covariance matrix of a gradient vector estimate g = ∇b k (for the Cartesian component b k of a vectorial observable b) is given by The mean square magnitude error is thus where ||Q|| is the Frobenius norm of the reciprocal matrix Q defined through . The mean square error of the divergence estimator has the same form as the square magnitude error: The error covariance matrix of the estimated curl vector c can be written in the form where E denotes the identity tensor (see Vogt et al., 2019). In the special case of a planar spacecraft array, ⟨ | c n | 2 ⟩ = ( b) 2 ||Q|| 2 for the mean square error of the curl component normal to the spacecraft plane (Vogt et al., 2019).

Implementation
The preceding derivation yields the following algorithm for spatial derivative estimators from the reciprocal vector family, derived from a local least squares principle.
1. Store the mesocentric spacecraft position vectors r , = 1, 2, … , S, as rows of the position matrix R. 2. Compute the SVD of the position matrix: R = USV t . 3. Using the nonzero elements of the diagonal matrix S, form its stable pseudoinverse S − . 4. Compute the reciprocal matrix Q, that is, the stable pseudoinverse of the position matrix as Q = US − V t .
The rows of Q are the (generalized) reciprocal vectors q , = 1, 2, … , S. 5. Assuming that the measurements b , = 1, 2, … , S, of a vectorial observable b form the rows of a data matrix B, the gradient matrix G = ∇b is estimated as G ≃ Q t B. The curl c = ∇ × b is estimated either using the skew-symmetric part of G or through c ≃ ∑ q × b , using the row vectors q of Q. 6. Isotropic and mutually uncorrelated instrumental errors b for each component b k of b produce a gradient vector (∇b k ) error covariance matrix of the form ( b) 2 Q t Q. The gradient magnitude error and, in the case of planar spacecraft arrays, the error of the normal curl component assume the form ||Q|| b where ||Q|| is the Frobenius norm of the reciprocal matrix Q.

Applications
The rasada algorithm is demonstrated using a jupyter notebook, available as Supporting Information S1. It is implemented in Python. The numerical part of the procedure rests solely on the numpy package, graphical display utilizes matplotlib, CDF files are read using CDFlib, and the pandas package helps to select data and to deal with data gaps.

Synthetic Data, Nondegenerate Satellite Arrays
The first section of the jupyter notebook introduces the rasada algorithm by means of two types of nondegenerate synthetic array orbit data.
In the first example, the reference four-point array is a regular tetrahedron with unit distances between its vertices, multiplied by a length scale L, and then propagated with constant velocity to produce artificial array orbit data. Following Steps 1-4 of the rasada procedure, the reciprocal matrix Q is constructed. In Step 5, the spatial derivatives are estimated for a linear noise-free vector field b = 1 2 c × r where c = ∇ × b is a constant vector.
The second example illustrates also Step 6, that is, error estimation. Since for the highly regular and symmetric array geometry of the first example all singular values are identical, implying that all potential correlations between the errors of the estimated curl components are zero, a less regular geometry is chosen where two vertices are shortened. The linear field is contaminated by white noise b to produce an ensemble of erroneous measurements of the form b(r) = 1 2 c × r + b(r). Ensemble statistics are compared with the error estimates produced through Step 6 of the rasada algorithm. Within statistical uncertainties, the results are fully consistent.

Synthetic Data, Partially Degenerate Satellite Arrays
The second section of the jupyter notebook makes use of partially degenerate synthetic satellite arrays to illustrate rasada's geometry control. Starting from a regular tetrahedron, only one spacecraft is allowed to move uniformly toward the plane spanned by the other three, temporarily producing a planar configuration. The effective array dimension parameter, computed as the number of singular values that are kept in rasada Step 3 when the pseudoinverse S − is constructed, jumps from three to two for the degenerate case, and the in-plane gradient component is obtained. In the same spirit, a transition from a regular quad to a linear configuration is simulated, also producing the expected results. Finally, the important case of normal curl estimation from virtual planar quads is considered as this is relevant for the determination of field-aligned currents (FACs) from the pair of SwA and SwC spacecraft (see Blaǎgǎu & Vogt, 2019, and section 3.4). Using ensembles of synthetic measurements, also, the error estimates are tested and confirmed.

Current Vectors Estimated From Cluster FGM Data
In the third section of the supplemented jupyter notebook, the rasada algorithm is applied to Cluster multisatellite fluxgate magnetometer (FGM) data measured on 11 June 2001 in the time interval 20:05-20:23 UT during a magnetopause crossing. The event was analyzed by Dunlop and Balogh (2005) and is also included on the multispacecraft webpage of the Cluster Science Archive (CSA; https://www.cosmos.esa.int/web/csa) as an application of the curlometer technique (Dunlop et al., 1988). CDF files from the CSA contain data gaps that are identified, flagged, and eliminated. The three components of the electric current vector j = −1 0 ∇ × b are computed and displayed (see Figure 1). The current components compare nicely with the results of Dunlop and Balogh (2005) and the plot published on page 18 of the CSA multispacecraft tutorial (ESDC-CSA-TN-001: The Curlometer Technique, Issue 1 Rev 1, 13 December 2016).

FAC Densities Estimated From Swarm VFM Data
The fourth section of the supplemented jupyter notebook deals with the estimation of auroral FAC densities using vector field magnetometer (VFM) data measured by ESA's Swarm spacecraft. To produce FAC density estimates from planar array data, the component of the curl vector normal to the spacecraft plane (i.e., the only meaningful component in this geometry) is extracted, the normal unit vector is obtained as the third row vector of the matrix V t , and the inclination of the ambient magnetic field with respect to the normal unit vector is computed. The procedure is summarized in a general-purpose function that can be applied to planar arrays formed by arbitrary numbers of satellites.
The first time interval chosen to demonstrate rasada FAC density estimation from planar satellite configurations is 08:49-08:52 UT on 23 May 2014. The third Swarm spacecraft SwB was sufficiently close to the lower pair SwA and SwC to yield nondegenerate three-point configurations and thus meaningful three-satellite FAC estimates. The resulting profile (displayed in the supplemented jupyter notebook) compares well with the FAC density estimate obtained by Blaǎgaǎu and Vogt (2019) who studied the same event using planar reciprocal vectors (Vogt et al., 2009). The only notable difference is a small vertical offset of a few tenths of μA/m 2 at the beginning of the time interval with the lower satellite pair close to the geographic pole, possibly related to nonlinear contributions from the ambient magnetic field present in the Level-1b magnetic field data (total field without model subtraction) that serve as input data to rasada.
To produce dual-spacecraft FAC estimates from VFM data measured by the two close Swarm satellites SwA and SwC, virtual four-point configurations are constructed using shifted positions and magnetic field data (Ritter & Lühr, 2006;Ritter et al., 2013;Vogt et al., 2013). Figure 2 shows rasada FAC density estimates in comparison with Level-2 FAC data products provided by ESA. The upper diagram shows the results for the time interval 19:35-19:40 UT on 11 July 2014, with stable and nondegenerate virtual four-point configurations. The Level-2 dual-s/c FAC data product is produced from filtered magnetic field data; hence, the overall profile is smoother than the rasada estimate that operates on the raw Level-1b magnetic field measurements without any filtering. Otherwise, the two dual-s/c FAC density estimates compare very well. This is remarkable because the rasada algorithm is applied to magnetic field data not only without any filtering but also without subtraction of a model magnetic field.
The lower diagram of Figure 2 shows a comparison of FAC density estimates for the time interval 11:19-11:23 UT on 27 September 2014. Since the chosen time interval includes the crossover point of the two satellites in the southern hemisphere, the virtual quad configuration temporarily degenerates to a quasi-linear array so that rasada's adaptability to the effective array dimension can be tested. To facilitate a comparison with the single-satellite (SwA and SwC) Level-2 FAC data products, the two time series are condensed into one profile after smoothing and averaging. While the ESA Level-2 dual-satellite FAC data product is not computed in the so-called exclusion zone within four degrees around the geographic pole due to errors caused by array degeneracy, the rasada dual-satellite FAC density estimate is available in the entire time interval and very consistent with average single-spacecraft FAC density (apart from a small vertical offset, possibly related to 10.1029/2019EA000953 differences in the input data). The rasada algorithm, through geometry adjustment based on the degeneracy parameter , provides proper dual-satellite estimates where the array can be considered two-dimensional and falls back to a single-spacecraft approach when the virtual quad degenerates into a quasi-linear array.

Conclusions and Outlook
In the context of spatial derivative estimation from multisatellite array data, the rasada algorithm introduced in this paper combines all earlier approaches of the reciprocal vector family (Chanteur, 1998;Vogt et al., 2008aVogt et al., , 2009Vogt et al., , 2013 in a unifying framework. By construction, the procedure adapts not only to the number of spacecraft in the array but also to its effective geometrical dimension and is robust in the sense that the errors resulting from degenerate configurations are constrained. The chosen implementation in Python is fully vectorized and thus efficient in terms of processing time and computational resources. Integration in Python-based analysis packages such as SpacePy (Morley et al., 2011), pysat , or HelioPy (Stansby et al., 2019) is expected to be straightforward.
Synthetic data are used to demonstrate the core qualities of the rasada procedure, namely, efficient gradient and curl estimation and array geometry control. The error formulas are shown to be consistent with the results of numerical Monte Carlo experiments. The procedure is validated through a reanalysis of Cluster FGM measurements of a magnetopause crossing and then illustrated further using Swarm VFM data to estimate the density of auroral FACs.
The versatile algorithm is not restricted to applications of the type discussed in this report. For instance, using dual-satellite measurements of the Swarm lower pair SwA and SwC, the number of points in the virtual configuration can be easily extended from four to arbitrary multiples of two by adding additional data along the track of each satellite, thus realizing data-adaptive smoothing instead of scale-based filtering as currently employed in the procedure to produce the official ESA dual-satellite FAC data product (Ritter et al., 2013). Also for Swarm three-satellite configurations, two or more position and magnetic field along-track data can be stacked to produce fully three-dimensional array configurations, so that for stationary situations, the full current vector can be obtained as done recently by Dunlop et al. (2015) using the curlometer technique.
In addition to multisatellite spatial derivative estimation, reciprocal vectors are employed beneficially, for example, in discontinuity analysis (Harvey, 1998;Vogt et al., 2011) and in wave vector estimation (Vogt et al., 2008b). It will be straightforward to generalize the rasada algorithm to those multispacecraft analysis tasks. Financial support from the Deutsche Forschungsgemeinschaft in the context of the SPP 1788 Dynamic Earth through Grant VO 855/4-1 is acknowledged. The work in Romania was supported by the ESA project MAGICS, PRODEX Contract 4000127660. Swarm Level 1b and Level 2 data are available from ESA (https:// www.esa.int/Swarm). Cluster FGM data are provided by the Cluster Science Archive (https://www.cosmos. esa.int/web/csa).