Particle Swarm Optimization for Magnetometer Calibration With Rotation Axis Fitting Using In-Orbit Data

This article demonstrates the performance of an improved particle swarm optimization (PSO) algorithm with scalar checking and rotation axis fitting objectives using in-orbit data, which is obtained from two CubeSats missions, Aalto-1 and ESTCube-1, as well as simulation as reference. The improved algorithm uses sequential objectives refinement process to combine the two optimization objectives. This improvement addresses some challenges of magnetometer calibration when using in-orbit data. First, the change in the magnetic field vector direction at different points in orbit which is uncorrelated to the rotation of the spacecraft itself. Second, the uncertainty of the rotation axis information used as the reference, e.g., from gyroscope noise. Third, the available data set is heavily affected by the rotation mode of the spacecraft, which imposes some limitation in the rotation axis information needed by the algorithm. The improved PSO algorithm is applied on simulated data in order to analyze the calibration performance under different spacecraft tumbling rates and noise levels. In ideal condition (varying rotation axis during measurements and sufficient sampling rate relative to the spin rate), the rotation axis fitting objective can reach ∼0.1° of correction accuracy.


I. INTRODUCTION
Attitude determination and control system performance of a spacecraft is a combination of its sensors' and actuators' technical performance and calibration, which ensures their reliability in acceptable range of error. One of the typical and critical attitude sensors in low earth orbit (LEO) spacecraft is a three-axis magnetometer. This sensor provides a very good long-term stability because of the well-defined knowledge on the Earth magnetic field vector model, such as the International Geomagnetic Reference Field model [1].
In practical applications, the influence of multiple error sources distorts the magnitude and direction of the magnetic field vector measured by the magnetometer. A calibration procedure is required to minimize these errors, and for LEO spacecraft application, most calibration processes are based on scalar checking algorithm, which only relies on the ambient magnetic field magnitude. Thus, in this case, the calibration procedure does not depend on attitude knowledge. However, this relaxed knowledge requirement implies that the usable information for calibration process is limited, resulting in an underdetermined estimation problem. This appears mathematically as a constraint in the calibration parameters, particularly in the calibration matrix: the matrix is assumed as a diagonal or triangular matrix (3-6 independent elements instead of 12 elements in a full 3 × 3 matrix) [2]- [9]. Physically, this means that the distortion which can be corrected is limited to nonorthogonality and nonuniform scaling of the axes, with the assumption that at least one axis is perfectly aligned with reference, i.e., the magnetic field scalar based calibration cannot resolve rotational error.
Several approaches have been discussed in different studies to resolve this rotational error, although the calibration parameters estimation is still sensitive to larger error in nondiagonal values of the total scaling matrix [10] or requires additional information from sensor unsuitable for space application [11]. Other approach requires an intermediary knowledge on the system attitude information and use a neural network to fit the known model of the ambient magnetic field with the measured magnetic field [12]. Riwanto et al. [13] have proposed a method to resolve this rotational error component using rotation axis information derived from the magnetometer reading and fit it into a reference rotation axis information, which, in spacecraft settings, can be read directly from gyroscopes or determined from other sensors. This rotation axis fitting method was used complementarily with the scalar checking method implemented inside a particle swarm optimization (PSO) algorithm. The testing method, however, only covered ground-based test under a specific rotation mode, where the spacecraft is rotated around one axis at a time, with at least two different rotation axes performed separately.
This article demonstrates the performance of an enhanced version of the rotation axis fitting algorithm, originally described in [13]. The algorithm is applied on simulated data that models the attitude dynamics in orbit as well as using real in-orbit flight data from Aalto-1 [14], [15] and ESTCube-1 [16], [17] satellites. This means that the algorithm robustness is tested under the changing magnetic field direction experienced by the spacecraft in orbit due to Earth's magnetic field vector variation with relation to the position in orbit, which will appear as if the magnetometer is rotating. The approach also takes into account the presence of noise in the reference rotation axis information, which, in this case, is represented by on-board gyroscope measurement noise. However, as the algorithm is using a PSO architecture, it inherits the disadvantage of the previous algorithm in computation cost and the nature of offline estimation. This means that it is not suitable for real-time calibration routine and is more suitable for establishing initial estimate for calibration using online estimation methods (e.g., filters) or in fault detection and isolation scenario.
In order to describe the design and validation methods of the enhanced algorithm, this article is organized as follows. Section II briefly presents the mathematical model of the magnetometer and gyroscope measurements used for simulations and calibration algorithm. Section III describes the implementation of PSO algorithm for optimizing the scalar checking and rotation fitting objectives. Section IV describes the test setups and results.

II. TEST ENVIRONMENT MODEL
The scalar-checking algorithm uses the mathematical model of the magnetometer and apply calibration parameters to fit the measurement data with the reference Earth magnetic field. Additionally, the rotation axis fitting algorithm uses the estimated rotation axis from magnetometer data to fit into the reference rotation axis; in this article, the reference rotation axis is gathered from gyroscope measurement. Both magnetometers and gyroscopes are affected by error sources, described in their models in the following subsections.

A. Magnetometer Model
The simplified mathematical model that describes the relationship between the measured magnetic field vectorb and the reference ambient magnetic field b iš where 3 × 3 matrix S m and 3 × 1 vector off m are the compounded calibration parameters from the combination of individual magnetometer error sources such as scaling factors, nonorthogonality, internal bias, soft iron, and hard iron errors described more extensively in [13], while the vector η m is the magnetometer measurement noise. Note that in magnetometer calibration, the calibration parameters need to transform the measured magnetic fieldb into the corrected magnetic fieldb, flipping (1a) intô where the ambient magnetic field b is now substituted by the corrected magnetic fieldb, while the measurement noise η m is assumed as zero-mean Gaussian random process. K m = S −1 m and k m = off m are the calibration parameters that need to be optimized by the calibration algorithm. With this definition, the scalar checking algorithm will serve to minimize the difference in the magnitude of the known reference field b and the corrected fieldb.

B. Gyroscope Model
A three-axis gyroscope is used as the source for reference axis information in the rotation axis fitting method. Gyroscope measurement can be described aš whereω is the measured angular rate, S g is the gyroscope total scale error, ω is the true angular rate, off g is the angular rate bias with a drifting ratė which is also known as the rate random walk (RRW), while η ARW is the measurement noise or angular random walk (ARW). Both η ARW and η RRW are assumed as zero-mean white noise. The variance of these noise sources can be determined, for example, using Allan Variance [18]- [20].

III. PSO ALGORITHM
This article implements a PSO algorithm based on the work in [13], where a rotation axis fitting objective was proposed to calibrate the rotation error factor that cannot be resolved by scalar checking method alone. In [13], however, the rotation axis fitting objective was only validated for ground-based test conditions, i.e., in the absence of varying magnetic field magnitude and with discrete separation between different rotation axis experienced by the magnetometer, where measurements from each rotation axis can be grouped into one measurement locus. The previous algorithm version also requires manual fine-tuning of the individual fitness value's weight for the weighted aggregate approach it used in combining the two objectives.
This article extends the algorithm validation with inorbit data and implements a different approach in the rotation axis fitting objective definition as well as the multiple objectives combination method. This improves the algorithm robustness in calibrating data from actual in-orbit conditions without the need to fine-tune the fitness values' weights.
To fully describe the implementation of the modified PSO algorithm, the following subsections are arranged as follows. Section III-A describes the summary of the barebone PSO algorithm, Section III-B describes the summary of the two objectives (scalar checking and rotation axis fitting) used for calibration, and Section III-C describes the architecture of the implemented PSO algorithm with the two different objectives.

A. Basic PSO Algorithm
The barebone PSO algorithm implemented in this article is identical with the previous work in [13]. A short summary on the key equations governing the swarm behavior is described in the following steps. 1) Swarm initialization: Initialize n p numbers of particles (the swarm size) randomly in n c -dimensional search space (the swarm dimension). Each particle is defined with two properties "position" p i (k) and "velocity" v i (k) where the subscript i is the particle index (i = [1, n p ] ∈ Z), subscript j is the particle component index ( j = [1, n c ] ∈ Z), and k is the iteration index.
In PSO, the number of optimized parameters also translates directly into the swarm dimension n c (i.e., the dimension of p i ). Each particle has its own local best position pbest i , which is the position with best fitness value in the search history of each particle i. Local best position is initialized as its initial position (pbest i (k = 0) ≡ p i (k = 0)). Additionally, a global best position gbest is defined as the best position in the whole swarm, i.e., the local best position of the particle with the best fitness. 2) Swarm positions update: The velocity and position of each particle are updated with the basic formula v i j (k + 1) = wv i j (k) + rand(0, 1)c 1 (pbest i j (k) − p i j (k)) + rand(0, 1)c 2 (gbest j (k) − p i j (k)) (3a) where rand(0, 1) is a random number in the range of [0 . . . 1] ∈ R, w is the inertia weight parameter, c 1 is the cognitive rate parameter, and c 2 is the social rate parameter. 3) Best position evaluation: pbest i for each particle i is updated when the fitness value during that iteration k is better than the last known best fitness. The magnetometer calibration is a minimization problem; thus, the best position at iteration k is defined as where pbest i and gbest are the local and global best positions which represent the known positions in the search space with the minimum fitness value for each particle i and among all particles, respectively. f () is the fitness function which returns the fitness value which will be evaluated in every iteration to generate the best local and global positions (F b i and F b g , respectively), which are defined as After the local best position for each particle is evaluated with (4a), the global best position is then reevaluated with (4b): If a new local best position with a fitness value lower than the current global best position exists, then the global best position will be updated with the new, better position-else, the global best position stays the same. 4) Iteration evaluation: The iteration ends when the required goals are met, e.g., global best fitness reached threshold value, number of iteration reached limit, or the swarm converged to the best global position in a certain range-else, the iteration will continue. This iteration termination marks the end of a single PSO run, where the algorithm then returns gbest as the optimization solution.

B. Optimization Objectives
In this article, instead of being optimized as a full 3 × 3 matrix, where all 12 components are optimized in a single run of multiobjective PSO as it was implemented in [13], the calibration matrix K m is decomposed with QR decomposition into orthogonal and upper triangular matrix. The upper triangular matrix K Rm represents distortion of the magnetometer measurement locus with the z-axis perfectly aligned with the reference z direction (e.g., the spacecraft body frame), while the orthogonal matrix K Qm represents the rotational error of the magnetometer measurement locus. Thus, K Qm is a direction cosine matrix describing the angular displacement of the magnetometer z-axis from the reference.
Scalar checking objective is only capable of solving bias, scaling, and skewing errors. Thus, the scalar checking method will only optimize the components of the upper triangular matrix K Rm and bias vector k m for every particle i, which are defined as Complementarily, rotation axis fitting objective will optimize the components of the orthogonal matrix K Qm . However, K Qm is not directly solved with PSO because from the nine elements in an orthogonal matrix, only three parameters are independent. To translate the actual physical rotation properties of K Qm into the PSO search space, K Qm could be represented in the form of quaternion (four parameters), Euler angle rotation sequence (three parameters), or axis-angle (four parameters), among other forms. In this article, K Qm is implemented in Euler rotation sequence form.
For Euler rotation sequence form, the PSO search parameters can be converted into the rotation matrix K Qm with the straightforward equation where s and c denote sine and cosine, and the subsets 1, 2, and 3 denote the Euler sequence rotation angle defined by K Qm1,i , K Qm2,i , and K Qm3,i , respectively, e.g., c 1 = cos(K Qm1,i ), s 3 = sin(K Qm3,i ), and so on. While it is possible to do this conversion in various sequence orders and notations, (6) is derived from Z-X-Z Euler sequence.
Conversions from different notations and representations (e.g., using quaternion search space) are available in literature [21]. For this article, the two objectives are optimized in separate PSO runs. Thus, from the components of K Rm , K Qm , and k m , we define two separate PSO search spaces Then, the final calibration matrix is calculated back from QR components with Note that (7b) shows K Qm defined with three parameters (e.g., euler angle sequence). For four-parameter representation (e.g., for quaternion or axis-angle form), an additional parameter K Qm4,i is required.
With this definition of the search space, we then define the PSO objectives.
1) Scalar Checking: Scalar checking objective minimizes the difference between the calibrated magnetic field vector magnitude |b| with the known reference magnetic field magnitude |b|. We can translate this mathematically as the fitness function or by substituting (1a) into (9a) where i = [1, n p ] ∈ Z and s = [1, n s ] ∈ Z are the index for swarm particles and measurement data, respectively. Note that the scalar checking objective only optimizes the upper triangular matrix part K Rm in (8), while the orthogonal rotation matrix part K Qm is optimized by the rotation axis fitting objective.
2) Rotation Axis Fitting: Instead of grouping multiple magnetic field measurements into a predetermined number of loci that correlate to a particular rotation plane as proposed in [13], which assumes a relatively constant rotation Algorithm 1: Algorithm for Rotation Axis Estimate With Direct Vector Cross Product. 1: procedure RotEstb,ω 2: for i = 1 . . . n p swarm particle i 3: for l = 1 . . . n l continuous measurement locus l 4: for sl = 1 . . . n sl measurement in the locus sl 5: calculate rotation axis vector with cross product: 6:n i,l,sl = norm(n i,l,sl ) normalize into unit vector 7: average and normalize reference rotation vector: ω l,sl = norm ω l,sl +ω l,sl+1 +ω l,sl+2 /3 8: end for 9: end for 10: end for 11: end procedure plane, this article updates the method for defining the rotation axis fitting objective in order to accommodate the magnetic field measurement locus of an in-orbit spacecraft, where the rotation mode of the spacecraft does not guarantee a perfectly stable rotation plane(s).
The measurement locus is defined by grouping each continuous measurements as a single locus l, which contains estimated calibrated datab i,l,sl , where sl = [1, n sl + 2] is the index of measurement data in that specific locus. In each locus, a rotation axis vector is estimated from every three magnetometer measurements, as defining a unique plane requires a minimum of three points in 3-D space-this means with n sl + 2 number of measurements in each locus, we can estimate n sl number of rotation axis. The calculated rotation axis is then compared with reference rotation axis direction using vector dot product in the fitness function where the rotation axis direction for every three measurement pointsn i,l,sl (subset sl indexes three-measurements pair in the locus) and the reference rotation axis direction ω l,sl (averaged over the three-measurements period) can be calculated directly with Algorithm 1.
It is important to note that since the rotation axis is estimated with a direct cross product, under ideal condition (zero noise and constant ambient magnetic field), the angular displacement between the two continuous magnetic field vectorsb i,l,sl andb i,l,sl+1 should not exceed 180°or the estimated rotation axis vector will be inverted. This also sets the condition of the continuous measurements that can be grouped as a single locus l: The magnetometer sampling rate needs to be at least twice the rotation rate, and delays in measurement interval that exceed half the rotation period should be discarded or grouped into separate locus.
For this article, the reference rotation axisω l,sl is obtained directly from gyroscope measurementω, which is sampled at the same time as the magnetometer. Alternatively, reference rotation axis can also be obtained from direct observation (e.g.,in preflight test) or estimated from vector sensors such as star trackers or other image processing techniques [22]- [25]. Thus, the rotation axis fitting function F 2,i represents the averaged sum of angular difference between the estimated rotation axis from magnetometer data and the reference rotation axis. Note that in (10), both the estimated and reference rotation axes are unit vectors, and, thus, F 2,i ranges from 0 to 2.

C. Multiobjective PSO Architecture
The more detailed PSO architecture used in a single iteration is mostly identical with the original study in [13]. This section will expand more on the updated components of the PSO architecture and parameters selection. 1) Initialization Condition: As the improved algorithm needs to deal with two separate swarms defined in (7), a set of swarm has to be initialized separately for each objective. The swarm positions and velocities initial states are determined by pmax and vmax, respectively. In this implementation, both parameters determine the search space boundary in which the swarm states are randomly generated during initialization. This process is described mathematically as where the subscript j is the swarm component index. The structure pmax j = [pmax s, j , pmax e, j ] defines the swarm position initialization range, while vmax j defines the maximum allowed velocity of the swarm. For scalar checking objective, the initial positions of K Rm and k m components are set to 0.1 unit around an identity matrix and 5000 nT around a zero bias vector respectively or using the same mathematical notation as (5) This represents an a priori estimate of the calibration parameters with no scaling and bias error, where the swarm is allowed to initially explore from that region until it finds the global optimum.
For rotation axis fitting objective, the initial positions of K Qm components cover the whole rotation range from 0 to 2π radians or in the same notation as (6) and (7b): which allows the swarm to explore every possible rotation sequence combinations possible from the beginning.
2) Boundary Condition: The main goal of boundary condition is to limit the swarm search dynamics from "exploding" in the search space. In the updated algorithm, only the boundary condition for the rotation axis fitting objective is different, while the rest of the implementation is identical with the previous version in [13]. This difference is because in rotation axis fitting objective, the search space of an Euler rotation sequence is a set of three rotation angles, and, thus, the search space is cyclic between 0 and 2π radians instead of expanding into infinity. This specific search space boundary characteristic is implemented in order to preserve continuity in the physical representation of the rotation transformation.
3) Sequential Objectives Refinement: This article implements a different approach in the algorithm that combines the two optimization objectives (represented by F 1 and F 2 , which translates to the best swarm gbest 1 and gbest 2 , respectively) and evaluates the termination condition when compared to the previous version in [13]: Instead of optimizing the weighted addition of the two fitness values in one PSO iteration, each PSO run, as summarized in Fig. 1, will optimize the two objectives alternatively and combine the results using a heuristic which will recognize stagnation in the fitness value as the final termination condition. This is implemented by adapting the refinement procedure from the previous version, where a PSO run is reinitialized with a smaller range of initial search space based on the gbest of the previous PSO run and p 0 that determines the range of initial search space around gbest. This will dictate the new pmax and vmax of the next PSO run. With this approach, fine-tuning the weight of each objective is not necessary for different conditions during measurements. This refinement procedure is described in Algorithm 2.
The first PSO run will optimize for F 1 (estimating the upper triangular matrix K Rm and bias vector k m ) while the rotational calibration matrix K Qm is initiated as identity matrix. This optimizes the refinement process because at the true global minimum, the scalar checking objective is not affected by rotational error. On the other hand, rotational calibration will try to compensate for the misalignment errors that is contained in the triangular matrix of the scalar vmax j = (pmax e, j − pmax s, j )/v lim define the velocity limit 5: for i = 1 . . . n p swarm particle i 6: initialize p i j and v i j using (11) 7: end for 8: end for 9: p i | i=randi(1,np ) = gbest Set one random particle in the swarm equal to previous global best (ensures at least the same fitness value with the previous run) 10: end procedure checking objective. The PSO run is then terminated when the maximum iteration k max is reached: Preliminary tests showed that gbest 2 requires less number of iterations to converge into minimum compared to gbest 1 . In this article, k max for F 1 and F 2 is set to 1000 and 200, respectively. Then, the PSO is reinitiated using the refinement procedure, and the fitness value of the subsequent PSO run is compared to the previous one: If the last fitness value is lower and the difference is smaller than a set threshold F t , then the objective is switched to the next one (from scalar checking to rotation axis fitting, and vice versa). The refinement and objective switching will stop and return the final solution for both objectives when the solutions of both objectives (gbest 1 and gbest 2 ) stagnate and do not compete with each other, i.e., the last refinement procedure does not show improvement larger than the set threshold for both objectives consecutively. The complete sequential objectives refinement procedure is summarized in Fig. 2.

IV. PSO CALIBRATION ALGORITHM VALIDATION
The PSO algorithm is tested with both simulation and real in-orbit flight data. The simulation data will emulate the condition of the real in-orbit data so that the result can be evaluated against the known simulated calibration error. The main evaluation metric for the simulation result is the accuracy of the rotation correction matrix K Qm , which can be calculated directly by converting the rotation matrix into axis-angle representation and taking the angular component [21] where S Qm is the orthogonal matrix from the QRdecomposition of S m (known magnetometer error matrix model; see (1). The accuracy of the scaling and misalignment matrix K Rm has been demonstrated in controlled test environment of the previous paper [13].

A. Calibration of Simulated Data
The simulation will emulate the orbit and attitude dynamics of Aalto-1 and ESTCube-1 in-orbit conditions, i.e., the simulation will replicate the orbit parameters and physical properties of each spacecraft, including the Earth magnetic model. Both satellites orbit around polar LEO, where the magnetic environments are practically identical. As for their physical geometries, Aalto-1 is a three-unit (3U) CubeSat [14], while ESTCube-1 is a one-unit (1U) CubeSat [16], both without any major deployable structures that will affect the attitude dynamics significantly during the time the data was collected. The difference in geometry does affect the rotation mode of the spacecraft as it naturally tumbles in orbit because 3U CubeSats have a distinct minor axis of inertia, while 1U CubeSats are closer to symmetrical in all axes. The simulation will run three sets of simulation model. Initial spin rate set on x-and z-axes, with small random initial rate around the y-axis. 3) Baseline model: The measurement data is collected from simulations with two initial rotation axes in (approximately) spin-stabilized condition, i.e., around the major and minor axes of inertia of a 3U CubeSat (geometrically identical with the 3U model). The measurements from this baseline condition is closer to the conditions of the ground test conducted in [13], as the spacecraft now have more than one distinct (initial) rotation axis. However, in-orbit factors such as the changing magnetic field along the orbit and nutations from spacecraft body dynamics still apply; this makes the baseline model close to ideal condition while still realistically emulating in-orbit data, especially when the spacecraft is capable of controlling its attitude to meet these conditions.
All three sets of simulations are run under varying initial spin rate (from 0 to 180°/s) and noise level, and the sensor data are sampled every 1 s. Noise level (in three standard deviation) for the magnetometer η m ranges from 0 up to 10 μT, while the gyroscope ARW η ARW and RRW η RRW ranges from 0 up to 0.5 and 0.05°/s, respectively. This gyroscope noise is relatively high compared to the typical noise variance of different types of gyroscopes in order to  [20].
The calibration algorithm is applied on all model sets with different combinations of noise levels. Table I shows the combinations of magnetometer and gyroscope noise levels used in the simulation, and the results are as follows. Fig. 3 shows an example of the simulated spin dynamics from one initial spin rate. Sample of 3-D plot of measured, calibrated, and modeled magnetic field vectors from each model with specific initial spin rate and noise level is shown in Fig. 4, while the rotation correction accuracy calculated from (14) for every model is shown in Fig. 5. Fig. 5(d) specifically showcases the rotational accuracy of the previous algorithm from [13] on every model with zero noise (index 1 in Table I). In order to adapt this algorithm to the dynamic rotation axis of the simulated data set, a single locus is defined from a group of every three measurements. The weight for each fitness value is set to a constant value (numbers taken from [13]) across all conditions.

1) Calibration of Aalto-1 In-Orbit Data:
Aalto-1 attitude is mostly in a tumbling state around its major or minor axis [15], [26]. The measurements were gathered at varying sampling period, mostly under 1 s, at an average of ≈ 284 ms, with a total of 5400 measurements that was taken on December 12, 2019 over a single ≈ 25 min session. The calibration algorithm is applied on the in-orbit data during a tumbling state with a spin rate of 24°/s around the z-axis, which is the spacecraft minor axis of inertia, with some   where K Qm for Aalto-1 corresponds to a 8.13°angular correction when converted to axis-angle form. The magnetic field measurement locus before and after calibration are shown in Fig. 8(a).

2) Calibration of ESTCube-1 In-Orbit Data:
ESTCube-1 attitude is largely dominated by the torque from residual magnetic moment of the spacecraft, resulting in its magnetometer reading mostly pointing on one side of the spacecraft. The measurement was taken over several measurement campaigns in May 2014, where several attitude control maneuvers were conducted, resulting in better data set for this calibration method [17], [24], [27]. The sensors sampling period varies from 400 to 700 ms, with a total of 2744 measurements. The gyroscope data can be seen in Fig. 7.
The estimated calibration parameters are where K Qm for ESTCube-1 corresponds to a 2.137°angular correction in axis-angle form. The magnetic field measurement locus before and after calibration are shown in Fig. 8(b).
C. Discussion 1) Simulated Data: The PSO-based calibration algorithm was applied to three different simulated data sets as described in Section IV-A. With the knowledge of the error models of the magnetometer, the calibration results can be validated against the reference model by examining the rotation correction errorθ , and this reveals several findings on the algorithm technical performance. a) Performance of the previous PSO algorithm. Fig. 5(d), compared to other figures in Fig. 5, shows that the performance of the previous PSO algorithm, even in noise-free data, is not as robust as the improved algorithm. The main cause of this is the competition between the two optimization objectives, which is affected by the selection of fitness function's weight for each objective in conjunction with different convergence rate from the initial estimates of the swarm. In this case, the previous algorithm has smaller window in initial spin rate for optimal convergence. For the 3U spin-stabilized model, the lack of rotation axis variation narrows this optimal window even further. Other factors such as noise will also affect the optimal weights for the fitness values. In a real data scenario where the "true" calibration parameters are not known, evaluation of the rotational correction element accuracy might be difficult, as in some cases of the simulations, the estimate for K Rm and k m still converges to global optimum while K Qm converges very far from global optimum.
b) In-depth performance of the improved PSO algorithm. The biggest advantage for the improved algorithm is the absence of a separate weight for each fitness value that needs to be manually fine-tuned. However, the improved algorithm is still affected by many factors such as spin rate, sampling rate relative to spin rate, magnetometer and reference rotation axis (gyroscope, in this article) noise, and the spin mode of the spacecraft.
Magnetometer data sampling rate in conjunction with the type of spin mode will affectθ on higher spin rate. This is in line with the requirement for a single continuous locus due to the rotation axis estimation in Algorithm 1, where the upper spin rate limit is half the sampling rate, i.e., 180°/s spin rate upper limit in the case of the simulated 1 Hz sampling rate [shown in Fig. 5(c) at 185°/s]. Measurement noise blurs this limit as the noise affects the 180 • angular displacement limit between the magnetic field vector.
Rotation correction errorθ will be larger for very low spin rate. This is because the magnetic field direction changes along the orbit, which will affect the rotation axis estimation in Algorithm 1 by adding rotation components that are uncorrelated to the spacecraft actual spin axis. While in noise-free condition the accuracy error is hitting minimum from around 40 to 50°/s spin rate, the presence of noise in general is a more dominant factor and the error already reaches minimum at a spin rate of 10-30°/s with higher magnitude of error.
As expected, the theoretically "ideal" condition of the baseline model [ Fig. 5(c)] gives the most consistent performance in minimizingθ across the spin rate limit.θ can be kept under 0.1 • magnitude until the highest gyroscope noise level (noise index 4 and 5 in Table I) is applied, where the gyroscope noise is more dominant than the magnetometer noise (except when the magnetometer noise is set very high as in index 6) and only higher spin rate can mask this presence of noise.
Spin mode with mostly one rotation axis with small nutation (3U spin-stabilized model) shows the weakest results in terms of calibration performance. This is expected, as more than one distinct rotation axis (ideally orthogonal to each other) is needed to resolve rotational ambiguity in all axes. The small nutation rate and the change of magnetic field direction along the orbit actually helps in adding variation to the available rotation axis [resulting in the locus coverage in Fig. 4(a)], but this model is very sensitive to noise: The difference inθ between different noise levels is much less pronounced compared to the 1U high nutation model and baseline models, whereθ is at several orders of magnitude higher even at low noise. Error in the order of 1°-10°can be expected even at low to moderate noise (noise index 2-3 in Table I).
Spin mode with more pronounced nutation shows earlier rise ofθ before the spin rate reaches the upper limit: The 1U high nutation model error starts to rise an order of magnitude at ≈50°/s initial spin rate, while 3U spin-stabilized and baseline model only starts to behave similarly at 100-120°/s initial spin rate (see Fig. 5). This is caused by the higher nutation rate [e.g., the oscillation of the spin rate around x-axis in Fig. 4(b)], where the estimated rotation axis that was averaged over three measurement points is undersampled compared to the speed the spin rate is changing.
2) In-Orbit Data: With this knowledge on algorithm performance with different simulated models, the calibration of real in-orbit flight data of Aalto-1 and ESTCube-1 can be evaluated.
The attitude mode of Aalto-1 is very similar to the 3U spin-stabilized model, although Aalto-1 has a higher nutation rate [see Fig. 6 compared to Fig. 3(a)]. Combined with the very high magnetometer noise level in Aalto-1 data (between noise index 3 and 6 in simulated model), the rotation correction accuracy of the calibration is likely in the order of ∼10°. The rotation correction factor of Aalto-1 calibration is around 8.13°, which is under the expected uncertainty itself. This means that there is no major misalignment between the magnetometer axes and the reference axes (in this case, the gyroscope axes), but no accuracy improvement can be expected from the magnetometer data besides the obvious bias correction.
ESTCube-1 in-orbit data, on the other hand, does not exactly fit any of the specific simulated models; while geometrically a 1U CubeSat, ESTCube-1 is affected by high residual magnetic moment and thus "locked" to the Earth magnetic field, resulting in most of the data showing a spin rate around one axis-visible in the mostly parallel measurement locus planes in Fig. 8(b). However, the data does include some variations in the rotation axis due to the attitude control attempt around the time of measurement, especially in the first part of the gyroscope reading in Fig. 7.
This variation in rotation axis is important as this allows a much better performance compared to the spin-stabilized model-closer to the high nutation or baseline model. Since ESTCube-1 also has relatively slow nutation rate compared to its data sampling rate, the rotation correction accuracy can be expected in the 0.1°-1°order. The rotation correction factor of ESTCube-1 calibration result is ≈2.14°, which shows possible misalignment between the magnetometer axes with the reference gyroscope axes.

V. CONCLUSION
An enhanced PSO-based magnetometer calibration using scalar checking and rotation axis fitting has been developed. This is based on an improvement of the rotation axis fitting method and a novel approach in the integration of the two optimization objectives using a sequential objectives refinement process. These improvements allow the algorithm to process measurements taken under different spacecraft dynamics without the restriction of a clearly defined rotation axis and the need to manually fine-tune the objectives' weight of the previous version [13].
Using sensors data with simulated in-orbit spacecraft dynamics, the accuracy of the rotation correction factor from the calibration result can be validated. In best cases, rotation accuracy in the order of ∼0.1°between the magnetometer and the reference sensor (in this case, gyroscope) can be achieved under reasonable noise level (up to noise index 4 in Table I). To achieve this, there are three main factors that determine the quality of the data for calibration: rotation mode of the spacecraft, spacecraft spin rate, and sampling rate of the sensors (magnetometer and the sensor for reference rotation axis).
The rotation mode and spin rate should ideally result in varying rotation axis during measurements, which can be achieved either from a significant angle of nutation or using attitude control to change the rotation axis. The sensor sampling rate defines the range of spin rate in which the calibration will perform best. Generally, the sampling rate lower limit is twice the rotation rate, i.e., a sampling rate of every 1 s if the spacecraft fully rotates every 2 s. For spacecraft dynamics with high nutation angle and rate, the sampling rate also needs to be higher than the oscillation of the spin rate. Faster sampling rate than the lower limit (or higher spin rate) is also needed to improve accuracy, especially in the presence of noise, as the natural change of magnetic field direction along the orbit lowers the accuracy in lower spin rate.
Applying the calibration on real in-orbit data of Aalto-1 and ESTCube-1 shows that the measurements do not always fulfill the ideal conditions desired from the spacecraft attitude dynamics. If a higher accuracy from the magnetometer is required, dedicated calibration session that actively utilizes the spacecraft attitude control can be used. This also depends on other uncertainties associated with the sensor used as the reference rotation axis for the rotation axis fitting-this article only focuses on random noise typically modeled for gyroscopes. Further improvements on the calibration process as a whole should take into account other error sources for the magnetometer and the calibration reference sources. Magnetometer model could include temperature dependency or time-varying magnetic moment from the spacecraft electronics. Reference for the scalar checking objective might contain error from the orbit estimation in combination with the Earth magnetic field model. The reference sensor used for rotation axis fitting could also include other uncertainty models.