Аналитическое решение задачи калибровки магнитометракосмического аппарата с помощью метода наименьших квадратов
- Авторы: Кириллов К.А.1, Кириллова С.В.1, Кузнецов А.А.1, Мелентьев Д.О.2,3, Сафонов К.В.1
-
Учреждения:
- Сибирский государственный университет науки и технологий имени академика М. Ф. Решетнева
- Сибирский федеральный университет
- Акционерное общество «Информационные спутниковые системы» имени академика М. Ф. Решетнёва»
- Выпуск: Том 26, № 2 (2025)
- Страницы: 181-194
- Раздел: Раздел 1. Информатика, вычислительная техника и управление
- Статья опубликована: 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
Цитировать
Полный текст
Аннотация
В настоящей работе предложен аналитический метод решения задачи калибровки магнитометра для модели, рассмотренной в [1]. При решении задачи определения калибровочных параметров блока магнитометра учитывается, что для измерений с любой пространственной ориентацией блока магнитометра величина измеряемого вектора магнитной индукции сохраняется и является известной модельной величиной. Вводится в рассмотрение штрафная функция 12 переменных, равная сумме квадратов невязок. Алгоритм решения задачи калибровки измерительных осей блока магнитометра сводится к поиску методом наименьших квадратов таких значений переменных этой функции, которые при заданном наборе векторов измерений магнитометра доставляют ей минимум. С этой целью указанная функция исследуется на условный экстремум при наличии трех уравнений связи. Составляется функция Лагранжа и, исходя из необходимого условия экстремума этой функции, формируется система 15 уравнений относительно 15 неизвестных. Доказывается, что эта система имеет четыре решения. Выведены формулы, позволяющие получить компоненты каждого из этих четырех решений. В качестве решения задачи калибровки магнитометра рекомендуется выбрать решение указанной системы, доставляющее минимум функции Лагранжа.
Полный текст
Введение
Магнитометры входят в состав системы ориентации и стабилизации низкоорбитальных малогабаритных космических аппаратов. Они широко используются благодаря тому, что являются легкими, недорогими и надежными. Однако вследствие физических особенностей чувствительного элемента современные магнитометры требуют проведения математической калибровки. На данный момент предложены различные методы калибровки магнитометров, этим методам посвящено немалое количество научных работ [1–12], в частности, статья [9], в которой приведен обзор различных способов проведения таких операций.
В вышеупомянутых работах задача калибровки магнитометра космического аппарата решалась с применением численных методов. В настоящей статье предложен аналитический метод решения этой задачи для модели, рассмотренной в [1].
1. Модель погрешностей измерений вектора магнитной индукции
Обозначим через h значение измеренного вектора магнитной индукции при некотором пространственном положении блока магнитометра (БМ). Воспользуемся упрощенной моделью измерений, рассмотренной в [1]:
(1)
В (1) использованы следующие обозначения [1]:
B = (B1, B2, B3)T – истинный вектор магнитной индукции;
b = (b1, b2, b3)T – постоянный вектор, отвечающий смещениям нуля для каждой из измерительных осей БМ;
n – случайный вектор, отвечающий некоррелированному шуму для каждой из измерительных осей;
Р – матрица, строки которой есть орты измерительных осей БМ, записанные в «базовой» системе БМ;
Q – диагональная матрица, содержащая на главной диагонали масштабные коэффициенты k1, k2, k3 для измерительных осей БМ, т. е.
Другими словами, матрица P описывает неортогональность измерительных осей БМ, а матрица Q отвечает масштабированию по этим осям.
Задача калибровки измерительных осей БМ сводится к нахождению элементов матриц P и Q, а также вектора смещений нуля b.
2. Разработка алгоритма определения калибровочных параметров БМ
При решении задачи определения калибровочных параметров БМ воспользуемся тем, что для измерений с любой пространственной ориентацией БМ величина измеряемого вектора магнитной индукции B сохраняется и является известной модельной величиной.
Будем считать, что в полёте в результате измерений магнитометра в дискретные моменты времени получен набор векторов h(l) = (h1(l), h2(l), h3(l))T, l = 1, …, N. При условии отсутствия шумов измерений из формулы (1) получаем:
(2)
где B(l) = (B1(l), B2(l), B3(l))T – истинный вектор магнитной индукции в той же точке пространства, что и измеренный вектор h(l), l = 1, …, N. Выразим векторы (l = 1, …, N) из равенств (2):
(3)
l = 1, …, N, где
(4)
Как отмечено в [1], без ограничения общности матрицу неортогональности Р можно представить с минимальным количеством неизвестных элементов следующим образом:
где ε1, ε2, ε3 – малые углы. Тогда обратная матрица P–1 примет вид
где αi = tg εi, βi = sec εi, i = 1, 2, 3, а так как
то в соответствии с (4) получаем:
(5)
где γi = ki – 1, i = 1, 2, 3.
Перепишем равенство (3) с учетом (5):
(6)
l = 1, …, N. Каждое из N векторных равенств (6) запишем в виде системы трех скалярных равенств:
(7)
l = 1, …, N. В силу первого равенства системы (7) второе равенство этой системы перепишем в виде
откуда получаем равенство
вследствие которого третье равенство системы (7) записывается следующим образом:
l = 1, …, N.
Введем в рассмотрение функцию двенадцати переменных αi, βi, γi, bi (i = 1, 2, 3)
где переменные αi и βi в силу тригонометрического тождества
связаны равенствами
(8)
i = 1, 2, 3.
Алгоритм решения задачи калибровки измерительных осей БМ сводится к поиску методом наименьших квадратов [13] с учетом (8) таких значений переменных αi, βi, γi, bi (i = 1, 2, 3), которые при заданном наборе векторов измерений {h(l)} (l = 1, …, N) доставляют минимум функции Ф. С этой целью требуется исследовать функцию Ф на условный экстремум [14] при наличии трех уравнений связи (8).
Составим функцию Лагранжа
, (9)
пятнадцати переменных αi, βi, γi, bi, λi (i = 1, 2, 3) и запишем необходимое условие локального экстремума этой функции:
Требуется найти стационарные точки функции F, т. е. решение системы уравнений (10).
Введем следующие обозначения:
Учитывая эти обозначения, для удобства разобьем систему уравнений (10) на три нелинейные системы – систему уравнений относительно двух неизвестных b1, γ1:
(11)
систему уравнений относительно пяти неизвестных b2, α1, β1, γ2, λ1:
(12)
и систему уравнений относительно восьми неизвестных b3, α2, α3, β2, β3, γ3, λ2, λ3:
(13)
В записи выражений βi через αi в последних уравнениях систем (12) и (13) мы воспользовались тем, что βi = sec εi > 0 вследствие очевидных неравенств –π/2 < εi < π/2, i = 1, 2, 3. Решим каждую из систем уравнений (11), (12), (13).
Сначала заметим, что справедливо неравенство
(14)
i = 1, 2, 3. Действительно, дискриминант квадратного уравнения – Fi + 2Hibi – Nbi² = 0 относительно неизвестной bi, будучи равным
отрицателен в силу очевидного следствия
неравенства Коши – Буняковского [15] (здесь поставлен знак « < » вместо « ≤ », так как хотя бы одно из слагаемых hi(1), …, hi(N) заведомо не равно нулю), поэтому рассматриваемое квадратное уравнение не имеет вещественных корней, i = 1, 2, 3.
Решим систему уравнений (11). Учитывая неравенство (14), выражаем γ1 через b1 из первого уравнения этой системы:
Подставляя вместо γ1 полученное выражение во второе уравнение системы (11) и принимая во внимание неравенство γ1 = k1– 1 ≠ 0, после простейших преобразований приходим к равенству
с учетом которого записанное выше равенство для γ1 приводится к виду
Таким образом, решение системы уравнений (11) найдено.
Рассмотрим систему уравнений (12). Поскольку γ2 = k2 – 1 ≠ 0, второе уравнение этой системы можно записать в виде
(15)
а так как β1 = sec ε1 ≠ 0, третье уравнение системы (12) равносильно следующему уравнению:
(16)
Из (15) с учетом (16) получаем, что λ1β1γ2– 1 = 0 и тогда в силу неравенства β1 ≠ 0 имеем
(17)
Имеет место неравенство D12 – C1b2 ≠ 0, так как непосредственная проверка показывает, что значение b2 = D12C1– 1, являющееся корнем уравнения D12 – C1b2 = 0, не удовлетворяет уравнениям системы (12).
Сначала рассмотрим случай H2 – Nb2 ≠ 0, т. е. b2 ≠ H2N– 1. Выразим β1γ2 из (16), а также из первого (с учетом равенства (17)) и четвертого (с учетом неравенства β1γ2 ≠ 0) уравнений системы (12):
(18)
По доказанному знаменатель дроби во втором из равенств (18) не равен нулю, знаменатель дроби в первом из равенств (18) отличен от нуля в силу (14). Из первого и третьего равенств (18) следует уравнение
из которого выражаем b2 через α1:
(19)
Аналогично, из второго и третьего равенств (18) следует уравнение
из которого также выражаем b2 через α1:
(20)
Приравнивая правые части равенств (19) и (20), после простейших преобразований получаем квадратное уравнение относительно α1:
(21)
Дискриминант уравнения (21) равен
а корни этого уравнения вычисляются по формулам
(22)
Заметим, что значения (22) неизвестной α1 не являются корнями знаменателей дробей в равенствах (19) и (20). Подставляя первое из значений (22) неизвестной α1 в четвертое уравнение системы (12), с учетом неравенства β1γ2 ≠ 0 приходим к равенству H2 – Nb2 = 0, которое противоречит рассматриваемому случаю.
Пусть теперь H2 – Nb2 = 0, т. е. b2 = H2N– 1. Принимая во внимание неравенство β1γ2 ≠ 0, из четвертого уравнения системы (12) получаем, что α1 = – C2C1– 1, т. е. найденное значение α1 совпадает с первым из корней (22) уравнения (21). Несложно показать, что равенства (19) и (20), доказанные в предположении b2 ≠ H2N– 1, справедливы и при b2 = H2N– 1. В случае b2 = H2N– 1 имеют место также и первые два равенства из (18), так как при их выводе никакие ограничения на значение неизвестной b2 не учитывались.
Из проведенных рассуждений вытекает, что система уравнений (12) имеет два решения. Компоненты обоих решений указанной системы получаем в следующем порядке:
- – значения неизвестной α1 находим по формулам (22);
- – значения β1, соответствующие найденным значениям α1, вычисляем по формуле, записанной в качестве последнего уравнения системы (12);
- – значения b2 – по любой из формул (19), (20);
- – значения γ2 – по любой из двух формул
которые вытекают соответственно из первого и второго равенств (18);
- – значение неизвестной λ1 в обоих решениях системы уравнений (12) в силу (17) принимаем равным нулю.
Перейдем к отысканию решений системы уравнений (13). В силу неравенств γ3 = k3– 1 ≠ 0, β2 = sec ε2 ≠ 0 и β3 = sec ε3 ≠ 0 второе и третье уравнения этой системы можно записать как
(23)
, (24)
соответственно. Из (23) с учетом (24) получаем: λ2β2(β3γ3)–1 = 0. Следовательно, справедливо равенство
(25)
Тогда первое уравнение системы (13) с учетом неравенства β3 ≠ 0 запишется в виде
(26)
Из шестого уравнения системы (13) с учетом (24) и (26) находим – λ3β3 = 0, откуда следует равенство
(27)
Имеют место неравенства H3 – Nb3 ≠ 0, Di3 – Cib3 ≠ 0, так как непосредственная проверка показывает, что ни корень b3 = H3N– 1 уравнения H3 – Nb3 = 0, ни корни b3 = Di3Ci– 1 уравнений Di3 – Cib3 = 0 не удовлетворяют уравнениям системы (13), i = 1, 2.
С учетом равенств (25), (27) и неравенств β2 ≠ 0, β3 ≠ 0, γ3 ≠ 0 выразим β2β3γ3 из первого, третьего, четвертого и пятого уравнений системы (13):
(28)
По доказанному знаменатели дробей в первом, третьем и четвертом из равенств (28) не равны нулю, знаменатель дроби во втором из равенств (28) отличен от нуля в силу (14). Из первого и третьего равенств (28) следует уравнение
из второго и третьего равенств (28) – уравнение
из третьего и четвертого равенств (28) – уравнение
Из последних трех уравнений выражаем α2β3 через α3 и b3:
(29)
(30)
(31)
Приравняв правую часть равенства (29) к правым частям равенств (30) и (31), после простейших преобразований получаем уравнения
из которых выражаем α3 через b3, предварительно разделив обе части каждого из них на (H3 – Nb3):
(32)
(33)
Приравняв правые части равенств (32) и (33), после простейших преобразований приходим к квадратному уравнению относительно b3
(34)
коэффициенты которого определяются следующим образом:
Дискриминант уравнения (34) равен
а корни этого уравнения вычисляются по формулам
(35)
Так как уравнение (34) имеет два решения, то два решения имеет и система уравнений (13). Компоненты обоих решений этой системы получаем в следующем порядке:
- – значения неизвестной b3 находим по формулам (35);
- – значения неизвестной α3, соответствующие найденным значениям b3, вычисляем по любой из формул (33), (34);
- – значения β3 – по формуле, записанной в качестве последнего уравнения системы (13);
- – значения α2 – по любой из трех формул
которые вытекают из равенств (29), (30) и (31) соответственно;
- – значения β2 – по формуле, записанной в качестве последнего уравнения системы (13);
- – значения γ3 – по любой из четырех формул
которые вытекают из (28);
- – значения неизвестных λ2 и λ3 в обоих решениях системы уравнений (13) в силу равенств (25) и (27) принимаем равными нулю.
Итак, система уравнений (11) имеет единственное решение, а каждая из систем (12), (13) – два решения. Следовательно, число решений исходной системы уравнений (10), определяемое как произведение числа решений систем (11), (12) и (13), равно четырем. Из полученных четырех решений
, (36)
(j = 1, 2, 3, 4) системы (10), которые также являются стационарными точками функции Лагранжа F, определенной равенством (9), нас интересует решение
удовлетворяющее равенству
Так как в каждом из решений (36) системы уравнений (10) последние три компоненты равны нулю, т. е.
j = 1, 2, 3, 4, то в силу (9) имеем:
j = 1, 2, 3, 4, поэтому искомое решение определяем по формуле
Искомые значения масштабных коэффициентов для измерительных осей БМ вычисляем по формулам
i = 1, 2, 3, искомые значения углов – по формулам
которые следуют из равенств
и очевидных неравенств i = 1, 2, 3.
ЗАКЛЮЧЕНИЕ
Итак, мы получили аналитическое решение задачи калибровки магнитометра космического аппарата для модели, рассмотренной в [1]. Процедура вычисления калибровочных параметров БМ по выведенным формулам обладает рядом очевидных преимуществ по сравнению с численными методами решения этой задачи:
- – существенно уменьшается число арифметических операций;
- – исчезает проблема возможной неустойчивости метода.
Об авторах
Кирилл Анатольевич Кириллов
Сибирский государственный университет науки и технологий имени академика М. Ф. Решетнева
Автор, ответственный за переписку.
Email: kkirillow@yandex.ru
ORCID iD: 0000-0002-3763-1303
доктор физико-математических наук, доцент, профессор кафедры прикладной математики
Россия, 660037, г. Красноярск, просп. им. газ. «Красноярский Рабочий», 31Светлана Владимировна Кириллова
Сибирский государственный университет науки и технологий имени академика М. Ф. Решетнева
Email: svkirillova2009@yandex.ru
ORCID iD: 0000-0003-3779-2825
кандидат технических наук, доцент, доцент кафедры прикладной математики и анализа данных
Россия, 660037, г. Красноярск, просп. им. газ. «Красноярский Рабочий», 31Александр Алексеевич Кузнецов
Сибирский государственный университет науки и технологий имени академика М. Ф. Решетнева
Email: alex_kuznetsov80@mail.ru
ORCID iD: 0000-0003-0944-1817
доктор физико-математических наук, профессор, директор НОЦ «Институт космических исследований и высоких технологий»
Россия, 660037, г. Красноярск, просп. им. газ. «Красноярский Рабочий», 31Денис Олегович Мелентьев
Сибирский федеральный университет; Акционерное общество «Информационные спутниковые системы» имени академика М. Ф. Решетнёва»
Email: denes.2000@mail.ru
ORCID iD: 0009-0009-6187-4098
аспирант, инженер
Россия, 660041, г. Красноярск, просп. Свободный, 79; 662972, г. Железногорск Красноярского края, ул. Ленина, 52Константин Владимирович Сафонов
Сибирский государственный университет науки и технологий имени академика М. Ф. Решетнева
Email: safonovkv@rambler.ru
доктор физико-математических наук, профессор, директор Института информатики и телекоммуникаций
Россия, 660037, г. Красноярск, просп. им. газ. «Красноярский Рабочий», 31Список литературы
- Калибровка магнитометра космического аппарата «Декарт» в полете / П. Е. Розин, А. В. Симонов, Е. С. Гордиенко, Ю. К. Зайко // Тр. МАИ. 2022. № 124. С. 1–20.
- Методика калибровки магнитометра на этапе наземной диагностики систем космического аппарата / И. О. Акимов, С. Н. Илюхин, Н. А. Ивлев, Г. Е. Колосов // Инженерный журнал: наука и инновации. 2018. № 5 (77). С. 1–18.
- Баллистическое обеспечение разработки и полётов межорбитального космического буксира «Фрегат» / И. М. Морской, А. В. Симонов, В. И. Лясковская, А. С. Ежов // Вестник «НПО им. С. А. Лавочкина». 2014. № 1. С. 10–15.
- 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.
- Verification and calibration of a commercial anisotropic magnetoresistive magnetometer by multivariate non-linear regression / N. Belsten, M. Knapp, R. Masterson et al. // Geosci. Instrum. Method. Data Syst. 2023. No. 12. P. 201–213.
- High precision magnetometer for geomagnetic exploration onboard of the China Seismo-Electromagnetic Satellite / B. J. Cheng et al. // 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.
- Calibration of strapdown magnetometers in magnetic field domain / D. Gebre-Egziabher, G. H. Elkaim, J. David Powell, B. W. Parkinson // 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.
- Machine learning-based calibration of the GOCE satellite platform magnetometers // K. Styp-Rekowski, I. Michaelis, C. Stolle et al. // Earth Planets Space. 2022. Vol. 74. P. 138.
- Бахвалов Н. С. Численные методы. М. : Наука, 1975. 632 с.
- Васильев Ф. П. Методы оптимизации. М. : Факториал Пресс, 2002. 824 с.
- Канторович Л. В., Акилов Г. П. Функциональный анализ. М. : Наука, 1984. 752 с.
Дополнительные файлы
