Обсудить в форуме Комментариев 12Редактировать в вики
Определение положения ИСЗ по орбитальным данным на заданное время по модели SGP4. Автоматизированное получение орбитальных данных с помощью API сервиса space-track.org. Пример реализации на языке Python.
Задачу определения положения того или иного искусственного спутника Земли в заданный момент времени (в прошлом или недалёком будущем) приходится решать для самых разнообразных целей, в том числе связанных с дистанционным зондированием Земли из космоса. Часть данных (например, многие продукты MODIS) распространяется без строгой географической привязки, а лишь с указанием времени непосредственного наблюдения территории для каждой сцены, — и для автоматизации поиска и загрузки таких данных требуется вычислять время пролёта спутника над исследуемыми объектами. Часто возникает и потребность определить время зондирования заданной территории в будущем - чаще всего для проведения подспутниковых наблюдений (в целях верификации, атмосферной коррекции и пр.).
В статье описывается подход к моделированию проекций орбит ИСЗ на поверхность Земли с использованием доступных средств: библиотек языка Python и API сервиса space-track.org.
Наиболее распространенной моделью для определения положения спутников на орбите является SGP (Simplified General Perturbations), различные модификации которой используются в оперативной работе по всему миру начиная с 70-х годов. Главная задача модели - вычислить скорость и геоцентрические координаты ИСЗ (X, Y, Z) на заданный момент времени, которые нетрудно пересчитать на поверхность эллипсоида, получив географические координаты проекции положения ИСЗ (широта, долгота). Сама модель достаточно сложна, хотя и сводится к линейным расчётам и удобна для алгоритмизации. Её описание и оригинальный FORTRAN-код можно найти в соответствующих документах [1,2].
В качестве входных параметров SGP использует данные телеметрии спутников в формате TLE (two-line element sets): это две линии по 69 символов, описывающие основные метаданные спутника и параметры телеметрии [3]. Содержание первой линии:
Номер | Положение | Содержание | Пример |
---|---|---|---|
1 | 01-01 | Номер строки | 1 |
2 | 03-07 | Номер спутника в базе данных NORAD | 25994 |
3 | 08-08 | Классификация (U=Unclassified — не секретный) | U |
4 | 10-11 | Международное обозначение (последние две цифры года запуска) | 99 |
5 | 12-14 | Международное обозначение (номер запуска в этом году) | 068 |
6 | 15-17 | Международное обозначение (часть запуска) | A |
7 | 19-20 | Год эпохи (последние две цифры) | 16 |
8 | 21-32 | Время эпохи (целая часть — номер дня в году, дробная — часть дня) | 052.07623983 |
9 | 34-43 | Первая производная от среднего движения (ускорение), деленная на два [виток/день^2] | .00001336 |
10 | 45-52 | Вторая производная от среднего движения, деленная на шесть (подразумевается, что число начинается с десятичного разделителя) [виток/день^3] | 00000-0 |
11 | 54-61 | Коэффициент торможения B* (подразумевается, что число начинается с десятичного разделителя) | 30635-3 |
12 | 63-63 | Изначально — типы эфемерид, сейчас — всегда число 0 | 0 |
13 | 65-68 | Номер (версия) элемента | 999 |
14 | 69-69 | Контрольная сумма по модулю 10 | 6 |
Собранный пример: 1 25994U 99068A 16052.07623983 .00001336 00000-0 30635-3 0 9996
Содержание второй линии:
Номер | Положение | Содержание | Пример |
---|---|---|---|
1 | 01-01 | Номер строки | 1 |
2 | 03-07 | Номер спутника в базе данных NORAD | 25994 |
3 | 09-16 | Наклонение в градусах | 98.1986 |
4 | 18-25 | Долгота восходящего узла в градусах | 128.0087 |
5 | 27-33 | Эксцентриситет (подразумевается, что число начинается с десятичного разделителя) | 0001485 |
6 | 35-42 | Аргумент перицентра в градусах | 109.3968 |
7 | 44-51 | Средняя аномалия в градусах | 250.7393 |
8 | 53-63 | Частота обращения (оборотов в день) (среднее движение) [виток/день] | 14.57136668 |
9 | 64-68 | Номер витка на момент эпохи | 86046 |
10 | 69-69 | Контрольная сумма по модулю 10 | 2 |
Собранный пример: 2 25994 98.1986 128.0087 0001485 109.3968 250.7393 14.57136668860462
Важно понимать, что такие эфемериды описывают мгновенное состояние ИСЗ, и, хотя описывают его поведение с довольно высокой точностью, при увеличении дальности прогноза (относительно данной эпохи) будут давать всё большую и большую ошибку.
Данные TLE сегодня публикуются многими поставщиками (например, последние данные TLE по ряду спутников ДЗЗ на сайте ScanEx), но нам нужно получать не только свежие данные, но и архивные, для моделирования положений спутников в прошлом.
Одним из лучших в сети ресурсов представляется портал space-track.org, предоставляющий доступ к обширной информации о спутниках различного назначения. Очень важно, что space-track имеет REST API, позволяющее получать нужные данные максимально удобно. Требуется авторизация (и для доступа к интерфейсу, и для программного обращения к API), регистрация при этом бесплатная и открытая. Забегая вперёд, скажем в пользу space-track ещё то, что для работы с его API существует открытая python-библиотека.
Непосредственно в интерфейсе сайта можно запрашивать данные TLE (в разделе Retrieve TLE Data by Satellite Catalog Number), заполнив небольшую форму с указанием названия или идентификатора спутника, а также интересующего вас периода времени. Для примера запросим данные TLE для спутника AQUA на середину мая 2012 года:
Результат вы получаете мгновенно. Примечательно, что сразу же при выдаче ответа сервис выводит команду API, соответствующую вашему запросу - это позволяет очень быстро разобраться в том, как оно организовано и как с ним работать.
https://www.space-track.org/basicspacedata/query/class/tle/EPOCH/2012-05-11--2012-05-12/NORAD_CAT_ID/27424/orderby/TLE_LINE1 ASC/format/tle
Выполнив этот запрос тут же в адресной строке браузера (т.е. реализовав простой HTTP-запрос), можно увидеть, что при работе с API данные представляются в незамысловатом текстовом виде, в котором их, учитывая строгую структуру формата, несложно интерпретировать программно.
В целом API подробно документировано. Для нашей задачи вполне достаточно рассмотреть тот пример, который был получен для майских приключений спутника AQUA. Изменяемыми в этой строке запроса будут всего два параметра:
Идентификатор нужного вам ИСЗ можно найти там же, на space-track, в разделе SATCAT, в удобном интерактивном интерфейсе. Нас интересует первая колонка таблицы результатов поиска. Например, поищем идентификаторы спутников программы Landsat:
Landsat 8 соотвествует номеру 39084. Попробуем найти актуальные TLE для этого спутника, заодно посмотрев, как изменится структура запроса при использовании не диапазона дат, а опции "Latest", т.е. "последние данные". Запрос:
https://www.space-track.org/basicspacedata/query/class/tle_latest/ORDINAL/1/NORAD_CAT_ID/39084/orderby/TLE_LINE1%20ASC/format/tle
и ответ:
1 39084U 13008A 16354.79369944 .00000065 00000-0 24444-4 0 9999 2 39084 98.2045 61.9547 0001318 94.3108 265.8241 14.57119154204938
Как видно, порядок аргументов в запросе изменился.
Открытая программная реализация модели SGP4 доступна для C++ и Python (библиотека pyorbital). Для примера будем использовать именно Python и pyorbital (есть и другая реализация на Python'e: python-sgp4). Для получения данных от API space-track.org доступна специальная библиотека. Чтобы представить результат в формате геоданных применим библиотеку pyshp. Поскольку за нас уже почти всё сделали, код очень прост. Разберём его по разделам.
Все библиотеки доступны в основном репозитории Python и устанавливаются очень просто
pip install pyorbital pip install spacetrack pip install pyshp
# Импортируем библиотеки
# Штатная библиотека для работы со временем
from datetime import datetime, date
# Собственно клиент для space-track
# Набор операторов для управления запросами. Отсюда нам понадобится время
import spacetrack.operators as op
# Главный класс для работы с space-track
from spacetrack import SpaceTrackClient
# Имя пользователя и пароль сейчас опишем как константы
USERNAME = <YOUR SPACE-TRACK USERNAME>
PASSWORD = <YOUR SPACE-TRACK PASSWORD>
# Для примера реализуем всё в виде одной простой функции
# На вход она потребует идентификатор спутника, диапазон дат, имя пользователя и пароль. Опциональный флаг для последних данных tle
def get_spacetrack_tle (sat_id, start_date, end_date, username, password, latest=False):
# Реализуем экземпляр класса SpaceTrackClient, инициализируя его именем пользователя и паролем
st = SpaceTrackClient(identity=username, password=password)
# Выполнение запроса для диапазона дат:
if not latest:
# Определяем диапазон дат через оператор библиотеки
daterange = op.inclusive_range(start_date, end_date)
# Собственно выполняем запрос через st.tle
data = st.tle(norad_cat_id=sat_id, orderby='epoch desc', limit=1, format='tle', epoch = daterange)
# Выполнение запроса для актуального состояния
else:
# Выполняем запрос через st.tle_latest
data = st.tle_latest(norad_cat_id=sat_id, orderby='epoch desc', limit=1, format='tle')
# Если данные недоступны
if not data:
return 0, 0
# Иначе возвращаем две строки
tle_1 = data[0:69]
tle_2 = data[70:139]
return tle_1, tle_2
Представлена очень простая функция, использующая клиентскую библиотеку space-track для получения одного (первого из запроса) набора tle. Пример её использования:
# Запросим данные о положении Landsat 8 11 мая 2016 года
# Обратите внимание, что даты указываем в формате date(y,m,d)
tle_1, tle_2 = get_spacetrack_tle (39084, date(2016,5,11), date(2016,5,12), USERNAME, PASSWORD)
print tle_1, tle_2
>>> 1 39084U 13008A 16132.92196421 +.00000109 +00000-0 +34320-4 0 9999
>>> 2 39084 098.2260 203.0765 0001471 094.1169 266.0197 14.57124417160799
# А теперь данные об актуальном положении
tle_1, tle_2 = get_spacetrack_tle (39084, None, None, USERNAME, PASSWORD, True)
print tle_1, tle_2
>>> 1 39084U 13008A 16354.79369944 .00000065 00000-0 24444-4 0 9999
>>> 2 39084 98.2045 61.9547 0001318 94.3108 265.8241 14.57119154204938
# Импортируем библиотеки
# Штатная библиотека для работы со временем
from datetime import datetime, date
# Ключевой класс библиотеки pyorbital
from pyorbital.orbital import Orbital
# Ещё одна простая функция, для демонстрации принципа.
# На вход она потребует две строки tle и время utc в формате datetime.datetime
def get_lat_lon_sgp (tle_1, tle_2, utc_time):
# Инициализируем экземпляр класса Orbital двумя строками TLE
orb = Orbital("N", line1=tle_1, line2=tle_2)
# Вычисляем географические координаты функцией get_lonlatalt, её аргумент - время в UTC.
lon, lat, alt = orb.get_lonlatalt(utc_time)
return lon, lat
Пример использования:
# Используем данные TLE полученные вручную на space-track.org для спутника Terra
tle_1 = '1 25994U 99068A 16355.18348138 .00000089 00000-0 29698-4 0 9992'
tle_2 = '2 25994 98.2045 66.7824 0000703 69.9253 290.2059 14.57115924904601'
# Нас интересует текущий момент времени
utc_time = datetime.utcnow()
# Обращаемся к фукнции и выводим результат
lon, lat = get_lat_lon_sgp (tle_1, tle_2, utc_time)
print lon, lat
>>> 175.589796941 -13.6408377148
Теперь объединим получение данных space-track и расчёт положения спутника, добавив создание точечного шейп-файла. Зададимся целью написать функцию, которая бы создавала точечный шейп-файл с положениями спутника в течение указанных суток с заданным шагом по времени, в атрибуты каждой точки записывая широту, долготу и время пролёта.
# Импортируем библиотеки - для начала оговоренные ранее
from datetime import datetime, date, timedelta
import spacetrack.operators as op
from spacetrack import SpaceTrackClient
from pyorbital.orbital import Orbital
# И pyshp, которая понадобится для создания шейп-файла
import shapefile
# Имя пользователя и пароль
USERNAME = <YOUR SPACE-TRACK USERNAME>
PASSWORD = <YOUR SPACE-TRACK PASSWORD>
# Уже описанная ранее функция get_spacetrack_tle может использоваться без изменений
def get_spacetrack_tle (sat_id, start_date, end_date, username, password, latest=False):
st = SpaceTrackClient(identity=username, password=password)
if not latest:
daterange = op.inclusive_range(start_date, end_date)
data = st.tle(norad_cat_id=sat_id, orderby='epoch desc', limit=1, format='tle', epoch = daterange)
else:
data = st.tle_latest(norad_cat_id=sat_id, orderby='epoch desc', limit=1, format='tle')
if not data:
return 0, 0
tle_1 = data[0:69]
tle_2 = data[70:139]
return tle_1, tle_2
# А вот функция get_lat_lon_sgp нам уже не пригодится в своём виде
# ведь создавать экземпляр класса Orbital для каждого момента времени
# не очень-то хочется
# На вход будем требовать идентификатор спутника, день (в формате date (y,m,d))
# шаг в минутах для определения положения спутника, путь для результирующего файла
def create_orbital_track_shapefile_for_day (sat_id, track_day, step_minutes, output_shapefile):
# Для начала получаем TLE
# Если запрошенная дата наступит в будущем, то запрашиваем самые последний набор TLE
if track_day > date.today():
tle_1, tle_2 = get_spacetrack_tle (sat_id, None, None, USERNAME, PASSWORD, True)
# Иначе на конкретный период, формируя запрос для указанной даты и дня после неё
else:
tle_1, tle_2 = get_spacetrack_tle (sat_id, track_day, track_day + timedelta(days = 1), USERNAME, PASSWORD, False)
# Если не получилось добыть
if not tle_1 or not tle_2:
print 'Impossible to retrieve TLE'
return
# Создаём экземляр класса Orbital
orb = Orbital("N", line1=tle_1, line2=tle_2)
# Создаём экземпляр класса Writer для создания шейп-файла, указываем тип геометрии
track_shape = shapefile.Writer(shapefile.POINT)
# Добавляем поля - идентификатор, время, широту и долготу
# N - целочисленный тип, C - строка, F - вещественное число
# Для времени придётся использовать строку, т.к. нет поддержки формата "дата и время"
track_shape.field('ID','N',40)
track_shape.field('TIME','C',40)
track_shape.field('LAT','F',40)
track_shape.field('LON','F',40)
# Объявляем счётчики, i для идентификаторов, minutes для времени
i = 0
minutes = 0
# Простой способ пройти сутки - с заданным в минутах шагом дойти до 1440 минут.
# Именно столько их в сутках!
while minutes < 1440:
# Расчитаем час, минуту, секунду (для текущего шага)
utc_hour = int(minutes // 60)
utc_minutes = int((minutes - (utc_hour*60)) // 1)
utc_seconds = int(round((minutes - (utc_hour*60) - utc_minutes)*60))
# Сформируем строку для атрибута
utc_string = str(utc_hour) + '-' + str(utc_minutes) + '-' + str(utc_seconds)
# И переменную с временем текущего шага в формате datetime
utc_time = datetime(track_day.year,track_day.month,track_day.day,utc_hour,utc_minutes,utc_seconds)
# Считаем положение спутника
lon, lat, alt = orb.get_lonlatalt(utc_time)
# Создаём в шейп-файле новый объект
# Определеяем геометрию
track_shape.point(lon,lat)
# и атрибуты
track_shape.record(i,utc_string,lat,lon)
# Не забываем про счётчики
i += 1
minutes += step_minutes
# Вне цикла нам осталось записать созданный шейп-файл на диск.
# Т.к. мы знаем, что координаты положений ИСЗ были получены в WGS84
# можно заодно создать файл .prj с нужным описанием
try:
# Создаем файл .prj с тем же именем, что и выходной .shp
prj = open("%s.prj" % output_shapefile.replace('.shp',''), "w")
# Создаем переменную с описанием EPSG:4326 (WGS84)
wgs84_wkt = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]'
# Записываем её в файл .prj
prj.write(wgs84_wkt)
# И закрываем его
prj.close()
# Функцией save также сохраняем и сам шейп.
track_shape.save(output_shapefile)
except:
# Вдруг нет прав на запись или вроде того...
print 'Unable to save shapefile'
return
Вот и всё, давайте проверим, как это работает, и проверим корректность! Данные о положении спутника Terra с частотой 5 минут публикуются на специальном сервисе Space Science and Engineering Data Center, с ним и сверимся. Смоделируем положения на 15 декабря 2016 года и визуализируем получившийся набор геоданных в QGIS.
create_orbital_track_shapefile_for_day(25994, date(2016,12,15), 5, '/home/silent/space/terra_15_12_2016_5min.shp')
Настроив в QGIS подписи и подложив OSM получим следующую картинку:
Найдём данные на тот же день у Space Science and Engineering Data Center:
Всё прекрасно сходится! Узнаем, где будет Landsat 8 в будущем? Например, 22 декабря 2016.
create_orbital_track_shapefile_for_day(39084, date(2016,12,22), 5, '/home/silent/space/landsat_8_22_12_2016_5min.shp')
Представляя себе полосу съемки, можно оценить охват снятой территории за 1 день.
В коде показан механизм, который несложно приспособить под собственные задачи. Таким образом можно осуществлять автоматизацию поиска и загрузки архивных данных, прогнозирование пролётов. Важно помнить, что в зависимости от особенностей аппаратуры, установленной на спутнике, соотношение между треком пролёта и снятой территорий будет сильно разниться. К примеру, Sentinel-1 оснащен радиолокатором бокового обзора, его наблюдения не надирные; полосы съемки MODIS (Terra) и ETM+ (Landsat) отличаются на порядки по степени охвата (хотя треки похожи); и так далее.
2. David A. Vallado, Paul Crawford. SGP4 Orbit Determination
Обсудить в форуме Комментариев 12Редактировать в вики
Последнее обновление: 2018-04-23 12:39
Дата создания:
Автор(ы): Эдуард Казаков
© GIS-Lab и авторы, 2002-2021. При использовании материалов сайта, ссылка на GIS-Lab и авторов обязательна. Содержание материалов - ответственность авторов. (подробнее).