Применение методов имитационного моделирования в задачах изучения движения околоземных космических аппаратов Текст научной статьи по специальности «Физика»

УДК 629.78

БОТ 10.18698/2308-6033-2016-07-1513

Применение методов имитационного моделирования в задачах изучения движения околоземных космических аппаратов

© Д. Г. Васильев1, В. В. Бетанов2

1МГТУ им. Н.Э. Баумана, Москва, 105005, Россия 2АО «Российские космические системы», Москва, 111250, Россия

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

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

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

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

времени, рассматривая различные сценарии его движения при изменении внешних и внутренних условий, применяют методы имитационного моделирования.

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

Для уточнения идентифицируемых параметров при наличии ошибок измерений предлагается использовать имитационно-моделирующий комплекс (ИМК).

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

Общие вопросы имитационного моделирования. Имитация дает возможность экспериментировать с большим числом переменных и правил принятия решений, с более сложными моделями. Можно выделить следующие преимущества этого метода [1].

1. Имитация позволяет исследовать сложные внутренние взаимодействия в рассматриваемой системе.

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

3. Детальное наблюдение за имитируемой системой позволяет лучше понять ее и разработать предложения по ее совершенствованию.

4. Имитация позволяет изучать динамические системы в реальном или приведенном времени.

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

• средства обучения;

• инструмента прогнозирования;

• средства постановки эксперимента.

В общем виде структуру модели математически можно представить так [2]:

Е = / ( X, У], ^ ),

где Е — результат действия системы; х1 — переменные и параметры, которыми можно управлять; У] — переменные и параметры, которыми управлять невозможно; И1 — случайные факторы, вносящие искажения в результаты моделирования; / — функциональная зависимость между параметрами х1 и У], определяющая величину Е.

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

• компонентов;

• переменных;

• параметров;

• функциональных зависимостей;

• ограничений;

• целевых функций.

В работе для создания модулей и блоков был использован математический пакет МЛТЬЛВ, его структура позволяет эффективно сочетать основные подходы к созданию модели: аналитический и имитационный.

Прогнозирование движения КА. Основные характеристики задачи определения параметров движения КА, прогнозирование их для заданных моментов времени непосредственно зависят от свойств математической модели движения (ММД), рационального построения ее функциональной структуры, своевременности получения и точности навигационно-баллистической информации.

Для определения параметров движения по измерениям текущих навигационных параметров необходимо иметь возможность вычисления расчетных значений вектора состояния в заданный момент времени. Для выполнения указанной процедуры используют ММД КА. Это совокупность дифференциальных уравнений движения КА с описанием математических моделей сил, действующих на КА в полете, и метода решения этой системы уравнений [3].

В ИМК используются две основные ММД КА:

1) ММД с системой дифференциальных уравнений, записанных в орбитальных (оскулирующих) кеплеровых элементах;

2) ММД в неособенных Х-переменных.

Здесь будет исследовано возмущенное движение КА. Причины возникновения основных возмущающих сил таковы:

• нецентральность поля тяготения, обусловленная сплюснутостью планеты и неравномерным распределением масс внутри нее (аномалии силы тяжести);

• сопротивление атмосферы;

• влияние Солнца и Луны;

• давление солнечного света.

Вид уравнений движения центра масс КА определяется выбранной системой координат и составом учитываемых действующих сил (на практике редко рассматривают все возмущающие силы одновременно). Для близких к Земле КА, как правило, учитываются нецентральность поля сил и сопротивление воздуха. В качестве характеристик возмущений будем рассматривать изменения кеплеровых элементов г, й, ш, е, р орбиты за один виток [4].

В соответствии с математической моделью КА движется над поверхностью Земли, имеющей форму эллипсоида вращения с экваториальным радиусом Яе = 6378,160 км и полярным радиусом Яр =

= 6356,863 км (эллипсоид Красовского) [5]. Аппарат движется относительно Земли под действием силы тяготения, полной аэродинамической силы и сил, обусловленных нецентральностью системы отсчета. Система дифференциальных уравнений, записанная в орбитальных (оскули-рующих) элементах, принимаемых за переменные, имеет вид [4]:

Лй = г

л ТмР

Лг

БШ и

Ж;

БШ г

л ТмР

Жсоб и;

ЛР = 2.1 РгТ; Лг

Л Х1 Лг

Л Х2 Лг

= - ¿соб и + Т

\+г •>

V

р)

Бт и Л— (Х1Т - Xг вт и) р

= . — ¿Бт и + Т

1+г ^

V

р)

соб и Л— (X2Т + Х1 Жс\% I Бт и) р

Лю = \р Лг М д

- 5 + Т

^ г ^ бш 1 1 + -

V

р

- ж—г и)

е р

Здесь в левых частях уравнений стоят производные от шести выбранных элементов орбиты: й, г, р, Х1, X2 и ш, в правых частях —

проекции 5, Т, Ж возмущающей силы, отнесенной к массе КА, на орбитальные оси г, Ь, п.

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

О = О0 + 80; р = р0 + 8р; Х2 = Х20 + 8Х2;

/ = /0 +8/; Х1 = Х10 + 8Х^ ш=ш0 + 8ш,

где 00, /0, р0, Х10, Х20, ш0 — значения оскулирующих элементов в начальный момент времени; 80, 8/, 8р, 8X1, 8X2, 8ш — изменения этих элементов, вычисленные интегрированием представленной выше системы уравнений.

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

йО

г у БШ и

йи цр Бт /

Ж;

й = г у йи цр

Жсоб и;

йР=21 г 3Т;

йи ц

й Х1 г 2 у

йи ц

й Х2 г 2 у

йг ц

(

- ¿соб и + Т

Л

1 + -

С

¿^п и + Т

1+г

Р )

\

Р )

Бт и +— (Х1Т - XЖсЩ / вт и) Р

соб и +— (X2Т + / Бт и)

Р

йю йи

г 2 у

ц

- 5

соб д

(

+ Т

Л

1 + -

Р )

Бт д ттг г

--Ж — ^ / бш и)

е р

Здесь

у = -

1

-; е =

+ X 2

+ Ае; г = ■

1--Ж ^ / Бт и

цр

1 + е соб и

- + Аг;

3

Т, Ж — проекции возмущающих ускорений на оси орбитальной системы координат с центром, совпадающим с центром масс КА: на радиус-вектор КА — на ось, перпендикулярную радиусу-вектору в плоскости орбиты, — Т; на бинормаль к плоскости оскулирующей орбиты — Ж [6],

£ = (з вт 2 * вт 2 и - 1);

£ 2 Т = —^т I б1П2м;

г

£

Ж = —4 8т2/ Бт и. г

На высоте более 150...200 км атмосфера сильно разрежена и поэтому оказывает малое сопротивление движущемуся КА [7]. Но сила сопротивления является постоянно действующей силой, поэтому, несмотря на малость, может значительно изменить элементы орбиты за достаточно большой промежуток времени.

Влияние сопротивления атмосферы на движение КА оценивают по характеру поведения и изменениям оскулирующих элементов орбиты [6]. Без учета атмосферы приближенные значения вековых возмущений некоторых элементов орбиты за один виток определяют по следующим зависимостям:

Дг = -4л£бРгср; Д ^ = ^бр^/^; Дв = 12^4; ДТ = -12л2 £б^/Т>,

где Дг, ДУу, Дв, ДТ — вековые возмущения соответственно среднего радиуса, продольной скорости, смещения вдоль орбиты и периода обращения КА.

Баллистический коэффициент КА массой т (площадь миделевого сечения БП1) определяют по формуле [4]

С £

о _ ха т

2m

где Cxa — коэффициент лобового сопротивления; р — плотность воздуха на рассматриваемой высоте полета КА; гср — средний радиус орбиты.

Под влиянием сопротивления атмосферы при движении КА по эллиптической орбите происходят вековые возмущения эксцентриситета в и фокального параметра p, при этом первоначальная орбита со

временем приближается к круговой. Период обращения монотонно уменьшается, средняя скорость полета возрастает. На рис. 1 представлен результат моделирования орбиты КА, а также приведены параметры полученной орбиты.

Рис. 1. Результаты моделирования формы и параметров орбиты

На рис. 2 показаны графические зависимости кеплеровых элементов от аргумента широты u, взятого в качестве параметра интегрирования. Полученные параметры соответствуют полету Международной космической станции.

На основе полученных элементов орбиты Q, i, p, ш, e определим координаты и составляющие скорости КА в инерциальной (абсолютной) системе координат в зависимости от истинной аномалии $ по следующим выражениям [6]:

х = r (cos u cos Q - sin u sin Q cos i); y = r (cos u sin Q + sin u cos Q cos i); z = r sin u sin i;

Vx = Vr (cos u cos Q - sin u sin Q cos i) - Vn (sin u cos Q + cos u sin Q cos i);

Vy = Vr (cos u sin Q + sin u cos Q cos i) - Vn (sin u sin Q - cos u cos Q cos i);

Vz = Vr sin u sin i + Vn cos u sin i.

Здесь Vr, Vn — радиальная и нормальная составляющие скорости, Vr = pe sin 9; Vn = p (l + e cos 9).

6792,4 6792,3 6792,2 6792,1 6792

\AJ\ <\S\f

О

10 u

15

20

2 1

св

oo n o "

8-1

-2

ir—^ ■Г"4 —

Г г

1

1 1

0

0,9008 0,9006 0,9004 0,9002 0,9000 0,8998

10 u

15

20

Рис. 2. Результаты моделирования зависимости кеплеровых элементов от аргумента широты

При использовании неособенных переменных вместо традиционных прямоугольных координат или кеплеровых элементов не возникает заметных осложнений, так как преобразование Х-переменных к кеплеровым элементам выполняется с помощью простых соотношений [3]. При обосновании спектра возмущающих сил (ускорений), который необходимо учитывать, описывая движение КА с помощью ММД, следует исходить из требований к точности прогнозирования движения тех или иных КА при одновременном выполнении требований к оперативности проведения баллистико-навигационных расчетов.

Неособенные переменные выражаются через кеплеровы элементы орбиты так:

X0 = a; X3 = sin i / 2 cos Q;

X = e cos (со + Q); X4 = sin i /2 sin Q;

X2 = e sin (со + Q); X5 = $ + ш + Q.

Здесь a — большая полуось орбиты; e — эксцентриситет; i — наклонение плоскости орбиты; Q — долгота восходящего узла; ш — аргумент перигея; $ — истинная аномалия.

Обратный переход от неособенных переменных к кеплеровым элементам орбиты или прямоугольным координатам осуществляется по формулам, приведенным в таблице.

Формулы перехода от неособенных Х-переменных к кеплеровым элементам и прямоугольным координатам

Кепле-ровы элементы Формула перевода Дополнительные параметры Прямоугольные координаты и скорости

a Xo P1 = X1 cos X5 + X2 sin X5 P2 = X1 sin X5 - X2 cos X5 P3 = X3 cos X5 + X4 sin X5 P4 = X3 sin X5 - X4 cos X5 P, = C 1 + P Ps = nXo n 4 03 C = V 1-X2 -X2 с = V 1-X22 -X4 У = P2 P6 r С у P6 t1 + P1 ) " С r = CP, a cos y1 = cos X 5 + 2X 4 P4 cos y2 = sin X5 - 2X3P4 cos y3 = 2cP4 x = r cos У!

e л/^2 + X 2 , = r cos у2

i 2arctg^X 2 + X 2 z = r cos у3

ю ( X 2X3 _ 4 ) arctg 1 1 ^ X1X3 + X2X4 J Vx =Vr cos yi - Vu cos yi

Q arctg fe ] V, = Vr cos у2 - Vu cos у2

3 arctg (P ] V = Vr cos уз - Vu cos уз

Система дифференциальных уравнений, описывающих движение КА в неособенных переменных, преобразованная к независимой переменной Х5, имеет следующий вид [8]:

dX0 2cP6

d Х5 nP5 у

SP2 +

Tс

Р

5

d = 1 d Х5 у

С

л

ScC2 ЕР^ + Т^с [Х1 + (2 + P1) соб Х5 ] - ЖСХ 2 P4 Р5

d

d Х5 у

л

- £сС2 -СО--1 + ТСс [ X 2 + (2 + Р1) мп Х5 ] - ЖСХ1Р4 —5

л

У

d= ЖС(соБ Х5 - Р3Х3);

d Х5

d Х4 d

2у

ЖС(б1п Х5 - Р3Х4) : 2у ;

=р6Сс d Х5 Р5 у

о П—6С

Здесь у = —6т

Р3

+ ЖР4С; п — среднее движение; Р1 - Р6, С, с — см. таб-

лицу.

Использование неособенных переменных позволяет описать систему дифференциальных уравнений таким образом, что правые части не вырождаются при значениях наклонения 0.90° и 90.180° и значениях эксцентриситета 0.1.

Применением в качестве независимой переменной Х5 (аналог аргумента широты) вместо времени достигается устойчивость получения численных решений интегрируемой системы дифференциальных уравнений для круговых и высокоэллиптических орбит. Возможности ММД позволяют осуществить баллистико-навигационное обеспечение всех типов КА с высотами орбит 120.40 000 км в штатном и нештатном режимах, включая падающие и спускаемые КА.

Структура ИМК. Комплекс состоит из следующих блоков (рис. 3):

• планирования модельного эксперимента;

• моделирования внешних воздействий;

• задания вида уравнений движения КА;

• расчета прогноза движения;

• расчета и построения трассы полета КА;

• задания координат наземных измерительных пунктов;

• расчета видимости КА наземными командно-измерительными пунктами;

• отображения зон видимости КА и светотеневой обстановки;

Планирование модельного эксперимента

Запуск исполняемого файла

Выбор вида дифференциальных уравнений движения

Ввод начальных условий

Прогноз движения КА

Модель внешних воздействий

Расчет трассы полета КА

Отображение трассы полета на развертке земной поверхности

V В оскулирующих элементах

V В неособенных ^.-переменных

V Атмосфера

\/ Гравитационное поле

V Солнце, Луна

V Давление солнечного света

Элементы орбиты

X

Задание координат наземных измерительных пунктов

Шумы, задержка сигнала

Расчет зон видимости КА наземными командно-измерительными пунктами

Отображение светотеневой обстановки на земной поверхности

Начало сеанса измерений навигационных параметров

Массив измерений навигационных параметров

Обработка измерений навигационных параметров

Отобрг на разверт поверх зон видии 1жение ке земной лости лости КА

Расчет времени сеанса связи

Построение орбиты, графические зависимости, сохранение результатов

Точность оценивания вектора состояния КА

Уточнение параметров

Рис. 3. Структурная схема разрабатываемого ИМК

• моделирования шумов измерений, а также задержки навигационного сигнала;

• накопления таблиц узловых значений;

• моделирования и обработки измерений навигационных параметров;

• точности оценивания параметров движения КА;

• уточнения параметров модели;

• вывода результатов.

Планирование модельного эксперимента. При планировании преследуют две основные цели [2]:

1) сокращение общего объема испытаний модели при соблюдении требований к достоверности и точности получаемых результатов;

2) повышение информативности каждого эксперимента в отдельности.

Поиск плана эксперимента производится в так называемом факторном пространстве, т. е. в множестве внешних и внутренних параметров модели, значения которых нужно контролировать в ходе подготовки и проведения модельного эксперимента.

План эксперимента строится относительно совокупности выходных параметров У, которые называются наблюдаемыми переменными. При этом предполагается, что значения наблюдаемых переменных, полученных в ходе эксперимента, складываются из двух составляющих:

у = Ф ( х ) + е ( х ),

где ф (х) — функция отклика (неслучайная функция факторов); е (х) — ошибка эксперимента (случайная величина); х — параметры модели.

Очевидно, что у является случайной переменной, так как зависит от случайной величины е (х). Дисперсия Б (у) наблюдаемой переменной, которая характеризует точность измерений, равна дисперсии ошибки эксперимента: Б (у ) = Б ( е ).

При таком рассмотрении можно выделить два основных варианта постановки задачи планирования имитационного эксперимента:

1) из всех допустимых выбрать такой план, который позволил бы получить наиболее достоверное значение функции отклика ф (х) при

фиксированном числе опытов;

2) выбрать такой допустимый план, при котором статистическая оценка функции отклика ф (х) может быть получена с заданной точностью при минимальном объеме испытаний.

В данном блоке также решается задача оптимального планирования измерений, которая включает:

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

• оптимизация временных программ измерений (определение общего числа измерений).

Моделирование внешних воздействий. В ходе изучения движения КА необходимо включать в процесс моделирования действие следующих внешних факторов:

• влияние атмосферы на движение КА;

• влияние гравитационного поля Земли;

• гравитационное действие Солнца и Луны;

• давление солнечного света.

На КА, движущиеся на высоте 150...1500 км, заметно влияет сопротивление атмосферы. Для учета торможения спутника в верхних слоях атмосферы Земли во внимание принимают вектор аэродинамических сил, направленный противоположно вектору относительной скорости спутника [9]:

Ра = - И )|

и и,

где 8Ь — баллистический коэффициент; р( И) — плотность воздуха

на высоте И над поверхностью Земли.

Приближенное значение баллистического коэффициента определяют по формуле

= 106 А

ч

т

2

где А3 — площадь поперечного сечения КА, м ; т8 — масса КА, кг.

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

„ 1Л-13 Г И - 200 ^ р(И) = 2-10 13 ехр^--^.

Текущую высоту полета (км) вычисляют по формуле

И = г - г0

2

1 -а — г

где г0 = 6378,14 км — экваториальный радиус Земли; а = 1 / 298,257 —

сжатие Земли.

Модуль расстояния (км)

г = л] х2 + у2 + г2.

Модуль вектора скорости (км/с)

V = у1 Ух2 + ГуУ + V2.

Затем вычисляют вектор ускорения Га, обусловленного торможением атмосферы, по следующим зависимостям:

(^) х = ; (^) у = -Бър^у;

(^ )г = -Бър^я.

Геопотенциал в земной системе координат имеет вид [9]

и = [1+ £ {Т Г рп ) (^ ф) \СпЪ ™ кX + ^ 8Ш кХ] [. Обозначим:

г I п=2 к=0^ г

и = ¿т.

и 0 - ; г

и1 = £ ТУпк;

п=2 к=0

ипк =

= >Гг0] р() (в1и ф) [Спк сов кX + Б* 81п кX].

г V г

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

. х ди

К = -.М—+-д-1;

г дх

г3 Зу

г Зи1

г г3 дг

•^тях п

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

Если высота орбиты менее 1000 км, то возмущения от Луны и Солнца можно не учитывать. Для более высоких орбит эти возмущения учитывать необходимо.

Вектор ускорения ¥м, обусловленного притяжением Луны, рассчитывают по следующему алгоритму [10].

1. Вычисляют прямоугольные координаты Луны для заданного момента времени Хм, Ум, 1м.

2. Вычисляют модуль расстояния от центра Земли до Луны:

гм = V ХМ + уМ + 'м.

3. Вычисляют модуль расстояния от Луны до КА:

г0 =^(Хм - *)2 + (м - У)2 + ('м - 2)2.

4. Вычисляют компоненты вектора ускорений Гм:

(Гм )х = - +тм Щг;

гм

3

г0

Ум - У

( 3

г0

'м - г

гм

(Гм ) = >м±м^ - №м\

г0 гм

Вектор ускорения Г*, обусловленного притяжением Солнца,

рассчитывают по следующему аналогичному алгоритму.

1. Вычисляют прямоугольные координаты Солнца для заданного момента времени X*, У*, .

2. Вычисляют модуль расстояния от центра Земли до Солнца:

г8 = ^X* + УI + 2*.

3. Вычисляют модуль расстояния от Солнца до КА:

= у1( - X)2 + ( - У)2 + ( - 2)2 .

г0

4. Вычисляют компоненты вектора ускорений Г|:

(Г| )*=- ;

* Г03 г1

(ъ) = - мА;

(Р ) = - !тг

Вектор ускорения ¥Р, обусловленного действием светового давления, рассчитывают по следующему алгоритму [10].

1. Вычисляют приближенное значение коэффициента эффективного отражения (км/с2):

-3 А8

Сгф = 10 Р0кг ,

т3

где кг — эмпирический коэффициент отражения (1 < кг < 2).

2. Вычисляют прямоугольные координаты Солнца для заданного момента времени Х5, У5, .

3. Вычисляют модуль расстояния от Солнца до КА Г).

4. Вычисляют компоненты вектора ускорения р:

(Р ) = с

(

гвА

1,4959787 -10й

(рр )у = с

гвА

(рр ) = с

х - X

5 .

1,4959787 -108

у - Ъ

5 .

гвА

1,4959787 -108

г -

Расчет трассы полета и зон радиовидимости КА наземными измерительными пунктами. Задача построения трассы, одна из наиболее важных, связана с обеспечением планирования полета. При решении этой задачи с определенным шагом по времени вычисляют высоту полета к и географические координаты — широту ф и долготу X, которые затем наносят на развертку земной поверхности.

По известному прогнозу движения КА расчет трассы (рис. 4) проводят по следующим формулам [7]:

г = ^ х2 + у2 + г2; г =>/ х2 + у2;

ф = агсБт

г

>/(1 -«)2 Г12 + г2

X = агй£—;

X

-90

Рис. 4. Пример расчета трассы полета КА

Задача расчета зон радиовидимости КА чрезвычайно важна для управления движением и планирования работы. При расчете на каждом витке орбиты для каждого измерительного пункта рассчитывают время входа КА в зону прямой радиовидимости и выхода из нее, а также суммарное время пребывания его в этой зоне. Кроме того, рассчитывают время достижения КА минимального расстояния от станции слежения до КА и угол места. Эти данные характеризуют условие радиовидимости [4].

Высота полета КА

И = Я

Бт (а' + Р) 1

бш а

Радиус зоны видимости

И

Р = агсБт

Я

+ 1 I Бт а'

-а' = 180°-а'-5',

5 .Г Я+и . '

5 = агоБт I-бш а

I Я

Наклонная дальность

ё = д/Я2 +(Я + И)2 -2(Я + И)Ясобр.

Зону видимости удобно представлять проекцией ее границы на поверхность Земли. В этом случае пересечение трассы КА с границей зоны видимости дает точки начала и конца этой зоны.

Для расчета границы зоны видимости должны быть заданы координаты измерительного пункта (широта и долгота) и известны элементы орбиты е, 1, р, ш.

Моделирование и обработка измерений навигационных параметров. В контексте разрабатываемого ИМК введем следующие ограничения на проводимые измерения:

• все навигационные измерения проводятся в виде навигационных сеансов стандартной длительности;

• частота измерений в течение каждого сеанса постоянна;

• навигационный сеанс реализуется в том случае, если выполняется условие радиовидимости КА в зоне действия соответствующего измерительного пункта;

• затраты на измерения определяются числом проводимых измерений текущих навигационных параметров;

• в качестве критерия, определяющего качество работы измерительного пункта, рассматривается апостериорная дисперсия параметра, характеризующего отклонение положения КА от истинного;

• информация, полученная с измерительных пунктов, обрабатывается с помощью метода наименьших квадратов и фильтра Калмана.

Модель измерений базируется также на методике накопления таблиц узловых значений (ТУЗ), неизменных относительно орбиты КА и представляющих собой совокупность шестимерных векторов, содержащих переменные и время в точках интегрирования уравнений движения [3]. Исходными данными для построения ТУЗ являются начальные условия движения КА, время начала и конца интервала ТУЗ. При построении ТУЗ для очередного интервала могут быть использованы как исходные начальные условия, так и параметры движения одной из узловых точек предыдущего интервала.

В каждом навигационном сеансе реализуется несколько навигационных измерений. Как правило, после каждого навигационного сеанса имеется пауза (период восстановления). Для упрощения модели введем некоторые ограничения на процесс проведения измерений:

• начало первого навигационного сеанса в каждой зоне радиовидимости совпадает с началом зоны;

• в отдельных случаях период восстановления может быть равен нулю, т. е. сеансы внутри зоны будут следовать один за другим;

• все навигационные измерения распределяются внутри одного навигационного сеанса равномерно.

Для того чтобы обеспечить заданную точность позиционирования КА, нужно управлять точностью процесса обработки поступающей информации, т. е. фильтрацией. Помимо этого необходимо контролировать процесс определения апостериорной матрицы расширенного вектора состояния КА, который включает не только координаты и скорости, но и дополнительные параметры, в том числе систематические ошибки навигационных измерений [11].

Ввиду дискретного характера навигационных измерений процесс

определения априорной (предсказанной) р и апостериорной (скор*

ректированной) р матриц ковариаций расширенного вектора состояния КА будет описываться с помощью соотношений дискретного фильтра Калмана [11]:

р* =( р-1 + нТ Б~1н-)-1-

* \ * I г|г I I '

р = А*—1р*—1АТ-1,

где А*-1 — фундаментальная матрица уравнений движения КА,

дУ*

* = 0, 1,.., Ы8; Н1 = ——, ] = 1, 2, ..., п — матрица частных произ-

дХ1

водных измеренных (навигационных) параметров по компонентам расширенного вектора состояния КА; — матрица ковариаций

случайных ошибок измерений.

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

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

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

у () = а ((г), Хтс ^), /, с (;)))(г),

где У () — вектор навигационных измерений (I х1); О (■) — вектор-функция ((х1), определяющая зависимости между измеряемыми величинами, компонентами расширенного вектора состояния КА и координатами наземного измерительного пункта; С (^) — вектор систематических ошибок измерений, который в дальнейшем подлежит

оцениванию // xl); ^/t) — вектор случайных ошибок измерений (/xl); X(t) — вектор состояния КА /nxl); XTS /1) — вектор, определяющий положение наземного измерительного пункта.

Компоненты вектора Y /1) включают:

• наклонную дальность;

• скорость изменения наклонной дальности;

• наклонение КА в топоцентрической системе координат;

• прямое восхождение КА в топоцентрической системе координат.

Как было сказано выше, задача в такой постановке интерпретируется как задача управления точностью калмановской фильтрации. Такой подход позволяет учесть стохастический характер движения КА, корреляцию ошибок измерений, ограничения на максимальную частоту измерений и др.

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

Возможны также и другие обобщения, касающиеся наблюдений:

• вырождение матрицы интенсивности шумов;

• отсутствие информации в канале наблюдения;

• повышение порядка формирующего фильтра.

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

Рассмотренная общая концепция создания ИМК и его основных элементов позволит в дальнейшем развивать и совершенствовать эти элементы, а также включать в состав программной части новые алгоритмы.

ЛИТЕРАТУРА

[1] Труды Всероссийской научной конференции «Проектирование научных и инженерных приложений в среде MATLAB». Секция 5. Моделирование в Simulink. Москва, ИПУ РАН, 2004, 1955 с.

[2] Гультяев А.К. MATLAB 5.2. Имитационное моделирование в среде WINDOWS: практическое пособие.

[3] Байрамов К.Р. Подход к разработке обобщенной технологической модели решения некорректных задач определения движения космических аппаратов по измерениям текущих навигационных параметров. Электроника, измерительная техника, радиотехника и связь. Доклады ТУСУР, 2010, № 2 (22), ч. 2, декабрь.

[4] Нариманов Г.С., Тихонравов М.К., ред. Основы теории полета космических аппаратов. Москва, Машиностроение, 1972, 608 с. (Серия «Основы теории полета и проектирования космических аппаратов»).

[5] Лазарев Ю.Н. Управление траекториями аэрокосмических аппаратов. Самара, Самар. науч. центр РАН, 2007, 274 с.

[6] Иванов Н.М., Лысенко Л.Н. Баллистика и навигация космических аппаратов. 2-е изд. Москва, Дрофа, 2004, 544 с.

[7] Бордовицына Т.В., Авдюшев В.А. Теория движения искусственных спутников Земли. Аналитические и численные методы. Томск: Изд-во Том. ун-та, 2007, 178 с.

[8] Лысенко Л.Н., Бетанов В.В., Звягин Ф.В. Теоретические основы балли-стико-навигационного обеспечения космических полетов. Москва, Изд-во МГТУ им. Н. Э. Баумана, 2014, 518 с.

[9] Аксенов Е.П., Чазов В.В. Модель движения ИСЗ. Главные проблемы. Основные алгоритмы. Москва, Наука, 2007, 188 с.

[10] Чазов В.В. Прогноз орбитального движения космического аппарата. Численная модель. Научно-технический отчет.

URL: http: //vadimchazov.narod.ru/text_pdf/comalg.pdf (дата обращения 05.06.2016).

[11] Малышев В.В., Красильщиков М.Н., Бобронников В.Т., Нестеренко О.П., Федоров А.В. Спутниковые системы мониторинга. Анализ, синтез и управление. В.В. Малышев, ред. Москва, Изд-во МАИ, 2000, 568 с.

Статья поступила в редакцию 18.05.2016

Ссылку на эту статью просим оформлять следующим образом:

Васильев Д.Г., Бетанов В.В. Применение методов имитационного моделирования в задачах изучения движения околоземных космических аппаратов. Инженерный журнал: наука и инновации, 2016, вып. 7.

http://dx.doi.org/10.18698/2308-6033-2016-07-1513

Статья подготовлена по материалам доклада, представленного на XL Академических чтениях по космонавтике, посвященных памяти академика С. П. Королева и других выдающихся отечественных ученых — пионеров освоения космического пространства, Москва, МГТУ им. Н.Э. Баумана,

26-29 января 2016 г.

Васильев Дмитрий Геннадьевич — аспирант заочной формы обучения кафедры «Динамика и управление полетом ракет и космических аппаратов» МГТУ им. Н. Э. Баумана. e-mail: vasiliev_dmitry83@list.ru

Бетанов Владимир Вадимович — д-р техн. наук, профессор, советник директора ОАО «Российские космические системы», профессор кафедры «Динамика и управление полетом ракет и космических аппаратов» МГТУ им. Н. Э. Баумана. Член-корреспондент Российской Академии ракетных и артиллерийских наук. Имеет более 300 научных работ и изобретений в области баллистики, динамики полета, управления движением ракет и космических аппаратов, пять монографий, два учебника, 30 учебных и научно-методических пособий, соавтор Военного энциклопедического словаря РВСН. e-mail: vlavab@mail.ru

Use of simulation modeling methods in problems of studying the near-Earth spacecraft motion

© D.G. Vasiliev1, V.V. Betanov2

'Bauman Moscow State Technical University, Moscow, 105005, Russia 2Joint-stock company Russian Space Systems, Russia

The paper presents the concept of creating a simulating and modeling complex and studies the structure of the algorithmic software. In the research we examine the interaction of the main elements of the complex. The modeling complex is designed for carrying out and processing the measurement sessions during the spacecraft orbital motion, as well as working out methods for determination and assessment of the vector of the spacecraft state according to the observations of its orbital motion using the distorted information. Moreover, we specified the model parameters according to the results of observations from the ground-based measuring stations. The concept of the modeling and simulating complex under consideration gives an insight into the work of its main subsystems and allows us to improve the algorithms of the individual elements and the whole system as well.

Keywords: spacecraft, simulation modeling, modeling complex, the prognosis of orbital motion, measurement processing of current navigation parameters.

REFERENCES

[1] Trudy vserossiyskoy nauchnoy konferentsii "Proektirovanie nauchnykh i izhe-nernykh prilozheniy v srede MITLAB". [Proceedings of all-Russian science conference Design of science and engineering applications in MITLAB medium]. Sektsiya 5. Modelirovanie v Simulink [Section 5. Modeling in Simulink]. 2004.

[2] Gultyaev A.K. MATLAB 5.2. Imitatsionnoe modelirovanie v srede Windows: prakticheskoe posobie [MATLAB 5.2. Simulation modeling in Windows medium: manual]. Moscow, Nauka Publ., 1990, 285 p.

[3] Bairamov K.R. Electronica, izmeritelnaya tekhnika, radiotekhnika i svyaz. Doklady TUSUR — Electronics, measurement technology, radio engineering and communications. Proceedings of TUSUR University, December 2010, no. 2 (22), part 2.

[4] Narimanov G.S., Tikhonravov M.K., ed. Osnovy teorii poleta kosmicheskikh apparatov. [Fundamdntals of spacecraft flight theory]. Seriya "Osnovy teorii poleta i proektirovaniya kosmicheskikh apparatov" [Ser. Fundamentals of space craft flight and design theory]. Moscow, Mashinostroenie Publ., 1972, 608 p.

[5] Lazarev Yu.N. Upravlenie traektoriyami aerokosmicheskikh apparatov [Spacecraft trajectory control]. Samara, Samara Research Centre RAS, 2007, 274 p.

[6] Ivanov N.M., Lysenko L.N. Ballistika i navigatsiya kosmicheskikh apparatov [Spacecraft ballistics and navigation]. 2nd ed. corr. and amend. Moscow, Dro-fa Publ., 2004, 544 p.

[7] Bordovitsyna T.V., Avdyushev V.A. Teoriya dvizheniya iskusstvenykh sput-nikov Zemli. Analiticheskie i chislennye metody [Theory of artificial satellites motion. Analytical and computational methods]. Tomsk, Tomsk University Publ., 2007, 178 p.

[8] Lysenko L.N., Betanov V.V., Zvyagin F.V. Teoreticheskie osnovy ballistiko-navigatsionnogo obespecheniya kosmicheskikh poletov [Theoretical fundamentals of ballistic and navigational equipment of space flights]. Moscow, BMSTU Publ., 2014, 518 p.

[9] Aksenov E.P., Chazov V.V. Model dvizheniya ISZ. Glavnye problemy. Osnov-nye algoritmy [Model of motion. Main problems. General algorithms]. Мoscow, 2007, 188 p.

[10] Chazov V.V. Prognoz orbitalnogo dvizheniya kosmicheskogo apparata. Chislennaya model [Prognosis of spacecraft orbital motion. Computational model]. Nauchno-tekhnicheskiy otchet [Scientific and technical report]. Available at: http://vadimchazov.narod.ru/text_pdf/comalg.pdf

[11] Malyshev V.V., Krasilschikov M.N., Bobronnikov V.T., Nesterenko O.P., Fe-dorov A.V. Sputnikovye sistemy monitoringa. Analiz, sintez i upravlenie [Satellite monitoring systems. Analysis, synthesis and control]. Moscow, MAI Publ., 2000, 568 p.

Vasiliev D.G., part-time post-graduate student of the Department of Dynamics and Flight Control of Rockets and Spacecraft, Bauman Moscow State Technical University. e-mail: vasiliev_dmitry83@list.ru

Betanov V.V., Dr. Sci. (Eng.), Professor, the director's advisor, Joint Stock company Russian Space Systems, Professor of the Department of Dynamics and Flight Control of rockets and spacecraft, Bauman Moscow State Technical University. Corresponding member of the Russian Academy of Rocket and Artillery Science. Author of more than 300 scientific papers and inventions in the field of ballistics, flight dynamics, rocket and spacecraft motion control, 5 monographs, 2 course books, 30 study guides methodological guidelines, co-author of the Strategic Missile Forces Military Encyclopedic Dictionary. e-mail: vlavab@mail.ru