CC BY
304
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2011 Математика и механика № 4(16)
УДК 523.24
Т.В. Бордовицына, А.Г. Александрова, И.Н. Чувашов ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ ОКОЛОЗЕМНЫХ КОСМИЧЕСКИХ ОБЪЕКТОВ ИСКУССТВЕННОГО ПРОИСХОЖДЕНИЯ С ИСПОЛЬЗОВАНИЕМ ПАРАЛЛЕЛЬНЫХ ВЫЧИСЛЕНИЙ1
Дается краткий обзор разработанного авторами обширного математического и программного обеспечения для исследования динамики больших совокупностей околоземных космических объектов искусственного происхождения. Представляемое программно-математическое обеспечение позволяет решать следующие задачи: исследовать одновременно орбитальную эволюцию большого числа искусственных спутников Земли (ИСЗ) и объектов космического мусора, включая анализ хаотичности движения; улучшать орбиты объектов по данным наблюдений; моделировать процесс образования и распределения космического мусора путем взрывов и столкновений; выявлять тесные сближения космических объектов и прогнозировать вероятность их столкновения
Ключевые слова: численные методы, искусственные спутники Земли, космический мусор, долговременная орбитальная эволюция, улучшение орбит, динамическая хаотичность, вероятность столкновения
Предполагается [1,2], что в результате деятельности человека в космосе на сегодня в околоземном пространстве находится около 14 000 объектов размером от 5 - 10 см и более, и только 4 % из них - работающие космические аппараты (КА). В совокупность неуправляемых объектов входят геодезические ИСЗ и космический мусор, состоящий из отработавших КА и верхних ступеней ракет-носителей, а также различных элементов конструкций КА, которые образуются вследствие разрушения КА под действием столкновений и взрывов.
По типу орбит все каталогизированные объекты делятся на следующие классы
или области:
LEO - low-Earth orbits, то есть низкоорбитальные объекты;
MEO - medium Earth orbits, объекты на орбитах между LEO и GEO;
GEO - geostationary orbits, объекты на геостационарных орбитах;
GTO - GEO transfer orbits, объекты на орбитах перехода в область GEO;
HEO - highly eccentric orbits, объекты с большими эксцентриситетами орбит. Последние два класса в значительной степени совпадают.
Прогнозирование движения и исследование орбитальной эволюции совокупности околоземных объектов требует создания разнообразного программноматематического обеспечения, ориентированного на решение следующих задач:
- высокоточное численное моделирование движения ИСЗ и представление лазерных, оптических и радиотехнических наблюдений;
- решение обратных задач динамики ИСЗ, то есть определение параметров движения и модели сил по данным измерений;
1 В статье обобщены результаты, полученные в рамках проектов № 2.1.1/2629, № 2.1.1/12782 по АВЦП «Развитие потенциала высшей школы» и госконтрактов № П1247 от 27 августа 2009 г., № П882 от 26 мая 2010 г. по ФЦП «Научные и научно-педагогические кадры инновационной России».
- моделирование процессов образования и распределения космического мусора;
- исследование динамики и долговременной орбитальной эволюции больших совокупностей околоземных объектов.
Для решения этих задач авторами настоящей работы создано с использованием параллельных вычислений обширное программно-математическое обеспечение, описанию которого посвящена данная статья.
«Численная модель движения систем ИСЗ» представляет собой методику и программу для высокоточного численного моделирования движения больших систем околоземных объектов с использованием параллельных вычислений. Данный программный комплекс имеет следующую структуру.
Уравнения движения объекта, рассматриваемого как материальная частица бесконечно малой массы, в поле тяготения центрального тела с массой М под действием сил, определенных потенциальными функциями V и Я, а также совокупности сил Р, не имеющих потенциала, представляются в виде
и интегрируются численно с помощью интегратор Гаусса - Эверхарта, модифицированного В.А. Авдюшевым [3]. Здесь У(Спт, £п,т) - потенциал притяжения Земли, Я = Ям + Я3, где Ям и Я3 - возмущающие функции, обусловленные соответственно притяжением Луны и Солнца; Р - возмущающее ускорение, создаваемое силами, не имеющими потенциала, такими, как световое давление и сопротивление атмосферы. О = О(/,И) - матрица перехода из вращающейся системы координат в инерциальную систему, где И - звездное время. Начальные параметры движения х0, х0, параметры модели сил Спт , $пт и И, а в случае необходимости и
сами силы Р, определяются из наблюдений в результате решения обратной задачи динамики ИСЗ.
Потенциал гравитационного поля Земли представлен в виде разложения по шаровым функциям в системе координат, жестко связанной с Землей. Шаровые функции и их частные производные вычисляются по рекуррентному алгоритму Каннингема [4]. В соответствии с рекомендациями Международной службы вращения Земли (ГЕКБ) все параметры разложения потенциала Земли берутся из модели геопотенциала БвМ2008, в которой определены коэффициенты гармоник до 2190 порядка и 2159 степени. Влияние приливных деформаций, происходящих в теле Земли под действием притяжения от Луны и Солнца, вводится в виде добавок [4] в мгновенные значения коэффициентов разложения гравитационного поля Земли. Учитываются: твердый прилив, модель Лява и отклонение модели Лява от модели Вара, полюсные и океанические приливы. При вычислении возмущений от Луны, Солнца и больших планет используются фонды координат больших планет: ББ405 - для высокоточных вычислений; и ББ406 - для исследования долговременной орбитальной эволюции околоземных космических объектов. При учете возмущений от светового давления вводится функция тени. Аналитические условия вхождения в тень, имеющую коническую форму, вычисляются через прямоугольные координаты спутника и Солнца [14].
Описание «Численная модель движения систем ИСЗ»
йґ йґ дх дх
(1)
с начальными условиями
Хо = х(ґ0), х о = х (ґо)
(2)
Кластер «Скиф СуЪепа» по структуре доступа к оперативной памяти относится к виду кластеров с распределенной оперативной памятью и позволяет задействовать в процессе обработки данных значительные ресурсы как оперативной памяти узла (до 8 Гб), так и процессорной памяти. Основной применяемый нами принцип распределения вычислений по ядрам кластера - это распределение по объектам. Разработанный нами программный комплекс позволяет отслеживать одновременно эволюцию орбит более 1000 объектов. Применение методов распараллеливания позволяет при запуске программного комплекса одновременно задействовать до 300 процессоров кластера и оптимально распределить между ними объекты. При таком распараллеливании быстродействие программного комплекса увеличивается в десятки раз по сравнению с одновременным интегрированием орбит 1000 объектов на одном процессоре. Второй важной особенностью кластера является возможность варьировать разрядную сетку от 32 до 128 бит. Это позволяет управлять точностью численной модели и ее быстродействием. При обращении к программе все параметры модели задаются во входном файле.
Использование 128-битной разрядной сетки позволяет существенно повысить точность при решении задач динамики ИСЗ по данным высокоточных, например лазерных, наблюдений, за счет учета слабых возмущений, использования метода численного интегрирования более высокого порядка и перевода ошибки округления в незначащие разряды. Распараллеливание вычислительного процесса оказывается эффективным при одновременном моделировании движения больших совокупностей околоземных объектов. Приведем примеры.
Нами были проведены численные эксперименты [4,5] по оценке порядка гармоник геопотенциала, которые можно учесть при работе на 64- и 128-битной сетках. Основные параметры орбит околоземных объектов, использованных в эксперименте, приведены в табл. 1.
Таблица 1
Параметры движения
Н, км а, км е и град Т, с
300 6678.13629 0.001 62.8 5431
500 6878.13629 0.001 62.8 5676
800 7178.13629 0.001 62.8 6052
1200 7578.13629 0.001 62.8 6565
1600 7978.13629 0.001 62.8 7091
2200 8578.13629 0.001 62.8 7906
3000 9378.13629 0.001 62.8 9038
4100 10738.13629 0.001 62.8 10674
5900(Лагеос) 12300.0 0.004 109.8 13680
20200(Эталон) 26600.0 0.01 63.4 43200
В качестве показателя влияния гармоник высокого порядка был принят вклад каждой следующей группы гармоник одного порядка относительно предыдущих. Были получены оценки влияний гармоник от 20 до 360 порядков. Оценки влияния гармоник сопоставлялись с величиной ошибки интегрирования. В качестве порядка гармоник, которые следует учитывать, выбирался порядок тех, влияние которых оказывалось больше ошибок интегрирования. В суммарном виде эти результаты приведены на рис. 1. На рис. 1, а даны оценки, полученные на 64-битной разрядной сетке, а на рис. 1, б - оценки, сделанные на 128-битной сетке. Эти
оценки показывают, во-первых, что использование большой разрядной сетки позволяет учитывать в процессе интегрирования значительно большее число гармоник геопотенциала и, во-вторых, что «Численная модель движения систем ИСЗ» позволяет подбирать оптимальный набор гармоник геопотенциала, обеспечивающий необходимую точность в процессе моделирования.
300
0 2000 4000 6000
Высота перигея Н, км
Высота перигея Н, км
Рис. 1. Зависимость числа значимых гармоник от высоты перигея спутника на различных разрядных сетках
Анализ влияния изменения коэффициентов Спт , Бпт разложения геопотенциала под действием приливных деформаций на точность прогнозирования движения также показывает [5], что в ряде случаев только использование 128-битной разрядной сетки позволяет отслеживать все особенности этого влияния.
Оценки накапливания методических ошибок численного моделирования на интервале 100 лет, полученные по данным прямого и обратного интегрирования уравнений движения спутников Эталон 1 и 2 с использованием 64-битной и 128битной разрядных сеток приведены на рис. 2.
Рис. 2. Оценки точности прогнозирования векторов положения объектов. На верхнем графике - для Эталона 1, на нижнем графике - для Эталона 2
Численное моделирование процесса образования орбитальной эволюции и распределения фрагментов космического мусора в околоземном пространстве
Математическая модель распада строится в предположении, что координаты каждого фрагмента совпадают с координатами родительского тела, а компоненты скорости фрагмента определяются по формулам
*10 = V +Av cos т,
x20 = v2 + Ду sin т cos ф,
*30 = v3 +Ду sin т sin ф,
где v1, v2, v3 - компоненты скорости родительского тела в геоцентрической системе координат. Параметры Ду, т, ф задают величину и направление вектора скорости фрагмента относительно родительского тела и рассматриваются как случайные величины, которые определяются с помощью метода обратных функций по заданным функциям плотности распределения. Дальнейшее движение фрагментов моделируется с помощью программы «Численная модель движения ИСЗ». Приведем сравнение с наблюдениями двух моделей распада: изотропного взрыва (рис. 3) и столкновения (рис. 4). Данные взяты из работы авторов [6].
Период, мин Период, мин
Рис. 3. Распределение фрагментов распада КА Космос 1275 через неделю после взрыва по данным каталога НАСА (слева) и по данным моделирования (справа)
Следует отметить, что высокая точность представления наблюдаемой картины разброса фрагментов достигается не в момент распада, а в момент наблюдения, который отстоит от момента распада на значительный интервал времени, что говорит о высоком качестве и математической модели распада и численной модели движения систем ИСЗ.
Период, мин Период, мин
Рис. 4. Распределение фрагментов распада USA 19 через день после столкновения по данным каталога NASA (слева) и по данным моделирования (справа)
MEGNO-анализ динамики околоземных объектов
Алгоритм выявления хаотичности в орбитальной эволюции объектов основан на MEGNO-анализе [7]. Параметр MEGNO (Mean Exponential Growth of Nearby Orbit) представляет собой оценку среднего экспоненциального расхождения двух близких орбит, которая характеризует уровень хаотичности движения. Причем известно, что для квазипериодических (регулярных) орбит осредненный параметр MEGNO Уф (t) осциллирует около 2, для круговых орбит Уф (t) всегда равно 2, а
для устойчивых орбит типа гармонического осциллятора Уф (t) = 0 . MEGNO-анализ долговременной орбитальной динамики объекта реализован нами с использованием программного комплекса «Численная модель движения систем ИСЗ». Комплекс дополнен интегрированием восьми уравнений для вычисления касательного вектора б = б (х,Х) и величин y(t) и w(t) [8], связанных с неосредненным
Y(t) и осредненным У (t) параметрами MEGNO соотношениями
У(t) = 2y(t)/1, У (t) = w(t)/1. (3)
Подробное описание алгоритма и результатов его тестирования дано в [9].
Алгоритм определения минимальных расстояний между объектами
Минимальное расстояние между искусственными спутниками Земли вычисляются путем нахождения локального минимума квадрата функции расстояния
[10], аппроксимируемого интерполяционным многочленом Лагранжа третьей степени относительно времени. Многочлен Лагранжа строится по узлам таблицы интегрирования уравнений движения
L3(t) =Е rf П t-^J~, j = 0Д (4)
i=0 j*i ti lj
где ri - расстояние между двумя спутниками в момент времени t.
Для нахождения минимума функции расстояния определим производную многочлена Лагранжа (4) и приравняем ее к нулю.
ёЬ3 (ґ)
X П (ґ - ґк)
’ = 2 г? 1 *гк^1------= 0, 1 = 0,3, к = 0,3. (5)
Ж г=0 г П ('г- - 0 )
1 *г
Решив квадратное уравнение (5), получим момент сближения (шт, и подставляя его в (4), определим минимальное расстояние до спутника гт1п = Ь3 ('т1п) . Вычисление минимального расстояния осуществляется на каждом шаге интегрирования.
Комплекс программ для решения обратных задач динамики ИСЗ
Обратная задача динамики ИСЗ и алгоритм ее решения могут быть сформулированы следующим образом. Пусть р = {р^)} есть вектор измеренных величин,
]=1,2,...Ы, а q={qi}- вектор определяемых параметров (г=1,2,...т), связь между которыми через решение уравнений (1) задается нелинейным соотношением
р = F(q),
причем 1 »г, то есть мы имеем избыточную нелинейную систему уравнений для определения неизвестных параметров. Решение этой системы осуществляется при условии существования минимума функции
N 2
Ф^) = 2[АРг (q)] = т1п, (6)
г=1
где Др/ = ро - рС - невязки, представляющие собой разность наблюденных и вычисленных значений измеряемых величин. Задачу минимизации (6) принято называть задачей наименьших квадратов (НК). Решение нелинейной задачи НК производится методом Гаусса - Ньютона путем нахождения дифференциальных поправок Дq в оценки q определяемых элементов q. Алгоритм вычисления дифференциальных поправок Дq можно записать в следующем виде:
дqп+1 = дqп + ((qп )т и йп) )-1 и йп )т Др(qп), (7)
(I п+1 = (|п +дq п+\
где Др(() - вектор невязок; Щ^) - матрица частных производных от измеряемых параметров по оцениваемым. Причем, как было сказано выше, N - число измерений, а т - число определяемых параметров.
Итерационный процесс считается завершенным при выполнении следующих условий:
\чГ1 -?Г| < в, (8)
где в - заданная точность вычислений.
Среднеквадратическая ошибка единицы веса определяется формулой
^0 =[ф(?)/^ - т)]1/2, а среднеквадратические ошибки определяемых параметров задаются как стг- =<з0л/^п, г = 1,2,...т, где - диагональный элемент обратной матрицы
п )т п)) 1. Обратная матрица вычисляется методами сингулярного анализа
[11].
Программный комплекс для решения обратных задач динамики ИСЗ тестировался на лазерных наблюдениях геодинамических ИСЗ Лагеос и Эталон. Обработка наблюдений для ИСЗ Лагеос проводилась на 7-суточной орбитальной дуге, а для ИСЗ Эталон - на 30-суточной орбитальной дуге, что в том и другом случаях соответствует 60 оборотам спутника. Данные, приведенные в работе [12] одного из авторов настоящей статьи, показывают, что при тщательной отбраковке наблюдений система параметров движения для спутника Лагеос имеет среднеквадратичную ошибку единицы веса равную нескольким сантиметрам, а для спутника Эталон 2 равную двум метрам, что хорошо согласуется с результатами других авторов.
Алгоритм исследования долговременной эволюции доверительных областей движения. Способ вероятностной оценки возможных столкновений объектов
Алгоритм исследования долговременной эволюции доверительных областей движения в виде граничной поверхности. Следуя [13], будем определять доверительную область движения объекта его граничной поверхностью по следующему алгоритму.
Как известно, начальная вероятностная область решения q для каждого объекта может быть определена относительно МНК-оценки q с использованием ковариационной матрицы Б по формуле
qг = Апг + q (г = 1,..., N), (9)
где п" - 6-мерный вектор случайных чисел, распределенных по нормальному закону, А - треугольная матрица, такая, что АТА = Б, а N - число рассматриваемых решений. Точки qI дают вероятностное распределение возможных значений q в фазовом пространстве определяемых параметров.
По количеству наблюдений, числу определяемых параметров т и заданной вероятности Р (в нашем случае мы полагали Р = 0,999) определяется верхняя квантиль функции Фишера Г . По известным значениям среднеквадратической ошиб-
2 7--*
ки единицы веса а0 и г легко наити правую часть уравнения:
(яг - q )То(^г - q ) = тст^* = у*, (Ю)
где о=кТи - так называемая нормальная матрица, а К - матрица изохронных производных.
Значение у определяет граничную поверхность доверительной области. Для
—I _ ^ *
всех точек q эллипсоидальной поверхности уг- = у .
По формуле (9) задается множество случайных точек q,, для каждой из которых вычисляется соответствующее значение у{:
Уг = ^ - q )Т о ^ - q).
С помощью коэффициентов и = у* / у,-, растягивающих вектор ^1 - q) до гра-
»_» * »_» —г
ничной поверхности у , определяется искомый вектор q :
(^г - q) = 1г (qг - q)
C помощью данного алгоритма методом Монте-Карло формируется доверительная область в виде множества случайных многомерных точек, заполняющих граничную поверхность доверительного эллипсоида. В отличие от классического алгоритма, основанного на генерировании случайных точек, заполняющих всю область, используемый нами алгоритм, как показано в работе [14], является более экономичным по числу моделируемых точек, и при одинаковом числе точек точнее определяет доверительные области.
Способ вероятностной оценки возможных столкновений объектов. Метод исследования долговременной эволюции доверительных областей движения в виде граничной поверхности, как уже отмечалось выше, является более экономичным по числу моделируемых точек, чем классический, однако он не применим для вероятностной оценки возможности столкновения исследуемых объектов. В случае пересечения доверительных областей движения КА, построенных в виде граничных поверхностей, предлагается строить для этих объектов повторно доверительные области классическим способом, с заполнением начальной вероятностной области методом Монте-Карло по формуле (9) и отслеживать их эволюцию до момента пересечения. А оценку вероятности столкновения объектов космического мусора определять как отношение числа траекторий, попавших в область пересечения, к их общему числу.
Исследование долговременной орбитальной эволюции всей совокупности объектов каталога ESA «Classification of geosynchronous objects»
MEGNO-анализ особенностей динамики объектов каталога. Как известно, все функционирующие КА зоны GEO располагаются над экватором Земли в различных его точках в зависимости от зоны обслуживания. После утраты управления отработавшие объекты переходят в состояние свободного полета, который определяется законами небесной механики. Особенности воздействия тессерального резонанса таковы, что дальнейшая динамика объекта будет определяться долготой X подспутноковой точки, с которой объект начал свое свободное движение.
На рис. 5 показана эволюция всей совокупности объектов каталога ESA на десятилетнем интервале времени. На графике хорошо видно, что достаточно большое число объектов находятся вблизи сепаратрис.
42000 -
41900
г
180
270
X,град
0
Рис. 5. Распределение КА каталога ESA в результате десятилетней эволюции
Таким образом, во всех случаях фракция неустойчивых объектов присутствует и может быть даже достаточно большой.
На рис. 6 показано значение усредненного параметра MEGNO для всех неуправляемых на 1 января 2009 г. КА геостационарной зоны, приведенных в каталоге ESA «Classification of geosynchronous objects», на интервалах времени 30, 100 и 200 лет. На графиках хорошо видно, что со временем все большее количество КА оказывается на неустойчивых орбитах. И если на интервале времени 30 лет неустойчивых объектов относительно немного, то через 100 и 200 лет их число возрастает в несколько раз.
a, км 42200
42150
42100 -
42050 1
-15 45 ' 105 165 225 285 X, град
42200
42150
42100 -
42050
-15 45
105
I I Г 165 225 285 X, град
a, км
a, км 42200
42150
42100
42050
-15 45 105 165 225 285 X, град
Рис. 6. Усредненный параметр MEGNO для объектов каталога ESA «Classification of geosynchronous objects» на интервале времени: а -30 лет; б - 100 лет; в - 200 лет
Исследование возможных сближений объектов каталога ESA «Classification of geosynchronous objects»
В данном разделе статьи приведены результаты исследования долговременной орбитальной эволюции всей совокупности объектов геостационарной зоны, приведенных в 11-м издании каталога ESA «Classification of geosynchronous objects». В этом каталоге орбитальные данные объектов на конкретный момент времени представлены без указания их точности и интервала времени, охваченного наблюдениями, по которым эти параметры были получены. Поэтому для построения доверительных областей нам пришлось моделировать наблюдения выбранных объектов и процесс улучшения их орбит.
Начальные значения параметров движения КА, которые имеются в каталоге, мы рассматривали как «точные». Затем вычислялись угловые положения КА на заданном интервале времени с определенным шагом. Мы выбрали интервал равный одному обороту КА с шагом 0,01 оборота. Далее путем внесения в угловые положения КА случайной ошибки из заданного интервала точностей мы формировали наблюдения. Решая задачу наименьших квадратов (НК) методом дифференциальных поправок, используя для этого сформированную выборку измерений, находили НК-оценки начальных параметров и их матрицу ковариации. В качестве начального приближения в итерационном методе использовали значения параметров движения КА, взятые из каталога. Затем по НК-оценке начальных параметров и матрице ковариации строили доверительную область движения КА, накрывающую с заданной вероятностью (P = 0,999) положение КА на начальный момент времени. Доверительная область (6-мерный эллипсоид) определялась методом Монте-Карло в виде множества 6-мерных точек. Эти точки, включая НК-оценку и значения параметров из каталога, рассматривались в дальнейшем в качестве начальных точек ансамбля траекторий, динамическая эволюция которых позволяет построить картину вероятностной эволюции движения КА.
Как уже отмечалось выше, области являются носителями информации о реальности наших знаний об орбитальном движении. Мы провели оценку влияния точности наблюдений на размеры доверительной области движения. Полученные оценки показывают зависимость размеров доверительной области от величины ошибки [14].
Для оценки сближения объектов геостационарной зоны между собой был выполнен прогноз движения на интервале времени 10 лет с шагом 1 с. Минимальные расстояния между объектами вычислялись по вышеизложенному алгоритму. Все выявленные сближения на расстояние менее 100 км, их оказалось 514556, приведены на рис. 7. Необходимо отметить, что произошло 30125 сближений на расстояние менее 20 км и 12274 - на расстояние менее 10 км, что считается крайне опасным. Данные о наиболее тесных сближениях (менее 100 м) приведены в табл. 2. Причем следует заметить, что данные объекты неоднократно имели тесные сближения.
В случае сближений менее 100 м нами были построены доверительные области для сближающихся объектов и прослежена их динамическая эволюция до момента наиболее тесного сближения. Для примера рассмотрим вариант сближения КА Cosmos 2133 и Raduga 10. На рис. 8 показаны начальные доверительные области для КА Cosmos 2133 и Raduga 10 (ст = 0,5") и на рис. 9 - в момент наиболее тесного сближения. Для КА Cosmos 2133 и Raduga 10 на момент сближения идет пересечение доверительных областей, что говорит о возможности столкновения.
AR, км —г
НО -
514556
40 -
3UJ25
12574
—1 | I | I | I | I |
2-15'IООП ^35000 2-156000 2-157000 ^5£000 10, сут
Рис. 7. Оценки сближения объектов геостационарной зоны между собой
Т аблица 2
Объекты, испытавшие тесные сближение
КА ДЯ, км Дата сближения
Gstar 3 OPS 9431 (DSCS II F-1) 0,15 13.05.2009
ACTS ASC 1 0,61 13.06.2009
LES 9 (RTGPP) 69101A Skynet 1A 0,16 14.05.2010
Raduga 1-3 Raduga 9 0,B3 2B.0B.2010
Cosmos 2133 Raduga 30 0,B3 20.01.2012
Cosmos 2133 Raduga 10 0,13 09.02.2012
N-Star 2 Eutelsat II F-3 0,64 21.11.2013
Raduga 32 Raduga 14 0,62 25.12.2013
ACTS Gstar 3 0,57 03.04.2014
ATS 1 Cosmos 1B97 0,4B 31.01.2015
Comstar 4 Cosmos 1546 0,9B 04.05.2015
Zhongxing б (B) Ekran 2 fragm. debris 0,26 24.09.2015
Raduga 12 Raduga 14 0,2B 10.11.2015
Raduga12 Raduga 14 0,76 22.11.2016
Cosmos 1961 Elektro 1 0,75 14.02.2017
Proton-K (B5070F) Proton-K (B0049F) 0,73 20.04.201B
IUS stage 2 Raduga 24 0,71 27.10.201B
Gorizont 27 Gorizont б 0,07 13.12.201B
На рис. 9 сплошными линиями показаны траектории номинальных орбит КА Cosmos 2133 (кривая 1) и Raduga 10 (кривая 2), точками - соответствующие им доверительные области движения на момент пересечения областей.
Вероятность риска столкновения рассчитывалась в виде отношения числа столкновительных траекторий к общему числу траекторий. В данном случае она оказалась низкой P = B-10-6.
5806.5 ■
5806.0 | 5805.5. * 5805.0
5804.5 ■
5804.0 -
IX
vm-'iffl
, ' ;#>; ‘ s'
■4'W’
3351.5 -5351.4 5351.4-S3JU-5351.35351.2-
5351.2' '.V -Vi'^05-2
. >.-1 J-H
W'
ЙЖЙУ ■■
rfiW.»,<■ ■■■■
' ...
6304
i.V
3-bq
•Cg?*
; ; 1" • • ' • ' .'Б.Ч
--.A' ;.4ii55.4fi
• •. • ' it355.92
-f-i"- --.':%-ii'355.sa ^ "■■'Й55.Я4 f,v
■4
Рис. 8. Начальные области возможных движений для КА Cosmos 2133 (слева) и Raduga 10 (справа)
■гькм *v
Рис. 9. Вероятностные области для КА Cosmos 2133 (кр. 1) и Raduga 10 (кр. 2) на момент сближения в разных масштабах
Заключение
Таким образом, в настоящей работе дано описание разработанного авторами алгоритмического и программного обеспечения, которое позволяет решать широкий круг задач динамики околоземных объектов. Все программные комплексы реализованы в среде параллельных вычислений и являются высокоэффективными как по точности, так и по быстродействию.
ЛИТЕРАТУРА
1. Рыхлова Л.В. Засоренность околоземного пространства объектами техногенного происхождения // Околоземная астрономия - 2003: тр. конф. Т.2. Терскол, сентябрь 2003 г. Институт астрономии РАН. СПб.: ВВМ, 2003. С. 11-19.
2. KlinkradH. Space debris. Springer, 2006. 430 p.
3. Авдюшев В.А. Интегратор Гаусса - Эверхарта // Вычисл. технологии. 2010. Т. 15. № 4. С. 31-47.
4. Бордовицына Т.В, Авдюшев В.А. Теория движения ИСЗ. Аналитические и численные методы. Томск: Изд-во Том. ун-та, 2007. 105 с.
5. Чувашов И.Н. Прогнозирование движения ИСЗ с использованием параллельных вычислений. Учет слабых возмущений // Изв. вузов. Физика. 2010. T. 53. № 8/2. С. 22-29.
6. Бордовицына Т.В., Александрова А.Г. Численное моделирование процесса образования орбитальной эволюции и распределения фрагментов космического мусора в околоземном пространстве // Астрон. вестн. 2010. Т. 44. С. 259-272.
7. Cincotta P.M.. Giordano C.M, Simob C. Phase space structure of multi-dimensional systems by meansof the mean exponential growth factor of nearby orbits // Physica D. 2003. V. 182. P. 151-178.
8. Valk S., Delsate N., Lemaitre A., Carletti T. Global dynamics of high area-to-mass ratios GEO space debris by means of the MEGNO indicator // Adv. Space Res. 2009. V. 43. P. 1509-1526.
9. Бордовицына Т.В., Александрова А.Г., Чувашов И.Н. Комплекс алгоритмов и программ для исследования хаотичности в динамике искусственных спутников Земли // Изв. вузов. Физика. 2010. T. 53. № 8/2. С. 14-21.
10. Шефер В.А. Регуляризирующие и стабилизирующие преобразования в задаче исследования движения особых малых планет и комет: автореф. дис. ... к.ф.-м.н. Казань, 1986. 13 с.
11. Форсайт Дж., Малькольм М., Моулер К. Машинные методы математических вычислений. М.: Мир, 1980. 279 с.
12. Чувашов И.Н. Программно-математическое обеспечение для решения обратных задач динамики ИСЗ с использованием параллельных вычислений // Изв. вузов. Физика. 2011. № 6/2. С. 5-12.
13. Сюсина О.М., Тамаров В.А., Черницов А.М. Новые алгоритмы построения методом Монте-Карло начальных доверительных областей движения малых тел // Изв. вузов. Физика. 2009. T. 52. № 10/2. С. 48-55.
14. Александрова А.Г., Бордовицына Т.В., Чувашов И.Н. Об исследовании долговременной эволюции доверительных областей движения объектов геостационарной зоны // Изв. вузов. Физика. 2009. T. 52. № 10/2. С. 20-25.
Статья поступила 21.11.2011 г.
Bordovitsyna T.V., Aleksandrova A.G., Chuvashov I.N. NUMERICAL SIMULATION OF NEAR EARTH ARTIFICIAL SPACE OBJECT DYNAMICS USING PARALLEL COMPUTATION. A short survey of algorithms and software developed by the authors for studying dynamics of large groups of near Earth artificial space object is given. The software permits one to solve the following problems: simultaneous investigation of orbital evolution of a large number of artificial satellites and objects of space debris, including the analysis of the chaotic state, improving orbits, simulating the process of space debris formation and distribution by means of explosions and collisions, discovering approaches of space objects, and forecasting the probability of their collisions.
Keywords: numerical methods, Earth artificial satellites, space debris, long-term orbital evolution, improvement of orbits, dynamical randomness, probability of collision.
BORDOVITSYNA Тatiana Valentinovna (Tomsk State University)
E-mail: tvbord@sibmail.com
ALEXANDROVA Anna Gennad’evna (Tomsk State University)
E-mail: chuvashov@sibmail.com
CHUVASHOV Ivan Nikolaevich (Tomsk State University)
E-mail: aleksann@sibmail.com
