Исследование Земли из Космоса, 2022, № 6, стр. 101-112

Географическая привязка данных дистанционных радиометрических измерений МТВЗА-ГЯ

И. Н. Садовский^a, *, Д. С. Сазонов^a

^a Институт космических исследований РАН
Москва, Россия

Поступила в редакцию 12.07.2022

EDN: LAUMDW
DOI: 10.31857/S0205961422060100

Аннотация

В работе представлено описание процедуры географической привязки данных радиометрических измерений микроволнового сканера/зондировщика МТВЗА-ГЯ, установленного на борту космического аппарата “Метеор-М” № 2-2. Географическая привязка выполняется на основе информации, содержащейся в файлах (в формате .hdf) с данными измерений этого инструмента (время наблюдения каждого элемента поверхности, единого для всех частотных каналов), а также на описании геометрии его сканирования, представленного в научно-периодической литературе (углы визирования, скорость сканирования и т.д.). В качестве дополнительной информации, необходимой для выполнения операции геопривязки (параметры орбиты и т.п.), используются данные, в свободном доступе присутствующие в сети Интернет. Процедура геопривязки состоит из нескольких последовательно выполняемых операций, включая: определение координат спутника в инерциальной системе координат, формирование векторов наблюдения в приборной системе координат, определение матриц прямых и обратных переходов для шести используемых координатных систем, определение координат точки наблюдения в инерциальной системе координат, их перевод в гринвичскую систему координат и т.д. В представленном виде имеется возможность гибкой настройки данной операции, учитывающей нестабильность положения носителя на орбите, вариации скорости сканирования прибора, нестабильность установки временных меток, ошибки монтажа прибора на носителе и т.д. Также описанная процедура геопривязки может быть использована в составе автоматизированных алгоритмов обработки данных измерений и процедурах поиска оптимальных значений корректирующих параметров (например, углов крена, тангажа и рыскания космического аппарата и/или прибора).

Ключевые слова: дистанционное зондирование, СВЧ-радиометрия, микроволновый сканер/зондировщик МТВЗА-ГЯ, спутник, радиометр, географическая привязка данных ДЗЗ

ВВЕДЕНИЕ

Микроволновый сканер температурно-влажностного зондирования атмосферы МТВЗА-ГЯ выступает единственным российским многофункциональным СВЧ-радиометром, функционирующим в настоящее время на орбите (Барсуков и др., 2021). По информационным характеристикам и области применения радиометр МТВЗА-ГЯ соответствует зарубежным аналогам – зондировщикам AMSU-A (англ. Advanced Microwave Sounding Unit), AMSU-B, ATMS (англ. Advanced Technology Microwave Sounder), сканерам AMSR-E (англ. Advanced Microwave Scanning Radiometer), AMSR2, и сканеру/зондировщику SSMIS (англ. Special Sensor Microwave Imager/Sounder).По данным МТВЗА-ГЯ могут быть восстановлены: паросодержание атмосферы, водозапас облаков, скорость приводного ветра, вертикальные профили температуры и влажности атмосферы и т.д. (Ермаков и др., 2021).

Прибор МТВЗА-ГЯ имеет круговую коническую систему сканирования с углом визирования 53.3° (относительно вертикальной строительной оси спутника-носителя – космического аппарата (КА) “Метеор-М” № 2-2), что обеспечивает угол встречи с Землей 65°. Сканирование производится в направлении кормы КА, направление сканирования – по часовой стрелке, рабочий сектор –105° (от –90 до +15° относительно плоскости орбиты КА), что соответствует ширине полосы обзора порядка 1500 км. Период сканирования составляет 2.5 с, за это время подспутниковая точка на поверхности Земли перемещается на 16 км, что соответствует пространственному разрешению высокочастотных каналов МТВЗА-ГЯ.

Данные измерений МТВЗА-ГЯ поступают в центры приема ФБГУ “НИЦ Планета” (http://planet.iitp.ru/) дважды в сутки и оттуда, по запросу, – в ЦКП (Центр коллективного пользования системами архивации, обработки и анализа данных спутниковых наблюдений Института космических исследований РАН (ИКИ РАН)) “ИКИ-Мониторинг”. Здесь результаты измерений проходят первичную обработку до уровня L1B (с использованием программного комплекса предварительной обработки измерений МТВЗА-ГЯ (созданного разработчиками аппаратуры)), включающую распаковку, фильтрацию, геопривязку и калибровку, после чего сохраняются в архивах и могут быть предоставлены потребителям.

Несмотря на высокий потенциал приборов серии МТВЗА и богатую историю функционирования на орбите (начиная с 2001 г., КА “Метеор-3М”, инструмент МТВЗА), практическое применение данных МТВЗА-ГЯ затруднено в силу ряда причин, часть из которых приведены в статье (Ермаков и др., 2021). Наиболее существенные из них – необходимость корректировки существующих алгоритмов географической привязки данных радиометрических измерений и их последующей калибровки.

Поиск возможных причин ошибок географической привязки, наблюдаемых при анализе данных измерений МТВЗА-ГЯ и достигающих в некоторых случаях 40–70 км, затруднен в силу следующих причин: во-первых, в литературе полностью отсутствует детальное техническое описание процедуры географической привязки данных МТВЗА-ГЯ; во-вторых, существующие программные продукты предварительной обработки данных МТВЗА-ГЯ не имеют возможности пошагового контроля этапов действующего алгоритма геопривязки.

Указанные причины побудили авторов настоящей публикации разработать свой, независимый алгоритм географической привязки данных измерений МТВЗА-ГЯ. Описание основных этапов представленного алгоритма приведено ниже. Основными требованиями при создании алгоритма были: прозрачность (детальное описание для возможности последующего сравнения с существующей реализацией), наличие широкого набора корректирующих параметров (для “тонкой” настройки, учитывающей максимальное число возмущающих факторов), а также использование в качестве входной информации только тех данных, которые содержатся в файлах данных измерений МТВЗА-ГЯ (в формате .hdf (англ. Hierarchical Data Format)).

ПОСТАНОВКА ЗАДАЧИ

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

Ранее было отмечено, что одним из основных требований к разрабатываемому алгоритму геопривязки было наличие всей необходимой информации в файлах измерений прибора. Файлы с результатами измерений МТВЗА-ГЯ представляются в формате .hdf и, помимо непосредственно результатов измерений яркостной температуры, содержат следующую информацию, необходимую для выполнения геопривязки радиометрических данных: время в формате UTC (англ. Coordinated Universal Time), соответствующее началу каждого скана МТВЗА-ГЯ, а также набор данных об элементах орбиты спутника-носителя “Метеор-М” № 2-2 в формате TLE (англ. Two-Line Element set). Следует отметить, что последняя информация в самих файлах данных зачастую отсутствует, поэтому требуется ее получение из сторонних источников. Так, TLE-данные для интересующего спутника за нужный период времени предоставляется в свободном доступе на сайте https://www.space-track.org.

Помимо указанной, выполнение геопривязки результатов измерений МТВЗА-ГЯ требует наличия информации о параметрах сканирования данного инструмента. Из официальных данных известно (Барсуков и др., 2021), что МТВЗА-ГЯ имеет круговую коническую систему сканирования с углом визирования 53.3° (относительно вертикальной строительной оси спутника-носителя КА “Метеор-М” № 2-2), направление скана – по часовой стрелке, период сканирования T = 2.5 с. Дополнительно разработчиками аппаратуры предоставлена информация, что измерения проводятся в секторе сканирования 145°, что составляет 200 пикселей, но из-за влияния элементов конструкции КА представление данных ограничено сектором 105°. Время формирования каждого пикселя составляет примерно 0.00505 с, а его угловой размер – 0.72714° (точные значений указанных величин немного отличаются от указанных, а соотношения для их расчета будут приведены ниже). Время i-го пикселя для каждого скана определяется по формуле: $t_{i}^{{{\text{UTC}}}} = {\text{БШВ}} + 0.95236 + 0.00505(i - 1),$ где БШВ – время начала скана, соответствующее бортовой шкале времени, которое фиксируется аппаратным датчиком прибора до начала рабочего сектора и записывается в файл с данными измерений. Общее количество пикселей в скане 1 ≤ i ≤ 200, однако реальная полоса обзора для КА “Метеор-М” № 2-2 –14 ≤ i ≤ 137. Описанная информация позволяет выполнить операцию геопривязки результатов измерений МТВЗА-ГЯ.

ИСПОЛЬЗУЕМЫЕ СИСТЕМЫ КООРДИНАТ

Для решения задачи геопривязки данных МТ-ВЗА-ГЯ в представленном подходе используются шесть систем координат.

1. Геоцентрическая инерциальная система координат (ИСК). Не вращающаяся прямоугольная система координат Ox^Иy^Иz^И. Начало координат O находится в центре масс Земли. Ось x^И лежит в экваториальной плоскости и направлена в истинную точку весеннего равноденствия. Ось z^И совпадает с осью вращения Земли и направлена на Северный полюс Земли. Ось y^И дополняет систему Ox^Иy^Иz^И до правой.

2. Гринвичская географическая система координат (ГСК). Прямоугольная система координат OXYZ, имеющая начало O в центре масс Земли, вращающаяся вместе с Землей. Ось Z направлена к так называемому международному условному началу (Северный полюс). Ось X лежит в плоскости среднего гринвичского меридиана и определяет положение нуль-пункта принятой системы отсчета долгот (пересечение среднего гринвичского меридиана с экватором). Ось Y дополняет систему OXYZ до правой.

3. Орбитальная скоростная система координат (ОСК). Прямоугольная система координат Ox^ОСy^ОСz^ОС, имеющая начало O, совпадающее в данный момент времени с положением космического аппарата в пространстве. Ось z^ОС совпадает по направлению с геоцентрическим радиусом (направлена от Земли). Ось x^ОС лежит в плоскости орбиты с положительным направлением, близким к направлению вектора мгновенной скорости орбитального движения КА. Ось y^ОС нормальна к плоскости орбиты и дополняет систему Ox^ОСy^ОСz^ОС до левой.

4. Система координат космического аппарата (СК КА). Прямоугольная система координат Ox^КАy^КАz^КА, реализуемая строительными осями космического аппарата. Предполагается, что СК КА совпадает с ОСК, если углы тангажа, рысканья и крена КА равны нулю.

5. Приборная система координат (ПСК). Прямоугольная система координат Ox^Пy^Пz^П, реализуемая строительными осями конкретного измерительного прибора (в частности МТВЗА-ГЯ). Предполагается, что ПСК совпадает с ОСК, если углы тангажа, рысканья и крена прибора равны нулю.

6. Система координат NEU (North-East-Up). Прямоугольная система координат Ox^NEUy^NEUz^NEU, имеющая начало O,совпадающее в данный момент времени с положением космического аппарата в пространстве. Ось z^NEU совпадает по направлению с геоцентрическим радиусом (направлена от Земли). Ось x^NEU перпендикулярна z^NEU и лежит в плоскости, образуемой осями z^NEU и Z. Ось y^NEU дополняет систему Ox^NEUy^NEUz^NEU до левой.

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

ПОСЛЕДОВАТЕЛЬНОСТЬ ВЫПОЛНЕНИЯ ГЕОГРАФИЧЕСКОЙ ПРИВЯЗКИ

В графическом виде структурная схема процедуры географической привязки данных МТВЗА-ГЯ (для одного выделенного скана), представлена на рис. 1. Серым цветом выделены блоки, в которых возможна замена параметров геопривязки для ее более точной настройки. Далее в тексте публикации приведено описание основных этапов реализации разработанной процедуры, а также используемых моделей и приближений. Следует отметить, что предложенная методика содержит общепринятые подходы к выполнению различных этапов географической привязки данных спутниковых измерений. Наиболее близким ее аналогом является подход, реализующий гепривязку данных радиометрических измерений SSM/I и SSMIS (Sapiano et al., 2010).

Рис. 1.

Общая схема реализации геопривязки.

ОПИСАНИЕ ОСНОВНЫХ ЭТАПОВ ПРОЦЕДУРЫ ГЕОПРИВЯЗКИ

Расчет времени наблюдения элемента на поверхности. В рамках предложенного подхода расчет времени наблюдения каждого элемента на поверхности выполняется на основе записанного в файле .hdf времени, соответствующего началу каждого цикла сканирования (БШВ). Таким образом, для каждого оборота МТВЗА-ГЯ (каждого скана) имеется лишь одна временная метка. Очевидно, что за время полного оборота прибора (2.5 с) КА меняет положение в пространстве и для достижения максимальной точности геопривязки следует определять положение КА для каждого измерения в отдельности. Ранее была приведена формула, предоставленная разработчиками, определяющая время измерения выделенного элемента наблюдения: $~t_{i}^{{{\text{UTC}}}} = {\text{БШВ}} + 0.95236 + 0.00505(i - 1).$ Данное соотношение требует корректировки для удовлетворения условия T = 2.5 с, а также оговоренного правила квантования: рабочий сектор 145° разбивается на 200 элементов наблюдения. Таким образом, соотношение для вычисления времени i‑го элемента скана будет иметь вид:

$t_{i}^{{{\text{UTC}}}} = {\text{БШВ}} + 0.95236 + \frac{{2.5}}{{360}}\frac{{145}}{{200 - 1}}(i - 1),$

если файлы измерений содержат 200 элементов в скане, т.е. 1 ≤ i ≤ 200; или:

$t_{i}^{{{\text{UTC}}}} = {\text{БШВ}} + 0.95236 + \frac{{2.5}}{{360}}\frac{{145}}{{200 - 1}}(i + 12),$

если файлы измерений содержат 123 элемента в скане, и 1 ≤ i ≤ 123 (в полном скане – это элементы 14 ≤ i ≤ 137). Вычисленное таким образом время выделенного измерения в скане используется в описываемой методике далее для определения координат КА “Метеор-М” № 2-2 в ИСК и последующего определения широты и долготы точки наблюдения в ГСК.

Определение координат спутника в ИСК. Информация о параметрах орбиты КА “Метеор-М” № 2-2 в формате TLE позволяет определить его координаты в ИСК для любого момента времени. Помимо служебной информации, включающей номер, международное обозначение, номер запуска КА и т.п., TLE-файлы содержат и кеплеровы элементы орбиты: наклонение орбиты, долготу восходящего узла, эксцентриситет, аргумент перицентра, среднюю аномалию.

Для определения вектора координат R и вектора мгновенной скорости V(в ИСК на момент проведения измерений $t_{i}^{{{\text{UTC}}}}$) используется модель SGP4 (англ. Simplified General Perturbations). Модель объединяет пять отдельных математических моделей (SGP, SGP4, SDP4 (англ. Simplified Deep Space Perturbations), SGP8 и SDP8), используемых для определения местоположения низкоорбитальных (SGP) объектов и взаимного расположения объектов, находящихся в глубоком космосе (SDP). Модель предполагает работу с TLE-данными и, несмотря на название “упрощенная”, позволяет определять координаты с ошибкой, не превышающей 1 км. Модель SGP4 распространяется бесплатно в виде программного кода на наиболее популярных языках программирования Python, Java, C++, Fortran и т.д. Также модель SGP4 уже интегрирована в популярные пакеты математического моделирования.

Входными параметрами расчета выступают расчетное время в формате UTC, для которого нужно определить расположение спутника ($t_{i}^{{{\text{UTC}}}}$), и TLE-информация для этого КА (актуальная на момент расчетного времени).

Определение координат вектора наблюдения в ПСК. Расчет координат вектора наблюдения в ПСК выполняется для каждого элемента наблюдения (для каждого i-го измерения или момента времени $t_{i}^{{{\text{UTC}}}}$). Если известны мгновенные значения углов θ_i и φ_i, вектор наблюдения k_i определяется следующим образом (рис. 2а):

${{{\mathbf{k}}}_{i}} = \left[ {\begin{array}{*{20}{c}} {\sin {{{{\theta }}}_{i}}\cos {{{{\varphi }}}_{i}}} \\ {\sin {{{{\theta }}}_{i}}\sin {{{{\varphi }}}_{i}}} \\ { - \cos {{{{\theta }}}_{i}}} \end{array}} \right],$

где θ_i = 53.3° для всех элементов, а азимутальный угол φ_i, должен быть рассчитан дополнительно, с учетом имеющейся информации об угловой скорости вращения МТВЗА-ГЯ и времени наблюдения конкретного элемента:

${{\varphi }_{i}} = \frac{{360}}{{2.5}}\left( {0.95236 + \frac{{2.5}}{{360}}\frac{{145}}{{200 - 1}}(i - 1)} \right).$

Рис. 2.

К вопросу об определении координат вектора наблюдения.

Однако представленное соотношение определяет поворот относительно аппаратного датчика (формирующего временную метку БШВ), расположение которого относительно плоскости орбиты неизвестно (известно лишь, что он находится вне зоны рабочего сектора). Таким образом, азимутальный угол относительно плоскости орбиты должен быть скорректирован (рис. 2, б):

${{\varphi }_{i}} = \frac{{360}}{{2.5}}\left( {0.95236 + \frac{{2.5}}{{360}}\frac{{145}}{{200 - 1}}(i - 1)} \right)~\,\, + {{\varphi }_{{{\text{кор}}}}},$

где φ_кор учитывает угловое смещение датчика формирования временных меток относительно плоскости орбиты. Как показали предварительные оценки, значение корректирующего угла составляет φ_кор = –25° (с учетом выбранных систем координат и положительного направления вращения углов). Приведенное соотношение справедливо при работе с полными наборами данных (1 ≤ i ≤ 200) и требует корректировки, если число элементов в скане составляет 123 (по аналогии с вычислением $t_{i}^{{{\text{UTC}}}}$).

Последовательность преобразования координат. Ранее было отмечено, что при определении трех систем координат (орбитальная скоростная система координат, система координат космического аппарата и приборная система координат) на первом этапе для упрощения можно исходить из приближения об отсутствии углов крена, тангажа и рыскания как самого космического аппарата, так и ошибки установки аппаратуры на нем. В данной ситуации соответствующие оси всех трех систем координат совпадают. В самом же общем случае это не так, а для перехода между этими системами координат требуется привлечение дополнительных матриц перехода. Тогда расчет координат мгновенного вектора наблюдение в ИСК K_i требует последовательного применения к вектору k_i (определенного ранее вектора наблюдения в ПСК) трех матриц преобразования координат: ${{{\mathbf{K}}}_{i}} = ~\,\,{\text{М}}_{i}^{{{\text{ИСК}} \leftarrow {\text{ОСК}}}}{\text{М}}_{i}^{{{\text{ОСК}} \leftarrow {\text{СК\;КА}}}}{\text{М}}_{i}^{{{\text{СК\;КА}} \leftarrow {\text{ПСК}}}}{{{\mathbf{k}}}_{i}}.$

Матрица перехода из СК КА в ОСК. В рамках настоящей работы сделано предположение об отсутствии углов крена, тангажа и рыскания КА “Метеор-М” № 2-2 (на основе анализа данных бортовых систем определения координат этого аппарата). Можно показать, что как средние значения этих углов, так и их вариации оказываются существенно меньше углового разрешения аппаратуры МТВЗА-ГЯ и, соответственно, ими можно пренебречь. При этом матрица перехода из ОСК в СК КА равна матрице перехода из СК КА в ОСК: $~{\text{М}}_{i}^{{{\text{ОСК}} \leftarrow {\text{СК\;КА}}}} = ~\,\,{\text{М}}_{i}^{{{\text{СК\;КА}} \leftarrow {\text{ОСК}}}}$ и является единичной. В силу приближения их постоянства во времени индекс i может быть опущен ${\text{М}}_{i}^{{{\text{ОСК}} \leftarrow {\text{СК\;КА}}}} = ~\,\,{\text{М}}_{{}}^{{{\text{ОСК}} \leftarrow {\text{СК\;КА}}}}$.

Следует отметить, что в дальнейшем для повышения точности геопривязки и при условии бесперебойного поступления качественных данных с бортовых систем определения координат КА “Метеор-М” № 2-2 (звездных датчиков) предлагаемая методика имеет возможность расчета матриц ${\text{М}}_{i}^{{{\text{ОСК}} \leftarrow {\text{СК\;КА}}}}$ для каждого измерения.

Матрица перехода из ПСК в СК КА. В рамках настоящей работы сделано предположение, что ошибки географической привязки данных измерений МТВЗА-ГЯ вызваны неточностью установки прибора на КА. Следствием этого становится смещение приборной системы координат Ox^Пy^Пz^П относительно системы координат космического аппарата Ox^КАy^КАz^КА. Данная трансформация может быть полностью описана посредством трех углов: рыскания α_р, крена α_к и тангажа α_т, а соответствующая матрица перехода от приборной системы координат к системе координат космического аппарата${\text{М}}_{i}^{{{\text{СК\;КА}} \leftarrow {\text{ПСК}}}}$ должна выполнять учет этих углов при расчете реальных углов сканирования поверхности. Поскольку ориентация МТВЗА-ГЯ относительно КА не меняется, индекс i также может быть опущен.

Положительный угол рыскания α_р соответствует повороту против часовой стрелки системы координат Ox^Пy^Пz^П вокруг оси z^П (если смотреть в направлении этой оси (рис. 3)). Положительный угол крена α_к соответствует повороту по часовой стрелке системы координат Ox^Пy^Пz^П вокруг оси x^П (если смотреть в направлении этой оси). Положительный угол тангажа α_т соответствует повороту против часовой стрелки системы координат Ox^Пy^Пz^П вокруг оси y^П (если смотреть в направлении этой оси).

Рис. 3.

Положительные направления вращения вводимых углов.

Последовательный поворот системы координат Ox^Пy^Пz^П вокруг каждой из трех осей обеспечивает ее полное совпадение с системой координат Ox^КАy^КАz^КА. Поскольку результат данной операции зависит от порядка выполнения поворотов, следует четко определить их порядок. В рамках настоящей работы принята следующая последовательность: Р-К-Т (рыскание – крен – тангаж).

Матрицы поворота вокруг каждой из трех осей имеют следующий вид:

$\begin{gathered} {{{\text{R}}}_{z}}({{{{\alpha }}}_{{\text{р}}}}) = \left[ {\begin{array}{*{20}{c}} {\cos {{{{\alpha }}}_{{\text{р}}}}}&{ - \sin {{{{\alpha }}}_{{\text{р}}}}}&0 \\ {\sin {{{{\alpha }}}_{{\text{р}}}}}&{\cos {{{{\alpha }}}_{{\text{р}}}}}&0 \\ 0&0&1 \end{array}} \right], \\ {{{\text{R}}}_{x}}({{{{\alpha }}}_{{\text{к}}}}) = \left[ {\begin{array}{*{20}{c}} 1&0&0 \\ 0&{\cos {{{{\alpha }}}_{{\text{к}}}}}&{\sin {{{{\alpha }}}_{{\text{к}}}}} \\ 0&{ - \sin {{{{\alpha }}}_{{\text{к}}}}}&{\cos {{{{\alpha }}}_{{\text{к}}}}} \end{array}} \right], \\ {{{\text{R}}}_{y}}({{{{\alpha }}}_{{\text{т}}}}) = \left[ {\begin{array}{*{20}{c}} {\cos {{{{\alpha }}}_{{\text{т}}}}}&0&{\sin {{{{\alpha }}}_{{\text{т}}}}} \\ 0&1&0 \\ { - \sin {{{{\alpha }}}_{{\text{т}}}}}&0&{\cos {{{{\alpha }}}_{{\text{т}}}}} \end{array}} \right]. \\ \end{gathered} $

Тогда матрица перехода из ПСК в СК КА будет иметь вид (с учетом оговоренной последовательности выполнения поворотов):

$\begin{gathered} {{{\text{M}}}^{{{\text{СК\;КА}} \leftarrow {\text{ПСК}}}}} = {{{\text{R}}}_{y}}({{{{\alpha }}}_{{\text{т}}}}){{{\text{R}}}_{x}}({{{{\alpha }}}_{{\text{к}}}}){{{\text{R}}}_{z}}({{{{\alpha }}}_{{\text{р}}}}) = \\ = \left[ {\begin{array}{*{20}{c}} {\cos {{{{\alpha }}}_{{\text{т}}}}}&0&{\sin {{{{\alpha }}}_{{\text{т}}}}} \\ 0&1&0 \\ { - \sin {{{{\alpha }}}_{{\text{т}}}}}&0&{\cos {{{{\alpha }}}_{{\text{т}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} 1&0&0 \\ 0&{\cos {{{{\alpha }}}_{{\text{к}}}}}&{\sin {{{{\alpha }}}_{{\text{к}}}}} \\ 0&{ - \sin {{{{\alpha }}}_{{\text{к}}}}}&{\cos {{{{\alpha }}}_{{\text{к}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\cos {{{{\alpha }}}_{{\text{р}}}}}&{ - \sin {{{{\alpha }}}_{{\text{р}}}}}&0 \\ {\sin {{{{\alpha }}}_{{\text{р}}}}}&{\cos {{{{\alpha }}}_{{\text{р}}}}}&0 \\ 0&0&1 \end{array}} \right] = \\ = \left[ {\begin{array}{*{20}{c}} {\cos {{{{\alpha }}}_{{\text{т}}}}\cos {{{{\alpha }}}_{{\text{р}}}} - \sin {{{{\alpha }}}_{{\text{р}}}}\sin {{{{\alpha }}}_{{\text{к}}}}\sin {{{{\alpha }}}_{{\text{т}}}}}&{ - \sin {{{{\alpha }}}_{{\text{р}}}}\cos {{{{\alpha }}}_{{\text{т}}}} - \sin {{{{\alpha }}}_{{\text{р}}}}\sin {{{{\alpha }}}_{{\text{к}}}}\sin {{{{\alpha }}}_{{\text{т}}}}}&{\sin {{{{\alpha }}}_{{\text{Т}}}}\cos {{{{\alpha }}}_{{\text{к}}}}} \\ {\cos {{{{\alpha }}}_{{\text{к}}}}\sin {{{{\alpha }}}_{{\text{р}}}}}&{\cos {{{{\alpha }}}_{{\text{к}}}}\cos {{{{\alpha }}}_{{\text{р}}}}}&{\sin {{{{\alpha }}}_{{\text{к}}}}} \\ { - \sin {{{{\alpha }}}_{{\text{т}}}}\cos {{{{\alpha }}}_{{\text{р}}}} - \cos {{{{\alpha }}}_{{\text{т}}}}\sin {{{{\alpha }}}_{{\text{к}}}}\sin {{{{\alpha }}}_{{\text{р}}}}}&{\sin {{{{\alpha }}}_{{\text{р}}}}\sin {{{{\alpha }}}_{{\text{т}}}} - \cos {{{{\alpha }}}_{{\text{р}}}}\sin {{{{\alpha }}}_{{\text{к}}}}\cos {{{{\alpha }}}_{{\text{т}}}}}&{\cos {{{{\alpha }}}_{{\text{к}}}}\cos {{{{\alpha }}}_{{\text{т}}}}} \end{array}} \right]. \\ \end{gathered} $

Матрица перехода из ОСК в ИСК. С учетом сделанных замечаний для завершения перехода от приборной системе координат к инерциальной необходимо определить матрицу перехода ${\text{М}}_{i}^{{{\text{ИСК}} \leftarrow {\text{ОСК}}}}.$ Для расчета элементов матрицы ${\text{М}}_{i}^{{{\text{ИСК}} \leftarrow {\text{ОСК}}}}$ необходимо определить координаты ортов координатных осей системы ОСК в инерциальной^{системе координат. Система координат ОСК постоянно меняет свое положение в ИСК (задавая
положение КА относительно “неподвижной” Земли), поэтому расчет должен выполняться
для каждого момента измерения $t_{i}^{{{\text{UTC}}}}.$ Координатные оси ОСК в ИСК определяются следующим образом (см. рис. 4):}

Рис. 4.

Определение координатных осей ОСК и СК NEU.

• вертикальная ось Z^ОС совпадает по направлению с радиус-вектором координат и имеет единичную длину: ${\mathbf{z}}_{i}^{{{\text{ОС}}}} = \,\,~{{{{{\mathbf{R}}}_{i}}} \mathord{\left/ {\vphantom {{{{{\mathbf{R}}}_{i}}} {\left| {{{{\mathbf{R}}}_{i}}} \right|}}} \right. \kern-0em} {\left| {{{{\mathbf{R}}}_{i}}} \right|}};$

• строительная ось y^ОС перпендикулярна плоскости орбиты (задаваемой векторами V_i и R_i): ${\mathbf{y}}_{i}^{{{\text{ОС}}}} = {{({{{\mathbf{V}}}_{i}} \times {{{\mathbf{R}}}_{i}})} \mathord{\left/ {\vphantom {{({{{\mathbf{V}}}_{i}} \times {{{\mathbf{R}}}_{i}})} {\left| {{{{\mathbf{V}}}_{i}} \times {{{\mathbf{R}}}_{i}}} \right|}}} \right. \kern-0em} {\left| {{{{\mathbf{V}}}_{i}} \times {{{\mathbf{R}}}_{i}}} \right|}}.$

Продольная строительная ось КА x^ОС дополнят систему Ox^ОСy^ОСz^ОС до левой: ${\mathbf{x}}_{i}^{{{\text{ОС}}}} = ~\,\,{\mathbf{z}}_{i}^{{{\text{ОС}}}} \times {\mathbf{y}}_{i}^{{{\text{ОС}}}}.$

Найденные значения направляющих косинусов для каждой из осей ОСК позволяют сформировать искомую матрицу перехода: ${\text{М}}_{i}^{{{\text{ИСК}} \leftarrow {\text{ОСК}}}} = \left[ {{\mathbf{x}}_{i}^{{{\text{ОС}}}}~{\mathbf{y}}_{i}^{{{\text{ОС}}}}~{\mathbf{z}}_{i}^{{{\text{ОС}}}}} \right].$

Расчет координат точки визирования на поверхности Земли. В качестве модели поверхности Земли выбран эллипсоид в системе координат WGS84 (англ. World Geodetic System 1984), который имеет следующие параметры: экваториальный радиус R_Э = 6378.137 км, полярный радиус R_П = 6356.755 км, коэффициент эллиптичности f = 1/298.257.

Координаты точки наблюдения на этой поверхности в ИСК задаются вектором (рис. 5):

${\mathbf{G}} = \left[ {\begin{array}{*{20}{c}} {{{{\text{R}}}_{{\text{X}}}} + d\,{{{\text{K}}}_{{\text{X}}}}} \\ {{{{\text{R}}}_{{\text{Y}}}} + d\,{{{\text{K}}}_{{\text{Y}}}}} \\ {{{{\text{R}}}_{{\text{Z}}}} + d\,{{{\text{K}}}_{{\text{Z}}}}} \end{array}} \right],$

где R = [R_X R_Y R_Z] – радиус-вектор координат КА в ИСК; K = [K_X K_Y K_Z] – вектор наблюдения в ИСК; d – наименьшее из реальных решений квадратного уравнения:

$a{{d}^{2}} + bd + c = 0,$

имеющее вид: $~d\, = \,\frac{{ - b\, \pm \,\sqrt {{{b}^{2}}\, - \,4ac} }}{{2a}},$ где $a\, = \,\frac{{{\text{K}}_{{\text{X}}}^{2}\, + \,{\text{K}}_{{\text{Y}}}^{2}}}{{R_{{\text{Э}}}^{2}}}\, + \,\frac{{{\text{K}}_{{\text{Z}}}^{2}}}{{R_{{\text{П}}}^{2}}}{\kern 1pt} ;$ $b = ~\frac{{2{{{\text{R}}}_{{\text{X}}}}{{{\text{K}}}_{{\text{X}}}} + 2{{{\text{R}}}_{{\text{Y}}}}{{{\text{K}}}_{{\text{Y}}}}}}{{R_{{\text{Э}}}^{2}}} + \frac{{2{{{\text{R}}}_{{\text{Z}}}}{{{\text{K}}}_{{\text{Z}}}}}}{{R_{{\text{П}}}^{2}}};$ $с = \frac{{{\text{R}}_{{\text{X}}}^{2} + {\text{R}}_{{\text{Y}}}^{2}}}{{R_{{\text{Э}}}^{2}}} + $ $\frac{{{\text{R}}_{{\text{Z}}}^{2}}}{{R_{{\text{П}}}^{2}}} - 1.$

Рис. 5.

К определению точки пересечения линии визирования с поверхностью эллипсоида.

Следует отметить, что в рамках данного подраздела опущены индексы i при векторах G, R и K (обозначающих их соответствие конкретному моменту измерений $t_{i}^{{{\text{UTC}}}}$), а также прочих вспомогательных переменных, для того чтобы не перегружать ими приведенные соотношения.

Знание координат вектора G позволяет определить геодезические координаты точки наблюдения в ИСК. В соответствии с работой (Patt, Gregg, 1994), геодезическая широта равна: ${{N}^{{{\text{ИСК}}}}} = \arctan \left( {{{{{{\text{G}}}_{{\text{Z}}}}} \mathord{\left/ {\vphantom {{{{{\text{G}}}_{{\text{Z}}}}} {\left( {(1 - {{f}^{2}})\sqrt {{\text{G}}_{{\text{X}}}^{2} + {\text{G}}_{{\text{Y}}}^{2}} } \right)}}} \right. \kern-0em} {\left( {(1 - {{f}^{2}})\sqrt {{\text{G}}_{{\text{X}}}^{2} + {\text{G}}_{{\text{Y}}}^{2}} } \right)}}} \right).$ Геодезическая долгота ${{E}^{{{\text{ИСК}}}}} = \arctan \left( {{{{{{\text{G}}}_{{\text{Y}}}}} \mathord{\left/ {\vphantom {{{{{\text{G}}}_{{\text{Y}}}}} {{{{\text{G}}}_{{\text{X}}}}}}} \right. \kern-0em} {{{{\text{G}}}_{{\text{X}}}}}}} \right).$

Следующим, завершающим этапом становится вычисление корректирующей добавки к найденному значению геодезической долготы, возникающей за счет вращения Земли и необходимой к учету при переходе к гринвичской системе координат. Найденное значение широты не изменяется. Таким образом: N^ГСК = N^ИСК и E^ГСК = = E^ИСК – θ_GAST, где θ_GAST – гринвичское истинное звездное время (в град) (англ. Greenwich Apparent Sidereal Time), соответствующее углу поворота ИСК вокруг оси Z (общей для ИСК и ГСК) на момент проведения измерений (рис. 6).

Рис. 6.

К вопросу перехода из ИСК в ГСК.

Определение гринвичского звездного времени. Переход из ИСК в ГСК представляет собой поворот системы ИСК вокруг оси Z (см. рис. 6) на угол θ_GAST. Следует отметить, что в самом общем случае, гринвичское истинное звездное время представляет собой сумму θ_GAST = θ_GMST + Δψ·cosε, где θ_GMST – гринвичское среднее звездное время (англ. Greenwich Mean Sidereal Time); Δψ – нутация по долготе; ε –средний наклон эклиптики к экватору.

В предлагаемом методе реализации геопривязки используются соотношения, применяемые в рамках модели SGP4 для расчета координат КА. В рамках этой модели углы θ_GAST и θ_GMST считаются равными, а эффектом коротко-периодических возмущений оси вращения Земли, нутацией, пренебрегают. Таким образом, расчет гринвичского истинного звездного времени θ_GAST упрощается и сводится к вычислению θ_GMST.

На первом шаге расчета θ_GMST необходимо скорректировать время, соответствующее моменту проведения измерений, в соответствии с соотношением $t_{i}^{{{\text{UT}}1}} = t_{i}^{{{\text{UTC}}}} + \Delta {\text{UT}}1,$ где$~t_{i}^{{{\text{UTC}}}}$ – время проведения измерений в формате UTC (в данном формате время записывается в файлах с результатами измерений МТВЗА-ГЯ); $t_{i}^{{{\text{UT}}1}}$ – время измерений в формате UT1 (современная реализация всемирного времени Universal Time), являющееся аргументом для последующих вычислений θ_GMST; ΔUT1 – поправка, учитывающая вариации вращения Земли, не поддающиеся точному прогнозированию. Приближенные значения поправок ΔUT1 можно найти, например, на сайте https://www.iers.org в разделе Earth Orientation Data.

На следующем шаге расчета θ_GMST определяется длина временного отрезка ΔD, предшествующего моменту измерения $t_{i}^{{{\text{UT}}1}}$ (здесь удобно использовать следующий формат представления времени: год-месяц-день часы:минуты:секунды. милисекунды, например, 2022-06-09 17:42:31.345), считая от стандартно используемой эпохи t₂₀₀ (2000-01-01 12:00:00.000): $\Delta D = {\text{\;}}t_{i}^{{{\text{UT}}1}} - {{t}_{{2000}}}.$ При этом результат должен быть представлен в виде действительного числа, в котором целая часть равна количеству полных дней, а дробная часть учитывает количество минут, секунд и миллисекунд найденной разницы. Далее определяется величина ΔD₁₀₀ = ΔD/36525 – значение ΔD в юлианских столетиях по 36 525 средних солнечных суток. Полученное значение используется для расчета итогового гринвичского среднего звездного времени θ_GMST:

$\begin{gathered} {{{{\theta }}}_{{{\text{GMST}}}}} = \,\,~{{{{\omega }}}_{{\text{З}}}}\left( {67310.54841 + \Delta {{D}_{{100}}} \times _{{}}^{{}}} \right. \\ \times \,\,\left( {876600 \times 3600 + {\text{\;}}8640184.812866 + _{{_{{}}^{{}}}}^{{}}} \right. \\ \left. {\left. { + \,\,\Delta {{D}_{{100}}}\left( {0.093104 - \Delta {{D}_{{100}}} \times 6.2 \times {{{10}}^{{ - 6}}}} \right)} \right)} \right), \\ \end{gathered} $

где ω_З –угловая скорость вращения Земли, составляющая ω_З = 360 град/86400 с = 1/240 град/с = = 7.2921158494 × 10^–5 град/с. Для упрощения последующих преобразований с использованием найденного значения θ_GMST следует работать не с его полной величиной, а с остатком от деления последней на 360°.

Расчет гринвичского звездного времени выполняется для каждого момента измерений $t_{i}^{{{\text{UTC}}}}$ (несмотря на то, что соответствующие индексы i при θ_GMST и θ_GAST опущены).

ОЦЕНКА КОРРЕКТНОСТИ РАБОТЫ ОПИСАННОГО ПОДХОДА

Для проверки правильности работы предложенного подхода были использованы данные измерений МТВЗА-ГЯ, полученные 16 февраля 2020 г. На рис. 7 приведен результат географической привязки данных радиометрических измерений в районе Австралии. Представленные изображения соответствуют частоте 31.5 ГГц, случай горизонтальной поляризации. Данные относятся к восходящим полувиткам – траектория движения КА “Метеор-М” № 2-2 проходит снизу вверх, направление сканирования МТВЗА-ГЯ – справа налево. Изображение рис. 7а получено для нулевых значений корректирующих значений крена, тангажа и рыскания (но с учетом поправки φ_кор = –25°). Отчетливо видно, что при общем правильном наложении радиометрического изображения на контурную карту (насыщенный красно-оранжевый цвет соответствует высокой яркостной температуре суши), наблюдается их смещение, наиболее заметное в области береговой линии (черная сплошная линия).

Рис. 7.

Смещение радиометрического портрета в результате ввода корректирующих углов. Комментарии даны в тексте.

Введение положительного угла рыскания (см. рис. 7, б) приводит, как и ожидалось, к “повороту” радиометрического изображения относительно географической береговой линии по часовой стрелке. Естественно, что траектория подспутниковой точки остается неизменной, поворачиваются лишь элементы наблюдения поверхности вокруг нее.

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

ОПРЕДЕЛЕНИЕ РЕАЛЬНЫХ УГЛОВ ВСТРЕЧИ С ЗЕМЛЕЙ

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

Для решения поставленной задачи сначала необходимо определить координаты ортов координатных осей системы NEU (стандартно применяемой для расчета азимутального угла относительно направления “на север”) в инерциальной системе координат (рис. 8). Используя соотношения работы (Patt, Gregg, 1994), при известных координатах точки наблюдения на поверхности Земли G, вектор нормали к поверхности, совпадающий с осью ${{{\mathbf{z}}}^{{{\text{NEU}}}}},$ имеет координаты:

${{{\mathbf{z}}}^{{{\text{NEU}}}}} = \left[ {\begin{array}{*{20}{c}} { - {{{\text{G}}}_{{\text{X}}}}\frac{{f{\kern 1pt} '}}{{\sqrt {{\text{G}}_{{\text{Z}}}^{2} + f{\kern 1pt} {{'}^{{~2}}}\left( {{\text{G}}_{{\text{X}}}^{2} + {\text{G}}_{{\text{Y}}}^{2}} \right)} }}} \\ { - {{{\text{G}}}_{{\text{Y}}}}\frac{{f{\kern 1pt} '}}{{\sqrt {{\text{G}}_{{\text{Z}}}^{2} + f{\kern 1pt} {{'}^{{~2}}}\left( {{\text{G}}_{{\text{X}}}^{2} + {\text{G}}_{{\text{Y}}}^{2}} \right)} }}} \\ { - {{{\text{G}}}_{{\text{Z}}}}\frac{1}{{\sqrt {{\text{G}}_{{\text{Z}}}^{2} + f{\kern 1pt} {{'}^{{~2}}}\left( {{\text{G}}_{{\text{X}}}^{2} + {\text{G}}_{{\text{Y}}}^{2}} \right)} }}} \end{array}} \right],$

где $f{\kern 1pt} ' = \,\,~{{\left( {1 - f} \right)}^{2}}$.

Рис. 8.

Определение осей координат NEU в ИСК.

Ось y^NEU определяется путем векторного перемножения оси Z ИСК и z^NEU: ${{{\mathbf{y}}}^{{{\text{NEU}}}}} = {{({\mathbf{Z}} \times {{{\mathbf{z}}}^{{{\text{NEU}}}}})} \mathord{\left/ {\vphantom {{({\mathbf{Z}} \times {{{\mathbf{z}}}^{{{\text{NEU}}}}})} {\left| {{\mathbf{Z}} \times {{{\mathbf{z}}}^{{{\text{NEU}}}}}} \right|}}} \right. \kern-0em} {\left| {{\mathbf{Z}} \times {{{\mathbf{z}}}^{{{\text{NEU}}}}}} \right|}}.$ Оставшаяся ось x^NEU дополняет систему Ox^NEUy^NEUz^NEU до левой: x^NEU = z^NEU×y^NEU.

Координаты вектора ${\mathbf{K}}{\kern 1pt} {\kern 1pt} ' = - {\mathbf{K}}$ (вектора, обратного вектору наблюдения K, определенного ранее в ИСК) в системе NEU определяются путем его скалярного перемножения на орты соответствующих осей: ${\text{K}}_{{\text{Z}}}^{'} = {\mathbf{K}}{\kern 1pt} '\,\, \cdot {{{\mathbf{z}}}^{{{\text{NEU}}}}},$ ${\text{K}}_{{\text{Y}}}^{'} = {\mathbf{K}}{\kern 1pt} '\,\, \cdot {{{\mathbf{y}}}^{{{\text{NEU}}}}}$ и ${\text{K}}_{{\text{X}}}^{'} = {\mathbf{K}}{\kern 1pt} '\,\, \cdot {{{\mathbf{x}}}^{{{\text{NEU}}}}}.$ Найденные значения позволяют определить искомые углы φ_EIA и θ_EIA (англ. Earth Incidence Angle) (рис. 9):

$\begin{gathered} {{{{\theta }}}_{{{\text{EIA}}}}} = \operatorname{arc} {\text{tg}}\left( {\frac{{\sqrt {{\text{K}}_{{\text{X}}}^{{'2}} + {\text{K}}_{{\text{Y}}}^{{'2}}} }}{{{\text{K}}_{{\text{Z}}}^{{'2}}}}} \right), \\ {{{{\varphi }}}_{{{\text{EIA}}}}} = \operatorname{arc} {\text{tg}}\left( {\frac{{{\text{K}}_{{\text{Y}}}^{'}}}{{{\text{K}}_{{\text{X}}}^{'}}}} \right). \\ \end{gathered} $

Рис. 9.

Определение углов встречи с Землей в СК NEU.

Расчет пары углов φ_EIA и θ_EIA также выполняется для каждого момента измерений $t_{i}^{{{\text{UTC}}}}$ (несмотря на то, что соответствующие индексы i при них опущены). В рамках реализации процедуры геопривязки, помимо изменения широты и долготы каждого из элементов наблюдения поверхности, в исходные файлы данных измерений МТ-ВЗА-ГЯ будет добавлена и информация об этой паре углов: φ_EIA и θ_EIA.

ЗАКЛЮЧЕНИЕ

Операция географической привязки данных спутниковых наблюдений – ключевой этап обработки, затрагивающий не только точность последующей локализации отслеживаемых процессов и явлений, но и кардинально влияющий на результаты восстановления всех параметров системы “подстилающая поверхность – атмосфера” по данным ДЗЗ за счет ошибочного учета (или игнорирования) возникающих угловых и поляризационных эффектов.

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

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

Представленный подход, как и сама публикация, должны стать основой необходимого информационного обеспечения прибора МТВЗА-ГЯ, что будет способствовать существенному росту востребованности данных измерений российского многофункционального радиометра и их научной и практической значимости, а также позволит продвинуться в направлении разработки сервисов по предоставлению готовых продуктов на основе измерений МТВЗА-ГЯ.

Список литературы

Барсуков И.А., Болдырев В.В., Гаврилов М.И., Евсеев Г.Е., Егоров А.Н., Ильгасов П.А., Панцов В.Ю., Стрельников Н.И., Стрельцов А.М., Черный И.В., Чернявский Г.М., Яковлев В.В. Спутниковая СВЧ-радиометрия для решения задач дистанционного зондирования Земли // Ракетно-косм. приборостроение и информац. системы. 2021. Т. 8. Вып. 1. С. 11–23. https://doi.org/10.30894/issn2409-0239.2021.8.1.11.23
Ермаков Д.М.,Кузьмин А.В., Мазуров А.А., Пашинов Е.В., Садовский И.Н., Сазонов Д.С., Стерлядкин В.В., Чернушич А.П., Черный И.В., Стрельцов А.М., Шарков Е.А., Екимов Н.С. Концепция потоковой обработки данных российских спутниковых СВЧ-радиометров серии МТВЗА на базе ЦКП “ИКИ-Мониторинг”// Современные проблемы дистанционного зондирования Земли из космоса. 2021. Т. 18. № 4. С. 298–303. https://doi.org/10.21046/2070-7401-2021-18-4-298-303
Patt F.S., Gregg W.W. Exact closed-form geolocation algorithm for Earth survey sensors// Intern. J. Remote Sensing. 1994. V. 15(18). P. 3719–3734. https://doi.org/10.1080/01431169408954354
Sapiano M.R.P., Bilanow S., Berg W.SSM/I and SSMIS Stewardship Code Geolocation Algorithm Theoretical Basis: CSU Technical Report. Colorado State University, 2010. 35 p. URL: http://rain.atmos.colostate.edu/FCDR/doc/CSU_FCDR_geolocation_tech_report.pdf.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Исследование Земли из Космоса

Архивы выпусков Информация о журнале Отправить рукопись в журнал