Analytical solution of the spacecraft magnetometer calibration problem using the method of least squares
- 作者: Kirillov K.A.1, Kirillova S.V.1, Kuznetsov A.A.1, Melent'ev D.O.2,3, Safonov K.V.1
-
隶属关系:
- Reshetnev Siberian State University of Science and Technology
- Siberian Federal University
- JSC “Information Satellite Systems” Academician M. F. Reshetnev Company”
- 期: 卷 26, 编号 2 (2025)
- 页面: 181-194
- 栏目: Section 1. Computer Science, Computer Engineering and Management
- ##submission.datePublished##: 30.06.2025
- URL: https://journals.eco-vector.com/2712-8970/article/view/686523
- DOI: https://doi.org/10.31772/2712-8970-2025-26-2-181-194
- ID: 686523
如何引用文章
详细
In this paper, an analytical method is proposed for solving the problem of magnetometer calibration for the model considered in [1]. When solving the problem of determining the calibration parameters of the magnetometer unit, it is taken into account that for measurements with any spatial orientation of the magnetometer unit, the value of the measured magnetic induction vector is preserved and is a known model value. A penalty function of 12 variables equal to the sum of the squares of the residuals is introduced into consideration. The algorithm for solving the problem of calibrating the measuring axes of the magnetometer unit is reduced to searching, by the method of least squares, for such values of the variables of this function that, for a given set of magnetometer measurement vectors, provide it with a minimum. For this purpose, the specified function is examined for a conditional extremum in the presence of three equality constraints. The Lagrangian function is compiled and, based on the necessary condition for the extremum of this function, the system of 15 equations in the 15 variables is formed. It is proved that the system has four solutions. Formulas are derived that make it possible to obtain the components of each of these four solutions. As a solution to the magnetometer calibration problem, it is recommended to choose a solution to the specified system that provides a minimum of the Lagrangian function.
全文:
Introduction
Magnetometers are part of the orientation and stabilization system of low-orbit small-sized spacecraft. They are widely used due to the fact that they are lightweight, inexpensive and reliable. However, due to the physical properties of the sensitive element, modern magnetometers require mathematical calibration. At present, various methods for calibrating magnetometers have been proposed, and a considerable number of scientific papers are devoted to these methods [1–12], in particular, article [9], which provides an overview of various methods for performing such operations.
In the above-mentioned works, the problem of calibrating the magnetometer of a spacecraft was solved using numerical methods. In this article, an analytical method for solving this problem for the model considered in [1] is proposed.
1. Model of errors in measurements of magnetic induction vector
We denote by h the value of the measured magnetic induction vector at a certain spatial position of the magnetometer unit (MU). We use the simplified measurement model considered in [1]:
(1)
The following notations [1] are used in(1):
B = (B1, B2, B3)T is a true magnetic induction vector;
b = (b1, b2, b3)T is a constant vector corresponding to zero offsets for each of the measuring axes of the MU;
n is a random vector corresponding to uncorrelated noise for each of the measurement axes;
Р is a matrix, the rows of which are the unit vectors of the measuring axes of the MU, written in the “base” MU system;
Q is a diagonal matrix containing on the main diagonal the scaling factors k1, k2, k3 for measuring axes MU i.e
.
In other words, the matrix P describes the non-orthogonality of the MU measurement axes, and the matrix Q corresponds to scaling along these axes.
The task of calibrating the MU measuring axes comes down to finding the elements of the matrices P and Q, as well as the zero offset vector b.
2. Development of an algorithm for determining the calibration parameters of the MU
When solving the problem of MU determining the calibration parameters, we will use the fact that for measurements with any MU spatial orientation, the value of the measured magnetic induction vector B is preserved and is a known model value.
We will assume that in flight, as a result of magnetometer measurements at discrete moments in time, a set of vectors is obtained h(l) = (h1(l), h2(l), h3(l))T, l = 1, …, N. Provided there is no measurement noise, from formula (1) we obtain:
(2)
where B(l) = (B1(l), B2(l), B3(l))T is the true vector of magnetic induction at the same point in space as the measured vector h(l), l = 1, …, N. We express vectors (l = 1, …, N) from equalities (2):
(3)
l = 1, …, N, где
(4)
As noted in [1], without loss of generality, the non-orthogonality matrix P can be represented with a minimum number of unknown elements as follows:
where ε1, ε2, ε3 are small angles. Then the inverse matrix P – 1 will take the form
where αi = tg εi, βi = sec εi, i = 1, 2, 3, and since
then in accordance with (4) we obtain:
(5)
where γi = ki – 1, i = 1, 2, 3.
We rewrite equality (3) taking into account (5):
(6)
l = 1, …, N. We write each of the N vector equalities (6) as a system of three scalar equalities:
(7)
l = 1, …, N. Due to the first equality of system (7), we rewrite the second equality of this system in the form
as a result of which the third equality of system (7) is written as follows:
l = 1, …, N.
We introduce into consideration a function of twelve variables αi, βi, γi, bi (i = 1, 2, 3)
where the variables αi и βi are related by the trigonometric identity
by the equalities
(8)
i = 1, 2, 3.
The algorithm for solving the problem of calibrating the MU measuring axes is reduced to searching by the least squares method [13] taking into account (8) such values of variables αi, βi, γi, bi (i = 1, 2, 3), which, for a given set of measurement vectors, {h(l)} (l = 1, …, N) provide a minimum of the function F. For this purpose, it is necessary to investigate the function F for a conditional extremum [14] in the presence of three constraint equations (8).
We compose the Lagrange function
, (9)
fifteen variables αi, βi, γi, bi, λi (i = 1, 2, 3) and we write down the necessary condition for the local extremum of this function:
It is required to find the stationary points of the function F, i.e. the solution of the system of equations (10). Let us introduce the following notations:
Taking into account these notations, for convenience we will divide the system of equations (10) into three nonlinear systems – a system of equations for two unknowns b1, γ1:
(11)
a system of equations for five unknowns b2, α1, β1, γ2, λ1:
(12)
and a system of equations for eight unknowns. b3, α2, α3, β2, β3, γ3, λ2, λ3:
(13)
неравенства Коши
(13)
In writing the expressions βi through αi in the last equations of systems (12) and (13), we used the fact that βi = sec εi > 0 due to the obvious inequalities–π/2 < εi < π/2, i = 1, 2, 3. We solve each of the systems of equations (11), (12), (13).
First, we note that the inequality is valid
(14)
i = 1, 2, 3. Indeed, the discriminant of the quadratic equation is – Fi + 2Hibi – Nbi2 = 0 with respect to the unknown bi, being equal to
is negative due to the obvious consequence
of the Cauchy–Bunyakovsky inequality [15] (here the sign “<” is used instead of “≤”, since at least one of the terms hi(1), …, hi(N) is obviously not equal to zero), therefore the quadratic equation under consideration has no real roots, i = 1, 2, 3.
We solve the system of equations (11). Taking into account inequality (14), we express γ1 through b1 from the first equation of this system:
Substituting the obtained expression for γ1 into the second equation of system (11) and taking into account the inequality γ1 = k1– 1 ≠ 0, after simple transformations we arrive at an equality,
taking into account which the equality written above for γ1 is reduced to the form
Thus, the solution to the system of equations (11) is found.
We consider the system of equations (12). Since γ2 = k2 – 1 ≠ 0, the second equation of this system can be written as
(15)
since β1 = sec ε1 ≠ 0, the third equation of system (12) is equivalent to the following equation:
(16)
From (15), taking into account (16), we obtain that λ1β1γ2– 1 = 0 and then, by virtue of the inequality β1 ≠ 0, we have
λ1 = 0 (17)
The inequality D12 – C1b2 ≠ 0 holds, since a direct check shows that the value b2 = D12C1– 1, which is the root of the equation D12 – C1b2 = 0, does not satisfy the equations of system (12)
First, we consider the case H2 – Nb2 ≠ 0, i.e. b2 ≠ H2N– 1. We express β1γ2 from (16), as well as from the first (taking into account equality (17)) and fourth (taking into account the inequality β1γ2 ≠ 0) equations of system (12):
(18)
As proven, the denominator of the fraction in the second of the equalities (18) is not equal to zero, the denominator of the fraction in the first of the equalities (18) is different from zero by virtue of (14). From the first and third equalities (18) follows the equation
we express b2 through α1:
(19)
Similarly, from the second and third equalities (18) follows the equation
we express b2 through α1:
(20)
Equating the right-hand sides of equalities (19) and (20), after simple transformations we obtain a quadratic equation with respect to α1:
(21)
The discriminant of equation (21) is equal to
and the roots of this equation are calculated using the formulae
(22)
We note that the values (22) of the unknown α1 are not the roots of the denominators of the fractions in equalities (19) and (20). Substituting the first of the values (22) of the unknown α1 into the fourth equation of system (12), taking into account the inequality β1γ2 ≠ 0, we arrive at the equality H2 – Nb2 = 0, which contradicts the case under consideration.
Now H2 – Nb2 = 0, i.e. b2 = H2N– 1. Taking into account the inequality β1γ2 ≠ 0, from the fourth equation of system (12) we obtain that α1 = – C2C1– 1, i.e. the found value of α1 coincides with the first of the roots (22) of equation (21). It is easy to show that equalities (19) and (20), proven under the assumption b2 ≠ H2N– 1, are also valid forb2 = H2N– 1. In the case b2 = H2N– 1, the first two equalities from (18) also hold, since no restrictions on the value of the unknown b2 were taken into account when deriving them.
From the above reasoning it follows that the system of equations (12) has two solutions. The components of both solutions of the specified system are obtained in the following order:
- – the values of the unknown α1 are found by formulae (22);
- – the values of β1 corresponding to the found values of α1 are calculated by the formula written as the last equation of the system (12);
- – the values of b2 – by any of the formulae (19), (20);
- – the values of γ2 – by any of the two formulae
which follow respectively from the first and second equalities (18);
– the value of the unknown λ1 in both solutions of the system of equations (12) due to (17) is taken to be equal to zero.
We move on to finding solutions to the system of equations (13). Due to the inequalities γ3 = k3– 1 ≠ 0, β2 = sec ε2 ≠ 0 and β3 = sec ε3 ≠ 0 the second and third equations of this system can be written as
(23)
(24)
respectively. From (23), taking into account (24), we obtain: λ2β2(β3γ3)–1 = 0. Consequently, the equality is valid
(25)
Then the first equation of system (13), taking into account the inequalityβ3 ≠ 0, can be written in the form
(26)
From the sixth equation of system (13), taking into account (24) and (26), we find – λ3β3 = 0, from which the equality follows
(27)
The inequalities H3 – Nb3 ≠ 0, Di3 – Cib3 ≠ 0 take place, since a direct check shows that neither the root b3 = H3N– 1 of the equation H3 – Nb3 = 0, nor the roots b3 = Di3Ci– 1 of the equations Di3 – Cib3 = 0 satisfy the equations of system (13), i = 1, 2.
Taking into account equalities (25), (27) and inequalities β2 ≠ 0, β3 ≠ 0, γ3 ≠ 0, we express β2β3γ3 from the first, third, fourth and fifth equations of system (13):
(28)
As proven, the denominators of the fractions in the first, third and fourth equalities (28) are not equal to zero, the denominator of the fraction in the second equalities (28) is different from zero by virtue of (14). From the first and third equalities (28) follows the equation
from the second and third equalities (28) follows the equation
from the third and thorthequalities (28) follows the equation
From the last three equations we express α2β3 through α3 and b3:
(29)
(30)
(31)
By equating the right-hand side of equality (29) to the right-hand sides of equalities (30) and (31), after simple transformations we obtain the equations
from which we express α3 through b3, having previously divided both parts of each of them by (H3 – Nb3):
(32)
(33)
Having equated the right-hand sides of equalities (32) and (33), after the simplest transformations we arrive at a quadratic equation with respect to b3
(34)
the coefficients of which are determined as follows:
The discriminant of equation (34) is equal to
and the roots of this equation are calculated using the formulae
(35)
Since equation (34) has two solutions, the system of equations (13) also has two solutions. The components of both solutions of this system are obtained in the following order:
- – the values of the unknown b3 are found using formulae (35);
- – the values of the unknown α3 corresponding to the found values of b3 are calculated using any of the formulae (33), (34);
- – the values of β3 – using the formula written as the last equation of system (13);
- – the values of α2 – using any of the three formulae
which follow from equalities (29), (30) and (31),
respectively;
- – the values of β2 – by the formula written as the last equation of system (13);
- – the values of γ3 – by any of the four formulae
which follow from (28);
- – the values of the unknowns λ2 and λ3 in both solutions of the system of equations (13) are taken to be equal to zero due to equalities (25) and (27)
Thus, the system of equations (11) has a single solution, and each of the systems (12), (13) has two solutions. Consequently, the number of solutions of the original system of equations (10), defined as the product of the number of solutions of systems (11), (12) and (13), is four. Of the four solutions obtained,
(36)
(j = 1, 2, 3, 4) the systems (10), which are also stationary points of the Lagrange function F, defined by equality (9), we are interested in the solution
satisfying the equality
Since in each of the solutions (36) of the system of equations (10) the last three components are equal to zero, i.e.
j = 1, 2, 3, 4, by virtue of (9) we have:
j = 1, 2, 3, 4, therefore, we determine the desired solution
using the formula
The required values of the scale factors for the measuring axes of the MU are calculated using the formulae
i = 1, 2, 3, required angle values are calculated using the formulae
which follow from the equalities
and obvious inequalities i = 1, 2, 3.
Conclusion
Thus, we have obtained an analytical solution to the problem of calibrating the magnetometer of a spacecraft for the model considered in [1]. The procedure for calculating the calibration parameters of the MU using the derived formulas has a number of obvious advantages over numerical methods for solving this problem:
- – the number of arithmetic operations is significantly reduced;
- – the problem of possible instability of the method disappears.
作者简介
Kirill Kirillov
Reshetnev Siberian State University of Science and Technology
编辑信件的主要联系方式.
Email: kkirillow@yandex.ru
ORCID iD: 0000-0002-3763-1303
Dr. Sc. (Phys. and Math.), Associate Professor, Professor of the Department of Applied Mathematics
俄罗斯联邦, 31, Krasnoyarskii Rabochii prospekt, Krasnoyarsk, 660037Svetlana Kirillova
Reshetnev Siberian State University of Science and Technology
Email: svkirillova2009@yandex.ru
ORCID iD: 0000-0003-3779-2825
Cand. Sc. (Technical Sciences), Associate Professor, Associate Professor of the Department of Applied Mathematics and Data Analysis
俄罗斯联邦, 31, Krasnoyarskii Rabochii prospekt, Krasnoyarsk, 660037Alexander Kuznetsov
Reshetnev Siberian State University of Science and Technology
Email: alex_kuznetsov80@mail.ru
ORCID iD: 0000-0003-0944-1817
Dr. Sc., Professor, Head of the Institute of Space Research and High Technologies
俄罗斯联邦, 31, Krasnoyarskii Rabochii prospekt, Krasnoyarsk, 660037Denis Melent'ev
Siberian Federal University; JSC “Information Satellite Systems” Academician M. F. Reshetnev Company”
Email: denes.2000@mail.ru
ORCID iD: 0009-0009-6187-4098
Graduate Student, Engineer
俄罗斯联邦, 79, Svobodny Av., Krasnoyarsk, 660041; 52, Lenin St., Zheleznogorsk, Krasnoyarsk region, 662972Konstantin Safonov
Reshetnev Siberian State University of Science and Technology
Email: safonovkv@rambler.ru
Dr. Sc. (Phys. and Math.), Professor, Head of the Institute of Informatics and Telecommunications
俄罗斯联邦, 31, Krasnoyarskii Rabochii prospekt, Krasnoyarsk, 660037参考
- Rozin P. E., Simonov A. V., Gordienko E. S., Zaiko Yu. K. [In-Flight Calibration of the “Dekart” Cubesat Magnetometer]. Trudy MAI. 2022, No. 124, P. 1–20 (In Russ.).
- Akimov I. O., Ilyukhin S. N., Ivlev N. A., Kolosov G. E. [Magnetometer calibration technique for the ground-based stage of spacecraft system diagnostics]. Inzhenernyy zhurnal: nauka i innovatsii. 2018, Vol. 77, No. 5, P. 1–18 (In Russ.).
- Morskoy I. M., Simonov A. V., Lyaskovskaya V. I., Ezhov A. S. [Ballistic support for the development and flights of the interorbital space tug “Fregat”]. Vestnik “NPO im. S. A. Lavochkina”. 2014, No. 1, P. 10–15 (In Russ.).
- Alonso R., Shuster M. D.: TWOSTEP: A Fast Robust Algorithm for Attitude-Independent Magnetometer Bias Determination. Journal of the Astronautical Sciences. 2002, Vol. 50, No. 4, P. 433–451.
- Belsten N., Knapp M., Masterson R., Payne C., Ammons K., Lind F.D., Cahoy K. Verification and calibration of a commercial anisotropic magnetoresistive magnetometer by multivariate non-linear regression. Geosci. Instrum. Method. Data Syst. 2023, No. 12, P. 201–213.
- Cheng B. J. et al. High precision magnetometer for geomagnetic exploration onboard of the China Seismo-Electromagnetic Satellite. Science China Technological Sciences. 2018, Vol. 61, No. 5, P. 659–668.
- Crassidis J. L., Kok-Lam Lai, Harman R.R. Real-Time Attitude-Independent Three-Axis Magnetometer Calibration. Journal of Guidance, Control, and Dynamics. 2005, Vol. 28, No. 1, P. 115–120.
- Gebre-Egziabher D., Elkaim G. H., David Powell J., Parkinson B. W. Calibration of strapdown magnetometers in magnetic field domain. J. Aerospace Eng. 2006, Vol. 19, P. 87–102.
- Soken H. E. A Survey of Calibration Algorithms for Small Satellite Magnetometers. Measurement. October 2017. doi: 10.1016/j.measurement.2017.10.017.
- Soken H. E., Sakai S. Attitude estimation and magnetometer calibration using reconfigurable TRIAD+ filtering approach. Aerospace Science and Technology. 2020, Vol. 99, 105754.
- Springmann J. C., Cutler J. W. Attitude-independent magnetometer calibration with time-varying bias. Journal of Guidance, Control, and Dynamics. 2012, Vol. 35, No. 4, P. 1080–1088.
- Styp-Rekowski K., Michaelis I., Stolle C., Baerenzung J., Korte M., Kao O. Machine learning-based calibration of the GOCE satellite platform magnetometers. Earth Planets Space. 2022, Vol. 74, P. 138.
- Bakhvalov N. S. Chislennye metody [Numerical methods]. Moscow, Nauka Publ., 1975, 632 p.
- Vasil'ev F. P. Metody optimizatsii [Optimization methods]. Moscow, Faktorial Press Publ., 2002, 824 p.
- Kantorovich L. V., Akilov G. P. Funktsional'nyy analiz [Functional analysis]. Moscow, Nauka Publ., 1984, 752 p.
补充文件
