СПОСОБ ОПРЕДЕЛЕНИЯ ПАРАМЕТРОВ ОРБИТЫ ИСКУССТВЕННОГО СПУТНИКА ЗЕМЛИ Российский патент 2018 года по МПК G01S5/00 

Изобретение относится к области космонавтики, а именно к технике выполнения траекторных измерений и определения параметров орбиты искусственного спутника Земли (ИСЗ), и может быть использовано на наземных и бортовых комплексах управления полетом ИСЗ для точного определения текущих параметров движения ИСЗ.

Известен способ определения параметров орбиты ИСЗ с использованием лазерных светодальномеров (ЛСД). Данный способ широко известен и описан, например, в статье [1]. Для реализации этого способа выполняют следующие этапы:

размещают несколько ЛСД на земной поверхности с известными координатами;

оборудуют ИСЗ уголковым световозвращателем;

излучают импульс лазера с помощью ЛСД в направлении на ИСЗ, запоминают его копию;

запоминают время излучения импульса лазера;

принимают с помощью ЛСД, отраженный от уголкового световозвращателя ИСЗ, импульс лазера;

измеряют и запоминают время приема импульса лазера;

рассчитывают по времени излучения и времени приема импульса лазера наклонную дальность от ЛСД до ИСЗ;

производят серию измерений наклонных дальностей от ИСЗ до одной либо нескольких ЛСД, используя которые определяют параметры орбиты ИСЗ.

Недостатками способа определения координат ИСЗ с использованием ЛСД [1] являются:

высокая суммарная стоимость устройств, реализующих способ;

низкая производительность способа, обусловленная необходимым условием работы только ночью, причем небо должно быть ясным над всеми ЛСД.

Известен способ определения координат навигационных ИСЗ [2], заключающийся в том, что:

выбирают группу из четырех навигационных ИСЗ, находящихся в зоне прямой видимости, состоящей из первой пары навигационных ИСЗ, находящихся на одной орбите, и второй пары навигационных ИСЗ, находящихся на другой орбите;

одновременно измеряют линейные межспутниковые расстояния между всеми четырьмя навигационными ИСЗ группы;

передают измеренные значения линейных расстояний от каждого ИСЗ группы к каждому;

вычисляют на каждом навигационном ИСЗ сферические расстояния между ним и координаты точки пересечения орбит выбранных навигационных ИСЗ;

определяют значения координат навигационных ИСЗ на основе значений сферических расстояний между ними и координат точки пересечения орбит.

К недостаткам способа определения координат навигационных ИСЗ [2] относят:

длительное время определения координат навигационных ИСЗ, связанное с необходимостью циркуляции информации между соседними ИСЗ;

высокую стоимость устройства, реализующего способ, обусловленную необходимостью размещения на борту навигационных ИСЗ дополнительных передатчиков, приемников, а также аппаратуры обработки информации.

Из известных способов наиболее близким аналогом (прототипом) предлагаемого способа по технической сущности является способ определения параметров орбиты геостационарного ИСЗ [3]. Для реализации данного способа последовательно выполняют следующие этапы:

размещают приемо-излучающую наземную радиотехническую станцию (НРТС) на позиции с известными координатами с известными координатами;

выбирают начальные (априорные) значения параметров орбиты геостационарного ИСЗ;

измеряют наклонную дальность от НРТС до геостационарного ИСЗ;

многократно, в течении мерного интервала, составляющего не менее двух суток, измеряют наклонную дальность от НРТС до геостационарного ИСЗ с интервалом 2…3 часа;

делят полученную совокупность измерений на две части: обучаемую и контрольную;

делят определяемые параметры орбиты ИСЗ на две группы: внутриплоскостные и внеплоскостные;

производят серию оценок внутриплоскостных параметров орбиты ИСЗ при различных фиксированных значениях внеплоскостных параметров;

выбирают параметры орбиты геостационарного ИСЗ на основе серии произведенных оценок.

Недостатком способа прототипа является большая продолжительность мерного интервала проведения измерений дальностей, которая должна составлять не менее двух последовательных суток.

Целью изобретения является разработка способа с меньшим временем, затрачиваемым на определение параметров орбиты ИСЗ на основе использования тестовых сигналов излучающих опорных реперных станции (ИОРС) на позициях с известными координатами.

Поставленная цель достигается тем, что в известном способе определения параметров орбиты ИСЗ, включающем размещение НРТС на позиции с известными координатами x_K, y_K, z_K, выбор начальных значений параметров орбиты ИСЗ, измерение наклонной дальности от ИСЗ до НРТС, выбор параметров орбиты ИСЗ на основе измеренной наклонной дальности, дополнительно устанавливают на земной поверхности N≥2 ИОРС на позициях с известными координатами x_In, y_In, z_In, где n=1…N - номер ИОРС. Перед измерением наклонной дальности R_SK от ИСЗ до НРТС предварительно, излучают в момент времени t₀ тестовые радиосигналы НРСТ на средней частоте ƒ₀ и N ИОРС на средних частотах ƒ_n соответственно, запоминают время излучения тестовых радиосигналов t₀, и номиналы частот тестовых радиосигналов ƒ₀ и ƒ_n, принимают в НРТС ретранслированные ИСЗ тестовые радиосигналы излученные НРТС и N ИОРС в моменты времени и , на средних частотах и соответственно. Далее передают на НРСТ номиналы частот тестовых сигналов ƒ_n от N ИОРС и определяют задержки во времени Δt₀, Δt_n, между излученными и принятыми тестовыми радиосигналами. Однократно измеряют наклонную дальность от ИСЗ до НРСТ R_SK, и дополнительно измеряют наклонные дальности от ИСЗ до N ИОРС . Вычисляют координаты ИСЗ x₀, y₀, z₀ по известным координатам НРСТ, известным координатам N ИОРС и измеренным наклонным дальностям R_SK, и . Вычисляют ортогональные составляющие вектора скорости ИСЗ по известным координатам НРСТ, известным координатам N ИОРС, вычисленным координатам ИСЗ x₀, y₀, z₀, запомненным номиналам средних частот переданных ƒ₀, ƒ_n и принятых тестовых радиосигналов, предварительно заданной частоте сдвига ƒ_G рабочей частоты ИСЗ. В качестве параметров орбиты ИСЗ принимают совокупность координат ИСЗ x₀, y₀, z₀ и ортогональных составляющих вектора его скорости в момент времени t₀.

Благодаря перечисленной новой совокупности существенных признаков, за счет использования тестовых радиосигналов N ИОРС на позициях с известными координатами достигается цель изобретения: снижение времени, затрачиваемого на определение параметров орбиты ИСЗ.

Заявленный способ поясняется чертежами, на которых показаны:

на фиг. 1 - положение орбиты ИСЗ,

на фиг. 2 - аномалии орбиты ИСЗ,

на фиг. 3 - схема алгоритма расчета аномалий ИСЗ,

на фиг. 4 - система трех возмущающих тел: Земли, Солнца и Луны,

на фиг. 5 - схема алгоритма прогноза координат ИСЗ,

на фиг. 6 - сферический треугольник Леонардо-Эйлера,

на фиг. 7 - структурная схема подсистемы определения координат ИСЗ,

на фиг. 8 - сферические поверхности положения ИСЗ,

на фиг. 9 - схема алгоритма расчета координат ИСЗ,

на фиг. 10 - структурная схема подсистемы определения ортогональных составляющих вектора скорости ИСЗ.

Теория полета ИСЗ, или, как ее еще называют, астродинамика, небесная механика, космическая баллистика, основана на законах И. Кеплера и законе всемирного тяготения И. Ньютона.

В первом приближении движение ИСЗ представляется как невозмущенное - такое движение, которое происходило бы только под влиянием силы притяжения Земли по закону Ньютона, т.е. точно соответствует задаче двух тел (Земля - ИСЗ) в небесной механике. Это движение называется движением по Кеплеровой орбите, так как подчиняется трем законам Кеплера [4]. Элементы Кеплеровой орбиты ИСЗ раскрыты в приложении А.

Достоинством Кеплеровой орбиты является простота вычисления координат и вектора скорости ИСЗ в прогнозируемый момент времени. Это предопределило широкое использование элементов Кеплеровой орбиты. В настоящем изобретении элементы Кеплеровой орбиты ИСЗ выступают в качестве априорных данных об ИСЗ. С помощью этих элементов ориентируют передающие антенны НРСТ и N ИОРС на ИСЗ, а также приемную антенну НРСТ на ИСЗ для организации передачи и приема тестовых радиосигналов. Кроме того, элементы Кеплеровой орбиты ИСЗ служат для устранения двузначности определения координат ИСЗ на основе тестовых сигналов НРСТ и двух ИОРС.

Недостатком Кеплеровой орбиты является относительно низкая точность определения координат и вектора скорости ИСЗ, которая недостаточна для выполнения ряда прикладных задач, например, для корректировки движения ИСЗ.

Более точно движение ИСЗ описывается с помощью возмущенной орбиты [5], которая задается:

каноническими параметрами ИСЗ, включающими координаты ИСЗ x₀, y₀, z₀ и ортогональные составляющие вектора его скорости в начальный момент времени t₀,

факторами, приводящими к отклонениям ИСЗ от идеальной (Кеплеровой) орбиты и называемыми возмущающими факторами.

Например, для геостационарного ИСЗ достаточно учитывать только три фактора, приводящие к отклонениям ИСЗ от идеальной (Кеплеровой) орбиты, - влияние Солнца, Луны и нецентральности гравитационного поля Земли. Схема трех возмущающих тел: Земли S_з, Солнца S_c и Луны S_л, отображена на фиг. 4. Точками S₀ (x₀, y₀, z₀) и S_пр (x_пр, y_пр, z_пр) показаны координаты, а стрелками и - векторы скоростей ИСЗ в начальный момент времени t₀ и в прогнозируемый момент времени t_пр соответственно в геоцентрической декартовой системе координат O_зXYZ.

Алгоритм прогноза координат геостационарного ИСЗ описан в приложении Б, при этом процедуры расчета координат Солнца и Луны вынесены в приложение В.

Прогноз канонических параметров ИСЗ x_пр, y_пр, z_пр, в нужный (прогнозируемый) момент времени t_пр осуществляют на основе исходных данных x₀, y₀, z₀, в момент t₀ и рассчитанных силовых функций (ускорений) U₁, U₂, U₃, действующих на ИСЗ со стороны Земли, Солнца и Луны соответственно на основе решения системы дифференциальных уравнений (Б.2).

Прогноз координат других (не геостационарных) ИСЗ проводят на основе алгоритма аналогичного рассмотренному, с той лишь разницей, что учитывают больше факторов, приводящих к отклонениям ИСЗ от идеальной (Кеплеровой) орбиты. В качестве таких факторов, например, для ИСЗ на низких орбитах выступают: влияние сопротивления атмосферы Земли, светового давления, притяжения планет.

Система определения канонических параметров ИСЗ S₀ содержит одну приемо-излучающую НРТС K и N≥2 ИОРС I_n, n=1…N (см. фиг. 7) на позициях с известными координатами, излучающих сигналы заданного вида. НРТС кроме того принимает и обрабатывает ретранслированные ИСЗ радиосигналы.

Для обеспечения синхронности работы передатчиков НРТС и N ИОРС используют метки времени высокостабильного генератора частот.

Для определения координат ИСЗ х₀, y₀, z₀ в момент времени t₀ используют временные задержки радиосигналов НРТС и N ИОРС. Для получения таких задержек предварительно излучают в момент времени to тестовые радиосигналы НРСТ и N ИОРС. Далее принимают в НРТС ретранслированные ИСЗ тестовые радиосигналы в моменты времени и . Определяют в НРСТ задержки во времени Δt₀, Δt_n между излученными и принятыми тестовыми радиосигналами:

На основе полученных задержек во времени Δt₀, Δt_n однократно измеряют наклонную дальность от ИСЗ до НРСТ R_SK, и дополнительно, измеряют наклонные дальности от ИСЗ до каждой из N ИОРС :

где с=3×10⁸ м/с - скорость света в вакууме.

Полученные наклонные дальности R_SK и являются координатно-информативными параметрами ИСЗ, а поверхностями положения ИСЗ будут сферы с центрами, совпадающими с известными координатами НРТС и N ИОРС и радиусами, равными R_SK и соответственно. Координаты точки пересечения полученных сфер соответствуют искомым координатам ИСЗ х₀, y₀, z₀ в момент времени t₀.

Для одномоментного определения координат ИСЗ необходимо наличие трех сферических поверхностей положения А₁, А₂ и А₃ (см. фиг. 8), с радиусами R_SK, и соответственно а, следовательно, двух ИОРС I₁ и I₂ (см. фиг. 8). Дальнейшее увеличение количества ИОРС будет приводить к повышению точности определения координат ИСЗ и сопутствующим увеличением суммарной стоимости устройства, реализующего способ. Двузначность определения координат ИСЗ, обусловленная тем, что три сферы пересекаются в двух точках S₀ и (см. фиг. 8), устраняется априорными данными о координатах ИСЗ, в качестве которых выступают известные элементы Кеплеровой орбиты ИСЗ.

Для определения координат ИСЗ возможно использование и одной ИОРС. Однако, в таком случае время, необходимое для определения координат ИСЗ, будет значительно возрастать.

В качестве примера в приложении Г представлена аналитическая интерпретация алгоритма построения трех поверхностей положения в виде сфер с центрами, совпадающими с известными координатами НРТС и двумя ИОРС с радиусами, равными R_SK, и соответственно, а также нахождения точки пересечения таких сфер. В качестве выходных результатов представленного алгоритма выступают координат ИСЗ х₀, y₀ и z₀ в момент времени t₀.

Для определения ортогональных составляющих вектора скорости ИСЗ в момент времени t₀ используют частотные сдвиги радиосигналов НРТС и N ИОРС. Для получения таких частотных сдвигов предварительно излучают в момент времени t₀ тестовые радиосигналы НРСТ на частоте ƒ₀ и N ИОРС на частотах ƒ_n соответственно. Далее принимают в НРТС излученные НРТС и N ИОРС и ретранслированные ИСЗ тестовые радиосигналы на средних частотах и соответственно. Запоминают номиналы средних частот переданных ƒ₀, ƒ_n и принятых тестовых радиосигналов.

Далее вычисляют ортогональные составляющие вектора скорости ИСЗ в момент времени t₀ по известным координатам НРСТ, известным координатам N ИОРС, вычисленным координатам ИСЗ x₀, y₀, z₀, запомненным номиналам частот переданных ƒ₀, ƒ_n и принятых тестовых радиосигналов, предварительно заданной частоте сдвига ƒ_G рабочей частоты ИСЗ.

Для одномоментного и однозначного определения ортогональных составляющих вектора скорости ИСЗ необходимо использование двух ИОРС. Дальнейшее увеличение количества ИОРС будет приводить к повышению точности определения ортогональных составляющих вектора скорости ИСЗ, но с сопутствующим увеличением суммарной стоимости устройства, реализующего способ.

Для определения ортогональных составляющих вектора скорости ИСЗ возможно использование и одной ИОРС. Однако в таком случае время, необходимое для определения ортогональных составляющих вектора скорости ИСЗ, неизбежно будет существенно возрастать.

В качестве примера в приложении Д представлена аналитическая интерпретация алгоритма определения ортогональных составляющих вектора скорости ИСЗ на основе частотных сдвигов между излученными и принятыми тестовыми радиосигналами НРТС и N ИОРС для частного случая, когда N=2.

Окончательно, в качестве параметров орбиты ИСЗ принимают координаты ИСЗ х₀, у₀, z₀ и ортогональные составляющие вектора его скорости в момент времени t₀.

Имитационное моделирование заявленного способа определения параметров орбиты ИСЗ показывает существенное снижение времени, затрачиваемого на определение параметров орбиты ИСЗ (с двух суток, заявленных в способе прототипе, до нескольких минут), при этом ошибка определения координат ИСЗ не превышает 200 м. Так, например, при использовании двух ИОРС, время определения параметров орбиты ИСЗ не превышает пяти минут.

Источники информации

1. Шаргородский В.Д., Косенко В.Е., Садовников М.А., Чубыкин А.А., Мокляк В.И. Роль лазерных средств в обеспечении точности системы ГЛОНАСС. Исследования наукограда. 2013. №3-4 (6). С. 17-23.

2. Патент RU №2615634, опубл. 06.04.2017. Бюл. №10.

3. Патент RU №2313104, опубл. 20.12.2007. Бюл. №8.

4. Абалакин В.К. Астрономический календарь. Постоянная часть. - М.: Наука, 1981. - 704 с.

5. Машбиц Л.М. Компьютерная картография и зоны спутниковой связи, - 2-е изд., перераб и доп. - М.: Горячая линия - Телеком, 2009. - 236 с.

6. Выгодский М.Я. Справочник по высшей математике. - М.: Астрель; ACT, 2002. - 992 с.

Приложение А

Элементы Кеплеровой орбиты искусственного спутника Земли

Элементы невозмущенной орбиты ИСЗ включают шесть величин [4], определяющих Кеплерово движение: a_S, e_S, i, ω, Ω, τ₀. Такие элементы орбиты, как большая полуось a_S (см. фиг. 1, фиг. 2) и эксцентриситет e_S, характеризуют размеры и форму орбиты ИСЗ (см. фиг. 1, фиг. 2), а наклонение i орбиты ИСЗ к плоскости эклиптики, аргумент перицентра ω и долгота восходящего узла Ω (см. фиг. 1) характеризуют ориентацию орбиты в пространстве. Для задания элементов невозмущенной орбиты ИСЗ предварительно задают направление на точку весеннего равноденствия γ - ось ОХ' (см. фиг. 1) и определяют линию узлов (см. фиг. 1), как линию пересечения плоскости эклиптики с плоскостью орбиты ИСЗ. Временной параметр τ₀ представляет собой момент прохода ИСЗ через перигей П_S (см. фиг. 1, фиг. 2).

Невозмущенная орбита ИСЗ представляет собой эллипс (в частном случае окружность). В одном из фокусов эллипса находится притягивающий центр земного эллипсоида О_З (см. фиг. 1, фиг. 2). Линия, по которой проходит ось орбитального эллипса, называется линией апсид (см. фиг. 2). Перигей П_S и апогей A_S (см. фиг. 1, фиг. 2) - ближняя и дальняя от притягивающего центра О_З точки пересечения орбиты с линией апсид соответственно.

Положение ИСЗ S (см. фиг. 1, фиг. 2), обращающегося по эллиптической орбите, часто определяют в геоцентрической декартовой орбитальной системе координат О_ЗξζСη, в которой начало О_З совмещено с центром Земли, ось ξ направлена по линии апсид, ось ζ - по нормали к плоскости орбиты ИСЗ, ось η дополняет систему координат до правой (см. фиг. 1). В орбитальной системе координатами ИСЗ S являются его проекции на оси О_Зξ, О_Зζ, О_Зη, обозначаемые ξ_S, О_Зζ, η_S.

Кроме того, положение ИСЗ S часто удобно задавать в декартовой геоцентрической системе координат O_ЗXYZ, в которой начало координат О_{З совмещено с центром Земли, ось X направлена на Гринвич, ось Z - на север, ось Y дополняет систему координат до правой. Координатами ИСЗ в такой системе координат являются его проекции x, у и z на оси ОЗХ, OЗY и ОЗZ соответственно.}

Известны соотношения для определения расстояния от центра эллипса до фокуса c_S (см. фиг. 1) и фокального параметра p_S (см. фиг. 1):

c_S=a_Se_S;

Для дальнейших расчетов вводят истинную Θ_S, эксцентрическую E_S и среднюю M_S аномалии ИСЗ (см. фиг. 2), а также вспомогательную окружность с расположенными на ней вспомогательными точками S_E и S_M (см. фиг. 2).

Схема алгоритма расчета аномалий ИСЗ представлена на фиг. 3. Исходными данными (этап 1) являются: текущее время t, время прохождения перигея t_П, большая полуось a_S, эксцентриситет e_S, произведение гравитационной постоянной на массу Земли μ=398602 км³/с², допустимая ошибка определения эксцентрической аномалии δ_E, начальное значение счетчика итераций m=1.

На этапе 2 определяют среднюю аномалию:

Средняя аномалия связана с эксцентрической аномалией соотношением, называемым уравнением Кеплера:

Решение уравнения (А. 1) относительно E_S находят методом Ньютона-Рафсона (этапы 3-6). Так, E_S на m-й итерации равна:

где значение функции ƒ(Е_Sm-1) в соответствии с равенством (А. 1) рассчитывается по формуле

Производная ƒ'(E_m-1) определяется как

С учетом выражений (А. 3) и (А. 4) равенство (А. 2) преобразуют:

В качестве начального приближения процедуры (А. 5) может быть принято значение E_S0=M_S или E_S0=0 (этап 3).

При выполнении условия этапа 5 итерационный алгоритм завершается, а значение E_S принимается равным значению E_Sm на последней итерации (этап 7).

Приближенное значение эксцентрической аномалии E_S можно также найти с помощью степенного ряда:

Эксцентрическую аномалию определяют на этапе 8 по формуле

Результатами алгоритма (этапе 9) являются значения аномалий ИСЗ M_S, E_S и Θ_S.

Значение истинной аномалии ИСЗ Θ_S позволяет рассчитать координаты ИСЗ в геоцентрической декартовой орбитальной системе координат:

где ρ_S - радиус-вектор, рассчитываемый как

ρ_S=a_S(1-e_ScosE_S).

Переход к декартовым геоцентрическим координатам осуществляют центроаффинным преобразованием с помощью девятичленной матрицы:

С помощью элементов Кеплеровой орбиты определяют также ортогональные составляющие вектора скорости ИСЗ

Приложение Б

Алгоритм прогноза координат геостационарного искусственного спутника Земли

Схема алгоритма прогноза координат геостационарного ИСЗ представлена на фиг. 5, которая включает следующие этапы:

На этапе 1 вводят исходные данные, которыми являются канонические параметры ИСЗ х₀, y₀, z₀, в момент времени t₀, время прогноза t_пр, время шага t_ш, начальное значение счетчика h=0.

На этапе 2 рассчитывают количество шагов алгоритма в соответствии с выражением

при этом время шага выбирается так, чтобы N_ш было целым числом.

На этапе 3 значение счетчика h увеличивают на единицу.

На этапе 4 рассчитывают текущее значение времени по формуле

t_h=t₀+t_шh.

На этапах 5 и 6 рассчитывают координаты Солнца x_c, y_с, z_c и Луны х_л, у_л, z_л для текущего момента времени. Алгоритмы таких расчетов исследованы в астродинамике [4], они представлены в приложении В.

На этапе 7 рассчитывают силовые функции (ускорения) U₀, U₁, U₂, обусловленные влиянием на ИСЗ в текущий момент времени Земли, Солнца и Луны соответственно.

На этапе 8 с учетом данных, полученных на этапах 4-7, рассчитывают текущие координаты ИСЗ.

Далее этапы 3-8 алгоритма повторяют до выполнения условия h>N_ш (этап 9). В качестве выходных результатов алгоритма выступают канонические параметры ИСЗ x_пр, y_пр, z_пр, в момент времени t_пр (этап 10).

Общее уравнение движения ИСЗ в гравитационном поле Земли, образуемом действием центрального тела (Земли) и других небесных тел (Солнца, Луны, планет и др.), в векторной форме записывают как

где Г - радиус-вектор ИСЗ; Ψ_сум - геометрическая сумма сил, действующих на ИСЗ.

С учетом того, что на ИСЗ действуют только три небесных тела (Земля, Солнце и Луна), а также того, что масса ИСЗ пренебрежимо мала по сравнению с массой любого из указанных тел, компоненты выражения (Б. 1) можно представить в инерциальной геоцентрической ДСК O_ЗXYZ:

где - силовая функция, действующая на ИСЗ со стороны k-го тела; U_xk, U_yk, U_zk - компоненты ускорения, сообщаемые ИСЗ k-м небесным телом, k=1…3 (k=1 - Земля, k=2 - Солнце, k=3 - Луна).

Система нелинейных дифференциальных уравнений (Б. 2) преобразуется в систему из шести уравнений первого порядка:

Для численного решения системы уравнений (Б. 3) необходимо задать шесть начальных условий и выразить в явном виде компоненты векторов ускорения U_xk, U_yk, U_zk. В качестве начальных условий выступают канонические параметры ИСЗ х₀, у₀, z₀, для момента времени t₀, рассчитанные на основе измерений временных задержек и частотных сдвигов радиосигналов ИОРС.

Силовая функция Земли Ψ₁ в настоящее время хорошо изучена [5]. Компоненты вектора ускорения применительно к модели Земли, учитывающие полярное сжатие, можно представить в следующем виде:

где g₁=398,602×10¹² м³/с²; g₂=66 043,987×10⁶ м²; g₃=-1,492×10¹⁵ м³; g₄=-2,193×10²¹ м⁴; x_пр, y_пр, z_пр - прогнозируемые координаты ИСЗ в геоцентрической системе координат; - расстояние между центрами масс ИСЗ и Земли.

Для определения возмущающих влияний Солнца и Луны пользуются общей формулой силовой функции небесного тела [5]:

где K_k - гравитационный параметр (для Солнца K₂=13,2318×10¹⁹ м³/с², для Луны K₃=4891,0532×10⁹ м³/с²), равный произведению гравитационной постоянной на массу возмущающего тела; - расстояние между центрами масс ИСЗ и возмущающего тела; x_k, y_k, z_k - координаты центра масс возмущающего тела; - расстояние между центрами масс Земли и возмущающего тела.

Возмущающие ускорения согласно системе уравнений (Б. 2) определяются как частные производные (Б. 5) по координатам ИСЗ:

Решая систему дифференциальных уравнений (Б. 3) одним из известных численных методов (метод Рунге-Кутта, метод Адамса и др. [6]), с учетом выражений (Б. 4) и (Б. 6), при предварительном определении координат Солнца и Луны в соответствии с выражениями (В. 1) и (В. 7), получают канонические параметры ИСЗ x_пр, y_пр, z_пр, в момент времени t_пр.

Приложение В

Расчет координат Солнца и Луны

Расчет координат Солнца и Луны производят на основе алгоритмов, хорошо известных в астрономии. Исходными данными являются декретное время и декретная добавка, а выходными результатами - декартовы координаты указанных небесных тел.

Текущие координаты Солнца рассчитывают в соответствии с алгоритмом:

На этапе 1 рассчитывают Юлианские даты начала года JD₀ и заданного момента времени JD:

JD₀=2440000+0,5+g_г,м;

где g_г,м - табличное значение для текущих «г» и «м» [4]; «г» - год; «м» - месяц; - год, часы, минуты, секунды местного декретного времени; t_д - декретная добавка (например, для Москвы t_д=3 ч).

Юлианская дата представляет собой текущее время, исчисляемое в днях (сутках с десятичными долями), отсчитываемых от полудня 1 января 4713 года до н. э.

На этапе 2 определяют прямое восхождение для среднего Солнца α_сс, которое согласно С. Ньюкому равно:

где Т_ю - эпоха в Юлианских столетиях от среднего гринвичского полудня (всемирного времени) 1 января 1900 года, которую определяют как

Т_ю=(JD-2415020)/36525.

На этапе 3 рассчитывают уравнения времени η_с и прямого восхождения истинного Солнца α_с в соответствии с выражениями

η_с=1,925sin(α_cc+78)-2,375sin(2α_cc);

α_c=α_сс+η_с.

На этапе 4 вычисляют склонение Солнца:

δ_c=arctg(sinα_ctgε_экл),

где ε_экл=23,4371722 - угол наклонения плоскости эклиптики к плоскости Экватора, рассчитанный для k=2016 по формуле

ε_экл=23,4457889-(k-1950)×1,31×10^-4.

На этапе 5 определяют истинную аномалию Θ_c и расстояние от Земли до Солнца r_с:

Θ_c=2π(JD₀-JD))/365,25636;

где а_с=149597870 км - большая полуось орбиты Земли; е_с=0,0167 - эксцентриситет орбиты Земли.

На этапе 6 рассчитывают линейные координаты Солнца в инерциальной геоцентрической системе координат:

Рассчитанные линейные координаты Солнца по формуле (В. 1) в инерциальной геоцентрической системе координат позволят вычислить силовую функцию - ускорение U₂, придаваемое Солнцем искусственному спутнику Земли в момент времени t₀.

Луна является спутником Земли, и в отношении нее справедливо все то, что касается движения спутников вообще. Однако вопрос о ее движении рассматривается особо, так как Луна - самое близкое к Земле небесное тело, движущееся непосредственно вокруг Земли, и ее движение изучается весьма детально. При этом отклонения Луны от эллиптического движения, т.е. возмущения, велики, они гораздо больше, чем возмущения планет.

Видимое движение Луны изучается с древних времен. Луна перемещается по небу, как и Солнце, все время с запада на восток и описывает за один месяц большой круг. За сутки перемещение Луны составляет около 12…13°.

При детальном рассмотрении обнаруживаются весьма сложные особенности движения Луны. Отклонения от равномерного перемещения обусловлены, во-первых, тем, что Луна движется не по окружности, а в первом приближении по эллипсу с эксцентриситетом, равным е_л=0,005, а во-вторых, возмущениями от Солнца и от сжатия Земли.

В теории движения Луны выводятся формулы непосредственно для геоцентрической эклиптической долготы λ_л, широты β_л и расстояния r_зл от Луны до Земли. В упрощенном виде (если ограничиться наибольшими членами) эти формулы записывают следующим образом:

где L_л - осредненная (т.е. освобожденная от периодических возмущений) средняя долгота Луны в орбите; - осредненная средняя аномалия Луны; D_л - разность средних долгот Луны и Солнца; - осредненная средняя аномалия Земли, увеличенная на 180°; F_л - средняя долгота Луны, отсчитываемая от восходящего узла орбиты.

В формулах (В. 2) и (В. 3) отброшены по сравнению с точными (содержащими каждая по нескольку сотен членов) членами, которые имеют коэффициенты менее 100'' (0,028°), а в формуле (Д. 4) отброшены члены с коэффициентами меньше 0,0003.

Углы D_л, , F_л, L_л называются основными аргументами в теории движения Луны. Они выведены по многолетним наблюдениям в сочетании с теорией и определяются как

В основе преобразований эклиптических координат в экваториальные лежит сферический треугольник Леонардо - Эйлера P_мS_лP_э, (см. фиг. 6). Вершинами этого треугольника являются полюс Мира Р_м, полюс эклиптики Р_э и небесное тело, в данном случае Луна S_л.

В соответствии с теоремой косинусов сферической тригонометрии справедливы соотношения, связывающие геоцентрические эклиптические координаты Луны λ_л, β_л с геоцентрическими экваториальными координатами Луны - прямым восхождением α_л и склонением δ_л:

После преобразования выражений (В. 5) и (В. 6) получают формулы для расчета прямого восхождения α_л и склонения Луны δ_л:

δ_л=arcsin(sinβ_лcosε_экл+cosβ_лsinε_эклsinλ_л);

Линейные координаты Луны х_л, у_л, z_л в инерциальной геоцентрической системе координат вычисляют по формулам

Рассчитанные по формулам (В. 7) линейные координаты Луны в инерциальной геоцентрической системе позволят вычислить силовую функцию - ускорение U₃, придаваемое Луной искусственному спутнику Земли в момент времени t₀.

Приложение Г

Алгоритм построения трех сферических поверхностей положения ИСЗ и нахождение точки этих поверхностей

Для расчета координат ИСЗ на основе построения трех сферических поверхностей положения ИСЗ разработан алгоритм, схема которого представлена на фиг. 9.

На этапе 1 производят ввод исходных данных, в качестве которых выступают координаты НРТС x_K, y_K, z_K; координаты двух ИОРС I₁ и I₂ ; временные задержки между излученными НРТС, двумя ИОРС и принятыми НРТС тестовыми радиосигналами Δt₀, Δt₁, и Δt₂, предварительно определенные в соответствии с формулами (1) и (2); порог точности δ₀ расчета координат ИСЗ.

На этапе 2 рассчитывают координатно-информативные параметры - расстояния R_SK, и по формулам (3) и (4).

На этапе 3 выбирают, на основе элементов Кеплеровой орбиты ИСЗ, координаты опорной точки

На этапе 4 рассчитывают расстояния и при условии равенства координат ИСЗ координатам опорной точки по формулам

На этапе 5 рассчитывают невязки n₁, n₂ и n₃ как разницы между определенными на этапе 2 расстояниями R_SK, и и координатно-информативными параметрами R_SK, и , определенными на этапе 4 соответственно

На этапе 6 формируют систему линейных уравнений при разложении в ряд Тейлора функций R_SK, и с точностью до первых членов, где в качестве переменных выступают поправки к координатам Δх₀, Δy₀, Δz₀:

где частные производные, в свою очередь, рассчитываются согласно выражениям

Решая систему линейных уравнений одним из известных методов, например методом Крамера [6], получают поправки к координатам ИСЗ Δх₀, Δу₀, Δz₀.

На этапе 7 рассчитывают координаты новой опорной точки

Этапы 4-7 в совокупности составляют первую итерацию. Далее итерации повторяют, используя каждый раз новую опорную точку, полученную на предыдущей итерации. Количество необходимых итераций зависит требуемой точности определения координат ИСЗ. С точностью определения координат ИСЗ напрямую связан шаг итерации d_ш.

На этапе 8 определяют шаг итерации d_ш как расстояние между текущей и предыдущей опорными точками:

На этапе 8 сравнивают d_ш с порогом δ₀, задаваемом на этапе 1.

Необходимое число итераций, как правило, составляет 2…4. В качестве координат ИСЗ х₀, у₀ и z₀ выбирают значения координат опорной точки на последней итерации, вывод которых осуществляют на этапе 10.

Приложение Д

Алгоритм определения ортогональных составляющих вектора скорости искусственного спутника Земли

Для одномоментного определения ортогональных составляющих вектора скорости ИСЗ возможно использование НРТС и двух ИОРС, обозначенных на фиг. 10, - K и I_n соответственно, где для рассматриваемого частного случая n=1…2. Введены следующие обозначения: Θ_K - угол между вектором и направлением на НРТС; - углы между вектором и направлениями на n-ю ИОРС; - радиальная скорость ИСЗ относительно НРТС; - радиальные скорости ИСЗ относительно n-й ИОРС;

Предполагается, что координаты ИОРС - известны, а координаты ИСЗ x₀, у₀, z₀ - рассчитаны в соответствии с алгоритмом, представленном в приложении Г.

Соотношения средних частот ƒ₀, ƒ_n тестовых радиосигналов излучаемых НРТС и тремя ИОРС, а также средних частот тех же радиосигналов принимаемых в НРТС после их ретрансляции ИСЗ, имеют вид

где и - доплеровские сдвиги частот на входе ИСЗ за счет его сближения (удаления) с (от) n-й ИОРС; ƒ_G - заданная частота сдвига рабочей частоты ИСЗ; и - доплеровские сдвиги частот на выходе ИСЗ за счет его сближения (удаления) с (от) K.

Предполагают, что нестабильность генератора частот ИСЗ известна и компенсируется. Влияние других эффектов на изменение частоты, например гравитационный и релятивистский эффекты в рамках рассматриваемой задачи, пренебрежимо мало и поэтому не учитывают.

Доплеровские сдвиги частот и на входе ИСЗ рассчитывают следующим образом:

Доплеровские сдвиги частот и на выходе ИСЗ рассчитывают следующим образом:

Согласно теореме о скалярном произведении векторов [6]

Модуль вектора скорости ИСЗ равен:

а расстояния от НРТС и n-х ИОРС до ИСЗ R_SK и рассчитывают как

Выражения (Д. 1), (Д. 2) с учетом уравнений (Д. 3)…(Д. 9) для частного случая, когда n=1…2, преобразуют в систему линейных уравнений:

где коэффициенты при переменных и свободные члены равны:

Систему из трех линейных уравнений с тремя неизвестными (Д. 12) решают одним из известных методов, например, методом Крамера [6]. Результатом решения системы уравнений (Д. 12) выступают ортогональные составляющие вектора скорости ИСЗ

В общем случае, когда количество ИОРС N≥3, алгоритм определения ортогональных составляющих вектора скорости ИСЗ остается прежним, с той лишь разницей, что система уравнений (Д. 12) будет содержать более четырех уравнений. Тогда такую систему уравнений решают, например, методом наименьших квадратов [6].

Изобретение относится к области космонавтики, а именно к технике выполнения траекторных измерений и определения параметров орбиты искусственного спутника Земли (ИСЗ), и может быть использовано на наземных и бортовых комплексах управления полетом ИСЗ для точного определения текущих параметров движения ИСЗ. Достигаемый технический результат - высокое быстродействие. Указанный результат достигается за счет использования наземной радиотехнической станции (НРТС), размещенной на позиции с известными координатами x_K, y_K, z_K, и N≥2 излучающих опорных реперных станций (ИОРС), размещенных на позициях с известными координатами где n=1…N - номер ИОРС, при этом способ основан на излучении тестовых радиосигналов НРТС на средней частоте ƒ₀ и N ИОРС на средних частотах ƒ_n в момент времени t₀ и последующего приема этих тестовых радиосигналов после ретрансляции ИЗС с помощью НРТС в моменты времени и на средних частотах и соответственно, а также основан на измерении наклонных дальностей от ИСЗ до НРСТ R_SK и до каждой из N ИОРС , на вычислении координат ИСЗ х₀, у₀, z₀ по известным координатам НРСТ, известным координатам N ИОРС и измеренным наклонным дальностям R_sk, и , вычислении ортогональных составляющих вектора скорости ИСЗ по известным координатам НРСТ, известным координатам N ИОРС, вычисленным координатам ИСЗ х₀, y₀, z₀, запомненным номиналам средних частот переданных ƒ₀, ƒ_n и принятых тестовых радиосигналов, предварительно заданной частоте сдвига ƒ_G рабочей частоты ИСЗ, причем в качестве параметров орбиты ИСЗ принимают совокупность координат ИСЗ х₀, у₀, z₀ и ортогональных составляющих вектора его скорости в момент времени t₀. 10 ил., 5 приложений.

Способ определения параметров орбиты искусственного спутника Земли (ИСЗ), заключающийся в том, что размещают приемо-передающую наземную радиотехническую станцию (НРТС) на позиции с известными координатами x_K, y_K, z_K, выбирают начальные значения параметров орбиты ИСЗ, измеряют наклонную дальность от ИСЗ до НРТС, выбирают параметры орбиты ИСЗ на основе измеренной наклонной дальности, отличающийся тем, что дополнительно устанавливают на земной поверхности N≥2 излучающих опорных реперных станции (ИОРС) на позициях с известными координатами x_In, y_In, z_In, где n=1…N - номер ИОРС, причем перед измерением наклонной дальности R_SK от ИСЗ до НРТС предварительно излучают в момент времени t₀ тестовые радиосигналы НРСТ на средней частоте и N ИОРС на средних частотах соответственно, запоминают время излучения тестовых радиосигналов t₀, и номиналы средних частот тестовых радиосигналов и принимают в НРТС в моменты времени и , на средних частотах и соответственно ретранслированные ИСЗ тестовые радиосигналы излученные НРТС и N ИОРС, передают на НРСТ номиналы средних частот тестовых сигналов от N ИОРС, определяют задержки во времени Δt₀, Δt_n между излученными и принятыми тестовыми радиосигналами, после чего однократно измеряют наклонную дальность от ИСЗ до НРСТ R_SK, и дополнительно, измеряют наклонные дальности от ИСЗ до каждой из N ИОРС вычисляют координаты ИСЗ x₀, y₀, z₀ по известным координатам НРСТ, известным координатам N ИОРС и измеренным наклонным дальностям R_SK^, и вычисляют ортогональные составляющие вектора скорости ИСЗ по известным координатам НРСТ, известным координатам N ИОРС, вычисленным координатам ИСЗ х₀, y₀, z₀, запомненным номиналам средних частот переданных и принятых тестовых радиосигналов, предварительно заданной частоте сдвига рабочей частоты ИСЗ, а в качестве параметров орбиты ИСЗ принимают совокупность координат ИСЗ x₀, y₀, z₀ и ортогональных составляющих вектора его скорости в момент времени t₀.