CC BY
139
УДК 62-83: 531.3 А.А. Карауш
Выбор численного метода интегрирования дифференциальных уравнений для задач спутниковых навигационных технологий
Проведен сравнительный анализ численных методов интегрирования уравнений движения навигационных спутников ГЛОНАСС. По показателям точности сравнивались численные методы интегрирования дифференциальных уравнений Рунге-Кутты, Адамса-Штёрмера, Адамса-Мультона-Коуэлла и Эверхарта. Необходимость такого сравнения следует из особенно высоких требований к точности, предъявляемых к численным схемам интегрирования уравнений движения НС в задачах сегмента эфемеридно-временного обеспечения спутниковой навигационной системы ГЛОНАСС. Целью исследований, проводимых в настоящей работе, является обоснованный выбор численных схем интегрирования для задач восстановления орбит навигационных спутников по данным беззапросных траекторных измерений, выполняемых наземными пунктами. Особое внимание в сравнительном анализе уделено тому, как рассматриваемые схемы ведут себя в условиях действия на правую часть уравнений движения скачкообразных возмущений. Они соответствуют воздействию на спутник радиационного давления солнечного излучения, которое претерпевает скачок в момент входа НС в тень Земли и выхода из неё. Приведены обоснования для использования метода Эверхарта в качестве наиболее эффективного метода численного интегрирования уравнений движения навигационных спутников ГЛОНАСС в рамках задачи эфемеридно-временного обеспечения системы.
Ключевые слова: навигационный спутник, расчет орбит, ГЛОНАСС, эфемеридно-временное обеспечение, численное интегрирование, метод Эверхарта.
При решении широкого круга задач спутниковых навигационных технологий возникает необходимость численного интегрирования уравнений движения навигационных спутников (НС) Земли [1, 2].
Задачи основного назначения спутниковых навигационных систем - определение координат и параметров движения объектов потребителя с помощью аппаратуры приёма навигационных сигналов. В этой аппаратуре, согласно [3], для численного интегрирования дифференциальных уравнений движения НС ГЛОНАСС применяются численные схемы Рунге-Кутты. В качестве начальных условий используются эфемериды, транслируемые с бортов НС в составе навигационного сигнала. Поскольку эфемеридная информация от каждого НС обновляется через 30 мин, погрешности интегрирования не успевают накапливаться и применение достаточно простых схем Рунге-Кутты в пользовательском сегменте считается вполне обоснованным.
Особенно высокие требования предъявляются к численным схемам интегрирования уравнений движения НС для задач, решаемых в сегменте эфемеридно-временного обеспечения спутниковой навигационной системы ГЛОНАСС. Здесь по результатам выполняемых с Земли траекторных измерений с помощью алгоритмов восстановления орбит строятся математические модели движения НС. Полученные таким образом математические модели используются для расчета орбит НС на интервалы времени длительностью в несколько суток. Для таких рассчётов требуются численные схемы, не накапливающие погрешностей интегрирования в специфичных условиях движения НС
ГЛОНАСС.
Большую важность в сегменте эфемеридно-временного обеспечения ГЛОНАСС играют задачи высокоточной синхронизации пространственно разнесенных часов по сигналам спутниковых навигационных систем [4]. В алгоритмах синхронизации используются вычисленные геометрические дальности от НС до антенны потребителя. Для расчета указанных геометрических дальностей также требуется высокоточное восстановление орбит НС путём интегрирования численными методами уравнений движения.
Приводимые в настоящей работе исследования направлены на обоснованный выбор численных схем интегрирования именно для задач восстановления орбит навигационных спутников. При этом основными требованиями к численным схемам являются:
- высокая точность,
- вычислительная эффективность,
- малый уровень погрешностей интегрирования, накапливающихся с течением времени.
Особое внимание в сравнительном анализе уделено тому, как рассматриваемые схемы ведут себя в условиях действия на правую часть уравнений движения скачкообразных возмущений. Они соответствуют воздействию на спутник радиационного давления солнечного излучения, которое претерпевает скачок в момент входа НС в тень Земли и выхода из неё.
Постановка задачи. Целью настоящей работы является сравнительный анализ численных методов, применяемых для интегрирования уравнений движения НС по околокруговым орбитам, и обоснованный выбор на этой основе наиболее эффективных численных схем.
Для указанного сравнительного анализа численных методов выбраны математические модели движения в классе дифференциальных уравнений [2, 6, 7]:
Ц 5
и (?) = —3--и(0 + Х 8/ (?), и(?о) = ио,й (?о) = и 0, (1)
р3(?) /=1
где иТ (?) = [х(?), у(?),¿(?)] - вектор текущих координат НС в квазиинерциальной системе координат і 2 2 2
(ИСК); р(?) = у]х (?) + у (?) + г (?) - текущий радиус орбиты НС; ц - геоцентрическая постоянная гравитационного поля Земли; 8/(?) - вектор действующих на НС возмущений; / = 1 - от несферич-ности гравитационного поля Земли; / = 2,3 - от гравитационного воздействия на НС Луны и Солнца; / = 4 - от радиационного давления на НС солнечного излучения.
Возмущения 8/(?), / =1,...,4 образуют группу моделируемых возмущений, которые с той или иной точностью могут быть представлены соответствующими математическими моделями и учтены при расчете орбитального движения НС.
/ = 5 - немоделируемые возмущения случайной природы, для которых в лучшем случае могут быть получены статистические характеристики.
Возмущения s/ (?) приводятся к центру масс НС в объектоцентрической системе координат и преобразуются в ИСК с помощью известных матричных преобразований [5, 7].
Наряду с общими требованиями точности к численным схемам интегрирования предъявляются дополнительные требования, порожденные особенностями уравнения движения НС вида (1):
1. Решение уравнения (1), не содержащего демпфирующих членов, находится на границе устойчивости. Выбранная численная схема не должна нарушать это условие и не накапливать погрешностей при интегрировании на больших интервалах времени.
2. Численная схема не должна терять точность интегрирования при скачкообразных возмущениях в правой части уравнения (1). Такое изменение возмущения 84 (?) возможно при вхождении НС в тень Земли, когда скачком снимается радиационное давление на НС солнечного излучения.
Анализ точности численных схем принято строить на основе сравнения аналитического решения уравнения движения и решений этого уравнения, полученных с помощью сравниваемых численных схем. Однако для общего случая, при условии действия на НС всех возмущений, уравнение
(1) в квадратурах не интегрируется [5, 6]. По этой причине для сравнительного анализа используется упрощенное уравнение
■¿(І) = -ю2 -г(?) + ^(?) + S4(t), ¿(?о) = го,¿(?о) = ¿0 , ®2 =“3— , (2)
р3(?)
решением которого является проекция на ось Т ИСК движения НС в плоскости (ХОТ) , для которого ^1 (?) имитирует возмущение от второй зональной гармоники амплитудой а, S4(t) - возмущение от радиационного давления на НС солнечного излучения в виде прямоугольного импульса 0(? - Т[) амплитудой Ь с моментом появления Т| и моментом снятия Т2 . Аналитическое решение уравнения
(2) представляется выражением [8]
2 2 a • sin (ю?) a• sin(ra?) • sin(3ra?) a• cos (ю?) a• cos(3ra?)• cos(ra?)
z (t)=---------2—+----------2------------------2—+-------------------2-+
2ю2 6ю2 2ю2 6ю2
b •Q(t— T¡) b •Q(t—T|) •sin(T|ra) •sin(ra?) b •Q(t— T¡) •cos(T]ra) •cos(raí)
+ 2 2 2
ю ю ю
b-9(?—T2) •sin2(юt) b ^9(t—T2) •sinT w) •sin(юt) b-9(?—T2) •cos2(юt)
2 + 2 2 +
ю ю ю
b ^9(t—T2) •cos(T2ю) •cos(юt) .
+---- ---—-------------- —- + c2 • sin(ю t) + q • cos(ю t).
ю2
Параметры q и C2 определяются из начальных условий z(?о) = ¿о,¿(?о) = ¿o.
Результаты. Сравнение численного интегрирования уравнения (2) с его аналитическим решением на интервале времени 105 с проводилось для схем:
1) Рунге-Кутты 4-го порядка [5, 6];
2) Адамса-Штёрмера 5-го порядка [5];
3) Эверхарта 15-го порядка [7];
4) Адамса-Мультона-Коуэлла 16-го порядка [6, 8].
Численные схемы 2, 3, 4 не требуют сведения уравнения (2) к системе уравнений 1-го порядка. Схемы 1, 3 одношаговые, для схем 2, 4 разгонная сетка решений получена на основе аналитического решения уравнения (2).
7 9 3 2
При моделировании выбирались: радиус орбиты р = 2,5-10 м, ц = 398600-10 м /с ,
—6 —7
а = ^,•10 , b = ^,•10 . Начальные условия ¿(?о) = ¿о, ¿(?о) = ¿о выбирались так, чтобы движение
НС начиналось от экватора.
В качестве конструктивного параметра численной схемы рассматривался шаг интегрирования (h = 100 с, h = 1000 с).
Результаты сравнительного анализа численных схем для различных значений шага интегрирования приведены в табл. 1-3. Для каждого из сравниваемых методов оценивалась погрешность интегрирования в метрах в виде:
- абсолютного значения погрешности на конец интервала интегрирования 8 -конечное, м;
- среднего арифметического на интервале интегрирования 8 -среднее арифм., м;
- максимального значения на интервале интегрирования 8 -max, м;
- времени выполнения интегрирования, с;
- количества вызовов функции вычисления правой части уравнения F_Eval.
Таблица 1
Результаты сравнения численных методов для h = 100 с, b = 0__________
Метод 8 -конечное, м 8 -среднее арифметическое, м 8 -max, м Время, с F_Eval
Рунге-Кутты 1,2000E-02 6,7600E-02 7,5570E-01 0,049 4000
Адамса-Штёрмера 2,6659E-06 -4,2711E-08 2,6934E-06 0,053 6006
Эверхарта 2,7940E-09 1,1191E-09 3,1665E-08 0,506 22051
Адамса-Мультона-Коуэлла 4,4694E-06 8,8574E-07 8,3372E-06 0,204 1991
Таблица 2
Результаты сравнения численных методов для h = 1000 с, b = 0_________
Метод 8 -конечное, м 8 -среднее арифметическое, м 8 -max, м Время, с F_Eval
Рунге-Кутты 1,8273E+03 1,3073E+01 1,8273E+03 0,005 400
Адамса-Штёрмера 1,8959E+00 -1,4830E-01 2,5660E+00 0,005 606
Эверхарта 2,6077E-08 -5,1868E-10 4,8429E-08 0,050 2251
Адамса-Мультона-Коуэлла 6,2514E+02 7,1383E+01 7,1766E+02 0,043 191
Таблица 3
Результаты сравнения численных методов для Н = 100 с со ступенчатым воздействием Ь ф 0
Метод 8 -конечное, м 8 -среднее арифметическое, м 8 -тах, м Время, с Р_Еуаі
Рунге-Кутты 3,1900Е-02 7,0400Е-02 8,0340Е-01 0,054 4000
Адамса-Штёрмера 1,0100Е-02 1,1000Е-03 2,0500Е-02 0,052 6006
Эверхарта 1,3970Е-08 -7,2660Е-10 5,2154Е-08 0,492 22051
Адамса-Мультона-Коуэлла 8,4000Е-03 1,2000Е-03 2,0500Е-02 0,276 1991
Анализ характера изменения абсолютных погрешностей интегрирования исследуемых схем показывает, что для схемы Рунге-Кутты погрешности нарастают на интервале интегрирования до 3 м, погрешности схемы Адамса-Штёрмера имеют приемлемый уровень, но существенно увеличиваются в случае скачкообразного возмущения правой части уравнения движения НС. Схема Эверхарта обеспечивает минимальный уровень погрешностей интегрирования.
В частности, на рис. 1 показан характер изменения абсолютных значений погрешностей интегрирования уравнений движения НС схемой Адамса-Штёрмера в режиме ступенчатого возмущения в правой части уравнения движения НС. Правая часть уравнения претерпевает скачок на 110-м шаге интегрирования, что приводит к видимому на графике гармоническому изменению абсолютной погрешности.
Погрешность, м Схема Адамса-Штёрмера
0,05
0,04
0,03
Рис. 1. Характер изменения абсолютных погрешностей интегрирования
Заключение
1. Схема Рунге-Кутты показала свою работоспособность при малом шаге интегрирования 100 с и неприемлемую точность при шаге 1000 с. Как правило, одношаговые схемы применяются в сочетании с многошаговыми схемами для заполнения с малым шагом разгонной сетки решений с целью запуска многошаговых схем.
2. Многошаговые методы Адамса-Штёрмера и Адамса-Мультона-Коуэлла обеспечивают приемлемую точность интегрирования при шаге 100 с и шаге 1000 с и требуют знания решения в узлах разгонной сетки. Количество узлов зависит от порядка применяемой схемы.
Это создает определенные неудобства, связанные с необходимостью сгущения сетки решений из-за риска потери точности интегрирования на участках вхождения НС в теневые зоны Земли, где возникают скачкообразные изменения радиационного давления солнечного излучения на НС.
3. Наилучшие по точности результаты показала одношаговая схема Эверхарта в условиях интегрирования с шагом 100 с и с шагом 1000 с гладкой функции и при скачкообразных возмущениях в правой части уравнения (2).
Дополнительный положительный эффект применения схемы Эверхарта заключается в том, что погрешность интегрирования оказывается знакопеременной функцией времени и обеспечивает тенденцию не накапливаться на интервале интегрирования. Об этом свидетельствует малость среднего арифметического погрешности по сравнению с выборочными значениями этой погрешности.
Определенным неудобством применения схемы Эверхарта, по сравнению с рассмотренными схемами, является большое число обращений к вычислению правой части уравнения (2), характеризующееся числом Р_Буа1 в табл. 1-3.
Литература
1. Сурнин Ю.В. Модели движения ИСЗ и точность численного прогнозирования / Ю.В. Сурнин, С.В. Кужелев // Геодезия и картография. - 1982. - № 10. - С. 8-13.
2. Бордовицына Т.В. Современные численные методы в задачах небесной механики. - М.: Наука, 1984. - 136 с.
3. Интерфейсный контрольный документ ГЛОНАСС. - М.: Российский научно-исследовательский институт космического приборостроения, 2008. - 74 с.
4. Толстиков А.С. Алгоритмы синхронизации пространственно разнесённых часов по сигналам спутниковых навигационных систем // Метрология: приложение к журналу «Измерительная техника». - 2009. - № 9. - С. 25-35.
5. Справочное руководство по небесной механике и астродинамике / В.К. Абалакин, Е.П. Аксенов, Е.А. Гребенников и др. - М.: Наука, 1976. - 864 с.
6. Чеботарев Г.А. Аналитические и численные методы небесной механики. - М.: Наука, 1965. -368 с.
7. Бордовицына Т.В. Теория движения искусственных спутников Земли / Т.В. Бордовицына, В. А. Авдюшев. - Томск: Изд-во Том. ун-та, 2007. - 178 с.
8. Камке Э. Справочник по обыкновенным дифференциальным уравнениям. - М.: Гос. изд-во физико-математической литературы, 1961. - 704 с.
Карауш Артем Андреевич
Аспирант Новосибирского государственного технического университета, науч. сотрудник ФГУП «СНИИМ»
Тел.: +7-923-603-30-39 Эл. почта: [email protected]
Karaush A.A.
Choice of numerical integration method of differential equations for global navigation satellite systems
Comparative analysis of numerical integration methods for GLONASS satellite movement equations is carried out in this paper. Four methods of differential equations numerical integration were compared in accuracy: Runge-Kutta, Adams-Stermer, Adams-Moulton-Cowell and Everhart methods. Problems of GLONASS ephem-eris and time supply require especially high accuracy in numerical integration of satellite orbits. It makes such comparison very necessary. Aim of this research is to make a well-reasoned choice of numerical integration method for recovering satellite orbits using trajectory observations carried out by terrestrial base stations. Methods behavior in case of step disturbance in the right part of satellite movement equation is specially considered. Step disturbance takes place in cases, when satellite comes into a shade of Earth and out of it. It is shown, that Everhart method is most effective for aims of numerical integration of satellites movement equations in a GLONASS ephemeris and time supply problem.
Keywords: navigation satellite, orbits calculation, GLONASS, ephemeris and time supply, numerical integration, Everhart method, mathematical movement model, trajectory observations.
