An orbit determination algorithm for small satellites based on the magnitude of the earth magnetic field

Autonomous attitude determination systems based on simple measurements of vector quantities such as magnetic field and the Sun direction are commonly used in very small satellites. However, those systems always require knowledge of the satellite position. This information can be either propagated from orbital elements periodically uplinked from the ground station or measured onboard by dedicated global positioning system (GPS) receiver. The former solution sacrifices satellite autonomy while the latter requires additional sensors which may represent a significant part of mass, volume, and power budget in case of pico- or nanosatellites. Hence, it is thought that a system for onboard satellite position determination without resorting to GPS receivers would be useful. In this paper, a novel algorithm for determining the satellite orbit semimajor-axis is presented. The methods exploit only the magnitude of the Earth magnetic field recorded onboard by magnetometers. This represents the first step toward an extended algorithm that can determine all orbital elements of the satellite. The method is validated by numerical analysis and real magnetic field measurements.

Autonomous attitude determination systems based on simple measurements of vector quantities such as magnetic ¦eld and the Sun direction are commonly used in very small satellites. However, those systems always require knowledge of the satellite position. This information can be either propagated from orbital elements periodically uplinked from the ground station or measured onboard by dedicated global positioning system (GPS) receiver. The former solution sacri¦ces satellite autonomy while the latter requires additional sensors which may represent a signi¦cant part of mass, volume, and power budget in case of pico-or nanosatellites. Hence, it is thought that a system for onboard satellite position determination without resorting to GPS receivers would be useful. In this paper, a novel algorithm for determining the satellite orbit semimajor-axis is presented. The methods exploit only the magnitude of the Earth magnetic ¦eld recorded onboard by magnetometers. This represents the ¦rst step toward an extended algorithm that can determine all orbital elements of the satellite. The method is validated by numerical analysis and real magnetic ¦eld measurements. at a¨ordable costs. Furthermore, nano-or picosatellites can be also used in scienti¦c swarm missions (see, for example, [1]) where the utilization of small platforms (see, for example, [2]) allows to launch multiple satellites with a single launch vehicle. Nonetheless, limitations coming from the small size pose a signi¦cant engineering challenge which may open new opportunities for reevaluation of a©rmed technologies performing essential tasks in larger satellites. For instance, the determination of the satellite attitude is usually estimated by measuring the Sun direction and magnetic ¦eld vector and comparing the measurements with geomagnetic and Sun position models. In turn, these models need the knowledge of the satellite£s position. Therefore, the determination of the satellite£s position becomes a requirement for attitude determination. Methods for measuring the spacecraft position include the utilization of specialized space-rated GPS receivers. Despite of the miniaturization and improving e©ciency of such devices, GPS receivers may still represent a signi¦cant part of mass, volume, and power budgets for nano-and picosatellites (see, for example, [37]). Another way for determining the satellite£s position consists of using ground or space-based measurements such as optical and radar tracking or laser ranging. The orbital elements are then calculated and uplinked to the satellite via a ground station network from time to time. Between two succeeding communications, the satellite position is calculated onboard through an orbital model. Main drawbacks of this solution are the lack of autonomy of the system and the progressive growth of the position estimate error.
From the discussion above, it turns out that it would be advantageous to develop a system for determining the position of the satellite which is both autonomous and light in terms of mass, volume, and energy. This work intends to be an initial step toward this goal proposing a novel approach for determining the semimajor axis of low-Earth orbits (LEO) class based only on the frequency analysis of the magnitude of the magnetic ¦eld measured onboard.
The paper is structured as follows. In section 2, theoretical principles of the algorithm are presented. Section 3 shows the performances of the method when applied to a simulated environment and to real magnetic ¦eld data measured by a SWARM satellite. The paper concludes with a short summary and gives perspectives for future work.

ORBITAL SEMIMAJOR AXIS DETERMINATION PROCEDURE
The idea behind the method is that the signal B(t) of the magnetic ¦eld magnitude measured by the satellite orbiting around the Earth presents clear periodicity which can bring information on the satellite orbital elements. Hence, it is expected that a Fourier transform of B(t) can lead to the calculation of the satellite orbiting frequency, f sat . Eventually, the semimajor axis of the satel-lite£s orbit, a, can be calculated from the satellite period through the well-known formula: In order to motivate the procedure, it is useful to recall here the expression of IGRF (International Geomagnetic Reference Field) magnetic ¦eld model that describes the magnetic ¦eld around the Earth. The model implements a spherical harmonic expression for the magnetic ¦eld potential V as follows: where R E is the reference radius of the Earth; r, θ, and φ are the geocentric coordinates, namely, altitude, latitude, and longitude, respectively; g m n and h m n are the coe©cients; and P m n are the Schmidt quasi-normalized associated Legendre functions. Since the magnetic ¦eld B(t) is directly calculated from the potential, an understanding of the potential spectrum gives much insight into the magnetic ¦eld as well. Here, the key point is to give an expression for θ and φ in terms of frequencies of the Earth and the orbiting satellite where θ and φ are to be intended as the geometric coordinates of the satellite with respect to a frame ¦xed to the Earth (Earth-centered rotational, ECR). In order to better justify the results in the following subsections, two particular orbits with inclinations of 0 • and 90 • are illustrated in the beginning. Afterwards, the general case with di¨erent inclination is presented.
Here, only two simple orbits will be studied in detail. First, a circular orbit with 0 degree inclination will be analyzed. Then, a circular orbit with 90 degree inclination will be considered. Eventually, short considerations are provided to justify the validity of the analysis in the general case.

Orbit with Inclination of 0 •
Since the magnetic ¦eld is connected with the Earth, onboard magnetic ¦eld measurements must refer to the relative motion of the satellite with respect to the Earth. For this case, it is not di©cult to prove that the angular velocity of the satellite with respect to the Earth is a constant vector parallel to the Earth rotation axis (Fig. 1) and with magnitude where ω sat and ω E are the magnitudes of the angular velocities of the satellite and the Earth, respectively. Since ω sat and ω E are constant, it is possible to introduce the frequencies f sat = ω sat /(2π) and f E = ω E /(2π). Finally, one can link the latitude and longitude of the relative motion of the satellite with respect to the Earth to f sat and f E as follows: Thus, the fundamental harmonics for θ and φ are: These results show that the potential in Eq. (2) is only function of the longitude. Moreover, its expression appears to be a linear combination of sinusoidal Therefore, the Fourier transform of B(t) is expected to presents peaks only at these frequencies. This behavior is pointed out in the |F (B(t))| spectrum of Fig. 2 where the Fast Fourier Transform (FFT) of B(t) is presented along with the original B(t) signal. Thus, the determination of the satellite frequency can be done by assessing f sat/E from the Fourier spectrum and summing f E to f sat/E . Eventually, the orbit radius can be calculated by Eq. (1).

Orbit with Inclination of 90 •
For a Keplerian circular orbit with 90 degree inclination, the analysis may proceed in a di¨erent manner. Referring to Fig. 3, parametric expressions for latitude and longitude take the form: and, consequently, Considering these results, ONE can conceive the FFT of the potential in Eq. (2) as a series of amplitude modulations where the carrier signals are the harmonics with frequency multiple of f θ and the modulating signal is composed of powers of harmonics with frequencies multiple of f φ . Thus, the B(t) spectrum

Orbit with Inclination Di¨erent from 0 • and 90 •
Due to the lack of orthogonality between ω sat and ω E , a visualization of the relative motion of the satellite is not straightforward. However, starting from the description of the motion with respect to Earth-centered inertial (ECI) coordinates and making a coordinate transformation, it is possible to write the parametric expressions of the satellite motion in ECR coordinates as follows: r x ′ = r(cos ω E t cos ω sat t + sin ω E t sin ω sat t cos i) ; r y ′ = r(− sin ω E t cos ω sat t + cos ω E t sin ω sat t cos i) ; r z ′ = r sin ω sat t sin i .
In the equations above, one can notice that r x ′ and r y ′ , which determine the satellite longitude, are the periodic functions with fundamental angular frequency ω sat − ω E while r z ′ , which determines the latitude, has an angular frequency equal to ω sat . Thus, the frequencies for latitude and longitude are: Using the same line of arguments as led to the explanation of the spectrum for the orbit with inclination of 90 • , it is possible to explain the frequency contents of the B(t) spectrum for this general case, too. The B(t) spectrum for an orbit with inclination of 30 • is given in Fig. 5. It is worth to notice that the two cases with inclinations of 0 • and 90 • can be also seen as particular cases of the general case presented here. Furthermore, Eqs. (3) are valid for eccentricity e di¨erent from 0, too. Indeed, the latitude and longitude periodicity are not the functions of e.

Procedure
The possibility to identify the satellite frequency from the B(t) spectrum enables one to develop a procedure for semimajor axis determination based on the analysis of the magnetic ¦eld magnitude. The procedure is shortly described in the following steps.
1. Record the B(t) signal. The signal length must be long enough to have a su©cient frequency resolution and to capture the lower frequency f E . It has been observed that an acquisition time of about 80 satellite orbits yields good results. Increasing the acquisition time does not yield signi¦cant improvements in the accuracy. On the other hand, for acquisition time lower than 70 orbits, the method considerably looses accuracy.
2. Multiply the time signal by a Blackman window and complete the missing samples to the nearest upper 2 N being number by 0-padding (with N a natural number). The sampling frequency is not a¨ecting the accuracy of the result in a sensible manner. In this work, a frequency of 0.01 Hz has been considered. 6. Determine the points de¦ning each lobe by a procedure which checks the curve gradient.
7. Calculate the weighted average frequencies of each lobe as follows: where J is the number of points forming the lobe; f j is the frequency of the jth point of the lobe; and B j is the magnitude of F (B(t)) at the frequency f j . For each lobe, calculate also the power p = J l j=1 B j .
8. Among all chosen lobes, select the lobe with the largest power. From this lobe, f sat is estimated by recalling its position in the spectrum.
9. Apply Eq. (1) for semimajor axis assessment. Figure 6 The set of selected lobes with datapoints marked with ¦lled and blank dots.
The lobe marked with black dots has the largest power and is used for determination of fsat The application of the procedure to the orbit presented in the previous subsection (see orbital data in Fig. 5) is presented in Fig. 6. Here, the set of selected lobes and their corresponding J points are marked with blank dots. The lobe with the largest power is the one with black markers and it is used for the determination of f sat .

APPLICATION CASES
The algorithm described in subsection 2.4 uses magnitudes of the magnetic ¦eld. For testing the algorithm performances, a simulating framework has been set up. The framework employs Scilab [8] with CelestLab library developed by the Centre National d£Etudes Spatiales (the French Space Agency) for simulating satellites orbital movement. CelestLab toolbox implements the IGRF11 geomagnetic ¦eld model. Orbit environment data generated in Scilab are exported to Matlab where the determination algorithm is implemented. Figure 7 presents the diagram of the simulation framework adopted in the work.
The algorithm has also been tested in a real world scenario. The magnetic ¦eld measurements gathered by SWARM [9] satellite constellation have been used for this purpose. SWARM mission was launched by the European Space Agency (ESA) to advance human understanding of how Earth magnetic ¦eld works. The three satellites (SWARM-A, SWARM-B, and SWARM-C) precisely measure the magnetic signals that stem from Earth£s core, mantle, crust, and oceans as well as its ionosphere and magnetosphere. Because of the very precise magnetometers present onboard, data gathered by SWARM are excellent for validating the orbit determination algorithm. The orbit of SWARM-A is almost circular with eccentricity 0.0003, inclination 87.4 • , and mean semimajor axis 6853 km. The magnitudes of the magnetic ¦eld measured by onboard magnetic sensor of about 5 days have been analyzed in this work. Since data provided by ESA are in common data format ¦les, it was necessary to develop a dedicated parser in Matlab to translate data into a magnetic ¦eld magnitude series. To verify the correctness of the extracted data, these latter have been compared with the ¦eld predicted by the IGRF model implemented in Scilab. The high accuracy of the model predictions is highlighted in the plots of Fig. 8 where the modeled and measured B(t) signals for a short time interval (Fig. 8a) and the error between the two signals for the whole available time (Fig. 8b) is given. Visual inspection suggests that the errors are few orders of magnitude smaller than the signal itself.

Simulated Data
The performances of the proposed procedure have been tested on a set of 110 randomly sampled orbits. The assumed ranges of variations for orbital elements for the 110 orbits are given in Table 1.  Figure 9 Algorithm accuracy resultant from 110 testing samples The procedure described in subsection 2.4 is applied to all orbits and the estimated semimajor axis is compared with the correct one. The accuracy of the estimated semimajor axes by the histogram is presented in Fig. 9. The procedure accuracy is remarkably within 0.2 km. Also, no signi¦cant correlation has been found between the error and orbital parameters.

SWARM Satellite
In contrast to the simulation data using Keplerian orbits without perturbations, SWARM satellites are a¨ected by orbit perturbations due to several sources. The most important one is the Earth oblateness whose two most relevant e¨ects are the nodal precession and the advance of perigee. Secular J 2 e¨ects are accounted for in the procedure by correcting f sat estimates with the following relationships: In Eqs. (4), ω nod , ω aps , and ω mna are the average rates of variations of the perturbed orbital elements calculated by the method of averaging [10]; e is the eccentricity; R E is the mean Earth radius; J 2 is the oblateness factor; and f sat is the estimated satellite frequency. The details of the enhanced procedure for secular orbits are skipped here and will be presented in future publications. Using real orbit eccentricity and inclination values yields a semimajor axis estimate of 716 m. Without correction, the estimate error is about 7 km.

CONCLUDING REMARKS
In this work, a novel procedure has been proposed for estimating the semimajor axis of a satellite orbit resorting only to the information of magnetic ¦eld magnitude measurements without the support of any additional analytic model. The method yields estimates of the semimajor axis with accuracy below 200 m for Keplerian orbits. The identi¦cation method has been also tested on real mission data. Estimate accuracy drops but still satisfying, with an error of less than 1 km. More work is to be done for extending the algorithm capability to the determination of inclination and eccentricity.